Merge branch 'brucefan1983:master' into master

README.md

@@ -66,6 +66,7 @@ There is a standalone C++ implementation of the neuroevolution potential (NEP) i
 | [15]                  | TNEP (tensorial NEP models of dipole and polarizability) |
 | [16]                  | MCMD (hybrid Monte Carlo and molecular dynamics simulations) |
 | [17]                  | PIMD/TRPMD (path-integral molecular dynamics/thermostatted ring-polymer molecular dynamics) |
+| [18]                  | NEMD and NPHug shock methods |
 
 ## References
 
@@ -119,3 +120,6 @@ arXiv:2404.13694 [cond-mat.mtrl-sci]
 
 [17] Penghua Ying, Wenjiang Zhou, Lucas Svensson, Esmée Berger, Erik Fransson, Fredrik Eriksson, Ke Xu, Ting Liang, Jianbin Xu, Bai Song, Shunda Chen, Paul Erhart, Zheyong Fan, [Highly efficient path-integral molecular dynamics simulations with GPUMD using neuroevolution potentials: Case studies on thermal properties of materials](https://arxiv.org/abs/2409.04430),
 arXiv:2409.04430 [cond-mat.mtrl-sci]
+
+[18] Shuning Pan, Jiuyang Shi, Zhixin Liang, Cong Liu, Junjie Wang, Yong Wang, Hui-Tian Wang, Dingyu Xing, and Jian Sun, [Shock compression pathways to pyrite silica from machine learning simulations](https://doi.org/10.1103/PhysRevB.110.224101),
+Phys. Rev. B **110**, 224101 (2024).

doc/bibliography.rst

@@ -20,7 +20,7 @@ Bibliography
    | Erik Bitzek, Pekka Koskinen, Franz Gähler, Michael Moseler, and Peter Gumbsch
    | *Structural Relaxation Made Simple*
    | Phys. Rev. Lett. **97**, 170201 (2006)
-   | DOI: `10.1103/PhysRevLett.97.170201 <https://doi.org/10.1103/PhysRevB.69.144113>`_
+   | DOI: `10.1103/PhysRevLett.97.170201 <https://doi.org/10.1103/PhysRevLett.97.170201>`_
 
 .. [Brorsson2021]
    | Joakim Brorsson, Arsalan Hashemi, Zheyong Fan, Erik Fransson, Fredrik Eriksson, Tapio Ala-Nissila, Arkady V. Krasheninnikov, Hannu-Pekka Komsa, and Paul Erhart

doc/gpumd/input_parameters/compute_adf.rst

@@ -0,0 +1,42 @@
+.. _kw_compute_adf:
+.. index::
+   single: compute_adf (keyword in run.in)
+
+:attr:`compute_adf`
+===================
+
+This keyword computes the angular distribution function (:term:`ADF`) for all atoms or specific triples of species. Each ADF is represented as a histogram, created by measuring the angles formed between a central atom and two neighboring atoms, and binning these angles into `num_bins` bins. Only neighbors with distances `rc_min < R < rc_max` are considered, where `rc_min` and `rc_max` are specified separately for the first and second neighbor atoms in each ADF calculation.
+Currently, this feature is only available for classical :term:`MD`.
+The results are written to the :ref:`adf.out <adf_out>` file.
+
+Syntax
+------
+
+For global ADF, the keyword is used as follows::
+
+  compute_adf <interval> <num_bins> <rc_min> <rc_max>
+
+This means the :term:`ADF` calculations will be performed every :attr:`interval` steps, with :attr:`num_bins` data points.
+
+For local ADF, the keyword is used as follows::
+
+  compute_adf <interval> <num_bins> <itype1> <jtype1> <ktype1> <rc_min_j1> <rc_max_j1> <rc_min_k1> <rc_max_k1> ...
+
+This means the :term:`ADF` calculations will be performed every :attr:`interval` steps, with :attr:`num_bins` data points. 
+The angle formed by the central atom I and neighboring atoms J and K is included in the ADF if the following conditions are met:
+
+- The distance between atoms I and J is between `rc_min_jN` and `rc_max_jN`.
+- The distance between atoms I and K is between `rc_min_kN` and `rc_max_kN`.
+- The type of atom I matches `itypeN`.
+- The type of atom J matches `jtypeN`.
+- The type of atom K matches `ktypeN`.
+- Atoms I, J, and K are distinct.
+
+The ADF value for a bin is computed by dividing the histogram count by the total number of triples that satisfy the criteria, ensuring that the integral of the ADF with respect to the angle is 1. In other words, the ADF is a probability density function.
+
+Example
+-------
+
+   compute_rdf 100 30 0.0 1.2  # Total ADF every 100 MD steps with 30 data points for bond lengths between 0.0 and 1.2
+   compute_rdf 500 50 0 1 1 0.0 1.2 0.0 1.3  # Calculate 0-1-1 ADF every 500 MD steps with 50 data points for I-J bond lengths between 0.0 and 1.2, and I-K bond lengths between 0.0 and 1.3
+   compute_rdf 500 50 0 1 1 0.0 1.2 0.0 1.3 1 0 1 0.0 1.2 0.0 1.3  # Calculate 0-1-1 and 1-0-1 ADF every 500 MD steps with 50 data points

doc/gpumd/input_parameters/index.rst

@@ -45,6 +45,7 @@ Below you can find a listing of keywords for the ``run.in`` input file.
    minimize
    run
    compute
+   compute_adf
    compute_cohesive
    compute_dos
    compute_elastic

doc/gpumd/output_files/adf_out.rst

@@ -0,0 +1,19 @@
+.. _adf_out:
+.. index::
+   single: adf.out (output file)
+
+``adf.out``
+===========
+
+This file contains the Angular Distribution Function (:term:`ADF`). 
+It is generated when the :ref:`compute_adf keyword <kw_compute_adf>` is used.
+
+File Format
+-------------
+For global ADF, the data in this file are organized as follows:
+* Column 1: Angles (in degrees)
+* Column 2: The (:term:`ADF`) for the entire system
+
+For local ADF, the data are organized as follows:
+* Column 1: Angles (in degrees)
+* The next Ntriples columns: The (:term:`ADF`) for each specific atom triple

doc/gpumd/output_files/compute_out.rst

@@ -19,7 +19,7 @@ Assuming that the system is divided into :math:`M` groups according to the group
       
       fx_1 ... fx_M fy_1 ... fy_M fz_1 ... fz_M
 
-  * if virial is computed, there are :math:`3M` columns of group virials (in units of eV) from left to right in a form similar to that for force
+  * if virial is computed, there are :math:`9M` columns of group virials (in units of eV) from left to right in the order of xx, xy, xz, yx, yy, yz, zx, zy, zz.
   * if the potential part of the heat current is computed, there are :math:`3M` columns of group potential heat currents (in units of eV\ :math:`^{3/2}` amu\ :math:`^{-1/2}`) from left to right in a form similar to that for force
   * if the kinetic part of the heat current is computed, there are :math:`3M` columns of group kinetic heat currents (in units of eV\ :math:`^{3/2}` amu\ :math:`^{-1/2}`) from left to right in a form similar to that for force
   * if momentum is computed, there are :math:`3M` columns of group momenta (in units of amu\ :math:`^{1/2}` eV\ :math:`^{1/2}`) from left to right in a form similar to that for force

doc/gpumd/output_files/index.rst

@@ -121,6 +121,10 @@ Output files
      - :ref:`compute_rdf <kw_compute_rdf>`
      - Radial distribution function (:term:`RDF`)
      - Append
+   * - :ref:`adf.out <adf_out>`
+     - :ref:`compute_adf <kw_compute_adf>`
+     - Angular distribution function (:term:`ADF`)
+     - Append
    * - :ref:`mcmd.out <mcmd_out>`
      - :ref:`mc <kw_mc>`
      - Acceptance ratio and species concentrations

src/main_gpumd/run.cu

@@ -449,6 +449,8 @@ void Run::parse_one_keyword(std::vector<std::string>& tokens)
     measure.msd.parse(param, num_param, group);
   } 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_adf") == 0) {
+    measure.adf.parse(param, num_param, box, number_of_types);
   } else if (strcmp(param[0], "compute_hac") == 0) {
     measure.hac.parse(param, num_param);
   } else if (strcmp(param[0], "compute_viscosity") == 0) {

src/main_nep/dataset.cu

@@ -134,13 +134,15 @@ void Dataset::initialize_gpu_data(Parameters& para)
 
   weight_cpu.resize(Nc);
   energy_ref_cpu.resize(Nc);
+  energy_weight_cpu.resize(Nc);
   virial_ref_cpu.resize(Nc * 6);
   force_ref_cpu.resize(N * 3);
   temperature_ref_cpu.resize(N);
 
   for (int n = 0; n < Nc; ++n) {
     weight_cpu[n] = structures[n].weight;
     energy_ref_cpu[n] = structures[n].energy;
+    energy_weight_cpu[n] = structures[n].energy_weight;
     for (int k = 0; k < 6; ++k) {
       virial_ref_cpu[k * Nc + n] = structures[n].virial[k];
     }
@@ -167,11 +169,13 @@ void Dataset::initialize_gpu_data(Parameters& para)
 
   type_weight_gpu.resize(NUM_ELEMENTS);
   energy_ref_gpu.resize(Nc);
+  energy_weight_gpu.resize(Nc);
   virial_ref_gpu.resize(Nc * 6);
   force_ref_gpu.resize(N * 3);
   temperature_ref_gpu.resize(N);
   type_weight_gpu.copy_from_host(para.type_weight_cpu.data());
   energy_ref_gpu.copy_from_host(energy_ref_cpu.data());
+  energy_weight_gpu.copy_from_host(energy_weight_cpu.data());
   virial_ref_gpu.copy_from_host(virial_ref_cpu.data());
   force_ref_gpu.copy_from_host(force_ref_cpu.data());
   temperature_ref_gpu.copy_from_host(temperature_ref_cpu.data());
@@ -436,7 +440,13 @@ std::vector<float> Dataset::get_rmse_force(Parameters& para, const bool use_weig
 }
 
 static __global__ void
-gpu_get_energy_shift(int* g_Na, int* g_Na_sum, float* g_pe, float* g_pe_ref, float* g_energy_shift)
+gpu_get_energy_shift(
+  int* g_Na, 
+  int* g_Na_sum, 
+  float* g_pe, 
+  float* g_pe_ref, 
+  float* g_pe_weight, 
+  float* g_energy_shift)
 {
   int tid = threadIdx.x;
   int bid = blockIdx.x;
@@ -460,12 +470,18 @@ gpu_get_energy_shift(int* g_Na, int* g_Na_sum, float* g_pe, float* g_pe_ref, flo
 
   if (tid == 0) {
     float diff = s_pe[0] / Na - g_pe_ref[bid];
-    g_energy_shift[bid] = diff;
+    g_energy_shift[bid] = diff * g_pe_weight[bid];
   }
 }
 
 static __global__ void gpu_sum_pe_error(
-  float energy_shift, int* g_Na, int* g_Na_sum, float* g_pe, float* g_pe_ref, float* error_gpu)
+  float energy_shift, 
+  int* g_Na, 
+  int* g_Na_sum, 
+  float* g_pe, 
+  float* g_pe_ref, 
+  float* g_pe_weight, 
+  float* error_gpu)
 {
   int tid = threadIdx.x;
   int bid = blockIdx.x;
@@ -489,7 +505,7 @@ static __global__ void gpu_sum_pe_error(
 
   if (tid == 0) {
     float diff = s_pe[0] / Na - g_pe_ref[bid] - energy_shift;
-    error_gpu[bid] = diff * diff;
+    error_gpu[bid] = diff * diff * g_pe_weight[bid];
   }
 }
 
@@ -508,12 +524,21 @@ std::vector<float> Dataset::get_rmse_energy(
 
   if (do_shift) {
     gpu_get_energy_shift<<<Nc, block_size, sizeof(float) * block_size>>>(
-      Na.data(), Na_sum.data(), energy.data(), energy_ref_gpu.data(), error_gpu.data());
+      Na.data(), 
+      Na_sum.data(), 
+      energy.data(), 
+      energy_ref_gpu.data(), 
+      energy_weight_gpu.data(), 
+      error_gpu.data());
     CHECK(gpuMemcpy(error_cpu.data(), error_gpu.data(), mem, gpuMemcpyDeviceToHost));
+    float Nc_with_weight = 0.0f;
     for (int n = 0; n < Nc; ++n) {
+      Nc_with_weight += energy_weight_cpu[n];
       energy_shift_per_structure += error_cpu[n];
     }
-    energy_shift_per_structure /= Nc;
+    if (Nc_with_weight > 0.0f) {
+      energy_shift_per_structure /= Nc_with_weight;
+    }
   }
 
   gpu_sum_pe_error<<<Nc, block_size, sizeof(float) * block_size>>>(
@@ -522,6 +547,7 @@ std::vector<float> Dataset::get_rmse_energy(
     Na_sum.data(),
     energy.data(),
     energy_ref_gpu.data(),
+    energy_weight_gpu.data(), 
     error_gpu.data());
   CHECK(gpuMemcpy(error_cpu.data(), error_gpu.data(), mem, gpuMemcpyDeviceToHost));
 

src/main_nep/dataset.cuh

@@ -46,10 +46,12 @@ public:
   std::vector<float> virial_cpu; // calculated virial in CPU
   std::vector<float> force_cpu;  // calculated force in CPU
 
+  GPU_Vector<float> energy_weight_gpu;    // energy weight in GPU
   GPU_Vector<float> energy_ref_gpu;       // reference energy in GPU
   GPU_Vector<float> virial_ref_gpu;       // reference virial in GPU
   GPU_Vector<float> force_ref_gpu;        // reference force in GPU
   GPU_Vector<float> temperature_ref_gpu;  // reference temperature in GPU
+  std::vector<float> energy_weight_cpu;   // energy weight in CPU
   std::vector<float> energy_ref_cpu;      // reference energy in CPU
   std::vector<float> virial_ref_cpu;      // reference virial in CPU
   std::vector<float> force_ref_cpu;       // reference force in CPU

src/main_nep/fitness.cu

@@ -278,13 +278,7 @@ void Fitness::write_nep_txt(FILE* fid_nep, Parameters& para, float* elite)
       } else {
         fprintf(fid_nep, "nep4 %d ", para.num_types);
       }
-    } else if (para.version == 5) {
-      if (para.enable_zbl) {
-        fprintf(fid_nep, "nep5_zbl %d ", para.num_types);
-      } else {
-        fprintf(fid_nep, "nep5 %d ", para.num_types);
-      }
-    }
+    } 
   } else if (para.train_mode == 1) { // dipole model
     if (para.version == 3) {
       fprintf(fid_nep, "nep3_dipole %d ", para.num_types);

src/main_nep/nep.cu

@@ -329,9 +329,6 @@ void NEP::update_potential(Parameters& para, float* parameters, ANN& ann)
     pointer += ann.num_neurons1;
     ann.w1[t] = pointer;
     pointer += ann.num_neurons1;
-    if (para.version == 5) {
-      pointer += 1; // one extra bias for NEP5 stored in ann.w1[t]
-    }
   }
   ann.b1 = pointer;
   pointer += 1;
@@ -399,29 +396,16 @@ static __global__ void apply_ann(
     // get energy and energy gradient
     float F = 0.0f, Fp[MAX_DIM] = {0.0f};
 
-    if (paramb.version == 5) {
-      apply_ann_one_layer_nep5(
-        annmb.dim,
-        annmb.num_neurons1,
-        annmb.w0[type],
-        annmb.b0[type],
-        annmb.w1[type],
-        annmb.b1,
-        q,
-        F,
-        Fp);
-    } else {
-      apply_ann_one_layer(
-        annmb.dim,
-        annmb.num_neurons1,
-        annmb.w0[type],
-        annmb.b0[type],
-        annmb.w1[type],
-        annmb.b1,
-        q,
-        F,
-        Fp);
-    }
+    apply_ann_one_layer(
+      annmb.dim,
+      annmb.num_neurons1,
+      annmb.w0[type],
+      annmb.b0[type],
+      annmb.w1[type],
+      annmb.b1,
+      q,
+      F,
+      Fp);
     g_pe[n1] = F;
 
     for (int d = 0; d < annmb.dim; ++d) {

src/main_nep/parameters.cu

@@ -158,10 +158,6 @@ void Parameters::read_zbl_in()
 
 void Parameters::calculate_parameters()
 {
-  if (version == 5 && train_mode != 0) {
-    PRINT_INPUT_ERROR("Can only use NEP5 for potential model.");
-  }
-
   if (train_mode != 0 && train_mode != 3) {
     // take virial as dipole or polarizability
     lambda_e = lambda_f = 0.0f;
@@ -193,8 +189,6 @@ void Parameters::calculate_parameters()
     number_of_variables_ann = (dim + 2) * num_neurons1 + 1;
   } else if (version == 4) {
     number_of_variables_ann = (dim + 2) * num_neurons1 * num_types + 1;
-  } else if (version == 5) {
-    number_of_variables_ann = ((dim + 2) * num_neurons1 + 1) * num_types + 1;
   }
 
   number_of_variables_descriptor =
@@ -541,8 +535,8 @@ void Parameters::parse_version(const char** param, int num_param)
   if (!is_valid_int(param[1], &version)) {
     PRINT_INPUT_ERROR("version should be an integer.\n");
   }
-  if (version < 3 || version > 5) {
-    PRINT_INPUT_ERROR("version should = 3 or 4 or 5.");
+  if (version < 3 || version > 4) {
+    PRINT_INPUT_ERROR("version should = 3 or 4.");
   }
 }
 

src/main_nep/snes.cu

@@ -131,13 +131,12 @@ void SNES::find_type_of_variable(Parameters& para)
   // NN part
   if (para.version != 3) {
     int num_ann = (para.train_mode == 2) ? 2 : 1;
-    int num_extra_bias = (para.version == 5) ? 1 : 0;
     for (int ann = 0; ann < num_ann; ++ann) {
       for (int t = 0; t < para.num_types; ++t) {
-        for (int n = 0; n < (para.dim + 2) * para.num_neurons1 + num_extra_bias; ++n) {
+        for (int n = 0; n < (para.dim + 2) * para.num_neurons1; ++n) {
           type_of_variable[n + offset] = t;
         }
-        offset += (para.dim + 2) * para.num_neurons1 + num_extra_bias;
+        offset += (para.dim + 2) * para.num_neurons1;
       }
       ++offset; // the bias
     }

src/main_nep/structure.cu

@@ -155,6 +155,15 @@ static void read_one_structure(
     PRINT_INPUT_ERROR("The second line for each frame should not be empty.");
   }
 
+  // get energy_weight (optional)
+  for (const auto& token : tokens) {
+    const std::string energy_weight_string = "energy_weight=";
+    if (token.substr(0, energy_weight_string.length()) == energy_weight_string) {
+      structure.energy_weight = get_double_from_token(
+        token.substr(energy_weight_string.length(), token.length()), xyz_filename.c_str(), line_number);
+    }
+  }
+
   bool has_energy_in_exyz = false;
   for (const auto& token : tokens) {
     const std::string energy_string = "energy=";
@@ -496,6 +505,7 @@ static void reorder(const int num_batches, std::vector<Structure>& structures)
     structures_copy[nc].weight = structures[nc].weight;
     structures_copy[nc].has_virial = structures[nc].has_virial;
     structures_copy[nc].energy = structures[nc].energy;
+    structures_copy[nc].energy_weight = structures[nc].energy_weight;
     structures_copy[nc].has_temperature = structures[nc].has_temperature;
     structures_copy[nc].temperature = structures[nc].temperature;
     structures_copy[nc].volume = structures[nc].volume;
@@ -534,6 +544,7 @@ static void reorder(const int num_batches, std::vector<Structure>& structures)
     structures[nc].weight = structures_copy[configuration_id[nc]].weight;
     structures[nc].has_virial = structures_copy[configuration_id[nc]].has_virial;
     structures[nc].energy = structures_copy[configuration_id[nc]].energy;
+    structures[nc].energy_weight = structures_copy[configuration_id[nc]].energy_weight;
     structures[nc].has_temperature = structures_copy[configuration_id[nc]].has_temperature;
     structures[nc].temperature = structures_copy[configuration_id[nc]].temperature;
     structures[nc].volume = structures_copy[configuration_id[nc]].volume;

src/main_nep/structure.cuh

@@ -25,6 +25,7 @@ struct Structure {
   int has_temperature;
   float weight;
   float energy = 0.0f;
+  float energy_weight = 1.0f;
   float virial[6];
   float box_original[9];
   float volume;