Merge pull request brucefan1983#662 from brucefan1983/add_efield

Add efield

src/main_gpumd/add_efield.cu

@@ -0,0 +1,168 @@
+/*
+    Copyright 2017 Zheyong Fan, Ville Vierimaa, Mikko Ervasti, and Ari Harju
+    This file is part of GPUMD.
+    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
+    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
+    along with GPUMD.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/*----------------------------------------------------------------------------80
+Add electric field to a group of atoms.
+------------------------------------------------------------------------------*/
+
+#include "add_efield.cuh"
+#include "model/atom.cuh"
+#include "model/group.cuh"
+#include "utilities/read_file.cuh"
+#include <iostream>
+#include <vector>
+
+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;
+  }
+}
+
+void Add_Efield::compute(const int step, const std::vector<Group>& groups, Atom& atom)
+{
+  for (int call = 0; call < num_calls_; ++call) {
+    const int step_mod_table_length = step % table_length_[call];
+    const double Ex = efield_table_[call][0 * table_length_[call] + step_mod_table_length];
+    const double Ey = efield_table_[call][1 * table_length_[call] + step_mod_table_length];
+    const double Ez = efield_table_[call][2 * table_length_[call] + step_mod_table_length];
+    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
+    );
+    CUDA_CHECK_KERNEL
+  }
+}
+
+void Add_Efield::parse(const char** param, int num_param, const std::vector<Group>& group)
+{
+  printf("Add electric field.\n");
+
+  // check the number of parameters
+  if (num_param != 6 && num_param != 4) {
+    PRINT_INPUT_ERROR("add_efield should have 5 or 3 parameters.\n");
+  }
+
+  // parse grouping method
+  if (!is_valid_int(param[1], &grouping_method_[num_calls_])) {
+    PRINT_INPUT_ERROR("grouping method should be an integer.\n");
+  }
+  if (grouping_method_[num_calls_] < 0) {
+    PRINT_INPUT_ERROR("grouping method should >= 0.\n");
+  }
+  if (grouping_method_[num_calls_] >= group.size()) {
+    PRINT_INPUT_ERROR("grouping method should < maximum number of grouping methods.\n");
+  }
+
+  // parse group id
+  if (!is_valid_int(param[2], &group_id_[num_calls_])) {
+    PRINT_INPUT_ERROR("group id should be an integer.\n");
+  }
+  if (group_id_[num_calls_] < 0) {
+    PRINT_INPUT_ERROR("group id should >= 0.\n");
+  }
+  if (group_id_[num_calls_] >= group[grouping_method_[num_calls_]].number) {
+    PRINT_INPUT_ERROR("group id should < maximum number of groups in the grouping method.\n");
+  }
+
+  printf(
+    "    for atoms in group %d of grouping method %d.\n", 
+    group_id_[num_calls_], 
+    grouping_method_[num_calls_]
+  );
+
+  if (num_param == 6) {
+    table_length_[num_calls_] = 1;
+    efield_table_[num_calls_].resize(table_length_[num_calls_] * 3);
+    if (!is_valid_real(param[3], &efield_table_[num_calls_][0])) {
+      PRINT_INPUT_ERROR("Ex should be a number.\n");
+    }
+    if (!is_valid_real(param[4], &efield_table_[num_calls_][1])) {
+      PRINT_INPUT_ERROR("Ey should be a number.\n");
+    }
+    if (!is_valid_real(param[5], &efield_table_[num_calls_][2])) {
+      PRINT_INPUT_ERROR("Ez should be a number.\n");
+    }
+    printf("    Ex = %g V/A.\n", efield_table_[num_calls_][0]);
+    printf("    Ey = %g V/A.\n", efield_table_[num_calls_][1]);
+    printf("    Ez = %g V/A.\n", efield_table_[num_calls_][2]);
+  } else {
+    std::ifstream input(param[3]);
+    if (!input.is_open()) {
+      printf("Failed to open %s.\n", param[3]);
+      exit(1);
+    }
+
+    std::vector<std::string> tokens = get_tokens(input);
+    if (tokens.size() != 1) {
+      PRINT_INPUT_ERROR("The first line of the add_efield file should have 1 value.");
+    }
+    table_length_[num_calls_] = get_int_from_token(tokens[0], __FILE__, __LINE__);
+    if (table_length_[num_calls_] < 2) {
+      PRINT_INPUT_ERROR("Number of steps in the add_efield file should >= 2.\n");
+    } else {
+      printf("    number of values in the add_efield file = %d.\n", table_length_[num_calls_]);
+    }
+
+    efield_table_[num_calls_].resize(table_length_[num_calls_] * 3);
+    for (int n = 0; n < table_length_[num_calls_]; ++n) {
+      std::vector<std::string> tokens = get_tokens(input);
+      if (tokens.size() != 3) {
+        PRINT_INPUT_ERROR("Number of electric field components at each step should be 3.");
+      }
+      for (int t = 0; t < 3; ++t) {
+        efield_table_[num_calls_][t * table_length_[num_calls_] + n] = get_double_from_token(tokens[t], __FILE__, __LINE__);
+      }
+    }
+  }
+
+  ++num_calls_;
+
+  if (num_calls_ > 10) {
+    PRINT_INPUT_ERROR("add_efield cannot be used more than 10 times in one run.");
+  }
+}
+
+void Add_Efield::finalize() 
+{ 
+  num_calls_ = 0;
+}

src/main_gpumd/add_efield.cuh

@@ -0,0 +1,38 @@
+/*
+    Copyright 2017 Zheyong Fan, Ville Vierimaa, Mikko Ervasti, and Ari Harju
+    This file is part of GPUMD.
+    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
+    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
+    along with GPUMD.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#pragma once
+
+#include <vector>
+
+class Atom;
+class Group;
+
+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 finalize();
+
+private:
+
+  int num_calls_ = 0;
+  int table_length_[10];
+  std::vector<double> efield_table_[10];
+  int grouping_method_[10];
+  int group_id_[10];
+};

src/main_gpumd/add_force.cu

@@ -24,7 +24,7 @@ Add force to a group of atoms.
 #include <iostream>
 #include <vector>
 
-void __global__
+static void __global__
 add_force(
   const int group_size,
   const int group_size_sum,
@@ -49,9 +49,9 @@ void Add_Force::compute(const int step, const std::vector<Group>& groups, Atom&
 {
   for (int call = 0; call < num_calls_; ++call) {
     const int step_mod_table_length = step % table_length_[call];
-    const float added_fx = force_table_[call][0 * table_length_[call] + step_mod_table_length];
-    const float added_fy = force_table_[call][1 * table_length_[call] + step_mod_table_length];
-    const float added_fz = force_table_[call][2 * table_length_[call] + step_mod_table_length];
+    const double added_fx = force_table_[call][0 * table_length_[call] + step_mod_table_length];
+    const double added_fy = force_table_[call][1 * table_length_[call] + step_mod_table_length];
+    const double added_fz = force_table_[call][2 * table_length_[call] + step_mod_table_length];
     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]];

src/main_gpumd/run.cu

@@ -17,6 +17,7 @@
 Run simulation according to the inputs in the run.in file.
 ------------------------------------------------------------------------------*/
 
+#include "add_efield.cuh"
 #include "add_force.cuh"
 #include "cohesive.cuh"
 #include "electron_stop.cuh"
@@ -259,6 +260,7 @@ void Run::perform_a_run()
 
     electron_stop.compute(time_step, atom);
     add_force.compute(step, group, atom);
+    add_efield.compute(step, group, atom);
 
     integrate.compute2(time_step, double(step) / number_of_steps, group, box, atom, thermo);
 
@@ -314,6 +316,7 @@ void Run::perform_a_run()
 
   electron_stop.finalize();
   add_force.finalize();
+  add_efield.finalize();
   integrate.finalize();
   mc.finalize();
   velocity.finalize();
@@ -466,6 +469,8 @@ void Run::parse_one_keyword(std::vector<std::string>& tokens)
     electron_stop.parse(param, num_param, atom.number_of_atoms, number_of_types);
   } else if (strcmp(param[0], "add_force") == 0) {
     add_force.parse(param, num_param, group);
+  } else if (strcmp(param[0], "add_efield") == 0) {
+    add_efield.parse(param, num_param, group);
   } else if (strcmp(param[0], "mc") == 0) {
     mc.parse_mc(param, num_param, group, atom);
   } else if (strcmp(param[0], "dftd3") == 0) {

src/main_gpumd/run.cuh

@@ -19,6 +19,7 @@ class Force;
 class Integrate;
 class Measure;
 
+#include "add_efield.cuh"
 #include "add_force.cuh"
 #include "electron_stop.cuh"
 #include "force/force.cuh"
@@ -70,4 +71,5 @@ private:
   Measure measure;
   Electron_Stop electron_stop;
   Add_Force add_force;
+  Add_Efield add_efield;
 };

src/model/atom.cuh

@@ -25,11 +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<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<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)