Merge pull request #176 from peastman/master

Added AmoebaMultipoleForce::getInducedDipoles()

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

@@ -301,6 +301,14 @@ public:
      */
     void setEwaldErrorTolerance(double tol);
 
+    /**
+     * Get the induced 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 getInducedDipoles(Context& context, std::vector<Vec3>& dipoles);
+    
     /**
      * Get the electrostatic potential.
      *

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

@@ -348,6 +348,8 @@ public:
      */
     virtual double execute(ContextImpl& context, bool includeForces, bool includeEnergy) = 0;
 
+    virtual void getInducedDipoles(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

@@ -82,6 +82,8 @@ public:
      */
     static void getCovalentDegree( const AmoebaMultipoleForce& force, std::vector<int>& covalentDegree );
 
+    void getInducedDipoles(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

@@ -226,6 +226,10 @@ void AmoebaMultipoleForce::getCovalentMaps(int index, std::vector< std::vector<i
     }
 }
 
+void AmoebaMultipoleForce::getInducedDipoles(Context& context, vector<Vec3>& dipoles) {
+    dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getInducedDipoles(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

@@ -183,6 +183,10 @@ void AmoebaMultipoleForceImpl::getCovalentDegree( const AmoebaMultipoleForce& fo
     return;
 }
 
+void AmoebaMultipoleForceImpl::getInducedDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
+    kernel.getAs<CalcAmoebaMultipoleForceKernel>().getInducedDipoles(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

@@ -1633,6 +1633,24 @@ void CudaCalcAmoebaMultipoleForceKernel::ensureMultipolesValid(ContextImpl& cont
         context.calcForcesAndEnergy(false, false, -1);
 }
 
+void CudaCalcAmoebaMultipoleForceKernel::getInducedDipoles(ContextImpl& context, vector<Vec3>& dipoles) {
+    ensureMultipolesValid(context);
+    int numParticles = cu.getNumAtoms();
+    dipoles.resize(numParticles);
+    if (cu.getUseDoublePrecision()) {
+        vector<double> d;
+        inducedDipole->download(d);
+        for (int i = 0; i < numParticles; i++)
+            dipoles[i] = Vec3(d[3*i], d[3*i+1], d[3*i+2]);
+    }
+    else {
+        vector<float> d;
+        inducedDipole->download(d);
+        for (int i = 0; i < numParticles; i++)
+            dipoles[i] = Vec3(d[3*i], d[3*i+1], d[3*i+2]);
+    }
+}
+
 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

@@ -327,6 +327,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Get the induced 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 getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
     /**
      * Execute the kernel to calculate the electrostatic potential
      *

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

@@ -2698,6 +2698,40 @@ static void testPMEMutualPolarizationLargeWater( FILE* log ) {
 
 }
 
+// test querying particle induced dipoles
+
+static void testParticleInducedDipoles() {
+    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->getInducedDipoles(context, dipole);
+    
+    // Compare to values calculated by TINKER.
+    
+    std::vector<Vec3> expectedDipole(numberOfParticles);
+    expectedDipole[0] = Vec3(0.0031710288, 9.3687453e-7, -0.0006919963);
+    expectedDipole[1] = Vec3(8.0279737504e-5, -0.000279376, 4.778060103e-5);
+    expectedDipole[2] = Vec3(0.000079322, 0.0002789804, 4.8696656126e-5);
+    expectedDipole[3] = Vec3(-0.0001407394, 1.540638116e-6, -0.0007077775);
+    expectedDipole[4] = Vec3(0.0019564439, -1.0409717e-7, 0.0007332188);
+    expectedDipole[5] = Vec3(0.0008213891, -0.0007749618, -0.0003883865);
+    expectedDipole[6] = Vec3(0.0046133992, -7.2868019e-7, 0.0002500622);
+    expectedDipole[7] = Vec3(0.0008204731, 0.0007772727, -0.0003856176);
+    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( FILE* log ) {
@@ -2892,6 +2926,10 @@ int main(int argc, char* argv[]) {
         testMultipoleIonsAndWaterPMEMutualPolarization( log );
         testMultipoleIonsAndWaterPMEDirectPolarization( log );
 
+        // test querying induced dipoles
+        
+        testParticleInducedDipoles();
+        
         // test computation of system multipole moments
 
         testSystemMultipoleMoments( log );

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

@@ -683,6 +683,25 @@ double ReferenceCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bo
     return static_cast<double>(energy);
 }
 
+void ReferenceCalcAmoebaMultipoleForceKernel::getInducedDipoles(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 induced dipoles.
+    
+    vector<RealVec> inducedDipoles;
+    amoebaReferenceMultipoleForce->calculateInducedDipoles(posData, charges, dipoles, quadrupoles, tholes,
+            dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs, multipoleAtomCovalentInfo, inducedDipoles);
+    for (int i = 0; i < numParticles; i++)
+        outputDipoles[i] = inducedDipoles[i];
+    delete amoebaReferenceMultipoleForce;
+}
+
 void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                                                         std::vector< double >& outputElectrostaticPotential ){
 
@@ -704,8 +723,6 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextI
     }
 
     delete amoebaReferenceMultipoleForce;
-
-    return;
 }
 
 void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, std::vector< double >& outputMultipoleMoments){
@@ -726,8 +743,6 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextI
                                                                           multipoleAtomCovalentInfo, outputMultipoleMoments );
 
     delete amoebaReferenceMultipoleForce;
-
-    return;
 }
 
 void ReferenceCalcAmoebaMultipoleForceKernel::copyParametersToContext(ContextImpl& context, const AmoebaMultipoleForce& force) {

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

@@ -366,6 +366,13 @@ public:
      * @return the potential energy due to the force
      */
     double execute(ContextImpl& context, bool includeForces, bool includeEnergy);
+    /**
+     * Get the induced 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 getInducedDipoles(ContextImpl& context, std::vector<Vec3>& dipoles);
     /** 
      * Calculate the electrostatic potential given vector of grid coordinates.
      *

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

@@ -1563,6 +1563,28 @@ RealOpenMM AmoebaReferenceMultipoleForce::calculateForceAndEnergy( const std::ve
     return energy;
 }
 
+void AmoebaReferenceMultipoleForce::calculateInducedDipoles(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< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
+                                                            std::vector<RealVec>& outputInducedDipoles) {
+    // setup, including calculating induced dipoles
+
+    std::vector<MultipoleParticleData> particleData;
+    setup( particlePositions, charges, dipoles, quadrupoles, tholes,
+           dampingFactors, polarity, axisTypes, multipoleAtomZs, multipoleAtomXs, multipoleAtomYs,
+           multipoleAtomCovalentInfo, particleData );
+    outputInducedDipoles = _inducedDipole;
+}
+
 void AmoebaReferenceMultipoleForce::calculateAmoebaSystemMultipoleMoments( const std::vector<RealOpenMM>& masses,
                                                                            const std::vector<RealVec>& particlePositions,
                                                                            const std::vector<RealOpenMM>& charges,

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

@@ -471,6 +471,38 @@ public:
                                         const std::vector< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
                                         std::vector<OpenMM::RealVec>& forces );
 
+    /**
+     * Calculate particle induced 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 calculateInducedDipoles(const std::vector<OpenMM::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< std::vector< std::vector<int> > >& multipoleAtomCovalentInfo,
+                                 std::vector<RealVec>& outputInducedDipoles);
+
     /**
      * Calculate system multipole moments.
      *

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

@@ -2605,6 +2605,40 @@ static void testPMEMutualPolarizationLargeWater( FILE* log ) {
 
 }
 
+// test querying particle induced dipoles
+
+static void testParticleInducedDipoles() {
+    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->getInducedDipoles(context, dipole);
+    
+    // Compare to values calculated by TINKER.
+    
+    std::vector<Vec3> expectedDipole(numberOfParticles);
+    expectedDipole[0] = Vec3(0.0031710288, 9.3687453e-7, -0.0006919963);
+    expectedDipole[1] = Vec3(8.0279737504e-5, -0.000279376, 4.778060103e-5);
+    expectedDipole[2] = Vec3(0.000079322, 0.0002789804, 4.8696656126e-5);
+    expectedDipole[3] = Vec3(-0.0001407394, 1.540638116e-6, -0.0007077775);
+    expectedDipole[4] = Vec3(0.0019564439, -1.0409717e-7, 0.0007332188);
+    expectedDipole[5] = Vec3(0.0008213891, -0.0007749618, -0.0003883865);
+    expectedDipole[6] = Vec3(0.0046133992, -7.2868019e-7, 0.0002500622);
+    expectedDipole[7] = Vec3(0.0008204731, 0.0007772727, -0.0003856176);
+    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( FILE* log ) {
@@ -2778,6 +2812,10 @@ int main( int numberOfArguments, char* argv[] ) {
 
         testMultipoleAmmoniaDirectPolarization( log );
 
+        // test querying induced dipoles
+        
+        testParticleInducedDipoles();
+
         // test mutual polarization, no cutoff
 
         testMultipoleAmmoniaMutualPolarization( log );