Merge branch 'master' of https://github.com/rafwiewiora/openmm

platforms/reference/src/SimTKReference/ReferenceLJCoulombIxn.cpp

@@ -395,8 +395,7 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<Vec3>& a
     double totalRealSpaceEwaldEnergy = 0.0f;
 
 
-    for (int i = 0; i < (int) neighborList->size(); i++) {
-        OpenMM::AtomPair pair = (*neighborList)[i];
+    for (auto& pair : *neighborList) {
         int ii = pair.first;
         int jj = pair.second;
 
@@ -489,10 +488,10 @@ void ReferenceLJCoulombIxn::calculateEwaldIxn(int numberOfAtoms, vector<Vec3>& a
     double totalExclusionEnergy = 0.0f;
     const double TWO_OVER_SQRT_PI = 2/sqrt(PI_M);
     for (int i = 0; i < numberOfAtoms; i++)
-        for (set<int>::const_iterator iter = exclusions[i].begin(); iter != exclusions[i].end(); ++iter) {
-            if (*iter > i) {
+        for (int exclusion : exclusions[i]) {
+            if (exclusion > i) {
                 int ii = i;
-                int jj = *iter;
+                int jj = exclusion;
 
                 double deltaR[2][ReferenceForce::LastDeltaRIndex];
                 ReferenceForce::getDeltaR(atomCoordinates[jj], atomCoordinates[ii], deltaR[0]);
@@ -581,11 +580,9 @@ void ReferenceLJCoulombIxn::calculatePairIxn(int numberOfAtoms, vector<Vec3>& at
     if (!includeDirect)
         return;
     if (cutoff) {
-        for (int i = 0; i < (int) neighborList->size(); i++) {
-            OpenMM::AtomPair pair = (*neighborList)[i];
+        for (auto& pair : *neighborList)
             calculateOneIxn(pair.first, pair.second, atomCoordinates, atomParameters, forces, energyByAtom, totalEnergy);
     }
-    }
     else {
         for (int ii = 0; ii < numberOfAtoms; ii++) {
             // loop over atom pairs

platforms/reference/src/SimTKReference/ReferenceMonteCarloBarostat.cpp

@@ -72,15 +72,15 @@ void ReferenceMonteCarloBarostat::applyBarostat(vector<Vec3>& atomPositions, con
 
     // Loop over molecules.
 
-    for (int i = 0; i < (int) molecules.size(); i++) {
+    for (auto& molecule : molecules) {
         // Find the molecule center.
 
         Vec3 pos(0, 0, 0);
-        for (int j = 0; j < (int) molecules[i].size(); j++) {
-            Vec3& atomPos = atomPositions[molecules[i][j]];
+        for (int atom : molecule) {
+            Vec3& atomPos = atomPositions[atom];
             pos += atomPos;
         }
-        pos /= molecules[i].size();
+        pos /= molecule.size();
 
         // Move it into the first periodic box.
 
@@ -95,8 +95,8 @@ void ReferenceMonteCarloBarostat::applyBarostat(vector<Vec3>& atomPositions, con
         newPos[1] *= scaleY;
         newPos[2] *= scaleZ;
         Vec3 offset = newPos-pos;
-        for (int j = 0; j < (int) molecules[i].size(); j++) {
-            Vec3& atomPos = atomPositions[molecules[i][j]];
+        for (int atom : molecule) {
+            Vec3& atomPos = atomPositions[atom];
             atomPos += offset;
         }
     }

platforms/reference/src/SimTKReference/ReferenceNeighborList.cpp

@@ -184,10 +184,10 @@ public:
                     const map<VoxelIndex, Voxel>::const_iterator voxelEntry = voxelMap.find(voxelIndex);
                     if (voxelEntry == voxelMap.end()) continue; // no such voxel; skip
                     const Voxel& voxel = voxelEntry->second;
-                    for (Voxel::const_iterator itemIter = voxel.begin(); itemIter != voxel.end(); ++itemIter)
+                    for (auto& item : voxel)
                     {
-                        const AtomIndex atomJ = itemIter->second;
-                        const Vec3& locationJ = *itemIter->first;
+                        const AtomIndex atomJ = item.second;
+                        const Vec3& locationJ = *item.first;
                         
                         // Ignore self hits
                         if (atomI == atomJ) continue;

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -52,6 +52,8 @@ AmoebaMultipoleForce::NonbondedMethod AmoebaMultipoleForce::getNonbondedMethod()
 }
 
 void AmoebaMultipoleForce::setNonbondedMethod(AmoebaMultipoleForce::NonbondedMethod method) {
+    if (method < 0 || method > 1)
+        throw OpenMMException("AmoebaMultipoleForce: Illegal value for nonbonded method");
     nonbondedMethod = method;
 }
 

plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp

@@ -123,6 +123,8 @@ AmoebaVdwForce::NonbondedMethod AmoebaVdwForce::getNonbondedMethod() const {
 }
 
 void AmoebaVdwForce::setNonbondedMethod(NonbondedMethod method) {
+    if (method < 0 || method > 1)
+        throw OpenMMException("AmoebaVdwForce: Illegal value for nonbonded method");
     nonbondedMethod = method;
 }
 

plugins/amoeba/openmmapi/src/AmoebaVdwForceImpl.cpp

@@ -161,14 +161,14 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
     // Double loop over different atom types.
     std::string sigmaCombiningRule = force.getSigmaCombiningRule();
     std::string epsilonCombiningRule = force.getEpsilonCombiningRule();
-    for (map<pair<double, double>, int>::const_iterator class1 = classCounts.begin(); class1 != classCounts.end(); ++class1) {
+    for (auto& class1 : classCounts) {
         k = 0;
-        for (map<pair<double, double>, int>::const_iterator class2 = classCounts.begin(); class2 != classCounts.end(); ++class2) { 
+        for (auto& class2 : classCounts) { 
             // AMOEBA combining rules, copied over from the CUDA code.
-            double iSigma = class1->first.first;
-            double jSigma = class2->first.first;
-            double iEpsilon = class1->first.second;
-            double jEpsilon = class2->first.second;
+            double iSigma = class1.first.first;
+            double jSigma = class2.first.first;
+            double iEpsilon = class1.first.second;
+            double jEpsilon = class2.first.second;
             // ARITHMETIC = 1
             // GEOMETRIC  = 2
             // CUBIC-MEAN = 3
@@ -207,7 +207,7 @@ double AmoebaVdwForceImpl::calcDispersionCorrection(const System& system, const
                 epsilon = 0.0;
               }
             }
-            int count = class1->second * class2->second;
+            int count = class1.second * class2.second;
             // Below is an exact copy of stuff from the previous block.
             double rv = sigma;
             double termik = 2.0 * M_PI * count; // termik is equivalent to 2 * pi * count.

plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp

@@ -973,8 +973,8 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
         molecularQuadrupolesVec.push_back((float) quadrupole[5]);
     }
     hasQuadrupoles = false;
-    for (int i = 0; i < (int) molecularQuadrupolesVec.size(); i++)
-        if (molecularQuadrupolesVec[i] != 0.0)
+    for (auto q : molecularQuadrupolesVec)
+        if (q != 0.0)
             hasQuadrupoles = true;
     int paddedNumAtoms = cu.getPaddedNumAtoms();
     for (int i = numMultipoles; i < paddedNumAtoms; i++) {
@@ -1049,15 +1049,15 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
         allAtoms.insert(atoms.begin(), atoms.end());
         force.getCovalentMap(i, AmoebaMultipoleForce::Covalent13, atoms);
         allAtoms.insert(atoms.begin(), atoms.end());
-        for (set<int>::const_iterator iter = allAtoms.begin(); iter != allAtoms.end(); ++iter)
-            covalentFlagValues.push_back(make_int3(i, *iter, 0));
+        for (int atom : allAtoms)
+            covalentFlagValues.push_back(make_int3(i, atom, 0));
         force.getCovalentMap(i, AmoebaMultipoleForce::Covalent14, atoms);
         allAtoms.insert(atoms.begin(), atoms.end());
-        for (int j = 0; j < (int) atoms.size(); j++)
-            covalentFlagValues.push_back(make_int3(i, atoms[j], 1));
+        for (int atom : atoms)
+            covalentFlagValues.push_back(make_int3(i, atom, 1));
         force.getCovalentMap(i, AmoebaMultipoleForce::Covalent15, atoms);
-        for (int j = 0; j < (int) atoms.size(); j++)
-            covalentFlagValues.push_back(make_int3(i, atoms[j], 2));
+        for (int atom : atoms)
+            covalentFlagValues.push_back(make_int3(i, atom, 2));
         allAtoms.insert(atoms.begin(), atoms.end());
         force.getCovalentMap(i, AmoebaMultipoleForce::PolarizationCovalent11, atoms);
         allAtoms.insert(atoms.begin(), atoms.end());
@@ -1068,15 +1068,14 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
 
         vector<int> atoms12;
         force.getCovalentMap(i, AmoebaMultipoleForce::PolarizationCovalent12, atoms12);
-        for (int j = 0; j < (int) atoms.size(); j++)
-            if (find(atoms12.begin(), atoms12.end(), atoms[j]) == atoms12.end())
-                polarizationFlagValues.push_back(make_int2(i, atoms[j]));
+        for (int atom : atoms)
+            if (find(atoms12.begin(), atoms12.end(), atom) == atoms12.end())
+                polarizationFlagValues.push_back(make_int2(i, atom));
     }
     set<pair<int, int> > tilesWithExclusions;
     for (int atom1 = 0; atom1 < (int) exclusions.size(); ++atom1) {
         int x = atom1/CudaContext::TileSize;
-        for (int j = 0; j < (int) exclusions[atom1].size(); ++j) {
-            int atom2 = exclusions[atom1][j];
+        for (int atom2 : exclusions[atom1]) {
             int y = atom2/CudaContext::TileSize;
             tilesWithExclusions.insert(make_pair(max(x, y), min(x, y)));
         }
@@ -1412,10 +1411,10 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
     }
     covalentFlags = CudaArray::create<uint2>(cu, nb.getExclusions().getSize(), "covalentFlags");
     vector<uint2> covalentFlagsVec(nb.getExclusions().getSize(), make_uint2(0, 0));
-    for (int i = 0; i < (int) covalentFlagValues.size(); i++) {
-        int atom1 = covalentFlagValues[i].x;
-        int atom2 = covalentFlagValues[i].y;
-        int value = covalentFlagValues[i].z;
+    for (int3 values : covalentFlagValues) {
+        int atom1 = values.x;
+        int atom2 = values.y;
+        int value = values.z;
         int x = atom1/CudaContext::TileSize;
         int offset1 = atom1-x*CudaContext::TileSize;
         int y = atom2/CudaContext::TileSize;
@@ -1446,9 +1445,9 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
     
     polarizationGroupFlags = CudaArray::create<unsigned int>(cu, nb.getExclusions().getSize(), "polarizationGroupFlags");
     vector<unsigned int> polarizationGroupFlagsVec(nb.getExclusions().getSize(), 0);
-    for (int i = 0; i < (int) polarizationFlagValues.size(); i++) {
-        int atom1 = polarizationFlagValues[i].x;
-        int atom2 = polarizationFlagValues[i].y;
+    for (int2 values : polarizationFlagValues) {
+        int atom1 = values.x;
+        int atom2 = values.y;
         int x = atom1/CudaContext::TileSize;
         int offset1 = atom1-x*CudaContext::TileSize;
         int y = atom2/CudaContext::TileSize;
@@ -1473,10 +1472,12 @@ void CudaCalcAmoebaMultipoleForceKernel::initializeScaleFactors() {
 double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool includeForces, bool includeEnergy) {
     if (!hasInitializedScaleFactors) {
         initializeScaleFactors();
-        for (int i = 0; i < (int) context.getForceImpls().size() && gkKernel == NULL; i++) {
-            AmoebaGeneralizedKirkwoodForceImpl* gkImpl = dynamic_cast<AmoebaGeneralizedKirkwoodForceImpl*>(context.getForceImpls()[i]);
-            if (gkImpl != NULL)
+        for (auto impl : context.getForceImpls()) {
+            AmoebaGeneralizedKirkwoodForceImpl* gkImpl = dynamic_cast<AmoebaGeneralizedKirkwoodForceImpl*>(impl);
+            if (gkImpl != NULL) {
                 gkKernel = dynamic_cast<CudaCalcAmoebaGeneralizedKirkwoodForceKernel*>(&gkImpl->getKernel().getImpl());
+                break;
+            }
         }
     }
     CudaNonbondedUtilities& nb = cu.getNonbondedUtilities();
@@ -2232,8 +2233,8 @@ void CudaCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& co
         molecularQuadrupolesVec.push_back((float) quadrupole[5]);
     }
     if (!hasQuadrupoles) {
-        for (int i = 0; i < (int) molecularQuadrupolesVec.size(); i++)
-            if (molecularQuadrupolesVec[i] != 0.0)
+        for (auto q : molecularQuadrupolesVec)
+            if (q != 0.0)
                 throw OpenMMException("updateParametersInContext: Cannot set a non-zero quadrupole moment, because quadrupoles were excluded from the kernel");
     }
     for (int i = force.getNumMultipoles(); i < cu.getPaddedNumAtoms(); i++) {

plugins/amoeba/platforms/cuda/src/kernels/multipoles.cu

@@ -27,12 +27,16 @@ extern "C" __global__ void computeLabFrameMoments(real4* __restrict__ posq, int4
         if (particles.z >= 0) {
             real4 thisParticlePos = posq[atom];
             real4 posZ = posq[particles.z];
-            real3 vectorZ = make_real3(posZ.x-thisParticlePos.x, posZ.y-thisParticlePos.y, posZ.z-thisParticlePos.z);
+            real3 vectorZ = normalize(make_real3(posZ.x-thisParticlePos.x, posZ.y-thisParticlePos.y, posZ.z-thisParticlePos.z));
             int axisType = particles.w; 
             real4 posX;
             real3 vectorX;
-            if (axisType >= 4)
-                vectorX = make_real3((real) 0.1f);
+            if (axisType >= 4) {
+                if (fabs(vectorZ.x) < 0.866)
+                    vectorX = make_real3(1, 0, 0);
+                else
+                    vectorX = make_real3(0, 1, 0);
+            }
             else {
                 posX = posq[particles.x];
                 vectorX = make_real3(posX.x-thisParticlePos.x, posX.y-thisParticlePos.y, posX.z-thisParticlePos.z);
@@ -81,8 +85,6 @@ extern "C" __global__ void computeLabFrameMoments(real4* __restrict__ posq, int4
         
             // branch based on axis type
                     
-            vectorZ = normalize(vectorZ);
-        
             if (axisType == 1) {
         
                 // bisector
@@ -362,8 +364,12 @@ extern "C" __global__ void mapTorqueToForce(unsigned long long* __restrict__ for
             norms[U] = normVector(vector[U]);
             if (axisType != 4 && particles.x >= 0)
                 vector[V] = atomPos - trimTo3(posq[particles.x]);
+            else {
+                if (fabs(vector[U].x/norms[U]) < 0.866)
+                    vector[V] = make_real3(1, 0, 0);
                 else
-                vector[V] = make_real3(0.1f);
+                    vector[V] = make_real3(0, 1, 0);
+            }
             norms[V] = normVector(vector[V]);
         
             // W = UxV
@@ -488,7 +494,7 @@ extern "C" __global__ void mapTorqueToForce(unsigned long long* __restrict__ for
             else if (axisType == 4) {
                 // z-only
         
-                forces[Z] = vector[UV]*dphi[V]/(norms[U]*angles[UV][1]);
+                forces[Z] = vector[UV]*dphi[V]/(norms[U]*angles[UV][1]) + vector[UW]*dphi[W]/norms[U];
                 forces[X] = make_real3(0);
                 forces[Y] = make_real3(0);
                 forces[I] = -forces[Z];

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaExtrapolatedPolarization.cpp

@@ -290,8 +290,8 @@ static void check_finite_differences(vector<Vec3> analytic_forces, Context &cont
 
     // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
     double norm = 0.0;
-    for (int i = 0; i < (int) analytic_forces.size(); ++i)
-        norm += analytic_forces[i].dot(analytic_forces[i]);
+    for (auto& f : analytic_forces)
+        norm += f.dot(f);
     norm = std::sqrt(norm);
     const double stepSize = 1e-3;
     double step = 0.5*stepSize/norm;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaGeneralizedKirkwoodForce.cpp

@@ -61,8 +61,8 @@ static void checkFiniteDifferences(vector<Vec3> forces, Context &context, vector
     // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
     double norm = 0.0;
-    for (int i = 0; i < (int) forces.size(); ++i)
-        norm += forces[i].dot(forces[i]);
+    for (auto& f : forces)
+        norm += f.dot(f);
     norm = std::sqrt(norm);
     const double stepSize = 1e-3;
     double step = 0.5*stepSize/norm;

plugins/amoeba/platforms/cuda/tests/TestCudaAmoebaMultipoleForce.cpp

@@ -2919,8 +2919,8 @@ static void testNoQuadrupoles(bool usePme) {
         int axisType, atomX, atomY, atomZ;
         vector<double> dipole, quadrupole;
         amoebaMultipoleForce->getMultipoleParameters(i, charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
-        for (int j = 0; j < (int) quadrupole.size(); j++)
-            quadrupole[j] = 0;
+        for (auto& q : quadrupole)
+            q = 0;
         amoebaMultipoleForce->setMultipoleParameters(i, charge, dipole, quadrupole, axisType, atomZ, atomX, atomY, thole, damping, polarity);
     }
     amoebaMultipoleForce->updateParametersInContext(context);
@@ -3223,6 +3223,55 @@ void testZBisect() {
     ASSERT_EQUAL_TOL(-84.1532, state.getPotentialEnergy(), 0.01);
 }
 
+void testZOnly() {
+    int numParticles = 3;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    AmoebaMultipoleForce* force = new AmoebaMultipoleForce();
+    system.addForce(force);
+    vector<double> d(3), q(9, 0.0);
+    d[0] = 0.05;
+    d[1] = -0.05;
+    d[2] = 0.1;
+    force->addMultipole(0.0, d, q, AmoebaMultipoleForce::ZOnly, 1, 0, 0, 0.39, 0.33, 0.001);
+    force->addMultipole(0.0, d, q, AmoebaMultipoleForce::Bisector, 0, 2, 0, 0.39, 0.33, 0.001);
+    force->addMultipole(0.0, d, q, AmoebaMultipoleForce::ZOnly, 1, 0, 0, 0.39, 0.33, 0.001);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0.2, 0, 0);
+    positions[2] = Vec3(0.2, 0.2, -0.05);
+    
+    // Evaluate the forces.
+    
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    context.setPositions(positions);
+    State state = context.getState(State::Forces);
+    double norm = 0.0;
+    for (Vec3 f : state.getForces())
+        norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+    norm = std::sqrt(norm);
+
+    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+    const double delta = 1e-3;
+    double step = 0.5*delta/norm;
+    vector<Vec3> positions2(numParticles), positions3(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = state.getForces()[i];
+        positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
+    }
+    context.setPositions(positions2);
+    State state2 = context.getState(State::Energy);
+    context.setPositions(positions3);
+    State state3 = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(state2.getPotentialEnergy(), state3.getPotentialEnergy()+norm*delta, 1e-3)
+}
+
+
 int main(int argc, char* argv[]) {
     try {
         std::cout << "TestCudaAmoebaMultipoleForce running test..." << std::endl;
@@ -3280,6 +3329,10 @@ int main(int argc, char* argv[]) {
         
         testZBisect();
         
+        // test the ZOnly axis type.
+        
+        testZOnly();
+
     } catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;
         std::cout << "FAIL - ERROR.  Test failed." << std::endl;

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceVdwForce.cpp

@@ -277,9 +277,8 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
  
         double sigmaI      = sigmas[ii];
         double epsilonI    = epsilons[ii];
-        for (std::set<int>::const_iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++) {
-            exclusions[*jj] = 1;
-        }
+        for (int jj : allExclusions[ii])
+            exclusions[jj] = 1;
 
         for (unsigned int jj = ii+1; jj < static_cast<unsigned int>(numParticles); jj++) {
             if (exclusions[jj] == 0) {
@@ -310,9 +309,8 @@ double AmoebaReferenceVdwForce::calculateForceAndEnergy(int numParticles,
             }
         }
 
-        for (std::set<int>::const_iterator jj = allExclusions[ii].begin(); jj != allExclusions[ii].end(); jj++) {
-            exclusions[*jj] = 0;
-        }
+        for (int jj : allExclusions[ii])
+            exclusions[jj] = 0;
     }
 
     return energy;

plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaExtrapolatedPolarization.cpp

@@ -287,8 +287,8 @@ static void check_finite_differences(vector<Vec3> analytic_forces, Context &cont
     // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
 
     double norm = 0.0;
-    for (int i = 0; i < (int) analytic_forces.size(); ++i)
-        norm += analytic_forces[i].dot(analytic_forces[i]);
+    for (auto& f : analytic_forces)
+        norm += f.dot(f);
     norm = std::sqrt(norm);
     const double stepSize = 1e-3;
     double step = 0.5*stepSize/norm;

plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaGeneralizedKirkwoodForce.cpp

@@ -61,8 +61,8 @@ static void checkFiniteDifferences(vector<Vec3> forces, Context &context, vector
     // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
     double norm = 0.0;
-    for (int i = 0; i < (int) forces.size(); ++i)
-        norm += forces[i].dot(forces[i]);
+    for (auto& f : forces)
+        norm += f.dot(f);
     norm = std::sqrt(norm);
     const double stepSize = 1e-3;
     double step = 0.5*stepSize/norm;

plugins/amoeba/platforms/reference/tests/TestReferenceAmoebaMultipoleForce.cpp

@@ -3030,6 +3030,54 @@ void testZBisect() {
     ASSERT_EQUAL_TOL(-84.1532, state.getPotentialEnergy(), 0.01);
 }
 
+void testZOnly() {
+    int numParticles = 3;
+    System system;
+    for (int i = 0; i < numParticles; i++)
+        system.addParticle(1.0);
+    AmoebaMultipoleForce* force = new AmoebaMultipoleForce();
+    system.addForce(force);
+    vector<double> d(3), q(9, 0.0);
+    d[0] = 0.05;
+    d[1] = -0.05;
+    d[2] = 0.1;
+    force->addMultipole(0.0, d, q, AmoebaMultipoleForce::ZOnly, 1, 0, 0, 0.39, 0.33, 0.001);
+    force->addMultipole(0.0, d, q, AmoebaMultipoleForce::Bisector, 0, 2, 0, 0.39, 0.33, 0.001);
+    force->addMultipole(0.0, d, q, AmoebaMultipoleForce::ZOnly, 1, 0, 0, 0.39, 0.33, 0.001);
+    vector<Vec3> positions(3);
+    positions[0] = Vec3(0, 0, 0);
+    positions[1] = Vec3(0.2, 0, 0);
+    positions[2] = Vec3(0.2, 0.2, -0.05);
+    
+    // Evaluate the forces.
+    
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    context.setPositions(positions);
+    State state = context.getState(State::Forces);
+    double norm = 0.0;
+    for (Vec3 f : state.getForces())
+        norm += f[0]*f[0] + f[1]*f[1] + f[2]*f[2];
+    norm = std::sqrt(norm);
+
+    // Take a small step in the direction of the energy gradient and see whether the potential energy changes by the expected amount.
+
+    const double delta = 1e-3;
+    double step = 0.5*delta/norm;
+    vector<Vec3> positions2(numParticles), positions3(numParticles);
+    for (int i = 0; i < numParticles; ++i) {
+        Vec3 p = positions[i];
+        Vec3 f = state.getForces()[i];
+        positions2[i] = Vec3(p[0]-f[0]*step, p[1]-f[1]*step, p[2]-f[2]*step);
+        positions3[i] = Vec3(p[0]+f[0]*step, p[1]+f[1]*step, p[2]+f[2]*step);
+    }
+    context.setPositions(positions2);
+    State state2 = context.getState(State::Energy);
+    context.setPositions(positions3);
+    State state3 = context.getState(State::Energy);
+    ASSERT_EQUAL_TOL(state2.getPotentialEnergy(), state3.getPotentialEnergy()+norm*delta, 1e-3)
+}
+
 int main(int numberOfArguments, char* argv[]) {
 
     try {
@@ -3083,6 +3131,10 @@ int main(int numberOfArguments, char* argv[]) {
         // test the ZBisect axis type.
         
         testZBisect();
+        
+        // test the ZOnly axis type.
+        
+        testZOnly();
     }
     catch(const std::exception& e) {
         std::cout << "exception: " << e.what() << std::endl;

plugins/amoeba/serialization/src/AmoebaBondForceProxy.cpp

@@ -72,11 +72,9 @@ void* AmoebaBondForceProxy::deserialize(const SerializationNode& node) const {
         force->setAmoebaGlobalBondCubic(node.getDoubleProperty("cubic"));
         force->setAmoebaGlobalBondQuartic(node.getDoubleProperty("quartic"));
         const SerializationNode& bonds = node.getChildNode("Bonds");
-        for (unsigned int ii = 0; ii < (int) bonds.getChildren().size(); ii++) {
-            const SerializationNode& bond = bonds.getChildren()[ii];
+        for (auto& bond : bonds.getChildren())
             force->addBond(bond.getIntProperty("p1"), bond.getIntProperty("p2"), bond.getDoubleProperty("d"), bond.getDoubleProperty("k"));
     }
-    }
     catch (...) {
         delete force;
         throw;

plugins/amoeba/serialization/src/AmoebaStretchBendForceProxy.cpp

@@ -67,8 +67,7 @@ void* AmoebaStretchBendForceProxy::deserialize(const SerializationNode& node) co
         if (version > 3)
             force->setUsesPeriodicBoundaryConditions(node.getBoolProperty("usesPeriodic"));
         const SerializationNode& bonds = node.getChildNode("StretchBendAngles");
-        for (unsigned int ii = 0; ii < (int) bonds.getChildren().size(); ii++) {
-            const SerializationNode& bond = bonds.getChildren()[ii];
+        for (auto& bond : bonds.getChildren()) {
             double k1, k2;
             if (version == 1)
                 k1 = k2 = bond.getDoubleProperty("k");

plugins/cpupme/src/CpuPmeKernels.cpp

@@ -515,8 +515,8 @@ CpuCalcPmeReciprocalForceKernel::~CpuCalcPmeReciprocalForceKernel() {
     pthread_mutex_destroy(&lock);
     pthread_cond_destroy(&startCondition);
     pthread_cond_destroy(&endCondition);
-    for (int i = 0; i < (int) tempGrid.size(); i++)
-        fftwf_free(tempGrid[i]);
+    for (auto grid : tempGrid)
+        fftwf_free(grid);
     if (complexGrid != NULL)
         fftwf_free(complexGrid);
     if (hasCreatedPlan) {
@@ -552,8 +552,8 @@ void CpuCalcPmeReciprocalForceKernel::runMainThread() {
         if (includeEnergy) {
             threads.resumeThreads(); // Signal threads to compute energy.
             threads.waitForThreads();
-            for (int i = 0; i < (int) threadEnergy.size(); i++)
-                energy += threadEnergy[i];
+            for (auto e : threadEnergy)
+                energy += e;
         }
         threads.resumeThreads(); // Signal threads to perform reciprocal convolution.
         threads.waitForThreads();
@@ -805,8 +805,8 @@ CpuCalcDispersionPmeReciprocalForceKernel::~CpuCalcDispersionPmeReciprocalForceK
     pthread_mutex_destroy(&lock);
     pthread_cond_destroy(&startCondition);
     pthread_cond_destroy(&endCondition);
-    for (int i = 0; i < (int) tempGrid.size(); i++)
-        fftwf_free(tempGrid[i]);
+    for (auto grid : tempGrid)
+        fftwf_free(grid);
     if (complexGrid != NULL)
         fftwf_free(complexGrid);
     if (hasCreatedPlan) {
@@ -843,8 +843,8 @@ void CpuCalcDispersionPmeReciprocalForceKernel::runMainThread() {
         if (includeEnergy) {
             threads.resumeThreads(); // Signal threads to compute energy.
             threads.waitForThreads();
-            for (int i = 0; i < (int) threadEnergy.size(); i++)
-                energy += threadEnergy[i];
+            for (auto e : threadEnergy)
+                energy += e;
         }
         threads.resumeThreads(); // Signal threads to perform reciprocal convolution.
         threads.waitForThreads();

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

@@ -267,8 +267,7 @@ void ReferenceIntegrateDrudeLangevinStepKernel::execute(ContextImpl& context, co
     const double fscale = (1-vscale)/integrator.getFriction();
     const double kT = BOLTZ*integrator.getTemperature();
     const double noisescale = sqrt(2*kT*integrator.getFriction())*sqrt(0.5*(1-vscale*vscale)/integrator.getFriction());
-    for (int i = 0; i < (int) normalParticles.size(); i++) {
-        int index = normalParticles[i];
+    for (int index : normalParticles) {
         double invMass = particleInvMass[index];
         if (invMass != 0.0) {
             double sqrtInvMass = sqrt(invMass);