From a5a2b4023ead6c4f9b87dc75897171116a775fe7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 8 Jul 2022 22:17:19 -0400 Subject: [PATCH] simplify and apply clang-format --- src/ML-PACE/compute_pace_extrapolation.cpp | 95 ++- src/ML-PACE/compute_pace_extrapolation.h | 29 +- src/ML-PACE/pair_pace_extrapolation.cpp | 704 ++++++++++----------- src/ML-PACE/pair_pace_extrapolation.h | 77 +-- 4 files changed, 439 insertions(+), 466 deletions(-) diff --git a/src/ML-PACE/compute_pace_extrapolation.cpp b/src/ML-PACE/compute_pace_extrapolation.cpp index 07f128f89b..ebbb7ccd19 100644 --- a/src/ML-PACE/compute_pace_extrapolation.cpp +++ b/src/ML-PACE/compute_pace_extrapolation.cpp @@ -16,8 +16,8 @@ // #include "compute_pace_extrapolation.h" -#include "pair_pace_extrapolation.h" +#include "pair_pace_extrapolation.h" #include "comm.h" #include "error.h" @@ -28,66 +28,63 @@ using namespace LAMMPS_NS; ComputePACEExtrapolation::ComputePACEExtrapolation(class LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) { - if (narg < 3) error->all(FLERR, "Illegal compute pace/extrapolation command"); - peratom_flag = 1; - size_peratom_cols = 0; - scalar_flag = 1; // compute max of gamma + Compute(lmp, narg, arg), pair_pace_extrapolation(nullptr) +{ + peratom_flag = 1; + size_peratom_cols = 0; + scalar_flag = 1; // get next timestep where gamma is available } -ComputePACEExtrapolation::~ComputePACEExtrapolation() { -} - -void ComputePACEExtrapolation::init() { - if (force->pair == nullptr) - error->all(FLERR, "Compute pace/extrapolation requires a pair style pace/extrapolation be defined"); - - if ((modify->get_compute_by_style("pace/extrapolation").size() > 1) && (comm->me == 0)) - error->warning(FLERR, "More than one instance of compute pace/atom"); - +void ComputePACEExtrapolation::init() +{ + if (force->pair) pair_pace_extrapolation = (PairPACEExtrapolation *) force->pair_match("pace/extrapolation", 1); - if (!pair_pace_extrapolation) - error->all(FLERR, "Compute pace/extrapolation requires a `pace/extrapolation` pair style"); + if (!pair_pace_extrapolation) + error->all(FLERR, "Compute pace/extrapolation requires a `pace/extrapolation` pair style"); + + if ((modify->get_compute_by_style("pace/extrapolation").size() > 1) && (comm->me == 0)) + error->warning(FLERR, "More than one instance of compute pace/atom"); } -void ComputePACEExtrapolation::invoke_compute_extrapolation_grades() { - bigint current_timestep = update->ntimestep; - pair_pace_extrapolation->bevaluator_timestep_shift = current_timestep; - int old_vflag_fdotr = pair_pace_extrapolation->vflag_fdotr; - pair_pace_extrapolation->vflag_fdotr = 0; - pair_pace_extrapolation->is_set_energies_forces = false; +void ComputePACEExtrapolation::invoke_compute_extrapolation_grades() +{ + bigint current_timestep = update->ntimestep; + pair_pace_extrapolation->bevaluator_timestep_shift = current_timestep; + int old_vflag_fdotr = pair_pace_extrapolation->vflag_fdotr; + pair_pace_extrapolation->vflag_fdotr = 0; + pair_pace_extrapolation->is_set_energies_forces = false; - pair_pace_extrapolation->compute(0, 0); + pair_pace_extrapolation->compute(0, 0); - pair_pace_extrapolation->is_set_energies_forces = true; - pair_pace_extrapolation->vflag_fdotr = old_vflag_fdotr; + pair_pace_extrapolation->is_set_energies_forces = true; + pair_pace_extrapolation->vflag_fdotr = old_vflag_fdotr; } -double ComputePACEExtrapolation::compute_scalar() { - invoked_scalar = update->ntimestep; +double ComputePACEExtrapolation::compute_scalar() +{ + invoked_scalar = update->ntimestep; - // check the coherence of bevaluator_timestep (when extrapolation grades are computed) and actual timestep - // if not coherent, change pair->bevaluator_timestep_shift to current timestep - // and call invoke_compute_extrapolation_grades without updating energies and forces - if (invoked_scalar != pair_pace_extrapolation->bevaluator_timestep) { - //utils::logmesg(lmp,"[ComputePaceAtom::compute_scalar] Reseting timestep shift to {} (pace timestep={}) and recomputing\n",invoked_scalar,pair->bevaluator_timestep); - invoke_compute_extrapolation_grades(); - } + // check the coherence of bevaluator_timestep (when extrapolation grades are computed) and actual timestep + // if not coherent, change pair->bevaluator_timestep_shift to current timestep + // and call invoke_compute_extrapolation_grades without updating energies and forces + if (invoked_scalar != pair_pace_extrapolation->bevaluator_timestep) { + invoke_compute_extrapolation_grades(); + } - scalar = pair_pace_extrapolation->max_gamma_grade_per_structure; - return scalar; + scalar = pair_pace_extrapolation->max_gamma_grade_per_structure; + return scalar; } -void ComputePACEExtrapolation::compute_peratom() { - invoked_peratom = update->ntimestep; +void ComputePACEExtrapolation::compute_peratom() +{ + invoked_peratom = update->ntimestep; - // check the coherence of bevaluator_timestep (when extrapolation grades are computed) and actual timestep - // if not coherent, change pair->bevaluator_timestep_shift to current timestep - // and call invoke_compute_extrapolation_grades without updating energies and forces - if (invoked_peratom != pair_pace_extrapolation->bevaluator_timestep) { - //utils::logmesg(lmp,"[ComputePaceAtom::compute_peratom] Reseting timestep shift to {} (pace timestep={}) and recomputing\n",invoked_peratom,pair->bevaluator_timestep); - invoke_compute_extrapolation_grades(); - } + // check the coherence of bevaluator_timestep (when extrapolation grades are computed) and actual timestep + // if not coherent, change pair->bevaluator_timestep_shift to current timestep + // and call invoke_compute_extrapolation_grades without updating energies and forces + if (invoked_peratom != pair_pace_extrapolation->bevaluator_timestep) { + invoke_compute_extrapolation_grades(); + } - vector_atom = pair_pace_extrapolation->extrapolation_grade_gamma; -} \ No newline at end of file + vector_atom = pair_pace_extrapolation->extrapolation_grade_gamma; +} diff --git a/src/ML-PACE/compute_pace_extrapolation.h b/src/ML-PACE/compute_pace_extrapolation.h index 1bc7d04189..0329379c45 100644 --- a/src/ML-PACE/compute_pace_extrapolation.h +++ b/src/ML-PACE/compute_pace_extrapolation.h @@ -24,25 +24,24 @@ ComputeStyle(pace/extrapolation,ComputePACEExtrapolation); #ifndef COMPUTE_PACE_H #define COMPUTE_PACE_H - #include "compute.h" #include "pair_pace_extrapolation.h" namespace LAMMPS_NS { - class PairPACEExtrapolation; +class PairPACEExtrapolation; - class ComputePACEExtrapolation : public Compute { - public: - ComputePACEExtrapolation(class LAMMPS *, int, char **); - ~ComputePACEExtrapolation() override; - void init() override; - double compute_scalar() override; - void compute_peratom() override; - private: - PairPACEExtrapolation *pair_pace_extrapolation; - void invoke_compute_extrapolation_grades(); - }; +class ComputePACEExtrapolation : public Compute { + public: + ComputePACEExtrapolation(class LAMMPS *, int, char **); + void init() override; + double compute_scalar() override; + void compute_peratom() override; + + private: + PairPACEExtrapolation *pair_pace_extrapolation; + void invoke_compute_extrapolation_grades(); +}; } // namespace LAMMPS_NS -#endif //COMPUTE_PACE_H -#endif \ No newline at end of file +#endif //COMPUTE_PACE_H +#endif diff --git a/src/ML-PACE/pair_pace_extrapolation.cpp b/src/ML-PACE/pair_pace_extrapolation.cpp index 3cfc182372..624cb69ff5 100644 --- a/src/ML-PACE/pair_pace_extrapolation.cpp +++ b/src/ML-PACE/pair_pace_extrapolation.cpp @@ -17,52 +17,49 @@ Copyright 2022 Yury Lysogorskiy^1, Anton Bochkarev^1, Matous Mrovec^1, Ralf Drau ^1: Ruhr-University Bochum, Bochum, Germany */ - // // Created by Lysogorskiy Yury on 2.01.22. // -#include -#include -#include -#include -#include #include "pair_pace_extrapolation.h" #include "atom.h" +#include "comm.h" #include "dump_custom.h" -#include "neighbor.h" +#include "error.h" +#include "force.h" +#include "math_const.h" +#include "memory.h" +#include "modify.h" #include "neigh_list.h" #include "neigh_request.h" -#include "force.h" -#include "comm.h" -#include "memory.h" -#include "error.h" +#include "neighbor.h" #include "update.h" -#include "modify.h" +#include +#include +#include -#include "math_const.h" - -#include "ace_b_evaluator.h" #include "ace_b_basis.h" +#include "ace_b_evaluator.h" #include "ace_recursive.h" #include "ace_version.h" #include "compute_pace_extrapolation.h" namespace LAMMPS_NS { - struct ACEALImpl { - ACEALImpl() : basis_set(nullptr), ace(nullptr), ctilde_basis_set(nullptr), rec_ace(nullptr) {} +struct ACEALImpl { + ACEALImpl() : basis_set(nullptr), ace(nullptr), ctilde_basis_set(nullptr), rec_ace(nullptr) {} - ~ACEALImpl() { - delete basis_set; - delete ace; - } + ~ACEALImpl() + { + delete basis_set; + delete ace; + } - ACEBBasisSet *basis_set; - ACEBEvaluator *ace; - ACECTildeBasisSet *ctilde_basis_set; - ACERecursiveEvaluator *rec_ace; - }; + ACEBBasisSet *basis_set; + ACEBEvaluator *ace; + ACECTildeBasisSet *ctilde_basis_set; + ACERecursiveEvaluator *rec_ace; +}; } // namespace LAMMPS_NS using namespace LAMMPS_NS; @@ -73,429 +70,416 @@ using namespace MathConst; //added YL - - int elements_num_pace_al = 104; -char const *const elements_pace_al[104] = {"X", "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", - "Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", - "Mn", - "Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", - "Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", - "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", - "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", - "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", - "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr" -}; +char const *const elements_pace_al[104] = { + "X", "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si", + "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", + "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru", + "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", + "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", + "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", + "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr"}; -int AtomicNumberByName_pace_al(char *elname) { - for (int i = 1; i < elements_num_pace_al; i++) - if (strcmp(elname, elements_pace_al[i]) == 0) - return i; - return -1; +int AtomicNumberByName_pace_al(char *elname) +{ + for (int i = 1; i < elements_num_pace_al; i++) + if (strcmp(elname, elements_pace_al[i]) == 0) return i; + return -1; } - /* ---------------------------------------------------------------------- */ -PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp) { - single_enable = 0; - restartinfo = 0; - one_coeff = 1; - manybody_flag = 1; +PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 0; + restartinfo = 0; + one_coeff = 1; + manybody_flag = 1; - natoms = 0; + natoms = 0; - aceimpl = new ACEALImpl; - scale = nullptr; - extrapolation_grade_gamma = nullptr; + aceimpl = new ACEALImpl; + scale = nullptr; + extrapolation_grade_gamma = nullptr; } /* ---------------------------------------------------------------------- check if allocated, since class can be destructed when incomplete ------------------------------------------------------------------------- */ -PairPACEExtrapolation::~PairPACEExtrapolation() { - if (copymode) return; +PairPACEExtrapolation::~PairPACEExtrapolation() +{ + if (copymode) return; - delete aceimpl; + delete aceimpl; - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - memory->destroy(scale); - memory->destroy(map); - memory->destroy(extrapolation_grade_gamma); - } + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(scale); + memory->destroy(map); + memory->destroy(extrapolation_grade_gamma); + } } /* ---------------------------------------------------------------------- */ -void PairPACEExtrapolation::compute(int eflag, int vflag) { - int i, j, ii, jj, inum, jnum; - double delx, dely, delz, evdwl; - double fij[3]; - int *ilist, *jlist, *numneigh, **firstneigh; - double max_gamma_grade = 0; - max_gamma_grade_per_structure = 0; - ev_init(eflag, vflag); +void PairPACEExtrapolation::compute(int eflag, int vflag) +{ + int i, j, ii, jj, inum, jnum; + double delx, dely, delz, evdwl; + double fij[3]; + int *ilist, *jlist, *numneigh, **firstneigh; + double max_gamma_grade = 0; + max_gamma_grade_per_structure = 0; + ev_init(eflag, vflag); - // downwards modified by YL + // downwards modified by YL - double **x = atom->x; - double **f = atom->f; - tagint *tag = atom->tag; - int *type = atom->type; - // number of atoms in cell - int nlocal = atom->nlocal; + double **x = atom->x; + double **f = atom->f; + tagint *tag = atom->tag; + int *type = atom->type; + // number of atoms in cell + int nlocal = atom->nlocal; - int newton_pair = force->newton_pair; + int newton_pair = force->newton_pair; - // number of atoms including ghost atoms - int nall = nlocal + atom->nghost; + // number of atoms including ghost atoms + int nall = nlocal + atom->nghost; - // inum: length of the neighborlists list - inum = list->inum; + // inum: length of the neighborlists list + inum = list->inum; - // ilist: list of "i" atoms for which neighbor lists exist - ilist = list->ilist; + // ilist: list of "i" atoms for which neighbor lists exist + ilist = list->ilist; - //numneigh: the length of each these neigbor list - numneigh = list->numneigh; + //numneigh: the length of each these neigbor list + numneigh = list->numneigh; - // the pointer to the list of neighbors of "i" - firstneigh = list->firstneigh; + // the pointer to the list of neighbors of "i" + firstneigh = list->firstneigh; - if (inum != nlocal) { - error->all(FLERR, "inum: %d nlocal: %d are different", inum, nlocal); + // this happens when used as substyle in pair style hybrid. + // So this check and error effectively disallows use with pair style hybrid. + if (inum != nlocal) { error->all(FLERR, "inum: {} nlocal: {} are different", inum, nlocal); } + + //grow extrapolation_grade_gamma array, that store per-atom extrapolation grades + if (atom->nlocal > natoms) { + memory->destroy(extrapolation_grade_gamma); + natoms = atom->nlocal; + memory->create(extrapolation_grade_gamma, natoms, "pace/atom:gamma"); + } + + //determine the maximum number of neighbours + int max_jnum = 0; + int nei = 0; + for (ii = 0; ii < list->inum; ii++) { + i = ilist[ii]; + jnum = numneigh[i]; + nei = nei + jnum; + if (jnum > max_jnum) max_jnum = jnum; + } + + bigint current_timestep = update->ntimestep; + bool is_bevaluator = (current_timestep - bevaluator_timestep_shift) % gamma_grade_eval_freq == 0; + + if (is_bevaluator) + aceimpl->ace->resize_neighbours_cache(max_jnum); + else + aceimpl->rec_ace->resize_neighbours_cache(max_jnum); + + //loop over atoms + for (ii = 0; ii < list->inum; ii++) { + i = list->ilist[ii]; + const int itype = type[i]; + + const double xtmp = x[i][0]; + const double ytmp = x[i][1]; + const double ztmp = x[i][2]; + + jlist = firstneigh[i]; + jnum = numneigh[i]; + + // checking if neighbours are actually within cutoff range is done inside compute_atom + // mapping from LAMMPS atom types ('type' array) to ACE species is done inside compute_atom + // by using 'ace->element_type_mapping' array + // x: [r0 ,r1, r2, ..., r100] + // i = 0 ,1 + // jnum(0) = 50 + // jlist(neigh ind of 0-atom) = [1,2,10,7,99,25, .. 50 element in total] + try { + if (is_bevaluator) + aceimpl->ace->compute_atom(i, x, type, jnum, jlist); + else + aceimpl->rec_ace->compute_atom(i, x, type, jnum, jlist); + } catch (std::exception &e) { + error->all(FLERR, e.what()); + exit(EXIT_FAILURE); } - - //grow extrapolation_grade_gamma array, that store per-atom extrapolation grades - if (atom->nlocal > natoms) { - memory->destroy(extrapolation_grade_gamma); - natoms = atom->nlocal; - memory->create(extrapolation_grade_gamma, natoms, "pace/atom:gamma"); - } - // Aidan Thompson told RD (26 July 2019) that practically always holds: - // inum = nlocal - // i = ilist(ii) < inum - // j = jlist(jj) < nall - // neighborlist contains neighbor atoms plus skin atoms, - // skin atoms can be removed by setting skin to zero but here - // they are disregarded anyway - - - //determine the maximum number of neighbours - int max_jnum = 0; - int nei = 0; - for (ii = 0; ii < list->inum; ii++) { - i = ilist[ii]; - jnum = numneigh[i]; - nei = nei + jnum; - if (jnum > max_jnum) - max_jnum = jnum; - } - - int current_timestep = update->ntimestep; - bool is_bevaluator = (current_timestep - bevaluator_timestep_shift) % gamma_grade_eval_freq == 0; - - if (is_bevaluator) - aceimpl->ace->resize_neighbours_cache(max_jnum); - else - aceimpl->rec_ace->resize_neighbours_cache(max_jnum); - - //loop over atoms - for (ii = 0; ii < list->inum; ii++) { - i = list->ilist[ii]; - const int itype = type[i]; - - const double xtmp = x[i][0]; - const double ytmp = x[i][1]; - const double ztmp = x[i][2]; - - jlist = firstneigh[i]; - jnum = numneigh[i]; - - // checking if neighbours are actually within cutoff range is done inside compute_atom - // mapping from LAMMPS atom types ('type' array) to ACE species is done inside compute_atom - // by using 'ace->element_type_mapping' array - // x: [r0 ,r1, r2, ..., r100] - // i = 0 ,1 - // jnum(0) = 50 - // jlist(neigh ind of 0-atom) = [1,2,10,7,99,25, .. 50 element in total] - try { - if (is_bevaluator) - aceimpl->ace->compute_atom(i, x, type, jnum, jlist); - else - aceimpl->rec_ace->compute_atom(i, x, type, jnum, jlist); - } catch (std::exception &e) { - error->all(FLERR, e.what()); - exit(EXIT_FAILURE); - } - // 'compute_atom' will update the `ace->e_atom` and `ace->neighbours_forces(jj, alpha)` arrays and max_gamma_grade - - if (is_bevaluator) { - double current_atom_gamma_grade = aceimpl->ace->max_gamma_grade; - if (max_gamma_grade < current_atom_gamma_grade) max_gamma_grade = current_atom_gamma_grade; - bevaluator_timestep = current_timestep; - extrapolation_grade_gamma[i] = current_atom_gamma_grade; - } - - Array2D &neighbours_forces = (is_bevaluator ? aceimpl->ace->neighbours_forces - : aceimpl->rec_ace->neighbours_forces); - //optionally assign global forces arrays - if (is_set_energies_forces) { - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - const int jtype = type[j]; - j &= NEIGHMASK; - delx = x[j][0] - xtmp; - dely = x[j][1] - ytmp; - delz = x[j][2] - ztmp; - - fij[0] = scale[itype][jtype] * neighbours_forces(jj, 0); - fij[1] = scale[itype][jtype] * neighbours_forces(jj, 1); - fij[2] = scale[itype][jtype] * neighbours_forces(jj, 2); - - - f[i][0] += fij[0]; - f[i][1] += fij[1]; - f[i][2] += fij[2]; - f[j][0] -= fij[0]; - f[j][1] -= fij[1]; - f[j][2] -= fij[2]; - - // tally per-atom virial contribution - if (vflag) - ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, - fij[0], fij[1], fij[2], - -delx, -dely, -delz); - } - } - - // tally energy contribution - if (eflag && is_set_energies_forces) { - // evdwl = energy of atom I - DOUBLE_TYPE e_atom; - if (is_bevaluator) - e_atom = aceimpl->ace->e_atom; - else - e_atom = aceimpl->rec_ace->e_atom; - evdwl = scale[itype][itype] * e_atom; - ev_tally_full(i, 2.0 * evdwl, 0.0, 0.0, 0.0, 0.0, 0.0); - } - } - - - if (vflag_fdotr) virial_fdotr_compute(); + // 'compute_atom' will update the `ace->e_atom` and `ace->neighbours_forces(jj, alpha)` arrays and max_gamma_grade if (is_bevaluator) { - //gather together max_gamma_grade_per_structure - MPI_Allreduce(&max_gamma_grade, &max_gamma_grade_per_structure, 1, MPI_DOUBLE, MPI_MAX, world); - - // check if gamma_upper_bound is exceeded - if (max_gamma_grade_per_structure > gamma_upper_bound) { - if (comm->me == 0) - error->all(FLERR, - "Extrapolation grade is too large (gamma={:.3f} > gamma_upper_bound={:.3f}, timestep={}), stopping...\n", - max_gamma_grade_per_structure, gamma_upper_bound, current_timestep); - - MPI_Barrier(world); - MPI_Abort(world, 1); //abort properly with error code '1' if not using many processes - exit(EXIT_FAILURE); - } + double current_atom_gamma_grade = aceimpl->ace->max_gamma_grade; + if (max_gamma_grade < current_atom_gamma_grade) max_gamma_grade = current_atom_gamma_grade; + bevaluator_timestep = current_timestep; + extrapolation_grade_gamma[i] = current_atom_gamma_grade; } -// end modifications YL + Array2D &neighbours_forces = + (is_bevaluator ? aceimpl->ace->neighbours_forces : aceimpl->rec_ace->neighbours_forces); + //optionally assign global forces arrays + if (is_set_energies_forces) { + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + const int jtype = type[j]; + j &= NEIGHMASK; + delx = x[j][0] - xtmp; + dely = x[j][1] - ytmp; + delz = x[j][2] - ztmp; + + fij[0] = scale[itype][jtype] * neighbours_forces(jj, 0); + fij[1] = scale[itype][jtype] * neighbours_forces(jj, 1); + fij[2] = scale[itype][jtype] * neighbours_forces(jj, 2); + + f[i][0] += fij[0]; + f[i][1] += fij[1]; + f[i][2] += fij[2]; + f[j][0] -= fij[0]; + f[j][1] -= fij[1]; + f[j][2] -= fij[2]; + + // tally per-atom virial contribution + if (vflag) + ev_tally_xyz(i, j, nlocal, newton_pair, 0.0, 0.0, fij[0], fij[1], fij[2], -delx, -dely, + -delz); + } + } + + // tally energy contribution + if (eflag && is_set_energies_forces) { + // evdwl = energy of atom I + DOUBLE_TYPE e_atom; + if (is_bevaluator) + e_atom = aceimpl->ace->e_atom; + else + e_atom = aceimpl->rec_ace->e_atom; + evdwl = scale[itype][itype] * e_atom; + ev_tally_full(i, 2.0 * evdwl, 0.0, 0.0, 0.0, 0.0, 0.0); + } + } + + if (vflag_fdotr) virial_fdotr_compute(); + + if (is_bevaluator) { + //gather together max_gamma_grade_per_structure + MPI_Allreduce(&max_gamma_grade, &max_gamma_grade_per_structure, 1, MPI_DOUBLE, MPI_MAX, world); + + // check if gamma_upper_bound is exceeded + if (max_gamma_grade_per_structure > gamma_upper_bound) { + if (comm->me == 0) + error->all(FLERR, + "Extrapolation grade is too large (gamma={:.3f} > gamma_upper_bound={:.3f}, " + "timestep={}), stopping...\n", + max_gamma_grade_per_structure, gamma_upper_bound, current_timestep); + + MPI_Barrier(world); + MPI_Abort(world, 1); //abort properly with error code '1' if not using many processes + exit(EXIT_FAILURE); + } + } + + // end modifications YL } /* ---------------------------------------------------------------------- */ -void PairPACEExtrapolation::allocate() { - allocated = 1; - int n = atom->ntypes; +void PairPACEExtrapolation::allocate() +{ + allocated = 1; + int n = atom->ntypes; - memory->create(setflag, n + 1, n + 1, "pair:setflag"); - memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); - memory->create(map, n + 1, "pair:map"); - memory->create(scale, n + 1, n + 1, "pair:scale"); + memory->create(setflag, n + 1, n + 1, "pair:setflag"); + memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); + memory->create(map, n + 1, "pair:map"); + memory->create(scale, n + 1, n + 1, "pair:scale"); } /* ---------------------------------------------------------------------- global settings ------------------------------------------------------------------------- */ -void PairPACEExtrapolation::settings(int narg, char **arg) { - if (narg > 3) { - error->all(FLERR, - "Illegal pair_style command. Correct form:\n\tpair_style pace/al [gamma_lower_bound] [gamma_upper_bound] [freq]"); - } +void PairPACEExtrapolation::settings(int narg, char **arg) +{ + if (narg > 3) { + error->all(FLERR, + "Illegal pair_style command. Correct form:\n\tpair_style pace/al " + "[gamma_lower_bound] [gamma_upper_bound] [freq]"); + } - if (narg > 0) { - double glb = atof(arg[0]); // gamma lower bound - if (glb < 1.0) - error->all(FLERR, - "Illegal gamma_lower_bound value: it should be real number >= 1.0"); - else - gamma_lower_bound = glb; - } + if (narg > 0) { + double glb = atof(arg[0]); // gamma lower bound + if (glb < 1.0) + error->all(FLERR, "Illegal gamma_lower_bound value: it should be real number >= 1.0"); + else + gamma_lower_bound = glb; + } - if (narg > 1) { - double gub = atof(arg[1]); // gamma upper bound - if (gub < gamma_lower_bound) - error->all(FLERR, - "Illegal gamma_upper_bound value: it should be real number >= gamma_lower_bound >= 1.0"); - else - gamma_upper_bound = gub; - } + if (narg > 1) { + double gub = atof(arg[1]); // gamma upper bound + if (gub < gamma_lower_bound) + error->all( + FLERR, + "Illegal gamma_upper_bound value: it should be real number >= gamma_lower_bound >= 1.0"); + else + gamma_upper_bound = gub; + } - if (narg > 2) { - gamma_grade_eval_freq = atoi(arg[2]); - if (gamma_grade_eval_freq < 1) - error->all(FLERR, - "Illegal gamma_grade_eval_freq value: it should be integer number >= 1"); - } + if (narg > 2) { + gamma_grade_eval_freq = atoi(arg[2]); + if (gamma_grade_eval_freq < 1) + error->all(FLERR, "Illegal gamma_grade_eval_freq value: it should be integer number >= 1"); + } - if (comm->me == 0) { - utils::logmesg(lmp, "ACE/AL version: {}.{}.{}\n", VERSION_YEAR, VERSION_MONTH, VERSION_DAY); - utils::logmesg(lmp, "Extrapolation grade thresholds (lower/upper): {}/{}\n", gamma_lower_bound, - gamma_upper_bound); - utils::logmesg(lmp, "Extrapolation grade evaluation frequency: {}\n", gamma_grade_eval_freq); - } + if (comm->me == 0) { + utils::logmesg(lmp, "ACE/AL version: {}.{}.{}\n", VERSION_YEAR, VERSION_MONTH, VERSION_DAY); + utils::logmesg(lmp, "Extrapolation grade thresholds (lower/upper): {}/{}\n", gamma_lower_bound, + gamma_upper_bound); + utils::logmesg(lmp, "Extrapolation grade evaluation frequency: {}\n", gamma_grade_eval_freq); + } } /* ---------------------------------------------------------------------- set coeffs for one or more type pairs ------------------------------------------------------------------------- */ -void PairPACEExtrapolation::coeff(int narg, char **arg) { +void PairPACEExtrapolation::coeff(int narg, char **arg) +{ - if (narg < 5) - error->all(FLERR, - "Incorrect args for pair coefficients. Correct form:\npair_coeff * * elem1 elem2 ..."); + if (narg < 5) + error->all(FLERR, + "Incorrect args for pair coefficients. Correct form:\npair_coeff * * " + " elem1 elem2 ..."); - if (!allocated) allocate(); + if (!allocated) allocate(); - map_element2type(narg - 4, arg + 4); + map_element2type(narg - 4, arg + 4); - auto potential_file_name = utils::get_potential_file_path(arg[2]); - auto active_set_inv_filename = utils::get_potential_file_path(arg[3]); - char **elemtypes = &arg[4]; + auto potential_file_name = utils::get_potential_file_path(arg[2]); + auto active_set_inv_filename = utils::get_potential_file_path(arg[3]); + char **elemtypes = &arg[4]; + delete aceimpl->basis_set; + delete aceimpl->ctilde_basis_set; - delete aceimpl->basis_set; - delete aceimpl->ctilde_basis_set; + //load potential file + aceimpl->basis_set = new ACEBBasisSet(); + if (comm->me == 0) utils::logmesg(lmp, "Loading {}\n", potential_file_name); + aceimpl->basis_set->load(potential_file_name); - //load potential file - aceimpl->basis_set = new ACEBBasisSet(); - if (comm->me == 0) utils::logmesg(lmp, "Loading {}\n", potential_file_name); - aceimpl->basis_set->load(potential_file_name); + //convert the basis set to CTilde format + aceimpl->ctilde_basis_set = new ACECTildeBasisSet(); + *aceimpl->ctilde_basis_set = aceimpl->basis_set->to_ACECTildeBasisSet(); - //convert the basis set to CTilde format - aceimpl->ctilde_basis_set = new ACECTildeBasisSet(); - *aceimpl->ctilde_basis_set = aceimpl->basis_set->to_ACECTildeBasisSet(); + if (comm->me == 0) { + utils::logmesg(lmp, "Total number of basis functions\n"); - if (comm->me == 0) { - utils::logmesg(lmp, "Total number of basis functions\n"); - - for (SPECIES_TYPE mu = 0; mu < aceimpl->basis_set->nelements; mu++) { - int n_r1 = aceimpl->basis_set->total_basis_size_rank1[mu]; - int n = aceimpl->basis_set->total_basis_size[mu]; - utils::logmesg(lmp, "\t{}: {} (r=1) {} (r>1)\n", aceimpl->basis_set->elements_name[mu], n_r1, - n); - } + for (SPECIES_TYPE mu = 0; mu < aceimpl->basis_set->nelements; mu++) { + int n_r1 = aceimpl->basis_set->total_basis_size_rank1[mu]; + int n = aceimpl->basis_set->total_basis_size[mu]; + utils::logmesg(lmp, "\t{}: {} (r=1) {} (r>1)\n", aceimpl->basis_set->elements_name[mu], n_r1, + n); } + } + // read args that map atom types to PACE elements + // map[i] = which element the Ith atom type is, -1 if not mapped + // map[0] is not used + delete aceimpl->ace; + delete aceimpl->rec_ace; - // read args that map atom types to PACE elements - // map[i] = which element the Ith atom type is, -1 if not mapped - // map[0] is not used - delete aceimpl->ace; - delete aceimpl->rec_ace; + aceimpl->ace = new ACEBEvaluator(); + aceimpl->ace->element_type_mapping.init(atom->ntypes + 1); - aceimpl->ace = new ACEBEvaluator(); - aceimpl->ace->element_type_mapping.init(atom->ntypes + 1); + aceimpl->rec_ace = new ACERecursiveEvaluator(); + aceimpl->rec_ace->set_recursive(true); + aceimpl->rec_ace->element_type_mapping.init(atom->ntypes + 1); + aceimpl->rec_ace->element_type_mapping.fill(-1); //-1 means atom not included into potential - aceimpl->rec_ace = new ACERecursiveEvaluator(); - aceimpl->rec_ace->set_recursive(true); - aceimpl->rec_ace->element_type_mapping.init(atom->ntypes + 1); - aceimpl->rec_ace->element_type_mapping.fill(-1); //-1 means atom not included into potential + FILE *species_type_file = nullptr; - FILE *species_type_file = nullptr; - - const int n = atom->ntypes; - element_names.resize(n); - for (int i = 1; i <= n; i++) { - char *elemname = elemtypes[i - 1]; - element_names[i - 1] = elemname; - if (strcmp(elemname, "NULL") == 0) { - // species_type=-1 value will not reach ACE Evaluator::compute_atom, - // but if it will ,then error will be thrown there - aceimpl->ace->element_type_mapping(i) = -1; - map[i] = -1; - if (comm->me == 0) utils::logmesg(lmp, "Skipping LAMMPS atom type #{}(NULL)\n", i); - } else { - // dump species types for reconstruction of atomic configurations - int atomic_number = AtomicNumberByName_pace_al(elemname); - if (atomic_number == -1) error->all(FLERR, "'{}' is not a valid element\n", elemname); - SPECIES_TYPE mu = aceimpl->basis_set->get_species_index_by_name(elemname); - if (mu != -1) { - if (comm->me == 0) - utils::logmesg(lmp, "Mapping LAMMPS atom type #{}({}) -> ACE species type #{}\n", i, - elemname, mu); - map[i] = mu; - // set up LAMMPS atom type to ACE species mapping for ace evaluators - aceimpl->ace->element_type_mapping(i) = mu; - aceimpl->rec_ace->element_type_mapping(i) = mu; - } else { - error->all(FLERR, "Element {} is not supported by ACE-potential from file {}", elemname, - potential_file_name); - } - } + const int n = atom->ntypes; + element_names.resize(n); + for (int i = 1; i <= n; i++) { + char *elemname = elemtypes[i - 1]; + element_names[i - 1] = elemname; + if (strcmp(elemname, "NULL") == 0) { + // species_type=-1 value will not reach ACE Evaluator::compute_atom, + // but if it will ,then error will be thrown there + aceimpl->ace->element_type_mapping(i) = -1; + map[i] = -1; + if (comm->me == 0) utils::logmesg(lmp, "Skipping LAMMPS atom type #{}(NULL)\n", i); + } else { + // dump species types for reconstruction of atomic configurations + int atomic_number = AtomicNumberByName_pace_al(elemname); + if (atomic_number == -1) error->all(FLERR, "'{}' is not a valid element\n", elemname); + SPECIES_TYPE mu = aceimpl->basis_set->get_species_index_by_name(elemname); + if (mu != -1) { + if (comm->me == 0) + utils::logmesg(lmp, "Mapping LAMMPS atom type #{}({}) -> ACE species type #{}\n", i, + elemname, mu); + map[i] = mu; + // set up LAMMPS atom type to ACE species mapping for ace evaluators + aceimpl->ace->element_type_mapping(i) = mu; + aceimpl->rec_ace->element_type_mapping(i) = mu; + } else { + error->all(FLERR, "Element {} is not supported by ACE-potential from file {}", elemname, + potential_file_name); + } } + } - aceimpl->ace->set_basis(*aceimpl->basis_set); - aceimpl->rec_ace->set_basis(*aceimpl->ctilde_basis_set); + aceimpl->ace->set_basis(*aceimpl->basis_set); + aceimpl->rec_ace->set_basis(*aceimpl->ctilde_basis_set); - if (comm->me == 0) utils::logmesg(lmp, "Loading ASI {}\n", active_set_inv_filename); - aceimpl->ace->load_active_set(active_set_inv_filename); - - // clear setflag since coeff() called once with I,J = * * - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - scale[i][j] = 1.0; + if (comm->me == 0) utils::logmesg(lmp, "Loading ASI {}\n", active_set_inv_filename); + aceimpl->ace->load_active_set(active_set_inv_filename); + // clear setflag since coeff() called once with I,J = * * + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) scale[i][j] = 1.0; } /* ---------------------------------------------------------------------- init specific to this pair style ------------------------------------------------------------------------- */ -void PairPACEExtrapolation::init_style() { - if (atom->tag_enable == 0) - error->all(FLERR, "Pair style PACE requires atom IDs"); - if (force->newton_pair == 0) - error->all(FLERR, "Pair style PACE requires newton pair on"); +void PairPACEExtrapolation::init_style() +{ + if (atom->tag_enable == 0) error->all(FLERR, "Pair style PACE requires atom IDs"); + if (force->newton_pair == 0) error->all(FLERR, "Pair style PACE requires newton pair on"); - // request a full neighbor list - neighbor->add_request(this, NeighConst::REQ_FULL); + // request a full neighbor list + neighbor->add_request(this, NeighConst::REQ_FULL); } /* ---------------------------------------------------------------------- init for one type pair i,j and corresponding j,i ------------------------------------------------------------------------- */ -double PairPACEExtrapolation::init_one(int i, int j) { - if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set"); - //cutoff from the basis set's radial functions settings - scale[j][i] = scale[i][j]; - return aceimpl->basis_set->radial_functions->cut(map[i], map[j]); +double PairPACEExtrapolation::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set"); + //cutoff from the basis set's radial functions settings + scale[j][i] = scale[i][j]; + return aceimpl->basis_set->radial_functions->cut(map[i], map[j]); } /* ---------------------------------------------------------------------- extract method for extracting value of scale variable ---------------------------------------------------------------------- */ -void *PairPACEExtrapolation::extract(const char *str, int &dim) { - dim = 2; - if (strcmp(str, "scale") == 0) return (void *) scale; - return nullptr; +void *PairPACEExtrapolation::extract(const char *str, int &dim) +{ + dim = 2; + if (strcmp(str, "scale") == 0) return (void *) scale; + return nullptr; } - diff --git a/src/ML-PACE/pair_pace_extrapolation.h b/src/ML-PACE/pair_pace_extrapolation.h index 568b4ce9b1..2ce0966c96 100644 --- a/src/ML-PACE/pair_pace_extrapolation.h +++ b/src/ML-PACE/pair_pace_extrapolation.h @@ -18,10 +18,9 @@ Copyright 2022 Yury Lysogorskiy^1, Anton Bochkarev^1, Matous Mrovec^1, Ralf Drau // Created by Lysogorskiy Yury on 1.01.22. // - #ifdef PAIR_CLASS // clang-format off -PairStyle(pace/extrapolation,PairPACEExtrapolation) +PairStyle(pace/extrapolation,PairPACEExtrapolation); // clang-format on #else @@ -33,53 +32,47 @@ PairStyle(pace/extrapolation,PairPACEExtrapolation) namespace LAMMPS_NS { - //forward declaration - class ComputePACEExtrapolation; - class DumpPACEExtrapolation; +//forward declaration +class ComputePACEExtrapolation; +class DumpPACEExtrapolation; - class PairPACEExtrapolation : public Pair { - friend class ComputePACEExtrapolation; - friend class DumpPACEExtrapolation; +class PairPACEExtrapolation : public Pair { + friend class ComputePACEExtrapolation; + friend class DumpPACEExtrapolation; - public: - PairPACEExtrapolation(class LAMMPS *); + public: + PairPACEExtrapolation(class LAMMPS *); + ~PairPACEExtrapolation() override; - ~PairPACEExtrapolation() override; + void compute(int, int) override; + void settings(int, char **) override; + void coeff(int, char **) override; + void init_style() override; + double init_one(int, int) override; + void *extract(const char *, int &) override; - void compute(int, int) override; + protected: + struct ACEALImpl *aceimpl; + bigint gamma_grade_eval_freq = 1; + bool is_set_energies_forces = true; // if set, then update forces and energies + int natoms; //total number of atoms - void settings(int, char **) override; + double gamma_lower_bound = 1.5; + double gamma_upper_bound = 10; + double max_gamma_grade_per_structure = 0; - void coeff(int, char **) override; + void allocate(); + std::vector element_names; // list of elements (used by dump pace/extrapolation) + double rcutmax; // max cutoff for all elements + int nelements; // # of unique elements + bigint bevaluator_timestep; // timestep, on which gamma grade were computed + bigint bevaluator_timestep_shift = 0; // + double *extrapolation_grade_gamma; //per-atom gamma value - void init_style() override; + double **scale; +}; - double init_one(int, int) override; - - void *extract(const char *, int &) override; - - protected: - struct ACEALImpl *aceimpl; - int gamma_grade_eval_freq = 1; - bool is_set_energies_forces = true; // if set, then update forces and energies - int natoms; //total number of atoms - - double gamma_lower_bound = 1.5; - double gamma_upper_bound = 10; - double max_gamma_grade_per_structure = 0; - - void allocate(); - std::vector element_names; // list of elements (used by dump pace/extrapolation) - double rcutmax; // max cutoff for all elements - int nelements; // # of unique elements - int bevaluator_timestep; // timestep, on which gamma grade were computed - int bevaluator_timestep_shift = 0; // - double *extrapolation_grade_gamma; //per-atom gamma value - - double **scale; - }; - -} +} // namespace LAMMPS_NS #endif -#endif \ No newline at end of file +#endif