v1.68 fixed the bug on the first and update charge and B0_prev in the CELL::Q model

Author: kylincaster

README.md

@@ -42,7 +42,8 @@ Note tat a NIL `compute pe/atom_GA` with the compute name of `__Electro_Control_
 
 **Otherwise**, copy all files in the `src_mod` folder to `src` and edit the`src\MAKE\OPTION\Makefile.XXX `. Then, type `make XXX` in the `src` folder. Here, **Makefile.intel_cpu_intelmpi** is also provided in `misc` folder as an example for installation configuration.
 
-Note that the default **single** precision of the FFTW setup is not higher enough for constant potential simulation that results in the order of error of 1e-6. To improve the accuracy, we suggests to switch to  **double** precision using the macro of `-DFFT_DOUBLE`. Thus, if the error of `undefined reference to 'LAMMPS_NS::FFT3d::compute(double*, double*, int)'` is encountered, type `make clean-XXX` to remove the existing **single** precision FFTW library and recompile it.
+Note that the default **single** precision of the FFTW setup is not higher enough for constant potential simulation that results in the order of error of 1e-6. To improve the accuracy, we suggests to switch to  **double** precision using the macro of `-DFFT_DOUBLE`. Thus, if the error of `undefined reference to 'LAMMPS_NS::FFT3d::compute(double*, double*, int)'` is encountered, type `make no-KSPACE;
+make yes-KSPACE` to recompile the existing *kspace* code compiled with **single** precision FFTW library.
 
 ## Syntax: pair_style lj/cut/point/long
 

makefile

@@ -15,17 +15,29 @@
 #include <cstdio>
 #include <unistd.h>
 
-namespace LAMMPS_NS {
+namespace LAMMPS_NS
+{
 
-void print_vec(double *arr, int N, const char *name) {
+void print_vec(double* arr, int N, const char* name)
+{
   printf("  %s = ", name);
   for (int i = 0; i < N; ++i) {
     printf("%12.10g ", arr[i]);
   }
   printf("\n");
 }
 
-double norm_mat(double *mat, int N) {
+void print_vec(int* arr, int N, const char* name)
+{
+  printf("  %s = ", name);
+  for (int i = 0; i < N; ++i) {
+    printf("%d ", arr[i]);
+  }
+  printf("\n");
+}
+
+double norm_mat(double* mat, int N)
+{
   double res = 0;
   for (int i = 0; i < N; ++i) {
     res += mat[i] * mat[i];
@@ -36,55 +48,61 @@ double norm_mat(double *mat, int N) {
     return sqrt(res);
 }
 
-double norm_mat_p(double *mat, int N, const char *info) {
+double norm_mat_p(double* mat, int N, const char* info)
+{
   double value = norm_mat(mat, N);
   printf("%3s %.16g with size = %d\n", info, value, N);
   return value;
 }
 
-void write_mat(const double *const ptr, int N, const char *name) {
-  FILE *outa = fopen(name, "wb");
+void write_mat(const double* const ptr, int N, const char* name)
+{
+  FILE* outa = fopen(name, "wb");
   printf("write as %s with size = %d\n", name, N);
-  if (ptr == NULL)
-    printf("the input ptr is NULL.\n");
+  if (ptr == NULL) printf("the input ptr is NULL.\n");
   fwrite(ptr, sizeof(double), N, outa);
   fclose(outa);
 }
 
-double hash_vec(double *mat, int N) {
+double hash_vec(double* mat, int N)
+{
   double res = 0, inv_N = 1 / double(N);
   for (int i = 0; i < N; ++i) {
     res += (mat[i] + i * inv_N) * (mat[i] + i * inv_N);
   }
   return sqrt(res);
 }
 
-double hash_vec(int *mat, int N) {
+double hash_vec(int* mat, int N)
+{
   double res = 0, inv_N = 1 / double(N);
   for (int i = 0; i < N; ++i) {
     res += (mat[i] + i * inv_N) * (mat[i] + i * inv_N);
   }
   return sqrt(res);
 }
 
-double *load_mat(int N, const char *name) {
-  double *ptr = (double *)malloc(sizeof(double) * N);
-  FILE *infile = fopen(name, "r");
-  int result = fread(ptr, sizeof(double), N, infile);
+double* load_mat(int N, const char* name)
+{
+  double* ptr  = (double*)malloc(sizeof(double) * N);
+  FILE* infile = fopen(name, "r");
+  int result   = fread(ptr, sizeof(double), N, infile);
   if (result != N) {
     printf("cannot read %d double from %s\n", N, name);
   }
   return ptr;
 }
 
-void print_hash(double *arr, int N, const char *name) {
+void print_hash(double* arr, int N, const char* name)
+{
   printf("  %s = ", name);
-  double hash = hash_vec(arr, N);
+  double hash   = hash_vec(arr, N);
   double normal = norm_mat(arr, N);
   printf("hash = %12.10g, Norm. = %12.10g Size = %d\n", hash, normal, N);
 }
 
-void eext(int n) {
+void eext(int n)
+{
   printf("LOG: exit with %d\n", n);
   fflush(stdout);
   sleep(n);

src_mod/debug_func.cpp

@@ -11,16 +11,18 @@
 #define CONP_GA_DEBUG_FUNC_H
 
 #ifndef CONP_NO_DEBUG
-namespace LAMMPS_NS {
+namespace LAMMPS_NS
+{
 
-void print_vec(double *arr, int N, const char *name);
-double norm_mat(double *mat, int N);
-double norm_mat_p(double *mat, int N, const char *info);
-void write_mat(const double *const ptr, int N, const char *name);
-double hash_vec(double *mat, int N);
-double hash_vec(int *mat, int N);
-double *load_mat(int N, const char *name);
-void print_hash(double *arr, int N, const char *name);
+void print_vec(double* arr, int N, const char* name);
+void print_vec(int* arr, int N, const char* name);
+double norm_mat(double* mat, int N);
+double norm_mat_p(double* mat, int N, const char* info);
+void write_mat(const double* const ptr, int N, const char* name);
+double hash_vec(double* mat, int N);
+double hash_vec(int* mat, int N);
+double* load_mat(int N, const char* name);
+void print_hash(double* arr, int N, const char* name);
 void eext(int n);
 
 } // namespace LAMMPS_NS

src_mod/debug_func.h

@@ -164,7 +164,7 @@ void Electro_Control_GA::setup()
   if (kspace) {
     kspace_conp_Flag = true;
     if (me == 0) {
-      utils::logmesg(lmp, "[ConpGA] using optimized electrostatic solver: %s\n", kspace->type_name());
+      utils::logmesg(lmp, "[ConpGA] using optimized electrostatic solver: {}\n", kspace->type_name());
     }
     kspace->electro         = this;
     kspace->DEBUG_LOG_LEVEL = DEBUG_LOG_LEVEL;
@@ -181,6 +181,8 @@ void Electro_Control_GA::setup()
 
 int* Electro_Control_GA::setup_comm(int* all_ghostGA_ind)
 {
+  if (me == 0) utils::logmesg(lmp, "[ConpGA] Setup inner communicator\n");
+
   int nprocs = comm->nprocs;
   std::map<int, int> tag_to_proc;
   const int* const localGA_tag = &FIX->localGA_tag.front();
@@ -490,22 +492,23 @@ void Electro_Control_GA::calc_GA_potential(STAT cl_stat_new, STAT ek_stat_new, S
   // OP: {0, NIL}, {1: PE->Potential}, {2: Potential->PE}
   // printf("ek_stat = %d\n", ek_stat);
   if (me == 0 && DEBUG_LOG_LEVEL > 1)
-    utils::logmesg(lmp, "[Debug] U<=>POT conv. {cl: {}=>{}} {ek: {}=>{}} sum_op {}\n", cl_stat, cl_stat_new, ek_stat, ek_stat_new, sum_op);
+    utils::logmesg(lmp, "[Debug] U<=>POT conv. {{cl: {}=>{}}} {{ek: {}=>{}}} sum_op {}\n", cl_stat, cl_stat_new, ek_stat, ek_stat_new, sum_op);
 
   if (sum_op != SUM_OP::BOTH) {
     int cl_op        = (3 + static_cast<int>(cl_stat_new) - static_cast<int>(cl_stat)) % 3;
     int ek_op        = (3 + static_cast<int>(ek_stat_new) - static_cast<int>(ek_stat)) % 3;
     bool sum_op_flag = sum_op == SUM_OP::ON;
     if (sum_op_flag && cl_stat_new != ek_stat_new) {
       if (me == 0)
-        utils::logmesg(lmp, "Sum between cl{{}: {}} vs ek{{}: {}}\n", cl_stat, name_stat[static_cast<int>(cl_stat)], ek_stat,
+        utils::logmesg(lmp, "Sum between cl{{{}: {}}} vs ek{{{}: {}}}\n", cl_stat, name_stat[static_cast<int>(cl_stat)], ek_stat,
                        name_stat[static_cast<int>(ek_stat)]);
       error->all(FLERR, "inconsistent sum type!");
     }
     conv_GA_potential(cl_op, ek_op, sum_op_flag);
     ek_stat = ek_stat_new;
     cl_stat = cl_stat_new;
   } else {
+    // SUM_OP::BOTH
     conv_GA_potential_both(cl_stat, ek_stat);
     cl_stat = STAT::E;
     ek_stat = STAT::POT;
@@ -846,4 +849,12 @@ Electro_Control_GA::~Electro_Control_GA()
     memory->sfree(sendlist);
     memory->sfree(recvlist);
   }
+  if (pair != NULL) {
+    pair->DEBUG_LOG_LEVEL = 0;
+    pair->electro         = nullptr;
+  }
+  if (kspace_conp_Flag == true && kspace != nullptr) {
+    kspace->DEBUG_LOG_LEVEL = 0;
+    kspace->electro         = nullptr;
+  }
 }

src_mod/electro_control_GA.cpp

@@ -95,6 +95,8 @@ class Electro_Control_GA : protected Pointers {
     } else
       return kspace->type_name();
   }
+  STAT get_ek_stat() { return ek_stat; }
+  STAT get_cl_stat() { return cl_stat; }
 
   /*
   virtual void compute_vector(bigint *, double *) = 0;