Skip to content

GitLab

Unverified Commit 71e02fa9 authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2130 from peastman/amoeba 

Minor optimizations to AMOEBA
parents 756e479c 99aa1207
Hide whitespace changes
Inline Side-by-side
Showing with 72 additions and 78 deletions
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.cpp
View file @ 71e02fa9
...@@ -1127,7 +1127,6 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const ...@@ -1127,7 +1127,6 @@ void CudaCalcAmoebaMultipoleForceKernel::initialize(const System& system, const
pmePhip.initialize(cu, 10*numMultipoles, elementSize, "pmePhip"); pmePhip.initialize(cu, 10*numMultipoles, elementSize, "pmePhip");
pmePhidp.initialize(cu, 20*numMultipoles, elementSize, "pmePhidp"); pmePhidp.initialize(cu, 20*numMultipoles, elementSize, "pmePhidp");
pmeCphi.initialize(cu, 10*numMultipoles, elementSize, "pmeCphi"); pmeCphi.initialize(cu, 10*numMultipoles, elementSize, "pmeCphi");
pmeAtomRange.initialize<int>(cu, gridSizeX*gridSizeY*gridSizeZ+1, "pmeAtomRange");
sort = new CudaSort(cu, new SortTrait(), cu.getNumAtoms()); sort = new CudaSort(cu, new SortTrait(), cu.getNumAtoms());
cufftResult result = cufftPlan3d(&fft, gridSizeX, gridSizeY, gridSizeZ, cu.getUseDoublePrecision() ? CUFFT_Z2Z : CUFFT_C2C); cufftResult result = cufftPlan3d(&fft, gridSizeX, gridSizeY, gridSizeZ, cu.getUseDoublePrecision() ? CUFFT_Z2Z : CUFFT_C2C);
if (result != CUFFT_SUCCESS) if (result != CUFFT_SUCCESS)
......
plugins/amoeba/platforms/cuda/src/AmoebaCudaKernels.h
View file @ 71e02fa9
...@@ -457,7 +457,6 @@ private: ...@@ -457,7 +457,6 @@ private:
CudaArray pmePhip; CudaArray pmePhip;
CudaArray pmePhidp; CudaArray pmePhidp;
CudaArray pmeCphi; CudaArray pmeCphi;
CudaArray pmeAtomRange;
CudaArray lastPositions; CudaArray lastPositions;
CudaSort* sort; CudaSort* sort;
cufftHandle fft; cufftHandle fft;
......
plugins/amoeba/platforms/cuda/src/kernels/multipoleInducedField.cu
View file @ 71e02fa9
...@@ -193,39 +193,39 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, real3 de ...@@ -193,39 +193,39 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, real3 de
real3 dipole = atom1.inducedDipole; real3 dipole = atom1.inducedDipole;
real muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z; real muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z;
atom2.fieldGradient[0] -= muDotR*deltaR.x*deltaR.x*scale3 - (2*dipole.x*deltaR.x + muDotR)*scale2; atom2.fieldGradient[0] -= (muDotR*scale3)*deltaR.x*deltaR.x - (2*dipole.x*deltaR.x + muDotR)*scale2;
atom2.fieldGradient[1] -= muDotR*deltaR.y*deltaR.y*scale3 - (2*dipole.y*deltaR.y + muDotR)*scale2; atom2.fieldGradient[1] -= (muDotR*scale3)*deltaR.y*deltaR.y - (2*dipole.y*deltaR.y + muDotR)*scale2;
atom2.fieldGradient[2] -= muDotR*deltaR.z*deltaR.z*scale3 - (2*dipole.z*deltaR.z + muDotR)*scale2; atom2.fieldGradient[2] -= (muDotR*scale3)*deltaR.z*deltaR.z - (2*dipole.z*deltaR.z + muDotR)*scale2;
atom2.fieldGradient[3] -= muDotR*deltaR.x*deltaR.y*scale3 - (dipole.x*deltaR.y + dipole.y*deltaR.x)*scale2; atom2.fieldGradient[3] -= (muDotR*scale3)*deltaR.x*deltaR.y - (dipole.x*deltaR.y + dipole.y*deltaR.x)*scale2;
atom2.fieldGradient[4] -= muDotR*deltaR.x*deltaR.z*scale3 - (dipole.x*deltaR.z + dipole.z*deltaR.x)*scale2; atom2.fieldGradient[4] -= (muDotR*scale3)*deltaR.x*deltaR.z - (dipole.x*deltaR.z + dipole.z*deltaR.x)*scale2;
atom2.fieldGradient[5] -= muDotR*deltaR.y*deltaR.z*scale3 - (dipole.y*deltaR.z + dipole.z*deltaR.y)*scale2; atom2.fieldGradient[5] -= (muDotR*scale3)*deltaR.y*deltaR.z - (dipole.y*deltaR.z + dipole.z*deltaR.y)*scale2;
dipole = atom1.inducedDipolePolar; dipole = atom1.inducedDipolePolar;
muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z; muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z;
atom2.fieldGradientPolar[0] -= muDotR*deltaR.x*deltaR.x*scale3 - (2*dipole.x*deltaR.x + muDotR)*scale2; atom2.fieldGradientPolar[0] -= (muDotR*scale3)*deltaR.x*deltaR.x - (2*dipole.x*deltaR.x + muDotR)*scale2;
atom2.fieldGradientPolar[1] -= muDotR*deltaR.y*deltaR.y*scale3 - (2*dipole.y*deltaR.y + muDotR)*scale2; atom2.fieldGradientPolar[1] -= (muDotR*scale3)*deltaR.y*deltaR.y - (2*dipole.y*deltaR.y + muDotR)*scale2;
atom2.fieldGradientPolar[2] -= muDotR*deltaR.z*deltaR.z*scale3 - (2*dipole.z*deltaR.z + muDotR)*scale2; atom2.fieldGradientPolar[2] -= (muDotR*scale3)*deltaR.z*deltaR.z - (2*dipole.z*deltaR.z + muDotR)*scale2;
atom2.fieldGradientPolar[3] -= muDotR*deltaR.x*deltaR.y*scale3 - (dipole.x*deltaR.y + dipole.y*deltaR.x)*scale2; atom2.fieldGradientPolar[3] -= (muDotR*scale3)*deltaR.x*deltaR.y - (dipole.x*deltaR.y + dipole.y*deltaR.x)*scale2;
atom2.fieldGradientPolar[4] -= muDotR*deltaR.x*deltaR.z*scale3 - (dipole.x*deltaR.z + dipole.z*deltaR.x)*scale2; atom2.fieldGradientPolar[4] -= (muDotR*scale3)*deltaR.x*deltaR.z - (dipole.x*deltaR.z + dipole.z*deltaR.x)*scale2;
atom2.fieldGradientPolar[5] -= muDotR*deltaR.y*deltaR.z*scale3 - (dipole.y*deltaR.z + dipole.z*deltaR.y)*scale2; atom2.fieldGradientPolar[5] -= (muDotR*scale3)*deltaR.y*deltaR.z - (dipole.y*deltaR.z + dipole.z*deltaR.y)*scale2;
dipole = atom2.inducedDipole; dipole = atom2.inducedDipole;
muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z; muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z;
atom1.fieldGradient[0] += muDotR*deltaR.x*deltaR.x*scale3 - (2*dipole.x*deltaR.x + muDotR)*scale2; atom1.fieldGradient[0] += (muDotR*scale3)*deltaR.x*deltaR.x - (2*dipole.x*deltaR.x + muDotR)*scale2;
atom1.fieldGradient[1] += muDotR*deltaR.y*deltaR.y*scale3 - (2*dipole.y*deltaR.y + muDotR)*scale2; atom1.fieldGradient[1] += (muDotR*scale3)*deltaR.y*deltaR.y - (2*dipole.y*deltaR.y + muDotR)*scale2;
atom1.fieldGradient[2] += muDotR*deltaR.z*deltaR.z*scale3 - (2*dipole.z*deltaR.z + muDotR)*scale2; atom1.fieldGradient[2] += (muDotR*scale3)*deltaR.z*deltaR.z - (2*dipole.z*deltaR.z + muDotR)*scale2;
atom1.fieldGradient[3] += muDotR*deltaR.x*deltaR.y*scale3 - (dipole.x*deltaR.y + dipole.y*deltaR.x)*scale2; atom1.fieldGradient[3] += (muDotR*scale3)*deltaR.x*deltaR.y - (dipole.x*deltaR.y + dipole.y*deltaR.x)*scale2;
atom1.fieldGradient[4] += muDotR*deltaR.x*deltaR.z*scale3 - (dipole.x*deltaR.z + dipole.z*deltaR.x)*scale2; atom1.fieldGradient[4] += (muDotR*scale3)*deltaR.x*deltaR.z - (dipole.x*deltaR.z + dipole.z*deltaR.x)*scale2;
atom1.fieldGradient[5] += muDotR*deltaR.y*deltaR.z*scale3 - (dipole.y*deltaR.z + dipole.z*deltaR.y)*scale2; atom1.fieldGradient[5] += (muDotR*scale3)*deltaR.y*deltaR.z - (dipole.y*deltaR.z + dipole.z*deltaR.y)*scale2;
dipole = atom2.inducedDipolePolar; dipole = atom2.inducedDipolePolar;
muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z; muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z;
atom1.fieldGradientPolar[0] += muDotR*deltaR.x*deltaR.x*scale3 - (2*dipole.x*deltaR.x + muDotR)*scale2; atom1.fieldGradientPolar[0] += (muDotR*scale3)*deltaR.x*deltaR.x - (2*dipole.x*deltaR.x + muDotR)*scale2;
atom1.fieldGradientPolar[1] += muDotR*deltaR.y*deltaR.y*scale3 - (2*dipole.y*deltaR.y + muDotR)*scale2; atom1.fieldGradientPolar[1] += (muDotR*scale3)*deltaR.y*deltaR.y - (2*dipole.y*deltaR.y + muDotR)*scale2;
atom1.fieldGradientPolar[2] += muDotR*deltaR.z*deltaR.z*scale3 - (2*dipole.z*deltaR.z + muDotR)*scale2; atom1.fieldGradientPolar[2] += (muDotR*scale3)*deltaR.z*deltaR.z - (2*dipole.z*deltaR.z + muDotR)*scale2;
atom1.fieldGradientPolar[3] += muDotR*deltaR.x*deltaR.y*scale3 - (dipole.x*deltaR.y + dipole.y*deltaR.x)*scale2; atom1.fieldGradientPolar[3] += (muDotR*scale3)*deltaR.x*deltaR.y - (dipole.x*deltaR.y + dipole.y*deltaR.x)*scale2;
atom1.fieldGradientPolar[4] += muDotR*deltaR.x*deltaR.z*scale3 - (dipole.x*deltaR.z + dipole.z*deltaR.x)*scale2; atom1.fieldGradientPolar[4] += (muDotR*scale3)*deltaR.x*deltaR.z - (dipole.x*deltaR.z + dipole.z*deltaR.x)*scale2;
atom1.fieldGradientPolar[5] += muDotR*deltaR.y*deltaR.z*scale3 - (dipole.y*deltaR.z + dipole.z*deltaR.y)*scale2; atom1.fieldGradientPolar[5] += (muDotR*scale3)*deltaR.y*deltaR.z - (dipole.y*deltaR.z + dipole.z*deltaR.y)*scale2;
#endif #endif
} }
#elif defined USE_GK #elif defined USE_GK
...@@ -315,39 +315,39 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, real3 de ...@@ -315,39 +315,39 @@ __device__ void computeOneInteraction(AtomData& atom1, AtomData& atom2, real3 de
real3 dipole = atom1.inducedDipole; real3 dipole = atom1.inducedDipole;
real muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z; real muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z;
atom2.fieldGradient[0] -= muDotR*deltaR.x*deltaR.x*rr7 - (2*dipole.x*deltaR.x + muDotR)*rr5; atom2.fieldGradient[0] -= (muDotR*rr7)*deltaR.x*deltaR.x - (2*dipole.x*deltaR.x + muDotR)*rr5;
atom2.fieldGradient[1] -= muDotR*deltaR.y*deltaR.y*rr7 - (2*dipole.y*deltaR.y + muDotR)*rr5; atom2.fieldGradient[1] -= (muDotR*rr7)*deltaR.y*deltaR.y - (2*dipole.y*deltaR.y + muDotR)*rr5;
atom2.fieldGradient[2] -= muDotR*deltaR.z*deltaR.z*rr7 - (2*dipole.z*deltaR.z + muDotR)*rr5; atom2.fieldGradient[2] -= (muDotR*rr7)*deltaR.z*deltaR.z - (2*dipole.z*deltaR.z + muDotR)*rr5;
atom2.fieldGradient[3] -= muDotR*deltaR.x*deltaR.y*rr7 - (dipole.x*deltaR.y + dipole.y*deltaR.x)*rr5; atom2.fieldGradient[3] -= (muDotR*rr7)*deltaR.x*deltaR.y - (dipole.x*deltaR.y + dipole.y*deltaR.x)*rr5;
atom2.fieldGradient[4] -= muDotR*deltaR.x*deltaR.z*rr7 - (dipole.x*deltaR.z + dipole.z*deltaR.x)*rr5; atom2.fieldGradient[4] -= (muDotR*rr7)*deltaR.x*deltaR.z - (dipole.x*deltaR.z + dipole.z*deltaR.x)*rr5;
atom2.fieldGradient[5] -= muDotR*deltaR.y*deltaR.z*rr7 - (dipole.y*deltaR.z + dipole.z*deltaR.y)*rr5; atom2.fieldGradient[5] -= (muDotR*rr7)*deltaR.y*deltaR.z - (dipole.y*deltaR.z + dipole.z*deltaR.y)*rr5;
dipole = atom1.inducedDipolePolar; dipole = atom1.inducedDipolePolar;
muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z; muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z;
atom2.fieldGradientPolar[0] -= muDotR*deltaR.x*deltaR.x*rr7 - (2*dipole.x*deltaR.x + muDotR)*rr5; atom2.fieldGradientPolar[0] -= (muDotR*rr7)*deltaR.x*deltaR.x - (2*dipole.x*deltaR.x + muDotR)*rr5;
atom2.fieldGradientPolar[1] -= muDotR*deltaR.y*deltaR.y*rr7 - (2*dipole.y*deltaR.y + muDotR)*rr5; atom2.fieldGradientPolar[1] -= (muDotR*rr7)*deltaR.y*deltaR.y - (2*dipole.y*deltaR.y + muDotR)*rr5;
atom2.fieldGradientPolar[2] -= muDotR*deltaR.z*deltaR.z*rr7 - (2*dipole.z*deltaR.z + muDotR)*rr5; atom2.fieldGradientPolar[2] -= (muDotR*rr7)*deltaR.z*deltaR.z - (2*dipole.z*deltaR.z + muDotR)*rr5;
atom2.fieldGradientPolar[3] -= muDotR*deltaR.x*deltaR.y*rr7 - (dipole.x*deltaR.y + dipole.y*deltaR.x)*rr5; atom2.fieldGradientPolar[3] -= (muDotR*rr7)*deltaR.x*deltaR.y - (dipole.x*deltaR.y + dipole.y*deltaR.x)*rr5;
atom2.fieldGradientPolar[4] -= muDotR*deltaR.x*deltaR.z*rr7 - (dipole.x*deltaR.z + dipole.z*deltaR.x)*rr5; atom2.fieldGradientPolar[4] -= (muDotR*rr7)*deltaR.x*deltaR.z - (dipole.x*deltaR.z + dipole.z*deltaR.x)*rr5;
atom2.fieldGradientPolar[5] -= muDotR*deltaR.y*deltaR.z*rr7 - (dipole.y*deltaR.z + dipole.z*deltaR.y)*rr5; atom2.fieldGradientPolar[5] -= (muDotR*rr7)*deltaR.y*deltaR.z - (dipole.y*deltaR.z + dipole.z*deltaR.y)*rr5;
dipole = atom2.inducedDipole; dipole = atom2.inducedDipole;
muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z; muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z;
atom1.fieldGradient[0] += muDotR*deltaR.x*deltaR.x*rr7 - (2*dipole.x*deltaR.x + muDotR)*rr5; atom1.fieldGradient[0] += (muDotR*rr7)*deltaR.x*deltaR.x - (2*dipole.x*deltaR.x + muDotR)*rr5;
atom1.fieldGradient[1] += muDotR*deltaR.y*deltaR.y*rr7 - (2*dipole.y*deltaR.y + muDotR)*rr5; atom1.fieldGradient[1] += (muDotR*rr7)*deltaR.y*deltaR.y - (2*dipole.y*deltaR.y + muDotR)*rr5;
atom1.fieldGradient[2] += muDotR*deltaR.z*deltaR.z*rr7 - (2*dipole.z*deltaR.z + muDotR)*rr5; atom1.fieldGradient[2] += (muDotR*rr7)*deltaR.z*deltaR.z - (2*dipole.z*deltaR.z + muDotR)*rr5;
atom1.fieldGradient[3] += muDotR*deltaR.x*deltaR.y*rr7 - (dipole.x*deltaR.y + dipole.y*deltaR.x)*rr5; atom1.fieldGradient[3] += (muDotR*rr7)*deltaR.x*deltaR.y - (dipole.x*deltaR.y + dipole.y*deltaR.x)*rr5;
atom1.fieldGradient[4] += muDotR*deltaR.x*deltaR.z*rr7 - (dipole.x*deltaR.z + dipole.z*deltaR.x)*rr5; atom1.fieldGradient[4] += (muDotR*rr7)*deltaR.x*deltaR.z - (dipole.x*deltaR.z + dipole.z*deltaR.x)*rr5;
atom1.fieldGradient[5] += muDotR*deltaR.y*deltaR.z*rr7 - (dipole.y*deltaR.z + dipole.z*deltaR.y)*rr5; atom1.fieldGradient[5] += (muDotR*rr7)*deltaR.y*deltaR.z - (dipole.y*deltaR.z + dipole.z*deltaR.y)*rr5;
dipole = atom2.inducedDipolePolar; dipole = atom2.inducedDipolePolar;
muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z; muDotR = dipole.x*deltaR.x + dipole.y*deltaR.y + dipole.z*deltaR.z;
atom1.fieldGradientPolar[0] += muDotR*deltaR.x*deltaR.x*rr7 - (2*dipole.x*deltaR.x + muDotR)*rr5; atom1.fieldGradientPolar[0] += (muDotR*rr7)*deltaR.x*deltaR.x - (2*dipole.x*deltaR.x + muDotR)*rr5;
atom1.fieldGradientPolar[1] += muDotR*deltaR.y*deltaR.y*rr7 - (2*dipole.y*deltaR.y + muDotR)*rr5; atom1.fieldGradientPolar[1] += (muDotR*rr7)*deltaR.y*deltaR.y - (2*dipole.y*deltaR.y + muDotR)*rr5;
atom1.fieldGradientPolar[2] += muDotR*deltaR.z*deltaR.z*rr7 - (2*dipole.z*deltaR.z + muDotR)*rr5; atom1.fieldGradientPolar[2] += (muDotR*rr7)*deltaR.z*deltaR.z - (2*dipole.z*deltaR.z + muDotR)*rr5;
atom1.fieldGradientPolar[3] += muDotR*deltaR.x*deltaR.y*rr7 - (dipole.x*deltaR.y + dipole.y*deltaR.x)*rr5; atom1.fieldGradientPolar[3] += (muDotR*rr7)*deltaR.x*deltaR.y - (dipole.x*deltaR.y + dipole.y*deltaR.x)*rr5;
atom1.fieldGradientPolar[4] += muDotR*deltaR.x*deltaR.z*rr7 - (dipole.x*deltaR.z + dipole.z*deltaR.x)*rr5; atom1.fieldGradientPolar[4] += (muDotR*rr7)*deltaR.x*deltaR.z - (dipole.x*deltaR.z + dipole.z*deltaR.x)*rr5;
atom1.fieldGradientPolar[5] += muDotR*deltaR.y*deltaR.z*rr7 - (dipole.y*deltaR.z + dipole.z*deltaR.y)*rr5; atom1.fieldGradientPolar[5] += (muDotR*rr7)*deltaR.y*deltaR.z - (dipole.y*deltaR.z + dipole.z*deltaR.y)*rr5;
#endif #endif
} }
#endif #endif
......
plugins/amoeba/platforms/cuda/src/kernels/multipolePme.cu
View file @ 71e02fa9
...@@ -235,17 +235,16 @@ extern "C" __global__ void gridSpreadFixedMultipoles(const real4* __restrict__ p ...@@ -235,17 +235,16 @@ extern "C" __global__ void gridSpreadFixedMultipoles(const real4* __restrict__ p
ybase -= (ybase >= GRID_SIZE_Y ? GRID_SIZE_Y : 0); ybase -= (ybase >= GRID_SIZE_Y ? GRID_SIZE_Y : 0);
ybase = xbase + ybase*GRID_SIZE_Z; ybase = xbase + ybase*GRID_SIZE_Z;
real4 u = theta2[iy]; real4 u = theta2[iy];
real term0 = atomCharge*t.x*u.x + atomDipoleY*t.x*u.y + atomQuadrupoleYY*t.x*u.z + atomDipoleX*t.y*u.x + atomQuadrupoleXY*t.y*u.y + atomQuadrupoleXX*t.z*u.x;
real term1 = atomDipoleZ*t.x*u.x + atomQuadrupoleYZ*t.x*u.y + atomQuadrupoleXZ*t.y*u.x;
real term2 = atomQuadrupoleZZ*t.x*u.x;
for (int iz = 0; iz < PME_ORDER; iz++) { for (int iz = 0; iz < PME_ORDER; iz++) {
int zindex = igrid3+iz; int zindex = igrid3+iz;
zindex -= (zindex >= GRID_SIZE_Z ? GRID_SIZE_Z : 0); zindex -= (zindex >= GRID_SIZE_Z ? GRID_SIZE_Z : 0);
int index = ybase + zindex; int index = ybase + zindex;
real4 v = theta3[iz]; real4 v = theta3[iz];
real add = term0*v.x + term1*v.y + term2*v.z;
real term0 = atomCharge*u.x*v.x + atomDipoleY*u.y*v.x + atomDipoleZ*u.x*v.y + atomQuadrupoleYY*u.z*v.x + atomQuadrupoleZZ*u.x*v.z + atomQuadrupoleYZ*u.y*v.y;
real term1 = atomDipoleX*u.x*v.x + atomQuadrupoleXY*u.y*v.x + atomQuadrupoleXZ*u.x*v.y;
real term2 = atomQuadrupoleXX * u.x * v.x;
real add = term0*t.x + term1*t.y + term2*t.z;
#ifdef USE_DOUBLE_PRECISION #ifdef USE_DOUBLE_PRECISION
unsigned long long * ulonglong_p = (unsigned long long *) pmeGrid; unsigned long long * ulonglong_p = (unsigned long long *) pmeGrid;
atomicAdd(&ulonglong_p[2*index], static_cast<unsigned long long>((long long) (add*0x100000000))); atomicAdd(&ulonglong_p[2*index], static_cast<unsigned long long>((long long) (add*0x100000000)));
...@@ -337,6 +336,10 @@ extern "C" __global__ void gridSpreadInducedDipoles(const real4* __restrict__ po ...@@ -337,6 +336,10 @@ extern "C" __global__ void gridSpreadInducedDipoles(const real4* __restrict__ po
ybase -= (ybase >= GRID_SIZE_Y ? GRID_SIZE_Y : 0); ybase -= (ybase >= GRID_SIZE_Y ? GRID_SIZE_Y : 0);
ybase = xbase + ybase*GRID_SIZE_Z; ybase = xbase + ybase*GRID_SIZE_Z;
real4 u = theta2[iy]; real4 u = theta2[iy];
real term01 = finducedDipole.y*t.x*u.y + finducedDipole.x*t.y*u.x;
real term11 = finducedDipole.z*t.x*u.x;
real term02 = finducedDipolePolar.y*t.x*u.y + finducedDipolePolar.x*t.y*u.x;
real term12 = finducedDipolePolar.z*t.x*u.x;
for (int iz = 0; iz < PME_ORDER; iz++) { for (int iz = 0; iz < PME_ORDER; iz++) {
int zindex = igrid3+iz; int zindex = igrid3+iz;
...@@ -344,12 +347,8 @@ extern "C" __global__ void gridSpreadInducedDipoles(const real4* __restrict__ po ...@@ -344,12 +347,8 @@ extern "C" __global__ void gridSpreadInducedDipoles(const real4* __restrict__ po
int index = ybase + zindex; int index = ybase + zindex;
real4 v = theta3[iz]; real4 v = theta3[iz];
real term01 = finducedDipole.y*u.y*v.x + finducedDipole.z*u.x*v.y; real add1 = term01*v.x + term11*v.y;
real term11 = finducedDipole.x*u.x*v.x; real add2 = term02*v.x + term12*v.y;
real term02 = finducedDipolePolar.y*u.y*v.x + finducedDipolePolar.z*u.x*v.y;
real term12 = finducedDipolePolar.x*u.x*v.x;
real add1 = term01*t.x + term11*t.y;
real add2 = term02*t.x + term12*t.y;
#ifdef USE_DOUBLE_PRECISION #ifdef USE_DOUBLE_PRECISION
unsigned long long * ulonglong_p = (unsigned long long *) pmeGrid; unsigned long long * ulonglong_p = (unsigned long long *) pmeGrid;
atomicAdd(&ulonglong_p[2*index], static_cast<unsigned long long>((long long) (add1*0x100000000))); atomicAdd(&ulonglong_p[2*index], static_cast<unsigned long long>((long long) (add1*0x100000000)));
......
plugins/amoeba/platforms/reference/src/SimTKReference/AmoebaReferenceMultipoleForce.cpp
View file @ 71e02fa9
...@@ -5477,18 +5477,18 @@ void AmoebaReferencePmeMultipoleForce::spreadFixedMultipolesOntoGrid(const vecto ...@@ -5477,18 +5477,18 @@ void AmoebaReferencePmeMultipoleForce::spreadFixedMultipolesOntoGrid(const vecto
IntVec& gridPoint = _iGrid[atomIndex]; IntVec& gridPoint = _iGrid[atomIndex];
for (int ix = 0; ix < AMOEBA_PME_ORDER; ix++) { for (int ix = 0; ix < AMOEBA_PME_ORDER; ix++) {
int x = (gridPoint[0]+ix) % _pmeGridDimensions[0]; int x = (gridPoint[0]+ix) % _pmeGridDimensions[0];
double4 t = _thetai[0][atomIndex*AMOEBA_PME_ORDER+ix];
for (int iy = 0; iy < AMOEBA_PME_ORDER; iy++) { for (int iy = 0; iy < AMOEBA_PME_ORDER; iy++) {
int y = (gridPoint[1]+iy) % _pmeGridDimensions[1]; int y = (gridPoint[1]+iy) % _pmeGridDimensions[1];
double4 u = _thetai[1][atomIndex*AMOEBA_PME_ORDER+iy];
double term0 = atomCharge*t[0]*u[0] + atomDipole[1]*t[0]*u[1] + atomQuadrupoleYY*t[0]*u[2] + atomDipole[0]*t[1]*u[0] + atomQuadrupoleXY*t[1]*u[1] + atomQuadrupoleXX*t[2]*u[0];
double term1 = atomDipole[2]*t[0]*u[0] + atomQuadrupoleYZ*t[0]*u[1] + atomQuadrupoleXZ*t[1]*u[0];
double term2 = atomQuadrupoleZZ*t[0]*u[0];
for (int iz = 0; iz < AMOEBA_PME_ORDER; iz++) { for (int iz = 0; iz < AMOEBA_PME_ORDER; iz++) {
int z = (gridPoint[2]+iz) % _pmeGridDimensions[2]; int z = (gridPoint[2]+iz) % _pmeGridDimensions[2];
double4 t = _thetai[0][atomIndex*AMOEBA_PME_ORDER+ix];
double4 u = _thetai[1][atomIndex*AMOEBA_PME_ORDER+iy];
double4 v = _thetai[2][atomIndex*AMOEBA_PME_ORDER+iz]; double4 v = _thetai[2][atomIndex*AMOEBA_PME_ORDER+iz];
double term0 = atomCharge*u[0]*v[0] + atomDipole[1]*u[1]*v[0] + atomDipole[2]*u[0]*v[1] + atomQuadrupoleYY*u[2]*v[0] + atomQuadrupoleZZ*u[0]*v[2] + atomQuadrupoleYZ*u[1]*v[1];
double term1 = atomDipole[0]*u[0]*v[0] + atomQuadrupoleXY*u[1]*v[0] + atomQuadrupoleXZ*u[0]*v[1];
double term2 = atomQuadrupoleXX * u[0] * v[0];
t_complex& gridValue = _pmeGrid[x*_pmeGridDimensions[1]*_pmeGridDimensions[2]+y*_pmeGridDimensions[2]+z]; t_complex& gridValue = _pmeGrid[x*_pmeGridDimensions[1]*_pmeGridDimensions[2]+y*_pmeGridDimensions[2]+z];
gridValue.re += term0*t[0] + term1*t[1] + term2*t[2]; gridValue.re += term0*v[0] + term1*v[1] + term2*v[2];
} }
} }
} }
...@@ -5668,23 +5668,20 @@ void AmoebaReferencePmeMultipoleForce::spreadInducedDipolesOnGrid(const vector<V ...@@ -5668,23 +5668,20 @@ void AmoebaReferencePmeMultipoleForce::spreadInducedDipolesOnGrid(const vector<V
IntVec& gridPoint = _iGrid[atomIndex]; IntVec& gridPoint = _iGrid[atomIndex];
for (int ix = 0; ix < AMOEBA_PME_ORDER; ix++) { for (int ix = 0; ix < AMOEBA_PME_ORDER; ix++) {
int x = (gridPoint[0]+ix) % _pmeGridDimensions[0]; int x = (gridPoint[0]+ix) % _pmeGridDimensions[0];
double4 t = _thetai[0][atomIndex*AMOEBA_PME_ORDER+ix];
for (int iy = 0; iy < AMOEBA_PME_ORDER; iy++) { for (int iy = 0; iy < AMOEBA_PME_ORDER; iy++) {
int y = (gridPoint[1]+iy) % _pmeGridDimensions[1]; int y = (gridPoint[1]+iy) % _pmeGridDimensions[1];
double4 u = _thetai[1][atomIndex*AMOEBA_PME_ORDER+iy];
double term01 = inducedDipole[1]*t[0]*u[1] + inducedDipole[0]*t[1]*u[0];
double term11 = inducedDipole[2]*t[0]*u[0];
double term02 = inducedDipolePolar[1]*t[0]*u[1] + inducedDipolePolar[0]*t[1]*u[0];
double term12 = inducedDipolePolar[2]*t[0]*u[0];
for (int iz = 0; iz < AMOEBA_PME_ORDER; iz++) { for (int iz = 0; iz < AMOEBA_PME_ORDER; iz++) {
int z = (gridPoint[2]+iz) % _pmeGridDimensions[2]; int z = (gridPoint[2]+iz) % _pmeGridDimensions[2];
double4 t = _thetai[0][atomIndex*AMOEBA_PME_ORDER+ix];
double4 u = _thetai[1][atomIndex*AMOEBA_PME_ORDER+iy];
double4 v = _thetai[2][atomIndex*AMOEBA_PME_ORDER+iz]; double4 v = _thetai[2][atomIndex*AMOEBA_PME_ORDER+iz];
double term01 = inducedDipole[1]*u[1]*v[0] + inducedDipole[2]*u[0]*v[1];
double term11 = inducedDipole[0]*u[0]*v[0];
double term02 = inducedDipolePolar[1]*u[1]*v[0] + inducedDipolePolar[2]*u[0]*v[1];
double term12 = inducedDipolePolar[0]*u[0]*v[0];
t_complex& gridValue = _pmeGrid[x*_pmeGridDimensions[1]*_pmeGridDimensions[2]+y*_pmeGridDimensions[2]+z]; t_complex& gridValue = _pmeGrid[x*_pmeGridDimensions[1]*_pmeGridDimensions[2]+y*_pmeGridDimensions[2]+z];
gridValue.re += term01*t[0] + term11*t[1]; gridValue.re += term01*v[0] + term11*v[1];
gridValue.im += term02*t[0] + term12*t[1]; gridValue.im += term02*v[0] + term12*v[1];
} }
} }
} }
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment