Merge pull request #2738 from peastman/drude

Added 2019 CHARMM polarizable force field

docs-source/usersguide/application.rst

@@ -650,8 +650,8 @@ such as :file:`charmm36/water.xml`, which specifies the default CHARMM water mod
 .. warning:: Drude polarizable sites and lone pairs are not yet supported
              by `ParmEd <https://github.com/parmed/parmed>`_ and the CHARMM36 forcefields
              that depend on these features are not included in this port.
-             To use the CHARMM 2013 polarizable force field\ :cite:`Lopes2013`,
-             include the single file :file:`charmm_polar_2013.xml`.
+             To use the CHARMM 2019 polarizable force field\ :cite:`Lopes2013`,
+             include the single file :file:`charmm_polar_2019.xml`.
 
 .. tip:: The solvent model XML files included under the :file:`charmm36/` directory
          include both water *and* ions compatible with that water model, so if you
@@ -712,17 +712,20 @@ recommended for most simulations.
 CHARMM Polarizable Force Field
 ------------------------------
 
-To use the CHARMM 2013 polarizable force field\ :cite:`Lopes2013`, include the
-single file :file:`charmm_polar_2013.xml`.  It includes parameters for proteins,
+To use the CHARMM 2019 polarizable force field\ :cite:`Lopes2013`, include the
+single file :file:`charmm_polar_2019.xml`.  It includes parameters for proteins, lipids,
 water, and ions.  When using this force field, remember to add extra particles to
 the :class:`Topology` as described in section :ref:`adding-or-removing-extra-particles`.
+This force field also requires that you use one of the special integrators that
+supports Drude particles.  The options are DrudeLangevinIntegrator, DrudeNoseHooverIntegrator,
+and DrudeSCFIntegrator.
 
-Older Amber Force Fields
-------------------------
+Older Force Fields
+------------------
 
-OpenMM includes several older Amber force fields as well.  For most simulations
-Amber14 is preferred over any of these, but they are still useful for reproducing
-older results.
+OpenMM includes several older force fields as well.  For most simulations, the
+newer force fields described above are preferred over any of these, but they are
+still useful for reproducing older results.
 
 .. tabularcolumns:: |l|L|
 
@@ -735,6 +738,7 @@ File                           Force Field
 :code:`amber99sbnmr.xml`       Amber99SB with modifications to fit NMR data\ :cite:`Li2010`
 :code:`amber03.xml`            Amber03\ :cite:`Duan2003`
 :code:`amber10.xml`            Amber10 (documented in the AmberTools_ manual as `ff10`)
+:code:`charmm_polar_2013.xml`  2013 version of the CHARMM polarizable force field\ :cite:`Lopes2013`
 =============================  ================================================================================
 
 Several of these force fields support implicit solvent.  To enable it, also

docs-source/usersguide/library.rst

@@ -3644,7 +3644,7 @@ The equations of motion can be integrated with three different methods:
 #. In the extended Lagrangian method, the positions of the Drude particles are
    treated as dynamical variables, just like any other particles.  A small amount
    of mass is transferred from the parent particles to the Drude particles,
-   allowing them to be integrated normally.  A dual Langevin integrator is used to
+   allowing them to be integrated normally.  A dual Langevin or Nose-Hoover integrator is used to
    maintain the center of mass of each Drude particle pair at the system
    temperature, while using a much lower temperature for their relative internal
    motion.  In practice, this produces dipole moments very close to those from the

wrappers/python/tests/TestForceField.py

@@ -1019,6 +1019,75 @@ END"""))
         # If the check is not done correctly, this will throw an exception.
         ff.createSystem(pdb.topology)
 
+    def test_CharmmPolar(self):
+        """Test the CHARMM polarizable force field."""
+        pdb = PDBFile('systems/ala_ala_ala_drude.pdb')
+        pdb.topology.setUnitCellDimensions(Vec3(3, 3, 3))
+        ff = ForceField('charmm_polar_2019.xml')
+        system = ff.createSystem(pdb.topology, nonbondedMethod=PME, nonbondedCutoff=1.2*nanometers)
+        for i,f in enumerate(system.getForces()):
+            f.setForceGroup(i)
+            if isinstance(f, NonbondedForce):
+                f.setPMEParameters(3.4, 64, 64, 64)
+        integrator = DrudeLangevinIntegrator(300, 1.0, 1.0, 10.0, 0.001)
+        context = Context(system, integrator, Platform.getPlatformByName('Reference'))
+        context.setPositions(pdb.positions)
+
+        # Compare the energy to values computed by CHARMM.  Here is what it outputs:
+
+        # ENER ENR:  Eval#     ENERgy      Delta-E         GRMS
+        # ENER INTERN:          BONDs       ANGLes       UREY-b    DIHEdrals    IMPRopers
+        # ENER CROSS:           CMAPs        PMF1D        PMF2D        PRIMO
+        # ENER EXTERN:        VDWaals         ELEC       HBONds          ASP         USER
+        # ENER EWALD:          EWKSum       EWSElf       EWEXcl       EWQCor       EWUTil
+        #  ----------       ---------    ---------    ---------    ---------    ---------
+        # ENER>        0    102.83992      0.00000     13.06415
+        # ENER INTERN>       54.72574     40.21459     11.61009     26.10373      0.14113
+        # ENER CROSS>        -3.37113      0.00000      0.00000      0.00000
+        # ENER EXTERN>       22.74761    -24.21667      0.00000      0.00000      0.00000
+        # ENER EWALD>        56.14258  -7279.07968   7197.82192      0.00000      0.00000
+        #  ----------       ---------    ---------    ---------    ---------    ---------
+
+        # First check the total energy.
+        
+        energy = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+        self.assertAlmostEqual(102.83992, energy, delta=energy*1e-3)
+
+        # Now check individual components.  CHARMM and OpenMM split them up a little differently.  I've tried to
+        # match things up, but I think there's still some inconsistency in where forces related to Drude particles
+        # are categorized.  That's why the Coulomb and bonds terms match less accurately than the other terms
+        # (and less accurately than the total energy, which agrees well).
+
+        coulomb = 0
+        vdw = 0
+        bonds = 0
+        angles = 0
+        propers = 0
+        impropers = 0
+        cmap = 0
+        for i,f in enumerate(system.getForces()):
+            energy = context.getState(getEnergy=True, groups={i}).getPotentialEnergy().value_in_unit(kilocalories_per_mole)
+            if isinstance(f, NonbondedForce):
+                coulomb += energy
+            elif isinstance(f, CustomNonbondedForce) or isinstance(f, CustomBondForce):
+                vdw += energy
+            elif isinstance(f, HarmonicBondForce) or isinstance(f, DrudeForce):
+                bonds += energy
+            elif isinstance(f, HarmonicAngleForce):
+                angles += energy
+            elif isinstance(f, PeriodicTorsionForce):
+                propers += energy
+            elif isinstance(f, CustomTorsionForce):
+                impropers += energy
+            elif isinstance(f, CMAPTorsionForce):
+                cmap += energy
+        self.assertAlmostEqual(-24.21667+56.14258-7279.07968+7197.82192, coulomb, delta=abs(coulomb)*5e-2) # ELEC+EWKSum+EWSElf+EWEXcl
+        self.assertAlmostEqual(22.74761, vdw, delta=vdw*1e-3) # VDWaals
+        self.assertAlmostEqual(54.72574+11.61009, bonds, delta=bonds*2e-2) # BONDs+UREY-b
+        self.assertAlmostEqual(40.21459, angles, delta=angles*1e-3) # ANGLes
+        self.assertAlmostEqual(26.10373, propers, delta=propers*1e-3) # DIHEdrals
+        self.assertAlmostEqual(0.14113, impropers, delta=impropers*1e-3) # IMPRopers
+
 class AmoebaTestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method with the AMOEBA forcefield."""
 

wrappers/python/tests/systems/ala_ala_ala_drude.pdb

+REMARK   1 CREATED WITH OPENMM 7.5, 2020-06-04
+ATOM      1  N   ALA A   1       0.072  -0.034  -0.113  1.00  0.00           N  
+ATOM      2 DN   ALA A   1       0.024  -0.103  -0.101  1.00  0.00          EP  
+ATOM      3  H   ALA A   1      -0.396   0.850  -0.349  1.00  0.00           H  
+ATOM      4  H2  ALA A   1      -0.547  -0.623   0.444  1.00  0.00           H  
+ATOM      5  H3  ALA A   1       0.343  -0.509  -0.974  1.00  0.00           H  
+ATOM      6  CA  ALA A   1       1.286   0.359   0.692  1.00  0.00           C  
+ATOM      7 DCA  ALA A   1       1.247   0.375   0.636  1.00  0.00          EP  
+ATOM      8  HA  ALA A   1       0.730   0.825   1.494  1.00  0.00           H  
+ATOM      9  C   ALA A   1       1.986   1.606   0.031  1.00  0.00           C  
+ATOM     10 DC   ALA A   1       1.956   1.579   0.036  1.00  0.00          EP  
+ATOM     11  O   ALA A   1       1.230   2.563  -0.208  1.00  0.00           O  
+ATOM     12 DO   ALA A   1       1.219   2.525  -0.201  1.00  0.00          EP  
+ATOM     13 LPOA ALA A   1       0.996   2.412  -0.097  1.00  0.00          EP  
+ATOM     14 LPOB ALA A   1       1.458   2.722  -0.321  1.00  0.00          EP  
+ATOM     15  CB  ALA A   1       2.135  -0.773   1.297  1.00  0.00           C  
+ATOM     16 DCB  ALA A   1       2.057  -0.772   1.289  1.00  0.00          EP  
+ATOM     17  HB1 ALA A   1       2.779  -1.237   0.521  1.00  0.00           H  
+ATOM     18  HB2 ALA A   1       1.483  -1.559   1.736  1.00  0.00           H  
+ATOM     19  HB3 ALA A   1       2.788  -0.373   2.102  1.00  0.00           H  
+ATOM     20  N   ALA A   2       3.340   1.626  -0.172  1.00  0.00           N  
+ATOM     21 DN   ALA A   2       3.289   1.631  -0.202  1.00  0.00          EP  
+ATOM     22  H   ALA A   2       3.913   0.773  -0.267  1.00  0.00           H  
+ATOM     23  CA  ALA A   2       3.977   2.932  -0.230  1.00  0.00           C  
+ATOM     24 DCA  ALA A   2       3.990   2.909  -0.215  1.00  0.00          EP  
+ATOM     25  HA  ALA A   2       3.675   3.404   0.715  1.00  0.00           H  
+ATOM     26  C   ALA A   2       5.518   2.682  -0.121  1.00  0.00           C  
+ATOM     27 DC   ALA A   2       5.487   2.654  -0.128  1.00  0.00          EP  
+ATOM     28  O   ALA A   2       5.947   1.505  -0.133  1.00  0.00           O  
+ATOM     29 DO   ALA A   2       5.889   1.489  -0.137  1.00  0.00          EP  
+ATOM     30 LPOA ALA A   2       5.668   1.399  -0.158  1.00  0.00          EP  
+ATOM     31 LPOB ALA A   2       6.229   1.603  -0.108  1.00  0.00          EP  
+ATOM     32  CB  ALA A   2       3.712   3.828  -1.442  1.00  0.00           C  
+ATOM     33 DCB  ALA A   2       3.662   3.802  -1.434  1.00  0.00          EP  
+ATOM     34  HB1 ALA A   2       2.636   4.055  -1.572  1.00  0.00           H  
+ATOM     35  HB2 ALA A   2       4.080   3.339  -2.370  1.00  0.00           H  
+ATOM     36  HB3 ALA A   2       4.238   4.801  -1.334  1.00  0.00           H  
+ATOM     37  N   ALA A   3       6.281   3.795  -0.031  1.00  0.00           N  
+ATOM     38 DN   ALA A   3       6.275   3.733  -0.037  1.00  0.00          EP  
+ATOM     39  H   ALA A   3       5.825   4.723   0.051  1.00  0.00           H  
+ATOM     40  CA  ALA A   3       7.692   3.798   0.060  1.00  0.00           C  
+ATOM     41 DCA  ALA A   3       7.707   3.802   0.068  1.00  0.00          EP  
+ATOM     42  HA  ALA A   3       8.015   3.302  -0.864  1.00  0.00           H  
+ATOM     43  C   ALA A   3       8.044   5.364   0.125  1.00  0.00           C  
+ATOM     44 DC   ALA A   3       8.018   5.323   0.136  1.00  0.00          EP  
+ATOM     45  CB  ALA A   3       8.274   3.163   1.304  1.00  0.00           C  
+ATOM     46 DCB  ALA A   3       8.233   3.093   1.333  1.00  0.00          EP  
+ATOM     47  HB1 ALA A   3       8.009   2.087   1.363  1.00  0.00           H  
+ATOM     48  HB2 ALA A   3       7.899   3.666   2.221  1.00  0.00           H  
+ATOM     49  HB3 ALA A   3       9.384   3.230   1.300  1.00  0.00           H  
+ATOM     50  O   ALA A   3       7.074   6.193   0.135  1.00  0.00           O  
+ATOM     51 DOT1 ALA A   3       7.032   6.119   0.127  1.00  0.00          EP  
+ATOM     52  OXT ALA A   3       9.258   5.664   0.187  1.00  0.00           O  
+ATOM     53 DOT2 ALA A   3       9.219   5.692   0.188  1.00  0.00          EP  
+ATOM     54 LPT1 ALA A   3       7.196   6.521   0.158  1.00  0.00          EP  
+ATOM     55 LPT2 ALA A   3       6.769   6.022   0.114  1.00  0.00          EP  
+ATOM     56 LPT3 ALA A   3       9.295   6.011   0.208  1.00  0.00          EP  
+ATOM     57 LPT4 ALA A   3       9.454   5.374   0.178  1.00  0.00          EP  
+TER      58      ALA A   3
+END