From c35ae8707e3e9be886084ec69658896d1af59d1d Mon Sep 17 00:00:00 2001 From: Brian Jimenez Date: Tue, 16 Apr 2024 18:15:43 +0200 Subject: [PATCH 1/5] Allow some Clippy lints --- Cargo.toml | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index d92298d..0d3b43f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "lightdock" -version = "0.3.0" +version = "0.3.1" authors = ["Brian Jimenez Garcia "] edition = "2021" @@ -13,3 +13,9 @@ lazy_static = "1.4.0" npyz = "0.8.3" log = "0.4.21" env_logger = "0.11.3" + +[lints.clippy] +borrowed_box = "allow" +needless_range_loop = "allow" +too_many_arguments = "allow" +new_ret_no_self = "allow" From 17f1231e33eb84115e8605e227bf251278947f9b Mon Sep 17 00:00:00 2001 From: Brian Jimenez Date: Tue, 16 Apr 2024 18:16:34 +0200 Subject: [PATCH 2/5] Use npyz without deprecated API --- src/bin/lightdock-rust.rs | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/src/bin/lightdock-rust.rs b/src/bin/lightdock-rust.rs index a9e6108..df4046f 100644 --- a/src/bin/lightdock-rust.rs +++ b/src/bin/lightdock-rust.rs @@ -13,11 +13,11 @@ use std::fs; use serde::{Serialize, Deserialize}; use std::error::Error; use std::fs::File; -use std::io::{Read, BufReader}; +use std::io::BufReader; use std::path::Path; use std::collections::HashMap; use std::thread; -use npyz::NpyData; +use npyz::NpyFile; // Use 8MB as binary stack const STACK_SIZE: usize = 8 * 1024 * 1024; @@ -161,7 +161,7 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step println!("Writing to swarm dir {:?}", swarm_directory); let positions = parse_input_coordinates(swarm_filename); - let receptor_filename = if simulation_path == "" { + let receptor_filename = if simulation_path.is_empty() { format!("{}{}", DEFAULT_LIGHTDOCK_PREFIX, setup.receptor_pdb) } else { format!("{}/{}{}", simulation_path, DEFAULT_LIGHTDOCK_PREFIX, setup.receptor_pdb) @@ -170,7 +170,7 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step println!("Reading receptor input structure: {}", receptor_filename); let (receptor, _errors) = pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Medium).unwrap(); - let ligand_filename = if simulation_path == "" { + let ligand_filename = if simulation_path.is_empty() { format!("{}{}", DEFAULT_LIGHTDOCK_PREFIX, setup.ligand_pdb) } else { format!("{}/{}{}", simulation_path, DEFAULT_LIGHTDOCK_PREFIX, setup.ligand_pdb) @@ -183,18 +183,28 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step let mut rec_nm: Vec = Vec::new(); let mut lig_nm: Vec = Vec::new(); if setup.use_anm { - let mut buf = vec![]; if setup.anm_rec > 0 { - std::fs::File::open(DEFAULT_REC_NM_FILE).unwrap().read_to_end(&mut buf).unwrap(); - rec_nm = NpyData::from_bytes(&buf).unwrap().to_vec(); + let bytes = match std::fs::read(DEFAULT_REC_NM_FILE) { + Ok(bytes) => bytes, + Err(e) => { + panic!("Error reading receptor ANM file [{:?}]: {:?}", DEFAULT_REC_NM_FILE, e.to_string()); + } + }; + let reader = NpyFile::new(&bytes[..]).unwrap(); + rec_nm = reader.into_vec::().unwrap(); if rec_nm.len() != receptor.atom_count() * 3 * setup.anm_rec { panic!("Number of read ANM in receptor does not correspond to the number of atoms"); } } if setup.anm_lig > 0 { - buf = vec![]; - std::fs::File::open(DEFAULT_LIG_NM_FILE).unwrap().read_to_end(&mut buf).unwrap(); - lig_nm = NpyData::from_bytes(&buf).unwrap().to_vec(); + let bytes = match std::fs::read(DEFAULT_LIG_NM_FILE) { + Ok(bytes) => bytes, + Err(e) => { + panic!("Error reading ligand ANM file [{:?}]: {:?}", DEFAULT_LIG_NM_FILE, e.to_string()); + } + }; + let reader = NpyFile::new(&bytes[..]).unwrap(); + lig_nm = reader.into_vec::().unwrap(); if lig_nm.len() != ligand.atom_count() * 3 * setup.anm_lig { panic!("Number of read ANM in ligand does not correspond to the number of atoms"); } From 13b882c6fde5ca2e0a19acf7d9654a8abe393425 Mon Sep 17 00:00:00 2001 From: Brian Jimenez Date: Tue, 16 Apr 2024 18:16:50 +0200 Subject: [PATCH 3/5] Removed unused lifetime --- src/dfire.rs | 2 +- src/dna.rs | 2 +- src/pydock.rs | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/dfire.rs b/src/dfire.rs index 13f12c9..4384618 100644 --- a/src/dfire.rs +++ b/src/dfire.rs @@ -229,7 +229,7 @@ impl<'a> DFIRE { } } -impl<'a> Score for DFIRE { +impl Score for DFIRE { fn energy(&self, translation: &[f64], rotation: &Quaternion, rec_nmodes: &[f64], lig_nmodes: &[f64]) -> f64 { diff --git a/src/dna.rs b/src/dna.rs index 258a31c..a710cca 100644 --- a/src/dna.rs +++ b/src/dna.rs @@ -374,7 +374,7 @@ impl<'a> DNA { } } -impl<'a> Score for DNA { +impl Score for DNA { fn energy(&self, translation: &[f64], rotation: &Quaternion, rec_nmodes: &[f64], lig_nmodes: &[f64]) -> f64 { diff --git a/src/pydock.rs b/src/pydock.rs index a29683b..cfcda53 100644 --- a/src/pydock.rs +++ b/src/pydock.rs @@ -387,7 +387,7 @@ impl<'a> PYDOCK { } } -impl<'a> Score for PYDOCK { +impl Score for PYDOCK { fn energy(&self, translation: &[f64], rotation: &Quaternion, rec_nmodes: &[f64], lig_nmodes: &[f64]) -> f64 { From 54a7f4a5fd6f8670d650e141398edbbc76d74c81 Mon Sep 17 00:00:00 2001 From: Brian Jimenez Date: Tue, 16 Apr 2024 18:43:43 +0200 Subject: [PATCH 4/5] Prepare Cargo.toml for publising crate --- Cargo.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/Cargo.toml b/Cargo.toml index 0d3b43f..4bd1fc9 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,6 +3,9 @@ name = "lightdock" version = "0.3.1" authors = ["Brian Jimenez Garcia "] edition = "2021" +license = "GPL-3.0-only" +description = "Macromolecular docking software based on the GSO algorithm" +homepage = "https://lightdock.org" [dependencies] rand = "0.7.3" From ce5bfbf425816005b908910de1ae14d77854427f Mon Sep 17 00:00:00 2001 From: Brian Jimenez Date: Wed, 17 Apr 2024 09:43:19 +0200 Subject: [PATCH 5/5] cargo fmt --- src/bin/lightdock-rust.rs | 189 ++++++++++++++++++++---------- src/constants.rs | 2 +- src/dfire.rs | 207 +++++++++++++++++++++------------ src/dna.rs | 231 +++++++++++++++++++++++-------------- src/glowworm.rs | 92 +++++++++------ src/lib.rs | 37 +++--- src/pydock.rs | 237 ++++++++++++++++++++++++-------------- src/qt.rs | 202 ++++++++++++++++++-------------- src/scoring.rs | 18 +-- src/swarm.rs | 154 +++++++++++++++---------- 10 files changed, 849 insertions(+), 520 deletions(-) diff --git a/src/bin/lightdock-rust.rs b/src/bin/lightdock-rust.rs index df4046f..9903e46 100644 --- a/src/bin/lightdock-rust.rs +++ b/src/bin/lightdock-rust.rs @@ -1,45 +1,47 @@ +extern crate npyz; extern crate serde; extern crate serde_json; -extern crate npyz; -use lightdock::GSO; -use lightdock::constants::{DEFAULT_LIGHTDOCK_PREFIX, DEFAULT_SEED, DEFAULT_REC_NM_FILE, DEFAULT_LIG_NM_FILE}; -use lightdock::scoring::{Score, Method}; +use lightdock::constants::{ + DEFAULT_LIGHTDOCK_PREFIX, DEFAULT_LIG_NM_FILE, DEFAULT_REC_NM_FILE, DEFAULT_SEED, +}; use lightdock::dfire::DFIRE; use lightdock::dna::DNA; use lightdock::pydock::PYDOCK; +use lightdock::scoring::{Method, Score}; +use lightdock::GSO; +use npyz::NpyFile; +use serde::{Deserialize, Serialize}; +use std::collections::HashMap; use std::env; -use std::fs; -use serde::{Serialize, Deserialize}; use std::error::Error; +use std::fs; use std::fs::File; use std::io::BufReader; use std::path::Path; -use std::collections::HashMap; use std::thread; -use npyz::NpyFile; // Use 8MB as binary stack const STACK_SIZE: usize = 8 * 1024 * 1024; #[derive(Serialize, Deserialize, Debug)] struct SetupFile { - seed: Option, - anm_seed: u64, - ftdock_file: Option, - noh: bool, - anm_rec: usize, - anm_lig: usize, - swarms: u32, - starting_points_seed: u32, - verbose_parser: bool, - noxt: bool, + seed: Option, + anm_seed: u64, + ftdock_file: Option, + noh: bool, + anm_rec: usize, + anm_lig: usize, + swarms: u32, + starting_points_seed: u32, + verbose_parser: bool, + noxt: bool, now: bool, - restraints: Option, - use_anm: bool, - glowworms: u32, - membrane: bool, - receptor_pdb: String, + restraints: Option, + use_anm: bool, + glowworms: u32, + membrane: bool, + receptor_pdb: String, ligand_pdb: String, receptor_restraints: Option>>, ligand_restraints: Option>>, @@ -57,15 +59,14 @@ fn read_setup_from_file>(path: P) -> Result Vec> { // Parse swarm filename content - let contents = fs::read_to_string(swarm_filename) - .expect("Error reading the input file"); + let contents = fs::read_to_string(swarm_filename).expect("Error reading the input file"); let mut positions: Vec> = Vec::new(); for s in contents.lines() { let vector_raw: String = String::from(s); let vector: Vec<&str> = vector_raw.split(' ').collect(); let mut position: Vec = Vec::new(); - for pos in vector.iter() { + for pos in vector.iter() { position.push(pos.trim().parse::().unwrap()); } positions.push(position); @@ -95,13 +96,11 @@ fn run() { let num_steps = &args[3]; // parse the number let steps: u32 = match num_steps.parse() { - Ok(n) => { - n - }, + Ok(n) => n, Err(_) => { eprintln!("Error: steps argument must be a number"); return; - }, + } }; let method_type = &args[4].to_lowercase(); // parse the type @@ -112,7 +111,7 @@ fn run() { _ => { eprintln!("Error: method not supported"); return; - }, + } }; // Load setup @@ -121,10 +120,19 @@ fn run() { // Simulation path let simulation_path = Path::new(setup_filename).parent().unwrap(); - simulate(simulation_path.to_str().unwrap(), &setup, swarm_filename, steps, method); + simulate( + simulation_path.to_str().unwrap(), + &setup, + swarm_filename, + steps, + method, + ); } _ => { - println!("Wrong command line. Usage: {} setup_filename swarm_filename steps method", args[0]); + println!( + "Wrong command line. Usage: {} setup_filename swarm_filename steps method", + args[0] + ); } } } @@ -137,15 +145,16 @@ fn parse_swarm_id(path: &Path) -> Option { .and_then(|s| s.parse::().ok()) } -fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, steps: u32, method: Method) { - - let seed:u64 = match setup.seed { - Some(seed) => { - seed - }, - None => { - DEFAULT_SEED - }, +fn simulate( + simulation_path: &str, + setup: &SetupFile, + swarm_filename: &str, + steps: u32, + method: Method, +) { + let seed: u64 = match setup.seed { + Some(seed) => seed, + None => DEFAULT_SEED, }; println!("Reading starting positions from {:?}", swarm_filename); @@ -154,7 +163,10 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step println!("Swarm ID {:?}", swarm_id); let swarm_directory = format!("swarm_{}", swarm_id); - if !fs::metadata(&swarm_directory).map(|m| m.is_dir()).unwrap_or(false) { + if !fs::metadata(&swarm_directory) + .map(|m| m.is_dir()) + .unwrap_or(false) + { panic!("Output directory does not exist for swarm {:?}", swarm_id); } @@ -164,20 +176,28 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step let receptor_filename = if simulation_path.is_empty() { format!("{}{}", DEFAULT_LIGHTDOCK_PREFIX, setup.receptor_pdb) } else { - format!("{}/{}{}", simulation_path, DEFAULT_LIGHTDOCK_PREFIX, setup.receptor_pdb) + format!( + "{}/{}{}", + simulation_path, DEFAULT_LIGHTDOCK_PREFIX, setup.receptor_pdb + ) }; // Parse receptor input PDB structure println!("Reading receptor input structure: {}", receptor_filename); - let (receptor, _errors) = pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Medium).unwrap(); + let (receptor, _errors) = + pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Medium).unwrap(); let ligand_filename = if simulation_path.is_empty() { format!("{}{}", DEFAULT_LIGHTDOCK_PREFIX, setup.ligand_pdb) } else { - format!("{}/{}{}", simulation_path, DEFAULT_LIGHTDOCK_PREFIX, setup.ligand_pdb) + format!( + "{}/{}{}", + simulation_path, DEFAULT_LIGHTDOCK_PREFIX, setup.ligand_pdb + ) }; // Parse ligand input PDB structure println!("Reading ligand input structure: {}", ligand_filename); - let (ligand, _errors) = pdbtbx::open(&ligand_filename, pdbtbx::StrictnessLevel::Medium).unwrap(); + let (ligand, _errors) = + pdbtbx::open(&ligand_filename, pdbtbx::StrictnessLevel::Medium).unwrap(); // Read ANM data if activated let mut rec_nm: Vec = Vec::new(); @@ -187,7 +207,11 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step let bytes = match std::fs::read(DEFAULT_REC_NM_FILE) { Ok(bytes) => bytes, Err(e) => { - panic!("Error reading receptor ANM file [{:?}]: {:?}", DEFAULT_REC_NM_FILE, e.to_string()); + panic!( + "Error reading receptor ANM file [{:?}]: {:?}", + DEFAULT_REC_NM_FILE, + e.to_string() + ); } }; let reader = NpyFile::new(&bytes[..]).unwrap(); @@ -200,7 +224,11 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step let bytes = match std::fs::read(DEFAULT_LIG_NM_FILE) { Ok(bytes) => bytes, Err(e) => { - panic!("Error reading ligand ANM file [{:?}]: {:?}", DEFAULT_LIG_NM_FILE, e.to_string()); + panic!( + "Error reading ligand ANM file [{:?}]: {:?}", + DEFAULT_LIG_NM_FILE, + e.to_string() + ); } }; let reader = NpyFile::new(&bytes[..]).unwrap(); @@ -213,31 +241,64 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step // Restraints let rec_active_restraints: Vec = match &setup.receptor_restraints { - Some(restraints) => { restraints["active"].clone() }, - None => { Vec::new() }, + Some(restraints) => restraints["active"].clone(), + None => Vec::new(), }; let rec_passive_restraints: Vec = match &setup.receptor_restraints { - Some(restraints) => { restraints["passive"].clone() }, - None => { Vec::new() }, + Some(restraints) => restraints["passive"].clone(), + None => Vec::new(), }; let lig_active_restraints: Vec = match &setup.ligand_restraints { - Some(restraints) => { restraints["active"].clone() }, - None => { Vec::new() }, + Some(restraints) => restraints["active"].clone(), + None => Vec::new(), }; let lig_passive_restraints: Vec = match &setup.ligand_restraints { - Some(restraints) => { restraints["passive"].clone() }, - None => { Vec::new() }, + Some(restraints) => restraints["passive"].clone(), + None => Vec::new(), }; // Scoring function println!("Loading {:?} scoring function", method); let scoring = match method { - Method::DFIRE => DFIRE::new(receptor, rec_active_restraints, rec_passive_restraints, rec_nm, setup.anm_rec, - ligand, lig_active_restraints, lig_passive_restraints, lig_nm, setup.anm_lig, setup.use_anm) as Box, - Method::DNA => DNA::new(receptor, rec_active_restraints, rec_passive_restraints, rec_nm, setup.anm_rec, - ligand, lig_active_restraints, lig_passive_restraints, lig_nm, setup.anm_lig, setup.use_anm) as Box, - Method::PYDOCK => PYDOCK::new(receptor, rec_active_restraints, rec_passive_restraints, rec_nm, setup.anm_rec, - ligand, lig_active_restraints, lig_passive_restraints, lig_nm, setup.anm_lig, setup.use_anm) as Box, + Method::DFIRE => DFIRE::new( + receptor, + rec_active_restraints, + rec_passive_restraints, + rec_nm, + setup.anm_rec, + ligand, + lig_active_restraints, + lig_passive_restraints, + lig_nm, + setup.anm_lig, + setup.use_anm, + ) as Box, + Method::DNA => DNA::new( + receptor, + rec_active_restraints, + rec_passive_restraints, + rec_nm, + setup.anm_rec, + ligand, + lig_active_restraints, + lig_passive_restraints, + lig_nm, + setup.anm_lig, + setup.use_anm, + ) as Box, + Method::PYDOCK => PYDOCK::new( + receptor, + rec_active_restraints, + rec_passive_restraints, + rec_nm, + setup.anm_rec, + ligand, + lig_active_restraints, + lig_passive_restraints, + lig_nm, + setup.anm_lig, + setup.use_anm, + ) as Box, }; // Glowworm Swarm Optimization algorithm @@ -249,7 +310,7 @@ fn simulate(simulation_path: &str, setup: &SetupFile, swarm_filename: &str, step setup.use_anm, setup.anm_rec, setup.anm_lig, - swarm_directory + swarm_directory, ); // Simulate for the given steps diff --git a/src/constants.rs b/src/constants.rs index 3c4cf4a..c1a2b1e 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -25,4 +25,4 @@ pub const DEFAULT_NMODES_STEP: f64 = 0.5; // 1D NumPy arrays containing calculated ANM from ProDy pub const DEFAULT_REC_NM_FILE: &str = "rec_nm.npy"; -pub const DEFAULT_LIG_NM_FILE: &str = "lig_nm.npy"; \ No newline at end of file +pub const DEFAULT_LIG_NM_FILE: &str = "lig_nm.npy"; diff --git a/src/dfire.rs b/src/dfire.rs index 4384618..f0e1d07 100644 --- a/src/dfire.rs +++ b/src/dfire.rs @@ -1,11 +1,11 @@ +use super::constants::{INTERFACE_CUTOFF, MEMBRANE_PENALTY_SCORE}; +use super::qt::Quaternion; +use super::scoring::{membrane_intersection, satisfied_restraints, Score}; +use pdbtbx::PDB; use std::collections::HashMap; +use std::env; use std::fs::File; use std::io::Read; -use std::env; -use pdbtbx::PDB; -use super::qt::Quaternion; -use super::constants::{INTERFACE_CUTOFF, MEMBRANE_PENALTY_SCORE}; -use super::scoring::{Score, satisfied_restraints, membrane_intersection}; macro_rules! hashmap { ($( $key: expr => $val: expr ),*) => {{ @@ -17,36 +17,40 @@ macro_rules! hashmap { pub fn r3_to_numerical(residue_name: &str) -> usize { match residue_name { - "ALA" => {0} - "CYS" => {1} - "ASP" => {2} - "GLU" => {3} - "PHE" => {4} - "GLY" => {5} - "HIS" => {6} - "ILE" => {7} - "LYS" => {8} - "LEU" => {9} - "MET" => {10} - "ASN" => {11} - "PRO" => {12} - "GLN" => {13} - "ARG" => {14} - "SER" => {15} - "THR" => {16} - "VAL" => {17} - "TRP" => {18} - "TYR" => {19} - "MMB" => {20} - "MMY" => {0} - _ => { panic!("Residue name not supported in DFIRE scoring function") } + "ALA" => 0, + "CYS" => 1, + "ASP" => 2, + "GLU" => 3, + "PHE" => 4, + "GLY" => 5, + "HIS" => 6, + "ILE" => 7, + "LYS" => 8, + "LEU" => 9, + "MET" => 10, + "ASN" => 11, + "PRO" => 12, + "GLN" => 13, + "ARG" => 14, + "SER" => 15, + "THR" => 16, + "VAL" => 17, + "TRP" => 18, + "TYR" => 19, + "MMB" => 20, + "MMY" => 0, + _ => { + panic!("Residue name not supported in DFIRE scoring function") + } } } // DFIRE only uses 20 distance bins -const DIST_TO_BINS: &[usize] = &[1, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 14, 15, 15, - 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, - 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 32]; +const DIST_TO_BINS: &[usize] = &[ + 1, 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, + 19, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, + 32, +]; lazy_static! { static ref ATOMNUMBER: HashMap<&'static str, usize> = hashmap![ @@ -97,7 +101,6 @@ lazy_static! { vec![74, 75, 76, 77, 78, 0, 0, 0, 0, 0, 0, 0, 0, 0]]; } - pub struct DFIREDockingModel { pub atoms: Vec, pub coordinates: Vec<[f64; 3]>, @@ -109,9 +112,13 @@ pub struct DFIREDockingModel { } impl<'a> DFIREDockingModel { - - fn new(structure: &'a PDB, active_restraints: &'a [String], passive_restraints: &'a [String], - nmodes: &[f64], num_anm: usize) -> DFIREDockingModel { + fn new( + structure: &'a PDB, + active_restraints: &'a [String], + passive_restraints: &'a [String], + nmodes: &[f64], + num_anm: usize, + ) -> DFIREDockingModel { let mut model = DFIREDockingModel { atoms: Vec::new(), coordinates: Vec::new(), @@ -145,10 +152,12 @@ impl<'a> DFIREDockingModel { match model.active_restraints.get_mut(&res_id) { Some(atom_indexes) => { atom_indexes.push(atom_index as usize); - }, + } None => { - model.active_restraints.insert(res_id.to_string(), vec![atom_index as usize]); - }, + model + .active_restraints + .insert(res_id.to_string(), vec![atom_index as usize]); + } } } @@ -156,10 +165,12 @@ impl<'a> DFIREDockingModel { match model.passive_restraints.get_mut(&res_id) { Some(atom_indexes) => { atom_indexes.push(atom_index as usize); - }, + } None => { - model.passive_restraints.insert(res_id.to_string(), vec![atom_index as usize]); - }, + model + .passive_restraints + .insert(res_id.to_string(), vec![atom_index as usize]); + } } } @@ -180,23 +191,42 @@ impl<'a> DFIREDockingModel { } pub struct DFIRE { - pub potential: Vec, + pub potential: Vec, pub receptor: DFIREDockingModel, pub ligand: DFIREDockingModel, pub use_anm: bool, } - impl<'a> DFIRE { - - pub fn new(receptor: PDB, rec_active_restraints: Vec, rec_passive_restraints: Vec, - rec_nmodes: Vec, rec_num_anm: usize, - ligand: PDB, lig_active_restraints: Vec, lig_passive_restraints: Vec, - lig_nmodes: Vec, lig_num_anm: usize, use_anm: bool) -> Box { + pub fn new( + receptor: PDB, + rec_active_restraints: Vec, + rec_passive_restraints: Vec, + rec_nmodes: Vec, + rec_num_anm: usize, + ligand: PDB, + lig_active_restraints: Vec, + lig_passive_restraints: Vec, + lig_nmodes: Vec, + lig_num_anm: usize, + use_anm: bool, + ) -> Box { let mut d = DFIRE { potential: Vec::with_capacity(169 * 169 * 20), - receptor: DFIREDockingModel::new(&receptor, &rec_active_restraints, &rec_passive_restraints, &rec_nmodes, rec_num_anm), - ligand: DFIREDockingModel::new(&ligand, &lig_active_restraints, &lig_passive_restraints, &lig_nmodes, lig_num_anm), + receptor: DFIREDockingModel::new( + &receptor, + &rec_active_restraints, + &rec_passive_restraints, + &rec_nmodes, + rec_num_anm, + ), + ligand: DFIREDockingModel::new( + &ligand, + &lig_active_restraints, + &lig_passive_restraints, + &lig_nmodes, + lig_num_anm, + ), use_anm, }; d.load_potentials(); @@ -213,13 +243,15 @@ impl<'a> DFIRE { let parameters_path: String = format!("{}/DCparams", data_folder); - File::open(parameters_path).expect("Unable to open DFIRE parameters") - .read_to_string(&mut raw_parameters).expect("Unable to read DFIRE parameters"); + File::open(parameters_path) + .expect("Unable to open DFIRE parameters") + .read_to_string(&mut raw_parameters) + .expect("Unable to read DFIRE parameters"); let split = raw_parameters.lines(); let params: Vec<&str> = split.collect(); - for param in params.iter().take(169*169*20) { + for param in params.iter().take(169 * 169 * 20) { self.potential.push(param.trim().parse::().unwrap()); } } @@ -230,9 +262,13 @@ impl<'a> DFIRE { } impl Score for DFIRE { - - fn energy(&self, translation: &[f64], rotation: &Quaternion, - rec_nmodes: &[f64], lig_nmodes: &[f64]) -> f64 { + fn energy( + &self, + translation: &[f64], + rotation: &Quaternion, + rec_nmodes: &[f64], + lig_nmodes: &[f64], + ) -> f64 { let mut score: f64 = 0.0; // Clone receptor coordinates @@ -255,9 +291,12 @@ impl Score for DFIRE { for i_nm in 0usize..self.ligand.num_anm { // (num_anm, num_atoms, 3) -> 1d // Endianness: i = i_nm * num_atoms * 3 + i_atom * 3 + coord - coordinate[0] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3] * lig_nmodes[i_nm]; - coordinate[1] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 1] * lig_nmodes[i_nm]; - coordinate[2] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 2] * lig_nmodes[i_nm]; + coordinate[0] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3] + * lig_nmodes[i_nm]; + coordinate[1] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 1] + * lig_nmodes[i_nm]; + coordinate[2] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 2] + * lig_nmodes[i_nm]; } } } @@ -268,9 +307,14 @@ impl Score for DFIRE { for i_nm in 0usize..self.receptor.num_anm { // (num_anm, num_atoms, 3) -> 1d // Endianness: i = i_nm * num_atoms * 3 + i_atom * 3 + coord - coordinate[0] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3] * rec_nmodes[i_nm]; - coordinate[1] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3 + 1] * rec_nmodes[i_nm]; - coordinate[2] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3 + 2] * rec_nmodes[i_nm]; + coordinate[0] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3] + * rec_nmodes[i_nm]; + coordinate[1] += self.receptor.nmodes + [i_nm * rec_num_atoms * 3 + i_atom * 3 + 1] + * rec_nmodes[i_nm]; + coordinate[2] += self.receptor.nmodes + [i_nm * rec_num_atoms * 3 + i_atom * 3 + 2] + * rec_nmodes[i_nm]; } } } @@ -284,12 +328,14 @@ impl Score for DFIRE { let z1 = ra[2]; let atoma = self.receptor.atoms[i]; for (j, la) in ligand_coordinates.iter().enumerate() { - let dist = (x1-la[0])*(x1-la[0]) + (y1-la[1])*(y1-la[1]) + (z1-la[2])*(z1-la[2]); + let dist = (x1 - la[0]) * (x1 - la[0]) + + (y1 - la[1]) * (y1 - la[1]) + + (z1 - la[2]) * (z1 - la[2]); if dist <= 225. { let atomb = self.ligand.atoms[j]; - let d = dist.sqrt()*2.0 - 1.0; + let d = dist.sqrt() * 2.0 - 1.0; let dfire_bin = DIST_TO_BINS[d as usize] - 1; - score += self.potential[atoma*169*20 + atomb*20 + dfire_bin]; + score += self.potential[atoma * 169 * 20 + atomb * 20 + dfire_bin]; if d <= INTERFACE_CUTOFF { interface_receptor[i] = 1; interface_ligand[j] = 1; @@ -301,10 +347,10 @@ impl Score for DFIRE { score = (score * 0.0157 - 4.7) * -1.0; // Bias the scoring depending on satisfied restraints - let perc_receptor_restraints: f64 = satisfied_restraints(&interface_receptor, - &self.receptor.active_restraints); - let perc_ligand_restraints: f64 = satisfied_restraints(&interface_ligand, - &self.ligand.active_restraints); + let perc_receptor_restraints: f64 = + satisfied_restraints(&interface_receptor, &self.receptor.active_restraints); + let perc_ligand_restraints: f64 = + satisfied_restraints(&interface_ligand, &self.ligand.active_restraints); // Take into account membrane intersection let mut membrane_penalty: f64 = 0.0; let intersection = membrane_intersection(&interface_receptor, &self.receptor.membrane); @@ -342,13 +388,26 @@ mod tests { let test_path: String = format!("{}/tests/2oob", cargo_path); let receptor_filename: String = format!("{}/2oob_receptor.pdb", test_path); - let (receptor, _errors) = pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); + let (receptor, _errors) = + pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); let ligand_filename: String = format!("{}/2oob_ligand.pdb", test_path); - let (ligand, _errors) = pdbtbx::open(&ligand_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); - - let scoring = DFIRE::new(receptor, Vec::new(), Vec::new(), Vec::new(), 0, - ligand, Vec::new(), Vec::new(), Vec::new(), 0, false); + let (ligand, _errors) = + pdbtbx::open(&ligand_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); + + let scoring = DFIRE::new( + receptor, + Vec::new(), + Vec::new(), + Vec::new(), + 0, + ligand, + Vec::new(), + Vec::new(), + Vec::new(), + 0, + false, + ); let translation = vec![0., 0., 0.]; let rotation = Quaternion::default(); diff --git a/src/dna.rs b/src/dna.rs index a710cca..2e63b2c 100644 --- a/src/dna.rs +++ b/src/dna.rs @@ -1,8 +1,8 @@ -use std::collections::HashMap; -use pdbtbx::PDB; -use super::qt::Quaternion; use super::constants::{INTERFACE_CUTOFF2, MEMBRANE_PENALTY_SCORE}; -use super::scoring::{Score, satisfied_restraints, membrane_intersection}; +use super::qt::Quaternion; +use super::scoring::{membrane_intersection, satisfied_restraints, Score}; +use pdbtbx::PDB; +use std::collections::HashMap; macro_rules! hashmap { ($( $key: expr => $val: expr ),*) => {{ @@ -18,36 +18,46 @@ const MAX_ES_CUTOFF: f64 = 1.0; const MIN_ES_CUTOFF: f64 = -1.0; const VDW_CUTOFF: f64 = 1.0; const ELEC_DIST_CUTOFF: f64 = 30.0; -const ELEC_DIST_CUTOFF2: f64 = ELEC_DIST_CUTOFF*ELEC_DIST_CUTOFF; +const ELEC_DIST_CUTOFF2: f64 = ELEC_DIST_CUTOFF * ELEC_DIST_CUTOFF; const VDW_DIST_CUTOFF: f64 = 10.0; -const VDW_DIST_CUTOFF2: f64 = VDW_DIST_CUTOFF*VDW_DIST_CUTOFF; -const ELEC_MAX_CUTOFF: f64 = MAX_ES_CUTOFF*EPSILON/FACTOR; -const ELEC_MIN_CUTOFF: f64 = MIN_ES_CUTOFF*EPSILON/FACTOR; +const VDW_DIST_CUTOFF2: f64 = VDW_DIST_CUTOFF * VDW_DIST_CUTOFF; +const ELEC_MAX_CUTOFF: f64 = MAX_ES_CUTOFF * EPSILON / FACTOR; +const ELEC_MIN_CUTOFF: f64 = MIN_ES_CUTOFF * EPSILON / FACTOR; pub fn atoms_in_residues(residue_name: &str) -> &'static [&'static str] { - match residue_name { - "ALA" => { &["N", "CA", "C", "O", "CB"] } - "CYS" => { &["N", "CA", "C", "O", "CB", "SG"] } - "ASP" => { &["N", "CA", "C", "O", "CB", "CG", "OD1", "OD2"] } - "GLU" => { &["N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "OE2"] } - "PHE" => { &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE1", "CE2", "CZ"] } - "GLY" => { &["N", "CA", "C", "O"] } - "HIS" => { &["N", "CA", "C", "O", "CB", "CG", "ND1", "CD2", "CE1", "NE2"] } - "ILE" => { &["N", "CA", "C", "O", "CB", "CG1", "CG2", "CD1"] } - "LYS" => { &["N", "CA", "C", "O", "CB", "CG", "CD", "CE", "NZ"] } - "LEU" => { &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2"] } - "MET" => { &["N", "CA", "C", "O", "CB", "CG", "SD", "CE"] } - "ASN" => { &["N", "CA", "C", "O", "CB", "CG", "OD1", "ND2"] } - "PRO" => { &["N", "CA", "C", "O", "CB", "CG", "CD"] } - "GLN" => { &["N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "NE2"] } - "ARG" => { &["N", "CA", "C", "O", "CB", "CG", "CD", "NE", "CZ", "NH1", "NH2"] } - "SER" => { &["N", "CA", "C", "O", "CB", "OG"] } - "THR" => { &["N", "CA", "C", "O", "CB", "OG1", "CG2"] } - "VAL" => { &["N", "CA", "C", "O", "CB", "CG1", "CG2"] } - "TRP" => { &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE2", "NE1", "CE3", "CZ3", "CH2", "CZ2"] } - "TYR" => { &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE1", "CE2", "CZ", "OH"] } - "MMB" => { &["BJ"] } - _ => { panic!("Residue name not supported in DNA scoring function") } + match residue_name { + "ALA" => &["N", "CA", "C", "O", "CB"], + "CYS" => &["N", "CA", "C", "O", "CB", "SG"], + "ASP" => &["N", "CA", "C", "O", "CB", "CG", "OD1", "OD2"], + "GLU" => &["N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "OE2"], + "PHE" => &[ + "N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE1", "CE2", "CZ", + ], + "GLY" => &["N", "CA", "C", "O"], + "HIS" => &["N", "CA", "C", "O", "CB", "CG", "ND1", "CD2", "CE1", "NE2"], + "ILE" => &["N", "CA", "C", "O", "CB", "CG1", "CG2", "CD1"], + "LYS" => &["N", "CA", "C", "O", "CB", "CG", "CD", "CE", "NZ"], + "LEU" => &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2"], + "MET" => &["N", "CA", "C", "O", "CB", "CG", "SD", "CE"], + "ASN" => &["N", "CA", "C", "O", "CB", "CG", "OD1", "ND2"], + "PRO" => &["N", "CA", "C", "O", "CB", "CG", "CD"], + "GLN" => &["N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "NE2"], + "ARG" => &[ + "N", "CA", "C", "O", "CB", "CG", "CD", "NE", "CZ", "NH1", "NH2", + ], + "SER" => &["N", "CA", "C", "O", "CB", "OG"], + "THR" => &["N", "CA", "C", "O", "CB", "OG1", "CG2"], + "VAL" => &["N", "CA", "C", "O", "CB", "CG1", "CG2"], + "TRP" => &[ + "N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE2", "NE1", "CE3", "CZ3", "CH2", "CZ2", + ], + "TYR" => &[ + "N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE1", "CE2", "CZ", "OH", + ], + "MMB" => &["BJ"], + _ => { + panic!("Residue name not supported in DNA scoring function") + } } } @@ -61,7 +71,6 @@ lazy_static! { "P" => 0.2, "S" => 0.25, "CR" => 0.086, "N2" => 0.17, "N3" => 0.17, "CW" => 0.086, "CV" => 0.086, "CT" => 0.1094, "MG" => 0.8947, "OH" => 0.2104, "H2" => 0.0157, "H3" => 0.0157, "H1" => 0.0157, "H4" => 0.015, "H5" => 0.015, "SH" => 0.25, "OW" => 0.152, "OS" => 0.17]; - static ref VDW_RADII: HashMap<&'static str, f64> = hashmap![ "IP" => 1.868, "HS" => 0.6, "HP" => 1.1, "Na" => 1.868, "N*" => 1.824, "Li" => 1.137, "HO" => 0.0001, "Rb" => 2.956, "HC" => 1.487, "HA" => 1.459, "O3" => 1.6612, "CQ" => 1.908, "C*" => 1.908, @@ -71,10 +80,8 @@ lazy_static! { "Zn" => 1.1, "O" => 1.6612, "N" => 1.824, "P" => 2.1, "S" => 2.0, "CR" => 1.908, "N2" => 1.824, "N3" => 1.875, "CW" => 1.908, "CV" => 1.908, "CT" => 1.908, "MG" => 0.7926, "OH" => 1.721, "H2" => 1.287, "H3" => 1.187, "H1" => 1.387, "H4" => 1.409, "H5" => 1.359, "SH" => 2.0, "OW" => 1.7683, "OS" => 1.6837]; - static ref RES_TO_TRANSLATE: HashMap<&'static str, &'static str> = hashmap![ "HIS" => "HID", "THY" => "DT", "ADE" => "DA", "CYT" => "DC", "GUA" => "DG"]; - static ref AMBER_TYPES: HashMap<&'static str, &'static str> = hashmap![ "ALA-C" => "C", "ALA-CA" => "CT", "ALA-CB" => "CT", "ALA-H" => "H", "ALA-HA" => "H1", "ALA-HB1" => "HC", "ALA-HB2" => "HC", "ALA-HB3" => "HC", "ALA-N" => "N", "ALA-O" => "O", "ARG-C" => "C", "ARG-CA" => "CT", "ARG-CB" => "CT", "ARG-CD" => "CT", "ARG-CG" => "CT", "ARG-CZ" => "CA", "ARG-H" => "H", "ARG-HA" => "H1", "ARG-HB2" => "HC", "ARG-HB3" => "HC", "ARG-HD2" => "H1", "ARG-HD3" => "H1", "ARG-HE" => "H", "ARG-HG2" => "HC", "ARG-HG3" => "HC", "ARG-HH11" => "H", "ARG-HH12" => "H", "ARG-HH21" => "H", "ARG-HH22" => "H", "ARG-N" => "N", "ARG-NE" => "N2", "ARG-NH1" => "N2", "ARG-NH2" => "N2", "ARG-O" => "O", @@ -136,7 +143,6 @@ lazy_static! { "TRP-C" => "C", "TRP-CA" => "CT", "TRP-CB" => "CT", "TRP-CD1" => "CW", "TRP-CD2" => "CB", "TRP-CE2" => "CN", "TRP-CE3" => "CA", "TRP-CG" => "C*", "TRP-CH2" => "CA", "TRP-CZ2" => "CA", "TRP-CZ3" => "CA", "TRP-H" => "H", "TRP-HA" => "H1", "TRP-HB2" => "HC", "TRP-HB3" => "HC", "TRP-HD1" => "H4", "TRP-HE1" => "H", "TRP-HE3" => "HA", "TRP-HH2" => "HA", "TRP-HZ2" => "HA", "TRP-HZ3" => "HA", "TRP-N" => "N", "TRP-NE1" => "NA", "TRP-O" => "O", "TYR-C" => "C", "TYR-CA" => "CT", "TYR-CB" => "CT", "TYR-CD1" => "CA", "TYR-CD2" => "CA", "TYR-CE1" => "CA", "TYR-CE2" => "CA", "TYR-CG" => "CA", "TYR-CZ" => "C", "TYR-H" => "H", "TYR-HA" => "H1", "TYR-HB2" => "HC", "TYR-HB3" => "HC", "TYR-HD1" => "HA", "TYR-HD2" => "HA", "TYR-HE1" => "HA", "TYR-HE2" => "HA", "TYR-HH" => "HO", "TYR-N" => "N", "TYR-O" => "O", "TYR-OH" => "OH", "VAL-C" => "C", "VAL-CA" => "CT", "VAL-CB" => "CT", "VAL-CG1" => "CT", "VAL-CG2" => "CT", "VAL-H" => "H", "VAL-HA" => "H1", "VAL-HB" => "HC", "VAL-HG11" => "HC", "VAL-HG12" => "HC", "VAL-HG13" => "HC", "VAL-HG21" => "HC", "VAL-HG22" => "HC", "VAL-HG23" => "HC", "VAL-N" => "N", "VAL-O" => "O"]; - static ref ELE_CHARGES: HashMap<&'static str, f64> = hashmap![ "ALA-C" => 0.5973, "ALA-CA" => 0.0337, "ALA-CB" => -0.1825, "ALA-H" => 0.2719, "ALA-HA" => 0.0823, "ALA-HB1" => 0.0603, "ALA-HB2" => 0.0603, "ALA-HB3" => 0.0603, "ALA-N" => -0.4157, "ALA-O" => -0.5679, "ARG-C" => 0.7341, "ARG-CA" => -0.2637, "ARG-CB" => -0.0007, "ARG-CD" => 0.0486, "ARG-CG" => 0.039, "ARG-CZ" => 0.8076, "ARG-H" => 0.2747, "ARG-HA" => 0.156, "ARG-HB2" => 0.0327, "ARG-HB3" => 0.0327, "ARG-HD2" => 0.0687, "ARG-HD3" => 0.0687, "ARG-HE" => 0.3456, "ARG-HG2" => 0.0285, "ARG-HG3" => 0.0285, "ARG-HH11" => 0.4478, "ARG-HH12" => 0.4478, "ARG-HH21" => 0.4478, "ARG-HH22" => 0.4478, "ARG-N" => -0.3479, "ARG-NE" => -0.5295, "ARG-NH1" => -0.8627, "ARG-NH2" => -0.8627, "ARG-O" => -0.5894, @@ -198,7 +204,6 @@ lazy_static! { "TRP-C" => 0.5973, "TRP-CA" => -0.0275, "TRP-CB" => -0.005, "TRP-CD1" => -0.1638, "TRP-CD2" => 0.1243, "TRP-CE2" => 0.138, "TRP-CE3" => -0.2387, "TRP-CG" => -0.1415, "TRP-CH2" => -0.1134, "TRP-CZ2" => -0.2601, "TRP-CZ3" => -0.1972, "TRP-H" => 0.2719, "TRP-HA" => 0.1123, "TRP-HB2" => 0.0339, "TRP-HB3" => 0.0339, "TRP-HD1" => 0.2062, "TRP-HE1" => 0.3412, "TRP-HE3" => 0.17, "TRP-HH2" => 0.1417, "TRP-HZ2" => 0.1572, "TRP-HZ3" => 0.1447, "TRP-N" => -0.4157, "TRP-NE1" => -0.3418, "TRP-O" => -0.5679, "TYR-C" => 0.5973, "TYR-CA" => -0.0014, "TYR-CB" => -0.0152, "TYR-CD1" => -0.1906, "TYR-CD2" => -0.1906, "TYR-CE1" => -0.2341, "TYR-CE2" => -0.2341, "TYR-CG" => -0.0011, "TYR-CZ" => 0.3226, "TYR-H" => 0.2719, "TYR-HA" => 0.0876, "TYR-HB2" => 0.0295, "TYR-HB3" => 0.0295, "TYR-HD1" => 0.1699, "TYR-HD2" => 0.1699, "TYR-HE1" => 0.1656, "TYR-HE2" => 0.1656, "TYR-HH" => 0.3992, "TYR-N" => -0.4157, "TYR-O" => -0.5679, "TYR-OH" => -0.5579, "VAL-C" => 0.5973, "VAL-CA" => -0.0875, "VAL-CB" => 0.2985, "VAL-CG1" => -0.3192, "VAL-CG2" => -0.3192, "VAL-H" => 0.2719, "VAL-HA" => 0.0969, "VAL-HB" => -0.0297, "VAL-HG11" => 0.0791, "VAL-HG12" => 0.0791, "VAL-HG13" => 0.0791, "VAL-HG21" => 0.0791, "VAL-HG22" => 0.0791, "VAL-HG23" => 0.0791, "VAL-N" => -0.4157, "VAL-O" => -0.5679]; - static ref NT_ELE_CHARGES: HashMap<&'static str, f64> = hashmap![ "ACE-C" => 0.5972, "ACE-CH3" => -0.3662, "ACE-HH31" => 0.1123, "ACE-HH32" => 0.1123, "ACE-HH33" => 0.1123, "ACE-O" => -0.5679, "ALA-C" => 0.6163, "ALA-CA" => 0.0962, "ALA-CB" => -0.0597, "ALA-H1" => 0.1997, "ALA-H2" => 0.1997, "ALA-H3" => 0.1997, "ALA-HA" => 0.0889, "ALA-HB1" => 0.03, "ALA-HB2" => 0.03, "ALA-HB3" => 0.03, "ALA-N" => 0.1414, "ALA-O" => -0.5722, @@ -241,9 +246,13 @@ pub struct DNADockingModel { } impl<'a> DNADockingModel { - - fn new(structure: &'a PDB, active_restraints: &'a [String], passive_restraints: &'a [String], - nmodes: &[f64], num_anm: usize) -> DNADockingModel { + fn new( + structure: &'a PDB, + active_restraints: &'a [String], + passive_restraints: &'a [String], + nmodes: &[f64], + num_anm: usize, + ) -> DNADockingModel { let mut model = DNADockingModel { atoms: Vec::new(), coordinates: Vec::new(), @@ -280,10 +289,12 @@ impl<'a> DNADockingModel { match model.active_restraints.get_mut(&res_id) { Some(atom_indexes) => { atom_indexes.push(atom_index as usize); - }, + } None => { - model.active_restraints.insert(res_id.to_string(), vec![atom_index as usize]); - }, + model + .active_restraints + .insert(res_id.to_string(), vec![atom_index as usize]); + } } } @@ -291,10 +302,12 @@ impl<'a> DNADockingModel { match model.passive_restraints.get_mut(&res_id) { Some(atom_indexes) => { atom_indexes.push(atom_index as usize); - }, + } None => { - model.passive_restraints.insert(res_id.to_string(), vec![atom_index as usize]); - }, + model + .passive_restraints + .insert(res_id.to_string(), vec![atom_index as usize]); + } } } @@ -314,18 +327,19 @@ impl<'a> DNADockingModel { } else { panic!("DNA Error: Atom [{:?}] not supported", atom_id); } - }, + } }; // Assign electrostatics charge let ele_charge = match ELE_CHARGES.get(&*atom_id) { Some(&charge) => charge, - _ => { - match NT_ELE_CHARGES.get(&*atom_id) { - Some(&charge) => charge, - _ => panic!("DNA Error: Atom [{:?}] electrostatics charge not found", atom_id), - } - } + _ => match NT_ELE_CHARGES.get(&*atom_id) { + Some(&charge) => charge, + _ => panic!( + "DNA Error: Atom [{:?}] electrostatics charge not found", + atom_id + ), + }, }; model.ele_charges.push(ele_charge); @@ -351,23 +365,42 @@ impl<'a> DNADockingModel { } pub struct DNA { - pub potential: Vec, + pub potential: Vec, pub receptor: DNADockingModel, pub ligand: DNADockingModel, pub use_anm: bool, } - impl<'a> DNA { - - pub fn new(receptor: PDB, rec_active_restraints: Vec, rec_passive_restraints: Vec, - rec_nmodes: Vec, rec_num_anm: usize, - ligand: PDB, lig_active_restraints: Vec, lig_passive_restraints: Vec, - lig_nmodes: Vec, lig_num_anm: usize, use_anm: bool) -> Box { + pub fn new( + receptor: PDB, + rec_active_restraints: Vec, + rec_passive_restraints: Vec, + rec_nmodes: Vec, + rec_num_anm: usize, + ligand: PDB, + lig_active_restraints: Vec, + lig_passive_restraints: Vec, + lig_nmodes: Vec, + lig_num_anm: usize, + use_anm: bool, + ) -> Box { let d = DNA { potential: Vec::with_capacity(168 * 168 * 20), - receptor: DNADockingModel::new(&receptor, &rec_active_restraints, &rec_passive_restraints, &rec_nmodes, rec_num_anm), - ligand: DNADockingModel::new(&ligand, &lig_active_restraints, &lig_passive_restraints, &lig_nmodes, lig_num_anm), + receptor: DNADockingModel::new( + &receptor, + &rec_active_restraints, + &rec_passive_restraints, + &rec_nmodes, + rec_num_anm, + ), + ligand: DNADockingModel::new( + &ligand, + &lig_active_restraints, + &lig_passive_restraints, + &lig_nmodes, + lig_num_anm, + ), use_anm, }; Box::new(d) @@ -375,10 +408,13 @@ impl<'a> DNA { } impl Score for DNA { - - fn energy(&self, translation: &[f64], rotation: &Quaternion, - rec_nmodes: &[f64], lig_nmodes: &[f64]) -> f64 { - + fn energy( + &self, + translation: &[f64], + rotation: &Quaternion, + rec_nmodes: &[f64], + lig_nmodes: &[f64], + ) -> f64 { // Clone receptor coordinates let mut receptor_coordinates: Vec<[f64; 3]> = self.receptor.coordinates.clone(); let rec_num_atoms = receptor_coordinates.len(); @@ -399,9 +435,12 @@ impl Score for DNA { for i_nm in 0usize..self.ligand.num_anm { // (num_anm, num_atoms, 3) -> 1d // Endianness: i = i_nm * num_atoms * 3 + i_atom * 3 + coord - coordinate[0] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3] * lig_nmodes[i_nm]; - coordinate[1] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 1] * lig_nmodes[i_nm]; - coordinate[2] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 2] * lig_nmodes[i_nm]; + coordinate[0] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3] + * lig_nmodes[i_nm]; + coordinate[1] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 1] + * lig_nmodes[i_nm]; + coordinate[2] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 2] + * lig_nmodes[i_nm]; } } } @@ -412,9 +451,14 @@ impl Score for DNA { for i_nm in 0usize..self.receptor.num_anm { // (num_anm, num_atoms, 3) -> 1d // Endianness: i = i_nm * num_atoms * 3 + i_atom * 3 + coord - coordinate[0] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3] * rec_nmodes[i_nm]; - coordinate[1] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3 + 1] * rec_nmodes[i_nm]; - coordinate[2] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3 + 2] * rec_nmodes[i_nm]; + coordinate[0] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3] + * rec_nmodes[i_nm]; + coordinate[1] += self.receptor.nmodes + [i_nm * rec_num_atoms * 3 + i_atom * 3 + 1] + * rec_nmodes[i_nm]; + coordinate[2] += self.receptor.nmodes + [i_nm * rec_num_atoms * 3 + i_atom * 3 + 2] + * rec_nmodes[i_nm]; } } } @@ -429,11 +473,14 @@ impl Score for DNA { let y1 = ra[1]; let z1 = ra[2]; for (j, la) in ligand_coordinates.iter().enumerate() { - let distance2 = (x1-la[0])*(x1-la[0]) + (y1-la[1])*(y1-la[1]) + (z1-la[2])*(z1-la[2]); + let distance2 = (x1 - la[0]) * (x1 - la[0]) + + (y1 - la[1]) * (y1 - la[1]) + + (z1 - la[2]) * (z1 - la[2]); // Electrostatics energy if distance2 <= ELEC_DIST_CUTOFF2 { - let mut atom_elec = self.receptor.ele_charges[i] * self.ligand.ele_charges[j] / distance2; + let mut atom_elec = + self.receptor.ele_charges[i] * self.ligand.ele_charges[j] / distance2; if atom_elec > ELEC_MAX_CUTOFF { atom_elec = ELEC_MAX_CUTOFF; } @@ -445,10 +492,11 @@ impl Score for DNA { // Van der Waals energy if distance2 <= VDW_DIST_CUTOFF2 { - let vdw_energy = (self.receptor.vdw_charges[i] * self.ligand.vdw_charges[j]).sqrt(); + let vdw_energy = + (self.receptor.vdw_charges[i] * self.ligand.vdw_charges[j]).sqrt(); let vdw_radius = self.receptor.vdw_radii[i] + self.ligand.vdw_radii[j]; let p6 = vdw_radius.powi(6) / distance2.powi(3); - let mut k = vdw_energy * (p6*p6 - 2.0 * p6); + let mut k = vdw_energy * (p6 * p6 - 2.0 * p6); if k > VDW_CUTOFF { k = VDW_CUTOFF; } @@ -466,10 +514,10 @@ impl Score for DNA { let score = (total_elec + total_vdw) * -1.0; // Bias the scoring depending on satisfied restraints - let perc_receptor_restraints: f64 = satisfied_restraints(&interface_receptor, - &self.receptor.active_restraints); - let perc_ligand_restraints: f64 = satisfied_restraints(&interface_ligand, - &self.ligand.active_restraints); + let perc_receptor_restraints: f64 = + satisfied_restraints(&interface_receptor, &self.receptor.active_restraints); + let perc_ligand_restraints: f64 = + satisfied_restraints(&interface_ligand, &self.ligand.active_restraints); // Take into account membrane intersection let mut membrane_penalty: f64 = 0.0; let intersection = membrane_intersection(&interface_receptor, &self.receptor.membrane); @@ -484,8 +532,8 @@ impl Score for DNA { #[cfg(test)] mod tests { use super::*; - use std::env; use crate::qt::Quaternion; + use std::env; #[test] fn test_1azp() { @@ -496,12 +544,26 @@ mod tests { let test_path: String = format!("{}/tests/1azp", cargo_path); let receptor_filename: String = format!("{}/1azp_receptor.pdb", test_path); - let (receptor, _errors) = pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); + let (receptor, _errors) = + pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); let ligand_filename: String = format!("{}/1azp_ligand.pdb", test_path); - let (ligand, _errors) = pdbtbx::open(&ligand_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); - - let scoring = DNA::new(receptor, Vec::new(), Vec::new(), Vec::new(), 0, ligand, Vec::new(), Vec::new(), Vec::new(), 0, false); + let (ligand, _errors) = + pdbtbx::open(&ligand_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); + + let scoring = DNA::new( + receptor, + Vec::new(), + Vec::new(), + Vec::new(), + 0, + ligand, + Vec::new(), + Vec::new(), + Vec::new(), + 0, + false, + ); let translation = vec![0., 0., 0.]; let rotation = Quaternion::default(); @@ -509,4 +571,3 @@ mod tests { assert_eq!(energy, -364.88126358158974); } } - diff --git a/src/glowworm.rs b/src/glowworm.rs index 2f0cbd8..d944ee8 100644 --- a/src/glowworm.rs +++ b/src/glowworm.rs @@ -1,7 +1,7 @@ -use std::f64; -use super::scoring::Score; +use super::constants::{DEFAULT_NMODES_STEP, DEFAULT_ROTATION_STEP, DEFAULT_TRANSLATION_STEP}; use super::qt::Quaternion; -use super::constants::{DEFAULT_TRANSLATION_STEP, DEFAULT_ROTATION_STEP, DEFAULT_NMODES_STEP}; +use super::scoring::Score; +use std::f64; pub struct Glowworm<'a> { pub id: u32, @@ -26,9 +26,15 @@ pub struct Glowworm<'a> { } impl<'a> Glowworm<'a> { - pub fn new(id: u32, translation:Vec, rotation:Quaternion, - rec_nmodes: Vec, lig_nmodes: Vec, scoring_function: &'a Box, - use_anm: bool) -> Self { + pub fn new( + id: u32, + translation: Vec, + rotation: Quaternion, + rec_nmodes: Vec, + lig_nmodes: Vec, + scoring_function: &'a Box, + use_anm: bool, + ) -> Self { Glowworm { id, translation, @@ -54,8 +60,12 @@ impl<'a> Glowworm<'a> { pub fn compute_luciferin(&mut self) { if self.moved || self.step == 0 { - self.scoring = self.scoring_function.energy(&self.translation, &self.rotation, - &self.rec_nmodes, &self.lig_nmodes); + self.scoring = self.scoring_function.energy( + &self.translation, + &self.rotation, + &self.rec_nmodes, + &self.lig_nmodes, + ); } self.luciferin = (1.0 - self.rho) * self.luciferin + self.gamma * self.scoring; self.step += 1; @@ -68,7 +78,7 @@ impl<'a> Glowworm<'a> { let y2 = other.translation[1]; let z1 = self.translation[2]; let z2 = other.translation[2]; - ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)).sqrt() + ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2)).sqrt() } pub fn is_neighbor(&mut self, other: &Glowworm) -> bool { @@ -79,9 +89,10 @@ impl<'a> Glowworm<'a> { } pub fn update_vision_range(&mut self) { - self.vision_range = (self.max_vision_range).min( - (0_f64).max(self.vision_range + self.beta * f64::from(self.max_neighbors as i32 - - (self.neighbors.len() as i32)))); + self.vision_range = (self.max_vision_range).min((0_f64).max( + self.vision_range + + self.beta * f64::from(self.max_neighbors as i32 - (self.neighbors.len() as i32)), + )); } pub fn compute_probability_moving_toward_neighbor(&mut self, luciferins: &[f64]) { @@ -95,72 +106,81 @@ impl<'a> Glowworm<'a> { total_sum += difference; } - for i in 0..self.neighbors.len(){ + for i in 0..self.neighbors.len() { self.probabilities[i] /= total_sum; } } - - pub fn select_random_neighbor(&mut self, random_number:f64) -> u32 { + pub fn select_random_neighbor(&mut self, random_number: f64) -> u32 { if self.neighbors.is_empty() { return self.id; } - let mut sum_probabilities:f64 = 0.0; - let mut i:usize = 0; + let mut sum_probabilities: f64 = 0.0; + let mut i: usize = 0; while sum_probabilities < random_number { sum_probabilities += self.probabilities[i]; i += 1; } - self.neighbors[i-1] + self.neighbors[i - 1] } - pub fn move_towards(&mut self, other_id:u32, other_position:&[f64], - other_rotation:&Quaternion, other_anm_rec:&[f64], other_anm_lig:&[f64]) { + pub fn move_towards( + &mut self, + other_id: u32, + other_position: &[f64], + other_rotation: &Quaternion, + other_anm_rec: &[f64], + other_anm_lig: &[f64], + ) { self.moved = self.id != other_id; if self.id != other_id { // Translation component - let mut delta_x:Vec = vec![other_position[0] - self.translation[0], - other_position[1] - self.translation[1], - other_position[2] - self.translation[2]]; - let norm:f64 = (delta_x[0]*delta_x[0] + delta_x[1]*delta_x[1] + delta_x[2]*delta_x[2]).sqrt(); - let coef:f64 = DEFAULT_TRANSLATION_STEP / norm; + let mut delta_x: Vec = vec![ + other_position[0] - self.translation[0], + other_position[1] - self.translation[1], + other_position[2] - self.translation[2], + ]; + let norm: f64 = + (delta_x[0] * delta_x[0] + delta_x[1] * delta_x[1] + delta_x[2] * delta_x[2]) + .sqrt(); + let coef: f64 = DEFAULT_TRANSLATION_STEP / norm; delta_x[0] *= coef; delta_x[1] *= coef; delta_x[2] *= coef; self.translation[0] += delta_x[0]; self.translation[1] += delta_x[1]; self.translation[2] += delta_x[2]; - + // Rotation component self.rotation = self.rotation.slerp(other_rotation, DEFAULT_ROTATION_STEP); - + // ANM component if self.use_anm && !self.rec_nmodes.is_empty() { - let mut delta_anm:Vec = Vec::new(); + let mut delta_anm: Vec = Vec::new(); let mut cum_norm: f64 = 0.0; for i in 0..self.rec_nmodes.len() { let diff = other_anm_rec[i] - self.rec_nmodes[i]; delta_anm.push(diff); - cum_norm += diff*diff + cum_norm += diff * diff } - let anm_rec_norm:f64 = cum_norm.sqrt(); - let anm_rec_coef:f64 = DEFAULT_NMODES_STEP / anm_rec_norm; + let anm_rec_norm: f64 = cum_norm.sqrt(); + let anm_rec_coef: f64 = DEFAULT_NMODES_STEP / anm_rec_norm; for i in 0..self.rec_nmodes.len() { delta_anm[i] *= anm_rec_coef; self.rec_nmodes[i] += delta_anm[i]; } } if self.use_anm && !self.lig_nmodes.is_empty() { - let mut delta_anm:Vec = Vec::new(); + let mut delta_anm: Vec = Vec::new(); let mut cum_norm: f64 = 0.0; for i in 0..self.lig_nmodes.len() { let diff = other_anm_lig[i] - self.lig_nmodes[i]; delta_anm.push(diff); - cum_norm += diff*diff + cum_norm += diff * diff } - let anm_lig_norm:f64 = cum_norm.sqrt(); - let anm_lig_coef:f64 = DEFAULT_NMODES_STEP / anm_lig_norm; + let anm_lig_norm: f64 = cum_norm.sqrt(); + let anm_lig_coef: f64 = DEFAULT_NMODES_STEP / anm_lig_norm; for i in 0..self.lig_nmodes.len() { delta_anm[i] *= anm_lig_coef; self.lig_nmodes[i] += delta_anm[i]; @@ -178,5 +198,5 @@ pub fn distance(one: &Glowworm, two: &Glowworm) -> f64 { let y2 = two.translation[1]; let z1 = one.translation[2]; let z2 = two.translation[2]; - ((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2)).sqrt() + ((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2)).sqrt() } diff --git a/src/lib.rs b/src/lib.rs index 3315664..fdcdab2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,45 +1,50 @@ #[macro_use] - extern crate lazy_static; extern crate rand; -pub mod glowworm; -pub mod swarm; -pub mod qt; pub mod constants; -pub mod scoring; pub mod dfire; pub mod dna; +pub mod glowworm; pub mod pydock; +pub mod qt; +pub mod scoring; +pub mod swarm; -use swarm::Swarm; -use rand::SeedableRng; +use log::info; use rand::rngs::StdRng; +use rand::SeedableRng; use scoring::Score; -use log::info; - +use swarm::Swarm; pub struct GSO<'a> { pub swarm: Swarm<'a>, pub rng: StdRng, - pub output_directory: String + pub output_directory: String, } - impl<'a> GSO<'a> { - pub fn new(positions: &[Vec], seed: u64, scoring: &'a Box, use_anm: bool, - rec_num_anm: usize, lig_num_anm: usize, output_directory: String) -> Self { + pub fn new( + positions: &[Vec], + seed: u64, + scoring: &'a Box, + use_anm: bool, + rec_num_anm: usize, + lig_num_anm: usize, + output_directory: String, + ) -> Self { let mut gso = GSO { swarm: Swarm::new(), rng: SeedableRng::seed_from_u64(seed), - output_directory + output_directory, }; - gso.swarm.add_glowworms(positions, scoring, use_anm, rec_num_anm, lig_num_anm); + gso.swarm + .add_glowworms(positions, scoring, use_anm, rec_num_anm, lig_num_anm); gso } pub fn run(&mut self, steps: u32) { - for step in 1..steps+1 { + for step in 1..steps + 1 { info!("Step {}", step); self.swarm.update_luciferin(); self.swarm.movement_phase(&mut self.rng); diff --git a/src/pydock.rs b/src/pydock.rs index cfcda53..7c77822 100644 --- a/src/pydock.rs +++ b/src/pydock.rs @@ -1,11 +1,10 @@ -use std::collections::HashMap; -use pdbtbx::PDB; -use super::qt::Quaternion; use super::constants::{INTERFACE_CUTOFF2, MEMBRANE_PENALTY_SCORE}; -use super::scoring::{Score, satisfied_restraints, membrane_intersection}; - -use log::{warn, info}; +use super::qt::Quaternion; +use super::scoring::{membrane_intersection, satisfied_restraints, Score}; +use pdbtbx::PDB; +use std::collections::HashMap; +use log::{info, warn}; macro_rules! hashmap { ($( $key: expr => $val: expr ),*) => {{ @@ -21,36 +20,46 @@ const MAX_ES_CUTOFF: f64 = 1.0; const MIN_ES_CUTOFF: f64 = -1.0; const VDW_CUTOFF: f64 = 1.0; const ELEC_DIST_CUTOFF: f64 = 30.0; -const ELEC_DIST_CUTOFF2: f64 = ELEC_DIST_CUTOFF*ELEC_DIST_CUTOFF; +const ELEC_DIST_CUTOFF2: f64 = ELEC_DIST_CUTOFF * ELEC_DIST_CUTOFF; const VDW_DIST_CUTOFF: f64 = 10.0; -const VDW_DIST_CUTOFF2: f64 = VDW_DIST_CUTOFF*VDW_DIST_CUTOFF; -const ELEC_MAX_CUTOFF: f64 = MAX_ES_CUTOFF*EPSILON/FACTOR; -const ELEC_MIN_CUTOFF: f64 = MIN_ES_CUTOFF*EPSILON/FACTOR; +const VDW_DIST_CUTOFF2: f64 = VDW_DIST_CUTOFF * VDW_DIST_CUTOFF; +const ELEC_MAX_CUTOFF: f64 = MAX_ES_CUTOFF * EPSILON / FACTOR; +const ELEC_MIN_CUTOFF: f64 = MIN_ES_CUTOFF * EPSILON / FACTOR; pub fn atoms_in_residues(residue_name: &str) -> &'static [&'static str] { - match residue_name { - "ALA" => { &["N", "CA", "C", "O", "CB"] } - "CYS" => { &["N", "CA", "C", "O", "CB", "SG"] } - "ASP" => { &["N", "CA", "C", "O", "CB", "CG", "OD1", "OD2"] } - "GLU" => { &["N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "OE2"] } - "PHE" => { &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE1", "CE2", "CZ"] } - "GLY" => { &["N", "CA", "C", "O"] } - "HIS" => { &["N", "CA", "C", "O", "CB", "CG", "ND1", "CD2", "CE1", "NE2"] } - "ILE" => { &["N", "CA", "C", "O", "CB", "CG1", "CG2", "CD1"] } - "LYS" => { &["N", "CA", "C", "O", "CB", "CG", "CD", "CE", "NZ"] } - "LEU" => { &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2"] } - "MET" => { &["N", "CA", "C", "O", "CB", "CG", "SD", "CE"] } - "ASN" => { &["N", "CA", "C", "O", "CB", "CG", "OD1", "ND2"] } - "PRO" => { &["N", "CA", "C", "O", "CB", "CG", "CD"] } - "GLN" => { &["N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "NE2"] } - "ARG" => { &["N", "CA", "C", "O", "CB", "CG", "CD", "NE", "CZ", "NH1", "NH2"] } - "SER" => { &["N", "CA", "C", "O", "CB", "OG"] } - "THR" => { &["N", "CA", "C", "O", "CB", "OG1", "CG2"] } - "VAL" => { &["N", "CA", "C", "O", "CB", "CG1", "CG2"] } - "TRP" => { &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE2", "NE1", "CE3", "CZ3", "CH2", "CZ2"] } - "TYR" => { &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE1", "CE2", "CZ", "OH"] } - "MMB" => { &["BJ"] } - _ => { panic!("Residue name not supported in PYDOCK scoring function") } + match residue_name { + "ALA" => &["N", "CA", "C", "O", "CB"], + "CYS" => &["N", "CA", "C", "O", "CB", "SG"], + "ASP" => &["N", "CA", "C", "O", "CB", "CG", "OD1", "OD2"], + "GLU" => &["N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "OE2"], + "PHE" => &[ + "N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE1", "CE2", "CZ", + ], + "GLY" => &["N", "CA", "C", "O"], + "HIS" => &["N", "CA", "C", "O", "CB", "CG", "ND1", "CD2", "CE1", "NE2"], + "ILE" => &["N", "CA", "C", "O", "CB", "CG1", "CG2", "CD1"], + "LYS" => &["N", "CA", "C", "O", "CB", "CG", "CD", "CE", "NZ"], + "LEU" => &["N", "CA", "C", "O", "CB", "CG", "CD1", "CD2"], + "MET" => &["N", "CA", "C", "O", "CB", "CG", "SD", "CE"], + "ASN" => &["N", "CA", "C", "O", "CB", "CG", "OD1", "ND2"], + "PRO" => &["N", "CA", "C", "O", "CB", "CG", "CD"], + "GLN" => &["N", "CA", "C", "O", "CB", "CG", "CD", "OE1", "NE2"], + "ARG" => &[ + "N", "CA", "C", "O", "CB", "CG", "CD", "NE", "CZ", "NH1", "NH2", + ], + "SER" => &["N", "CA", "C", "O", "CB", "OG"], + "THR" => &["N", "CA", "C", "O", "CB", "OG1", "CG2"], + "VAL" => &["N", "CA", "C", "O", "CB", "CG1", "CG2"], + "TRP" => &[ + "N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE2", "NE1", "CE3", "CZ3", "CH2", "CZ2", + ], + "TYR" => &[ + "N", "CA", "C", "O", "CB", "CG", "CD1", "CD2", "CE1", "CE2", "CZ", "OH", + ], + "MMB" => &["BJ"], + _ => { + panic!("Residue name not supported in PYDOCK scoring function") + } } } @@ -64,7 +73,6 @@ lazy_static! { "P" => 0.2, "S" => 0.25, "CR" => 0.086, "N2" => 0.17, "N3" => 0.17, "CW" => 0.086, "CV" => 0.086, "CT" => 0.1094, "MG" => 0.8947, "OH" => 0.2104, "H2" => 0.0157, "H3" => 0.0157, "H1" => 0.0157, "H4" => 0.015, "H5" => 0.015, "SH" => 0.25, "OW" => 0.152, "OS" => 0.17]; - static ref VDW_RADII: HashMap<&'static str, f64> = hashmap![ "IP" => 1.868, "HS" => 0.6, "HP" => 1.1, "Na" => 1.868, "N*" => 1.824, "Li" => 1.137, "HO" => 0.0001, "Rb" => 2.956, "HC" => 1.487, "HA" => 1.459, "O3" => 1.6612, "CQ" => 1.908, "C*" => 1.908, @@ -74,10 +82,8 @@ lazy_static! { "Zn" => 1.1, "O" => 1.6612, "N" => 1.824, "P" => 2.1, "S" => 2.0, "CR" => 1.908, "N2" => 1.824, "N3" => 1.875, "CW" => 1.908, "CV" => 1.908, "CT" => 1.908, "MG" => 0.7926, "OH" => 1.721, "H2" => 1.287, "H3" => 1.187, "H1" => 1.387, "H4" => 1.409, "H5" => 1.359, "SH" => 2.0, "OW" => 1.7683, "OS" => 1.6837]; - static ref RES_TO_TRANSLATE: HashMap<&'static str, &'static str> = hashmap![ "HIS" => "HID", "THY" => "DT", "ADE" => "DA", "CYT" => "DC", "GUA" => "DG"]; - static ref AMBER_TYPES: HashMap<&'static str, &'static str> = hashmap![ "ALA-C" => "C", "ALA-CA" => "CT", "ALA-CB" => "CT", "ALA-H" => "H", "ALA-HA" => "H1", "ALA-HB1" => "HC", "ALA-HB2" => "HC", "ALA-HB3" => "HC", "ALA-N" => "N", "ALA-O" => "O", "ARG-C" => "C", "ARG-CA" => "CT", "ARG-CB" => "CT", "ARG-CD" => "CT", "ARG-CG" => "CT", "ARG-CZ" => "CA", "ARG-H" => "H", "ARG-HA" => "H1", "ARG-HB2" => "HC", "ARG-HB3" => "HC", "ARG-HD2" => "H1", "ARG-HD3" => "H1", "ARG-HE" => "H", "ARG-HG2" => "HC", "ARG-HG3" => "HC", "ARG-HH11" => "H", "ARG-HH12" => "H", "ARG-HH21" => "H", "ARG-HH22" => "H", "ARG-N" => "N", "ARG-NE" => "N2", "ARG-NH1" => "N2", "ARG-NH2" => "N2", "ARG-O" => "O", @@ -140,7 +146,6 @@ lazy_static! { "TYR-C" => "C", "TYR-CA" => "CT", "TYR-CB" => "CT", "TYR-CD1" => "CA", "TYR-CD2" => "CA", "TYR-CE1" => "CA", "TYR-CE2" => "CA", "TYR-CG" => "CA", "TYR-CZ" => "C", "TYR-H" => "H", "TYR-HA" => "H1", "TYR-HB2" => "HC", "TYR-HB3" => "HC", "TYR-HD1" => "HA", "TYR-HD2" => "HA", "TYR-HE1" => "HA", "TYR-HE2" => "HA", "TYR-HH" => "HO", "TYR-N" => "N", "TYR-O" => "O", "TYR-OH" => "OH", "VAL-C" => "C", "VAL-CA" => "CT", "VAL-CB" => "CT", "VAL-CG1" => "CT", "VAL-CG2" => "CT", "VAL-H" => "H", "VAL-HA" => "H1", "VAL-HB" => "HC", "VAL-HG11" => "HC", "VAL-HG12" => "HC", "VAL-HG13" => "HC", "VAL-HG21" => "HC", "VAL-HG22" => "HC", "VAL-HG23" => "HC", "VAL-N" => "N", "VAL-O" => "O", "*-C" => "C", "*-H" => "H", "*-N" => "N", "*-O" => "O", "*-S" => "S", "*-F" => "F"]; - static ref ELE_CHARGES: HashMap<&'static str, f64> = hashmap![ "ALA-C" => 0.5973, "ALA-CA" => 0.0337, "ALA-CB" => -0.1825, "ALA-H" => 0.2719, "ALA-HA" => 0.0823, "ALA-HB1" => 0.0603, "ALA-HB2" => 0.0603, "ALA-HB3" => 0.0603, "ALA-N" => -0.4157, "ALA-O" => -0.5679, "ARG-C" => 0.7341, "ARG-CA" => -0.2637, "ARG-CB" => -0.0007, "ARG-CD" => 0.0486, "ARG-CG" => 0.039, "ARG-CZ" => 0.8076, "ARG-H" => 0.2747, "ARG-HA" => 0.156, "ARG-HB2" => 0.0327, "ARG-HB3" => 0.0327, "ARG-HD2" => 0.0687, "ARG-HD3" => 0.0687, "ARG-HE" => 0.3456, "ARG-HG2" => 0.0285, "ARG-HG3" => 0.0285, "ARG-HH11" => 0.4478, "ARG-HH12" => 0.4478, "ARG-HH21" => 0.4478, "ARG-HH22" => 0.4478, "ARG-N" => -0.3479, "ARG-NE" => -0.5295, "ARG-NH1" => -0.8627, "ARG-NH2" => -0.8627, "ARG-O" => -0.5894, @@ -203,7 +208,6 @@ lazy_static! { "TYR-C" => 0.5973, "TYR-CA" => -0.0014, "TYR-CB" => -0.0152, "TYR-CD1" => -0.1906, "TYR-CD2" => -0.1906, "TYR-CE1" => -0.2341, "TYR-CE2" => -0.2341, "TYR-CG" => -0.0011, "TYR-CZ" => 0.3226, "TYR-H" => 0.2719, "TYR-HA" => 0.0876, "TYR-HB2" => 0.0295, "TYR-HB3" => 0.0295, "TYR-HD1" => 0.1699, "TYR-HD2" => 0.1699, "TYR-HE1" => 0.1656, "TYR-HE2" => 0.1656, "TYR-HH" => 0.3992, "TYR-N" => -0.4157, "TYR-O" => -0.5679, "TYR-OH" => -0.5579, "VAL-C" => 0.5973, "VAL-CA" => -0.0875, "VAL-CB" => 0.2985, "VAL-CG1" => -0.3192, "VAL-CG2" => -0.3192, "VAL-H" => 0.2719, "VAL-HA" => 0.0969, "VAL-HB" => -0.0297, "VAL-HG11" => 0.0791, "VAL-HG12" => 0.0791, "VAL-HG13" => 0.0791, "VAL-HG21" => 0.0791, "VAL-HG22" => 0.0791, "VAL-HG23" => 0.0791, "VAL-N" => -0.4157, "VAL-O" => -0.5679, "*-C" => 0.5973, "*-H" => 0.2719, "*-N" => -0.4157, "*-O" => -0.5679, "*-S" => -0.2737, "*-F" => -0.342]; - static ref NT_ELE_CHARGES: HashMap<&'static str, f64> = hashmap![ "ACE-C" => 0.5972, "ACE-CH3" => -0.3662, "ACE-HH31" => 0.1123, "ACE-HH32" => 0.1123, "ACE-HH33" => 0.1123, "ACE-O" => -0.5679, "ALA-C" => 0.6163, "ALA-CA" => 0.0962, "ALA-CB" => -0.0597, "ALA-H1" => 0.1997, "ALA-H2" => 0.1997, "ALA-H3" => 0.1997, "ALA-HA" => 0.0889, "ALA-HB1" => 0.03, "ALA-HB2" => 0.03, "ALA-HB3" => 0.03, "ALA-N" => 0.1414, "ALA-O" => -0.5722, @@ -246,9 +250,13 @@ pub struct PYDOCKDockingModel { } impl<'a> PYDOCKDockingModel { - - fn new(structure: &'a PDB, active_restraints: &'a [String], passive_restraints: &'a [String], - nmodes: &[f64], num_anm: usize) -> PYDOCKDockingModel { + fn new( + structure: &'a PDB, + active_restraints: &'a [String], + passive_restraints: &'a [String], + nmodes: &[f64], + num_anm: usize, + ) -> PYDOCKDockingModel { let mut model = PYDOCKDockingModel { atoms: Vec::new(), coordinates: Vec::new(), @@ -285,10 +293,12 @@ impl<'a> PYDOCKDockingModel { match model.active_restraints.get_mut(&res_id) { Some(atom_indexes) => { atom_indexes.push(atom_index as usize); - }, + } None => { - model.active_restraints.insert(res_id.to_string(), vec![atom_index as usize]); - }, + model + .active_restraints + .insert(res_id.to_string(), vec![atom_index as usize]); + } } } @@ -296,10 +306,12 @@ impl<'a> PYDOCKDockingModel { match model.passive_restraints.get_mut(&res_id) { Some(atom_indexes) => { atom_indexes.push(atom_index as usize); - }, + } None => { - model.passive_restraints.insert(res_id.to_string(), vec![atom_index as usize]); - }, + model + .passive_restraints + .insert(res_id.to_string(), vec![atom_index as usize]); + } } } @@ -317,7 +329,10 @@ impl<'a> PYDOCKDockingModel { _ => panic!("PYDOCK Error: Atom [{:?}] not supported", atom_id), } } else { - warn!("PYDOCK Warning: Atom [{:?}] not supported, trying generic", atom_id); + warn!( + "PYDOCK Warning: Atom [{:?}] not supported, trying generic", + atom_id + ); let atom_element = match atom_name.chars().nth(0) { Some(element) => element, _ => panic!("PYDOCK Error: Atom element could not be guessed from [{:?}]", atom_name), @@ -328,18 +343,19 @@ impl<'a> PYDOCKDockingModel { _ => panic!("PYDOCK Error: Atom [{:?}] not supported", atom_id), } } - }, + } }; // Assign electrostatics charge let ele_charge = match ELE_CHARGES.get(&*atom_id) { Some(&charge) => charge, - _ => { - match NT_ELE_CHARGES.get(&*atom_id) { - Some(&charge) => charge, - _ => panic!("PYDOCK Error: Atom [{:?}] electrostatics charge not found", atom_id), - } - } + _ => match NT_ELE_CHARGES.get(&*atom_id) { + Some(&charge) => charge, + _ => panic!( + "PYDOCK Error: Atom [{:?}] electrostatics charge not found", + atom_id + ), + }, }; model.ele_charges.push(ele_charge); @@ -371,16 +387,35 @@ pub struct PYDOCK { pub use_anm: bool, } - impl<'a> PYDOCK { - - pub fn new(receptor: PDB, rec_active_restraints: Vec, rec_passive_restraints: Vec, - rec_nmodes: Vec, rec_num_anm: usize, - ligand: PDB, lig_active_restraints: Vec, lig_passive_restraints: Vec, - lig_nmodes: Vec, lig_num_anm: usize, use_anm: bool) -> Box { + pub fn new( + receptor: PDB, + rec_active_restraints: Vec, + rec_passive_restraints: Vec, + rec_nmodes: Vec, + rec_num_anm: usize, + ligand: PDB, + lig_active_restraints: Vec, + lig_passive_restraints: Vec, + lig_nmodes: Vec, + lig_num_anm: usize, + use_anm: bool, + ) -> Box { let d = PYDOCK { - receptor: PYDOCKDockingModel::new(&receptor, &rec_active_restraints, &rec_passive_restraints, &rec_nmodes, rec_num_anm), - ligand: PYDOCKDockingModel::new(&ligand, &lig_active_restraints, &lig_passive_restraints, &lig_nmodes, lig_num_anm), + receptor: PYDOCKDockingModel::new( + &receptor, + &rec_active_restraints, + &rec_passive_restraints, + &rec_nmodes, + rec_num_anm, + ), + ligand: PYDOCKDockingModel::new( + &ligand, + &lig_active_restraints, + &lig_passive_restraints, + &lig_nmodes, + lig_num_anm, + ), use_anm, }; Box::new(d) @@ -388,10 +423,13 @@ impl<'a> PYDOCK { } impl Score for PYDOCK { - - fn energy(&self, translation: &[f64], rotation: &Quaternion, - rec_nmodes: &[f64], lig_nmodes: &[f64]) -> f64 { - + fn energy( + &self, + translation: &[f64], + rotation: &Quaternion, + rec_nmodes: &[f64], + lig_nmodes: &[f64], + ) -> f64 { // Clone receptor coordinates let mut receptor_coordinates: Vec<[f64; 3]> = self.receptor.coordinates.clone(); let rec_num_atoms = receptor_coordinates.len(); @@ -412,9 +450,12 @@ impl Score for PYDOCK { for i_nm in 0usize..self.ligand.num_anm { // (num_anm, num_atoms, 3) -> 1d // Endianness: i = i_nm * num_atoms * 3 + i_atom * 3 + coord - coordinate[0] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3] * lig_nmodes[i_nm]; - coordinate[1] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 1] * lig_nmodes[i_nm]; - coordinate[2] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 2] * lig_nmodes[i_nm]; + coordinate[0] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3] + * lig_nmodes[i_nm]; + coordinate[1] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 1] + * lig_nmodes[i_nm]; + coordinate[2] += self.ligand.nmodes[i_nm * lig_num_atoms * 3 + i_atom * 3 + 2] + * lig_nmodes[i_nm]; } } } @@ -425,9 +466,14 @@ impl Score for PYDOCK { for i_nm in 0usize..self.receptor.num_anm { // (num_anm, num_atoms, 3) -> 1d // Endianness: i = i_nm * num_atoms * 3 + i_atom * 3 + coord - coordinate[0] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3] * rec_nmodes[i_nm]; - coordinate[1] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3 + 1] * rec_nmodes[i_nm]; - coordinate[2] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3 + 2] * rec_nmodes[i_nm]; + coordinate[0] += self.receptor.nmodes[i_nm * rec_num_atoms * 3 + i_atom * 3] + * rec_nmodes[i_nm]; + coordinate[1] += self.receptor.nmodes + [i_nm * rec_num_atoms * 3 + i_atom * 3 + 1] + * rec_nmodes[i_nm]; + coordinate[2] += self.receptor.nmodes + [i_nm * rec_num_atoms * 3 + i_atom * 3 + 2] + * rec_nmodes[i_nm]; } } } @@ -442,11 +488,14 @@ impl Score for PYDOCK { let y1 = ra[1]; let z1 = ra[2]; for (j, la) in ligand_coordinates.iter().enumerate() { - let distance2 = (x1-la[0])*(x1-la[0]) + (y1-la[1])*(y1-la[1]) + (z1-la[2])*(z1-la[2]); + let distance2 = (x1 - la[0]) * (x1 - la[0]) + + (y1 - la[1]) * (y1 - la[1]) + + (z1 - la[2]) * (z1 - la[2]); // Electrostatics energy if distance2 <= ELEC_DIST_CUTOFF2 { - let mut atom_elec = self.receptor.ele_charges[i] * self.ligand.ele_charges[j] / distance2; + let mut atom_elec = + self.receptor.ele_charges[i] * self.ligand.ele_charges[j] / distance2; if atom_elec > ELEC_MAX_CUTOFF { atom_elec = ELEC_MAX_CUTOFF; } @@ -458,10 +507,11 @@ impl Score for PYDOCK { // Van der Waals energy if distance2 <= VDW_DIST_CUTOFF2 { - let vdw_energy = (self.receptor.vdw_charges[i] * self.ligand.vdw_charges[j]).sqrt(); + let vdw_energy = + (self.receptor.vdw_charges[i] * self.ligand.vdw_charges[j]).sqrt(); let vdw_radius = self.receptor.vdw_radii[i] + self.ligand.vdw_radii[j]; let p6 = vdw_radius.powi(6) / distance2.powi(3); - let mut k = vdw_energy * (p6*p6 - 2.0 * p6); + let mut k = vdw_energy * (p6 * p6 - 2.0 * p6); if k > VDW_CUTOFF { k = VDW_CUTOFF; } @@ -479,10 +529,10 @@ impl Score for PYDOCK { let score = (total_elec + total_vdw) * -1.0; // Bias the scoring depending on satisfied restraints - let perc_receptor_restraints: f64 = satisfied_restraints(&interface_receptor, - &self.receptor.active_restraints); - let perc_ligand_restraints: f64 = satisfied_restraints(&interface_ligand, - &self.ligand.active_restraints); + let perc_receptor_restraints: f64 = + satisfied_restraints(&interface_receptor, &self.receptor.active_restraints); + let perc_ligand_restraints: f64 = + satisfied_restraints(&interface_ligand, &self.ligand.active_restraints); // Take into account membrane intersection let mut membrane_penalty: f64 = 0.0; let intersection = membrane_intersection(&interface_receptor, &self.receptor.membrane); @@ -497,8 +547,8 @@ impl Score for PYDOCK { #[cfg(test)] mod tests { use super::*; - use std::env; use crate::qt::Quaternion; + use std::env; #[test] fn test_1azp() { @@ -509,12 +559,26 @@ mod tests { let test_path: String = format!("{}/tests/1azp", cargo_path); let receptor_filename: String = format!("{}/1azp_receptor.pdb", test_path); - let (receptor, _errors) = pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); + let (receptor, _errors) = + pdbtbx::open(&receptor_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); let ligand_filename: String = format!("{}/1azp_ligand.pdb", test_path); - let (ligand, _errors) = pdbtbx::open(&ligand_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); - - let scoring = PYDOCK::new(receptor, Vec::new(), Vec::new(), Vec::new(), 0, ligand, Vec::new(), Vec::new(), Vec::new(), 0, false); + let (ligand, _errors) = + pdbtbx::open(&ligand_filename, pdbtbx::StrictnessLevel::Strict).unwrap(); + + let scoring = PYDOCK::new( + receptor, + Vec::new(), + Vec::new(), + Vec::new(), + 0, + ligand, + Vec::new(), + Vec::new(), + Vec::new(), + 0, + false, + ); let translation = vec![0., 0., 0.]; let rotation = Quaternion::default(); @@ -522,4 +586,3 @@ mod tests { assert_eq!(energy, -364.88126358158974); } } - diff --git a/src/qt.rs b/src/qt.rs index efe4f21..ea44a1f 100644 --- a/src/qt.rs +++ b/src/qt.rs @@ -1,14 +1,13 @@ -use std::ops; -use std::f64; +use super::constants::LINEAR_THRESHOLD; use rand::Rng; +use std::f64; use std::f64::consts::PI; -use super::constants::LINEAR_THRESHOLD; +use std::ops; -fn float_equals(x:f64, y:f64) -> bool { +fn float_equals(x: f64, y: f64) -> bool { (x - y).abs() < f64::EPSILON } - #[derive(Debug, Copy, Clone)] pub struct Quaternion { pub w: f64, @@ -17,48 +16,42 @@ pub struct Quaternion { pub z: f64, } - impl Quaternion { - pub fn new(w:f64, x:f64, y:f64, z:f64) -> Quaternion { - Quaternion { - w, - x, - y, - z - } + pub fn new(w: f64, x: f64, y: f64, z: f64) -> Quaternion { + Quaternion { w, x, y, z } } pub fn conjugate(&self) -> Quaternion { - Quaternion::new(self.w, -self.x, -self.y, -self.z) - } + Quaternion::new(self.w, -self.x, -self.y, -self.z) + } - pub fn dot(&self, other: Quaternion) -> f64 { - self.w*other.w + self.x*other.x + self.y*other.y + self.z*other.z - } + pub fn dot(&self, other: Quaternion) -> f64 { + self.w * other.w + self.x * other.x + self.y * other.y + self.z * other.z + } - pub fn norm2(&self) -> f64 { - self.w*self.w + self.x*self.x + self.y*self.y + self.z*self.z - } + pub fn norm2(&self) -> f64 { + self.w * self.w + self.x * self.x + self.y * self.y + self.z * self.z + } - pub fn norm(&self) -> f64 { - (self.w*self.w + self.x*self.x + self.y*self.y + self.z*self.z).sqrt() - } + pub fn norm(&self) -> f64 { + (self.w * self.w + self.x * self.x + self.y * self.y + self.z * self.z).sqrt() + } pub fn normalize(&mut self) { - let norm = self.norm(); - self.w /= norm; + let norm = self.norm(); + self.w /= norm; self.x /= norm; self.y /= norm; self.z /= norm; } - pub fn inverse(&self) -> Quaternion{ - self.conjugate() / self.norm2() + pub fn inverse(&self) -> Quaternion { + self.conjugate() / self.norm2() } pub fn distance(&self, other: Quaternion) -> f64 { - let dot = self.dot(other); - 1.0 - dot*dot + let dot = self.dot(other); + 1.0 - dot * dot } pub fn rotate(&self, vec3: Vec) -> Vec { @@ -71,10 +64,10 @@ impl Quaternion { *self * (1.0 - t) + other * t } - pub fn slerp(&self, other: &Quaternion, t: f64) -> Quaternion { - let mut q1 = *self; - let mut q2 = *other; - q1.normalize(); + pub fn slerp(&self, other: &Quaternion, t: f64) -> Quaternion { + let mut q1 = *self; + let mut q2 = *other; + q1.normalize(); q2.normalize(); let mut q_dot = q1.dot(q2); @@ -83,17 +76,17 @@ impl Quaternion { q1 = -q1; q_dot *= -1.0; } - + if q_dot > LINEAR_THRESHOLD { // Linear interpolation if quaternions are too close - let mut result = q1 + (q2-q1) * t; + let mut result = q1 + (q2 - q1) * t; result.normalize(); result } else { - q_dot =((q_dot).min(1.0)).max(-1.0); + q_dot = ((q_dot).min(1.0)).max(-1.0); let omega = q_dot.acos(); let so = omega.sin(); - q1 * (((1.0-t)*omega).sin() / so) + q2 * ((t*omega).sin()/so) + q1 * (((1.0 - t) * omega).sin() / so) + q2 * ((t * omega).sin() / so) } } @@ -102,15 +95,15 @@ impl Quaternion { let u2 = rng.gen::(); let u3 = rng.gen::(); Quaternion::new( - (1.0-u1).sqrt() * (2.0 * PI * u2).sin(), - (1.0-u1).sqrt() * (2.0 * PI * u2).cos(), + (1.0 - u1).sqrt() * (2.0 * PI * u2).sin(), + (1.0 - u1).sqrt() * (2.0 * PI * u2).cos(), u1.sqrt() * (2.0 * PI * u3).sin(), - u1.sqrt() * (2.0 * PI * u3).cos() + u1.sqrt() * (2.0 * PI * u3).cos(), ) } - } +} - impl Default for Quaternion { +impl Default for Quaternion { fn default() -> Quaternion { Quaternion { w: 1.0, @@ -125,7 +118,12 @@ impl ops::Sub for Quaternion { type Output = Self; fn sub(self, other: Quaternion) -> Self::Output { - Quaternion::new(self.w-other.w, self.x-other.x, self.y-other.y, self.z-other.z) + Quaternion::new( + self.w - other.w, + self.x - other.x, + self.y - other.y, + self.z - other.z, + ) } } @@ -133,16 +131,21 @@ impl ops::Add for Quaternion { type Output = Self; fn add(self, other: Quaternion) -> Self::Output { - Quaternion::new(self.w+other.w, self.x+other.x, self.y+other.y, self.z+other.z) + Quaternion::new( + self.w + other.w, + self.x + other.x, + self.y + other.y, + self.z + other.z, + ) } } impl PartialEq for Quaternion { fn eq(&self, other: &Self) -> bool { - float_equals(self.w, other.w) && - float_equals(self.x, other.x) && - float_equals(self.y, other.y) && - float_equals(self.z, other.z) + float_equals(self.w, other.w) + && float_equals(self.x, other.x) + && float_equals(self.y, other.y) + && float_equals(self.z, other.z) } } impl Eq for Quaternion {} @@ -159,39 +162,48 @@ impl ops::Mul for Quaternion { type Output = Self; fn mul(self, scalar: f64) -> Self::Output { - Quaternion::new(scalar*self.w, scalar*self.x, scalar*self.y, scalar*self.z) + Quaternion::new( + scalar * self.w, + scalar * self.x, + scalar * self.y, + scalar * self.z, + ) } } impl ops::Mul for Quaternion { - type Output = Self; - - fn mul(self, other: Quaternion) -> Self::Output { - Quaternion::new( - self.w * other.w - self.x * other.x - self.y * other.y - self.z * other.z, - self.w * other.x + self.x * other.w + self.y * other.z - self.z * other.y, - self.w * other.y - self.x * other.z + self.y * other.w + self.z * other.x, - self.w * other.z + self.x * other.y - self.y * other.x + self.z * other.w - ) - } + type Output = Self; + + fn mul(self, other: Quaternion) -> Self::Output { + Quaternion::new( + self.w * other.w - self.x * other.x - self.y * other.y - self.z * other.z, + self.w * other.x + self.x * other.w + self.y * other.z - self.z * other.y, + self.w * other.y - self.x * other.z + self.y * other.w + self.z * other.x, + self.w * other.z + self.x * other.y - self.y * other.x + self.z * other.w, + ) + } } impl ops::Div for Quaternion { type Output = Self; fn div(self, scalar: f64) -> Self::Output { - Quaternion::new(self.w/scalar, self.x/scalar, self.y/scalar, self.z/scalar) + Quaternion::new( + self.w / scalar, + self.x / scalar, + self.y / scalar, + self.z / scalar, + ) } } - #[cfg(test)] mod tests { use super::*; #[test] fn quaternion_default() { - let q1:Quaternion = Default::default(); + let q1: Quaternion = Default::default(); assert_eq!(q1.w, 1.0); assert_eq!(q1.x, 0.0); assert_eq!(q1.y, 0.0); @@ -222,14 +234,14 @@ mod tests { #[test] fn quaternion_eq() { - let q1:Quaternion = Default::default(); - let q2:Quaternion = Default::default(); + let q1: Quaternion = Default::default(); + let q2: Quaternion = Default::default(); assert!(q1 == q2); - let q3:Quaternion = Default::default(); - let q4:Quaternion = Quaternion::new(1.000000000000001, 0.0, 0.0, 0.0); + let q3: Quaternion = Default::default(); + let q4: Quaternion = Quaternion::new(1.000000000000001, 0.0, 0.0, 0.0); assert!(q3 != q4); - let q5:Quaternion = Default::default(); - let q6:Quaternion = Quaternion::new(1.0000000000000001, 0.0, 0.0, 0.0); + let q5: Quaternion = Default::default(); + let q6: Quaternion = Quaternion::new(1.0000000000000001, 0.0, 0.0, 0.0); assert!(q5 == q6); } @@ -259,18 +271,18 @@ mod tests { let q1 = Quaternion::new(1.0, 0.0, 0.0, 2.0); let q2 = Quaternion::new(3.0, -1.0, 4.0, 3.0); let expected = Quaternion::new(-3.0, -9.0, 2.0, 9.0); - assert!(expected == q1*q2); + assert!(expected == q1 * q2); let q1 = Quaternion::new(1.0, 0.0, 0.0, 2.0); let q2 = Quaternion::new(3.0, -1.0, 4.0, 3.0); let expected = Quaternion::new(-3.0, 7.0, 6.0, 9.0); - assert!(expected == q2*q1); + assert!(expected == q2 * q1); let q1 = Quaternion::new(1.0, 0.0, 0.0, 2.0); let q2 = Quaternion::new(3.0, -1.0, 4.0, 3.0); - let q3 = Quaternion::new(1.0/2.0, -3.0, 2.0, 9.0); - let expected = Quaternion::new(-147.0/2.0, 97.0/2.0, -93.0, 19.0/2.0); - assert!(expected == q2*q1*q3); + let q3 = Quaternion::new(1.0 / 2.0, -3.0, 2.0, 9.0); + let expected = Quaternion::new(-147.0 / 2.0, 97.0 / 2.0, -93.0, 19.0 / 2.0); + assert!(expected == q2 * q1 * q3); } #[test] @@ -278,13 +290,13 @@ mod tests { let q1 = Quaternion::new(1.0, 0.0, 0.0, 2.0); let q2 = Quaternion::new(3.0, -1.0, 4.0, 3.0); let expected = Quaternion::new(35.0, 0.0, 0.0, 0.0); - assert!( (q1*q2).conjugate() == q2.conjugate() * q1.conjugate() ); - assert!(expected == q2.conjugate()*q2); + assert!((q1 * q2).conjugate() == q2.conjugate() * q1.conjugate()); + assert!(expected == q2.conjugate() * q2); } #[test] fn test_dot_product() { - let q = Quaternion::new(2_f64.sqrt()/2.0, 0.0, 2_f64.sqrt()/2.0, 0.0); + let q = Quaternion::new(2_f64.sqrt() / 2.0, 0.0, 2_f64.sqrt() / 2.0, 0.0); assert_eq!(1.0000000000000002, q.dot(q)); } @@ -293,13 +305,18 @@ mod tests { let q1 = Quaternion::new(1.0, -3.0, 4.0, 3.0); let q2 = Quaternion::new(3.0, -1.0, 4.0, 3.0); assert_eq!(5.916079783099616, q1.norm()); - assert_eq!((q1*q2).norm(), q1.norm()*q2.norm()); + assert_eq!((q1 * q2).norm(), q1.norm() * q2.norm()); } #[test] fn test_normalize() { let mut q1 = Quaternion::new(1.0, -3.0, 4.0, 3.0); - let expected = Quaternion::new(0.1690308509457033, -0.50709255283711, 0.6761234037828132, 0.50709255283711); + let expected = Quaternion::new( + 0.1690308509457033, + -0.50709255283711, + 0.6761234037828132, + 0.50709255283711, + ); q1.normalize(); assert!(expected == q1); } @@ -308,8 +325,8 @@ mod tests { fn test_inverse() { let q1 = Quaternion::new(1.0, 0.0, 0.0, 2.0); let q2 = Quaternion::new(3.0, -1.0, 4.0, 3.0); - let expected = Quaternion::new(-3.0/175.0, 9.0/175.0, -2.0/175.0, -9.0/175.0); - assert!(expected == (q1*q2).inverse()); + let expected = Quaternion::new(-3.0 / 175.0, 9.0 / 175.0, -2.0 / 175.0, -9.0 / 175.0); + assert!(expected == (q1 * q2).inverse()); } #[test] @@ -332,7 +349,7 @@ mod tests { assert_eq!(0.5000000002638181, q1.distance(q2)); } - #[test] + #[test] fn test_distance_composite_rotation() { let q1 = Quaternion::new(1.0, 0.0, 0.0, 0.0); let q2 = Quaternion::new(0.5, 0.5, 0.5, 0.5); @@ -386,7 +403,12 @@ mod tests { fn test_slerp_t_1() { let q1 = Quaternion::new(1.0, 0.0, 0.0, 2.0); let q2 = Quaternion::new(3.0, -1.0, 4.0, 3.0); - let expected = Quaternion::new(0.50709255283711, -0.1690308509457033, 0.6761234037828132, 0.50709255283711); + let expected = Quaternion::new( + 0.50709255283711, + -0.1690308509457033, + 0.6761234037828132, + 0.50709255283711, + ); let s = q1.slerp(&q2, 1.0); @@ -427,12 +449,16 @@ mod tests { #[test] fn test_random_quaternion() { - use rand::SeedableRng; - let mut rng = SeedableRng::seed_from_u64(324324324); + use rand::SeedableRng; + let mut rng = SeedableRng::seed_from_u64(324324324); let q = Quaternion::random(&mut rng); - let expected = Quaternion::new(0.31924330894562036, -0.5980633213833059, - 0.5444724265858514, 0.49391674399349367); + let expected = Quaternion::new( + 0.31924330894562036, + -0.5980633213833059, + 0.5444724265858514, + 0.49391674399349367, + ); assert!(expected == q); } -} \ No newline at end of file +} diff --git a/src/scoring.rs b/src/scoring.rs index 45ed8f8..47372cb 100644 --- a/src/scoring.rs +++ b/src/scoring.rs @@ -1,6 +1,5 @@ -use std::collections::HashMap; use super::qt::Quaternion; - +use std::collections::HashMap; #[derive(Debug)] pub enum Method { @@ -10,14 +9,19 @@ pub enum Method { } pub trait Score { - fn energy(&self, translation: &[f64], rotation: &Quaternion, - rec_nmodes: &[f64], lig_nmodes: &[f64]) -> f64; + fn energy( + &self, + translation: &[f64], + rotation: &Quaternion, + rec_nmodes: &[f64], + lig_nmodes: &[f64], + ) -> f64; } pub fn satisfied_restraints(interface: &[usize], restraints: &HashMap>) -> f64 { // Calculate the percentage of satisfied restraints if restraints.is_empty() { - return 0.0 + return 0.0; } let mut num_residues = 0; for (_k, atom_indexes) in restraints.iter() { @@ -33,11 +37,11 @@ pub fn satisfied_restraints(interface: &[usize], restraints: &HashMap f64 { if membrane.is_empty() { - return 0.0 + return 0.0; } let mut num_beads = 0; for &i_bead in membrane.iter() { num_beads += interface[i_bead]; } num_beads as f64 / membrane.len() as f64 -} \ No newline at end of file +} diff --git a/src/swarm.rs b/src/swarm.rs index 2b13846..9ac7b49 100644 --- a/src/swarm.rs +++ b/src/swarm.rs @@ -1,14 +1,13 @@ -use super::glowworm::Glowworm; use super::glowworm::distance; -use super::scoring::Score; +use super::glowworm::Glowworm; use super::qt::Quaternion; +use super::scoring::Score; use rand::Rng; use std::fs::File; -use std::io::{Write, Error}; - +use std::io::{Error, Write}; pub struct Swarm<'a> { - pub glowworms: Vec>, + pub glowworms: Vec>, } impl<'a> Default for Swarm<'a> { @@ -19,14 +18,20 @@ impl<'a> Default for Swarm<'a> { impl<'a> Swarm<'a> { pub fn new() -> Self { - Swarm { - glowworms: Vec::new(), - } + Swarm { + glowworms: Vec::new(), + } } - pub fn add_glowworms(&mut self, positions: &[Vec], scoring: &'a Box, use_anm: bool, - rec_num_anm: usize, lig_num_anm: usize) { - for (i, position) in positions.iter().enumerate() { + pub fn add_glowworms( + &mut self, + positions: &[Vec], + scoring: &'a Box, + use_anm: bool, + rec_num_anm: usize, + lig_num_anm: usize, + ) { + for (i, position) in positions.iter().enumerate() { // Translation component let translation = vec![position[0], position[1], position[2]]; // Rotation component @@ -34,74 +39,82 @@ impl<'a> Swarm<'a> { // ANM for receptor let mut rec_nmodes: Vec = Vec::new(); if use_anm && rec_num_anm > 0 { - for j in 7..7+rec_num_anm { + for j in 7..7 + rec_num_anm { rec_nmodes.push(positions[i][j]); } } // ANM for ligand let mut lig_nmodes: Vec = Vec::new(); if use_anm && lig_num_anm > 0 { - for j in 7+rec_num_anm..positions[i].len() { + for j in 7 + rec_num_anm..positions[i].len() { lig_nmodes.push(positions[i][j]); } } - let glowworm = Glowworm::new(i as u32, translation, rotation, rec_nmodes, lig_nmodes, scoring, use_anm); - self.glowworms.push(glowworm); - } + let glowworm = Glowworm::new( + i as u32, + translation, + rotation, + rec_nmodes, + lig_nmodes, + scoring, + use_anm, + ); + self.glowworms.push(glowworm); + } } pub fn update_luciferin(&mut self) { - for glowworm in self.glowworms.iter_mut() { - glowworm.compute_luciferin(); - } + for glowworm in self.glowworms.iter_mut() { + glowworm.compute_luciferin(); + } } pub fn movement_phase(&mut self, rng: &mut rand::prelude::StdRng) { - // Save original positions - let mut positions: Vec> = Vec::new(); + // Save original positions + let mut positions: Vec> = Vec::new(); let mut rotations: Vec = Vec::new(); let mut anm_recs: Vec> = Vec::new(); let mut anm_ligs: Vec> = Vec::new(); - for glowworm in self.glowworms.iter(){ - positions.push(glowworm.translation.clone()); + for glowworm in self.glowworms.iter() { + positions.push(glowworm.translation.clone()); rotations.push(glowworm.rotation); anm_recs.push(glowworm.rec_nmodes.clone()); anm_ligs.push(glowworm.lig_nmodes.clone()); - } + } + + // First search for each glowworm's neighbors + let mut neighbors: Vec> = Vec::new(); + for i in 0..self.glowworms.len() { + let mut this_neighbors = Vec::new(); + let g1 = &self.glowworms[i]; + for j in 0..self.glowworms.len() { + if i != j { + let g2 = &self.glowworms[j]; + if g1.luciferin < g2.luciferin { + let distance = distance(g1, g2); + if distance < g1.vision_range { + this_neighbors.push(g2.id); + } + } + } + } + neighbors.push(this_neighbors); + } + + // Second compute probability moving towards the neighbor + let mut luciferins = Vec::new(); + for glowworm in self.glowworms.iter_mut() { + luciferins.push(glowworm.luciferin); + } + for i in 0..self.glowworms.len() { + let glowworm = &mut self.glowworms[i]; + glowworm.neighbors = neighbors[i].clone(); + glowworm.compute_probability_moving_toward_neighbor(&luciferins); + } - // First search for each glowworm's neighbors - let mut neighbors: Vec> = Vec::new(); - for i in 0..self.glowworms.len() { - let mut this_neighbors = Vec::new(); - let g1 = &self.glowworms[i]; - for j in 0..self.glowworms.len() { - if i != j { - let g2 = &self.glowworms[j]; - if g1.luciferin < g2.luciferin { - let distance = distance(g1, g2); - if distance < g1.vision_range { - this_neighbors.push(g2.id); - } - } - } - } - neighbors.push(this_neighbors); - } - - // Second compute probability moving towards the neighbor - let mut luciferins = Vec::new(); - for glowworm in self.glowworms.iter_mut() { - luciferins.push(glowworm.luciferin); - } - for i in 0..self.glowworms.len() { - let glowworm = &mut self.glowworms[i]; - glowworm.neighbors = neighbors[i].clone(); - glowworm.compute_probability_moving_toward_neighbor(&luciferins); - } - - // Finally move to the selected position + // Finally move to the selected position for i in 0..self.glowworms.len() { - let glowworm = &mut self.glowworms[i]; + let glowworm = &mut self.glowworms[i]; let neighbor_id = glowworm.select_random_neighbor(rng.gen::()); let position = &positions[neighbor_id as usize]; let rotation = &rotations[neighbor_id as usize]; @@ -115,11 +128,22 @@ impl<'a> Swarm<'a> { pub fn save(&mut self, step: u32, output_directory: &str) -> Result<(), Error> { let path = format!("{}/gso_{:?}.out", output_directory, step); let mut output = File::create(path)?; - writeln!(output, "#Coordinates RecID LigID Luciferin Neighbor's number Vision Range Scoring")?; + writeln!( + output, + "#Coordinates RecID LigID Luciferin Neighbor's number Vision Range Scoring" + )?; for glowworm in self.glowworms.iter() { - write!(output, "({:.7}, {:.7}, {:.7}, {:.7}, {:.7}, {:.7}, {:.7}", - glowworm.translation[0], glowworm.translation[1], glowworm.translation[2], - glowworm.rotation.w, glowworm.rotation.x, glowworm.rotation.y, glowworm.rotation.z)?; + write!( + output, + "({:.7}, {:.7}, {:.7}, {:.7}, {:.7}, {:.7}, {:.7}", + glowworm.translation[0], + glowworm.translation[1], + glowworm.translation[2], + glowworm.rotation.w, + glowworm.rotation.x, + glowworm.rotation.y, + glowworm.rotation.z + )?; if glowworm.use_anm && !glowworm.rec_nmodes.is_empty() { for i in 0..glowworm.rec_nmodes.len() { write!(output, ", {:.7}", glowworm.rec_nmodes[i])?; @@ -130,8 +154,14 @@ impl<'a> Swarm<'a> { write!(output, ", {:.7}", glowworm.lig_nmodes[i])?; } } - writeln!(output, ") 0 0 {:.8} {:?} {:.3} {:.8}", - glowworm.luciferin, glowworm.neighbors.len(), glowworm.vision_range, glowworm.scoring)?; + writeln!( + output, + ") 0 0 {:.8} {:?} {:.3} {:.8}", + glowworm.luciferin, + glowworm.neighbors.len(), + glowworm.vision_range, + glowworm.scoring + )?; } Ok(()) }