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

modernize and consolidate style across all three temper command variants

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-07-15 12:37:12 -04:00
parent 77bdcb3e19
commit fd0a72eab5
5 changed files with 82 additions and 91 deletions
Download Patch File Download Diff File Expand all files Collapse all files

58
src/REPLICA/temper.cpp
View File

@ -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

3
src/REPLICA/temper.h
View File

@ -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

75
src/REPLICA/temper_grem.cpp
View File

@ -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<FixGrem*>(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<FixGrem*>(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<FixNH *>(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

4
src/REPLICA/temper_grem.h
View File

@ -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;
};

33
src/REPLICA/temper_npt.cpp
View File

@ -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