Merge branch 'brucefan1983:master' into multi-layers

potentials/nep/readme.md

@@ -5,7 +5,7 @@
 
 | NEP model file   | Paper for the training data      | Paper for the NEP model   | Comments  |
 | --------- | -------------------------------------- | -------------------------------------- | ------------------- |
-| Si_2022_NEP2.txt   | Albert P. Bartók et al., [Machine Learning a General-Purpose Interatomic Potential for Silicon](https://doi.org/10.1103/PhysRevX.8.041048), Phys. Rev. X **8**, 041048 (2018).   | Zheyong Fan et al., [Neuroevolution machine learning potentials: Combining high accuracy and low cost in atomistic simulations and application to heat transport](https://doi.org/10.1103/PhysRevB.104.104309), Phys. Rev. B. **104**, 104309 (2021). | General-purpose potential for silicon |
+| Si_2022_NEP3_3body.txt; Si_2022_NEP3_4body.txt; Si_2022_NEP3_5body.txt   | Albert P. Bartók et al., [Machine Learning a General-Purpose Interatomic Potential for Silicon](https://doi.org/10.1103/PhysRevX.8.041048), Phys. Rev. X **8**, 041048 (2018).   | Zheyong Fan et al., [GPUMD: A package for constructing accurate machine-learned potentials and performing highly efficient atomistic simulations](https://aip.scitation.org/doi/10.1063/5.0106617), The Journal of Chemical Physics **157**, 114801 (2022). | General-purpose potential for silicon |
 | C_2022_NEP3.txt  | Volker L. Deringer et al., [Machine learning based interatomic potential for amorphous carbon](https://doi.org/10.1103/PhysRevB.95.094203), Phys. Rev. B **95**, 094203 (2017).  | Zheyong Fan et al., [GPUMD: A package for constructing accurate machine-learned potentials and performing highly efficient atomistic simulations](https://aip.scitation.org/doi/10.1063/5.0106617), The Journal of Chemical Physics **157**, 114801 (2022).  | Mainly for amorphous carbon |
 | C_2024_NEP4.txt  | Patrick Rowe et al., [An accurate and transferable machine learning potential for carbon](https://doi.org/10.1063/5.0005084), J. Chem. Phys. **153**, 034702 (2020);  | Zheyong Fan et al., [Combining linear-scaling quantum transport and machine-learning molecular dynamics to study thermal and electronic transports in complex materials](https://doi.org/10.1088/1361-648X/ad31c2). Journal of Physics: Condensed Matter **36**, 245901 (2024)  | General-purpose potential for carbon |
 | Song-2024-UNEP-v1-AgAlAuCrCuMgMoNiPbPdPtTaTiVWZr.txt  | Keke Song et al., [General-purpose machine-learned potential for 16 elemental metals and their alloys](https://doi.org/10.48550/arXiv.2311.04732), arXiv:2311.04732 [cond-mat.mtrl-sci]  | Keke Song et al., [General-purpose machine-learned potential for 16 elemental metals and their alloys](https://doi.org/10.48550/arXiv.2311.04732), arXiv:2311.04732 [cond-mat.mtrl-sci]  | General-purpose potential for 16 metals and their arbitrary alloys |

src/main_nep/structure.cu

@@ -87,7 +87,9 @@ static void read_force(
   const int force_offset,
   std::ifstream& input,
   const Parameters& para,
-  Structure& structure)
+  Structure& structure,
+  std::string& xyz_filename,
+  int& line_number)
 {
   structure.type.resize(structure.num_atom);
   structure.x.resize(structure.num_atom);
@@ -99,17 +101,19 @@ static void read_force(
 
   for (int na = 0; na < structure.num_atom; ++na) {
     std::vector<std::string> tokens = get_tokens(input);
+    line_number++;
+
     if (tokens.size() != num_columns) {
       PRINT_INPUT_ERROR("Number of items for an atom line mismatches properties.");
     }
     std::string atom_symbol(tokens[0 + species_offset]);
-    structure.x[na] = get_float_from_token(tokens[0 + pos_offset], __FILE__, __LINE__);
-    structure.y[na] = get_float_from_token(tokens[1 + pos_offset], __FILE__, __LINE__);
-    structure.z[na] = get_float_from_token(tokens[2 + pos_offset], __FILE__, __LINE__);
+    structure.x[na] = get_float_from_token(tokens[0 + pos_offset], xyz_filename.c_str(), line_number);
+    structure.y[na] = get_float_from_token(tokens[1 + pos_offset], xyz_filename.c_str(), line_number);
+    structure.z[na] = get_float_from_token(tokens[2 + pos_offset], xyz_filename.c_str(), line_number);
     if (num_columns > 4) {
-      structure.fx[na] = get_float_from_token(tokens[0 + force_offset], __FILE__, __LINE__);
-      structure.fy[na] = get_float_from_token(tokens[1 + force_offset], __FILE__, __LINE__);
-      structure.fz[na] = get_float_from_token(tokens[2 + force_offset], __FILE__, __LINE__);
+      structure.fx[na] = get_float_from_token(tokens[0 + force_offset], xyz_filename.c_str(), line_number);
+      structure.fy[na] = get_float_from_token(tokens[1 + force_offset], xyz_filename.c_str(), line_number);
+      structure.fz[na] = get_float_from_token(tokens[2 + force_offset], xyz_filename.c_str(), line_number);
     }
 
     bool is_allowed_element = false;
@@ -125,9 +129,16 @@ static void read_force(
   }
 }
 
-static void read_one_structure(const Parameters& para, std::ifstream& input, Structure& structure)
+static void read_one_structure(
+  const Parameters& para, 
+  std::ifstream& input, 
+  Structure& structure,
+  std::string& xyz_filename,
+  int& line_number)
 {
   std::vector<std::string> tokens = get_tokens_without_unwanted_spaces(input);
+  line_number++;
+
   for (auto& token : tokens) {
     std::transform(
       token.begin(), token.end(), token.begin(), [](unsigned char c) { return std::tolower(c); });
@@ -143,7 +154,7 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
     if (token.substr(0, energy_string.length()) == energy_string) {
       has_energy_in_exyz = true;
       structure.energy = get_float_from_token(
-        token.substr(energy_string.length(), token.length()), __FILE__, __LINE__);
+        token.substr(energy_string.length(), token.length()), xyz_filename.c_str(), line_number);
       structure.energy /= structure.num_atom;
     }
   }
@@ -157,7 +168,7 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
     if (token.substr(0, temperature_string.length()) == temperature_string) {
       structure.has_temperature = true;
       structure.temperature = get_float_from_token(
-        token.substr(temperature_string.length(), token.length()), __FILE__, __LINE__);
+        token.substr(temperature_string.length(), token.length()), xyz_filename.c_str(), line_number);
     }
   }
   if (para.train_mode == 3 && !structure.has_temperature) {
@@ -172,7 +183,7 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
     const std::string weight_string = "weight=";
     if (token.substr(0, weight_string.length()) == weight_string) {
       structure.weight = get_float_from_token(
-        token.substr(weight_string.length(), token.length()), __FILE__, __LINE__);
+        token.substr(weight_string.length(), token.length()), xyz_filename.c_str(), line_number);
       if (structure.weight <= 0.0f || structure.weight > 100.0f) {
         PRINT_INPUT_ERROR("Configuration weight should > 0 and <= 100.");
       }
@@ -190,8 +201,7 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
           tokens[n + m].substr(
             (m == 0) ? (lattice_string.length() + 1) : 0,
             (m == 8) ? (tokens[n + m].length() - 1) : tokens[n + m].length()),
-          __FILE__,
-          __LINE__);
+          xyz_filename.c_str(), line_number);
       }
       change_box(para, structure);
     }
@@ -211,8 +221,7 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
           tokens[n + m].substr(
             (m == 0) ? (virial_string.length() + 1) : 0,
             (m == 8) ? (tokens[n + m].length() - 1) : tokens[n + m].length()),
-          __FILE__,
-          __LINE__);
+          xyz_filename.c_str(), line_number);
         structure.virial[reduced_index[m]] /= structure.num_atom;
       }
     }
@@ -231,8 +240,7 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
           tokens[n + m].substr(
             (m == 0) ? (stress_string.length() + 1) : 0,
             (m == 8) ? (tokens[n + m].length() - 1) : tokens[n + m].length()),
-          __FILE__,
-          __LINE__);
+          xyz_filename.c_str(), line_number);
         virials_from_stress[reduced_index[m]] *= -volume / structure.num_atom;
       }
     }
@@ -279,8 +287,7 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
             tokens[n + m].substr(
               (m == 0) ? (dipole_string.length() + 1) : 0,
               (m == 2) ? (tokens[n + m].length() - 1) : tokens[n + m].length()),
-            __FILE__,
-            __LINE__);
+            xyz_filename.c_str(), line_number);
           structure.virial[m] /= structure.num_atom;
         }
       }
@@ -309,8 +316,7 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
             tokens[n + m].substr(
               (m == 0) ? (pol_string.length() + 1) : 0,
               (m == 8) ? (tokens[n + m].length() - 1) : tokens[n + m].length()),
-            __FILE__,
-            __LINE__);
+            xyz_filename.c_str(), line_number);
           structure.virial[reduced_index[m]] /= structure.num_atom;
         }
       }
@@ -365,39 +371,46 @@ static void read_one_structure(const Parameters& para, std::ifstream& input, Str
       }
       for (int k = 0; k < sub_tokens.size() / 3; ++k) {
         if (k < species_position) {
-          species_offset += get_int_from_token(sub_tokens[k * 3 + 2], __FILE__, __LINE__);
+          species_offset += get_int_from_token(sub_tokens[k * 3 + 2], xyz_filename.c_str(), line_number);
         }
         if (k < pos_position) {
-          pos_offset += get_int_from_token(sub_tokens[k * 3 + 2], __FILE__, __LINE__);
+          pos_offset += get_int_from_token(sub_tokens[k * 3 + 2], xyz_filename.c_str(), line_number);
         }
         if (k < force_position) {
-          force_offset += get_int_from_token(sub_tokens[k * 3 + 2], __FILE__, __LINE__);
+          force_offset += get_int_from_token(sub_tokens[k * 3 + 2], xyz_filename.c_str(), line_number);
         }
-        num_columns += get_int_from_token(sub_tokens[k * 3 + 2], __FILE__, __LINE__);
+        num_columns += get_int_from_token(sub_tokens[k * 3 + 2], xyz_filename.c_str(), line_number);
       }
     }
   }
 
-  read_force(num_columns, species_offset, pos_offset, force_offset, input, para, structure);
+  read_force(num_columns, species_offset, pos_offset, force_offset, input, para, structure, xyz_filename, line_number);
 }
 
 static void
-read_exyz(const Parameters& para, std::ifstream& input, std::vector<Structure>& structures)
+read_exyz(
+  const Parameters& para, 
+  std::ifstream& input, 
+  std::vector<Structure>& structures, 
+  std::string& xyz_filename)
 {
+  int line_number = 0;
   int Nc = 0;
   while (true) {
     std::vector<std::string> tokens = get_tokens(input);
+    line_number++;
+
     if (tokens.size() == 0) {
       break;
     } else if (tokens.size() > 1) {
       PRINT_INPUT_ERROR("The first line for each frame should have one value.");
     }
     Structure structure;
-    structure.num_atom = get_int_from_token(tokens[0], __FILE__, __LINE__);
+    structure.num_atom = get_int_from_token(tokens[0], xyz_filename.c_str(), line_number);
     if (structure.num_atom < 1) {
       PRINT_INPUT_ERROR("Number of atoms for each frame should >= 1.");
     }
-    read_one_structure(para, input, structure);
+    read_one_structure(para, input, structure, xyz_filename, line_number);
     structures.emplace_back(structure);
     ++Nc;
   }
@@ -545,9 +558,10 @@ bool read_structures(bool is_train, Parameters& para, std::vector<Structure>& st
     }
   } else {
     print_line_1();
-    is_train ? printf("Started reading train.xyz.\n") : printf("Started reading test.xyz.\n");
+    std::string xyz_filename = is_train ? "train.xyz" : "test.xyz";
+    std::cout << "Started reading " << xyz_filename << std::endl;
     print_line_2();
-    read_exyz(para, input, structures);
+    read_exyz(para, input, structures, xyz_filename);
     input.close();
   }
 

tools/shift_energy_to_zero/shift_energy_to_zero.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 # coding: utf-8
 '''
-Goal:
+Purpose:
     Shift total energy to ~0 by substracting atomic energies using the least square method.
 Run example:
     $ python shift_energy_to_zero.py your.xyz
@@ -34,17 +34,56 @@ def SVD_A(A,b):
 all_elements = sorted(set(flatten_comprehension([i.get_chemical_symbols() for i in all_frames])))
 coeff_matrix = np.zeros((len(all_frames),len(all_elements)))
 energy_matrix = np.zeros((len(all_frames),1))
-for i in range(len(all_frames)):
-    for j in range(len(all_elements)):
-        coeff_matrix[i][j]= all_frames[i].get_chemical_symbols().count(all_elements[j])
-    energy_matrix[i][0] = all_frames[i].get_potential_energy()
+
+#test whether get_potential_energy works
+default_energy_store = True
+try:
+    all_frames[0].get_potential_energy()
+except RuntimeError:
+    tmp = all_frames[0].info
+    tmp =list(tmp.keys())
+    tmp1 = [i.upper() for i in tmp]
+    if "ENERGY" not in tmp1:
+        raise SystemError("Ensure your xyz file has energy.")
+    keyword = tmp[tmp1.index("ENERGY")]
+    default_energy_store= False
+
+if default_energy_store:
+    for i in range(len(all_frames)):
+        for j in range(len(all_elements)):
+            coeff_matrix[i][j]= all_frames[i].get_chemical_symbols().count(all_elements[j])
+        energy_matrix[i][0] = all_frames[i].get_potential_energy()
+else:
+    for i in range(len(all_frames)):
+        for j in range(len(all_elements)):
+            coeff_matrix[i][j]= all_frames[i].get_chemical_symbols().count(all_elements[j])
+        energy_matrix[i][0] = all_frames[i].info[keyword]    
+print('Normalizing energy....')
+#check if matrix is underdetermined
+
+if np.linalg.matrix_rank(coeff_matrix) < len(all_elements):
+    print("Warning! The coeff_matrix is underdetermined, adding constrains....")
+    number_of_constrains = len(all_elements) - np.linalg.matrix_rank(coeff_matrix)
+    import itertools
+    to_add_constrain_pairs = []
+    for i in itertools.combinations(range(len(all_elements)), 2):
+        additional_matrix = np.zeros(len(all_elements))
+        additional_matrix[i[0]] = 1
+        additional_matrix[i[1]] = -1
+        additional_energy = np.zeros(1)
+        coeff_matrix = np.r_[coeff_matrix,[additional_matrix]]
+        energy_matrix = np.r_[energy_matrix,[additional_energy]]
+
 print('Atomic energy---->')
 atomic_shifted_energy = SVD_A(coeff_matrix,energy_matrix)
 for i in range(len(all_elements)):
     print("%s:%f" %(all_elements[i],atomic_shifted_energy[i][0]))
-print('Normalized energy---->')
+    
 shifted_energy = (energy_matrix-np.matmul(coeff_matrix,atomic_shifted_energy)).flatten()
-print(shifted_energy)
+
+print("Averaged energies now: %f eV. " %shifted_energy[:len(all_frames)].mean())
+print("Absolute maximum energy now: %f eV." %max(abs(shifted_energy[:len(all_frames)])))
+
 for i in range(len(all_frames)):
     try:
         forces = all_frames[i].get_forces()
@@ -54,4 +93,4 @@ def SVD_A(A,b):
     all_frames[i].calc = None
     all_frames[i].info['energy'] = shifted_energy[i]
 write("shifted.xyz",all_frames)
-print("Done!")
+print("Done! New shifted.xyz is generated.")