/* ############################################################################################## # # # sdilicon.c # # # # Guillermo Moyna et al. # # # # (c) Department of Chemistry & Biochemistry, University # # of the Sciences in Philadelphia, 2002, 2003 # # # # DESCRIPTION: # # # # sdilicon.c is the standalone version of SDiLiCon, originally implemented in Sybyl as a # # prototype. It is not only the combination of all the functions used to compute chemical # # shifts, but also a wrapper that reads in data from the protein-ligand complex, and has a # # built-in force field that considers vdW/electrostatic interactions, as well as different # # optimization algorithms to do the shift-directed optimization. # # # # The program reads the details of the run from a sdilicon control file with the same format # # as the Sybyl prototype. It also reads either .mol2 or .pdb files to get information on the # # protein-ligand complex. Once the data is read in, the information of the .sdl file is used # # to extract the appropriate XYZ, charge, and vdW data from the .mol2 or .pdb files and make # # sub-arrays with the ligands, rings, anisotropic groups, etc. These arrays are the ones # # employed in the calculations. # # # # Two different methods can be used in the minimization. One is the 'ragdoll' minimizer that # # was implemented in the Sybyl prototype. The other one is a random-incremental pulse search # # (RIPS) to cover a large region of the accessible conformational space. # # # # Once the calculation is complete a new .mol2 or .pdb (depending on the type of input) is # # created with the optimized protein-ligand complex. # # # # The program is called with a series of command line options wich are described within the # # 'main' section. By the way, SDiLiCon means Shift Directed Ligand Conformation. Have fun! # # # # Send comments and suggestions to the author via e-mail at g.moyna@usip.edu # # # ############################################################################################## # # # LICENSE: # # # # This program is free software; you can redistribute it and/or modify it under the terms of # # the GNU General Public License as published by the Free Software Foundation; either # # version 2 of the License, or (at your option) any later version. # # # # This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; # # without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. # # See the GNU General Public License for more details. # # # # You should have received a copy of the GNU General Public License along with this program; # # if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, # # MA 02111-1307 USA # # # ############################################################################################## */ /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Includes Go Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */ #include #include #include #include #include #include #include /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Functions Go Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */ /* This are the vdW assignment function. The vdW radius function uses a 'small' hydrogen atom, but everything else is Tripos... */ void vdW( char AT[6], double &vdWr, double &vdWe, double &CovRad ) { int i; char AtomType[6]; strcpy( AtomType, AT ); for( i = 0; i <= strlen( AtomType ) - 1; ++i ) AtomType[i] = toupper( AtomType[i] ); if( ! strcmp( AtomType, "C.3" ) ) { vdWr = 1.7; vdWe = 0.107; CovRad = 0.762; } else if( ! strcmp( AtomType, "C.2" ) ) { vdWr = 1.7; vdWe = 0.107; CovRad = 0.762; } else if( ! strcmp( AtomType, "C.CAT" ) ) { vdWr = 1.7; vdWe = 0.107; CovRad = 0.762; } else if( ! strcmp( AtomType, "C.AR" ) ) { vdWr = 1.7; vdWe = 0.107; CovRad = 0.762; } else if( ! strcmp( AtomType, "C.1" ) ) { vdWr = 1.7; vdWe = 0.107; CovRad = 0.762; } else if( ! strcmp( AtomType, "C" ) ) { vdWr = 1.7; vdWe = 0.107; CovRad = 0.762; } else if( ! strcmp( AtomType, "N.3" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } else if( ! strcmp( AtomType, "N.2" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } else if( ! strcmp( AtomType, "N.1" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } else if( ! strcmp( AtomType, "N.AR" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } else if( ! strcmp( AtomType, "N.4" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } else if( ! strcmp( AtomType, "N.PL3" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } else if( ! strcmp( AtomType, "N.AM" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } else if( ! strcmp( AtomType, "N" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } else if( ! strcmp( AtomType, "O.3" ) ) { vdWr = 1.52; vdWe = 0.116; CovRad = 0.640; } else if( ! strcmp( AtomType, "O.2" ) ) { vdWr = 1.52; vdWe = 0.116; CovRad = 0.640; } else if( ! strcmp( AtomType, "O.CO2" ) ) { vdWr = 1.52; vdWe = 0.116; CovRad = 0.640; } else if( ! strcmp( AtomType, "O.T3P" ) ) { vdWr = 1.76825; vdWe = 0.15207; CovRad = 0.640; } else if( ! strcmp( AtomType, "O.SPC" ) ) { vdWr = 1.7766; vdWe = 0.1554; CovRad = 0.640; } else if( ! strcmp( AtomType, "O" ) ) { vdWr = 1.52; vdWe = 0.116; CovRad = 0.640; } else if( ! strcmp( AtomType, "P.3" ) ) { vdWr = 1.8; vdWe = 0.314; CovRad = 1.094; } else if( ! strcmp( AtomType, "P" ) ) { vdWr = 1.8; vdWe = 0.314; CovRad = 1.094; } else if( ! strcmp( AtomType, "H" ) ) { vdWr = 1.0;vdWe = 0.042; CovRad = 0.352; } else if( ! strcmp( AtomType, "H.T3P" ) ) { vdWr = 1.008; vdWe = 0; CovRad = 0.352; } else if( ! strcmp( AtomType, "H.SPC" ) ) { vdWr = 1.008; vdWe = 0; CovRad = 0.352; } else if( ! strcmp( AtomType, "S.3" ) ) { vdWr = 1.8; vdWe = 0.314; CovRad = 1.053; } else if( ! strcmp( AtomType, "S.2" ) ) { vdWr = 1.8; vdWe = 0.314; CovRad = 1.053; } else if( ! strcmp( AtomType, "S.O" ) ) { vdWr = 1.7; vdWe = 0.314; CovRad = 1.053; } else if( ! strcmp( AtomType, "S.O2" ) ) { vdWr = 1.7; vdWe = 0.314; CovRad = 1.053; } else if( ! strcmp( AtomType, "S" ) ) { vdWr = 1.75; vdWe = 0.314; CovRad = 1.053; } else if( ! strcmp( AtomType, "BR" ) ) { vdWr = 1.85; vdWe = 0.434; CovRad = 1.187; } else if( ! strcmp( AtomType, "CL" ) ) { vdWr = 1.75; vdWe = 0.314; CovRad = 1.033; } else if( ! strcmp( AtomType, "F" ) ) { vdWr = 1.47; vdWe = 0.109; CovRad = 0.630; } else if( ! strcmp( AtomType, "I" ) ) { vdWr = 1.98; vdWe = 0.623; CovRad = 1.387; } else if( ! strcmp( AtomType, "LP" ) ) { vdWr = 1e-05; vdWe = 1e-05; CovRad = 0.750; } else if( ! strcmp( AtomType, "NA" ) ) { vdWr = 1.2; vdWe = 0.4; CovRad = 1.540; } else if( ! strcmp( AtomType, "K" ) ) { vdWr = 1.2; vdWe = 0.4; CovRad = 2.030; } else if( ! strcmp( AtomType, "CA" ) ) { vdWr = 1.2; vdWe = 0.6; CovRad = 1.740; } else if( ! strcmp( AtomType, "LI" ) ) { vdWr = 1.2; vdWe = 0.4; CovRad = 1.230; } else if( ! strcmp( AtomType, "AL" ) ) { vdWr = 1.2; vdWe = 0.042; CovRad = 1.180; } else if( ! strcmp( AtomType, "DU" ) ) { vdWr = 0; vdWe = 0; CovRad = 0.750; } else if( ! strcmp( AtomType, "DU.C" ) ) { vdWr = 0; vdWe = 0; CovRad = 0.750; } else if( ! strcmp( AtomType, "SI" ) ) { vdWr = 2.1; vdWe = 0.042; CovRad = 1.118; } else if( ! strcmp( AtomType, "HEV" ) ) { vdWr = 1.7; vdWe = 0.107; CovRad = 0.762; } else if( ! strcmp( AtomType, "ANY" ) ) { vdWr = 1.7; vdWe = 0.107; CovRad = 0.762; } else if( ! strcmp( AtomType, "HAL" ) ) { vdWr = 1.75; vdWe = 0.314; CovRad = 1.033; } else if( ! strcmp( AtomType, "HET" ) ) { vdWr = 1.55; vdWe = 0.095; CovRad = 0.676; } return; } /* This is the torsional parameter assignment function. Lifted from tripos... */ void gettorpar( char atX[8], char atA[8], char atB[8], char atY[8], double &pot, int &phase ) { int i; char atlX[8], atlA[8], atlB[8], atlY[8]; strcpy( atlX, atX ); strcpy( atlA, atA ); strcpy( atlB, atB ); strcpy( atlY, atY ); for( i = 0; i <= strlen( atlX ) - 1; ++i ) atlX[i] = toupper( atlX[i] ); for( i = 0; i <= strlen( atlA ) - 1; ++i ) atlA[i] = toupper( atlA[i] ); for( i = 0; i <= strlen( atlB ) - 1; ++i ) atlB[i] = toupper( atlB[i] ); for( i = 0; i <= strlen( atlY ) - 1; ++i ) atlY[i] = toupper( atlY[i] ); if( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) { if( ( ! strcmp( atlX, "C.3" ) ) && ( ! strcmp( atlY, "C.3" ) ) ) { pot = 0.5; phase = 3; } else if( ( ! strcmp( atlX, "C.2" ) ) && ( ! strcmp( atlY, "C.2" ) ) ) { pot = 0.04; phase = 3; } else if( ( ( ! strcmp( atlX, "C.3" ) ) && ( ! strcmp( atlY, "C.2" ) ) ) || ( ( ! strcmp( atlX, "C.2" ) ) && ( ! strcmp( atlY, "C.3" ) ) ) ) { pot = 0.126; phase = 3; } else if( ( ! strcmp( atlX, "H" ) ) || ( ! strcmp( atlY, "H" ) ) ) { pot = 0.32; phase = 3; } else { pot = 0.2; phase = 3; } } else if( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) { pot = 1.424; phase = -2; } else if( ( ! strcmp( atlA, "C.1" ) ) && ( ! strcmp( atlB, "C.1" ) ) ) { pot = 0.0; phase = 1; } else if( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) { pot = 0.6; phase = -2; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) || ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { if( ( ! strcmp( atlX, "C.3" ) ) && ( ! strcmp( atlY, "C.3" ) ) ) { pot = 0.126; phase = 3; } else if( ( ! strcmp( atlX, "C.2" ) ) && ( ! strcmp( atlY, "C.2" ) ) ) { pot = 0.126; phase = -3; } else if( ( ! strcmp( atlX, "H" ) ) && ( ! strcmp( atlY, "H" ) ) ) { pot = 0.274; phase = 3; } else if( ( ( ! strcmp( atlX, "C.2" ) ) && ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlY, "H" ) ) ) || ( ( ! strcmp( atlX, "H" ) ) && ( ! strcmp( atlB, "C.2" ) ) && ( ! strcmp( atlY, "C.2" ) ) ) ) { pot = 0.274; phase = -3; } else if( ( ( ! strcmp( atlX, "H" ) ) && ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlY, "C.2" ) ) ) || ( ( ! strcmp( atlX, "C.2" ) ) && ( ! strcmp( atlB, "C.2" ) ) && ( ! strcmp( atlY, "H" ) ) ) ) { pot = 0.274; phase = 3; } else if( ( ( ! strcmp( atlX, "C.3" ) ) && ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlY, "C.2" ) ) ) || ( ( ! strcmp( atlX, "C.2" ) ) && ( ! strcmp( atlB, "C.2" ) ) && ( ! strcmp( atlY, "C.3" ) ) ) ) { pot = 0.126; phase = 3; } else if( ( ( ! strcmp( atlX, "C.2" ) ) && ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlY, "C.3" ) ) ) || ( ( ! strcmp( atlX, "C.3" ) ) && ( ! strcmp( atlB, "C.2" ) ) && ( ! strcmp( atlY, "C.2" ) ) ) ) { pot = 0.126; phase = -3; } else if( ( ( ! strcmp( atlX, "O.2" ) ) && ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlY, "C.3" ) ) ) || ( ( ! strcmp( atlX, "C.3" ) ) && ( ! strcmp( atlB, "C.2" ) ) && ( ! strcmp( atlY, "O.2" ) ) ) ) { pot = 0.7; phase = -3; } else if( ( ( ! strcmp( atlX, "O.CO2" ) ) && ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlY, "C.3" ) ) ) || ( ( ! strcmp( atlX, "C.3" ) ) && ( ! strcmp( atlB, "C.2" ) ) && ( ! strcmp( atlY, "O.CO2" ) ) ) ) { pot = 0.7; phase = -3; } else if( ( ( ! strcmp( atlX, "C.3" ) ) && ( ! strcmp( atlY, "H" ) ) ) || ( ( ! strcmp( atlX, "H" ) ) && ( ! strcmp( atlY, "C.3" ) ) ) ) { pot = 0.274; phase = 3; } else if( ( ! strcmp( atlX, "C.2" ) ) || ( ! strcmp( atlY, "C.2" ) ) ) { pot = 0.274; phase = 3; } else if( ( ! strcmp( atlX, "C.3" ) ) || ( ! strcmp( atlY, "C.3" ) ) ) { pot = 0.274; phase = 3; } else if( ( ! strcmp( atlX, "H" ) ) || ( ! strcmp( atlY, "H" ) ) ) { pot = 0.274; phase = 3; } else { pot = 0.12; phase = -3; } } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "C.1" ) ) ) || ( ( ! strcmp( atlA, "C.1" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 0.0; phase = 1; } else if( ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "C.1" ) ) ) || ( ( ! strcmp( atlA, "C.1" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) ) { pot = 0.0; phase = 1; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) || ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 0.12; phase = -3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) || ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 1.6; phase = -2; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "C.1" ) ) ) || ( ( ! strcmp( atlA, "C.1" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 0.0; phase = 1; } else if( ( ( ! strcmp( atlA, "C.1" ) ) && ( ! strcmp( atlB, "N.2" ) ) ) || ( ( ! strcmp( atlA, "N.2" ) ) && ( ! strcmp( atlB, "C.1" ) ) ) ) { pot = 0.0; phase = 1; } else if( ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "N.2" ) ) ) || ( ( ! strcmp( atlA, "N.2" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) ) { pot = 12.0; phase = -2; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "N.2" ) ) ) || ( ( ! strcmp( atlA, "N.2" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 0.4; phase = -3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "N.2" ) ) ) || ( ( ! strcmp( atlA, "N.2" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 1.6; phase = -2; } else if( ( ! strcmp( atlA, "N.2" ) ) && ( ! strcmp( atlB, "N.2" ) ) ) { pot = 1.6; phase = -2; } else if( ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "N.3" ) ) ) || ( ( ! strcmp( atlA, "N.3" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) ) { pot = 0.12; phase = -3; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "N.3" ) ) ) || ( ( ! strcmp( atlA, "N.3" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 0.2; phase = 3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "N.3" ) ) ) || ( ( ! strcmp( atlA, "N.3" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 0.12; phase = -3; } else if( ( ! strcmp( atlA, "N.3" ) ) && ( ! strcmp( atlB, "N.3" ) ) ) { pot = 0.2; phase = 3; } else if( ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "N.AM" ) ) ) || ( ( ! strcmp( atlA, "N.AM" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) ) { pot = 6.46; phase = -2; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "N.AM" ) ) ) || ( ( ! strcmp( atlA, "N.AM" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 0.2; phase = 3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "N.AM" ) ) ) || ( ( ! strcmp( atlA, "N.AM" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 1.6; phase = -2; } else if( ( ( ! strcmp( atlA, "N.2" ) ) && ( ! strcmp( atlB, "N.AM" ) ) ) || ( ( ! strcmp( atlA, "N.AM" ) ) && ( ! strcmp( atlB, "N.2" ) ) ) ) { pot = 1.6; phase = -2; } else if( ( ( ! strcmp( atlA, "N.3" ) ) && ( ! strcmp( atlB, "N.AM" ) ) ) || ( ( ! strcmp( atlA, "N.AM" ) ) && ( ! strcmp( atlB, "N.3" ) ) ) ) { pot = 0.12; phase = -3; } else if( ( ! strcmp( atlA, "N.AM" ) ) && ( ! strcmp( atlB, "N.AM" ) ) ) { pot = 1.6; phase = -2; } else if( ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "N.PL3" ) ) ) || ( ( ! strcmp( atlA, "N.PL3" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) ) { pot = 12.0; phase = -2; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "N.PL3" ) ) ) || ( ( ! strcmp( atlA, "N.PL3" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 0.4; phase = -3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "N.PL3" ) ) ) || ( ( ! strcmp( atlA, "N.PL3" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 1.6; phase = -2; } else if( ( ( ! strcmp( atlA, "N.2" ) ) && ( ! strcmp( atlB, "N.PL3" ) ) ) || ( ( ! strcmp( atlA, "N.PL3" ) ) && ( ! strcmp( atlB, "N.2" ) ) ) ) { pot = 1.6; phase = -2; } else if( ( ! strcmp( atlA, "N.PL3" ) ) && ( ! strcmp( atlB, "N.PL3" ) ) ) { pot = 1.6; phase = -2; } else if( ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "O.3" ) ) ) || ( ( ! strcmp( atlA, "O.3" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) ) { pot = 5.8; phase = -2; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "O.3" ) ) ) || ( ( ! strcmp( atlA, "O.3" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 1.2; phase = 3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "O.3" ) ) ) || ( ( ! strcmp( atlA, "O.3" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 1.2; phase = -2; } else if( ( ( ! strcmp( atlA, "N.2" ) ) && ( ! strcmp( atlB, "O.3" ) ) ) || ( ( ! strcmp( atlA, "O.3" ) ) && ( ! strcmp( atlB, "N.2" ) ) ) ) { pot = 1.0; phase = 2; } else if( ( ( ! strcmp( atlA, "N.3" ) ) && ( ! strcmp( atlB, "O.3" ) ) ) || ( ( ! strcmp( atlA, "O.3" ) ) && ( ! strcmp( atlB, "N.3" ) ) ) ) { pot = 0.2; phase = 3; } else if( ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "P.3" ) ) ) || ( ( ! strcmp( atlA, "P.3" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) ) { pot = 1.0; phase = -2; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "P.3" ) ) ) || ( ( ! strcmp( atlA, "P.3" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 0.4; phase = 3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "P.3" ) ) ) || ( ( ! strcmp( atlA, "P.3" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 1.0; phase = 3; } else if( ( ( ! strcmp( atlA, "O.3" ) ) && ( ! strcmp( atlB, "P.3" ) ) ) || ( ( ! strcmp( atlA, "P.3" ) ) && ( ! strcmp( atlB, "O.3" ) ) ) ) { pot = 0.4; phase = 3; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "S.2" ) ) ) || ( ( ! strcmp( atlA, "S.2" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 0.4; phase = 3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "S.2" ) ) ) || ( ( ! strcmp( atlA, "S.2" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 1.0; phase = 3; } else if( ( ( ! strcmp( atlA, "N.3" ) ) && ( ! strcmp( atlB, "S.2" ) ) ) || ( ( ! strcmp( atlA, "S.2" ) ) && ( ! strcmp( atlB, "N.3" ) ) ) ) { pot = 0.4; phase = 3; } else if( ( ( ! strcmp( atlA, "C.2" ) ) && ( ! strcmp( atlB, "S.3" ) ) ) || ( ( ! strcmp( atlA, "S.3" ) ) && ( ! strcmp( atlB, "C.2" ) ) ) ) { pot = 1.0; phase = -2; } else if( ( ( ! strcmp( atlA, "C.3" ) ) && ( ! strcmp( atlB, "S.3" ) ) ) || ( ( ! strcmp( atlA, "S.3" ) ) && ( ! strcmp( atlB, "C.3" ) ) ) ) { pot = 0.4; phase = 3; } else if( ( ( ! strcmp( atlA, "C.AR" ) ) && ( ! strcmp( atlB, "S.3" ) ) ) || ( ( ! strcmp( atlA, "S.3" ) ) && ( ! strcmp( atlB, "C.AR" ) ) ) ) { pot = 1.0; phase = 3; } else if( ( ! strcmp( atlA, "S.3" ) ) && ( ! strcmp( atlB, "S.3" ) ) ) { pot = 4.0; phase = 3; } else { pot = 0.2; phase = 3; } return; } /* This is the translation function. Excedingly simple... */ void translate( double **LigArray, int nLigAt, double delta, int coord ) { int i; for( i = 0; i <= nLigAt - 1; ++i ) LigArray[i][coord] += delta; return; } /* This is the rotation function - Not so simple. For future reference, here's the rotation matrix: | cos(phi) -sin(phi) 0 | around z: | sin(phi) cos(phi) 0 | | 0 0 1 | | cos(phi) 0 -sin(phi) | around y: | 0 1 0 | | sin(phi) 0 cos(phi) | | 1 0 0 | around x: | 0 cos(phi) -sin(phi) | | 0 sin(phi) cos(phi) | */ void rotate( double **LigArray, int nLigAt, double theta, int coord ) { double relorg_x = 0, relorg_y = 0, relorg_z = 0; double relpos_x, relpos_y, relpos_z, rotated_x, rotated_y, rotated_z; int i; theta *= ( M_PI / 180.0 ); for( i = 0; i <= nLigAt - 1; ++i ) { relorg_x += LigArray[i][0]; relorg_y += LigArray[i][1]; relorg_z += LigArray[i][2]; } relorg_x /= nLigAt; relorg_y /= nLigAt; relorg_z /= nLigAt; if( coord == 2 ) { for( i = 0; i <= nLigAt - 1; ++i ) { relpos_x = LigArray[i][0] - relorg_x; relpos_y = LigArray[i][1] - relorg_y; rotated_x = cos( theta ) * relpos_x - sin( theta ) * relpos_y; rotated_y = sin( theta ) * relpos_x + cos( theta ) * relpos_y; LigArray[i][0] = rotated_x + relorg_x; LigArray[i][1] = rotated_y + relorg_y; } } else if( coord == 1 ) { for( i = 0; i <= nLigAt - 1; ++i ) { relpos_x = LigArray[i][0] - relorg_x; relpos_z = LigArray[i][2] - relorg_z; rotated_x = cos( theta ) * relpos_x - sin( theta ) * relpos_z; rotated_z = sin( theta ) * relpos_x + cos( theta ) * relpos_z; LigArray[i][0] = rotated_x + relorg_x; LigArray[i][2] = rotated_z + relorg_z; } } else if( coord == 0 ) { for( i = 0; i <= nLigAt - 1; ++i ) { relpos_y = LigArray[i][1] - relorg_y; relpos_z = LigArray[i][2] - relorg_z; rotated_y = cos( theta ) * relpos_y - sin( theta ) * relpos_z; rotated_z = sin( theta ) * relpos_y + cos( theta ) * relpos_z; LigArray[i][1] = rotated_y + relorg_y; LigArray[i][2] = rotated_z + relorg_z; } } return; } /* Define the functions for the elliptics... */ double ellint( double qc, double p, double a, double b ) { double e, f, g, em, q; bool go = true; qc = fabs( qc ); e = qc; em = 1.0; if( p > 0.0 ) { p = sqrt( p ); b /= p; } else { f = qc * qc; q = 1.0 - f; g = 1.0 - p; f -= p; q *= b - a * p; p = sqrt( f / g ); a = ( a - b ) / g; b = - q / ( g * g * p ) + a * p; } while( go ) { f = a; a += b / p; g = e / p; b = 2.0 * ( b + f * g ); p += g; g = em; em += qc; if( fabs( g - qc ) <= g * 0.0003 ) break; qc = sqrt( e ); qc *= 2.0; e = qc * em; } return 0.5 * M_PI * ( b + a * em ) / ( em * ( em + p ) ); } double ellpk( double k ) { return ellint( sqrt( 1.0 - k * k ), 1.0, 1.0, 1.0 ); } double ellpe( double k ) { double kc; kc = sqrt( 1.0 - k * k ); return ellint( kc, 1.0, 1.0, kc * kc ); } /* Define the anisotropy shift function... */ double anishift( double H[3], double **AniArray, int nAniAt, double chiiso, double chiani, char bondtype[4] ) { double porb_x, porb_y, porb_z, porb_norm, grp_normx, grp_normy, grp_normz, grp_nmag, r_H, cos_phi, cos_theta; double shift = 0, ani_centx = 0, ani_centy = 0, ani_centz = 0; int i; /* Now loop over the three atoms in the group, and get their coordinates, and get the center of the group/bond... */ for( i = 0; i <= nAniAt - 1; ++i ) { ani_centx += AniArray[i][0]; ani_centy += AniArray[i][1]; ani_centz += AniArray[i][2]; } ani_centx /= nAniAt; ani_centy /= nAniAt; ani_centz /= nAniAt; /* If this was a originally a double bond, move the center back to the center of the bond... */ if( ( ! strcmp( bondtype, "BND" ) ) && ( nAniAt == 3 ) ) { ani_centx = 3 * ( ani_centx - AniArray[2][0] / 3 ) / 2; ani_centy = 3 * ( ani_centy - AniArray[2][1] / 3 ) / 2; ani_centz = 3 * ( ani_centz - AniArray[2][2] / 3 ) / 2; } /* Now we should get the distance from the proton to the center of the group/bond... */ r_H = sqrt( ( H[0] - ani_centx ) * ( H[0] - ani_centx ) + ( H[1] - ani_centy ) * ( H[1] - ani_centy ) + ( H[2] - ani_centz ) * ( H[2] - ani_centz ) ); /* If we have three atoms, we have a group or double bond. Calculate the perpendicular to the group's plane (i.e., the cross-product of the different vectors). The plane coordinates are called porb, because they can also be used as bond vector coordinates (i.e., plane or bond)... */ if( nAniAt == 3 ) { porb_x = ( AniArray[1][1] - AniArray[0][1] ) * ( AniArray[2][2] - AniArray[0][2] ) - ( AniArray[1][2] - AniArray[0][2] ) * ( AniArray[2][1] - AniArray[0][1] ); porb_y = ( AniArray[1][2] - AniArray[0][2] ) * ( AniArray[2][0] - AniArray[0][0] ) - ( AniArray[1][0] - AniArray[0][0] ) * ( AniArray[2][2] - AniArray[0][2] ); porb_z = ( AniArray[1][0] - AniArray[0][0] ) * ( AniArray[2][1] - AniArray[0][1] ) - ( AniArray[1][1] - AniArray[0][1] ) * ( AniArray[2][0] - AniArray[0][0] ); porb_norm = sqrt( porb_x * porb_x + porb_y * porb_y + porb_z * porb_z ); cos_phi = ( porb_x * ( H[0] - ani_centx ) + porb_y * ( H[1] - ani_centy ) + porb_z * ( H[2] - ani_centz ) ) / ( porb_norm * r_H ); /* Depending on the selection of symetric or non-symetric bond/group, we calculate the anisotropy. For symetric bond/groups, we have all we need... */ shift = chiiso * ( 1 - 3 * cos_phi * cos_phi ) / ( 3 * r_H * r_H * r_H ); /* For non-symetric groups/bonds, we still have to calculate cos_theta... */ if( chiani > 1.0e-6 ) { /* Get the normal to the plane made by the group's normal and the bond vector... */ grp_normx = porb_y * ( ani_centz - AniArray[0][2] ) - porb_z * ( ani_centy - AniArray[0][1] ); grp_normy = porb_z * ( ani_centx - AniArray[0][0] ) - porb_x * ( ani_centz - AniArray[0][2] ); grp_normz = porb_x * ( ani_centy - AniArray[0][1] ) - porb_y * ( ani_centx - AniArray[0][0] ); grp_nmag = sqrt( grp_normx * grp_normx + grp_normy * grp_normy + grp_normz * grp_normz ); cos_theta = ( grp_normx * ( H[0] - ani_centx ) + grp_normy * ( H[1] - ani_centy ) + grp_normz * ( H[2] - ani_centz ) ) / ( r_H * grp_nmag ); shift = shift + chiani * ( 1 - 3 * cos_theta * cos_theta ) / ( 3 * r_H * r_H * r_H ); } } /* Now, if we passed only two atoms as arguments, we have a single/triple bond. There's not much to do but calculate the dot product between the bond vector and the radial vector, and then use the regular McConnell equation... */ else { porb_x = AniArray[0][0] - ani_centx; porb_y = AniArray[0][1] - ani_centy; porb_z = AniArray[0][2] - ani_centz; porb_norm = sqrt( porb_x * porb_x + porb_y * porb_y + porb_z * porb_z ); cos_phi = ( porb_x * ( H[0] - ani_centx ) + porb_y * ( H[1] - ani_centy ) + porb_z * ( H[2] - ani_centz ) ) / ( porb_norm * r_H ); shift = chiiso * ( 1 - 3 * cos_phi * cos_phi ) / ( 3 * r_H * r_H * r_H ); } return shift; } /* Define the ring current shift function... */ double ringshift( double H[3], double **RingArray, int nRingAt, double ifactor, char method[3] ) { double **tempRA; double B, plane_x, plane_y, plane_z, plane_norm, plnn_x, plnn_y, plnn_z, p, modulo, r_H, z_to_plane, cos_phi, z_updown, ro_Hjb, xt, yt, zt, vpx, vpy, vpz, r1_x, r1_y, r1_z, r2_x, r2_y, r2_z, s12, r1sq, r2sq, d1, d2; double rc_x = 0, rc_y = 0, rc_z = 0, ring_radius = 0, shift = 0; int i, j, k, at1, at2; for( i = 0; i <= nRingAt - 1; ++i ) { rc_x += RingArray[i][0]; rc_y += RingArray[i][1]; rc_z += RingArray[i][2]; } rc_x /= nRingAt; rc_y /= nRingAt; rc_z /= nRingAt; /* Now we should get the distance from the center of the ring to the proton, find the equation for the normal to the plane, and the distance from the proton to the plane. We will use three of the atoms in the ring. Use atoms 1, 3, and 5... */ r_H = sqrt( ( H[0] - rc_x ) * ( H[0] - rc_x ) + ( H[1] - rc_y ) * ( H[1] - rc_y ) + ( H[2] - rc_z ) * ( H[2] - rc_z ) ); /* Get the normal vector to the plane defined by 3 of the 5 or 6 atoms, get the distance to plane, and get the distance from the center of the ring to the projection of the proton onto the ring's plane... */ plane_x = ( RingArray[1][1] - RingArray[0][1] ) * ( RingArray[2][2] - RingArray[0][2] ) - ( RingArray[1][2] - RingArray[0][2] ) * ( RingArray[2][1] - RingArray[0][1] ); plane_y = ( RingArray[1][2] - RingArray[0][2] ) * ( RingArray[2][0] - RingArray[0][0] ) - ( RingArray[1][0] - RingArray[0][0] ) * ( RingArray[2][2] - RingArray[0][2] ); plane_z = ( RingArray[1][0] - RingArray[0][0] ) * ( RingArray[2][1] - RingArray[0][1] ) - ( RingArray[1][1] - RingArray[0][1] ) * ( RingArray[2][0] - RingArray[0][0] ); plane_norm = sqrt( plane_x * plane_x + plane_y * plane_y + plane_z * plane_z ); plnn_x = plane_x / plane_norm; plnn_y = plane_y / plane_norm; plnn_z = plane_z / plane_norm; z_to_plane = plnn_x * ( H[0] - rc_x ) + plnn_y * ( H[1] - rc_y ) + plnn_z * ( H[2] - rc_z ); ro_Hjb = sqrt( r_H * r_H - z_to_plane * z_to_plane ); /* Decide which method to use... */ if( ! strcmp( method, "PD" ) ) /* Get the contribution to the shift using the Pople point-dipole aproximation. Use the B factor from Moyna et al., JCICS, 1998... */ { B = 22.34 * ifactor; cos_phi = z_to_plane / r_H; shift = B * ( 1 - 3 * cos_phi * cos_phi ) / ( r_H * r_H * r_H ); } else if( ! strcmp( method, "JB" ) ) /* Get the contribution to the shift using the Johnson-Bovey method. There are currents above and below the ring. These rings are at 0.64 (q) / ring_radius from the plane. The B factor here is '[ 1x1014 * enumber * ( echarge )^2 ] / [ 6 * pi * emass * ( speedoflight )^2 ]'... */ { /* Calculate the ring radius... */ for( i = 0; i <= nRingAt - 1; ++i ) { ring_radius += sqrt( ( RingArray[i][0] - rc_x ) * ( RingArray[i][0] - rc_x ) + ( RingArray[i][1] - rc_y ) * ( RingArray[i][1] - rc_y ) + ( RingArray[i][2] - rc_z ) * ( RingArray[i][2] - rc_z ) ); } ring_radius /= nRingAt; B = 4.485758 * ifactor / ring_radius; p = ro_Hjb / ring_radius; for( k = 0; k <= 1; ++k ) { z_updown = ( z_to_plane + 0.64 * ( 2 * k - 1 ) ) / ring_radius; modulo = sqrt( 4.0 * p / ( ( 1 + p ) * ( 1 + p ) + z_updown * z_updown ) ); shift -= B * 1 / sqrt( ( 1 + p ) * ( 1 + p ) + z_updown * z_updown ) * ( ellpk( modulo ) + ( 1 - p * p - z_updown * z_updown ) / ( ( 1 - p ) * ( 1 - p ) + z_updown * z_updown ) * ellpe( modulo ) ); } } else if( ! strcmp( method, "HM" ) ) { /* Store the pointer in a new pointer, otherwise it will get shuffled... */ tempRA = ( double** ) malloc( nRingAt * sizeof( double* ) ); for( i = 0; i <= nRingAt - 1; ++i ) tempRA[i] = ( double* ) malloc( 3 * sizeof( double ) ); for( i = 0; i <= nRingAt - 1; ++i ) { for( j = 0; j <= 2; ++j ) tempRA[i][j] = RingArray[i][j]; } /* Before going any further, do a bubble sort to order the atoms in the ring... */ for( i = 0; i <= nRingAt - 1; ++i ) { j = i + 1; while( j <= nRingAt - 1 ) { if( ( ( tempRA[i][0] - tempRA[j][0] ) * ( tempRA[i][0] - tempRA[j][0] ) + ( tempRA[i][1] - tempRA[j][1] ) * ( tempRA[i][1] - tempRA[j][1] ) + ( tempRA[i][2] - tempRA[j][2] ) * ( tempRA[i][2] - tempRA[j][2] ) ) < 2.50 ) { /* After we find the 'jth' atom next to the 'ith' atom, swap them and modify the 'j' index so that we bail out of the loop and avoid going all the way around the ring... */ xt = tempRA[i + 1][0]; yt = tempRA[i + 1][1]; zt = tempRA[i + 1][2]; tempRA[i + 1][0] = tempRA[j][0]; tempRA[i + 1][1] = tempRA[j][1]; tempRA[i + 1][2] = tempRA[j][2]; tempRA[j][0] = xt; tempRA[j][1] = yt; tempRA[j][2] = zt; j = nRingAt - 1; } ++j; } } vpx = ( tempRA[0][1] - rc_y ) * ( tempRA[2][2] - rc_z ) - ( tempRA[0][2] - rc_z ) * ( tempRA[2][1] - rc_y ); vpy = ( tempRA[0][2] - rc_z ) * ( tempRA[2][0] - rc_x ) - ( tempRA[0][0] - rc_x ) * ( tempRA[2][2] - rc_z ); vpz = ( tempRA[0][0] - rc_x ) * ( tempRA[2][1] - rc_y ) - ( tempRA[0][1] - rc_y ) * ( tempRA[2][0] - rc_x ); B = 2.7274 * ifactor; if( ( vpz * plnn_z + vpy * plnn_y + vpx * plnn_x ) > 0 ) B *= -1.0; for( at1 = 0; at1 <= nRingAt - 1; ++at1 ) { at2 = at1 + 1; if( at2 > nRingAt - 1 ) at2 = 0; r1_x = tempRA[at1][0] - H[0]; r1_y = tempRA[at1][1] - H[1]; r1_z = tempRA[at1][2] - H[2]; r2_x = tempRA[at2][0] - H[0]; r2_y = tempRA[at2][1] - H[1]; r2_z = tempRA[at2][2] - H[2]; s12 = r1_x * ( r2_y * plnn_z - r2_z * plnn_y ) + r1_y * ( r2_z * plnn_x - r2_x * plnn_z ) + r1_z * ( r2_x * plnn_y - r2_y * plnn_x ); r1sq = r1_x * r1_x + r1_y * r1_y + r1_z * r1_z; d1 = sqrt( r1sq ); r2sq = r2_x * r2_x + r2_y * r2_y + r2_z * r2_z; d2 = sqrt( r2sq ); shift += B * s12 * ( 1 / ( r1sq * d1 ) + 1 / ( r2sq * d2 ) ); } for( i = 0; i <= nRingAt - 1; ++i ) free( tempRA[i] ); free( tempRA ); } return shift; } /* Define the electrostatic shift function... */ double electroshift( double H[3], double C[3], double **ElecAtoms, int nEleAt, double *charges, double eleA, double eleB, double diel, char dieldep[3] ) { double CHx, CHy, CHz, XHx, XHy, XHz, rCH, rXH, dotCHXH, shift; int i; shift = 0; CHx = H[0] - C[0]; CHy = H[1] - C[1]; CHz = H[2] - C[2]; rCH = sqrt( CHx * CHx + CHy * CHy + CHz * CHz ); for( i = 0; i <= nEleAt - 1; ++i ) { XHx = ElecAtoms[i][0] - H[0]; XHy = ElecAtoms[i][1] - H[1]; XHz = ElecAtoms[i][2] - H[2]; rXH = sqrt( XHx * XHx + XHy * XHy + XHz * XHz ); dotCHXH = ( CHx * XHx + CHy * XHy + CHz * XHz ) / ( rXH * rCH ); if( ! strcmp( dieldep, "CD" ) ) shift += charges[i] * ( eleA * dotCHXH + eleB ) / ( rXH * rXH * diel ); else if( ! strcmp( dieldep, "RD" ) ) shift += charges[i] * ( eleA * dotCHXH + eleB ) / ( rXH * rXH * rXH * diel ); } return shift; } /* Define the non-bonded potential functions. They use the Tripos 5.2 function... */ double nonbonded( double ***LigArray, int *nLigAt, int nLigands, double **ProtArray, int nProtAt, double cutoff, char useelectro[4], double diel, char dieldep[3], char rmix[4] ) { double energy = 0; double rLAPA, rLALA, aij, aij6; int i, j, k, l; for( i = 0; i <= nLigands - 1; ++i ) { for( j = 0; j <= nLigAt[i] - 1; ++j ) { for( k = 0; k <= nProtAt - 1; ++k ) { rLAPA = sqrt( ( LigArray[i][j][0] - ProtArray[k][0] ) * ( LigArray[i][j][0] - ProtArray[k][0] ) + ( LigArray[i][j][1] - ProtArray[k][1] ) * ( LigArray[i][j][1] - ProtArray[k][1] ) + ( LigArray[i][j][2] - ProtArray[k][2] ) * ( LigArray[i][j][2] - ProtArray[k][2] ) ); if( rLAPA <= cutoff ) { aij = rLAPA / ( LigArray[i][j][4] + ProtArray[k][4] ); aij6 = aij * aij * aij * aij * aij * aij; energy += sqrt( LigArray[i][j][5] * ProtArray[k][5] ) * ( 1.0 / ( aij6 * aij6 ) - 2.0 / aij6 ); if( ! strcmp( useelectro, "YES" ) ) { if( ! strcmp( dieldep, "CD" ) ) energy += 332.17 * LigArray[i][j][3] * ProtArray[k][3] / ( diel * rLAPA ); else if( ! strcmp( dieldep, "RD" ) ) energy += 332.17 * LigArray[i][j][3] * ProtArray[k][3] / ( diel * rLAPA * rLAPA ); } } } } } /* If there is more than one ligand, do the ligand-ligand interactions... */ if( ( nLigands > 1 ) || ( strcmp( rmix, "NO" ) ) ) { for( i = 0; i <= nLigands - 1; ++i ) { for( j = i + 1; j <= nLigands - 1; ++j ) { for( k = 0; k <= nLigAt[i] - 1; ++k ) { for( l = 0; l <= nLigAt[j] - 1; ++l ) { rLALA = sqrt( ( LigArray[i][k][0] - LigArray[j][l][0] ) * ( LigArray[i][k][0] - LigArray[j][l][0] ) + ( LigArray[i][k][1] - LigArray[j][l][1] ) * ( LigArray[i][k][1] - LigArray[j][l][1] ) + ( LigArray[i][k][2] - LigArray[j][l][2] ) * ( LigArray[i][k][2] - LigArray[j][l][2] ) ); if( rLALA <= cutoff ) { aij = rLALA / ( LigArray[i][k][4] + LigArray[j][l][4] ); aij6 = aij * aij * aij * aij * aij * aij; energy += sqrt( LigArray[i][k][5] * LigArray[j][l][5] ) * ( 1.0 / ( aij6 * aij6 ) - 2.0 / aij6 ); if( ! strcmp( useelectro, "YES" ) ) { if( ! strcmp( dieldep, "CD" ) ) energy += 332.17 * LigArray[i][k][3] * LigArray[j][l][3] / ( diel * rLALA ); else if( ! strcmp( dieldep, "RD" ) ) energy += 332.17 * LigArray[i][k][3] * LigArray[j][l][3] / ( diel * rLALA * rLALA ); } } } } } } } if( ! strcmp( rmix, "YES" ) ) energy /= nLigands; return energy; } double slfligene( double ***LigArray, int *nLigAt, int nLigands, int ***LigNBList, int **nLigNBList, double cutoff, char useelectro[4], double diel, char dieldep[3] ) { double energy = 0; double rLALA, aij, aij6; int k, i, j, nbi; for( i = 0; i <= nLigands - 1; ++i ) { for( j = 0; j <= nLigAt[i] - 1; ++j ) { for( k = 0; k <= nLigNBList[i][j] - 1; ++k ) { nbi = LigNBList[i][j][k]; if( nbi > j ) { rLALA = sqrt( ( LigArray[i][j][0] - LigArray[i][nbi][0] ) * ( LigArray[i][j][0] - LigArray[i][nbi][0] ) + ( LigArray[i][j][1] - LigArray[i][nbi][1] ) * ( LigArray[i][j][1] - LigArray[i][nbi][1] ) + ( LigArray[i][j][2] - LigArray[i][nbi][2] ) * ( LigArray[i][j][2] - LigArray[i][nbi][2] ) ); if( rLALA <= cutoff ) { aij = rLALA / ( LigArray[i][j][4] + LigArray[i][nbi][4] ); aij6 = aij * aij * aij * aij * aij * aij; energy += sqrt( LigArray[i][j][5] * LigArray[i][nbi][5] ) * ( 1.0 / ( aij6 * aij6 ) - 2.0 / aij6 ); if( ! strcmp( useelectro, "YES" ) ) { if( ! strcmp( dieldep, "CD" ) ) energy += 332.17 * LigArray[i][j][3] * LigArray[i][nbi][3] / ( diel * rLALA ); else if( ! strcmp( dieldep, "RD" ) ) energy += 332.17 * LigArray[i][j][3] * LigArray[i][nbi][3] / ( diel * rLALA * rLALA ); } } } } } } return energy; } double centpot( int nLigands, double ***LigArray, int *nLigAt, double **initcent, double radius, double force ) { double energy; double dist; double currentcent[3]; int i, j, k; energy = 0; for( i = 0; i <= nLigands - 1; ++i ) { for( j = 0; j <= 2; ++j ) currentcent[j] = 0; for( k = 0; k <= nLigAt[i] - 1; ++k ) { for( j = 0; j <= 2; ++ j ) currentcent[j] += LigArray[i][k][j]; } for( j = 0; j <= 2; ++j ) currentcent[j] /= nLigAt[i]; dist = sqrt( ( initcent[i][0] - currentcent[0] ) * ( initcent[i][0] - currentcent[0] ) + ( initcent[i][1] - currentcent[1] ) * ( initcent[i][1] - currentcent[1] ) + ( initcent[i][2] - currentcent[2] ) * ( initcent[i][2] - currentcent[2] ) ); if( dist > radius ) energy += force * ( dist - radius ) * ( dist - radius ); } return energy; } double noepot( double **nOeAtLig, int *nOeAtLigID, double **nOeAtProt, int *nOeAtProtID, int nnOe, double **nOeRange, double *nOePenalty, char verbose[3] ) { double dist, distdiff, energy, dVAvg, dVSD, dV2Avg, dVAvg2; int i, nViolations; energy = 0; dVAvg = 0; dVSD = 0; dV2Avg = 0; dVAvg2 = 0; nViolations = 0; if( ( ! strcmp( verbose, "YES" ) ) || ( ! strcmp( verbose, "INI" ) ) || ( ! strcmp( verbose, "END" ) ) ) printf( "\n" ); for( i = 0; i <= nnOe - 1; ++i ) { dist = sqrt( ( nOeAtLig[i][0] - nOeAtProt[i][0] ) * ( nOeAtLig[i][0] - nOeAtProt[i][0] ) + ( nOeAtLig[i][1] - nOeAtProt[i][1] ) * ( nOeAtLig[i][1] - nOeAtProt[i][1] ) + ( nOeAtLig[i][2] - nOeAtProt[i][2] ) * ( nOeAtLig[i][2] - nOeAtProt[i][2] ) ); if( ( ! strcmp( verbose, "YES" ) ) || ( ! strcmp( verbose, "INI" ) ) || ( ! strcmp( verbose, "END" ) ) ) printf( "Proton pair [%d - %d] - Experimental range = [%f...%f] - Actual distance = %f\n", nOeAtLigID[i], nOeAtProtID[i], nOeRange[i][0], nOeRange[i][1], dist ); if( dist <= nOeRange[i][0] ) { distdiff = nOeRange[i][0] - dist; energy += nOePenalty[i] * distdiff * distdiff; dVAvg += distdiff; dV2Avg += distdiff * distdiff; ++nViolations; } else if( dist >= nOeRange[i][1] ) { distdiff = dist - nOeRange[i][1]; energy += nOePenalty[i] * distdiff * distdiff; dVAvg += distdiff; dV2Avg += distdiff * distdiff; ++nViolations; } else energy += 0; } dVAvg /= nViolations; dV2Avg /= nViolations; dVAvg2 = dVAvg * dVAvg; dVSD = sqrt( fabs( dV2Avg - dVAvg2 ) ); if( ! strcmp( verbose, "INI" ) ) { printf( "\n" ); printf( "Initial number of nOe violations = %d - Total distance violation = %f +/- %f\n", nViolations, dVAvg, dVSD ); } if( ! strcmp( verbose, "END" ) ) { printf( "\n" ); printf( "Final number of nOe violations = %d - Total distance violation = %f +/- %f\n", nViolations, dVAvg, dVSD ); } return energy; } double gettorpot( int RotBonds, int *nTor, double *pot, int *phase, double ****TorAts ) { double torpot, x21, y21, z21, x32, y32, z32, x43, y43, z43, xt, yt, zt, xu, yu, zu, rtru, cosine; int i, j, k; torpot = 0; for( i = 0; i <= RotBonds - 1; ++i ) { for( j = 0; j <= nTor[i] - 1; ++j ) { x21 = TorAts[i][j][1][0] - TorAts[i][j][0][0]; y21 = TorAts[i][j][1][1] - TorAts[i][j][0][1]; z21 = TorAts[i][j][1][2] - TorAts[i][j][0][2]; x32 = TorAts[i][j][2][0] - TorAts[i][j][1][0]; y32 = TorAts[i][j][2][1] - TorAts[i][j][1][1]; z32 = TorAts[i][j][2][2] - TorAts[i][j][1][2]; x43 = TorAts[i][j][3][0] - TorAts[i][j][2][0]; y43 = TorAts[i][j][3][1] - TorAts[i][j][2][1]; z43 = TorAts[i][j][3][2] - TorAts[i][j][2][2]; xt = y21 * z32 - y32 * z21; yt = z21 * x32 - z32 * x21; zt = x21 * y32 - x32 * y21; xu = y32 * z43 - y43 * z32; yu = z32 * x43 - z43 * x32; zu = x32 * y43 - x43 * y32; rtru = sqrt( ( xt * xt + yt * yt + zt * zt ) * ( xu * xu + yu * yu + zu * zu ) ); cosine = ( xt * xu + yt * yu + zt * zu ) / rtru; torpot += pot[i] * ( 1 + phase[i] / abs( phase[i] ) * cos( phase[i] * acos( cosine ) ) ); } } return ( 0.5 * torpot ); } /* This recursive fucntion finds all the atoms one side of a bond defined as rotatable... */ void markneighbor( int k, short *hitlist, int *nLigBondedList, int **LigBondedList ) { int i; hitlist[k] = 1; for( i = 0; i <= nLigBondedList[k] - 1; ++i ) { if( hitlist[LigBondedList[k][i]] ) continue; else markneighbor( LigBondedList[k][i], hitlist, nLigBondedList, LigBondedList ); } return; } /* This is the bond rotation function... */ void rotbond( int nFragAt, double **FragAt, double theta ) { double cos1, sin1, cos2, sin2, tempx, tempy, tempz, lenght1, lenght2, relpos[3]; int i, j; theta *= ( M_PI / 180.0 ); for( i = 0; i <= 2; ++i ) relpos[i] = FragAt[0][i]; /* First, translate the whole fragment so that the first atoms in the array is at the origin... */ for( i = 0; i <= nFragAt - 1; ++i ) { for( j = 0; j <= 2; ++j ) FragAt[i][j] -= relpos[j]; } lenght1 = sqrt( FragAt[1][0] * FragAt[1][0] + FragAt[1][2] * FragAt[1][2] ); if( lenght1 > 0 ) { cos1 = FragAt[1][2] / lenght1; sin1 = FragAt[1][0] / lenght1; for( i = 0; i <= nFragAt - 1; ++i ) { tempx = FragAt[i][0] * cos1 - FragAt[i][2] * sin1; tempz = FragAt[i][0] * sin1 + FragAt[i][2] * cos1; FragAt[i][0] = tempx; FragAt[i][2] = tempz; } } lenght2 = sqrt( FragAt[1][1] * FragAt[1][1] + FragAt[1][2] * FragAt[1][2] ); cos2 = FragAt[1][2] / lenght2; sin2 = FragAt[1][1] / lenght2; for( i = 0; i <= nFragAt - 1; ++i ) { tempy = FragAt[i][1] * cos2 - FragAt[i][2] * sin2; tempz = FragAt[i][1] * sin2 + FragAt[i][2] * cos2; FragAt[i][1] = tempy; FragAt[i][2] = tempz; } for( i = 0; i <= nFragAt - 1; ++i ) { tempx = FragAt[i][0] * cos( theta ) - FragAt[i][1] * sin( theta ); tempy = FragAt[i][0] * sin( theta ) + FragAt[i][1] * cos( theta ); FragAt[i][0] = tempx; FragAt[i][1] = tempy; } for( i = 0; i <= nFragAt - 1; ++i ) { tempy = FragAt[i][1] * cos2 + FragAt[i][2] * sin2; tempz = FragAt[i][2] * cos2 - FragAt[i][1] * sin2; FragAt[i][1] = tempy; FragAt[i][2] = tempz; } if( lenght1 > 0 ) { for( i = 0; i <= nFragAt - 1; ++i ) { tempx = FragAt[i][0] * cos1 + FragAt[i][2] * sin1; tempz = FragAt[i][2] * cos1 - FragAt[i][0] * sin1; FragAt[i][0] = tempx; FragAt[i][2] = tempz; } } /* After we are done, translate everything back to its original position... */ for( i = 0; i <= nFragAt - 1; ++i ) { for( j = 0; j <= 2; ++j ) FragAt[i][j] += relpos[j]; } return; } /* This is the H-bond shift contribution. Uses iffy equation by Wishart, and does not allow for bifurcated H-bonds... */ double hbshift( double H[3], double N[3], double HBSlope, double HBConst, int nHBAcc, double **HBAccAtoms ) { double NHx, NHy, NHz, DHx, DHy, DHz, rNH, rDH, dotNHDH, shift, besthbradius; bool ishbonded; int i; shift = 0; besthbradius = INFINITY; ishbonded = false; NHx = H[0] - N[0]; NHy = H[1] - N[1]; NHz = H[2] - N[2]; rNH = sqrt( NHx * NHx + NHy * NHy + NHz * NHz ); for( i = 0; i <= nHBAcc - 1; ++i ) { DHx = H[0] - HBAccAtoms[i][0]; DHy = H[1] - HBAccAtoms[i][1]; DHz = H[2] - HBAccAtoms[i][2]; rDH = sqrt( DHx * DHx + DHy * DHy + DHz * DHz ); dotNHDH = ( NHx * DHx + NHy * DHy + NHz * DHz ) / ( rNH * rDH ); if( rDH <= 2.50 ) { if( acos( dotNHDH ) >= 0.75 * M_PI ) { ishbonded = true; if( rDH < besthbradius ) besthbradius = rDH; } } } if( ishbonded ) shift += ( HBSlope / ( besthbradius * besthbradius * besthbradius ) + HBConst ); return shift; } /* Define the word count function (blatantly lifted from the web)... */ int parseString( char *line, char ***args, int delimiter ) { char *buffer; int nargs; buffer = ( char* ) malloc( ( strlen( line ) + 1 ) * sizeof( char ) ); strcpy( buffer, line ); /* Allow for as many arguments as characters in the 'line' string. This will always be more than the actual number of arguments in the line, so we are safe... */ ( *args ) = ( char** ) malloc( strlen( line ) * sizeof( char** ) ); nargs = 0; if( delimiter == 1 ) { ( *args )[nargs++] = strtok( buffer, " \t" ); while( ( ( ( *args )[nargs] = strtok( NULL, " \t\n" ) ) != NULL ) && ( nargs < strlen( line ) ) ) ++nargs; } else if( delimiter == 2 ) { ( *args )[nargs++] = strtok( buffer, "," ); while( ( ( ( *args )[nargs] = strtok( NULL, ",\n" ) ) != NULL ) && ( nargs < strlen( line ) ) ) ++nargs; } else if( delimiter == 3 ) { ( *args )[nargs++] = strtok( buffer, ":" ); while( ( ( ( *args )[nargs] = strtok( NULL, ":\n" ) ) != NULL ) && ( nargs < strlen( line ) ) ) ++nargs; } return nargs; } int writemol( char ft[5], char of[256], int hc, char **hj, int ac, char **rad, int nl, int *nla, int **laID, double ***la, int fc, char **fj, double sdE ) { int i,j,k, ligj, ligk; char **sl; bool tip; FILE *OutFile; OutFile = fopen( of, "w" ); if( ! OutFile ) return 1; /* Write out the SDILICON energy for further ref... */ if( ! strcmp( ft, "MOL2" ) ) fprintf( OutFile, "#\tSDILICON Energy:\t%f Kcal/mol\n", sdE ); else if( ! strcmp( ft, "PDB" ) ) fprintf( OutFile, "REMARK 0 SDILICON Energy = %f Kcal/mol\n", sdE ); /* Now save the stuff that was in the header. Ommit SDILICON refences already there... */ if( hc > 0 ) { for( i = 0; i <= hc; ++i ) { if( ! strstr( hj[i], "SDILICON Energy" ) ) fputs( hj[i], OutFile ); } } /* Now we should save the modified atom info in mol2/pdb format... */ for( i = 0; i <= ac - 1; ++i ) { tip = true; for( j = 0; j <= nl - 1; ++j ) { for( k = 0; k <= nla[j] - 1; ++k ) { if( ( laID[j][k] - 1 ) == i ) { tip = false; ligj = j; ligk = k; } } } if( tip ) fputs( rad[i], OutFile ); else { parseString( rad[i], &sl, 1 ); if( ! strcmp( ft, "MOL2" ) ) { fprintf( OutFile, "%7s %-8s", sl[0], sl[1] ); fprintf( OutFile, "%10.4f%10.4f%10.4f", la[ligj][ligk][0], la[ligj][ligk][1], la[ligj][ligk][2] ); fprintf( OutFile, " %-7s%4s %-8s%10.4f ", sl[5], sl[6], sl[7], atof( sl[8] ) ); if( sl[9] ) fprintf( OutFile, "%s", sl[9] ); fprintf( OutFile, "\n" ); } else if( ! strcmp( ft, "PDB" ) ) { fprintf( OutFile, "%-6s%5s %-4s%4s%6s", sl[0], sl[1], sl[2], sl[3], sl[4] ); fprintf( OutFile, "%12.3f%8.3f%8.3f", la[ligj][ligk][0], la[ligj][ligk][1], la[ligj][ligk][2] ); fprintf( OutFile, "%6s%6s%12s\n", sl[8], sl[9], sl[10] ); } } } /* And finally save the contents of the footer... */ if( fc > 0 ) { for( i = 0; i <= fc - 1; ++i ) fputs( fj[i], OutFile ); } fclose( OutFile ); return 0; } /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Program Starts Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */ main( int argc, char *argv[] ) { /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Variables Go Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */ double ***RingArray, ***AniArray, ***LigArray, **randtrans, **randrot, *randbrot, ****TorArray, ***FragArray, **ElecAtoms, **protons, **initcent, **ProtArray, **protoneigh, **nOeAtLig, **nOeAtProt,**nOeRange, *dexp, *shift, *shftavg, *force, *ifactor, *chiiso, *chiani, *charges, *nOePenalty, *pertfact, *TorPot, *maxbrot, ****transgene, ****transchld, ****rotgene, ****rotchld, ***brotgene, ***brotchld, **gaE, **gascore, **pstart, **pend, **LigInicoord, **HBAccAtoms, *scoresum, *popE, *temperature, *currshft, *shftwavg; int *protonid, *nRingAt, *nAniAt, *nLigAt, **LigAtID, *ProtAtID, *nOeAtLigID, *nOeAtProtID, *nTorsions, *TorPhase, ***LigBondedList, **nLigBondedList, ***LigNB13List, **nLigNB13List, ***LigNBList, **nLigNBList, *nFragAt, *TorAtA, *TorAtB, **rankvec, *ligICidx, *popsize, *diversity, *pii, *pin, *prank; char **ligName, **SplitLine, **HeaderJunk, **RawAtomData, **Protattype, ***Ligattype, **pName, **FooterJunk, **bndtyp, **tempEleString, **protonFullName, **ProtAtName, **ProtAtResName; short *hitlist; double delta, theta, convtrans, convrot, diel, shiftval, rPRPA, ripsE, egrad, nOeE, shiftE, ffE, selfE, torE, gaEsum, mindistance, radiuscp, forcecp, cutoff, upE, downE, iniE, eleA, eleB, econv, AvgDeltaDiff, bondrotstep, brconvavg, AvgDelta2Diff, AvgDeltaDiff2, AvgDeltaSD, jthresh, ProtH[3], ProtHneigh[3], LigAtDist, anchorE, convbrot, geneconv, mingaE, maxgaE, pselect, mutfreq, finalE, genediff, HBSlope, HBConst, pfact, bestislandE, migfreqf, tempgene, tlkelvin, thkelvin, expfact, gasconst, acceptratio, transwin, rotwin, torwin, tpivot, tau, constweight; int i, j, k, l, m, n, x, ncla, maxiter, nRings, iter, neighindex, rseed, pstep, nnOe, nRotBonds, nInicoord, nProtons, nAni, nEleAt, nLigands, nProtAt, ligj, ligk, npstep, eterms, HeaderCount, LigLine, activegavar, AtomCount, FooterCount, TerCount, axes, maxphits, TorArgs, rgdlfctr, xpoint, a1, a2, nHBAcc, nIslands, matingA, matingB, whattomutate, mutatethis, sdlctr, maxgaconf, temprank, keepbest, nargpline, bestisland, dIsland, mIsland, travellers, smallestpop, migcounter, migfreqi, poptotal, sastep, sastpvt, savectr, mdcycles; char verboutput[4], potential[6], method[3], constcent[4], ftype[5], useel[4], attype[8], optmethod[5], pRes[32], pResNo[32], dotrans[4], dorot[4], dielfunc[3], infile[256], jptype[4], outfile[256], sdlfile[256], line[4096], tvlabel[3][3], rvlabel[3][3], calcsdl[256], multconf[4], conffile[256], scorefunc[5], useinicoord[4], dobonds[4], racemate[4], migmethod[5], freqoritvl[5], adjuststep[4], annealfunc[5], runningavg[4]; bool up, oozing, NonBonded, ThisIsProtein, singleorg, noletters, norescore, didmigration, domigration, nonproductive, dosave; FILE *InFile, *SdlFile, *CalcSdl; /* >>>>>>>>>>>>>>>>>>>>>>>>>>>> Functions Declarations Go Here <<<<<<<<<<<<<<<<<<<<<<<<<<< */ /* Define the prototypes for the van der Waals type functions... */ void vdW( char[6], double &, double &, double & ); /* Define the prototypes of the elliptic integral functions... */ double ellint( double, double, double, double ); double ellpk( double ); double ellpe( double ); /* Define the prototypes for the shift functions... */ double anishift( double[3], double **, int, double, double, char[4] ); double electroshift( double[3], double[3], double **, int, double *, double, double, double, char[3] ); double ringshift( double[3], double **, int, double, char[3] ); double hbshift( double[3], double[3], double, double, int, double ** ); /* Define the prototypes for the translation, rotation, and bond rotation functions... */ void translate( double **, int, double, int ); void rotate( double **, int, double, int ); void rotbond( int, double **, double ); /* Define the prototype for the non-bonded and center constraint potential functions... */ double nonbonded( double ***, int *, int, double **, int, double, char[4], double, char[3], char[4] ); double slfligene( double ***, int *, int, int ***, int **, double, char[4], double, char[3] ); double centpot( int, double ***, int*, double **, double, double ); /* Define the prototype for the nOe potential... */ double noepot( double **, int *, double **, int *, int, double **, double *, char[3] ); /* Define the prototype for the torsional potential parameters function... */ void gettorpar( char[8], char[8], char[8], char[8], double &, int & ); /* Define the prototype for the torsional potential function... */ double gettorpot( int, int *, double *, int *, double **** ); /* Define the prototype fragment location recursive functions... */ void markneighbor( int, short *, int *, int ** ); /* Define the prototype for word count function... */ int parseString( char *, char ***, int ); int whitespace = 1, comma = 2, colon = 3; /* Define the prototype of the writemol function... */ int writemol( char[5], char[256], int, char **, int, char **, int, int *, int **, double ***, int, char **, double ); /* Define default variables and do the argument stuff here... */ strcpy( tvlabel[0], "XT" ); strcpy( tvlabel[1], "YT" ); strcpy( tvlabel[2], "ZT" ); strcpy( rvlabel[0], "XR" ); strcpy( rvlabel[1], "YR" ); strcpy( rvlabel[2], "ZR" ); strcpy( infile, "" ); strcpy( outfile, "" ); strcpy( sdlfile, "" ); strcpy( calcsdl, "" ); delta = 1.0; theta = 30.0; bondrotstep = 30.0; convtrans = 0.01; convrot = 0.05; convbrot = 0.05; transwin = 1.0; rotwin = 1.0; torwin = 1.0; maxiter = 100; diel = 1.0; cutoff = 8.0; rseed = time( NULL ); maxphits = 0; egrad = 0.0; jthresh = 0.1; nIslands = 1; popsize = ( int* ) malloc( nIslands * sizeof( int ) ); popsize[0] = 500; poptotal = 500; mutfreq = 0.0; maxgaconf = 0; keepbest = 0; travellers = 5; migfreqi = 10; migfreqf = 0.2; tlkelvin = 298.0; thkelvin = 2000.0; mdcycles = 0; gasconst = 0.001987; acceptratio = 0.5; tpivot = 0.5; tau = 1.0; constweight = 1.0; strcpy( optmethod, "RGDL" ); strcpy( ftype, "MOL2" ); strcpy( verboutput, "NO" ); strcpy( potential, "FFNMR" ); strcpy( method, "HM" ); strcpy( constcent, "NO" ); strcpy( dotrans, "YES" ); strcpy( dorot, "YES" ); strcpy( dobonds, "YES" ); strcpy( dielfunc, "RD" ); strcpy( useel, "NO" ); strcpy( multconf, "NO" ); strcpy( scorefunc, "RANK" ); strcpy( useinicoord, "NO" ); strcpy( jptype, "ALL" ); strcpy( racemate, "NO" ); strcpy( migmethod, "SWAP" ); strcpy( freqoritvl, "ITVL" ); strcpy( adjuststep, "NO" ); strcpy( annealfunc, "EXP" ); strcpy( runningavg, "NO" ); /* We should now check for the correct number of arguments. See if the quotient is an integer.. */ if( ! ( ( argc - 1 ) % 2 == 0 ) ) { printf( "Incorrect number of arguments. Check and try again...\n" ); return 0; } ncla = ( argc - 1 ) / 2; /* The following are the valid keywords that can be passed to sdilicon, together with a short explanation of the expceted values for each control flag: General Keywords ---------------- -if - File to use as input (whith directory path) -of - File to use as output (whith directory path) -cf - The name of the control file (whith directory path) -ft - The file type - Sybyl mol2 - Brookhaven PDB -it - Initial step for translations -ir - Initial step for rotations -br - Initial step for bond rotations -ct - Convergence criterion for translations -cr - Convergence criterion for rotations -cb - Convergence criterion for bond rotations -mi - Maximum number of iterations -ff - Potential to use - NMR potential - NMR and FF potential - FF potential only -rc - Ring-current model to use - Haigh-Mallion - Johnson-Bovey - Pople point-dipole -vo - for verbose updates of the shifts after each cycle -cc - Entering maintains ligand center around sphere of radius with force -dt - will turn off translations during optimization -dr - will turn off rotations during optimization -db - will turn of bond rotations during optimization -dc - Dielectric constant for all electrostatic calculations -df - Dielectric function - Constant dielectric - Distance-dependent dielectric -co - The non-bonded cutoff in Angstroms for FF calculations -ue - to use compute non-bonded electrostatics interactions -om - Optimization method - Ragdoll optimization - Random incremental pulse search - Metropolis molecular dynamics - Metropolis simulated annealing - Simple genetic algorithm - Print calculated shifts only -mc - to save conformers as optimization progresses -rs - Random-seed for RIPS and GENA optimizations. Default is current time in seconds -ic - If initial coordinates for the ligands are present in the sdilicon control file, move the ligands to those coordinates before starting the optimization -rm - Considers ligands as independent, and computes the shift as the average from the perturbations of each ligand. -cw - Weigh all the chemical shift constraint forces by the factor entered RIPS/MCMD/MCSA Specific Keywords -------------------------------- -ph - Maximum number of RIPS optimization productive hits -eg - Energy gradient to use for RIPS optimization convergence -tl - Temperature, in K, to use for Metropolis dynamics or as low temperature in Metropolis simulated annealing -th - Temperature, in K, to use as high temperature in Metropolis simulated annealing -nc - The number of cycles in the MCSA run, or the number of cycles between saving structures to file in MCMD -as - to adjust the step of the translations, rotations, and torsions to have a ratio of productive versus non-productive hits, specified with '-ar' -ar - The ratio between productive and total steps to use for adjusting the step -af - The annealing function to use - Exponential annealing - Linear annealing - Square-step annealing -ra - to compute running averages of the chemical shifts during MCMD or MCSA runs. Additionally, will use the weighed running average in the penalty function following the method of Torda -tu - When using the weighed running average in the penalty functions, this will be the factor to use in the history 'memory'. The value representes the decay factor, in iteration units. Only works with MSMD or MCSA calculations if running averages are turned to -tp - In MCSA, the point in the cycle were the temperature starts decreasing Genetic Algorithm Specific Keywords ----------------------------------- -ps - Population size. Use a comma-separated list to specify number of islands and their populations -kb - The number of best scoring population members to to pass straight into the next generation. In other words, elitism -mf - Mutation frequency for genes in the off-spring. <0.0> will give no mutations and <1.0> will make all the genes to be mutated -nt - The number of individuals to use in migrations between islands -mm - Migration method - Swap best conformers from two islands - Replace worst conformers of one island with best of another -mp - The probability of having a migration. <0.0> will result in no migrations between islands, and <1.0> will force migrations between islands after every generation -nc - Maximum number of conformer to save from the final population -sf - The function for the fitness evaluation - Use only energies - Use the log10 of the energies - Use the logn of the energies - Uses the square root of the energies - Use the ranking for scoring -tw - The window to use in mutations of translation genes. Between 0 and 1 -rw - The window to use in mutations of rotation genes. Between 0 and 1 -bw - The window to use in mutations of bond torsion genes. Between 0 and 1 Calculated J-Surface Specific Keywords ------------------------------------- -jc - Print out a sdilicon control file with calculated shifts for all protein protons to generate theoretical j-surfaces -jt - Threshold in ppm for calculated shifts in calculated sdilicon control file -jp - Specifies the proton type for which calculated shifts will be included in the calculated sdilicon control file New flags can be easily added to the code later on... */ /* First, find which optimization method we are using. Depending on that, a number of defaults are modified... */ for( i = 1; i <= ncla; ++i ) { for( j = 0; j <= strlen( argv[2*i - 1] ) - 1; ++j ) argv[2*i - 1][j] = tolower( argv[2*i - 1][j] ); if( strcmp( argv[2*i - 1], "-if" ) && strcmp( argv[2*i - 1], "-cf" ) && strcmp( argv[2*i - 1], "-of" ) && strcmp( argv[2*i - 1], "-jc" ) ) { for( j = 0; j <= strlen( argv[2*i] ) - 1; ++j ) argv[2*i][j] = toupper( argv[2*i][j] ); } if( ! strcmp( argv[2*i - 1], "-om" ) ) strcpy( optmethod, argv[2*i] ); } /* Now re-set defaults... */ if( ! strcmp( optmethod, "MCMD" ) ) { maxphits = 1000; strcpy( adjuststep, "YES" ); } else if( ! strcmp( optmethod, "RIPS" ) ) { maxiter = 10000; } else if( ! strcmp( optmethod, "MCSA" ) ) { mdcycles = 10; maxphits = 1000; tlkelvin = 200; strcpy( adjuststep, "YES" ); } else if( ! strcmp( optmethod, "GENA" ) ) { delta = 2.0; theta = 180.0; bondrotstep = 180.0; } else if( ( strcmp( optmethod, "RGDL" ) ) && ( strcmp( optmethod, "RPRT" ) ) ) { printf( "Unknown optimization method %s. Check and try again...\n", optmethod ); return 0; } for( i = 1; i <= ncla; ++i ) { for( j = 0; j <= strlen( argv[2*i - 1] ) - 1; ++j ) argv[2*i - 1][j] = tolower( argv[2*i - 1][j] ); if( strcmp( argv[2*i - 1], "-if" ) && strcmp( argv[2*i - 1], "-cf" ) && strcmp( argv[2*i - 1], "-of" ) && strcmp( argv[2*i - 1], "-jc" ) ) { for( j = 0; j <= strlen( argv[2*i] ) - 1; ++j ) argv[2*i][j] = toupper( argv[2*i][j] ); } if( ! strcmp( argv[2*i - 1], "-if" ) ) strcpy( infile, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-of" ) ) strcpy( outfile, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-cf" ) ) strcpy( sdlfile, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-it" ) ) delta = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-ir" ) ) theta = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-br" ) ) bondrotstep = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-ct" ) ) convtrans = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-cr" ) ) convrot = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-cb" ) ) convbrot = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-mi" ) ) { if( ( ! strcmp( optmethod, "MCMD" ) ) || ( ! strcmp( optmethod, "MCSA" ) ) ) maxphits = atoi( argv[2*i] ); else maxiter = atoi( argv[2*i] ); } else if( ! strcmp( argv[2*i - 1], "-ra" ) ) strcpy( runningavg, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-mc" ) ) strcpy( multconf, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-ff" ) ) strcpy( potential, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-co" ) ) cutoff = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-rc" ) ) strcpy( method, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-vo" ) ) strcpy( verboutput, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-dt" ) ) strcpy( dotrans, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-dr" ) ) strcpy( dorot, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-dc" ) ) diel = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-df" ) ) strcpy( dielfunc, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-ft" ) ) strcpy( ftype, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-mm" ) ) strcpy( migmethod, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-rs" ) ) rseed = atoi( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-ps" ) ) { nIslands = parseString( argv[2*i], &SplitLine, comma ); popsize = ( int* ) realloc( popsize, nIslands * sizeof( int ) ); poptotal = 0; for( j = 0; j <= nIslands - 1; ++j ) { popsize[j] = atoi( SplitLine[j] ); poptotal += popsize[j]; } } else if( ! strcmp( argv[2*i - 1], "-ph" ) ) { if( ( ! strcmp( optmethod, "MCMD" ) ) || ( ! strcmp( optmethod, "MCSA" ) ) ) continue; else maxphits = atoi( argv[2*i] ); } else if( ! strcmp( argv[2*i - 1], "-eg" ) ) egrad = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-ue" ) ) strcpy( useel, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-jc" ) ) strcpy( calcsdl, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-jt" ) ) jthresh = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-jp" ) ) strcpy( jptype, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-kb" ) ) keepbest = atoi( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-mf" ) ) mutfreq = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-tw" ) ) transwin = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-rw" ) ) rotwin = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-bw" ) ) torwin = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-nc" ) ) { if( ! strcmp( optmethod, "GENA" ) ) maxgaconf = atoi( argv[2*i] ); else if( ( ! strcmp( optmethod, "MCMD" ) ) || ( ! strcmp( optmethod, "MCSA" ) ) ) mdcycles = atoi( argv[2*i] ); } else if( ! strcmp( argv[2*i - 1], "-nt" ) ) travellers = atoi( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-mp" ) ) { if( atof( argv[2*i] ) > 1 ) { migfreqi = atoi( argv[2*i] ); strcpy( freqoritvl, "ITVL" ); } else { migfreqf = atof( argv[2*i] ); strcpy( freqoritvl, "FREQ" ); } } else if( ! strcmp( argv[2*i - 1], "-sf" ) ) strcpy( scorefunc, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-ic" ) ) strcpy( useinicoord, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-db" ) ) strcpy( dobonds, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-rm" ) ) strcpy( racemate, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-cc" ) ) { strcpy( constcent, "YES" ); parseString( argv[2*i], &SplitLine, comma ); radiuscp = atof( SplitLine[0] ); forcecp = atof( SplitLine[1] ); } else if( ! strcmp( argv[2*i - 1], "-tl" ) ) tlkelvin = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-th" ) ) thkelvin = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-ar" ) ) acceptratio = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-as" ) ) strcpy( adjuststep, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-af" ) ) strcpy( annealfunc, argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-tp" ) ) tpivot = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-tu" ) ) tau = atof( argv[2*i] ); else if( ! strcmp( argv[2*i - 1], "-cw" ) ) constweight = atof( argv[2*i] ); else if( strcmp( argv[2*i - 1], "-om" ) ) { printf( "Unknown command line argument %s. Check and try again...\n", argv[2*i - 1] ); return 0; } } /* Bail out if the user forgot the input or control files... */ if( ( ! strcmp( infile, "" ) ) || ( ! strcmp( sdlfile, "" ) ) ) { printf( "Input molecule file and sdl files are required. Check and try again...\n" ); return 0; } /* If no output file was specified, use the input file to store the result... */ if( ! strcmp( outfile, "" ) ) strcpy( outfile, infile ); /* Now we should have all the arguments. Check if the mol2/pdb and sdl files exist. If they do, open the mol2/pdb file first and count the number of lines of header, atoms, and footer. Then open the sdl file and get the number of ligands, protons, rings, anisotropic groups. Don't worry about charged groups - We'll get them later if they are there... */ HeaderCount = 0; AtomCount = 0; FooterCount = 0; TerCount = 0; InFile = fopen( infile, "r" ); if( ! InFile ) { printf( "Input file %s not found or corrupt. Check and try again...\n", infile ); return 0; } if( ! strcmp( ftype, "MOL2" ) ) { while( ! strstr( fgets( line, sizeof( line ), InFile ), "@ATOM" ) ) ++HeaderCount; while( ! strstr( fgets( line, sizeof( line ), InFile ), "@BOND" ) ) ++AtomCount; while( ! feof( InFile ) ) { fgets( line, sizeof( line ), InFile ); ++FooterCount; } } else if( ! strcmp( ftype, "PDB" ) ) { while( ! strstr( fgets( line, sizeof( line ), InFile ), "ATOM" ) ) ++HeaderCount; if( HeaderCount > 0 ) --HeaderCount; while( ! ( strstr( line, "MASTER" ) || strstr( line, "CONECT" ) || strstr( line, "END" ) ) ) { ++AtomCount; if( strstr( line, "TER " ) ) ++TerCount; fgets( line, sizeof( line ), InFile ); } while( ! feof( InFile ) ) { ++FooterCount; fgets( line, sizeof( line ), InFile ); } } /* Now we have this. Allocate memory for the header, footer, and raw atom data arrays... */ HeaderJunk = ( char** ) malloc( ( HeaderCount + 1 ) * sizeof( char* ) ); for( j = 0; j <= HeaderCount; ++j ) HeaderJunk[j] = ( char* ) malloc( 1024 * sizeof( char ) ); RawAtomData = ( char** ) malloc( AtomCount * sizeof( char* ) ); for( j = 0; j <= AtomCount - 1; ++j ) RawAtomData[j] = ( char* ) malloc( 1024 * sizeof( char ) ); FooterJunk = ( char** ) malloc( FooterCount * sizeof( char* ) ); for( j = 0; j <= FooterCount - 1; ++j ) FooterJunk[j] = ( char* ) malloc( 1024 * sizeof( char ) ); /* Done. Now rewind the file and read all the stuff in the mol2/pdb file. The reading is the same for both... */ rewind( InFile ); if( HeaderCount > 0 ) { for( i = 0; i <= HeaderCount; ++i ) { fgets( line, sizeof( line ), InFile ); strcpy( HeaderJunk[i], line ); } } for( i = 0; i <= AtomCount - 1; ++i ) { fgets( line, sizeof( line ), InFile ); strcpy( RawAtomData[i], line ); } if( FooterCount > 0 ) { for( i = 0; i <= FooterCount - 1; ++i ) { fgets( line, sizeof( line ), InFile ); strcpy( FooterJunk[i], line ); } } fclose( InFile ); /* Now collect info from the sdl file... */ nLigands = 0; sdlctr = 0; nRotBonds = 0; nProtAt = 0; nProtons = 0; nnOe = 0; nRings = 0; nAni = 0; nEleAt = 0; nInicoord = 0; nHBAcc = 0; SdlFile = fopen( sdlfile, "r" ); if( ! SdlFile ) { printf( "Input sdl file %s not found or corrupt. Check and try again...\n", sdlfile ); return 0; } fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! feof( SdlFile ) ) { if( strstr( line, "@LIGAND" ) ) { LigLine = sdlctr; fgets( line, sizeof( line ), SdlFile ); while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { ++nLigands; fgets( line, sizeof( line ), SdlFile ); } } else if( strstr( line, "@ROTBONDS" ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { ++nRotBonds; fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } else if( strstr( line, "@PROTONS" ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { ++nProtons; fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } else if( strstr( line, "@NOE" ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { ++nnOe; fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } else if( strstr( line, "@AROMRINGS" ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { ++nRings; fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } else if( strstr( line, "@ANIGROUPS" ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { ++nAni; fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } else if( strstr( line, "@INICOORDS" ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { ++nInicoord; fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } else if( strstr( line, "@ELECTRO" ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } else if( strstr( line, "@HBACCEPT" ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; while( ! ( strstr( line, "@" ) || feof( SdlFile ) ) ) { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } else { fgets( line, sizeof( line ), SdlFile ); ++sdlctr; } } /* Done, the dimenssions of all the raw data is in. Now we have to arrange it into the appropriate arrays. We will do this at on the fly as we read the information form the sdl file. We also do the allocation here to use memory wisely. In particular, to allocate ligands, rotatable bonds, 5- and 6-membered rings, and 2- and 3-atom anisotropic groups with the right ammount of memory. In the case of the ligand we also have to keep map of the atom id to sdilicon atom index so that we can save it after the calculation is over... */ rewind( SdlFile ); /* The very first thing to do is to go fish the ligands out, as many of the other arrays depend on the ligand atoms. This way the items in the sdl file can be in any order... */ for( sdlctr = 0; sdlctr <= LigLine - 1; ++sdlctr ) fgets( line, sizeof( line ), SdlFile ); nLigAt = ( int* ) malloc( nLigands * sizeof( int ) ); LigArray = ( double*** ) malloc( nLigands * sizeof( double** ) ); ligName = ( char** ) malloc( nLigands * sizeof( char* ) ); LigAtID = ( int** ) malloc( nLigands * sizeof( int* ) ); Ligattype = ( char*** ) malloc( nLigands * sizeof( char** ) ); for( i = 0; i <= nLigands - 1; ++i ) { nLigAt[i] = parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ) - 1; LigArray[i] = ( double** ) malloc( nLigAt[i] * sizeof( double* ) ); ligName[i] = ( char* ) malloc( 256 * sizeof( char ) ); LigAtID[i] = ( int* ) malloc( nLigAt[i] * sizeof( int ) ); Ligattype[i] = ( char** ) malloc( nLigAt[i] * sizeof( char* ) ); strcpy( ligName[i], SplitLine[0] ); nProtAt += nLigAt[i]; for( j = 0; j <= nLigAt[i] - 1; ++j ) { LigArray[i][j] = ( double* ) malloc( 7 * sizeof( double ) ); Ligattype[i][j] = ( char* ) malloc( 8 * sizeof( char ) ); LigAtID[i][j] = atoi( SplitLine[j + 1] ); if( ! strcmp( ftype, "MOL2" ) ) sscanf( RawAtomData[LigAtID[i][j] - 1], "%*d%*s%lf%lf%lf%s%*d%*s%lf", &LigArray[i][j][0], &LigArray[i][j][1], &LigArray[i][j][2], Ligattype[i][j], &LigArray[i][j][3] ); else if( ! strcmp( ftype, "PDB" ) ) { sscanf( RawAtomData[LigAtID[i][j] - 1], "%*s%*d%*s%*s%*d%lf%lf%lf%*lf%*lf%s", &LigArray[i][j][0], &LigArray[i][j][1], &LigArray[i][j][2], Ligattype[i][j] ); LigArray[i][j][3] = 0; // PDB's have no charge... } vdW( Ligattype[i][j], LigArray[i][j][4], LigArray[i][j][5], LigArray[i][j][6] ); } } nProtAt = AtomCount - nProtAt - TerCount; /* With this info, we can get the bonded/non-bonded list for the ligands. However, do it only if there are rotatable bonds defined in the sdl file... */ if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { LigBondedList = ( int*** ) malloc( nLigands * sizeof( int** ) ); nLigBondedList = ( int** ) malloc( nLigands * sizeof( int* ) ); for( i = 0; i <= nLigands - 1; ++i ) { LigBondedList[i] = ( int** ) malloc( nLigAt[i] * sizeof( int* ) ); nLigBondedList[i] = ( int* ) malloc( nLigAt[i] * sizeof( int ) ); for( j = 0; j <= nLigAt[i] - 1; ++j ) { l = 0; LigBondedList[i][j] = ( int* ) malloc( 1 * sizeof( int ) ); // Start with 1, then realloc memory... for( k = 0; k <= nLigAt[i] - 1; ++k ) { if( k != j ) { LigAtDist = sqrt( ( LigArray[i][j][0] - LigArray[i][k][0] ) * ( LigArray[i][j][0] - LigArray[i][k][0] ) + ( LigArray[i][j][1] - LigArray[i][k][1] ) * ( LigArray[i][j][1] - LigArray[i][k][1] ) + ( LigArray[i][j][2] - LigArray[i][k][2] ) * ( LigArray[i][j][2] - LigArray[i][k][2] ) ); if( LigAtDist < 1.2 * ( LigArray[i][j][6] + LigArray[i][k][6] ) ) { LigBondedList[i][j][l] = k; LigBondedList[i][j] = ( int* ) realloc( LigBondedList[i][j], ( l + 2 ) * sizeof( int ) ); ++l; } } } nLigBondedList[i][j] = l; } } /* Get the 1-3 bonded list... */ LigNB13List = ( int*** ) malloc( nLigands * sizeof( int** ) ); nLigNB13List = ( int** ) malloc( nLigands * sizeof( int* ) ); for( i = 0; i <= nLigands - 1; ++i ) { LigNB13List[i] = ( int** ) malloc( nLigAt[i] * sizeof( int* ) ); nLigNB13List[i] = ( int* ) malloc( nLigAt[i] * sizeof( int ) ); for( j = 0; j <= nLigAt[i] - 1; ++j ) { l = 0; LigNB13List[i][j] = ( int* ) malloc( 1 * sizeof( int ) ); // Start with 1, then realloc memory... for( k = 0; k <= nLigBondedList[i][j] - 1; ++k ) { for( m = 0; m <= nLigBondedList[i][LigBondedList[i][j][k]] - 1; ++m ) { if( LigBondedList[i][LigBondedList[i][j][k]][m] != j ) { LigNB13List[i][j][l] = LigBondedList[i][LigBondedList[i][j][k]][m]; LigNB13List[i][j] = ( int* ) realloc( LigNB13List[i][j], ( l + 2 ) * sizeof( int ) ); ++l; } } } nLigNB13List[i][j] = l; } } /* Now we have the 1-2 and 1-3 bonded lists. Get the non-bonded list... */ LigNBList = ( int*** ) malloc( nLigands * sizeof( int** ) ); nLigNBList = ( int** ) malloc( nLigands * sizeof( int* ) ); for( i = 0; i <= nLigands - 1; ++i ) { LigNBList[i] = ( int** ) malloc( nLigAt[i] * sizeof( int* ) ); nLigNBList[i] = ( int* ) malloc( nLigAt[i] * sizeof( int ) ); for( j = 0; j <= nLigAt[i] - 1; ++j ) { m = 0; LigNBList[i][j] = ( int* ) malloc( 1 * sizeof( int ) ); // Start with 1, then realloc memory... for( k = 0; k <= nLigAt[i] - 1; ++k ) { NonBonded = true; if( k != j ) { for( l = 0; l <= nLigBondedList[i][j] - 1; ++l ) { if( k == LigBondedList[i][j][l] ) NonBonded = false; } if( NonBonded ) { for( l = 0; l <= nLigNB13List[i][j] - 1; ++l ) { if( k == LigNB13List[i][j][l] ) NonBonded = false; } } if( NonBonded ) { LigNBList[i][j][m] = k; LigNBList[i][j] = ( int* ) realloc( LigNBList[i][j], ( m + 2 ) * sizeof( int ) ); ++m; } } } nLigNBList[i][j] = m; } } /* Remember to 'free' the 13 array - It is not used for anything else... */ for( i = 0; i <= nLigands - 1; ++i ) { for( j = 0; j <= nLigAt[i] - 1; ++j ) free( LigNB13List[i][j] ); free( LigNB13List[i] ); } free( LigNB13List ); } /* Also, now that we have all the ligands, fish out the protein atoms. Allocate memory for it. We also need a counter to conunt the number of hits and index them. Use 'l'.. */ ProtArray = ( double ** ) malloc( nProtAt * sizeof( double * ) ); Protattype = ( char ** ) malloc( nProtAt * sizeof( char * ) ); ProtAtName = ( char ** ) malloc( nProtAt * sizeof( char * ) ); ProtAtResName = ( char ** ) malloc( nProtAt * sizeof( char * ) ); for( i = 0; i <= nProtAt - 1; ++i ) { ProtArray[i] = ( double * ) malloc( 7 * sizeof( double ) ); Protattype[i] = ( char * ) malloc( 8 * sizeof( char ) ); ProtAtName[i] = ( char * ) malloc( 16 * sizeof( char ) ); ProtAtResName[i] = ( char * ) malloc( 32 * sizeof( char ) ); } ProtAtID = ( int* ) malloc( nProtAt * sizeof( int ) ); l = 0; for( i = 0; i <= AtomCount - 1; ++i ) { ThisIsProtein = true; for( j = 0; j <= nLigands - 1; ++j ) { for( k = 0; k <= nLigAt[j] - 1; ++k ) { if( ( LigAtID[j][k] - 1 ) == i ) { ThisIsProtein = false; break; } } if( ! ThisIsProtein ) break; } if( ThisIsProtein ) { if( ! strcmp( ftype, "MOL2" ) ) { sscanf( RawAtomData[i], "%d%s%lf%lf%lf%s%*d%s%lf", &ProtAtID[l], ProtAtName[l], &ProtArray[l][0], &ProtArray[l][1], &ProtArray[l][2], Protattype[l], ProtAtResName[l], &ProtArray[l][3] ); vdW( Protattype[l], ProtArray[l][4], ProtArray[l][5], ProtArray[l][6] ); ++l; } else if( ( ! strcmp( ftype, "PDB" ) ) && ( ! strstr( RawAtomData[i], "TER " ) ) ) { sscanf( RawAtomData[i], "%*s%d%s%s%s%lf%lf%lf%*lf%*lf%s", &ProtAtID[l], &ProtAtName[l], ProtAtResName[l], pResNo, &ProtArray[l][0], &ProtArray[l][1], &ProtArray[l][2], Protattype[l] ); strcat( ProtAtResName[l], pResNo ); ProtArray[l][3] = 0; // PDB's have no charge... vdW( Protattype[l], ProtArray[l][4], ProtArray[l][5], ProtArray[l][6] ); ++l; } } } /* After we read the ligand, we can go on with everything else. Rewind the file... */ rewind( SdlFile ); fgets( line, sizeof( line ), SdlFile ); while( ! feof( SdlFile ) ) { /* See if there are rotatable bonds defined, and create all the arrays needed for the calculation of torsional potential and bond rotations... */ if( strstr( line, "@ROTBONDS" ) ) { if( ! strcmp( dobonds, "YES" ) ) { nTorsions = ( int* ) malloc( nRotBonds * sizeof( int ) ); TorPot = ( double* ) malloc( nRotBonds * sizeof( double ) ); TorAtA = ( int* ) malloc( nRotBonds * sizeof( int ) ); TorAtB = ( int* ) malloc( nRotBonds * sizeof( int ) ); TorPhase = ( int* ) malloc( nRotBonds * sizeof( int ) ); maxbrot = ( double* ) malloc( nRotBonds * sizeof( double ) ); TorArray = ( double**** ) malloc( nRotBonds * sizeof( double*** ) ); FragArray = ( double*** ) malloc( nRotBonds * sizeof( double** ) ); nFragAt = ( int* ) malloc( nRotBonds * sizeof( int ) ); for( i = 0; i <= nRotBonds - 1; ++i ) { TorArgs = parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ); TorAtA[i] = atoi( SplitLine[0] ); TorAtB[i] = atoi( SplitLine[1] ); nTorsions[i] = 0; /* Find all the atoms involved in the torsion... */ for( j = 0; j <= nLigands - 1; ++j ) { for( k = 0; k <= nLigAt[j] - 1; ++k ) { if( ( LigAtID[j][k] ) == TorAtA[i] ) { nTorsions[i] = ( nLigBondedList[j][k] - 1 ); for( l = 0; l <= nLigAt[j] - 1; ++l ) { if( ( LigAtID[j][l] ) == TorAtB[i] ) { nTorsions[i] *= ( nLigBondedList[j][l] - 1 ); TorArray[i] = ( double*** ) malloc( nTorsions[i] * sizeof( double** ) ); for( m = 0; m <= nTorsions[i] - 1; ++m ) { TorArray[i][m] = ( double** ) malloc( 4 * sizeof( double* ) ); for( n = 0; n <= 3; ++n ) TorArray[i][m][n] = ( double* ) malloc( 3 * sizeof( double ) ); } x = 0; for( m = 0; m <= nLigBondedList[j][l] - 1; ++m ) { if( LigBondedList[j][l][m] != k ) { for( n = 0; n <= nLigBondedList[j][k] - 1; ++n ) { if( LigBondedList[j][k][n] != l ) { TorArray[i][x][0] = LigArray[j][LigBondedList[j][k][n]]; TorArray[i][x][1] = LigArray[j][k]; TorArray[i][x][2] = LigArray[j][l]; TorArray[i][x][3] = LigArray[j][LigBondedList[j][l][m]]; gettorpar( Ligattype[j][LigBondedList[j][k][n]], Ligattype[j][k], Ligattype[j][l], Ligattype[j][LigBondedList[j][l][m]], TorPot[i], TorPhase[i] ); ++x; } } } } maxbrot[i] = bondrotstep; if( TorArgs == 3 ) maxbrot[i] = atof( SplitLine[2] ); else if( TorArgs == 4 ) { TorPot[i] = atof( SplitLine[2] ); TorPhase[i] = atoi( SplitLine[3] ); } else if( TorArgs == 5 ) { TorPot[i] = atof( SplitLine[2] ); TorPhase[i] = atoi( SplitLine[3] ); maxbrot[i] = atof( SplitLine[4] ); } /* Now find the fragments to each side of the bond. Uses a recursive function to do it... */ hitlist = ( short* ) calloc( nLigAt[j], sizeof( short ) ); hitlist[k] = 1; nFragAt[i] = 0; markneighbor( l, hitlist, nLigBondedList[j], LigBondedList[j] ); for( m = 0; m <= nLigAt[j] - 1; ++m ) nFragAt[i] += hitlist[m]; FragArray[i] = ( double** ) malloc( nFragAt[i] * sizeof( double* ) ); for( m = 0; m <= nFragAt[i] - 1; ++m ) FragArray[i][m] = ( double* ) malloc( 3 * sizeof( double ) ); FragArray[i][0] = LigArray[j][k]; FragArray[i][1] = LigArray[j][l]; x = 2; for( m = 0; m <= nLigAt[j] - 1; ++m ) { if( ( m != k ) && ( m != l ) && ( hitlist[m] ) ) { FragArray[i][x] = LigArray[j][m]; ++x; } } free( hitlist ); break; } } break; } } } } } else { for( i = 0; i <= nRotBonds - 1; ++i ) fgets( line, sizeof( line ), SdlFile ); } fgets( line, sizeof( line ), SdlFile ); } /* Once we get protons, allocate the memory and read in the data. Fish the coordinates from the raw atom data... */ /* ratios: H-alkyl = 1.00, H = 1.38, N|C|CB = 0.195 , CA = 0.29 */ else if( strstr( line, "@PROTONS" ) ) { shift = ( double* ) malloc( nProtons * sizeof( double ) ); dexp = ( double* ) malloc( nProtons * sizeof( double ) ); force = ( double* ) malloc( nProtons * sizeof( double ) ); pertfact = ( double* ) malloc( nProtons * sizeof( double ) ); protonid = ( int* ) malloc( nProtons * sizeof( int ) ); protons = ( double** ) malloc( nProtons * sizeof( double* ) ); pName = ( char** ) malloc( nProtons * sizeof( char* ) ); protonFullName = ( char** ) malloc( nProtons * sizeof( char* ) ); for( i = 0; i <= nProtons - 1; ++i ) { protons[i] = ( double * ) malloc( 3 * sizeof( double ) ); pName[i] = ( char* ) malloc( 16 * sizeof( char ) ); protonFullName[i] = ( char* ) malloc( 32 * sizeof( char ) ); nargpline = parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ); protonid[i] = atoi( SplitLine[0] ); dexp[i] = atof( SplitLine[1] ); force[i] = constweight * atof( SplitLine[2] ); if( ! strcmp( ftype, "MOL2" ) ) sscanf( RawAtomData[protonid[i] - 1], "%*d%s%lf%lf%lf%*s%*d%s", pName[i], &protons[i][0], &protons[i][1], &protons[i][2], pRes ); else if( ! strcmp( ftype, "PDB" ) ) { sscanf( RawAtomData[protonid[i] - 1], "%*s%*d%s%s%s%lf%lf%lf", pName[i], pRes, pResNo, &protons[i][0], &protons[i][1], &protons[i][2] ); strcat( pRes, pResNo ); } strcpy( protonFullName[i], strcat( strcat( pRes, "." ), pName[i] ) ); if( nargpline > 3 ) { if( ! strcmp( SplitLine[3], "*" ) ) { if( ! strcmp( pName[i], "H" ) ) pertfact[i] = 1.376; else if( ( ! strcmp( pName[i], "N" ) ) || ( ! strcmp( pName[i], "C" ) ) || ( ! strcmp( pName[i], "CB" ) ) ) pertfact[i] = 0.195; else if( ! strcmp( pName[i], "CA" ) ) pertfact[i] = 0.292; else pertfact[i] = 1.000; } else { noletters = true; for( k = 0; k <= strlen( SplitLine[3] ) - 1; ++k ) { if( isalpha( SplitLine[3][k] ) ) { noletters = false; break; } } if( noletters ) pertfact[i] = atof( SplitLine[3] ); else pertfact[i] = 1.000; } } else pertfact[i] = 1.000; } fgets( line, sizeof( line ), SdlFile ); } /* Now see if we have any nOe data.... */ else if( strstr( line, "@NOE" ) ) { nOeAtLig = ( double** ) malloc( nnOe * sizeof( double* ) ); nOeAtProt = ( double** ) malloc( nnOe * sizeof( double* ) ); nOeAtLigID = ( int* ) malloc( nnOe * sizeof( int ) ); nOeAtProtID = ( int* ) malloc( nnOe * sizeof( int ) ); nOeRange = ( double** ) malloc( nnOe * sizeof( double* ) ); nOePenalty = ( double* ) malloc( nnOe * sizeof( double ) ); for( i = 0; i <= nnOe - 1; ++i ) { nOeRange[i] = ( double* ) malloc( 2 * sizeof( double ) ); nOeAtLig[i] = ( double* ) malloc( 3 * sizeof( double ) ); nOeAtProt[i] = ( double* ) malloc( 3 * sizeof( double ) ); parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ); nOeAtLigID[i] = atoi( SplitLine[0] ); nOeAtProtID[i] = atoi( SplitLine[1] ); nOeRange[i][0] = atof( SplitLine[2] ); nOeRange[i][1] = atof( SplitLine[3] ); nOePenalty[i] = atof( SplitLine[4] ); /* Look for the protons from the ligands with nOe on the arrays for the ligands. Done this way they will 'move' when the ligands move... */ for( j = 0; j <= nLigands - 1; ++j ) { for( k = 0; k <= nLigAt[j] - 1; ++k ) { if( ( LigAtID[j][k] ) == nOeAtLigID[i] ) nOeAtLig[i] = LigArray[j][k]; } } /* Now look for the protons with nOe data in the protein... */ for( j = 0; j <= nProtAt - 1; ++j ) { if( ( ProtAtID[j] ) == nOeAtProtID[i] ) nOeAtProt[i] = ProtArray[j]; } } fgets( line, sizeof( line ), SdlFile ); } /* Aromatic ring and anisotropic group arrays are allocated/filled in a similar fashion, but using the internal atom arrays - The same for mol2, pdb, or whatever... */ else if( strstr( line, "@AROMRINGS" ) ) { ifactor = ( double* ) malloc( nRings * sizeof( double ) ); nRingAt = ( int* ) malloc( nRings * sizeof( int ) ); RingArray = ( double*** ) malloc( nRings * sizeof( double** ) ); for( i = 0; i <= nRings - 1; ++i ) { nRingAt[i] = parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ) - 1; ifactor[i] = atof( SplitLine[0] ); RingArray[i] = ( double** ) malloc( nRingAt[i] * sizeof( double* ) ); for( j = 0; j <= nRingAt[i] - 1; ++j ) { RingArray[i][j] = ( double* ) malloc( 3 * sizeof( double ) ); for( k = 0; k <= nLigands - 1; ++k ) { for( l = 0; l <= nLigAt[k] - 1; ++l ) { if( ( LigAtID[k][l] ) == atoi( SplitLine[j + 1] ) ) RingArray[i][j] = LigArray[k][l]; } } } } fgets( line, sizeof( line ), SdlFile ); } else if( strstr( line, "@ANIGROUPS" ) ) { chiiso = ( double* ) malloc( nAni * sizeof( double ) ); chiani = ( double* ) malloc( nAni * sizeof( double ) ); nAniAt = ( int* ) malloc( nAni * sizeof( int ) ); AniArray = ( double*** ) malloc( nAni * sizeof( double** ) ); bndtyp = ( char** ) malloc( nAni * sizeof( char* ) ); for( i = 0; i <= nAni - 1; ++i ) { bndtyp[i] = ( char* ) malloc( 4 * sizeof( char ) ); nAniAt[i] = parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ) - 3; strcpy( bndtyp[i], SplitLine[0] ); chiiso[i] = atof( SplitLine[1] ); chiani[i] = atof( SplitLine[2] ); AniArray[i] = ( double** ) malloc( nAniAt[i] * sizeof( double* ) ); for( j = 0; j <= nAniAt[i] - 1; ++j ) { AniArray[i][j] = ( double* ) malloc( 3 * sizeof( double ) ); for( k = 0; k <= nLigands - 1; ++k ) { for( l = 0; l <= nLigAt[k] - 1; ++l ) { if( ( LigAtID[k][l] ) == atoi( SplitLine[j + 3] ) ) AniArray[i][j] = LigArray[k][l]; } } } } fgets( line, sizeof( line ), SdlFile ); } else if( strstr( line, "@INICOORDS" ) ) { if( nInicoord <= nLigands ) x = nInicoord; else x = nLigands; LigInicoord = ( double** ) malloc( x * sizeof( double* ) ); ligICidx = ( int* ) malloc( x * sizeof( int ) ); for( i = 0; i <= nInicoord - 1; ++i ) { if( i <= nLigands ) { LigInicoord[i] = ( double* ) malloc( 3 * sizeof( double ) ); if( parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ) > 3 ) { for( j = 0; j <= nLigands - 1; ++j ) { if( ( ! strcmp( ligName[j], SplitLine[0] ) ) || ( ( atoi( SplitLine[0] ) - 1 ) == j ) ) { ligICidx[i] = j; LigInicoord[i][0] = atof( SplitLine[1] ); LigInicoord[i][1] = atof( SplitLine[2] ); LigInicoord[i][2] = atof( SplitLine[3] ); break; } } } else { ligICidx[i] = i; LigInicoord[i][0] = atof( SplitLine[0] ); LigInicoord[i][1] = atof( SplitLine[1] ); LigInicoord[i][2] = atof( SplitLine[2] ); } } else fgets( line, sizeof( line ), SdlFile ); } fgets( line, sizeof( line ), SdlFile ); } /* Charged groups are easy. The only difference is that we have to split across the ":" to get the charge... */ else if( strstr( line, "@ELECTRO" ) ) { nEleAt = parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ) - 2; eleA = atof( SplitLine[0] ); eleB = atof( SplitLine[1] ); ElecAtoms = ( double** ) malloc( nEleAt * sizeof( double* ) ); charges = ( double* ) malloc( nEleAt * sizeof( double ) ); for( i = 0; i <= nEleAt - 1; ++i ) { ElecAtoms[i] = ( double* ) malloc( 3 * sizeof( double ) ); parseString( SplitLine[i + 2], &tempEleString, colon ); for( k = 0; k <= nLigands - 1; ++k ) { for( l = 0; l <= nLigAt[k] - 1; ++l ) { if( ( LigAtID[k][l] ) == atoi( tempEleString[0] ) ) { ElecAtoms[i] = LigArray[k][l]; if( ! strcmp( tempEleString[1], "sb" ) ) charges[i] = LigArray[k][l][3]; else charges[i] = atof( tempEleString[1] ); } } } } fgets( line, sizeof( line ), SdlFile ); } /* If there are H-bond acceptors defined, read them in... */ else if( strstr( line, "@HBACCEPT" ) ) { nHBAcc = parseString( fgets( line, sizeof( line ), SdlFile ), &SplitLine, whitespace ) - 2; HBSlope = atof( SplitLine[0] ); HBConst = atof( SplitLine[1] ); HBAccAtoms = ( double** ) malloc( nHBAcc * sizeof( double* ) ); for( i = 0; i <= nHBAcc - 1; ++i ) { HBAccAtoms[i] = ( double* ) malloc( 3 * sizeof( double ) ); for( k = 0; k <= nLigands - 1; ++k ) { for( l = 0; l <= nLigAt[k] - 1; ++l ) { if( ( LigAtID[k][l] ) == atoi( SplitLine[i+2] ) ) HBAccAtoms[i] = LigArray[k][l]; } } } fgets( line, sizeof( line ), SdlFile ); } else fgets( line, sizeof( line ), SdlFile ); } fclose( SdlFile ); /* Once we have the proton find their neighbors. Do it only if we are doing electrostatics or we have H-bonding acceptors defined in the ligands... */ if( ( nEleAt > 0 ) || ( nHBAcc > 0 ) ) { protoneigh = ( double** ) malloc( nProtons * sizeof( double* ) ); for( i = 0; i <= nProtons - 1; ++i ) { protoneigh[i] = ( double* ) malloc( 3 * sizeof( double ) ); mindistance = INFINITY; for( j = 0; j <= nProtAt - 1; ++j ) { if( protonid[i] != ProtAtID[j] ) { rPRPA = sqrt( ( protons[i][0] - ProtArray[j][0] ) * ( protons[i][0] - ProtArray[j][0] ) + ( protons[i][1] - ProtArray[j][1] ) * ( protons[i][1] - ProtArray[j][1] ) + ( protons[i][2] - ProtArray[j][2] ) * ( protons[i][2] - ProtArray[j][2] ) ); if( rPRPA < mindistance ) { neighindex = j; mindistance = rPRPA; } } } protoneigh[i] = ProtArray[neighindex]; } } /* Finally, move ligands to the positions indicated in the sdl file if requested. Also, see if we want to constraint the center of the ligands. Start by calculating the current centers of the ligands... */ if( ( ! strcmp( constcent, "YES" ) ) || ( ( nInicoord > 0 ) && ( ! strcmp( useinicoord, "YES" ) ) ) ) { initcent = ( double** ) malloc( nLigands * sizeof( double* ) ); for( i = 0; i <= nLigands - 1; ++i ) { initcent[i] = ( double* ) calloc( 3, sizeof( double ) ); for( j = 0; j <= nLigAt[i] - 1; ++j ) { for( k = 0; k <= 2; ++k ) initcent[i][k] += LigArray[i][j][k]; } for( k = 0; k <= 2; ++k ) initcent[i][k] /= nLigAt[i]; } if( ( nInicoord > 0 ) && ( ! strcmp( useinicoord, "YES" ) ) ) { for( i = 0; i <= nInicoord - 1; ++i ) { for( k = 0; k <= 2; ++k ) { translate( LigArray[ligICidx[i]], nLigAt[ligICidx[i]], ( LigInicoord[i][k] - initcent[ligICidx[i]][k] ), k ); initcent[ligICidx[i]][k] = LigInicoord[i][k]; } } for( i = 0; i <= nInicoord - 1; ++i ) free( LigInicoord[i] ); free( LigInicoord ); free( ligICidx ); } if( ! strcmp( constcent, "NO" ) ) { for( i = 0; i <= nLigands - 1; ++i ) free( initcent[i] ); free( initcent ); } } /* OK, so here we go. All the data is in place, and we have to start with the optimization. Set the initial iteration and zero enegy and statistics variables... */ up = true; iter = 1; AvgDeltaDiff = 0; AvgDelta2Diff = 0; AvgDeltaDiff2 = 0; AvgDeltaSD = 0; ffE = 0; shiftE = 0; nOeE = 0; selfE = 0; torE = 0; anchorE = 0; eterms = 0; /* First calculate the initial energy contributions from each term and let the user know how far the calculated shifts are from the experimental ones. Calculate the average difference between observed and calculated shifts and its standard deviation... */ if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { ffE = nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); ++eterms; } if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) ++eterms; for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) shiftE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); printf( "Proton %d: %s - Experimental shift = %f - Calculated shift = %f\n", protonid[j], protonFullName[j], dexp[j], shift[j] ); /* Compute statistics... */ AvgDeltaDiff += fabs( shift[j] - dexp[j] ); AvgDelta2Diff += ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } AvgDeltaDiff /= nProtons; AvgDelta2Diff /= nProtons; AvgDeltaDiff2 = AvgDeltaDiff * AvgDeltaDiff; AvgDeltaSD = sqrt( fabs( AvgDelta2Diff - AvgDeltaDiff2 ) ); printf( "\n" ); printf( "Initial shift deviation = %f +/- %f\n", AvgDeltaDiff, AvgDeltaSD ); /* If we have nOe data in the .sdl file, call the function here... */ if( nnOe > 0 ) { nOeE = noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "INI" ); ++eterms; } /* If there are rotatable bonds, get the appropriate energy terms... */ if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { selfE = slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); torE = gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); eterms += 2; } /* If we defined an anchor for the ligand centroid, increase the eterm counter, but don't print the violation, as initially there should be none... */ if( ! strcmp( constcent, "YES" ) ) ++eterms; /* Get the total energy... */ iniE = ffE + shiftE + nOeE + selfE + torE; /* If we have more than one energy term contributing to the total energy, print our a breakdown of the energy contributions... */ if( eterms > 1 ) { printf( "\n" ); if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) printf( "Initial force-field energy contribution = %f Kcal/mol\n", ffE ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) printf( "Initial shift penalty energy contribution = %f Kcal/mol\n", shiftE ); if( nnOe > 0 ) printf( "Initial nOe penalty energy contribution = %f Kcal/mol\n", nOeE ); if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { printf( "Initial self-energy contribution = %f Kcal/mol\n", selfE ); printf( "Initial torsional energy contribution = %f Kcal/mol\n", torE ); } if( ! strcmp( constcent, "YES" ) ) printf( "Anchor potential activated\n" ); } /* If we want to do the simple 'ragdoll' line-search, go for it... */ if( ! strcmp( optmethod, "RGDL" ) ) { printf( "\n" ); printf( "Enetering SDiLiCon ragdoll geometry optimization routine\n" ); printf( "\n" ); printf( "Starting optimization cylce - Iteration %d - Current energy = %f\n", iter, iniE ); rgdlfctr = 0; while( up ) { /* Define the 'oozing' variable, which lets the program decide if we are just oozing around a certain point with no improvements. Every time we improve in one cycle we may require to go over all the variables using that step again, so if there is an improvement, the variable is reset to FALSE... */ oozing = true; /* Claculate the self and torsional energy here, as they may have changed after bond rotations... */ if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { selfE = slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); torE = gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); } /* We have to loop over all the defined ligands... */ for( l = 0; l <= nLigands - 1; ++ l ) { if( nLigands > 1 ) { printf( "\n" ); printf( "Working on ligand %s...\n", ligName[l] ); printf( "\n" ); } /* Treat translations and rotations separately. First, translations, but check that we really want to do translations here... */ if( ! strcmp( dotrans, "YES" ) ) { for( axes = 0; axes <= 2; ++axes ) { /* Move the ligand one way and compute the energy... */ printf( "Trying %f Angstrom translations along %s\n", delta, tvlabel[axes] ); upE = selfE + torE; translate( LigArray[l], nLigAt[l], delta, axes ); if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) upE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; upE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } /* If we have nOe data, get the energy... */ if( nnOe > 0 ) upE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } /* If we want to keep the molecule within a certain distance from its original position, calculate the constraint here... */ if( ! strcmp( constcent, "YES" ) ) upE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); /* If this was productive, keep going this way... */ if( upE < iniE ) { printf( "Initial energy = %f - Energy after +%s translation = %f\n", iniE, tvlabel[axes], upE ); while( upE < iniE ) { oozing = false; iniE = upE; upE = selfE + torE; translate( LigArray[l], nLigAt[l], delta, axes ); if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) upE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; upE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) upE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ! strcmp( constcent, "YES" ) ) upE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); printf( "Initial energy = %f - Energy after +%s translation = %f\n", iniE, tvlabel[axes], upE ); } /* Now, if this direction stops decreasing, reset the position to the previous one... */ translate( LigArray[l], nLigAt[l], -1.0 * delta, axes ); printf( "No more progress with +%s translation...\n", tvlabel[axes] ); } /* If our initial move did not improve the energy of the system, move it back to the original position and go the other way... */ else { downE = selfE + torE; translate( LigArray[l], nLigAt[l], -2.0 * delta, axes ); if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) downE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni -1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; downE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) downE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ! strcmp( constcent, "YES" ) ) downE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); /* If this way turns out to be good, keep at it... */ if( downE < iniE ) { printf( "Initial energy = %f - Energy after -%s translation = %f\n", iniE, tvlabel[axes], downE ); while( downE < iniE ) { oozing = false; iniE = downE; downE = selfE + torE; translate( LigArray[l], nLigAt[l], -1.0 * delta, axes ); if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) downE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; downE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) downE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ! strcmp( constcent, "YES" ) ) downE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); printf( "Initial energy = %f - Energy after -%s translation = %f\n", iniE, tvlabel[axes], downE ); } /* As before, when we go past the minimum, bail out and go back to the previous conformation... */ translate( LigArray[l], nLigAt[l], delta, axes ); printf( "No more progress with -%s translation...\n", tvlabel[axes] ); } /* If nothing happens this way either, move it back to the starting point and repeat the cycle... */ else translate( LigArray[l], nLigAt[l], delta, axes ); } } } if( ! strcmp( dorot, "YES" ) ) { /* Calculate the anchor energy once, as it won't change during rigid-body rotations... */ if( ! strcmp( constcent, "YES" ) ) anchorE = centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); for( axes = 0; axes <= 2; ++axes ) { printf( "Trying %f degree rotations about %s\n", theta, rvlabel[axes] ); rotate( LigArray[l], nLigAt[l], theta, axes ); upE = selfE + torE + anchorE; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) upE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; upE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) upE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( upE < iniE ) { printf( "Initial energy = %f - Energy after +%s rotation = %f\n", iniE, rvlabel[axes], upE ); while( upE < iniE ) { oozing = false; iniE = upE; rotate( LigArray[l], nLigAt[l], theta, axes ); upE = selfE + torE + anchorE; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) upE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; upE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) upE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } printf( "Initial energy = %f - Energy after +%s rotation = %f\n", iniE, rvlabel[axes], upE ); } rotate( LigArray[l], nLigAt[l], -1.0 * theta, axes ); printf( "No more progress with +%s rotation...\n", rvlabel[axes] ); } else { rotate( LigArray[l], nLigAt[l], -2.0 * theta, axes ); downE = selfE + torE + anchorE; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) downE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; downE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) downE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( downE < iniE ) { printf( "Initial energy = %f - Energy after -%s rotation = %f\n", iniE, rvlabel[axes], downE ); while( downE < iniE ) { oozing = false; iniE = downE; rotate( LigArray[l], nLigAt[l], -1.0 * theta, axes ); downE = selfE + torE + anchorE; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) downE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; downE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) downE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } printf( "Initial energy = %f - Energy after -%s rotation = %f\n", iniE, rvlabel[axes], downE ); } rotate( LigArray[l], nLigAt[l], theta, axes ); printf( "No more progress with -%s rotation...\n", rvlabel[axes] ); } else rotate( LigArray[l], nLigAt[l], theta, axes ); } } } } /* Bond rotations would have to be done here, irrespective of what ligand they belong to... */ if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { printf( "\n" ); printf( "Working on %d rotatable bonds...\n", nRotBonds ); printf( "\n" ); for( l = 0; l <= nRotBonds - 1; ++l ) { printf( "Trying %f degree rotations about %d-%d bond\n", maxbrot[l], TorAtA[l], TorAtB[l] ); rotbond( nFragAt[l], FragArray[l], maxbrot[l] ); selfE = slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); torE = gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); upE = selfE + torE; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) upE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; upE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) upE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ! strcmp( constcent, "YES" ) ) upE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); if( upE < iniE ) { printf( "Initial energy = %f - Energy after +%d-%d bond rotation = %f\n", iniE, TorAtA[l], TorAtB[l], upE ); while( upE < iniE ) { oozing = false; iniE = upE; rotbond( nFragAt[l], FragArray[l], maxbrot[l] ); selfE = slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); torE = gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); upE = selfE + torE; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) upE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; upE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) upE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ! strcmp( constcent, "YES" ) ) upE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); printf( "Initial energy = %f - Energy after +%d-%d bond rotation = %f\n", iniE, TorAtA[l], TorAtB[l], upE ); } rotbond( nFragAt[l], FragArray[l], -1.0 * maxbrot[l] ); printf( "No more progress with +%d-%d bond rotation...\n", TorAtA[l], TorAtB[l] ); } else { rotbond( nFragAt[l], FragArray[l], -2.0 * maxbrot[l] ); selfE = slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); torE = gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); downE = selfE + torE; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) downE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; downE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) downE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ! strcmp( constcent, "YES" ) ) downE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); if( downE < iniE ) { printf( "Initial energy = %f - Energy after -%d-%d bond rotation = %f\n", iniE, TorAtA[l], TorAtB[l], downE ); while( downE < iniE ) { oozing = false; iniE = downE; rotbond( nFragAt[l], FragArray[l], -1.0 * maxbrot[l] ); selfE = slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); torE = gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); downE = selfE + torE; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) downE += nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; downE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) downE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ! strcmp( constcent, "YES" ) ) downE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); printf( "Initial energy = %f - Energy after -%d-%d bond rotation = %f\n", iniE, TorAtA[l], TorAtB[l], downE ); } rotbond( nFragAt[l], FragArray[l], maxbrot[l] ); printf( "No more progress with -%d-%d bond rotation...\n", TorAtA[l], TorAtB[l] ); } else rotbond( nFragAt[l], FragArray[l], maxbrot[l] ); } } } /* End of grand loop over all ligands and bonds. If none of the translations/rotations improved things, 'oozing' will be TRUE. We can therefore half the step to improve the resolution of the search... */ if( oozing ) { printf( "\n" ); printf( "Reducing Step...\n" ); if( ! strcmp( dotrans, "YES" ) ) delta /= 2; if( ! strcmp( dorot, "YES" ) ) theta /= 2; if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { brconvavg = 0; for( l = 0; l <= nRotBonds - 1; ++l ) { maxbrot[l] /= 2; brconvavg += maxbrot[l] / nRotBonds; } } /* If the resolution is smaller than the desired one, bail out of the search... */ if( ( ( delta < convtrans ) && ( ! strcmp( dotrans, "YES" ) ) ) || ( ( theta < convrot ) && ( ! strcmp( dorot, "YES" ) ) ) || ( ( brconvavg < convbrot ) && ( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) ) ) { printf( "\n" ); printf( "Convergnce reached!\n" ); printf( "\n" ); up = false; } /* if the user selected to save multiple conformers, save those from successive steps... */ if( ! strcmp( multconf, "YES" ) ) { ++rgdlfctr; sprintf( conffile, "%s%s%03d", outfile, ".", rgdlfctr ); if( writemol( ftype, conffile, HeaderCount, HeaderJunk, AtomCount, RawAtomData, nLigands, nLigAt, LigAtID, LigArray, FooterCount, FooterJunk, iniE ) ) { printf( "Error opening output file %s. Check file permissions and try again...\n", conffile ); return 0; } } } if( ( iter >= maxiter ) && ( up ) ) { printf( "\n" ); printf( "Maximum number of iterations reached!\n" ); printf( "\n" ); up = false; } if( up ) { ++iter; if( ! strcmp( verboutput, "YES" ) ) { printf( "\n" ); for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; printf( "Proton %d: %s - Experimental shift = %f - Calculated shift = %f\n", protonid[j], protonFullName[j], dexp[j], shift[j] ); } if( nnOe > 0 ) noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "YES" ); } printf( "\n" ); printf( "Re-starting optimization cycle - Iteration %d - Current energy = %f\n", iter, iniE ); } } } /* If, on the other hand, we want to do a stochastic RIPS search, do it here... */ else if( ( ! strcmp( optmethod, "RIPS" ) ) || ( ! strcmp( optmethod, "MCMD" ) ) || ( ! strcmp( optmethod, "MCSA" ) ) ) { printf( "\n" ); printf( "Enetering SDiLiCon %s geometry optimization routine\n", optmethod ); printf( "\n" ); /* Initialize the randtrans, randrot, and randbrot arrays and Seed the random number generator... */ if( ! strcmp( dotrans, "YES" ) ) { randtrans = ( double** ) malloc( nLigands * sizeof( double* ) ); for( i = 0; i <= nLigands - 1; ++i ) randtrans[i] = ( double* ) malloc( 3 * sizeof( double ) ); } if( ! strcmp( dorot, "YES" ) ) { randrot = ( double** ) malloc( nLigands * sizeof( double* ) ); for( i = 0; i <= nLigands - 1; ++i ) randrot[i] = ( double* ) malloc( 3 * sizeof( double ) ); } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) randbrot = ( double* ) malloc( nRotBonds * sizeof( double ) ); srandom( rseed ); /* Initialize the productive and non-productive step counters... */ pstep = 0; npstep = 0; savectr = 1; dosave = false; /* See if we want to do Metropolis MD or SA. If so, initialize the temperature array... */ if( strcmp( optmethod, "RIPS" ) ) { if( ( ! strcmp( runningavg, "YES" ) ) || ( ! strcmp( runningavg, "OPT" ) ) ) { currshft = ( double* ) malloc( nProtons * sizeof( double ) ); shftavg = ( double* ) calloc( nProtons, sizeof( double ) ); if( ! strcmp( runningavg, "OPT" ) ) shftwavg = ( double* ) calloc( nProtons, sizeof( double ) ); for( j = 0; j <= nProtons - 1; ++j ) currshft[j] = shift[j]; } temperature = ( double* ) malloc( maxphits * sizeof( double ) ); for( i = 0; i <= maxphits - 1; ++i ) temperature[i] = tlkelvin; if( ! strcmp( optmethod, "MCSA" ) ) { sastep = maxphits / mdcycles; sastpvt = int( tpivot * sastep ); if( ! strcmp( multconf, "YES" ) ) { dosave = true; savectr = sastep; } /* This exponential factor will get within 1% of the low temperature... */ expfact = ( 1 + sastpvt - sastep ) / ( log( 0.01 * tlkelvin / ( thkelvin - tlkelvin ) ) ); for( i = 0; i <= mdcycles - 1; ++i ) { for( j = 0; j <= sastpvt - 1; ++j ) temperature[i*sastep+j] = thkelvin; for( j = 0; j <= ( sastep - sastpvt ) - 1; ++j ) { if( ! strcmp( annealfunc, "LIN" ) ) temperature[i*sastep+sastpvt+j] = thkelvin - ( thkelvin - tlkelvin ) * j / ( sastep - sastpvt - 1 ); else if( ! strcmp( annealfunc, "STEP" ) ) temperature[i*sastep+sastpvt+j] = tlkelvin; else temperature[i*sastep+sastpvt+j] = tlkelvin + ( thkelvin - tlkelvin ) * exp( -1 * j / expfact ); } } } else { if( ( ! strcmp( multconf, "YES" ) ) && ( mdcycles == 0 ) ) mdcycles = 10; if( mdcycles > 0 ) { dosave = true; savectr = mdcycles; } } } else { if( ! strcmp( multconf, "YES" ) ) dosave = true; } while( up ) { /* Initialize the energy... */ ripsE = 0; shiftE = 0; /* Loop through the ligands, translating/rotating them randomly... */ if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( axes = 0; axes <= 2; ++axes ) { randtrans[l][axes] = delta * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; translate( LigArray[l], nLigAt[l], randtrans[l][axes], axes ); } } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( axes = 0; axes <= 2; ++axes ) { randrot[l][axes] = theta * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; rotate( LigArray[l], nLigAt[l], randrot[l][axes], axes ); } } } /* If there are rotatable bonds defined, rotate them here... */ if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( l = 0; l <= nRotBonds - 1; ++l ) { randbrot[l] = maxbrot[l] * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; rotbond( nFragAt[l], FragArray[l], randbrot[l] ); } } /* Compute/re-compute energy... */ if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) ripsE = nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; if( ! strcmp( runningavg, "OPT" ) ) shiftE += force[j] * ( ( shift[j] / ( tau * ( 1.0 - exp( -1.0 * iter / tau ) ) ) + exp( -1.0 / tau ) * shftwavg[j] ) - dexp[j] ) * ( ( shift[j] / ( tau * ( 1.0 - exp( -1.0 * iter / tau ) ) ) + exp( -1.0 / tau ) * shftwavg[j] ) - dexp[j] ); else shiftE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { ripsE += shiftE; if( nnOe > 0 ) ripsE += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { ripsE += slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); ripsE += gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); } if( ! strcmp( constcent, "YES" ) ) ripsE += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); /* If this random step was non-productive, reset to the original position. When we do the rotations backwards, make them in the reverse order. The rotation operator does not commute when applied to more than one axis. Also be careful, as the backward bond rotations also have to be done first... */ nonproductive = true; if( ! strcmp( optmethod, "RIPS" ) ) { if( ripsE < iniE ) nonproductive = false; } else { if( ripsE < iniE ) nonproductive = false; else if( ( ( double ) random( ) ) / RAND_MAX < exp( ( iniE - ripsE ) / ( gasconst * temperature[pstep] ) ) ) nonproductive = false; } if( nonproductive ) { if( ( ! strcmp( runningavg, "YES" ) ) || ( ! strcmp( runningavg, "OPT" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shftavg[j] = ( ( iter - 1 ) * shftavg[j] + currshft[j] ) / iter; if( ! strcmp( runningavg, "OPT" ) ) shftwavg[j] = currshft[j] / ( tau * ( 1 - exp( -1.0 * iter / tau ) ) ) + exp( -1.0 / tau ) * shftwavg[j]; } } ++npstep; if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( l = 0; l <= nRotBonds - 1; ++l ) rotbond( nFragAt[l], FragArray[l], -1.0 * randbrot[l] ); } if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( axes = 0; axes <= 2; ++axes ) translate( LigArray[l], nLigAt[l], -1.0 * randtrans[l][axes], axes ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( axes = 2; axes >= 0; --axes ) rotate( LigArray[l], nLigAt[l], -1.0 * randrot[l][axes], axes ); } } if( ! strcmp( verboutput, "YES" ) ) { if( ! strcmp( optmethod, "RIPS" ) ) printf( "Non-productive step %d - Iteration %d - Initial energy = %g - Current energy = %g\n", npstep, iter, iniE, ripsE ); else printf( "Non-productive step %d - Iteration %d - Temperature = %.2f - Energy = %g\n", npstep, iter, temperature[pstep], ripsE ); } } /* If this step was productive, use these as the starting conformations of the ligands for the next pulse... */ else { if( ( ! strcmp( runningavg, "YES" ) ) || ( ! strcmp( runningavg, "OPT" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { currshft[j] = shift[j]; shftavg[j] = ( ( iter - 1 ) * shftavg[j] + currshft[j] ) / iter; if( ! strcmp( runningavg, "OPT" ) ) shftwavg[j] = currshft[j] / ( tau * ( 1 - exp( -1.0 * iter / tau ) ) ) + exp( -1.0 / tau ) * shftwavg[j]; } } ++pstep; if( ! strcmp( optmethod, "RIPS" ) ) printf( "Productive step %d - Iteration %d - Initial Energy = %g - Current energy = %g\n", pstep, iter, iniE, ripsE ); else printf( "Productive step %d - Iteration %d - Temperature = %.2f - Energy = %g\n", pstep, iter, temperature[pstep-1], ripsE ); if( ! strcmp( verboutput, "YES" ) ) { printf( "\n" ); if( ( ! strcmp( runningavg, "YES" ) ) || ( ! strcmp( runningavg, "OPT" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) printf( "Proton %d: %s - Exp. shift = %f - Calc. shift = %f - Avg. shift = %f\n", protonid[j], protonFullName[j], dexp[j], shift[j], shftavg[j] ); } else { for( j = 0; j <= nProtons - 1; ++j ) printf( "Proton %d: %s - Experimental shift = %f - Calculated shift = %f\n", protonid[j], protonFullName[j], dexp[j], shift[j] ); } if( nnOe > 0 ) noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "YES" ); printf( "\n" ); } /* If we have the multiple confomrer flag 'ON', save the conformer here... */ if( ( dosave ) && ( pstep % savectr == 0 ) ) { sprintf( conffile, "%s%s%03d", outfile, ".", pstep ); if( writemol( ftype, conffile, HeaderCount, HeaderJunk, AtomCount, RawAtomData, nLigands, nLigAt, LigAtID, LigArray, FooterCount, FooterJunk, ripsE ) ) { printf( "Error opening output file %s. Check file permissions and try again...\n", conffile ); return 0; } } /* Compute the energy gradient between productive runs... */ econv = fabs( iniE - ripsE ); iniE = ripsE; } ++iter; /* If the user decided to use the enrgy gradient between runs to stop the search, check at this point... */ if( ( ! strcmp( optmethod, "RIPS" ) ) && ( egrad > 0 ) && ( econv <= egrad ) ) { printf( "\n" ); printf( "Energy convergence criterion reached!\n" ); printf( "\n" ); up = false; } /* Increment or decrement the maximum step sizes to get 50/50 productive versus non-productive hits. We do this multiplying the steps by 5%... */ if( ! strcmp( adjuststep, "YES" ) ) { if( ( double ) pstep / ( npstep + pstep ) < acceptratio ) { delta *= 0.95; theta *= 0.95; for( i = 0; i <= nRotBonds - 1; ++i ) maxbrot[i] *= 0.95; } else { delta *= 1.05; theta *= 1.05; for( i = 0; i <= nRotBonds - 1; ++i ) maxbrot[i] *= 1.05; } } /* If a maximum number of productive steps was entered and we reach it, bail out... */ if( ( maxphits > 0 ) && ( pstep >= maxphits ) ) { if( strcmp( verboutput, "YES" ) ) printf( "\n" ); printf( "Maximum number of productive steps reached!\n" ); printf( "\n" ); up = false; } /* Finally, if we go over the maximum number of steps, halt the search... */ if( ( ! strcmp( optmethod, "RIPS" ) ) && ( iter >= maxiter ) ) { if( strcmp( verboutput, "YES" ) ) printf( "\n" ); printf( "Maximum number of iterations reached!\n" ); printf( "\n" ); up = false; } } printf( "Total iterations completed: %d\n", iter ); printf( "Productive steps completed: %d\n", pstep ); printf( "Non-productive steps completed: %d\n", npstep ); printf( "\n" ); if( ( ! strcmp( runningavg, "YES" ) ) || ( ! strcmp( runningavg, "OPT" ) ) ) { AvgDeltaDiff = 0; AvgDelta2Diff = 0; for( j = 0; j <= nProtons - 1; ++j ) { printf( "Proton %d: %s - Experimental shift = %f - Avgerage shift = %f\n", protonid[j], protonFullName[j], dexp[j], shftavg[j] ); AvgDeltaDiff += fabs( shftavg[j] - dexp[j] ); AvgDelta2Diff += ( shftavg[j] - dexp[j] ) * ( shftavg[j] - dexp[j] ); } AvgDeltaDiff /= nProtons; AvgDelta2Diff /= nProtons; AvgDeltaDiff2 = AvgDeltaDiff * AvgDeltaDiff; AvgDeltaSD = sqrt( fabs( AvgDelta2Diff - AvgDeltaDiff2 ) ); printf( "\n" ); printf( "Final running average shift deviation = %f +/- %f\n", AvgDeltaDiff, AvgDeltaSD ); printf( "\n" ); } } /* Finally, we may have chosen a genetic algorithm... */ else if( ! strcmp( optmethod, "GENA" ) ) { printf( "\n" ); printf( "Enetering SDiLiCon GA geometry optimization routine\n" ); printf( "\n" ); /* Compute the number of variables in each gene as well as the 'average' threshold for checking gene convergence (i.e., duplicate checking)... */ activegavar = 0; geneconv = 0; if( ! strcmp( dotrans, "YES" ) ) { activegavar += 3; geneconv += 3 * convtrans; } if( ! strcmp( dorot, "YES" ) ) { activegavar += 3; geneconv += 3 * convrot; } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { activegavar += nRotBonds; geneconv += nRotBonds * convbrot; } /* Figure out the smallest population... */ smallestpop = INT_MAX; for( k = 0; k <= nIslands - 1; ++k ) { if( popsize[k] < smallestpop ) smallestpop = popsize[k]; } if( keepbest > smallestpop ) keepbest = smallestpop; if( travellers > smallestpop ) travellers = smallestpop; /* Seed the random number generator, allocate memory for the genes of parents and offspring, and fill with the initial population(s)... */ migcounter = 1; srandom( rseed ); gaE = ( double** ) malloc( nIslands * sizeof( double* ) ); gascore = ( double** ) malloc( nIslands * sizeof( double* ) ); pstart = ( double** ) malloc( nIslands * sizeof( double* ) ); pend = ( double** ) malloc( nIslands * sizeof( double* ) ); rankvec = ( int** ) malloc( nIslands * sizeof( int* ) ); diversity = ( int* ) malloc( nIslands * sizeof( int ) ); scoresum = ( double* ) malloc( nIslands * sizeof( double ) ); transgene = ( double**** ) malloc( nIslands * sizeof( double*** ) ); transchld = ( double**** ) malloc( nIslands * sizeof( double*** ) ); rotgene = ( double**** ) malloc( nIslands * sizeof( double*** ) ); rotchld = ( double**** ) malloc( nIslands * sizeof( double*** ) ); brotgene = ( double*** ) malloc( nIslands * sizeof( double** ) ); brotchld = ( double*** ) malloc( nIslands * sizeof( double** ) ); for( k = 0; k <= nIslands - 1; ++ k ) { gaE[k] = ( double* ) malloc( popsize[k] * sizeof( double ) ); gascore[k] = ( double* ) malloc( popsize[k] * sizeof( double ) ); pstart[k] = ( double* ) malloc( popsize[k] * sizeof( double ) ); pend[k] = ( double* ) malloc( popsize[k] * sizeof( double ) ); rankvec[k] = ( int* ) malloc( popsize[k] * sizeof( int ) ); if( ! strcmp( dotrans, "YES" ) ) { transgene[k] = ( double*** ) malloc( nLigands * sizeof( double** ) ); transchld[k] = ( double*** ) malloc( nLigands * sizeof( double** ) ); for( l = 0; l <= nLigands - 1; ++l ) { transgene[k][l] = ( double** ) malloc( popsize[k] * sizeof( double* ) ); transchld[k][l] = ( double** ) malloc( popsize[k] * sizeof( double* ) ); for( i = 0; i <= popsize[k] - 1; ++i ) { transgene[k][l][i] = ( double* ) malloc( 3 * sizeof( double ) ); transchld[k][l][i] = ( double* ) malloc( 3 * sizeof( double ) ); for( j = 0; j <= 2; ++j ) transgene[k][l][i][j] = delta * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; } } } if( ! strcmp( dorot, "YES" ) ) { rotgene[k] = ( double*** ) malloc( nLigands * sizeof( double** ) ); rotchld[k] = ( double*** ) malloc( nLigands * sizeof( double** ) ); for( l = 0; l <= nLigands - 1; ++l ) { rotgene[k][l] = ( double** ) malloc( popsize[k] * sizeof( double* ) ); rotchld[k][l] = ( double** ) malloc( popsize[k] * sizeof( double* ) ); for( i = 0; i <= popsize[k] - 1; ++i ) { rotgene[k][l][i] = ( double* ) malloc( 3 * sizeof( double ) ); rotchld[k][l][i] = ( double* ) malloc( 3 * sizeof( double ) ); for( j = 0; j <= 2; ++j ) rotgene[k][l][i][j] = theta * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; } } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { brotgene[k] = ( double** ) malloc( popsize[k] * sizeof( double* ) ); brotchld[k] = ( double** ) malloc( popsize[k] * sizeof( double* ) ); for( i = 0; i <= popsize[k] - 1; ++i ) { brotgene[k][i] = ( double* ) malloc( nRotBonds * sizeof( double ) ); brotchld[k][i] = ( double* ) malloc( nRotBonds * sizeof( double ) ); for( j = 0; j <= nRotBonds - 1; ++j ) brotgene[k][i][j] = maxbrot[j] * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; } } } /* Now we have the initial population(s) set up. Here we should have a continuous 'while'... */ while( up ) { norescore = false; didmigration = false; /* This loop is needed in case migrations take place. If they do, the populations have to be re-ranked... */ while( up ) { for( k = 0; k <= nIslands - 1; ++k ) { for( i = 0; i <= popsize[k] - 1; ++i ) rankvec[k][i]= i; /* Now loop over all the members of the population, or in other words, get the 'gene products' of each gene by rotating/translating/etc... */ for( i = 0; i <= popsize[k] - 1; ++i ) { if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) translate( LigArray[l], nLigAt[l], transgene[k][l][i][j], j ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) rotate( LigArray[l], nLigAt[l], rotgene[k][l][i][j], j ); } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) rotbond( nFragAt[j], FragArray[j], brotgene[k][i][j] ); } /* Now compute the energy... */ gaE[k][i] = 0; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) gaE[k][i] = nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( m = 0; m <= nRings - 1; ++m ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[m], nRingAt[m], ifactor[m], method ); } if( nAni > 0 ) { for( m = 0; m <= nAni - 1; ++m ) shift[j] += anishift( protons[j], AniArray[m], nAniAt[m], chiiso[m], chiani[m], bndtyp[m] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; gaE[k][i] += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } if( nnOe > 0 ) gaE[k][i] += noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "NO" ); } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { gaE[k][i] += slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); gaE[k][i] += gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); } if( ! strcmp( constcent, "YES" ) ) gaE[k][i] += centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); /* Finally, move it back to the origin. A bit ineffective, but simple... */ if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) rotbond( nFragAt[j], FragArray[j], -1.0 * brotgene[k][i][j] ); } if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) translate( LigArray[l], nLigAt[l], -1.0 * transgene[k][l][i][j], j ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 2; j >= 0; --j ) rotate( LigArray[l], nLigAt[l], -1.0 * rotgene[k][l][i][j], j ); } } } /* Rank the energies using a bubble-sort... */ for( j = 0; j <= popsize[k] - 1; ++j ) { for( i = 0; i <= popsize[k] - 2 - j; ++i ) { if( gaE[k][rankvec[k][i+1]] < gaE[k][rankvec[k][i]] ) { temprank = rankvec[k][i]; rankvec[k][i] = rankvec[k][i+1]; rankvec[k][i+1] = temprank; } } } /* Score the members of the current population. ... */ scoresum[k] = 0; if( ! strcmp( scorefunc, "RANK" ) ) { for( i = 0; i <= popsize[k] - 1; ++i ) { gascore[k][rankvec[k][i]] = ( double ) ( popsize[k] - i ); scoresum[k] += gascore[k][rankvec[k][i]]; } } else if( ! strcmp( scorefunc, "LOGT" ) ) { for( i = 0; i <= popsize[k] - 1; ++i ) { gascore[k][i] = log10( 1 + fabs( gaE[k][i] - gaE[k][rankvec[k][popsize[k]-1]] ) ); scoresum[k] += gascore[k][i]; } } else if( ! strcmp( scorefunc, "LOGN" ) ) { for( i = 0; i <= popsize[k] - 1; ++i ) { gascore[k][i] = log( 1 + fabs( gaE[k][i] - gaE[k][rankvec[k][popsize[k]-1]] ) ); scoresum[k] += gascore[k][i]; } } else if( ! strcmp( scorefunc, "SQRT" ) ) { for( i = 0; i <= popsize[k] - 1; ++i ) { gascore[k][i] = sqrt( fabs( gaE[k][i] - gaE[k][rankvec[k][popsize[k]-1]] ) ); scoresum[k] += gascore[k][i]; } } else { for( i = 0; i <= popsize[k] - 1; ++i ) { gascore[k][i] = fabs( gaE[k][i] - gaE[k][rankvec[k][popsize[k]-1]] ); scoresum[k] += gascore[k][i]; } } } /* Here check for the migration frequency, and if we have to do a migration, pick the 'origin' and 'destination' islands at random. Once selected, take the number of individuals specified to migrate from one, and use them to replace the worst individuals from the destination island... */ if( ( nIslands > 1 ) && ( ! didmigration ) && ( iter > 1 ) ) { domigration = false; if( ! strcmp( freqoritvl, "FREQ" ) ) { if( ( ( ( double ) random( ) ) / RAND_MAX ) < migfreqf ) domigration = true; } else { if( migcounter == migfreqi ) { domigration = true; migcounter = 1; } else ++migcounter; } if( domigration ) { didmigration = true; mIsland = int( ( nIslands * ( double ) random( ) ) / RAND_MAX ); dIsland = int( ( nIslands * ( double ) random( ) ) / RAND_MAX ); while( mIsland == dIsland ) dIsland = int( ( nIslands * ( double ) random( ) ) / RAND_MAX ); if( ! strcmp( migmethod, "SWAP" ) ) { for( i = 0; i <= travellers - 1; ++i ) { if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) { tempgene = transgene[dIsland][l][rankvec[dIsland][i]][j]; transgene[dIsland][l][rankvec[dIsland][i]][j] = transgene[mIsland][l][rankvec[mIsland][i]][j]; transgene[mIsland][l][rankvec[mIsland][i]][j] = tempgene; } } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) { tempgene = rotgene[dIsland][l][rankvec[dIsland][i]][j]; rotgene[dIsland][l][rankvec[dIsland][i]][j] = rotgene[mIsland][l][rankvec[mIsland][i]][j]; rotgene[mIsland][l][rankvec[mIsland][i]][j] = tempgene; } } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) { tempgene = brotgene[dIsland][rankvec[dIsland][i]][j]; brotgene[dIsland][rankvec[dIsland][i]][j] = brotgene[mIsland][rankvec[mIsland][i]][j]; brotgene[mIsland][rankvec[mIsland][i]][j] = tempgene; } } } } else if( ! strcmp( migmethod, "COPY" ) ) { for( i = 0; i <= travellers - 1; ++i ) { if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) transgene[dIsland][l][rankvec[dIsland][popsize[dIsland]-i-1]][j] = transgene[mIsland][l][rankvec[mIsland][i]][j]; } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) rotgene[dIsland][l][rankvec[dIsland][popsize[dIsland]-i-1]][j] = rotgene[mIsland][l][rankvec[mIsland][i]][j]; } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) brotgene[dIsland][rankvec[dIsland][popsize[dIsland]-i-1]][j] = brotgene[mIsland][rankvec[mIsland][i]][j]; } } } } else norescore = true; } else norescore = true; if( norescore ) break; } /* Get the index of the best island... */ bestislandE = INFINITY; for( k = 0; k <= nIslands - 1; ++k ) { if( gaE[k][rankvec[k][0]] < bestislandE ) { bestislandE = gaE[k][rankvec[k][0]]; bestisland = k; } } /* Now that we got a whole population, print some info, etc., etc... */ printf( "Generation of population %d completed\n", iter ); if( didmigration ) printf( "There was a migration from island %d to %d\n", mIsland + 1, dIsland + 1 ); /* See if we have a single individual overtaking the population(s)... */ singleorg = true; for( k = 0; k <= nIslands - 1; ++k ) { diversity[k] = 1; for( i = 1; i <= popsize[k] - 1; ++i ) { genediff = 0; if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) genediff += fabs( transgene[k][l][rankvec[k][i-1]][j] - transgene[k][l][rankvec[k][i]][j] ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) genediff += fabs( rotgene[k][l][rankvec[k][i-1]][j] - rotgene[k][l][rankvec[k][i]][j] ); } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) genediff += fabs( brotgene[k][rankvec[k][i-1]][j] - brotgene[k][rankvec[k][i]][j] ); } if( genediff >= geneconv ) { singleorg = false; ++diversity[k]; } } } /* If we have islands with identical individuals, just check if the islands are identical by comparing one of the individuals of one island to those in others. Once we find two different ones, break out.. */ if( ( singleorg ) && ( nIslands > 1 ) ) { for( k = 0; k <= nIslands - 2; ++k ) { genediff = 0; if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) genediff += fabs( transgene[k][l][0][j] - transgene[k+1][l][0][j] ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) genediff += fabs( rotgene[k][l][0][j] - rotgene[k+1][l][0][j] ); } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) genediff += fabs( brotgene[k][0][j] - brotgene[k+1][0][j] ); } if( genediff >= geneconv ) { singleorg = false; break; } } } if( singleorg ) { if( nIslands > 1 ) { printf( "Minimum energy (island) = " ); for( k = 0; k <= nIslands - 1; ++k ) { printf( "%g (%d)", gaE[k][rankvec[k][popsize[k]-1]], k + 1 ); if( k < nIslands - 1 ) printf( ", " ); } printf( "\n" ); printf( "Minumum energy (island) = " ); for( k = 0; k <= nIslands - 1; ++k ) { printf( "%g (%d)", gaE[k][rankvec[k][0]], k + 1 ); if( k < nIslands - 1 ) printf( ", " ); } printf( "\n" ); } else { printf( "Maximum energy = %g\n", gaE[0][rankvec[0][popsize[0]-1]] ); printf( "Minimum energy = %g\n", gaE[0][rankvec[0][0]] ); } printf( "\n" ); printf( "The population consists of a single confomer!\n" ); printf( "\n" ); break; } else { if( nIslands > 1 ) { printf( "Minimum energy (island) = " ); for( k = 0; k <= nIslands - 1; ++k ) { printf( "%g (%d)", gaE[k][rankvec[k][popsize[k]-1]], k + 1 ); if( k < nIslands - 1 ) printf( ", " ); } printf( "\n" ); printf( "Minumum energy (island) = " ); for( k = 0; k <= nIslands - 1; ++k ) { printf( "%g (%d)", gaE[k][rankvec[k][0]], k + 1 ); if( k < nIslands - 1 ) printf( ", " ); } printf( "\n" ); printf( "Unique conformers (island) = " ); for( k = 0; k <= nIslands - 1; ++k ) { printf( "%d (%d)", diversity[k], k + 1 ); if( k < nIslands - 1 ) printf( ", " ); } printf( "\n" ); } else { printf( "Maximum energy = %g\n", gaE[0][rankvec[0][popsize[0]-1]] ); printf( "Minimum energy = %g\n", gaE[0][rankvec[0][0]] ); printf( "Unique conformers = %d\n", diversity[0] ); } printf( "\n" ); } /* Check if we have to break out of the endless loop here, otherwise we do an additional mating... */ if( iter >= maxiter ) { printf( "Maximum number of generations reached!\n" ); printf( "\n" ); break; } /* If we want the chemical shift data updated after every generation, do it here... */ if( ! strcmp( verboutput, "YES" ) ) { if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) translate( LigArray[l], nLigAt[l], transgene[bestisland][l][rankvec[bestisland][0]][j], j ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) rotate( LigArray[l], nLigAt[l], rotgene[bestisland][l][rankvec[bestisland][0]][j], j ); } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) rotbond( nFragAt[j], FragArray[j], brotgene[bestisland][rankvec[bestisland][0]][j] ); } /* Now calculate and print the chemical shifts... */ for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( m = 0; m <= nRings - 1; ++m ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[m], nRingAt[m], ifactor[m], method ); } if( nAni > 0 ) { for( m = 0; m <= nAni - 1; ++m ) shift[j] += anishift( protons[j], AniArray[m], nAniAt[m], chiiso[m], chiani[m], bndtyp[m] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; printf( "Proton %d: %s - Experimental shift = %f - Calculated shift = %f\n", protonid[j], protonFullName[j], dexp[j], shift[j] ); } if( nnOe > 0 ) noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "YES" ); printf( "\n" ); /* Now move things back... */ if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) rotbond( nFragAt[j], FragArray[j], -1.0 * brotgene[bestisland][rankvec[bestisland][0]][j] ); } if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) translate( LigArray[l], nLigAt[l], -1.0 * transgene[bestisland][l][rankvec[bestisland][0]][j], j ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 2; j >= 0; --j ) rotate( LigArray[l], nLigAt[l], -1.0 * rotgene[bestisland][l][rankvec[bestisland][0]][j], j ); } } } /* If not there yet, increase the counter... */ ++iter; /* Generate the roulette wheel based on the fitness function... */ for( k = 0; k <= nIslands - 1; ++k ) { pstart[k][0] = 0; pend[k][popsize[k]-1] = 1; for( i = 0; i <= popsize[k] - 2; ++i ) { pstart[k][i+1] = gascore[k][i] / scoresum[k] + pstart[k][i]; pend[k][i] = pstart[k][i+1]; } /* If we want to keep the best individuals, do so here... */ if( keepbest > 0 ) { for( i = 0; i <= keepbest - 1; ++i ) { if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) transchld[k][l][i][j] = transgene[k][l][rankvec[k][i]][j]; } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) rotchld[k][l][i][j] = rotgene[k][l][rankvec[k][i]][j]; } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) brotchld[k][i][j] = brotgene[k][rankvec[k][i]][j]; } } } for( i = keepbest; i <= popsize[k] - 1; ++i ) { /* Select two genes for mating... */ pselect = ( ( double ) random( ) ) / RAND_MAX; for( j = 0; j <= popsize[k] - 1; ++j ) { if( ( pselect >= pstart[k][j] ) && ( pselect < pend[k][j] ) ) { matingA = j; break; } } pselect = ( ( double ) random( ) ) / RAND_MAX; for( j = 0; j <= popsize[k] - 1; ++j ) { if( ( pselect >= pstart[k][j] ) && ( pselect < pend[k][j] ) ) { matingB = j; break; } } /* Do not allow the same organism to mate with itself.... */ while( matingA == matingB ) { pselect = ( ( double ) random( ) ) / RAND_MAX; for( j = 0; j <= popsize[k] - 1; ++j ) { if( ( pselect >= pstart[k][j] ) && ( pselect < pend[k][j] ) ) { matingB = j; break; } } } /* Allow the two selected members of the current population to mate and generate a member of the new population. A simple single-point cross-over is used... */ if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { xpoint = int( ( 4 * ( double ) random( ) ) / RAND_MAX ); for( j = 0; j <= 2; ++j ) { if( j < xpoint ) transchld[k][l][i][j] = transgene[k][l][matingA][j]; else transchld[k][l][i][j] = transgene[k][l][matingB][j]; } } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { xpoint = int( ( 4 * ( double ) random( ) ) / RAND_MAX ); for( j = 0; j <= 2; ++j ) { if( j < xpoint ) rotchld[k][l][i][j] = rotgene[k][l][matingA][j]; else rotchld[k][l][i][j] = rotgene[k][l][matingB][j]; } } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { xpoint = int( ( ( nRotBonds + 1 ) * ( double ) random( ) ) / RAND_MAX ); for( j = 0; j <= nRotBonds - 1; ++j ) { if( j < xpoint ) brotchld[k][i][j] = brotgene[k][matingA][j]; else brotchld[k][i][j] = brotgene[k][matingB][j]; } } /* If we have a mutation frequency defined, check if a mutation is in order and do it... */ if( ( ( ( double ) random( ) ) / RAND_MAX ) < mutfreq ) { mutatethis = 0; whattomutate = int( ( activegavar * ( double ) random( ) ) / RAND_MAX ); if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) { if( whattomutate == mutatethis ) transchld[k][l][i][j] = transwin * delta * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; ++mutatethis; } } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) { if( whattomutate == mutatethis ) rotchld[k][l][i][j] = rotwin * theta * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; ++mutatethis; } } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) { if( whattomutate == mutatethis ) brotchld[k][i][j] = torwin * maxbrot[j] * ( ( double ) ( RAND_MAX - 2 * random( ) ) ) / RAND_MAX; ++mutatethis; } } } } /* Now we have offspring. Put their genes back in the original gene arrays... */ for( i = 0; i <= popsize[k] - 1; ++i ) { if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) transgene[k][l][i][j] = transchld[k][l][i][j]; } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) rotgene[k][l][i][j] = rotchld[k][l][i][j]; } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) brotgene[k][i][j] = brotchld[k][i][j]; } } } } /* Done with all the generations. If we want to save the best ones, pluck them out from the last, which should be the most fit, population... */ if( ( ! strcmp( multconf, "YES" ) ) && ( maxgaconf == 0 ) ) maxgaconf = 10; if( maxgaconf > 0 ) { /* Since we can have multiple islands, the easist thing to do is to make 'index' arrays with the the whole population, and then sort those. Once we have that, go back to the correspongin islands and get the genes, etc., etc... */ prank = ( int* ) malloc( poptotal * sizeof( int ) ); pii = ( int* ) malloc( poptotal * sizeof( int ) ); pin = ( int* ) malloc( poptotal * sizeof( int ) ); popE = ( double* ) malloc( poptotal * sizeof( double ) ); x = 0; for( k = 0; k <= nIslands - 1; ++ k ) { for( i = 0; i <= popsize[k] - 1; ++i ) { prank[x] = x; pii[x] = i; pin[x] = k; popE[x] = gaE[k][i]; ++x; } } /* Again, do a bubble-sort to get the ranks... */ for( j = 0; j <= poptotal - 1; ++j ) { for( i = 0; i <= poptotal - 2 - j; ++i ) { if( popE[prank[i+1]] < popE[prank[i]] ) { temprank = prank[i]; prank[i] = prank[i+1]; prank[i+1] = temprank; } } } n = 1; x = 1; while( ( n <= maxgaconf ) && ( x <= poptotal - 1 ) ) { /* Don't save anything if they are all the same. Compare the genes... */ genediff = 0; if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) genediff += fabs( transgene[pin[prank[x-1]]][l][pii[prank[x-1]]][j] - transgene[pin[prank[x]]][l][pii[prank[x]]][j] ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) genediff += fabs( rotgene[pin[prank[x-1]]][l][pii[prank[x-1]]][j] - rotgene[pin[prank[x]]][l][pii[prank[x]]][j] ); } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) genediff += fabs( brotgene[pin[prank[x-1]]][pii[prank[x-1]]][j] - brotgene[pin[prank[x]]][pii[prank[x]]][j] ); } if( genediff >= geneconv ) { /* If we find a unique conformer, save it. First get the 'gene product'... */ if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) translate( LigArray[l], nLigAt[l], transgene[pin[prank[x-1]]][l][pii[prank[x-1]]][j], j ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) rotate( LigArray[l], nLigAt[l], rotgene[pin[prank[x-1]]][l][pii[prank[x-1]]][j], j ); } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) rotbond( nFragAt[j], FragArray[j], brotgene[pin[prank[x-1]]][pii[prank[x-1]]][j] ); } /* Now save it... */ sprintf( conffile, "%s%s%03d", outfile, ".", n ); if( writemol( ftype, conffile, HeaderCount, HeaderJunk, AtomCount, RawAtomData, nLigands, nLigAt, LigAtID, LigArray, FooterCount, FooterJunk, popE[prank[x-1]] ) ) { printf( "Error opening output file %s. Check file permissions and try again...\n", conffile ); return 0; } /* And finally move things back... */ if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) rotbond( nFragAt[j], FragArray[j], -1.0 * brotgene[pin[prank[x-1]]][pii[prank[x-1]]][j] ); } if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) translate( LigArray[l], nLigAt[l], -1.0 * transgene[pin[prank[x-1]]][l][pii[prank[x-1]]][j], j ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 2; j >= 0; --j ) rotate( LigArray[l], nLigAt[l], -1.0 * rotgene[pin[prank[x-1]]][l][pii[prank[x-1]]][j], j ); } } ++n; } ++x; } } /* Finally, save the best conformer from the last generation in the ligand array... */ if( ! strcmp( dotrans, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) translate( LigArray[l], nLigAt[l], transgene[bestisland][l][rankvec[bestisland][0]][j], j ); } } if( ! strcmp( dorot, "YES" ) ) { for( l = 0; l <= nLigands - 1; ++l ) { for( j = 0; j <= 2; ++j ) rotate( LigArray[l], nLigAt[l], rotgene[bestisland][l][rankvec[bestisland][0]][j], j ); } } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { for( j = 0; j <= nRotBonds - 1; ++j ) rotbond( nFragAt[j], FragArray[j], brotgene[bestisland][rankvec[bestisland][0]][j] ); } } /* Once we are done, print out a table with the final chemical shifts and statistics, unless we just wanted a report... */ if( strcmp( optmethod, "RPRT" ) ) { AvgDeltaDiff = 0; AvgDelta2Diff = 0; AvgDeltaDiff2 = 0; AvgDeltaSD = 0; ffE = 0; shiftE = 0; nOeE = 0; selfE = 0; torE = 0; anchorE = 0; eterms = 0; if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) { ffE = nonbonded( LigArray, nLigAt, nLigands, ProtArray, nProtAt, cutoff, useel, diel, dielfunc, racemate ); ++eterms; } if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) ++eterms; for( j = 0; j <= nProtons - 1; ++j ) { shift[j] = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shift[j] += pertfact[j] * ringshift( protons[j], RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shift[j] += anishift( protons[j], AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( nEleAt > 0 ) shift[j] += electroshift( protons[j], protoneigh[j], ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( pName[j], "H" ) ) ) shift[j] += hbshift( protons[j], protoneigh[j], HBSlope, HBConst, nHBAcc, HBAccAtoms ); if( ! strcmp( racemate, "YES" ) ) shift[j] /= nLigands; if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) shiftE += force[j] * ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); printf( "Proton %d: %s - Experimental shift = %f - Calculated shift = %f\n", protonid[j], protonFullName[j], dexp[j], shift[j] ); /* Compute statistics... */ AvgDeltaDiff += fabs( shift[j] - dexp[j] ); AvgDelta2Diff += ( shift[j] - dexp[j] ) * ( shift[j] - dexp[j] ); } AvgDeltaDiff /= nProtons; AvgDelta2Diff /= nProtons; AvgDeltaDiff2 = AvgDeltaDiff * AvgDeltaDiff; AvgDeltaSD = sqrt( fabs( AvgDelta2Diff - AvgDeltaDiff2 ) ); printf( "\n" ); printf( "Final shift deviation = %f +/- %f\n", AvgDeltaDiff, AvgDeltaSD ); if( nnOe > 0 ) { nOeE = noepot( nOeAtLig, nOeAtLigID, nOeAtProt, nOeAtProtID, nnOe, nOeRange, nOePenalty, "END" ); ++eterms; } if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { selfE = slfligene( LigArray, nLigAt, nLigands, LigNBList, nLigNBList, cutoff, useel, diel, dielfunc ); torE = gettorpot( nRotBonds, nTorsions, TorPot, TorPhase, TorArray ); eterms +=2; } if( ! strcmp( constcent, "YES" ) ) { anchorE = centpot( nLigands, LigArray, nLigAt, initcent, radiuscp, forcecp ); ++eterms; } finalE = nOeE + shiftE + ffE + selfE + torE + anchorE; /* If we have more than one energy term contributing to the total energy, print out a breakdown of the energy contributions... */ if( eterms > 1 ) { printf( "\n" ); if( ( ! strcmp( potential, "FF" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) printf( "Final force-field energy contribution = %f Kcal/mol\n", ffE ); if( ( ! strcmp( potential, "NMR" ) ) || ( ! strcmp( potential, "FFNMR" ) ) ) printf( "Final shift penalty energy contribution = %f Kcal/mol\n", shiftE ); if( nnOe > 0 ) printf( "Final nOe penalty energy contribution = %f Kcal/mol\n", nOeE ); if( ( nRotBonds > 0 ) && ( ! strcmp( dobonds, "YES" ) ) ) { printf( "Final self-energy contribution = %f Kcal/mol\n", selfE ); printf( "Final torsional energy contribution = %f Kcal/mol\n", torE ); } if( ! strcmp( constcent, "YES" ) ) printf( "Final anchor potential contribution = %f Kcal/mol\n", anchorE ); } printf( "\n" ); printf( "Final energy of the system = %f Kcal/mol\n", ( nOeE + shiftE + ffE + selfE + torE + anchorE ) ); /* We are done. After all this crap we should save the new mol2/pdb file. First delete the output file, no questions asked... */ if( writemol( ftype, outfile, HeaderCount, HeaderJunk, AtomCount, RawAtomData, nLigands, nLigAt, LigAtID, LigArray, FooterCount, FooterJunk, finalE ) ) { printf( "Error opening output file %s. Check file permissions and try again...\n", outfile ); return 0; } } /* If we chose to get a 'calculated' .sdl file for use with jsurf, do it here for all protons of the protein. If we are doing electrostatics we have to recalculate the neighbor table, because the one we have is only for the protons from the observed .sdl file. This is done only if we are taking electrostatic effects into account... */ if( strcmp( calcsdl, "" ) ) { CalcSdl = fopen( calcsdl, "w" ); if( ! CalcSdl ) { printf( "Error opening output file %s. Check file permissions and try again...\n", calcsdl ); return 0; } /* We need to print the ligand section of the .sdl file in case we want to filter by vdW radii during the j-surface calculation... */ fprintf( CalcSdl, "@LIGAND\n" ); for( i = 0; i <= nLigands - 1; ++i ) { fprintf( CalcSdl, "%s ", ligName[i] ); for( j = 0; j <= nLigAt[i] - 1; ++j ) { fprintf( CalcSdl, "%d", LigAtID[i][j] ); if( j < ( nLigAt[i] - 1 ) ) fprintf( CalcSdl, " " ); } fprintf( CalcSdl, "\n" ); } fprintf( CalcSdl, "@PROTONS\n" ); for( k = 0; k <= nProtAt - 1; ++k ) { if( ( ! strcmp( Protattype[k], "H" ) ) || ( ! strcmp( ProtAtName[k], "C" ) ) || ( ! strcmp( ProtAtName[k], "N" ) ) || ( ! strcmp( ProtAtName[k], "CA" ) ) || ( ! strcmp( ProtAtName[k], "CB" ) ) ) { if( ! strcmp( ProtAtName[k], "H" ) ) pfact = 1.376; else if( ( ! strcmp( ProtAtName[k], "N" ) ) || ( ! strcmp( ProtAtName[k], "C" ) ) || ( ! strcmp( ProtAtName[k], "CB" ) ) ) pfact = 0.195; else if( ! strcmp( ProtAtName[k], "CA" ) ) pfact = 0.292; else pfact = 1.000; for( j = 0; j <= 2; ++j) ProtH[j] = ProtArray[k][j]; shiftval = 0; if( nRings > 0 ) { for( i = 0; i <= nRings - 1; ++i ) shiftval += pfact * ringshift( ProtH, RingArray[i], nRingAt[i], ifactor[i], method ); } if( nAni > 0 ) { for( i = 0; i <= nAni - 1; ++i ) shiftval += anishift( ProtH, AniArray[i], nAniAt[i], chiiso[i], chiani[i], bndtyp[i] ); } if( ( nEleAt > 0 ) || ( nHBAcc > 0 ) ) { mindistance = INFINITY; for( j = 0; j <= nProtAt - 1; ++j ) { if( j != k ) { rPRPA = sqrt( ( ProtH[0] - ProtArray[j][0] ) * ( ProtH[0] - ProtArray[j][0] ) + ( ProtH[1] - ProtArray[j][1] ) * ( ProtH[1] - ProtArray[j][1] ) + ( ProtH[2] - ProtArray[j][2] ) * ( ProtH[2] - ProtArray[j][2] ) ); if( rPRPA < mindistance ) { neighindex = j; mindistance = rPRPA; } } } for( j = 0; j <= 2; ++j) ProtHneigh[j] = ProtArray[neighindex][j]; if( nEleAt > 0 ) shiftval += electroshift( ProtH, ProtHneigh, ElecAtoms, nEleAt, charges, eleA, eleB, diel, dielfunc ); if( ( nHBAcc > 0 ) && ( ! strcmp( ProtAtName[k], "H" ) ) ) shiftval += hbshift( ProtH, ProtHneigh, HBSlope, HBConst, nHBAcc, HBAccAtoms ); } if( ! strcmp( racemate, "YES" ) ) shiftval /= nLigands; if( fabs( shiftval ) >= jthresh ) { if( ! strcmp( jptype, "ALLA" ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); else if( ( ! strcmp( jptype, "ALLH" ) ) && ( ! strcmp( Protattype[k], "H" ) ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); else if( ( ! strcmp( jptype, "BCKA" ) ) && ( ( ! strcmp( ProtAtName[k], "H" ) ) || ( ! strcmp( ProtAtName[k], "HA" ) ) || ( ! strcmp( ProtAtName[k], "N" ) ) || ( ! strcmp( ProtAtName[k], "CA" ) ) ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); else if( ( ! strcmp( jptype, "BCKH" ) ) && ( ( ! strcmp( ProtAtName[k], "H" ) ) || ( ! strcmp( ProtAtName[k], "HA" ) ) ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); else if( ( ! strcmp( jptype, "ALKH" ) ) && ( ( ! strcmp( Protattype[k], "H" ) ) && ( strcmp( ProtAtName[k], "H" ) ) ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); else if( ( ! strcmp( jptype, "AMDH" ) ) && ( ! strcmp( ProtAtName[k], "H" ) ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); else if( ( ! strcmp( jptype, "NATM" ) ) && ( ! strcmp( ProtAtName[k], "N" ) ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); else if( ( ! strcmp( jptype, "CAAT" ) ) && ( ! strcmp( ProtAtName[k], "CA" ) ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); else if( ( ! strcmp( jptype, "CBAT" ) ) && ( ! strcmp( ProtAtName[k], "CB" ) ) ) fprintf( CalcSdl, "%d\t%f\t%f\t%f\t%s\t%s\n", ProtAtID[k], shiftval, 0.00, pfact, ProtAtName[k], ProtAtResName[k] ); } } } } printf( "\n" ); } /* >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Program Ends Here <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< */