Merge pull request #1646 from peastman/zonly

Fixed errors with multipoles set to ZOnly mode

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

@@ -24,13 +24,19 @@ extern "C" __global__ void computeLabFrameMoments(real4* __restrict__ posq, int4
         // code common to ZThenX and Bisector
         
         int4 particles = multipoleParticles[atom];
-        if (particles.x >= 0 && particles.z >= 0) {
+        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);
-            real4 posX = posq[particles.x];
-            real3 vectorX = make_real3(posX.x-thisParticlePos.x, posX.y-thisParticlePos.y, posX.z-thisParticlePos.z);
             int axisType = particles.w; 
+            real4 posX;
+            real3 vectorX;
+            if (axisType >= 4)
+                vectorX = make_real3((real) 0.1f);
+            else {
+                posX = posq[particles.x];
+                vectorX = make_real3(posX.x-thisParticlePos.x, posX.y-thisParticlePos.y, posX.z-thisParticlePos.z);
+            }
     
             /*
                 z-only
@@ -108,8 +114,6 @@ extern "C" __global__ void computeLabFrameMoments(real4* __restrict__ posq, int4
                 }
          
             }
-            else if (axisType >= 4)
-                vectorX = make_real3((real) 0.1f);
             
             // x = x - (x.z)z
         

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

@@ -397,8 +397,8 @@ void AmoebaReferenceMultipoleForce::checkChiral(vector<MultipoleParticleData>& p
 }
 
 void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle(      MultipoleParticleData& particleI,
-                                                                  const MultipoleParticleData& particleZ,
-                                                                  const MultipoleParticleData& particleX,
+                                                                  const MultipoleParticleData* particleZ,
+                                                                  const MultipoleParticleData* particleX,
                                                                         MultipoleParticleData* particleY,
                                                                         int axisType) const
 {
@@ -410,14 +410,20 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle(      Multipol
     // this atom and the axis atom
 
 
-    RealVec vectorY;
-    RealVec vectorZ = particleZ.position - particleI.position;
-    RealVec vectorX = particleX.position - particleI.position;
-
+    RealVec vectorX, vectorY;
+    RealVec vectorZ = particleZ->position - particleI.position;
     normalizeRealVec(vectorZ);
 
     // branch based on axis type
 
+    if (axisType == AmoebaMultipoleForce::ZOnly) {
+
+        // z-only
+
+        vectorX = RealVec(0.1, 0.1, 0.1);
+    }
+    else {
+        vectorX = particleX->position - particleI.position;
         if (axisType == AmoebaMultipoleForce::Bisector) {
 
             // bisector
@@ -427,8 +433,8 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle(      Multipol
             normalizeRealVec(vectorX);
             vectorZ      += vectorX;
             normalizeRealVec(vectorZ);
-
-    } else if (axisType == AmoebaMultipoleForce::ZBisect) {
+        }
+        else if (axisType == AmoebaMultipoleForce::ZBisect) {
 
             // z-bisect
 
@@ -441,8 +447,8 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle(      Multipol
 
             vectorX += vectorY;
             normalizeRealVec(vectorX);
-
-    } else if (axisType == AmoebaMultipoleForce::ThreeFold) {
+        }
+        else if (axisType == AmoebaMultipoleForce::ThreeFold) {
 
             // 3-fold
 
@@ -455,13 +461,7 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrixToParticle(      Multipol
 
             vectorZ  += vectorX +  vectorY;
             normalizeRealVec(vectorZ);
-
-    } else if (axisType == AmoebaMultipoleForce::ZOnly) {
-
-        // z-only
-
-        vectorX = RealVec(0.1, 0.1, 0.1);
-
+        }
     }
 
     RealOpenMM dot      = vectorZ.dot(vectorX);
@@ -652,8 +652,9 @@ void AmoebaReferenceMultipoleForce::applyRotationMatrix(vector<MultipoleParticle
 {
 
     for (unsigned int ii = 0; ii < _numParticles; ii++) {
-        if (multipoleAtomZs[ii] >= 0 && multipoleAtomXs[ii] >= 0) {
-            applyRotationMatrixToParticle(particleData[ii], particleData[multipoleAtomZs[ii]], particleData[multipoleAtomXs[ii]],
+        if (multipoleAtomZs[ii] >= 0) {
+            applyRotationMatrixToParticle(particleData[ii], &particleData[multipoleAtomZs[ii]],
+                                          multipoleAtomXs[ii] > -1 ? &particleData[multipoleAtomXs[ii]] : NULL,
                                           multipoleAtomYs[ii] > -1 ? &particleData[multipoleAtomYs[ii]] : NULL, axisTypes[ii]);
         }
     }

plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.h

@@ -900,8 +900,8 @@ protected:
      * @param  axisType             axis type
      */
     void applyRotationMatrixToParticle(      MultipoleParticleData& particleI,
-                                       const MultipoleParticleData& particleZ,
-                                       const MultipoleParticleData& particleX,
+                                       const MultipoleParticleData* particleZ,
+                                       const MultipoleParticleData* particleX,
                                              MultipoleParticleData* particleY, int axisType) const;
 
     /**

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

@@ -4907,7 +4907,7 @@ class AmoebaMultipoleGenerator(object):
                         bondedAtomZ = data.atoms[bondedAtomZIndex]
 
                         if (kx == 0 and kz == bondedAtomZType):
-                            kz = bondedAtomZIndex
+                            zaxis = bondedAtomZIndex
                             savedMultipoleDict = multipoleDict
                             hit = 5