Merge pull request #1534 from saurabhbelsare/master

Added functions to get per atom lab frame permanent dipoles and per atom total dipole moments for AMOEBA

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

@@ -370,7 +370,13 @@ public:
      * This can be overridden by explicitly setting an alpha parameter and grid dimensions to use.
      */
     void setEwaldErrorTolerance(double tol);
-
+    /**
+     * Get the fixed dipole moments of all particles in the global reference frame.
+     *
+     * @param context         the Context for which to get the fixed dipoles
+     * @param[out] dipoles    the fixed dipole moment of particle i is stored into the i'th element
+     */
+    void getLabFramePermanentDipoles(Context& context, std::vector<Vec3>& dipoles);
     /**
      * Get the induced dipole moments of all particles.
      *
@@ -379,6 +385,14 @@ public:
      */
     void getInducedDipoles(Context& context, std::vector<Vec3>& dipoles);
 
+    /**
+     * Get the total dipole moments (fixed plus induced) of all particles.
+     *
+     * @param context         the Context for which to get the total dipoles
+     * @param[out] dipoles    the total 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

@@ -59,7 +59,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the AmoebaBondForce this kernel will be used for
      */
@@ -99,7 +99,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the AmoebaAngleForce this kernel will be used for
      */
@@ -139,7 +139,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the AmoebaInPlaneAngleForce this kernel will be used for
      */
@@ -179,7 +179,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the PiTorsionForce this kernel will be used for
      */
@@ -219,7 +219,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the StretchBendForce this kernel will be used for
      */
@@ -259,7 +259,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the OutOfPlaneBendForce this kernel will be used for
      */
@@ -299,7 +299,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the TorsionTorsionForce this kernel will be used for
      */
@@ -332,7 +332,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the MultipoleForce this kernel will be used for
      */
@@ -348,7 +348,9 @@ public:
      */
     virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
 
+    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;
@@ -364,7 +366,7 @@ public:
 
     /**
      * Get the parameters being used for PME.
-     * 
+     *
      * @param alpha   the separation parameter
      * @param nx      the number of grid points along the X axis
      * @param ny      the number of grid points along the Y axis
@@ -389,7 +391,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the GBSAOBCForce this kernel will be used for
      */
@@ -429,7 +431,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the GBSAOBCForce this kernel will be used for
      */
@@ -469,7 +471,7 @@ public:
 
     /**
      * Initialize the kernel.
-     * 
+     *
      * @param system     the System this kernel will be applied to
      * @param force      the GBSAOBCForce this kernel will be used for
      */

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

@@ -64,7 +64,7 @@ public:
 
     /**
      * Get the CovalentMap for an atom
-     * 
+     *
      * @param force                AmoebaMultipoleForce force reference
      * @param index                the index of the atom for which to set parameters
      * @param minCovalentIndex     minimum covalent index
@@ -76,13 +76,14 @@ public:
 
     /**
      * Get the covalent degree for the  CovalentEnd lists
-     * 
+     *
      * @param force                AmoebaMultipoleForce force reference
      * @param covalentDegree      covalent degrees for the CovalentEnd lists
      */
     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);
@@ -90,7 +91,7 @@ public:
     void getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments);
     void updateParametersInContext(ContextImpl& context);
     void getPMEParameters(double& alpha, int& nx, int& ny, int& nz) const;
- 
+
 
 private:
     const AmoebaMultipoleForce& owner;

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -147,13 +147,13 @@ void AmoebaMultipoleForce::setEwaldErrorTolerance(double tol) {
     ewaldErrorTol = tol;
 }
 
-int AmoebaMultipoleForce::addMultipole(double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole, int axisType, 
+int AmoebaMultipoleForce::addMultipole(double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole, int axisType,
                                        int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) {
     multipoles.push_back(MultipoleInfo(charge, molecularDipole, molecularQuadrupole,  axisType, multipoleAtomZ,  multipoleAtomX, multipoleAtomY, thole, dampingFactor, polarity));
     return multipoles.size()-1;
 }
 
-void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole, 
+void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std::vector<double>& molecularDipole, std::vector<double>& molecularQuadrupole,
                                                   int& axisType, int& multipoleAtomZ, int& multipoleAtomX, int& multipoleAtomY, double& thole, double& dampingFactor, double& polarity) const {
     charge                      = multipoles[index].charge;
 
@@ -183,7 +183,7 @@ void AmoebaMultipoleForce::getMultipoleParameters(int index, double& charge, std
     polarity                    = multipoles[index].polarity;
 }
 
-void AmoebaMultipoleForce::setMultipoleParameters(int index, double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole, 
+void AmoebaMultipoleForce::setMultipoleParameters(int index, double charge, const std::vector<double>& molecularDipole, const std::vector<double>& molecularQuadrupole,
                                                   int axisType, int multipoleAtomZ, int multipoleAtomX, int multipoleAtomY, double thole, double dampingFactor, double polarity) {
 
     multipoles[index].charge                      = charge;
@@ -250,6 +250,14 @@ void AmoebaMultipoleForce::getInducedDipoles(Context& context, vector<Vec3>& dip
     dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getInducedDipoles(getContextImpl(context), dipoles);
 }
 
+void AmoebaMultipoleForce::getLabFramePermanentDipoles(Context& context, vector<Vec3>& dipoles) {
+    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

@@ -61,7 +61,7 @@ void AmoebaMultipoleForceImpl::initialize(ContextImpl& context) {
         double cutoff = owner.getCutoffDistance();
         if (cutoff > 0.5*boxVectors[0][0] || cutoff > 0.5*boxVectors[1][1] || cutoff > 0.5*boxVectors[2][2])
             throw OpenMMException("AmoebaMultipoleForce: The cutoff distance cannot be greater than half the periodic box size.");
-    }   
+    }
 
     double quadrupoleValidationTolerance = 1.0e-05;
     for (int ii = 0; ii < system.getNumParticles(); ii++) {
@@ -170,7 +170,7 @@ void AmoebaMultipoleForceImpl::getCovalentRange(const AmoebaMultipoleForce& forc
                *maxCovalentIndex = covalentList[ii];
             }
         }
-    }   
+    }
     return;
 }
 
@@ -179,14 +179,22 @@ void AmoebaMultipoleForceImpl::getCovalentDegree(const AmoebaMultipoleForce& for
     const int* CovalentDegrees = AmoebaMultipoleForceImpl::getCovalentDegrees();
     for (unsigned int kk = 0; kk < AmoebaMultipoleForce::CovalentEnd; kk++) {
         covalentDegree[kk] = CovalentDegrees[kk];
-    }   
+    }
     return;
 }
 
+void AmoebaMultipoleForceImpl::getLabFramePermanentDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
+    kernel.getAs<CalcAmoebaMultipoleForceKernel>().getLabFramePermanentDipoles(context, dipoles);
+}
+
 void AmoebaMultipoleForceImpl::getInducedDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
     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

@@ -1992,6 +1992,27 @@ void CudaCalcAmoebaMultipoleForceKernel::ensureMultipolesValid(ContextImpl& cont
         context.calcForcesAndEnergy(false, false, -1);
 }
 
+
+void CudaCalcAmoebaMultipoleForceKernel::getLabFramePermanentDipoles(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> labDipoleVec;
+        labFrameDipoles->download(labDipoleVec);
+        for (int i = 0; i < numParticles; i++)
+            dipoles[order[i]] = Vec3(labDipoleVec[3*i], labDipoleVec[3*i+1], labDipoleVec[3*i+2]);
+    }
+    else {
+        vector<float> labDipoleVec;
+        labFrameDipoles->download(labDipoleVec);
+        for (int i = 0; i < numParticles; i++)
+            dipoles[order[i]] = Vec3(labDipoleVec[3*i], labDipoleVec[3*i+1], labDipoleVec[3*i+2]);
+    }
+}
+
+
 void CudaCalcAmoebaMultipoleForceKernel::getInducedDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
     ensureMultipolesValid(context);
     int numParticles = cu.getNumAtoms();
@@ -2011,6 +2032,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<double4> posqVec;
+        vector<double> labDipoleVec;
+        vector<double> inducedDipoleVec;
+        double totalDipoleVecX;
+        double totalDipoleVecY;
+        double totalDipoleVecZ;
+        inducedDipole->download(inducedDipoleVec);
+        labFrameDipoles->download(labDipoleVec);
+        cu.getPosq().download(posqVec);
+        for (int i = 0; i < numParticles; i++) {
+            totalDipoleVecX = labDipoleVec[3*i] + inducedDipoleVec[3*i];
+            totalDipoleVecY = labDipoleVec[3*i+1] + inducedDipoleVec[3*i+1];
+            totalDipoleVecZ = labDipoleVec[3*i+2] + inducedDipoleVec[3*i+2];
+            dipoles[order[i]] = Vec3(totalDipoleVecX, totalDipoleVecY, totalDipoleVecZ);
+        }
+    }
+    else {
+        vector<float4> posqVec;
+        vector<float> labDipoleVec;
+        vector<float> inducedDipoleVec;
+        float totalDipoleVecX;
+        float totalDipoleVecY;
+        float totalDipoleVecZ;
+        inducedDipole->download(inducedDipoleVec);
+        labFrameDipoles->download(labDipoleVec);
+        cu.getPosq().download(posqVec);
+        for (int i = 0; i < numParticles; i++) {
+            totalDipoleVecX = labDipoleVec[3*i] + inducedDipoleVec[3*i];
+            totalDipoleVecY = labDipoleVec[3*i+1] + inducedDipoleVec[3*i+1];
+            totalDipoleVecZ = labDipoleVec[3*i+2] + inducedDipoleVec[3*i+2];
+            dipoles[order[i]] = Vec3(totalDipoleVecX, totalDipoleVecY, totalDipoleVecZ);
+        }
+    }
+}
+
 void CudaCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl& context, const vector<Vec3>& inputGrid, vector<double>& outputElectrostaticPotential) {
     ensureMultipolesValid(context);
     int numPoints = inputGrid.size();
@@ -2164,6 +2227,7 @@ void CudaCalcAmoebaMultipoleForceKernel::computeSystemMultipoleMoments(ContextIm
     outputMultipoleMoments[12] = 100.0*zzqdp*debye;
 }
 
+
 void CudaCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, vector<double>& outputMultipoleMoments) {
     ensureMultipolesValid(context);
     if (cu.getUseDoublePrecision())

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

@@ -328,6 +328,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+     /**
+     * Get the LabFrame 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 getLabFramePermanentDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
     /**
      * Get the induced dipole moments of all particles.
      * 
@@ -335,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/cuda/tests/TestCudaAmoebaMultipoleForce.cpp

@@ -2666,6 +2666,76 @@ static void testParticleInducedDipoles() {
         ASSERT_EQUAL_VEC(expectedDipole[i], dipole[i], 1e-4);
 }
 
+// test querying particle lab frame permanent dipoles
+
+static void testParticleLabFramePermanentDipoles() {
+    int numberOfParticles     = 8;
+    int inputPmeGridDimension = 0;
+    double cutoff             = 9000000.0;
+    std::vector<Vec3> forces;
+    double energy;
+
+    System system;
+    AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
+    setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual, 
+                                             cutoff, inputPmeGridDimension);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    getForcesEnergyMultipoleAmmonia(context, forces, energy);
+    std::vector<Vec3> dipole;
+    amoebaMultipoleForce->getLabFramePermanentDipoles(context, dipole);
+    
+    // Compare to values calculated by TINKER.
+    
+    std::vector<Vec3> expectedDipole(numberOfParticles);
+    expectedDipole[0] = Vec3(0.00876454250, -2.04310718E-06, -0.00227593519);
+    expectedDipole[1] = Vec3(0.000780382180, -0.00432882849, 0.00236926761);
+    expectedDipole[2] = Vec3(0.000801345883, 0.00431830946, 0.00238143437);
+    expectedDipole[3] = Vec3(-0.00109746996, 1.16087953e-5, -0.00487407492);
+    expectedDipole[4] = Vec3(0.00203814102, -2.26554196e-5, 0.00882284298);
+    expectedDipole[5] = Vec3(-0.00239443187, 0.00432388648, 0.000729303209);
+    expectedDipole[6] = Vec3(0.00491086743, 2.86430963e-6, -0.000918996348);
+    expectedDipole[7] = Vec3(-0.00239301946, -0.00432743976, 0.000712674115);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(expectedDipole[i], dipole[i], 1e-4);
+}
+
+// test querying particle total dipoles (fixed + induced)
+
+static void testParticleTotalDipoles() {
+    int numberOfParticles     = 8;
+    int inputPmeGridDimension = 0;
+    double cutoff             = 9000000.0;
+    std::vector<Vec3> forces;
+    double energy;
+
+    System system;
+    AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
+    setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual, 
+                                             cutoff, inputPmeGridDimension);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("CUDA"));
+    getForcesEnergyMultipoleAmmonia(context, forces, energy);
+    std::vector<Vec3> dipole;
+    amoebaMultipoleForce->getTotalDipoles(context, dipole);
+    
+    // Compare to values calculated by TINKER.
+    
+    std::vector<Vec3> expectedDipole(numberOfParticles);
+    expectedDipole[0] = Vec3(0.0119356307, -1.11302433e-6, -0.00296793872);
+    expectedDipole[1] = Vec3(8.60636211e-4, -0.00460821816, 0.00241705344);
+    expectedDipole[2] = Vec3(8.80646403e-4, 0.00459728769, 0.00243013245);
+    expectedDipole[3] = Vec3(-0.00123822377, 1.31555550e-5, -0.00558185336);
+    expectedDipole[4] = Vec3(0.00399455556, -2.27511931e-5, 0.00955607952);
+    expectedDipole[5] = Vec3(-0.00157302682, 0.00354892386, 3.40921137e-4);
+    expectedDipole[6] = Vec3(0.00952428069, 2.14171505e-6, -6.68945865e-4);
+    expectedDipole[7] = Vec3(-0.00157252460, -0.00355015528, 3.27055162e-4);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(expectedDipole[i], dipole[i], 1e-4);
+}
+
+
+
 // test computation of system multipole moments
 
 static void testSystemMultipoleMoments() {

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

@@ -735,6 +735,48 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getInducedDipoles(ContextImpl& con
     delete amoebaReferenceMultipoleForce;
 }
 
+void ReferenceCalcAmoebaMultipoleForceKernel::getLabFramePermanentDipoles(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> labFramePermanentDipoles;
+    amoebaReferenceMultipoleForce->calculateLabFramePermanentDipoles(posData, charges, dipoles, quadrupoles, tholes, 
+            dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, labFramePermanentDipoles);
+    for (int i = 0; i < numParticles; i++)
+        outputDipoles[i] = labFramePermanentDipoles[i];
+    delete amoebaReferenceMultipoleForce;
+}
+
+
+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, totalDipoles);
+
+    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,20 @@ 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 fixed 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 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.h

@@ -533,6 +533,75 @@ 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
+     * @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 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.
      *

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

@@ -2587,6 +2587,76 @@ static void testParticleInducedDipoles() {
         ASSERT_EQUAL_VEC(expectedDipole[i], dipole[i], 1e-4);
 }
 
+// test querying particle lab frame permanent dipoles
+
+static void testParticleLabFramePermanentDipoles() {
+    int numberOfParticles     = 8;
+    int inputPmeGridDimension = 0;
+    double cutoff             = 9000000.0;
+    std::vector<Vec3> forces;
+    double energy;
+
+    System system;
+    AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
+    setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual, 
+                                             cutoff, inputPmeGridDimension);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    getForcesEnergyMultipoleAmmonia(context, forces, energy);
+    std::vector<Vec3> dipole;
+    amoebaMultipoleForce->getLabFramePermanentDipoles(context, dipole);
+    
+    // Compare to values calculated by TINKER.
+    
+    std::vector<Vec3> expectedDipole(numberOfParticles);
+    expectedDipole[0] = Vec3(0.00876454250, -2.04310718E-06, -0.00227593519);
+    expectedDipole[1] = Vec3(0.000780382180, -0.00432882849, 0.00236926761);
+    expectedDipole[2] = Vec3(0.000801345883, 0.00431830946, 0.00238143437);
+    expectedDipole[3] = Vec3(-0.00109746996, 1.16087953e-5, -0.00487407492);
+    expectedDipole[4] = Vec3(0.00203814102, -2.26554196e-5, 0.00882284298);
+    expectedDipole[5] = Vec3(-0.00239443187, 0.00432388648, 0.000729303209);
+    expectedDipole[6] = Vec3(0.00491086743, 2.86430963e-6, -0.000918996348);
+    expectedDipole[7] = Vec3(-0.00239301946, -0.00432743976, 0.000712674115);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(expectedDipole[i], dipole[i], 1e-4);
+}
+
+
+// test querying particle total dipoles (fixed + induced)
+
+static void testParticleTotalDipoles() {
+    int numberOfParticles     = 8;
+    int inputPmeGridDimension = 0;
+    double cutoff             = 9000000.0;
+    std::vector<Vec3> forces;
+    double energy;
+
+    System system;
+    AmoebaMultipoleForce* amoebaMultipoleForce = new AmoebaMultipoleForce();;
+    setupMultipoleAmmonia(system, amoebaMultipoleForce, AmoebaMultipoleForce::NoCutoff, AmoebaMultipoleForce::Mutual, 
+                                             cutoff, inputPmeGridDimension);
+    LangevinIntegrator integrator(0.0, 0.1, 0.01);
+    Context context(system, integrator, Platform::getPlatformByName("Reference"));
+    getForcesEnergyMultipoleAmmonia(context, forces, energy);
+    std::vector<Vec3> dipole;
+    amoebaMultipoleForce->getTotalDipoles(context, dipole);
+    
+    // Compare to values calculated by TINKER.
+    
+    std::vector<Vec3> expectedDipole(numberOfParticles);
+    expectedDipole[0] = Vec3(0.0119356307, -1.11302433e-6, -0.00296793872);
+    expectedDipole[1] = Vec3(8.60636211e-4, -0.00460821816, 0.00241705344);
+    expectedDipole[2] = Vec3(8.80646403e-4, 0.00459728769, 0.00243013245);
+    expectedDipole[3] = Vec3(-0.00123822377, 1.31555550e-5, -0.00558185336);
+    expectedDipole[4] = Vec3(0.00399455556, -2.27511931e-5, 0.00955607952);
+    expectedDipole[5] = Vec3(-0.00157302682, 0.00354892386, 3.40921137e-4);
+    expectedDipole[6] = Vec3(0.00952428069, 2.14171505e-6, -6.68945865e-4);
+    expectedDipole[7] = Vec3(-0.00157252460, -0.00355015528, 3.27055162e-4);
+    for (int i = 0; i < numberOfParticles; i++)
+        ASSERT_EQUAL_VEC(expectedDipole[i], dipole[i], 1e-4);
+}
+
+
 // test computation of system multipole moments
 
 static void testSystemMultipoleMoments() {

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -124,6 +124,8 @@ NO_OUTPUT_ARGS = [('LocalEnergyMinimizer', 'minimize', 'context'),
                   ('AmoebaMultipoleForce', 'setCovalentMap', 'covalentAtoms'),
                   ('AmoebaMultipoleForce', 'getElectrostaticPotential', 'context'),
                   ('AmoebaMultipoleForce', 'getInducedDipoles', 'context'),
+                  ('AmoebaMultipoleForce', 'getLabFramePermanentDipoles', 'context'),
+                  ('AmoebaMultipoleForce', 'getTotalDipoles', 'context'),
 ]
 
 # SWIG assumes the target language shadow class owns the C++ class
@@ -274,6 +276,8 @@ UNITS = {
 #("AmoebaMultipoleForce",                 "getElectrostaticPotential")                     :  ( ('unit.kilojoule_per_mole'), ()),
 ("AmoebaMultipoleForce",                 "getElectrostaticPotential")                     :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getInducedDipoles")                             :  ( None, ()),
+("AmoebaMultipoleForce",                 "getLabFramePermanentDipoles")                   :  ( None, ()),
+("AmoebaMultipoleForce",                 "getTotalDipoles")                               :  ( None, ()),
 ("AmoebaMultipoleForce",                 "getSystemMultipoleMoments")                     :  ( None, ()),
 
 ("AmoebaOutOfPlaneBendForce",            "getNumOutOfPlaneBends")                         :  ( None, ()),