Merge pull request lammps#3962 from ndtrung81/amoeba-pppm-accuracy

AMOEBA/HIPPO RMS force accuracy estimate

examples/amoeba/amoeba_ubiquitin.key

@@ -12,7 +12,8 @@ ewald
 ewald-alpha 0.4
 pewald-alpha 0.5
 ewald-cutoff                    7.0
-#pme-grid                   60 45 45
 pme-grid                   60 48 48
-pme-order                         5
 polar-eps                   0.00001
+#pme-grid                   15 12 12
+#polar-eps                   0.0002
+pme-order                         5

src/AMOEBA/amoeba_multipole.cpp

@@ -18,6 +18,7 @@
 #include "atom.h"
 #include "comm.h"
 #include "domain.h"
+#include "force.h"
 #include "math_const.h"
 #include "math_special.h"
 #include "neigh_list.h"
@@ -1000,3 +1001,61 @@ void PairAmoeba::damppole(double r, int rorder, double alphai, double alphak,
     }
   }
 }
+
+/* ----------------------------------------------------------------------
+   estimate the accuracy of m_kspace solver based on the monopoles
+   based on Ewald
+------------------------------------------------------------------------- */
+
+double PairAmoeba::final_accuracy_mpole()
+{
+  const int nlocal = atom->nlocal;
+  double qsqsum_local(0.0), qsqsum;
+  for (int i = 0; i < nlocal; i++) {
+    qsqsum_local += rpole[i][0]*rpole[i][0];
+  }
+  MPI_Allreduce(&qsqsum_local,&qsqsum,1,MPI_DOUBLE,MPI_SUM,world);
+  double q2 = qsqsum * force->qqrd2e;
+
+  const double * const prd = domain->prd;
+  const double xprd = prd[0];
+  const double yprd = prd[1];
+  const double zprd = prd[2];
+  const double slab_volfactor = 1.0;
+  const double zprd_slab = zprd*slab_volfactor;
+  bigint natoms = atom->natoms;
+
+  int nx_fft = m_kspace->nx;
+  int ny_fft = m_kspace->ny;
+  int nz_fft = m_kspace->nz;
+  double cutoff = mpolecut;
+
+  double lprx = rms(nx_fft,xprd,natoms,aeewald,q2);
+  double lpry = rms(ny_fft,yprd,natoms,aeewald,q2);
+  double lprz = rms(nz_fft,zprd_slab,natoms,aeewald,q2);
+  double lpr = sqrt(lprx*lprx + lpry*lpry + lprz*lprz) / sqrt(3.0);
+  double q2_over_sqrt = q2 / sqrt(natoms*cutoff*xprd*yprd*zprd_slab);
+  double spr = 2.0 *q2_over_sqrt * exp(-aeewald*aeewald*cutoff*cutoff);
+  double tpr = 0;
+  double estimated_accuracy = sqrt(lpr*lpr + spr*spr + tpr*tpr);
+
+  two_charge_force = force->qqr2e *
+    (force->qelectron * force->qelectron) /
+    (force->angstrom * force->angstrom);
+
+  return estimated_accuracy;
+}
+
+/* ----------------------------------------------------------------------
+   compute RMS accuracy for a dimension
+------------------------------------------------------------------------- */
+
+double PairAmoeba::rms(int km, double prd, bigint natoms, double g_ewald, double q2)
+{
+  if (natoms == 0) natoms = 1;   // avoid division by zero
+  double value = 2.0*q2*g_ewald/prd *
+    sqrt(1.0/(MY_PI*km*natoms)) *
+    exp(-MY_PI*MY_PI*km*km/(g_ewald*g_ewald*prd*prd));
+
+  return value;
+}

src/AMOEBA/pair_amoeba.cpp

@@ -347,8 +347,6 @@ void PairAmoeba::compute(int eflag, int vflag)
     }
   }
 
-  first_flag_compute = 0;
-
   // -------------------------------------------------------------------
   // end of one-time initializations
   // -------------------------------------------------------------------
@@ -428,6 +426,12 @@ void PairAmoeba::compute(int eflag, int vflag)
   else cfstyle = SETUP_HIPPO;
   comm->forward_comm(this);
 
+  // output FF settings to screen and logfile
+  //   delay until here because RMS force accuracy is computed based on rpole
+
+  if (first_flag_compute && (comm->me == 0)) print_settings();
+  first_flag_compute = 0;
+
   if (amoeba) pbc_xred();
   time1 = platform::walltime();
 
@@ -978,10 +982,6 @@ void PairAmoeba::init_style()
     for (int i = 0; i < nlocal; i++) pval[i] = 0.0;
   }
 
-  // output FF settings to screen and logfile
-
-  if (first_flag && (comm->me == 0)) print_settings();
-
   // all done with one-time initializations
 
   first_flag = 0;
@@ -1098,9 +1098,13 @@ void PairAmoeba::print_settings()
 
   if (use_ewald) {
     choose(MPOLE_LONG);
-    mesg += fmt::format("  multipole: cut {} aewald {} bsorder {} FFT {} {} {} "
+    double estimated_accuracy = final_accuracy_mpole();
+    mesg += fmt::format("  multipole: cut {} aewald {} bsorder {} FFT {} {} {}; "
+                        "estimated absolute RMS force accuracy = {:.8g}; "
+                        "estimated relative RMS force accuracy = {:.8g}; "
                         "mscale {} {} {} {}\n",
                         sqrt(off2),aewald,bseorder,nefft1,nefft2,nefft3,
+                        estimated_accuracy,estimated_accuracy/two_charge_force,
                         special_mpole[1],special_mpole[2],special_mpole[3],special_mpole[4]);
   } else {
     choose(MPOLE);

src/AMOEBA/pair_amoeba.h

@@ -419,6 +419,10 @@ class PairAmoeba : public Pair {
   double ewaldcof(double);
   int factorable(int);
 
+  double final_accuracy_mpole();
+  double rms(int km, double prd, bigint natoms, double g_ewald, double q2);
+  double two_charge_force;
+
   // debug methods
 
   FILE *fp_uind;