Added extrapolated polarization support to Python and to manual

docs-source/usersguide/application.rst

@@ -922,29 +922,33 @@ If :code:`vdwCutoff` is not specified, then the value of
 Specifying the Polarization Method
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-OpenMM allows the setting of several other parameters particular to the AMOEBA
-force field.  The :code:`mutualInducedTargetEpsilon` option allows you to
-specify the accuracy to which the induced dipoles are calculated at each time
-step; the default value is 0.01.  The :code:`polarization` setting
-determines whether the calculation of the induced dipoles is continued until the
-dipoles are self-consistent to within the tolerance specified by
-:code:`mutualInducedTargetEpsilon` or whether a quick estimate of the induced
-dipoles is used instead.  The first option corresponds to the
-:code:`polarization='mutual'` setting and is the default; the quick estimate
-option is given by :code:`polarization='direct'` and in this case,
-:code:`mutualInducedTargetEpsilon` is ignored, if provided.  Simulations using
-:code:`polarization='direct'` will be significantly faster than those with
-:code:`polarization='mutual'`\ , but less accurate.  Examples using the two
-options are given below:
-::
-
-    system = forcefield.createSystem(nonbondedMethod=PME,
-        nonbondedCutoff=1*nanometer,ewaldErrorTolerance=0.00001,
-        vdwCutoff=1.2*nanometer, mutualInducedTargetEpsilon=0.01)
-
-    system = forcefield.createSystem(nonbondedMethod=PME,
-        nonbondedCutoff=1*nanometer,ewaldErrorTolerance=0.00001,
-        vdwCutoff=1.2*nanometer, polarization ='direct')
+When using the AMOEBA force field, OpenMM allows the induced dipoles to be
+calculated in any of three different ways.  The slowest but potentially most
+accurate method is to iterate the calculation until the dipoles converge to a
+specified tolerance.  To select this, specify :code:`polarization='mutual'`.
+Use the :code:`mutualInducedTargetEpsilon` option to select the tolerance; for
+most situations, a value of 0.00001 works well.  Alternatively you can specify
+:code:`polarization='extrapolated'`.  This uses an analytic approximation
+:cite:`Simmonett2015` to estimate what the fully converged dipoles will be without
+actually continuing the calculation to convergence.  In many cases this can be
+significantly faster with only a small loss in accuracy.  Finally, you can
+specify :code:`polarization='direct'` to use the direct polarization
+approximation, in which induced dipoles depend only on the fixed multipoles, not
+on other induced dipoles.  This is even faster, but it produces very different
+forces from mutual polarization, so it should only be used with force fields
+that have been specifically parameterized for use with this approximation.
+
+Here are examples of using each method:
+::
+
+    system = forcefield.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
+        vdwCutoff=1.2*nanometer, polarization='mutual', mutualInducedTargetEpsilon=0.00001)
+
+    system = forcefield.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
+        vdwCutoff=1.2*nanometer, polarization='extrapolated')
+
+    system = forcefield.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
+        vdwCutoff=1.2*nanometer, polarization='direct')
 
 
 Implicit Solvent and Solute Dielectrics

docs-source/usersguide/references.bib

@@ -389,6 +389,17 @@
    type = {Journal Article}
 }
 
+@article{Simmonett2015
+   author = {Simmonett, Andrew C. and Pickard, Frank C. and Shao, Yihan and Cheatham, Thomas E. and Brooks, Bernard R.},
+   title = {Efficient treatment of induced dipoles},
+   journal = {Journal of Chemical Physics},
+   year = {2015},
+   volume = {143},
+   number = {7},
+   pages = {074115},
+   type = {Journal Article}
+}
+
 @article{Sindhikara2009,
   author =   {Sindhikara, Daniel J. and Kim, Seonah and Voter,
                   Arthur F. and Roitberg, Adrian E.},

examples/benchmark.py

@@ -36,22 +36,18 @@ def runOneTest(testName, options):
     if amoeba:
         constraints = None
         epsilon = float(options.epsilon)
-        if epsilon == 0:
-            polarization = 'direct'
-        else:
-            polarization = 'mutual'
         if explicit:
             ff = app.ForceField('amoeba2009.xml')
             pdb = app.PDBFile('5dfr_solv-cube_equil.pdb')
             cutoff = 0.7*unit.nanometers
             vdwCutoff = 0.9*unit.nanometers
-            system = ff.createSystem(pdb.topology, nonbondedMethod=app.PME, nonbondedCutoff=cutoff, vdwCutoff=vdwCutoff, constraints=constraints, ewaldErrorTolerance=0.00075, mutualInducedTargetEpsilon=epsilon, polarization=polarization)
+            system = ff.createSystem(pdb.topology, nonbondedMethod=app.PME, nonbondedCutoff=cutoff, vdwCutoff=vdwCutoff, constraints=constraints, ewaldErrorTolerance=0.00075, mutualInducedTargetEpsilon=epsilon, polarization=options.polarization)
         else:
             ff = app.ForceField('amoeba2009.xml', 'amoeba2009_gk.xml')
             pdb = app.PDBFile('5dfr_minimized.pdb')
             cutoff = 2.0*unit.nanometers
             vdwCutoff = 1.2*unit.nanometers
-            system = ff.createSystem(pdb.topology, nonbondedMethod=app.NoCutoff, constraints=constraints, mutualInducedTargetEpsilon=epsilon, polarization=polarization)
+            system = ff.createSystem(pdb.topology, nonbondedMethod=app.NoCutoff, constraints=constraints, mutualInducedTargetEpsilon=epsilon, polarization=options.polarization)
         for f in system.getForces():
             if isinstance(f, mm.AmoebaMultipoleForce) or isinstance(f, mm.AmoebaVdwForce) or isinstance(f, mm.AmoebaGeneralizedKirkwoodForce) or isinstance(f, mm.AmoebaWcaDispersionForce):
                 f.setForceGroup(1)
@@ -127,7 +123,8 @@ parser.add_option('--platform', dest='platform', choices=platformNames, help='na
 parser.add_option('--test', dest='test', choices=('gbsa', 'rf', 'pme', 'amoebagk', 'amoebapme'), help='the test to perform: gbsa, rf, pme, amoebagk, or amoebapme [default: all]')
 parser.add_option('--pme-cutoff', default='0.9', dest='cutoff', type='float', help='direct space cutoff for PME in nm [default: 0.9]')
 parser.add_option('--seconds', default='60', dest='seconds', type='float', help='target simulation length in seconds [default: 60]')
-parser.add_option('--mutual-epsilon', default='1e-4', dest='epsilon', type='float', help='mutual induced epsilon for AMOEBA [default: 1e-4]')
+parser.add_option('--polarization', default='mutual', dest='polarization', choices=('direct', 'extrapolated', 'mutual'), help='the polarization method for AMOEBA: direct, extrapolated, or mutual [default: mutual]')
+parser.add_option('--mutual-epsilon', default='1e-5', dest='epsilon', type='float', help='mutual induced epsilon for AMOEBA [default: 1e-5]')
 parser.add_option('--heavy-hydrogens', action='store_true', default=False, dest='heavy', help='repartition mass to allow a larger time step')
 parser.add_option('--device', default=None, dest='device', help='device index for CUDA or OpenCL')
 parser.add_option('--precision', default='single', dest='precision', choices=('single', 'mixed', 'double'), help='precision mode for CUDA or OpenCL: single, mixed, or double [default: single]')

wrappers/python/simtk/openmm/app/forcefield.py

@@ -3601,6 +3601,8 @@ class AmoebaMultipoleGenerator:
                 polarizationType = args['polarization']
                 if (polarizationType.lower() == 'direct'):
                     force.setPolarizationType(mm.AmoebaMultipoleForce.Direct)
+                elif (polarizationType.lower() == 'extrapolated'):
+                    force.setPolarizationType(mm.AmoebaMultipoleForce.Extrapolated)
                 else:
                     force.setPolarizationType(mm.AmoebaMultipoleForce.Mutual)