Mods to output of forces

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernelFactory.cpp

@@ -37,18 +37,20 @@ extern "C" void OPENMM_EXPORT registerKernelFactories() {
     for( int ii = 0; ii < Platform::getNumPlatforms(); ii++ ){
         Platform& platform = Platform::getPlatform(ii);
         if( platform.getName() == "Reference" ){
+
              AmoebaReferenceKernelFactory* factory = new AmoebaReferenceKernelFactory();
-             platform.registerKernelFactory(CalcAmoebaHarmonicBondForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaUreyBradleyForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaHarmonicAngleForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaPiTorsionForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaOutOfPlaneBendForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(), factory);
-             platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(), factory);
+
+             platform.registerKernelFactory(CalcAmoebaHarmonicBondForceKernel::Name(),          factory);
+             platform.registerKernelFactory(CalcAmoebaUreyBradleyForceKernel::Name(),           factory);
+             platform.registerKernelFactory(CalcAmoebaHarmonicAngleForceKernel::Name(),         factory);
+             platform.registerKernelFactory(CalcAmoebaHarmonicInPlaneAngleForceKernel::Name(),  factory);
+             platform.registerKernelFactory(CalcAmoebaTorsionForceKernel::Name(),               factory);
+             platform.registerKernelFactory(CalcAmoebaPiTorsionForceKernel::Name(),             factory);
+             platform.registerKernelFactory(CalcAmoebaStretchBendForceKernel::Name(),           factory);
+             platform.registerKernelFactory(CalcAmoebaOutOfPlaneBendForceKernel::Name(),        factory);
+             platform.registerKernelFactory(CalcAmoebaTorsionTorsionForceKernel::Name(),        factory);
+             platform.registerKernelFactory(CalcAmoebaVdwForceKernel::Name(),                   factory);
+             platform.registerKernelFactory(CalcAmoebaMultipoleForceKernel::Name(),             factory);
 /*
              platform.registerKernelFactory(CalcAmoebaGeneralizedKirkwoodForceKernel::Name(), factory);
 */

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

@@ -175,9 +175,9 @@ RealOpenMM AmoebaReferenceHarmonicAngleForce::calculateAngleIxn( const RealVec&
     // accumulate forces
  
     for( int jj = 0; jj < 3; jj++ ){
-        forces[jj][0] += deltaCrossP[jj][0];
-        forces[jj][1] += deltaCrossP[jj][1];
-        forces[jj][2] += deltaCrossP[jj][2];
+        forces[jj][0] = deltaCrossP[jj][0];
+        forces[jj][1] = deltaCrossP[jj][1];
+        forces[jj][2] = deltaCrossP[jj][2];
     }
  
     return energy;
@@ -202,11 +202,14 @@ RealOpenMM AmoebaReferenceHarmonicAngleForce::calculateForceAndEnergy( int numAn
         RealOpenMM idealAngle   = angle[ii];
         RealOpenMM angleK       = kQuadratic[ii];
         RealVec forces[3];
-        forces[0]               = forceData[particle1Index];
-        forces[1]               = forceData[particle2Index];
-        forces[2]               = forceData[particle3Index];
         energy                 += calculateAngleIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index],
                                                      idealAngle, angleK, angleCubic, angleQuartic, anglePentic, angleSextic, forces );
+
+        for( unsigned int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] += forces[0][jj];
+            forceData[particle2Index][jj] += forces[1][jj];
+            forceData[particle3Index][jj] += forces[2][jj];
+        }
     }   
     return energy;
 }

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

@@ -75,13 +75,14 @@ RealOpenMM AmoebaReferenceHarmonicBondForce::calculateBondIxn( const RealVec& po
    dEdR                      *= two*bondK*deltaIdeal;
    dEdR                       = r > zero ? (dEdR/r) : zero;
 
-   forces[0][0]              += dEdR*deltaR[0];
-   forces[0][1]              += dEdR*deltaR[1];
-   forces[0][2]              += dEdR*deltaR[2];
+   forces[0][0]               = dEdR*deltaR[0];
+   forces[0][1]               = dEdR*deltaR[1];
+   forces[0][2]               = dEdR*deltaR[2];
 
-   forces[1][0]              -= dEdR*deltaR[0];
-   forces[1][1]              -= dEdR*deltaR[1];
-   forces[1][2]              -= dEdR*deltaR[2];
+   dEdR                      *= -1.0;
+   forces[1][0]               = dEdR*deltaR[0];
+   forces[1][1]               = dEdR*deltaR[1];
+   forces[1][2]               = dEdR*deltaR[2];
 
    RealOpenMM energy          = bondK*deltaIdeal2*( one + bondCubic*deltaIdeal + bondQuartic*deltaIdeal2 );
    return energy;
@@ -103,11 +104,16 @@ RealOpenMM AmoebaReferenceHarmonicBondForce::calculateForceAndEnergy( int numBon
         RealOpenMM bondLength   = length[ii];
         RealOpenMM bondK        = kQuadratic[ii];
         RealVec forces[2];
-        forces[0]               = forceData[particle1Index];
-        forces[1]               = forceData[particle2Index];
+
         energy                 += calculateBondIxn( particlePositions[particle1Index], particlePositions[particle2Index],
                                                     bondLength, bondK, globalHarmonicBondCubic, globalHarmonicBondQuartic,
                                                     forces );
+
+        for( int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] += forces[0][jj];
+            forceData[particle2Index][jj] += forces[1][jj];
+        }
+
     }   
     return energy;
 }

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

@@ -236,9 +236,9 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateAngleIxn( const Re
     // accumulate forces
  
     for( int jj = 0; jj < 4; jj++ ){
-        forces[jj][0] -= forceTerm[jj][0];
-        forces[jj][1] -= forceTerm[jj][1];
-        forces[jj][2] -= forceTerm[jj][2];
+        forces[jj][0] = forceTerm[jj][0];
+        forces[jj][1] = forceTerm[jj][1];
+        forces[jj][2] = forceTerm[jj][2];
     }
  
     return energy;
@@ -266,12 +266,18 @@ RealOpenMM AmoebaReferenceHarmonicInPlaneAngleForce::calculateForceAndEnergy( in
         RealOpenMM idealAngle   = angle[ii];
         RealOpenMM angleK       = kQuadratic[ii];
         RealVec forces[4];
-        forces[0]               = forceData[particle1Index];
-        forces[1]               = forceData[particle2Index];
-        forces[2]               = forceData[particle3Index];
-        forces[3]               = forceData[particle4Index];
         energy                 += calculateAngleIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
                                                      idealAngle, angleK, angleCubic, angleQuartic, anglePentic, angleSextic, forces );
+
+        // accumulate forces
+     
+        for( int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] -= forces[0][jj];
+            forceData[particle2Index][jj] -= forces[1][jj];
+            forceData[particle3Index][jj] -= forces[2][jj];
+            forceData[particle4Index][jj] -= forces[3][jj];
+        }
+ 
     }   
     return energy;
 }

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

@@ -182,7 +182,7 @@ RealOpenMM AmoebaReferenceOutOfPlaneBendForce::calculateOutOfPlaneBendIxn( const
  
     for( int jj = 0; jj < LastAtomIndex; jj++ ){
        for( int ii = 0; ii < 3; ii++ ){
-          forces[jj][ii] -= subForce[jj][ii];
+          forces[jj][ii] = subForce[jj][ii];
        }
     }
  
@@ -215,12 +215,15 @@ RealOpenMM AmoebaReferenceOutOfPlaneBendForce::calculateForceAndEnergy( int numO
         int particle4Index      = particle4[ii];
         RealOpenMM kAngle       = kQuadratic[ii];
         RealVec forces[4];
-        forces[0]               = forceData[particle1Index];
-        forces[1]               = forceData[particle2Index];
-        forces[2]               = forceData[particle3Index];
-        forces[3]               = forceData[particle4Index];
         energy                 += calculateOutOfPlaneBendIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
                                                               kAngle, angleCubic, angleQuartic, anglePentic, angleSextic, forces );
+        for( int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] -= forces[0][jj];
+            forceData[particle2Index][jj] -= forces[1][jj];
+            forceData[particle3Index][jj] -= forces[2][jj];
+            forceData[particle4Index][jj] -= forces[3][jj];
+        }
+
     }   
     return energy;
 }

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

@@ -154,29 +154,29 @@ RealOpenMM AmoebaReferencePiTorsionForce::calculatePiTorsionIxn( const RealVec&
  
     // add in forces
  
-    forces[0][0] -= d[0][0];
-    forces[0][1] -= d[0][1];
-    forces[0][2] -= d[0][2];
+    forces[0][0] = d[0][0];
+    forces[0][1] = d[0][1];
+    forces[0][2] = d[0][2];
  
-    forces[1][0] -= d[1][0];
-    forces[1][1] -= d[1][1];
-    forces[1][2] -= d[1][2];
+    forces[1][0] = d[1][0];
+    forces[1][1] = d[1][1];
+    forces[1][2] = d[1][2];
  
-    forces[2][0] -= d[2][0];
-    forces[2][1] -= d[2][1];
-    forces[2][2] -= d[2][2];
+    forces[2][0] = d[2][0];
+    forces[2][1] = d[2][1];
+    forces[2][2] = d[2][2];
  
-    forces[3][0] -= d[3][0];
-    forces[3][1] -= d[3][1];
-    forces[3][2] -= d[3][2];
+    forces[3][0] = d[3][0];
+    forces[3][1] = d[3][1];
+    forces[3][2] = d[3][2];
  
-    forces[4][0] -= d[4][0];
-    forces[4][1] -= d[4][1];
-    forces[4][2] -= d[4][2];
+    forces[4][0] = d[4][0];
+    forces[4][1] = d[4][1];
+    forces[4][2] = d[4][2];
  
-    forces[5][0] -= d[5][0];
-    forces[5][1] -= d[5][1];
-    forces[5][2] -= d[5][2];
+    forces[5][0] = d[5][0];
+    forces[5][1] = d[5][1];
+    forces[5][2] = d[5][2];
  
     // ---------------------------------------------------------------------------------------
  
@@ -206,16 +206,21 @@ RealOpenMM AmoebaReferencePiTorsionForce::calculateForceAndEnergy( int numPiTors
         int particle6Index      = particle6[ii];
 
         RealVec forces[6];
-        forces[0]               = forceData[particle1Index];
-        forces[1]               = forceData[particle2Index];
-        forces[2]               = forceData[particle3Index];
-        forces[3]               = forceData[particle4Index];
-        forces[4]               = forceData[particle5Index];
-        forces[5]               = forceData[particle6Index];
         energy                 += calculatePiTorsionIxn( posData[particle1Index], posData[particle2Index],
                                                          posData[particle3Index], posData[particle4Index],
                                                          posData[particle5Index], posData[particle6Index],
                                                          kTorsion[ii], forces );
+        // accumulate forces
+     
+        for( int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] -= forces[0][jj];
+            forceData[particle2Index][jj] -= forces[1][jj];
+            forceData[particle3Index][jj] -= forces[2][jj];
+            forceData[particle4Index][jj] -= forces[3][jj];
+            forceData[particle5Index][jj] -= forces[4][jj];
+            forceData[particle6Index][jj] -= forces[5][jj];
+        }   
+
     }   
     return energy;
 }

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

@@ -137,9 +137,9 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateStretchBendIxn( const RealV
     // add in forces
  
     for( int jj = 0; jj < LastAtomIndex; jj++ ){
-        forces[jj][0] -= subForce[jj][0];
-        forces[jj][1] -= subForce[jj][1];
-        forces[jj][2] -= subForce[jj][2];
+        forces[jj][0] = subForce[jj][0];
+        forces[jj][1] = subForce[jj][1];
+        forces[jj][2] = subForce[jj][2];
     }
  
     // ---------------------------------------------------------------------------------------
@@ -166,11 +166,16 @@ RealOpenMM AmoebaReferenceStretchBendForce::calculateForceAndEnergy( int numStre
         RealOpenMM idealAngle   = angle[ii];
         RealOpenMM angleK       = kQuadratic[ii];
         RealVec forces[3];
-        forces[0]               = forceData[particle1Index];
-        forces[1]               = forceData[particle2Index];
-        forces[2]               = forceData[particle3Index];
         energy                 += calculateStretchBendIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index],
                                                            abLength, cbLength, idealAngle, angleK, forces );
+        // accumulate forces
+    
+        for( int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] -= forces[0][jj];
+            forceData[particle2Index][jj] -= forces[1][jj];
+            forceData[particle3Index][jj] -= forces[2][jj];
+        }
+
     }   
     return energy;
 }

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

@@ -185,21 +185,21 @@ RealOpenMM AmoebaReferenceTorsionForce::calculateTorsionIxn( const RealVec& posi
      
     // forces
   
-    forces[0][0] -= tempVector[0][0];
-    forces[0][1] -= tempVector[0][1];
-    forces[0][2] -= tempVector[0][2];
+    forces[0][0] = tempVector[0][0];
+    forces[0][1] = tempVector[0][1];
+    forces[0][2] = tempVector[0][2];
  
-    forces[1][0] -= tempVector[1][0];
-    forces[1][1] -= tempVector[1][1];
-    forces[1][2] -= tempVector[1][2];
+    forces[1][0] = tempVector[1][0];
+    forces[1][1] = tempVector[1][1];
+    forces[1][2] = tempVector[1][2];
  
-    forces[2][0] -= tempVector[2][0];
-    forces[2][1] -= tempVector[2][1];
-    forces[2][2] -= tempVector[2][2];
+    forces[2][0] = tempVector[2][0];
+    forces[2][1] = tempVector[2][1];
+    forces[2][2] = tempVector[2][2];
  
-    forces[3][0] -= tempVector[3][0];
-    forces[3][1] -= tempVector[3][1];
-    forces[3][2] -= tempVector[3][2];
+    forces[3][0] = tempVector[3][0];
+    forces[3][1] = tempVector[3][1];
+    forces[3][2] = tempVector[3][2];
  
     // ---------------------------------------------------------------------------------------
  
@@ -228,12 +228,17 @@ RealOpenMM AmoebaReferenceTorsionForce::calculateForceAndEnergy( int numTorsions
         int particle3Index      = particle3[ii];
         int particle4Index      = particle4[ii];
         RealVec forces[4];
-        forces[0]               = forceData[particle1Index];
-        forces[1]               = forceData[particle2Index];
-        forces[2]               = forceData[particle3Index];
-        forces[3]               = forceData[particle4Index];
         energy                 += calculateTorsionIxn( posData[particle1Index], posData[particle2Index], posData[particle3Index], posData[particle4Index],
                                                        torsionParameters1[ii], torsionParameters2[ii], torsionParameters3[ii], forces );
+        // accumulate forces
+    
+        for( int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] -= forces[0][jj];
+            forceData[particle2Index][jj] -= forces[1][jj];
+            forceData[particle3Index][jj] -= forces[2][jj];
+            forceData[particle4Index][jj] -= forces[3][jj];
+        }
+
     }   
     return energy;
 }

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

@@ -542,9 +542,9 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateTorsionTorsionIxn( const
     // add in forces
  
     for( int jj = 0; jj < LastAtomIndex; jj++ ){
-        forces[jj][0] -= d[jj][0];
-        forces[jj][1] -= d[jj][1];
-        forces[jj][2] -= d[jj][2];
+        forces[jj][0] = d[jj][0];
+        forces[jj][1] = d[jj][1];
+        forces[jj][2] = d[jj][2];
     }
  
     // ---------------------------------------------------------------------------------------
@@ -576,12 +576,6 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy( int numT
         int gridIndex            = gridIndices[ii];
 
         RealVec forces[5];
-        forces[0]                = forceData[particle1Index];
-        forces[1]                = forceData[particle2Index];
-        forces[2]                = forceData[particle3Index];
-        forces[3]                = forceData[particle4Index];
-        forces[4]                = forceData[particle5Index];
-
         RealVec* chiralCheckAtom;
         if( chiralCheckAtomIndex >= 0 ){
             chiralCheckAtom = &posData[chiralCheckAtomIndex];
@@ -593,6 +587,16 @@ RealOpenMM AmoebaReferenceTorsionTorsionForce::calculateForceAndEnergy( int numT
                                                               posData[particle5Index], chiralCheckAtom, torsionTorsionGrids[gridIndex],
                                                               forces );
 
+        // accumulate forces
+    
+        for( int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] -= forces[0][jj];
+            forceData[particle2Index][jj] -= forces[1][jj];
+            forceData[particle3Index][jj] -= forces[2][jj];
+            forceData[particle4Index][jj] -= forces[3][jj];
+            forceData[particle5Index][jj] -= forces[4][jj];
+        }
+
     }   
     return energy;
 }

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

@@ -75,13 +75,13 @@ RealOpenMM AmoebaReferenceUreyBradleyForce::calculateUreyBradleyIxn( const RealV
    dEdR                      *= two*kForceConstant*deltaIdeal;
    dEdR                       = r > zero ? (dEdR/r) : zero;
 
-   forces[0][0]              += dEdR*deltaR[0];
-   forces[0][1]              += dEdR*deltaR[1];
-   forces[0][2]              += dEdR*deltaR[2];
+   forces[0][0]              = dEdR*deltaR[0];
+   forces[0][1]              = dEdR*deltaR[1];
+   forces[0][2]              = dEdR*deltaR[2];
 
-   forces[1][0]              -= dEdR*deltaR[0];
-   forces[1][1]              -= dEdR*deltaR[1];
-   forces[1][2]              -= dEdR*deltaR[2];
+   forces[1][0]              = dEdR*deltaR[0];
+   forces[1][1]              = dEdR*deltaR[1];
+   forces[1][2]              = dEdR*deltaR[2];
 
    RealOpenMM energy          = kForceConstant*deltaIdeal2*( one + cubicK*deltaIdeal + quarticK*deltaIdeal2 );
    return energy;
@@ -103,11 +103,16 @@ RealOpenMM AmoebaReferenceUreyBradleyForce::calculateForceAndEnergy( int numIxns
         RealOpenMM idealLength        = length[ii];
         RealOpenMM kForceConstant     = kQuadratic[ii];
         RealVec forces[2];
-        forces[0]                     = forceData[particle1Index];
-        forces[1]                     = forceData[particle2Index];
         energy                       += calculateUreyBradleyIxn( particlePositions[particle1Index], particlePositions[particle2Index],
                                                                  idealLength, kForceConstant, globalUreyBradleyCubic, globalUreyBradleyQuartic,
                                                                  forces );
+        // accumulate forces
+    
+        for( int jj = 0; jj < 3; jj++ ){
+            forceData[particle1Index][jj] += forces[0][jj];
+            forceData[particle2Index][jj] -= forces[1][jj];
+        }
+
     }   
     return energy;
 }