From 3b9df291c1b27017c74e557c58785ae80bf20d14 Mon Sep 17 00:00:00 2001 From: Eric Lindgren Date: Tue, 13 Dec 2022 15:19:59 +0100 Subject: [PATCH] Add functionality for running MD on an average PES --- src/force/force.cu | 118 +++++++++++++++++++++++++++++++--- src/force/force.cuh | 6 ++ src/main_gpumd/run.cu | 2 +- src/measure/dump_exyz.cu | 2 +- src/measure/dump_exyz.cuh | 2 +- src/measure/dump_observer.cu | 42 ++++++++---- src/measure/dump_observer.cuh | 8 ++- src/measure/measure.cu | 5 +- src/measure/measure.cuh | 2 +- 9 files changed, 157 insertions(+), 30 deletions(-) diff --git a/src/force/force.cu b/src/force/force.cu index 10ed3caf2..45e1b0c72 100644 --- a/src/force/force.cu +++ b/src/force/force.cu @@ -4,9 +4,7 @@ GPUMD is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - GPUMD is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of + (at your option) any later version. GPUMD is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License @@ -52,7 +50,7 @@ void Force::parse_potential( PRINT_INPUT_ERROR("reading error for potential file."); } int num_types = get_number_of_types(fid_potential); - + number_of_atoms_ = number_of_atoms; bool is_nep = false; // determine the potential if (strcmp(potential_name, "tersoff_1989") == 0) { @@ -384,6 +382,37 @@ static __global__ void gpu_apply_pbc(int N, Box box, double* g_x, double* g_y, d } } +static __global__ void gpu_add_vectors(int length, double* a, double* b) +{ + // Add vector b onto a + int n = blockIdx.x * blockDim.x + threadIdx.x; + if (n < length) { + a[n] += b[n]; + } +} + + +static __global__ void gpu_divide_and_copy_vector(int length, double* source, double* target, double denominator) +{ + // Divide source vector by denominator, and copy into target + int n = blockIdx.x * blockDim.x + threadIdx.x; + if (n < length) { + source[n] /= denominator; + target[n] = source[n]; + } +} + + +void Force::set_multiple_potentials_mode(const char *mode) +{ + multiple_potentials_mode_ = mode; + if(strcmp(multiple_potentials_mode_, "average") == 0){ + average_potential_per_atom_.resize(number_of_atoms_); + average_force_per_atom_.resize(number_of_atoms_ * 3); + average_virial_per_atom_.resize(number_of_atoms_ * 9); + } +} + void Force::compute( Box& box, GPU_Vector& position_per_atom, @@ -404,9 +433,45 @@ void Force::compute( number_of_atoms, force_per_atom.data(), force_per_atom.data() + number_of_atoms, force_per_atom.data() + number_of_atoms * 2, potential_per_atom.data(), virial_per_atom.data()); CUDA_CHECK_KERNEL - - potentials[0]->compute( - box, type, position_per_atom, potential_per_atom, force_per_atom, virial_per_atom); + + if(strcmp(multiple_potentials_mode_, "observe") == 0) + { + // If observing, calculate using main potential only + potentials[0]->compute( + box, type, position_per_atom, potential_per_atom, force_per_atom, virial_per_atom); + } + else if(strcmp(multiple_potentials_mode_, "average") == 0) + { + // Calculate average potential, force and virial per atom. + for (int i = 0; i < potentials.size(); i++){ + potentials[i]->compute( + box, type, position_per_atom, potential_per_atom, force_per_atom, virial_per_atom); + // Save properties to average vectors + gpu_add_vectors<<<(number_of_atoms - 1) / 128 + 1, 128>>>( + number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data()); + CUDA_CHECK_KERNEL + gpu_add_vectors<<<(3*number_of_atoms - 1) / 128 + 1, 128>>>( + 3*number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data()); + CUDA_CHECK_KERNEL + gpu_add_vectors<<<(9*number_of_atoms - 1) / 128 + 1, 128>>>( + 9*number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data()); + CUDA_CHECK_KERNEL + } + // Compute average and copy properties back into original vectors. + gpu_divide_and_copy_vector<<<(number_of_atoms - 1) / 128 + 1, 128>>>( + number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data(), potentials.size()); + CUDA_CHECK_KERNEL + gpu_divide_and_copy_vector<<<(number_of_atoms - 1) / 128 + 1, 128>>>( + 3*number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data(), potentials.size()); + CUDA_CHECK_KERNEL + gpu_divide_and_copy_vector<<<(number_of_atoms - 1) / 128 + 1, 128>>>( + 9*number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data(), potentials.size()); + CUDA_CHECK_KERNEL + } + else { + PRINT_INPUT_ERROR("Invalid mode for multiple potentials.\n"); + } + if (compute_hnemd_) { // the virial tensor: @@ -621,8 +686,43 @@ void Force::compute( force_per_atom.data() + number_of_atoms * 2, potential_per_atom.data(), virial_per_atom.data()); CUDA_CHECK_KERNEL - potentials[0]->compute( - box, type, position_per_atom, potential_per_atom, force_per_atom, virial_per_atom); + if(strcmp(multiple_potentials_mode_, "observe") == 0) + { + // If observing, calculate using main potential only + potentials[0]->compute( + box, type, position_per_atom, potential_per_atom, force_per_atom, virial_per_atom); + } + else if(strcmp(multiple_potentials_mode_, "average") == 0) + { + // Calculate average potential, force and virial per atom. + for (int i = 0; i < potentials.size(); i++){ + potentials[i]->compute( + box, type, position_per_atom, potential_per_atom, force_per_atom, virial_per_atom); + // Save properties to average vectors + gpu_add_vectors<<<(number_of_atoms - 1) / 128 + 1, 128>>>( + number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data()); + CUDA_CHECK_KERNEL + gpu_add_vectors<<<(3*number_of_atoms - 1) / 128 + 1, 128>>>( + 3*number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data()); + CUDA_CHECK_KERNEL + gpu_add_vectors<<<(9*number_of_atoms - 1) / 128 + 1, 128>>>( + 9*number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data()); + CUDA_CHECK_KERNEL + } + // Compute average and copy properties back into original vectors. + gpu_divide_and_copy_vector<<<(number_of_atoms - 1) / 128 + 1, 128>>>( + number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data(), potentials.size()); + CUDA_CHECK_KERNEL + gpu_divide_and_copy_vector<<<(number_of_atoms - 1) / 128 + 1, 128>>>( + 3*number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data(), potentials.size()); + CUDA_CHECK_KERNEL + gpu_divide_and_copy_vector<<<(number_of_atoms - 1) / 128 + 1, 128>>>( + 9*number_of_atoms, average_potential_per_atom_.data(), potential_per_atom.data(), potentials.size()); + CUDA_CHECK_KERNEL + } + else { + PRINT_INPUT_ERROR("Invalid mode for multiple potentials.\n"); + } if (compute_hnemd_) { // the virial tensor: diff --git a/src/force/force.cuh b/src/force/force.cuh index bf63be028..4b7b756d0 100644 --- a/src/force/force.cuh +++ b/src/force/force.cuh @@ -61,6 +61,7 @@ public: const std::vector& type, const std::vector& type_size, const double T); + void set_multiple_potentials_mode(const char* mode); bool compute_hnemd_ = false; int compute_hnemdec_ = -1; @@ -70,6 +71,11 @@ public: std::vector> potentials; private: + int number_of_atoms_ = -1; bool is_fcp = false; bool has_non_nep = false; + const char* multiple_potentials_mode_ = "observe"; // "observe" or "average" + GPU_Vector average_potential_per_atom_; + GPU_Vector average_force_per_atom_; + GPU_Vector average_virial_per_atom_; }; diff --git a/src/main_gpumd/run.cu b/src/main_gpumd/run.cu index 313799e00..26053bc9c 100644 --- a/src/main_gpumd/run.cu +++ b/src/main_gpumd/run.cu @@ -154,7 +154,7 @@ void Run::execute_run_in() void Run::perform_a_run() { integrate.initialize(N, time_step, group); - measure.initialize(number_of_steps, time_step, box, group, atom, force.potentials.size()); + measure.initialize(number_of_steps, time_step, box, group, atom, force); #ifdef USE_PLUMED if (measure.plmd.use_plumed == 1) { diff --git a/src/measure/dump_exyz.cu b/src/measure/dump_exyz.cu index 169cc95fb..a13d59903 100644 --- a/src/measure/dump_exyz.cu +++ b/src/measure/dump_exyz.cu @@ -219,7 +219,7 @@ void Dump_EXYZ::process( fflush(fid_); } -void Dump_EXYZ::setupObserverDump( +void Dump_EXYZ::setup_observer_dump( bool dump, int dump_interval, std::string file_label, diff --git a/src/measure/dump_exyz.cuh b/src/measure/dump_exyz.cuh index 93e19dc5f..16d88af45 100644 --- a/src/measure/dump_exyz.cuh +++ b/src/measure/dump_exyz.cuh @@ -40,7 +40,7 @@ public: GPU_Vector& gpu_thermo, const int file_index); void postprocess(); - void setupObserverDump( + void setup_observer_dump( bool dump, int dump_interval, std::string file_label, diff --git a/src/measure/dump_observer.cu b/src/measure/dump_observer.cu index a5ea12f80..a5dfe4298 100644 --- a/src/measure/dump_observer.cu +++ b/src/measure/dump_observer.cu @@ -25,33 +25,36 @@ Dump energy/force/virial with all loaded potentials at a given interval. #include "utilities/read_file.cuh" #include #include -#include -#include -#include void Dump_Observer::parse(const char** param, int num_param) { dump_ = true; printf("Dump observer.\n"); - if (num_param != 2) { + if (num_param != 3) { PRINT_INPUT_ERROR("dump_observer should have 2 parameters."); } - if (!is_valid_int(param[1], &dump_interval_)) { + mode_ = param[1]; + if (strcmp(mode_, "observe") != 0 && strcmp(mode_, "average") != 0) { + PRINT_INPUT_ERROR("observer mode should be 'observe' or 'average'"); + } + if (!is_valid_int(param[2], &dump_interval_)) { PRINT_INPUT_ERROR("observer dump interval should be an integer."); } if (dump_interval_ <= 0) { PRINT_INPUT_ERROR("observer dump interval should > 0."); } + printf(" every %d steps.\n", dump_interval_); } -void Dump_Observer::preprocess(const int number_of_atoms, const int number_of_potentials) +void Dump_Observer::preprocess(const int number_of_atoms, const int number_of_potentials, Force& force) { // Setup a dump_exyz with the dump_interval for dump_observer. - dump_exyz_.setupObserverDump(dump_, dump_interval_, "observer", 1, 1); + dump_exyz_.setup_observer_dump(dump_, dump_interval_, "observer", 1, 1); dump_exyz_.preprocess(number_of_atoms, number_of_potentials); + force.set_multiple_potentials_mode(mode_); } void Dump_Observer::process( @@ -68,14 +71,29 @@ void Dump_Observer::process( return; if ((step + 1) % dump_interval_ != 0) return; - const int number_of_potentials = force.potentials.size(); - for (int potential_index = 0; potential_index < number_of_potentials; potential_index++) { - force.potentials[potential_index]->compute(box, atom.type, atom.position_per_atom, - atom.potential_per_atom, atom.force_per_atom, atom.virial_per_atom); + if(strcmp(mode_, "observe") == 0) + { + // If observing, calculate properties with all potentials. + const int number_of_potentials = force.potentials.size(); + for (int potential_index = 0; potential_index < number_of_potentials; potential_index++) { + force.potentials[potential_index]->compute(box, atom.type, atom.position_per_atom, + atom.potential_per_atom, atom.force_per_atom, atom.virial_per_atom); + dump_exyz_.process( + step, global_time, box, atom.cpu_atom_symbol, atom.cpu_type, atom.position_per_atom, + atom.cpu_position_per_atom, atom.velocity_per_atom, atom.cpu_velocity_per_atom, + atom.force_per_atom, atom.virial_per_atom, thermo, potential_index); + } + } + else if(strcmp(mode_, "average") == 0) + { + // If average, dump already computed properties to file. dump_exyz_.process( step, global_time, box, atom.cpu_atom_symbol, atom.cpu_type, atom.position_per_atom, atom.cpu_position_per_atom, atom.velocity_per_atom, atom.cpu_velocity_per_atom, - atom.force_per_atom, atom.virial_per_atom, thermo, potential_index); + atom.force_per_atom, atom.virial_per_atom, thermo, 0); + } + else { + PRINT_INPUT_ERROR("Invalid observer mode.\n"); } } diff --git a/src/measure/dump_observer.cuh b/src/measure/dump_observer.cuh index 1112da99f..b481bb41b 100644 --- a/src/measure/dump_observer.cuh +++ b/src/measure/dump_observer.cuh @@ -28,7 +28,10 @@ class Dump_Observer { public: void parse(const char** param, int num_param); - void preprocess(const int number_of_atoms, const int number_of_files); + void preprocess( + const int number_of_atoms, + const int number_of_files, + Force& force); void process( int step, const double global_time, @@ -38,10 +41,9 @@ public: GPU_Vector& thermo); void postprocess(); - private: bool dump_ = false; int dump_interval_ = 1; - std::string mode_ = "observe"; // observe or average + const char* mode_ = "observe"; // observe or average Dump_EXYZ dump_exyz_; // Local member of dump_exyz to dump at a possibly different interval }; diff --git a/src/measure/measure.cu b/src/measure/measure.cu index 2f61b2ca0..19100b753 100644 --- a/src/measure/measure.cu +++ b/src/measure/measure.cu @@ -29,9 +29,10 @@ void Measure::initialize( Box& box, std::vector& group, Atom& atom, - const int number_of_potentials) + Force& force) { const int number_of_atoms = atom.mass.size(); + const int number_of_potentials = force.potentials.size(); dos.preprocess(time_step, group, atom.mass); sdc.preprocess(number_of_atoms, time_step, group); msd.preprocess(number_of_atoms, time_step, group); @@ -48,7 +49,7 @@ void Measure::initialize( dump_thermo.preprocess(); dump_force.preprocess(number_of_atoms, group); dump_exyz.preprocess(number_of_atoms, 1); - dump_observer.preprocess(number_of_atoms, number_of_potentials); + dump_observer.preprocess(number_of_atoms, number_of_potentials, force); #ifdef USE_NETCDF dump_netcdf.preprocess(number_of_atoms); #endif diff --git a/src/measure/measure.cuh b/src/measure/measure.cuh index 82cd4cfd8..0ea274055 100644 --- a/src/measure/measure.cuh +++ b/src/measure/measure.cuh @@ -54,7 +54,7 @@ public: Box& box, std::vector& group, Atom& atom, - const int number_of_potentials); + Force& force); void finalize( const int number_of_steps,