Derive the default PPPM spacing of qNEP models from the radial cutoff

doc/gpumd-replica/input_parameters/kspace.rst

@@ -10,14 +10,30 @@ Syntax
 
 ::
 
-  kspace <method>
+  kspace ewald
+  kspace pppm [spacing]
 
 Select ``pppm`` (the default) or ``ewald`` for charge NEP electrostatics.
 This setting has no effect for ordinary NEP models.
 
+The optional ``spacing`` parameter sets the maximum PPPM mesh spacing in Angstrom
+(allowed range: 0.2 to 2.0, inclusive).
+The default spacing is :math:`\min(1, r_c^R/6)` Angstrom,
+where :math:`r_c^R` is the radial cutoff of the model in Angstrom.
+The mesh grows as needed when the box changes.
+
+At the first force evaluation, ``gpumd_replica`` prints the mesh
+and an estimate of the root-mean-square error of the PPPM forces relative to the Ewald sum.
+The estimate assumes randomly distributed charges.
+For a crystal, the actual error can be smaller by an order of magnitude or more.
+
 Example
 -------
 
 ::
 
   kspace pppm
+
+To request a PPPM spacing of at most 0.8 Angstrom::
+
+  kspace pppm 0.8

doc/gpumd-replica/output_files/restart.rst

@@ -25,7 +25,7 @@ Caveats
 
 * Restart files with an incompatible format are rejected.
 * ``run.in``, ``model.xyz``, and the potential file are still required.
-  Preserve atom count, species, masses, box, potential content, replica count, time step, thermostat, and method parameters.
+  Preserve atom count, species, masses, box, potential content, replica count, time step, thermostat, ``kspace`` setting, and method parameters.
 * REMD resume cannot request positive ``equilibrate``.
 * Restarts are written at normal completion, not periodically.
   Appended diagnostics and trajectories should correspond to the saved checkpoint.

doc/gpumd/input_parameters/kspace.rst

@@ -17,13 +17,25 @@ Syntax
 
 The default method is ``pppm`` (particle-particle particle-mesh).
 Its optional ``spacing`` parameter sets the maximum mesh spacing in Angstrom
-(default: 1.0; allowed range: 0.2 to 2.0, inclusive).
+(allowed range: 0.2 to 2.0, inclusive).
 Smaller values give finer meshes at a higher computational cost.
 
+The default spacing is :math:`\min(1, r_c^R/6)` Angstrom,
+where :math:`r_c^R` is the radial cutoff of the model in Angstrom.
+The Ewald splitting parameter is :math:`\alpha = \pi/r_c^R`.
+The PPPM error grows with the product of :math:`\alpha` and the spacing.
+The default spacing keeps this product at or below :math:`\pi/6`.
+
 The mesh is chosen automatically and grows as needed when the box changes.
 It does not shrink. The actual spacing may be smaller than the requested value.
 
-Specify ``kspace`` at most once in ``run.in``, before the first ``run``.
+At the first force evaluation, GPUMD prints the mesh
+and an estimate of the root-mean-square error of the PPPM forces relative to the Ewald sum [Deserno1998]_.
+The estimate uses the charges of that step and assumes that they are randomly distributed.
+For a crystal, the actual error can be smaller by an order of magnitude or more.
+
+Specify ``kspace`` at most once in ``run.in``.
+The setting applies to every ``run``, wherever the line appears in ``run.in``.
 
 Example
 -------
@@ -35,3 +47,8 @@ To use Ewald::
 To request a PPPM spacing of at most 1.5 Angstrom use::
 
    kspace pppm 1.5
+
+References
+----------
+
+.. [Deserno1998] M. Deserno and C. Holm, *How to mesh up Ewald sums. II. An accurate error estimate for the particle-particle-particle-mesh algorithm*, J. Chem. Phys. **109**, 7694 (1998).

src/force/nep_charge.cu

@@ -49,6 +49,8 @@ const std::string ELEMENTS[NUM_ELEMENTS] = {
 
 void NEP_Charge::check_ewald_pppm(const RunInput& run_input)
 {
+  // With alpha = pi / rc_radial, the default spacing keeps alpha * spacing at or below pi / 6.
+  pppm_spacing = std::fmin(1.0, paramb.rc_radial / 6.0);
   for (const auto& line : run_input.lines()) {
     const std::vector<std::string>& tokens = line.tokens;
     if (!tokens.empty() && tokens[0] == "kspace") {

src/force/nep_charge.cuh

@@ -183,7 +183,7 @@ private:
   bool need_bec = false;
   void check_need_bec(const RunInput& run_input);
   bool use_pppm = true; // use PPPM by default
-  double pppm_spacing = 1.0;
+  double pppm_spacing;
   void check_ewald_pppm(const RunInput& run_input);
   bool has_dftd3 = false;
   void initialize_dftd3(const RunInput& run_input);

src/force/pppm.cu

@@ -24,6 +24,7 @@ The k-space part of the PPPM method.
 #include <cmath>
 #include <cstdio>
 #include <iostream>
+#include <vector>
 
 namespace{
 
@@ -623,7 +624,6 @@ void PPPM::find_para(const int N, const Box& box)
   const float two_pi = 6.2831853f;
   const double volume = box.get_volume();
   para.two_pi_over_V = two_pi / volume;
-  const bool first_mesh = !plan_initialized;
   for (int d = 0; d < 3; ++d) {
     const double required = volume / box.get_area(d) / mesh_spacing;
     if (required > para.K[d]) {
@@ -641,14 +641,6 @@ void PPPM::find_para(const int N, const Box& box)
     para.K0K1K2 = static_cast<int>(number_of_points);
     allocate_memory();
   }
-  if (first_mesh) {
-    printf(
-      "PPPM mesh: %d x %d x %d (target spacing %.17g A; actual spacing %.17g %.17g %.17g A).\n",
-      para.K[0], para.K[1], para.K[2], mesh_spacing,
-      volume / box.get_area(0) / para.K[0],
-      volume / box.get_area(1) / para.K[1],
-      volume / box.get_area(2) / para.K[2]);
-  }
   para.potential_factor = K_C_SP / N;
   for (int d = 0; d < 3; ++d) {
     para.b[0][d] = two_pi * (float)box.cpu_h[9 + d];
@@ -657,6 +649,54 @@ void PPPM::find_para(const int N, const Box& box)
   }
 }
 
+void PPPM::print_mesh(const int N1, const int N2, const Box& box, const GPU_Vector<float>& charge)
+{
+  const double volume = box.get_volume();
+  printf(
+    "PPPM mesh: %d x %d x %d (target spacing %.17g A; actual spacing %.17g %.17g %.17g A).\n",
+    para.K[0],
+    para.K[1],
+    para.K[2],
+    mesh_spacing,
+    volume / box.get_area(0) / para.K[0],
+    volume / box.get_area(1) / para.K[1],
+    volume / box.get_area(2) / para.K[2]);
+
+  // Root-mean-square force error for ik differentiation and five-point charge assignment
+  // with the optimal influence function, from M. Deserno and C. Holm, JCP 109, 7694 (1998).
+  // The coefficients are those of PPPM::estimate_ik_error in LAMMPS.
+  const double coefficients[5] = {
+    1.0 / 23232.0,
+    7601.0 / 13628160.0,
+    143.0 / 69120.0,
+    517231.0 / 106536960.0,
+    106640677.0 / 11737571328.0};
+  const int number_of_charges = N2 - N1;
+  std::vector<float> q(number_of_charges);
+  CHECK(gpuMemcpy(
+    q.data(), charge.data() + N1, sizeof(float) * number_of_charges, gpuMemcpyDeviceToHost));
+  double sum_of_squared_charges = 0.0;
+  for (const float q_i : q) {
+    sum_of_squared_charges += double(q_i) * q_i;
+  }
+  const double alpha = para.alpha;
+  double squared_error = 0.0;
+  for (int d = 0; d < 3; ++d) {
+    const double thickness = volume / box.get_area(d);
+    const double h_alpha = thickness / para.K[d] * alpha;
+    double sum = 0.0;
+    for (int m = 0; m < 5; ++m) {
+      sum += coefficients[m] * std::pow(h_alpha, 2 * m);
+    }
+    const double error =
+      K_C_SP * sum_of_squared_charges * std::pow(h_alpha, 5) *
+      std::sqrt(alpha * thickness * std::sqrt(2.0 * PI) * sum / number_of_charges) /
+      (thickness * thickness);
+    squared_error += error * error;
+  }
+  printf("PPPM force error estimate: %.3e eV/A.\n", std::sqrt(squared_error / 3.0));
+}
+
 void PPPM::find_force(
   const int N,
   const int N1,
@@ -669,7 +709,11 @@ void PPPM::find_force(
   GPU_Vector<double>& virial_per_atom,
   GPU_Vector<double>& potential_per_atom)
 {
+  const bool first_mesh = !plan_initialized;
   find_para(N, box);
+  if (first_mesh) {
+    print_mesh(N1, N2, box, charge);
+  }
 
   find_k_and_G_opt<<<(para.K0K1K2 - 1) / 64 + 1, 64>>>(
     para, 

src/force/pppm.cuh

@@ -70,6 +70,7 @@ private:
   void destroy_plans();
   void allocate_memory();
   void find_para(const int N, const Box& box);
+  void print_mesh(const int N1, const int N2, const Box& box, const GPU_Vector<float>& charge);
 
   bool need_peratom_virial = false;
   GPU_Vector<gpufftComplex> mesh_virial;

src/main_replica/force/force.cu

@@ -116,7 +116,7 @@ Force::~Force() = default;
 
 void Force::parse_potential(
   const char** param, const int num_param, const Box& box,
-  const int number_of_atoms, const bool use_pppm)
+  const int number_of_atoms, const Kspace_Setting& kspace)
 {
   if (num_param != 2)
     PRINT_INPUT_ERROR("potential should have one parameter.\n");
@@ -126,7 +126,7 @@ void Force::parse_potential(
   const std::string potential_name = read_potential_name(param[1]);
   potential_file_ = param[1];
   number_of_atoms_ = number_of_atoms;
-  use_pppm_ = use_pppm;
+  kspace_ = kspace;
   if (is_energy_nep(potential_name.c_str())) {
     model_type_ = Model_Type::energy_nep;
     potential_.reset(new NEP(param[1], number_of_atoms, false));
@@ -166,7 +166,8 @@ void Force::prepare_stream(const gpuStream_t stream)
       potential_file_.c_str(),
       number_of_atoms_,
       stream,
-      use_pppm_,
+      kspace_.method == "pppm",
+      kspace_.pppm_spacing,
       charge_potentials_.empty()));
     potential->N1 = 0;
     potential->N2 = number_of_atoms_;

src/main_replica/force/force.cuh

@@ -26,6 +26,11 @@
 class NEP;
 class NEP_Charge;
 
+struct Kspace_Setting {
+  std::string method = "pppm"; // pppm or ewald, and none for models without charges
+  double pppm_spacing = 0.0;   // requested PPPM spacing in A; 0 selects min(1, rc_radial / 6)
+};
+
 class Force
 {
 public:
@@ -34,7 +39,7 @@ public:
 
   void parse_potential(
     const char** param, int num_param, const Box& box,
-    int number_of_atoms, bool use_pppm);
+    int number_of_atoms, const Kspace_Setting& kspace);
 
   void clear();
   bool has_potential() const;
@@ -68,5 +73,5 @@ private:
   std::map<gpuStream_t, std::unique_ptr<NEP_Charge>> charge_potentials_;
   std::string potential_file_;
   int number_of_atoms_ = 0;
-  bool use_pppm_ = true;
+  Kspace_Setting kspace_;
 };

src/main_replica/force/nep_charge.cu

@@ -54,6 +54,7 @@ NEP_Charge::NEP_Charge(
   const int num_atoms,
   const gpuStream_t stream,
   const bool use_pppm_input,
+  const double pppm_spacing,
   const bool verbose)
   : use_pppm(use_pppm_input), stream_(stream), verbose_(verbose)
 {
@@ -335,7 +336,11 @@ NEP_Charge::NEP_Charge(
   // charge related parameters and data
   charge_para.alpha = float(PI) / paramb.rc_radial; // a good value
   if (use_pppm) {
-    pppm.initialize(charge_para.alpha, stream_);
+    // With alpha = pi / rc_radial, the default spacing keeps alpha * spacing at or below pi / 6.
+    pppm.initialize(
+      charge_para.alpha,
+      pppm_spacing > 0.0 ? pppm_spacing : std::fmin(1.0, paramb.rc_radial / 6.0),
+      verbose_);
   } else {
     ewald.initialize(charge_para.alpha);
   }

src/main_replica/force/nep_charge.cuh

@@ -125,6 +125,7 @@ public:
     const int num_atoms,
     const gpuStream_t stream,
     const bool use_pppm,
+    const double pppm_spacing,
     const bool verbose);
   virtual ~NEP_Charge(void);
   virtual void compute(