diff --git a/src/REPLICA/temper_npt.cpp b/src/REPLICA/temper_npt.cpp deleted file mode 100644 index a043c6140..000000000 --- a/src/REPLICA/temper_npt.cpp +++ /dev/null @@ -1,381 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing Authors: Amulya K. Pervaje and Cody K. Addington - Contact Email: amulyapervaje@gmail.com - temper_npt is a modification of temper that is applicable to the NPT ensemble - uses the npt acceptance criteria for parallel tempering (replica exchange) as given in - Mori, Y .; Okamoto, Y . Generalized-Ensemble Algorithms for the Isobaric–Isothermal Ensemble. J. Phys. Soc. Japan 2010, 79, 74003. - - temper_npt N M temp fix-ID seed1 seed2 pressure index(optional) - refer to documentation for temper, only difference with temper_npt is that the pressure is specified as the 7th argument, the 8th argument is the same optional index argument used in temper -------------------------------------------------------------------------- */ - -#include -#include -#include -#include "temper_npt.h" -#include "universe.h" -#include "domain.h" -#include "atom.h" -#include "update.h" -#include "integrate.h" -#include "modify.h" -#include "compute.h" -#include "force.h" -#include "output.h" -#include "thermo.h" -#include "fix.h" -#include "random_park.h" -#include "finish.h" -#include "timer.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -#define TEMPER_DEBUG 0 - -/* ---------------------------------------------------------------------- */ - -TemperNpt::TemperNpt(LAMMPS *lmp) : Pointers(lmp) {} - -/* ---------------------------------------------------------------------- */ - -TemperNpt::~TemperNpt() -{ - MPI_Comm_free(&roots); - if (ranswap) delete ranswap; - delete ranboltz; - delete [] set_temp; - delete [] temp2world; - delete [] world2temp; - delete [] world2root; -} - -/* ---------------------------------------------------------------------- - perform tempering with inter-world swaps -------------------------------------------------------------------------- */ - -void TemperNpt::command(int narg, char **arg) -{ - if (universe->nworlds == 1) - 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 command"); - - int nsteps = force->inumeric(FLERR,arg[0]); - nevery = force->inumeric(FLERR,arg[1]); - double temp = force->numeric(FLERR,arg[2]); - double press_set = force->numeric(FLERR,arg[6]); - - 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"); - - seed_swap = force->inumeric(FLERR,arg[4]); - seed_boltz = force->inumeric(FLERR,arg[5]); - - my_set_temp = universe->iworld; - if (narg == 8) my_set_temp = force->inumeric(FLERR,arg[6]); - - // swap frequency must evenly divide total # of timesteps - - if (nevery <= 0) - error->universe_all(FLERR,"Invalid frequency in temper command"); - nswaps = nsteps/nevery; - if (nswaps*nevery != nsteps) - error->universe_all(FLERR,"Non integer # of swaps in temper command"); - - // 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 ((strcmp(modify->fix[whichfix]->style,"npt") != 0)) error->universe_all(FLERR,"Tempering temperature fix is not supported"); - - // setup for long tempering run - - update->whichflag = 1; - 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"); - - lmp->init(); - - // local storage - - me_universe = universe->me; - MPI_Comm_rank(world,&me); - nworlds = universe->nworlds; - iworld = universe->iworld; - boltz = force->boltz; - nktv2p = force->nktv2p; - - // 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]; - pe_compute->addstep(update->ntimestep + nevery); - - // create MPI communicator for root proc from each world - - int color; - if (me == 0) color = 0; - else color = 1; - MPI_Comm_split(universe->uworld,color,0,&roots); - - // RNGs for swaps and Boltzmann test - // warm up Boltzmann RNG - - if (seed_swap) ranswap = new RanPark(lmp,seed_swap); - else ranswap = NULL; - ranboltz = new RanPark(lmp,seed_boltz + me_universe); - for (int i = 0; i < 100; i++) ranboltz->uniform(); - - // world2root[i] = global proc that is root proc of world i - - world2root = new int[nworlds]; - if (me == 0) - MPI_Allgather(&me_universe,1,MPI_INT,world2root,1,MPI_INT,roots); - MPI_Bcast(world2root,nworlds,MPI_INT,0,world); - - // create static list of set temperatures - // allgather tempering arg "temp" across root procs - // bcast from each root to other procs in world - - set_temp = new double[nworlds]; - if (me == 0) MPI_Allgather(&temp,1,MPI_DOUBLE,set_temp,1,MPI_DOUBLE,roots); - MPI_Bcast(set_temp,nworlds,MPI_DOUBLE,0,world); - - // create world2temp only on root procs from my_set_temp - // create temp2world on root procs from world2temp, - // then bcast to all procs within world - - world2temp = new int[nworlds]; - temp2world = new int[nworlds]; - if (me == 0) { - MPI_Allgather(&my_set_temp,1,MPI_INT,world2temp,1,MPI_INT,roots); - for (int i = 0; i < nworlds; i++) temp2world[world2temp[i]] = i; - } - MPI_Bcast(temp2world,nworlds,MPI_INT,0,world); - - // if restarting tempering, reset temp target of Fix to current my_set_temp - - if (narg == 8) { - double new_temp = set_temp[my_set_temp]; - modify->fix[whichfix]->reset_target(new_temp); - } - - // setup tempering runs - - int i,which,partner,swap,partner_set_temp,partner_world; - double pe,pe_partner, delr,boltz_factor,new_temp, press_units; - - if (me_universe == 0 && universe->uscreen) - fprintf(universe->uscreen,"Setting up tempering ...\n"); - - update->integrate->setup(); - - if (me_universe == 0) { - if (universe->uscreen) { - fprintf(universe->uscreen,"Step"); - for (int i = 0; i < nworlds; i++) - fprintf(universe->uscreen," T%d",i); - fprintf(universe->uscreen,"\n"); - } - if (universe->ulogfile) { - fprintf(universe->ulogfile,"Step"); - for (int i = 0; i < nworlds; i++) - fprintf(universe->ulogfile," T%d",i); - fprintf(universe->ulogfile,"\n"); - } - print_status(); - } - - timer->init(); - timer->barrier_start(); - - for (int iswap = 0; iswap < nswaps; iswap++) { - - // run for nevery timesteps - - update->integrate->run(nevery); - - // compute PE - // notify compute it will be called at next swap - - pe = pe_compute->compute_scalar(); - pe_compute->addstep(update->ntimestep + nevery); - double boxlox=domain->boxlo[0]; - double boxhix=domain->boxhi[0]; - double boxloy=domain->boxlo[1]; - double boxhiy=domain->boxhi[1]; - double boxloz=domain->boxlo[2]; - double boxhiz=domain->boxhi[2]; - double vol = (boxhix - boxlox)*(boxhiy - boxloy)*(boxhiz - boxloz); - double vol_partner = vol; - // which = which of 2 kinds of swaps to do (0,1) - - if (!ranswap) which = iswap % 2; - else if (ranswap->uniform() < 0.5) which = 0; - else which = 1; - - // partner_set_temp = which set temp I am partnering with for this swap - - if (which == 0) { - if (my_set_temp % 2 == 0) partner_set_temp = my_set_temp + 1; - else partner_set_temp = my_set_temp - 1; - } else { - if (my_set_temp % 2 == 1) partner_set_temp = my_set_temp + 1; - else partner_set_temp = my_set_temp - 1; - } - - // partner = proc ID to swap with - // if partner = -1, then I am not a proc that swaps - - partner = -1; - if (me == 0 && partner_set_temp >= 0 && partner_set_temp < nworlds) { - partner_world = temp2world[partner_set_temp]; - partner = world2root[partner_world]; - } - - // swap with a partner, only root procs in each world participate - // hi proc sends PE to low proc - // lo proc make Boltzmann decision on whether to swap - // lo proc communicates decision back to hi proc - - swap = 0; - if (partner != -1) { - if (me_universe > partner) { - MPI_Send(&pe,1,MPI_DOUBLE,partner,0,universe->uworld); - } - else { - MPI_Recv(&pe_partner,1,MPI_DOUBLE,partner,0,universe->uworld,MPI_STATUS_IGNORE); - } - if (me_universe > partner) { - MPI_Send(&vol,1, MPI_DOUBLE,partner,0,universe->uworld); - } - else { - MPI_Recv(&vol_partner,1,MPI_DOUBLE,partner,0,universe->uworld,MPI_STATUS_IGNORE); - } - // Acceptance criteria changed 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); - boltz_factor = -delr; - if (boltz_factor >= 0.0) swap = 1; - else if (ranboltz->uniform() < exp(boltz_factor)) swap = 1; - } - - if (me_universe < partner) - 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 (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 - - MPI_Bcast(&swap,1,MPI_INT,0,world); - - // rescale kinetic energy via velocities if move is accepted - - if (swap) scale_velocities(partner_set_temp,my_set_temp); - - // if my world swapped, all procs in world reset temp target of Fix - - if (swap) { - new_temp = set_temp[partner_set_temp]; - modify->fix[whichfix]->reset_target(new_temp); - } - - // update my_set_temp and temp2world on every proc - // root procs update their value if swap took place - // allgather across root procs - // bcast within my world - - if (swap) my_set_temp = partner_set_temp; - if (me == 0) { - MPI_Allgather(&my_set_temp,1,MPI_INT,world2temp,1,MPI_INT,roots); - for (i = 0; i < nworlds; i++) temp2world[world2temp[i]] = i; - } - MPI_Bcast(temp2world,nworlds,MPI_INT,0,world); - - // print out current swap status - - if (me_universe == 0) print_status(); - } - - timer->barrier_stop(); - - update->integrate->cleanup(); - - Finish finish(lmp); - finish.end(1); - - update->whichflag = 0; - update->firststep = update->laststep = 0; - update->beginstep = update->endstep = 0; -} - -/* ---------------------------------------------------------------------- - scale kinetic energy via velocities a la Sugita -------------------------------------------------------------------------- */ - -void TemperNpt::scale_velocities(int t_partner, int t_me) -{ - double sfactor = sqrt(set_temp[t_partner]/set_temp[t_me]); - - double **v = atom->v; - int nlocal = atom->nlocal; - - for (int i = 0; i < nlocal; i++) { - v[i][0] = v[i][0]*sfactor; - v[i][1] = v[i][1]*sfactor; - v[i][2] = v[i][2]*sfactor; - } -} - -/* ---------------------------------------------------------------------- - proc 0 prints current tempering status -------------------------------------------------------------------------- */ - -void TemperNpt::print_status() -{ - if (universe->uscreen) { - fprintf(universe->uscreen,BIGINT_FORMAT,update->ntimestep); - for (int i = 0; i < nworlds; i++) - fprintf(universe->uscreen," %d",world2temp[i]); - fprintf(universe->uscreen,"\n"); - } - if (universe->ulogfile) { - fprintf(universe->ulogfile,BIGINT_FORMAT,update->ntimestep); - for (int i = 0; i < nworlds; i++) - fprintf(universe->ulogfile," %d",world2temp[i]); - fprintf(universe->ulogfile,"\n"); - fflush(universe->ulogfile); - } -} diff --git a/src/REPLICA/temper_npt.h b/src/REPLICA/temper_npt.h deleted file mode 100644 index 984f2839b..000000000 --- a/src/REPLICA/temper_npt.h +++ /dev/null @@ -1,118 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. - - Contributing Authors: Amulya K. Pervaje and Cody K. Addington - Contact Email: amulyapervaje@gmail.com - temper_npt is a modification of temper that is applicable to the NPT ensemble - uses the npt acceptance criteria for parallel tempering (replica exchange) as given in - Mori, Y .; Okamoto, Y . Generalized-Ensemble Algorithms for the Isobaric–Isothermal Ensemble. J. Phys. Soc. Japan 2010, 79, 74003. - - temper_npt N M temp fix-ID seed1 seed2 pressure index(optional) - refer to documentation for temper, only difference with temper_npt is that the pressure is s -pecified as the 7th argument, the 8th argument is the same optional index argument used in temp -er -------------------------------------------------------------------------- */ - -#ifdef COMMAND_CLASS - -CommandStyle(temper_npt,TemperNpt) - -#else - -#ifndef LMP_TEMPERNPT_H -#define LMP_TEMPERNPT_H - -#include "pointers.h" - -namespace LAMMPS_NS { - -class TemperNpt : protected Pointers { - public: - TemperNpt(class LAMMPS *); - ~TemperNpt(); - void command(int, char **); - - private: - int me,me_universe; // my proc ID in world and universe - int iworld,nworlds; // world info - double boltz; // copy from output->boltz - double nktv2p; - MPI_Comm roots; // MPI comm with 1 root proc from each world - class RanPark *ranswap,*ranboltz; // RNGs for swapping and Boltz factor - int nevery; // # of timesteps between swaps - 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 - - int my_set_temp; // which set temp I am simulating - double *set_temp; // static list of replica set temperatures - int *temp2world; // temp2world[i] = world simulating set temp i - int *world2temp; // world2temp[i] = temp simulated by world i - int *world2root; // world2root[i] = root proc of world i - - void scale_velocities(int, int); - void print_status(); -}; - -} - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Must have more than one processor partition to temper - -Cannot use the temper command with only one processor partition. Use -the -partition command-line option. - -E: TemperNpt command before simulation box is defined - -The temper command cannot be used before a read_data, read_restart, or -create_box command. - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Tempering fix ID is not defined - -The fix ID specified by the temper command does not exist. - -E: Invalid frequency in temper command - -Nevery must be > 0. - -E: Non integer # of swaps in temper command - -Swap frequency in temper command must evenly divide the total # of -timesteps. - -E: Tempering temperature fix is not valid - -The fix specified by the temper command is not one that controls -temperature (nvt or langevin). - -E: Too many timesteps - -The cummulative timesteps must fit in a 64-bit integer. - -E: Tempering could not find thermo_pe compute - -This compute is created by the thermo command. It must have been -explicitly deleted by a uncompute command. - -*/