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

new accessor APIs for fixes and computes in Modify plus a few applications

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-08-08 17:25:06 -04:00
parent ef04f6ca69
commit 5b40e4cb38
20 changed files with 397 additions and 388 deletions
Download Patch File Download Diff File Expand all files Collapse all files

55
src/GRANULAR/fix_pour.cpp
View File

@ -41,7 +41,6 @@ using namespace MathConst;
enum{ATOM,MOLECULE}; enum{ATOM,MOLECULE};
enum{ONE,RANGE,POLY}; enum{ONE,RANGE,POLY};
enum{CONSTANT,EQUAL}; // same as FixGravity
#define EPSILON 0.001 #define EPSILON 0.001
#define SMALL 1.0e-10 #define SMALL 1.0e-10
@ -186,10 +185,12 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
// grav = gravity in distance/time^2 units // grav = gravity in distance/time^2 units
// assume grav = -magnitude at this point, enforce in init() // assume grav = -magnitude at this point, enforce in init()
int ifix = modify->find_fix_by_style("^gravity"); auto fixlist = modify->get_fix_by_style("^gravity");
if (ifix == -1) if (fixlist.size() != 1)
error->all(FLERR,"No fix gravity defined for fix pour"); error->all(FLERR,"There must be exactly one fix gravity defined for fix pour");
grav = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v; auto fixgrav = (FixGravity *)fixlist.front();
grav = -fixgrav->magnitude * force->ftm2v;
// nfreq = timesteps between insertions // nfreq = timesteps between insertions
// should be time for a particle to fall from top of insertion region // should be time for a particle to fall from top of insertion region
@ -208,8 +209,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) :
v_relative = vy - rate; v_relative = vy - rate;
delta = yhi - ylo; delta = yhi - ylo;
} }
double t = double t = (-v_relative - sqrt(v_relative*v_relative - 2.0*grav*delta)) / grav;
(-v_relative - sqrt(v_relative*v_relative - 2.0*grav*delta)) / grav;
nfreq = static_cast<int>(t/update->dt + 0.5); nfreq = static_cast<int>(t/update->dt + 0.5);
// 1st insertion on next timestep // 1st insertion on next timestep
@ -316,17 +316,16 @@ void FixPour::init()
// for 3d must point in -z, for 2d must point in -y // for 3d must point in -z, for 2d must point in -y
// else insertion cannot work // else insertion cannot work
int ifix = modify->find_fix_by_style("^gravity"); auto fixlist = modify->get_fix_by_style("^gravity");
if (ifix == -1) if (fixlist.size() != 1)
error->all(FLERR,"No fix gravity defined for fix pour"); error->all(FLERR,"There must be exactly one fix gravity defined for fix pour");
auto fixgrav = (FixGravity *)fixlist.front();
int varflag = ((FixGravity *) modify->fix[ifix])->varflag; if (fixgrav->varflag != FixGravity::CONSTANT)
if (varflag != CONSTANT)
error->all(FLERR,"Fix gravity for fix pour must be constant"); error->all(FLERR,"Fix gravity for fix pour must be constant");
double xgrav = ((FixGravity *) modify->fix[ifix])->xgrav; double xgrav = fixgrav->xgrav;
double ygrav = ((FixGravity *) modify->fix[ifix])->ygrav; double ygrav = fixgrav->ygrav;
double zgrav = ((FixGravity *) modify->fix[ifix])->zgrav; double zgrav = fixgrav->zgrav;
if (domain->dimension == 3) { if (domain->dimension == 3) {
if (fabs(xgrav) > EPSILON || fabs(ygrav) > EPSILON || if (fabs(xgrav) > EPSILON || fabs(ygrav) > EPSILON ||
@ -338,37 +337,29 @@ void FixPour::init()
error->all(FLERR,"Gravity must point in -y to use with fix pour in 2d"); error->all(FLERR,"Gravity must point in -y to use with fix pour in 2d");
} }
double gnew = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v; double gnew = -fixgrav->magnitude * force->ftm2v;
if (gnew != grav) if (gnew != grav) error->all(FLERR,"Gravity changed since fix pour was created");
error->all(FLERR,"Gravity changed since fix pour was created");
// if rigidflag defined, check for rigid/small fix // if rigidflag defined, check for rigid/small fix
// its molecule template must be same as this one // its molecule template must be same as this one
fixrigid = nullptr; fixrigid = modify->get_fix_by_id(idrigid);
if (rigidflag) { if (rigidflag) {
int ifix = modify->find_fix(idrigid); if (!fixrigid) error->all(FLERR,"Fix pour rigid fix does not exist");
if (ifix < 0) error->all(FLERR,"Fix pour rigid fix does not exist");
fixrigid = modify->fix[ifix];
int tmp; int tmp;
if (onemols != (Molecule **) fixrigid->extract("onemol",tmp)) if (onemols != (Molecule **) fixrigid->extract("onemol",tmp))
error->all(FLERR, error->all(FLERR,"Fix pour and fix rigid/small not using same molecule template ID");
"Fix pour and fix rigid/small not using "
"same molecule template ID");
} }
// if shakeflag defined, check for SHAKE fix // if shakeflag defined, check for SHAKE fix
// its molecule template must be same as this one // its molecule template must be same as this one
fixshake = nullptr; fixshake = modify->get_fix_by_id(idshake);
if (shakeflag) { if (shakeflag) {
int ifix = modify->find_fix(idshake); if (!fixshake) error->all(FLERR,"Fix pour shake fix does not exist");
if (ifix < 0) error->all(FLERR,"Fix pour shake fix does not exist");
fixshake = modify->fix[ifix];
int tmp; int tmp;
if (onemols != (Molecule **) fixshake->extract("onemol",tmp)) if (onemols != (Molecule **) fixshake->extract("onemol",tmp))
error->all(FLERR,"Fix pour and fix shake not using " error->all(FLERR,"Fix pour and fix shake not using same molecule template ID");
"same molecule template ID");
} }
} }

2
src/KOKKOS/fix_gravity_kokkos.cpp
View File

@ -24,8 +24,6 @@
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
enum{CONSTANT,EQUAL};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
template<class DeviceType> template<class DeviceType>

1
src/OPENMP/fix_gravity_omp.cpp
View File

@ -30,7 +30,6 @@ using namespace LAMMPS_NS;
using namespace FixConst; using namespace FixConst;
enum{CHUTE,SPHERICAL,GRADIENT,VECTOR}; enum{CHUTE,SPHERICAL,GRADIENT,VECTOR};
enum{CONSTANT,EQUAL};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

11
src/atom.cpp
View File

@ -1700,8 +1700,7 @@ void Atom::data_bodies(int n, char *buf, AtomVec *avec_body, tagint id_offset)
void Atom::data_fix_compute_variable(int nprev, int nnew) void Atom::data_fix_compute_variable(int nprev, int nnew)
{ {
for (int m = 0; m < modify->nfix; m++) { for (const auto &fix : modify->get_fix_list()) {
Fix *fix = modify->fix[m];
if (fix->create_attribute) if (fix->create_attribute)
for (int i = nprev; i < nnew; i++) for (int i = nprev; i < nnew; i++)
fix->set_arrays(i); fix->set_arrays(i);
@ -2238,15 +2237,13 @@ void Atom::setup_sort_bins()
#ifdef LMP_GPU #ifdef LMP_GPU
if (userbinsize == 0.0) { if (userbinsize == 0.0) {
int ifix = modify->find_fix("package_gpu"); FixGPU *fix = (FixGPU *)modify->get_fix_by_id("package_gpu");
if (ifix >= 0) { if (fix) {
const double subx = domain->subhi[0] - domain->sublo[0]; const double subx = domain->subhi[0] - domain->sublo[0];
const double suby = domain->subhi[1] - domain->sublo[1]; const double suby = domain->subhi[1] - domain->sublo[1];
const double subz = domain->subhi[2] - domain->sublo[2]; const double subz = domain->subhi[2] - domain->sublo[2];
FixGPU *fix = static_cast<FixGPU *>(modify->fix[ifix]); binsize = fix->binsize(subx, suby, subz, atom->nlocal,neighbor->cutneighmax);
binsize = fix->binsize(subx, suby, subz, atom->nlocal,
neighbor->cutneighmax);
bininv = 1.0 / binsize; bininv = 1.0 / binsize;
nbinx = static_cast<int> (ceil(subx * bininv)); nbinx = static_cast<int> (ceil(subx * bininv));

7
src/balance.cpp
View File

@ -496,11 +496,8 @@ void Balance::weight_storage(char *prefix)
if (prefix) cmd = prefix; if (prefix) cmd = prefix;
cmd += "IMBALANCE_WEIGHTS"; cmd += "IMBALANCE_WEIGHTS";
int ifix = modify->find_fix(cmd); fixstore = (FixStore *) modify->get_fix_by_id(cmd);
if (ifix < 1) { if (!fixstore) fixstore = (FixStore *) modify->add_fix(cmd + " all STORE peratom 0 1");
cmd += " all STORE peratom 0 1";
fixstore = (FixStore *) modify->add_fix(cmd);
} else fixstore = (FixStore *) modify->fix[ifix];
// do not carry weights with atoms during normal atom migration // do not carry weights with atoms during normal atom migration

25
src/change_box.cpp
View File

@ -291,16 +291,12 @@ void ChangeBox::command(int narg, char **arg)
} else if (ops[m].style == ORTHO) { } else if (ops[m].style == ORTHO) {
if (domain->xy != 0.0 || domain->yz != 0.0 || domain->xz != 0.0) if (domain->xy != 0.0 || domain->yz != 0.0 || domain->xz != 0.0)
error->all(FLERR, error->all(FLERR,"Cannot change box to orthogonal when tilt is non-zero");
"Cannot change box to orthogonal when tilt is non-zero");
if (output->ndump) if (output->ndump)
error->all(FLERR, error->all(FLERR,"Cannot change box ortho/triclinic with dumps defined");
"Cannot change box ortho/triclinic with dumps defined"); for (const auto &fix : modify->get_fix_list())
for (i = 0; i < modify->nfix; i++) if (fix->no_change_box)
if (modify->fix[i]->no_change_box) error->all(FLERR,"Cannot change box ortho/triclinic with certain fixes defined");
error->all(FLERR,
"Cannot change box ortho/triclinic with "
"certain fixes defined");
domain->triclinic = 0; domain->triclinic = 0;
domain->set_initial_box(); domain->set_initial_box();
domain->set_global_box(); domain->set_global_box();
@ -309,13 +305,10 @@ void ChangeBox::command(int narg, char **arg)
} else if (ops[m].style == TRICLINIC) { } else if (ops[m].style == TRICLINIC) {
if (output->ndump) if (output->ndump)
error->all(FLERR, error->all(FLERR,"Cannot change box ortho/triclinic with dumps defined");
"Cannot change box ortho/triclinic with dumps defined"); for (const auto &fix : modify->get_fix_list())
for (i = 0; i < modify->nfix; i++) if (fix->no_change_box)
if (modify->fix[i]->no_change_box) error->all(FLERR,"Cannot change box ortho/triclinic with certain fixes defined");
error->all(FLERR,
"Cannot change box ortho/triclinic with "
"certain fixes defined");
domain->triclinic = 1; domain->triclinic = 1;
domain->set_lamda_box(); domain->set_lamda_box();
domain->set_initial_box(); domain->set_initial_box();

27
src/compute_angmom_chunk.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -29,8 +28,8 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeAngmomChunk::ComputeAngmomChunk(LAMMPS *lmp, int narg, char **arg) : ComputeAngmomChunk::ComputeAngmomChunk(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg), idchunk(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr),
idchunk(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr), comall(nullptr), angmom(nullptr), angmomall(nullptr) comall(nullptr), angmom(nullptr), angmomall(nullptr)
{ {
if (narg != 4) error->all(FLERR, "Illegal compute angmom/chunk command"); if (narg != 4) error->all(FLERR, "Illegal compute angmom/chunk command");
@ -70,11 +69,8 @@ ComputeAngmomChunk::~ComputeAngmomChunk()
void ComputeAngmomChunk::init() void ComputeAngmomChunk::init()
{ {
int icompute = modify->find_compute(idchunk); cchunk = (ComputeChunkAtom *) modify->get_compute_by_id(idchunk);
if (icompute < 0) if (!cchunk) error->all(FLERR, "Chunk/atom compute does not exist for compute angmom/chunk");
error->all(FLERR,"Chunk/atom compute does not exist for "
"compute angmom/chunk");
cchunk = (ComputeChunkAtom *) modify->compute[icompute];
if (strcmp(cchunk->style, "chunk/atom") != 0) if (strcmp(cchunk->style, "chunk/atom") != 0)
error->all(FLERR, "Compute angmom/chunk does not use chunk/atom compute"); error->all(FLERR, "Compute angmom/chunk does not use chunk/atom compute");
} }
@ -122,8 +118,10 @@ void ComputeAngmomChunk::compute_array()
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
index = ichunk[i] - 1; index = ichunk[i] - 1;
if (index < 0) continue; if (index < 0) continue;
if (rmass) massone = rmass[i]; if (rmass)
else massone = mass[type[i]]; massone = rmass[i];
else
massone = mass[type[i]];
domain->unmap(x[i], image[i], unwrap); domain->unmap(x[i], image[i], unwrap);
massproc[index] += massone; massproc[index] += massone;
com[index][0] += unwrap[0] * massone; com[index][0] += unwrap[0] * massone;
@ -154,15 +152,16 @@ void ComputeAngmomChunk::compute_array()
dx = unwrap[0] - comall[index][0]; dx = unwrap[0] - comall[index][0];
dy = unwrap[1] - comall[index][1]; dy = unwrap[1] - comall[index][1];
dz = unwrap[2] - comall[index][2]; dz = unwrap[2] - comall[index][2];
if (rmass) massone = rmass[i]; if (rmass)
else massone = mass[type[i]]; massone = rmass[i];
else
massone = mass[type[i]];
angmom[index][0] += massone * (dy * v[i][2] - dz * v[i][1]); angmom[index][0] += massone * (dy * v[i][2] - dz * v[i][1]);
angmom[index][1] += massone * (dz * v[i][0] - dx * v[i][2]); angmom[index][1] += massone * (dz * v[i][0] - dx * v[i][2]);
angmom[index][2] += massone * (dx * v[i][1] - dy * v[i][0]); angmom[index][2] += massone * (dx * v[i][1] - dy * v[i][0]);
} }
MPI_Allreduce(&angmom[0][0],&angmomall[0][0],3*nchunk, MPI_Allreduce(&angmom[0][0], &angmomall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world);
MPI_DOUBLE,MPI_SUM,world);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------

129
src/compute_centroid_stress_atom.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -13,21 +12,23 @@
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "compute_centroid_stress_atom.h" #include "compute_centroid_stress_atom.h"
#include <cstring>
#include "atom.h"
#include "update.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h" #include "angle.h"
#include "atom.h"
#include "bond.h"
#include "comm.h"
#include "dihedral.h" #include "dihedral.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "improper.h" #include "improper.h"
#include "kspace.h" #include "kspace.h"
#include "modify.h"
#include "fix.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "modify.h"
#include "pair.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
@ -36,8 +37,7 @@ enum{NOBIAS,BIAS};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeCentroidStressAtom::ComputeCentroidStressAtom(LAMMPS *lmp, int narg, char **arg) : ComputeCentroidStressAtom::ComputeCentroidStressAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg), id_temp(nullptr), stress(nullptr)
id_temp(nullptr), stress(nullptr)
{ {
if (narg < 4) error->all(FLERR, "Illegal compute centroid/stress/atom command"); if (narg < 4) error->all(FLERR, "Illegal compute centroid/stress/atom command");
@ -50,17 +50,16 @@ ComputeCentroidStressAtom::ComputeCentroidStressAtom(LAMMPS *lmp, int narg, char
// store temperature ID used by stress computation // store temperature ID used by stress computation
// insure it is valid for temperature computation // insure it is valid for temperature computation
if (strcmp(arg[3],"NULL") == 0) id_temp = nullptr; if (strcmp(arg[3], "NULL") == 0)
id_temp = nullptr;
else { else {
id_temp = utils::strdup(arg[3]); id_temp = utils::strdup(arg[3]);
int icompute = modify->find_compute(id_temp); auto compute = modify->get_compute_by_id(id_temp);
if (icompute < 0) if (!compute)
error->all(FLERR,"Could not find compute centroid/stress/atom temperature ID"); error->all(FLERR, "Could not find compute centroid/stress/atom temperature ID {}", id_temp);
if (modify->compute[icompute]->tempflag == 0) if (compute->tempflag == 0)
error->all(FLERR, error->all(FLERR, "Compute centroid/stress/atom temperature ID does not compute temperature");
"Compute centroid/stress/atom temperature ID does not "
"compute temperature");
} }
// process optional args // process optional args
@ -79,19 +78,28 @@ ComputeCentroidStressAtom::ComputeCentroidStressAtom(LAMMPS *lmp, int narg, char
fixflag = 0; fixflag = 0;
int iarg = 4; int iarg = 4;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"ke") == 0) keflag = 1; if (strcmp(arg[iarg], "ke") == 0)
else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1; keflag = 1;
else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1; else if (strcmp(arg[iarg], "pair") == 0)
else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1; pairflag = 1;
else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1; else if (strcmp(arg[iarg], "bond") == 0)
else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1; bondflag = 1;
else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1; else if (strcmp(arg[iarg], "angle") == 0)
else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1; angleflag = 1;
else if (strcmp(arg[iarg], "dihedral") == 0)
dihedralflag = 1;
else if (strcmp(arg[iarg], "improper") == 0)
improperflag = 1;
else if (strcmp(arg[iarg], "kspace") == 0)
kspaceflag = 1;
else if (strcmp(arg[iarg], "fix") == 0)
fixflag = 1;
else if (strcmp(arg[iarg], "virial") == 0) { else if (strcmp(arg[iarg], "virial") == 0) {
pairflag = 1; pairflag = 1;
bondflag = angleflag = dihedralflag = improperflag = 1; bondflag = angleflag = dihedralflag = improperflag = 1;
kspaceflag = fixflag = 1; kspaceflag = fixflag = 1;
} else error->all(FLERR,"Illegal compute centroid/stress/atom command"); } else
error->all(FLERR, "Illegal compute centroid/stress/atom command");
iarg++; iarg++;
} }
} }
@ -115,13 +123,15 @@ void ComputeCentroidStressAtom::init()
// fixes could have changed or compute_modify could have changed it // fixes could have changed or compute_modify could have changed it
if (id_temp) { if (id_temp) {
int icompute = modify->find_compute(id_temp); temperature = modify->get_compute_by_id(id_temp);
if (icompute < 0) if (!temperature)
error->all(FLERR,"Could not find compute centroid/stress/atom temperature ID"); error->all(FLERR, "Could not find compute centroid/stress/atom temperature ID {}",id_temp);
temperature = modify->compute[icompute]; if (temperature->tempbias)
if (temperature->tempbias) biasflag = BIAS; biasflag = BIAS;
else biasflag = NOBIAS; else
} else biasflag = NOBIAS; biasflag = NOBIAS;
} else
biasflag = NOBIAS;
// check if force components and fixes support centroid atom stress // check if force components and fixes support centroid atom stress
// all bond styles support it as CENTROID_SAME // all bond styles support it as CENTROID_SAME
@ -194,8 +204,7 @@ void ComputeCentroidStressAtom::compute_peratom()
// clear local stress array // clear local stress array
for (i = 0; i < ntotal; i++) for (i = 0; i < ntotal; i++)
for (j = 0; j < 9; j++) for (j = 0; j < 9; j++) stress[i][j] = 0.0;
stress[i][j] = 0.0;
// add in per-atom contributions from all force components and fixes // add in per-atom contributions from all force components and fixes
@ -205,15 +214,12 @@ void ComputeCentroidStressAtom::compute_peratom()
if (force->pair->centroidstressflag == CENTROID_AVAIL) { if (force->pair->centroidstressflag == CENTROID_AVAIL) {
double **cvatom = force->pair->cvatom; double **cvatom = force->pair->cvatom;
for (i = 0; i < npair; i++) for (i = 0; i < npair; i++)
for (j = 0; j < 9; j++) for (j = 0; j < 9; j++) stress[i][j] += cvatom[i][j];
stress[i][j] += cvatom[i][j];
} else { } else {
double **vatom = force->pair->vatom; double **vatom = force->pair->vatom;
for (i = 0; i < npair; i++) { for (i = 0; i < npair; i++) {
for (j = 0; j < 6; j++) for (j = 0; j < 6; j++) stress[i][j] += vatom[i][j];
stress[i][j] += vatom[i][j]; for (j = 6; j < 9; j++) stress[i][j] += vatom[i][j - 3];
for (j = 6; j < 9; j++)
stress[i][j] += vatom[i][j-3];
} }
} }
} }
@ -226,41 +232,34 @@ void ComputeCentroidStressAtom::compute_peratom()
if (bondflag && force->bond) { if (bondflag && force->bond) {
double **vatom = force->bond->vatom; double **vatom = force->bond->vatom;
for (i = 0; i < nbond; i++) { for (i = 0; i < nbond; i++) {
for (j = 0; j < 6; j++) for (j = 0; j < 6; j++) stress[i][j] += vatom[i][j];
stress[i][j] += vatom[i][j]; for (j = 6; j < 9; j++) stress[i][j] += vatom[i][j - 3];
for (j = 6; j < 9; j++)
stress[i][j] += vatom[i][j-3];
} }
} }
if (angleflag && force->angle) { if (angleflag && force->angle) {
double **cvatom = force->angle->cvatom; double **cvatom = force->angle->cvatom;
for (i = 0; i < nbond; i++) for (i = 0; i < nbond; i++)
for (j = 0; j < 9; j++) for (j = 0; j < 9; j++) stress[i][j] += cvatom[i][j];
stress[i][j] += cvatom[i][j];
} }
if (dihedralflag && force->dihedral) { if (dihedralflag && force->dihedral) {
double **cvatom = force->dihedral->cvatom; double **cvatom = force->dihedral->cvatom;
for (i = 0; i < nbond; i++) for (i = 0; i < nbond; i++)
for (j = 0; j < 9; j++) for (j = 0; j < 9; j++) stress[i][j] += cvatom[i][j];
stress[i][j] += cvatom[i][j];
} }
if (improperflag && force->improper) { if (improperflag && force->improper) {
double **cvatom = force->improper->cvatom; double **cvatom = force->improper->cvatom;
for (i = 0; i < nbond; i++) for (i = 0; i < nbond; i++)
for (j = 0; j < 9; j++) for (j = 0; j < 9; j++) stress[i][j] += cvatom[i][j];
stress[i][j] += cvatom[i][j];
} }
if (kspaceflag && force->kspace && force->kspace->compute_flag) { if (kspaceflag && force->kspace && force->kspace->compute_flag) {
double **vatom = force->kspace->vatom; double **vatom = force->kspace->vatom;
for (i = 0; i < nkspace; i++) { for (i = 0; i < nkspace; i++) {
for (j = 0; j < 6; j++) for (j = 0; j < 6; j++) stress[i][j] += vatom[i][j];
stress[i][j] += vatom[i][j]; for (j = 6; j < 9; j++) stress[i][j] += vatom[i][j - 3];
for (j = 6; j < 9; j++)
stress[i][j] += vatom[i][j-3];
} }
} }
@ -279,18 +278,15 @@ void ComputeCentroidStressAtom::compute_peratom()
double **vatom = fix[ifix]->vatom; double **vatom = fix[ifix]->vatom;
if (vatom) if (vatom)
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
for (j = 0; j < 6; j++) for (j = 0; j < 6; j++) stress[i][j] += vatom[i][j];
stress[i][j] += vatom[i][j]; for (j = 6; j < 9; j++) stress[i][j] += vatom[i][j - 3];
for (j = 6; j < 9; j++)
stress[i][j] += vatom[i][j-3];
} }
} }
} }
// communicate ghost virials between neighbor procs // communicate ghost virials between neighbor procs
if (force->newton || if (force->newton || (force->kspace && force->kspace->tip4pflag && force->kspace->compute_flag))
(force->kspace && force->kspace->tip4pflag && force->kspace->compute_flag))
comm->reverse_comm_compute(this); comm->reverse_comm_compute(this);
// zero virial of atoms not in group // zero virial of atoms not in group
@ -359,8 +355,7 @@ void ComputeCentroidStressAtom::compute_peratom()
// invoke temperature if it hasn't been already // invoke temperature if it hasn't been already
// this insures bias factor is pre-computed // this insures bias factor is pre-computed
if (keflag && temperature->invoked_scalar != update->ntimestep) if (keflag && temperature->invoked_scalar != update->ntimestep) temperature->compute_scalar();
temperature->compute_scalar();
if (rmass) { if (rmass) {
for (i = 0; i < nlocal; i++) for (i = 0; i < nlocal; i++)

72
src/compute_chunk_atom.cpp
View File

@ -316,40 +316,30 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
} }
if (which == ArgInfo::COMPUTE) { if (which == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(cfvid); cchunk = modify->get_compute_by_id(cfvid);
if (icompute < 0) if (!cchunk)
error->all(FLERR,"Compute ID for compute chunk /atom does not exist"); error->all(FLERR,"Compute ID {} for compute chunk /atom does not exist",cfvid);
if (modify->compute[icompute]->peratom_flag == 0) if (cchunk->peratom_flag == 0)
error->all(FLERR, error->all(FLERR,"Compute chunk/atom compute does not calculate per-atom values");
"Compute chunk/atom compute does not calculate " if ((argindex == 0) && (cchunk->size_peratom_cols != 0))
"per-atom values"); error->all(FLERR,"Compute chunk/atom compute does not calculate a per-atom vector");
if (argindex == 0 && if (argindex && (cchunk->size_peratom_cols == 0))
modify->compute[icompute]->size_peratom_cols != 0) error->all(FLERR,"Compute chunk/atom compute does not calculate a per-atom array");
error->all(FLERR,"Compute chunk/atom compute does not " if (argindex && argindex > cchunk->size_peratom_cols)
"calculate a per-atom vector"); error->all(FLERR,"Compute chunk/atom compute array is accessed out-of-range");
if (argindex && modify->compute[icompute]->size_peratom_cols == 0)
error->all(FLERR,"Compute chunk/atom compute does not "
"calculate a per-atom array");
if (argindex &&
argindex > modify->compute[icompute]->size_peratom_cols)
error->all(FLERR,"Compute chunk/atom compute array is "
"accessed out-of-range");
} }
if (which == ArgInfo::FIX) { if (which == ArgInfo::FIX) {
int ifix = modify->find_fix(cfvid); fchunk = modify->get_fix_by_id(cfvid);
if (ifix < 0) if (!fchunk)
error->all(FLERR,"Fix ID for compute chunk/atom does not exist"); error->all(FLERR,"Fix ID {} for compute chunk/atom does not exist",cfvid);
if (modify->fix[ifix]->peratom_flag == 0) if (fchunk->peratom_flag == 0)
error->all(FLERR,"Compute chunk/atom fix does not calculate " error->all(FLERR,"Compute chunk/atom fix does not calculate per-atom values");
"per-atom values"); if (argindex == 0 && fchunk->size_peratom_cols != 0)
if (argindex == 0 && modify->fix[ifix]->size_peratom_cols != 0) error->all(FLERR,"Compute chunk/atom fix does not calculate a per-atom vector");
error->all(FLERR, if (argindex && fchunk->size_peratom_cols == 0)
"Compute chunk/atom fix does not calculate a per-atom vector"); error->all(FLERR,"Compute chunk/atom fix does not calculate a per-atom array");
if (argindex && modify->fix[ifix]->size_peratom_cols == 0) if (argindex && argindex > fchunk->size_peratom_cols)
error->all(FLERR,
"Compute chunk/atom fix does not calculate a per-atom array");
if (argindex && argindex > modify->fix[ifix]->size_peratom_cols)
error->all(FLERR,"Compute chunk/atom fix array is accessed out-of-range"); error->all(FLERR,"Compute chunk/atom fix array is accessed out-of-range");
} }
@ -358,16 +348,14 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) :
if (ivariable < 0) if (ivariable < 0)
error->all(FLERR,"Variable name for compute chunk/atom does not exist"); error->all(FLERR,"Variable name for compute chunk/atom does not exist");
if (input->variable->atomstyle(ivariable) == 0) if (input->variable->atomstyle(ivariable) == 0)
error->all(FLERR,"Compute chunk/atom variable is not " error->all(FLERR,"Compute chunk/atom variable is not atom-style variable");
"atom-style variable");
} }
// setup scaling // setup scaling
if (binflag) { if (binflag) {
if (domain->triclinic == 1 && scaleflag != REDUCED) if (domain->triclinic == 1 && scaleflag != REDUCED)
error->all(FLERR,"Compute chunk/atom for triclinic boxes " error->all(FLERR,"Compute chunk/atom for triclinic boxes requires units reduced");
"requires units reduced");
} }
if (scaleflag == LATTICE) { if (scaleflag == LATTICE) {
@ -501,15 +489,13 @@ void ComputeChunkAtom::init()
// set compute,fix,variable // set compute,fix,variable
if (which == ArgInfo::COMPUTE) { if (which == ArgInfo::COMPUTE) {
int icompute = modify->find_compute(cfvid); cchunk = modify->get_compute_by_id(cfvid);
if (icompute < 0) if (!cchunk)
error->all(FLERR,"Compute ID for compute chunk/atom does not exist"); error->all(FLERR,"Compute ID {} for compute chunk/atom does not exist",cfvid);
cchunk = modify->compute[icompute];
} else if (which == ArgInfo::FIX) { } else if (which == ArgInfo::FIX) {
int ifix = modify->find_fix(cfvid); fchunk = modify->get_fix_by_id(cfvid);
if (ifix < 0) if (!fchunk)
error->all(FLERR,"Fix ID for compute chunk/atom does not exist"); error->all(FLERR,"Fix ID {} for compute chunk/atom does not exist",cfvid);
fchunk = modify->fix[ifix];
} else if (which == ArgInfo::VARIABLE) { } else if (which == ArgInfo::VARIABLE) {
int ivariable = input->variable->find(cfvid); int ivariable = input->variable->find(cfvid);
if (ivariable < 0) if (ivariable < 0)

30
src/domain.cpp
View File

@ -140,19 +140,19 @@ void Domain::init()
int box_change_x=0, box_change_y=0, box_change_z=0; int box_change_x=0, box_change_y=0, box_change_z=0;
int box_change_yz=0, box_change_xz=0, box_change_xy=0; int box_change_yz=0, box_change_xz=0, box_change_xy=0;
Fix **fixes = modify->fix; const auto &fixes = modify->get_fix_list();
if (nonperiodic == 2) box_change_size = 1; if (nonperiodic == 2) box_change_size = 1;
for (int i = 0; i < modify->nfix; i++) { for (const auto &fix : fixes) {
if (fixes[i]->box_change & Fix::BOX_CHANGE_SIZE) box_change_size = 1; if (fix->box_change & Fix::BOX_CHANGE_SIZE) box_change_size = 1;
if (fixes[i]->box_change & Fix::BOX_CHANGE_SHAPE) box_change_shape = 1; if (fix->box_change & Fix::BOX_CHANGE_SHAPE) box_change_shape = 1;
if (fixes[i]->box_change & Fix::BOX_CHANGE_DOMAIN) box_change_domain = 1; if (fix->box_change & Fix::BOX_CHANGE_DOMAIN) box_change_domain = 1;
if (fixes[i]->box_change & Fix::BOX_CHANGE_X) box_change_x++; if (fix->box_change & Fix::BOX_CHANGE_X) box_change_x++;
if (fixes[i]->box_change & Fix::BOX_CHANGE_Y) box_change_y++; if (fix->box_change & Fix::BOX_CHANGE_Y) box_change_y++;
if (fixes[i]->box_change & Fix::BOX_CHANGE_Z) box_change_z++; if (fix->box_change & Fix::BOX_CHANGE_Z) box_change_z++;
if (fixes[i]->box_change & Fix::BOX_CHANGE_YZ) box_change_yz++; if (fix->box_change & Fix::BOX_CHANGE_YZ) box_change_yz++;
if (fixes[i]->box_change & Fix::BOX_CHANGE_XZ) box_change_xz++; if (fix->box_change & Fix::BOX_CHANGE_XZ) box_change_xz++;
if (fixes[i]->box_change & Fix::BOX_CHANGE_XY) box_change_xy++; if (fix->box_change & Fix::BOX_CHANGE_XY) box_change_xy++;
} }
std::string mesg = "Must not have multiple fixes change box parameter "; std::string mesg = "Must not have multiple fixes change box parameter ";
@ -174,12 +174,12 @@ void Domain::init()
// check for fix deform // check for fix deform
deform_flag = deform_vremap = deform_groupbit = 0; deform_flag = deform_vremap = deform_groupbit = 0;
for (int i = 0; i < modify->nfix; i++) for (const auto &fix : fixes)
if (utils::strmatch(modify->fix[i]->style,"^deform")) { if (utils::strmatch(fix->style,"^deform")) {
deform_flag = 1; deform_flag = 1;
if (((FixDeform *) modify->fix[i])->remapflag == Domain::V_REMAP) { if (((FixDeform *) fix)->remapflag == Domain::V_REMAP) {
deform_vremap = 1; deform_vremap = 1;
deform_groupbit = modify->fix[i]->groupbit; deform_groupbit = fix->groupbit;
} }
} }

4
src/dump.cpp
View File

@ -248,8 +248,8 @@ void Dump::init()
reorderflag = 0; reorderflag = 0;
int gcmcflag = 0; int gcmcflag = 0;
for (int i = 0; i < modify->nfix; i++) for (const auto &fix : modify->get_fix_list())
if ((strcmp(modify->fix[i]->style,"gcmc") == 0)) if (utils::strmatch(fix->style,"^gcmc"))
gcmcflag = 1; gcmcflag = 1;
if (sortcol == 0 && atom->tag_consecutive() && !gcmcflag) { if (sortcol == 0 && atom->tag_consecutive() && !gcmcflag) {

5
src/finish.cpp
View File

@ -364,14 +364,13 @@ void Finish::end(int flag)
} }
#ifdef LMP_OPENMP #ifdef LMP_OPENMP
int ifix = modify->find_fix("package_omp"); FixOMP *fixomp = (FixOMP *) modify->get_fix_by_id("package_omp");
// print thread breakdown only with full timer detail // print thread breakdown only with full timer detail
if ((ifix >= 0) && timer->has_full() && me == 0) { if (fixomp && timer->has_full() && me == 0) {
double thr_total = 0.0; double thr_total = 0.0;
ThrData *td; ThrData *td;
FixOMP *fixomp = static_cast<FixOMP *>(lmp->modify->fix[ifix]);
for (i=0; i < nthreads; ++i) { for (i=0; i < nthreads; ++i) {
td = fixomp->get_thr(i); td = fixomp->get_thr(i);
thr_total += td->get_time(Timer::ALL); thr_total += td->get_time(Timer::ALL);

1
src/fix_gravity.cpp
View File

@ -32,7 +32,6 @@ using namespace FixConst;
using namespace MathConst; using namespace MathConst;
enum{CHUTE,SPHERICAL,VECTOR}; enum{CHUTE,SPHERICAL,VECTOR};
enum{CONSTANT,EQUAL}; // same as FixPour
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */

1
src/fix_gravity.h
View File

@ -37,6 +37,7 @@ class FixGravity : public Fix {
virtual void post_force_respa(int, int, int); virtual void post_force_respa(int, int, int);
double compute_scalar(); double compute_scalar();
void *extract(const char *, int &); void *extract(const char *, int &);
enum { CONSTANT, EQUAL };
protected: protected:
int style, disable; int style, disable;

78
src/modify.cpp
View File

@ -1103,6 +1103,45 @@ int Modify::find_fix_by_style(const char *style)
return -1; return -1;
} }
/* ----------------------------------------------------------------------
look up pointer to Fix class by fix-ID
return null pointer if ID not found
------------------------------------------------------------------------- */
Fix *Modify::get_fix_by_id(const std::string &id) const
{
if (id.empty()) return nullptr;
for (int ifix = 0; ifix < nfix; ifix++)
if (id == fix[ifix]->id) return fix[ifix];
return nullptr;
}
/* ----------------------------------------------------------------------
look up pointer to fixes by fix style name
return vector of matching pointers
------------------------------------------------------------------------- */
const std::vector<Fix *> Modify::get_fix_by_style(const std::string &style) const
{
std::vector<Fix *> matches;
if (style.empty()) return matches;
for (int ifix = 0; ifix < nfix; ifix++)
if (utils::strmatch(fix[ifix]->style,style)) matches.push_back(fix[ifix]);
return matches;
}
/* ----------------------------------------------------------------------
return list of fixes as vector
------------------------------------------------------------------------- */
const std::vector<Fix *> &Modify::get_fix_list()
{
fix_list = std::vector<Fix *>(fix, fix+nfix);
return fix_list;
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
check for fix associated with package name in compiled list check for fix associated with package name in compiled list
return 1 if found else 0 return 1 if found else 0
@ -1359,6 +1398,45 @@ int Modify::find_compute_by_style(const char *style)
return -1; return -1;
} }
/* ----------------------------------------------------------------------
look up pointer to Compute class by compute-ID
return null pointer if ID not found
------------------------------------------------------------------------- */
Compute *Modify::get_compute_by_id(const std::string &id) const
{
if (id.empty()) return nullptr;
for (int icompute = 0; icompute < ncompute; icompute++)
if (id == compute[icompute]->id) return compute[icompute];
return nullptr;
}
/* ----------------------------------------------------------------------
look up pointer to compute by compute style name
return null pointer if style not used
------------------------------------------------------------------------- */
const std::vector<Compute *> Modify::get_compute_by_style(const std::string &style) const
{
std::vector<Compute *> matches;
if (style.empty()) return matches;
for (int icompute = 0; icompute < ncompute; icompute++)
if (utils::strmatch(compute[icompute]->style,style)) matches.push_back(compute[icompute]);
return matches;
}
/* ----------------------------------------------------------------------
return vector with Computes
------------------------------------------------------------------------- */
const std::vector<Compute *> &Modify::get_compute_list()
{
compute_list = std::vector<Compute *>(compute, compute+ncompute);
return compute_list;
}
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
clear invoked flag of all computes clear invoked flag of all computes
called everywhere that computes are used, before computes are invoked called everywhere that computes are used, before computes are invoked

17
src/modify.h
View File

@ -17,6 +17,7 @@
#include "pointers.h" #include "pointers.h"
#include <map> #include <map>
#include <vector>
namespace LAMMPS_NS { namespace LAMMPS_NS {
@ -107,16 +108,28 @@ class Modify : protected Pointers {
void modify_fix(int, char **); void modify_fix(int, char **);
void delete_fix(const std::string &); void delete_fix(const std::string &);
void delete_fix(int); void delete_fix(int);
// deprecated API
int find_fix(const std::string &); int find_fix(const std::string &);
int find_fix_by_style(const char *); int find_fix_by_style(const char *);
// new API
Fix *get_fix_by_id(const std::string &) const;
const std::vector<Fix *> get_fix_by_style(const std::string &) const;
const std::vector<Fix *> &get_fix_list();
Compute *add_compute(int, char **, int trysuffix = 1); Compute *add_compute(int, char **, int trysuffix = 1);
Compute *add_compute(const std::string &, int trysuffix = 1); Compute *add_compute(const std::string &, int trysuffix = 1);
void modify_compute(int, char **); void modify_compute(int, char **);
void delete_compute(const std::string &); void delete_compute(const std::string &);
void delete_compute(int); void delete_compute(int);
// deprecated API
int find_compute(const std::string &); int find_compute(const std::string &);
int find_compute_by_style(const char *); int find_compute_by_style(const char *);
// new API
Compute *get_compute_by_id(const std::string &) const;
const std::vector<Compute *> get_compute_by_style(const std::string &) const;
const std::vector<Compute *> &get_compute_list();
void clearstep_compute(); void clearstep_compute();
void addstep_compute(bigint); void addstep_compute(bigint);
@ -165,6 +178,10 @@ class Modify : protected Pointers {
int index_permanent; // fix/compute index returned to library call int index_permanent; // fix/compute index returned to library call
// vectors to be used for new-API accessors as wrapper
std::vector<Fix *>fix_list;
std::vector<Compute *>compute_list;
void list_init(int, int &, int *&); void list_init(int, int &, int *&);
void list_init_end_of_step(int, int &, int *&); void list_init_end_of_step(int, int &, int *&);
void list_init_energy_couple(int &, int *&); void list_init_energy_couple(int &, int *&);

122
src/variable.cpp
View File

@ -963,8 +963,7 @@ double Variable::compute_equal(int ivar)
else if (style[ivar] == PYTHON) { else if (style[ivar] == PYTHON) {
int ifunc = python->find(data[ivar][0]); int ifunc = python->find(data[ivar][0]);
if (ifunc < 0) if (ifunc < 0)
print_var_error(FLERR,fmt::format("cannot find python function {}", print_var_error(FLERR,fmt::format("cannot find python function {}",data[ivar][0]),ivar);
data[ivar][0]),ivar);
python->invoke_function(ifunc,data[ivar][1]); python->invoke_function(ifunc,data[ivar][1]);
value = atof(data[ivar][1]); value = atof(data[ivar][1]);
} }
@ -1324,10 +1323,9 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
// compute // compute
// ---------------- // ----------------
if (strncmp(word,"c_",2) == 0 || strncmp(word,"C_",2) == 0) { if (utils::strmatch(word,"^[Cc]_")) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
print_var_error(FLERR,"Variable evaluation before " print_var_error(FLERR,"Variable evaluation before simulation box is defined",ivar);
"simulation box is defined",ivar);
// uppercase used to force access of // uppercase used to force access of
// global vector vs global scalar, and global array vs global vector // global vector vs global scalar, and global array vs global vector
@ -1335,14 +1333,10 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
int lowercase = 1; int lowercase = 1;
if (word[0] == 'C') lowercase = 0; if (word[0] == 'C') lowercase = 0;
int icompute = modify->find_compute(word+2); Compute *compute = modify->get_compute_by_id(word+2);
if (icompute < 0) { if (!compute)
std::string mesg = "Invalid compute ID '"; print_var_error(FLERR,fmt::format("Invalid compute ID '{}' in variable formula", word+2),ivar);
mesg += (word+2);
mesg += "' in variable formula";
print_var_error(FLERR,mesg,ivar);
}
Compute *compute = modify->compute[icompute];
// parse zero or one or two trailing brackets // parse zero or one or two trailing brackets
// point i beyond last bracket // point i beyond last bracket
@ -1371,8 +1365,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_scalar != update->ntimestep) if (compute->invoked_scalar != update->ntimestep)
print_var_error(FLERR,"Compute used in variable between " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"runs is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) { } else if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) {
compute->compute_scalar(); compute->compute_scalar();
compute->invoked_flag |= Compute::INVOKED_SCALAR; compute->invoked_flag |= Compute::INVOKED_SCALAR;
@ -1392,12 +1385,10 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
if (index1 > compute->size_vector && if (index1 > compute->size_vector &&
compute->size_vector_variable == 0) compute->size_vector_variable == 0)
print_var_error(FLERR,"Variable formula compute vector " print_var_error(FLERR,"Variable formula compute vector is accessed out-of-range",ivar,0);
"is accessed out-of-range",ivar,0);
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_vector != update->ntimestep) if (compute->invoked_vector != update->ntimestep)
print_var_error(FLERR,"Compute used in variable between runs " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
compute->compute_vector(); compute->compute_vector();
compute->invoked_flag |= Compute::INVOKED_VECTOR; compute->invoked_flag |= Compute::INVOKED_VECTOR;
@ -1419,15 +1410,12 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
if (index1 > compute->size_array_rows && if (index1 > compute->size_array_rows &&
compute->size_array_rows_variable == 0) compute->size_array_rows_variable == 0)
print_var_error(FLERR,"Variable formula compute array " print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0);
"is accessed out-of-range",ivar,0);
if (index2 > compute->size_array_cols) if (index2 > compute->size_array_cols)
print_var_error(FLERR,"Variable formula compute array " print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0);
"is accessed out-of-range",ivar,0);
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_array != update->ntimestep) if (compute->invoked_array != update->ntimestep)
print_var_error(FLERR,"Compute used in variable between runs " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) {
compute->compute_array(); compute->compute_array();
compute->invoked_flag |= Compute::INVOKED_ARRAY; compute->invoked_flag |= Compute::INVOKED_ARRAY;
@ -1448,18 +1436,14 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
} else if (nbracket == 0 && compute->vector_flag) { } else if (nbracket == 0 && compute->vector_flag) {
if (tree == nullptr) if (tree == nullptr)
print_var_error(FLERR,"Compute global vector in " print_var_error(FLERR,"Compute global vector in equal-style variable formula",ivar);
"equal-style variable formula",ivar);
if (treetype == ATOM) if (treetype == ATOM)
print_var_error(FLERR,"Compute global vector in " print_var_error(FLERR,"Compute global vector in atom-style variable formula",ivar);
"atom-style variable formula",ivar);
if (compute->size_vector == 0) if (compute->size_vector == 0)
print_var_error(FLERR,"Variable formula compute " print_var_error(FLERR,"Variable formula compute vector is zero length",ivar);
"vector is zero length",ivar);
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_vector != update->ntimestep) if (compute->invoked_vector != update->ntimestep)
print_var_error(FLERR,"Compute used in variable between " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"runs is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
compute->compute_vector(); compute->compute_vector();
compute->invoked_flag |= Compute::INVOKED_VECTOR; compute->invoked_flag |= Compute::INVOKED_VECTOR;
@ -1477,18 +1461,14 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
} else if (nbracket == 1 && compute->array_flag) { } else if (nbracket == 1 && compute->array_flag) {
if (tree == nullptr) if (tree == nullptr)
print_var_error(FLERR,"Compute global vector in " print_var_error(FLERR,"Compute global vector in equal-style variable formula",ivar);
"equal-style variable formula",ivar);
if (treetype == ATOM) if (treetype == ATOM)
print_var_error(FLERR,"Compute global vector in " print_var_error(FLERR,"Compute global vector in atom-style variable formula",ivar);
"atom-style variable formula",ivar);
if (compute->size_array_rows == 0) if (compute->size_array_rows == 0)
print_var_error(FLERR,"Variable formula compute array " print_var_error(FLERR,"Variable formula compute array is zero length",ivar);
"is zero length",ivar);
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_array != update->ntimestep) if (compute->invoked_array != update->ntimestep)
print_var_error(FLERR,"Compute used in variable between " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"runs is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) {
compute->compute_array(); compute->compute_array();
compute->invoked_flag |= Compute::INVOKED_ARRAY; compute->invoked_flag |= Compute::INVOKED_ARRAY;
@ -1508,8 +1488,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep) if (compute->invoked_peratom != update->ntimestep)
print_var_error(FLERR,"Compute used in variable " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"between runs is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom(); compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM; compute->invoked_flag |= Compute::INVOKED_PERATOM;
@ -1550,15 +1529,12 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
compute->size_peratom_cols == 0) { compute->size_peratom_cols == 0) {
if (tree == nullptr) if (tree == nullptr)
print_var_error(FLERR,"Per-atom compute in " print_var_error(FLERR,"Per-atom compute in equal-style variable formula",ivar);
"equal-style variable formula",ivar);
if (treetype == VECTOR) if (treetype == VECTOR)
print_var_error(FLERR,"Per-atom compute in " print_var_error(FLERR,"Per-atom compute in vector-style variable formula",ivar);
"vector-style variable formula",ivar);
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep) if (compute->invoked_peratom != update->ntimestep)
print_var_error(FLERR,"Compute used in variable " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"between runs is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom(); compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM; compute->invoked_flag |= Compute::INVOKED_PERATOM;
@ -1576,18 +1552,14 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
compute->size_peratom_cols > 0) { compute->size_peratom_cols > 0) {
if (tree == nullptr) if (tree == nullptr)
print_var_error(FLERR,"Per-atom compute in " print_var_error(FLERR,"Per-atom compute in equal-style variable formula",ivar);
"equal-style variable formula",ivar);
if (treetype == VECTOR) if (treetype == VECTOR)
print_var_error(FLERR,"Per-atom compute in " print_var_error(FLERR,"Per-atom compute in vector-style variable formula",ivar);
"vector-style variable formula",ivar);
if (index1 > compute->size_peratom_cols) if (index1 > compute->size_peratom_cols)
print_var_error(FLERR,"Variable formula compute array " print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0);
"is accessed out-of-range",ivar,0);
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_peratom != update->ntimestep) if (compute->invoked_peratom != update->ntimestep)
print_var_error(FLERR,"Compute used in variable " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"between runs is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) {
compute->compute_peratom(); compute->compute_peratom();
compute->invoked_flag |= Compute::INVOKED_PERATOM; compute->invoked_flag |= Compute::INVOKED_PERATOM;
@ -1602,17 +1574,15 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
} else if (nbracket == 1 && compute->local_flag) { } else if (nbracket == 1 && compute->local_flag) {
print_var_error(FLERR,"Cannot access local data via indexing",ivar); print_var_error(FLERR,"Cannot access local data via indexing",ivar);
} else print_var_error(FLERR, } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar);
"Mismatched compute in variable formula",ivar);
// ---------------- // ----------------
// fix // fix
// ---------------- // ----------------
} else if (strncmp(word,"f_",2) == 0 || strncmp(word,"F_",2) == 0) { } else if (utils::strmatch(word,"^[fF]_")) {
if (domain->box_exist == 0) if (domain->box_exist == 0)
print_var_error(FLERR,"Variable evaluation before " print_var_error(FLERR,"Variable evaluation before simulation box is defined",ivar);
"simulation box is defined",ivar);
// uppercase used to force access of // uppercase used to force access of
// global vector vs global scalar, and global array vs global vector // global vector vs global scalar, and global array vs global vector
@ -1620,11 +1590,10 @@ double Variable::evaluate(char *str, Tree **tree, int ivar)
int lowercase = 1; int lowercase = 1;
if (word[0] == 'F') lowercase = 0; if (word[0] == 'F') lowercase = 0;
int ifix = modify->find_fix(word+2); Fix *fix = modify->get_fix_by_id(word+2);
if (ifix < 0) if (!fix)
print_var_error(FLERR,fmt::format("Invalid fix ID '{}' in variable" print_var_error(FLERR,fmt::format("Invalid fix ID '{}' in variable formula",word+2),ivar);
" formula",word+2),ivar);
Fix *fix = modify->fix[ifix];
// parse zero or one or two trailing brackets // parse zero or one or two trailing brackets
// point i beyond last bracket // point i beyond last bracket
@ -4057,19 +4026,17 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
*ptr1 = '\0'; *ptr1 = '\0';
} else index = 0; } else index = 0;
int icompute = modify->find_compute(&args[0][2]); compute = modify->get_compute_by_id(&args[0][2]);
if (icompute < 0) { if (!compute) {
std::string mesg = "Invalid compute ID '"; std::string mesg = "Invalid compute ID '";
mesg += (args[0]+2); mesg += (args[0]+2);
mesg += "' in variable formula"; mesg += "' in variable formula";
print_var_error(FLERR,mesg,ivar); print_var_error(FLERR,mesg,ivar);
} }
compute = modify->compute[icompute];
if (index == 0 && compute->vector_flag) { if (index == 0 && compute->vector_flag) {
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_vector != update->ntimestep) if (compute->invoked_vector != update->ntimestep)
print_var_error(FLERR,"Compute used in variable between runs " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) {
compute->compute_vector(); compute->compute_vector();
compute->invoked_flag |= Compute::INVOKED_VECTOR; compute->invoked_flag |= Compute::INVOKED_VECTOR;
@ -4078,12 +4045,10 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
nstride = 1; nstride = 1;
} else if (index && compute->array_flag) { } else if (index && compute->array_flag) {
if (index > compute->size_array_cols) if (index > compute->size_array_cols)
print_var_error(FLERR,"Variable formula compute array " print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0);
"is accessed out-of-range",ivar,0);
if (update->whichflag == 0) { if (update->whichflag == 0) {
if (compute->invoked_array != update->ntimestep) if (compute->invoked_array != update->ntimestep)
print_var_error(FLERR,"Compute used in variable between runs " print_var_error(FLERR,"Compute used in variable between runs is not current",ivar);
"is not current",ivar);
} else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) {
compute->compute_array(); compute->compute_array();
compute->invoked_flag |= Compute::INVOKED_ARRAY; compute->invoked_flag |= Compute::INVOKED_ARRAY;
@ -4102,14 +4067,13 @@ int Variable::special_function(char *word, char *contents, Tree **tree,
*ptr1 = '\0'; *ptr1 = '\0';
} else index = 0; } else index = 0;
int ifix = modify->find_fix(&args[0][2]); fix = modify->get_fix_by_id(&args[0][2]);
if (ifix < 0) { if (!fix) {
std::string mesg = "Invalid fix ID '"; std::string mesg = "Invalid fix ID '";
mesg += (args[0]+2); mesg += (args[0]+2);
mesg += "' in variable formula"; mesg += "' in variable formula";
print_var_error(FLERR,mesg,ivar); print_var_error(FLERR,mesg,ivar);
} }
fix = modify->fix[ifix];
if (index == 0 && fix->vector_flag) { if (index == 0 && fix->vector_flag) {
if (update->whichflag > 0 && update->ntimestep % fix->global_freq) { if (update->whichflag > 0 && update->ntimestep % fix->global_freq) {
std::string mesg = "Fix with ID '"; std::string mesg = "Fix with ID '";

42
src/velocity.cpp
View File

@ -43,7 +43,7 @@ enum{NONE,CONSTANT,EQUAL,ATOM};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
Velocity::Velocity(LAMMPS *lmp) : Command(lmp) {} Velocity::Velocity(LAMMPS *lmp) : Command(lmp), rigid_fix(nullptr), temperature(nullptr) {}
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -91,7 +91,6 @@ void Velocity::command(int narg, char **arg)
bias_flag = 0; bias_flag = 0;
loop_flag = ALL; loop_flag = ALL;
scale_flag = 1; scale_flag = 1;
rfix = -1;
// read options from end of input line // read options from end of input line
// change defaults as options specify // change defaults as options specify
@ -108,8 +107,7 @@ void Velocity::command(int narg, char **arg)
// b/c methods invoked in the compute/fix perform forward/reverse comm // b/c methods invoked in the compute/fix perform forward/reverse comm
int initcomm = 0; int initcomm = 0;
if (style == ZERO && rfix >= 0 && if (style == ZERO && rigid_fix && utils::strmatch(rigid_fix->style,"^rigid.*/small.*")) initcomm = 1;
utils::strmatch(modify->fix[rfix]->style,"^rigid.*/small.*")) initcomm = 1;
if ((style == CREATE || style == SET) && temperature && if ((style == CREATE || style == SET) && temperature &&
strcmp(temperature->style,"temp/cs") == 0) initcomm = 1; strcmp(temperature->style,"temp/cs") == 0) initcomm = 1;
@ -154,7 +152,6 @@ void Velocity::init_external(const char *extgroup)
momentum_flag = 1; momentum_flag = 1;
rotation_flag = 0; rotation_flag = 0;
loop_flag = ALL; loop_flag = ALL;
rfix = -1;
scale_flag = 1; scale_flag = 1;
bias_flag = 0; bias_flag = 0;
} }
@ -692,21 +689,21 @@ void Velocity::ramp(int /*narg*/, char **arg)
void Velocity::zero(int /*narg*/, char **arg) void Velocity::zero(int /*narg*/, char **arg)
{ {
if (strcmp(arg[0],"linear") == 0) { if (strcmp(arg[0],"linear") == 0) {
if (rfix < 0) zero_momentum(); if (!rigid_fix) zero_momentum();
else if (utils::strmatch(modify->fix[rfix]->style,"^rigid/small")) { else if (utils::strmatch(rigid_fix->style,"^rigid/small")) {
modify->fix[rfix]->setup_pre_neighbor(); rigid_fix->setup_pre_neighbor();
modify->fix[rfix]->zero_momentum(); rigid_fix->zero_momentum();
} else if (utils::strmatch(modify->fix[rfix]->style,"^rigid")) { } else if (utils::strmatch(rigid_fix->style,"^rigid")) {
modify->fix[rfix]->zero_momentum(); rigid_fix->zero_momentum();
} else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID"); } else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID");
} else if (strcmp(arg[0],"angular") == 0) { } else if (strcmp(arg[0],"angular") == 0) {
if (rfix < 0) zero_rotation(); if (!rigid_fix) zero_rotation();
else if (utils::strmatch(modify->fix[rfix]->style,"^rigid/small")) { else if (utils::strmatch(rigid_fix->style,"^rigid/small")) {
modify->fix[rfix]->setup_pre_neighbor(); rigid_fix->setup_pre_neighbor();
modify->fix[rfix]->zero_rotation(); rigid_fix->zero_rotation();
} else if (utils::strmatch(modify->fix[rfix]->style,"^rigid")) { } else if (utils::strmatch(rigid_fix->style,"^rigid")) {
modify->fix[rfix]->zero_rotation(); rigid_fix->zero_rotation();
} else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID"); } else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID");
} else error->all(FLERR,"Illegal velocity command"); } else error->all(FLERR,"Illegal velocity command");
@ -843,11 +840,10 @@ void Velocity::options(int narg, char **arg)
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"temp") == 0) { } else if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
int icompute = modify->find_compute(arg[iarg+1]); temperature = modify->get_compute_by_id(arg[iarg+1]);
if (icompute < 0) error->all(FLERR,"Could not find velocity temperature ID"); if (!temperature) error->all(FLERR,"Could not find velocity temperature compute ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0) if (temperature->tempflag == 0)
error->all(FLERR,"Velocity temperature ID does not compute temperature"); error->all(FLERR,"Velocity temperature compute {} does not compute temperature", arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"bias") == 0) { } else if (strcmp(arg[iarg],"bias") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
@ -862,8 +858,8 @@ void Velocity::options(int narg, char **arg)
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"rigid") == 0) { } else if (strcmp(arg[iarg],"rigid") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");
rfix = modify->find_fix(arg[iarg+1]); rigid_fix = modify->get_fix_by_id(arg[iarg+1]);
if (rfix < 0) error->all(FLERR,"Fix ID for velocity does not exist"); if (!rigid_fix) error->all(FLERR,"Fix ID {} for velocity does not exist", arg[iarg+1]);
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"units") == 0) { } else if (strcmp(arg[iarg],"units") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command");

3
src/velocity.h
View File

@ -36,8 +36,9 @@ class Velocity : public Command {
int igroup, groupbit; int igroup, groupbit;
int style; int style;
int dist_flag, sum_flag, momentum_flag, rotation_flag; int dist_flag, sum_flag, momentum_flag, rotation_flag;
int bias_flag, loop_flag, scale_flag, rfix; int bias_flag, loop_flag, scale_flag;
double xscale, yscale, zscale; double xscale, yscale, zscale;
class Fix *rigid_fix;
class Compute *temperature; class Compute *temperature;
void set(int, char **); void set(int, char **);

3
src/verlet.cpp
View File

@ -68,8 +68,7 @@ void Verlet::init()
// detect if fix omp is present for clearing force arrays // detect if fix omp is present for clearing force arrays
int ifix = modify->find_fix("package_omp"); if (modify->get_fix_by_id("package_omp")) external_force_clear = 1;
if (ifix >= 0) external_force_clear = 1;
// set flags for arrays to clear in force_clear() // set flags for arrays to clear in force_clear()