Modified the AmoebaMultipoleForce.getSystemMultipoleMoments method with two...

Modified the AmoebaMultipoleForce.getSystemMultipoleMoments method with two extra arguments; toggle whether to evaluate energy/force and select maximum multipole order.

plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.h

@@ -319,16 +319,21 @@ public:
      * <i>lowest nonvanishing moment</i> has a well defined value.  This means that if the system has a net
      * nonzero charge, the dipole and quadrupole moments are not well defined and should be ignored.  If the
      * net charge is zero, the dipole moment is well defined (and really represents a dipole density), but
-     * the quadrupole moment is still undefined and should be ignored.
+     * the quadrupole moment is still undefined and should be ignored.  Specify "compute = True" if you want
+     * to do an energy and force calculation prior to getting the moments; if compute = False, it will be faster
+     * but we will lose accuracy.  Specify "order" for the maximum order computed (order = 2 means we go up 
+     * through quadrupoles.)
      *
      * @param context      context
+     * @param compute      Whether to compute energies and forces first
+     * @param order        Maximum order to go through, above 2 for quadrupole has no effect
      * @param outputMultipoleMonents (charge,
                                       dipole_x, dipole_y, dipole_z,
                                       quadrupole_xx, quadrupole_xy, quadrupole_xz,
                                       quadrupole_yx, quadrupole_yy, quadrupole_yz,
                                       quadrupole_zx, quadrupole_zy, quadrupole_zz)
      */
-    void getSystemMultipoleMoments(Context& context, std::vector< double >& outputMultipoleMoments);
+    void getSystemMultipoleMoments(Context& context, bool compute, int order, std::vector< double >& outputMultipoleMoments);
     /**
      * Update the multipole parameters in a Context to match those stored in this Force object.  This method
      * provides an efficient method to update certain parameters in an existing Context without needing to reinitialize it.

plugins/amoeba/openmmapi/include/openmm/amoebaKernels.h

@@ -351,7 +351,7 @@ public:
     virtual void getElectrostaticPotential( ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                             std::vector< double >& outputElectrostaticPotential ) = 0;
 
-    virtual void getSystemMultipoleMoments( ContextImpl& context, std::vector< double >& outputMultipoleMonents ) = 0;
+    virtual void getSystemMultipoleMoments( ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMonents ) = 0;
     /**
      * Copy changed parameters over to a context.
      *

plugins/amoeba/openmmapi/include/openmm/internal/AmoebaMultipoleForceImpl.h

@@ -85,7 +85,7 @@ public:
     void getElectrostaticPotential( ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                     std::vector< double >& outputElectrostaticPotential );
 
-    void getSystemMultipoleMoments( ContextImpl& context, std::vector< double >& outputMultipoleMonents );
+    void getSystemMultipoleMoments( ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMonents );
     void updateParametersInContext(ContextImpl& context);
  
 

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -230,8 +230,8 @@ void AmoebaMultipoleForce::getElectrostaticPotential( const std::vector< Vec3 >&
     dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getElectrostaticPotential(getContextImpl(context), inputGrid, outputElectrostaticPotential);
 }
 
-void AmoebaMultipoleForce::getSystemMultipoleMoments(Context& context, std::vector< double >& outputMultipoleMonents ){
-    dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getSystemMultipoleMoments(getContextImpl(context), outputMultipoleMonents);
+void AmoebaMultipoleForce::getSystemMultipoleMoments(Context& context, bool compute, int order, std::vector< double >& outputMultipoleMonents ){
+    dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getSystemMultipoleMoments(getContextImpl(context), compute, order, outputMultipoleMonents);
 }
 
 ForceImpl* AmoebaMultipoleForce::createImpl() {

plugins/amoeba/openmmapi/src/AmoebaMultipoleForceImpl.cpp

@@ -188,8 +188,8 @@ void AmoebaMultipoleForceImpl::getElectrostaticPotential( ContextImpl& context,
     kernel.getAs<CalcAmoebaMultipoleForceKernel>().getElectrostaticPotential(context, inputGrid, outputElectrostaticPotential);
 }
 
-void AmoebaMultipoleForceImpl::getSystemMultipoleMoments( ContextImpl& context, std::vector< double >& outputMultipoleMonents ){
-    kernel.getAs<CalcAmoebaMultipoleForceKernel>().getSystemMultipoleMoments(context, outputMultipoleMonents);
+void AmoebaMultipoleForceImpl::getSystemMultipoleMoments( ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMonents ){
+    kernel.getAs<CalcAmoebaMultipoleForceKernel>().getSystemMultipoleMoments(context, compute, order, outputMultipoleMonents);
 }
 
 void AmoebaMultipoleForceImpl::updateParametersInContext(ContextImpl& context) {

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

@@ -1596,7 +1596,7 @@ void CudaCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl&
 }
 
 template <class T, class T4, class M4>
-void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextImpl& context, vector<double>& outputMultipoleMoments) {
+void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextImpl& context, int order, vector<double>& outputMultipoleMoments) {
     // Compute the local coordinates relative to the center of mass.
     int numAtoms = cu.getNumAtoms();
     vector<T4> posq;
@@ -1641,56 +1641,66 @@ void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextIm
     double zyqdp = 0.0;
     double zzqdp = 0.0;
     vector<T> labDipoleVec, inducedDipoleVec, quadrupoleVec;
-    labFrameDipoles->download(labDipoleVec);
-    inducedDipole->download(inducedDipoleVec);
-    labFrameQuadrupoles->download(quadrupoleVec);
+    if (order >= 1) {
+        labFrameDipoles->download(labDipoleVec);
+        inducedDipole->download(inducedDipoleVec);
+    }
+    if (order >= 2) {
+        labFrameQuadrupoles->download(quadrupoleVec);
+    }
     for (int i = 0; i < numAtoms; i++) {
         totalCharge += posqLocal[i].w;
-        double netDipoleX = (labDipoleVec[3*i]    + inducedDipoleVec[3*i]);
-        double netDipoleY = (labDipoleVec[3*i+1]  + inducedDipoleVec[3*i+1]);
-        double netDipoleZ = (labDipoleVec[3*i+2]  + inducedDipoleVec[3*i+2]);
-        xdpl += posqLocal[i].x*posqLocal[i].w + netDipoleX;
-        ydpl += posqLocal[i].y*posqLocal[i].w + netDipoleY;
-        zdpl += posqLocal[i].z*posqLocal[i].w + netDipoleZ;
-        xxqdp += posqLocal[i].x*posqLocal[i].x*posqLocal[i].w + 2*posqLocal[i].x*netDipoleX;
-        xyqdp += posqLocal[i].x*posqLocal[i].y*posqLocal[i].w + posqLocal[i].x*netDipoleY + posqLocal[i].y*netDipoleX;
-        xzqdp += posqLocal[i].x*posqLocal[i].z*posqLocal[i].w + posqLocal[i].x*netDipoleZ + posqLocal[i].z*netDipoleX;
-        yxqdp += posqLocal[i].y*posqLocal[i].x*posqLocal[i].w + posqLocal[i].y*netDipoleX + posqLocal[i].x*netDipoleY;
-        yyqdp += posqLocal[i].y*posqLocal[i].y*posqLocal[i].w + 2*posqLocal[i].y*netDipoleY;
-        yzqdp += posqLocal[i].y*posqLocal[i].z*posqLocal[i].w + posqLocal[i].y*netDipoleZ + posqLocal[i].z*netDipoleY;
-        zxqdp += posqLocal[i].z*posqLocal[i].x*posqLocal[i].w + posqLocal[i].z*netDipoleX + posqLocal[i].x*netDipoleZ;
-        zyqdp += posqLocal[i].z*posqLocal[i].y*posqLocal[i].w + posqLocal[i].z*netDipoleY + posqLocal[i].y*netDipoleZ;
-        zzqdp += posqLocal[i].z*posqLocal[i].z*posqLocal[i].w + 2*posqLocal[i].z*netDipoleZ;
+	double netDipoleX = (labDipoleVec[3*i]    + inducedDipoleVec[3*i]);
+	double netDipoleY = (labDipoleVec[3*i+1]  + inducedDipoleVec[3*i+1]);
+	double netDipoleZ = (labDipoleVec[3*i+2]  + inducedDipoleVec[3*i+2]);
+	if (order >= 1) {
+            xdpl += posqLocal[i].x*posqLocal[i].w + netDipoleX;
+            ydpl += posqLocal[i].y*posqLocal[i].w + netDipoleY;
+            zdpl += posqLocal[i].z*posqLocal[i].w + netDipoleZ;
+	}
+	if (order >= 2) {
+            xxqdp += posqLocal[i].x*posqLocal[i].x*posqLocal[i].w + 2*posqLocal[i].x*netDipoleX;
+            xyqdp += posqLocal[i].x*posqLocal[i].y*posqLocal[i].w + posqLocal[i].x*netDipoleY + posqLocal[i].y*netDipoleX;
+            xzqdp += posqLocal[i].x*posqLocal[i].z*posqLocal[i].w + posqLocal[i].x*netDipoleZ + posqLocal[i].z*netDipoleX;
+            yxqdp += posqLocal[i].y*posqLocal[i].x*posqLocal[i].w + posqLocal[i].y*netDipoleX + posqLocal[i].x*netDipoleY;
+            yyqdp += posqLocal[i].y*posqLocal[i].y*posqLocal[i].w + 2*posqLocal[i].y*netDipoleY;
+            yzqdp += posqLocal[i].y*posqLocal[i].z*posqLocal[i].w + posqLocal[i].y*netDipoleZ + posqLocal[i].z*netDipoleY;
+            zxqdp += posqLocal[i].z*posqLocal[i].x*posqLocal[i].w + posqLocal[i].z*netDipoleX + posqLocal[i].x*netDipoleZ;
+            zyqdp += posqLocal[i].z*posqLocal[i].y*posqLocal[i].w + posqLocal[i].z*netDipoleY + posqLocal[i].y*netDipoleZ;
+            zzqdp += posqLocal[i].z*posqLocal[i].z*posqLocal[i].w + 2*posqLocal[i].z*netDipoleZ;
+	}
     }
 
     // Convert the quadrupole from traced to traceless form.
  
-    double qave = (xxqdp + yyqdp + zzqdp)/3;
-    xxqdp = 1.5*(xxqdp-qave);
-    xyqdp = 1.5*xyqdp;
-    xzqdp = 1.5*xzqdp;
-    yxqdp = 1.5*yxqdp;
-    yyqdp = 1.5*(yyqdp-qave);
-    yzqdp = 1.5*yzqdp;
-    zxqdp = 1.5*zxqdp;
-    zyqdp = 1.5*zyqdp;
-    zzqdp = 1.5*(zzqdp-qave);
+    if (order >= 2) {
+        double qave = (xxqdp + yyqdp + zzqdp)/3;
+        xxqdp = 1.5*(xxqdp-qave);
+        xyqdp = 1.5*xyqdp;
+        xzqdp = 1.5*xzqdp;
+        yxqdp = 1.5*yxqdp;
+        yyqdp = 1.5*(yyqdp-qave);
+        yzqdp = 1.5*yzqdp;
+        zxqdp = 1.5*zxqdp;
+        zyqdp = 1.5*zyqdp;
+        zzqdp = 1.5*(zzqdp-qave);
 
     // Add the traceless atomic quadrupoles to the total quadrupole moment.
 
-    for (int i = 0; i < numAtoms; i++) {
-        xxqdp = xxqdp + 3*quadrupoleVec[5*i];
-        xyqdp = xyqdp + 3*quadrupoleVec[5*i+1];
-        xzqdp = xzqdp + 3*quadrupoleVec[5*i+2];
-        yxqdp = yxqdp + 3*quadrupoleVec[5*i+1];
-        yyqdp = yyqdp + 3*quadrupoleVec[5*i+3];
-        yzqdp = yzqdp + 3*quadrupoleVec[5*i+4];
-        zxqdp = zxqdp + 3*quadrupoleVec[5*i+2];
-        zyqdp = zyqdp + 3*quadrupoleVec[5*i+4];
-        zzqdp = zzqdp + -3*(quadrupoleVec[5*i]+quadrupoleVec[5*i+3]);
-    }
- 
-    double debye = 4.80321;
+        for (int i = 0; i < numAtoms; i++) {
+            xxqdp = xxqdp + 3*quadrupoleVec[5*i];
+            xyqdp = xyqdp + 3*quadrupoleVec[5*i+1];
+            xzqdp = xzqdp + 3*quadrupoleVec[5*i+2];
+            yxqdp = yxqdp + 3*quadrupoleVec[5*i+1];
+            yyqdp = yyqdp + 3*quadrupoleVec[5*i+3];
+            yzqdp = yzqdp + 3*quadrupoleVec[5*i+4];
+            zxqdp = zxqdp + 3*quadrupoleVec[5*i+2];
+            zyqdp = zyqdp + 3*quadrupoleVec[5*i+4];
+            zzqdp = zzqdp + -3*(quadrupoleVec[5*i]+quadrupoleVec[5*i+3]);
+        }
+    } 
+
+    double debye = 48.0321; // If we use 4.80321 then the result is ten times too small.
     outputMultipoleMoments.resize(13);
     outputMultipoleMoments[0] = totalCharge;
     outputMultipoleMoments[1] = xdpl*debye;
@@ -1707,14 +1717,15 @@ void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextIm
     outputMultipoleMoments[12] = zzqdp*debye;
 }
 
-void CudaCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, vector<double>& outputMultipoleMoments) {
-    context.calcForcesAndEnergy(false, false, -1);
+void CudaCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, bool compute, int order, vector<double>& outputMultipoleMoments) {
+    if (compute)
+        context.calcForcesAndEnergy(false, false, -1);
     if (cu.getUseDoublePrecision())
-        computeSystemMultipoleMoments<double, double4, double4>(context, outputMultipoleMoments);
+        computeSystemMultipoleMoments<double, double4, double4>(context, order, outputMultipoleMoments);
     else if (cu.getUseMixedPrecision())
-        computeSystemMultipoleMoments<float, float4, double4>(context, outputMultipoleMoments);
+        computeSystemMultipoleMoments<float, float4, double4>(context, order, outputMultipoleMoments);
     else
-        computeSystemMultipoleMoments<float, float4, float4>(context, outputMultipoleMoments);
+        computeSystemMultipoleMoments<float, float4, float4>(context, order, outputMultipoleMoments);
 }
 
 void CudaCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force) {

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

@@ -341,13 +341,15 @@ public:
      * Get the system multipole moments
      *
      * @param context      context
+     * @param compute      Whether to compute energies and forces first
+     * @param order        Maximum order to go through, above 2 for quadrupole has no effect
      * @param outputMultipoleMoments (charge,
      *                                dipole_x, dipole_y, dipole_z,
      *                                quadrupole_xx, quadrupole_xy, quadrupole_xz,
      *                                quadrupole_yx, quadrupole_yy, quadrupole_yz,
      *                                quadrupole_zx, quadrupole_zy, quadrupole_zz )
      */
-    void getSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
+    void getSystemMultipoleMoments(ContextImpl& context, bool compute, int order, std::vector<double>& outputMultipoleMoments);
     /**
      * Copy changed parameters over to a context.
      *
@@ -368,7 +370,7 @@ private:
         const char* getSortKey() const {return "value.y";}
     };
     void initializeScaleFactors();
-    template <class T, class T4, class M4> void computeSystemMultipoleMoments(ContextImpl& context, std::vector<double>& outputMultipoleMoments);
+    template <class T, class T4, class M4> void computeSystemMultipoleMoments(ContextImpl& context, int order, std::vector<double>& outputMultipoleMoments);
     int numMultipoles, maxInducedIterations;
     int fixedFieldThreads, inducedFieldThreads, electrostaticsThreads;
     double inducedEpsilon;

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

@@ -2009,7 +2009,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::No
     context.setPositions(positions);
 
     if( testName == "testSystemMultipoleMoments" ){
-        amoebaMultipoleForce->getSystemMultipoleMoments(context, outputMultipoleMoments );
+        amoebaMultipoleForce->getSystemMultipoleMoments(context, 1, 2, outputMultipoleMoments );
     } else if( testName == "testMultipoleGridPotential" ){
         amoebaMultipoleForce->getElectrostaticPotential( inputGrid, context, outputGridPotential );
     } else {

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

@@ -708,7 +708,7 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextI
     return;
 }
 
-void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments){
+void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMoments){
 
     // retrieve masses
 

plugins/amoeba/platforms/reference/src/AmoebaReferenceKernels.h

@@ -380,6 +380,8 @@ public:
      * Get the system multipole moments.
      *
      * @param context                context 
+     * @param compute                Whether to compute energies and forces first
+     * @param order                  Maximum order to go through, above 2 for quadrupole has no effect
      * @param outputMultipoleMonents vector of multipole moments:
                                      (charge,
                                       dipole_x, dipole_y, dipole_z,
@@ -387,7 +389,7 @@ public:
                                       quadrupole_yx, quadrupole_yy, quadrupole_yz,
                                       quadrupole_zx, quadrupole_zy, quadrupole_zz )
      */
-    void getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments);
+    void getSystemMultipoleMoments(ContextImpl& context, bool compute, int order, std::vector< double >& outputMultipoleMoments);
     /**
      * Copy changed parameters over to a context.
      *

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

@@ -1915,7 +1915,7 @@ static void setupAndGetForcesEnergyMultipoleLargeWater( AmoebaMultipoleForce::No
     context.setPositions(positions);
 
     if( testName == "testSystemMultipoleMoments" ){
-        amoebaMultipoleForce->getSystemMultipoleMoments( context, outputMultipoleMoments );
+        amoebaMultipoleForce->getSystemMultipoleMoments( context, 1, 2, outputMultipoleMoments );
     } else if( testName == "testMultipoleGridPotential" ){
         amoebaMultipoleForce->getElectrostaticPotential( inputGrid, context, outputGridPotential );
     } else {