diff --git a/src/min_fire.cpp b/src/min_fire.cpp new file mode 100644 index 0000000000..5b54be5ec0 --- /dev/null +++ b/src/min_fire.cpp @@ -0,0 +1,192 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "min_fire.h" +#include "atom.h" +#include "force.h" +#include "update.h" +#include "output.h" +#include "timer.h" +#include "error.h" + +using namespace LAMMPS_NS; + +// same as in other min classes + +enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD}; + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +#define DELAYSTEP 5 +#define DT_GROW 1.1 +#define DT_SHRINK 0.5 +#define ALPHA0 0.1 +#define ALPHA_SHRINK 0.99 +#define TMAX 10.0 + +/* ---------------------------------------------------------------------- */ + +MinFire::MinFire(LAMMPS *lmp) : Min(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void MinFire::init_style() +{ + dt = update->dt; +} + +/* ---------------------------------------------------------------------- */ + +void MinFire::setup_style() +{ + double **v = atom->v; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + v[i][0] = v[i][1] = v[i][2] = 0.0; +} + +/* ---------------------------------------------------------------------- + set current vector lengths and pointers + called after atoms have migrated +------------------------------------------------------------------------- */ + +void MinFire::reset_vectors() +{ + // atomic dof + + nvec = 3 * atom->nlocal; + if (nvec) xvec = atom->x[0]; + if (nvec) fvec = atom->f[0]; +} + +/* ---------------------------------------------------------------------- */ + +int MinFire::iterate(int maxiter) +{ + int ntimestep; + double vmax,vdotf,vdotfall,vdotv,vdotvall,fdotf,fdotfall; + double scale1,scale2; + double dtvone,dtv,dtfm; + + alpha_final = 0.0; + double alpha = ALPHA0; + double dtmax = TMAX * dt; + int last_negative = update->ntimestep; + + for (int iter = 0; iter < maxiter; iter++) { + ntimestep = ++update->ntimestep; + niter++; + + // vdotfall = v dot f + + double **v = atom->v; + double **f = atom->f; + int nlocal = atom->nlocal; + + vdotf = 0.0; + for (int i = 0; i < nlocal; i++) + vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2]; + MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world); + + // if (v dot f) > 0: + // v = (1-alpha) v + alpha |v| Fhat + // |v| = length of v, Fhat = unit f + // if more than DELAYSTEP since v dot f was negative: + // increase timestep and decrease alpha + + if (vdotfall > 0.0) { + vdotv = 0.0; + for (int i = 0; i < nlocal; i++) + vdotv += v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; + MPI_Allreduce(&vdotv,&vdotvall,1,MPI_DOUBLE,MPI_SUM,world); + fdotf = 0.0; + for (int i = 0; i < nlocal; i++) + fdotf += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2]; + MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world); + + scale1 = 1.0 - alpha; + if (fdotfall == 0.0) scale2 = 0.0; + else scale2 = alpha * sqrt(vdotvall/fdotfall); + for (int i = 0; i < nlocal; i++) { + v[i][0] = scale1*v[i][0] + scale2*f[i][0]; + v[i][1] = scale1*v[i][1] + scale2*f[i][1]; + v[i][2] = scale1*v[i][2] + scale2*f[i][2]; + } + + if (ntimestep - last_negative > DELAYSTEP) { + dt = MIN(dt*DT_GROW,dtmax); + alpha *= ALPHA_SHRINK; + } + + // else (v dot f) <= 0: + // decrease timestep, reset alpha, set v = 0 + + } else { + last_negative = ntimestep; + dt *= DT_SHRINK; + alpha = ALPHA0; + for (int i = 0; i < nlocal; i++) + v[i][0] = v[i][1] = v[i][2] = 0.0; + } + + // limit timestep so no particle moves further than dmax + + double *mass = atom->mass; + int *type = atom->type; + + dtvone = dt; + + for (int i = 0; i < nlocal; i++) { + vmax = MAX(fabs(v[i][0]),fabs(v[i][1])); + vmax = MAX(vmax,fabs(v[i][2])); + if (dtvone*vmax > dmax) dtvone = dmax/vmax; + } + MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world); + + // Euler integration step + + double **x = atom->x; + + for (int i = 0; i < nlocal; i++) { + dtfm = dtv / mass[type[i]]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + } + + eprevious = ecurrent; + ecurrent = energy_force(0); + neval++; + + // force tolerance criterion + + //fdotf = fnorm_sqr(); + //if (fdotf < update->ftol*update->ftol) return FTOL; + + // output for thermo, dump, restart files + + if (output->next == ntimestep) { + timer->stamp(); + output->write(ntimestep); + timer->stamp(TIME_OUTPUT); + } + } + + return MAXITER; +} diff --git a/src/min_fire.h b/src/min_fire.h new file mode 100644 index 0000000000..0160ab4afb --- /dev/null +++ b/src/min_fire.h @@ -0,0 +1,43 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef MINIMIZE_CLASS + +MinimizeStyle(fire,MinFire) + +#else + +#ifndef LMP_MIN_FIRE_H +#define LMP_MIN_FIRE_H + +#include "min.h" + +namespace LAMMPS_NS { + +class MinFire : public Min { + public: + MinFire(class LAMMPS *); + ~MinFire() {} + void init_style(); + void setup_style(); + void reset_vectors(); + int iterate(int); + + private: + double dt; +}; + +} + +#endif +#endif diff --git a/src/min_quickmin.cpp b/src/min_quickmin.cpp new file mode 100644 index 0000000000..dd4fdeb8b2 --- /dev/null +++ b/src/min_quickmin.cpp @@ -0,0 +1,159 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#include "math.h" +#include "min_quickmin.h" +#include "atom.h" +#include "force.h" +#include "update.h" +#include "output.h" +#include "timer.h" +#include "error.h" + +using namespace LAMMPS_NS; + +// same as in other min classes + +enum{MAXITER,MAXEVAL,ETOL,FTOL,DOWNHILL,ZEROALPHA,ZEROFORCE,ZEROQUAD}; + +#define MIN(A,B) ((A) < (B)) ? (A) : (B) +#define MAX(A,B) ((A) > (B)) ? (A) : (B) + +/* ---------------------------------------------------------------------- */ + +MinQuickmin::MinQuickmin(LAMMPS *lmp) : Min(lmp) {} + +/* ---------------------------------------------------------------------- */ + +void MinQuickmin::init_style() +{ + dt = update->dt; +} + +/* ---------------------------------------------------------------------- */ + +void MinQuickmin::setup_style() +{ + double **v = atom->v; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + v[i][0] = v[i][1] = v[i][2] = 0.0; +} + +/* ---------------------------------------------------------------------- + set current vector lengths and pointers + called after atoms have migrated +------------------------------------------------------------------------- */ + +void MinQuickmin::reset_vectors() +{ + // atomic dof + + nvec = 3 * atom->nlocal; + if (nvec) xvec = atom->x[0]; + if (nvec) fvec = atom->f[0]; +} + +/* ---------------------------------------------------------------------- */ + +int MinQuickmin::iterate(int maxiter) +{ + int ntimestep; + double vmax,vdotf,vdotfall,fdotf,fdotfall,scale; + double dtvone,dtv,dtfm; + + alpha_final = 0.0; + + for (int iter = 0; iter < maxiter; iter++) { + ntimestep = ++update->ntimestep; + niter++; + + // zero velocity if anti-parallel to force + // else project velocity in direction of force + + double **v = atom->v; + double **f = atom->f; + int nlocal = atom->nlocal; + + vdotf = 0.0; + for (int i = 0; i < nlocal; i++) + vdotf += v[i][0]*f[i][0] + v[i][1]*f[i][1] + v[i][2]*f[i][2]; + MPI_Allreduce(&vdotf,&vdotfall,1,MPI_DOUBLE,MPI_SUM,world); + + if (vdotfall < 0.0) { + for (int i = 0; i < nlocal; i++) + v[i][0] = v[i][1] = v[i][2] = 0.0; + + } else { + fdotf = 0.0; + for (int i = 0; i < nlocal; i++) + fdotf += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2]; + MPI_Allreduce(&fdotf,&fdotfall,1,MPI_DOUBLE,MPI_SUM,world); + if (fdotfall == 0.0) scale = 0.0; + else scale = vdotfall/fdotfall; + for (int i = 0; i < nlocal; i++) { + v[i][0] = scale * f[i][0]; + v[i][1] = scale * f[i][1]; + v[i][2] = scale * f[i][2]; + } + } + + // limit timestep so no particle moves further than dmax + + double *mass = atom->mass; + int *type = atom->type; + + dtvone = dt; + + for (int i = 0; i < nlocal; i++) { + vmax = MAX(fabs(v[i][0]),fabs(v[i][1])); + vmax = MAX(vmax,fabs(v[i][2])); + if (dtvone*vmax > dmax) dtvone = dmax/vmax; + } + MPI_Allreduce(&dtvone,&dtv,1,MPI_DOUBLE,MPI_MIN,world); + + // Euler integration step + + double **x = atom->x; + + for (int i = 0; i < nlocal; i++) { + dtfm = dtv / mass[type[i]]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + } + + eprevious = ecurrent; + ecurrent = energy_force(0); + neval++; + + // force tolerance criterion + + //fdotf = fnorm_sqr(); + //if (fdotf < update->ftol*update->ftol) return FTOL; + + // output for thermo, dump, restart files + + if (output->next == ntimestep) { + timer->stamp(); + output->write(ntimestep); + timer->stamp(TIME_OUTPUT); + } + } + + return MAXITER; +} diff --git a/src/min_quickmin.h b/src/min_quickmin.h new file mode 100644 index 0000000000..b6362922fa --- /dev/null +++ b/src/min_quickmin.h @@ -0,0 +1,43 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef MINIMIZE_CLASS + +MinimizeStyle(quickmin,MinQuickmin) + +#else + +#ifndef LMP_MIN_QUICKMIN_H +#define LMP_MIN_QUICKMIN_H + +#include "min.h" + +namespace LAMMPS_NS { + +class MinQuickmin : public Min { + public: + MinQuickmin(class LAMMPS *); + ~MinQuickmin() {} + void init_style(); + void setup_style(); + void reset_vectors(); + int iterate(int); + + private: + double dt; +}; + +} + +#endif +#endif diff --git a/src/temper.cpp b/src/temper.cpp deleted file mode 100644 index 592e012f3e..0000000000 --- a/src/temper.cpp +++ /dev/null @@ -1,357 +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 author: Mark Sears (SNL) -------------------------------------------------------------------------- */ - -#include "math.h" -#include "stdlib.h" -#include "string.h" -#include "temper.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 1 - -/* ---------------------------------------------------------------------- */ - -Temper::Temper(LAMMPS *lmp) : Pointers(lmp) {} - -/* ---------------------------------------------------------------------- */ - -Temper::~Temper() -{ - 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 Temper::command(int narg, char **arg) -{ - if (universe->nworlds == 1) - error->all("Must have more than one processor partition to temper"); - if (domain->box_exist == 0) - error->all("Temper command before simulation box is defined"); - if (narg != 6 && narg != 7) error->universe_all("Illegal temper command"); - - int nsteps = atoi(arg[0]); - nevery = atoi(arg[1]); - double temp = atof(arg[2]); - - for (whichfix = 0; whichfix < modify->nfix; whichfix++) - if (strcmp(arg[3],modify->fix[whichfix]->id) == 0) break; - if (whichfix == modify->nfix) - error->universe_all("Tempering fix ID is not defined"); - - seed_swap = atoi(arg[4]); - seed_boltz = atoi(arg[5]); - - my_set_temp = universe->iworld; - if (narg == 7) my_set_temp = atoi(arg[6]); - - // swap frequency must evenly divide total # of timesteps - - if (nevery == 0) error->universe_all("Invalid frequency in temper command"); - nswaps = nsteps/nevery; - if (nswaps*nevery != nsteps) - error->universe_all("Non integer # of swaps in temper command"); - - // fix style must be appropriate for temperature control - - if ((strcmp(modify->fix[whichfix]->style,"nvt") != 0) && - (strcmp(modify->fix[whichfix]->style,"langevin") != 0) && - (strcmp(modify->fix[whichfix]->style,"temp/berendsen") != 0) && - (strcmp(modify->fix[whichfix]->style,"temp/rescale") != 0)) - error->universe_all("Tempering temperature fix is not valid"); - - // setup for long tempering run - - update->whichflag = 1; - update->nsteps = nsteps; - update->beginstep = update->firststep = update->ntimestep; - update->endstep = update->laststep = update->firststep + nsteps; - - lmp->init(); - - // local storage - - me_universe = universe->me; - MPI_Comm_rank(world,&me); - nworlds = universe->nworlds; - iworld = universe->iworld; - boltz = force->boltz; - - // 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("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 == 7) { - 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,boltz_factor,new_temp; - MPI_Status status; - - 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->barrier_start(TIME_LOOP); - - 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); - - // 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,&status); - - if (me_universe < partner) { - boltz_factor = (pe - pe_partner) * - (1.0/(boltz*set_temp[my_set_temp]) - - 1.0/(boltz*set_temp[partner_set_temp])); - 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,&status); - -#ifdef 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)); -#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(TIME_LOOP); - - 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 Temper::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 Temper::print_status() -{ - if (universe->uscreen) { - fprintf(universe->uscreen,"%d ",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,"%d ",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/temper.h b/src/temper.h deleted file mode 100644 index 58c3a7aa7e..0000000000 --- a/src/temper.h +++ /dev/null @@ -1,59 +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. -------------------------------------------------------------------------- */ - -#ifdef COMMAND_CLASS - -CommandStyle(temper,Temper) - -#else - -#ifndef LMP_TEMPER_H -#define LMP_TEMPER_H - -#include "pointers.h" - -namespace LAMMPS_NS { - -class Temper : protected Pointers { - public: - Temper(class LAMMPS *); - ~Temper(); - 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 - 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