diff --git a/doc/gpumd-replica/input_parameters/kspace.rst b/doc/gpumd-replica/input_parameters/kspace.rst index e65315ba7..06420fa51 100644 --- a/doc/gpumd-replica/input_parameters/kspace.rst +++ b/doc/gpumd-replica/input_parameters/kspace.rst @@ -10,14 +10,30 @@ Syntax :: - kspace + 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 diff --git a/doc/gpumd-replica/output_files/restart.rst b/doc/gpumd-replica/output_files/restart.rst index bdc90c4c7..8d5bd298f 100644 --- a/doc/gpumd-replica/output_files/restart.rst +++ b/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. diff --git a/doc/gpumd/input_parameters/kspace.rst b/doc/gpumd/input_parameters/kspace.rst index a466d2f9a..ef860736b 100644 --- a/doc/gpumd/input_parameters/kspace.rst +++ b/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). diff --git a/src/force/nep_charge.cu b/src/force/nep_charge.cu index cb150d3d9..6fbfa419f 100644 --- a/src/force/nep_charge.cu +++ b/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& tokens = line.tokens; if (!tokens.empty() && tokens[0] == "kspace") { diff --git a/src/force/nep_charge.cuh b/src/force/nep_charge.cuh index 72534ce32..dd71fc054 100644 --- a/src/force/nep_charge.cuh +++ b/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); diff --git a/src/force/pppm.cu b/src/force/pppm.cu index 4ba05a52c..b173cfc28 100644 --- a/src/force/pppm.cu +++ b/src/force/pppm.cu @@ -24,6 +24,7 @@ The k-space part of the PPPM method. #include #include #include +#include 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(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& 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 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& virial_per_atom, GPU_Vector& 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, diff --git a/src/force/pppm.cuh b/src/force/pppm.cuh index afe1c3755..6824abe57 100644 --- a/src/force/pppm.cuh +++ b/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& charge); bool need_peratom_virial = false; GPU_Vector mesh_virial; diff --git a/src/main_replica/force/force.cu b/src/main_replica/force/force.cu index 1e75e6e2e..e0ab057c7 100644 --- a/src/main_replica/force/force.cu +++ b/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_; diff --git a/src/main_replica/force/force.cuh b/src/main_replica/force/force.cuh index 0c8684762..d7332c6f2 100644 --- a/src/main_replica/force/force.cuh +++ b/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> charge_potentials_; std::string potential_file_; int number_of_atoms_ = 0; - bool use_pppm_ = true; + Kspace_Setting kspace_; }; diff --git a/src/main_replica/force/nep_charge.cu b/src/main_replica/force/nep_charge.cu index d6263d3eb..518f99139 100644 --- a/src/main_replica/force/nep_charge.cu +++ b/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); } diff --git a/src/main_replica/force/nep_charge.cuh b/src/main_replica/force/nep_charge.cuh index 0a807e12b..8f4ec95d2 100644 --- a/src/main_replica/force/nep_charge.cuh +++ b/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( diff --git a/src/main_replica/force/pppm.cu b/src/main_replica/force/pppm.cu index 3f4e4df2f..12c999e35 100644 --- a/src/main_replica/force/pppm.cu +++ b/src/main_replica/force/pppm.cu @@ -19,18 +19,39 @@ The k-space part of the PPPM method. #include "pppm.cuh" #include "utilities/common.cuh" +#include "utilities/error.cuh" #include "utilities/gpu_macro.cuh" #include -#include +#include #include +#include namespace{ -int get_best_K(const int m) +constexpr int max_mesh_points = 512 * 512 * 512; + +bool is_good_K(int n) +{ + const int primes[4] = {2, 3, 5, 7}; + for (int p : primes) { + while (n % p == 0) { + n /= p; + } + } + return n == 1; +} + +int get_best_K(const double required) { - int n = 16; - while (n < m) { - n *= 2; + int n = static_cast(std::ceil(std::fmin(required, max_mesh_points))); + if (n < 16) { + n = 16; + } + if (n % 2 != 0) { + ++n; + } + while (!is_good_K(n)) { + n += 2; } return n; } @@ -589,38 +610,37 @@ void PPPM::allocate_memory(const gpuStream_t stream) } } -void PPPM::initialize(const float alpha_input, const gpuStream_t stream) +void PPPM::initialize( + const float alpha_input, const double mesh_spacing_input, const bool verbose_input) { need_peratom_virial = false; + para = {}; + mesh_spacing = mesh_spacing_input; + verbose = verbose_input; para.alpha = alpha_input; para.alpha_factor = 0.25f / (para.alpha * para.alpha); - para.K[0] = 16; - para.K[1] = 16; - para.K[2] = 16; - para.K0K1K2 = para.K[0] * para.K[1] * para.K[2]; - allocate_memory(stream); } void PPPM::find_para(const int N, const Box& box, const gpuStream_t stream) { const float two_pi = 6.2831853f; - const double mesh_spacing = 1.0; // Is this good enough? const double volume = box.get_volume(); para.two_pi_over_V = two_pi / volume; - int K[3] = {0}; for (int d = 0; d < 3; ++d) { - const double box_thickness = volume / box.get_area(d); - K[d] = box_thickness / mesh_spacing; - K[d] = get_best_K(K[d]); - para.K_half[d] = K[d] / 2; - para.two_pi_over_K[d] = two_pi / K[d]; - } - para.K0K1 = K[0] * K[1]; - para.K0K1K2 = para.K0K1 * K[2]; - if (K[0] != para.K[0] || K[1] != para.K[1] || K[2] != para.K[2]) { - para.K[0] = K[0]; - para.K[1] = K[1]; - para.K[2] = K[2]; + const double required = volume / box.get_area(d) / mesh_spacing; + if (required > para.K[d]) { + para.K[d] = get_best_K(required); + } + para.K_half[d] = para.K[d] / 2; + para.two_pi_over_K[d] = two_pi / para.K[d]; + } + const double number_of_points = double(para.K[0]) * para.K[1] * para.K[2]; + if (number_of_points > max_mesh_points) { + PRINT_INPUT_ERROR("PPPM mesh is too large; increase spacing or reduce the box size."); + } + if (number_of_points != para.K0K1K2) { + para.K0K1 = para.K[0] * para.K[1]; + para.K0K1K2 = static_cast(number_of_points); allocate_memory(stream); } para.potential_factor = K_C_SP / N; @@ -631,6 +651,64 @@ void PPPM::find_para(const int N, const Box& box, const gpuStream_t stream) } } +void PPPM::print_mesh( + const int N1, + const int N2, + const Box& box, + const GPU_Vector& charge, + const gpuStream_t stream) +{ + 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 q(number_of_charges); + CHECK(gpuMemcpyAsync( + q.data(), + charge.data() + N1, + sizeof(float) * number_of_charges, + gpuMemcpyDeviceToHost, + stream)); + CHECK(gpuStreamSynchronize(stream)); + 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, @@ -644,7 +722,11 @@ void PPPM::find_force( GPU_Vector& potential_per_atom, const gpuStream_t stream) { + const bool first_mesh = !plan_initialized; find_para(N, box, stream); + if (first_mesh && verbose) { + print_mesh(N1, N2, box, charge, stream); + } find_k_and_G_opt<<<(para.K0K1K2 - 1) / 64 + 1, 64, 0, stream>>>( para, diff --git a/src/main_replica/force/pppm.cuh b/src/main_replica/force/pppm.cuh index 12b28e6b0..3fea4a75e 100644 --- a/src/main_replica/force/pppm.cuh +++ b/src/main_replica/force/pppm.cuh @@ -27,7 +27,8 @@ class PPPM public: PPPM(); ~PPPM(); - void initialize(const float alpha_input, const gpuStream_t stream); + void + initialize(const float alpha_input, const double mesh_spacing_input, const bool verbose_input); void find_force( const int N, const int N1, @@ -53,7 +54,9 @@ public: float two_pi_over_K[3]; // 2 * pi ./ K }; private: - Para para; + Para para = {}; + double mesh_spacing = 1.0; + bool verbose = true; GPU_Vector kx; GPU_Vector ky; GPU_Vector kz; @@ -67,6 +70,12 @@ private: bool plan_initialized = false; void allocate_memory(const gpuStream_t stream); void find_para(const int N, const Box& box, const gpuStream_t stream); + void print_mesh( + const int N1, + const int N2, + const Box& box, + const GPU_Vector& charge, + const gpuStream_t stream); void find_k_and_G(const double* box); bool need_peratom_virial = false; diff --git a/src/main_replica/main.cu b/src/main_replica/main.cu index 238577e4e..1cee95b64 100644 --- a/src/main_replica/main.cu +++ b/src/main_replica/main.cu @@ -13,11 +13,12 @@ along with GPUMD. If not, see . */ -#include "replica.cuh" +#include "force/force.cuh" #include "model/atom.cuh" #include "model/box.cuh" #include "model/group.cuh" #include "model/read_xyz.cuh" +#include "replica.cuh" #include "utilities/common.cuh" #include "utilities/error.cuh" #include "utilities/gpu_macro.cuh" @@ -37,7 +38,8 @@ struct Replica_Input { std::vector ensemble_command; std::vector multi_replica_command; std::vector dump_xyz_command; - std::vector kspace_command; + Kspace_Setting kspace; + bool has_kspace = false; double time_step = 1.0 / TIME_UNIT_CONVERSION; int number_of_steps = -1; }; @@ -101,10 +103,21 @@ Replica_Input parse_run_input( PRINT_INPUT_ERROR("Only one dump_xyz command is allowed in a replica run."); input.dump_xyz_command = command; } else if (keyword == "kspace") { - if (!input.kspace_command.empty() || command.size() != 2 || - (command[1] != "ewald" && command[1] != "pppm")) - PRINT_INPUT_ERROR("kspace should be specified once as ewald or pppm."); - input.kspace_command = command; + if (input.has_kspace) + PRINT_INPUT_ERROR("kspace can only appear once."); + if (command.size() < 2 || command.size() > 3) + PRINT_INPUT_ERROR("kspace requires ewald or pppm [spacing]."); + if (command[1] != "ewald" && command[1] != "pppm") + PRINT_INPUT_ERROR("kspace method can only be ewald or pppm."); + if (command[1] == "ewald" && command.size() == 3) + PRINT_INPUT_ERROR("kspace ewald does not accept spacing."); + if ( + command.size() == 3 && + (!is_valid_real(command[2].c_str(), &input.kspace.pppm_spacing) || + !(input.kspace.pppm_spacing >= 0.2 && input.kspace.pppm_spacing <= 2.0))) + PRINT_INPUT_ERROR("PPPM spacing must be a finite number between 0.2 and 2.0 A."); + input.kspace.method = command[1]; + input.has_kspace = true; } else if (keyword == "run") { if (input.number_of_steps >= 0 || command.size() != 2 || !is_valid_int(command[1].c_str(), &input.number_of_steps) || @@ -171,7 +184,7 @@ int main() input.ensemble_command, input.multi_replica_command, input.dump_xyz_command, - input.kspace_command.empty() ? "pppm" : input.kspace_command[1], + input.kspace, atom, box, group, diff --git a/src/main_replica/replica.cu b/src/main_replica/replica.cu index 6ee8f5bf9..d4ee8db69 100644 --- a/src/main_replica/replica.cu +++ b/src/main_replica/replica.cu @@ -23,7 +23,7 @@ void run_replica_md( const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -40,7 +40,7 @@ void run_replica_md( ensemble_command, multi_replica_command, dump_xyz_command, - kspace_method, + kspace, source_atom, source_box, source_group, @@ -52,7 +52,7 @@ void run_replica_md( ensemble_command, multi_replica_command, dump_xyz_command, - kspace_method, + kspace, source_atom, source_box, source_group, diff --git a/src/main_replica/replica.cuh b/src/main_replica/replica.cuh index 0ff7221de..131786168 100644 --- a/src/main_replica/replica.cuh +++ b/src/main_replica/replica.cuh @@ -20,6 +20,7 @@ class Atom; class Box; +struct Kspace_Setting; class Group; namespace replica @@ -29,7 +30,7 @@ void run_remd( const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -41,7 +42,7 @@ void run_prd( const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -54,7 +55,7 @@ void run_replica_md( const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, diff --git a/src/main_replica/replica_common.cu b/src/main_replica/replica_common.cu index 6233205c8..df8b8a3a3 100644 --- a/src/main_replica/replica_common.cu +++ b/src/main_replica/replica_common.cu @@ -277,13 +277,13 @@ Device_Context::~Device_Context() Replica_Runtime::Replica_Runtime( const std::vector>& potential_commands, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, const double time_step) : potential_commands_(potential_commands) - , kspace_method_(kspace_method) + , kspace_(kspace) , source_atom_(source_atom) , source_box_(source_box) , source_group_(source_group) @@ -327,7 +327,7 @@ void Replica_Runtime::initialize_device_force(Device_Context& device) static_cast(parameters.size()), source_box_, source_atom_.number_of_atoms, - kspace_method_ == "pppm"); + kspace_); if (!device.force.has_potential()) PRINT_INPUT_ERROR( "The selected potential is not an explicit-stream NEP energy model."); diff --git a/src/main_replica/replica_common.cuh b/src/main_replica/replica_common.cuh index 0da3845f4..41ce7f824 100644 --- a/src/main_replica/replica_common.cuh +++ b/src/main_replica/replica_common.cuh @@ -132,7 +132,7 @@ public: Replica_Runtime( const std::vector>& potential_commands, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -157,7 +157,7 @@ public: private: const std::vector>& potential_commands_; - std::string kspace_method_; + Kspace_Setting kspace_; Atom& source_atom_; Box& source_box_; std::vector& source_group_; diff --git a/src/main_replica/replica_prd.cu b/src/main_replica/replica_prd.cu index 79631cd01..14ffc60e0 100644 --- a/src/main_replica/replica_prd.cu +++ b/src/main_replica/replica_prd.cu @@ -60,7 +60,7 @@ using replica::values_match; constexpr const char* event_output_marker = "# GPUMD PRD event output format 1"; -constexpr int restart_format_version = 6; +constexpr int restart_format_version = 7; void validate_append_event_file(const std::string& filename) { @@ -383,7 +383,7 @@ public: const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -397,7 +397,7 @@ private: PRD_Config config_; const std::vector>& potential_commands_; const std::vector& ensemble_command_; - std::string kspace_method_; + Kspace_Setting kspace_; Atom& source_atom_; Box& source_box_; std::vector& source_group_; @@ -464,7 +464,7 @@ PRD_Driver::PRD_Driver( const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -472,7 +472,7 @@ PRD_Driver::PRD_Driver( const int number_of_steps) : potential_commands_(potential_commands) , ensemble_command_(ensemble_command) - , kspace_method_(kspace_method) + , kspace_(kspace) , source_atom_(source_atom) , source_box_(source_box) , source_group_(source_group) @@ -480,7 +480,7 @@ PRD_Driver::PRD_Driver( , number_of_steps_(number_of_steps) , runtime_( potential_commands, - kspace_method, + kspace, source_atom, source_box, source_group, @@ -551,7 +551,7 @@ void PRD_Driver::validate_input() (!source_box_.pbc_x || !source_box_.pbc_y || !source_box_.pbc_z)) PRINT_INPUT_ERROR("qNEP PRD requires periodic boundaries in all directions."); if (!is_supported_charge_nep_name(potential_name)) - kspace_method_ = "none"; + kspace_ = {"none", 0.0}; potential_hash_ = hash_file(potential_commands_[0][1]); double final_temperature = 0.0; @@ -1286,7 +1286,8 @@ void PRD_Driver::write_restart_metadata() const "temperature_coupling %.17g\n", temperature_coupling_); fprintf(file, "potential_hash %llu\n", potential_hash_); - fprintf(file, "kspace %s\n", kspace_method_.c_str()); + fprintf(file, "kspace %s\n", kspace_.method.c_str()); + fprintf(file, "pppm_spacing %.17g\n", kspace_.pppm_spacing); fprintf(file, "event_interval %d\n", config_.event_interval); fprintf(file, "dephase_steps %d\n", config_.dephase_steps); fprintf(file, "correlation_steps %d\n", config_.correlation_steps); @@ -1505,8 +1506,11 @@ void PRD_Driver::read_restart_metadata() PRINT_INPUT_ERROR("PRD restart potential file does not match the input."); std::string stored_kspace; read_restart_value(input, "kspace", stored_kspace); - if (stored_kspace != kspace_method_) + if (stored_kspace != kspace_.method) PRINT_INPUT_ERROR("PRD restart kspace does not match the input."); + read_restart_value(input, "pppm_spacing", value); + if (!values_match(value, kspace_.pppm_spacing)) + PRINT_INPUT_ERROR("PRD restart pppm_spacing does not match the input."); int integer_value = 0; read_restart_value(input, "event_interval", integer_value); @@ -1668,7 +1672,7 @@ void run_prd( const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -1680,7 +1684,7 @@ void run_prd( ensemble_command, multi_replica_command, dump_xyz_command, - kspace_method, + kspace, source_atom, source_box, source_group, diff --git a/src/main_replica/replica_remd.cu b/src/main_replica/replica_remd.cu index a9a63965b..82d0bc045 100644 --- a/src/main_replica/replica_remd.cu +++ b/src/main_replica/replica_remd.cu @@ -56,7 +56,7 @@ using replica::validate_xyz_box; using replica::values_match; constexpr double target_acceptance = 0.2; -constexpr int restart_format_version = 5; +constexpr int restart_format_version = 6; constexpr const char* exchange_output_marker = "# GPUMD REMD exchange output format 1"; constexpr const char* status_output_marker = @@ -270,7 +270,7 @@ public: const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -284,7 +284,7 @@ private: REMD_Config config_; const std::vector>& potential_commands_; const std::vector& ensemble_command_; - std::string kspace_method_; + Kspace_Setting kspace_; Atom& source_atom_; Box& source_box_; std::vector& source_group_; @@ -334,7 +334,7 @@ REMD_Driver::REMD_Driver( const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -342,13 +342,13 @@ REMD_Driver::REMD_Driver( const int number_of_steps) : potential_commands_(potential_commands) , ensemble_command_(ensemble_command) - , kspace_method_(kspace_method) + , kspace_(kspace) , source_atom_(source_atom) , source_box_(source_box) , source_group_(source_group) , time_step_(time_step) , number_of_steps_(number_of_steps) - , runtime_(potential_commands, kspace_method, source_atom, source_box, source_group, time_step) + , runtime_(potential_commands, kspace, source_atom, source_box, source_group, time_step) { config_.parse(multi_replica_command, dump_xyz_command); } @@ -375,7 +375,7 @@ void REMD_Driver::validate_input() (!source_box_.pbc_x || !source_box_.pbc_y || !source_box_.pbc_z)) PRINT_INPUT_ERROR("qNEP REMD requires periodic boundaries in all directions."); if (!is_supported_charge_nep_name(potential_name)) - kspace_method_ = "none"; + kspace_ = {"none", 0.0}; potential_hash_ = hash_file(potential_commands_[0][1]); if (ensemble_command_.size() != 5 || ensemble_command_[0] != "ensemble" || @@ -858,7 +858,8 @@ void REMD_Driver::write_restart_metadata() const fprintf(file, "time_step %.17g\n", time_step_); fprintf(file, "temperature_coupling %.17g\n", temperature_coupling_); fprintf(file, "potential_hash %llu\n", potential_hash_); - fprintf(file, "kspace %s\n", kspace_method_.c_str()); + fprintf(file, "kspace %s\n", kspace_.method.c_str()); + fprintf(file, "pppm_spacing %.17g\n", kspace_.pppm_spacing); fprintf(file, "temperatures"); for (const double temperature : config_.temperatures) fprintf(file, " %.17g", temperature); @@ -934,6 +935,7 @@ void REMD_Driver::read_restart_metadata() bool has_coupling = false; bool has_hash = false; bool has_kspace = false; + bool has_pppm_spacing = false; bool has_temperatures = false; bool has_replica_mapping = false; bool has_temperature_mapping = false; @@ -992,9 +994,14 @@ void REMD_Driver::read_restart_metadata() has_hash = true; } else if (key == "kspace") { std::string value; - if (has_kspace || !(input >> value) || value != kspace_method_) + if (has_kspace || !(input >> value) || value != kspace_.method) PRINT_INPUT_ERROR("REMD restart kspace does not match the input."); has_kspace = true; + } else if (key == "pppm_spacing") { + double value = 0.0; + if (has_pppm_spacing || !(input >> value) || !values_match(value, kspace_.pppm_spacing)) + PRINT_INPUT_ERROR("REMD restart pppm_spacing does not match the input."); + has_pppm_spacing = true; } else if (key == "temperatures") { if (has_temperatures) PRINT_INPUT_ERROR("REMD restart repeats temperatures."); @@ -1068,12 +1075,12 @@ void REMD_Driver::read_restart_metadata() } } - if (!has_version || !has_mode || !has_replicas || !has_completed_steps || - !has_exchange_count || !has_exchange_interval || !has_time_step || - !has_coupling || !has_hash || !has_kspace || !has_temperatures || !has_replica_mapping || - !has_temperature_mapping || !has_attempts || !has_accepts || - !has_probability || !has_exchange_rng || - std::find(has_bdp.begin(), has_bdp.end(), false) != has_bdp.end()) + if ( + !has_version || !has_mode || !has_replicas || !has_completed_steps || !has_exchange_count || + !has_exchange_interval || !has_time_step || !has_coupling || !has_hash || !has_kspace || + !has_pppm_spacing || !has_temperatures || !has_replica_mapping || !has_temperature_mapping || + !has_attempts || !has_accepts || !has_probability || !has_exchange_rng || + std::find(has_bdp.begin(), has_bdp.end(), false) != has_bdp.end()) PRINT_INPUT_ERROR("REMD restart metadata is incomplete."); if (exchange_count_ != completed_steps_ / config_.exchange_interval) PRINT_INPUT_ERROR("REMD restart exchange_count is inconsistent with completed_steps."); @@ -1282,7 +1289,7 @@ void replica::run_remd( const std::vector& ensemble_command, const std::vector& multi_replica_command, const std::vector& dump_xyz_command, - const std::string& kspace_method, + const Kspace_Setting& kspace, Atom& source_atom, Box& source_box, std::vector& source_group, @@ -1294,7 +1301,7 @@ void replica::run_remd( ensemble_command, multi_replica_command, dump_xyz_command, - kspace_method, + kspace, source_atom, source_box, source_group, diff --git a/tests_pytest/test_kspace_consistency.py b/tests_pytest/test_kspace_consistency.py index 1706c0ad4..7b78e453e 100644 --- a/tests_pytest/test_kspace_consistency.py +++ b/tests_pytest/test_kspace_consistency.py @@ -14,9 +14,13 @@ existing bulk_water and bulk_perovskite fixtures, plus a larger supercell built by tiling bulk_perovskite further (no new fixture file needed). """ +import re +import subprocess + import numpy as np import pytest from calorine.calculators import GPUNEP +from calorine.gpumd import write_xyz from conftest import MODELS_DIR, approx_tol, make_bulk_perovskite, make_bulk_water @@ -72,3 +76,74 @@ def test_pppm_agrees_with_ewald(structure_name, model_type, gpumd_command): assert energy_pppm == approx_tol(energy_ewald, PPPM_VS_EWALD_ENERGY_TOLERANCE) assert np.allclose(forces_pppm, forces_ewald, **PPPM_VS_EWALD_FORCE_TOLERANCE) + + +def _run_gpumd(directory, model_path, atoms, kspace_before_run, kspace_after_run, gpumd_command): + write_xyz(str(directory / 'model.xyz'), atoms) + lines = [f'potential {model_path}'] + lines += [f'kspace {kspace_before_run}'] if kspace_before_run else [] + lines += ['ensemble nve', 'time_step 1', 'run 1'] + lines += [f'kspace {kspace_after_run}'] if kspace_after_run else [] + (directory / 'run.in').write_text('\n'.join(lines) + '\n') + result = subprocess.run( + [gpumd_command], cwd=directory, capture_output=True, text=True, check=False + ) + assert result.returncode == 0, result.stdout + result.stderr + return result.stdout + + +def _pppm_target_spacing(directory, rc_radial, kspace_before_run, kspace_after_run, gpumd_command): + model_text = (MODELS_DIR / 'qnep_mode1_water.txt').read_text() + model_path = directory / 'nep.txt' + model_path.write_text(model_text.replace('cutoff 4 4 ', f'cutoff {rc_radial} 4 ', 1)) + stdout = _run_gpumd( + directory, model_path, make_bulk_water(), kspace_before_run, kspace_after_run, + gpumd_command + ) + return float(re.search(r'PPPM mesh: .*\(target spacing (\S+) A', stdout).group(1)) + + +# The default spacing is min(1, rc_radial / 6) A. +# With alpha = pi / rc_radial, it keeps the product of alpha and the spacing at or below pi / 6. +@pytest.mark.parametrize( + 'rc_radial, kspace_before_run, kspace_after_run, expected_spacing', + [ + (4, None, None, 4 / 6), + (4, 'pppm', None, 4 / 6), + (6, 'pppm', None, 1.0), + (8, 'pppm', None, 1.0), + (4, 'pppm 0.9', None, 0.9), + (4, None, 'pppm 0.9', 0.9), + ], + ids=['absent_rc4', 'default_rc4', 'default_rc6', 'default_rc8', 'explicit', 'explicit_late'], +) +def test_pppm_spacing( + tmp_path, gpumd_command, rc_radial, kspace_before_run, kspace_after_run, expected_spacing +): + spacing = _pppm_target_spacing( + tmp_path, rc_radial, kspace_before_run, kspace_after_run, gpumd_command + ) + assert spacing == pytest.approx(expected_spacing, rel=1e-12) + + +# The ratio of the estimate to the root-mean-square force difference is 0.97, 1.13 and 1.20 +# at spacings of 0.5, 0.7 and 1.0 A for this water box of 1701 atoms. +PPPM_FORCE_ERROR_ESTIMATE_TOLERANCE = 0.25 + + +@pytest.mark.parametrize('spacing', [0.5, 0.7, 1.0]) +def test_pppm_force_error_estimate(tmp_path, gpumd_command, spacing): + atoms = make_bulk_water().repeat(3) + model_path = MODELS_DIR / 'qnep_mode1_water.txt' + stdout = _run_gpumd(tmp_path, model_path, atoms, f'pppm {spacing}', None, gpumd_command) + estimate = float(re.search(r'PPPM force error estimate: (\S+) eV/A', stdout).group(1)) + + forces = {} + for kspace in [f'pppm {spacing}', 'ewald']: + calc = GPUNEP(str(model_path), command=gpumd_command) + calc.single_point_parameters = calc.single_point_parameters + [('kspace', kspace)] + atoms.calc = calc + forces[kspace] = atoms.get_forces() + difference = forces[f'pppm {spacing}'] - forces['ewald'] + rms_difference = np.sqrt(np.mean(np.sum(difference**2, axis=1))) + assert estimate == pytest.approx(rms_difference, rel=PPPM_FORCE_ERROR_ESTIMATE_TOLERANCE) diff --git a/tests_regression/pppm_checks.py b/tests_regression/pppm_checks.py index e9890f9e7..c799c341b 100644 --- a/tests_regression/pppm_checks.py +++ b/tests_regression/pppm_checks.py @@ -12,6 +12,7 @@ rf"PPPM mesh: (\d+) x (\d+) x (\d+) \(target spacing ({NUMBER}) A; " rf"actual spacing ({NUMBER}) ({NUMBER}) ({NUMBER}) A\)\." ) +ESTIMATE = re.compile(rf"PPPM force error estimate: ({NUMBER}) eV/A\.") def validate_spec(spec): @@ -44,12 +45,17 @@ def good_size(n): def check(spec, result): text = Path(result["stdout_path"]).read_text(encoding="utf-8") initial = [] + estimates = [] for line in text.splitlines(): if line.startswith("PPPM "): match = INITIAL.fullmatch(line) - if not match: + estimate = ESTIMATE.fullmatch(line) + if match: + initial.append(match.groups()) + elif estimate and len(initial) == 1 and not estimates: + estimates.append(float(estimate.group(1))) + else: raise ValueError("unexpected or malformed PPPM mesh output") - initial.append(match.groups()) expected = spec["initial_mesh"] if expected is None: if initial: @@ -57,6 +63,8 @@ def check(spec, result): return {"initial_mesh": None} if len(initial) != 1: raise ValueError(f"expected one initial mesh line, got {len(initial)}") + if not all(math.isfinite(e) and e >= 0 for e in estimates): + raise ValueError("invalid PPPM force error estimate") mesh = [int(x) for x in initial[0][:3]] spacing = float(initial[0][3]) actual = [float(x) for x in initial[0][4:]] diff --git a/tests_regression/replica/README.md b/tests_regression/replica/README.md index 387b43e0d..8ddea21ff 100644 --- a/tests_regression/replica/README.md +++ b/tests_regression/replica/README.md @@ -1,6 +1,6 @@ # Replica regression tests -These 12 basic integration tests exercise the standalone `gpumd_replica` +These 15 basic integration tests exercise the standalone `gpumd_replica` executable. They run independently of the parent directory's 221-case GPUMD/NEP differential suite and do not require a historical executable. @@ -32,6 +32,8 @@ The tests check: every replica's thermostat RNG state; - rejection of unsupported restart versions and changed exchange/dephasing settings; - qNEP Ewald/PPPM selection, inline comments, default solver, and restart validation; +- the default and explicit PPPM spacing, its input validation, its restart validation, + and the PPPM mesh and force error estimate lines; - single-segment dephasing syntax and the retry limit; - rejection of the unsupported `verbose` and `verbose_output` options in both modes. diff --git a/tests_regression/replica/test_replica.py b/tests_regression/replica/test_replica.py index 0a286bf4b..137ad58a6 100644 --- a/tests_regression/replica/test_replica.py +++ b/tests_regression/replica/test_replica.py @@ -6,6 +6,7 @@ import os from pathlib import Path +import re import shutil import subprocess import tempfile @@ -226,6 +227,53 @@ def test_qnep_comments_and_restart_solver_validation(self): # Omitting kspace must resolve to the same default. self.run_md(directory, mode, 1, options, resume=True) + def test_qnep_pppm_spacing_and_force_error_estimate(self): + options = "replicas 2 exchange 2 temp 300 400" + # The default spacing is min(1, rc_radial / 6) A. + cases = [(8, "", 1.0), (8, "pppm", 1.0), (4, "pppm", 4 / 6), (4, "pppm 0.9", 0.9)] + for rc_radial, kspace, expected in cases: + with self.subTest(rc_radial=rc_radial, kspace=kspace): + directory = self.prepare(f"spacing_{rc_radial}_{kspace.replace(' ', '_')}", + charge=True) + model = directory / "nep.txt" + model.write_text(model.read_text().replace( + "cutoff 8 4 ", f"cutoff {rc_radial} 4 ", 1)) + output = self.run_md(directory, "remd", 1, options, kspace=kspace).stdout + spacings = re.findall(r"PPPM mesh: .*\(target spacing (\S+) A", output) + self.assertEqual(len(spacings), 1, output) + self.assertAlmostEqual(float(spacings[0]), expected, delta=1e-12) + estimates = re.findall(r"PPPM force error estimate: (\S+) eV/A\.", output) + self.assertEqual(len(estimates), 1, output) + self.assertGreater(float(estimates[0]), 0.0) + + def test_qnep_pppm_spacing_input_validation(self): + options = "replicas 2 exchange 2 temp 300 400" + cases = [ + ("pppm 0.1", "PPPM spacing must be a finite number between 0.2 and 2.0 A."), + ("pppm 1.0oops", "PPPM spacing must be a finite number between 0.2 and 2.0 A."), + ("ewald 1.0", "kspace ewald does not accept spacing."), + ("pppm 1.0 2.0", "kspace requires ewald or pppm [spacing]."), + ] + for kspace, error in cases: + with self.subTest(kspace=kspace): + directory = self.prepare(f"invalid_{kspace.replace(' ', '_')}", charge=True) + self.run_md(directory, "remd", 1, options, kspace=kspace, expected_error=error) + + def test_restart_pppm_spacing_validation(self): + for mode in ("remd", "prd"): + options = ("replicas 2 exchange 2 temp 300 400" if mode == "remd" else + "replicas 2 event 2 dephase 0 correlate 5 distance 0.5 quench 0.001 2000") + with self.subTest(mode=mode): + directory = self.prepare(f"{mode}_spacing", charge=True) + self.run_md(directory, mode, 1, options, kspace="pppm 0.9") + self.assertEqual(float(metadata(directory, mode)["pppm_spacing"]), 0.9) + self.run_md(directory, mode, 1, options, kspace="pppm 0.9", resume=True) + before = (directory / f"{mode}_restart.meta").read_bytes() + for kspace in ("pppm 0.8", "pppm", ""): + self.run_md(directory, mode, 1, options, kspace=kspace, resume=True, + expected_error=f"{mode.upper()} restart pppm_spacing does not match") + self.assertEqual((directory / f"{mode}_restart.meta").read_bytes(), before) + def test_split_runs_preserve_state(self): for mode in ("remd", "prd"): with self.subTest(mode=mode): diff --git a/tests_regression/test_pppm_checks.py b/tests_regression/test_pppm_checks.py index bf0f9367b..cdfabccbe 100644 --- a/tests_regression/test_pppm_checks.py +++ b/tests_regression/test_pppm_checks.py @@ -6,6 +6,7 @@ INITIAL = "PPPM mesh: 16 x 16 x 16 (target spacing 1 A; actual spacing 0.99375 0.5 0.5 A).\n" +ESTIMATE = "PPPM force error estimate: 1.234e-04 eV/A.\n" class PPPMCheckTests(unittest.TestCase): @@ -23,6 +24,15 @@ def check(self, text): def test_initial_mesh(self): self.assertEqual(self.check(INITIAL), {"initial_mesh": [16, 16, 16]}) + def test_force_error_estimate(self): + self.assertEqual(self.check(INITIAL + ESTIMATE), {"initial_mesh": [16, 16, 16]}) + with self.assertRaisesRegex(ValueError, "unexpected"): + self.check(ESTIMATE + INITIAL) + with self.assertRaisesRegex(ValueError, "unexpected"): + self.check(INITIAL + ESTIMATE + ESTIMATE) + with self.assertRaisesRegex(ValueError, "invalid PPPM force error"): + self.check(INITIAL + ESTIMATE.replace("1.234e-04", "-1.234e-04")) + def test_unexpected_intermediate_mesh_output_fails(self): with self.assertRaisesRegex(ValueError, "unexpected"): self.check(INITIAL + "PPPM mesh update: 18 x 16 x 16\n")