From 5b40e4cb38bc5e008a5b8c427eaa52311cb64cfc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 8 Aug 2021 17:25:06 -0400 Subject: [PATCH] new accessor APIs for fixes and computes in Modify plus a few applications --- src/GRANULAR/fix_pour.cpp | 57 +++---- src/KOKKOS/fix_gravity_kokkos.cpp | 2 - src/OPENMP/fix_gravity_omp.cpp | 1 - src/atom.cpp | 11 +- src/balance.cpp | 7 +- src/change_box.cpp | 25 ++- src/compute_angmom_chunk.cpp | 75 +++++---- src/compute_centroid_stress_atom.cpp | 229 +++++++++++++-------------- src/compute_chunk_atom.cpp | 72 ++++----- src/domain.cpp | 30 ++-- src/dump.cpp | 4 +- src/finish.cpp | 5 +- src/fix_gravity.cpp | 1 - src/fix_gravity.h | 1 + src/modify.cpp | 78 +++++++++ src/modify.h | 17 ++ src/variable.cpp | 122 +++++--------- src/velocity.cpp | 42 +++-- src/velocity.h | 3 +- src/verlet.cpp | 3 +- 20 files changed, 397 insertions(+), 388 deletions(-) diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index a2f4d91840..063efda351 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -41,7 +41,6 @@ using namespace MathConst; enum{ATOM,MOLECULE}; enum{ONE,RANGE,POLY}; -enum{CONSTANT,EQUAL}; // same as FixGravity #define EPSILON 0.001 #define SMALL 1.0e-10 @@ -186,10 +185,12 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : // grav = gravity in distance/time^2 units // assume grav = -magnitude at this point, enforce in init() - int ifix = modify->find_fix_by_style("^gravity"); - if (ifix == -1) - error->all(FLERR,"No fix gravity defined for fix pour"); - grav = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v; + auto fixlist = modify->get_fix_by_style("^gravity"); + if (fixlist.size() != 1) + error->all(FLERR,"There must be exactly one fix gravity defined for fix pour"); + auto fixgrav = (FixGravity *)fixlist.front(); + + grav = -fixgrav->magnitude * force->ftm2v; // nfreq = timesteps between insertions // should be time for a particle to fall from top of insertion region @@ -208,9 +209,8 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : v_relative = vy - rate; delta = yhi - ylo; } - double t = - (-v_relative - sqrt(v_relative*v_relative - 2.0*grav*delta)) / grav; - nfreq = static_cast (t/update->dt + 0.5); + double t = (-v_relative - sqrt(v_relative*v_relative - 2.0*grav*delta)) / grav; + nfreq = static_cast(t/update->dt + 0.5); // 1st insertion on next timestep @@ -316,17 +316,16 @@ void FixPour::init() // for 3d must point in -z, for 2d must point in -y // else insertion cannot work - int ifix = modify->find_fix_by_style("^gravity"); - if (ifix == -1) - error->all(FLERR,"No fix gravity defined for fix pour"); - - int varflag = ((FixGravity *) modify->fix[ifix])->varflag; - if (varflag != CONSTANT) + auto fixlist = modify->get_fix_by_style("^gravity"); + if (fixlist.size() != 1) + error->all(FLERR,"There must be exactly one fix gravity defined for fix pour"); + auto fixgrav = (FixGravity *)fixlist.front(); + if (fixgrav->varflag != FixGravity::CONSTANT) error->all(FLERR,"Fix gravity for fix pour must be constant"); - double xgrav = ((FixGravity *) modify->fix[ifix])->xgrav; - double ygrav = ((FixGravity *) modify->fix[ifix])->ygrav; - double zgrav = ((FixGravity *) modify->fix[ifix])->zgrav; + double xgrav = fixgrav->xgrav; + double ygrav = fixgrav->ygrav; + double zgrav = fixgrav->zgrav; if (domain->dimension == 3) { 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"); } - double gnew = - ((FixGravity *) modify->fix[ifix])->magnitude * force->ftm2v; - if (gnew != grav) - error->all(FLERR,"Gravity changed since fix pour was created"); + double gnew = -fixgrav->magnitude * force->ftm2v; + if (gnew != grav) error->all(FLERR,"Gravity changed since fix pour was created"); // if rigidflag defined, check for rigid/small fix // its molecule template must be same as this one - fixrigid = nullptr; + fixrigid = modify->get_fix_by_id(idrigid); if (rigidflag) { - int ifix = modify->find_fix(idrigid); - if (ifix < 0) error->all(FLERR,"Fix pour rigid fix does not exist"); - fixrigid = modify->fix[ifix]; + if (!fixrigid) error->all(FLERR,"Fix pour rigid fix does not exist"); int tmp; if (onemols != (Molecule **) fixrigid->extract("onemol",tmp)) - error->all(FLERR, - "Fix pour and fix rigid/small not using " - "same molecule template ID"); + error->all(FLERR,"Fix pour and fix rigid/small not using same molecule template ID"); } // if shakeflag defined, check for SHAKE fix // its molecule template must be same as this one - fixshake = nullptr; + fixshake = modify->get_fix_by_id(idshake); if (shakeflag) { - int ifix = modify->find_fix(idshake); - if (ifix < 0) error->all(FLERR,"Fix pour shake fix does not exist"); - fixshake = modify->fix[ifix]; + if (!fixshake) error->all(FLERR,"Fix pour shake fix does not exist"); int tmp; if (onemols != (Molecule **) fixshake->extract("onemol",tmp)) - error->all(FLERR,"Fix pour and fix shake not using " - "same molecule template ID"); + error->all(FLERR,"Fix pour and fix shake not using same molecule template ID"); } } diff --git a/src/KOKKOS/fix_gravity_kokkos.cpp b/src/KOKKOS/fix_gravity_kokkos.cpp index ec6afbc133..63ac87a786 100644 --- a/src/KOKKOS/fix_gravity_kokkos.cpp +++ b/src/KOKKOS/fix_gravity_kokkos.cpp @@ -24,8 +24,6 @@ using namespace LAMMPS_NS; -enum{CONSTANT,EQUAL}; - /* ---------------------------------------------------------------------- */ template diff --git a/src/OPENMP/fix_gravity_omp.cpp b/src/OPENMP/fix_gravity_omp.cpp index 7354ec2aa4..72d2fdc59b 100644 --- a/src/OPENMP/fix_gravity_omp.cpp +++ b/src/OPENMP/fix_gravity_omp.cpp @@ -30,7 +30,6 @@ using namespace LAMMPS_NS; using namespace FixConst; enum{CHUTE,SPHERICAL,GRADIENT,VECTOR}; -enum{CONSTANT,EQUAL}; /* ---------------------------------------------------------------------- */ diff --git a/src/atom.cpp b/src/atom.cpp index 796c0ba156..68eedd6eaa 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -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) { - for (int m = 0; m < modify->nfix; m++) { - Fix *fix = modify->fix[m]; + for (const auto &fix : modify->get_fix_list()) { if (fix->create_attribute) for (int i = nprev; i < nnew; i++) fix->set_arrays(i); @@ -2238,15 +2237,13 @@ void Atom::setup_sort_bins() #ifdef LMP_GPU if (userbinsize == 0.0) { - int ifix = modify->find_fix("package_gpu"); - if (ifix >= 0) { + FixGPU *fix = (FixGPU *)modify->get_fix_by_id("package_gpu"); + if (fix) { const double subx = domain->subhi[0] - domain->sublo[0]; const double suby = domain->subhi[1] - domain->sublo[1]; const double subz = domain->subhi[2] - domain->sublo[2]; - FixGPU *fix = static_cast(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; nbinx = static_cast (ceil(subx * bininv)); diff --git a/src/balance.cpp b/src/balance.cpp index bd3ba007ef..0543de9971 100644 --- a/src/balance.cpp +++ b/src/balance.cpp @@ -496,11 +496,8 @@ void Balance::weight_storage(char *prefix) if (prefix) cmd = prefix; cmd += "IMBALANCE_WEIGHTS"; - int ifix = modify->find_fix(cmd); - if (ifix < 1) { - cmd += " all STORE peratom 0 1"; - fixstore = (FixStore *) modify->add_fix(cmd); - } else fixstore = (FixStore *) modify->fix[ifix]; + fixstore = (FixStore *) modify->get_fix_by_id(cmd); + if (!fixstore) fixstore = (FixStore *) modify->add_fix(cmd + " all STORE peratom 0 1"); // do not carry weights with atoms during normal atom migration diff --git a/src/change_box.cpp b/src/change_box.cpp index 91d65bfcad..bbac78ab3d 100644 --- a/src/change_box.cpp +++ b/src/change_box.cpp @@ -291,16 +291,12 @@ void ChangeBox::command(int narg, char **arg) } else if (ops[m].style == ORTHO) { if (domain->xy != 0.0 || domain->yz != 0.0 || domain->xz != 0.0) - error->all(FLERR, - "Cannot change box to orthogonal when tilt is non-zero"); + error->all(FLERR,"Cannot change box to orthogonal when tilt is non-zero"); if (output->ndump) - error->all(FLERR, - "Cannot change box ortho/triclinic with dumps defined"); - for (i = 0; i < modify->nfix; i++) - 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 dumps defined"); + for (const auto &fix : modify->get_fix_list()) + if (fix->no_change_box) + error->all(FLERR,"Cannot change box ortho/triclinic with certain fixes defined"); domain->triclinic = 0; domain->set_initial_box(); domain->set_global_box(); @@ -309,13 +305,10 @@ void ChangeBox::command(int narg, char **arg) } else if (ops[m].style == TRICLINIC) { if (output->ndump) - error->all(FLERR, - "Cannot change box ortho/triclinic with dumps defined"); - for (i = 0; i < modify->nfix; i++) - 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 dumps defined"); + for (const auto &fix : modify->get_fix_list()) + if (fix->no_change_box) + error->all(FLERR,"Cannot change box ortho/triclinic with certain fixes defined"); domain->triclinic = 1; domain->set_lamda_box(); domain->set_initial_box(); diff --git a/src/compute_angmom_chunk.cpp b/src/compute_angmom_chunk.cpp index 787dc43fc6..1351d379b7 100644 --- a/src/compute_angmom_chunk.cpp +++ b/src/compute_angmom_chunk.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -29,10 +28,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeAngmomChunk::ComputeAngmomChunk(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - idchunk(nullptr), massproc(nullptr), masstotal(nullptr), com(nullptr), comall(nullptr), angmom(nullptr), angmomall(nullptr) + Compute(lmp, narg, arg), idchunk(nullptr), massproc(nullptr), masstotal(nullptr), com(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"); array_flag = 1; size_array_cols = 3; @@ -57,7 +56,7 @@ ComputeAngmomChunk::ComputeAngmomChunk(LAMMPS *lmp, int narg, char **arg) : ComputeAngmomChunk::~ComputeAngmomChunk() { - delete [] idchunk; + delete[] idchunk; memory->destroy(massproc); memory->destroy(masstotal); memory->destroy(com); @@ -70,21 +69,18 @@ ComputeAngmomChunk::~ComputeAngmomChunk() void ComputeAngmomChunk::init() { - int icompute = modify->find_compute(idchunk); - if (icompute < 0) - 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) - error->all(FLERR,"Compute angmom/chunk does not use chunk/atom compute"); + cchunk = (ComputeChunkAtom *) modify->get_compute_by_id(idchunk); + if (!cchunk) error->all(FLERR, "Chunk/atom compute does not exist for compute angmom/chunk"); + if (strcmp(cchunk->style, "chunk/atom") != 0) + error->all(FLERR, "Compute angmom/chunk does not use chunk/atom compute"); } /* ---------------------------------------------------------------------- */ void ComputeAngmomChunk::compute_array() { - int i,index; - double dx,dy,dz,massone; + int i, index; + double dx, dy, dz, massone; double unwrap[3]; invoked_array = update->ntimestep; @@ -120,19 +116,21 @@ void ComputeAngmomChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - domain->unmap(x[i],image[i],unwrap); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + domain->unmap(x[i], image[i], unwrap); massproc[index] += massone; com[index][0] += unwrap[0] * massone; com[index][1] += unwrap[1] * massone; com[index][2] += unwrap[2] * massone; } - MPI_Allreduce(massproc,masstotal,nchunk,MPI_DOUBLE,MPI_SUM,world); - MPI_Allreduce(&com[0][0],&comall[0][0],3*nchunk,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(massproc, masstotal, nchunk, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(&com[0][0], &comall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); for (i = 0; i < nchunk; i++) { if (masstotal[i] > 0.0) { @@ -148,21 +146,22 @@ void ComputeAngmomChunk::compute_array() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - index = ichunk[i]-1; + index = ichunk[i] - 1; if (index < 0) continue; - domain->unmap(x[i],image[i],unwrap); + domain->unmap(x[i], image[i], unwrap); dx = unwrap[0] - comall[index][0]; dy = unwrap[1] - comall[index][1]; dz = unwrap[2] - comall[index][2]; - if (rmass) massone = rmass[i]; - else massone = mass[type[i]]; - 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][2] += massone * (dx*v[i][1] - dy*v[i][0]); + if (rmass) + massone = rmass[i]; + else + massone = mass[type[i]]; + 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][2] += massone * (dx * v[i][1] - dy * v[i][0]); } - MPI_Allreduce(&angmom[0][0],&angmomall[0][0],3*nchunk, - MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&angmom[0][0], &angmomall[0][0], 3 * nchunk, MPI_DOUBLE, MPI_SUM, world); } /* ---------------------------------------------------------------------- @@ -209,7 +208,7 @@ int ComputeAngmomChunk::lock_length() void ComputeAngmomChunk::lock(Fix *fixptr, bigint startstep, bigint stopstep) { - cchunk->lock(fixptr,startstep,stopstep); + cchunk->lock(fixptr, startstep, stopstep); } /* ---------------------------------------------------------------------- @@ -234,12 +233,12 @@ void ComputeAngmomChunk::allocate() memory->destroy(angmom); memory->destroy(angmomall); maxchunk = nchunk; - memory->create(massproc,maxchunk,"angmom/chunk:massproc"); - memory->create(masstotal,maxchunk,"angmom/chunk:masstotal"); - memory->create(com,maxchunk,3,"angmom/chunk:com"); - memory->create(comall,maxchunk,3,"angmom/chunk:comall"); - memory->create(angmom,maxchunk,3,"angmom/chunk:angmom"); - memory->create(angmomall,maxchunk,3,"angmom/chunk:angmomall"); + memory->create(massproc, maxchunk, "angmom/chunk:massproc"); + memory->create(masstotal, maxchunk, "angmom/chunk:masstotal"); + memory->create(com, maxchunk, 3, "angmom/chunk:com"); + memory->create(comall, maxchunk, 3, "angmom/chunk:comall"); + memory->create(angmom, maxchunk, 3, "angmom/chunk:angmom"); + memory->create(angmomall, maxchunk, 3, "angmom/chunk:angmomall"); array = angmomall; } @@ -250,7 +249,7 @@ void ComputeAngmomChunk::allocate() double ComputeAngmomChunk::memory_usage() { double bytes = (bigint) maxchunk * 2 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); - bytes += (double) maxchunk * 2*3 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); + bytes += (double) maxchunk * 2 * 3 * sizeof(double); return bytes; } diff --git a/src/compute_centroid_stress_atom.cpp b/src/compute_centroid_stress_atom.cpp index c35b6a63e7..a050c8bb6a 100644 --- a/src/compute_centroid_stress_atom.cpp +++ b/src/compute_centroid_stress_atom.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,33 +12,34 @@ ------------------------------------------------------------------------- */ #include "compute_centroid_stress_atom.h" -#include -#include "atom.h" -#include "update.h" -#include "comm.h" -#include "force.h" -#include "pair.h" -#include "bond.h" + #include "angle.h" +#include "atom.h" +#include "bond.h" +#include "comm.h" #include "dihedral.h" +#include "error.h" +#include "fix.h" +#include "force.h" #include "improper.h" #include "kspace.h" -#include "modify.h" -#include "fix.h" #include "memory.h" -#include "error.h" +#include "modify.h" +#include "pair.h" +#include "update.h" + +#include using namespace LAMMPS_NS; -enum{NOBIAS,BIAS}; +enum { NOBIAS, BIAS }; /* ---------------------------------------------------------------------- */ ComputeCentroidStressAtom::ComputeCentroidStressAtom(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - id_temp(nullptr), stress(nullptr) + Compute(lmp, narg, arg), 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"); peratom_flag = 1; size_peratom_cols = 9; @@ -50,17 +50,16 @@ ComputeCentroidStressAtom::ComputeCentroidStressAtom(LAMMPS *lmp, int narg, char // store temperature ID used by stress 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 { id_temp = utils::strdup(arg[3]); - int icompute = modify->find_compute(id_temp); - if (icompute < 0) - error->all(FLERR,"Could not find compute centroid/stress/atom temperature ID"); - if (modify->compute[icompute]->tempflag == 0) - error->all(FLERR, - "Compute centroid/stress/atom temperature ID does not " - "compute temperature"); + auto compute = modify->get_compute_by_id(id_temp); + if (!compute) + error->all(FLERR, "Could not find compute centroid/stress/atom temperature ID {}", id_temp); + if (compute->tempflag == 0) + error->all(FLERR, "Compute centroid/stress/atom temperature ID does not compute temperature"); } // process optional args @@ -79,19 +78,28 @@ ComputeCentroidStressAtom::ComputeCentroidStressAtom(LAMMPS *lmp, int narg, char fixflag = 0; int iarg = 4; while (iarg < narg) { - if (strcmp(arg[iarg],"ke") == 0) keflag = 1; - else if (strcmp(arg[iarg],"pair") == 0) pairflag = 1; - else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1; - else if (strcmp(arg[iarg],"angle") == 0) 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) { + if (strcmp(arg[iarg], "ke") == 0) + keflag = 1; + else if (strcmp(arg[iarg], "pair") == 0) + pairflag = 1; + else if (strcmp(arg[iarg], "bond") == 0) + bondflag = 1; + else if (strcmp(arg[iarg], "angle") == 0) + 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) { pairflag = 1; bondflag = angleflag = dihedralflag = improperflag = 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++; } } @@ -103,7 +111,7 @@ ComputeCentroidStressAtom::ComputeCentroidStressAtom(LAMMPS *lmp, int narg, char ComputeCentroidStressAtom::~ComputeCentroidStressAtom() { - delete [] id_temp; + delete[] id_temp; memory->destroy(stress); } @@ -115,13 +123,15 @@ void ComputeCentroidStressAtom::init() // fixes could have changed or compute_modify could have changed it if (id_temp) { - int icompute = modify->find_compute(id_temp); - if (icompute < 0) - error->all(FLERR,"Could not find compute centroid/stress/atom temperature ID"); - temperature = modify->compute[icompute]; - if (temperature->tempbias) biasflag = BIAS; - else biasflag = NOBIAS; - } else biasflag = NOBIAS; + temperature = modify->get_compute_by_id(id_temp); + if (!temperature) + error->all(FLERR, "Could not find compute centroid/stress/atom temperature ID {}",id_temp); + if (temperature->tempbias) + biasflag = BIAS; + else + biasflag = NOBIAS; + } else + biasflag = NOBIAS; // check if force components and fixes support centroid atom stress // all bond styles support it as CENTROID_SAME @@ -158,12 +168,12 @@ void ComputeCentroidStressAtom::init() void ComputeCentroidStressAtom::compute_peratom() { - int i,j; + int i, j; double onemass; invoked_peratom = update->ntimestep; if (update->vflag_atom != invoked_peratom) - error->all(FLERR,"Per-atom virial was not tallied on needed timestep"); + error->all(FLERR, "Per-atom virial was not tallied on needed timestep"); // grow local stress array if necessary // needs to be atom->nmax in length @@ -171,7 +181,7 @@ void ComputeCentroidStressAtom::compute_peratom() if (atom->nmax > nmax) { memory->destroy(stress); nmax = atom->nmax; - memory->create(stress,nmax,9,"centroid/stress/atom:stress"); + memory->create(stress, nmax, 9, "centroid/stress/atom:stress"); array_atom = stress; } @@ -194,8 +204,7 @@ void ComputeCentroidStressAtom::compute_peratom() // clear local stress array for (i = 0; i < ntotal; i++) - for (j = 0; j < 9; j++) - stress[i][j] = 0.0; + for (j = 0; j < 9; j++) stress[i][j] = 0.0; // 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) { double **cvatom = force->pair->cvatom; for (i = 0; i < npair; i++) - for (j = 0; j < 9; j++) - stress[i][j] += cvatom[i][j]; + for (j = 0; j < 9; j++) stress[i][j] += cvatom[i][j]; } else { double **vatom = force->pair->vatom; for (i = 0; i < npair; i++) { - for (j = 0; j < 6; j++) - stress[i][j] += vatom[i][j]; - for (j = 6; j < 9; j++) - stress[i][j] += vatom[i][j-3]; + for (j = 0; j < 6; j++) stress[i][j] += vatom[i][j]; + 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) { double **vatom = force->bond->vatom; for (i = 0; i < nbond; i++) { - for (j = 0; j < 6; j++) - stress[i][j] += vatom[i][j]; - for (j = 6; j < 9; j++) - stress[i][j] += vatom[i][j-3]; + for (j = 0; j < 6; j++) stress[i][j] += vatom[i][j]; + for (j = 6; j < 9; j++) stress[i][j] += vatom[i][j - 3]; } } if (angleflag && force->angle) { double **cvatom = force->angle->cvatom; for (i = 0; i < nbond; i++) - for (j = 0; j < 9; j++) - stress[i][j] += cvatom[i][j]; + for (j = 0; j < 9; j++) stress[i][j] += cvatom[i][j]; } if (dihedralflag && force->dihedral) { double **cvatom = force->dihedral->cvatom; for (i = 0; i < nbond; i++) - for (j = 0; j < 9; j++) - stress[i][j] += cvatom[i][j]; + for (j = 0; j < 9; j++) stress[i][j] += cvatom[i][j]; } if (improperflag && force->improper) { double **cvatom = force->improper->cvatom; for (i = 0; i < nbond; i++) - for (j = 0; j < 9; j++) - stress[i][j] += cvatom[i][j]; + for (j = 0; j < 9; j++) stress[i][j] += cvatom[i][j]; } if (kspaceflag && force->kspace && force->kspace->compute_flag) { double **vatom = force->kspace->vatom; for (i = 0; i < nkspace; i++) { - for (j = 0; j < 6; j++) - stress[i][j] += vatom[i][j]; - for (j = 6; j < 9; j++) - stress[i][j] += vatom[i][j-3]; + for (j = 0; j < 6; j++) stress[i][j] += vatom[i][j]; + 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; if (vatom) for (i = 0; i < nlocal; i++) { - for (j = 0; j < 6; j++) - stress[i][j] += vatom[i][j]; - for (j = 6; j < 9; j++) - stress[i][j] += vatom[i][j-3]; + for (j = 0; j < 6; j++) stress[i][j] += vatom[i][j]; + for (j = 6; j < 9; j++) stress[i][j] += vatom[i][j - 3]; } } } // communicate ghost virials between neighbor procs - if (force->newton || - (force->kspace && force->kspace->tip4pflag && force->kspace->compute_flag)) + if (force->newton || (force->kspace && force->kspace->tip4pflag && force->kspace->compute_flag)) comm->reverse_comm_compute(this); // zero virial of atoms not in group @@ -327,30 +323,30 @@ void ComputeCentroidStressAtom::compute_peratom() for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { onemass = mvv2e * rmass[i]; - stress[i][0] += onemass*v[i][0]*v[i][0]; - stress[i][1] += onemass*v[i][1]*v[i][1]; - stress[i][2] += onemass*v[i][2]*v[i][2]; - stress[i][3] += onemass*v[i][0]*v[i][1]; - stress[i][4] += onemass*v[i][0]*v[i][2]; - stress[i][5] += onemass*v[i][1]*v[i][2]; - stress[i][6] += onemass*v[i][1]*v[i][0]; - stress[i][7] += onemass*v[i][2]*v[i][0]; - stress[i][8] += onemass*v[i][2]*v[i][1]; + stress[i][0] += onemass * v[i][0] * v[i][0]; + stress[i][1] += onemass * v[i][1] * v[i][1]; + stress[i][2] += onemass * v[i][2] * v[i][2]; + stress[i][3] += onemass * v[i][0] * v[i][1]; + stress[i][4] += onemass * v[i][0] * v[i][2]; + stress[i][5] += onemass * v[i][1] * v[i][2]; + stress[i][6] += onemass * v[i][1] * v[i][0]; + stress[i][7] += onemass * v[i][2] * v[i][0]; + stress[i][8] += onemass * v[i][2] * v[i][1]; } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { onemass = mvv2e * mass[type[i]]; - stress[i][0] += onemass*v[i][0]*v[i][0]; - stress[i][1] += onemass*v[i][1]*v[i][1]; - stress[i][2] += onemass*v[i][2]*v[i][2]; - stress[i][3] += onemass*v[i][0]*v[i][1]; - stress[i][4] += onemass*v[i][0]*v[i][2]; - stress[i][5] += onemass*v[i][1]*v[i][2]; - stress[i][6] += onemass*v[i][1]*v[i][0]; - stress[i][7] += onemass*v[i][2]*v[i][0]; - stress[i][8] += onemass*v[i][2]*v[i][1]; + stress[i][0] += onemass * v[i][0] * v[i][0]; + stress[i][1] += onemass * v[i][1] * v[i][1]; + stress[i][2] += onemass * v[i][2] * v[i][2]; + stress[i][3] += onemass * v[i][0] * v[i][1]; + stress[i][4] += onemass * v[i][0] * v[i][2]; + stress[i][5] += onemass * v[i][1] * v[i][2]; + stress[i][6] += onemass * v[i][1] * v[i][0]; + stress[i][7] += onemass * v[i][2] * v[i][0]; + stress[i][8] += onemass * v[i][2] * v[i][1]; } } @@ -359,41 +355,40 @@ void ComputeCentroidStressAtom::compute_peratom() // invoke temperature if it hasn't been already // this insures bias factor is pre-computed - if (keflag && temperature->invoked_scalar != update->ntimestep) - temperature->compute_scalar(); + if (keflag && temperature->invoked_scalar != update->ntimestep) temperature->compute_scalar(); if (rmass) { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); + temperature->remove_bias(i, v[i]); onemass = mvv2e * rmass[i]; - stress[i][0] += onemass*v[i][0]*v[i][0]; - stress[i][1] += onemass*v[i][1]*v[i][1]; - stress[i][2] += onemass*v[i][2]*v[i][2]; - stress[i][3] += onemass*v[i][0]*v[i][1]; - stress[i][4] += onemass*v[i][0]*v[i][2]; - stress[i][5] += onemass*v[i][1]*v[i][2]; - stress[i][6] += onemass*v[i][1]*v[i][0]; - stress[i][7] += onemass*v[i][2]*v[i][0]; - stress[i][8] += onemass*v[i][2]*v[i][1]; - temperature->restore_bias(i,v[i]); + stress[i][0] += onemass * v[i][0] * v[i][0]; + stress[i][1] += onemass * v[i][1] * v[i][1]; + stress[i][2] += onemass * v[i][2] * v[i][2]; + stress[i][3] += onemass * v[i][0] * v[i][1]; + stress[i][4] += onemass * v[i][0] * v[i][2]; + stress[i][5] += onemass * v[i][1] * v[i][2]; + stress[i][6] += onemass * v[i][1] * v[i][0]; + stress[i][7] += onemass * v[i][2] * v[i][0]; + stress[i][8] += onemass * v[i][2] * v[i][1]; + temperature->restore_bias(i, v[i]); } } else { for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - temperature->remove_bias(i,v[i]); + temperature->remove_bias(i, v[i]); onemass = mvv2e * mass[type[i]]; - stress[i][0] += onemass*v[i][0]*v[i][0]; - stress[i][1] += onemass*v[i][1]*v[i][1]; - stress[i][2] += onemass*v[i][2]*v[i][2]; - stress[i][3] += onemass*v[i][0]*v[i][1]; - stress[i][4] += onemass*v[i][0]*v[i][2]; - stress[i][5] += onemass*v[i][1]*v[i][2]; - stress[i][6] += onemass*v[i][1]*v[i][0]; - stress[i][7] += onemass*v[i][2]*v[i][0]; - stress[i][8] += onemass*v[i][2]*v[i][1]; - temperature->restore_bias(i,v[i]); + stress[i][0] += onemass * v[i][0] * v[i][0]; + stress[i][1] += onemass * v[i][1] * v[i][1]; + stress[i][2] += onemass * v[i][2] * v[i][2]; + stress[i][3] += onemass * v[i][0] * v[i][1]; + stress[i][4] += onemass * v[i][0] * v[i][2]; + stress[i][5] += onemass * v[i][1] * v[i][2]; + stress[i][6] += onemass * v[i][1] * v[i][0]; + stress[i][7] += onemass * v[i][2] * v[i][0]; + stress[i][8] += onemass * v[i][2] * v[i][1]; + temperature->restore_bias(i, v[i]); } } } @@ -420,7 +415,7 @@ void ComputeCentroidStressAtom::compute_peratom() int ComputeCentroidStressAtom::pack_reverse_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; m = 0; last = first + n; @@ -442,7 +437,7 @@ int ComputeCentroidStressAtom::pack_reverse_comm(int n, int first, double *buf) void ComputeCentroidStressAtom::unpack_reverse_comm(int n, int *list, double *buf) { - int i,j,m; + int i, j, m; m = 0; for (i = 0; i < n; i++) { @@ -465,6 +460,6 @@ void ComputeCentroidStressAtom::unpack_reverse_comm(int n, int *list, double *bu double ComputeCentroidStressAtom::memory_usage() { - double bytes = (double)nmax*9 * sizeof(double); + double bytes = (double) nmax * 9 * sizeof(double); return bytes; } diff --git a/src/compute_chunk_atom.cpp b/src/compute_chunk_atom.cpp index fc4f17e53a..daa8fc64e7 100644 --- a/src/compute_chunk_atom.cpp +++ b/src/compute_chunk_atom.cpp @@ -316,40 +316,30 @@ ComputeChunkAtom::ComputeChunkAtom(LAMMPS *lmp, int narg, char **arg) : } if (which == ArgInfo::COMPUTE) { - int icompute = modify->find_compute(cfvid); - if (icompute < 0) - error->all(FLERR,"Compute ID for compute chunk /atom does not exist"); - if (modify->compute[icompute]->peratom_flag == 0) - error->all(FLERR, - "Compute chunk/atom compute does not calculate " - "per-atom values"); - if (argindex == 0 && - modify->compute[icompute]->size_peratom_cols != 0) - error->all(FLERR,"Compute chunk/atom compute does not " - "calculate a per-atom vector"); - 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"); + cchunk = modify->get_compute_by_id(cfvid); + if (!cchunk) + error->all(FLERR,"Compute ID {} for compute chunk /atom does not exist",cfvid); + if (cchunk->peratom_flag == 0) + error->all(FLERR,"Compute chunk/atom compute does not calculate per-atom values"); + if ((argindex == 0) && (cchunk->size_peratom_cols != 0)) + error->all(FLERR,"Compute chunk/atom compute does not calculate a per-atom vector"); + if (argindex && (cchunk->size_peratom_cols == 0)) + error->all(FLERR,"Compute chunk/atom compute does not calculate a per-atom array"); + if (argindex && argindex > cchunk->size_peratom_cols) + error->all(FLERR,"Compute chunk/atom compute array is accessed out-of-range"); } if (which == ArgInfo::FIX) { - int ifix = modify->find_fix(cfvid); - if (ifix < 0) - error->all(FLERR,"Fix ID for compute chunk/atom does not exist"); - if (modify->fix[ifix]->peratom_flag == 0) - error->all(FLERR,"Compute chunk/atom fix does not calculate " - "per-atom values"); - if (argindex == 0 && modify->fix[ifix]->size_peratom_cols != 0) - error->all(FLERR, - "Compute chunk/atom fix does not calculate a per-atom vector"); - if (argindex && modify->fix[ifix]->size_peratom_cols == 0) - error->all(FLERR, - "Compute chunk/atom fix does not calculate a per-atom array"); - if (argindex && argindex > modify->fix[ifix]->size_peratom_cols) + fchunk = modify->get_fix_by_id(cfvid); + if (!fchunk) + error->all(FLERR,"Fix ID {} for compute chunk/atom does not exist",cfvid); + if (fchunk->peratom_flag == 0) + error->all(FLERR,"Compute chunk/atom fix does not calculate per-atom values"); + if (argindex == 0 && fchunk->size_peratom_cols != 0) + error->all(FLERR,"Compute chunk/atom fix does not calculate a per-atom vector"); + if (argindex && fchunk->size_peratom_cols == 0) + error->all(FLERR,"Compute chunk/atom fix does not calculate a per-atom array"); + if (argindex && argindex > fchunk->size_peratom_cols) 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) error->all(FLERR,"Variable name for compute chunk/atom does not exist"); if (input->variable->atomstyle(ivariable) == 0) - error->all(FLERR,"Compute chunk/atom variable is not " - "atom-style variable"); + error->all(FLERR,"Compute chunk/atom variable is not atom-style variable"); } // setup scaling if (binflag) { if (domain->triclinic == 1 && scaleflag != REDUCED) - error->all(FLERR,"Compute chunk/atom for triclinic boxes " - "requires units reduced"); + error->all(FLERR,"Compute chunk/atom for triclinic boxes requires units reduced"); } if (scaleflag == LATTICE) { @@ -501,15 +489,13 @@ void ComputeChunkAtom::init() // set compute,fix,variable if (which == ArgInfo::COMPUTE) { - int icompute = modify->find_compute(cfvid); - if (icompute < 0) - error->all(FLERR,"Compute ID for compute chunk/atom does not exist"); - cchunk = modify->compute[icompute]; + cchunk = modify->get_compute_by_id(cfvid); + if (!cchunk) + error->all(FLERR,"Compute ID {} for compute chunk/atom does not exist",cfvid); } else if (which == ArgInfo::FIX) { - int ifix = modify->find_fix(cfvid); - if (ifix < 0) - error->all(FLERR,"Fix ID for compute chunk/atom does not exist"); - fchunk = modify->fix[ifix]; + fchunk = modify->get_fix_by_id(cfvid); + if (!fchunk) + error->all(FLERR,"Fix ID {} for compute chunk/atom does not exist",cfvid); } else if (which == ArgInfo::VARIABLE) { int ivariable = input->variable->find(cfvid); if (ivariable < 0) diff --git a/src/domain.cpp b/src/domain.cpp index daaf41338f..4759490ae8 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -140,19 +140,19 @@ void Domain::init() 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; - Fix **fixes = modify->fix; + const auto &fixes = modify->get_fix_list(); if (nonperiodic == 2) box_change_size = 1; - for (int i = 0; i < modify->nfix; i++) { - if (fixes[i]->box_change & Fix::BOX_CHANGE_SIZE) box_change_size = 1; - if (fixes[i]->box_change & Fix::BOX_CHANGE_SHAPE) box_change_shape = 1; - if (fixes[i]->box_change & Fix::BOX_CHANGE_DOMAIN) box_change_domain = 1; - if (fixes[i]->box_change & Fix::BOX_CHANGE_X) box_change_x++; - if (fixes[i]->box_change & Fix::BOX_CHANGE_Y) box_change_y++; - if (fixes[i]->box_change & Fix::BOX_CHANGE_Z) box_change_z++; - if (fixes[i]->box_change & Fix::BOX_CHANGE_YZ) box_change_yz++; - if (fixes[i]->box_change & Fix::BOX_CHANGE_XZ) box_change_xz++; - if (fixes[i]->box_change & Fix::BOX_CHANGE_XY) box_change_xy++; + for (const auto &fix : fixes) { + if (fix->box_change & Fix::BOX_CHANGE_SIZE) box_change_size = 1; + if (fix->box_change & Fix::BOX_CHANGE_SHAPE) box_change_shape = 1; + if (fix->box_change & Fix::BOX_CHANGE_DOMAIN) box_change_domain = 1; + if (fix->box_change & Fix::BOX_CHANGE_X) box_change_x++; + if (fix->box_change & Fix::BOX_CHANGE_Y) box_change_y++; + if (fix->box_change & Fix::BOX_CHANGE_Z) box_change_z++; + if (fix->box_change & Fix::BOX_CHANGE_YZ) box_change_yz++; + if (fix->box_change & Fix::BOX_CHANGE_XZ) box_change_xz++; + if (fix->box_change & Fix::BOX_CHANGE_XY) box_change_xy++; } std::string mesg = "Must not have multiple fixes change box parameter "; @@ -174,12 +174,12 @@ void Domain::init() // check for fix deform deform_flag = deform_vremap = deform_groupbit = 0; - for (int i = 0; i < modify->nfix; i++) - if (utils::strmatch(modify->fix[i]->style,"^deform")) { + for (const auto &fix : fixes) + if (utils::strmatch(fix->style,"^deform")) { deform_flag = 1; - if (((FixDeform *) modify->fix[i])->remapflag == Domain::V_REMAP) { + if (((FixDeform *) fix)->remapflag == Domain::V_REMAP) { deform_vremap = 1; - deform_groupbit = modify->fix[i]->groupbit; + deform_groupbit = fix->groupbit; } } diff --git a/src/dump.cpp b/src/dump.cpp index d00c42086d..4c02aa3070 100644 --- a/src/dump.cpp +++ b/src/dump.cpp @@ -248,8 +248,8 @@ void Dump::init() reorderflag = 0; int gcmcflag = 0; - for (int i = 0; i < modify->nfix; i++) - if ((strcmp(modify->fix[i]->style,"gcmc") == 0)) + for (const auto &fix : modify->get_fix_list()) + if (utils::strmatch(fix->style,"^gcmc")) gcmcflag = 1; if (sortcol == 0 && atom->tag_consecutive() && !gcmcflag) { diff --git a/src/finish.cpp b/src/finish.cpp index 81ba59c1df..f3dbc83498 100644 --- a/src/finish.cpp +++ b/src/finish.cpp @@ -364,14 +364,13 @@ void Finish::end(int flag) } #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 - if ((ifix >= 0) && timer->has_full() && me == 0) { + if (fixomp && timer->has_full() && me == 0) { double thr_total = 0.0; ThrData *td; - FixOMP *fixomp = static_cast(lmp->modify->fix[ifix]); for (i=0; i < nthreads; ++i) { td = fixomp->get_thr(i); thr_total += td->get_time(Timer::ALL); diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp index 743fc705d2..16cd97469d 100644 --- a/src/fix_gravity.cpp +++ b/src/fix_gravity.cpp @@ -32,7 +32,6 @@ using namespace FixConst; using namespace MathConst; enum{CHUTE,SPHERICAL,VECTOR}; -enum{CONSTANT,EQUAL}; // same as FixPour /* ---------------------------------------------------------------------- */ diff --git a/src/fix_gravity.h b/src/fix_gravity.h index 31b6eb3478..389c5a1af1 100644 --- a/src/fix_gravity.h +++ b/src/fix_gravity.h @@ -37,6 +37,7 @@ class FixGravity : public Fix { virtual void post_force_respa(int, int, int); double compute_scalar(); void *extract(const char *, int &); + enum { CONSTANT, EQUAL }; protected: int style, disable; diff --git a/src/modify.cpp b/src/modify.cpp index 89a26025af..07301c4779 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -1103,6 +1103,45 @@ int Modify::find_fix_by_style(const char *style) 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 Modify::get_fix_by_style(const std::string &style) const +{ + std::vector 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 &Modify::get_fix_list() +{ + fix_list = std::vector(fix, fix+nfix); + return fix_list; +} + /* ---------------------------------------------------------------------- check for fix associated with package name in compiled list return 1 if found else 0 @@ -1359,6 +1398,45 @@ int Modify::find_compute_by_style(const char *style) 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 Modify::get_compute_by_style(const std::string &style) const +{ + std::vector 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 &Modify::get_compute_list() +{ + compute_list = std::vector(compute, compute+ncompute); + return compute_list; +} + /* ---------------------------------------------------------------------- clear invoked flag of all computes called everywhere that computes are used, before computes are invoked diff --git a/src/modify.h b/src/modify.h index 9446f285e7..0ea3bc2c7c 100644 --- a/src/modify.h +++ b/src/modify.h @@ -17,6 +17,7 @@ #include "pointers.h" #include +#include namespace LAMMPS_NS { @@ -107,16 +108,28 @@ class Modify : protected Pointers { void modify_fix(int, char **); void delete_fix(const std::string &); void delete_fix(int); + + // deprecated API int find_fix(const std::string &); int find_fix_by_style(const char *); + // new API + Fix *get_fix_by_id(const std::string &) const; + const std::vector get_fix_by_style(const std::string &) const; + const std::vector &get_fix_list(); Compute *add_compute(int, char **, int trysuffix = 1); Compute *add_compute(const std::string &, int trysuffix = 1); void modify_compute(int, char **); void delete_compute(const std::string &); void delete_compute(int); + + // deprecated API int find_compute(const std::string &); int find_compute_by_style(const char *); + // new API + Compute *get_compute_by_id(const std::string &) const; + const std::vector get_compute_by_style(const std::string &) const; + const std::vector &get_compute_list(); void clearstep_compute(); void addstep_compute(bigint); @@ -165,6 +178,10 @@ class Modify : protected Pointers { int index_permanent; // fix/compute index returned to library call + // vectors to be used for new-API accessors as wrapper + std::vectorfix_list; + std::vectorcompute_list; + void list_init(int, int &, int *&); void list_init_end_of_step(int, int &, int *&); void list_init_energy_couple(int &, int *&); diff --git a/src/variable.cpp b/src/variable.cpp index 22f792b5ec..125793a2d7 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -963,8 +963,7 @@ double Variable::compute_equal(int ivar) else if (style[ivar] == PYTHON) { int ifunc = python->find(data[ivar][0]); if (ifunc < 0) - print_var_error(FLERR,fmt::format("cannot find python function {}", - data[ivar][0]),ivar); + print_var_error(FLERR,fmt::format("cannot find python function {}",data[ivar][0]),ivar); python->invoke_function(ifunc,data[ivar][1]); value = atof(data[ivar][1]); } @@ -1324,10 +1323,9 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) // compute // ---------------- - if (strncmp(word,"c_",2) == 0 || strncmp(word,"C_",2) == 0) { + if (utils::strmatch(word,"^[Cc]_")) { if (domain->box_exist == 0) - print_var_error(FLERR,"Variable evaluation before " - "simulation box is defined",ivar); + print_var_error(FLERR,"Variable evaluation before simulation box is defined",ivar); // uppercase used to force access of // 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; if (word[0] == 'C') lowercase = 0; - int icompute = modify->find_compute(word+2); - if (icompute < 0) { - std::string mesg = "Invalid compute ID '"; - mesg += (word+2); - mesg += "' in variable formula"; - print_var_error(FLERR,mesg,ivar); - } - Compute *compute = modify->compute[icompute]; + Compute *compute = modify->get_compute_by_id(word+2); + if (!compute) + print_var_error(FLERR,fmt::format("Invalid compute ID '{}' in variable formula", word+2),ivar); + // parse zero or one or two trailing brackets // point i beyond last bracket @@ -1371,8 +1365,7 @@ double Variable::evaluate(char *str, Tree **tree, int ivar) if (update->whichflag == 0) { if (compute->invoked_scalar != update->ntimestep) - print_var_error(FLERR,"Compute used in variable between " - "runs is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_SCALAR)) { compute->compute_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 && compute->size_vector_variable == 0) - print_var_error(FLERR,"Variable formula compute vector " - "is accessed out-of-range",ivar,0); + print_var_error(FLERR,"Variable formula compute vector is accessed out-of-range",ivar,0); if (update->whichflag == 0) { if (compute->invoked_vector != update->ntimestep) - print_var_error(FLERR,"Compute used in variable between runs " - "is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { compute->compute_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 && compute->size_array_rows_variable == 0) - print_var_error(FLERR,"Variable formula compute array " - "is accessed out-of-range",ivar,0); + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); if (index2 > compute->size_array_cols) - print_var_error(FLERR,"Variable formula compute array " - "is accessed out-of-range",ivar,0); + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); if (update->whichflag == 0) { if (compute->invoked_array != update->ntimestep) - print_var_error(FLERR,"Compute used in variable between runs " - "is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { compute->compute_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) { if (tree == nullptr) - print_var_error(FLERR,"Compute global vector in " - "equal-style variable formula",ivar); + print_var_error(FLERR,"Compute global vector in equal-style variable formula",ivar); if (treetype == ATOM) - print_var_error(FLERR,"Compute global vector in " - "atom-style variable formula",ivar); + print_var_error(FLERR,"Compute global vector in atom-style variable formula",ivar); if (compute->size_vector == 0) - print_var_error(FLERR,"Variable formula compute " - "vector is zero length",ivar); + print_var_error(FLERR,"Variable formula compute vector is zero length",ivar); if (update->whichflag == 0) { if (compute->invoked_vector != update->ntimestep) - print_var_error(FLERR,"Compute used in variable between " - "runs is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { compute->compute_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) { if (tree == nullptr) - print_var_error(FLERR,"Compute global vector in " - "equal-style variable formula",ivar); + print_var_error(FLERR,"Compute global vector in equal-style variable formula",ivar); if (treetype == ATOM) - print_var_error(FLERR,"Compute global vector in " - "atom-style variable formula",ivar); + print_var_error(FLERR,"Compute global vector in atom-style variable formula",ivar); if (compute->size_array_rows == 0) - print_var_error(FLERR,"Variable formula compute array " - "is zero length",ivar); + print_var_error(FLERR,"Variable formula compute array is zero length",ivar); if (update->whichflag == 0) { if (compute->invoked_array != update->ntimestep) - print_var_error(FLERR,"Compute used in variable between " - "runs is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { compute->compute_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 (compute->invoked_peratom != update->ntimestep) - print_var_error(FLERR,"Compute used in variable " - "between runs is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { compute->compute_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) { if (tree == nullptr) - print_var_error(FLERR,"Per-atom compute in " - "equal-style variable formula",ivar); + print_var_error(FLERR,"Per-atom compute in equal-style variable formula",ivar); if (treetype == VECTOR) - print_var_error(FLERR,"Per-atom compute in " - "vector-style variable formula",ivar); + print_var_error(FLERR,"Per-atom compute in vector-style variable formula",ivar); if (update->whichflag == 0) { if (compute->invoked_peratom != update->ntimestep) - print_var_error(FLERR,"Compute used in variable " - "between runs is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { compute->compute_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) { if (tree == nullptr) - print_var_error(FLERR,"Per-atom compute in " - "equal-style variable formula",ivar); + print_var_error(FLERR,"Per-atom compute in equal-style variable formula",ivar); if (treetype == VECTOR) - print_var_error(FLERR,"Per-atom compute in " - "vector-style variable formula",ivar); + print_var_error(FLERR,"Per-atom compute in vector-style variable formula",ivar); if (index1 > compute->size_peratom_cols) - print_var_error(FLERR,"Variable formula compute array " - "is accessed out-of-range",ivar,0); + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); if (update->whichflag == 0) { if (compute->invoked_peratom != update->ntimestep) - print_var_error(FLERR,"Compute used in variable " - "between runs is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_PERATOM)) { compute->compute_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) { print_var_error(FLERR,"Cannot access local data via indexing",ivar); - } else print_var_error(FLERR, - "Mismatched compute in variable formula",ivar); + } else print_var_error(FLERR,"Mismatched compute in variable formula",ivar); // ---------------- // fix // ---------------- - } else if (strncmp(word,"f_",2) == 0 || strncmp(word,"F_",2) == 0) { + } else if (utils::strmatch(word,"^[fF]_")) { if (domain->box_exist == 0) - print_var_error(FLERR,"Variable evaluation before " - "simulation box is defined",ivar); + print_var_error(FLERR,"Variable evaluation before simulation box is defined",ivar); // uppercase used to force access of // 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; if (word[0] == 'F') lowercase = 0; - int ifix = modify->find_fix(word+2); - if (ifix < 0) - print_var_error(FLERR,fmt::format("Invalid fix ID '{}' in variable" - " formula",word+2),ivar); - Fix *fix = modify->fix[ifix]; + Fix *fix = modify->get_fix_by_id(word+2); + if (!fix) + print_var_error(FLERR,fmt::format("Invalid fix ID '{}' in variable formula",word+2),ivar); + // parse zero or one or two trailing brackets // point i beyond last bracket @@ -4057,19 +4026,17 @@ int Variable::special_function(char *word, char *contents, Tree **tree, *ptr1 = '\0'; } else index = 0; - int icompute = modify->find_compute(&args[0][2]); - if (icompute < 0) { + compute = modify->get_compute_by_id(&args[0][2]); + if (!compute) { std::string mesg = "Invalid compute ID '"; mesg += (args[0]+2); mesg += "' in variable formula"; print_var_error(FLERR,mesg,ivar); } - compute = modify->compute[icompute]; if (index == 0 && compute->vector_flag) { if (update->whichflag == 0) { if (compute->invoked_vector != update->ntimestep) - print_var_error(FLERR,"Compute used in variable between runs " - "is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_VECTOR)) { compute->compute_vector(); compute->invoked_flag |= Compute::INVOKED_VECTOR; @@ -4078,12 +4045,10 @@ int Variable::special_function(char *word, char *contents, Tree **tree, nstride = 1; } else if (index && compute->array_flag) { if (index > compute->size_array_cols) - print_var_error(FLERR,"Variable formula compute array " - "is accessed out-of-range",ivar,0); + print_var_error(FLERR,"Variable formula compute array is accessed out-of-range",ivar,0); if (update->whichflag == 0) { if (compute->invoked_array != update->ntimestep) - print_var_error(FLERR,"Compute used in variable between runs " - "is not current",ivar); + print_var_error(FLERR,"Compute used in variable between runs is not current",ivar); } else if (!(compute->invoked_flag & Compute::INVOKED_ARRAY)) { compute->compute_array(); compute->invoked_flag |= Compute::INVOKED_ARRAY; @@ -4102,14 +4067,13 @@ int Variable::special_function(char *word, char *contents, Tree **tree, *ptr1 = '\0'; } else index = 0; - int ifix = modify->find_fix(&args[0][2]); - if (ifix < 0) { + fix = modify->get_fix_by_id(&args[0][2]); + if (!fix) { std::string mesg = "Invalid fix ID '"; mesg += (args[0]+2); mesg += "' in variable formula"; print_var_error(FLERR,mesg,ivar); } - fix = modify->fix[ifix]; if (index == 0 && fix->vector_flag) { if (update->whichflag > 0 && update->ntimestep % fix->global_freq) { std::string mesg = "Fix with ID '"; diff --git a/src/velocity.cpp b/src/velocity.cpp index 4a449d00e7..34b8d5a700 100644 --- a/src/velocity.cpp +++ b/src/velocity.cpp @@ -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; loop_flag = ALL; scale_flag = 1; - rfix = -1; // read options from end of input line // 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 int initcomm = 0; - if (style == ZERO && rfix >= 0 && - utils::strmatch(modify->fix[rfix]->style,"^rigid.*/small.*")) initcomm = 1; + if (style == ZERO && rigid_fix && utils::strmatch(rigid_fix->style,"^rigid.*/small.*")) initcomm = 1; if ((style == CREATE || style == SET) && temperature && strcmp(temperature->style,"temp/cs") == 0) initcomm = 1; @@ -154,7 +152,6 @@ void Velocity::init_external(const char *extgroup) momentum_flag = 1; rotation_flag = 0; loop_flag = ALL; - rfix = -1; scale_flag = 1; bias_flag = 0; } @@ -692,21 +689,21 @@ void Velocity::ramp(int /*narg*/, char **arg) void Velocity::zero(int /*narg*/, char **arg) { if (strcmp(arg[0],"linear") == 0) { - if (rfix < 0) zero_momentum(); - else if (utils::strmatch(modify->fix[rfix]->style,"^rigid/small")) { - modify->fix[rfix]->setup_pre_neighbor(); - modify->fix[rfix]->zero_momentum(); - } else if (utils::strmatch(modify->fix[rfix]->style,"^rigid")) { - modify->fix[rfix]->zero_momentum(); + if (!rigid_fix) zero_momentum(); + else if (utils::strmatch(rigid_fix->style,"^rigid/small")) { + rigid_fix->setup_pre_neighbor(); + rigid_fix->zero_momentum(); + } else if (utils::strmatch(rigid_fix->style,"^rigid")) { + rigid_fix->zero_momentum(); } else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID"); } else if (strcmp(arg[0],"angular") == 0) { - if (rfix < 0) zero_rotation(); - else if (utils::strmatch(modify->fix[rfix]->style,"^rigid/small")) { - modify->fix[rfix]->setup_pre_neighbor(); - modify->fix[rfix]->zero_rotation(); - } else if (utils::strmatch(modify->fix[rfix]->style,"^rigid")) { - modify->fix[rfix]->zero_rotation(); + if (!rigid_fix) zero_rotation(); + else if (utils::strmatch(rigid_fix->style,"^rigid/small")) { + rigid_fix->setup_pre_neighbor(); + rigid_fix->zero_rotation(); + } else if (utils::strmatch(rigid_fix->style,"^rigid")) { + rigid_fix->zero_rotation(); } else error->all(FLERR,"Velocity rigid used with non-rigid fix-ID"); } else error->all(FLERR,"Illegal velocity command"); @@ -843,11 +840,10 @@ void Velocity::options(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"temp") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); - int icompute = modify->find_compute(arg[iarg+1]); - if (icompute < 0) error->all(FLERR,"Could not find velocity temperature ID"); - temperature = modify->compute[icompute]; + temperature = modify->get_compute_by_id(arg[iarg+1]); + if (!temperature) error->all(FLERR,"Could not find velocity temperature compute ID"); 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; } else if (strcmp(arg[iarg],"bias") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); @@ -862,8 +858,8 @@ void Velocity::options(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"rigid") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); - rfix = modify->find_fix(arg[iarg+1]); - if (rfix < 0) error->all(FLERR,"Fix ID for velocity does not exist"); + rigid_fix = modify->get_fix_by_id(arg[iarg+1]); + if (!rigid_fix) error->all(FLERR,"Fix ID {} for velocity does not exist", arg[iarg+1]); iarg += 2; } else if (strcmp(arg[iarg],"units") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal velocity command"); diff --git a/src/velocity.h b/src/velocity.h index f0a756350c..cce8508021 100644 --- a/src/velocity.h +++ b/src/velocity.h @@ -36,8 +36,9 @@ class Velocity : public Command { int igroup, groupbit; int style; 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; + class Fix *rigid_fix; class Compute *temperature; void set(int, char **); diff --git a/src/verlet.cpp b/src/verlet.cpp index b9b0b392e1..3ea9d72a9e 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -68,8 +68,7 @@ void Verlet::init() // detect if fix omp is present for clearing force arrays - int ifix = modify->find_fix("package_omp"); - if (ifix >= 0) external_force_clear = 1; + if (modify->get_fix_by_id("package_omp")) external_force_clear = 1; // set flags for arrays to clear in force_clear()