Skip to content

GitLab

Commit 5c83e423 authored by peastman's avatar peastman
Browse files

Added 2019 CHARMM polarizable force field

parent e4c1f73c
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 19922 additions and 10 deletions
docs-source/usersguide/application.rst
View file @ 5c83e423
......@@ -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
View file @ 5c83e423
......@@ -3644,7 +3644,7 @@ The equations of motion can be integrated with two 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/simtk/openmm/app/data/charmm_polar_2019.xml 0 → 100644
View file @ 5c83e423
This diff is collapsed.
wrappers/python/tests/TestForceField.py
View file @ 5c83e423
......@@ -1019,6 +1019,72 @@ 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.46525 0.00000 13.22189
# 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 -23.94811 0.00000 0.00000 0.00000
# ENER EWALD> 55.93962 -7324.67582 7242.97779 0.00000 0.00000
# ---------- --------- --------- --------- --------- ---------
# First check the total energy.
energy = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilocalories_per_mole)
self.assertAlmostEqual(102.46525, energy, delta=energy*5e-3)
# Now check individual components. CHARMM and OpenMM split them up a little differently.
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(-23.94811+55.93962-7324.67582+7242.97779, coulomb, delta=abs(coulomb)*5e-2) # ELEC+EWKSum+EWSElf+EWEXcl
self.assertAlmostEqual(22.74761, vdw, delta=vdw*5e-3) # VDWaals
self.assertAlmostEqual(54.72574+11.61009, bonds, delta=bonds*2e-2) # BONDs+UREY-b
self.assertAlmostEqual(40.21459, angles, delta=angles*5e-3) # ANGLes
self.assertAlmostEqual(26.10373, propers, delta=propers*5e-3) # DIHEdrals
self.assertAlmostEqual(0.14113, impropers, delta=impropers*5e-3) # IMPRopers
class AmoebaTestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method with the AMOEBA forcefield."""
......
wrappers/python/tests/systems/ala_ala_ala_drude.pdb 0 → 100644
View file @ 5c83e423
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
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment