Add functionality for running MD on an average PES

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<double>& 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:

src/force/force.cuh

@@ -61,6 +61,7 @@ public:
     const std::vector<int>& type,
     const std::vector<int>& 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<std::unique_ptr<Potential>> 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<double> average_potential_per_atom_;
+  GPU_Vector<double> average_force_per_atom_;
+  GPU_Vector<double> average_virial_per_atom_;
 };

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) {

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, 

src/measure/dump_exyz.cuh

@@ -40,7 +40,7 @@ public:
     GPU_Vector<double>& gpu_thermo,
     const int file_index);
   void postprocess();
-  void setupObserverDump(
+  void setup_observer_dump(
     bool dump, 
     int dump_interval, 
     std::string file_label, 

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 <vector>
 #include <iostream>
-#include <execinfo.h>
-#include <stdio.h>
-#include <stdlib.h>
 
 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");
   }
 }
 

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<double>& 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
 };

src/measure/measure.cu

@@ -29,9 +29,10 @@ void Measure::initialize(
   Box& box,
   std::vector<Group>& 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

src/measure/measure.cuh

@@ -54,7 +54,7 @@ public:
     Box& box,
     std::vector<Group>& group,
     Atom& atom,
-    const int number_of_potentials);
+    Force& force);
 
   void finalize(
     const int number_of_steps,