From aa4787f6046ae9071f56d9f55a83d7643a8768ce Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 16 Apr 2022 00:12:43 -0400 Subject: [PATCH] complete region handling refactor --- lib/atc/LammpsInterface.cpp | 49 +- src/ADIOS/dump_custom_adios.cpp | 2 +- src/EFF/compute_temp_region_eff.cpp | 95 ++-- src/EFF/compute_temp_region_eff.h | 2 +- src/EXTRA-FIX/fix_electron_stopping.cpp | 102 ++-- src/EXTRA-FIX/fix_electron_stopping.h | 5 +- src/EXTRA-FIX/fix_oneway.cpp | 58 ++- src/EXTRA-FIX/fix_oneway.h | 4 +- src/EXTRA-FIX/fix_wall_region_ees.cpp | 197 ++++---- src/EXTRA-FIX/fix_wall_region_ees.h | 2 +- src/GRANULAR/fix_pour.cpp | 623 ++++++++++++------------ src/GRANULAR/fix_pour.h | 4 +- src/GRANULAR/fix_wall_gran_region.cpp | 196 ++++---- src/MACHDYN/fix_smd_setvel.cpp | 490 +++++++++---------- src/MACHDYN/fix_smd_setvel.h | 3 +- src/MC/fix_atom_swap.cpp | 270 +++++----- src/MC/fix_atom_swap.h | 3 +- src/MC/fix_gcmc.cpp | 101 ++-- src/MC/fix_gcmc.h | 21 +- src/MC/fix_widom.cpp | 82 ++-- src/MC/fix_widom.h | 17 +- src/PLUGIN/plugin.cpp | 2 +- src/REPLICA/hyper.cpp | 9 +- src/REPLICA/prd.cpp | 9 +- src/RIGID/fix_ehex.cpp | 372 +++++++------- src/RIGID/fix_ehex.h | 2 +- src/VTK/dump_vtk.cpp | 24 +- src/domain.cpp | 97 ++-- src/domain.h | 10 +- src/group.cpp | 77 +-- src/group.h | 20 +- src/info.cpp | 26 +- src/library.cpp | 16 +- src/math_const.h | 1 + src/modify.cpp | 4 +- src/set.cpp | 8 +- src/variable.cpp | 144 +++--- src/variable.h | 10 +- unittest/commands/test_groups.cpp | 11 +- 39 files changed, 1504 insertions(+), 1664 deletions(-) diff --git a/lib/atc/LammpsInterface.cpp b/lib/atc/LammpsInterface.cpp index 8f7d20361c..9e2df3e46c 100644 --- a/lib/atc/LammpsInterface.cpp +++ b/lib/atc/LammpsInterface.cpp @@ -387,7 +387,7 @@ double LammpsInterface::atom_quantity_conversion(FundamentalAtomQuantity quantit int LammpsInterface::dimension() const { return lammps_->domain->dimension; } -int LammpsInterface::nregion() const { return lammps_->domain->nregion; } +int LammpsInterface::nregion() const { return lammps_->domain->get_region_list().size(); } void LammpsInterface::box_bounds(double & boxxlo, double & boxxhi, double & boxylo, double & boxyhi, @@ -527,14 +527,15 @@ void LammpsInterface::box_periodicity(int & xperiodic, zperiodic = lammps_->domain->zperiodic; } -int LammpsInterface::region_id(const char * regionName) const { - int nregion = this->nregion(); - for (int iregion = 0; iregion < nregion; iregion++) { - if (strcmp(regionName, region_name(iregion)) == 0) { +int LammpsInterface::region_id(const char *regionName) const { + auto regions = lammps_->domain->get_region_list(); + int iregion = 0; + for (auto reg : regions) { + if (strcmp(regionName, reg->id) == 0) { return iregion; } + ++iregion; } - throw ATC_Error("Region has not been defined"); return -1; } @@ -1322,61 +1323,73 @@ int** LammpsInterface::bond_list() const { return lammps_->neighbor->bondlist; char * LammpsInterface::region_name(int iRegion) const { - return lammps_->domain->regions[iRegion]->id; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->id; } char * LammpsInterface::region_style(int iRegion) const { - return lammps_->domain->regions[iRegion]->style; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->style; } double LammpsInterface::region_xlo(int iRegion) const { - return lammps_->domain->regions[iRegion]->extent_xlo; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->extent_xlo; } double LammpsInterface::region_xhi(int iRegion) const { - return lammps_->domain->regions[iRegion]->extent_xhi; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->extent_xhi; } double LammpsInterface::region_ylo(int iRegion) const { - return lammps_->domain->regions[iRegion]->extent_ylo; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->extent_ylo; } double LammpsInterface::region_yhi(int iRegion) const { - return lammps_->domain->regions[iRegion]->extent_yhi; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->extent_yhi; } double LammpsInterface::region_zlo(int iRegion) const { - return lammps_->domain->regions[iRegion]->extent_zlo; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->extent_zlo; } double LammpsInterface::region_zhi(int iRegion) const { - return lammps_->domain->regions[iRegion]->extent_zhi; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->extent_zhi; } double LammpsInterface::region_xscale(int iRegion) const { - return lammps_->domain->regions[iRegion]->xscale; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->xscale; } double LammpsInterface::region_yscale(int iRegion) const { - return lammps_->domain->regions[iRegion]->yscale; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->yscale; } double LammpsInterface::region_zscale(int iRegion) const { - return lammps_->domain->regions[iRegion]->zscale; + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->zscale; } int LammpsInterface::region_match(int iRegion, double x, double y, double z) const { - return lammps_->domain->regions[iRegion]->match(x,y,z); + auto regions = lammps_->domain->get_region_list(); + return regions[iRegion]->match(x,y,z); } // ----------------------------------------------------------------- diff --git a/src/ADIOS/dump_custom_adios.cpp b/src/ADIOS/dump_custom_adios.cpp index 5249021ed4..0ffe612f2c 100644 --- a/src/ADIOS/dump_custom_adios.cpp +++ b/src/ADIOS/dump_custom_adios.cpp @@ -273,7 +273,7 @@ void DumpCustomADIOS::init_style() } // set index and check validity of region - if (idregion && !domain->find_region(idregion)) + if (idregion && !domain->get_region_by_id(idregion)) error->all(FLERR, "Region {} for dump custom/adios does not exist", idregion); /* Define the group of variables for the atom style here since it's a fixed diff --git a/src/EFF/compute_temp_region_eff.cpp b/src/EFF/compute_temp_region_eff.cpp index 3b5147b900..de61abc4b7 100644 --- a/src/EFF/compute_temp_region_eff.cpp +++ b/src/EFF/compute_temp_region_eff.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -16,7 +15,6 @@ Contributing author: Andres Jaramillo-Botero (Caltech) ------------------------------------------------------------------------- */ - #include "compute_temp_region_eff.h" #include "atom.h" @@ -33,16 +31,15 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) + Compute(lmp, narg, arg), region(nullptr), idregion(nullptr) { if (!atom->electron_flag) - error->all(FLERR,"Compute temp/region/eff requires atom style electron"); + error->all(FLERR, "Compute temp/region/eff requires atom style electron"); - if (narg != 4) error->all(FLERR,"Illegal compute temp/region/eff command"); + if (narg != 4) error->all(FLERR, "Illegal compute temp/region/eff command"); - iregion = domain->find_region(arg[3]); - if (iregion == -1) - error->all(FLERR,"Region ID for compute temp/region/eff does not exist"); + region = domain->get_region_by_id(arg[3]); + if (!region) error->all(FLERR, "Region {} for compute temp/region/eff does not exist", arg[3]); idregion = utils::strdup(arg[3]); scalar_flag = vector_flag = 1; @@ -61,9 +58,9 @@ ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) : ComputeTempRegionEff::~ComputeTempRegionEff() { - delete [] idregion; + delete[] idregion; memory->destroy(vbiasall); - delete [] vector; + delete[] vector; } /* ---------------------------------------------------------------------- */ @@ -72,9 +69,8 @@ void ComputeTempRegionEff::init() { // set index and check validity of region - iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for compute temp/region/eff does not exist"); + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for compute temp/region/eff does not exist", idregion); } /* ---------------------------------------------------------------------- */ @@ -90,7 +86,7 @@ void ComputeTempRegionEff::setup() void ComputeTempRegionEff::dof_remove_pre() { - domain->regions[iregion]->prematch(); + region->prematch(); } /* ---------------------------------------------------------------------- */ @@ -98,7 +94,7 @@ void ComputeTempRegionEff::dof_remove_pre() int ComputeTempRegionEff::dof_remove(int i) { double *x = atom->x[i]; - if (domain->regions[iregion]->match(x[0],x[1],x[2])) return 0; + if (region->match(x[0], x[1], x[2])) return 0; return 1; } @@ -116,9 +112,8 @@ double ComputeTempRegionEff::compute_scalar() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - double mefactor = domain->dimension/4.0; + double mefactor = domain->dimension / 4.0; - Region *region = domain->regions[iregion]; region->prematch(); int count = 0; @@ -127,34 +122,35 @@ double ComputeTempRegionEff::compute_scalar() if (mass) { for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { + if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) { count++; - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) * - mass[type[i]]; - if (abs(spin[i])==1) { - t += mefactor*mass[type[i]]*ervel[i]*ervel[i]; + t += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]) * mass[type[i]]; + if (abs(spin[i]) == 1) { + t += mefactor * mass[type[i]] * ervel[i] * ervel[i]; ecount++; } } } - double tarray[2],tarray_all[2]; + double tarray[2], tarray_all[2]; // Assume 3/2 k T per nucleus - tarray[0] = count-ecount; + tarray[0] = count - ecount; tarray[1] = t; - MPI_Allreduce(tarray,tarray_all,2,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(tarray, tarray_all, 2, MPI_DOUBLE, MPI_SUM, world); dof = domain->dimension * tarray_all[0] - extra_dof; if (dof < 0.0 && tarray_all[0] > 0.0) - error->all(FLERR,"Temperature compute degrees of freedom < 0"); + error->all(FLERR, "Temperature compute degrees of freedom < 0"); int one = 0; for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { - if (abs(spin[i])==1) one++; + if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) { + if (abs(spin[i]) == 1) one++; } - if (dof > 0.0) scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz); - else scalar = 0.0; + if (dof > 0.0) + scalar = force->mvv2e * tarray_all[1] / (dof * force->boltz); + else + scalar = 0.0; return scalar; } @@ -174,33 +170,32 @@ void ComputeTempRegionEff::compute_vector() int *type = atom->type; int *mask = atom->mask; int nlocal = atom->nlocal; - double mefactor = domain->dimension/4.0; + double mefactor = domain->dimension / 4.0; - Region *region = domain->regions[iregion]; region->prematch(); - double massone,t[6]; + double massone, t[6]; for (i = 0; i < 6; i++) t[i] = 0.0; for (i = 0; i < nlocal; i++) - if (mask[i] & groupbit && region->match(x[i][0],x[i][1],x[i][2])) { + if (mask[i] & groupbit && region->match(x[i][0], x[i][1], x[i][2])) { massone = mass[type[i]]; - t[0] += massone * v[i][0]*v[i][0]; - t[1] += massone * v[i][1]*v[i][1]; - t[2] += massone * v[i][2]*v[i][2]; - t[3] += massone * v[i][0]*v[i][1]; - t[4] += massone * v[i][0]*v[i][2]; - t[5] += massone * v[i][1]*v[i][2]; + t[0] += massone * v[i][0] * v[i][0]; + t[1] += massone * v[i][1] * v[i][1]; + t[2] += massone * v[i][2] * v[i][2]; + t[3] += massone * v[i][0] * v[i][1]; + t[4] += massone * v[i][0] * v[i][2]; + t[5] += massone * v[i][1] * v[i][2]; - if (abs(spin[i])==1) { - t[0] += mefactor * massone * ervel[i]*ervel[i]; - t[1] += mefactor * massone * ervel[i]*ervel[i]; - t[2] += mefactor * massone * ervel[i]*ervel[i]; + if (abs(spin[i]) == 1) { + t[0] += mefactor * massone * ervel[i] * ervel[i]; + t[1] += mefactor * massone * ervel[i] * ervel[i]; + t[2] += mefactor * massone * ervel[i] * ervel[i]; } } - MPI_Allreduce(t,vector,6,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(t, vector, 6, MPI_DOUBLE, MPI_SUM, world); for (i = 0; i < 6; i++) vector[i] *= force->mvv2e; } @@ -212,7 +207,7 @@ void ComputeTempRegionEff::compute_vector() void ComputeTempRegionEff::remove_bias(int i, double *v) { double *x = atom->x[i]; - if (domain->regions[iregion]->match(x[0],x[1],x[2])) + if (region->match(x[0], x[1], x[2])) vbias[0] = vbias[1] = vbias[2] = 0.0; else { vbias[0] = v[0]; @@ -236,14 +231,12 @@ void ComputeTempRegionEff::remove_bias_all() if (atom->nmax > maxbias) { memory->destroy(vbiasall); maxbias = atom->nmax; - memory->create(vbiasall,maxbias,3,"temp/region:vbiasall"); + memory->create(vbiasall, maxbias, 3, "temp/region:vbiasall"); } - Region *region = domain->regions[iregion]; - for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (region->match(x[i][0],x[i][1],x[i][2])) + if (region->match(x[i][0], x[i][1], x[i][2])) vbiasall[i][0] = vbiasall[i][1] = vbiasall[i][2] = 0.0; else { vbiasall[i][0] = v[i][0]; @@ -289,6 +282,6 @@ void ComputeTempRegionEff::restore_bias_all() double ComputeTempRegionEff::memory_usage() { - double bytes = (double)maxbias * sizeof(double); + double bytes = (double) maxbias * sizeof(double); return bytes; } diff --git a/src/EFF/compute_temp_region_eff.h b/src/EFF/compute_temp_region_eff.h index bc8f930374..a3dd7e5dd4 100644 --- a/src/EFF/compute_temp_region_eff.h +++ b/src/EFF/compute_temp_region_eff.h @@ -42,7 +42,7 @@ class ComputeTempRegionEff : public Compute { double memory_usage() override; protected: - int iregion; + class Region *region; char *idregion; }; diff --git a/src/EXTRA-FIX/fix_electron_stopping.cpp b/src/EXTRA-FIX/fix_electron_stopping.cpp index a61cedcf27..e9cc501744 100644 --- a/src/EXTRA-FIX/fix_electron_stopping.cpp +++ b/src/EXTRA-FIX/fix_electron_stopping.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -39,63 +38,53 @@ using namespace LAMMPS_NS; using namespace FixConst; -#define MAXLINE 1024 - /* ---------------------------------------------------------------------- */ FixElectronStopping::FixElectronStopping(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) + Fix(lmp, narg, arg), elstop_ranges(nullptr), idregion(nullptr), region(nullptr), list(nullptr) { - scalar_flag = 1; // Has compute_scalar - global_freq = 1; // SeLoss computed every step - extscalar = 0; // SeLoss compute_scalar is intensive - nevery = 1; // Run fix every step - + scalar_flag = 1; // Has compute_scalar + global_freq = 1; // SeLoss computed every step + extscalar = 0; // SeLoss compute_scalar is intensive + nevery = 1; // Run fix every step // args: 0 = fix ID, 1 = group ID, 2 = "electron/stopping" // 3 = Ecut, 4 = file path // optional rest: "region" // "minneigh" - if (narg < 5) error->all(FLERR, - "Illegal fix electron/stopping command: too few arguments"); + if (narg < 5) error->all(FLERR, "Illegal fix electron/stopping command: too few arguments"); - Ecut = utils::numeric(FLERR, arg[3],false,lmp); - if (Ecut <= 0.0) error->all(FLERR, - "Illegal fix electron/stopping command: Ecut <= 0"); + Ecut = utils::numeric(FLERR, arg[3], false, lmp); + if (Ecut <= 0.0) error->all(FLERR, "Illegal fix electron/stopping command: Ecut <= 0"); int iarg = 5; - iregion = -1; minneigh = 1; bool minneighflag = false; while (iarg < narg) { if (strcmp(arg[iarg], "region") == 0) { - if (iregion >= 0) error->all(FLERR, - "Illegal fix electron/stopping command: region given twice"); - if (iarg+2 > narg) error->all(FLERR, - "Illegal fix electron/stopping command: region name missing"); - iregion = domain->find_region(arg[iarg+1]); - if (iregion < 0) error->all(FLERR, - "Region ID for fix electron/stopping does not exist"); + if (region) error->all(FLERR, "Illegal fix electron/stopping command: region given twice"); + if (iarg + 2 > narg) + error->all(FLERR, "Illegal fix electron/stopping command: region name missing"); + region = domain->get_region_by_id(arg[iarg + 1]); + if (!region) + error->all(FLERR, "Region {} for fix electron/stopping does not exist", arg[iarg + 1]); + idregion = utils::strdup(arg[iarg + 1]); iarg += 2; - } - else if (strcmp(arg[iarg], "minneigh") == 0) { - if (minneighflag) error->all(FLERR, - "Illegal fix electron/stopping command: minneigh given twice"); + } else if (strcmp(arg[iarg], "minneigh") == 0) { + if (minneighflag) + error->all(FLERR, "Illegal fix electron/stopping command: minneigh given twice"); minneighflag = true; - if (iarg+2 > narg) error->all(FLERR, - "Illegal fix electron/stopping command: minneigh number missing"); - minneigh = utils::inumeric(FLERR, arg[iarg+1],false,lmp); - if (minneigh < 0) error->all(FLERR, - "Illegal fix electron/stopping command: minneigh < 0"); + if (iarg + 2 > narg) + error->all(FLERR, "Illegal fix electron/stopping command: minneigh number missing"); + minneigh = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (minneigh < 0) error->all(FLERR, "Illegal fix electron/stopping command: minneigh < 0"); iarg += 2; - } - else error->all(FLERR, - "Illegal fix electron/stopping command: unknown argument"); + } else + error->all(FLERR, "Illegal fix electron/stopping command: unknown argument"); } - // Read the input file for energy ranges and stopping powers. // First proc 0 reads the file, then bcast to others. const int ncol = atom->ntypes + 1; @@ -105,13 +94,12 @@ FixElectronStopping::FixElectronStopping(LAMMPS *lmp, int narg, char **arg) : read_table(arg[4]); } - MPI_Bcast(&maxlines, 1 , MPI_INT, 0, world); - MPI_Bcast(&table_entries, 1 , MPI_INT, 0, world); + MPI_Bcast(&maxlines, 1, MPI_INT, 0, world); + MPI_Bcast(&table_entries, 1, MPI_INT, 0, world); - if (comm->me != 0) - memory->create(elstop_ranges, ncol, maxlines, "electron/stopping:table"); + if (comm->me != 0) memory->create(elstop_ranges, ncol, maxlines, "electron/stopping:table"); - MPI_Bcast(&elstop_ranges[0][0], ncol*maxlines, MPI_DOUBLE, 0, world); + MPI_Bcast(&elstop_ranges[0][0], ncol * maxlines, MPI_DOUBLE, 0, world); } /* ---------------------------------------------------------------------- */ @@ -136,6 +124,10 @@ void FixElectronStopping::init() { SeLoss_sync_flag = 0; SeLoss = 0.0; + if (idregion) { + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for fix electron/stopping does not exist", idregion); + } // need an occasional full neighbor list neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); @@ -176,18 +168,17 @@ void FixElectronStopping::post_force(int /*vflag*/) int itype = type[i]; double massone = (atom->rmass) ? atom->rmass[i] : atom->mass[itype]; - double v2 = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; + double v2 = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]; double energy = 0.5 * force->mvv2e * massone * v2; if (energy < Ecut) continue; if (energy < elstop_ranges[0][0]) continue; - if (energy > elstop_ranges[0][table_entries - 1]) error->one(FLERR, - "Atom kinetic energy too high for fix electron/stopping"); + if (energy > elstop_ranges[0][table_entries - 1]) + error->one(FLERR, "Atom kinetic energy too high for fix electron/stopping"); - if (iregion >= 0) { + if (region) { // Only apply in the given region - if (domain->regions[iregion]->match(x[i][0], x[i][1], x[i][2]) != 1) - continue; + if (region->match(x[i][0], x[i][1], x[i][2]) != 1) continue; } // Binary search to find correct energy range @@ -196,8 +187,10 @@ void FixElectronStopping::post_force(int /*vflag*/) while (true) { int ihalf = idown + (iup - idown) / 2; if (ihalf == idown) break; - if (elstop_ranges[0][ihalf] < energy) idown = ihalf; - else iup = ihalf; + if (elstop_ranges[0][ihalf] < energy) + idown = ihalf; + else + iup = ihalf; } double Se_lo = elstop_ranges[itype][idown]; @@ -215,7 +208,7 @@ void FixElectronStopping::post_force(int /*vflag*/) f[i][1] += v[i][1] * factor; f[i][2] += v[i][2] * factor; - SeLoss += Se * vabs * dt; // very rough approx + SeLoss += Se * vabs * dt; // very rough approx } } @@ -254,19 +247,17 @@ void FixElectronStopping::read_table(const char *file) ValueTokenizer values(line); elstop_ranges[0][nlines] = values.next_double(); if (elstop_ranges[0][nlines] <= oldvalue) - throw TokenizerException("energy values must be positive and in ascending order",line); + throw TokenizerException("energy values must be positive and in ascending order", line); oldvalue = elstop_ranges[0][nlines]; - for (int i = 1; i < ncol; ++i) - elstop_ranges[i][nlines] = values.next_double(); + for (int i = 1; i < ncol; ++i) elstop_ranges[i][nlines] = values.next_double(); ++nlines; } } catch (std::exception &e) { error->one(FLERR, "Problem parsing electron stopping data: {}", e.what()); } - if (nlines == 0) - error->one(FLERR, "Did not find any data in electron/stopping table file"); + if (nlines == 0) error->one(FLERR, "Did not find any data in electron/stopping table file"); table_entries = nlines; } @@ -281,8 +272,7 @@ void FixElectronStopping::grow_table() double **new_array; memory->create(new_array, ncol, new_maxlines, "electron/stopping:table"); - for (int i = 0; i < ncol; i++) - memcpy(new_array[i], elstop_ranges[i], maxlines*sizeof(double)); + for (int i = 0; i < ncol; i++) memcpy(new_array[i], elstop_ranges[i], maxlines * sizeof(double)); memory->destroy(elstop_ranges); elstop_ranges = new_array; diff --git a/src/EXTRA-FIX/fix_electron_stopping.h b/src/EXTRA-FIX/fix_electron_stopping.h index 015d444996..762435e966 100644 --- a/src/EXTRA-FIX/fix_electron_stopping.h +++ b/src/EXTRA-FIX/fix_electron_stopping.h @@ -52,8 +52,9 @@ class FixElectronStopping : public Fix { double **elstop_ranges; // [ 0][i]: energies // [>0][i]: stopping powers per type - int iregion; // region index if used, else -1 - int minneigh; // minimum number of neighbors + char *idregion; // region id + class Region *region; // region pointer if used, else NULL + int minneigh; // minimum number of neighbors class NeighList *list; }; diff --git a/src/EXTRA-FIX/fix_oneway.cpp b/src/EXTRA-FIX/fix_oneway.cpp index 465813ef57..f4a29f4437 100644 --- a/src/EXTRA-FIX/fix_oneway.cpp +++ b/src/EXTRA-FIX/fix_oneway.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -28,35 +27,36 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{NONE=-1,X=0,Y=1,Z=2,XYZMASK=3,MINUS=4,PLUS=0}; +enum { NONE = -1, X = 0, Y = 1, Z = 2, XYZMASK = 3, MINUS = 4, PLUS = 0 }; /* ---------------------------------------------------------------------- */ -FixOneWay::FixOneWay(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) +FixOneWay::FixOneWay(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), region(nullptr), idregion(nullptr) { direction = NONE; - regionidx = 0; - regionstr = nullptr; - if (narg < 6) error->all(FLERR,"Illegal fix oneway command"); + if (narg < 6) error->all(FLERR, "Illegal fix oneway command"); - nevery = utils::inumeric(FLERR,arg[3],false,lmp); - if (nevery < 1) error->all(FLERR,"Illegal fix oneway command"); + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + if (nevery < 1) error->all(FLERR, "Illegal fix oneway command"); - regionstr = utils::strdup(arg[4]); + idregion = utils::strdup(arg[4]); + if (!domain->get_region_by_id(idregion)) + error->all(FLERR, "Region {} for fix oneway does not exist", idregion); - if (strcmp(arg[5], "x") == 0) direction = X|PLUS; - if (strcmp(arg[5], "X") == 0) direction = X|PLUS; - if (strcmp(arg[5], "y") == 0) direction = Y|PLUS; - if (strcmp(arg[5], "Y") == 0) direction = Y|PLUS; - if (strcmp(arg[5], "z") == 0) direction = Z|PLUS; - if (strcmp(arg[5], "Z") == 0) direction = Z|PLUS; - if (strcmp(arg[5],"-x") == 0) direction = X|MINUS; - if (strcmp(arg[5],"-X") == 0) direction = X|MINUS; - if (strcmp(arg[5],"-y") == 0) direction = Y|MINUS; - if (strcmp(arg[5],"-Y") == 0) direction = Y|MINUS; - if (strcmp(arg[5],"-z") == 0) direction = Z|MINUS; - if (strcmp(arg[5],"-Z") == 0) direction = Z|MINUS; + if (strcmp(arg[5], "x") == 0) direction = X | PLUS; + if (strcmp(arg[5], "X") == 0) direction = X | PLUS; + if (strcmp(arg[5], "y") == 0) direction = Y | PLUS; + if (strcmp(arg[5], "Y") == 0) direction = Y | PLUS; + if (strcmp(arg[5], "z") == 0) direction = Z | PLUS; + if (strcmp(arg[5], "Z") == 0) direction = Z | PLUS; + if (strcmp(arg[5], "-x") == 0) direction = X | MINUS; + if (strcmp(arg[5], "-X") == 0) direction = X | MINUS; + if (strcmp(arg[5], "-y") == 0) direction = Y | MINUS; + if (strcmp(arg[5], "-Y") == 0) direction = Y | MINUS; + if (strcmp(arg[5], "-z") == 0) direction = Z | MINUS; + if (strcmp(arg[5], "-Z") == 0) direction = Z | MINUS; global_freq = nevery; } @@ -65,7 +65,7 @@ FixOneWay::FixOneWay(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) FixOneWay::~FixOneWay() { - delete[] regionstr; + delete[] idregion; } /* ---------------------------------------------------------------------- */ @@ -79,26 +79,25 @@ int FixOneWay::setmask() void FixOneWay::init() { - regionidx = domain->find_region(regionstr); - if (regionidx < 0) - error->all(FLERR,"Region for fix oneway does not exist"); + region = domain->get_region_by_id(idregion); + if (!region) + error->all(FLERR, "Region {} for fix oneway does not exist", idregion); } /* ---------------------------------------------------------------------- */ void FixOneWay::end_of_step() { - Region *region = domain->regions[regionidx]; region->prematch(); const int idx = direction & XYZMASK; - const double * const * const x = atom->x; - double * const * const v = atom->v; + const double *const *const x = atom->x; + double *const *const v = atom->v; const int *mask = atom->mask; const int nlocal = atom->nlocal; for (int i = 0; i < nlocal; ++i) { - if ((mask[i] & groupbit) && region->match(x[i][0],x[i][1],x[i][2])) { + if ((mask[i] & groupbit) && region->match(x[i][0], x[i][1], x[i][2])) { if (direction & MINUS) { if (v[i][idx] > 0.0) v[i][idx] = -v[i][idx]; } else { @@ -107,4 +106,3 @@ void FixOneWay::end_of_step() } } } - diff --git a/src/EXTRA-FIX/fix_oneway.h b/src/EXTRA-FIX/fix_oneway.h index ad6a5dc106..c9017aad52 100644 --- a/src/EXTRA-FIX/fix_oneway.h +++ b/src/EXTRA-FIX/fix_oneway.h @@ -34,8 +34,8 @@ class FixOneWay : public Fix { protected: int direction; - int regionidx; - char *regionstr; + class Region *region; + char *idregion; }; } // namespace LAMMPS_NS diff --git a/src/EXTRA-FIX/fix_wall_region_ees.cpp b/src/EXTRA-FIX/fix_wall_region_ees.cpp index eb1ede5c4a..a9f6205373 100644 --- a/src/EXTRA-FIX/fix_wall_region_ees.cpp +++ b/src/EXTRA-FIX/fix_wall_region_ees.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -23,6 +22,7 @@ #include "domain.h" #include "error.h" #include "math_extra.h" +#include "math_special.h" #include "region.h" #include "respa.h" #include "update.h" @@ -31,13 +31,14 @@ using namespace LAMMPS_NS; using namespace FixConst; +using MathSpecial::powint; /* ---------------------------------------------------------------------- */ FixWallRegionEES::FixWallRegionEES(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) + Fix(lmp, narg, arg), idregion(nullptr), region(nullptr) { - if (narg != 7) error->all(FLERR,"Illegal fix wall/region/ees command"); + if (narg != 7) error->all(FLERR, "Illegal fix wall/region/ees command"); scalar_flag = 1; vector_flag = 1; @@ -49,15 +50,14 @@ FixWallRegionEES::FixWallRegionEES(LAMMPS *lmp, int narg, char **arg) : // parse args - iregion = domain->find_region(arg[3]); - if (iregion == -1) - error->all(FLERR,"Region ID for fix wall/region/ees does not exist"); + region = domain->get_region_by_id(arg[3]); + if (!region) error->all(FLERR, "Region {} for fix wall/region/ees does not exist", arg[3]); idregion = utils::strdup(arg[3]); - epsilon = utils::numeric(FLERR,arg[4],false,lmp); - sigma = utils::numeric(FLERR,arg[5],false,lmp); - cutoff = utils::numeric(FLERR,arg[6],false,lmp); + epsilon = utils::numeric(FLERR, arg[4], false, lmp); + sigma = utils::numeric(FLERR, arg[5], false, lmp); + cutoff = utils::numeric(FLERR, arg[6], false, lmp); - if (cutoff <= 0.0) error->all(FLERR,"Fix wall/region/ees cutoff <= 0.0"); + if (cutoff <= 0.0) error->all(FLERR, "Fix wall/region/ees cutoff <= 0.0"); eflag = 0; ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0; @@ -67,7 +67,7 @@ FixWallRegionEES::FixWallRegionEES(LAMMPS *lmp, int narg, char **arg) : FixWallRegionEES::~FixWallRegionEES() { - delete [] idregion; + delete[] idregion; } /* ---------------------------------------------------------------------- */ @@ -87,13 +87,11 @@ void FixWallRegionEES::init() { // set index and check validity of region - iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for fix wall/region/ees does not exist"); + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for fix wall/region/ees does not exist", idregion); - avec = dynamic_cast( atom->style_match("ellipsoid")); - if (!avec) - error->all(FLERR,"Fix wall/region/ees requires atom style ellipsoid"); + avec = dynamic_cast(atom->style_match("ellipsoid")); + if (!avec) error->all(FLERR, "Fix wall/region/ees requires atom style ellipsoid"); // check that all particles are finite-size ellipsoids // no point particles allowed, spherical is OK @@ -105,33 +103,33 @@ void FixWallRegionEES::init() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (ellipsoid[i] < 0) - error->one(FLERR,"Fix wall/region/ees requires extended particles"); + error->one(FLERR, "Fix wall/region/ees requires only extended particles"); // setup coefficients - coeff1 = ( 2. / 4725. ) * epsilon * pow(sigma,12.0); - coeff2 = ( 1. / 24. ) * epsilon * pow(sigma,6.0); - coeff3 = ( 2. / 315. ) * epsilon * pow(sigma,12.0); - coeff4 = ( 1. / 3. ) * epsilon * pow(sigma,6.0); - coeff5 = ( 4. / 315. ) * epsilon * pow(sigma,12.0); - coeff6 = ( 1. / 12. ) * epsilon * pow(sigma,6.0); + coeff1 = (2.0 / 4725.0) * epsilon * powint(sigma, 12); + coeff2 = (1.0 / 24.0) * epsilon * powint(sigma, 6); + coeff3 = (2.0 / 315.0) * epsilon * powint(sigma, 12); + coeff4 = (1.0 / 3.0) * epsilon * powint(sigma, 6); + coeff5 = (4.0 / 315.0) * epsilon * powint(sigma, 12); + coeff6 = (1.0 / 12.0) * epsilon * powint(sigma, 6); offset = 0; - - if (utils::strmatch(update->integrate_style,"^respa")) - nlevels_respa = (dynamic_cast( update->integrate))->nlevels; + if (utils::strmatch(update->integrate_style, "^respa")) + nlevels_respa = (dynamic_cast(update->integrate))->nlevels; } /* ---------------------------------------------------------------------- */ void FixWallRegionEES::setup(int vflag) { - if (utils::strmatch(update->integrate_style,"^verlet")) + if (utils::strmatch(update->integrate_style, "^respa")) { + auto respa = dynamic_cast(update->integrate); + respa->copy_flevel_f(nlevels_respa - 1); + post_force_respa(vflag, nlevels_respa - 1, 0); + respa->copy_f_flevel(nlevels_respa - 1); + } else { post_force(vflag); - else { - (dynamic_cast( update->integrate))->copy_flevel_f(nlevels_respa-1); - post_force_respa(vflag,nlevels_respa-1,0); - (dynamic_cast( update->integrate))->copy_f_flevel(nlevels_respa-1); } } @@ -139,7 +137,7 @@ void FixWallRegionEES::setup(int vflag) void FixWallRegionEES::min_setup(int vflag) { - post_force(vflag); + post_force(vflag); } /* ---------------------------------------------------------------------- */ @@ -149,8 +147,8 @@ void FixWallRegionEES::post_force(int /*vflag*/) //sth is needed here, but I dont know what //that is calculation of sn - int i,m,n; - double rinv,fx,fy,fz,sn,tooclose[3]; + int i, m, n; + double rinv, fx, fy, fz, sn, tooclose[3]; eflag = 0; ewall[0] = ewall[1] = ewall[2] = ewall[3] = 0.0; @@ -165,7 +163,6 @@ void FixWallRegionEES::post_force(int /*vflag*/) int *mask = atom->mask; int nlocal = atom->nlocal; - Region *region = domain->regions[iregion]; region->prematch(); int onflag = 0; @@ -176,33 +173,34 @@ void FixWallRegionEES::post_force(int /*vflag*/) for (i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - if (!region->match(x[i][0],x[i][1],x[i][2])) { + if (!region->match(x[i][0], x[i][1], x[i][2])) { onflag = 1; continue; } - double A[3][3] = {{0,0,0},{0,0,0},{0,0,0}}; - double tempvec[3]= {0,0,0}; + double A[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}; + double tempvec[3] = {0, 0, 0}; double sn2 = 0.0; - double nhat[3] = {0,0,0}; - double* shape = bonus[ellipsoid[i]].shape;; - MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A); + double nhat[3] = {0, 0, 0}; + double *shape = bonus[ellipsoid[i]].shape; + ; + MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat, A); - for (int which = 0 ; which < 3; which ++) {//me - nhat[which]=1; - nhat[(which+1)%3] = 0 ; - nhat[(which+2)%3] = 0 ; - sn2 = 0 ; - MathExtra::transpose_matvec(A,nhat,tempvec); - for (int k = 0; k<3; k++) { + for (int which = 0; which < 3; which++) { //me + nhat[which] = 1; + nhat[(which + 1) % 3] = 0; + nhat[(which + 2) % 3] = 0; + sn2 = 0; + MathExtra::transpose_matvec(A, nhat, tempvec); + for (int k = 0; k < 3; k++) { tempvec[k] *= shape[k]; - sn2 += tempvec[k]*tempvec[k]; + sn2 += tempvec[k] * tempvec[k]; } sn = sqrt(sn2); tooclose[which] = sn; } - n = region->surface(x[i][0],x[i][1],x[i][2],cutoff); + n = region->surface(x[i][0], x[i][1], x[i][2], cutoff); for (m = 0; m < n; m++) { @@ -212,12 +210,13 @@ void FixWallRegionEES::post_force(int /*vflag*/) } else if (region->contact[m].dely != 0 && region->contact[m].r <= tooclose[1]) { onflag = 1; continue; - } else if (region->contact[m].delz !=0 && region->contact[m].r <= tooclose[2]) { + } else if (region->contact[m].delz != 0 && region->contact[m].r <= tooclose[2]) { onflag = 1; continue; - } else rinv = 1.0/region->contact[m].r; + } else + rinv = 1.0 / region->contact[m].r; - ees(m,i); + ees(m, i); ewall[0] += eng; fx = fwall * region->contact[m].delx * rinv; @@ -237,15 +236,15 @@ void FixWallRegionEES::post_force(int /*vflag*/) } } - if (onflag) error->one(FLERR,"Particle on or inside surface of region " - "used in fix wall/region/ees"); + if (onflag) + error->one(FLERR, "Particle on or inside surface of region used in fix wall/region/ees"); } /* ---------------------------------------------------------------------- */ void FixWallRegionEES::post_force_respa(int vflag, int ilevel, int /*iloop*/) { - if (ilevel == nlevels_respa-1) post_force(vflag); + if (ilevel == nlevels_respa - 1) post_force(vflag); } /* ---------------------------------------------------------------------- */ @@ -264,7 +263,7 @@ double FixWallRegionEES::compute_scalar() // only sum across procs one time if (eflag == 0) { - MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(ewall, ewall_all, 4, MPI_DOUBLE, MPI_SUM, world); eflag = 1; } return ewall_all[0]; @@ -279,10 +278,10 @@ double FixWallRegionEES::compute_vector(int n) // only sum across procs one time if (eflag == 0) { - MPI_Allreduce(ewall,ewall_all,4,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(ewall, ewall_all, 4, MPI_DOUBLE, MPI_SUM, world); eflag = 1; } - return ewall_all[n+1]; + return ewall_all[n + 1]; } /* ---------------------------------------------------------------------- @@ -292,24 +291,23 @@ double FixWallRegionEES::compute_vector(int n) void FixWallRegionEES::ees(int m, int i) { - Region *region = domain->regions[iregion]; region->prematch(); double delta, delta2, delta3, delta4, delta5, delta6; - double sigman, sigman2 , sigman3, sigman4, sigman5, sigman6; - double hhss, hhss2, hhss4, hhss7, hhss8; //h^2 - s_n^2 - double hps; //h+s_n - double hms; //h-s_n + double sigman, sigman2, sigman3, sigman4, sigman5, sigman6; + double hhss, hhss2, hhss4, hhss7, hhss8; //h^2 - s_n^2 + double hps; //h+s_n + double hms; //h-s_n double twall; double A[3][3], nhat[3], SAn[3], that[3]; - double tempvec[3]= {0,0,0}; - double tempvec2[3]= {0,0,0}; + double tempvec[3] = {0, 0, 0}; + double tempvec2[3] = {0, 0, 0}; - double Lx[3][3] = {{0,0,0},{0,0,-1},{0,1,0}}; - double Ly[3][3] = {{0,0,1},{0,0,0},{-1,0,0}}; - double Lz[3][3] = {{0,-1,0},{1,0,0},{0,0,0}}; + double Lx[3][3] = {{0, 0, 0}, {0, 0, -1}, {0, 1, 0}}; + double Ly[3][3] = {{0, 0, 1}, {0, 0, 0}, {-1, 0, 0}}; + double Lz[3][3] = {{0, -1, 0}, {1, 0, 0}, {0, 0, 0}}; nhat[0] = region->contact[m].delx / region->contact[m].r; nhat[1] = region->contact[m].dely / region->contact[m].r; @@ -318,14 +316,15 @@ void FixWallRegionEES::ees(int m, int i) AtomVecEllipsoid::Bonus *bonus = avec->bonus; int *ellipsoid = atom->ellipsoid; - double* shape = bonus[ellipsoid[i]].shape;; - MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat,A); + double *shape = bonus[ellipsoid[i]].shape; + ; + MathExtra::quat_to_mat(bonus[ellipsoid[i]].quat, A); sigman2 = 0.0; - MathExtra::transpose_matvec(A,nhat,tempvec); - for (int k = 0; k<3; k++) { + MathExtra::transpose_matvec(A, nhat, tempvec); + for (int k = 0; k < 3; k++) { tempvec[k] *= shape[k]; - sigman2 += tempvec[k]*tempvec[k]; + sigman2 += tempvec[k] * tempvec[k]; SAn[k] = tempvec[k]; } @@ -337,14 +336,14 @@ void FixWallRegionEES::ees(int m, int i) sigman5 = sigman4 * sigman; sigman6 = sigman3 * sigman3; - delta2 = delta * delta; + delta2 = delta * delta; delta3 = delta2 * delta; delta4 = delta2 * delta2; delta5 = delta3 * delta2; delta6 = delta3 * delta3; hhss = delta2 - sigman2; - hhss2 = hhss * hhss; + hhss2 = hhss * hhss; hhss4 = hhss2 * hhss2; hhss8 = hhss4 * hhss4; hhss7 = hhss4 * hhss2 * hhss; @@ -352,31 +351,31 @@ void FixWallRegionEES::ees(int m, int i) hps = delta + sigman; hms = delta - sigman; - fwall = -1*coeff4/hhss2 + coeff3 - * (21*delta6 + 63*delta4*sigman2 + 27*delta2*sigman4 + sigman6) / hhss8; + fwall = -1 * coeff4 / hhss2 + + coeff3 * (21 * delta6 + 63 * delta4 * sigman2 + 27 * delta2 * sigman4 + sigman6) / hhss8; - eng = -1*coeff2 * (4*delta/sigman2/hhss + 2*log(hms/hps)/sigman3) + - coeff1 * (35*delta5 + 70*delta3*sigman2 + 15*delta*sigman4) / hhss7; + eng = -1 * coeff2 * (4 * delta / sigman2 / hhss + 2 * log(hms / hps) / sigman3) + + coeff1 * (35 * delta5 + 70 * delta3 * sigman2 + 15 * delta * sigman4) / hhss7; - twall = coeff6 * (6*delta3/sigman4/hhss2 - 10*delta/sigman2/hhss2 - + 3*log(hms/hps)/sigman5) - + coeff5 * (21.*delta5 + 30.*delta3*sigman2 + 5.*delta*sigman4) / hhss8; + twall = coeff6 * + (6 * delta3 / sigman4 / hhss2 - 10 * delta / sigman2 / hhss2 + + 3 * log(hms / hps) / sigman5) + + coeff5 * (21. * delta5 + 30. * delta3 * sigman2 + 5. * delta * sigman4) / hhss8; - MathExtra::matvec(Lx,nhat,tempvec); - MathExtra::transpose_matvec(A,tempvec,tempvec2); - for (int k = 0; k<3; k++) tempvec2[k] *= shape[k]; - that[0] = MathExtra::dot3(SAn,tempvec2); - - MathExtra::matvec(Ly,nhat,tempvec); - MathExtra::transpose_matvec(A,tempvec,tempvec2); - for (int k = 0; k<3; k++) tempvec2[k] *= shape[k]; - that[1] = MathExtra::dot3(SAn,tempvec2); - - MathExtra::matvec(Lz,nhat,tempvec); - MathExtra::transpose_matvec(A,tempvec,tempvec2); + MathExtra::matvec(Lx, nhat, tempvec); + MathExtra::transpose_matvec(A, tempvec, tempvec2); for (int k = 0; k < 3; k++) tempvec2[k] *= shape[k]; - that[2] = MathExtra::dot3(SAn,tempvec2); + that[0] = MathExtra::dot3(SAn, tempvec2); - for (int j = 0; j<3 ; j++) - torque[j] = twall * that[j]; + MathExtra::matvec(Ly, nhat, tempvec); + MathExtra::transpose_matvec(A, tempvec, tempvec2); + for (int k = 0; k < 3; k++) tempvec2[k] *= shape[k]; + that[1] = MathExtra::dot3(SAn, tempvec2); + + MathExtra::matvec(Lz, nhat, tempvec); + MathExtra::transpose_matvec(A, tempvec, tempvec2); + for (int k = 0; k < 3; k++) tempvec2[k] *= shape[k]; + that[2] = MathExtra::dot3(SAn, tempvec2); + + for (int j = 0; j < 3; j++) torque[j] = twall * that[j]; } diff --git a/src/EXTRA-FIX/fix_wall_region_ees.h b/src/EXTRA-FIX/fix_wall_region_ees.h index 5163d99e90..ee5932e959 100644 --- a/src/EXTRA-FIX/fix_wall_region_ees.h +++ b/src/EXTRA-FIX/fix_wall_region_ees.h @@ -41,12 +41,12 @@ class FixWallRegionEES : public Fix { private: class AtomVecEllipsoid *avec; - int iregion; double epsilon, sigma, cutoff; int eflag; double ewall[4], ewall_all[4]; int nlevels_respa; char *idregion; + class Region *region; double coeff1, coeff2, coeff3, coeff4, offset; double coeff5, coeff6; diff --git a/src/GRANULAR/fix_pour.cpp b/src/GRANULAR/fix_pour.cpp index 2a25efc6f2..7a8bf630df 100644 --- a/src/GRANULAR/fix_pour.cpp +++ b/src/GRANULAR/fix_pour.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -37,106 +36,101 @@ using namespace LAMMPS_NS; using namespace FixConst; -using namespace MathConst; +using MathConst::MY_2PI; +using MathConst::MY_4PI3; +using MathConst::MY_PI; -enum{ATOM,MOLECULE}; -enum{ONE,RANGE,POLY}; +enum { ATOM, MOLECULE }; +enum { ONE, RANGE, POLY }; -#define EPSILON 0.001 -#define SMALL 1.0e-10 +static constexpr double EPSILON = 0.001; +static constexpr double SMALL = 1.0e-10; /* ---------------------------------------------------------------------- */ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), radius_poly(nullptr), frac_poly(nullptr), - idrigid(nullptr), idshake(nullptr), onemols(nullptr), molfrac(nullptr), coords(nullptr), - imageflags(nullptr), fixrigid(nullptr), fixshake(nullptr), recvcounts(nullptr), - displs(nullptr), random(nullptr), random2(nullptr) + Fix(lmp, narg, arg), radius_poly(nullptr), frac_poly(nullptr), idrigid(nullptr), + idshake(nullptr), idregion(nullptr), region(nullptr), onemols(nullptr), molfrac(nullptr), + coords(nullptr), imageflags(nullptr), fixrigid(nullptr), fixshake(nullptr), recvcounts(nullptr), + displs(nullptr), random(nullptr), random2(nullptr) { - if (narg < 6) error->all(FLERR,"Illegal fix pour command"); + if (narg < 6) error->all(FLERR, "Illegal fix pour command"); - if (lmp->kokkos) - error->all(FLERR,"Cannot yet use fix pour with the KOKKOS package"); + if (lmp->kokkos) error->all(FLERR, "Cannot yet use fix pour with the KOKKOS package"); time_depend = 1; if (!atom->radius_flag || !atom->rmass_flag) - error->all(FLERR,"Fix pour requires atom attributes radius, rmass"); + error->all(FLERR, "Fix pour requires atom attributes radius, rmass"); // required args - ninsert = utils::inumeric(FLERR,arg[3],false,lmp); - ntype = utils::inumeric(FLERR,arg[4],false,lmp); - seed = utils::inumeric(FLERR,arg[5],false,lmp); + ninsert = utils::inumeric(FLERR, arg[3], false, lmp); + ntype = utils::inumeric(FLERR, arg[4], false, lmp); + seed = utils::inumeric(FLERR, arg[5], false, lmp); - if (seed <= 0) error->all(FLERR,"Illegal fix pour command"); + if (seed <= 0) error->all(FLERR, "Illegal fix pour command"); // read options from end of input line - options(narg-6,&arg[6]); + options(narg - 6, &arg[6]); // error check on type if (mode == ATOM && (ntype <= 0 || ntype > atom->ntypes)) - error->all(FLERR,"Invalid atom type in fix pour command"); + error->all(FLERR, "Invalid atom type in fix pour command"); // error checks on region and its extent being inside simulation box - if (iregion == -1) error->all(FLERR,"Must specify a region in fix pour"); - if (domain->regions[iregion]->bboxflag == 0) - error->all(FLERR,"Fix pour region does not support a bounding box"); - if (domain->regions[iregion]->dynamic_check()) - error->all(FLERR,"Fix pour region cannot be dynamic"); + if (!region) error->all(FLERR, "Must specify a region in fix pour"); + if (region->bboxflag == 0) + error->all(FLERR, "Fix pour region {} does not support a bounding box", idregion); + if (region->dynamic_check()) error->all(FLERR, "Fix pour region {} cannot be dynamic", idregion); - if (strcmp(domain->regions[iregion]->style,"block") == 0) { + if (strcmp(region->style, "block") == 0) { + auto block = dynamic_cast(region); region_style = 1; - xlo = (dynamic_cast( domain->regions[iregion]))->xlo; - xhi = (dynamic_cast( domain->regions[iregion]))->xhi; - ylo = (dynamic_cast( domain->regions[iregion]))->ylo; - yhi = (dynamic_cast( domain->regions[iregion]))->yhi; - zlo = (dynamic_cast( domain->regions[iregion]))->zlo; - zhi = (dynamic_cast( domain->regions[iregion]))->zhi; - if (xlo < domain->boxlo[0] || xhi > domain->boxhi[0] || - ylo < domain->boxlo[1] || yhi > domain->boxhi[1] || - zlo < domain->boxlo[2] || zhi > domain->boxhi[2]) - error->all(FLERR,"Insertion region extends outside simulation box"); - } else if (strcmp(domain->regions[iregion]->style,"cylinder") == 0) { + xlo = block->xlo; + xhi = block->xhi; + ylo = block->ylo; + yhi = block->yhi; + zlo = block->zlo; + zhi = block->zhi; + if (xlo < domain->boxlo[0] || xhi > domain->boxhi[0] || ylo < domain->boxlo[1] || + yhi > domain->boxhi[1] || zlo < domain->boxlo[2] || zhi > domain->boxhi[2]) + error->all(FLERR, "Insertion region extends outside simulation box"); + } else if (strcmp(region->style, "cylinder") == 0) { + auto cylinder = dynamic_cast(region); region_style = 2; - char axis = (dynamic_cast( domain->regions[iregion]))->axis; - xc = (dynamic_cast( domain->regions[iregion]))->c1; - yc = (dynamic_cast( domain->regions[iregion]))->c2; - rc = (dynamic_cast( domain->regions[iregion]))->radius; - zlo = (dynamic_cast( domain->regions[iregion]))->lo; - zhi = (dynamic_cast( domain->regions[iregion]))->hi; - if (axis != 'z') - error->all(FLERR,"Must use a z-axis cylinder region with fix pour"); - if (xc-rc < domain->boxlo[0] || xc+rc > domain->boxhi[0] || - yc-rc < domain->boxlo[1] || yc+rc > domain->boxhi[1] || - zlo < domain->boxlo[2] || zhi > domain->boxhi[2]) - error->all(FLERR,"Insertion region extends outside simulation box"); - } else error->all(FLERR,"Must use a block or cylinder region with fix pour"); + char axis = cylinder->axis; + xc = cylinder->c1; + yc = cylinder->c2; + rc = cylinder->radius; + zlo = cylinder->lo; + zhi = cylinder->hi; + if (axis != 'z') error->all(FLERR, "Must use a z-axis cylinder region with fix pour"); + if (xc - rc < domain->boxlo[0] || xc + rc > domain->boxhi[0] || yc - rc < domain->boxlo[1] || + yc + rc > domain->boxhi[1] || zlo < domain->boxlo[2] || zhi > domain->boxhi[2]) + error->all(FLERR, "Insertion region extends outside simulation box"); + } else + error->all(FLERR, "Must use a block or cylinder region with fix pour"); if (region_style == 2 && domain->dimension == 2) - error->all(FLERR, - "Must use a block region with fix pour for 2d simulations"); + error->all(FLERR, "Must use a block region with fix pour for 2d simulations"); // error check and further setup for mode = MOLECULE - if (atom->tag_enable == 0) - error->all(FLERR,"Cannot use fix_pour unless atoms have IDs"); + if (atom->tag_enable == 0) error->all(FLERR, "Cannot use fix_pour unless atoms have IDs"); if (mode == MOLECULE) { for (int i = 0; i < nmol; i++) { - if (onemols[i]->xflag == 0) - error->all(FLERR,"Fix pour molecule must have coordinates"); - if (onemols[i]->typeflag == 0) - error->all(FLERR,"Fix pour molecule must have atom types"); - if (ntype+onemols[i]->ntypes <= 0 || - ntype+onemols[i]->ntypes > atom->ntypes) - error->all(FLERR,"Invalid atom type in fix pour mol command"); + if (onemols[i]->xflag == 0) error->all(FLERR, "Fix pour molecule must have coordinates"); + if (onemols[i]->typeflag == 0) error->all(FLERR, "Fix pour molecule must have atom types"); + if (ntype + onemols[i]->ntypes <= 0 || ntype + onemols[i]->ntypes > atom->ntypes) + error->all(FLERR, "Invalid atom type in fix pour mol command"); if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols) - error->all(FLERR,"Fix pour molecule template ID must be same as atom style template ID"); + error->all(FLERR, "Fix pour molecule template ID must be same as atom style template ID"); onemols[i]->check_attributes(0); // fix pour uses geoemetric center of molecule for insertion @@ -145,23 +139,20 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : } } - if (rigidflag && mode == ATOM) - error->all(FLERR,"Cannot use fix pour rigid and not molecule"); - if (shakeflag && mode == ATOM) - error->all(FLERR,"Cannot use fix pour shake and not molecule"); - if (rigidflag && shakeflag) - error->all(FLERR,"Cannot use fix pour rigid and shake"); + if (rigidflag && mode == ATOM) error->all(FLERR, "Cannot use fix pour rigid and not molecule"); + if (shakeflag && mode == ATOM) error->all(FLERR, "Cannot use fix pour shake and not molecule"); + if (rigidflag && shakeflag) error->all(FLERR, "Cannot use fix pour rigid and shake"); // setup of coords and imageflags array - if (mode == ATOM) natom_max = 1; + if (mode == ATOM) + natom_max = 1; else { natom_max = 0; - for (int i = 0; i < nmol; i++) - natom_max = MAX(natom_max,onemols[i]->natoms); + for (int i = 0; i < nmol; i++) natom_max = MAX(natom_max, onemols[i]->natoms); } - memory->create(coords,natom_max,4,"pour:coords"); - memory->create(imageflags,natom_max,"pour:imageflags"); + memory->create(coords, natom_max, 4, "pour:coords"); + memory->create(imageflags, natom_max, "pour:imageflags"); // find max atom and molecule IDs just once @@ -171,13 +162,13 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : // warm up the generator 30x to avoid correlations in first-particle // positions if runs are repeated with consecutive seeds - random = new RanPark(lmp,seed); - for (int ii=0; ii < 30; ii++) random->uniform(); + random = new RanPark(lmp, seed); + for (int ii = 0; ii < 30; ii++) random->uniform(); // allgather arrays - MPI_Comm_rank(world,&me); - MPI_Comm_size(world,&nprocs); + MPI_Comm_rank(world, &me); + MPI_Comm_size(world, &nprocs); recvcounts = new int[nprocs]; displs = new int[nprocs]; @@ -186,7 +177,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : 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"); + error->all(FLERR, "There must be exactly one fix gravity defined for fix pour"); auto fixgrav = dynamic_cast(fixlist.front()); grav = -fixgrav->magnitude * force->ftm2v; @@ -200,7 +191,7 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : // gives t = [-(vz-rate) - sqrt((vz-rate)^2 - 2*grav*(zhi-zlo))] / grav // where zhi-zlo > 0, grav < 0, and vz & rate can be either > 0 or < 0 - double v_relative,delta; + double v_relative, delta; if (domain->dimension == 3) { v_relative = vz - rate; delta = zhi - zlo; @@ -208,8 +199,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 @@ -224,52 +215,51 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : // volume_one = volume of inserted particle (with max possible radius) // in 3d, insure dy >= 1, for quasi-2d simulations - double volume,volume_one=1.0; + double volume, volume_one = 1.0; molradius_max = 0.0; if (mode == MOLECULE) { - for (int i = 0; i < nmol; i++) - molradius_max = MAX(molradius_max,onemols[i]->molradius); + for (int i = 0; i < nmol; i++) molradius_max = MAX(molradius_max, onemols[i]->molradius); } if (domain->dimension == 3) { if (region_style == 1) { double dy = yhi - ylo; if (dy < 1.0) dy = 1.0; - volume = (xhi-xlo) * dy * (zhi-zlo); - } else volume = MY_PI*rc*rc * (zhi-zlo); + volume = (xhi - xlo) * dy * (zhi - zlo); + } else + volume = MY_PI * rc * rc * (zhi - zlo); if (mode == MOLECULE) { - volume_one = 4.0/3.0 * MY_PI * molradius_max*molradius_max*molradius_max; + volume_one = MY_4PI3 * molradius_max * molradius_max * molradius_max; } else if (dstyle == ONE || dstyle == RANGE) { - volume_one = 4.0/3.0 * MY_PI * radius_max*radius_max*radius_max; + volume_one = MY_4PI3 * radius_max * radius_max * radius_max; } else if (dstyle == POLY) { volume_one = 0.0; for (int i = 0; i < npoly; i++) - volume_one += (4.0/3.0 * MY_PI * - radius_poly[i]*radius_poly[i]*radius_poly[i]) * frac_poly[i]; + volume_one += (MY_4PI3 * radius_poly[i] * radius_poly[i] * radius_poly[i]) * frac_poly[i]; } } else { - volume = (xhi-xlo) * (yhi-ylo); + volume = (xhi - xlo) * (yhi - ylo); if (mode == MOLECULE) { - volume_one = MY_PI * molradius_max*molradius_max; + volume_one = MY_PI * molradius_max * molradius_max; } else if (dstyle == ONE || dstyle == RANGE) { - volume_one = MY_PI * radius_max*radius_max; + volume_one = MY_PI * radius_max * radius_max; } else if (dstyle == POLY) { volume_one = 0.0; for (int i = 0; i < npoly; i++) - volume_one += (MY_PI * radius_poly[i]*radius_poly[i]) * frac_poly[i]; + volume_one += (MY_PI * radius_poly[i] * radius_poly[i]) * frac_poly[i]; } } - nper = static_cast (volfrac*volume/volume_one); - if (nper == 0) error->all(FLERR,"Fix pour insertion count per timestep is 0"); - int nfinal = update->ntimestep + 1 + ((bigint)ninsert-1)/nper * nfreq; + nper = static_cast(volfrac * volume / volume_one); + if (nper == 0) error->all(FLERR, "Fix pour insertion count per timestep is 0"); + int nfinal = update->ntimestep + 1 + ((bigint) ninsert - 1) / nper * nfreq; // print stats if (me == 0) - utils::logmesg(lmp, "Particle insertion: {} every {} steps, {} by step {}\n", - nper,nfreq,ninsert,nfinal); + utils::logmesg(lmp, "Particle insertion: {} every {} steps, {} by step {}\n", nper, nfreq, + ninsert, nfinal); } /* ---------------------------------------------------------------------- */ @@ -277,15 +267,15 @@ FixPour::FixPour(LAMMPS *lmp, int narg, char **arg) : FixPour::~FixPour() { delete random; - delete [] molfrac; - delete [] idrigid; - delete [] idshake; - delete [] radius_poly; - delete [] frac_poly; + delete[] molfrac; + delete[] idrigid; + delete[] idshake; + delete[] radius_poly; + delete[] frac_poly; memory->destroy(coords); memory->destroy(imageflags); - delete [] recvcounts; - delete [] displs; + delete[] recvcounts; + delete[] displs; } /* ---------------------------------------------------------------------- */ @@ -301,8 +291,10 @@ int FixPour::setmask() void FixPour::init() { - if (domain->triclinic) - error->all(FLERR,"Cannot use fix pour with triclinic box"); + if (domain->triclinic) error->all(FLERR, "Cannot use fix pour with triclinic box"); + + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Fix pour region {} does not exist", idregion); // insure gravity fix (still) exists // for 3d must point in -z, for 2d must point in -y @@ -310,37 +302,35 @@ void FixPour::init() 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"); + error->all(FLERR, "There must be exactly one fix gravity defined for fix pour"); auto fixgrav = dynamic_cast(fixlist.front()); if (fixgrav->varflag != FixGravity::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 = fixgrav->xgrav; double ygrav = fixgrav->ygrav; double zgrav = fixgrav->zgrav; if (domain->dimension == 3) { - if (fabs(xgrav) > EPSILON || fabs(ygrav) > EPSILON || - fabs(zgrav+1.0) > EPSILON) - error->all(FLERR,"Gravity must point in -z to use with fix pour in 3d"); + if (fabs(xgrav) > EPSILON || fabs(ygrav) > EPSILON || fabs(zgrav + 1.0) > EPSILON) + error->all(FLERR, "Gravity must point in -z to use with fix pour in 3d"); } else { - if (fabs(xgrav) > EPSILON || fabs(ygrav+1.0) > EPSILON || - fabs(zgrav) > EPSILON) - error->all(FLERR,"Gravity must point in -y to use with fix pour in 2d"); + if (fabs(xgrav) > EPSILON || fabs(ygrav + 1.0) > EPSILON || fabs(zgrav) > EPSILON) + error->all(FLERR, "Gravity must point in -y to use with fix pour in 2d"); } double gnew = -fixgrav->magnitude * force->ftm2v; - if (gnew != grav) error->all(FLERR,"Gravity changed since fix pour was created"); + 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 if (rigidflag) { fixrigid = modify->get_fix_by_id(idrigid); - if (!fixrigid) error->all(FLERR,"Fix pour rigid fix does not exist"); + if (!fixrigid) error->all(FLERR, "Fix pour rigid fix {} does not exist", idrigid); int tmp; - if (onemols != (Molecule **) fixrigid->extract("onemol",tmp)) - error->all(FLERR,"Fix pour and fix rigid/small not using same molecule template ID"); + if (onemols != (Molecule **) fixrigid->extract("onemol", tmp)) + error->all(FLERR, "Fix pour and fix rigid/small not using same molecule template ID"); } // if shakeflag defined, check for SHAKE fix @@ -348,10 +338,10 @@ void FixPour::init() if (shakeflag) { fixshake = modify->get_fix_by_id(idshake); - if (!fixshake) error->all(FLERR,"Fix pour shake fix does not exist"); + if (!fixshake) error->all(FLERR, "Fix pour shake fix {} does not exist", idshake); int tmp; - if (onemols != (Molecule **) fixshake->extract("onemol",tmp)) - error->all(FLERR,"Fix pour and fix shake not using same molecule template ID"); + if (onemols != (Molecule **) fixshake->extract("onemol", tmp)) + error->all(FLERR, "Fix pour and fix shake not using same molecule template ID"); } } @@ -359,8 +349,10 @@ void FixPour::init() void FixPour::setup_pre_exchange() { - if (ninserted < ninsert) next_reneighbor = update->ntimestep + 1; - else next_reneighbor = 0; + if (ninserted < ninsert) + next_reneighbor = update->ntimestep + 1; + else + next_reneighbor = 0; } /* ---------------------------------------------------------------------- @@ -369,8 +361,8 @@ void FixPour::setup_pre_exchange() void FixPour::pre_exchange() { - int i,m,flag,nlocalprev,imol,natom; - double r[3],rotmat[3][3],quat[4],vnew[3]; + int i, m, flag, nlocalprev, imol, natom; + double r[3], rotmat[3][3], quat[4], vnew[3]; double *newcoord; // just return if should not be called on this timestep @@ -414,24 +406,24 @@ void FixPour::pre_exchange() if (overlap(i)) ncount++; int nprevious; - MPI_Allreduce(&ncount,&nprevious,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&ncount, &nprevious, 1, MPI_INT, MPI_SUM, world); // xmine is for my atoms // xnear is for atoms from all procs + atoms to be inserted - double **xmine,**xnear; - memory->create(xmine,ncount,4,"fix_pour:xmine"); - memory->create(xnear,nprevious+nnew*natom_max,4,"fix_pour:xnear"); + double **xmine, **xnear; + memory->create(xmine, ncount, 4, "fix_pour:xmine"); + memory->create(xnear, nprevious + nnew * natom_max, 4, "fix_pour:xnear"); int nnear = nprevious; // setup for allgatherv - int n = 4*ncount; - MPI_Allgather(&n,1,MPI_INT,recvcounts,1,MPI_INT,world); + int n = 4 * ncount; + MPI_Allgather(&n, 1, MPI_INT, recvcounts, 1, MPI_INT, world); displs[0] = 0; for (int iproc = 1; iproc < nprocs; iproc++) - displs[iproc] = displs[iproc-1] + recvcounts[iproc-1]; + displs[iproc] = displs[iproc - 1] + recvcounts[iproc - 1]; // load up xmine array @@ -452,8 +444,7 @@ void FixPour::pre_exchange() double *ptr = nullptr; if (ncount) ptr = xmine[0]; - MPI_Allgatherv(ptr,4*ncount,MPI_DOUBLE, - xnear[0],recvcounts,displs,MPI_DOUBLE,world); + MPI_Allgatherv(ptr, 4 * ncount, MPI_DOUBLE, xnear[0], recvcounts, displs, MPI_DOUBLE, world); // insert new particles into xnear list, one by one // check against all nearby atoms and previously inserted ones @@ -468,7 +459,7 @@ void FixPour::pre_exchange() // store image flag modified due to PBC int success; - double radtmp,delx,dely,delz,rsq,radsum,rn,h; + double radtmp, delx, dely, delz, rsq, radsum, rn, h; double coord[3]; double denstmp; @@ -481,13 +472,13 @@ void FixPour::pre_exchange() while (nsuccess < nnew) { rn = random->uniform(); - h = hi_current - rn*rn * (hi_current-lo_current); + h = hi_current - rn * rn * (hi_current - lo_current); if (mode == ATOM) radtmp = radius_sample(); success = 0; while (attempt < maxiter) { attempt++; - xyz_random(h,coord); + xyz_random(h, coord); if (mode == ATOM) { natom = 1; @@ -495,8 +486,7 @@ void FixPour::pre_exchange() coords[0][1] = coord[1]; coords[0][2] = coord[2]; coords[0][3] = radtmp; - imageflags[0] = ((imageint) IMGMAX << IMG2BITS) | - ((imageint) IMGMAX << IMGBITS) | IMGMAX; + imageflags[0] = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; } else { double rng = random->uniform(); imol = 0; @@ -512,10 +502,10 @@ void FixPour::pre_exchange() } double theta = random->uniform() * MY_2PI; MathExtra::norm3(r); - MathExtra::axisangle_to_quat(r,theta,quat); - MathExtra::quat_to_mat(quat,rotmat); + MathExtra::axisangle_to_quat(r, theta, quat); + MathExtra::quat_to_mat(quat, rotmat); for (i = 0; i < natom; i++) { - MathExtra::matvec(rotmat,onemols[imol]->dx[i],coords[i]); + MathExtra::matvec(rotmat, onemols[imol]->dx[i], coords[i]); coords[i][0] += coord[0]; coords[i][1] += coord[1]; coords[i][2] += coord[2]; @@ -526,11 +516,11 @@ void FixPour::pre_exchange() if (onemols[imol]->radiusflag) coords[i][3] = onemols[imol]->radius[i]; - else coords[i][3] = 0.5; + else + coords[i][3] = 0.5; - imageflags[i] = ((imageint) IMGMAX << IMG2BITS) | - ((imageint) IMGMAX << IMGBITS) | IMGMAX; - domain->remap(coords[i],imageflags[i]); + imageflags[i] = ((imageint) IMGMAX << IMG2BITS) | ((imageint) IMGMAX << IMGBITS) | IMGMAX; + domain->remap(coords[i], imageflags[i]); } } @@ -542,10 +532,10 @@ void FixPour::pre_exchange() delx = coords[m][0] - xnear[i][0]; dely = coords[m][1] - xnear[i][1]; delz = coords[m][2] - xnear[i][2]; - domain->minimum_image(delx,dely,delz); - rsq = delx*delx + dely*dely + delz*delz; + domain->minimum_image(delx, dely, delz); + rsq = delx * delx + dely * dely + delz * delz; radsum = coords[m][3] + xnear[i][3]; - if (rsq <= radsum*radsum) break; + if (rsq <= radsum * radsum) break; } if (i < nnear) break; } @@ -582,12 +572,12 @@ void FixPour::pre_exchange() // coord[2] = hi_current + vz*t + 1/2 grav t^2 if (dimension == 3) { - vnew[0] = vxlo + random->uniform() * (vxhi-vxlo); - vnew[1] = vylo + random->uniform() * (vyhi-vylo); - vnew[2] = -sqrt(vz*vz + 2.0*grav*(coord[2]-hi_current)); + vnew[0] = vxlo + random->uniform() * (vxhi - vxlo); + vnew[1] = vylo + random->uniform() * (vyhi - vylo); + vnew[2] = -sqrt(vz * vz + 2.0 * grav * (coord[2] - hi_current)); } else { - vnew[0] = vxlo + random->uniform() * (vxhi-vxlo); - vnew[1] = -sqrt(vy*vy + 2.0*grav*(coord[1]-hi_current)); + vnew[0] = vxlo + random->uniform() * (vxhi - vxlo); + vnew[1] = -sqrt(vy * vy + 2.0 * grav * (coord[1] - hi_current)); vnew[2] = 0.0; } @@ -597,45 +587,47 @@ void FixPour::pre_exchange() // set group mask to "all" plus fix group for (m = 0; m < natom; m++) { - if (mode == ATOM) - denstmp = density_lo + random->uniform() * (density_hi-density_lo); + if (mode == ATOM) denstmp = density_lo + random->uniform() * (density_hi - density_lo); newcoord = coords[m]; flag = 0; - if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && - newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] && - newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1; + if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && newcoord[1] >= sublo[1] && + newcoord[1] < subhi[1] && newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) + flag = 1; else if (dimension == 3 && newcoord[2] >= domain->boxhi[2]) { if (comm->layout != Comm::LAYOUT_TILED) { - if (comm->myloc[2] == comm->procgrid[2]-1 && - newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && - newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; + if (comm->myloc[2] == comm->procgrid[2] - 1 && newcoord[0] >= sublo[0] && + newcoord[0] < subhi[0] && newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) + flag = 1; } else { - if (comm->mysplit[2][1] == 1.0 && - newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && - newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1; + if (comm->mysplit[2][1] == 1.0 && newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] && + newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) + flag = 1; } } else if (dimension == 2 && newcoord[1] >= domain->boxhi[1]) { if (comm->layout != Comm::LAYOUT_TILED) { - if (comm->myloc[1] == comm->procgrid[1]-1 && - newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; + if (comm->myloc[1] == comm->procgrid[1] - 1 && newcoord[0] >= sublo[0] && + newcoord[0] < subhi[0]) + flag = 1; } else { - if (comm->mysplit[1][1] == 1.0 && - newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1; + if (comm->mysplit[1][1] == 1.0 && newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) + flag = 1; } } if (flag) { - if (mode == ATOM) atom->avec->create_atom(ntype,coords[m]); - else atom->avec->create_atom(ntype+onemols[imol]->type[m],coords[m]); + if (mode == ATOM) + atom->avec->create_atom(ntype, coords[m]); + else + atom->avec->create_atom(ntype + onemols[imol]->type[m], coords[m]); int n = atom->nlocal - 1; - atom->tag[n] = maxtag_all + m+1; + atom->tag[n] = maxtag_all + m + 1; if (mode == MOLECULE) { if (atom->molecule_flag) { if (onemols[imol]->moleculeflag) { atom->molecule[n] = maxmol_all + onemols[imol]->molecule[m]; } else { - atom->molecule[n] = maxmol_all+1; + atom->molecule[n] = maxmol_all + 1; } } if (atom->molecular == Atom::TEMPLATE) { @@ -651,10 +643,10 @@ void FixPour::pre_exchange() if (mode == ATOM) { radtmp = newcoord[3]; atom->radius[n] = radtmp; - atom->rmass[n] = 4.0*MY_PI/3.0 * radtmp*radtmp*radtmp * denstmp; + atom->rmass[n] = 4.0 * MY_PI / 3.0 * radtmp * radtmp * radtmp * denstmp; } else { onemols[imol]->quat_external = quat; - atom->add_molecule_atom(onemols[imol],m,n,maxtag_all); + atom->add_molecule_atom(onemols[imol], m, n, maxtag_all); } modify->create_attribute(n); @@ -666,9 +658,9 @@ void FixPour::pre_exchange() // FixShake::set_molecule stores shake info for molecule if (rigidflag) - fixrigid->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat); + fixrigid->set_molecule(nlocalprev, maxtag_all, imol, coord, vnew, quat); else if (shakeflag) - fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat); + fixshake->set_molecule(nlocalprev, maxtag_all, imol, coord, vnew, quat); maxtag_all += natom; if (mode == MOLECULE && atom->molecule_flag) { @@ -685,8 +677,7 @@ void FixPour::pre_exchange() int ninserted_atoms = nnear - nprevious; int ninserted_mols = ninserted_atoms / natom; ninserted += ninserted_mols; - if (ninserted_mols < nnew && me == 0) - error->warning(FLERR,"Less insertions than requested"); + if (ninserted_mols < nnew && me == 0) error->warning(FLERR, "Less insertions than requested"); // reset global natoms,nbonds,etc // increment maxtag_all and maxmol_all if necessary @@ -696,16 +687,14 @@ void FixPour::pre_exchange() if (ninserted_atoms) { atom->natoms += ninserted_atoms; - if (atom->natoms < 0) - error->all(FLERR,"Too many total atoms"); + if (atom->natoms < 0) error->all(FLERR, "Too many total atoms"); if (mode == MOLECULE) { - atom->nbonds += (bigint)onemols[imol]->nbonds * ninserted_mols; - atom->nangles += (bigint)onemols[imol]->nangles * ninserted_mols; - atom->ndihedrals += (bigint)onemols[imol]->ndihedrals * ninserted_mols; - atom->nimpropers += (bigint)onemols[imol]->nimpropers * ninserted_mols; + atom->nbonds += (bigint) onemols[imol]->nbonds * ninserted_mols; + atom->nangles += (bigint) onemols[imol]->nangles * ninserted_mols; + atom->ndihedrals += (bigint) onemols[imol]->ndihedrals * ninserted_mols; + atom->nimpropers += (bigint) onemols[imol]->nimpropers * ninserted_mols; } - if (maxtag_all >= MAXTAGINT) - error->all(FLERR,"New atom IDs exceed maximum allowed ID"); + if (maxtag_all >= MAXTAGINT) error->all(FLERR, "New atom IDs exceed maximum allowed ID"); } // rebuild atom map @@ -722,8 +711,10 @@ void FixPour::pre_exchange() // next timestep to insert - if (ninserted < ninsert) next_reneighbor += nfreq; - else next_reneighbor = 0; + if (ninserted < ninsert) + next_reneighbor += nfreq; + else + next_reneighbor = 0; } /* ---------------------------------------------------------------------- @@ -738,13 +729,13 @@ void FixPour::find_maxid() int nlocal = atom->nlocal; tagint max = 0; - for (int i = 0; i < nlocal; i++) max = MAX(max,tag[i]); - MPI_Allreduce(&max,&maxtag_all,1,MPI_LMP_TAGINT,MPI_MAX,world); + for (int i = 0; i < nlocal; i++) max = MAX(max, tag[i]); + MPI_Allreduce(&max, &maxtag_all, 1, MPI_LMP_TAGINT, MPI_MAX, world); if (mode == MOLECULE && molecule) { max = 0; - for (int i = 0; i < nlocal; i++) max = MAX(max,molecule[i]); - MPI_Allreduce(&max,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world); + for (int i = 0; i < nlocal; i++) max = MAX(max, molecule[i]); + MPI_Allreduce(&max, &maxmol_all, 1, MPI_LMP_TAGINT, MPI_MAX, world); } } @@ -766,29 +757,31 @@ int FixPour::overlap(int i) if (ignoretri && atom->tri[i] >= 0) return 0; } - if (mode == ATOM) delta = atom->radius[i] + radius_max; - else delta = atom->radius[i] + molradius_max; + if (mode == ATOM) + delta = atom->radius[i] + radius_max; + else + delta = atom->radius[i] + molradius_max; double *x = atom->x[i]; if (domain->dimension == 3) { if (region_style == 1) { - if (outside(0,x[0],xlo-delta,xhi+delta)) return 0; - if (outside(1,x[1],ylo-delta,yhi+delta)) return 0; - if (outside(2,x[2],lo_current-delta,hi_current+delta)) return 0; + if (outside(0, x[0], xlo - delta, xhi + delta)) return 0; + if (outside(1, x[1], ylo - delta, yhi + delta)) return 0; + if (outside(2, x[2], lo_current - delta, hi_current + delta)) return 0; } else { double delx = x[0] - xc; double dely = x[1] - yc; double delz = 0.0; - domain->minimum_image(delx,dely,delz); - double rsq = delx*delx + dely*dely; + domain->minimum_image(delx, dely, delz); + double rsq = delx * delx + dely * dely; double r = rc + delta; - if (rsq > r*r) return 0; - if (outside(2,x[2],lo_current-delta,hi_current+delta)) return 0; + if (rsq > r * r) return 0; + if (outside(2, x[2], lo_current - delta, hi_current + delta)) return 0; } } else { - if (outside(0,x[0],xlo-delta,xhi+delta)) return 0; - if (outside(1,x[1],lo_current-delta,hi_current+delta)) return 0; + if (outside(0, x[0], xlo - delta, xhi + delta)) return 0; + if (outside(1, x[1], lo_current - delta, hi_current + delta)) return 0; } return 1; @@ -836,22 +829,22 @@ void FixPour::xyz_random(double h, double *coord) { if (domain->dimension == 3) { if (region_style == 1) { - coord[0] = xlo + random->uniform() * (xhi-xlo); - coord[1] = ylo + random->uniform() * (yhi-ylo); + coord[0] = xlo + random->uniform() * (xhi - xlo); + coord[1] = ylo + random->uniform() * (yhi - ylo); coord[2] = h; } else { - double r1,r2; + double r1, r2; while (true) { r1 = random->uniform() - 0.5; r2 = random->uniform() - 0.5; - if (r1*r1 + r2*r2 < 0.25) break; + if (r1 * r1 + r2 * r2 < 0.25) break; } - coord[0] = xc + 2.0*r1*rc; - coord[1] = yc + 2.0*r2*rc; + coord[0] = xc + 2.0 * r1 * rc; + coord[1] = yc + 2.0 * r2 * rc; coord[2] = h; } } else { - coord[0] = xlo + random->uniform() * (xhi-xlo); + coord[0] = xlo + random->uniform() * (xhi - xlo); coord[1] = h; coord[2] = 0.0; } @@ -862,8 +855,7 @@ void FixPour::xyz_random(double h, double *coord) double FixPour::radius_sample() { if (dstyle == ONE) return radius_one; - if (dstyle == RANGE) return radius_lo + - random->uniform()*(radius_hi-radius_lo); + if (dstyle == RANGE) return radius_lo + random->uniform() * (radius_hi - radius_lo); double value = random->uniform(); @@ -873,7 +865,7 @@ double FixPour::radius_sample() sum += frac_poly[i]; i++; } - return radius_poly[i-1]; + return radius_poly[i - 1]; } /* ---------------------------------------------------------------------- @@ -884,18 +876,13 @@ void FixPour::options(int narg, char **arg) { // defaults - iregion = -1; mode = ATOM; - molfrac = nullptr; rigidflag = 0; - idrigid = nullptr; shakeflag = 0; - idshake = nullptr; idnext = 0; ignoreflag = ignoreline = ignoretri = 0; dstyle = ONE; radius_max = radius_one = 0.5; - radius_poly = frac_poly = nullptr; density_lo = density_hi = 1.0; volfrac = 0.25; maxattempt = 50; @@ -904,139 +891,141 @@ void FixPour::options(int narg, char **arg) int iarg = 0; while (iarg < narg) { - if (strcmp(arg[iarg],"region") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); - iregion = domain->find_region(arg[iarg+1]); - if (iregion == -1) error->all(FLERR,"Fix pour region ID does not exist"); + if (strcmp(arg[iarg], "region") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command"); + region = domain->get_region_by_id(arg[iarg + 1]); + if (!region) error->all(FLERR, "Fix pour region {} does not exist", arg[iarg + 1]); + idregion = utils::strdup(arg[iarg + 1]); iarg += 2; - - } else if (strcmp(arg[iarg],"mol") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); - int imol = atom->find_molecule(arg[iarg+1]); - if (imol == -1) - error->all(FLERR,"Molecule template ID for fix pour does not exist"); + } else if (strcmp(arg[iarg], "mol") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command"); + int imol = atom->find_molecule(arg[iarg + 1]); + if (imol == -1) error->all(FLERR, "Molecule template ID for fix pour does not exist"); mode = MOLECULE; onemols = &atom->molecules[imol]; nmol = onemols[0]->nset; - delete [] molfrac; + delete[] molfrac; molfrac = new double[nmol]; - molfrac[0] = 1.0/nmol; - for (int i = 1; i < nmol-1; i++) molfrac[i] = molfrac[i-1] + 1.0/nmol; - molfrac[nmol-1] = 1.0; + molfrac[0] = 1.0 / nmol; + for (int i = 1; i < nmol - 1; i++) molfrac[i] = molfrac[i - 1] + 1.0 / nmol; + molfrac[nmol - 1] = 1.0; iarg += 2; - } else if (strcmp(arg[iarg],"molfrac") == 0) { - if (mode != MOLECULE) error->all(FLERR,"Illegal fix pour command"); - if (iarg+nmol+1 > narg) error->all(FLERR,"Illegal fix pour command"); - molfrac[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "molfrac") == 0) { + if (mode != MOLECULE) error->all(FLERR, "Illegal fix pour command"); + if (iarg + nmol + 1 > narg) error->all(FLERR, "Illegal fix pour command"); + molfrac[0] = utils::numeric(FLERR, arg[iarg + 1], false, lmp); for (int i = 1; i < nmol; i++) - molfrac[i] = molfrac[i-1] + utils::numeric(FLERR,arg[iarg+i+1],false,lmp); - if (molfrac[nmol-1] < 1.0-EPSILON || molfrac[nmol-1] > 1.0+EPSILON) - error->all(FLERR,"Illegal fix pour command"); - molfrac[nmol-1] = 1.0; - iarg += nmol+1; + molfrac[i] = molfrac[i - 1] + utils::numeric(FLERR, arg[iarg + i + 1], false, lmp); + if (molfrac[nmol - 1] < 1.0 - EPSILON || molfrac[nmol - 1] > 1.0 + EPSILON) + error->all(FLERR, "Illegal fix pour command"); + molfrac[nmol - 1] = 1.0; + iarg += nmol + 1; - } else if (strcmp(arg[iarg],"rigid") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); - delete [] idrigid; - idrigid = utils::strdup(arg[iarg+1]); + } else if (strcmp(arg[iarg], "rigid") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command"); + delete[] idrigid; + idrigid = utils::strdup(arg[iarg + 1]); rigidflag = 1; iarg += 2; - } else if (strcmp(arg[iarg],"shake") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); - delete [] idshake; - idshake = utils::strdup(arg[iarg+1]); + } else if (strcmp(arg[iarg], "shake") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command"); + delete[] idshake; + idshake = utils::strdup(arg[iarg + 1]); shakeflag = 1; iarg += 2; - } else if (strcmp(arg[iarg],"id") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); - if (strcmp(arg[iarg+1],"max") == 0) idnext = 0; - else if (strcmp(arg[iarg+1],"next") == 0) idnext = 1; - else error->all(FLERR,"Illegal fix pour command"); + } else if (strcmp(arg[iarg], "id") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command"); + if (strcmp(arg[iarg + 1], "max") == 0) + idnext = 0; + else if (strcmp(arg[iarg + 1], "next") == 0) + idnext = 1; + else + error->all(FLERR, "Illegal fix pour command"); iarg += 2; - } else if (strcmp(arg[iarg],"ignore") == 0) { + } else if (strcmp(arg[iarg], "ignore") == 0) { if (atom->line_flag) ignoreline = 1; if (atom->tri_flag) ignoretri = 1; if (ignoreline || ignoretri) ignoreflag = 1; iarg += 1; - } else if (strcmp(arg[iarg],"diam") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); - if (strcmp(arg[iarg+1],"one") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command"); + } else if (strcmp(arg[iarg], "diam") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command"); + if (strcmp(arg[iarg + 1], "one") == 0) { + if (iarg + 3 > narg) error->all(FLERR, "Illegal fix pour command"); dstyle = ONE; - radius_one = 0.5 * utils::numeric(FLERR,arg[iarg+2],false,lmp); + radius_one = 0.5 * utils::numeric(FLERR, arg[iarg + 2], false, lmp); radius_max = radius_one; iarg += 3; - } else if (strcmp(arg[iarg+1],"range") == 0) { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix pour command"); + } else if (strcmp(arg[iarg + 1], "range") == 0) { + if (iarg + 4 > narg) error->all(FLERR, "Illegal fix pour command"); dstyle = RANGE; - radius_lo = 0.5 * utils::numeric(FLERR,arg[iarg+2],false,lmp); - radius_hi = 0.5 * utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (radius_lo > radius_hi) error->all(FLERR,"Illegal fix pour command"); + radius_lo = 0.5 * utils::numeric(FLERR, arg[iarg + 2], false, lmp); + radius_hi = 0.5 * utils::numeric(FLERR, arg[iarg + 3], false, lmp); + if (radius_lo > radius_hi) error->all(FLERR, "Illegal fix pour command"); radius_max = radius_hi; iarg += 4; - } else if (strcmp(arg[iarg+1],"poly") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command"); + } else if (strcmp(arg[iarg + 1], "poly") == 0) { + if (iarg + 3 > narg) error->all(FLERR, "Illegal fix pour command"); dstyle = POLY; - npoly = utils::inumeric(FLERR,arg[iarg+2],false,lmp); - if (npoly <= 0) error->all(FLERR,"Illegal fix pour command"); - if (iarg+3 + 2*npoly > narg) - error->all(FLERR,"Illegal fix pour command"); + npoly = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); + if (npoly <= 0) error->all(FLERR, "Illegal fix pour command"); + if (iarg + 3 + 2 * npoly > narg) error->all(FLERR, "Illegal fix pour command"); radius_poly = new double[npoly]; frac_poly = new double[npoly]; iarg += 3; radius_max = 0.0; for (int i = 0; i < npoly; i++) { - radius_poly[i] = 0.5 * utils::numeric(FLERR,arg[iarg++],false,lmp); - frac_poly[i] = utils::numeric(FLERR,arg[iarg++],false,lmp); + radius_poly[i] = 0.5 * utils::numeric(FLERR, arg[iarg++], false, lmp); + frac_poly[i] = utils::numeric(FLERR, arg[iarg++], false, lmp); if (radius_poly[i] <= 0.0 || frac_poly[i] < 0.0) - error->all(FLERR,"Illegal fix pour command"); - radius_max = MAX(radius_max,radius_poly[i]); + error->all(FLERR, "Illegal fix pour command"); + radius_max = MAX(radius_max, radius_poly[i]); } double sum = 0.0; for (int i = 0; i < npoly; i++) sum += frac_poly[i]; if (fabs(sum - 1.0) > SMALL) - error->all(FLERR,"Fix pour polydisperse fractions do not sum to 1.0"); - } else error->all(FLERR,"Illegal fix pour command"); + error->all(FLERR, "Fix pour polydisperse fractions do not sum to 1.0"); + } else + error->all(FLERR, "Illegal fix pour command"); - } else if (strcmp(arg[iarg],"dens") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command"); - density_lo = utils::numeric(FLERR,arg[iarg+1],false,lmp); - density_hi = utils::numeric(FLERR,arg[iarg+2],false,lmp); - if (density_lo > density_hi) error->all(FLERR,"Illegal fix pour command"); + } else if (strcmp(arg[iarg], "dens") == 0) { + if (iarg + 3 > narg) error->all(FLERR, "Illegal fix pour command"); + density_lo = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + density_hi = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + if (density_lo > density_hi) error->all(FLERR, "Illegal fix pour command"); iarg += 3; - } else if (strcmp(arg[iarg],"vol") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix pour command"); - volfrac = utils::numeric(FLERR,arg[iarg+1],false,lmp); - maxattempt = utils::inumeric(FLERR,arg[iarg+2],false,lmp); + } else if (strcmp(arg[iarg], "vol") == 0) { + if (iarg + 3 > narg) error->all(FLERR, "Illegal fix pour command"); + volfrac = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + maxattempt = utils::inumeric(FLERR, arg[iarg + 2], false, lmp); iarg += 3; - } else if (strcmp(arg[iarg],"rate") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command"); - rate = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "rate") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix pour command"); + rate = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"vel") == 0) { + } else if (strcmp(arg[iarg], "vel") == 0) { if (domain->dimension == 3) { - if (iarg+6 > narg) error->all(FLERR,"Illegal fix pour command"); - vxlo = utils::numeric(FLERR,arg[iarg+1],false,lmp); - vxhi = utils::numeric(FLERR,arg[iarg+2],false,lmp); - vylo = utils::numeric(FLERR,arg[iarg+3],false,lmp); - vyhi = utils::numeric(FLERR,arg[iarg+4],false,lmp); - if (vxlo > vxhi || vylo > vyhi) - error->all(FLERR,"Illegal fix pour command"); - vz = utils::numeric(FLERR,arg[iarg+5],false,lmp); + if (iarg + 6 > narg) error->all(FLERR, "Illegal fix pour command"); + vxlo = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + vxhi = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + vylo = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + vyhi = utils::numeric(FLERR, arg[iarg + 4], false, lmp); + if (vxlo > vxhi || vylo > vyhi) error->all(FLERR, "Illegal fix pour command"); + vz = utils::numeric(FLERR, arg[iarg + 5], false, lmp); iarg += 6; } else { - if (iarg+4 > narg) error->all(FLERR,"Illegal fix pour command"); - vxlo = utils::numeric(FLERR,arg[iarg+1],false,lmp); - vxhi = utils::numeric(FLERR,arg[iarg+2],false,lmp); - vy = utils::numeric(FLERR,arg[iarg+3],false,lmp); + if (iarg + 4 > narg) error->all(FLERR, "Illegal fix pour command"); + vxlo = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + vxhi = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + vy = utils::numeric(FLERR, arg[iarg + 3], false, lmp); vz = 0.0; - if (vxlo > vxhi) error->all(FLERR,"Illegal fix pour command"); + if (vxlo > vxhi) error->all(FLERR, "Illegal fix pour command"); iarg += 4; } - } else error->all(FLERR,"Illegal fix pour command"); + } else + error->all(FLERR, "Illegal fix pour command"); } } @@ -1044,7 +1033,7 @@ void FixPour::options(int narg, char **arg) void FixPour::reset_dt() { - error->all(FLERR,"Cannot change timestep with fix pour"); + error->all(FLERR, "Cannot change timestep with fix pour"); } /* ---------------------------------------------------------------------- @@ -1053,10 +1042,12 @@ void FixPour::reset_dt() void *FixPour::extract(const char *str, int &itype) { - if (strcmp(str,"radius") == 0) { + if (strcmp(str, "radius") == 0) { if (mode == ATOM) { - if (itype == ntype) oneradius = radius_max; - else oneradius = 0.0; + if (itype == ntype) + oneradius = radius_max; + else + oneradius = 0.0; } else { @@ -1065,7 +1056,7 @@ void *FixPour::extract(const char *str, int &itype) oneradius = 0.0; for (int i = 0; i < nmol; i++) { - if (itype > ntype+onemols[i]->ntypes) continue; + if (itype > ntype + onemols[i]->ntypes) continue; double *radius = onemols[i]->radius; int *type = onemols[i]->type; int natoms = onemols[i]->natoms; @@ -1075,9 +1066,11 @@ void *FixPour::extract(const char *str, int &itype) // same as atom->avec->create_atom(), invoked in pre_exchange() for (int i = 0; i < natoms; i++) - if (type[i]+ntype == itype) { - if (radius) oneradius = MAX(oneradius,radius[i]); - else oneradius = MAX(oneradius,0.5); + if (type[i] + ntype == itype) { + if (radius) + oneradius = MAX(oneradius, radius[i]); + else + oneradius = MAX(oneradius, 0.5); } } } diff --git a/src/GRANULAR/fix_pour.h b/src/GRANULAR/fix_pour.h index 16ef3b7f32..d22b6e4b1e 100644 --- a/src/GRANULAR/fix_pour.h +++ b/src/GRANULAR/fix_pour.h @@ -37,7 +37,7 @@ class FixPour : public Fix { private: int ninsert, ntype, seed; - int iregion, mode, idnext, dstyle, npoly, rigidflag, shakeflag; + int mode, idnext, dstyle, npoly, rigidflag, shakeflag; int ignoreflag, ignoreline, ignoretri; double radius_one, radius_max; double radius_lo, radius_hi; @@ -52,6 +52,8 @@ class FixPour : public Fix { double xc, yc, rc; double grav; char *idrigid, *idshake; + char *idregion; + class Region *region; class Molecule **onemols; int nmol, natom_max; diff --git a/src/GRANULAR/fix_wall_gran_region.cpp b/src/GRANULAR/fix_wall_gran_region.cpp index f04b03b340..dc2bfce142 100644 --- a/src/GRANULAR/fix_wall_gran_region.cpp +++ b/src/GRANULAR/fix_wall_gran_region.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -35,21 +34,16 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) : - FixWallGran(lmp, narg, arg), region(nullptr), region_style(nullptr), - ncontact(nullptr), - walls(nullptr), history_many(nullptr), c2r(nullptr) + FixWallGran(lmp, narg, arg), region(nullptr), ncontact(nullptr), walls(nullptr), + history_many(nullptr), c2r(nullptr) { restart_global = 1; motion_resetflag = 0; - int iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for fix wall/gran/region does not exist"); - region = domain->regions[iregion]; - region_style = utils::strdup(region->style); + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for fix wall/gran/region does not exist", idregion); nregion = region->nregion; - - tmax = domain->regions[iregion]->tmax; + tmax = region->tmax; c2r = new int[tmax]; // re-allocate atom-based arrays with nshear @@ -67,8 +61,7 @@ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) : if (use_history) { int nlocal = atom->nlocal; - for (int i = 0; i < nlocal; i++) - ncontact[i] = 0; + for (int i = 0; i < nlocal; i++) ncontact[i] = 0; } } @@ -76,8 +69,8 @@ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) : FixWallGranRegion::~FixWallGranRegion() { - delete [] c2r; - delete [] region_style; + delete[] c2r; + delete[] region_style; memory->destroy(ncontact); memory->destroy(walls); @@ -90,25 +83,32 @@ void FixWallGranRegion::init() { FixWallGran::init(); - int iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for fix wall/gran/region does not exist"); - region = domain->regions[iregion]; + auto newregion = domain->get_region_by_id(idregion); + if (!newregion) error->all(FLERR, "Region {} for fix wall/gran/region does not exist", idregion); // check if region properties changed between runs // reset if restart info was inconsistent - if ((strcmp(idregion,region->id) != 0) - || (strcmp(region_style,region->style) != 0) - || (nregion != region->nregion)) { - error->warning(FLERR,"Region properties for region {} changed between " - "runs, resetting its motion",idregion); + if (newregion != region) { + region = newregion; + if (comm->me == 0) + error->warning(FLERR, + "Region properties for region {} changed between runs, resetting its motion", + idregion); + nregion = region->nregion; + tmax = region->tmax; + delete[] c2r; + c2r = new int[tmax]; + region = newregion; region->reset_vel(); } if (motion_resetflag) { - error->warning(FLERR,"Region properties for region {} are inconsistent " - "with restart file, resetting its motion",idregion); + if (comm->me == 0) + error->warning(FLERR, + "Region properties for region {} are inconsistent with restart file, " + "resetting its motion", + idregion); region->reset_vel(); } } @@ -117,8 +117,8 @@ void FixWallGranRegion::init() void FixWallGranRegion::post_force(int /*vflag*/) { - int i,m,nc,iwall; - double dx,dy,dz,rsq,meff; + int i, m, nc, iwall; + double dx, dy, dz, rsq, meff; double vwall[3]; // do not update shear history during setup @@ -133,17 +133,19 @@ void FixWallGranRegion::post_force(int /*vflag*/) if (neighbor->ago == 0 && fix_rigid) { int tmp; - int *body = (int *) fix_rigid->extract("body",tmp); - auto mass_body = (double *) fix_rigid->extract("masstotal",tmp); + int *body = (int *) fix_rigid->extract("body", tmp); + auto mass_body = (double *) fix_rigid->extract("masstotal", tmp); if (atom->nmax > nmax) { memory->destroy(mass_rigid); nmax = atom->nmax; - memory->create(mass_rigid,nmax,"wall/gran:mass_rigid"); + memory->create(mass_rigid, nmax, "wall/gran:mass_rigid"); } int nlocal = atom->nlocal; for (i = 0; i < nlocal; i++) { - if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; - else mass_rigid[i] = 0.0; + if (body[i] >= 0) + mass_rigid[i] = mass_body[body[i]]; + else + mass_rigid[i] = 0.0; } } @@ -169,23 +171,18 @@ void FixWallGranRegion::post_force(int /*vflag*/) region->set_velocity(); } - if (peratom_flag) { - clear_stored_contacts(); - } + if (peratom_flag) { clear_stored_contacts(); } for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (!region->match(x[i][0],x[i][1],x[i][2])) continue; + if (!region->match(x[i][0], x[i][1], x[i][2])) continue; if (pairstyle == FixWallGran::GRANULAR && normal_model == FixWallGran::JKR) { - nc = region->surface(x[i][0],x[i][1],x[i][2], - radius[i]+pulloff_distance(radius[i])); + nc = region->surface(x[i][0], x[i][1], x[i][2], radius[i] + pulloff_distance(radius[i])); + } else { + nc = region->surface(x[i][0], x[i][1], x[i][2], radius[i]); } - else{ - nc = region->surface(x[i][0],x[i][1],x[i][2],radius[i]); - } - if (nc > tmax) - error->one(FLERR,"Too many wall/gran/region contacts for one particle"); + if (nc > tmax) error->one(FLERR, "Too many wall/gran/region contacts for one particle"); // shear history maintenance // update ncontact,walls,shear2many for particle I @@ -204,11 +201,11 @@ void FixWallGranRegion::post_force(int /*vflag*/) if (ncontact[i] == 0) { ncontact[i] = 1; walls[i][0] = iwall; - for (m = 0; m < size_history; m++) - history_many[i][0][m] = 0.0; + for (m = 0; m < size_history; m++) history_many[i][0][m] = 0.0; } else if (ncontact[i] > 1 || iwall != walls[i][0]) - update_contacts(i,nc); - } else update_contacts(i,nc); + update_contacts(i, nc); + } else + update_contacts(i, nc); } // process current contacts @@ -217,12 +214,11 @@ void FixWallGranRegion::post_force(int /*vflag*/) // rsq = squared contact distance // xc = contact point - rsq = region->contact[ic].r*region->contact[ic].r; + rsq = region->contact[ic].r * region->contact[ic].r; if (pairstyle == FixWallGran::GRANULAR && normal_model == FixWallGran::JKR) { - if (history_many[i][c2r[ic]][0] == 0.0 && rsq > radius[i]*radius[i]) { - for (m = 0; m < size_history; m++) - history_many[i][0][m] = 0.0; + if (history_many[i][c2r[ic]][0] == 0.0 && rsq > radius[i] * radius[i]) { + for (m = 0; m < size_history; m++) history_many[i][0][m] = 0.0; continue; } } @@ -256,20 +252,16 @@ void FixWallGranRegion::post_force(int /*vflag*/) contact = nullptr; if (pairstyle == FixWallGran::HOOKE) - hooke(rsq,dx,dy,dz,vwall,v[i],f[i], - omega[i],torque[i],radius[i],meff, contact); + hooke(rsq, dx, dy, dz, vwall, v[i], f[i], omega[i], torque[i], radius[i], meff, contact); else if (pairstyle == FixWallGran::HOOKE_HISTORY) - hooke_history(rsq,dx,dy,dz,vwall,v[i],f[i], - omega[i],torque[i],radius[i],meff, - history_many[i][c2r[ic]], contact); + hooke_history(rsq, dx, dy, dz, vwall, v[i], f[i], omega[i], torque[i], radius[i], meff, + history_many[i][c2r[ic]], contact); else if (pairstyle == FixWallGran::HERTZ_HISTORY) - hertz_history(rsq,dx,dy,dz,vwall,region->contact[ic].radius, - v[i],f[i],omega[i],torque[i], - radius[i],meff,history_many[i][c2r[ic]], contact); + hertz_history(rsq, dx, dy, dz, vwall, region->contact[ic].radius, v[i], f[i], omega[i], + torque[i], radius[i], meff, history_many[i][c2r[ic]], contact); else if (pairstyle == FixWallGran::GRANULAR) - granular(rsq,dx,dy,dz,vwall,region->contact[ic].radius, - v[i],f[i],omega[i],torque[i], - radius[i],meff,history_many[i][c2r[ic]],contact); + granular(rsq, dx, dy, dz, vwall, region->contact[ic].radius, v[i], f[i], omega[i], + torque[i], radius[i], meff, history_many[i][c2r[ic]], contact); } } } @@ -285,7 +277,7 @@ void FixWallGranRegion::post_force(int /*vflag*/) void FixWallGranRegion::update_contacts(int i, int nc) { - int j,m,iold,nold,ilast,inew,iadd,iwall; + int j, m, iold, nold, ilast, inew, iadd, iwall; // loop over old contacts // if not in new contact list: @@ -296,12 +288,12 @@ void FixWallGranRegion::update_contacts(int i, int nc) for (m = 0; m < nc; m++) if (region->contact[m].iwall == walls[i][iold]) break; if (m >= nc) { - ilast = ncontact[i]-1; - for (j = 0; j < size_history; j++) - history_many[i][iold][j] = history_many[i][ilast][j]; + ilast = ncontact[i] - 1; + for (j = 0; j < size_history; j++) history_many[i][iold][j] = history_many[i][ilast][j]; walls[i][iold] = walls[i][ilast]; ncontact[i]--; - } else iold++; + } else + iold++; } // loop over new contacts @@ -315,13 +307,13 @@ void FixWallGranRegion::update_contacts(int i, int nc) iwall = region->contact[inew].iwall; for (m = 0; m < nold; m++) if (walls[i][m] == iwall) break; - if (m < nold) c2r[m] = inew; + if (m < nold) + c2r[m] = inew; else { iadd = ncontact[i]; c2r[iadd] = inew; - for (j = 0; j < size_history; j++) - history_many[i][iadd][j] = 0.0; + for (j = 0; j < size_history; j++) history_many[i][iadd][j] = 0.0; walls[i][iadd] = iwall; ncontact[i]++; } @@ -336,12 +328,12 @@ double FixWallGranRegion::memory_usage() { int nmax = atom->nmax; double bytes = 0.0; - if (use_history) { // shear history - bytes += (double)nmax * sizeof(int); // ncontact - bytes += (double)nmax*tmax * sizeof(int); // walls - bytes += (double)nmax*tmax*size_history * sizeof(double); // history_many + if (use_history) { // shear history + bytes += (double) nmax * sizeof(int); // ncontact + bytes += (double) nmax * tmax * sizeof(int); // walls + bytes += (double) nmax * tmax * size_history * sizeof(double); // history_many } - if (fix_rigid) bytes += (double)nmax * sizeof(int); // mass_rigid + if (fix_rigid) bytes += (double) nmax * sizeof(int); // mass_rigid return bytes; } @@ -352,12 +344,11 @@ double FixWallGranRegion::memory_usage() void FixWallGranRegion::grow_arrays(int nmax) { if (use_history) { - memory->grow(ncontact,nmax,"fix_wall_gran:ncontact"); - memory->grow(walls,nmax,tmax,"fix_wall_gran:walls"); - memory->grow(history_many,nmax,tmax,size_history,"fix_wall_gran:history_many"); + memory->grow(ncontact, nmax, "fix_wall_gran:ncontact"); + memory->grow(walls, nmax, tmax, "fix_wall_gran:walls"); + memory->grow(history_many, nmax, tmax, size_history, "fix_wall_gran:history_many"); } - if (peratom_flag) - memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom"); + if (peratom_flag) memory->grow(array_atom, nmax, size_peratom_cols, "fix_wall_gran:array_atom"); } /* ---------------------------------------------------------------------- @@ -366,21 +357,19 @@ void FixWallGranRegion::grow_arrays(int nmax) void FixWallGranRegion::copy_arrays(int i, int j, int /*delflag*/) { - int m,n,iwall; + int m, n, iwall; if (use_history) { n = ncontact[i]; for (iwall = 0; iwall < n; iwall++) { walls[j][iwall] = walls[i][iwall]; - for (m = 0; m < size_history; m++) - history_many[j][iwall][m] = history_many[i][iwall][m]; + for (m = 0; m < size_history; m++) history_many[j][iwall][m] = history_many[i][iwall][m]; } ncontact[j] = ncontact[i]; } if (peratom_flag) { - for (int m = 0; m < size_peratom_cols; m++) - array_atom[j][m] = array_atom[i][m]; + for (int m = 0; m < size_peratom_cols; m++) array_atom[j][m] = array_atom[i][m]; } } @@ -390,11 +379,9 @@ void FixWallGranRegion::copy_arrays(int i, int j, int /*delflag*/) void FixWallGranRegion::set_arrays(int i) { - if (use_history) - ncontact[i] = 0; + if (use_history) ncontact[i] = 0; if (peratom_flag) { - for (int m = 0; m < size_peratom_cols; m++) - array_atom[i][m] = 0; + for (int m = 0; m < size_peratom_cols; m++) array_atom[i][m] = 0; } } @@ -412,13 +399,11 @@ int FixWallGranRegion::pack_exchange(int i, double *buf) buf[n++] = ubuf(count).d; for (int iwall = 0; iwall < count; iwall++) { buf[n++] = ubuf(walls[i][iwall]).d; - for (m = 0; m < size_history; m++) - buf[n++] = history_many[i][iwall][m]; + for (m = 0; m < size_history; m++) buf[n++] = history_many[i][iwall][m]; } } if (peratom_flag) { - for (int m = 0; m < size_peratom_cols; m++) - buf[n++] = array_atom[i][m]; + for (int m = 0; m < size_peratom_cols; m++) buf[n++] = array_atom[i][m]; } return n; @@ -432,19 +417,16 @@ int FixWallGranRegion::unpack_exchange(int nlocal, double *buf) { int m; - int n = 0; if (use_history) { int count = ncontact[nlocal] = (int) ubuf(buf[n++]).i; for (int iwall = 0; iwall < count; iwall++) { walls[nlocal][iwall] = (int) ubuf(buf[n++]).i; - for (m = 0; m < size_history; m++) - history_many[nlocal][iwall][m] = buf[n++]; + for (m = 0; m < size_history; m++) history_many[nlocal][iwall][m] = buf[n++]; } } if (peratom_flag) { - for (int m = 0; m < size_peratom_cols; m++) - array_atom[nlocal][m] = buf[n++]; + for (int m = 0; m < size_peratom_cols; m++) array_atom[nlocal][m] = buf[n++]; } return n; @@ -466,8 +448,7 @@ int FixWallGranRegion::pack_restart(int i, double *buf) buf[n++] = ubuf(count).d; for (int iwall = 0; iwall < count; iwall++) { buf[n++] = ubuf(walls[i][iwall]).d; - for (m = 0; m < size_history; m++) - buf[n++] = history_many[i][iwall][m]; + for (m = 0; m < size_history; m++) buf[n++] = history_many[i][iwall][m]; } // pack buf[0] this way because other fixes unpack it buf[0] = n; @@ -490,14 +471,13 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth) // unpack the Nth first values this way because other fixes pack them int m = 0; - for (int i = 0; i < nth; i++) m += static_cast (extra[nlocal][m]); + for (int i = 0; i < nth; i++) m += static_cast(extra[nlocal][m]); m++; int count = ncontact[nlocal] = (int) ubuf(extra[nlocal][m++]).i; for (int iwall = 0; iwall < count; iwall++) { walls[nlocal][iwall] = (int) ubuf(extra[nlocal][m++]).i; - for (k = 0; k < size_history; k++) - history_many[nlocal][iwall][k] = extra[nlocal][m++]; + for (k = 0; k < size_history; k++) history_many[nlocal][iwall][k] = extra[nlocal][m++]; } } @@ -508,7 +488,7 @@ void FixWallGranRegion::unpack_restart(int nlocal, int nth) int FixWallGranRegion::maxsize_restart() { if (!use_history) return 0; - return 2 + tmax*(size_history+1); + return 2 + tmax * (size_history + 1); } /* ---------------------------------------------------------------------- @@ -518,7 +498,7 @@ int FixWallGranRegion::maxsize_restart() int FixWallGranRegion::size_restart(int nlocal) { if (!use_history) return 0; - return 2 + ncontact[nlocal]*(size_history+1); + return 2 + ncontact[nlocal] * (size_history + 1); } /* ---------------------------------------------------------------------- @@ -530,7 +510,7 @@ void FixWallGranRegion::write_restart(FILE *fp) if (comm->me) return; int len = 0; region->length_restart_string(len); - fwrite(&len, sizeof(int),1,fp); + fwrite(&len, sizeof(int), 1, fp); region->write_restart(fp); } @@ -541,5 +521,5 @@ void FixWallGranRegion::write_restart(FILE *fp) void FixWallGranRegion::restart(char *buf) { int n = 0; - if (!region->restart(buf,n)) motion_resetflag = 1; + if (!region->restart(buf, n)) motion_resetflag = 1; } diff --git a/src/MACHDYN/fix_smd_setvel.cpp b/src/MACHDYN/fix_smd_setvel.cpp index 7b76e08ec6..1b5578c7e0 100644 --- a/src/MACHDYN/fix_smd_setvel.cpp +++ b/src/MACHDYN/fix_smd_setvel.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- * * *** Smooth Mach Dynamics *** @@ -41,320 +40,297 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum { - NONE, CONSTANT, EQUAL, ATOM -}; +enum { NONE, CONSTANT, EQUAL, ATOM }; /* ---------------------------------------------------------------------- */ FixSMDSetVel::FixSMDSetVel(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) { - if (narg < 6) - error->all(FLERR, "Illegal fix setvelocity command"); + Fix(lmp, narg, arg), xstr(nullptr), ystr(nullptr), zstr(nullptr), idregion(nullptr), + region(nullptr), sforce(nullptr) +{ + if (narg < 6) error->all(FLERR, "Illegal fix setvelocity command"); - dynamic_group_allow = 1; - vector_flag = 1; - size_vector = 3; - global_freq = 1; - extvector = 1; + dynamic_group_allow = 1; + vector_flag = 1; + size_vector = 3; + global_freq = 1; + extvector = 1; - xstr = ystr = zstr = nullptr; + if (strstr(arg[3], "v_") == arg[3]) { + xstr = utils::strdup(&arg[3][2]); + } else if (strcmp(arg[3], "NULL") == 0) { + xstyle = NONE; + } else { + xvalue = utils::numeric(FLERR, arg[3], false, lmp); + xstyle = CONSTANT; + } + if (strstr(arg[4], "v_") == arg[4]) { + ystr = utils::strdup(&arg[4][2]); + } else if (strcmp(arg[4], "NULL") == 0) { + ystyle = NONE; + } else { + yvalue = utils::numeric(FLERR, arg[4], false, lmp); + ystyle = CONSTANT; + } + if (strstr(arg[5], "v_") == arg[5]) { + zstr = utils::strdup(&arg[5][2]); + } else if (strcmp(arg[5], "NULL") == 0) { + zstyle = NONE; + } else { + zvalue = utils::numeric(FLERR, arg[5], false, lmp); + zstyle = CONSTANT; + } - if (strstr(arg[3], "v_") == arg[3]) { - xstr = utils::strdup( &arg[3][2]); - } else if (strcmp(arg[3], "NULL") == 0) { - xstyle = NONE; - } else { - xvalue = utils::numeric(FLERR, arg[3],false,lmp); - xstyle = CONSTANT; - } - if (strstr(arg[4], "v_") == arg[4]) { - ystr = utils::strdup( &arg[4][2]); - } else if (strcmp(arg[4], "NULL") == 0) { - ystyle = NONE; - } else { - yvalue = utils::numeric(FLERR, arg[4],false,lmp); - ystyle = CONSTANT; - } - if (strstr(arg[5], "v_") == arg[5]) { - zstr = utils::strdup( &arg[5][2]); - } else if (strcmp(arg[5], "NULL") == 0) { - zstyle = NONE; - } else { - zvalue = utils::numeric(FLERR, arg[5],false,lmp); - zstyle = CONSTANT; - } + // optional args - // optional args + int iarg = 6; + while (iarg < narg) { + if (strcmp(arg[iarg], "region") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix setvelocity command"); + region = domain->get_region_by_id(arg[iarg + 1]); + if (!region) error->all(FLERR, "Region {} for fix setvelocity does not exist", arg[iarg + 1]); + idregion = utils::strdup(arg[iarg + 1]); + iarg += 2; + } else + error->all(FLERR, "Illegal fix setvelocity command"); + } - iregion = -1; - idregion = nullptr; + force_flag = 0; + foriginal[0] = foriginal[1] = foriginal[2] = 0.0; - int iarg = 6; - while (iarg < narg) { - if (strcmp(arg[iarg], "region") == 0) { - if (iarg + 2 > narg) - error->all(FLERR, "Illegal fix setvelocity command"); - iregion = domain->find_region(arg[iarg + 1]); - if (iregion == -1) - error->all(FLERR, "Region ID for fix setvelocity does not exist"); - idregion = utils::strdup( arg[iarg + 1]); - iarg += 2; - } else - error->all(FLERR, "Illegal fix setvelocity command"); - } - - force_flag = 0; - foriginal[0] = foriginal[1] = foriginal[2] = 0.0; - - maxatom = atom->nmax; - memory->create(sforce, maxatom, 3, "setvelocity:sforce"); + maxatom = atom->nmax; + memory->create(sforce, maxatom, 3, "setvelocity:sforce"); } /* ---------------------------------------------------------------------- */ -FixSMDSetVel::~FixSMDSetVel() { - delete[] xstr; - delete[] ystr; - delete[] zstr; - delete[] idregion; - memory->destroy(sforce); +FixSMDSetVel::~FixSMDSetVel() +{ + delete[] xstr; + delete[] ystr; + delete[] zstr; + delete[] idregion; + memory->destroy(sforce); } /* ---------------------------------------------------------------------- */ -int FixSMDSetVel::setmask() { - int mask = 0; - //mask |= INITIAL_INTEGRATE; - mask |= POST_FORCE; - return mask; +int FixSMDSetVel::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + return mask; } /* ---------------------------------------------------------------------- */ -void FixSMDSetVel::init() { - // check variables +void FixSMDSetVel::init() +{ + // check variables - if (xstr) { - xvar = input->variable->find(xstr); - if (xvar < 0) - error->all(FLERR, "Variable name for fix setvelocity does not exist"); - if (input->variable->equalstyle(xvar)) - xstyle = EQUAL; - else if (input->variable->atomstyle(xvar)) - xstyle = ATOM; - else - error->all(FLERR, "Variable for fix setvelocity is invalid style"); - } - if (ystr) { - yvar = input->variable->find(ystr); - if (yvar < 0) - error->all(FLERR, "Variable name for fix setvelocity does not exist"); - if (input->variable->equalstyle(yvar)) - ystyle = EQUAL; - else if (input->variable->atomstyle(yvar)) - ystyle = ATOM; - else - error->all(FLERR, "Variable for fix setvelocity is invalid style"); - } - if (zstr) { - zvar = input->variable->find(zstr); - if (zvar < 0) - error->all(FLERR, "Variable name for fix setvelocity does not exist"); - if (input->variable->equalstyle(zvar)) - zstyle = EQUAL; - else if (input->variable->atomstyle(zvar)) - zstyle = ATOM; - else - error->all(FLERR, "Variable for fix setvelocity is invalid style"); - } + if (xstr) { + xvar = input->variable->find(xstr); + if (xvar < 0) error->all(FLERR, "Variable name for fix setvelocity does not exist"); + if (input->variable->equalstyle(xvar)) + xstyle = EQUAL; + else if (input->variable->atomstyle(xvar)) + xstyle = ATOM; + else + error->all(FLERR, "Variable for fix setvelocity is invalid style"); + } + if (ystr) { + yvar = input->variable->find(ystr); + if (yvar < 0) error->all(FLERR, "Variable name for fix setvelocity does not exist"); + if (input->variable->equalstyle(yvar)) + ystyle = EQUAL; + else if (input->variable->atomstyle(yvar)) + ystyle = ATOM; + else + error->all(FLERR, "Variable for fix setvelocity is invalid style"); + } + if (zstr) { + zvar = input->variable->find(zstr); + if (zvar < 0) error->all(FLERR, "Variable name for fix setvelocity does not exist"); + if (input->variable->equalstyle(zvar)) + zstyle = EQUAL; + else if (input->variable->atomstyle(zvar)) + zstyle = ATOM; + else + error->all(FLERR, "Variable for fix setvelocity is invalid style"); + } - // set index and check validity of region + // set index and check validity of region - if (iregion >= 0) { - iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR, "Region ID for fix setvelocity does not exist"); - } + if (idregion) { + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for fix setvelocity does not exist", idregion); + } - if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM) - varflag = ATOM; - else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL) - varflag = EQUAL; - else - varflag = CONSTANT; + if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM) + varflag = ATOM; + else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL) + varflag = EQUAL; + else + varflag = CONSTANT; - // cannot use non-zero forces for a minimization since no energy is integrated - // use fix addforce instead + // cannot use non-zero forces for a minimization since no energy is integrated + // use fix addforce instead - int flag = 0; - if (update->whichflag == 2) { - if (xstyle == EQUAL || xstyle == ATOM) - flag = 1; - if (ystyle == EQUAL || ystyle == ATOM) - flag = 1; - if (zstyle == EQUAL || zstyle == ATOM) - flag = 1; - if (xstyle == CONSTANT && xvalue != 0.0) - flag = 1; - if (ystyle == CONSTANT && yvalue != 0.0) - flag = 1; - if (zstyle == CONSTANT && zvalue != 0.0) - flag = 1; - } - if (flag) - error->all(FLERR, "Cannot use non-zero forces in an energy minimization"); + int flag = 0; + if (update->whichflag == 2) { + if (xstyle == EQUAL || xstyle == ATOM) flag = 1; + if (ystyle == EQUAL || ystyle == ATOM) flag = 1; + if (zstyle == EQUAL || zstyle == ATOM) flag = 1; + if (xstyle == CONSTANT && xvalue != 0.0) flag = 1; + if (ystyle == CONSTANT && yvalue != 0.0) flag = 1; + if (zstyle == CONSTANT && zvalue != 0.0) flag = 1; + } + if (flag) error->all(FLERR, "Cannot use non-zero forces in an energy minimization"); } /* ---------------------------------------------------------------------- */ -void FixSMDSetVel::setup(int vflag) { - if (utils::strmatch(update->integrate_style,"^verlet")) - post_force(vflag); - else - error->all(FLERR,"Fix smd/setvel does not support RESPA"); +void FixSMDSetVel::setup(int vflag) +{ + if (utils::strmatch(update->integrate_style, "^verlet")) + post_force(vflag); + else + error->all(FLERR, "Fix smd/setvel does not support RESPA"); } /* ---------------------------------------------------------------------- */ -void FixSMDSetVel::min_setup(int vflag) { - post_force(vflag); +void FixSMDSetVel::min_setup(int vflag) +{ + post_force(vflag); } /* ---------------------------------------------------------------------- */ -//void FixSMDSetVel::initial_integrate(int vflag) { -void FixSMDSetVel::post_force(int /*vflag*/) { - double **x = atom->x; - double **f = atom->f; - double **v = atom->v; - double **vest = atom->vest; - int *mask = atom->mask; - int nlocal = atom->nlocal; +void FixSMDSetVel::post_force(int /*vflag*/) +{ + double **x = atom->x; + double **f = atom->f; + double **v = atom->v; + double **vest = atom->vest; + int *mask = atom->mask; + int nlocal = atom->nlocal; - // update region if necessary + // update region if necessary - Region *region = nullptr; - if (iregion >= 0) { - region = domain->regions[iregion]; - region->prematch(); + if (region) region->prematch(); + + // reallocate sforce array if necessary + + if (varflag == ATOM && atom->nmax > maxatom) { + maxatom = atom->nmax; + memory->destroy(sforce); + memory->create(sforce, maxatom, 3, "setvelocity:sforce"); + } + + foriginal[0] = foriginal[1] = foriginal[2] = 0.0; + force_flag = 0; + + if (varflag == CONSTANT) { + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue; + foriginal[0] += f[i][0]; + foriginal[1] += f[i][1]; + foriginal[2] += f[i][2]; + if (xstyle) { + v[i][0] = xvalue; + vest[i][0] = xvalue; + f[i][0] = 0.0; + } + if (ystyle) { + v[i][1] = yvalue; + vest[i][1] = yvalue; + f[i][1] = 0.0; + } + if (zstyle) { + v[i][2] = zvalue; + vest[i][2] = zvalue; + f[i][2] = 0.0; + } + } + + // variable force, wrap with clear/add + + } else { + + modify->clearstep_compute(); + + if (xstyle == EQUAL) + xvalue = input->variable->compute_equal(xvar); + else if (xstyle == ATOM) + input->variable->compute_atom(xvar, igroup, &sforce[0][0], 3, 0); + if (ystyle == EQUAL) + yvalue = input->variable->compute_equal(yvar); + else if (ystyle == ATOM) + input->variable->compute_atom(yvar, igroup, &sforce[0][1], 3, 0); + if (zstyle == EQUAL) + zvalue = input->variable->compute_equal(zvar); + else if (zstyle == ATOM) + input->variable->compute_atom(zvar, igroup, &sforce[0][2], 3, 0); + + modify->addstep_compute(update->ntimestep + 1); + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + if (region && !region->match(x[i][0], x[i][1], x[i][2])) continue; + foriginal[0] += f[i][0]; + foriginal[1] += f[i][1]; + foriginal[2] += f[i][2]; + if (xstyle == ATOM) { + vest[i][0] = v[i][0] = sforce[i][0]; + f[i][0] = 0.0; + } else if (xstyle) { + vest[i][0] = v[i][0] = xvalue; + f[i][0] = 0.0; } - // reallocate sforce array if necessary - - if (varflag == ATOM && atom->nmax > maxatom) { - maxatom = atom->nmax; - memory->destroy(sforce); - memory->create(sforce, maxatom, 3, "setvelocity:sforce"); + if (ystyle == ATOM) { + vest[i][1] = v[i][1] = sforce[i][1]; + f[i][1] = 0.0; + } else if (ystyle) { + vest[i][1] = v[i][1] = yvalue; + f[i][1] = 0.0; } - foriginal[0] = foriginal[1] = foriginal[2] = 0.0; - force_flag = 0; - - if (varflag == CONSTANT) { - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (region && !region->match(x[i][0], x[i][1], x[i][2])) - continue; - foriginal[0] += f[i][0]; - foriginal[1] += f[i][1]; - foriginal[2] += f[i][2]; - if (xstyle) { - v[i][0] = xvalue; - vest[i][0] = xvalue; - f[i][0] = 0.0; - } - if (ystyle) { - v[i][1] = yvalue; - vest[i][1] = yvalue; - f[i][1] = 0.0; - } - if (zstyle) { - v[i][2] = zvalue; - vest[i][2] = zvalue; - f[i][2] = 0.0; - } - } - - // variable force, wrap with clear/add - - } else { - - modify->clearstep_compute(); - - if (xstyle == EQUAL) - xvalue = input->variable->compute_equal(xvar); - else if (xstyle == ATOM) - input->variable->compute_atom(xvar, igroup, &sforce[0][0], 3, 0); - if (ystyle == EQUAL) - yvalue = input->variable->compute_equal(yvar); - else if (ystyle == ATOM) - input->variable->compute_atom(yvar, igroup, &sforce[0][1], 3, 0); - if (zstyle == EQUAL) - zvalue = input->variable->compute_equal(zvar); - else if (zstyle == ATOM) - input->variable->compute_atom(zvar, igroup, &sforce[0][2], 3, 0); - - modify->addstep_compute(update->ntimestep + 1); - - //printf("setting velocity at timestep %d\n", update->ntimestep); - - for (int i = 0; i < nlocal; i++) - if (mask[i] & groupbit) { - if (region && !region->match(x[i][0], x[i][1], x[i][2])) - continue; - foriginal[0] += f[i][0]; - foriginal[1] += f[i][1]; - foriginal[2] += f[i][2]; - if (xstyle == ATOM) { - vest[i][0] = v[i][0] = sforce[i][0]; - f[i][0] = 0.0; - } else if (xstyle) { - vest[i][0] = v[i][0] = xvalue; - f[i][0] = 0.0; - } - - if (ystyle == ATOM) { - vest[i][1] = v[i][1] = sforce[i][1]; - f[i][1] = 0.0; - } else if (ystyle) { - vest[i][1] = v[i][1] = yvalue; - f[i][1] = 0.0; - } - - if (zstyle == ATOM) { - vest[i][2] = v[i][2] = sforce[i][2]; - f[i][2] = 0.0; - } else if (zstyle) { - vest[i][2] = v[i][2] = zvalue; - f[i][2] = 0.0; - } - - } + if (zstyle == ATOM) { + vest[i][2] = v[i][2] = sforce[i][2]; + f[i][2] = 0.0; + } else if (zstyle) { + vest[i][2] = v[i][2] = zvalue; + f[i][2] = 0.0; } + } + } } /* ---------------------------------------------------------------------- return components of total force on fix group before force was changed ------------------------------------------------------------------------- */ -double FixSMDSetVel::compute_vector(int n) { -// only sum across procs one time +double FixSMDSetVel::compute_vector(int n) +{ + // only sum across procs one time - if (force_flag == 0) { - MPI_Allreduce(foriginal, foriginal_all, 3, MPI_DOUBLE, MPI_SUM, world); - force_flag = 1; - } - return foriginal_all[n]; + if (force_flag == 0) { + MPI_Allreduce(foriginal, foriginal_all, 3, MPI_DOUBLE, MPI_SUM, world); + force_flag = 1; + } + return foriginal_all[n]; } /* ---------------------------------------------------------------------- memory usage of local atom-based array ------------------------------------------------------------------------- */ -double FixSMDSetVel::memory_usage() { - double bytes = 0.0; - if (varflag == ATOM) - bytes = atom->nmax * 3 * sizeof(double); - return bytes; +double FixSMDSetVel::memory_usage() +{ + double bytes = 0.0; + if (varflag == ATOM) bytes = atom->nmax * 3 * sizeof(double); + return bytes; } diff --git a/src/MACHDYN/fix_smd_setvel.h b/src/MACHDYN/fix_smd_setvel.h index 30071d0875..5fd0d2637c 100644 --- a/src/MACHDYN/fix_smd_setvel.h +++ b/src/MACHDYN/fix_smd_setvel.h @@ -50,9 +50,10 @@ class FixSMDSetVel : public Fix { private: double xvalue, yvalue, zvalue; - int varflag, iregion; + int varflag; char *xstr, *ystr, *zstr; char *idregion; + class Region *region; int xvar, yvar, zvar, xstyle, ystyle, zstyle; double foriginal[3], foriginal_all[3]; int force_flag; diff --git a/src/MC/fix_atom_swap.cpp b/src/MC/fix_atom_swap.cpp index 690e628fe3..2f8cd64169 100644 --- a/src/MC/fix_atom_swap.cpp +++ b/src/MC/fix_atom_swap.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -40,9 +39,9 @@ #include "region.h" #include "update.h" -#include #include #include +#include #include using namespace LAMMPS_NS; @@ -51,13 +50,12 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - idregion(nullptr), type_list(nullptr), mu(nullptr), qtype(nullptr), - sqrt_mass_ratio(nullptr), local_swap_iatom_list(nullptr), - local_swap_jatom_list(nullptr), local_swap_atom_list(nullptr), - random_equal(nullptr), random_unequal(nullptr), c_pe(nullptr) + Fix(lmp, narg, arg), region(nullptr), idregion(nullptr), type_list(nullptr), mu(nullptr), + qtype(nullptr), sqrt_mass_ratio(nullptr), local_swap_iatom_list(nullptr), + local_swap_jatom_list(nullptr), local_swap_atom_list(nullptr), random_equal(nullptr), + random_unequal(nullptr), c_pe(nullptr) { - if (narg < 10) error->all(FLERR,"Illegal fix atom/swap command"); + if (narg < 10) error->all(FLERR, "Illegal fix atom/swap command"); dynamic_group_allow = 1; @@ -70,33 +68,33 @@ FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) : // required args - nevery = utils::inumeric(FLERR,arg[3],false,lmp); - ncycles = utils::inumeric(FLERR,arg[4],false,lmp); - seed = utils::inumeric(FLERR,arg[5],false,lmp); - double temperature = utils::numeric(FLERR,arg[6],false,lmp); + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + ncycles = utils::inumeric(FLERR, arg[4], false, lmp); + seed = utils::inumeric(FLERR, arg[5], false, lmp); + double temperature = utils::numeric(FLERR, arg[6], false, lmp); - if (nevery <= 0) error->all(FLERR,"Illegal fix atom/swap command"); - if (ncycles < 0) error->all(FLERR,"Illegal fix atom/swap command"); - if (seed <= 0) error->all(FLERR,"Illegal fix atom/swap command"); - if (temperature <= 0.0) error->all(FLERR,"Illegal fix atom/swap command"); + if (nevery <= 0) error->all(FLERR, "Illegal fix atom/swap command"); + if (ncycles < 0) error->all(FLERR, "Illegal fix atom/swap command"); + if (seed <= 0) error->all(FLERR, "Illegal fix atom/swap command"); + if (temperature <= 0.0) error->all(FLERR, "Illegal fix atom/swap command"); - beta = 1.0/(force->boltz*temperature); + beta = 1.0 / (force->boltz * temperature); - memory->create(type_list,atom->ntypes,"atom/swap:type_list"); - memory->create(mu,atom->ntypes+1,"atom/swap:mu"); + memory->create(type_list, atom->ntypes, "atom/swap:type_list"); + memory->create(mu, atom->ntypes + 1, "atom/swap:mu"); for (int i = 1; i <= atom->ntypes; i++) mu[i] = 0.0; // read options from end of input line - options(narg-7,&arg[7]); + options(narg - 7, &arg[7]); // random number generator, same for all procs - random_equal = new RanPark(lmp,seed); + random_equal = new RanPark(lmp, seed); // random number generator, not the same for all procs - random_unequal = new RanPark(lmp,seed); + random_unequal = new RanPark(lmp, seed); // set up reneighboring @@ -115,9 +113,10 @@ FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) : // set comm size needed by this Fix - if (atom->q_flag) comm_forward = 2; - else comm_forward = 1; - + if (atom->q_flag) + comm_forward = 2; + else + comm_forward = 1; } /* ---------------------------------------------------------------------- */ @@ -130,7 +129,7 @@ FixAtomSwap::~FixAtomSwap() memory->destroy(sqrt_mass_ratio); memory->destroy(local_swap_iatom_list); memory->destroy(local_swap_jatom_list); - if (regionflag) delete [] idregion; + delete[] idregion; delete random_equal; delete random_unequal; } @@ -141,54 +140,51 @@ FixAtomSwap::~FixAtomSwap() void FixAtomSwap::options(int narg, char **arg) { - if (narg < 0) error->all(FLERR,"Illegal fix atom/swap command"); + if (narg < 0) error->all(FLERR, "Illegal fix atom/swap command"); - regionflag = 0; ke_flag = 1; semi_grand_flag = 0; nswaptypes = 0; nmutypes = 0; - iregion = -1; int iarg = 0; while (iarg < narg) { - if (strcmp(arg[iarg],"region") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix atom/swap command"); - iregion = domain->find_region(arg[iarg+1]); - if (iregion == -1) - error->all(FLERR,"Region ID for fix atom/swap does not exist"); - idregion = utils::strdup(arg[iarg+1]); - regionflag = 1; + if (strcmp(arg[iarg], "region") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix atom/swap command"); + region = domain->get_region_by_id(arg[iarg + 1]); + if (!region) error->all(FLERR, "Region {} for fix atom/swap does not exist", arg[iarg + 1]); + idregion = utils::strdup(arg[iarg + 1]); iarg += 2; - } else if (strcmp(arg[iarg],"ke") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix atom/swap command"); - ke_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "ke") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix atom/swap command"); + ke_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"semi-grand") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix atom/swap command"); - semi_grand_flag = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "semi-grand") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix atom/swap command"); + semi_grand_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"types") == 0) { - if (iarg+3 > narg) error->all(FLERR,"Illegal fix atom/swap command"); + } else if (strcmp(arg[iarg], "types") == 0) { + if (iarg + 3 > narg) error->all(FLERR, "Illegal fix atom/swap command"); iarg++; while (iarg < narg) { if (isalpha(arg[iarg][0])) break; - if (nswaptypes >= atom->ntypes) error->all(FLERR,"Illegal fix atom/swap command"); - type_list[nswaptypes] = utils::numeric(FLERR,arg[iarg],false,lmp); + if (nswaptypes >= atom->ntypes) error->all(FLERR, "Illegal fix atom/swap command"); + type_list[nswaptypes] = utils::numeric(FLERR, arg[iarg], false, lmp); nswaptypes++; iarg++; } - } else if (strcmp(arg[iarg],"mu") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix atom/swap command"); + } else if (strcmp(arg[iarg], "mu") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix atom/swap command"); iarg++; while (iarg < narg) { if (isalpha(arg[iarg][0])) break; nmutypes++; - if (nmutypes > atom->ntypes) error->all(FLERR,"Illegal fix atom/swap command"); - mu[nmutypes] = utils::numeric(FLERR,arg[iarg],false,lmp); + if (nmutypes > atom->ntypes) error->all(FLERR, "Illegal fix atom/swap command"); + mu[nmutypes] = utils::numeric(FLERR, arg[iarg], false, lmp); iarg++; } - } else error->all(FLERR,"Illegal fix atom/swap command"); + } else + error->all(FLERR, "Illegal fix atom/swap command"); } } @@ -205,36 +201,40 @@ int FixAtomSwap::setmask() void FixAtomSwap::init() { - auto id_pe = (char *) "thermo_pe"; - int ipe = modify->find_compute(id_pe); - c_pe = modify->compute[ipe]; + c_pe = modify->get_compute_by_id("thermo_pe"); int *type = atom->type; - if (nswaptypes < 2) - error->all(FLERR,"Must specify at least 2 types in fix atom/swap command"); + if (nswaptypes < 2) error->all(FLERR, "Must specify at least 2 types in fix atom/swap command"); if (semi_grand_flag) { if (nswaptypes != nmutypes) - error->all(FLERR,"Need nswaptypes mu values in fix atom/swap command"); + error->all(FLERR, "Need nswaptypes mu values in fix atom/swap command"); } else { if (nswaptypes != 2) - error->all(FLERR,"Only 2 types allowed when not using semi-grand in fix atom/swap command"); + error->all(FLERR, "Only 2 types allowed when not using semi-grand in fix atom/swap command"); if (nmutypes != 0) - error->all(FLERR,"Mu not allowed when not using semi-grand in fix atom/swap command"); + error->all(FLERR, "Mu not allowed when not using semi-grand in fix atom/swap command"); + } + + // set index and check validity of region + + if (idregion) { + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for fix setforce does not exist", idregion); } for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++) if (type_list[iswaptype] <= 0 || type_list[iswaptype] > atom->ntypes) - error->all(FLERR,"Invalid atom type in fix atom/swap command"); + error->all(FLERR, "Invalid atom type in fix atom/swap command"); // this is only required for non-semi-grand // in which case, nswaptypes = 2 if (atom->q_flag && !semi_grand_flag) { - double qmax,qmin; - int firstall,first; - memory->create(qtype,nswaptypes,"atom/swap:qtype"); + double qmax, qmin; + int firstall, first; + memory->create(qtype, nswaptypes, "atom/swap:qtype"); for (int iswaptype = 0; iswaptype < nswaptypes; iswaptype++) { first = 1; for (int i = 0; i < atom->nlocal; i++) { @@ -244,24 +244,26 @@ void FixAtomSwap::init() qtype[iswaptype] = atom->q[i]; first = 0; } else if (qtype[iswaptype] != atom->q[i]) - error->one(FLERR,"All atoms of a swapped type must have the same charge."); + error->one(FLERR, "All atoms of a swapped type must have the same charge."); } } } - MPI_Allreduce(&first,&firstall,1,MPI_INT,MPI_MIN,world); - if (firstall) error->all(FLERR,"At least one atom of each swapped type must be present to define charges."); + MPI_Allreduce(&first, &firstall, 1, MPI_INT, MPI_MIN, world); + if (firstall) + error->all(FLERR, + "At least one atom of each swapped type must be present to define charges."); if (first) qtype[iswaptype] = -DBL_MAX; - MPI_Allreduce(&qtype[iswaptype],&qmax,1,MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(&qtype[iswaptype], &qmax, 1, MPI_DOUBLE, MPI_MAX, world); if (first) qtype[iswaptype] = DBL_MAX; - MPI_Allreduce(&qtype[iswaptype],&qmin,1,MPI_DOUBLE,MPI_MIN,world); - if (qmax != qmin) error->all(FLERR,"All atoms of a swapped type must have same charge."); + MPI_Allreduce(&qtype[iswaptype], &qmin, 1, MPI_DOUBLE, MPI_MIN, world); + if (qmax != qmin) error->all(FLERR, "All atoms of a swapped type must have same charge."); } } - memory->create(sqrt_mass_ratio,atom->ntypes+1,atom->ntypes+1,"atom/swap:sqrt_mass_ratio"); + memory->create(sqrt_mass_ratio, atom->ntypes + 1, atom->ntypes + 1, "atom/swap:sqrt_mass_ratio"); for (int itype = 1; itype <= atom->ntypes; itype++) for (int jtype = 1; jtype <= atom->ntypes; jtype++) - sqrt_mass_ratio[itype][jtype] = sqrt(atom->mass[itype]/atom->mass[jtype]); + sqrt_mass_ratio[itype][jtype] = sqrt(atom->mass[itype] / atom->mass[jtype]); // check to see if itype and jtype cutoffs are the same // if not, reneighboring will be needed between swaps @@ -286,10 +288,9 @@ void FixAtomSwap::init() if ((mask[i] == groupbit) && (mask[i] && firstgroupbit)) flag = 1; int flagall; - MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&flag, &flagall, 1, MPI_INT, MPI_SUM, world); - if (flagall) - error->all(FLERR,"Cannot do atom/swap on atoms in atom_modify first group"); + if (flagall) error->all(FLERR, "Cannot do atom/swap on atoms in atom_modify first group"); } } @@ -309,7 +310,7 @@ void FixAtomSwap::pre_exchange() domain->pbc(); comm->exchange(); comm->borders(); - if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(1); @@ -353,14 +354,14 @@ int FixAtomSwap::attempt_semi_grand() // pick a random atom and perform swap - int itype,jtype,jswaptype; + int itype, jtype, jswaptype; int i = pick_semi_grand_atom(); if (i >= 0) { - jswaptype = static_cast (nswaptypes*random_unequal->uniform()); + jswaptype = static_cast(nswaptypes * random_unequal->uniform()); jtype = type_list[jswaptype]; itype = atom->type[i]; while (itype == jtype) { - jswaptype = static_cast (nswaptypes*random_unequal->uniform()); + jswaptype = static_cast(nswaptypes * random_unequal->uniform()); jtype = type_list[jswaptype]; } atom->type[i] = jtype; @@ -374,7 +375,7 @@ int FixAtomSwap::attempt_semi_grand() if (domain->triclinic) domain->x2lamda(atom->nlocal); comm->exchange(); comm->borders(); - if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(1); } else { @@ -389,11 +390,11 @@ int FixAtomSwap::attempt_semi_grand() int success = 0; if (i >= 0) if (random_unequal->uniform() < - exp(beta*(energy_before - energy_after - + mu[jtype] - mu[itype]))) success = 1; + exp(beta * (energy_before - energy_after + mu[jtype] - mu[itype]))) + success = 1; int success_all = 0; - MPI_Allreduce(&success,&success_all,1,MPI_INT,MPI_MAX,world); + MPI_Allreduce(&success, &success_all, 1, MPI_INT, MPI_MAX, world); // swap accepted, return 1 @@ -460,7 +461,7 @@ int FixAtomSwap::attempt_swap() domain->pbc(); comm->exchange(); comm->borders(); - if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost); + if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost); if (modify->n_pre_neighbor) modify->pre_neighbor(); neighbor->build(1); } else { @@ -474,8 +475,7 @@ int FixAtomSwap::attempt_swap() // swap accepted, return 1 // if ke_flag, rescale atom velocities - if (random_equal->uniform() < - exp(beta*(energy_before - energy_after))) { + if (random_equal->uniform() < exp(beta * (energy_before - energy_after))) { update_swap_atoms_list(); if (ke_flag) { if (i >= 0) { @@ -521,16 +521,16 @@ double FixAtomSwap::energy_full() if (modify->n_pre_force) modify->pre_force(vflag); - if (force->pair) force->pair->compute(eflag,vflag); + if (force->pair) force->pair->compute(eflag, vflag); if (atom->molecular != Atom::ATOMIC) { - if (force->bond) force->bond->compute(eflag,vflag); - if (force->angle) force->angle->compute(eflag,vflag); - if (force->dihedral) force->dihedral->compute(eflag,vflag); - if (force->improper) force->improper->compute(eflag,vflag); + if (force->bond) force->bond->compute(eflag, vflag); + if (force->angle) force->angle->compute(eflag, vflag); + if (force->dihedral) force->dihedral->compute(eflag, vflag); + if (force->improper) force->improper->compute(eflag, vflag); } - if (force->kspace) force->kspace->compute(eflag,vflag); + if (force->kspace) force->kspace->compute(eflag, vflag); if (modify->n_post_force) modify->post_force(vflag); @@ -546,9 +546,8 @@ double FixAtomSwap::energy_full() int FixAtomSwap::pick_semi_grand_atom() { int i = -1; - int iwhichglobal = static_cast (nswap*random_equal->uniform()); - if ((iwhichglobal >= nswap_before) && - (iwhichglobal < nswap_before + nswap_local)) { + int iwhichglobal = static_cast(nswap * random_equal->uniform()); + if ((iwhichglobal >= nswap_before) && (iwhichglobal < nswap_before + nswap_local)) { int iwhichlocal = iwhichglobal - nswap_before; i = local_swap_atom_list[iwhichlocal]; } @@ -562,9 +561,8 @@ int FixAtomSwap::pick_semi_grand_atom() int FixAtomSwap::pick_i_swap_atom() { int i = -1; - int iwhichglobal = static_cast (niswap*random_equal->uniform()); - if ((iwhichglobal >= niswap_before) && - (iwhichglobal < niswap_before + niswap_local)) { + int iwhichglobal = static_cast(niswap * random_equal->uniform()); + if ((iwhichglobal >= niswap_before) && (iwhichglobal < niswap_before + niswap_local)) { int iwhichlocal = iwhichglobal - niswap_before; i = local_swap_iatom_list[iwhichlocal]; } @@ -578,9 +576,8 @@ int FixAtomSwap::pick_i_swap_atom() int FixAtomSwap::pick_j_swap_atom() { int j = -1; - int jwhichglobal = static_cast (njswap*random_equal->uniform()); - if ((jwhichglobal >= njswap_before) && - (jwhichglobal < njswap_before + njswap_local)) { + int jwhichglobal = static_cast(njswap * random_equal->uniform()); + if ((jwhichglobal >= njswap_before) && (jwhichglobal < njswap_before + njswap_local)) { int jwhichlocal = jwhichglobal - njswap_before; j = local_swap_jatom_list[jwhichlocal]; } @@ -600,16 +597,15 @@ void FixAtomSwap::update_semi_grand_atoms_list() if (atom->nmax > atom_swap_nmax) { memory->sfree(local_swap_atom_list); atom_swap_nmax = atom->nmax; - local_swap_atom_list = (int *) memory->smalloc(atom_swap_nmax*sizeof(int), - "MCSWAP:local_swap_atom_list"); + local_swap_atom_list = + (int *) memory->smalloc(atom_swap_nmax * sizeof(int), "MCSWAP:local_swap_atom_list"); } nswap_local = 0; - if (regionflag) { - + if (region) { for (int i = 0; i < nlocal; i++) { - if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) == 1) { + if (region->match(x[i][0], x[i][1], x[i][2]) == 1) { if (atom->mask[i] & groupbit) { int itype = atom->type[i]; int iswaptype; @@ -625,23 +621,22 @@ void FixAtomSwap::update_semi_grand_atoms_list() } else { for (int i = 0; i < nlocal; i++) { if (atom->mask[i] & groupbit) { - int itype = atom->type[i]; - int iswaptype; - for (iswaptype = 0; iswaptype < nswaptypes; iswaptype++) - if (itype == type_list[iswaptype]) break; - if (iswaptype == nswaptypes) continue; + int itype = atom->type[i]; + int iswaptype; + for (iswaptype = 0; iswaptype < nswaptypes; iswaptype++) + if (itype == type_list[iswaptype]) break; + if (iswaptype == nswaptypes) continue; local_swap_atom_list[nswap_local] = i; nswap_local++; } } } - MPI_Allreduce(&nswap_local,&nswap,1,MPI_INT,MPI_SUM,world); - MPI_Scan(&nswap_local,&nswap_before,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&nswap_local, &nswap, 1, MPI_INT, MPI_SUM, world); + MPI_Scan(&nswap_local, &nswap_before, 1, MPI_INT, MPI_SUM, world); nswap_before -= nswap_local; } - /* ---------------------------------------------------------------------- update the list of gas atoms ------------------------------------------------------------------------- */ @@ -656,24 +651,24 @@ void FixAtomSwap::update_swap_atoms_list() memory->sfree(local_swap_iatom_list); memory->sfree(local_swap_jatom_list); atom_swap_nmax = atom->nmax; - local_swap_iatom_list = (int *) memory->smalloc(atom_swap_nmax*sizeof(int), - "MCSWAP:local_swap_iatom_list"); - local_swap_jatom_list = (int *) memory->smalloc(atom_swap_nmax*sizeof(int), - "MCSWAP:local_swap_jatom_list"); + local_swap_iatom_list = + (int *) memory->smalloc(atom_swap_nmax * sizeof(int), "MCSWAP:local_swap_iatom_list"); + local_swap_jatom_list = + (int *) memory->smalloc(atom_swap_nmax * sizeof(int), "MCSWAP:local_swap_jatom_list"); } niswap_local = 0; njswap_local = 0; - if (regionflag) { + if (region) { for (int i = 0; i < nlocal; i++) { - if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) == 1) { + if (region->match(x[i][0], x[i][1], x[i][2]) == 1) { if (atom->mask[i] & groupbit) { - if (type[i] == type_list[0]) { + if (type[i] == type_list[0]) { local_swap_iatom_list[niswap_local] = i; niswap_local++; - } else if (type[i] == type_list[1]) { + } else if (type[i] == type_list[1]) { local_swap_jatom_list[njswap_local] = i; njswap_local++; } @@ -684,10 +679,10 @@ void FixAtomSwap::update_swap_atoms_list() } else { for (int i = 0; i < nlocal; i++) { if (atom->mask[i] & groupbit) { - if (type[i] == type_list[0]) { + if (type[i] == type_list[0]) { local_swap_iatom_list[niswap_local] = i; niswap_local++; - } else if (type[i] == type_list[1]) { + } else if (type[i] == type_list[1]) { local_swap_jatom_list[njswap_local] = i; njswap_local++; } @@ -695,12 +690,12 @@ void FixAtomSwap::update_swap_atoms_list() } } - MPI_Allreduce(&niswap_local,&niswap,1,MPI_INT,MPI_SUM,world); - MPI_Scan(&niswap_local,&niswap_before,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&niswap_local, &niswap, 1, MPI_INT, MPI_SUM, world); + MPI_Scan(&niswap_local, &niswap_before, 1, MPI_INT, MPI_SUM, world); niswap_before -= niswap_local; - MPI_Allreduce(&njswap_local,&njswap,1,MPI_INT,MPI_SUM,world); - MPI_Scan(&njswap_local,&njswap_before,1,MPI_INT,MPI_SUM,world); + MPI_Allreduce(&njswap_local, &njswap, 1, MPI_INT, MPI_SUM, world); + MPI_Scan(&njswap_local, &njswap_before, 1, MPI_INT, MPI_SUM, world); njswap_before -= njswap_local; } @@ -708,7 +703,7 @@ void FixAtomSwap::update_swap_atoms_list() int FixAtomSwap::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/) { - int i,j,m; + int i, j, m; int *type = atom->type; double *q = atom->q; @@ -735,7 +730,7 @@ int FixAtomSwap::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag void FixAtomSwap::unpack_forward_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; int *type = atom->type; double *q = atom->q; @@ -745,12 +740,11 @@ void FixAtomSwap::unpack_forward_comm(int n, int first, double *buf) if (atom->q_flag) { for (i = first; i < last; i++) { - type[i] = static_cast (buf[m++]); + type[i] = static_cast(buf[m++]); q[i] = buf[m++]; } } else { - for (i = first; i < last; i++) - type[i] = static_cast (buf[m++]); + for (i = first; i < last; i++) type[i] = static_cast(buf[m++]); } } @@ -771,7 +765,7 @@ double FixAtomSwap::compute_vector(int n) double FixAtomSwap::memory_usage() { - double bytes = (double)atom_swap_nmax * sizeof(int); + double bytes = (double) atom_swap_nmax * sizeof(int); return bytes; } @@ -792,8 +786,8 @@ void FixAtomSwap::write_restart(FILE *fp) if (comm->me == 0) { int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(list,sizeof(double),n,fp); + fwrite(&size, sizeof(int), 1, fp); + fwrite(list, sizeof(double), n, fp); } } @@ -806,10 +800,10 @@ void FixAtomSwap::restart(char *buf) int n = 0; auto list = (double *) buf; - seed = static_cast (list[n++]); + seed = static_cast(list[n++]); random_equal->reset(seed); - seed = static_cast (list[n++]); + seed = static_cast(list[n++]); random_unequal->reset(seed); next_reneighbor = (bigint) ubuf(list[n++]).i; @@ -819,5 +813,5 @@ void FixAtomSwap::restart(char *buf) bigint ntimestep_restart = (bigint) ubuf(list[n++]).i; if (ntimestep_restart != update->ntimestep) - error->all(FLERR,"Must not reset timestep when restarting fix atom/swap"); + error->all(FLERR, "Must not reset timestep when restarting fix atom/swap"); } diff --git a/src/MC/fix_atom_swap.h b/src/MC/fix_atom_swap.h index e66b383031..c4e280cf06 100644 --- a/src/MC/fix_atom_swap.h +++ b/src/MC/fix_atom_swap.h @@ -49,8 +49,7 @@ class FixAtomSwap : public Fix { int nswap; // # of swap atoms on all procs int nswap_local; // # of swap atoms on this proc int nswap_before; // # of swap atoms on procs < this proc - int regionflag; // 0 = anywhere in box, 1 = specific region - int iregion; // swap region + class Region *region; // swap region char *idregion; // swap region id int nswaptypes, nmutypes; diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index 6937b90202..f6ed234c07 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -68,7 +68,7 @@ enum{NONE,MOVEATOM,MOVEMOL}; // movemode FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - idregion(nullptr), full_flag(false), ngroups(0), groupstrings(nullptr), ngrouptypes(0), + region(nullptr), idregion(nullptr), full_flag(false), groupstrings(nullptr), grouptypestrings(nullptr), grouptypebits(nullptr), grouptypes(nullptr), local_gas_list(nullptr), molcoords(nullptr), molq(nullptr), molimage(nullptr), random_equal(nullptr), random_unequal(nullptr), fixrigid(nullptr), fixshake(nullptr), idrigid(nullptr), idshake(nullptr) @@ -87,6 +87,9 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) : restart_global = 1; time_depend = 1; + ngroups = 0; + ngrouptypes = 0; + // required args nevery = utils::inumeric(FLERR,arg[3],false,lmp); @@ -122,18 +125,18 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) : region_xlo = region_xhi = region_ylo = region_yhi = region_zlo = region_zhi = 0.0; - if (regionflag) { - if (domain->regions[iregion]->bboxflag == 0) + if (region) { + if (region->bboxflag == 0) error->all(FLERR,"Fix gcmc region does not support a bounding box"); - if (domain->regions[iregion]->dynamic_check()) + if (region->dynamic_check()) error->all(FLERR,"Fix gcmc region cannot be dynamic"); - region_xlo = domain->regions[iregion]->extent_xlo; - region_xhi = domain->regions[iregion]->extent_xhi; - region_ylo = domain->regions[iregion]->extent_ylo; - region_yhi = domain->regions[iregion]->extent_yhi; - region_zlo = domain->regions[iregion]->extent_zlo; - region_zhi = domain->regions[iregion]->extent_zhi; + region_xlo = region->extent_xlo; + region_xhi = region->extent_xhi; + region_ylo = region->extent_ylo; + region_yhi = region->extent_yhi; + region_zlo = region->extent_zlo; + region_zhi = region->extent_zhi; if (region_xlo < domain->boxlo[0] || region_xhi > domain->boxhi[0] || region_ylo < domain->boxlo[1] || region_yhi > domain->boxhi[1] || @@ -149,15 +152,14 @@ FixGCMC::FixGCMC(LAMMPS *lmp, int narg, char **arg) : coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - if (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) != 0) + if (region->match(coord[0],coord[1],coord[2]) != 0) inside++; } - double max_region_volume = (region_xhi - region_xlo)* - (region_yhi - region_ylo)*(region_zhi - region_zlo); + double max_region_volume = (region_xhi - region_xlo) * + (region_yhi - region_ylo) * (region_zhi - region_zlo); - region_volume = max_region_volume*static_cast (inside)/ - static_cast (attempts); + region_volume = max_region_volume * static_cast(inside) / static_cast(attempts); } // error check and further setup for exchmode = EXCHMOL @@ -241,8 +243,6 @@ void FixGCMC::options(int narg, char **arg) pmolrotate = 0.0; pmctot = 0.0; max_rotation_angle = 10*MY_PI/180; - regionflag = 0; - iregion = -1; region_volume = 0; max_region_attempts = 1000; molecule_group = 0; @@ -300,11 +300,10 @@ void FixGCMC::options(int narg, char **arg) iarg += 4; } else if (strcmp(arg[iarg],"region") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command"); - iregion = domain->find_region(arg[iarg+1]); - if (iregion == -1) - error->all(FLERR,"Region ID for fix gcmc does not exist"); + region = domain->get_region_by_id(arg[iarg+1]); + if (!region) + error->all(FLERR,"Region {} for fix gcmc does not exist",arg[iarg+1]); idregion = utils::strdup(arg[iarg+1]); - regionflag = 1; iarg += 2; } else if (strcmp(arg[iarg],"maxangle") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command"); @@ -397,7 +396,7 @@ void FixGCMC::options(int narg, char **arg) FixGCMC::~FixGCMC() { - if (regionflag) delete [] idregion; + delete[] idregion; delete random_equal; delete random_unequal; @@ -406,12 +405,12 @@ FixGCMC::~FixGCMC() memory->destroy(molq); memory->destroy(molimage); - delete [] idrigid; - delete [] idshake; + delete[] idrigid; + delete[] idshake; if (ngroups > 0) { for (int igroup = 0; igroup < ngroups; igroup++) - delete [] groupstrings[igroup]; + delete[] groupstrings[igroup]; memory->sfree(groupstrings); } @@ -419,7 +418,7 @@ FixGCMC::~FixGCMC() memory->destroy(grouptypes); memory->destroy(grouptypebits); for (int igroup = 0; igroup < ngrouptypes; igroup++) - delete [] grouptypestrings[igroup]; + delete[] grouptypestrings[igroup]; memory->sfree(grouptypestrings); } if (full_flag && group) { @@ -443,6 +442,13 @@ int FixGCMC::setmask() void FixGCMC::init() { + // set index and check validity of region + + if (idregion) { + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for fix gcmc does not exist", idregion); + } + triclinic = domain->triclinic; // set probabilities for MC moves @@ -719,7 +725,7 @@ void FixGCMC::pre_exchange() subhi = domain->subhi; } - if (regionflag) volume = region_volume; + if (region) volume = region_volume; else volume = domain->xprd * domain->yprd * domain->zprd; if (triclinic) domain->x2lamda(atom->nlocal); @@ -801,8 +807,7 @@ void FixGCMC::attempt_atomic_translation() double **x = atom->x; double energy_before = energy(i,ngcmc_type,-1,x[i]); if (overlap_flag && energy_before > MAXENERGYTEST) - error->warning(FLERR,"Energy of old configuration in " - "fix gcmc is > MAXENERGYTEST."); + error->warning(FLERR,"Energy of old configuration in fix gcmc is > MAXENERGYTEST."); double rsq = 1.1; double rx,ry,rz; rx = ry = rz = 0.0; @@ -816,8 +821,8 @@ void FixGCMC::attempt_atomic_translation() coord[0] = x[i][0] + displace*rx; coord[1] = x[i][1] + displace*ry; coord[2] = x[i][2] + displace*rz; - if (regionflag) { - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + if (region) { + while (region->match(coord[0],coord[1],coord[2]) == 0) { rsq = 1.1; while (rsq > 1.0) { rx = 2*random_unequal->uniform() - 1.0; @@ -913,12 +918,12 @@ void FixGCMC::attempt_atomic_insertion() // pick coordinates for insertion point double coord[3]; - if (regionflag) { + if (region) { int region_attempt = 0; coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + while (region->match(coord[0],coord[1],coord[2]) == 0) { coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); @@ -1043,7 +1048,7 @@ void FixGCMC::attempt_molecule_translation() com_displace[1] = displace*ry; com_displace[2] = displace*rz; - if (regionflag) { + if (region) { int *mask = atom->mask; for (int i = 0; i < atom->nlocal; i++) { if (atom->molecule[i] == translation_molecule) { @@ -1058,7 +1063,7 @@ void FixGCMC::attempt_molecule_translation() coord[0] = com[0] + displace*rx; coord[1] = com[1] + displace*ry; coord[2] = com[2] + displace*rz; - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + while (region->match(coord[0],coord[1],coord[2]) == 0) { rsq = 1.1; while (rsq > 1.0) { rx = 2*random_equal->uniform() - 1.0; @@ -1266,7 +1271,7 @@ void FixGCMC::attempt_molecule_insertion() if (ngas >= max_ngas) return; double com_coord[3]; - if (regionflag) { + if (region) { int region_attempt = 0; com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); @@ -1274,7 +1279,7 @@ void FixGCMC::attempt_molecule_insertion() (region_yhi-region_ylo); com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - while (domain->regions[iregion]->match(com_coord[0],com_coord[1], + while (region->match(com_coord[0],com_coord[1], com_coord[2]) == 0) { com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); @@ -1485,8 +1490,8 @@ void FixGCMC::attempt_atomic_translation_full() coord[0] = x[i][0] + displace*rx; coord[1] = x[i][1] + displace*ry; coord[2] = x[i][2] + displace*rz; - if (regionflag) { - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + if (region) { + while (region->match(coord[0],coord[1],coord[2]) == 0) { rsq = 1.1; while (rsq > 1.0) { rx = 2*random_unequal->uniform() - 1.0; @@ -1602,12 +1607,12 @@ void FixGCMC::attempt_atomic_insertion_full() double energy_before = energy_stored; double coord[3]; - if (regionflag) { + if (region) { int region_attempt = 0; coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + while (region->match(coord[0],coord[1],coord[2]) == 0) { coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); @@ -1726,7 +1731,7 @@ void FixGCMC::attempt_molecule_translation_full() com_displace[1] = displace*ry; com_displace[2] = displace*rz; - if (regionflag) { + if (region) { int *mask = atom->mask; for (int i = 0; i < atom->nlocal; i++) { if (atom->molecule[i] == translation_molecule) { @@ -1741,7 +1746,7 @@ void FixGCMC::attempt_molecule_translation_full() coord[0] = com[0] + displace*rx; coord[1] = com[1] + displace*ry; coord[2] = com[2] + displace*rz; - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + while (region->match(coord[0],coord[1],coord[2]) == 0) { rsq = 1.1; while (rsq > 1.0) { rx = 2*random_equal->uniform() - 1.0; @@ -1998,7 +2003,7 @@ void FixGCMC::attempt_molecule_insertion_full() int nlocalprev = atom->nlocal; double com_coord[3]; - if (regionflag) { + if (region) { int region_attempt = 0; com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); @@ -2006,7 +2011,7 @@ void FixGCMC::attempt_molecule_insertion_full() (region_yhi-region_ylo); com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - while (domain->regions[iregion]->match(com_coord[0],com_coord[1], + while (region->match(com_coord[0],com_coord[1], com_coord[2]) == 0) { com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); @@ -2408,7 +2413,7 @@ void FixGCMC::update_gas_atoms_list() ngas_local = 0; - if (regionflag) { + if (region) { if (exchmode == EXCHMOL || movemode == MOVEMOL) { @@ -2441,7 +2446,7 @@ void FixGCMC::update_gas_atoms_list() for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (domain->regions[iregion]->match(comx[molecule[i]], + if (region->match(comx[molecule[i]], comy[molecule[i]],comz[molecule[i]]) == 1) { local_gas_list[ngas_local] = i; ngas_local++; @@ -2454,7 +2459,7 @@ void FixGCMC::update_gas_atoms_list() } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) == 1) { + if (region->match(x[i][0],x[i][1],x[i][2]) == 1) { local_gas_list[ngas_local] = i; ngas_local++; } diff --git a/src/MC/fix_gcmc.h b/src/MC/fix_gcmc.h index 79ce515b20..bddf0c425a 100644 --- a/src/MC/fix_gcmc.h +++ b/src/MC/fix_gcmc.h @@ -65,17 +65,16 @@ class FixGCMC : public Fix { int ngcmc_type, nevery, seed; int ncycles, nexchanges, nmcmoves; double patomtrans, pmoltrans, pmolrotate, pmctot; - int ngas; // # of gas atoms on all procs - int ngas_local; // # of gas atoms on this proc - int ngas_before; // # of gas atoms on procs < this proc - int exchmode; // exchange ATOM or MOLECULE - int movemode; // move ATOM or MOLECULE - int regionflag; // 0 = anywhere in box, 1 = specific region - int iregion; // gcmc region - char *idregion; // gcmc region id - bool pressure_flag; // true if user specified reservoir pressure - bool charge_flag; // true if user specified atomic charge - bool full_flag; // true if doing full system energy calculations + int ngas; // # of gas atoms on all procs + int ngas_local; // # of gas atoms on this proc + int ngas_before; // # of gas atoms on procs < this proc + int exchmode; // exchange ATOM or MOLECULE + int movemode; // move ATOM or MOLECULE + class Region *region; // gcmc region + char *idregion; // gcmc region id + bool pressure_flag; // true if user specified reservoir pressure + bool charge_flag; // true if user specified atomic charge + bool full_flag; // true if doing full system energy calculations int natoms_per_molecule; // number of atoms in each inserted molecule int nmaxmolatoms; // number of atoms allocated for molecule arrays diff --git a/src/MC/fix_widom.cpp b/src/MC/fix_widom.cpp index bb20313ea8..0adabe5eae 100644 --- a/src/MC/fix_widom.cpp +++ b/src/MC/fix_widom.cpp @@ -46,7 +46,6 @@ #include #include -using namespace std; using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; @@ -59,8 +58,8 @@ enum{EXCHATOM,EXCHMOL}; // exchmode FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - idregion(nullptr), full_flag(false), local_gas_list(nullptr), molcoords(nullptr), - molq(nullptr), molimage(nullptr), random_equal(nullptr) + region(nullptr), idregion(nullptr), full_flag(false), local_gas_list(nullptr), + molcoords(nullptr),molq(nullptr), molimage(nullptr), random_equal(nullptr) { if (narg < 8) error->all(FLERR,"Illegal fix widom command"); @@ -76,8 +75,6 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) : restart_global = 1; time_depend = 1; - //ave_widom_chemical_potential = 0; - // required args nevery = utils::inumeric(FLERR,arg[3],false,lmp); @@ -104,18 +101,18 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) : region_xlo = region_xhi = region_ylo = region_yhi = region_zlo = region_zhi = 0.0; - if (regionflag) { - if (domain->regions[iregion]->bboxflag == 0) + if (region) { + if (region->bboxflag == 0) error->all(FLERR,"Fix widom region does not support a bounding box"); - if (domain->regions[iregion]->dynamic_check()) + if (region->dynamic_check()) error->all(FLERR,"Fix widom region cannot be dynamic"); - region_xlo = domain->regions[iregion]->extent_xlo; - region_xhi = domain->regions[iregion]->extent_xhi; - region_ylo = domain->regions[iregion]->extent_ylo; - region_yhi = domain->regions[iregion]->extent_yhi; - region_zlo = domain->regions[iregion]->extent_zlo; - region_zhi = domain->regions[iregion]->extent_zhi; + region_xlo = region->extent_xlo; + region_xhi = region->extent_xhi; + region_ylo = region->extent_ylo; + region_yhi = region->extent_yhi; + region_zlo = region->extent_zlo; + region_zhi = region->extent_zhi; if (region_xlo < domain->boxlo[0] || region_xhi > domain->boxhi[0] || region_ylo < domain->boxlo[1] || region_yhi > domain->boxhi[1] || @@ -131,15 +128,14 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) : coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - if (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) != 0) + if (region->match(coord[0],coord[1],coord[2]) != 0) inside++; } - double max_region_volume = (region_xhi - region_xlo)* - (region_yhi - region_ylo)*(region_zhi - region_zlo); + double max_region_volume = (region_xhi - region_xlo) * + (region_yhi - region_ylo) * (region_zhi - region_zlo); - region_volume = max_region_volume*static_cast (inside)/ - static_cast (attempts); + region_volume = max_region_volume * static_cast(inside) / static_cast(attempts); } // error check and further setup for exchmode = EXCHMOL @@ -191,8 +187,6 @@ void FixWidom::options(int narg, char **arg) // defaults exchmode = EXCHATOM; - regionflag = 0; - iregion = -1; region_volume = 0; max_region_attempts = 1000; molecule_group = 0; @@ -221,11 +215,10 @@ void FixWidom::options(int narg, char **arg) iarg += 2; } else if (strcmp(arg[iarg],"region") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix widom command"); - iregion = domain->find_region(arg[iarg+1]); - if (iregion == -1) - error->all(FLERR,"Region ID for fix widom does not exist"); + region = domain->get_region_by_id(arg[iarg+1]); + if (!region) + error->all(FLERR,"Region {} for fix widom does not exist",arg[iarg+1]); idregion = utils::strdup(arg[iarg+1]); - regionflag = 1; iarg += 2; } else if (strcmp(arg[iarg],"charge") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix widom command"); @@ -247,7 +240,7 @@ void FixWidom::options(int narg, char **arg) FixWidom::~FixWidom() { - if (regionflag) delete [] idregion; + delete[] idregion; delete random_equal; memory->destroy(local_gas_list); @@ -271,11 +264,18 @@ int FixWidom::setmask() void FixWidom::init() { + // set index and check validity of region + + if (idregion) { + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for fix widom does not exist", idregion); + } + triclinic = domain->triclinic; ave_widom_chemical_potential = 0; - if (regionflag) volume = region_volume; + if (region) volume = region_volume; else volume = domain->xprd * domain->yprd * domain->zprd; // decide whether to switch to the full_energy option @@ -283,8 +283,8 @@ void FixWidom::init() if ((force->kspace) || (force->pair == nullptr) || (force->pair->single_enable == 0) || - (force->pair_match("hybrid",0)) || - (force->pair_match("eam",0)) || + (force->pair_match("^hybrid",0)) || + (force->pair_match("^eam",0)) || (force->pair->tail_flag)) { full_flag = true; if (comm->me == 0) @@ -434,7 +434,7 @@ void FixWidom::pre_exchange() subhi = domain->subhi; } - if (regionflag) volume = region_volume; + if (region) volume = region_volume; else volume = domain->xprd * domain->yprd * domain->zprd; if (triclinic) domain->x2lamda(atom->nlocal); @@ -486,12 +486,12 @@ void FixWidom::attempt_atomic_insertion() // pick coordinates for insertion point - if (regionflag) { + if (region) { int region_attempt = 0; coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + while (region->match(coord[0],coord[1],coord[2]) == 0) { coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); @@ -562,7 +562,7 @@ void FixWidom::attempt_molecule_insertion() for (int imove = 0; imove < ninsertions; imove++) { - if (regionflag) { + if (region) { int region_attempt = 0; com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); @@ -570,7 +570,7 @@ void FixWidom::attempt_molecule_insertion() (region_yhi-region_ylo); com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - while (domain->regions[iregion]->match(com_coord[0],com_coord[1], + while (region->match(com_coord[0],com_coord[1], com_coord[2]) == 0) { com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); @@ -688,12 +688,12 @@ void FixWidom::attempt_atomic_insertion_full() for (int imove = 0; imove < ninsertions; imove++) { - if (regionflag) { + if (region) { int region_attempt = 0; coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) { + while (region->match(coord[0],coord[1],coord[2]) == 0) { coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo); coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); @@ -801,7 +801,7 @@ void FixWidom::attempt_molecule_insertion_full() for (int imove = 0; imove < ninsertions; imove++) { double com_coord[3]; - if (regionflag) { + if (region) { int region_attempt = 0; com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); @@ -809,7 +809,7 @@ void FixWidom::attempt_molecule_insertion_full() (region_yhi-region_ylo); com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo); - while (domain->regions[iregion]->match(com_coord[0],com_coord[1], + while (region->match(com_coord[0],com_coord[1], com_coord[2]) == 0) { com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo); @@ -1081,7 +1081,7 @@ void FixWidom::update_gas_atoms_list() ngas_local = 0; - if (regionflag) { + if (region) { if (exchmode == EXCHMOL) { @@ -1114,7 +1114,7 @@ void FixWidom::update_gas_atoms_list() for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (domain->regions[iregion]->match(comx[molecule[i]], + if (region->match(comx[molecule[i]], comy[molecule[i]],comz[molecule[i]]) == 1) { local_gas_list[ngas_local] = i; ngas_local++; @@ -1127,7 +1127,7 @@ void FixWidom::update_gas_atoms_list() } else { for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) == 1) { + if (region->match(x[i][0],x[i][1],x[i][2]) == 1) { local_gas_list[ngas_local] = i; ngas_local++; } diff --git a/src/MC/fix_widom.h b/src/MC/fix_widom.h index d2240fd6c8..ec273f73a7 100644 --- a/src/MC/fix_widom.h +++ b/src/MC/fix_widom.h @@ -53,15 +53,14 @@ class FixWidom : public Fix { int exclusion_group, exclusion_group_bit; int nwidom_type, nevery, seed; int ninsertions; - int ngas; // # of gas atoms on all procs - int ngas_local; // # of gas atoms on this proc - int exchmode; // exchange ATOM or MOLECULE - int movemode; // move ATOM or MOLECULE - int regionflag; // 0 = anywhere in box, 1 = specific region - int iregion; // widom region - char *idregion; // widom region id - bool charge_flag; // true if user specified atomic charge - bool full_flag; // true if doing full system energy calculations + int ngas; // # of gas atoms on all procs + int ngas_local; // # of gas atoms on this proc + int exchmode; // exchange ATOM or MOLECULE + int movemode; // move ATOM or MOLECULE + class Region *region; // widom region + char *idregion; // widom region id + bool charge_flag; // true if user specified atomic charge + bool full_flag; // true if doing full system energy calculations int natoms_per_molecule; // number of atoms in each inserted molecule int nmaxmolatoms; // number of atoms allocated for molecule arrays diff --git a/src/PLUGIN/plugin.cpp b/src/PLUGIN/plugin.cpp index f8dfb8af22..4d608304a0 100644 --- a/src/PLUGIN/plugin.cpp +++ b/src/PLUGIN/plugin.cpp @@ -389,7 +389,7 @@ void plugin_unload(const char *style, const char *name, LAMMPS *lmp) if (found != region_map->end()) region_map->erase(name); for (auto iregion : lmp->domain->get_region_by_style(name)) - lmp->domain->delete_region(iregion->id); + lmp->domain->delete_region(iregion); } else if (pstyle == "command") { diff --git a/src/REPLICA/hyper.cpp b/src/REPLICA/hyper.cpp index 00c594e7a3..9fa4910473 100644 --- a/src/REPLICA/hyper.cpp +++ b/src/REPLICA/hyper.cpp @@ -155,12 +155,11 @@ void Hyper::command(int narg, char **arg) // cannot use hyper with time-dependent fixes or regions - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->time_depend) - error->all(FLERR,"Cannot use hyper with a time-dependent fix defined"); + for (auto ifix : modify->get_fix_list()) + if (ifix->time_depend) error->all(FLERR,"Cannot use hyper with a time-dependent fix defined"); - for (int i = 0; i < domain->nregion; i++) - if (domain->regions[i]->dynamic_check()) + for (auto reg : domain->get_region_list()) + if (reg->dynamic_check()) error->all(FLERR,"Cannot use hyper with a time-dependent region defined"); // perform hyperdynamics simulation diff --git a/src/REPLICA/prd.cpp b/src/REPLICA/prd.cpp index 671bb49b64..79d09080ec 100644 --- a/src/REPLICA/prd.cpp +++ b/src/REPLICA/prd.cpp @@ -229,12 +229,11 @@ void PRD::command(int narg, char **arg) // cannot use PRD with time-dependent fixes or regions - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->time_depend) - error->all(FLERR,"Cannot use PRD with a time-dependent fix defined"); + for (auto ifix : modify->get_fix_list()) + if (ifix->time_depend) error->all(FLERR,"Cannot use PRD with a time-dependent fix defined"); - for (int i = 0; i < domain->nregion; i++) - if (domain->regions[i]->dynamic_check()) + for (auto reg : domain->get_region_list()) + if (reg->dynamic_check()) error->all(FLERR,"Cannot use PRD with a time-dependent region defined"); // perform PRD simulation diff --git a/src/RIGID/fix_ehex.cpp b/src/RIGID/fix_ehex.cpp index 01fe76757b..4dff4c6a2d 100644 --- a/src/RIGID/fix_ehex.cpp +++ b/src/RIGID/fix_ehex.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -43,18 +42,18 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{CONSTANT,EQUAL,ATOM}; +enum { CONSTANT, EQUAL, ATOM }; /* ---------------------------------------------------------------------- */ -FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - idregion(nullptr), x(nullptr), f(nullptr), v(nullptr), - mass(nullptr), rmass(nullptr), type(nullptr), scalingmask(nullptr) +FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), region(nullptr), idregion(nullptr), x(nullptr), f(nullptr), v(nullptr), + mass(nullptr), rmass(nullptr), type(nullptr), scalingmask(nullptr) { MPI_Comm_rank(world, &me); // check - if (narg < 4) error->all(FLERR,"Illegal fix ehex command: wrong number of parameters "); + if (narg < 4) error->all(FLERR, "Illegal fix ehex command: wrong number of parameters "); scalar_flag = 1; global_freq = 1; @@ -62,18 +61,16 @@ FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), // apply fix every nevery timesteps - nevery = utils::inumeric(FLERR,arg[3],false,lmp); + nevery = utils::inumeric(FLERR, arg[3], false, lmp); - if (nevery <= 0) error->all(FLERR,"Illegal fix ehex command"); + if (nevery <= 0) error->all(FLERR, "Illegal fix ehex command"); // heat flux into the reservoir - heat_input = utils::numeric(FLERR,arg[4],false,lmp); + heat_input = utils::numeric(FLERR, arg[4], false, lmp); // optional args - iregion = -1; - // NOTE: constraints are deactivated by default constraints = 0; @@ -89,12 +86,12 @@ FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), int iarg = 5; while (iarg < narg) { - if (strcmp(arg[iarg],"region") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix ehex command: wrong number of parameters "); - iregion = domain->find_region(arg[iarg+1]); - if (iregion == -1) - error->all(FLERR,"Region ID for fix ehex does not exist"); - idregion = utils::strdup(arg[iarg+1]); + if (strcmp(arg[iarg], "region") == 0) { + if (iarg + 2 > narg) + error->all(FLERR, "Illegal fix ehex command: wrong number of parameters "); + region = domain->get_region_by_id(arg[iarg + 1]); + if (!region) error->all(FLERR, "Region {} for fix ehex does not exist", arg[iarg + 1]); + idregion = utils::strdup(arg[iarg + 1]); iarg += 2; } @@ -115,10 +112,9 @@ FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), // don't apply a coordinate correction if this keyword is specified else if (strcmp(arg[iarg], "hex") == 0) { - hex = 1; - iarg+= 1; - } - else + hex = 1; + iarg += 1; + } else error->all(FLERR, "Illegal fix ehex keyword "); } @@ -128,28 +124,25 @@ FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), error->all(FLERR, "You can only use the keyword 'com' together with the keyword 'constrain' "); scale = 1.0; - scalingmask = nullptr; + scalingmask = nullptr; FixEHEX::grow_arrays(atom->nmax); atom->add_callback(Atom::GROW); - } - /* ---------------------------------------------------------------------- */ -void FixEHEX::grow_arrays(int nmax) { - memory->grow(scalingmask, nmax,"ehex:scalingmask"); +void FixEHEX::grow_arrays(int nmax) +{ + memory->grow(scalingmask, nmax, "ehex:scalingmask"); } - /* ---------------------------------------------------------------------- */ FixEHEX::~FixEHEX() { - atom->delete_callback(id,Atom::GROW); - delete [] idregion; + atom->delete_callback(id, Atom::GROW); + delete[] idregion; memory->destroy(scalingmask); - } /* ---------------------------------------------------------------------- */ @@ -167,16 +160,14 @@ void FixEHEX::init() { // set index and check validity of region - if (iregion >= 0) { - iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for fix ehex does not exist"); + if (idregion) { + region = domain->get_region_by_id(idregion); + if (!region) error->all(FLERR, "Region {} for fix ehex does not exist",idregion); } // cannot have 0 atoms in group - if (group->count(igroup) == 0) - error->all(FLERR,"Fix ehex group has no atoms"); + if (group->count(igroup) == 0) error->all(FLERR, "Fix ehex group has no atoms"); fshake = nullptr; if (constraints) { @@ -194,30 +185,29 @@ void FixEHEX::init() } if (cnt_shake > 1) - error->all(FLERR,"Multiple instances of fix shake/rattle detected (not supported yet)"); - else if (cnt_shake == 1) { - fshake = (dynamic_cast( modify->fix[id_shake])); - } - else if (cnt_shake == 0) - error->all(FLERR, "Fix ehex was configured with keyword constrain, but shake/rattle was not defined"); + error->all(FLERR, "Multiple instances of fix shake/rattle detected (not supported yet)"); + else if (cnt_shake == 1) { + fshake = (dynamic_cast(modify->fix[id_shake])); + } else if (cnt_shake == 0) + error->all( + FLERR, + "Fix ehex was configured with keyword constrain, but shake/rattle was not defined"); } } - - /* ---------------------------------------------------------------------- */ - -void FixEHEX::end_of_step() { +void FixEHEX::end_of_step() +{ // store local pointers - x = atom->x; - f = atom->f; - v = atom->v; - mass = atom->mass; - rmass = atom->rmass; - type = atom->type; - nlocal = atom->nlocal; + x = atom->x; + f = atom->f; + v = atom->v; + mass = atom->mass; + rmass = atom->rmass; + type = atom->type; + nlocal = atom->nlocal; // determine which sites are to be rescaled @@ -229,20 +219,18 @@ void FixEHEX::end_of_step() { // if required use shake/rattle to correct coordinates and velocities - if (constraints && fshake) - fshake->shake_end_of_step(0); + if (constraints && fshake) fshake->shake_end_of_step(0); } - - /* ---------------------------------------------------------------------- Iterate over all atoms, rescale the velocities and apply coordinate corrections. ------------------------------------------------------------------------- */ -void FixEHEX::rescale() { +void FixEHEX::rescale() +{ double Kr, Ke, escale; - double vsub[3],vcm[3], sfr[3]; + double vsub[3], vcm[3], sfr[3]; double mi; double dt; double F, mr, epsr_ik, sfvr, eta_ik; @@ -255,54 +243,54 @@ void FixEHEX::rescale() { // heat flux into the reservoir - F = heat_input * force->ftm2v * nevery; + F = heat_input * force->ftm2v * nevery; // total mass - mr = masstotal; + mr = masstotal; // energy scaling factor - escale = 1. + (F*dt)/Kr; + escale = 1. + (F * dt) / Kr; // safety check for kinetic energy - if (escale < 0.0) error->all(FLERR,"Fix ehex kinetic energy went negative"); + if (escale < 0.0) error->all(FLERR, "Fix ehex kinetic energy went negative"); scale = sqrt(escale); - vsub[0] = (scale-1.0) * vcm[0]; - vsub[1] = (scale-1.0) * vcm[1]; - vsub[2] = (scale-1.0) * vcm[2]; + vsub[0] = (scale - 1.0) * vcm[0]; + vsub[1] = (scale - 1.0) * vcm[1]; + vsub[2] = (scale - 1.0) * vcm[2]; for (int i = 0; i < nlocal; i++) { if (scalingmask[i]) { - mi = (rmass) ? rmass[i] : mass[type[i]]; + mi = (rmass) ? rmass[i] : mass[type[i]]; - for (int k=0; k<3; k++) { + for (int k = 0; k < 3; k++) { // apply coordinate correction unless running in hex mode if (!hex) { - // epsr_ik implements Eq. (20) in the paper + // epsr_ik implements Eq. (20) in the paper - eta_ik = mi * F/(2.*Kr) * (v[i][k] - vcm[k]); - epsr_ik = eta_ik / (mi*Kr) * (F/48. + sfvr/6.*force->ftm2v) - F/(12.*Kr) * (f[i][k]/mi - sfr[k]/mr)*force->ftm2v; + eta_ik = mi * F / (2. * Kr) * (v[i][k] - vcm[k]); + epsr_ik = eta_ik / (mi * Kr) * (F / 48. + sfvr / 6. * force->ftm2v) - + F / (12. * Kr) * (f[i][k] / mi - sfr[k] / mr) * force->ftm2v; - x[i][k] -= dt*dt*dt * epsr_ik; + x[i][k] -= dt * dt * dt * epsr_ik; } // rescale the velocity - v[i][k] = scale*v[i][k] - vsub[k]; + v[i][k] = scale * v[i][k] - vsub[k]; } } } } - /* ---------------------------------------------------------------------- */ double FixEHEX::compute_scalar() @@ -317,17 +305,17 @@ double FixEHEX::compute_scalar() double FixEHEX::memory_usage() { double bytes = 0.0; - bytes += (double)atom->nmax * sizeof(double); + bytes += (double) atom->nmax * sizeof(double); return bytes; } - /* ---------------------------------------------------------------------- Update the array scalingmask depending on which individual atoms will be rescaled or not. ------------------------------------------------------------------------- */ -void FixEHEX::update_scalingmask() { +void FixEHEX::update_scalingmask() +{ int m; int lid; bool stat; @@ -335,11 +323,7 @@ void FixEHEX::update_scalingmask() { // prematch region - Region *region = nullptr; - if (iregion >= 0) { - region = domain->regions[iregion]; - region->prematch(); - } + if (region) region->prematch(); // only rescale molecules whose center of mass if fully contained in the region @@ -347,28 +331,33 @@ void FixEHEX::update_scalingmask() { // loop over all clusters - for (int i=0; i < fshake->nlist; i++) { + for (int i = 0; i < fshake->nlist; i++) { // cluster id - m = fshake->list[i]; + m = fshake->list[i]; // check if the centre of mass of the cluster is inside the region // if region == nullptr, just check the group information of all sites - if (fshake->shake_flag[m] == 1) nsites = 3; - else if (fshake->shake_flag[m] == 2) nsites = 2; - else if (fshake->shake_flag[m] == 3) nsites = 3; - else if (fshake->shake_flag[m] == 4) nsites = 4; - else nsites = 0; + if (fshake->shake_flag[m] == 1) + nsites = 3; + else if (fshake->shake_flag[m] == 2) + nsites = 2; + else if (fshake->shake_flag[m] == 3) + nsites = 3; + else if (fshake->shake_flag[m] == 4) + nsites = 4; + else + nsites = 0; if (nsites == 0) { - error->all(FLERR,"Internal error: shake_flag[m] has to be between 1 and 4 for m in nlist"); + error->all(FLERR, "Internal error: shake_flag[m] has to be between 1 and 4 for m in nlist"); } stat = check_cluster(fshake->shake_atom[m], nsites, region); - for (int l=0; l < nsites; l++) { + for (int l = 0; l < nsites; l++) { lid = atom->map(fshake->shake_atom[m][l]); scalingmask[lid] = stat; } @@ -376,9 +365,8 @@ void FixEHEX::update_scalingmask() { // check atoms that do not belong to any cluster - for (int i=0; inlocal; i++) { - if (fshake->shake_flag[i] == 0) - scalingmask[i] = rescale_atom(i,region); + for (int i = 0; i < atom->nlocal; i++) { + if (fshake->shake_flag[i] == 0) scalingmask[i] = rescale_atom(i, region); } } @@ -386,41 +374,39 @@ void FixEHEX::update_scalingmask() { // no clusters, just individual sites (e.g. monatomic system or flexible molecules) else { - for (int i=0; inlocal; i++) - scalingmask[i] = rescale_atom(i,region); + for (int i = 0; i < atom->nlocal; i++) scalingmask[i] = rescale_atom(i, region); } - } - /* ---------------------------------------------------------------------- Check if the centre of mass of the cluster to be constrained is inside the region. ------------------------------------------------------------------------- */ -bool FixEHEX::check_cluster(tagint *shake_atom, int n, Region * region) { +bool FixEHEX::check_cluster(tagint *shake_atom, int n, Region *region) +{ // IMPORTANT NOTE: If any site of the cluster belongs to a group // which should not be rescaled than all of the sites // will be ignored! - double **x = atom->x; - double * rmass = atom->rmass; - double * mass = atom->mass; - int * type = atom->type; - int * mask = atom->mask; - double xcom[3], xtemp[3]; - double mcluster, mi; - bool stat; - int lid[4]; + double **x = atom->x; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + int *mask = atom->mask; + double xcom[3], xtemp[3]; + double mcluster, mi; + bool stat; + int lid[4]; // accumulate mass and centre of mass position - stat = true; - xcom[0] = 0.; - xcom[1] = 0.; - xcom[2] = 0.; - mcluster = 0; + stat = true; + xcom[0] = 0.; + xcom[1] = 0.; + xcom[2] = 0.; + mcluster = 0; for (int i = 0; i < n; i++) { @@ -432,26 +418,24 @@ bool FixEHEX::check_cluster(tagint *shake_atom, int n, Region * region) { stat = stat && (mask[lid[i]] & groupbit); - if (region && stat) { + if (region && stat) { // check if reduced mass is used - mi = (rmass) ? rmass[lid[i]] : mass[type[lid[i]]]; + mi = (rmass) ? rmass[lid[i]] : mass[type[lid[i]]]; mcluster += mi; // accumulate centre of mass // NOTE: you can either use unwrapped coordinates or take site x[lid[0]] as reference, // i.e. reconstruct the molecule around this site and calculate the com. - for (int k=0; k<3; k++) - xtemp[k] = x[lid[i]][k] - x[lid[0]][k]; + for (int k = 0; k < 3; k++) xtemp[k] = x[lid[i]][k] - x[lid[0]][k]; // take into account pbc domain->minimum_image(xtemp); - for (int k=0; k<3; k++) - xcom[k] += mi * (x[lid[0]][k] + xtemp[k]) ; + for (int k = 0; k < 3; k++) xcom[k] += mi * (x[lid[0]][k] + xtemp[k]); } } @@ -461,14 +445,11 @@ bool FixEHEX::check_cluster(tagint *shake_atom, int n, Region * region) { // check mass - if (mcluster < 1.e-14) { - error->all(FLERR, "Fix ehex shake cluster has almost zero mass."); - } + if (mcluster < 1.e-14) { error->all(FLERR, "Fix ehex shake cluster has almost zero mass."); } // divide by total mass - for (int k=0; k<3; k++) - xcom[k] = xcom[k]/mcluster; + for (int k = 0; k < 3; k++) xcom[k] = xcom[k] / mcluster; // apply periodic boundary conditions (centre of mass could be outside the box) // and check if molecule is inside the region @@ -480,12 +461,12 @@ bool FixEHEX::check_cluster(tagint *shake_atom, int n, Region * region) { return stat; } - /* ---------------------------------------------------------------------- Check if atom i has the correct group and is inside the region. ------------------------------------------------------------------------- */ -bool FixEHEX::rescale_atom(int i, Region*region) { +bool FixEHEX::rescale_atom(int i, Region *region) +{ bool stat; double x_r[3]; @@ -505,7 +486,7 @@ bool FixEHEX::rescale_atom(int i, Region*region) { // check if the atom is in the group/region - stat = stat && region->match(x_r[0],x_r[1],x_r[2]); + stat = stat && region->match(x_r[0], x_r[1], x_r[2]); } return stat; @@ -516,102 +497,101 @@ bool FixEHEX::rescale_atom(int i, Region*region) { (e.g. com velocity, kinetic energy, total mass,...) ------------------------------------------------------------------------- */ -void FixEHEX::com_properties(double * vr, double * sfr, double *sfvr, double *K, double *Kr, double *mr) { - double ** f = atom->f; - double ** v = atom->v; - int nlocal = atom->nlocal; - double *rmass= atom->rmass; - double *mass = atom->mass; - int *type = atom->type; - double l_vr[3]; - double l_mr; - double l_sfr[3]; - double l_sfvr; - double l_K; - double mi; - double l_buf[9]; - double buf[9]; +void FixEHEX::com_properties(double *vr, double *sfr, double *sfvr, double *K, double *Kr, + double *mr) +{ + double **f = atom->f; + double **v = atom->v; + int nlocal = atom->nlocal; + double *rmass = atom->rmass; + double *mass = atom->mass; + int *type = atom->type; + double l_vr[3]; + double l_mr; + double l_sfr[3]; + double l_sfvr; + double l_K; + double mi; + double l_buf[9]; + double buf[9]; - // calculate partial sums + // calculate partial sums - l_vr[0] = l_vr[1] = l_vr[2] = 0; - l_sfr[0] = l_sfr[1] = l_sfr[2] = 0; - l_sfvr = 0; - l_mr = 0; - l_K = 0; + l_vr[0] = l_vr[1] = l_vr[2] = 0; + l_sfr[0] = l_sfr[1] = l_sfr[2] = 0; + l_sfvr = 0; + l_mr = 0; + l_K = 0; - for (int i = 0; i < nlocal; i++) { - if (scalingmask[i]) { + for (int i = 0; i < nlocal; i++) { + if (scalingmask[i]) { - // check if reduced mass is used + // check if reduced mass is used - mi = (rmass) ? rmass[i] : mass[type[i]]; + mi = (rmass) ? rmass[i] : mass[type[i]]; - // accumulate total mass + // accumulate total mass - l_mr += mi; + l_mr += mi; - // accumulate kinetic energy + // accumulate kinetic energy - l_K += mi/2. * (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]); + l_K += mi / 2. * (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]); - // sum_j f_j * v_j + // sum_j f_j * v_j - l_sfvr += f[i][0]*v[i][0] + f[i][1]*v[i][1] + f[i][2]*v[i][2]; + l_sfvr += f[i][0] * v[i][0] + f[i][1] * v[i][1] + f[i][2] * v[i][2]; - // accumulate com velocity and sum of forces + // accumulate com velocity and sum of forces - for (int k=0; k<3; k++) { - l_vr[k] += mi * v[i][k]; - l_sfr[k]+= f[i][k]; - } - } - } + for (int k = 0; k < 3; k++) { + l_vr[k] += mi * v[i][k]; + l_sfr[k] += f[i][k]; + } + } + } - // reduce sums + // reduce sums - l_buf[0] = l_vr[0]; - l_buf[1] = l_vr[1]; - l_buf[2] = l_vr[2]; - l_buf[3] = l_K; - l_buf[4] = l_mr; - l_buf[5] = l_sfr[0]; - l_buf[6] = l_sfr[1]; - l_buf[7] = l_sfr[2]; - l_buf[8] = l_sfvr; + l_buf[0] = l_vr[0]; + l_buf[1] = l_vr[1]; + l_buf[2] = l_vr[2]; + l_buf[3] = l_K; + l_buf[4] = l_mr; + l_buf[5] = l_sfr[0]; + l_buf[6] = l_sfr[1]; + l_buf[7] = l_sfr[2]; + l_buf[8] = l_sfvr; - MPI_Allreduce(l_buf, buf, 9, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(l_buf, buf, 9, MPI_DOUBLE, MPI_SUM, world); - // total mass of region + // total mass of region - *mr = buf[4]; + *mr = buf[4]; - if (*mr < 1.e-14) { - error->all(FLERR, "Fix ehex error mass of region is close to zero"); - } + if (*mr < 1.e-14) { error->all(FLERR, "Fix ehex error mass of region is close to zero"); } - // total kinetic energy of region + // total kinetic energy of region - *K = buf[3]; + *K = buf[3]; - // centre of mass velocity of region + // centre of mass velocity of region - vr[0] = buf[0]/(*mr); - vr[1] = buf[1]/(*mr); - vr[2] = buf[2]/(*mr); + vr[0] = buf[0] / (*mr); + vr[1] = buf[1] / (*mr); + vr[2] = buf[2] / (*mr); - // sum of forces + // sum of forces - sfr[0] = buf[5]; - sfr[1] = buf[6]; - sfr[2] = buf[7]; + sfr[0] = buf[5]; + sfr[1] = buf[6]; + sfr[2] = buf[7]; - // calculate non-translational kinetic energy + // calculate non-translational kinetic energy - *Kr = *K - 0.5* (*mr) * (vr[0]*vr[0]+vr[1]*vr[1]+vr[2]*vr[2]); + *Kr = *K - 0.5 * (*mr) * (vr[0] * vr[0] + vr[1] * vr[1] + vr[2] * vr[2]); - // calculate sum_j f_j * (v_j - v_r) = sum_j f_j * v_j - v_r * sum_j f_j + // calculate sum_j f_j * (v_j - v_r) = sum_j f_j * v_j - v_r * sum_j f_j - *sfvr = buf[8] - (vr[0]*sfr[0] + vr[1]*sfr[1] + vr[2]*sfr[2]); + *sfvr = buf[8] - (vr[0] * sfr[0] + vr[1] * sfr[1] + vr[2] * sfr[2]); } - diff --git a/src/RIGID/fix_ehex.h b/src/RIGID/fix_ehex.h index f4ea872533..7b3b9a5ff8 100644 --- a/src/RIGID/fix_ehex.h +++ b/src/RIGID/fix_ehex.h @@ -43,10 +43,10 @@ class FixEHEX : public Fix { bool check_cluster(tagint *shake_atom, int n, class Region *region); private: - int iregion; double heat_input; double masstotal; double scale; + class Region *region; char *idregion; int me; diff --git a/src/VTK/dump_vtk.cpp b/src/VTK/dump_vtk.cpp index fd7f4b2c2b..51c6632b78 100644 --- a/src/VTK/dump_vtk.cpp +++ b/src/VTK/dump_vtk.cpp @@ -246,12 +246,11 @@ void DumpVTK::init_style() else if (flag && cols) custom_flag[i] = DARRAY; } - // set index and check validity of region + // check validity of region - if (iregion >= 0) { - iregion = domain->find_region(idregion); - if (iregion == -1) - error->all(FLERR,"Region ID for dump vtk does not exist"); + if (idregion) { + if (!domain->get_region_by_id(idregion)) + error->all(FLERR,"Region {} for dump vtk does not exist",idregion); } } @@ -335,8 +334,8 @@ int DumpVTK::count() // un-choose if not in region - if (iregion >= 0) { - Region *region = domain->regions[iregion]; + auto region = domain->get_region_by_id(idregion); + if (region) { region->prematch(); double **x = atom->x; for (i = 0; i < nlocal; i++) @@ -2054,11 +2053,12 @@ int DumpVTK::modify_param(int narg, char **arg) { if (strcmp(arg[0],"region") == 0) { if (narg < 2) error->all(FLERR,"Illegal dump_modify command"); - if (strcmp(arg[1],"none") == 0) iregion = -1; - else { - iregion = domain->find_region(arg[1]); - if (iregion == -1) - error->all(FLERR,"Dump_modify region ID {} does not exist",arg[1]); + if (strcmp(arg[1],"none") == 0) { + delete[] idregion; + idregion = nullptr; + } else { + if (!domain->get_region_by_id(arg[1])) + error->all(FLERR,"Dump_modify region {} does not exist",arg[1]); delete[] idregion; idregion = utils::strdup(arg[1]); } diff --git a/src/domain.cpp b/src/domain.cpp index d3804d75bd..3abe351e36 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -106,9 +106,6 @@ Domain::Domain(LAMMPS *lmp) : Pointers(lmp) set_lattice(2,args); delete [] args; - nregion = maxregion = 0; - regions = nullptr; - copymode = 0; region_map = new RegionCreatorMap(); @@ -128,10 +125,8 @@ Domain::~Domain() { if (copymode) return; + region_list.clear(); delete lattice; - for (int i = 0; i < nregion; i++) delete regions[i]; - memory->sfree(regions); - delete region_map; } @@ -194,7 +189,7 @@ void Domain::init() // region inits - for (int i = 0; i < nregion; i++) regions[i]->init(); + for (auto reg : region_list) reg->init(); } /* ---------------------------------------------------------------------- @@ -1758,88 +1753,67 @@ void Domain::add_region(int narg, char **arg) if (strcmp(arg[1],"none") == 0) error->all(FLERR,"Unrecognized region style 'none'"); - if (find_region(arg[0]) >= 0) error->all(FLERR,"Reuse of region ID"); - - // extend Region list if necessary - - if (nregion == maxregion) { - maxregion += DELTAREGION; - regions = (Region **) - memory->srealloc(regions,maxregion*sizeof(Region *),"domain:regions"); - } + if (get_region_by_id(arg[0])) error->all(FLERR,"Reuse of region ID {}", arg[0]); // create the Region + Region *newregion = nullptr; if (lmp->suffix_enable) { if (lmp->suffix) { std::string estyle = std::string(arg[1]) + "/" + lmp->suffix; if (region_map->find(estyle) != region_map->end()) { RegionCreator ®ion_creator = (*region_map)[estyle]; - regions[nregion] = region_creator(lmp, narg, arg); - regions[nregion]->init(); - nregion++; - return; + newregion = region_creator(lmp, narg, arg); } } - if (lmp->suffix2) { + if (!newregion && lmp->suffix2) { std::string estyle = std::string(arg[1]) + "/" + lmp->suffix2; if (region_map->find(estyle) != region_map->end()) { RegionCreator ®ion_creator = (*region_map)[estyle]; - regions[nregion] = region_creator(lmp, narg, arg); - regions[nregion]->init(); - nregion++; - return; + newregion = region_creator(lmp, narg, arg); } } } - if (region_map->find(arg[1]) != region_map->end()) { + if (!newregion && (region_map->find(arg[1]) != region_map->end())) { RegionCreator ®ion_creator = (*region_map)[arg[1]]; - regions[nregion] = region_creator(lmp, narg, arg); - } else error->all(FLERR,utils::check_packages_for_style("region",arg[1],lmp)); + newregion = region_creator(lmp, narg, arg); + } + + if (!newregion) + error->all(FLERR,utils::check_packages_for_style("region",arg[1],lmp)); // initialize any region variables via init() // in case region is used between runs, e.g. to print a variable - regions[nregion]->init(); - nregion++; + newregion->init(); + region_list.push_back(newregion); } /* ---------------------------------------------------------------------- delete a region ------------------------------------------------------------------------- */ -void Domain::delete_region(int iregion) +void Domain::delete_region(Region *reg) { - if ((iregion < 0) || (iregion >= nregion)) return; + if (!reg) return; // delete and move other Regions down in list one slot - - delete regions[iregion]; - for (int i = iregion+1; i < nregion; ++i) - regions[i-1] = regions[i]; - nregion--; + bool match = false; + for (std::size_t i = 0; i < region_list.size(); ++i) { + if (match) region_list[i-1] = region_list[i]; + if (reg == region_list[i]) match = true; + } + delete reg; + region_list.resize(region_list.size() - 1); } void Domain::delete_region(const std::string &id) { - int iregion = find_region(id); - if (iregion == -1) error->all(FLERR,"Delete region ID does not exist"); - - delete_region(iregion); -} - -/* ---------------------------------------------------------------------- - return region index if name matches existing region ID - return -1 if no such region -------------------------------------------------------------------------- */ - -int Domain::find_region(const std::string &name) const -{ - for (int iregion = 0; iregion < nregion; iregion++) - if (name == regions[iregion]->id) return iregion; - return -1; + auto reg = get_region_by_id(id); + if (!reg) error->all(FLERR,"Delete region {} does not exist", id); + delete_region(reg); } /* ---------------------------------------------------------------------- @@ -1849,8 +1823,8 @@ int Domain::find_region(const std::string &name) const Region *Domain::get_region_by_id(const std::string &name) const { - for (int iregion = 0; iregion < nregion; iregion++) - if (name == regions[iregion]->id) return regions[iregion]; + for (auto ® : region_list) + if (name == reg->id) return reg; return nullptr; } @@ -1864,12 +1838,21 @@ const std::vector Domain::get_region_by_style(const std::string &name) std::vector matches; if (name.empty()) return matches; - for (int iregion = 0; iregion < nregion; iregion++) - if (name == regions[iregion]->style) matches.push_back(regions[iregion]); + for (auto ® : region_list) + if (name == reg->style) matches.push_back(reg); return matches; } +/* ---------------------------------------------------------------------- + return list of fixes as vector +------------------------------------------------------------------------- */ + +const std::vector &Domain::get_region_list() +{ + return region_list; +} + /* ---------------------------------------------------------------------- (re)set boundary settings flag = 0, called from the input script diff --git a/src/domain.h b/src/domain.h index 68ff5aece5..07ca199236 100644 --- a/src/domain.h +++ b/src/domain.h @@ -17,6 +17,7 @@ #include "pointers.h" #include +#include #include namespace LAMMPS_NS { @@ -98,10 +99,6 @@ class Domain : protected Pointers { class Lattice *lattice; // user-defined lattice - int nregion; // # of defined Regions - int maxregion; // max # list can hold - Region **regions; // list of defined Regions - int copymode; enum { NO_REMAP, X_REMAP, V_REMAP }; @@ -137,11 +134,11 @@ class Domain : protected Pointers { void set_lattice(int, char **); void add_region(int, char **); - void delete_region(int); + void delete_region(Region *); void delete_region(const std::string &); - int find_region(const std::string &) const; Region *get_region_by_id(const std::string &) const; const std::vector get_region_by_style(const std::string &) const; + const std::vector &get_region_list(); void set_boundary(int, char **, int); void set_box(int, char **); void print_box(const std::string &); @@ -175,6 +172,7 @@ class Domain : protected Pointers { protected: double small[3]; // fractions of box lengths + std::vector region_list; }; } // namespace LAMMPS_NS diff --git a/src/group.cpp b/src/group.cpp index d678f2de86..41c24b8f48 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -40,8 +40,8 @@ using namespace LAMMPS_NS; -#define MAX_GROUP 32 -#define EPSILON 1.0e-6 +static constexpr int MAX_GROUP = 32; +static constexpr double EPSILON = 1.0e-6; enum{NONE,TYPE,MOLECULE,ID}; enum{LT,LE,GT,GE,EQ,NEQ,BETWEEN}; @@ -176,13 +176,13 @@ void Group::assign(int narg, char **arg) if (narg != 3) error->all(FLERR,"Illegal group command"); - int iregion = domain->find_region(arg[2]); - if (iregion == -1) error->all(FLERR,"Group region ID does not exist"); - domain->regions[iregion]->init(); - domain->regions[iregion]->prematch(); + auto region = domain->get_region_by_id(arg[2]); + if (!region) error->all(FLERR,"Group region {} does not exist", arg[2]); + region->init(); + region->prematch(); for (i = 0; i < nlocal; i++) - if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) + if (region->match(x[i][0],x[i][1],x[i][2])) mask[i] |= bit; // create an empty group @@ -813,15 +813,6 @@ bigint Group::count(int igroup, Region *region) return nall; } -/* ---------------------------------------------------------------------- - count atoms in group and region -------------------------------------------------------------------------- */ - -bigint Group::count(int igroup, int iregion) -{ - return count(igroup, domain->regions[iregion]); -} - /* ---------------------------------------------------------------------- compute the total mass of group of atoms use either per-type mass or per-atom rmass @@ -886,16 +877,6 @@ double Group::mass(int igroup, Region *region) return all; } -/* ---------------------------------------------------------------------- - compute the total mass of group of atoms in region - use either per-type mass or per-atom rmass -------------------------------------------------------------------------- */ - -double Group::mass(int igroup, int iregion) -{ - return mass(igroup, domain->regions[iregion]); -} - /* ---------------------------------------------------------------------- compute the total charge of group of atoms ------------------------------------------------------------------------- */ @@ -921,10 +902,9 @@ double Group::charge(int igroup) compute the total charge of group of atoms in region ------------------------------------------------------------------------- */ -double Group::charge(int igroup, int iregion) +double Group::charge(int igroup, Region *region) { int groupbit = bitmask[igroup]; - Region *region = domain->regions[iregion]; region->prematch(); double **x = atom->x; @@ -990,10 +970,9 @@ void Group::bounds(int igroup, double *minmax) periodic images are not considered, so atoms are NOT unwrapped ------------------------------------------------------------------------- */ -void Group::bounds(int igroup, double *minmax, int iregion) +void Group::bounds(int igroup, double *minmax, Region *region) { int groupbit = bitmask[igroup]; - Region *region = domain->regions[iregion]; region->prematch(); double extent[6]; @@ -1090,10 +1069,9 @@ void Group::xcm(int igroup, double masstotal, double *cm) must unwrap atoms to compute center-of-mass correctly ------------------------------------------------------------------------- */ -void Group::xcm(int igroup, double masstotal, double *cm, int iregion) +void Group::xcm(int igroup, double masstotal, double *cm, Region *region) { int groupbit = bitmask[igroup]; - Region *region = domain->regions[iregion]; region->prematch(); double **x = atom->x; @@ -1232,17 +1210,6 @@ void Group::vcm(int igroup, double masstotal, double *cm, Region *region) } } -/* ---------------------------------------------------------------------- - compute the center-of-mass velocity of group of atoms in region - masstotal = total mass - return center-of-mass velocity in cm[] -------------------------------------------------------------------------- */ - -void Group::vcm(int igroup, double masstotal, double *cm, int iregion) -{ - vcm(igroup, masstotal, cm, domain->regions[iregion]); -} - /* ---------------------------------------------------------------------- compute the total force on group of atoms ------------------------------------------------------------------------- */ @@ -1272,10 +1239,9 @@ void Group::fcm(int igroup, double *cm) compute the total force on group of atoms in region ------------------------------------------------------------------------- */ -void Group::fcm(int igroup, double *cm, int iregion) +void Group::fcm(int igroup, double *cm, Region *region) { int groupbit = bitmask[igroup]; - Region *region = domain->regions[iregion]; region->prematch(); double **x = atom->x; @@ -1368,15 +1334,6 @@ double Group::ke(int igroup, Region *region) return all; } -/* ---------------------------------------------------------------------- - compute the total kinetic energy of group of atoms in region and return it -------------------------------------------------------------------------- */ - -double Group::ke(int igroup, int iregion) -{ - return ke(igroup, domain->regions[iregion]); -} - /* ---------------------------------------------------------------------- compute the radius-of-gyration of group of atoms around center-of-mass cm @@ -1422,10 +1379,9 @@ double Group::gyration(int igroup, double masstotal, double *cm) must unwrap atoms to compute Rg correctly ------------------------------------------------------------------------- */ -double Group::gyration(int igroup, double masstotal, double *cm, int iregion) +double Group::gyration(int igroup, double masstotal, double *cm, Region *region) { int groupbit = bitmask[igroup]; - Region *region = domain->regions[iregion]; region->prematch(); double **x = atom->x; @@ -1504,10 +1460,9 @@ void Group::angmom(int igroup, double *cm, double *lmom) must unwrap atoms to compute L correctly ------------------------------------------------------------------------- */ -void Group::angmom(int igroup, double *cm, double *lmom, int iregion) +void Group::angmom(int igroup, double *cm, double *lmom, Region *region) { int groupbit = bitmask[igroup]; - Region *region = domain->regions[iregion]; region->prematch(); double **x = atom->x; @@ -1583,10 +1538,9 @@ void Group::torque(int igroup, double *cm, double *tq) must unwrap atoms to compute T correctly ------------------------------------------------------------------------- */ -void Group::torque(int igroup, double *cm, double *tq, int iregion) +void Group::torque(int igroup, double *cm, double *tq, Region *region) { int groupbit = bitmask[igroup]; - Region *region = domain->regions[iregion]; region->prematch(); double **x = atom->x; @@ -1669,12 +1623,11 @@ void Group::inertia(int igroup, double *cm, double itensor[3][3]) must unwrap atoms to compute itensor correctly ------------------------------------------------------------------------- */ -void Group::inertia(int igroup, double *cm, double itensor[3][3], int iregion) +void Group::inertia(int igroup, double *cm, double itensor[3][3], Region *region) { int i,j; int groupbit = bitmask[igroup]; - Region *region = domain->regions[iregion]; region->prematch(); double **x = atom->x; diff --git a/src/group.h b/src/group.h index 6b075bb4f8..66c5e737f9 100644 --- a/src/group.h +++ b/src/group.h @@ -41,33 +41,29 @@ class Group : protected Pointers { bigint count_all(); // count atoms in group all bigint count(int); // count atoms in group - bigint count(int, int); // count atoms in group & region bigint count(int, Region *); // count atoms in group & region double mass(int); // total mass of atoms in group - double mass(int, int); double mass(int, Region *); double charge(int); // total charge of atoms in group - double charge(int, int); + double charge(int, Region *); void bounds(int, double *); // bounds of atoms in group - void bounds(int, double *, int); + void bounds(int, double *, Region *); void xcm(int, double, double *); // center-of-mass coords of group - void xcm(int, double, double *, int); + void xcm(int, double, double *, Region *); void vcm(int, double, double *); // center-of-mass velocity of group - void vcm(int, double, double *, int); void vcm(int, double, double *, Region *); void fcm(int, double *); // total force on group - void fcm(int, double *, int); + void fcm(int, double *, Region *); double ke(int); // kinetic energy of group - double ke(int, int); double ke(int, Region *); double gyration(int, double, double *); // radius-of-gyration of group - double gyration(int, double, double *, int); + double gyration(int, double, double *, Region *); void angmom(int, double *, double *); // angular momentum of group - void angmom(int, double *, double *, int); + void angmom(int, double *, double *, Region *); void torque(int, double *, double *); // torque on group - void torque(int, double *, double *, int); + void torque(int, double *, double *, Region *); void inertia(int, double *, double[3][3]); // inertia tensor - void inertia(int, double *, double[3][3], int); + void inertia(int, double *, double[3][3], Region *); void omega(double *, double[3][3], double *); // angular velocity private: diff --git a/src/info.cpp b/src/info.cpp index 5e53864253..d50c45457b 100644 --- a/src/info.cpp +++ b/src/info.cpp @@ -546,19 +546,17 @@ void Info::command(int narg, char **arg) } if (flags & REGIONS) { - int nreg = domain->nregion; - Region **regs = domain->regions; fputs("\nRegion information:\n",out); - for (int i=0; i < nreg; ++i) { + int i=0; + for (auto ® : domain->get_region_list()) { fmt::print(out,"Region[{:3d}]: {:16} style = {:16} side = {}\n", - i, std::string(regs[i]->id)+',', - std::string(regs[i]->style)+',', - regs[i]->interior ? "in" : "out"); - if (regs[i]->bboxflag) + i, std::string(reg->id)+',', std::string(reg->style)+',', + reg->interior ? "in" : "out"); + if (reg->bboxflag) fmt::print(out," Boundary: lo {:.8} {:.8} {:.8} hi {:.8} {:.8} {:.8}\n", - regs[i]->extent_xlo, regs[i]->extent_ylo, - regs[i]->extent_zlo, regs[i]->extent_xhi, - regs[i]->extent_yhi, regs[i]->extent_zhi); + reg->extent_xlo, reg->extent_ylo, + reg->extent_zlo, reg->extent_xhi, + reg->extent_yhi, reg->extent_zhi); else fputs(" No Boundary\n",out); } } @@ -916,12 +914,8 @@ bool Info::is_defined(const char *category, const char *name) return true; } } else if (strcmp(category,"region") == 0) { - int nreg = domain->nregion; - Region **regs = domain->regions; - for (int i=0; i < nreg; ++i) { - if (strcmp(regs[i]->id,name) == 0) - return true; - } + for (auto ® : domain->get_region_list()) + if (strcmp(reg->id,name) == 0) return true; } else if (strcmp(category,"variable") == 0) { int nvar = input->variable->nvar; char **names = input->variable->names; diff --git a/src/library.cpp b/src/library.cpp index 8ab97321ba..793c817e9a 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -4772,10 +4772,8 @@ int lammps_has_id(void *handle, const char *category, const char *name) { if (strcmp(name,molecule[i]->id) == 0) return 1; } } else if (strcmp(category,"region") == 0) { - int nregion = lmp->domain->nregion; - Region **region = lmp->domain->regions; - for (int i=0; i < nregion; ++i) { - if (strcmp(name,region[i]->id) == 0) return 1; + for (auto ® : lmp->domain->get_region_list()) { + if (strcmp(name,reg->id) == 0) return 1; } } else if (strcmp(category,"variable") == 0) { int nvariable = lmp->input->variable->nvar; @@ -4818,7 +4816,7 @@ int lammps_id_count(void *handle, const char *category) { } else if (strcmp(category,"molecule") == 0) { return lmp->atom->nmolecule; } else if (strcmp(category,"region") == 0) { - return lmp->domain->nregion; + return lmp->domain->get_region_list().size(); } else if (strcmp(category,"variable") == 0) { return lmp->input->variable->nvar; } @@ -4848,8 +4846,7 @@ set to an empty string, otherwise 1. * \param buf_size size of the provided string buffer * \return 1 if successful, otherwise 0 */ -int lammps_id_name(void *handle, const char *category, int idx, - char *buffer, int buf_size) { +int lammps_id_name(void *handle, const char *category, int idx, char *buffer, int buf_size) { auto lmp = (LAMMPS *) handle; if (strcmp(category,"compute") == 0) { @@ -4878,8 +4875,9 @@ int lammps_id_name(void *handle, const char *category, int idx, return 1; } } else if (strcmp(category,"region") == 0) { - if ((idx >=0) && (idx < lmp->domain->nregion)) { - strncpy(buffer, lmp->domain->regions[idx]->id, buf_size); + auto regions = lmp->domain->get_region_list(); + if ((idx >=0) && (idx < (int) regions.size())) { + strncpy(buffer, regions[idx]->id, buf_size); return 1; } } else if (strcmp(category,"variable") == 0) { diff --git a/src/math_const.h b/src/math_const.h index 1ad0b3e717..1a6b5c3b47 100644 --- a/src/math_const.h +++ b/src/math_const.h @@ -23,6 +23,7 @@ namespace MathConst { static constexpr double MY_2PI = 6.28318530717958647692; // 2pi static constexpr double MY_3PI = 9.42477796076937971538; // 3pi static constexpr double MY_4PI = 12.56637061435917295384; // 4pi + static constexpr double MY_4PI3 = 4.18879020478639098461; // 4/3pi static constexpr double MY_PI2 = 1.57079632679489661923; // pi/2 static constexpr double MY_PI4 = 0.78539816339744830962; // pi/4 static constexpr double MY_PIS = 1.77245385090551602729; // sqrt(pi) diff --git a/src/modify.cpp b/src/modify.cpp index 7554079e2a..fa4c8fd01b 100644 --- a/src/modify.cpp +++ b/src/modify.cpp @@ -859,7 +859,7 @@ Fix *Modify::add_fix(int narg, char **arg, int trysuffix) fix[ifix]->style = utils::strdup(estyle); } } - if (fix[ifix] == nullptr && lmp->suffix2) { + if ((fix[ifix] == nullptr) && lmp->suffix2) { std::string estyle = arg[2] + std::string("/") + lmp->suffix2; if (fix_map->find(estyle) != fix_map->end()) { FixCreator &fix_creator = (*fix_map)[estyle]; @@ -870,7 +870,7 @@ Fix *Modify::add_fix(int narg, char **arg, int trysuffix) } } - if (fix[ifix] == nullptr && fix_map->find(arg[2]) != fix_map->end()) { + if ((fix[ifix] == nullptr) && (fix_map->find(arg[2]) != fix_map->end())) { FixCreator &fix_creator = (*fix_map)[arg[2]]; fix[ifix] = fix_creator(lmp, narg, arg); } diff --git a/src/set.cpp b/src/set.cpp index 1abd6c758c..d644abe780 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -700,13 +700,13 @@ void Set::selection(int n) else select[i] = 0; } else if (style == REGION_SELECT) { - int iregion = domain->find_region(id); - if (iregion == -1) error->all(FLERR,"Set region ID does not exist"); - domain->regions[iregion]->prematch(); + auto region = domain->get_region_by_id(id); + if (!region) error->all(FLERR,"Set region {} does not exist", id); + region->prematch(); double **x = atom->x; for (int i = 0; i < n; i++) - if (domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2])) + if (region->match(x[i][0],x[i][1],x[i][2])) select[i] = 1; else select[i] = 0; } diff --git a/src/variable.cpp b/src/variable.cpp index 42afda5cd8..8455227ad7 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -3043,18 +3043,18 @@ double Variable::eval_tree(Tree *tree, int i) } if (tree->type == GMASK) { - if (atom->mask[i] & tree->ivalue1) return 1.0; + if (atom->mask[i] & tree->ivalue) return 1.0; else return 0.0; } if (tree->type == RMASK) { - if (domain->regions[tree->ivalue1]->match(atom->x[i][0], atom->x[i][1], atom->x[i][2])) return 1.0; + if (tree->region->match(atom->x[i][0], atom->x[i][1], atom->x[i][2])) return 1.0; else return 0.0; } if (tree->type == GRMASK) { - if ((atom->mask[i] & tree->ivalue1) && - (domain->regions[tree->ivalue2]->match(atom->x[i][0], atom->x[i][1], atom->x[i][2]))) return 1.0; + if ((atom->mask[i] & tree->ivalue) && + (tree->region->match(atom->x[i][0], atom->x[i][1], atom->x[i][2]))) return 1.0; else return 0.0; } @@ -3664,32 +3664,31 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre int igroup = group->find(args[0]); if (igroup == -1) { - std::string mesg = "Group ID '"; - mesg += args[0]; - mesg += "' in variable formula does not exist"; - print_var_error(FLERR,mesg,ivar); + const auto errmesg = fmt::format("Group {} in variable formula does not exist", args[0]); + print_var_error(FLERR, errmesg, ivar); } // match word to group function double value = 0.0; + const auto group_errmesg = fmt::format("Invalid {}() function in variable formula", word); if (strcmp(word,"count") == 0) { if (narg == 1) value = group->count(igroup); else if (narg == 2) value = group->count(igroup,region_function(args[1],ivar)); - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"mass") == 0) { if (narg == 1) value = group->mass(igroup); else if (narg == 2) value = group->mass(igroup,region_function(args[1],ivar)); - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"charge") == 0) { if (narg == 1) value = group->charge(igroup); else if (narg == 2) value = group->charge(igroup,region_function(args[1],ivar)); - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"xcm") == 0) { atom->check_mass(FLERR); @@ -3698,14 +3697,14 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre double masstotal = group->mass(igroup); group->xcm(igroup,masstotal,xcm); } else if (narg == 3) { - int iregion = region_function(args[2],ivar); - double masstotal = group->mass(igroup,iregion); - group->xcm(igroup,masstotal,xcm,iregion); - } else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + auto region = region_function(args[2],ivar); + double masstotal = group->mass(igroup,region); + group->xcm(igroup,masstotal,xcm,region); + } else print_var_error(FLERR,group_errmesg,ivar); if (strcmp(args[1],"x") == 0) value = xcm[0]; else if (strcmp(args[1],"y") == 0) value = xcm[1]; else if (strcmp(args[1],"z") == 0) value = xcm[2]; - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"vcm") == 0) { atom->check_mass(FLERR); @@ -3714,38 +3713,38 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre double masstotal = group->mass(igroup); group->vcm(igroup,masstotal,vcm); } else if (narg == 3) { - int iregion = region_function(args[2],ivar); - double masstotal = group->mass(igroup,iregion); - group->vcm(igroup,masstotal,vcm,iregion); - } else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + auto region = region_function(args[2],ivar); + double masstotal = group->mass(igroup,region); + group->vcm(igroup,masstotal,vcm,region); + } else print_var_error(FLERR,group_errmesg,ivar); if (strcmp(args[1],"x") == 0) value = vcm[0]; else if (strcmp(args[1],"y") == 0) value = vcm[1]; else if (strcmp(args[1],"z") == 0) value = vcm[2]; - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"fcm") == 0) { double fcm[3]; if (narg == 2) group->fcm(igroup,fcm); else if (narg == 3) group->fcm(igroup,fcm,region_function(args[2],ivar)); - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); if (strcmp(args[1],"x") == 0) value = fcm[0]; else if (strcmp(args[1],"y") == 0) value = fcm[1]; else if (strcmp(args[1],"z") == 0) value = fcm[2]; - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"bound") == 0) { double minmax[6]; if (narg == 2) group->bounds(igroup,minmax); else if (narg == 3) group->bounds(igroup,minmax,region_function(args[2],ivar)); - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); if (strcmp(args[1],"xmin") == 0) value = minmax[0]; else if (strcmp(args[1],"xmax") == 0) value = minmax[1]; else if (strcmp(args[1],"ymin") == 0) value = minmax[2]; else if (strcmp(args[1],"ymax") == 0) value = minmax[3]; else if (strcmp(args[1],"zmin") == 0) value = minmax[4]; else if (strcmp(args[1],"zmax") == 0) value = minmax[5]; - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"gyration") == 0) { atom->check_mass(FLERR); @@ -3755,16 +3754,16 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre group->xcm(igroup,masstotal,xcm); value = group->gyration(igroup,masstotal,xcm); } else if (narg == 2) { - int iregion = region_function(args[1],ivar); - double masstotal = group->mass(igroup,iregion); - group->xcm(igroup,masstotal,xcm,iregion); - value = group->gyration(igroup,masstotal,xcm,iregion); - } else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + auto region = region_function(args[1],ivar); + double masstotal = group->mass(igroup,region); + group->xcm(igroup,masstotal,xcm,region); + value = group->gyration(igroup,masstotal,xcm,region); + } else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"ke") == 0) { if (narg == 1) value = group->ke(igroup); else if (narg == 2) value = group->ke(igroup,region_function(args[1],ivar)); - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"angmom") == 0) { atom->check_mass(FLERR); @@ -3774,15 +3773,15 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre group->xcm(igroup,masstotal,xcm); group->angmom(igroup,xcm,lmom); } else if (narg == 3) { - int iregion = region_function(args[2],ivar); - double masstotal = group->mass(igroup,iregion); - group->xcm(igroup,masstotal,xcm,iregion); - group->angmom(igroup,xcm,lmom,iregion); - } else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + auto region = region_function(args[2],ivar); + double masstotal = group->mass(igroup,region); + group->xcm(igroup,masstotal,xcm,region); + group->angmom(igroup,xcm,lmom,region); + } else print_var_error(FLERR,group_errmesg,ivar); if (strcmp(args[1],"x") == 0) value = lmom[0]; else if (strcmp(args[1],"y") == 0) value = lmom[1]; else if (strcmp(args[1],"z") == 0) value = lmom[2]; - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"torque") == 0) { atom->check_mass(FLERR); @@ -3792,15 +3791,15 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre group->xcm(igroup,masstotal,xcm); group->torque(igroup,xcm,tq); } else if (narg == 3) { - int iregion = region_function(args[2],ivar); - double masstotal = group->mass(igroup,iregion); - group->xcm(igroup,masstotal,xcm,iregion); - group->torque(igroup,xcm,tq,iregion); - } else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + auto region = region_function(args[2],ivar); + double masstotal = group->mass(igroup,region); + group->xcm(igroup,masstotal,xcm,region); + group->torque(igroup,xcm,tq,region); + } else print_var_error(FLERR,group_errmesg,ivar); if (strcmp(args[1],"x") == 0) value = tq[0]; else if (strcmp(args[1],"y") == 0) value = tq[1]; else if (strcmp(args[1],"z") == 0) value = tq[2]; - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"inertia") == 0) { atom->check_mass(FLERR); @@ -3810,18 +3809,18 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre group->xcm(igroup,masstotal,xcm); group->inertia(igroup,xcm,inertia); } else if (narg == 3) { - int iregion = region_function(args[2],ivar); - double masstotal = group->mass(igroup,iregion); - group->xcm(igroup,masstotal,xcm,iregion); - group->inertia(igroup,xcm,inertia,iregion); - } else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + auto region = region_function(args[2],ivar); + double masstotal = group->mass(igroup,region); + group->xcm(igroup,masstotal,xcm,region); + group->inertia(igroup,xcm,inertia,region); + } else print_var_error(FLERR,group_errmesg,ivar); if (strcmp(args[1],"xx") == 0) value = inertia[0][0]; else if (strcmp(args[1],"yy") == 0) value = inertia[1][1]; else if (strcmp(args[1],"zz") == 0) value = inertia[2][2]; else if (strcmp(args[1],"xy") == 0) value = inertia[0][1]; else if (strcmp(args[1],"yz") == 0) value = inertia[1][2]; else if (strcmp(args[1],"xz") == 0) value = inertia[0][2]; - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } else if (strcmp(word,"omega") == 0) { atom->check_mass(FLERR); @@ -3833,17 +3832,17 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre group->inertia(igroup,xcm,inertia); group->omega(angmom,inertia,omega); } else if (narg == 3) { - int iregion = region_function(args[2],ivar); - double masstotal = group->mass(igroup,iregion); - group->xcm(igroup,masstotal,xcm,iregion); - group->angmom(igroup,xcm,angmom,iregion); - group->inertia(igroup,xcm,inertia,iregion); + auto region = region_function(args[2],ivar); + double masstotal = group->mass(igroup,region); + group->xcm(igroup,masstotal,xcm,region); + group->angmom(igroup,xcm,angmom,region); + group->inertia(igroup,xcm,inertia,region); group->omega(angmom,inertia,omega); - } else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + } else print_var_error(FLERR,group_errmesg,ivar); if (strcmp(args[1],"x") == 0) value = omega[0]; else if (strcmp(args[1],"y") == 0) value = omega[1]; else if (strcmp(args[1],"z") == 0) value = omega[2]; - else print_var_error(FLERR,"Invalid group function in variable formula",ivar); + else print_var_error(FLERR,group_errmesg,ivar); } // delete stored args @@ -3864,21 +3863,16 @@ int Variable::group_function(char *word, char *contents, Tree **tree, Tree **tre /* ---------------------------------------------------------------------- */ -int Variable::region_function(char *id, int ivar) +Region *Variable::region_function(char *id, int ivar) { - int iregion = domain->find_region(id); - if (iregion == -1) { - std::string mesg = "Region ID '"; - mesg += id; - mesg += "' in variable formula does not exist"; - print_var_error(FLERR,mesg,ivar); - } + auto region = domain->get_region_by_id(id); + if (!region) + print_var_error(FLERR, fmt::format("Region {} in variable formula does not exist", id), ivar); // init region in case sub-regions have been deleted - domain->regions[iregion]->init(); - - return iregion; + region->init(); + return region; } /* ---------------------------------------------------------------------- @@ -4149,7 +4143,7 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t auto newtree = new Tree(); newtree->type = GMASK; - newtree->ivalue1 = group->bitmask[igroup]; + newtree->ivalue = group->bitmask[igroup]; treestack[ntreestack++] = newtree; } else if (strcmp(word,"rmask") == 0) { @@ -4158,12 +4152,12 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t if (narg != 1) print_var_error(FLERR,"Invalid special function in variable formula",ivar); - int iregion = region_function(args[0],ivar); - domain->regions[iregion]->prematch(); + auto region = region_function(args[0],ivar); + region->prematch(); auto newtree = new Tree(); newtree->type = RMASK; - newtree->ivalue1 = iregion; + newtree->region = region; treestack[ntreestack++] = newtree; } else if (strcmp(word,"grmask") == 0) { @@ -4175,13 +4169,13 @@ int Variable::special_function(char *word, char *contents, Tree **tree, Tree **t int igroup = group->find(args[0]); if (igroup == -1) print_var_error(FLERR,"Group ID in variable formula does not exist",ivar); - int iregion = region_function(args[1],ivar); - domain->regions[iregion]->prematch(); + auto region = region_function(args[1],ivar); + region->prematch(); auto newtree = new Tree(); newtree->type = GRMASK; - newtree->ivalue1 = group->bitmask[igroup]; - newtree->ivalue2 = iregion; + newtree->ivalue = group->bitmask[igroup]; + newtree->region = region; treestack[ntreestack++] = newtree; // special function for file-style or atomfile-style variables diff --git a/src/variable.h b/src/variable.h index 7b51b45f38..909c3f8af9 100644 --- a/src/variable.h +++ b/src/variable.h @@ -17,6 +17,7 @@ #include "pointers.h" namespace LAMMPS_NS { +class Region; class Variable : protected Pointers { friend class Info; @@ -110,14 +111,15 @@ class Variable : protected Pointers { int nvector; // length of array for vector-style variable int nstride; // stride between atoms if array is a 2d array int selfalloc; // 1 if array is allocated here, else 0 - int ivalue1, ivalue2; // extra values needed for gmask,rmask,grmask + int ivalue; // extra value needed for gmask, grmask int nextra; // # of additional args beyond first 2 + Region *region; // region pointer for rmask, grmask Tree *first, *second; // ptrs further down tree for first 2 args Tree **extra; // ptrs further down tree for nextra args Tree() : - array(nullptr), iarray(nullptr), barray(nullptr), selfalloc(0), ivalue1(0), ivalue2(0), - nextra(0), first(nullptr), second(nullptr), extra(nullptr) + array(nullptr), iarray(nullptr), barray(nullptr), selfalloc(0), ivalue(0), + nextra(0), region(nullptr), first(nullptr), second(nullptr), extra(nullptr) { } }; @@ -135,7 +137,7 @@ class Variable : protected Pointers { int find_matching_paren(char *, int, char *&, int); int math_function(char *, char *, Tree **, Tree **, int &, double *, int &, int); int group_function(char *, char *, Tree **, Tree **, int &, double *, int &, int); - int region_function(char *, int); + Region *region_function(char *, int); int special_function(char *, char *, Tree **, Tree **, int &, double *, int &, int); void peratom2global(int, char *, double *, int, tagint, Tree **, Tree **, int &, double *, int &); int is_atom_vector(char *); diff --git a/unittest/commands/test_groups.cpp b/unittest/commands/test_groups.cpp index a356a02cca..b0706ea775 100644 --- a/unittest/commands/test_groups.cpp +++ b/unittest/commands/test_groups.cpp @@ -151,7 +151,8 @@ TEST_F(GroupTest, RegionClear) ASSERT_EQ(group->count_all(), lmp->atom->natoms); TEST_FAILURE(".*ERROR: Illegal group command.*", command("group three region left xxx");); - TEST_FAILURE(".*ERROR: Group region ID does not exist.*", command("group four region dummy");); + TEST_FAILURE(".*ERROR: Group region dummy does not exist.*", + command("group four region dummy");); BEGIN_HIDE_OUTPUT(); command("group one clear"); @@ -196,8 +197,8 @@ TEST_F(GroupTest, SelectRestart) ASSERT_EQ(group->count(group->find("four")), 32); ASSERT_EQ(group->count(group->find("five")), 16); ASSERT_EQ(group->count(group->find("six")), 8); - ASSERT_EQ(group->count(group->find("half"), domain->find_region("top")), 8); - ASSERT_DOUBLE_EQ(group->mass(group->find("half"), domain->find_region("top")), 8.0); + ASSERT_EQ(group->count(group->find("half"), domain->get_region_by_id("top")), 8); + ASSERT_DOUBLE_EQ(group->mass(group->find("half"), domain->get_region_by_id("top")), 8.0); BEGIN_HIDE_OUTPUT(); command("write_restart group.restart"); @@ -243,9 +244,9 @@ TEST_F(GroupTest, Molecular) ASSERT_EQ(group->count(group->find("two")), 16); ASSERT_EQ(group->count(group->find("three")), 15); ASSERT_DOUBLE_EQ(group->mass(group->find("half")), 40); - ASSERT_DOUBLE_EQ(group->mass(group->find("half"), domain->find_region("top")), 10); + ASSERT_DOUBLE_EQ(group->mass(group->find("half"), domain->get_region_by_id("top")), 10); ASSERT_NEAR(group->charge(group->find("top")), 0, 1.0e-14); - ASSERT_NEAR(group->charge(group->find("right"), domain->find_region("top")), 0, 1.0e-14); + ASSERT_NEAR(group->charge(group->find("right"), domain->get_region_by_id("top")), 0, 1.0e-14); TEST_FAILURE(".*ERROR: Illegal group command.*", command("group three include xxx");); }