Pass more forcefield parameters to the GPU for consistent values

GOMC-WSU · Dec 13, 2024 · f451add · f451add
1 parent 713622a
commit f451add
Show file tree

Hide file tree

Showing 10 changed files with 97 additions and 84 deletions.
diff --git a/src/Ewald.cpp b/src/Ewald.cpp
@@ -1530,8 +1530,7 @@ void Ewald::BoxForceReciprocal(XYZArray const &molCoords,
     CallBoxForceReciprocalGPU(
         ff.particles->getCUDAVars(), atomForceRec, molForceRec, particleCharge,
         particleMol, particleHasNoCharge, particleUsed, startMol, lengthMol,
-        ff.alpha[box], ff.alphaSq[box], constValue, imageSizeRef[box],
-        molCoords, currentAxes, box);
+        constValue, imageSizeRef[box], molCoords, currentAxes, box);
     delete[] particleUsed;
 #else
     // molecule iterator

diff --git a/src/FFParticle.cpp b/src/FFParticle.cpp
@@ -88,9 +88,10 @@ void FFParticle::Init(ff_setup::Particle const &mie,
   double diElectric_1 = 1.0 / forcefield.dielectric;
   InitGPUForceField(*varCUDA, sigmaSq, epsilon_cn, n, forcefield.vdwKind,
                     forcefield.isMartini, count, forcefield.rCut,
-                    forcefield.rCutCoulomb, forcefield.rCutLow,
-                    forcefield.rswitch, forcefield.alpha, forcefield.ewald,
-                    diElectric_1);
+                    forcefield.rCutSq, forcefield.rCutCoulomb,
+                    forcefield.rCutCoulombSq, forcefield.rCutLow,
+                    forcefield.rswitch, forcefield.alpha, forcefield.alphaSq,
+                    forcefield.ewald, diElectric_1);
 #endif
 }
 

diff --git a/src/Forcefield.h b/src/Forcefield.h
@@ -54,12 +54,12 @@ class Forcefield {
   double tolerance;                // Ewald sum terms
   double rswitch;                  // Switch distance
   double dielectric;               // dielectric for martini
-  double scaling_14; //!< Scaling factor for 1-4 pairs' ewald interactions
+  double scaling_14; //!< Scaling factor for 1-4 pairs' Ewald interactions
   double sc_alpha;   // Free energy parameter
   double sc_sigma, sc_sigma_6; // Free energy parameter
 
   bool OneThree, OneFour, OneN; // To include 1-3, 1-4 and more interaction
-  bool electrostatic, ewald;    // To consider columb interaction
+  bool electrostatic, ewald;    // To consider coulomb interaction
   bool vdwGeometricSigma;       // For sigma combining rule
   bool isMartini;
   bool exp6;

diff --git a/src/GPU/CalculateEwaldCUDAKernel.cu b/src/GPU/CalculateEwaldCUDAKernel.cu
@@ -450,8 +450,8 @@ void CallBoxForceReciprocalGPU(
     const std::vector<int> &particleMol,
     const std::vector<bool> &particleHasNoCharge, const bool *particleUsed,
     const std::vector<int> &startMol, const std::vector<int> &lengthMol,
-    double alpha, double alphaSq, double constValue, uint imageSize,
-    XYZArray const &molCoords, BoxDimensions const &boxAxes, int box) {
+    double constValue, uint imageSize, XYZArray const &molCoords,
+    BoxDimensions const &boxAxes, int box) {
   int atomCount = atomForceRec.Count();
   int molCount = molForceRec.Count();
   double *gpu_particleCharge;
@@ -518,13 +518,13 @@ void CallBoxForceReciprocalGPU(
       vars->gpu_aForceRecx, vars->gpu_aForceRecy, vars->gpu_aForceRecz,
       vars->gpu_mForceRecx, vars->gpu_mForceRecy, vars->gpu_mForceRecz,
       gpu_particleCharge, gpu_particleMol, gpu_particleHasNoCharge,
-      gpu_particleUsed, gpu_startMol, gpu_lengthMol, alpha, alphaSq, constValue,
-      imageSize, vars->gpu_kxRef[box], vars->gpu_kyRef[box],
-      vars->gpu_kzRef[box], vars->gpu_x, vars->gpu_y, vars->gpu_z,
-      vars->gpu_prefactRef[box], vars->gpu_sumRnew[box], vars->gpu_sumInew[box],
-      vars->gpu_isFraction, vars->gpu_molIndex, vars->gpu_lambdaCoulomb,
-      vars->gpu_cell_x[box], vars->gpu_cell_y[box], vars->gpu_cell_z[box],
-      vars->gpu_Invcell_x[box], vars->gpu_Invcell_y[box],
+      gpu_particleUsed, gpu_startMol, gpu_lengthMol, vars->gpu_alpha,
+      vars->gpu_alphaSq, constValue, imageSize, vars->gpu_kxRef[box],
+      vars->gpu_kyRef[box], vars->gpu_kzRef[box], vars->gpu_x, vars->gpu_y,
+      vars->gpu_z, vars->gpu_prefactRef[box], vars->gpu_sumRnew[box],
+      vars->gpu_sumInew[box], vars->gpu_isFraction, vars->gpu_molIndex,
+      vars->gpu_lambdaCoulomb, vars->gpu_cell_x[box], vars->gpu_cell_y[box],
+      vars->gpu_cell_z[box], vars->gpu_Invcell_x[box], vars->gpu_Invcell_y[box],
       vars->gpu_Invcell_z[box], vars->gpu_nonOrth, boxAxes.GetAxis(box).x,
       boxAxes.GetAxis(box).y, boxAxes.GetAxis(box).z, box, atomCount);
   cudaDeviceSynchronize();
@@ -558,7 +558,7 @@ __global__ void BoxForceReciprocalGPU(
     double *gpu_mForceRecx, double *gpu_mForceRecy, double *gpu_mForceRecz,
     double *gpu_particleCharge, int *gpu_particleMol,
     bool *gpu_particleHasNoCharge, bool *gpu_particleUsed, int *gpu_startMol,
-    int *gpu_lengthMol, double alpha, double alphaSq, double constValue,
+    int *gpu_lengthMol, double *gpu_alpha, double *gpu_alphaSq, double constValue,
     int imageSize, double *gpu_kx, double *gpu_ky, double *gpu_kz,
     double *gpu_x, double *gpu_y, double *gpu_z, double *gpu_prefact,
     double *gpu_sumRnew, double *gpu_sumInew, bool *gpu_isFraction,
@@ -627,11 +627,11 @@ __global__ void BoxForceReciprocalGPU(
                      gpu_Invcell_z);
         dist = sqrt(distSq);
 
-        double expConstValue = exp(-1.0 * alphaSq * distSq);
+        double expConstValue = exp(-1.0 * gpu_alphaSq[box] * distSq);
         double qiqj = gpu_particleCharge[particleID] *
                       gpu_particleCharge[otherParticle] * qqFactGPU;
         intraForce = qiqj * lambdaCoef * lambdaCoef / distSq;
-        intraForce *= ((erf(alpha * dist) / dist) - constValue * expConstValue);
+        intraForce *= ((erf(gpu_alpha[box] * dist) / dist) - constValue * expConstValue);
         forceX -= intraForce * distVect.x;
         forceY -= intraForce * distVect.y;
         forceZ -= intraForce * distVect.z;

diff --git a/src/GPU/CalculateEwaldCUDAKernel.cuh b/src/GPU/CalculateEwaldCUDAKernel.cuh
@@ -23,8 +23,6 @@ void CallBoxForceReciprocalGPU(VariablesCUDA *vars,
                                const bool *particleUsed,
                                const std::vector<int> &startMol,
                                const std::vector<int> &lengthMol,
-                               double alpha,
-                               double alphaSq,
                                double constValue,
                                uint imageSize,
                                XYZArray const &molCoords,
@@ -103,8 +101,8 @@ __global__ void BoxForceReciprocalGPU(double *gpu_aForceRecx,
                                        bool *gpu_particleUsed,
                                      int *gpu_startMol,
                                       int *gpu_lengthMol,
-                                      double alpha,
-                                      double alphaSq,
+                                      double *gpu_alpha,
+                                      double *gpu_alphaSq,
                                       double constValue,
                                       int imageSize,
                                       double *gpu_kx,