Added calculateTotalDipoles function

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

@@ -331,6 +331,14 @@ public:
      */
     void getInducedDipoles(Context& context, std::vector<Vec3>& dipoles);
 
+    /**
+     * Get the total dipole moments of all particles.
+     *
+     * @param context         the Context for which to get the induced dipoles
+     * @param[out] dipoles    the induced dipole moment of particle i is stored into the i'th element
+     */
+    void getTotalDipoles(Context& context, std::vector<Vec3>& dipoles);
+
     /**
      * Get the electrostatic potential.
      *

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

@@ -350,6 +350,7 @@ public:
 
     virtual void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles) = 0;
     virtual void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles) = 0;
+    virtual void getTotalDipoles(ContextImpl& context, std::vector<Vec3>& dipoles) = 0;
 
     virtual void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                            std::vector< double >& outputElectrostaticPotential) = 0;

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

@@ -83,6 +83,7 @@ public:
     static void getCovalentDegree(const AmoebaMultipoleForce& force, std::vector<int>& covalentDegree);
     void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
     void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
+    void getTotalDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
 
     void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                    std::vector< double >& outputElectrostaticPotential);

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -238,6 +238,10 @@ void AmoebaMultipoleForce::getLabFramePermanentDipoles(Context& context, vector<
     dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getLabFramePermanentDipoles(getContextImpl(context), dipoles);
 }
 
+void AmoebaMultipoleForce::getTotalDipoles(Context& context, vector<Vec3>& dipoles) {
+    dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getTotalDipoles(getContextImpl(context), dipoles);
+}
+
 void AmoebaMultipoleForce::getElectrostaticPotential(const std::vector< Vec3 >& inputGrid, Context& context, std::vector< double >& outputElectrostaticPotential) {
     dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getElectrostaticPotential(getContextImpl(context), inputGrid, outputElectrostaticPotential);
 }

plugins/amoeba/openmmapi/src/AmoebaMultipoleForceImpl.cpp

@@ -191,6 +191,10 @@ void AmoebaMultipoleForceImpl::getInducedDipoles(ContextImpl& context, vector<Ve
     kernel.getAs<CalcAmoebaMultipoleForceKernel>().getInducedDipoles(context, dipoles);
 }
 
+void AmoebaMultipoleForceImpl::getTotalDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
+    kernel.getAs<CalcAmoebaMultipoleForceKernel>().getTotalDipoles(context, dipoles);
+}
+
 void AmoebaMultipoleForceImpl::getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                                          std::vector< double >& outputElectrostaticPotential) {
     kernel.getAs<CalcAmoebaMultipoleForceKernel>().getElectrostaticPotential(context, inputGrid, outputElectrostaticPotential);

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

@@ -1812,6 +1812,48 @@ void CudaCalcAmoebaMultipoleForceKernel::getInducedDipoles(ContextImpl& context,
     }
 }
 
+
+void CudaCalcAmoebaMultipoleForceKernel::getTotalDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
+    ensureMultipolesValid(context);
+    int numParticles = cu.getNumAtoms();
+    dipoles.resize(numParticles);
+    const vector<int>& order = cu.getAtomIndex();
+    if (cu.getUseDoublePrecision()) {
+        vector<double> posqVec;
+        vector<double> labDipoleVec;
+        vector<double> inducedDipoleVec;
+        vector<double> totalDipoleVecX(numParticles);
+        vector<double> totalDipoleVecY(numParticles);
+        vector<double> totalDipoleVecZ(numParticles);
+        inducedDipole->download(inducedDipoleVec);
+        labFrameDipoles->download(labDipoleVec);
+        cu.getPosq().download(posqVec);
+        for (int i = 0; i < numParticles; i++) {
+            totalDipoleVecX[i] = posqVec[i].x * posqVec[i].w + labDipoleVec[3*i] + inducedDipoleVec[3*i];
+            totalDipoleVecY[i] = posqVec[i].y * posqVec[i].w + labDipoleVec[3*i+1] + inducedDipoleVec[3*i+1];
+            totalDipoleVecZ[i] = posqVec[i].z * posqVec[i].w + labDipoleVec[3*i+2] + inducedDipoleVec[3*i+2];
+            dipoles[order[i]] = Vec3(totalDipoleVecX[i], totalDipoleVecY[i], totalDipoleVecZ[i]);
+        }
+    }
+    else {
+        vector<float> posqVec;
+        vector<float> labFramePermanentDipolesVec;
+        vector<float> inducedDipolesVec;
+        vector<float> totalDipoleVecX(numParticles);
+        vector<float> totalDipoleVecY(numParticles);
+        vector<float> totalDipoleVecZ(numParticles);
+        inducedDipole->download(inducedDipoleVec);
+        labFrameDipoles->download(labDipoleVec);
+        cu.getPosq().download(posqVec);
+        for (int i = 0; i < numParticles; i++) {
+            totalDipoleVecX[i] = posqVec[i].x * posqVec[i].w + labDipoleVec[3*i] + inducedDipoleVec[3*i];
+            totalDipoleVecY[i] = posqVec[i].y * posqVec[i].w + labDipoleVec[3*i+1] + inducedDipoleVec[3*i+1];
+            totalDipoleVecZ[i] = posqVec[i].z * posqVec[i].w + labDipoleVec[3*i+2] + inducedDipoleVec[3*i+2];
+            dipoles[order[i]] = Vec3(totalDipoleVecX[i], totalDipoleVecY[i], totalDipoleVecZ[i]);
+        }
+    }
+}
+
 void CudaCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl& context, const vector<Vec3>& inputGrid, vector<double>& outputElectrostaticPotential) {
     ensureMultipolesValid(context);
     int numPoints = inputGrid.size();

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

@@ -342,6 +342,13 @@ public:
      * @param dipoles    the induced dipole moment of particle i is stored into the i'th element
      */
     void getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
+    /**
+     * Get the total dipole moments of all particles.
+     * 
+     * @param context    the Context for which to get the induced dipoles
+     * @param dipoles    the induced dipole moment of particle i is stored into the i'th element
+     */
+    void getTotalDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
     /**
      * Execute the kernel to calculate the electrostatic potential
      *

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

@@ -730,6 +730,27 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getLabFramePermanentDipoles(Contex
 }
 
 
+void ReferenceCalcAmoebaMultipoleForceKernel::getTotalDipoles(ContextImpl& context, vector<Vec3>& outputDipoles) {
+    int numParticles = context.getSystem().getNumParticles();
+    outputDipoles.resize(numParticles);
+
+    // Create an AmoebaReferenceMultipoleForce to do the calculation.
+    
+    AmoebaReferenceMultipoleForce* amoebaReferenceMultipoleForce = setupAmoebaReferenceMultipoleForce(context);
+    vector<RealVec>& posData = extractPositions(context);
+    
+    // Retrieve the permanent dipoles in the lab frame.
+    
+    vector<RealVec> totalDipoles;
+    amoebaReferenceMultipoleForce->calculateTotalDipoles(posData, charges, dipoles, quadrupoles, tholes,
+            dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, inducedDipoles);
+
+    for (int i = 0; i < numParticles; i++)
+        outputDipoles[i] = totalDipoles[i];
+    delete amoebaReferenceMultipoleForce;
+}
+
+
 
 void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                                                         std::vector< double >& outputElectrostaticPotential) {

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

@@ -381,6 +381,13 @@ public:
      * @param dipoles    the fixed dipole moment of particle i is stored into the i'th element
      */
     void getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
+    /**
+     * Get the total dipole moments of all particles in the global reference frame.
+     * 
+     * @param context    the Context for which to get the fixed dipoles
+     * @param dipoles    the fixed dipole moment of particle i is stored into the i'th element
+     */
+    void getTotalDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
     /** 
      * Calculate the electrostatic potential given vector of grid coordinates.
      *

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

@@ -1819,6 +1819,8 @@ void AmoebaReferenceMultipoleForce::calculateInducedDipoles(const vector<RealVec
 }
 
 
+
+
 void AmoebaReferenceMultipoleForce::calculateLabFramePermanentDipoles(const vector<RealVec>& particlePositions,
                                                                       const vector<RealOpenMM>& charges,
                                                                       const vector<RealOpenMM>& dipoles,
@@ -1845,6 +1847,35 @@ void AmoebaReferenceMultipoleForce::calculateLabFramePermanentDipoles(const vect
       }
 }
 
+void AmoebaReferenceMultipoleForce::calculateTotalDipoles(const vector<RealVec>& particlePositions,
+                                                          const vector<RealOpenMM>& charges,
+                                                          const vector<RealOpenMM>& dipoles,
+                                                          const vector<RealOpenMM>& quadrupoles,
+                                                          const vector<RealOpenMM>& tholes,
+                                                          const vector<RealOpenMM>& dampingFactors,
+                                                          const vector<RealOpenMM>& polarity,
+                                                          const vector<int>& axisTypes,
+                                                          const vector<int>& multipoleAtomZs,
+                                                          const vector<int>& multipoleAtomXs,
+                                                          const vector<int>& multipoleAtomYs,
+                                                          const vector< vector< vector<int> > >& multipoleAtomCovalentInfo,
+                                                          vector<RealVec>& outputTotalDipoles) {
+    // setup, including calculating permanent dipoles
+
+    vector<MultipoleParticleData> particleData;
+    setup(particlePositions, charges, dipoles, quadrupoles, tholes,
+           dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
+           multipoleAtomCovalentInfo, particleData);
+    outputTotalDipoles.resize(_numParticles);
+    for (int i = 0; i < _numParticles; i++)
+      {
+      for (int j = 0; j < 3; j++)
+        {
+        outputTotalDipoles[i][j] = particleData[i].position[j] * particleData[i].charge + particleData[i].dipole[j] + _inducedDipole[j];
+        }
+      }
+}
+
 void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments(const vector<RealOpenMM>& masses,
                                                                           const vector<RealVec>& particlePositions,
                                                                           const vector<RealOpenMM>& charges,

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

@@ -519,7 +519,7 @@ public:
                                  const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
                                  std::vector<RealVec>& outputInducedDipoles);
 
-    /*
+    /**
      * Calculate particle permanent dipoles rotated in the lab frame.
      *
      * @param masses                    particle masses
@@ -538,19 +538,55 @@ public:
      * @param outputMultipoleMoments    output multipole moments
      */
 
-    void calculateLabFramePermanentDipoles(const vector<RealVec>& particlePositions,
-                                           const vector<RealOpenMM>& charges,
-                                           const vector<RealOpenMM>& dipoles,
-                                           const vector<RealOpenMM>& quadrupoles,
-                                           const vector<RealOpenMM>& tholes,
-                                           const vector<RealOpenMM>& dampingFactors,
-                                           const vector<RealOpenMM>& polarity,
-                                           const vector<int>& axisTypes,
-                                           const vector<int>& multipoleAtomZs,
-                                           const vector<int>& multipoleAtomXs,
-                                           const vector<int>& multipoleAtomYs,
-                                           const vector< vector< vector<int> > >& multipoleAtomCovalentInfo,
-                                           vector<RealVec>& outputRotatedPermanentDipoles);
+    void calculateLabFramePermanentDipoles(const std::vector<RealVec>& particlePositions,
+                                           const std::vector<RealOpenMM>& charges,
+                                           const std::vector<RealOpenMM>& dipoles,
+                                           const std::vector<RealOpenMM>& quadrupoles,
+                                           const std::vector<RealOpenMM>& tholes,
+                                           const std::vector<RealOpenMM>& dampingFactors,
+                                           const std::vector<RealOpenMM>& polarity,
+                                           const std::vector<int>& axisTypes,
+                                           const std::vector<int>& multipoleAtomZs,
+                                           const std::vector<int>& multipoleAtomXs,
+                                           const std::vector<int>& multipoleAtomYs,
+                                           const std::vector< vector< vector<int> > >& multipoleAtomCovalentInfo,
+                                           std::vector<RealVec>& outputRotatedPermanentDipoles);
+
+    /**
+     * Calculate particle total dipoles.
+     *
+     * @param masses                    particle masses
+     * @param particlePositions         Cartesian coordinates of particles
+     * @param charges                   scalar charges for each particle
+     * @param dipoles                   molecular frame dipoles for each particle
+     * @param quadrupoles               molecular frame quadrupoles for each particle
+     * @param tholes                    Thole factors for each particle
+     * @param dampingFactors            dampling factors for each particle
+     * @param polarity                  polarity for each particle
+     * @param axisTypes                 axis type (Z-then-X, ...) for each particle
+     * @param multipoleAtomZs           indicies of particle specifying the molecular frame z-axis for each particle
+     * @param multipoleAtomXs           indicies of particle specifying the molecular frame x-axis for each particle
+     * @param multipoleAtomYs           indicies of particle specifying the molecular frame y-axis for each particle
+     * @param multipoleAtomCovalentInfo covalent info needed to set scaling factors
+     * @param outputMultipoleMoments    output multipole moments
+     */
+
+
+    void calculateTotalDipoles(const std::vector<RealVec>& particlePositions,
+                                           const std::vector<RealOpenMM>& charges,
+                                           const std::vector<RealOpenMM>& dipoles,
+                                           const std::vector<RealOpenMM>& quadrupoles,
+                                           const std::vector<RealOpenMM>& tholes,
+                                           const std::vector<RealOpenMM>& dampingFactors,
+                                           const std::vector<RealOpenMM>& polarity,
+                                           const std::vector<int>& axisTypes,
+                                           const std::vector<int>& multipoleAtomZs,
+                                           const std::vector<int>& multipoleAtomXs,
+                                           const std::vector<int>& multipoleAtomYs,
+                                           const std::vector< vector< vector<int> > >& multipoleAtomCovalentInfo,
+                                           std::vector<RealVec>& outputRotatedPermanentDipoles);
+
+
 
     /**
      * Calculate system multipole moments.
@@ -641,6 +677,10 @@ protected:
             RealOpenMM thole;
             RealOpenMM dampingFactor;
             RealOpenMM polarity;
+
+            std::vector<RealOpenMM> xDipoleTraj;
+            std::vector<RealOpenMM> yDipoleTraj;
+            std::vector<RealOpenMM> zDipoleTraj;
     };
     
     /**

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -125,6 +125,7 @@ NO_OUTPUT_ARGS = [('LocalEnergyMinimizer', 'minimize', 'context'),
                   ('AmoebaMultipoleForce', 'getElectrostaticPotential', 'context'),
                   ('AmoebaMultipoleForce', 'getInducedDipoles', 'context'),
                   ('AmoebaMultipoleForce', 'getLabFramePermanentDipoles', 'context'),
+                  ('AmoebaMultipoleForce', 'getTotalDipoles', 'context'),
 ]
 
 # SWIG assumes the target language shadow class owns the C++ class
@@ -275,6 +276,7 @@ UNITS = {
 ("AmoebaMultipoleForce",                 "getElectrostaticPotential")                     :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getInducedDipoles")                             :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getLabFramePermanentDipoles")                   :  ( None, ()),
+("AmoebaMultipoleForce",                 "getTotalDipoles")                               :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getSystemMultipoleMoments")                     :  ( None, ()),
 
 ("AmoebaOutOfPlaneBendForce",            "getNumOutOfPlaneBends")                         :  ( None, ()),