ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

simplify and apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-07-08 22:17:19 -04:00
parent dcf9082591
commit a5a2b4023e
4 changed files with 439 additions and 466 deletions
Download Patch File Download Diff File Expand all files Collapse all files

95
src/ML-PACE/compute_pace_extrapolation.cpp
View File

@ -16,8 +16,8 @@
// //
#include "compute_pace_extrapolation.h" #include "compute_pace_extrapolation.h"
#include "pair_pace_extrapolation.h"
#include "pair_pace_extrapolation.h"
#include "comm.h" #include "comm.h"
#include "error.h" #include "error.h"
@ -28,66 +28,63 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
ComputePACEExtrapolation::ComputePACEExtrapolation(class LAMMPS *lmp, int narg, char **arg) : ComputePACEExtrapolation::ComputePACEExtrapolation(class LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg) { Compute(lmp, narg, arg), pair_pace_extrapolation(nullptr)
if (narg < 3) error->all(FLERR, "Illegal compute pace/extrapolation command"); {
peratom_flag = 1; peratom_flag = 1;
size_peratom_cols = 0; size_peratom_cols = 0;
scalar_flag = 1; // compute max of gamma scalar_flag = 1; // get next timestep where gamma is available
} }
ComputePACEExtrapolation::~ComputePACEExtrapolation() { void ComputePACEExtrapolation::init()
} {
if (force->pair)
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");
pair_pace_extrapolation = (PairPACEExtrapolation *) force->pair_match("pace/extrapolation", 1); pair_pace_extrapolation = (PairPACEExtrapolation *) force->pair_match("pace/extrapolation", 1);
if (!pair_pace_extrapolation) if (!pair_pace_extrapolation)
error->all(FLERR, "Compute pace/extrapolation requires a `pace/extrapolation` pair style"); 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() { void ComputePACEExtrapolation::invoke_compute_extrapolation_grades()
bigint current_timestep = update->ntimestep; {
pair_pace_extrapolation->bevaluator_timestep_shift = current_timestep; bigint current_timestep = update->ntimestep;
int old_vflag_fdotr = pair_pace_extrapolation->vflag_fdotr; pair_pace_extrapolation->bevaluator_timestep_shift = current_timestep;
pair_pace_extrapolation->vflag_fdotr = 0; int old_vflag_fdotr = pair_pace_extrapolation->vflag_fdotr;
pair_pace_extrapolation->is_set_energies_forces = false; 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->is_set_energies_forces = true;
pair_pace_extrapolation->vflag_fdotr = old_vflag_fdotr; pair_pace_extrapolation->vflag_fdotr = old_vflag_fdotr;
} }
double ComputePACEExtrapolation::compute_scalar() { double ComputePACEExtrapolation::compute_scalar()
invoked_scalar = update->ntimestep; {
invoked_scalar = update->ntimestep;
// check the coherence of bevaluator_timestep (when extrapolation grades are computed) and actual timestep // 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 // if not coherent, change pair->bevaluator_timestep_shift to current timestep
// and call invoke_compute_extrapolation_grades without updating energies and forces // and call invoke_compute_extrapolation_grades without updating energies and forces
if (invoked_scalar != pair_pace_extrapolation->bevaluator_timestep) { 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();
invoke_compute_extrapolation_grades(); }
}
scalar = pair_pace_extrapolation->max_gamma_grade_per_structure; scalar = pair_pace_extrapolation->max_gamma_grade_per_structure;
return scalar; return scalar;
} }
void ComputePACEExtrapolation::compute_peratom() { void ComputePACEExtrapolation::compute_peratom()
invoked_peratom = update->ntimestep; {
invoked_peratom = update->ntimestep;
// check the coherence of bevaluator_timestep (when extrapolation grades are computed) and actual timestep // 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 // if not coherent, change pair->bevaluator_timestep_shift to current timestep
// and call invoke_compute_extrapolation_grades without updating energies and forces // and call invoke_compute_extrapolation_grades without updating energies and forces
if (invoked_peratom != pair_pace_extrapolation->bevaluator_timestep) { 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();
invoke_compute_extrapolation_grades(); }
}
vector_atom = pair_pace_extrapolation->extrapolation_grade_gamma; vector_atom = pair_pace_extrapolation->extrapolation_grade_gamma;
} }

29
src/ML-PACE/compute_pace_extrapolation.h
View File

@ -24,25 +24,24 @@ ComputeStyle(pace/extrapolation,ComputePACEExtrapolation);
#ifndef COMPUTE_PACE_H #ifndef COMPUTE_PACE_H
#define COMPUTE_PACE_H #define COMPUTE_PACE_H
#include "compute.h" #include "compute.h"
#include "pair_pace_extrapolation.h" #include "pair_pace_extrapolation.h"
namespace LAMMPS_NS { namespace LAMMPS_NS {
class PairPACEExtrapolation; class PairPACEExtrapolation;
class ComputePACEExtrapolation : public Compute { class ComputePACEExtrapolation : public Compute {
public: public:
ComputePACEExtrapolation(class LAMMPS *, int, char **); ComputePACEExtrapolation(class LAMMPS *, int, char **);
~ComputePACEExtrapolation() override; void init() override;
void init() override; double compute_scalar() override;
double compute_scalar() override; void compute_peratom() override;
void compute_peratom() override;
private: private:
PairPACEExtrapolation *pair_pace_extrapolation; PairPACEExtrapolation *pair_pace_extrapolation;
void invoke_compute_extrapolation_grades(); void invoke_compute_extrapolation_grades();
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS
#endif //COMPUTE_PACE_H #endif //COMPUTE_PACE_H
#endif #endif

704
src/ML-PACE/pair_pace_extrapolation.cpp
View File

@ -17,52 +17,49 @@ Copyright 2022 Yury Lysogorskiy^1, Anton Bochkarev^1, Matous Mrovec^1, Ralf Drau
^1: Ruhr-University Bochum, Bochum, Germany ^1: Ruhr-University Bochum, Bochum, Germany
*/ */
// //
// Created by Lysogorskiy Yury on 2.01.22. // Created by Lysogorskiy Yury on 2.01.22.
// //
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <mpi.h>
#include "pair_pace_extrapolation.h" #include "pair_pace_extrapolation.h"
#include "atom.h" #include "atom.h"
#include "comm.h"
#include "dump_custom.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_list.h"
#include "neigh_request.h" #include "neigh_request.h"
#include "force.h" #include "neighbor.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "update.h" #include "update.h"
#include "modify.h"
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "math_const.h"
#include "ace_b_evaluator.h"
#include "ace_b_basis.h" #include "ace_b_basis.h"
#include "ace_b_evaluator.h"
#include "ace_recursive.h" #include "ace_recursive.h"
#include "ace_version.h" #include "ace_version.h"
#include "compute_pace_extrapolation.h" #include "compute_pace_extrapolation.h"
namespace LAMMPS_NS { namespace LAMMPS_NS {
struct ACEALImpl { struct ACEALImpl {
ACEALImpl() : basis_set(nullptr), ace(nullptr), ctilde_basis_set(nullptr), rec_ace(nullptr) {} ACEALImpl() : basis_set(nullptr), ace(nullptr), ctilde_basis_set(nullptr), rec_ace(nullptr) {}
~ACEALImpl() { ~ACEALImpl()
delete basis_set; {
delete ace; delete basis_set;
} delete ace;
}
ACEBBasisSet *basis_set; ACEBBasisSet *basis_set;
ACEBEvaluator *ace; ACEBEvaluator *ace;
ACECTildeBasisSet *ctilde_basis_set; ACECTildeBasisSet *ctilde_basis_set;
ACERecursiveEvaluator *rec_ace; ACERecursiveEvaluator *rec_ace;
}; };
} // namespace LAMMPS_NS } // namespace LAMMPS_NS
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -73,429 +70,416 @@ using namespace MathConst;
//added YL //added YL
int elements_num_pace_al = 104; 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", char const *const elements_pace_al[104] = {
"Mg", "Al", "Si", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", "X", "H", "He", "Li", "Be", "B", "C", "N", "O", "F", "Ne", "Na", "Mg", "Al", "Si",
"Mn", "P", "S", "Cl", "Ar", "K", "Ca", "Sc", "Ti", "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu",
"Fe", "Co", "Ni", "Cu", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Zn", "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr", "Y", "Zr", "Nb", "Mo", "Tc", "Ru",
"Y", "Zr", "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Rh", "Pd", "Ag", "Cd", "In", "Sn", "Sb", "Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr",
"Te", "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd", "Pm", "Sm", "Eu", "Gd", "Nd", "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W",
"Tb", "Dy", "Ho", "Er", "Tm", "Yb", "Lu", "Hf", "Ta", "W", "Re", "Os", "Ir", "Re", "Os", "Ir", "Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac",
"Pt", "Au", "Hg", "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th", "Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr"};
"Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr"
};
int AtomicNumberByName_pace_al(char *elname) { int AtomicNumberByName_pace_al(char *elname)
for (int i = 1; i < elements_num_pace_al; i++) {
if (strcmp(elname, elements_pace_al[i]) == 0) for (int i = 1; i < elements_num_pace_al; i++)
return i; if (strcmp(elname, elements_pace_al[i]) == 0) return i;
return -1; return -1;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp) { PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp)
single_enable = 0; {
restartinfo = 0; single_enable = 0;
one_coeff = 1; restartinfo = 0;
manybody_flag = 1; one_coeff = 1;
manybody_flag = 1;
natoms = 0; natoms = 0;
aceimpl = new ACEALImpl; aceimpl = new ACEALImpl;
scale = nullptr; scale = nullptr;
extrapolation_grade_gamma = nullptr; extrapolation_grade_gamma = nullptr;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
PairPACEExtrapolation::~PairPACEExtrapolation() { PairPACEExtrapolation::~PairPACEExtrapolation()
if (copymode) return; {
if (copymode) return;
delete aceimpl; delete aceimpl;
if (allocated) { if (allocated) {
memory->destroy(setflag); memory->destroy(setflag);
memory->destroy(cutsq); memory->destroy(cutsq);
memory->destroy(scale); memory->destroy(scale);
memory->destroy(map); memory->destroy(map);
memory->destroy(extrapolation_grade_gamma); memory->destroy(extrapolation_grade_gamma);
} }
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
void PairPACEExtrapolation::compute(int eflag, int vflag) { void PairPACEExtrapolation::compute(int eflag, int vflag)
int i, j, ii, jj, inum, jnum; {
double delx, dely, delz, evdwl; int i, j, ii, jj, inum, jnum;
double fij[3]; double delx, dely, delz, evdwl;
int *ilist, *jlist, *numneigh, **firstneigh; double fij[3];
double max_gamma_grade = 0; int *ilist, *jlist, *numneigh, **firstneigh;
max_gamma_grade_per_structure = 0; double max_gamma_grade = 0;
ev_init(eflag, vflag); max_gamma_grade_per_structure = 0;
ev_init(eflag, vflag);
// downwards modified by YL // downwards modified by YL
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
tagint *tag = atom->tag; tagint *tag = atom->tag;
int *type = atom->type; int *type = atom->type;
// number of atoms in cell // number of atoms in cell
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
int newton_pair = force->newton_pair; int newton_pair = force->newton_pair;
// number of atoms including ghost atoms // number of atoms including ghost atoms
int nall = nlocal + atom->nghost; int nall = nlocal + atom->nghost;
// inum: length of the neighborlists list // inum: length of the neighborlists list
inum = list->inum; inum = list->inum;
// ilist: list of "i" atoms for which neighbor lists exist // ilist: list of "i" atoms for which neighbor lists exist
ilist = list->ilist; ilist = list->ilist;
//numneigh: the length of each these neigbor list //numneigh: the length of each these neigbor list
numneigh = list->numneigh; numneigh = list->numneigh;
// the pointer to the list of neighbors of "i" // the pointer to the list of neighbors of "i"
firstneigh = list->firstneigh; firstneigh = list->firstneigh;
if (inum != nlocal) { // this happens when used as substyle in pair style hybrid.
error->all(FLERR, "inum: %d nlocal: %d are different", inum, nlocal); // 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);
} }
// 'compute_atom' will update the `ace->e_atom` and `ace->neighbours_forces(jj, alpha)` arrays and max_gamma_grade
//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<DOUBLE_TYPE> &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) { if (is_bevaluator) {
//gather together max_gamma_grade_per_structure double current_atom_gamma_grade = aceimpl->ace->max_gamma_grade;
MPI_Allreduce(&max_gamma_grade, &max_gamma_grade_per_structure, 1, MPI_DOUBLE, MPI_MAX, world); if (max_gamma_grade < current_atom_gamma_grade) max_gamma_grade = current_atom_gamma_grade;
bevaluator_timestep = current_timestep;
// check if gamma_upper_bound is exceeded extrapolation_grade_gamma[i] = current_atom_gamma_grade;
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 Array2D<DOUBLE_TYPE> &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() { void PairPACEExtrapolation::allocate()
allocated = 1; {
int n = atom->ntypes; allocated = 1;
int n = atom->ntypes;
memory->create(setflag, n + 1, n + 1, "pair:setflag"); memory->create(setflag, n + 1, n + 1, "pair:setflag");
memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); memory->create(cutsq, n + 1, n + 1, "pair:cutsq");
memory->create(map, n + 1, "pair:map"); memory->create(map, n + 1, "pair:map");
memory->create(scale, n + 1, n + 1, "pair:scale"); memory->create(scale, n + 1, n + 1, "pair:scale");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
global settings global settings
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairPACEExtrapolation::settings(int narg, char **arg) { void PairPACEExtrapolation::settings(int narg, char **arg)
if (narg > 3) { {
error->all(FLERR, if (narg > 3) {
"Illegal pair_style command. Correct form:\n\tpair_style pace/al [gamma_lower_bound] [gamma_upper_bound] [freq]"); error->all(FLERR,
} "Illegal pair_style command. Correct form:\n\tpair_style pace/al "
"[gamma_lower_bound] [gamma_upper_bound] [freq]");
}
if (narg > 0) { if (narg > 0) {
double glb = atof(arg[0]); // gamma lower bound double glb = atof(arg[0]); // gamma lower bound
if (glb < 1.0) if (glb < 1.0)
error->all(FLERR, error->all(FLERR, "Illegal gamma_lower_bound value: it should be real number >= 1.0");
"Illegal gamma_lower_bound value: it should be real number >= 1.0"); else
else gamma_lower_bound = glb;
gamma_lower_bound = glb; }
}
if (narg > 1) { if (narg > 1) {
double gub = atof(arg[1]); // gamma upper bound double gub = atof(arg[1]); // gamma upper bound
if (gub < gamma_lower_bound) if (gub < gamma_lower_bound)
error->all(FLERR, error->all(
"Illegal gamma_upper_bound value: it should be real number >= gamma_lower_bound >= 1.0"); FLERR,
else "Illegal gamma_upper_bound value: it should be real number >= gamma_lower_bound >= 1.0");
gamma_upper_bound = gub; else
} gamma_upper_bound = gub;
}
if (narg > 2) { if (narg > 2) {
gamma_grade_eval_freq = atoi(arg[2]); gamma_grade_eval_freq = atoi(arg[2]);
if (gamma_grade_eval_freq < 1) if (gamma_grade_eval_freq < 1)
error->all(FLERR, error->all(FLERR, "Illegal gamma_grade_eval_freq value: it should be integer number >= 1");
"Illegal gamma_grade_eval_freq value: it should be integer number >= 1"); }
}
if (comm->me == 0) { if (comm->me == 0) {
utils::logmesg(lmp, "ACE/AL version: {}.{}.{}\n", VERSION_YEAR, VERSION_MONTH, VERSION_DAY); 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, utils::logmesg(lmp, "Extrapolation grade thresholds (lower/upper): {}/{}\n", gamma_lower_bound,
gamma_upper_bound); gamma_upper_bound);
utils::logmesg(lmp, "Extrapolation grade evaluation frequency: {}\n", gamma_grade_eval_freq); utils::logmesg(lmp, "Extrapolation grade evaluation frequency: {}\n", gamma_grade_eval_freq);
} }
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
set coeffs for one or more type pairs 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) if (narg < 5)
error->all(FLERR, error->all(FLERR,
"Incorrect args for pair coefficients. Correct form:\npair_coeff * * <potential.yaml> <potential.asi> elem1 elem2 ..."); "Incorrect args for pair coefficients. Correct form:\npair_coeff * * "
"<potential.yaml> <potential.asi> 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 potential_file_name = utils::get_potential_file_path(arg[2]);
auto active_set_inv_filename = utils::get_potential_file_path(arg[3]); auto active_set_inv_filename = utils::get_potential_file_path(arg[3]);
char **elemtypes = &arg[4]; char **elemtypes = &arg[4];
delete aceimpl->basis_set;
delete aceimpl->ctilde_basis_set;
delete aceimpl->basis_set; //load potential file
delete aceimpl->ctilde_basis_set; 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 //convert the basis set to CTilde format
aceimpl->basis_set = new ACEBBasisSet(); aceimpl->ctilde_basis_set = new ACECTildeBasisSet();
if (comm->me == 0) utils::logmesg(lmp, "Loading {}\n", potential_file_name); *aceimpl->ctilde_basis_set = aceimpl->basis_set->to_ACECTildeBasisSet();
aceimpl->basis_set->load(potential_file_name);
//convert the basis set to CTilde format if (comm->me == 0) {
aceimpl->ctilde_basis_set = new ACECTildeBasisSet(); utils::logmesg(lmp, "Total number of basis functions\n");
*aceimpl->ctilde_basis_set = aceimpl->basis_set->to_ACECTildeBasisSet();
if (comm->me == 0) { for (SPECIES_TYPE mu = 0; mu < aceimpl->basis_set->nelements; mu++) {
utils::logmesg(lmp, "Total number of basis functions\n"); int n_r1 = aceimpl->basis_set->total_basis_size_rank1[mu];
int n = aceimpl->basis_set->total_basis_size[mu];
for (SPECIES_TYPE mu = 0; mu < aceimpl->basis_set->nelements; mu++) { utils::logmesg(lmp, "\t{}: {} (r=1) {} (r>1)\n", aceimpl->basis_set->elements_name[mu], n_r1,
int n_r1 = aceimpl->basis_set->total_basis_size_rank1[mu]; n);
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 aceimpl->ace = new ACEBEvaluator();
// map[i] = which element the Ith atom type is, -1 if not mapped aceimpl->ace->element_type_mapping.init(atom->ntypes + 1);
// map[0] is not used
delete aceimpl->ace;
delete aceimpl->rec_ace;
aceimpl->ace = new ACEBEvaluator(); aceimpl->rec_ace = new ACERecursiveEvaluator();
aceimpl->ace->element_type_mapping.init(atom->ntypes + 1); 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(); FILE *species_type_file = nullptr;
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; const int n = atom->ntypes;
element_names.resize(n);
const int n = atom->ntypes; for (int i = 1; i <= n; i++) {
element_names.resize(n); char *elemname = elemtypes[i - 1];
for (int i = 1; i <= n; i++) { element_names[i - 1] = elemname;
char *elemname = elemtypes[i - 1]; if (strcmp(elemname, "NULL") == 0) {
element_names[i - 1] = elemname; // species_type=-1 value will not reach ACE Evaluator::compute_atom,
if (strcmp(elemname, "NULL") == 0) { // but if it will ,then error will be thrown there
// species_type=-1 value will not reach ACE Evaluator::compute_atom, aceimpl->ace->element_type_mapping(i) = -1;
// but if it will ,then error will be thrown there map[i] = -1;
aceimpl->ace->element_type_mapping(i) = -1; if (comm->me == 0) utils::logmesg(lmp, "Skipping LAMMPS atom type #{}(NULL)\n", i);
map[i] = -1; } else {
if (comm->me == 0) utils::logmesg(lmp, "Skipping LAMMPS atom type #{}(NULL)\n", i); // dump species types for reconstruction of atomic configurations
} else { int atomic_number = AtomicNumberByName_pace_al(elemname);
// dump species types for reconstruction of atomic configurations if (atomic_number == -1) error->all(FLERR, "'{}' is not a valid element\n", elemname);
int atomic_number = AtomicNumberByName_pace_al(elemname); SPECIES_TYPE mu = aceimpl->basis_set->get_species_index_by_name(elemname);
if (atomic_number == -1) error->all(FLERR, "'{}' is not a valid element\n", elemname); if (mu != -1) {
SPECIES_TYPE mu = aceimpl->basis_set->get_species_index_by_name(elemname); if (comm->me == 0)
if (mu != -1) { utils::logmesg(lmp, "Mapping LAMMPS atom type #{}({}) -> ACE species type #{}\n", i,
if (comm->me == 0) elemname, mu);
utils::logmesg(lmp, "Mapping LAMMPS atom type #{}({}) -> ACE species type #{}\n", i, map[i] = mu;
elemname, mu); // set up LAMMPS atom type to ACE species mapping for ace evaluators
map[i] = mu; aceimpl->ace->element_type_mapping(i) = mu;
// set up LAMMPS atom type to ACE species mapping for ace evaluators aceimpl->rec_ace->element_type_mapping(i) = mu;
aceimpl->ace->element_type_mapping(i) = mu; } else {
aceimpl->rec_ace->element_type_mapping(i) = mu; error->all(FLERR, "Element {} is not supported by ACE-potential from file {}", elemname,
} else { potential_file_name);
error->all(FLERR, "Element {} is not supported by ACE-potential from file {}", elemname, }
potential_file_name);
}
}
} }
}
aceimpl->ace->set_basis(*aceimpl->basis_set); aceimpl->ace->set_basis(*aceimpl->basis_set);
aceimpl->rec_ace->set_basis(*aceimpl->ctilde_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); if (comm->me == 0) utils::logmesg(lmp, "Loading ASI {}\n", active_set_inv_filename);
aceimpl->ace->load_active_set(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;
// 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 init specific to this pair style
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
void PairPACEExtrapolation::init_style() { void PairPACEExtrapolation::init_style()
if (atom->tag_enable == 0) {
error->all(FLERR, "Pair style PACE requires atom IDs"); if (atom->tag_enable == 0) error->all(FLERR, "Pair style PACE requires atom IDs");
if (force->newton_pair == 0) if (force->newton_pair == 0) error->all(FLERR, "Pair style PACE requires newton pair on");
error->all(FLERR, "Pair style PACE requires newton pair on");
// request a full neighbor list // request a full neighbor list
neighbor->add_request(this, NeighConst::REQ_FULL); neighbor->add_request(this, NeighConst::REQ_FULL);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
double PairPACEExtrapolation::init_one(int i, int 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 if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
scale[j][i] = scale[i][j]; //cutoff from the basis set's radial functions settings
return aceimpl->basis_set->radial_functions->cut(map[i], map[j]); 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 extract method for extracting value of scale variable
---------------------------------------------------------------------- */ ---------------------------------------------------------------------- */
void *PairPACEExtrapolation::extract(const char *str, int &dim) { void *PairPACEExtrapolation::extract(const char *str, int &dim)
dim = 2; {
if (strcmp(str, "scale") == 0) return (void *) scale; dim = 2;
return nullptr; if (strcmp(str, "scale") == 0) return (void *) scale;
return nullptr;
} }

77
src/ML-PACE/pair_pace_extrapolation.h
View File

@ -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. // Created by Lysogorskiy Yury on 1.01.22.
// //
#ifdef PAIR_CLASS #ifdef PAIR_CLASS
// clang-format off // clang-format off
PairStyle(pace/extrapolation,PairPACEExtrapolation) PairStyle(pace/extrapolation,PairPACEExtrapolation);
// clang-format on // clang-format on
#else #else
@ -33,53 +32,47 @@ PairStyle(pace/extrapolation,PairPACEExtrapolation)
namespace LAMMPS_NS { namespace LAMMPS_NS {
//forward declaration //forward declaration
class ComputePACEExtrapolation; class ComputePACEExtrapolation;
class DumpPACEExtrapolation; class DumpPACEExtrapolation;
class PairPACEExtrapolation : public Pair { class PairPACEExtrapolation : public Pair {
friend class ComputePACEExtrapolation; friend class ComputePACEExtrapolation;
friend class DumpPACEExtrapolation; friend class DumpPACEExtrapolation;
public: public:
PairPACEExtrapolation(class LAMMPS *); 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<std::string> 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; } // namespace LAMMPS_NS
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<std::string> 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;
};
}
#endif #endif
#endif #endif