Commits

Rio Yokota committed 6d3b7a6

Added exclusion.

  • Participants
  • Parent commits 4668267

Comments (0)

Files changed (2)

File wrappers/charmm.cxx

   int ic = 0;
   for (int i=0; i<*nglobal; i++) {
     if (icpumap[i] == 1) {
-      double pp = 0, fx = 0, fy = 0, fz = 0;
       for (int jc=0; jc<numex[i]; jc++, ic++) {
 	int j = natex[ic]-1;
 	float dx = x[3*i+0] - x[3*j+0];
 	float dy = x[3*i+1] - x[3*j+1];
 	float dz = x[3*i+2] - x[3*j+2];
+	if (dx < -*cycle/2) dx += *cycle;
+	if (dy < -*cycle/2) dy += *cycle;
+	if (dz < -*cycle/2) dz += *cycle;
+	if (dx >  *cycle/2) dx -= *cycle;
+	if (dy >  *cycle/2) dy -= *cycle;
+	if (dz >  *cycle/2) dz -= *cycle;
 	float R2 = dx * dx + dy * dy + dz * dz;
 	float invR = 1 / std::sqrt(R2);
 	if (R2 == 0) invR = 0;
-	float invR3 = q[j] * invR * invR * invR;
-	pp -= q[j] * invR;
-	fx += dx * invR3;
-	fy += dy * invR3;
-	fz += dz * invR3;
+	float invR3 = invR * invR * invR;
+        float qinvR3 = q[j] * invR3;
+	p[i] -= q[j] * invR;
+	f[3*i+0] += dx * qinvR3;
+	f[3*i+1] += dy * qinvR3;
+	f[3*i+2] += dz * qinvR3;
+	p[j] -= q[i] * invR;
+	qinvR3 = q[i] * invR3;
+	f[3*j+0] -= dx * qinvR3;
+	f[3*j+1] -= dy * qinvR3;
+	f[3*j+2] -= dz * qinvR3;
       }
-      p[i] -= pp;
-      f[3*i+0] -= fx;
-      f[3*i+1] -= fy;
-      f[3*i+2] -= fz;
     } else {
       ic += numex[i];
     }

File wrappers/test_charmm.f90

   call mpi_comm_size(mpi_comm_world, mpisize, ierr)
   call mpi_comm_rank(mpi_comm_world, mpirank, ierr)
   nglobal = 1000
-  images = 1
+  images = 3
   ksize = 11
   pcycle = 2 * pi
   sigma = .25 / pi
   end do
   call fmm_init(images)
   call fmm_partition(nglobal, icpumap, x, q, pcycle)
-  call fmm(nglobal, icpumap, x, q, p, f, pcycle, natex, numex)
+  call fmm(nglobal, icpumap, x, q, p, f, pcycle)
+  call exclusion(nglobal, icpumap, x, q, p, f, pcycle, numex, natex)
   do i = 1, nglobal
      x2(3*i-2) = x(3*i-2)
      x2(3*i-1) = x(3*i-1)
      f2(3*i-1) = 0
      f2(3*i-0) = 0
   end do
-#if 0
+#if 1
   call ewald(nglobal, icpumap, x2, q2, p2, f2, ksize, alpha, sigma, cutoff, pcycle)
 #else
   prange = 0
      f2(3*i-0) = f2(3*i-0) - coef * dipole(3)
   end do
 #endif
+  call exclusion(nglobal, icpumap, x2, q2, p2, f2, pcycle, numex, natex)
   potSum = 0
   potSum2 = 0
   accDif = 0