Added AmoebaMultipoleForce::getSystemMultipoleMoments() to get system...

Added AmoebaMultipoleForce::getSystemMultipoleMoments() to get system multipole moments; based on TINKER subroutine moments()

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

@@ -330,6 +330,20 @@ public:
                                     Context& context, std::vector< double >& outputElectrostaticPotential );
 
 
+    /**
+     * Get the system multipole moments
+     *
+     * @param origin       origin
+     * @param context      context
+     * @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( const Vec3& origin, Context& context, std::vector< double >& outputMultipoleMonents );
+
 protected:
     ForceImpl* createImpl();
 private:

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

@@ -374,6 +374,8 @@ public:
     virtual void getElectrostaticPotential( ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                             std::vector< double >& outputElectrostaticPotential ) = 0;
 
+    virtual void getSystemMultipoleMoments( ContextImpl& context, const Vec3& origin, std::vector< double >& outputMultipoleMonents ) = 0;
+
 };
 
 /**

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

@@ -85,6 +85,8 @@ public:
     void getElectrostaticPotential( ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                     std::vector< double >& outputElectrostaticPotential );
 
+    void getSystemMultipoleMoments( ContextImpl& context, const Vec3& origin, std::vector< double >& outputMultipoleMonents );
+ 
 
 private:
     AmoebaMultipoleForce& owner;

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -243,6 +243,10 @@ void AmoebaMultipoleForce::getElectrostaticPotential( const std::vector< Vec3 >&
     dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getElectrostaticPotential(getContextImpl(context), inputGrid, outputElectrostaticPotential);
 }
 
+void AmoebaMultipoleForce::getSystemMultipoleMoments( const Vec3& origin, Context& context, std::vector< double >& outputMultipoleMonents ){
+    dynamic_cast<AmoebaMultipoleForceImpl&>(getImplInContext(context)).getSystemMultipoleMoments(getContextImpl(context), origin, outputMultipoleMonents);
+}
+
 ForceImpl* AmoebaMultipoleForce::createImpl() {
     return new AmoebaMultipoleForceImpl(*this);
 }

plugins/amoeba/openmmapi/src/AmoebaMultipoleForceImpl.cpp

@@ -144,3 +144,7 @@ void AmoebaMultipoleForceImpl::getElectrostaticPotential( ContextImpl& context,
     kernel.getAs<CalcAmoebaMultipoleForceKernel>().getElectrostaticPotential(context, inputGrid, outputElectrostaticPotential);
 }
 
+void AmoebaMultipoleForceImpl::getSystemMultipoleMoments( ContextImpl& context, const Vec3& origin, std::vector< double >& outputMultipoleMonents ){
+    kernel.getAs<CalcAmoebaMultipoleForceKernel>().getSystemMultipoleMoments(context, origin, outputMultipoleMonents);
+}
+

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

@@ -860,6 +860,17 @@ static void computeAmoebaMultipolePotential( AmoebaCudaData& data, const std::ve
     }
 }
 
+static void computeAmoebaSystemMultipoleMoments( AmoebaCudaData& data, const Vec3& origin,
+                                                 std::vector< double >& outputMultipoleMonents) {
+
+    amoebaGpuContext gpu = data.getAmoebaGpu();
+
+    data.setGpuInitialized( false );
+    data.initializeGpu();
+    kCalculateAmoebaSystemMultipoleMoments( gpu, origin, outputMultipoleMonents );
+
+}
+
 class CudaCalcAmoebaMultipoleForceKernel::ForceInfo : public CudaForceInfo {
 public:
     ForceInfo(const AmoebaMultipoleForce& force) : force(force) {
@@ -1075,6 +1086,12 @@ void CudaCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextImpl&
     return;
 }
 
+void CudaCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context,  const Vec3& origin,
+                                                                   std::vector< double >& outputMultipoleMonents) {
+    computeAmoebaSystemMultipoleMoments( data, origin, outputMultipoleMonents);
+    return;
+}
+
 /* -------------------------------------------------------------------------- *
  *                       AmoebaGeneralizedKirkwood                            *
  * -------------------------------------------------------------------------- */

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

@@ -344,6 +344,22 @@ public:
      */
     void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                    std::vector< double >& outputElectrostaticPotential );
+
+   /** 
+     * Get the system multipole moments
+     *
+     * @param origin       origin
+     * @param context      context
+     * @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( ContextImpl& context, const Vec3& origin, std::vector< double >& outputMultipoleMonents );
+
+
 private:
     class ForceInfo;
     int numMultipoles;

plugins/amoeba/platforms/cuda/src/kernels/amoebaCudaKernels.h

@@ -28,6 +28,7 @@
  * -------------------------------------------------------------------------- */
 
 #include "amoebaGpuTypes.h"
+#include "openmm/Vec3.h"
 
 #include <string>
 #include <vector>
@@ -53,6 +54,7 @@ extern void SetCalculateAmoebaMultipoleForcesSim(amoebaGpuContext gpu);
 extern void GetCalculateAmoebaMultipoleForcesSim(amoebaGpuContext gpu);
 extern void kCalculateAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool performGk );
 extern void kCalculateAmoebaMultipolePotential(amoebaGpuContext amoebaGpu );
+extern void kCalculateAmoebaSystemMultipoleMoments(amoebaGpuContext amoebaGpu, const OpenMM::Vec3& origin, std::vector< double >& outputMultipoleMonents );
 
 // vdw
 

plugins/amoeba/platforms/cuda/src/kernels/kCalculateAmoebaCudaRotateFrame.cu

@@ -422,9 +422,9 @@ void cudaComputeAmoebaLabFrameMoments( amoebaGpuContext amoebaGpu )
 
 }
 
-void kCalculateAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool hasAmoebaGeneralizedKirkwood ) 
+static void kSetupAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool hasAmoebaGeneralizedKirkwood ) 
 {
-    std::string methodName = "kCalculateAmoebaMultipoleForces";
+    std::string methodName  = "kSetupAmoebaMultipoleForces";
 
     // compute lab frame moments
 
@@ -484,6 +484,13 @@ void kCalculateAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool hasAmoebaG
        throw OpenMM::OpenMMException("Induced dipole calculation did not converge" );
     }
 
+}
+
+void kCalculateAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool hasAmoebaGeneralizedKirkwood ) 
+{
+
+    kSetupAmoebaMultipoleForces(amoebaGpu, hasAmoebaGeneralizedKirkwood ); 
+
     // calculate electrostatic forces
 
     if( amoebaGpu->multipoleNonbondedMethod  == AMOEBA_NO_CUTOFF ){
@@ -496,62 +503,156 @@ void kCalculateAmoebaMultipoleForces(amoebaGpuContext amoebaGpu, bool hasAmoebaG
 
 void kCalculateAmoebaMultipolePotential(amoebaGpuContext amoebaGpu ) 
 {
-    std::string methodName = "kCalculateAmoebaMultipolePotential";
 
-    // compute lab frame moments
+    // setup
 
-    cudaComputeAmoebaLabFrameMoments( amoebaGpu );
+    kSetupAmoebaMultipoleForces(amoebaGpu, false ); 
+
+    // calculate electrostatic potential
+
+    cudaComputeAmoebaElectrostaticPotential( amoebaGpu );
+
+}
+
+void kCalculateAmoebaSystemMultipoleMoments(amoebaGpuContext amoebaGpu, const OpenMM::Vec3& origin, std::vector< double >& outputMultipoleMoments ) 
+{
+
+    // setup
+
+    kSetupAmoebaMultipoleForces(amoebaGpu, false ); 
 
-    if( 0 ){
     gpuContext gpu         = amoebaGpu->gpuContext;
-        std::vector<int> fileId;
-        //fileId.push_back( 0 );
-        VectorOfDoubleVectors outputVector;
-        //cudaLoadCudaFloat4Array( gpu->natoms, 3, gpu->psPosq4,              outputVector, gpu->psAtomIndex->_pSysData, 1.0f );
-        cudaLoadCudaFloatArray( gpu->natoms,  3, amoebaGpu->psLabFrameDipole,     outputVector, gpu->psAtomIndex->_pSysData, 1.0f );
-        cudaLoadCudaFloatArray( gpu->natoms,  9, amoebaGpu->psLabFrameQuadrupole, outputVector, gpu->psAtomIndex->_pSysData, 1.0f );
-        cudaWriteVectorOfDoubleVectorsToFile( "CudaLabMoments", fileId, outputVector );
+
+    gpu->psPosq4->Download();
+    gpu->psVelm4->Download();
+    float4* posq           = gpu->psPosq4->_pSysData;    
+    float4* velm           = gpu->psVelm4->_pSysData;    
+    float totalMass        = 0.0f;
+    float centerOfMass[3]  = { 0.0f, 0.0f, 0.0f };
+    for( unsigned int ii  = 0; ii < gpu->natoms; ii++ ){
+        float mass;
+        if( velm->w > 0.0f ){
+            mass        = 1.0f/velm[ii].w;
+        } else {
+            mass        = 0.0f;
+        }
+        totalMass        += mass;
+        centerOfMass[0]  += mass*posq[ii].x;
+        centerOfMass[1]  += mass*posq[ii].y;
+        centerOfMass[2]  += mass*posq[ii].z;
     }
 
-    // compute fixed E-field and mutual induced field 
+    std::vector<float4>  posqLocal(gpu->natoms);
+    if( totalMass > 0.0f ){
+        centerOfMass[0]  /= totalMass;
+        centerOfMass[1]  /= totalMass;
+        centerOfMass[2]  /= totalMass;
+    }
+    for( unsigned int ii  = 0; ii < gpu->natoms; ii++ ){
+        posqLocal[ii].x = posq[ii].x  - centerOfMass[0];
+        posqLocal[ii].y = posq[ii].y  - centerOfMass[1];
+        posqLocal[ii].z = posq[ii].z  - centerOfMass[2];
+        posqLocal[ii].w = posq[ii].w;
+    }
 
-    if( amoebaGpu->multipoleNonbondedMethod == AMOEBA_NO_CUTOFF ){
+    float netchg  = 0.0f;
 
-        cudaComputeAmoebaFixedEField( amoebaGpu );
-        cudaComputeAmoebaMutualInducedField( amoebaGpu );
+    float xdpl    = 0.0f;
+    float ydpl    = 0.0f;
+    float zdpl    = 0.0f;
 
-    } else {
+    float xxqdp   = 0.0f;
+    float xyqdp   = 0.0f;
+    float xzqdp   = 0.0f;
 
-        gpuContext gpu = amoebaGpu->gpuContext;
-        kFindBlockBoundsPeriodic_kernel<<<(gpu->psGridBoundingBox->_length+63)/64, 64>>>();
-        LAUNCHERROR("kFindBlockBoundsPeriodic");
-        kFindBlocksWithInteractionsPeriodic_kernel<<<gpu->sim.interaction_blocks, gpu->sim.interaction_threads_per_block>>>();
-        LAUNCHERROR("kFindBlocksWithInteractionsPeriodic");
+    float yxqdp   = 0.0f;
+    float yyqdp   = 0.0f;
+    float yzqdp   = 0.0f;
 
-        compactStream(gpu->compactPlan, gpu->sim.pInteractingWorkUnit, gpu->sim.pWorkUnit, gpu->sim.pInteractionFlag, gpu->sim.workUnits, gpu->sim.pInteractionCount);
+    float zxqdp   = 0.0f;
+    float zyqdp   = 0.0f;
+    float zzqdp   = 0.0f;
 
-        //compactStream( gpu->compactPlan, 
-        //               gpu->sim.pInteractingWorkUnit, unsigned int* dOut
-        //               amoebaGpu->psWorkUnit->_pDevData, const unsigned int* dIn
-        //               gpu->sim.pInteractionFlag,        const unsigned int* dValid
-        //               gpu->sim.workUnits,               gpu
-        //               gpu->sim.pInteractionCount);
-        kFindInteractionsWithinBlocksPeriodic_kernel<<<gpu->sim.nonbond_blocks, gpu->sim.nonbond_threads_per_block,
-                sizeof(unsigned int)*gpu->sim.nonbond_threads_per_block>>>(gpu->sim.pInteractingWorkUnit);
-        LAUNCHERROR("kFindInteractionsWithinBlocksPeriodic");
+    amoebaGpu->psLabFrameDipole->Download();
+    float* labFrameDipole      = amoebaGpu->psLabFrameDipole->_pSysData;    
 
-        cudaComputeAmoebaPmeFixedEField( amoebaGpu );
-        cudaComputeAmoebaPmeMutualInducedField( amoebaGpu );
-    }
+    amoebaGpu->psInducedDipole->Download();
+    float* inducedDipole       = amoebaGpu->psInducedDipole->_pSysData;    
 
-    // check if induce dipole calculation converged -- abort if it did not
+    amoebaGpu->psLabFrameQuadrupole->Download();
+    float* labFrameQuadrupole  = amoebaGpu->psLabFrameQuadrupole->_pSysData;    
+    for( unsigned int ii  = 0; ii < gpu->natoms; ii++ ){
+
+        netchg              += posqLocal[ii].w;
+
+        float netDipoleX     = (labFrameDipole[3*ii]    + inducedDipole[3*ii]);
+        float netDipoleY     = (labFrameDipole[3*ii+1]  + inducedDipole[3*ii+1]);
+        float netDipoleZ     = (labFrameDipole[3*ii+2]  + inducedDipole[3*ii+2]);
+
+        xdpl    += posqLocal[ii].x*posqLocal[ii].w + netDipoleX;
+        ydpl    += posqLocal[ii].y*posqLocal[ii].w + netDipoleY;
+        zdpl    += posqLocal[ii].z*posqLocal[ii].w + netDipoleZ;
+
+        xxqdp   += posqLocal[ii].x*posqLocal[ii].x*posqLocal[ii].w + 2.0f*posqLocal[ii].x*netDipoleX;
+        xyqdp   += posqLocal[ii].x*posqLocal[ii].y*posqLocal[ii].w + posqLocal[ii].x*netDipoleY + posqLocal[ii].y*netDipoleX;
+        xzqdp   += posqLocal[ii].x*posqLocal[ii].z*posqLocal[ii].w + posqLocal[ii].x*netDipoleZ + posqLocal[ii].z*netDipoleX;
+
+        yxqdp   += posqLocal[ii].y*posqLocal[ii].x*posqLocal[ii].w + posqLocal[ii].y*netDipoleX + posqLocal[ii].x*netDipoleY;
+        yyqdp   += posqLocal[ii].y*posqLocal[ii].y*posqLocal[ii].w + 2.0f*posqLocal[ii].y*netDipoleY;
+        yzqdp   += posqLocal[ii].y*posqLocal[ii].z*posqLocal[ii].w + posqLocal[ii].y*netDipoleZ + posqLocal[ii].z*netDipoleY;
+
+        zxqdp   += posqLocal[ii].z*posqLocal[ii].x*posqLocal[ii].w + posqLocal[ii].z*netDipoleX + posqLocal[ii].x*netDipoleZ;
+        zyqdp   += posqLocal[ii].z*posqLocal[ii].y*posqLocal[ii].w + posqLocal[ii].z*netDipoleY + posqLocal[ii].y*netDipoleZ;
+        zzqdp   += posqLocal[ii].z*posqLocal[ii].z*posqLocal[ii].w + 2.0f*posqLocal[ii].z*netDipoleZ;
 
-    if( amoebaGpu->mutualInducedDone == 0 ){
-       throw OpenMM::OpenMMException("Induced dipole calculation did not converge" );
     }
 
-    // calculate electrostatic potential
+//  convert the quadrupole from traced to traceless form
+ 
+    float qave   = (xxqdp + yyqdp + zzqdp)/3.0f;
+          xxqdp  = 1.5f*(xxqdp-qave);
+          xyqdp  = 1.5f*xyqdp;
+          xzqdp  = 1.5f*xzqdp;
+          yxqdp  = 1.5f*yxqdp;
+          yyqdp  = 1.5f*(yyqdp-qave);
+          yzqdp  = 1.5f*yzqdp;
+          zxqdp  = 1.5f*zxqdp;
+          zyqdp  = 1.5f*zyqdp;
+          zzqdp  = 1.5f*(zzqdp-qave);
+
+//  add the traceless atomic quadrupoles to total quadrupole
+
+    for( unsigned int ii  = 0; ii < gpu->natoms; ii++ ){
+        xxqdp  = xxqdp + 3.0f*labFrameQuadrupole[9*ii];
+        xyqdp  = xyqdp + 3.0f*labFrameQuadrupole[9*ii+1];
+        xzqdp  = xzqdp + 3.0f*labFrameQuadrupole[9*ii+2];
+        yxqdp  = yxqdp + 3.0f*labFrameQuadrupole[9*ii+3];
+        yyqdp  = yyqdp + 3.0f*labFrameQuadrupole[9*ii+4];
+        yzqdp  = yzqdp + 3.0f*labFrameQuadrupole[9*ii+5];
+        zxqdp  = zxqdp + 3.0f*labFrameQuadrupole[9*ii+6];
+        zyqdp  = zyqdp + 3.0f*labFrameQuadrupole[9*ii+7];
+        zzqdp  = zzqdp + 3.0f*labFrameQuadrupole[9*ii+8];
+    }
  
-    cudaComputeAmoebaElectrostaticPotential( amoebaGpu );
+    float debye                = 4.80321f;
+    outputMultipoleMoments.resize( 13 );
+
+    outputMultipoleMoments[0]  = netchg;
+
+    outputMultipoleMoments[1]  = xdpl*debye;
+    outputMultipoleMoments[2]  = ydpl*debye;
+    outputMultipoleMoments[3]  = zdpl*debye;
+    
+    outputMultipoleMoments[4]  = xxqdp*debye;
+    outputMultipoleMoments[5]  = xyqdp*debye;
+    outputMultipoleMoments[6]  = xzqdp*debye;
+
+    outputMultipoleMoments[7]  = yxqdp*debye;
+    outputMultipoleMoments[8]  = yyqdp*debye;
+    outputMultipoleMoments[9]  = yzqdp*debye;
+
+    outputMultipoleMoments[10] = zxqdp*debye;
+    outputMultipoleMoments[11] = zyqdp*debye;
+    outputMultipoleMoments[12] = zzqdp*debye;
 
 }

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

@@ -568,6 +568,11 @@ void ReferenceCalcAmoebaMultipoleForceKernel::getElectrostaticPotential(ContextI
     return;
 }
 
+void ReferenceCalcAmoebaMultipoleForceKernel::getSystemMultipoleMoments(ContextImpl& context, const Vec3& origin,
+                                                                        std::vector< double >& outputMultipoleMonents){
+    return;
+}
+
 ///* -------------------------------------------------------------------------- *
 // *                       AmoebaGeneralizedKirkwood                            *
 // * -------------------------------------------------------------------------- */

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

@@ -393,6 +393,20 @@ public:
      */
     void getElectrostaticPotential(ContextImpl& context, const std::vector< Vec3 >& inputGrid,
                                    std::vector< double >& outputElectrostaticPotential );
+
+    /**
+     * Get the system multipole moments
+     *
+     * @param origin       origin
+     * @param context      context
+     * @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(ContextImpl& context, const Vec3& origin, std::vector< double >& outputMultipoleMonents);
+
 private:
     int numMultipoles;
     int polarizationType;

wrappers/python/src/swig_doxygen/swigInputConfig.py

@@ -290,6 +290,7 @@ UNITS = {
 #("AmoebaMultipoleForce",                 "getElectrostaticPotential")                     :  ( None, ('unit.kilojoule_per_mole')),
 #("AmoebaMultipoleForce",                 "getElectrostaticPotential")                     :  ( ('unit.kilojoule_per_mole'), ()),
 ("AmoebaMultipoleForce",                 "getElectrostaticPotential")                     :  ( None, ()),
+("AmoebaMultipoleForce",                 "getSystemMultipoleMoments")                     :  ( None, ()),
 
 ("AmoebaOutOfPlaneBendForce",            "getNumOutOfPlaneBends")                         :  ( None, ()),
 ("AmoebaOutOfPlaneBendForce",            "getAmoebaGlobalOutOfPlaneBendCubic")            :  ( None,()),