PairPACEExtrapolation: rely on usage "fix pair"

This commit is contained in:
Yury Lysogorskiy
2022-09-01 11:00:20 +02:00
parent 4e5604daef
commit c63df4ac87
6 changed files with 343 additions and 573 deletions

View File

@ -1,90 +0,0 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
//
// Created by Yury Lysogorskiy on 23.06.22.
//
#include "compute_pace_extrapolation.h"
#include "pair_pace_extrapolation.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "modify.h"
#include "update.h"
using namespace LAMMPS_NS;
ComputePACEExtrapolation::ComputePACEExtrapolation(class LAMMPS *lmp, int narg, char **arg) :
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
}
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 ((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;
pair_pace_extrapolation->compute(0, 0);
pair_pace_extrapolation->is_set_energies_forces = true;
pair_pace_extrapolation->vflag_fdotr = old_vflag_fdotr;
}
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) {
invoke_compute_extrapolation_grades();
}
scalar = pair_pace_extrapolation->max_gamma_grade_per_structure;
return scalar;
}
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) {
invoke_compute_extrapolation_grades();
}
vector_atom = pair_pace_extrapolation->extrapolation_grade_gamma;
}

View File

@ -1,47 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
//
// Created by Yury Lysogorskiy on 23.06.22.
//
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(pace/extrapolation,ComputePACEExtrapolation);
// clang-format on
#else
#ifndef COMPUTE_PACE_H
#define COMPUTE_PACE_H
#include "compute.h"
#include "pair_pace_extrapolation.h"
namespace LAMMPS_NS {
class PairPACEExtrapolation;
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

View File

@ -1,56 +0,0 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "dump_pace_extrapolation.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "update.h"
#include "modify.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
DumpPACEExtrapolation::DumpPACEExtrapolation(struct LAMMPS *lmp, int nargs, char **argv) : DumpCustom(lmp, nargs,
argv) {
pairPaceExtrapolation = (PairPACEExtrapolation *) force->pair_match("pace/extrapolation", 1);
if (!pairPaceExtrapolation)
error->all(FLERR, "Dump pace/extrapolation requires a `pace/extrapolation` pair style");
// Save atomtypes to elements mapping into SPECIES_TYPE_FNAME to be used later in Python for loading extrapolative structures
FILE *species_type_file = fopen(SPECIES_TYPE_FNAME.c_str(), "w");
const int n = atom->ntypes;
for (int i = 0; i < n; i++) {
auto elemname = pairPaceExtrapolation->element_names[i].c_str();
fprintf(species_type_file, "%s ", elemname);
}
fclose(species_type_file);
}
void DumpPACEExtrapolation::write() {
int current_time_step = update->ntimestep;
// dump only if
// 1) extrapolation grades were computed on current timestep AND
// 2) max extrapolation grade > gamma_lower_bound
if (current_time_step == pairPaceExtrapolation->bevaluator_timestep) {
if (pairPaceExtrapolation->max_gamma_grade_per_structure > pairPaceExtrapolation->gamma_lower_bound) {
DumpCustom::write();
MPI_Barrier(world);
}
}
}

View File

@ -1,43 +0,0 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef DUMP_CLASS
// clang-format off
DumpStyle(pace/extrapolation,DumpPACEExtrapolation);
// clang-format on
#else
#ifndef LMP_DUMP_PACE_AL_H
#define LMP_DUMP_PACE_AL_H
#include "pair.h"
#include "dump_custom.h"
#include "pair_pace_extrapolation.h"
namespace LAMMPS_NS {
// forward declaration
class PairPACEExtrapolation;
class DumpPACEExtrapolation : public DumpCustom {
const std::string SPECIES_TYPE_FNAME = "species_types.dat";
PairPACEExtrapolation *pairPaceExtrapolation;
public:
DumpPACEExtrapolation(class LAMMPS *lmp, int nargs, char **argv);
void write() override;
};
} // namespace LAMMPS_NS
#endif
#endif

View File

@ -44,14 +44,13 @@ Copyright 2022 Yury Lysogorskiy^1, Anton Bochkarev^1, Matous Mrovec^1, Ralf Drau
#include "ace_b_evaluator.h"
#include "ace_recursive.h"
#include "ace_version.h"
#include "compute_pace_extrapolation.h"
//#include "compute_pace_extrapolation.h"
namespace LAMMPS_NS {
struct ACEALImpl {
ACEALImpl() : basis_set(nullptr), ace(nullptr), ctilde_basis_set(nullptr), rec_ace(nullptr) {}
~ACEALImpl()
{
~ACEALImpl() {
delete basis_set;
delete ace;
}
@ -76,16 +75,14 @@ static char const *const elements_pace_al[] = {
"Th", "Pa", "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm", "Md", "No", "Lr"};
static constexpr int elements_num_pace_al = sizeof(elements_pace_al) / sizeof(const char *);
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) return i;
return -1;
}
/* ---------------------------------------------------------------------- */
PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp)
{
PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp) {
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
@ -102,8 +99,7 @@ PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp)
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairPACEExtrapolation::~PairPACEExtrapolation()
{
PairPACEExtrapolation::~PairPACEExtrapolation() {
if (copymode) return;
delete aceimpl;
@ -119,8 +115,7 @@ PairPACEExtrapolation::~PairPACEExtrapolation()
/* ---------------------------------------------------------------------- */
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;
double fij[3];
@ -177,9 +172,9 @@ void PairPACEExtrapolation::compute(int eflag, int vflag)
}
bigint current_timestep = update->ntimestep;
bool is_bevaluator = (current_timestep - bevaluator_timestep_shift) % gamma_grade_eval_freq == 0;
// flag_compute_extrapolation_grade = (current_timestep - bevaluator_timestep_shift) % gamma_grade_eval_freq == 0;
if (is_bevaluator)
if (flag_compute_extrapolation_grade)
aceimpl->ace->resize_neighbours_cache(max_jnum);
else
aceimpl->rec_ace->resize_neighbours_cache(max_jnum);
@ -204,7 +199,7 @@ void PairPACEExtrapolation::compute(int eflag, int vflag)
// jnum(0) = 50
// jlist(neigh ind of 0-atom) = [1,2,10,7,99,25, .. 50 element in total]
try {
if (is_bevaluator)
if (flag_compute_extrapolation_grade)
aceimpl->ace->compute_atom(i, x, type, jnum, jlist);
else
aceimpl->rec_ace->compute_atom(i, x, type, jnum, jlist);
@ -213,17 +208,18 @@ void PairPACEExtrapolation::compute(int eflag, int vflag)
}
// 'compute_atom' will update the `ace->e_atom` and `ace->neighbours_forces(jj, alpha)` arrays and max_gamma_grade
if (is_bevaluator) {
if (flag_compute_extrapolation_grade) {
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;
// 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);
(flag_compute_extrapolation_grade ? 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];
@ -248,13 +244,13 @@ void PairPACEExtrapolation::compute(int eflag, int 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) {
if (eflag) {
// evdwl = energy of atom I
DOUBLE_TYPE e_atom;
if (is_bevaluator)
if (flag_compute_extrapolation_grade)
e_atom = aceimpl->ace->e_atom;
else
e_atom = aceimpl->rec_ace->e_atom;
@ -265,10 +261,11 @@ void PairPACEExtrapolation::compute(int eflag, int vflag)
if (vflag_fdotr) virial_fdotr_compute();
if (is_bevaluator) {
if (flag_compute_extrapolation_grade) {
//gather together max_gamma_grade_per_structure
MPI_Allreduce(&max_gamma_grade, &max_gamma_grade_per_structure, 1, MPI_DOUBLE, MPI_MAX, world);
// TODO: check, whether to stop here or externally in LAMMPS
// check if gamma_upper_bound is exceeded
if (max_gamma_grade_per_structure > gamma_upper_bound) {
if (comm->me == 0)
@ -288,8 +285,7 @@ void PairPACEExtrapolation::compute(int eflag, int vflag)
/* ---------------------------------------------------------------------- */
void PairPACEExtrapolation::allocate()
{
void PairPACEExtrapolation::allocate() {
allocated = 1;
int n = atom->ntypes;
@ -303,8 +299,7 @@ void PairPACEExtrapolation::allocate()
global settings
------------------------------------------------------------------------- */
void PairPACEExtrapolation::settings(int narg, char **arg)
{
void PairPACEExtrapolation::settings(int narg, char **arg) {
if (narg > 3) {
error->all(FLERR,
"Illegal pair_style command. Correct form:\n\tpair_style pace/al "
@ -329,17 +324,17 @@ void PairPACEExtrapolation::settings(int narg, char **arg)
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);
// utils::logmesg(lmp, "Extrapolation grade evaluation frequency: {}\n", gamma_grade_eval_freq);
}
}
@ -347,8 +342,7 @@ void PairPACEExtrapolation::settings(int narg, char **arg)
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,
@ -448,8 +442,7 @@ void PairPACEExtrapolation::coeff(int narg, char **arg)
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 (force->newton_pair == 0) error->all(FLERR, "Pair style PACE requires newton pair on");
@ -461,8 +454,7 @@ void PairPACEExtrapolation::init_style()
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
scale[j][i] = scale[i][j];
@ -472,9 +464,29 @@ double PairPACEExtrapolation::init_one(int i, int j)
/* ----------------------------------------------------------------------
extract method for extracting value of scale variable
---------------------------------------------------------------------- */
void *PairPACEExtrapolation::extract(const char *str, int &dim)
{
void *PairPACEExtrapolation::extract(const char *str, int &dim) {
//check if str=="gamma_flag" then compute extrapolation grades on this iteration
dim = 0;
if (strcmp(str, "gamma_flag")) return (void *) &flag_compute_extrapolation_grade;
dim = 2;
if (strcmp(str, "scale") == 0) return (void *) scale;
return nullptr;
}
/* ----------------------------------------------------------------------
peratom requests from FixPair
return ptr to requested data
also return ncol = # of quantites per atom
0 = per-atom vector
1 or more = # of columns in per-atom array
return NULL if str is not recognized
---------------------------------------------------------------------- */
void *PairPACEExtrapolation::extract_peratom(const char *str, int &ncol) {
if (strcmp(str, "gamma") == 0) {
ncol = 0;
return (void *) extrapolation_grade_gamma;
}
return nullptr;
}

View File

@ -20,7 +20,7 @@ Copyright 2022 Yury Lysogorskiy^1, Anton Bochkarev^1, Matous Mrovec^1, Ralf Drau
#ifdef PAIR_CLASS
// clang-format off
PairStyle(pace/extrapolation,PairPACEExtrapolation);
PairStyle(pace/extrapolation,PairPACEExtrapolation)
// clang-format on
#else
@ -32,14 +32,7 @@ PairStyle(pace/extrapolation,PairPACEExtrapolation);
namespace LAMMPS_NS {
//forward declaration
class ComputePACEExtrapolation;
class DumpPACEExtrapolation;
class PairPACEExtrapolation : public Pair {
friend class ComputePACEExtrapolation;
friend class DumpPACEExtrapolation;
public:
PairPACEExtrapolation(class LAMMPS *);
~PairPACEExtrapolation() override;
@ -50,11 +43,12 @@ class PairPACEExtrapolation : public Pair {
void init_style() override;
double init_one(int, int) override;
void *extract(const char *, int &) override;
void *extract_peratom(const char *, 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
// bigint gamma_grade_eval_freq = 1;
// bool is_set_energies_forces = true; // if set, then update forces and energies
int nmax;
double gamma_lower_bound = 1.5;
@ -65,10 +59,10 @@ class PairPACEExtrapolation : public Pair {
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
int flag_compute_extrapolation_grade;
double **scale;
};