Began CUDA implementation of spherical harmonics for multipoles

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

@@ -806,8 +806,8 @@ private:
 
 CudaCalcAmoebaMultipoleForceKernel::CudaCalcAmoebaMultipoleForceKernel(std::string name, const Platform& platform, CudaContext& cu, const System& system) : 
         CalcAmoebaMultipoleForceKernel(name, platform), cu(cu), system(system), hasInitializedScaleFactors(false), hasInitializedFFT(false), multipolesAreValid(false),
-        multipoleParticles(NULL), molecularDipoles(NULL), molecularQuadrupoles(NULL), labFrameDipoles(NULL), labFrameQuadrupoles(NULL), fracDipoles(NULL),
-        fracQuadrupoles(NULL), field(NULL), fieldPolar(NULL), inducedField(NULL), inducedFieldPolar(NULL), torque(NULL), dampingAndThole(NULL), inducedDipole(NULL),
+        multipoleParticles(NULL), molecularDipoles(NULL), molecularQuadrupoles(NULL), labFrameDipoles(NULL), labFrameQuadrupoles(NULL), sphericalDipoles(NULL), sphericalQuadrupoles(NULL),
+        fracDipoles(NULL), fracQuadrupoles(NULL), field(NULL), fieldPolar(NULL), inducedField(NULL), inducedFieldPolar(NULL), torque(NULL), dampingAndThole(NULL), inducedDipole(NULL),
         diisCoefficients(NULL), inducedDipolePolar(NULL), inducedDipoleErrors(NULL), prevDipoles(NULL), prevDipolesPolar(NULL), prevDipolesGk(NULL),
         prevDipolesGkPolar(NULL), prevErrors(NULL), diisMatrix(NULL), polarizability(NULL), covalentFlags(NULL), polarizationGroupFlags(NULL),
         pmeGrid(NULL), pmeBsplineModuliX(NULL), pmeBsplineModuliY(NULL), pmeBsplineModuliZ(NULL), pmeIgrid(NULL), pmePhi(NULL),
@@ -826,6 +826,10 @@ CudaCalcAmoebaMultipoleForceKernel::~CudaCalcAmoebaMultipoleForceKernel() {
         delete labFrameDipoles;
     if (labFrameQuadrupoles != NULL)
         delete labFrameQuadrupoles;
+    if (sphericalDipoles != NULL)
+        delete sphericalDipoles;
+    if (sphericalQuadrupoles != NULL)
+        delete sphericalQuadrupoles;
     if (fracDipoles != NULL)
         delete fracDipoles;
     if (fracQuadrupoles != NULL)
@@ -985,6 +989,8 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
     int elementSize = (cu.getUseDoublePrecision() ? sizeof(double) : sizeof(float));
     labFrameDipoles = new CudaArray(cu, 3*paddedNumAtoms, elementSize, "labFrameDipoles");
     labFrameQuadrupoles = new CudaArray(cu, 5*paddedNumAtoms, elementSize, "labFrameQuadrupoles");
+    sphericalDipoles = new CudaArray(cu, 3*paddedNumAtoms, elementSize, "sphericalDipoles");
+    sphericalQuadrupoles = new CudaArray(cu, 5*paddedNumAtoms, elementSize, "sphericalQuadrupoles");
     fracDipoles = new CudaArray(cu, 3*paddedNumAtoms, elementSize, "fracDipoles");
     fracQuadrupoles = new CudaArray(cu, 6*paddedNumAtoms, elementSize, "fracQuadrupoles");
     field = new CudaArray(cu, 3*paddedNumAtoms, sizeof(long long), "field");
@@ -1433,7 +1439,8 @@ double CudaCalcAmoebaMultipoleForceKernel::execute(ContextImpl& context, bool in
 
     void* computeMomentsArgs[] = {&cu.getPosq().getDevicePointer(), &multipoleParticles->getDevicePointer(),
         &molecularDipoles->getDevicePointer(), &molecularQuadrupoles->getDevicePointer(),
-        &labFrameDipoles->getDevicePointer(), &labFrameQuadrupoles->getDevicePointer()};
+        &labFrameDipoles->getDevicePointer(), &labFrameQuadrupoles->getDevicePointer(),
+        &sphericalDipoles->getDevicePointer(), &sphericalQuadrupoles->getDevicePointer()};
     cu.executeKernel(computeMomentsKernel, computeMomentsArgs, cu.getNumAtoms());
     int startTileIndex = nb.getStartTileIndex();
     int numTileIndices = nb.getNumTiles();

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

@@ -392,6 +392,8 @@ private:
     CudaArray* molecularQuadrupoles;
     CudaArray* labFrameDipoles;
     CudaArray* labFrameQuadrupoles;
+    CudaArray* sphericalDipoles;
+    CudaArray* sphericalQuadrupoles;
     CudaArray* fracDipoles;
     CudaArray* fracQuadrupoles;
     CudaArray* field;

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

 extern "C" __global__ void computeLabFrameMoments(real4* __restrict__ posq, int4* __restrict__ multipoleParticles, float* __restrict__ molecularDipoles,
-        float* __restrict__ molecularQuadrupoles, real* __restrict__ labFrameDipoles, real* __restrict__ labFrameQuadrupoles) {
-    // get coordinates of this atom and the z & x axis atoms
-    // compute the vector between the atoms and 1/sqrt(d2), d2 is distance between
-    // this atom and the axis atom
- 
-    // this atom is referred to as the k-atom in notes below
- 
-    // code common to ZThenX and Bisector
-    
+        float* __restrict__ molecularQuadrupoles, real* __restrict__ labFrameDipoles, real* __restrict__ labFrameQuadrupoles,
+        real* __restrict__ sphericalDipoles, real* __restrict__ sphericalQuadrupoles) {
     for (int atom = blockIdx.x*blockDim.x+threadIdx.x; atom < NUM_ATOMS; atom += gridDim.x*blockDim.x) {
+        // Load the spherical multipoles.
+        
+        int offset = 3*atom;
+        sphericalDipoles[offset+0] = molecularDipoles[offset+2]; // z -> Q_10
+        sphericalDipoles[offset+1] = molecularDipoles[offset+0]; // x -> Q_11c
+        sphericalDipoles[offset+2] = molecularDipoles[offset+1]; // y -> Q_11s
+        offset = 5*atom;
+        sphericalQuadrupoles[offset+0] = -3.0f*(molecularQuadrupoles[offset+0]+molecularQuadrupoles[offset+3]); // zz -> Q_20
+        sphericalQuadrupoles[offset+1] = (2*SQRT(3))*molecularQuadrupoles[offset+2]; // xz -> Q_21c
+        sphericalQuadrupoles[offset+2] = (2*SQRT(3))*molecularQuadrupoles[offset+4]; // yz -> Q_21s
+        sphericalQuadrupoles[offset+3] = SQRT(3)*(molecularQuadrupoles[offset+0]-molecularQuadrupoles[offset+3]); // xx-yy -> Q_22c
+        sphericalQuadrupoles[offset+4] = (2*SQRT(3))*molecularQuadrupoles[offset+1]; // xy -> Q_22s
+        
+        // get coordinates of this atom and the z & x axis atoms
+        // compute the vector between the atoms and 1/sqrt(d2), d2 is distance between
+        // this atom and the axis atom
+
+        // this atom is referred to as the k-atom in notes below
+
+        // code common to ZThenX and Bisector
+        
         int4 particles = multipoleParticles[atom];
         if (particles.x >= 0 && particles.z >= 0) {
             real4 thisParticlePos = posq[atom];
@@ -149,7 +163,7 @@ extern "C" __global__ void computeLabFrameMoments(real4* __restrict__ posq, int4
         
             // Transform the dipole
             
-            unsigned int offset = 3*atom;
+            offset = 3*atom;
             real molDipole[3];
             molDipole[0] = molecularDipoles[offset];
             molDipole[1] = molecularDipoles[offset+1];
@@ -192,6 +206,67 @@ extern "C" __global__ void computeLabFrameMoments(real4* __restrict__ posq, int4
             labFrameQuadrupoles[offset+4] = vectorX.y*(vectorX.z*mPoleXX + vectorY.z*mPoleXY + vectorZ.z*mPoleXZ)
                                         + vectorY.y*(vectorX.z*mPoleXY + vectorY.z*mPoleYY + vectorZ.z*mPoleYZ)
                                         + vectorZ.y*(vectorX.z*mPoleXZ + vectorY.z*mPoleYZ + vectorZ.z*mPoleZZ);
+            
+            // ---------------------------------------------------------------------------------------
+            
+            // Now transform the spherical multipoles.  First do the dipoles.
+
+            offset = 3*atom;
+            real sphericalDipole[3];
+            sphericalDipole[0] = sphericalDipoles[offset];
+            sphericalDipole[1] = sphericalDipoles[offset+1];
+            sphericalDipole[2] = sphericalDipoles[offset+2];
+            if (reverse)
+                sphericalDipole[2] *= -1;
+            sphericalDipoles[offset] = sphericalDipole[0]*vectorZ.z + sphericalDipole[1]*vectorX.z + sphericalDipole[2]*vectorY.z;
+            sphericalDipoles[offset+1] = sphericalDipole[0]*vectorZ.x + sphericalDipole[1]*vectorX.x + sphericalDipole[2]*vectorY.x;
+            sphericalDipoles[offset+2] = sphericalDipole[0]*vectorZ.y + sphericalDipole[1]*vectorX.y + sphericalDipole[2]*vectorY.y;
+            
+            // Now the quadrupoles.
+
+            offset = 5*atom;
+            real sphericalQuadrupole[5];
+            sphericalQuadrupole[0] = sphericalQuadrupoles[offset];
+            sphericalQuadrupole[1] = sphericalQuadrupoles[offset+1];
+            sphericalQuadrupole[2] = sphericalQuadrupoles[offset+2];
+            sphericalQuadrupole[3] = sphericalQuadrupoles[offset+3];
+            sphericalQuadrupole[4] = sphericalQuadrupoles[offset+4];
+            if (reverse) {
+                sphericalQuadrupole[2] *= -1;
+                sphericalQuadrupole[4] *= -1;
+            }
+            real rotatedQuadrupole[5] = {0, 0, 0, 0, 0};
+            real sqrtThree = SQRT(3);
+            rotatedQuadrupole[0] += sphericalQuadrupole[0]*0.5f*(3.0f*vectorZ.z*vectorZ.z - 1.0f);
+            rotatedQuadrupole[1] += sphericalQuadrupole[0]*sqrtThree*vectorZ.z*vectorZ.x;
+            rotatedQuadrupole[2] += sphericalQuadrupole[0]*sqrtThree*vectorZ.z*vectorZ.y;
+            rotatedQuadrupole[3] += sphericalQuadrupole[0]*0.5f*sqrtThree*(vectorZ.x*vectorZ.x - vectorZ.y*vectorZ.y);
+            rotatedQuadrupole[4] += sphericalQuadrupole[0]*sqrtThree*vectorZ.x*vectorZ.y;
+            rotatedQuadrupole[0] += sphericalQuadrupole[1]*sqrtThree*vectorZ.z*vectorX.z;
+            rotatedQuadrupole[1] += sphericalQuadrupole[1]*(vectorZ.x*vectorX.z + vectorZ.z*vectorX.x);
+            rotatedQuadrupole[2] += sphericalQuadrupole[1]*(vectorZ.y*vectorX.z + vectorZ.z*vectorX.y);
+            rotatedQuadrupole[3] += sphericalQuadrupole[1]*(vectorZ.x*vectorX.x - vectorZ.y*vectorX.y);
+            rotatedQuadrupole[4] += sphericalQuadrupole[1]*(vectorZ.y*vectorX.x + vectorZ.x*vectorX.y);
+            rotatedQuadrupole[0] += sphericalQuadrupole[2]*sqrtThree*vectorZ.z*vectorY.z;
+            rotatedQuadrupole[1] += sphericalQuadrupole[2]*(vectorZ.x*vectorY.z + vectorZ.z*vectorY.x);
+            rotatedQuadrupole[2] += sphericalQuadrupole[2]*(vectorZ.y*vectorY.z + vectorZ.z*vectorY.y);
+            rotatedQuadrupole[3] += sphericalQuadrupole[2]*(vectorZ.x*vectorY.x - vectorZ.y*vectorY.y);
+            rotatedQuadrupole[4] += sphericalQuadrupole[2]*(vectorZ.y*vectorY.x + vectorZ.x*vectorY.y);
+            rotatedQuadrupole[0] += sphericalQuadrupole[3]*0.5f*sqrtThree*(vectorX.z*vectorX.z - vectorY.z*vectorY.z);
+            rotatedQuadrupole[1] += sphericalQuadrupole[3]*(vectorX.z*vectorX.x - vectorY.z*vectorY.x);
+            rotatedQuadrupole[2] += sphericalQuadrupole[3]*(vectorX.z*vectorX.y - vectorY.z*vectorY.y);
+            rotatedQuadrupole[3] += sphericalQuadrupole[3]*0.5f*(vectorX.x*vectorX.x - vectorX.y*vectorX.y - vectorY.x*vectorY.x + vectorY.y*vectorY.y);
+            rotatedQuadrupole[4] += sphericalQuadrupole[3]*(vectorX.x*vectorX.y - vectorY.x*vectorY.y);
+            rotatedQuadrupole[0] += sphericalQuadrupole[4]*sqrtThree*vectorX.z*vectorY.z;
+            rotatedQuadrupole[1] += sphericalQuadrupole[4]*(vectorX.x*vectorY.z + vectorX.z*vectorY.x);
+            rotatedQuadrupole[2] += sphericalQuadrupole[4]*(vectorX.y*vectorY.z + vectorX.z*vectorY.y);
+            rotatedQuadrupole[3] += sphericalQuadrupole[4]*(vectorX.x*vectorY.x - vectorX.y*vectorY.y);
+            rotatedQuadrupole[4] += sphericalQuadrupole[4]*(vectorX.y*vectorY.x + vectorX.x*vectorY.y);
+            sphericalQuadrupoles[offset] = rotatedQuadrupole[0];
+            sphericalQuadrupoles[offset+1] = rotatedQuadrupole[1];
+            sphericalQuadrupoles[offset+2] = rotatedQuadrupole[2];
+            sphericalQuadrupoles[offset+3] = rotatedQuadrupole[3];
+            sphericalQuadrupoles[offset+4] = rotatedQuadrupole[4];
         }
         else {
             labFrameDipoles[3*atom] = molecularDipoles[3*atom];

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

@@ -25,9 +25,6 @@
 #include "AmoebaReferenceMultipoleForce.h"
 #include "jama_svd.h"
 #include <algorithm>
-#include <iostream>
-#include <iomanip>
-#define FMT(x) std::setw(16) << std::setprecision(10) << (x)
 
 // In case we're using some primitive version of Visual Studio this will
 // make sure that erf() and erfc() are defined.
@@ -362,6 +359,9 @@ void AmoebaReferenceMultipoleForce::checkChiralCenterAtParticle(MultipoleParticl
     RealOpenMM volume = deltaC.dot(deltaAD);
 
     if (volume < 0.0) {
+        particleI.dipole[1] *= -1.0; // pole(3,i)
+        particleI.quadrupole[QXY] *= -1.0; // pole(6,i) && pole(8,i)
+        particleI.quadrupole[QYZ] *= -1.0; // pole(10,i) && pole(12,i)
         particleI.sphericalDipole[2]     *= -1.0;   // y
         particleI.sphericalQuadrupole[2] *= -1.0;   // yz
         particleI.sphericalQuadrupole[4] *= -1.0;   // xy