Commits

Rio Yokota committed b6309c9

Update wrapper.

Comments (0)

Files changed (3)

 test_parallel_coulombVdW: test_parallel_coulombVdW.cxx
 	rm -f libparallel_coulombVdW.a
 	make libparallel_coulombVdW.a
-	$(CXX) $? -L. -lparallel_coulombVdW -lmd3withg80vg $(LFLAGS)
+	$(CXX) $? -L. -lparallel_coulombVdW $(LFLAGS)
 	mpirun -np 2 ./a.out
 
 ewald: ewald.cxx $(MR3)

wrapper/serial_coulombVdW.cxx

 #include "serialfmm.h"
-typedef real              realTyp; // Real Type per Kerenl setup. real variable overwritten in md.h  
 extern "C" {
 #include "md.h"
 #include "vtgrapeproto.h"
 }
 
-#define realTyp float
-
-extern "C" void FMMcalccoulomb_ij(int ni, realTyp* xi, realTyp* qi, realTyp* fi,
-  int nj, realTyp* xj, realTyp* qj, realTyp, int tblno, realTyp size, int periodicflag,
-  int ksize, realTyp alpha) {
-  std::cout << "NEW tblno: " << tblno << std::endl;
-  IMAGES = ((periodicflag & 0x1) == 0) ? 0 : 5;
+extern "C" void FMMcalccoulomb_ij(int ni, double* xi, double* qi, double* fi,
+  int nj, double* xj, double* qj, double, int tblno, double size, int periodicflag) {
+  std::cout << "tblno: " << tblno << std::endl;
+  IMAGES = ((periodicflag & 0x1) == 0) ? 0 : 3;
   THETA = .5;
   Bodies bodies(ni),jbodies(nj);
   Cells cells,jcells;
   FMM.finalize();
 
 #if 1
-  // dipole correction.
-  realTyp dipole[3]={0.0,0.0,0.0};
-  for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) {
-    for( int d=0; d!=3; ++d ) {
-      dipole[d] += (B->X[d] - size/2 ) * B->SRC;
-    }
-  }
-  realTyp coef = 4*M_PI / (3 * size  * size  * size );
-  realTyp coefP = coef*(dipole[0] * dipole[0] + dipole[1] * dipole[1] + dipole[2] * dipole[2])/bodies.size();
-  realTyp coefF[3] = {coef*dipole[0], coef*dipole[1], coef*dipole[2]};
-
-  // sending the data back.
   for( B_iter B=bodies.begin(); B!=bodies.end(); ++B ) {
     int i = B-bodies.begin();
+//    xi[3*i+0] = B->X[0];
+//    xi[3*i+1] = B->X[1];
+//    xi[3*i+2] = B->X[2];
+//    qi[i]     = B->SRC;
     switch (tblno) {
     case 0 :
-      fi[3*i+0] = -B->SRC * (B->TRG[1] - coefF[0]);
-      fi[3*i+1] = -B->SRC * (B->TRG[2] - coefF[1]);
-      fi[3*i+2] = -B->SRC * (B->TRG[3] - coefF[2]);
+      fi[3*i+0] = -B->SRC * B->TRG[1];
+      fi[3*i+1] = -B->SRC * B->TRG[2];
+      fi[3*i+2] = -B->SRC * B->TRG[3];
       break;
     case 1 :
-      fi[3*i+0] = 0.5 * (B->SRC * B->TRG[0] - coefP);
+      fi[3*i+0] = 0.5 * B->SRC * B->TRG[0];
       break;
     }
   }
   switch (tblno) {
   case 0 :
     for( int i=0; i!=ni; ++i ) {
-      realTyp Fx = 0, Fy = 0, Fz = 0;
+      double Fx = 0, Fy = 0, Fz = 0;
       for( int j=0; j!=nj; ++j ) {
-        realTyp dx = xi[3*i+0] - xj[3*j+0];
-        realTyp dy = xi[3*i+1] - xj[3*j+1];
-        realTyp dz = xi[3*i+2] - xj[3*j+2];
-        realTyp R2 = dx * dx + dy * dy + dz * dz;
-        realTyp invR = 1 / std::sqrt(R2);
+        double dx = xi[3*i+0] - xj[3*j+0];
+        double dy = xi[3*i+1] - xj[3*j+1];
+        double dz = xi[3*i+2] - xj[3*j+2];
+        double R2 = dx * dx + dy * dy + dz * dz;
+        double invR = 1 / std::sqrt(R2);
         if( R2 == 0 ) invR = 0;
-        realTyp invR3 = qj[j] * invR * invR * invR;
+        double invR3 = qj[j] * invR * invR * invR;
         Fx += dx * invR3;
         Fy += dy * invR3;
         Fz += dz * invR3;
     break;
   case 1:
     for( int i=0; i!=ni; ++i ) {
-      realTyp Po = 0;
+      double Po = 0;
       for( int j=0; j!=nj; ++j ) {
-        realTyp dx = xi[3*i+0] - xj[3*j+0];
-        realTyp dy = xi[3*i+1] - xj[3*j+1];
-        realTyp dz = xi[3*i+2] - xj[3*j+2];
-        realTyp R2 = dx * dx + dy * dy + dz * dz;
-        realTyp invR = 1 / std::sqrt(R2);
+        double dx = xi[3*i+0] - xj[3*j+0];
+        double dy = xi[3*i+1] - xj[3*j+1];
+        double dz = xi[3*i+2] - xj[3*j+2];
+        double R2 = dx * dx + dy * dy + dz * dz;
+        double invR = 1 / std::sqrt(R2);
         if( R2 == 0 ) invR = 0;
         Po += qj[j] * invR;
       }
     }
     break;
   }
-  // This is the correction factor from FMM to MD Ewald.
-  realTyp fc[3];
-  for( int d=0; d!=3; ++d ) fc[d]=0;
-  for( int i=0; i!=ni; ++i ) { 
-    for( int d=0; d!=3; ++d ) { 
-      fc[d] += qi[i] * xi[3*i+d];
-    }   
-  }
-  if( tblno == 0 ) { 
-    for( int i=0; i!=ni; ++i ) { 
-      for( int d=0; d!=3; ++d ) {
-        fi[3*i+d] -= 4.0 * M_PI * qi[i] * fc[d] / (3.0 * size * size * size);
-      }
-    }   
-  } else {
-    for( int i=0; i!=ni; ++i ) { 
-      fi[3*i+0] += M_PI / (3.0 * size * size * size)
-                * (fc[0] * fc[0] + fc[1] * fc[1] + fc[2] * fc[2]) / ni;
-    }   
-  }
 #endif
 }
 
-/*
-extern "C" void FMMcalcvdw_ij(int ni, realTyp* xi, int* atypei, realTyp* fi,
-  int nj, realTyp* xj, int* atypej, int nat, realTyp* gscale, realTyp* rscale,
-  int tblno, realTyp size, int periodicflag) {
+extern "C" void FMMcalcvdw_ij(int ni, double* xi, int* atypei, double* fi,
+  int nj, double* xj, int* atypej, int nat, double* gscale, double* rscale,
+  int tblno, double size, int periodicflag) {
   std::cout << "tblno: " << tblno << std::endl;
   IMAGES = ((periodicflag & 0x1) == 0) ? 0 : 3;
   THETA = .5;
   switch (tblno) {
   case 2 :
     for( int i=0; i!=ni; ++i ) {
-      realTyp Fx = 0, Fy = 0, Fz = 0;
+      double Fx = 0, Fy = 0, Fz = 0;
       for( int j=0; j!=nj; ++j ) {
-        realTyp dx = xi[3*i+0] - xj[3*j+0];
-        realTyp dy = xi[3*i+1] - xj[3*j+1];
-        realTyp dz = xi[3*i+2] - xj[3*j+2];
-        realTyp R2 = dx * dx + dy * dy + dz * dz;
+        double dx = xi[3*i+0] - xj[3*j+0];
+        double dy = xi[3*i+1] - xj[3*j+1];
+        double dz = xi[3*i+2] - xj[3*j+2];
+        double R2 = dx * dx + dy * dy + dz * dz;
         if( R2 != 0 ) {
-          realTyp rs = rscale[atypei[i]*nat+atypej[j]];
-          realTyp gs = gscale[atypei[i]*nat+atypej[j]];
-          realTyp R2s = R2 * rs;
+          double rs = rscale[atypei[i]*nat+atypej[j]];
+          double gs = gscale[atypei[i]*nat+atypej[j]];
+          double R2s = R2 * rs;
           if( R2MIN <= R2s && R2s < R2MAX ) {
-            realTyp invR2 = 1.0 / R2s;
-            realTyp invR6 = invR2 * invR2 * invR2;
-            realTyp dtmp = gs * invR6 * invR2 * (2.0 * invR6 - 1.0);
+            double invR2 = 1.0 / R2s;
+            double invR6 = invR2 * invR2 * invR2;
+            double dtmp = gs * invR6 * invR2 * (2.0 * invR6 - 1.0);
             Fx += dx * dtmp;
             Fy += dy * dtmp;
             Fz += dz * dtmp;
     break;
   case 3:
     for( int i=0; i!=ni; ++i ) {
-      realTyp Po = 0;
+      double Po = 0;
       for( int j=0; j!=nj; ++j ) {
-        realTyp dx = xi[3*i+0] - xj[3*j+0];
-        realTyp dy = xi[3*i+1] - xj[3*j+1];
-        realTyp dz = xi[3*i+2] - xj[3*j+2];
-        realTyp R2 = dx * dx + dy * dy + dz * dz;
+        double dx = xi[3*i+0] - xj[3*j+0];
+        double dy = xi[3*i+1] - xj[3*j+1];
+        double dz = xi[3*i+2] - xj[3*j+2];
+        double R2 = dx * dx + dy * dy + dz * dz;
         if( R2 != 0 ) {
-          realTyp rs = rscale[atypei[i]*nat+atypej[j]];
-          realTyp gs = gscale[atypei[i]*nat+atypej[j]];
-          realTyp R2s = R2 * rs;
+          double rs = rscale[atypei[i]*nat+atypej[j]];
+          double gs = gscale[atypei[i]*nat+atypej[j]];
+          double R2s = R2 * rs;
           if( R2MIN <= R2s && R2s < R2MAX ) {
-            realTyp invR2 = 1.0 / R2s;
-            realTyp invR6 = invR2 * invR2 * invR2;
+            double invR2 = 1.0 / R2s;
+            double invR6 = invR2 * invR2 * invR2;
             Po += gs * invR6 * (invR6 - 1.0);
           }
         }
 #endif
 }
 
-extern "C" void fmmcalccoulomb_ij_(int *ni, realTyp* xi, realTyp* qi, realTyp* fi,
-  int *nj, realTyp* xj, realTyp* qj, realTyp *rscale, int *tblno, realTyp *size, int *periodicflag) {
+extern "C" void fmmcalccoulomb_ij_(int *ni, double* xi, double* qi, double* fi,
+  int *nj, double* xj, double* qj, double *rscale, int *tblno, double *size, int *periodicflag) {
   std::cout << "Starting FMM" << std::endl;
   FMMcalccoulomb_ij(*ni,xi,qi,fi,*nj,xj,qj,*rscale,*tblno-6,*size,*periodicflag);
 }
 
-extern "C" void fmmcalccoulomb_ij_exlist_(int *ni, realTyp* xi, realTyp* qi, realTyp* fi,
-  int *nj, realTyp* xj, realTyp* qj, realTyp *rscale, int *tblno, realTyp *size, int *periodicflag,
+extern "C" void fmmcalccoulomb_ij_exlist_(int *ni, double* xi, double* qi, double* fi,
+  int *nj, double* xj, double* qj, double *rscale, int *tblno, double *size, int *periodicflag,
   int *numex, int* natex) {
   std::cout << "Starting FMM" << std::endl;
   FMMcalccoulomb_ij(*ni,xi,qi,fi,*nj,xj,qj,*rscale,*tblno-6,*size,*periodicflag);
   switch (*tblno-6) {
   case 0 :
-#if 0
-    for( int i=0,ic=0; i<*ni; i++ ) {
-      for( int j=0; j<numex[i]; j++,ic++ ) natex[ic]--;
-    }
-//    MR3calccoulomb_nlist_ij_host(*ni,xi,qi,fi,*nj,xj,qj,
-//                                *rscale,*tblno-6,*size,*periodicflag&3,numex,natex,-1.0);
-//    MR3calccoulomb_nlist_ij_emu(*ni,xi,qi,fi,*nj,xj,qj,
-//			        *rscale,*tblno-6,*size,*periodicflag&3,numex,natex,-1.0);
-    for( int i=0,ic=0; i<*ni; i++ ) {
-      for( int j=0; j<numex[i]; j++,ic++ ) natex[ic]++;
-    }
-#else
     for( int i=0,ic=0; i!=*ni; ++i ) {
-      realTyp Fx = 0, Fy = 0, Fz = 0;
+      double Fx = 0, Fy = 0, Fz = 0;
       for( int in=0; in!=numex[i]; in++,ic++ ) {
         int j = natex[ic]-1;
-        realTyp dx = xi[3*i+0] - xj[3*j+0];
-        realTyp dy = xi[3*i+1] - xj[3*j+1];
-        realTyp dz = xi[3*i+2] - xj[3*j+2];
-        realTyp R2 = dx * dx + dy * dy + dz * dz;
-        realTyp invR = 1 / std::sqrt(R2);
+        double dx = xi[3*i+0] - xj[3*j+0];
+        double dy = xi[3*i+1] - xj[3*j+1];
+        double dz = xi[3*i+2] - xj[3*j+2];
+        double R2 = dx * dx + dy * dy + dz * dz;
+        double invR = 1 / std::sqrt(R2);
         if( R2 == 0 ) invR = 0;
-        realTyp invR3 = qj[j] * invR * invR * invR;
+        double invR3 = qj[j] * invR * invR * invR;
         Fx += dx * invR3;
         Fy += dy * invR3;
         Fz += dz * invR3;
       fi[3*i+1] -= qi[i] * Fy;
       fi[3*i+2] -= qi[i] * Fz;
     }
-#endif
     break;
   case 1:
-    realTyp Total = 0;
+    double Total = 0;
     int ic = 0;
     for( int i=0; i!=*ni; ++i ) {
-      realTyp Po = 0;
+      double Po = 0;
       for( int in=0; in!=numex[i]; in++,ic++ ) {
         int j = natex[ic]-1;
-        realTyp dx = xi[3*i+0] - xj[3*j+0];
-        realTyp dy = xi[3*i+1] - xj[3*j+1];
-        realTyp dz = xi[3*i+2] - xj[3*j+2];
-        realTyp R2 = dx * dx + dy * dy + dz * dz;
-        realTyp invR = 1 / std::sqrt(R2);
+        double dx = xi[3*i+0] - xj[3*j+0];
+        double dy = xi[3*i+1] - xj[3*j+1];
+        double dz = xi[3*i+2] - xj[3*j+2];
+        double R2 = dx * dx + dy * dy + dz * dz;
+        double invR = 1 / std::sqrt(R2);
         assert( i != j );
         if( R2 == 0 ) invR = 0;
         Po += qj[j] * invR;
 }
 
 
-extern "C" void fmmcalcvdw_ij_(int *ni, realTyp* xi, int* atypei, realTyp* fi,
-  int *nj, realTyp* xj, int* atypej, int *nat, realTyp* gscale, realTyp* rscale,
-  int *tblno, realTyp *size, int *periodicflag) {
+extern "C" void fmmcalcvdw_ij_(int *ni, double* xi, int* atypei, double* fi,
+  int *nj, double* xj, int* atypej, int *nat, double* gscale, double* rscale,
+  int *tblno, double *size, int *periodicflag) {
   std::cout << "Starting FMM" << std::endl;
   for( int i=0; i<*ni; i++ ) atypei[i]--;
   for( int i=0; i<*nj; i++ ) atypej[i]--;
   for( int i=0; i<*nj; i++ ) atypej[i]++;
 }
 
-extern "C" void fmmcalcvdw_ij_exlist_(int *ni, realTyp* xi, int* atypei, realTyp* fi,
-  int *nj, realTyp* xj, int* atypej, int *nat, realTyp* gscale, realTyp* rscale,
-  int *tblno, realTyp *size, int *periodicflag, int* numex, int* natex) {
+extern "C" void fmmcalcvdw_ij_exlist_(int *ni, double* xi, int* atypei, double* fi,
+  int *nj, double* xj, int* atypej, int *nat, double* gscale, double* rscale,
+  int *tblno, double *size, int *periodicflag, int* numex, int* natex) {
   std::cout << "Starting FMM" << std::endl;
   for( int i=0; i<*ni; i++ ) atypei[i]--;
   for( int i=0; i<*nj; i++ ) atypej[i]--;
   for( int i=0,ic=0; i<*ni; i++ ) {
     for( int j=0; j<numex[i]; j++,ic++ ) natex[ic]--;
   }
-#if 0
-  MR3calcvdw_ij(*ni,xi,atypei,fi,*nj,xj,atypej,*nat,gscale,rscale,*tblno,*size,*periodicflag);
-  MR3calcvdw_nlist_ij_emu(*ni,xi,atypei,fi,*nj,xj,atypej,*nat,gscale,rscale,*tblno,
-                          *size,(*periodicflag & 3),numex,natex,-1.0);
-#else
   FMMcalcvdw_ij(*ni,xi,atypei,fi,*nj,xj,atypej,*nat,gscale,rscale,*tblno,*size,*periodicflag);
-//  MR3calcvdw_nlist_ij_host(*ni,xi,atypei,fi,*nj,xj,atypej,*nat,gscale,rscale,*tblno,
-//                          *size,(*periodicflag & 3),numex,natex,-1.0);
-#endif
+  //MR3calcvdw_nlist_ij_host(*ni,xi,atypei,fi,*nj,xj,atypej,*nat,gscale,rscale,*tblno,
+  //                        *size,(*periodicflag & 3),numex,natex,-1.0);
   for( int i=0; i<*ni; i++ ) atypei[i]++;
   for( int i=0; i<*nj; i++ ) atypej[i]++;
   for( int i=0,ic=0; i<*ni; i++ ) {
     for( int j=0; j<numex[i]; j++,ic++ ) natex[ic]++;
   }
 }
-*/
+

wrapper/test_serial_coulombVdW.cxx

 const float R2MIN    = 0.0001;                                   //!< Minimum value for L-J R^2
 const float R2MAX    = 100.0;                                    //!< Maximum value for L-J R^2
 extern "C" {
-//#include "md.h"
-//#include "vtgrapeproto.h"
+#include "md.h"
+#include "vtgrapeproto.h"
 }
 //#define DATAFILE
 
 #else
   const int N = 10000;
   const int nat = 16;
-  const double size = 2;
+  const double size = 20;
 #endif
   double *xi     = new double [3*N];
   double *qi     = new double [N];
   }
   fid.close();
   for( int i=0; i<N; i++ ) {
-//    xi[3*i+0] = drand48() * size - size/2;
-//    xi[3*i+1] = drand48() * size - size/2;
-//    xi[3*i+2] = drand48() * size - size/2;
     xj[3*i+0] = xi[3*i+0];
     xj[3*i+1] = xi[3*i+1];
     xj[3*i+2] = xi[3*i+2];
     fd[3*i+2] += qi[i] * Fz;
   }
 #endif
-  double fc[3];
-  for( int d=0; d<3; ++d ) fc[d]=0;
-  for( int i=0; i<N; ++i ) {
-    for( int d=0; d<3; ++d ) {
-      fc[d] += qi[i] * xi[3*i+d];
-    }
-  }
-  for( int i=0; i<N; ++i ) {
-    pd[3*i+0] += 2.0 * M_PI / (3.0 * size * size * size)
-              * (fc[0] * fc[0] + fc[1] * fc[1] + fc[2] * fc[2]) / N;
-  }
-  for( int i=0; i<N; ++i ) {
-    for( int d=0; d<3; ++d ) {
-      fd[3*i+d] -= 4.0 * M_PI * qi[i] * fc[d] / (3.0 * size * size * size);
-    }
-  }
   double Pd = 0, Pn = 0, Fd = 0, Fn = 0;
   for( int i=0; i<N; i++ ) {
     pd[3*i+0] *= 0.5;