Merge pull request brucefan1983#1158 from brucefan1983/add-efield-nep-charge

Add efield nep charge

Author: brucefan1983

src/force/nep_charge.cu

@@ -2124,3 +2124,5 @@ void NEP_Charge::compute(
 const GPU_Vector<int>& NEP_Charge::get_NN_radial_ptr() { return nep_data.NN_radial; }
 
 const GPU_Vector<int>& NEP_Charge::get_NL_radial_ptr() { return nep_data.NL_radial; }
+
+GPU_Vector<float>& NEP_Charge::get_charge_reference() { return nep_data.charge; }

src/force/nep_charge.cuh

@@ -137,6 +137,8 @@ public:
 
   const GPU_Vector<int>& get_NL_radial_ptr();
 
+  GPU_Vector<float>& get_charge_reference();
+
 private:
   ParaMB paramb;
   ANN annmb;

src/force/potential.cuh

@@ -75,6 +75,12 @@ public:
     return dummy_NL; // Return the const reference to NL_radial
   }
 
+  virtual GPU_Vector<float>& get_charge_reference()
+  {
+    static GPU_Vector<float> dummy_charge;
+    return dummy_charge;
+  }
+
 protected:
   void find_properties_many_body(
     Box& box,

src/main_gpumd/add_efield.cu

@@ -18,8 +18,10 @@ Add electric field to a group of atoms.
 ------------------------------------------------------------------------------*/
 
 #include "add_efield.cuh"
+#include "force/force.cuh"
 #include "model/atom.cuh"
 #include "model/group.cuh"
+#include "utilities/gpu_vector.cuh"
 #include "utilities/gpu_macro.cuh"
 #include "utilities/read_file.cuh"
 #include <iostream>
@@ -48,7 +50,30 @@ static void __global__ add_efield(
   }
 }
 
-void Add_Efield::compute(const int step, const std::vector<Group>& groups, Atom& atom)
+// for NEP-charge
+static void __global__ add_efield(
+  const int group_size,
+  const int group_size_sum,
+  const int* g_group_contents,
+  const double Ex,
+  const double Ey,
+  const double Ez,
+  const float* g_charge,
+  double* g_fx,
+  double* g_fy,
+  double* g_fz)
+{
+  const int tid = blockIdx.x * blockDim.x + threadIdx.x;
+  if (tid < group_size) {
+    const int atom_id = g_group_contents[group_size_sum + tid];
+    const double charge = g_charge[atom_id];
+    g_fx[atom_id] += charge * Ex;
+    g_fy[atom_id] += charge * Ey;
+    g_fz[atom_id] += charge * Ez;
+  }
+}
+
+void Add_Efield::compute(const int step, const std::vector<Group>& groups, Atom& atom, Force& force)
 {
   for (int call = 0; call < num_calls_; ++call) {
     const int step_mod_table_length = step % table_length_[call];
@@ -58,17 +83,33 @@ void Add_Efield::compute(const int step, const std::vector<Group>& groups, Atom&
     const int num_atoms_total = atom.force_per_atom.size() / 3;
     const int group_size = groups[grouping_method_[call]].cpu_size[group_id_[call]];
     const int group_size_sum = groups[grouping_method_[call]].cpu_size_sum[group_id_[call]];
-    add_efield<<<(group_size - 1) / 64 + 1, 64>>>(
-      group_size,
-      group_size_sum,
-      groups[grouping_method_[call]].contents.data(),
-      Ex,
-      Ey,
-      Ez,
-      atom.charge.data(),
-      atom.force_per_atom.data(),
-      atom.force_per_atom.data() + num_atoms_total,
-      atom.force_per_atom.data() + num_atoms_total * 2);
+    if (is_nep_charge) {
+      GPU_Vector<float>& nep_charge = force.potentials[0]->get_charge_reference();
+      add_efield<<<(group_size - 1) / 64 + 1, 64>>>(
+        group_size,
+        group_size_sum,
+        groups[grouping_method_[call]].contents.data(),
+        Ex,
+        Ey,
+        Ez,
+        nep_charge.data(),
+        atom.force_per_atom.data(),
+        atom.force_per_atom.data() + num_atoms_total,
+        atom.force_per_atom.data() + num_atoms_total * 2);
+    }
+    else {
+      add_efield<<<(group_size - 1) / 64 + 1, 64>>>(
+        group_size,
+        group_size_sum,
+        groups[grouping_method_[call]].contents.data(),
+        Ex,
+        Ey,
+        Ez,
+        atom.charge.data(),
+        atom.force_per_atom.data(),
+        atom.force_per_atom.data() + num_atoms_total,
+        atom.force_per_atom.data() + num_atoms_total * 2);
+    }
     GPU_CHECK_KERNEL
   }
 }
@@ -160,6 +201,14 @@ void Add_Efield::parse(const char** param, int num_param, const std::vector<Grou
   if (num_calls_ > 10) {
     PRINT_INPUT_ERROR("add_efield cannot be used more than 10 times in one run.");
   }
+
+  is_nep_charge = check_is_nep_charge();
+  if (is_nep_charge) {
+    printf("    using the charge values predicted by the NEP-Charge model.\n");
+  } else {
+    printf("    using the charge values specified in model.xyz.\n");
+  }
+
 }
 
 void Add_Efield::finalize() { num_calls_ = 0; }

src/main_gpumd/add_efield.cuh

@@ -19,12 +19,13 @@
 
 class Atom;
 class Group;
+class Force;
 
 class Add_Efield
 {
 public:
   void parse(const char** param, int num_param, const std::vector<Group>& group);
-  void compute(const int step, const std::vector<Group>& groups, Atom& atom);
+  void compute(const int step, const std::vector<Group>& groups, Atom& atom, Force& force);
   void finalize();
 
 private:
@@ -33,4 +34,5 @@ private:
   std::vector<double> efield_table_[10];
   int grouping_method_[10];
   int group_id_[10];
+  bool is_nep_charge = false;
 };

src/main_gpumd/run.cu

@@ -293,7 +293,7 @@ void Run::perform_a_run()
     electron_stop.compute(time_step, atom);
     add_force.compute(step, group, atom);
     add_random_force.compute(step, atom);
-    add_efield.compute(step, group, atom);
+    add_efield.compute(step, group, atom, force);
 
     integrate.compute2(time_step, double(step) / number_of_steps, group, box, atom, thermo, force);
 

src/utilities/read_file.cu

@@ -56,3 +56,53 @@ int is_valid_real(const char* s, double* result)
     return 1;
   }
 }
+
+static std::string get_potential_file_name()
+{
+  std::ifstream input_run("run.in");
+  if (!input_run.is_open()) {
+    PRINT_INPUT_ERROR("Cannot open run.in.");
+  }
+  std::string potential_file_name;
+  std::string line;
+  while (std::getline(input_run, line)) {
+    std::vector<std::string> tokens = get_tokens(line);
+    if (tokens.size() != 0) {
+      if (tokens[0] == "potential") {
+        potential_file_name = tokens[1];
+        break;
+      }
+    }
+  }
+
+  input_run.close();
+  return potential_file_name;
+}
+
+bool check_is_nep_charge()
+{
+  bool is_nep_charge = false;
+  std::string potential_file_name = get_potential_file_name();
+
+  std::ifstream input_potential(potential_file_name);
+  if (!input_potential.is_open()) {
+    PRINT_INPUT_ERROR("Cannot open potential file.");
+  }
+  std::string line;
+  std::getline(input_potential, line);
+  std::vector<std::string> tokens = get_tokens(line);
+  if (tokens[0].size() >= 12) {
+    if (tokens[0].substr(0, 11) == "nep4_charge") {
+      is_nep_charge = true;
+    }
+  } 
+  if (tokens[0].size() >= 16) {
+    if (tokens[0].substr(0, 15) == "nep4_zbl_charge") {
+      is_nep_charge = true;
+    }
+  }
+  input_potential.close();
+
+  return is_nep_charge;
+}
+

src/utilities/read_file.cuh

@@ -17,3 +17,4 @@
 
 int is_valid_int(const char*, int*);
 int is_valid_real(const char*, double*);
+bool check_is_nep_charge();