Merge pull request brucefan1983#1166 from brucefan1983/dump_charge

dump charge

Author: brucefan1983

src/main_gpumd/add_efield.cu

@@ -28,29 +28,6 @@ Add electric field to a group of atoms.
 #include <vector>
 #include <cstring>
 
-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 double* 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;
-  }
-}
-
-// for NEP-charge
 static void __global__ add_efield(
   const int group_size,
   const int group_size_sum,

src/measure/dump_xyz.cu

@@ -18,6 +18,7 @@ Dump per-atom data to user-specified file(s) in the extended XYZ format
 --------------------------------------------------------------------------------------------------*/
 
 #include "dump_xyz.cuh"
+#include "force/force.cuh"
 #include "model/atom.cuh"
 #include "model/box.cuh"
 #include "utilities/common.cuh"
@@ -52,6 +53,8 @@ static __global__ void gpu_sum(const int N, const double* g_data, double* g_data
 
 Dump_XYZ::Dump_XYZ(const char** param, int num_param, const std::vector<Group>& groups, Atom& atom) 
 {
+  is_nep_charge = check_is_nep_charge();
+
   parse(param, num_param, groups);
   if (atom.unwrapped_position.size() < atom.number_of_atoms * 3) {
     atom.unwrapped_position.resize(atom.number_of_atoms * 3);
@@ -139,6 +142,15 @@ void Dump_XYZ::parse(const char** param, int num_param, const std::vector<Group>
       quantities.has_mass_ = true;
       printf("    has mass.\n");
     }
+    if (strcmp(param[m], "charge") == 0) {
+      quantities.has_charge_ = true;
+      if (is_nep_charge){
+        printf("    has charge predicted by NEP-charge.\n");
+      } else {
+        printf("    has charge specified in model.xyz.\n");
+      }
+      
+    }
     if (strcmp(param[m], "virial") == 0) {
       quantities.has_virial_ = true;
       printf("    has virial.\n");
@@ -247,6 +259,9 @@ void Dump_XYZ::output_line2(
   if (quantities.has_mass_) {
     fprintf(fid_, ":mass:R:1");
   }
+  if (quantities.has_charge_) {
+    fprintf(fid_, ":charge:R:1");
+  }
   if (quantities.has_velocity_) {
     fprintf(fid_, ":vel:R:3");
   }
@@ -297,6 +312,14 @@ void Dump_XYZ::process(
   if (quantities.has_mass_) {
     atom.mass.copy_to_host(atom.cpu_mass.data());
   }
+  if (quantities.has_charge_) {
+    if (is_nep_charge) {
+      GPU_Vector<float>& nep_charge = force.potentials[0]->get_charge_reference();
+      nep_charge.copy_to_host(atom.cpu_charge.data());
+    } else {
+      atom.charge.copy_to_host(atom.cpu_charge.data());
+    }
+  }
   if (quantities.has_velocity_) {
     atom.velocity_per_atom.copy_to_host(atom.cpu_velocity_per_atom.data());
   }
@@ -340,6 +363,9 @@ void Dump_XYZ::process(
     if (quantities.has_mass_) {
       fprintf(fid_, " %.8f", atom.cpu_mass[m]);
     }
+    if (quantities.has_charge_) {
+      fprintf(fid_, " %.8f", atom.cpu_charge[m]);
+    }
     if (quantities.has_velocity_) {
       const double natural_to_A_per_fs = 1.0 / TIME_UNIT_CONVERSION;
       for (int d = 0; d < 3; ++d) {

src/measure/dump_xyz.cuh

@@ -64,12 +64,14 @@ public:
     bool has_potential_ = false;
     bool has_unwrapped_position_ = false;
     bool has_mass_ = false;
+    bool has_charge_ = false;
     bool has_virial_ = false;
     bool has_group_ = false;
   };
 
 private:
 
+  bool is_nep_charge = false;
   int grouping_method_ = -1;
   int group_id_ = -1;
   int dump_interval_ = 1;

src/model/atom.cuh

@@ -25,13 +25,13 @@ public:
   std::vector<int> cpu_type;
   std::vector<int> cpu_type_size;
   std::vector<double> cpu_mass;
-  std::vector<double> cpu_charge;
+  std::vector<float> cpu_charge; // float is sufficient for charge
   std::vector<double> cpu_position_per_atom;
   std::vector<double> cpu_velocity_per_atom;
   std::vector<std::string> cpu_atom_symbol;
   GPU_Vector<int> type;                  // per-atom type (1 component)
   GPU_Vector<double> mass;               // per-atom mass (1 component)
-  GPU_Vector<double> charge;             // per-atom charge (1 component)
+  GPU_Vector<float> charge;              // per-atom charge (1 component)
   GPU_Vector<double> position_per_atom;  // per-atom position (3 components)
   GPU_Vector<double> position_temp;      // used to calculated unwrapped_position
   GPU_Vector<double> unwrapped_position; // unwrapped per-atom position (3 components)

src/model/read_xyz.cu

@@ -322,7 +322,7 @@ void read_xyz_in_line_3(
   std::vector<std::string>& cpu_atom_symbol,
   std::vector<int>& cpu_type,
   std::vector<double>& cpu_mass,
-  std::vector<double>& cpu_charge,
+  std::vector<float>& cpu_charge,
   std::vector<double>& cpu_position_per_atom,
   std::vector<double>& cpu_velocity_per_atom,
   std::vector<Group>& group)