|
| 1 | +/** |
| 2 | + * @file esolver_nep.cpp |
| 3 | +#include "source_io/module_parameter/parameter.h" |
| 4 | + * @brief Implementation of ESolver_NEP class for neuroevolution potential (NEP). |
| 5 | + * |
| 6 | + * This file contains the implementation of the ESolver_NEP class, which is used for solving the energy and forces in a |
| 7 | + * NEP simulation. |
| 8 | + * NEP is a method for training deep neural networks to accurately predict the potential energy surface of a |
| 9 | + * molecular system. |
| 10 | + * |
| 11 | + * For more information about NEP, see the following reference: |
| 12 | + * 1. https://gpumd.org/potentials/nep.html |
| 13 | + * 2. https://doi.org/10.1002/mgea.70028 |
| 14 | + * |
| 15 | + * @author MoseyQAQ |
| 16 | + * @date 2025-10-10 |
| 17 | + */ |
| 18 | +#include "esolver_nep.h" |
| 19 | + |
| 20 | +#include "source_base/parallel_common.h" |
| 21 | +#include "source_base/timer.h" |
| 22 | +#include "source_io/output_log.h" |
| 23 | +#include "source_io/cif_io.h" |
| 24 | + |
| 25 | +#include <numeric> |
| 26 | +#include <unordered_map> |
| 27 | + |
| 28 | +using namespace ModuleESolver; |
| 29 | + |
| 30 | +void ESolver_NEP::before_all_runners(UnitCell& ucell, const Input_para& inp) |
| 31 | +{ |
| 32 | + nep_potential = 0.0; |
| 33 | + nep_force.create(ucell.nat, 3); |
| 34 | + nep_virial.create(3, 3); |
| 35 | + atype.resize(ucell.nat); |
| 36 | + _e.resize(ucell.nat); |
| 37 | + _f.resize(3 * ucell.nat); |
| 38 | + _v.resize(9 * ucell.nat); |
| 39 | + |
| 40 | + ModuleIO::CifParser::write(PARAM.globalv.global_out_dir + "STRU.cif", |
| 41 | + ucell, |
| 42 | + "# Generated by ABACUS ModuleIO::CifParser", |
| 43 | + "data_?"); |
| 44 | + |
| 45 | +#ifdef __NEP |
| 46 | + /// determine the type map from STRU to NEP model |
| 47 | + type_map(ucell); |
| 48 | +#endif |
| 49 | +} |
| 50 | + |
| 51 | +void ESolver_NEP::runner(UnitCell& ucell, const int istep) |
| 52 | +{ |
| 53 | + ModuleBase::TITLE("ESolver_NEP", "runner"); |
| 54 | + ModuleBase::timer::tick("ESolver_NEP", "runner"); |
| 55 | + |
| 56 | + // note that NEP are column major, thus a transpose is needed |
| 57 | + // cell |
| 58 | + std::vector<double> cell(9, 0.0); |
| 59 | + cell[0] = ucell.latvec.e11 * ucell.lat0_angstrom; |
| 60 | + cell[1] = ucell.latvec.e21 * ucell.lat0_angstrom; |
| 61 | + cell[2] = ucell.latvec.e31 * ucell.lat0_angstrom; |
| 62 | + cell[3] = ucell.latvec.e12 * ucell.lat0_angstrom; |
| 63 | + cell[4] = ucell.latvec.e22 * ucell.lat0_angstrom; |
| 64 | + cell[5] = ucell.latvec.e32 * ucell.lat0_angstrom; |
| 65 | + cell[6] = ucell.latvec.e13 * ucell.lat0_angstrom; |
| 66 | + cell[7] = ucell.latvec.e23 * ucell.lat0_angstrom; |
| 67 | + cell[8] = ucell.latvec.e33 * ucell.lat0_angstrom; |
| 68 | + |
| 69 | + // coord |
| 70 | + std::vector<double> coord(3 * ucell.nat, 0.0); |
| 71 | + int iat = 0; |
| 72 | + const int nat = ucell.nat; |
| 73 | + for (int it = 0; it < ucell.ntype; ++it) |
| 74 | + { |
| 75 | + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) |
| 76 | + { |
| 77 | + coord[iat] = ucell.atoms[it].tau[ia].x * ucell.lat0_angstrom; |
| 78 | + coord[iat + nat] = ucell.atoms[it].tau[ia].y * ucell.lat0_angstrom; |
| 79 | + coord[iat + 2 * nat] = ucell.atoms[it].tau[ia].z * ucell.lat0_angstrom; |
| 80 | + iat++; |
| 81 | + } |
| 82 | + } |
| 83 | + assert(ucell.nat == iat); |
| 84 | + |
| 85 | +#ifdef __NEP |
| 86 | + nep_potential = 0.0; |
| 87 | + nep_force.zero_out(); |
| 88 | + nep_virial.zero_out(); |
| 89 | + |
| 90 | + nep.compute(atype, cell, coord, _e, _f, _v); |
| 91 | + |
| 92 | + // unit conversion |
| 93 | + const double fact_e = 1.0 / ModuleBase::Ry_to_eV; |
| 94 | + const double fact_f = 1.0 / (ModuleBase::Ry_to_eV * ModuleBase::ANGSTROM_AU); |
| 95 | + const double fact_v = 1.0 / (ucell.omega * ModuleBase::Ry_to_eV); |
| 96 | + |
| 97 | + |
| 98 | + // potential energy |
| 99 | + nep_potential = fact_e * std::accumulate(_e.begin(), _e.end(), 0.0) ; |
| 100 | + GlobalV::ofs_running << " #TOTAL ENERGY# " << std::setprecision(11) << nep_potential * ModuleBase::Ry_to_eV << " eV" |
| 101 | + << std::endl; |
| 102 | + |
| 103 | + // forces |
| 104 | + for (int i = 0; i < nat; ++i) |
| 105 | + { |
| 106 | + nep_force(i, 0) = _f[i] * fact_f; |
| 107 | + nep_force(i, 1) = _f[i + nat] * fact_f; |
| 108 | + nep_force(i, 2) = _f[i + 2 * nat] * fact_f; |
| 109 | + } |
| 110 | + |
| 111 | + // virial |
| 112 | + std::vector<double> v_sum(9, 0.0); |
| 113 | + for (int j = 0; j < 9; ++j) |
| 114 | + { |
| 115 | + for (int i = 0; i < nat; ++i) |
| 116 | + { |
| 117 | + int index = j * nat + i; |
| 118 | + v_sum[j] += _v[index]; |
| 119 | + } |
| 120 | + } |
| 121 | + |
| 122 | + // virial -> stress |
| 123 | + for (int i = 0; i < 3; ++i) |
| 124 | + { |
| 125 | + for (int j = 0; j < 3; ++j) |
| 126 | + { |
| 127 | + nep_virial(i, j) = v_sum[3 * i + j] * fact_v; |
| 128 | + } |
| 129 | + } |
| 130 | +#else |
| 131 | + ModuleBase::WARNING_QUIT("ESolver_NEP", "Please recompile with -D__NEP"); |
| 132 | +#endif |
| 133 | + ModuleBase::timer::tick("ESolver_NEP", "runner"); |
| 134 | +} |
| 135 | + |
| 136 | +double ESolver_NEP::cal_energy() |
| 137 | +{ |
| 138 | + return nep_potential; |
| 139 | +} |
| 140 | + |
| 141 | +void ESolver_NEP::cal_force(UnitCell& ucell, ModuleBase::matrix& force) |
| 142 | +{ |
| 143 | + force = nep_force; |
| 144 | + ModuleIO::print_force(GlobalV::ofs_running, ucell, "TOTAL-FORCE (eV/Angstrom)", force, false); |
| 145 | +} |
| 146 | + |
| 147 | +void ESolver_NEP::cal_stress(UnitCell& ucell, ModuleBase::matrix& stress) |
| 148 | +{ |
| 149 | + stress = nep_virial; |
| 150 | + ModuleIO::print_stress("TOTAL-STRESS", stress, true, false, GlobalV::ofs_running); |
| 151 | + |
| 152 | + // external stress |
| 153 | + double unit_transform = ModuleBase::RYDBERG_SI / pow(ModuleBase::BOHR_RADIUS_SI, 3) * 1.0e-8; |
| 154 | + double external_stress[3] = {PARAM.inp.press1, PARAM.inp.press2, PARAM.inp.press3}; |
| 155 | + for (int i = 0; i < 3; i++) |
| 156 | + { |
| 157 | + stress(i, i) -= external_stress[i] / unit_transform; |
| 158 | + } |
| 159 | +} |
| 160 | + |
| 161 | +void ESolver_NEP::after_all_runners(UnitCell& ucell) |
| 162 | +{ |
| 163 | + GlobalV::ofs_running << "\n --------------------------------------------" << std::endl; |
| 164 | + GlobalV::ofs_running << std::setprecision(16); |
| 165 | + GlobalV::ofs_running << " !FINAL_ETOT_IS " << nep_potential * ModuleBase::Ry_to_eV << " eV" << std::endl; |
| 166 | + GlobalV::ofs_running << " --------------------------------------------\n\n" << std::endl; |
| 167 | +} |
| 168 | + |
| 169 | +#ifdef __NEP |
| 170 | +void ESolver_NEP::type_map(const UnitCell& ucell) |
| 171 | +{ |
| 172 | + // parse the element list from NEP model file |
| 173 | + std::unordered_map<std::string, int> label; |
| 174 | + std::string temp; |
| 175 | + for (int i = 0; i < nep.element_list.size(); ++i) |
| 176 | + { |
| 177 | + label[nep.element_list[i]] = i; //> label: map from element string to index int. |
| 178 | + } |
| 179 | + |
| 180 | + std::cout << "\n Element list of model file " << nep_file << " " << std::endl; |
| 181 | + std::cout << " ----------------------------------------------------------------"; |
| 182 | + int count = 0; |
| 183 | + for (auto it = label.begin(); it != label.end(); ++it) |
| 184 | + { |
| 185 | + if (count % 5 == 0) |
| 186 | + { |
| 187 | + std::cout << std::endl; |
| 188 | + std::cout << " "; |
| 189 | + } |
| 190 | + count++; |
| 191 | + temp = it->first + ": " + std::to_string(it->second); |
| 192 | + std::cout << std::left << std::setw(10) << temp; |
| 193 | + } |
| 194 | + std::cout << "\n -----------------------------------------------------------------" << std::endl; |
| 195 | + |
| 196 | + // parse the atype based on the element list |
| 197 | + int iat = 0; |
| 198 | + for (int it = 0; it < ucell.ntype; ++it) |
| 199 | + { |
| 200 | + for (int ia = 0; ia < ucell.atoms[it].na; ++ia) |
| 201 | + { |
| 202 | + if (label.find(ucell.atoms[it].label) == label.end()) |
| 203 | + { |
| 204 | + ModuleBase::WARNING_QUIT("ESolver_NEP", |
| 205 | + "The label " + ucell.atoms[it].label + " is not found in the type map."); |
| 206 | + } |
| 207 | + atype[iat] = label[ucell.atoms[it].label]; |
| 208 | + iat++; |
| 209 | + } |
| 210 | + } |
| 211 | + assert(ucell.nat == iat); |
| 212 | +} |
| 213 | +#endif |
0 commit comments