Merge pull request #531 from peastman/charmm

Implemented CHARMM polarizable force field

docs/usersguide/application.rst

@@ -651,19 +651,20 @@ For the main force field, OpenMM provides the following options:
 
 .. tabularcolumns:: |l|L|
 
-=================  ================================================================================
-File               Force Field                                                                     
-=================  ================================================================================
-amber96.xml        AMBER96\ :cite:`Kollman1997`
-amber99sb.xml      AMBER99\ :cite:`Wang2000` with modified backbone torsions\ :cite:`Hornak2006`
-amber99sbildn.xml  AMBER99SB plus improved side chain torsions\ :cite:`Lindorff-Larsen2010`
-amber99sbnmr.xml   AMBER99SB with modifications to fit NMR data\ :cite:`Li2010`
-amber03.xml        AMBER03\ :cite:`Duan2003`
-amber10.xml        AMBER10
-amoeba2009.xml     AMOEBA 2009\ :cite:`Ren2002`.  This force field is deprecated.  It is 
-                   recommended to use AMOEBA 2013 instead.
-amoeba2013.xml     AMOEBA 2013\ :cite:`Shi2013`
-=================  ================================================================================
+=====================  ================================================================================
+File                   Force Field                                                                     
+=====================  ================================================================================
+amber96.xml            AMBER96\ :cite:`Kollman1997`
+amber99sb.xml          AMBER99\ :cite:`Wang2000` with modified backbone torsions\ :cite:`Hornak2006`
+amber99sbildn.xml      AMBER99SB plus improved side chain torsions\ :cite:`Lindorff-Larsen2010`
+amber99sbnmr.xml       AMBER99SB with modifications to fit NMR data\ :cite:`Li2010`
+amber03.xml            AMBER03\ :cite:`Duan2003`
+amber10.xml            AMBER10
+amoeba2009.xml         AMOEBA 2009\ :cite:`Ren2002`.  This force field is deprecated.  It is 
+                       recommended to use AMOEBA 2013 instead.
+amoeba2013.xml         AMOEBA 2013\ :cite:`Shi2013`
+charmm_polar_2013.xml  CHARMM 2013 polarizable force field\ :cite:`Lopes2013`
+=====================  ================================================================================
 
 
 The AMBER files do not include parameters for water molecules.  This allows you
@@ -686,10 +687,10 @@ swm4ndp.xml  SWM4-NDP water model\ :cite:`Lamoureux2006`
 ===========  ============================================
 
 
-For the AMOEBA force field, only one explicit water model is currently available
-and the water parameters are included in the file :code:`amoeba2009.xml`\ .
-Also the AMOEBA force field file only includes the parameters for amino acids
-and ions; nucleic acids will be included in a future release.
+For the polarizable force fields (AMOEBA and CHARMM), only one explicit water model
+is currently available and the water parameters are included in the same file as
+the macromolecule parameters.  Also, the polarizable force fields only include
+parameters for amino acids and ions, not for nucleic acids.
 
 If you want to include an implicit solvation model, you can also specify one of
 the following files:

docs/usersguide/references.bib

@@ -215,6 +215,17 @@
    type = {Journal Article}
 }
 
+@article{Lopes2013,
+    author = {Lopes, Pedro E. M. and Huang, Jing and Shim, Jihyun and Luo, Yun and Li, Hui and Roux, Benoît and MacKerell, Alexander D.},
+    title = {Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator},
+    journal = {Journal of Chemical Theory and Computation},
+    volume = {9},
+    number = {12},
+    pages = {5430-5449},
+    year = {2013},
+    type = {Journal Article}
+}
+
 @article{Mahoney2000
    author = {Mahoney, Michael W. and Jorgensen, William L.},
    title = {A five-site model for liquid water and the reproduction of the density anomaly by rigid, nonpolarizable potential functions},

plugins/drude/platforms/cuda/src/CudaDrudeKernels.cpp

@@ -126,9 +126,9 @@ void CudaCalcDrudeForceKernel::initialize(const System& system, const DrudeForce
             double a1 = (atoms[i][2] == -1 ? 1 : aniso12);
             double a2 = (atoms[i][3] == -1 || atoms[i][4] == -1 ? 1 : aniso34);
             double a3 = 3-a1-a2;
-            double k3 = charge*charge/(polarizability*a3);
-            double k1 = charge*charge/(polarizability*a1) - k3;
-            double k2 = charge*charge/(polarizability*a2) - k3;
+            double k3 = ONE_4PI_EPS0*charge*charge/(polarizability*a3);
+            double k1 = ONE_4PI_EPS0*charge*charge/(polarizability*a1) - k3;
+            double k2 = ONE_4PI_EPS0*charge*charge/(polarizability*a2) - k3;
             if (atoms[i][2] == -1) {
                 atoms[i][2] = 0;
                 k1 = 0;

plugins/drude/platforms/cuda/src/kernels/drudePairForce.cu

@@ -13,7 +13,7 @@ real u = drudeParams.x*r;
 real screening = 1-(1+0.5f*u)*EXP(-u);
 real pairEnergy = drudeParams.y*screening*rInv;
 energy += pairEnergy;
-real3 f = delta*(pairEnergy*rInv*rInv);
+real3 f = delta*(drudeParams.y*rInv*rInv*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x));
 force1 += f;
 force3 -= f;
 
@@ -26,7 +26,7 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = -drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(-drudeParams.y*rInv*rInv*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x));
 force1 += f;
 force4 -= f;
 
@@ -39,7 +39,7 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = -drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(-drudeParams.y*rInv*rInv*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x));
 force2 += f;
 force3 -= f;
 
@@ -52,6 +52,6 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(drudeParams.y*rInv*rInv*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x));
 force2 += f;
 force4 -= f;

plugins/drude/platforms/cuda/tests/TestCudaDrudeForce.cpp

@@ -76,9 +76,9 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg
 }
 
 void testSingleParticle() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     System system;
     system.addParticle(1.0);
     system.addParticle(1.0);
@@ -92,9 +92,9 @@ void testSingleParticle() {
 }
 
 void testAnisotropicParticle() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double a1 = 0.8;
     const double a2 = 1.1;
     const double k1 = k/a1;
@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) {
 }
 
 void testThole() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double thole = 2.5;
     System system;
     system.addParticle(1.0);
@@ -157,9 +157,9 @@ void testThole() {
 }
 
 void testChangingParameters() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     Platform& platform = Platform::getPlatformByName("CUDA");
     
     // Create the system.
@@ -184,9 +184,9 @@ void testChangingParameters() {
     
     // Modify the parameters.
     
-    const double k2 = 2.2;
+    const double k2 = ONE_4PI_EPS0*2.2;
     const double charge2 = 0.3;
-    const double alpha2 = charge2*charge2/k2;
+    const double alpha2 = ONE_4PI_EPS0*charge2*charge2/k2;
     drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
     drude->updateParametersInContext(context);
     state = context.getState(State::Energy);

plugins/drude/platforms/cuda/tests/TestCudaDrudeLangevinIntegrator.cpp

@@ -53,9 +53,9 @@ extern "C" OPENMM_EXPORT void registerDrudeCudaKernelFactories();
 void testSinglePair() {
     const double temperature = 300.0;
     const double temperatureDrude = 10.0;
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double mass1 = 1.0;
     const double mass2 = 0.1;
     const double totalMass = mass1+mass2;
@@ -132,7 +132,7 @@ void testWater() {
         system.addConstraint(startIndex, startIndex+3, 0.09572);
         system.addConstraint(startIndex+2, startIndex+3, 0.15139);
         system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
-        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, ONE_4PI_EPS0*1.71636*1.71636/(100000*4.184), 1, 1);
     }
     vector<Vec3> positions;
     for (int i = 0; i < gridSize; i++)

plugins/drude/platforms/cuda/tests/TestCudaDrudeSCFIntegrator.cpp

@@ -85,7 +85,7 @@ void testWater() {
         system.addConstraint(startIndex, startIndex+3, 0.09572);
         system.addConstraint(startIndex+2, startIndex+3, 0.15139);
         system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
-        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, ONE_4PI_EPS0*1.71636*1.71636/(100000*4.184), 1, 1);
     }
     vector<Vec3> positions;
     for (int i = 0; i < gridSize; i++)

plugins/drude/platforms/opencl/src/OpenCLDrudeKernels.cpp

@@ -131,9 +131,9 @@ void OpenCLCalcDrudeForceKernel::initialize(const System& system, const DrudeFor
             double a1 = (atoms[i][2] == -1 ? 1 : aniso12);
             double a2 = (atoms[i][3] == -1 || atoms[i][4] == -1 ? 1 : aniso34);
             double a3 = 3-a1-a2;
-            double k3 = charge*charge/(polarizability*a3);
-            double k1 = charge*charge/(polarizability*a1) - k3;
-            double k2 = charge*charge/(polarizability*a2) - k3;
+            double k3 = ONE_4PI_EPS0*charge*charge/(polarizability*a3);
+            double k1 = ONE_4PI_EPS0*charge*charge/(polarizability*a1) - k3;
+            double k2 = ONE_4PI_EPS0*charge*charge/(polarizability*a2) - k3;
             if (atoms[i][2] == -1) {
                 atoms[i][2] = 0;
                 k1 = 0;

plugins/drude/platforms/opencl/src/kernels/drudePairForce.cl

@@ -13,7 +13,7 @@ real u = drudeParams.x*r;
 real screening = 1-(1+0.5f*u)*EXP(-u);
 real pairEnergy = drudeParams.y*screening*rInv;
 energy += pairEnergy;
-real4 f = delta*(pairEnergy*rInv*rInv);
+real4 f = delta*(drudeParams.y*rInv*rInv)*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x);
 force1 += f;
 force3 -= f;
 
@@ -26,7 +26,7 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = -drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(-drudeParams.y*rInv*rInv)*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x);
 force1 += f;
 force4 -= f;
 
@@ -39,7 +39,7 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = -drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(-drudeParams.y*rInv*rInv)*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x);
 force2 += f;
 force3 -= f;
 
@@ -52,6 +52,6 @@ u = drudeParams.x*r;
 screening = 1-(1+0.5f*u)*EXP(-u);
 pairEnergy = drudeParams.y*screening*rInv;
 energy += pairEnergy;
-f = delta*(pairEnergy*rInv*rInv);
+f = delta*(drudeParams.y*rInv*rInv)*(screening*rInv-0.5f*(1+u)*EXP(-u)*drudeParams.x);
 force2 += f;
 force4 -= f;

plugins/drude/platforms/opencl/tests/TestOpenCLDrudeForce.cpp

@@ -76,9 +76,9 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg
 }
 
 void testSingleParticle() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     System system;
     system.addParticle(1.0);
     system.addParticle(1.0);
@@ -92,9 +92,9 @@ void testSingleParticle() {
 }
 
 void testAnisotropicParticle() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double a1 = 0.8;
     const double a2 = 1.1;
     const double k1 = k/a1;
@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) {
 }
 
 void testThole() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double thole = 2.5;
     System system;
     system.addParticle(1.0);
@@ -157,9 +157,9 @@ void testThole() {
 }
 
 void testChangingParameters() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     Platform& platform = Platform::getPlatformByName("OpenCL");
     
     // Create the system.
@@ -184,9 +184,9 @@ void testChangingParameters() {
     
     // Modify the parameters.
     
-    const double k2 = 2.2;
+    const double k2 = ONE_4PI_EPS0*2.2;
     const double charge2 = 0.3;
-    const double alpha2 = charge2*charge2/k2;
+    const double alpha2 = ONE_4PI_EPS0*charge2*charge2/k2;
     drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
     drude->updateParametersInContext(context);
     state = context.getState(State::Energy);

plugins/drude/platforms/opencl/tests/TestOpenCLDrudeLangevinIntegrator.cpp

@@ -53,9 +53,9 @@ extern "C" OPENMM_EXPORT void registerDrudeOpenCLKernelFactories();
 void testSinglePair() {
     const double temperature = 300.0;
     const double temperatureDrude = 10.0;
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double mass1 = 1.0;
     const double mass2 = 0.1;
     const double totalMass = mass1+mass2;
@@ -132,7 +132,7 @@ void testWater() {
         system.addConstraint(startIndex, startIndex+3, 0.09572);
         system.addConstraint(startIndex+2, startIndex+3, 0.15139);
         system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
-        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, ONE_4PI_EPS0*1.71636*1.71636/(100000*4.184), 1, 1);
     }
     vector<Vec3> positions;
     for (int i = 0; i < gridSize; i++)

plugins/drude/platforms/opencl/tests/TestOpenCLDrudeSCFIntegrator.cpp

@@ -85,7 +85,7 @@ void testWater() {
         system.addConstraint(startIndex, startIndex+3, 0.09572);
         system.addConstraint(startIndex+2, startIndex+3, 0.15139);
         system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
-        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, ONE_4PI_EPS0*1.71636*1.71636/(100000*4.184), 1, 1);
     }
     vector<Vec3> positions;
     for (int i = 0; i < gridSize; i++)

plugins/drude/platforms/reference/src/ReferenceDrudeKernels.cpp

@@ -6,7 +6,7 @@
  * Biological Structures at Stanford, funded under the NIH Roadmap for        *
  * Medical Research, grant U54 GM072970. See https://simtk.org.               *
  *                                                                            *
- * Portions copyright (c) 2011-2013 Stanford University and the Authors.      *
+ * Portions copyright (c) 2011-2014 Stanford University and the Authors.      *
  * Authors: Peter Eastman                                                     *
  * Contributors:                                                              *
  *                                                                            *
@@ -136,9 +136,9 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
         RealOpenMM a1 = (p2 == -1 ? 1 : aniso12[i]);
         RealOpenMM a2 = (p3 == -1 || p4 == -1 ? 1 : aniso34[i]);
         RealOpenMM a3 = 3-a1-a2;
-        RealOpenMM k3 = charge[i]*charge[i]/(polarizability[i]*a3);
-        RealOpenMM k1 = charge[i]*charge[i]/(polarizability[i]*a1) - k3;
-        RealOpenMM k2 = charge[i]*charge[i]/(polarizability[i]*a2) - k3;
+        RealOpenMM k3 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a3);
+        RealOpenMM k1 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a1) - k3;
+        RealOpenMM k2 = ONE_4PI_EPS0*charge[i]*charge[i]/(polarizability[i]*a2) - k3;
         
         // Compute the isotropic force.
         
@@ -188,6 +188,7 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
         int dipole2 = pair2[i];
         int dipole1Particles[] = {particle[dipole1], particle1[dipole1]};
         int dipole2Particles[] = {particle[dipole2], particle1[dipole2]};
+        RealOpenMM uscale = pairThole[i]/pow(polarizability[dipole1]*polarizability[dipole2], 1.0/6.0);
         for (int j = 0; j < 2; j++)
             for (int k = 0; k < 2; k++) {
                 int p1 = dipole1Particles[j];
@@ -195,10 +196,10 @@ double ReferenceCalcDrudeForceKernel::execute(ContextImpl& context, bool include
                 RealOpenMM chargeProduct = charge[dipole1]*charge[dipole2]*(j == k ? 1 : -1);
                 RealVec delta = pos[p1]-pos[p2];
                 RealOpenMM r = sqrt(delta.dot(delta));
-                RealOpenMM u = r*pairThole[i]/pow(polarizability[dipole1]*polarizability[dipole2], 1.0/6.0);
+                RealOpenMM u = r*uscale;
                 RealOpenMM screening = 1.0 - (1.0+0.5*u)*exp(-u);
                 energy += ONE_4PI_EPS0*chargeProduct*screening/r;
-                RealVec f = delta*(ONE_4PI_EPS0*chargeProduct*screening/(r*r*r));
+                RealVec f = delta*(ONE_4PI_EPS0*chargeProduct/(r*r))*(screening/r-0.5*(1+u)*exp(-u)*uscale);
                 force[p1] += f;
                 force[p2] -= f;
             }
@@ -461,4 +462,4 @@ void ReferenceIntegrateDrudeSCFStepKernel::minimize(ContextImpl& context, double
     lbfgsfloatval_t fx;
     MinimizerData data(context, drudeParticles);
     lbfgs(numDrudeParticles*3, minimizerPos, &fx, evaluate, NULL, &data, &minimizerParams);
-}
+}
\ No newline at end of file

plugins/drude/platforms/reference/tests/TestReferenceDrudeForce.cpp

@@ -71,14 +71,14 @@ void validateForce(System& system, vector<Vec3>& positions, double expectedEnerg
             offsetPos[i][j] = positions[i][j]+offset;
             context.setPositions(offsetPos);
             double e2 = context.getState(State::Energy | State::Forces).getPotentialEnergy();
-            ASSERT_EQUAL_TOL(state.getForces()[i][j], (e1-e2)/(2*offset), 1e-3);
+            ASSERT_EQUAL_TOL(state.getForces()[i][j], (e1-e2)/(2*offset), 1e-5);
         }
 }
 
 void testSingleParticle() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     System system;
     system.addParticle(1.0);
     system.addParticle(1.0);
@@ -92,9 +92,9 @@ void testSingleParticle() {
 }
 
 void testAnisotropicParticle() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double a1 = 0.8;
     const double a2 = 1.1;
     const double k1 = k/a1;
@@ -124,9 +124,9 @@ double computeScreening(double r, double thole, double alpha1, double alpha2) {
 }
 
 void testThole() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double thole = 2.5;
     System system;
     system.addParticle(1.0);
@@ -157,9 +157,9 @@ void testThole() {
 }
 
 void testChangingParameters() {
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     Platform& platform = Platform::getPlatformByName("Reference");
     
     // Create the system.
@@ -184,9 +184,9 @@ void testChangingParameters() {
     
     // Modify the parameters.
     
-    const double k2 = 2.2;
+    const double k2 = ONE_4PI_EPS0*2.2;
     const double charge2 = 0.3;
-    const double alpha2 = charge2*charge2/k2;
+    const double alpha2 = ONE_4PI_EPS0*charge2*charge2/k2;
     drude->setParticleParameters(0, 1, 0, -1, -1, -1, charge2, alpha2, 1, 1);
     drude->updateParametersInContext(context);
     state = context.getState(State::Energy);

plugins/drude/platforms/reference/tests/TestReferenceDrudeLangevinIntegrator.cpp

@@ -53,9 +53,9 @@ extern "C" OPENMM_EXPORT void registerDrudeReferenceKernelFactories();
 void testSinglePair() {
     const double temperature = 300.0;
     const double temperatureDrude = 10.0;
-    const double k = 1.5;
+    const double k = ONE_4PI_EPS0*1.5;
     const double charge = 0.1;
-    const double alpha = charge*charge/k;
+    const double alpha = ONE_4PI_EPS0*charge*charge/k;
     const double mass1 = 1.0;
     const double mass2 = 0.1;
     const double totalMass = mass1+mass2;
@@ -132,7 +132,7 @@ void testWater() {
         system.addConstraint(startIndex, startIndex+3, 0.09572);
         system.addConstraint(startIndex+2, startIndex+3, 0.15139);
         system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
-        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, ONE_4PI_EPS0*1.71636*1.71636/(100000*4.184), 1, 1);
     }
     vector<Vec3> positions;
     for (int i = 0; i < gridSize; i++)

plugins/drude/platforms/reference/tests/TestReferenceDrudeSCFIntegrator.cpp

@@ -85,7 +85,7 @@ void testWater() {
         system.addConstraint(startIndex, startIndex+3, 0.09572);
         system.addConstraint(startIndex+2, startIndex+3, 0.15139);
         system.setVirtualSite(startIndex+4, new ThreeParticleAverageSite(startIndex, startIndex+2, startIndex+3, 0.786646558, 0.106676721, 0.106676721));
-        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, 1.71636*1.71636/(100000*4.184), 1, 1);
+        drude->addParticle(startIndex+1, startIndex, -1, -1, -1, -1.71636, ONE_4PI_EPS0*1.71636*1.71636/(100000*4.184), 1, 1);
     }
     vector<Vec3> positions;
     for (int i = 0; i < gridSize; i++)

wrappers/python/simtk/openmm/app/data/swm4ndp.xml

@@ -30,7 +30,7 @@
   <Atom type="swm4ndp-OD" charge="-1.71636" sigma="1" epsilon="0"/>
  </NonbondedForce>
  <DrudeForce>
-  <Particle type1="swm4ndp-OD" type2="swm4ndp-O" charge="-1.71636" polarizability="7.040850e-6" thole="1.3"/>
+  <Particle type1="swm4ndp-OD" type2="swm4ndp-O" charge="-1.71636" polarizability="0.000978253" thole="1.3"/>
  </DrudeForce>
 </ForceField>
 

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

@@ -240,6 +240,7 @@ class ForceField(object):
         def __init__(self):
             self.atomType = {}
             self.atoms = []
+            self.excludeAtomWith = []
             self.virtualSites = {}
             self.bonds = []
             self.angles = []
@@ -302,6 +303,10 @@ class ForceField(object):
                 self.localPos = [float(attrib['p1']), float(attrib['p2']), float(attrib['p3'])]
             else:
                 raise ValueError('Unknown virtual site type: %s' % self.type)
+            if 'excludeWith' in attrib:
+                self.excludeWith = int(attrib['excludeWith'])
+            else:
+                self.excludeWith = self.atoms[0]
 
     def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
                      constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, **args):
@@ -324,6 +329,8 @@ class ForceField(object):
         """
         data = ForceField._SystemData()
         data.atoms = list(topology.atoms())
+        for atom in data.atoms:
+            data.excludeAtomWith.append([])
 
         # Make a list of all bonds
 
@@ -363,7 +370,7 @@ class ForceField(object):
                     data.atomType[atom] = template.atoms[match].type
                     for site in template.virtualSites:
                         if match == site.index:
-                            data.virtualSites[atom] = (site, [matchAtoms[i].index for i in site.atoms])
+                            data.virtualSites[atom] = (site, [matchAtoms[i].index for i in site.atoms], matchAtoms[site.excludeWith].index)
 
         # Create the System and add atoms
 
@@ -478,8 +485,9 @@ class ForceField(object):
         # Add virtual sites
 
         for atom in data.virtualSites:
-            (site, atoms) = data.virtualSites[atom]
+            (site, atoms, excludeWith) = data.virtualSites[atom]
             index = atom.index
+            data.excludeAtomWith[excludeWith].append(index)
             if site.type == 'average2':
                 sys.setVirtualSite(index, mm.TwoParticleAverageSite(atoms[0], atoms[1], site.weights[0], site.weights[1]))
             elif site.type == 'average3':
@@ -895,7 +903,7 @@ class PeriodicTorsionGenerator:
                     for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                         if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
                             # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
-                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
+                            # impropers, which leaves the ordering ambiguous.  It then follows some bizarre rules
                             # to pick the order.
                             a1 = torsion[t2[1]]
                             a2 = torsion[t3[1]]
@@ -991,18 +999,22 @@ class RBTorsionGenerator:
                 if type1 in types1:
                     for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                         if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
-                            # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
-                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
-                            # to pick the order.
-                            a1 = torsion[t2[1]]
-                            a2 = torsion[t3[1]]
-                            e1 = data.atoms[a1].element
-                            e2 = data.atoms[a2].element
-                            if e1 == e2 and a1 > a2:
-                                (a1, a2) = (a2, a1)
-                            elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
-                                (a1, a2) = (a2, a1)
-                            force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
+                            if wildcard in (types1, types2, types3, types4):
+                                # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
+                                # impropers, which leaves the ordering ambiguous.  It then follows some bizarre rules
+                                # to pick the order.
+                                a1 = torsion[t2[1]]
+                                a2 = torsion[t3[1]]
+                                e1 = data.atoms[a1].element
+                                e2 = data.atoms[a2].element
+                                if e1 == e2 and a1 > a2:
+                                    (a1, a2) = (a2, a1)
+                                elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
+                                    (a1, a2) = (a2, a1)
+                                force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
+                            else:
+                                # There are no wildcards, so the order is unambiguous.
+                                force.addTorsion(torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
                             done = True
                             break
 
@@ -1157,33 +1169,33 @@ class NonbondedGenerator:
         sys.addForce(force)
 
     def postprocessSystem(self, sys, data, args):
-        # Create exceptions based on bonds, virtual sites, and Drude particles.
+        # Create exceptions based on bonds.
+        
         bondIndices = []
         for bond in data.bonds:
             bondIndices.append((bond.atom1, bond.atom2))
+
+        # If a virtual site does *not* share exclusions with another atom, add a bond between it and its first parent atom.
+
         for i in range(sys.getNumParticles()):
             if sys.isVirtualSite(i):
-                site = sys.getVirtualSite(i)
-                for j in range(site.getNumParticles()):
-                    bondIndices.append((i, site.getParticle(j)))
-        drude = [f for f in sys.getForces() if isinstance(f, mm.DrudeForce)]
-        if len(drude) > 0:
-            drude = drude[0]
-            # For purposes of creating exceptions, a Drude particle is "bonded" to anything
-            # its parent atom is bonded to.
-            drudeMap = {}
-            for i in range(drude.getNumParticles()):
-                params = drude.getParticleParameters(i)
-                drudeMap[params[1]] = params[0]
-            for atom1, atom2 in bondIndices:
-                drude1 = drudeMap[atom1] if atom1 in drudeMap else None
-                drude2 = drudeMap[atom2] if atom2 in drudeMap else None
-                if drude1 is not None:
-                    bondIndices.append((drude1, atom2))
-                    if drude2 is not None:
-                        bondIndices.append((drude1, drude2))
-                if drude2 is not None:
-                    bondIndices.append((atom1, drude2))
+                (site, atoms, excludeWith) = data.virtualSites[data.atoms[i]]
+                if excludeWith is None:
+                    bondIndices.append((i, site.getParticle(0)))
+        
+        # Certain particles, such as lone pairs and Drude particles, share exclusions with a parent atom.
+        # If the parent atom does not interact with an atom, the child particle does not either.
+        
+        for atom1, atom2 in bondIndices:
+            for child1 in data.excludeAtomWith[atom1]:
+                bondIndices.append((child1, atom2))
+                for child2 in data.excludeAtomWith[atom2]:
+                    bondIndices.append((child1, child2))
+            for child2 in data.excludeAtomWith[atom2]:
+                bondIndices.append((atom1, child2))
+
+        # Create the exceptions.
+
         nonbonded = [f for f in sys.getForces() if isinstance(f, mm.NonbondedForce)][0]
         nonbonded.createExceptionsFromBonds(bondIndices, self.coulomb14scale, self.lj14scale)
 
@@ -1519,18 +1531,22 @@ class CustomTorsionGenerator:
                 if type1 in types1:
                     for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
                         if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
-                            # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
-                            # impropers, which leaves the ordering ambigous.  It then follows some bizarre rules
-                            # to pick the order.
-                            a1 = torsion[t2[1]]
-                            a2 = torsion[t3[1]]
-                            e1 = data.atoms[a1].element
-                            e2 = data.atoms[a2].element
-                            if e1 == e2 and a1 > a2:
-                                (a1, a2) = (a2, a1)
-                            elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
-                                (a1, a2) = (a2, a1)
-                            force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.paramValues)
+                            if wildcard in (types1, types2, types3, types4):
+                                # Workaround to be more consistent with AMBER.  It uses wildcards to define most of its
+                                # impropers, which leaves the ordering ambiguous.  It then follows some bizarre rules
+                                # to pick the order.
+                                a1 = torsion[t2[1]]
+                                a2 = torsion[t3[1]]
+                                e1 = data.atoms[a1].element
+                                e2 = data.atoms[a2].element
+                                if e1 == e2 and a1 > a2:
+                                    (a1, a2) = (a2, a1)
+                                elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
+                                    (a1, a2) = (a2, a1)
+                                force.addTorsion(a1, a2, torsion[0], torsion[t4[1]], tordef.paramValues)
+                            else:
+                                # There are no wildcards, so the order is unambiguous.
+                                force.addTorsion(torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef.paramValues)
                             done = True
                             break
 
@@ -4290,8 +4306,6 @@ class DrudeGenerator:
 
         # Add Drude particles.
 
-        drudeMap = {}
-        parentMap = {}
         for atom in data.atoms:
             t = data.atomType[atom]
             if t in self.typeMap:
@@ -4308,9 +4322,8 @@ class DrudeGenerator:
                         p[2] = atom2.index
                     elif values[3] is not None and type2 in values[3]:
                         p[3] = atom2.index
-                drudeIndex = force.addParticle(atom.index, p[0], p[1], p[2], p[3], values[4], values[5], values[6], values[7])
-                drudeMap[atom.index] = p[0]
-                parentMap[p[0]] = (atom.index, drudeIndex)
+                force.addParticle(atom.index, p[0], p[1], p[2], p[3], values[4], values[5], values[6], values[7])
+                data.excludeAtomWith[p[0]].append(atom.index)
         sys.addForce(force)
 
     def postprocessSystem(self, sys, data, args):
@@ -4323,14 +4336,14 @@ class DrudeGenerator:
             particleMap[drude.getParticleParameters(i)[0]] = i
         for i in range(nonbonded.getNumExceptions()):
             (particle1, particle2, charge, sigma, epsilon) = nonbonded.getExceptionParameters(i)
-            if charge == 0 and epsilon == 0:
+            if charge._value == 0 and epsilon._value == 0:
                 # This is an exclusion.
                 if particle1 in particleMap and particle2 in particleMap:
                     # It connects two Drude particles, so add a screened pair.
                     drude1 = particleMap[particle1]
                     drude2 = particleMap[particle2]
-                    type1 = data.atomType[data.atoms[drude1]]
-                    type2 = data.atomType[data.atoms[drude2]]
+                    type1 = data.atomType[data.atoms[particle1]]
+                    type2 = data.atomType[data.atoms[particle2]]
                     thole1 = self.typeMap[type1][8]
                     thole2 = self.typeMap[type2][8]
                     drude.addScreenedPair(drude1, drude2, thole1+thole2)