Merge pull request #1 from amilumas/temper_npt

# Solvated 5-mer peptide

# Demonstrating temper_npt

units		real

atom_style	full

pair_style	lj/charmm/coul/long 8.0 10.0 10.0

bond_style      harmonic

angle_style     charmm

dihedral_style  charmm

improper_style  harmonic

kspace_style	pppm 0.0001

read_data	data.peptide

neighbor	2.0 bin

neigh_modify	delay 5

timestep	2.0

thermo_style	custom step temp epair emol etotal press density

thermo		50

variable	temper_T world 275 280 285 290 295 300 305 310

variable 	rep world 0 1 2 3 4 5 6 7 

fix		myfix all npt temp ${temper_T} ${temper_T} 100.0 iso 1 1 1000

run		500

temper_npt	2000 100 ${temper_T} myfix 0 58728 1 

fix		2 all shake 0.0001 10 100 b 4 6 8 10 12 14 18 a 31

group		peptide type <= 12

/* ----------------------------------------------------------------------

   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 <math.h>

#include <stdlib.h>

#include <string.h>

#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);

  }

}

/* -*- 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.

*/