From 02d189fb23d7b61f1cacb3bbf38bc5e157ad8739 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 13 Jul 2023 10:53:56 -0400 Subject: [PATCH 1/9] avoid false positive with static code analysis --- src/compute_count_type.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/compute_count_type.cpp b/src/compute_count_type.cpp index 3cb2e6bc3a..3d4815f9ff 100644 --- a/src/compute_count_type.cpp +++ b/src/compute_count_type.cpp @@ -159,6 +159,8 @@ void ComputeCountType::compute_vector() nvec = count_dihedrals(); else if (mode == IMPROPER) nvec = count_impropers(); + else + nvec = 0; // sum across procs as bigint, then convert to double // correct for multiple counting if newton_bond off From 78d3d4948f281c62e57e4f5ec200afcca5145ef2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 14 Jul 2023 20:06:31 -0400 Subject: [PATCH 2/9] fix off-by-one bug in argument parsing --- src/REPLICA/temper_npt.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/REPLICA/temper_npt.cpp b/src/REPLICA/temper_npt.cpp index 6567cb83e1..6aa895e813 100644 --- a/src/REPLICA/temper_npt.cpp +++ b/src/REPLICA/temper_npt.cpp @@ -89,7 +89,9 @@ void TemperNPT::command(int narg, char **arg) seed_boltz = utils::inumeric(FLERR,arg[5],false,lmp); my_set_temp = universe->iworld; - if (narg == 8) my_set_temp = utils::inumeric(FLERR,arg[6],false,lmp); + if (narg == 8) my_set_temp = utils::inumeric(FLERR,arg[7],false,lmp); + if ((my_set_temp < 0) || (my_set_temp > 7)) + error->universe_all(FLERR,"Invalid partition number for temperature index keyword"); // swap frequency must evenly divide total # of timesteps From e6873bb7c816cea6e5f46ad47f8b650b4b70ec79 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 14 Jul 2023 20:07:09 -0400 Subject: [PATCH 3/9] modernize code --- src/REPLICA/temper_npt.cpp | 35 ++++++++++++++++------------------- src/REPLICA/temper_npt.h | 3 +-- 2 files changed, 17 insertions(+), 21 deletions(-) diff --git a/src/REPLICA/temper_npt.cpp b/src/REPLICA/temper_npt.cpp index 6aa895e813..303eff59b3 100644 --- a/src/REPLICA/temper_npt.cpp +++ b/src/REPLICA/temper_npt.cpp @@ -52,10 +52,10 @@ TemperNPT::~TemperNPT() MPI_Comm_free(&roots); if (ranswap) delete ranswap; delete ranboltz; - delete [] set_temp; - delete [] temp2world; - delete [] world2temp; - delete [] world2root; + delete[] set_temp; + delete[] temp2world; + delete[] world2temp; + delete[] world2root; } /* ---------------------------------------------------------------------- @@ -68,8 +68,7 @@ void TemperNPT::command(int narg, char **arg) error->all(FLERR,"Must have more than one processor partition to temper"); if (domain->box_exist == 0) error->all(FLERR,"temper/npt command before simulation box is defined"); - if (narg != 7 && narg != 8) - error->universe_all(FLERR,"Illegal temper/npt command"); + if (narg != 7 && narg != 8) error->universe_all(FLERR,"Illegal temper/npt command"); int nsteps = utils::inumeric(FLERR,arg[0],false,lmp); nevery = utils::inumeric(FLERR,arg[1],false,lmp); @@ -80,10 +79,9 @@ void TemperNPT::command(int narg, char **arg) if (timer->is_timeout()) return; - for (whichfix = 0; whichfix < modify->nfix; whichfix++) - if (strcmp(arg[3],modify->fix[whichfix]->id) == 0) break; - if (whichfix == modify->nfix) - error->universe_all(FLERR,"Tempering fix ID is not defined"); + whichfix = modify->get_fix_by_id(arg[3]); + if (!whichfix) + error->universe_all(FLERR,fmt::format("Tempering fix ID {} is not defined", arg[3])); seed_swap = utils::inumeric(FLERR,arg[4],false,lmp); seed_boltz = utils::inumeric(FLERR,arg[5],false,lmp); @@ -106,8 +104,8 @@ void TemperNPT::command(int narg, char **arg) // change the volume. This currently only applies to fix npt and // fix rigid/npt variants - if ( (!utils::strmatch(modify->fix[whichfix]->style,"^npt")) && - (!utils::strmatch(modify->fix[whichfix]->style,"^rigid/npt")) ) + if ( (!utils::strmatch(whichfix->style,"^npt")) && + (!utils::strmatch(whichfix->style,"^rigid/npt")) ) error->universe_all(FLERR,"Tempering temperature and pressure fix is not supported"); // setup for long tempering run @@ -118,8 +116,7 @@ void TemperNPT::command(int narg, char **arg) update->nsteps = nsteps; update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + nsteps; - if (update->laststep < 0) - error->all(FLERR,"Too many timesteps"); + if (update->laststep < 0) error->all(FLERR,"Too many timesteps"); lmp->init(); @@ -135,9 +132,9 @@ void TemperNPT::command(int narg, char **arg) // pe_compute = ptr to thermo_pe compute // notify compute it will be called at first swap - int id = modify->find_compute("thermo_pe"); - if (id < 0) error->all(FLERR,"Tempering could not find thermo_pe compute"); - Compute *pe_compute = modify->compute[id]; + Compute *pe_compute = modify->get_compute_by_id("thermo_pe"); + if (!pe_compute) error->all(FLERR,"Tempering could not find thermo_pe compute"); + pe_compute->addstep(update->ntimestep + nevery); // create MPI communicator for root proc from each world @@ -186,7 +183,7 @@ void TemperNPT::command(int narg, char **arg) if (narg == 8) { double new_temp = set_temp[my_set_temp]; - modify->fix[whichfix]->reset_target(new_temp); + whichfix->reset_target(new_temp); } // setup tempering runs @@ -327,7 +324,7 @@ void TemperNPT::command(int narg, char **arg) if (swap) { new_temp = set_temp[partner_set_temp]; - modify->fix[whichfix]->reset_target(new_temp); + whichfix->reset_target(new_temp); } // update my_set_temp and temp2world on every proc diff --git a/src/REPLICA/temper_npt.h b/src/REPLICA/temper_npt.h index 3a51fc74e4..2e831c55ef 100644 --- a/src/REPLICA/temper_npt.h +++ b/src/REPLICA/temper_npt.h @@ -42,8 +42,7 @@ class TemperNPT : public Command { int nswaps; // # of tempering swaps to perform int seed_swap; // 0 = toggle swaps, n = RNG for swap direction int seed_boltz; // seed for Boltz factor comparison - int whichfix; // index of temperature fix to use - int fixstyle; // what kind of temperature fix is used + class Fix *whichfix; // index of temperature fix to use int my_set_temp; // which set temp I am simulating double *set_temp; // static list of replica set temperatures From baac049aed55d644783289d424ed0528c81e0ad3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 15 Jul 2023 12:25:03 -0400 Subject: [PATCH 4/9] update LAMMPS input file syntax highlighting for recent changes --- tools/vim/lammps.vim | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/tools/vim/lammps.vim b/tools/vim/lammps.vim index 008ae1ee3f..9565d0a86e 100644 --- a/tools/vim/lammps.vim +++ b/tools/vim/lammps.vim @@ -1,23 +1,27 @@ " Vim syntax file -" Language: Lammps Simulation Script File +" Language: Lammps Simulation Script File " Maintainer: Gerolf Ziegenhain " Updates: Axel Kohlmeyer , Sam Bateman , Daniel Möller Montull , Eryk Skalinski -" Latest Revision: 2022-08-17 +" Latest Revision: 2023-07-15 syn clear -syn keyword lammpsOutput log write_data write_dump info shell write_restart restart dump undump thermo thermo_modify thermo_style print timer +" Add '/' to list of valid keyword characters +set iskeyword+=/ + +syn keyword lammpsOutput log write_data write_dump write_coeff info shell write_restart restart dump undump thermo thermo_modify +syn keyword lammpsOutput thermo_style print timer syn keyword lammpsRead include read_restart read_data read_dump molecule syn keyword lammpsLattice boundary units atom_style lattice region create_box create_atoms dielectric syn keyword lammpsLattice delete_atoms displace_atoms change_box dimension replicate -syn keyword lammpsParticle pair_coeff pair_style pair_modify pair_write mass velocity angle_coeff angle_style +syn keyword lammpsParticle pair_coeff pair_style pair_modify pair_write mass velocity angle_coeff angle_style angle_write syn keyword lammpsParticle atom_modify atom_style bond_coeff bond_style bond_write create_bonds delete_bonds kspace_style -syn keyword lammpsParticle kspace_modify dihedral_style dihedral_coeff improper_style improper_coeff -syn keyword lammpsSetup min_style fix_modify run_style timestep neighbor neigh_modify fix unfix suffix special_bonds -syn keyword lammpsSetup balance box clear comm_modify comm_style newton package processors reset_ids reset_timestep -syn keyword lammpsRun minimize run rerun tad neb prd quit server temper temper/grem temper/npt -syn keyword lammpsRun min/spin message hyper dynamical_matrix -syn keyword lammpsDefine variable group compute python set uncompute kim_query +syn keyword lammpsParticle kspace_modify dihedral_style dihedral_coeff dihedral_write improper_style improper_coeff labelmap +syn keyword lammpsSetup min_style min_modify fix_modify run_style timestep neighbor neigh_modify fix unfix suffix special_bonds +syn keyword lammpsSetup balance box clear comm_modify comm_style newton package processors reset_atoms reset_ids reset_timestep +syn keyword lammpsRun minimize minimize/kk run rerun tad neb neb/spin prd quit server temper/npt temper/grem temper +syn keyword lammpsRun message hyper dynamical_matrix dynamical_matrix/kk third_order third_order/kk fitpod +syn keyword lammpsDefine variable group compute python set uncompute kim_query kim group2ndx ndx2group mdi syn keyword lammpsRepeat jump next loop From 77bdcb3e19679ae43ed1b4b969e5b818939be5f1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 15 Jul 2023 12:36:41 -0400 Subject: [PATCH 5/9] small doc style update --- doc/src/dynamical_matrix.rst | 2 +- doc/src/minimize.rst | 4 ++-- doc/src/third_order.rst | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/dynamical_matrix.rst b/doc/src/dynamical_matrix.rst index 93e4c2a1aa..0243fb49db 100644 --- a/doc/src/dynamical_matrix.rst +++ b/doc/src/dynamical_matrix.rst @@ -4,7 +4,7 @@ dynamical_matrix command ======================== -Accelerator Variants: dynamical_matrix/kk +Accelerator Variant: dynamical_matrix/kk Syntax """""" diff --git a/doc/src/minimize.rst b/doc/src/minimize.rst index c71c8e3632..56efa12f44 100644 --- a/doc/src/minimize.rst +++ b/doc/src/minimize.rst @@ -1,10 +1,10 @@ .. index:: minimize +.. index:: minimize/kk minimize command ================ -minimize/kk command -=================== +Accelerator Variant: minimize/kk Syntax """""" diff --git a/doc/src/third_order.rst b/doc/src/third_order.rst index 2635c92f84..bbbf424c44 100644 --- a/doc/src/third_order.rst +++ b/doc/src/third_order.rst @@ -4,7 +4,7 @@ third_order command =================== -Accelerator Variants: third_order/kk +Accelerator Variant: third_order/kk Syntax """""" From fd0a72eab53889af860c33d406b1390fd59ce609 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 15 Jul 2023 12:37:12 -0400 Subject: [PATCH 6/9] modernize and consolidate style across all three temper command variants --- src/REPLICA/temper.cpp | 58 +++++++++++++--------------- src/REPLICA/temper.h | 3 +- src/REPLICA/temper_grem.cpp | 75 ++++++++++++++++++------------------- src/REPLICA/temper_grem.h | 4 +- src/REPLICA/temper_npt.cpp | 33 ++++++++-------- 5 files changed, 82 insertions(+), 91 deletions(-) diff --git a/src/REPLICA/temper.cpp b/src/REPLICA/temper.cpp index 9e7e67d79d..adbdb4d742 100644 --- a/src/REPLICA/temper.cpp +++ b/src/REPLICA/temper.cpp @@ -37,7 +37,7 @@ using namespace LAMMPS_NS; -// #define TEMPER_DEBUG 1 +#define TEMPER_DEBUG 0 /* ---------------------------------------------------------------------- */ @@ -50,10 +50,10 @@ Temper::~Temper() MPI_Comm_free(&roots); if (ranswap) delete ranswap; delete ranboltz; - delete [] set_temp; - delete [] temp2world; - delete [] world2temp; - delete [] world2root; + delete[] set_temp; + delete[] temp2world; + delete[] world2temp; + delete[] world2root; } /* ---------------------------------------------------------------------- @@ -63,11 +63,10 @@ Temper::~Temper() void Temper::command(int narg, char **arg) { if (universe->nworlds == 1) - error->all(FLERR,"Must have more than one processor partition to temper"); + error->universe_all(FLERR,"More than one processor partition required for temper command"); if (domain->box_exist == 0) - error->all(FLERR,"Temper command before simulation box is defined"); - if (narg != 6 && narg != 7) - error->universe_all(FLERR,"Illegal temper command"); + error->universe_all(FLERR,"Temper command before simulation box is defined"); + if (narg != 6 && narg != 7) error->universe_all(FLERR,"Illegal temper command"); int nsteps = utils::inumeric(FLERR,arg[0],false,lmp); nevery = utils::inumeric(FLERR,arg[1],false,lmp); @@ -77,10 +76,9 @@ void Temper::command(int narg, char **arg) if (timer->is_timeout()) return; - for (whichfix = 0; whichfix < modify->nfix; whichfix++) - if (strcmp(arg[3],modify->fix[whichfix]->id) == 0) break; - if (whichfix == modify->nfix) - error->universe_all(FLERR,"Tempering fix ID is not defined"); + whichfix = modify->get_fix_by_id(arg[3]); + if (!whichfix) + error->universe_all(FLERR,fmt::format("Tempering fix ID {} is not defined", arg[3])); seed_swap = utils::inumeric(FLERR,arg[4],false,lmp); seed_boltz = utils::inumeric(FLERR,arg[5],false,lmp); @@ -88,7 +86,7 @@ void Temper::command(int narg, char **arg) my_set_temp = universe->iworld; if (narg == 7) my_set_temp = utils::inumeric(FLERR,arg[6],false,lmp); if ((my_set_temp < 0) || (my_set_temp >= universe->nworlds)) - error->universe_one(FLERR,"Illegal temperature index"); + error->universe_one(FLERR,"Invalid temperature index value"); // swap frequency must evenly divide total # of timesteps @@ -101,11 +99,11 @@ void Temper::command(int narg, char **arg) // fix style must be appropriate for temperature control, i.e. it needs // to provide a working Fix::reset_target() and must not change the volume. - if ((!utils::strmatch(modify->fix[whichfix]->style,"^nvt")) && - (!utils::strmatch(modify->fix[whichfix]->style,"^langevin")) && - (!utils::strmatch(modify->fix[whichfix]->style,"^gl[de]$")) && - (!utils::strmatch(modify->fix[whichfix]->style,"^rigid/nvt")) && - (!utils::strmatch(modify->fix[whichfix]->style,"^temp/"))) + if ((!utils::strmatch(whichfix->style,"^nvt")) && + (!utils::strmatch(whichfix->style,"^langevin")) && + (!utils::strmatch(whichfix->style,"^gl[de]$")) && + (!utils::strmatch(whichfix->style,"^rigid/nvt")) && + (!utils::strmatch(whichfix->style,"^temp/"))) error->universe_all(FLERR,"Tempering temperature fix is not supported"); // setup for long tempering run @@ -116,8 +114,7 @@ void Temper::command(int narg, char **arg) update->nsteps = nsteps; update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + nsteps; - if (update->laststep < 0) - error->all(FLERR,"Too many timesteps"); + if (update->laststep < 0) error->all(FLERR,"Too many timesteps"); lmp->init(); @@ -132,9 +129,9 @@ void Temper::command(int narg, char **arg) // pe_compute = ptr to thermo_pe compute // notify compute it will be called at first swap - int id = modify->find_compute("thermo_pe"); - if (id < 0) error->all(FLERR,"Tempering could not find thermo_pe compute"); - Compute *pe_compute = modify->compute[id]; + Compute *pe_compute = modify->get_compute_by_id("thermo_pe"); + if (!pe_compute) error->all(FLERR,"Tempering could not find thermo_pe compute"); + pe_compute->addstep(update->ntimestep + nevery); // create MPI communicator for root proc from each world @@ -183,7 +180,7 @@ void Temper::command(int narg, char **arg) if (narg == 7) { double new_temp = set_temp[my_set_temp]; - modify->fix[whichfix]->reset_target(new_temp); + whichfix->reset_target(new_temp); } // setup tempering runs @@ -289,13 +286,12 @@ void Temper::command(int narg, char **arg) else MPI_Recv(&swap,1,MPI_INT,partner,0,universe->uworld,MPI_STATUS_IGNORE); -#ifdef TEMPER_DEBUG +#if TEMPER_DEBUG if (me_universe < partner) - printf("SWAP %d & %d: yes = %d,Ts = %d %d, PEs = %g %g, Bz = %g %g\n", - me_universe,partner,swap,my_set_temp,partner_set_temp, - pe,pe_partner,boltz_factor,exp(boltz_factor)); + fprintf(universe->uscreen,"SWAP %d & %d: yes = %d,Ts = %d %d, PEs = %g %g, Bz = %g %g\n", + me_universe,partner,swap,my_set_temp,partner_set_temp, + pe,pe_partner,boltz_factor,exp(boltz_factor)); #endif - } // bcast swap result to other procs in my world @@ -310,7 +306,7 @@ void Temper::command(int narg, char **arg) if (swap) { new_temp = set_temp[partner_set_temp]; - modify->fix[whichfix]->reset_target(new_temp); + whichfix->reset_target(new_temp); } // update my_set_temp and temp2world on every proc diff --git a/src/REPLICA/temper.h b/src/REPLICA/temper.h index b557e8ddd8..d959892587 100644 --- a/src/REPLICA/temper.h +++ b/src/REPLICA/temper.h @@ -40,8 +40,7 @@ class Temper : public Command { int nswaps; // # of tempering swaps to perform int seed_swap; // 0 = toggle swaps, n = RNG for swap direction int seed_boltz; // seed for Boltz factor comparison - int whichfix; // index of temperature fix to use - int fixstyle; // what kind of temperature fix is used + class Fix *whichfix; // index of temperature fix to use int my_set_temp; // which set temp I am simulating double *set_temp; // static list of replica set temperatures diff --git a/src/REPLICA/temper_grem.cpp b/src/REPLICA/temper_grem.cpp index 68fbb4eef6..303f502309 100644 --- a/src/REPLICA/temper_grem.cpp +++ b/src/REPLICA/temper_grem.cpp @@ -24,6 +24,7 @@ #include "finish.h" #include "fix.h" #include "fix_grem.h" +#include "fix_nh.h" #include "force.h" #include "integrate.h" #include "modify.h" @@ -37,7 +38,7 @@ using namespace LAMMPS_NS; -//#define TEMPER_DEBUG 1 +#define TEMPER_DEBUG 0 /* ---------------------------------------------------------------------- */ @@ -50,11 +51,10 @@ TemperGrem::~TemperGrem() MPI_Comm_free(&roots); if (ranswap) delete ranswap; delete ranboltz; - delete [] set_lambda; - delete [] lambda2world; - delete [] world2lambda; - delete [] world2root; - delete [] id_nh; + delete[] set_lambda; + delete[] lambda2world; + delete[] world2lambda; + delete[] world2root; } /* ---------------------------------------------------------------------- @@ -64,46 +64,50 @@ TemperGrem::~TemperGrem() void TemperGrem::command(int narg, char **arg) { if (universe->nworlds == 1) - error->all(FLERR,"Must have more than one processor partition to temper"); + error->universe_all(FLERR,"More than one processor partition required for temper/grem command"); if (domain->box_exist == 0) - error->all(FLERR,"Temper/gREM command before simulation box is defined"); - if (narg != 7 && narg != 8) - error->universe_all(FLERR,"Illegal temper command"); + error->universe_all(FLERR,"Temper/grem command before simulation box is defined"); + if (narg != 7 && narg != 8) error->universe_all(FLERR,"Illegal temper/grem command"); int nsteps = utils::inumeric(FLERR,arg[0],false,lmp); nevery = utils::inumeric(FLERR,arg[1],false,lmp); double lambda = utils::numeric(FLERR,arg[2],false,lmp); // ignore temper command, if walltime limit was already reached + if (timer->is_timeout()) return; - // Get and check if gREM fix exists - for (whichfix = 0; whichfix < modify->nfix; whichfix++) - if (strcmp(arg[3],modify->fix[whichfix]->id) == 0) break; - if (whichfix == modify->nfix) - error->universe_all(FLERR,"Tempering fix ID is not defined"); - fix_grem = dynamic_cast(modify->fix[whichfix]); + // Get and check if gREM fix exists and is correct style + + auto ifix = modify->get_fix_by_id(arg[3]); + if (!ifix) error->universe_all(FLERR,fmt::format("Tempering fix ID {} is not defined", arg[3])); + + fix_grem = dynamic_cast(ifix); + if (!fix_grem || (strcmp(ifix->style,"grem") != 0)) + error->universe_all(FLERR,"Tempering temperature fix is of incorrect style"); // Check input values lambdas should be equal, assign other gREM values if (lambda != fix_grem->lambda) - error->universe_all(FLERR,"Lambda from tempering and fix in the same world" - " must be the same"); + error->universe_all(FLERR,"Lambda from tempering and fix in the same world must be the same"); double eta = fix_grem->eta; double h0 = fix_grem->h0; double pressref = 0; // Get and check for nh fix - id_nh = utils::strdup(arg[4]); - int ifix = modify->find_fix(id_nh); - if (ifix < 0) - error->all(FLERR,"Fix id for nvt or npt fix does not exist"); - Fix *nh = modify->fix[ifix]; + + FixNH *nh = dynamic_cast(modify->get_fix_by_id(arg[4])); + if (!nh) + error->all(FLERR,fmt::format("Fix {} for Nose-Hoover fix does not exist", arg[4])); // get result from nvt vs npt check from fix_grem + int pressflag = fix_grem->pressflag; + // fix_grem does all the checking... + if (pressflag) { - auto p_start = (double *) nh->extract("p_start",ifix); + int dummy; + auto p_start = (double *) nh->extract("p_start",dummy); pressref = p_start[0]; } @@ -123,11 +127,6 @@ void TemperGrem::command(int narg, char **arg) if (nswaps*nevery != nsteps) error->universe_all(FLERR,"Non integer # of swaps in temper command"); - // Must be used with fix_grem - - if (strcmp(modify->fix[whichfix]->style,"grem") != 0) - error->universe_all(FLERR,"Tempering temperature fix is not supported"); - // setup for long tempering run update->whichflag = 1; @@ -136,8 +135,7 @@ void TemperGrem::command(int narg, char **arg) update->nsteps = nsteps; update->beginstep = update->firststep = update->ntimestep; update->endstep = update->laststep = update->firststep + nsteps; - if (update->laststep < 0) - error->all(FLERR,"Too many timesteps"); + if (update->laststep < 0) error->all(FLERR,"Too many timesteps"); lmp->init(); @@ -152,9 +150,9 @@ void TemperGrem::command(int narg, char **arg) // pe_compute = ptr to thermo_pe compute // notify compute it will be called at first swap - int id = modify->find_compute("thermo_pe"); - if (id < 0) error->all(FLERR,"Tempering could not find thermo_pe compute"); - Compute *pe_compute = modify->compute[id]; + Compute *pe_compute = modify->get_compute_by_id("thermo_pe"); + if (!pe_compute) error->all(FLERR,"Tempering could not find thermo_pe compute"); + pe_compute->addstep(update->ntimestep + nevery); // create MPI communicator for root proc from each world @@ -319,13 +317,12 @@ void TemperGrem::command(int narg, char **arg) else MPI_Recv(&swap,1,MPI_INT,partner,0,universe->uworld,MPI_STATUS_IGNORE); -#ifdef TEMPER_DEBUG +#if TEMPER_DEBUG if (me_universe < partner) - printf("SWAP %d & %d: yes = %d,Ts = %d %d, PEs = %g %g, Bz = %g %g\n", - me_universe,partner,swap,my_set_lambda,partner_set_lambda, - weight,weight_partner,boltz_factor,exp(boltz_factor)); + fprintf(universe->uscreen,"SWAP %d & %d: yes = %d,Ts = %d %d, PEs = %g %g, Bz = %g %g\n", + me_universe,partner,swap,my_set_lambda,partner_set_lambda, + weight,weight_partner,boltz_factor,exp(boltz_factor)); #endif - } // bcast swap result to other procs in my world diff --git a/src/REPLICA/temper_grem.h b/src/REPLICA/temper_grem.h index 3c41f9b082..c66239a794 100644 --- a/src/REPLICA/temper_grem.h +++ b/src/REPLICA/temper_grem.h @@ -40,8 +40,7 @@ class TemperGrem : public Command { int nswaps; // # of tempering swaps to perform int seed_swap; // 0 = toggle swaps, n = RNG for swap direction int seed_boltz; // seed for Boltz factor comparison - int whichfix; // index of temperature fix to use - int fixstyle; // what kind of temperature fix is used + class Fix *whichfix; // index of temperature fix to use int my_set_lambda; // which set lambda I am simulating double *set_lambda; // static list of replica set lambdas @@ -54,7 +53,6 @@ class TemperGrem : public Command { class FixGrem *fix_grem; protected: - char *id_nh; int pressflag; }; diff --git a/src/REPLICA/temper_npt.cpp b/src/REPLICA/temper_npt.cpp index 303eff59b3..d814bf6725 100644 --- a/src/REPLICA/temper_npt.cpp +++ b/src/REPLICA/temper_npt.cpp @@ -65,9 +65,9 @@ TemperNPT::~TemperNPT() void TemperNPT::command(int narg, char **arg) { if (universe->nworlds == 1) - error->all(FLERR,"Must have more than one processor partition to temper"); + error->universe_all(FLERR,"More than one processor partition required for temper/npt command"); if (domain->box_exist == 0) - error->all(FLERR,"temper/npt command before simulation box is defined"); + error->universe_all(FLERR,"Temper/npt command before simulation box is defined"); if (narg != 7 && narg != 8) error->universe_all(FLERR,"Illegal temper/npt command"); int nsteps = utils::inumeric(FLERR,arg[0],false,lmp); @@ -88,8 +88,8 @@ void TemperNPT::command(int narg, char **arg) my_set_temp = universe->iworld; if (narg == 8) my_set_temp = utils::inumeric(FLERR,arg[7],false,lmp); - if ((my_set_temp < 0) || (my_set_temp > 7)) - error->universe_all(FLERR,"Invalid partition number for temperature index keyword"); + if ((my_set_temp < 0) || (my_set_temp >= universe->nworlds)) + error->universe_one(FLERR,"Invalid temperature index value"); // swap frequency must evenly divide total # of timesteps @@ -278,22 +278,23 @@ void TemperNPT::command(int narg, char **arg) swap = 0; if (partner != -1) { - if (me_universe > partner) { + if (me_universe > partner) MPI_Send(&pe,1,MPI_DOUBLE,partner,0,universe->uworld); - } - else { + else MPI_Recv(&pe_partner,1,MPI_DOUBLE,partner,0,universe->uworld,MPI_STATUS_IGNORE); - } - if (me_universe > partner) { + + if (me_universe > partner) MPI_Send(&vol,1, MPI_DOUBLE,partner,0,universe->uworld); - } - else { + else MPI_Recv(&vol_partner,1,MPI_DOUBLE,partner,0,universe->uworld,MPI_STATUS_IGNORE); - } - // Acceptance criteria changed for NPT ensemble + + // Acceptance criteria changed versus temper command for NPT ensemble if (me_universe < partner) { press_units = press_set/nktv2p; - delr = (pe_partner - pe)*(1.0/(boltz*set_temp[my_set_temp]) - 1.0/(boltz*set_temp[partner_set_temp])) + press_units*(1.0/(boltz*set_temp[my_set_temp]) - 1.0/(boltz*set_temp[partner_set_temp]))*(vol_partner - vol); + delr = (pe_partner - pe)*(1.0/(boltz*set_temp[my_set_temp]) + - 1.0/(boltz*set_temp[partner_set_temp])) + + press_units*(1.0/(boltz*set_temp[my_set_temp]) + - 1.0/(boltz*set_temp[partner_set_temp]))*(vol_partner - vol); boltz_factor = -delr; if (boltz_factor >= 0.0) swap = 1; else if (ranboltz->uniform() < exp(boltz_factor)) swap = 1; @@ -303,13 +304,13 @@ void TemperNPT::command(int narg, char **arg) MPI_Send(&swap,1,MPI_INT,partner,0,universe->uworld); else MPI_Recv(&swap,1,MPI_INT,partner,0,universe->uworld,MPI_STATUS_IGNORE); -#ifdef TEMPER_DEBUG + +#if TEMPER_DEBUG if (me_universe < partner) fprintf(universe->uscreen,"SWAP %d & %d: yes = %d,Ts = %d %d, PEs = %g %g, Bz = %g %g, vol = %g %g\n", me_universe,partner,swap,my_set_temp,partner_set_temp, pe,pe_partner,boltz_factor,exp(boltz_factor), vol, vol_partner); #endif - } // bcast swap result to other procs in my world From 89d82fde22eb625bcf42cdba8d8b5d489ffddba2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 16 Jul 2023 15:20:36 -0400 Subject: [PATCH 7/9] modernize access to list of fixes --- src/fix_wall_reflect.cpp | 9 ++++----- src/reset_atoms_id.cpp | 6 +++--- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/fix_wall_reflect.cpp b/src/fix_wall_reflect.cpp index 60a055fbd6..00ef968828 100644 --- a/src/fix_wall_reflect.cpp +++ b/src/fix_wall_reflect.cpp @@ -176,12 +176,11 @@ void FixWallReflect::init() } int nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) nrigid++; + for (const auto &ifix : modify->get_fix_list()) + if (ifix->rigid_flag) nrigid++; - if (nrigid && comm->me == 0) - error->warning(FLERR,"Should not allow rigid bodies to bounce off " - "relecting walls"); + if (nrigid && (comm->me == 0)) + error->warning(FLERR,"Should not use reflecting walls with rigid bodies"); } /* ---------------------------------------------------------------------- */ diff --git a/src/reset_atoms_id.cpp b/src/reset_atoms_id.cpp index c0fce95326..9e8ba3630d 100644 --- a/src/reset_atoms_id.cpp +++ b/src/reset_atoms_id.cpp @@ -55,9 +55,9 @@ void ResetAtomsID::command(int narg, char **arg) error->all(FLERR, "Reset_atoms id command before simulation box is defined"); if (atom->tag_enable == 0) error->all(FLERR, "Cannot use reset_atoms id unless atoms have IDs"); - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->stores_ids) - error->all(FLERR, "Cannot use reset_atoms id when a fix exists that stores atom IDs"); + for (const auto &ifix : modify->get_fix_list()) + if (ifix->stores_ids) + error->all(FLERR, "Cannot use reset_atoms id with a fix {} storing atom IDs", ifix->style); if (comm->me == 0) utils::logmesg(lmp, "Resetting atom IDs ...\n"); From 27aa6898f8b137d3db9e6504774f1a83f64d82fd Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 16 Jul 2023 20:05:52 -0400 Subject: [PATCH 8/9] modernize access for fixes and computes --- src/atom.cpp | 13 ++++++------- src/atom_vec_sphere.cpp | 15 +++++++++++---- src/compute_aggregate_atom.cpp | 6 ++---- src/compute_centro_atom.cpp | 8 +++----- src/compute_cluster_atom.cpp | 6 ++---- src/compute_cna_atom.cpp | 14 +++++--------- src/compute_coord_atom.cpp | 16 ++++++++++------ src/compute_erotate_sphere_atom.cpp | 7 ++----- src/compute_fragment_atom.cpp | 7 ++----- src/compute_ke_atom.cpp | 6 ++---- 10 files changed, 45 insertions(+), 53 deletions(-) diff --git a/src/atom.cpp b/src/atom.cpp index 6b014fc9c8..c4521a244e 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -1818,17 +1818,16 @@ void Atom::data_bodies(int n, char *buf, AtomVec *avec_body, tagint id_offset) void Atom::data_fix_compute_variable(int nprev, int nnew) { - for (const auto &fix : modify->get_fix_list()) { - if (fix->create_attribute) + for (const auto &ifix : modify->get_fix_list()) { + if (ifix->create_attribute) for (int i = nprev; i < nnew; i++) - fix->set_arrays(i); + ifix->set_arrays(i); } - for (int m = 0; m < modify->ncompute; m++) { - Compute *compute = modify->compute[m]; - if (compute->create_attribute) + for (const auto &icompute : modify->get_compute_list()) { + if (icompute->create_attribute) for (int i = nprev; i < nnew; i++) - compute->set_arrays(i); + icompute->set_arrays(i); } for (int i = nprev; i < nnew; i++) diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp index d2f9c51880..27ea4b6fef 100644 --- a/src/atom_vec_sphere.cpp +++ b/src/atom_vec_sphere.cpp @@ -17,6 +17,7 @@ #include "error.h" #include "fix.h" #include "fix_adapt.h" +#include "fix_adapt_fep.h" #include "math_const.h" #include "modify.h" @@ -88,12 +89,18 @@ void AtomVecSphere::init() // check if optional radvary setting should have been set to 1 - for (int i = 0; i < modify->nfix; i++) - if (strcmp(modify->fix[i]->style, "adapt") == 0) { - auto fix = dynamic_cast(modify->fix[i]); - if (fix->diamflag && radvary == 0) + for (auto &ifix : modify->get_fix_by_style("^adapt")) { + if (utils::strmatch(ifix->style, "^adapt$")) { + auto fix = dynamic_cast(ifix); + if (fix && fix->diamflag && radvary == 0) error->all(FLERR, "Fix adapt changes particle radii but atom_style sphere is not dynamic"); } + if (utils::strmatch(ifix->style, "^adapt/fep$")) { + auto fix = dynamic_cast(ifix); + if (fix && fix->diamflag && radvary == 0) + error->all(FLERR, "Fix adapt/fep changes particle radii but atom_style sphere is not dynamic"); + } + } } /* ---------------------------------------------------------------------- diff --git a/src/compute_aggregate_atom.cpp b/src/compute_aggregate_atom.cpp index 05526494c0..e68062b73e 100644 --- a/src/compute_aggregate_atom.cpp +++ b/src/compute_aggregate_atom.cpp @@ -82,10 +82,8 @@ void ComputeAggregateAtom::init() neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style, "aggregate/atom") == 0) count++; - if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute aggregate/atom"); + if (modify->get_compute_by_style(style).size() > 1) + if (comm->me == 0) error->warning(FLERR, "More than one compute {}", style); } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_centro_atom.cpp b/src/compute_centro_atom.cpp index f726a95594..10d2a79c45 100644 --- a/src/compute_centro_atom.cpp +++ b/src/compute_centro_atom.cpp @@ -93,14 +93,12 @@ void ComputeCentroAtom::init() if (force->pair == nullptr) error->all(FLERR, "Compute centro/atom requires a pair style be defined"); - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style, "centro/atom") == 0) count++; - if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute centro/atom"); - // need an occasional full neighbor list neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); + + if (modify->get_compute_by_style(style).size() > 1) + if (comm->me == 0) error->warning(FLERR, "More than one compute {}", style); } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_cluster_atom.cpp b/src/compute_cluster_atom.cpp index 0a0a36debd..ae44fbcd37 100644 --- a/src/compute_cluster_atom.cpp +++ b/src/compute_cluster_atom.cpp @@ -69,10 +69,8 @@ void ComputeClusterAtom::init() neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style, "cluster/atom") == 0) count++; - if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute cluster/atom"); + if (modify->get_compute_by_style(style).size() > 1) + if (comm->me == 0) error->warning(FLERR, "More than one compute {}", style); } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_cna_atom.cpp b/src/compute_cna_atom.cpp index ccb4183707..a09a671c07 100644 --- a/src/compute_cna_atom.cpp +++ b/src/compute_cna_atom.cpp @@ -76,19 +76,15 @@ void ComputeCNAAtom::init() // cannot use neighbor->cutneighmax b/c neighbor has not yet been init - if (2.0 * sqrt(cutsq) > force->pair->cutforce + neighbor->skin && comm->me == 0) - error->warning(FLERR, - "Compute cna/atom cutoff may be too large to find " - "ghost atom neighbors"); - - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style, "cna/atom") == 0) count++; - if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute cna/atom defined"); + if ((2.0 * sqrt(cutsq)) > (force->pair->cutforce + neighbor->skin) && (comm->me == 0)) + error->warning(FLERR, "Compute cna/atom cutoff may be too large to find ghost atom neighbors"); // need an occasional full neighbor list neighbor->add_request(this, NeighConst::REQ_FULL | NeighConst::REQ_OCCASIONAL); + + if (modify->get_compute_by_style(style).size() > 1) + if (comm->me == 0) error->warning(FLERR, "More than one compute {}", style); } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp index 0d56ea6573..ea19f01b91 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -82,10 +82,11 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : id_orientorder = utils::strdup(arg[4]); - int iorientorder = modify->find_compute(id_orientorder); - if (iorientorder < 0) error->all(FLERR, "Could not find compute coord/atom compute ID"); - if (!utils::strmatch(modify->compute[iorientorder]->style, "^orientorder/atom")) - error->all(FLERR, "Compute coord/atom compute ID is not orientorder/atom"); + auto iorientorder = modify->get_compute_by_id(id_orientorder); + if (!iorientorder) + error->all(FLERR, "Could not find compute coord/atom compute ID {}", id_orientorder); + if (!utils::strmatch(iorientorder->style, "^orientorder/atom")) + error->all(FLERR, "Compute coord/atom compute ID {} is not orientorder/atom", id_orientorder); threshold = utils::numeric(FLERR, arg[5], false, lmp); if (threshold <= -1.0 || threshold >= 1.0) @@ -128,8 +129,11 @@ ComputeCoordAtom::~ComputeCoordAtom() void ComputeCoordAtom::init() { if (cstyle == ORIENT) { - int iorientorder = modify->find_compute(id_orientorder); - c_orientorder = dynamic_cast(modify->compute[iorientorder]); + c_orientorder = + dynamic_cast(modify->get_compute_by_id(id_orientorder)); + if (!c_orientorder) + error->all(FLERR, "Could not find compute coord/atom compute ID {}", id_orientorder); + cutsq = c_orientorder->cutsq; l = c_orientorder->qlcomp; // communicate real and imaginary 2*l+1 components of the normalized vector diff --git a/src/compute_erotate_sphere_atom.cpp b/src/compute_erotate_sphere_atom.cpp index 9bfa0afaea..3ec0f402a8 100644 --- a/src/compute_erotate_sphere_atom.cpp +++ b/src/compute_erotate_sphere_atom.cpp @@ -57,11 +57,8 @@ ComputeErotateSphereAtom::~ComputeErotateSphereAtom() void ComputeErotateSphereAtom::init() { - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"erotate/sphere/atom") == 0) count++; - if (count > 1 && comm->me == 0) - error->warning(FLERR,"More than one compute erotate/sphere/atom"); + if (modify->get_compute_by_style(style).size() > 1) + if (comm->me == 0) error->warning(FLERR, "More than one compute {}", style); pfactor = 0.5 * force->mvv2e * INERTIA; } diff --git a/src/compute_fragment_atom.cpp b/src/compute_fragment_atom.cpp index 810ffa4d66..035f554c8d 100644 --- a/src/compute_fragment_atom.cpp +++ b/src/compute_fragment_atom.cpp @@ -84,11 +84,8 @@ void ComputeFragmentAtom::init() if (atom->molecular != Atom::MOLECULAR) error->all(FLERR,"Compute fragment/atom requires a molecular system"); - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style,"fragment/atom") == 0) count++; - if (count > 1 && comm->me == 0) - error->warning(FLERR,"More than one compute fragment/atom"); + if (modify->get_compute_by_style(style).size() > 1) + if (comm->me == 0) error->warning(FLERR, "More than one compute {}", style); } /* ---------------------------------------------------------------------- */ diff --git a/src/compute_ke_atom.cpp b/src/compute_ke_atom.cpp index aab687755e..9a329232b3 100644 --- a/src/compute_ke_atom.cpp +++ b/src/compute_ke_atom.cpp @@ -47,10 +47,8 @@ ComputeKEAtom::~ComputeKEAtom() void ComputeKEAtom::init() { - int count = 0; - for (int i = 0; i < modify->ncompute; i++) - if (strcmp(modify->compute[i]->style, "ke/atom") == 0) count++; - if (count > 1 && comm->me == 0) error->warning(FLERR, "More than one compute ke/atom"); + if (modify->get_compute_by_style(style).size() > 1) + if (comm->me == 0) error->warning(FLERR, "More than one compute {}", style); } /* ---------------------------------------------------------------------- */ From 3568cced4bd6b5c42841a9c61ee658ae3031d325 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 17 Jul 2023 09:52:44 -0400 Subject: [PATCH 9/9] update comment --- src/REPLICA/temper.h | 2 +- src/REPLICA/temper_grem.h | 2 +- src/REPLICA/temper_npt.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/REPLICA/temper.h b/src/REPLICA/temper.h index d959892587..46a6666f07 100644 --- a/src/REPLICA/temper.h +++ b/src/REPLICA/temper.h @@ -40,7 +40,7 @@ class Temper : public Command { int nswaps; // # of tempering swaps to perform int seed_swap; // 0 = toggle swaps, n = RNG for swap direction int seed_boltz; // seed for Boltz factor comparison - class Fix *whichfix; // index of temperature fix to use + class Fix *whichfix; // temperature fix to use int my_set_temp; // which set temp I am simulating double *set_temp; // static list of replica set temperatures diff --git a/src/REPLICA/temper_grem.h b/src/REPLICA/temper_grem.h index c66239a794..8fe49ae804 100644 --- a/src/REPLICA/temper_grem.h +++ b/src/REPLICA/temper_grem.h @@ -40,7 +40,7 @@ class TemperGrem : public Command { int nswaps; // # of tempering swaps to perform int seed_swap; // 0 = toggle swaps, n = RNG for swap direction int seed_boltz; // seed for Boltz factor comparison - class Fix *whichfix; // index of temperature fix to use + class Fix *whichfix; // temperature fix to use int my_set_lambda; // which set lambda I am simulating double *set_lambda; // static list of replica set lambdas diff --git a/src/REPLICA/temper_npt.h b/src/REPLICA/temper_npt.h index 2e831c55ef..3ea29a21f6 100644 --- a/src/REPLICA/temper_npt.h +++ b/src/REPLICA/temper_npt.h @@ -42,7 +42,7 @@ class TemperNPT : public Command { int nswaps; // # of tempering swaps to perform int seed_swap; // 0 = toggle swaps, n = RNG for swap direction int seed_boltz; // seed for Boltz factor comparison - class Fix *whichfix; // index of temperature fix to use + class Fix *whichfix; // temperature fix to use int my_set_temp; // which set temp I am simulating double *set_temp; // static list of replica set temperatures