lj speed up

Author: brucefan1983

src/force/lj.cu

@@ -26,6 +26,8 @@ The class dealing with the Lennard-Jones (LJ) pairwise potentials.
 // best block size here: 128
 #define BLOCK_SIZE_FORCE 128
 
+#define USE_FIXED_NEIGHBOR
+
 LJ::LJ(FILE* fid, int num_types, int num_atoms)
 {
   printf("Use %d-element LJ potential with elements:\n", num_types);
@@ -69,11 +71,11 @@ LJ::~LJ(void)
 
 // get U_ij and (d U_ij / d r_ij) / r_ij (the LJ potential)
 static __device__ void
-find_p2_and_f2(double s6e4, double s12e4, double d12sq, double& p2, double& f2)
+find_p2_and_f2(float s6e4, float s12e4, float d12sq, float& p2, float& f2)
 {
-  double d12inv2 = 1.0 / d12sq;
-  double d12inv6 = d12inv2 * d12inv2 * d12inv2;
-  f2 = 6.0 * (s6e4 * d12inv6 - s12e4 * 2.0 * d12inv6 * d12inv6) * d12inv2;
+  float d12inv2 = 1.0f / d12sq;
+  float d12inv6 = d12inv2 * d12inv2 * d12inv2;
+  f2 = 6.0f * (s6e4 * d12inv6 - s12e4 * 2.0f * d12inv6 * d12inv6) * d12inv2;
   p2 = s12e4 * d12inv6 * d12inv6 - s6e4 * d12inv6;
 }
 
@@ -97,19 +99,19 @@ static __global__ void gpu_find_force(
   double* g_potential)
 {
   int n1 = blockIdx.x * blockDim.x + threadIdx.x + N1; // particle index
-  double s_fx = 0.0;                                   // force_x
-  double s_fy = 0.0;                                   // force_y
-  double s_fz = 0.0;                                   // force_z
-  double s_pe = 0.0;                                   // potential energy
-  double s_sxx = 0.0;                                  // virial_stress_xx
-  double s_sxy = 0.0;                                  // virial_stress_xy
-  double s_sxz = 0.0;                                  // virial_stress_xz
-  double s_syx = 0.0;                                  // virial_stress_yx
-  double s_syy = 0.0;                                  // virial_stress_yy
-  double s_syz = 0.0;                                  // virial_stress_yz
-  double s_szx = 0.0;                                  // virial_stress_zx
-  double s_szy = 0.0;                                  // virial_stress_zy
-  double s_szz = 0.0;                                  // virial_stress_zz
+  float s_fx = 0.0f;                                   // force_x
+  float s_fy = 0.0f;                                   // force_y
+  float s_fz = 0.0f;                                   // force_z
+  float s_pe = 0.0f;                                   // potential energy
+  float s_sxx = 0.0f;                                  // virial_stress_xx
+  float s_sxy = 0.0f;                                  // virial_stress_xy
+  float s_sxz = 0.0f;                                  // virial_stress_xz
+  float s_syx = 0.0f;                                  // virial_stress_yx
+  float s_syy = 0.0f;                                  // virial_stress_yy
+  float s_syz = 0.0f;                                  // virial_stress_yz
+  float s_szx = 0.0f;                                  // virial_stress_zx
+  float s_szy = 0.0f;                                  // virial_stress_zy
+  float s_szz = 0.0f;                                  // virial_stress_zz
 
   if (n1 < N2) {
     int neighbor_number = g_neighbor_number[n1];
@@ -122,33 +124,33 @@ static __global__ void gpu_find_force(
       int n2 = g_neighbor_list[n1 + number_of_particles * i1];
       int type2 = g_type[n2];
 
-      double x12 = g_x[n2] - x1;
-      double y12 = g_y[n2] - y1;
-      double z12 = g_z[n2] - z1;
+      float x12 = g_x[n2] - x1;
+      float y12 = g_y[n2] - y1;
+      float z12 = g_z[n2] - z1;
       apply_mic(box, x12, y12, z12);
-      double d12sq = x12 * x12 + y12 * y12 + z12 * z12;
+      float d12sq = x12 * x12 + y12 * y12 + z12 * z12;
 
-      double p2, f2;
+      float p2, f2;
       if (d12sq >= lj.cutoff_square[type1][type2]) {
         continue;
       }
       find_p2_and_f2(lj.s6e4[type1][type2], lj.s12e4[type1][type2], d12sq, p2, f2);
 
       // treat two-body potential in the same way as many-body potential
-      double f12x = f2 * x12 * 0.5;
-      double f12y = f2 * y12 * 0.5;
-      double f12z = f2 * z12 * 0.5;
-      double f21x = -f12x;
-      double f21y = -f12y;
-      double f21z = -f12z;
+      float f12x = f2 * x12 * 0.5f;
+      float f12y = f2 * y12 * 0.5f;
+      float f12z = f2 * z12 * 0.5f;
+      float f21x = -f12x;
+      float f21y = -f12y;
+      float f21z = -f12z;
 
       // accumulate force
       s_fx += f12x - f21x;
       s_fy += f12y - f21y;
       s_fz += f12z - f21z;
 
       // accumulate potential energy and virial
-      s_pe += p2 * 0.5; // two-body potential
+      s_pe += p2 * 0.5f; // two-body potential
       s_sxx += x12 * f21x;
       s_sxy += x12 * f21y;
       s_sxz += x12 * f21z;

src/force/lj.cuh

@@ -22,9 +22,9 @@
 #define MAX_TYPE 10 // == max number of potentials
 
 struct LJ_Para {
-  double s6e4[MAX_TYPE][MAX_TYPE];
-  double s12e4[MAX_TYPE][MAX_TYPE];
-  double cutoff_square[MAX_TYPE][MAX_TYPE];
+  float s6e4[MAX_TYPE][MAX_TYPE];
+  float s12e4[MAX_TYPE][MAX_TYPE];
+  float cutoff_square[MAX_TYPE][MAX_TYPE];
 };
 
 struct LJ_Data {