mcmd

src/main_gpumd/run.cu

@@ -158,6 +158,7 @@ void Run::execute_run_in()
 void Run::perform_a_run()
 {
   integrate.initialize(N, time_step, group, atom);
+  mc.initialize();
   measure.initialize(number_of_steps, time_step, integrate, group, atom, force);
 
 #ifdef USE_PLUMED
@@ -214,6 +215,8 @@ void Run::perform_a_run()
 
     integrate.compute2(time_step, double(step) / number_of_steps, group, box, atom, thermo);
 
+    mc.compute(step, atom, box);
+
     measure.process(
       number_of_steps,
       step,
@@ -252,10 +255,17 @@ void Run::perform_a_run()
   printf("Speed of this run = %g atom*step/second.\n", run_speed);
   print_line_2();
 
-  measure.finalize(integrate, number_of_steps, time_step, integrate.temperature2, box.get_volume(),atom.number_of_beads);
+  measure.finalize(
+    integrate,
+    number_of_steps,
+    time_step,
+    integrate.temperature2,
+    box.get_volume(),
+    atom.number_of_beads);
 
   electron_stop.finalize();
   integrate.finalize();
+  mc.finalize();
   velocity.finalize();
   max_distance_per_step = 0.0;
 }
@@ -366,7 +376,7 @@ void Run::parse_one_keyword(std::vector<std::string>& tokens)
     measure.sdc.parse(param, num_param, group);
   } else if (strcmp(param[0], "compute_msd") == 0) {
     measure.msd.parse(param, num_param, group);
-  }  else if (strcmp(param[0], "compute_rdf") == 0) {
+  } else if (strcmp(param[0], "compute_rdf") == 0) {
     measure.rdf.parse(param, num_param, box, number_of_types, number_of_steps);
   } else if (strcmp(param[0], "compute_hac") == 0) {
     measure.hac.parse(param, num_param);
@@ -392,6 +402,8 @@ void Run::parse_one_keyword(std::vector<std::string>& tokens)
     integrate.parse_move(param, num_param, group);
   } else if (strcmp(param[0], "electron_stop") == 0) {
     electron_stop.parse(param, num_param, atom.number_of_atoms, number_of_types);
+  } else if (strcmp(param[0], "cmc") == 0) {
+    mc.parse_cmc(param, num_param);
   } else if (strcmp(param[0], "run") == 0) {
     parse_run(param, num_param);
   } else {

src/main_gpumd/run.cuh

@@ -22,6 +22,7 @@ class Measure;
 #include "electron_stop.cuh"
 #include "force/force.cuh"
 #include "integrate/integrate.cuh"
+#include "mc/mc.cuh"
 #include "measure/measure.cuh"
 #include "model/atom.cuh"
 #include "model/box.cuh"
@@ -65,6 +66,7 @@ private:
 
   Force force;
   Integrate integrate;
+  MC mc;
   Measure measure;
   Electron_Stop electron_stop;
 };

src/makefile

@@ -32,8 +32,9 @@ LIBS = -lcublas -lcusolver
 SOURCES_GPUMD =                   \
 	$(wildcard main_gpumd/*.cu)   \
 	$(wildcard minimize/*.cu)     \
-	$(wildcard phonon/*.cu)  \
+	$(wildcard phonon/*.cu)       \
 	$(wildcard integrate/*.cu)    \
+	$(wildcard mc/*.cu)           \
 	$(wildcard force/*.cu)        \
 	$(wildcard measure/*.cu)      \
 	$(wildcard model/*.cu)        \
@@ -62,11 +63,12 @@ HEADERS =                         \
 	$(wildcard utilities/*.cuh)   \
 	$(wildcard main_gpumd/*.cuh)  \
 	$(wildcard integrate/*.cuh)   \
+	$(wildcard mc/*.cuh)          \
 	$(wildcard minimize/*.cuh)    \
 	$(wildcard force/*.cuh)       \
 	$(wildcard measure/*.cuh)     \
 	$(wildcard model/*.cuh)       \
-	$(wildcard phonon/*.cuh) \
+	$(wildcard phonon/*.cuh)      \
 	$(wildcard main_nep/*.cuh)
 
 
@@ -86,6 +88,8 @@ nep: $(OBJ_NEP)
 ifdef OS # for Windows
 integrate/%.obj: integrate/%.cu $(HEADERS)
 	$(CC) $(CFLAGS) $(INC) -c $< -o $@
+mc/%.obj: mc/%.cu $(HEADERS)
+	$(CC) $(CFLAGS) $(INC) -c $< -o $@
 minimize/%.obj: minimize/%.cu $(HEADERS)
 	$(CC) $(CFLAGS) $(INC) -c $< -o $@
 force/%.obj: force/%.cu $(HEADERS)
@@ -105,6 +109,8 @@ main_nep/%.obj: main_nep/%.cu $(HEADERS)
 else # for Linux
 integrate/%.o: integrate/%.cu $(HEADERS)
 	$(CC) $(CFLAGS) $(INC) -c $< -o $@
+mc/%.o: mc/%.cu $(HEADERS)
+	$(CC) $(CFLAGS) $(INC) -c $< -o $@
 minimize/%.o: minimize/%.cu $(HEADERS)
 	$(CC) $(CFLAGS) $(INC) -c $< -o $@
 force/%.o: force/%.cu $(HEADERS)

src/mc/mc.cu

@@ -0,0 +1,77 @@
+/*
+    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
+The driver class for the various MC ensembles.
+------------------------------------------------------------------------------*/
+
+#include "mc.cuh"
+#include "mc_ensemble_canonical.cuh"
+#include "model/atom.cuh"
+#include "utilities/common.cuh"
+#include "utilities/read_file.cuh"
+
+void MC::initialize(void)
+{
+  // todo
+}
+
+void MC::finalize(void) { do_mcmd = false; }
+
+void MC::compute(int step, Atom& atom, Box& box)
+{
+  if (do_mcmd) {
+    if ((step + 1) % num_steps_md == 0) {
+      printf("Now do MC after %d MD steps:\n", step + 1);
+      mc_ensemble->compute(atom, box);
+    }
+  }
+}
+
+void MC::parse_cmc(const char** param, int num_param)
+{
+  if (num_param != 4) {
+    PRINT_INPUT_ERROR("cmc should have 3 parameter.\n");
+  }
+
+  if (!is_valid_int(param[1], &num_steps_md)) {
+    PRINT_INPUT_ERROR("number of MD steps for cmc should be an integer.\n");
+  }
+  if (num_steps_md <= 0) {
+    PRINT_INPUT_ERROR("number of MD steps for cmc should be positive.\n");
+  }
+
+  if (!is_valid_int(param[2], &num_steps_mc)) {
+    PRINT_INPUT_ERROR("number of MC steps for cmc should be an integer.\n");
+  }
+  if (num_steps_mc <= 0) {
+    PRINT_INPUT_ERROR("number of MC steps for cmc should be positive.\n");
+  }
+
+  if (!is_valid_real(param[3], &temperature)) {
+    PRINT_INPUT_ERROR("temperature for cmc should be a number.\n");
+  }
+  if (temperature <= 0) {
+    PRINT_INPUT_ERROR("temperature for cmc should be positive.\n");
+  }
+
+  printf("Perform canonical MC:\n");
+  printf("    after every %d MD steps, do %d MC trials.\n", num_steps_md, num_steps_mc);
+  printf("    with a temperature of %g K.\n", temperature);
+
+  mc_ensemble.reset(new MC_Ensemble_Canonical(num_steps_mc, temperature));
+
+  do_mcmd = true;
+}

src/mc/mc.cuh

@@ -0,0 +1,42 @@
+/*
+    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 "mc_ensemble.cuh"
+#include "model/box.cuh"
+#include "model/group.cuh"
+#include <memory>
+#include <vector>
+
+class Atom;
+
+class MC
+{
+public:
+  std::unique_ptr<MC_Ensemble> mc_ensemble;
+
+  void initialize(void);
+  void finalize(void);
+  void compute(int step, Atom& atom, Box& box);
+
+  void parse_cmc(const char** param, int num_param);
+
+private:
+  bool do_mcmd = false;
+  int num_steps_md = 0;
+  int num_steps_mc = 0;
+  double temperature = 0.0;
+};

src/mc/mc_ensemble.cu

@@ -0,0 +1,58 @@
+/*
+    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
+The abstract base class (ABC) for the MC_Ensemble classes.
+------------------------------------------------------------------------------*/
+
+#include "mc_ensemble.cuh"
+#include "utilities/common.cuh"
+#include <chrono>
+
+MC_Ensemble::MC_Ensemble(void)
+{
+  const int n_max = 1000;
+  const int m_max = 1000;
+  NN_radial.resize(n_max);
+  NN_angular.resize(n_max);
+  local_type_before.resize(n_max);
+  local_type_after.resize(n_max);
+  t2_radial_before.resize(n_max * m_max);
+  t2_radial_after.resize(n_max * m_max);
+  t2_angular_before.resize(n_max * m_max);
+  t2_angular_after.resize(n_max * m_max);
+  x12_radial.resize(n_max * m_max);
+  y12_radial.resize(n_max * m_max);
+  z12_radial.resize(n_max * m_max);
+  x12_angular.resize(n_max * m_max);
+  y12_angular.resize(n_max * m_max);
+  z12_angular.resize(n_max * m_max);
+  pe_before.resize(n_max);
+  pe_after.resize(n_max);
+
+  // TODO
+  nep_energy.initialize("../examples/11_NEP_potential_PbTe/nep.txt");
+
+#ifdef DEBUG
+  rng = std::mt19937(13579);
+#else
+  rng = std::mt19937(std::chrono::system_clock::now().time_since_epoch().count());
+#endif
+}
+
+MC_Ensemble::~MC_Ensemble(void)
+{
+  // nothing now
+}

src/mc/mc_ensemble.cuh

@@ -0,0 +1,58 @@
+/*
+    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 "model/atom.cuh"
+#include "model/box.cuh"
+#include "model/group.cuh"
+#include "nep_energy.cuh"
+#include "utilities/gpu_vector.cuh"
+#include <random>
+#include <vector>
+
+class MC_Ensemble
+{
+public:
+  MC_Ensemble(void);
+  virtual ~MC_Ensemble(void);
+
+  virtual void compute(Atom& atom, Box& box) = 0;
+
+protected:
+  int num_steps_mc = 0;
+  double temperature = 0.0;
+  std::mt19937 rng;
+
+  GPU_Vector<int> NN_radial;
+  GPU_Vector<int> NN_angular;
+  GPU_Vector<int> type_before;
+  GPU_Vector<int> type_after;
+  GPU_Vector<int> local_type_before;
+  GPU_Vector<int> local_type_after;
+  GPU_Vector<int> t2_radial_before;
+  GPU_Vector<int> t2_radial_after;
+  GPU_Vector<int> t2_angular_before;
+  GPU_Vector<int> t2_angular_after;
+  GPU_Vector<float> x12_radial;
+  GPU_Vector<float> y12_radial;
+  GPU_Vector<float> z12_radial;
+  GPU_Vector<float> x12_angular;
+  GPU_Vector<float> y12_angular;
+  GPU_Vector<float> z12_angular;
+  GPU_Vector<float> pe_before;
+  GPU_Vector<float> pe_after;
+
+  NEP_Energy nep_energy;
+};