Merge pull request #2022 from ghutchis/ff-dielectric-const

Add force field support for dielectric constants in charge terms.
openbabel · Sep 14, 2019 · c195713 · c195713
2 parents 62bc9b4 + 166b7bd
commit c195713
Show file tree

Hide file tree

Showing 16 changed files with 132 additions and 27 deletions.
diff --git a/include/openbabel/forcefield.h b/include/openbabel/forcefield.h
@@ -554,6 +554,7 @@ namespace OpenBabel
     bool 	_cutoff; //!< true = cut-off enabled
     double 	_rvdw; //!< VDW cut-off distance
     double 	_rele; //!< Electrostatic cut-off distance
+    double _epsilon; //!< Dielectric constant for electrostatics
     OBBitVec	_vdwpairs; //!< VDW pairs that should be calculated
     OBBitVec	_elepairs; //!< Electrostatic pairs that should be calculated
     int 	_pairfreq; //!< The frequence to update non-bonded pairs
@@ -817,6 +818,20 @@ namespace OpenBabel
     {
       return _rele;
     }
+    /*! Set the dielectric constant for electrostatic SetupCalculations
+     * \param epsilon The relative permittivity to use (default = 1.0)
+     */
+     void SetDielectricConstant(double epsilon)
+     {
+       _epsilon = epsilon;
+     }
+     /* Get the dielectric permittivity used for electrostatic calculations
+     * \rreturn The current relative permittivity
+     */
+     double GetDielectricConstant()
+     {
+       return _epsilon;
+     }
     /*! Set the frequency by which non-bonded pairs are updated. Values from 10 to 20
      *  are recommended. Too low will decrease performance, too high will cause
      *  non-bonded interactions within cut-off not to be calculated.

diff --git a/src/forcefields/forcefieldgaff.cpp b/src/forcefields/forcefieldgaff.cpp
@@ -1040,7 +1040,7 @@ namespace OpenBabel
           continue;
       }
 
-      elecalc.qq = KCAL_TO_KJ * 332.17 * a->GetPartialCharge() * b->GetPartialCharge();
+      elecalc.qq = KCAL_TO_KJ * 332.17 * a->GetPartialCharge() * b->GetPartialCharge() / _epsilon;
 
       if (elecalc.qq) {
         elecalc.a = &*a;

diff --git a/src/forcefields/forcefieldgaff.h b/src/forcefields/forcefieldgaff.h
@@ -126,6 +126,7 @@ namespace OpenBabel
         _init = false;
         _rvdw = 7.0;
         _rele = 15.0;
+        _epsilon = 1.0;
         _pairfreq = 10;
         _cutoff = false;
         _linesearch = LineSearchType::Newton2Num;
@@ -209,4 +210,3 @@ namespace OpenBabel
 
 //! \file forcefieldgaff.h
 //! \brief Gaff force field
-
diff --git a/src/forcefields/forcefieldghemical.cpp b/src/forcefields/forcefieldghemical.cpp
@@ -815,7 +815,7 @@ namespace OpenBabel
           continue;
       }
 
-      elecalc.qq = KCAL_TO_KJ * 332.17 * a->GetPartialCharge() * b->GetPartialCharge();
+      elecalc.qq = KCAL_TO_KJ * 332.17 * a->GetPartialCharge() * b->GetPartialCharge() / _epsilon;
 
       if (elecalc.qq) {
         elecalc.a = &*a;

diff --git a/src/forcefields/forcefieldghemical.h b/src/forcefields/forcefieldghemical.h
@@ -117,6 +117,7 @@ namespace OpenBabel
         _init = false;
         _rvdw = 7.0;
         _rele = 15.0;
+        _epsilon = 1.0;
         _pairfreq = 10;
         _cutoff = false;
         _linesearch = LineSearchType::Newton2Num;
@@ -194,4 +195,3 @@ namespace OpenBabel
 
 //! \file forcefieldghemical.h
 //! \brief Ghemical force field
-
diff --git a/src/forcefields/forcefieldmm2.h b/src/forcefields/forcefieldmm2.h
@@ -98,7 +98,7 @@ namespace OpenBabel
 
 
       virtual const char* Description()
-	{ return "MM2 force field.";};
+      { return "MM2 force field.";};
 
 
       //! Destructor
@@ -129,4 +129,3 @@ namespace OpenBabel
 
 //! \file forcefieldmm2.h
 //! \brief MM2 force field
-
diff --git a/src/forcefields/forcefieldmmff94.cpp b/src/forcefields/forcefieldmmff94.cpp
@@ -3656,7 +3656,7 @@ namespace OpenBabel
           continue;
       }
 
-      elecalc.qq = 332.0716 * a->GetPartialCharge() * b->GetPartialCharge();
+      elecalc.qq = 332.0716 * a->GetPartialCharge() * b->GetPartialCharge() / _epsilon;
 
       if (elecalc.qq) {
         elecalc.a = &*a;

diff --git a/src/forcefields/forcefieldmmff94.h b/src/forcefields/forcefieldmmff94.h
@@ -235,6 +235,7 @@ namespace OpenBabel
         _init = false;
         _rvdw = 7.0;
         _rele = 15.0;
+        _epsilon = 1.0; // default electrostatics
         _pairfreq = 15;
         _cutoff = false;
         _linesearch = LineSearchType::Newton2Num;
@@ -335,4 +336,3 @@ namespace OpenBabel
 
 //! \file forcefieldmmff94.h
 //! \brief MMFF94 force field
-
diff --git a/src/forcefields/forcefielduff.h b/src/forcefields/forcefielduff.h
@@ -127,6 +127,7 @@ namespace OpenBabel
       _init = false;
       _rvdw = 7.0;
       _rele = 15.0;
+      _epsilon = 1.0; // electrostatics not used
       _pairfreq = 10;
       _cutoff = false;
       _linesearch = LineSearchType::Newton2Num;

diff --git a/src/ops/forcefield.cpp b/src/ops/forcefield.cpp
@@ -54,7 +54,8 @@ namespace OpenBabel
           "Typical usage: obabel infile.xxx -O outfile.yy --energy --log\n"
           " options:         description\n"
           " --log        output a log of the energies (default = no log)\n"
-          " --ff #       select a forcefield (default = Ghemical)\n"
+          " --epsilon #  set the dielectric constant for electrostatics\n"
+          " --ff #       select a forcefield (default = MMFF94)\n"
           " The hydrogens are always made explicit before energy evaluation.\n"
           " The energy is put in an OBPairData object \"Energy\" which is\n"
           "   accessible via an SDF or CML property or --append (to title).\n"
@@ -82,10 +83,14 @@ namespace OpenBabel
     bool log = false;
 
     string ff = "MMFF94";
+    double epsilon = 1.0;
     OpMap::const_iterator iter = pmap->find("ff");
     if(iter!=pmap->end())
       ff = iter->second;
     OBForceField* pFF = OBForceField::FindForceField(ff);
+    iter = pmap->find("epsilon");
+    if (iter!=pmap->end())
+      epsilon = atof(iter->second.c_str());
 
     iter = pmap->find("log");
     if(iter!=pmap->end())
@@ -95,6 +100,7 @@ namespace OpenBabel
     pFF->SetLogFile(&clog);
     pFF->SetLogLevel(log ? OBFF_LOGLVL_MEDIUM : OBFF_LOGLVL_NONE);
 
+    pFF->SetDielectricConstant(epsilon);
     if (!pFF->Setup(*pmol)) {
       cerr  << "Could not setup force field." << endl;
       return false;
@@ -135,6 +141,7 @@ namespace OpenBabel
           " --ff #       select a forcefield (default = Ghemical)\n"
           " --steps #    specify the maximum number of steps (default = 2500)\n"
           " --cut        use cut-off (default = don't use cut-off)\n"
+          " --epsilon #  relative dielectric constant (default = 1.0)\n"
           " --rvdw #     specify the VDW cut-off distance (default = 6.0)\n"
           " --rele #     specify the Electrostatic cut-off distance (default = 10.0)\n"
           " --freq #     specify the frequency to update the non-bonded pairs (default = 10)\n"
@@ -167,6 +174,7 @@ namespace OpenBabel
     bool sd = false;
     bool cut = false;
     bool newton = true;
+    double epsilon = 1.0;
     double rvdw = 6.0;
     double rele = 10.0;
     int freq = 10;
@@ -198,6 +206,10 @@ namespace OpenBabel
     if(iter!=pmap->end())
       steps = atoi(iter->second.c_str());
 
+    iter = pmap->find("epsilon");
+    if(iter!=pmap->end())
+      epsilon = atof(iter->second.c_str());
+
     iter = pmap->find("rvdw");
     if(iter!=pmap->end())
       rvdw = atof(iter->second.c_str());
@@ -226,6 +238,7 @@ namespace OpenBabel
     pFF->SetVDWCutOff(rvdw);
     pFF->SetElectrostaticCutOff(rele);
     pFF->SetUpdateFrequency(freq);
+    pFF->SetDielectricConstant(epsilon);
     pFF->EnableCutOff(cut);
 
     if (!pFF->Setup(*pmol)) {

diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
@@ -71,7 +71,7 @@ set(origtests
     unitcell
     )
 set (atom_parts 1 2 3 4)
-set (ffmmff94_parts 1 2)
+set (ffmmff94_parts 1 2 3 4 5 6)
 set (math_parts 1 2 3 4)
 set (pdbreadfile_parts 1 2 3 4)
 

diff --git a/test/ffmmff94.cpp b/test/ffmmff94.cpp
@@ -5,11 +5,11 @@ This file is part of the Open Babel project.
 For more information, see <http://openbabel.org/>
 
 Some portions Copyright (C) 2008 Geoffrey R. Hutchison
- 
+
 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 version 2 of the License.
- 
+
 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
@@ -40,7 +40,7 @@ using namespace OpenBabel;
 
   int currentTest = 0;
 
-void GenerateEnergies(string molecules_file, string results_file)
+void GenerateEnergies(string molecules_file, string results_file, string method, double epsilon = 1.0)
 {
   std::ifstream ifs;
   if (!SafeOpen(ifs, molecules_file.c_str()))
@@ -53,18 +53,19 @@ void GenerateEnergies(string molecules_file, string results_file)
   OBMol mol;
   OBConversion conv(&ifs, &cout);
   char buffer[BUFF_SIZE];
-  
+
   if(! conv.SetInAndOutFormats("SDF","SDF"))
     {
       cerr << "SDF format is not loaded" << endl;
       return;
     }
 
-  OBForceField* pFF = OBForceField::FindForceField("MMFF94");
+  OBForceField* pFF = OBForceField::FindForceField(method);
   OB_REQUIRE(pFF != NULL);
 
   pFF->SetLogFile(&cout);
   pFF->SetLogLevel(OBFF_LOGLVL_NONE);
+  pFF->SetDielectricConstant(epsilon);
 
   for (;ifs;)
     {
@@ -77,7 +78,7 @@ void GenerateEnergies(string molecules_file, string results_file)
         cerr << "Could not setup force field on molecule: " << mol.GetTitle() << endl;
         return;
       }
-      
+
       // Don't compute gradients
       sprintf(buffer, "%15.5f\n", pFF->Energy(false));
       ofs << buffer;
@@ -87,7 +88,7 @@ void GenerateEnergies(string molecules_file, string results_file)
   return;
 }
 
-void TestFile(string filename, string results_file)
+void TestFile(string filename, string results_file, string method, double epsilon = 1.0)
 {
   std::ifstream mifs;
   if (!SafeOpen(mifs, filename.c_str()))
@@ -113,12 +114,13 @@ void TestFile(string filename, string results_file)
       cout << "Bail out! SDF format is not loaded" << endl;
       return;
     }
-    
-  OBForceField* pFF = OBForceField::FindForceField("MMFF94");
+
+  OBForceField* pFF = OBForceField::FindForceField(method);
   OB_REQUIRE(pFF != NULL);
 
   pFF->SetLogFile(&cout);
   pFF->SetLogLevel(OBFF_LOGLVL_NONE);
+  pFF->SetDielectricConstant(epsilon);
 
   double energy;
   while(mifs)
@@ -132,7 +134,7 @@ void TestFile(string filename, string results_file)
           cout << "Bail out! error reading reference data" << endl;
           return;
         }
-        
+
       if (!pFF->Setup(mol)) {
         cout << "Bail out! could not setup force field on " << mol.GetTitle() << endl;
         return;
@@ -164,7 +166,7 @@ void TestFile(string filename, string results_file)
 int ffmmff94(int argc, char* argv[])
 {
   int defaultchoice = 1;
-  
+
   int choice = defaultchoice;
 
   if (argc > 1) {
@@ -174,22 +176,46 @@ int ffmmff94(int argc, char* argv[])
     }
   }
 
+  string testdatadir = TESTDATADIR;
+
+  if (choice == 99)
+    {
+      GenerateEnergies(testdatadir + "forcefield.sdf", testdatadir + "mmff94results.txt", "MMFF94");
+      GenerateEnergies(testdatadir + "more-mmff94.sdf", testdatadir + "more-mmff94results.txt", "MMFF94"); // provided by Paolo Tosco
+      GenerateEnergies(testdatadir + "forcefield.sdf", testdatadir + "mmff94sresults.txt", "MMFF94s");
+      GenerateEnergies(testdatadir + "more-mmff94.sdf", testdatadir + "more-mmff94sresults.txt", "MMFF94s"); // ditto
+      GenerateEnergies(testdatadir + "forcefield.sdf", testdatadir + "mmff94e4results.txt", "MMFF94", 4.0);
+      GenerateEnergies(testdatadir + "more-mmff94.sdf", testdatadir + "more-mmff94e4results.txt", "MMFF94", 4.0); // provided by Paolo Tosco
+
+      return 0;
+    }
+
   // Define location of file formats for testing
   #ifdef FORMATDIR
     char env[BUFF_SIZE];
     snprintf(env, BUFF_SIZE, "BABEL_LIBDIR=%s", FORMATDIR);
     putenv(env);
-  #endif  
-
-  string testdatadir = TESTDATADIR;
+  #endif
 
   cout << "# Testing MMFF94 Force Field..." << endl;
   switch(choice) {
   case 1:
-    TestFile(testdatadir + "forcefield.sdf", testdatadir + "mmff94results.txt");
+    TestFile(testdatadir + "forcefield.sdf", testdatadir + "mmff94results.txt", "MMFF94");
     break;
   case 2:
-    TestFile(testdatadir + "more-mmff94.sdf", testdatadir + "more-mmff94results.txt"); // provided by Paolo Tosco
+    TestFile(testdatadir + "more-mmff94.sdf", testdatadir + "more-mmff94results.txt", "MMFF94"); // provided by Paolo Tosco
+    break;
+  case 3:
+    TestFile(testdatadir + "forcefield.sdf", testdatadir + "mmff94sresults.txt", "MMFF94s");
+    break;
+  case 4:
+    TestFile(testdatadir + "more-mmff94.sdf", testdatadir + "more-mmff94sresults.txt", "MMFF94s"); // ditto
+    break;
+  case 5:
+    TestFile(testdatadir + "forcefield.sdf", testdatadir + "mmff94e4results.txt", "MMFF94", 4.0);
+    break;
+  case 6:
+    TestFile(testdatadir + "more-mmff94.sdf", testdatadir + "more-mmff94e4sresults.txt", "MMFF94", 4.0);
     break;
   default:
     cout << "Test number " << choice << " does not exist!\n";
@@ -199,4 +225,3 @@ int ffmmff94(int argc, char* argv[])
   // Passed tests
   return 0;
 }
-
diff --git a/test/files/mmff94e4results.txt b/test/files/mmff94e4results.txt
@@ -0,0 +1,18 @@
+       81.07932
+        2.91283
+       -3.69683
+        4.24811
+       21.22882
+      118.74484
+      105.47746
+       20.86046
+     1207.63670
+       -0.14958
+      100.53756
+       82.09361
+       64.68310
+       81.07932
+      316.34690
+       35.71481
+       44.77390
+       72.62709
diff --git a/test/files/mmff94sresults.txt b/test/files/mmff94sresults.txt
@@ -0,0 +1,18 @@
+       68.11491
+        8.64766
+      -28.35351
+       -2.84895
+       21.22882
+      118.74484
+       79.03067
+       18.04199
+     1207.63670
+     -113.74578
+      112.23600
+      107.75248
+       93.43627
+       68.11491
+      312.57114
+       19.17321
+       39.99337
+       14.10039
diff --git a/test/files/more-mmff94e4results.txt b/test/files/more-mmff94e4results.txt
@@ -0,0 +1,8 @@
+       69.34483
+       11.26353
+        3.80210
+       -2.68512
+       16.18134
+       35.18203
+       -2.57316
+       62.69012