1 /********************************************************************** 2 forcefield.cpp - A OBOp to calculate and minimize the energy using a 3 forcefield (re-wrap of obminimize and obenergy) 4 5 Copyright (C) 2006-2007 by Tim Vandermeersch 6 (C) 2009 by Chris Morley 7 8 This file is part of the Open Babel project. 9 For more information, see <http://openbabel.org/> 10 11 This program is free software; you can redistribute it and/or modify 12 it under the terms of the GNU General Public License as published by 13 the Free Software Foundation version 2 of the License. 14 15 This program is distributed in the hope that it will be useful, 16 but WITHOUT ANY WARRANTY; without even the implied warranty of 17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18 GNU General Public License for more details. 19 ***********************************************************************/ 20 21 /****************************************************************************** 22 **** CURRENTLY ONLY SUITABLE FOR USE WITH THE OBABEL COMMANDLINE INTERFACE **** 23 This allows options to have parameters, e.g. --ff Ghemical 24 Compile with tools/obabel.cpp rather than tools/babel.cpp 25 26 *******************************************************************************/ 27 28 #include <openbabel/babelconfig.h> 29 #include <iostream> 30 #include <cstdlib> 31 #include<openbabel/op.h> 32 #include<openbabel/mol.h> 33 #include<openbabel/forcefield.h> 34 #include<openbabel/generic.h> 35 36 namespace OpenBabel 37 { 38 using namespace std; 39 40 ////////////////////////////////////////////////////////// 41 // 42 // OpEnergy 43 // 44 ////////////////////////////////////////////////////////// 45 46 class OpEnergy : public OBOp 47 { 48 public: OpEnergy(const char * ID)49 OpEnergy(const char *ID) : OBOp(ID, false) {} 50 Description()51 const char* Description() 52 { 53 return "ForceField Energy Evaluation (not displayed in GUI)\n" 54 "Typical usage: obabel infile.xxx -O outfile.yy --energy --log\n" 55 " options: description\n" 56 " --log output a log of the energies (default = no log)\n" 57 " --epsilon # set the dielectric constant for electrostatics\n" 58 " --noh don't add explicit hydrogens (default = make explicit)\n" 59 " --ff # select a forcefield (default = MMFF94)\n" 60 " The hydrogens are made explicit by default before energy evaluation.\n" 61 " The energy is put in an OBPairData object \"Energy\" which is\n" 62 " accessible via an SDF or CML property or --append (to title).\n" 63 ; 64 } 65 WorksWith(OBBase * pOb) const66 virtual bool WorksWith(OBBase* pOb) const 67 { 68 return dynamic_cast<OBMol*>(pOb) != nullptr; 69 } 70 virtual bool Do(OBBase* pOb, const char* OptionText, OpMap* pmap, OBConversion*); 71 }; 72 73 ////////////////////////////////////////////////////////// 74 OpEnergy theOpEnergy("energy"); //Global instance 75 76 ////////////////////////////////////////////////////////// Do(OBBase * pOb,const char * OptionText,OpMap * pmap,OBConversion *)77 bool OpEnergy::Do(OBBase* pOb, const char* OptionText, OpMap* pmap, OBConversion*) 78 { 79 OBMol* pmol = dynamic_cast<OBMol*>(pOb); 80 if(!pmol) 81 return false; 82 83 bool log = false; 84 bool addh = true; 85 86 string ff = "MMFF94"; 87 double epsilon = 1.0; 88 OpMap::const_iterator iter = pmap->find("ff"); 89 if(iter!=pmap->end()) 90 ff = iter->second; 91 OBForceField* pFF = OBForceField::FindForceField(ff); 92 iter = pmap->find("epsilon"); 93 if (iter!=pmap->end()) 94 epsilon = atof(iter->second.c_str()); 95 96 iter = pmap->find("log"); 97 if(iter!=pmap->end()) 98 log=true; 99 100 iter = pmap->find("noh"); 101 if(iter!=pmap->end()) 102 addh=false; 103 104 if (addh) 105 pmol->AddHydrogens(false, false); 106 107 // set some force field variables 108 pFF->SetLogFile(&clog); 109 pFF->SetLogLevel(log ? OBFF_LOGLVL_MEDIUM : OBFF_LOGLVL_NONE); 110 111 pFF->SetDielectricConstant(epsilon); 112 if (!pFF->Setup(*pmol)) { 113 cerr << "Could not setup force field." << endl; 114 return false; 115 } 116 117 //Put the energy in a OBPairData object 118 OBPairData *dp = new OBPairData; 119 dp->SetAttribute("Energy"); 120 stringstream ss; 121 ss << pFF->Energy(false); 122 dp->SetValue(ss.str()); 123 dp->SetOrigin(fileformatInput); 124 pmol->SetData(dp); 125 126 return true; 127 } 128 129 ////////////////////////////////////////////////////////// 130 // 131 // OpMinimize 132 // 133 ////////////////////////////////////////////////////////// 134 135 class OpMinimize : public OBOp 136 { 137 public: OpMinimize(const char * ID)138 OpMinimize(const char* ID) : OBOp(ID, false) {} 139 Description()140 const char* Description() 141 { 142 return "ForceField Energy Minimization (not displayed in GUI)\n" 143 "Typical usage: obabel infile.xxx -O outfile.yyy --minimize --steps 1500 --sd\n" 144 " options: description:\n" 145 " --log output a log of the minimization process(default= no log)\n" 146 " --crit # set convergence criteria (default=1e-6)\n" 147 " --sd use steepest descent algorithm (default = conjugate gradient)\n" 148 " --newton use Newton2Num linesearch (default = Simple)\n" 149 " --ff # select a forcefield (default = Ghemical)\n" 150 " --steps # specify the maximum number of steps (default = 2500)\n" 151 " --cut use cut-off (default = don't use cut-off)\n" 152 " --noh don't add explicit hydrogens (default = make explicit)\n" 153 " --epsilon # relative dielectric constant (default = 1.0)\n" 154 " --rvdw # specify the VDW cut-off distance (default = 6.0)\n" 155 " --rele # specify the Electrostatic cut-off distance (default = 10.0)\n" 156 " --freq # specify the frequency to update the non-bonded pairs (default = 10)\n" 157 " The hydrogens are made explicit before minimization by default.\n" 158 " The energy is put in an OBPairData object \"Energy\" which is\n" 159 " accessible via an SDF or CML property or --append (to title).\n" 160 ; 161 } 162 WorksWith(OBBase * pOb) const163 virtual bool WorksWith(OBBase* pOb) const 164 { 165 return dynamic_cast<OBMol*>(pOb) != nullptr; 166 } 167 virtual bool Do(OBBase* pOb, const char* OptionText, OpMap* pmap, OBConversion*); 168 }; 169 170 ////////////////////////////////////////////////////////// 171 OpMinimize theOpMinimize("minimize"); //Global instance 172 173 ////////////////////////////////////////////////////////// Do(OBBase * pOb,const char * OptionText,OpMap * pmap,OBConversion *)174 bool OpMinimize::Do(OBBase* pOb, const char* OptionText, OpMap* pmap, OBConversion*) 175 { 176 OBMol* pmol = dynamic_cast<OBMol*>(pOb); 177 if(!pmol) 178 return false; 179 180 int steps = 2500; 181 double crit = 1e-6; 182 bool sd = false; 183 bool cut = false; 184 bool addh = true; 185 bool newton = true; 186 double epsilon = 1.0; 187 double rvdw = 6.0; 188 double rele = 10.0; 189 int freq = 10; 190 bool log = false; 191 192 string ff = "MMFF94"; 193 OpMap::const_iterator iter = pmap->find("ff"); 194 if(iter!=pmap->end()) 195 ff = iter->second; 196 OBForceField* pFF = OBForceField::FindForceField(ff); 197 198 iter = pmap->find("sd"); 199 if(iter!=pmap->end()) 200 sd=true; 201 202 iter = pmap->find("newton"); 203 if(iter!=pmap->end()) 204 newton=true; 205 206 iter = pmap->find("cut"); 207 if(iter!=pmap->end()) 208 cut=true; 209 210 iter = pmap->find("noh"); 211 if(iter!=pmap->end()) 212 addh=false; 213 214 iter = pmap->find("crit"); 215 if(iter!=pmap->end()) 216 crit = atof(iter->second.c_str()); 217 218 iter = pmap->find("steps"); 219 if(iter!=pmap->end()) 220 steps = atoi(iter->second.c_str()); 221 222 iter = pmap->find("epsilon"); 223 if(iter!=pmap->end()) 224 epsilon = atof(iter->second.c_str()); 225 226 iter = pmap->find("rvdw"); 227 if(iter!=pmap->end()) 228 rvdw = atof(iter->second.c_str()); 229 230 iter = pmap->find("rele"); 231 if(iter!=pmap->end()) 232 rele = atof(iter->second.c_str()); 233 234 iter = pmap->find("pf"); 235 if(iter!=pmap->end()) { 236 freq = atoi(iter->second.c_str()); 237 if (freq < 1) 238 freq = 10; // don't divide by zero 239 } 240 241 iter = pmap->find("log"); 242 if(iter!=pmap->end()) 243 log=true; 244 245 if (newton) 246 pFF->SetLineSearchType(LineSearchType::Newton2Num); 247 248 // set some force field variables 249 pFF->SetLogFile(&clog); 250 pFF->SetLogLevel(log ? OBFF_LOGLVL_LOW : OBFF_LOGLVL_NONE); 251 pFF->SetVDWCutOff(rvdw); 252 pFF->SetElectrostaticCutOff(rele); 253 pFF->SetUpdateFrequency(freq); 254 pFF->SetDielectricConstant(epsilon); 255 pFF->EnableCutOff(cut); 256 257 if (addh) 258 pmol->AddHydrogens(false, false); 259 260 if (!pFF->Setup(*pmol)) { 261 cerr << "Could not setup force field." << endl; 262 return false; 263 } 264 265 bool done = true; 266 if (sd) 267 pFF->SteepestDescent(steps, crit); 268 else 269 pFF->ConjugateGradients(steps, crit); 270 271 pFF->GetCoordinates(*pmol); 272 273 //Put the energy in a OBPairData object 274 OBPairData *dp = new OBPairData; 275 dp->SetAttribute("Energy"); 276 stringstream ss; 277 ss << pFF->Energy(false); 278 dp->SetValue(ss.str()); 279 dp->SetOrigin(fileformatInput); 280 pmol->SetData(dp); 281 282 return true; 283 } 284 285 }//namespace 286