clean up ewald, remove dependency on SimulationRegion

source/api_cc/tests/test_dipolecharge.cc

@@ -6,7 +6,7 @@
 #include "DeepTensor.h"
 #include "DataModifier.h"
 #include "SimulationRegion.h"
-#include "Ewald.h"
+#include "ewald.h"
 #include "neighbor_list.h"
 #include "test_utils.h"
 
@@ -174,12 +174,12 @@ TEST_F(TestDipoleCharge, cpu_lmp_nlist)
   // compute the recp part of the ele interaction
   double eener;
   std::vector<double> eforce, evirial;
-  SimulationRegion<double> region;
-  region.reinitBox(&box[0]);
+  Region<double> region;
+  init_region_cpu(region, &box[0]);
   EwaldParameters<double> eparam;
   eparam.beta = 0.2;
   eparam.spacing = 4;
-  EwaldReciprocal(eener, eforce, evirial, coord, charge, region, eparam);
+  ewald_recp(eener, eforce, evirial, coord, charge, region, eparam);
 
   EXPECT_LT(fabs(eener - expected_e[0]), 1e-6);
   EXPECT_EQ(eforce.size(), coord.size());

source/api_cc/tests/test_ewald.cc

@@ -5,7 +5,7 @@
 #include <vector>
 #include "neighbor_list.h"
 #include "test_utils.h"
-#include "Ewald.h"
+#include "ewald.h"
 
 class TestInferEwald : public ::testing::Test
 {  
@@ -48,9 +48,9 @@ TEST_F(TestInferEwald, cpu_numfv)
 	std::vector<double> &	virial,
 	const std::vector<double> & coord,
 	const std::vector<double> & box) {
-      SimulationRegion<double> region;
-      region.reinitBox(&box[0]);
-      EwaldReciprocal(ener, force, virial, coord, charge, region, eparam);
+      Region<double> region;
+      init_region_cpu(region, &box[0]);
+      ewald_recp(ener, force, virial, coord, charge, region, eparam);
     }
   };
   MyModel model(charge);

source/lib/include/ewald.h

@@ -0,0 +1,34 @@
+#pragma once
+
+#include<algorithm>
+#include<cassert>
+#include<omp.h>
+
+#include "utilities.h"
+#include "region.h"
+
+// 8.988e9 / pc.electron_volt / pc.angstrom * (1.602e-19)**2
+const double ElectrostaticConvertion = 14.39964535475696995031;
+
+template <typename VALUETYPE>
+struct EwaldParameters 
+{
+  VALUETYPE rcut = 6.0;
+  VALUETYPE beta = 2;
+  VALUETYPE spacing = 4;
+};
+
+// compute the reciprocal part of the Ewald sum.
+// outputs: energy force virial
+// inputs: coordinates charges region
+template <typename VALUETYPE>
+void 
+ewald_recp(
+    VALUETYPE &				ener, 
+    std::vector<VALUETYPE> &		force,
+    std::vector<VALUETYPE> &		virial,
+    const std::vector<VALUETYPE>&	coord,
+    const std::vector<VALUETYPE>&	charge,
+    const Region<VALUETYPE>&		region, 
+    const EwaldParameters<VALUETYPE>&	param);
+

source/lib/include/Ewald.h → source/lib/src/ewald.cc

@@ -1,23 +1,6 @@
-#pragma once
-
-#include<algorithm>
-#include<cassert>
-#include<omp.h>
-
-#include "utilities.h"
+#include "ewald.h"
 #include "SimulationRegion.h"
 
-// 8.988e9 / pc.electron_volt / pc.angstrom * (1.602e-19)**2
-const double ElectrostaticConvertion = 14.39964535475696995031;
-
-template <typename VALUETYPE>
-struct EwaldParameters 
-{
-  VALUETYPE rcut = 6.0;
-  VALUETYPE beta = 2;
-  VALUETYPE spacing = 4;
-};
-
 template<typename VALUETYPE> 
 VALUETYPE
 dir_err_esti(const VALUETYPE & test_q,
@@ -82,14 +65,9 @@ rec_err_esti(const VALUETYPE & test_q,
 template <typename VALUETYPE>
 void
 cmpt_k(std::vector<int> & KK,
-       const SimulationRegion<double>&		region, 
+       const VALUETYPE *		boxt, 
        const EwaldParameters<VALUETYPE>&	param)
 {
-  const double * boxt_ = region.getBoxTensor();
-  VALUETYPE boxt[9];
-  for (int dd = 0; dd < 9; ++dd){
-    boxt[dd] = static_cast<VALUETYPE>(boxt_[dd]);
-  }  
   KK.resize(3);
   for (int dd = 0; dd < 3; ++dd){
     VALUETYPE ll = sqrt(dot3(boxt+dd*3, boxt+dd*3));
@@ -108,13 +86,14 @@ cmpt_k(std::vector<int> & KK,
 // inputs: coordinates charges region
 template <typename VALUETYPE>
 void 
-EwaldReciprocal(VALUETYPE &				ener, 
-		std::vector<VALUETYPE> &		force,
-		std::vector<VALUETYPE> &		virial,
-		const std::vector<VALUETYPE>&		coord,
-		const std::vector<VALUETYPE>&		charge,
-		const SimulationRegion<double>&		region, 
-		const EwaldParameters<VALUETYPE>&	param)
+ewald_recp(
+    VALUETYPE &				ener, 
+    std::vector<VALUETYPE> &		force,
+    std::vector<VALUETYPE> &		virial,
+    const std::vector<VALUETYPE>&	coord,
+    const std::vector<VALUETYPE>&	charge,
+    const Region<VALUETYPE>&		region, 
+    const EwaldParameters<VALUETYPE>&	param)
 {
   // natoms
   int natoms = charge.size();
@@ -137,7 +116,7 @@ EwaldReciprocal(VALUETYPE &				ener,
   // K grid
   std::vector<int> KK(3);
   int totK = 1;
-  cmpt_k<VALUETYPE>(KK, region, param);
+  cmpt_k<VALUETYPE>(KK, region.boxt, param);
   for (int dd = 0; dd < 3; ++dd){
     totK *= (KK[dd]+1);
   }  
@@ -154,11 +133,12 @@ EwaldReciprocal(VALUETYPE &				ener,
 #pragma omp parallel for num_threads(nthreads)
   for (int ii = 0; ii < natoms; ++ii){
     int thread_id = omp_get_thread_num();
-    double ir[3];
-    double tmpcoord[3] = {coord[ii*3], coord[ii*3+1], coord[ii*3+2]};
-    region.phys2Inter(ir, tmpcoord);
+    VALUETYPE ir[3];
+    VALUETYPE tmpcoord[3] = {coord[ii*3], coord[ii*3+1], coord[ii*3+2]};
+    convert_to_inter_cpu(ir, region, tmpcoord);
+    // region.phys2Inter(ir, tmpcoord);
     for (int mm0 = -KK[0]/2; mm0 <= KK[0]/2; ++mm0){
-      double mr[3];
+      VALUETYPE mr[3];
       mr[0] = ir[0] * mm0;      
       int shift0 = (mm0 + KK[0]/2) * stride[1] * stride[2];
       for (int mm1 = -KK[1]/2; mm1 <= KK[1]/2; ++mm1){
@@ -168,7 +148,7 @@ EwaldReciprocal(VALUETYPE &				ener,
 	  if (mm0 == 0 && mm1 == 0 && mm2 == 0) continue;
 	  int mc = shift0 + shift1 + mm2 + KK[2]/2;
 	  mr[2] = ir[2] * mm2;
-	  double mdotr = 2. * M_PI * (mr[0]+mr[1]+mr[2]);
+	  VALUETYPE mdotr = 2. * M_PI * (mr[0]+mr[1]+mr[2]);
 	  thread_sqr[thread_id][mc] += charge[ii] * cos(mdotr);
 	  thread_sqi[thread_id][mc] += charge[ii] * sin(mdotr);
 	}
@@ -187,11 +167,7 @@ EwaldReciprocal(VALUETYPE &				ener,
   }  
 
   // get rbox
-  VALUETYPE rec_box[9];
-  const double * rec_box_ = region.getRecBoxTensor();
-  for (int ii = 0; ii < 9; ++ii){
-    rec_box[ii] = static_cast<VALUETYPE>(rec_box_[ii]);
-  }
+  const VALUETYPE * rec_box = region.rec_boxt;
 
   std::vector<VALUETYPE> thread_ener(nthreads, 0.);
   std::vector<std::vector<VALUETYPE> > thread_force(nthreads);
@@ -272,7 +248,7 @@ EwaldReciprocal(VALUETYPE &				ener,
     }
   }
 
-  VALUETYPE vol = static_cast<VALUETYPE>(region.getVolume());
+  VALUETYPE vol = volume_cpu(region);
   ener /= 2 * M_PI * vol;
   ener *= ElectrostaticConvertion;
   for (int ii = 0; ii < 3*natoms; ++ii){
@@ -287,3 +263,25 @@ EwaldReciprocal(VALUETYPE &				ener,
   delete[]sqi;
 }
 
+
+template
+void 
+ewald_recp<float>(
+    float &				ener, 
+    std::vector<float> &		force,
+    std::vector<float> &		virial,
+    const std::vector<float>&		coord,
+    const std::vector<float>&		charge,
+    const Region<float>&		region, 
+    const EwaldParameters<float>&	param);
+
+template
+void 
+ewald_recp<double>(
+    double &				ener, 
+    std::vector<double> &		force,
+    std::vector<double> &		virial,
+    const std::vector<double>&		coord,
+    const std::vector<double>&		charge,
+    const Region<double>&		region, 
+    const EwaldParameters<double>&	param);

source/lib/tests/test_ewald.cc

@@ -1,7 +1,7 @@
 #include <gtest/gtest.h>
 #include <cmath>
 #include <algorithm>
-#include "Ewald.h"
+#include "ewald.h"
 
 class TestEwald : public ::testing::Test
 {

source/op/ewald_recp.cc

@@ -1,5 +1,5 @@
 #include "custom_op.h"
-#include "Ewald.h"
+#include "ewald.h"
 
 typedef double boxtensor_t ;
 typedef double compute_t;
@@ -79,29 +79,19 @@ class EwaldRecpOp : public OpKernel {
       int coord_iter = kk * nloc * 3;
       int charge_iter = kk * nloc;
       // set region
-      boxtensor_t boxt [9] = {0};
-      for (int dd = 0; dd < 9; ++dd) {
-	boxt[dd] = box(box_iter + dd);
-      }
-      SimulationRegion<boxtensor_t > region;
-      region.reinitBox (boxt);
+      Region<FPTYPE> region;
+      init_region_cpu(region, &box(box_iter));
 
       // set & normalize coord
-      std::vector<boxtensor_t > d_coord3_ (nloc*3);
+      std::vector<FPTYPE > d_coord3 (nloc*3);
       for (int ii = 0; ii < nloc; ++ii){
-	for (int dd = 0; dd < 3; ++dd){
-	  d_coord3_[ii*3+dd] = coord(coord_iter + ii*3+dd);
-	}
-	double inter[3];
-	region.phys2Inter (inter, &d_coord3_[3*ii]);
+	FPTYPE inter[3];
+	convert_to_inter_cpu(inter, region, &coord(coord_iter + ii*3));
 	for (int dd = 0; dd < 3; ++dd){
 	  if      (inter[dd] < 0 ) inter[dd] += 1.;
 	  else if (inter[dd] >= 1) inter[dd] -= 1.;
 	}
-      }
-      std::vector<FPTYPE > d_coord3 (nloc*3);
-      for (int ii = 0; ii < nloc * 3; ++ii) {
-	d_coord3[ii] = d_coord3_[ii];
+	convert_to_phys_cpu(&d_coord3[ii*3], region, inter);
       }
 
       // set charge
@@ -114,7 +104,7 @@ class EwaldRecpOp : public OpKernel {
       std::vector<FPTYPE> d_virial(9);
 
       // compute
-      EwaldReciprocal(d_ener, d_force, d_virial, d_coord3, d_charge, region, ep);
+      ewald_recp(d_ener, d_force, d_virial, d_coord3, d_charge, region, ep);
 
       // copy output
       energy(kk) = d_ener;