diff --git a/src/PRD/fix_event.cpp b/src/PRD/fix_event.cpp
index 30e5bfaf8..9f8e48b50 100644
--- a/src/PRD/fix_event.cpp
+++ b/src/PRD/fix_event.cpp
@@ -1,249 +1,249 @@
/* ----------------------------------------------------------------------
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: Mike Brown (SNL)
------------------------------------------------------------------------- */
#include "stdlib.h"
#include "string.h"
#include "fix_event.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "neighbor.h"
#include "comm.h"
#include "universe.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixEvent::FixEvent(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg != 3) error->all("Illegal fix event command");
restart_global = 1;
// perform initial allocation of atom-based array
// register with Atom class
xevent = NULL;
xold = NULL;
imageold = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
event_number = 0;
event_timestep = update->ntimestep;
clock = 0;
}
/* ---------------------------------------------------------------------- */
FixEvent::~FixEvent()
{
// unregister callbacks to this fix from Atom class
atom->delete_callback(id,0);
// delete locally stored array
memory->destroy_2d_double_array(xevent);
memory->destroy_2d_double_array(xold);
memory->sfree(imageold);
}
/* ---------------------------------------------------------------------- */
int FixEvent::setmask()
{
return 0;
}
/* ----------------------------------------------------------------------
save current atom coords as an event
called when an event occurs in some replica
set event_timestep = when event occurred in a particular replica
update clock = elapsed time since last event, across all replicas
------------------------------------------------------------------------- */
void FixEvent::store_event(int timestep, int delta_clock)
{
double **x = atom->x;
int *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
domain->unmap(x[i],image[i],xevent[i]);
event_timestep = timestep;
clock += delta_clock;
event_number++;
}
/* ----------------------------------------------------------------------
store state of all atoms
called before quench and subsequent check for event
so can later restore pre-quench state if no event occurs
------------------------------------------------------------------------- */
void FixEvent::store_state()
{
double **x = atom->x;
double **f = atom->f;
int *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
xold[i][0] = x[i][0];
xold[i][1] = x[i][1];
xold[i][2] = x[i][2];
imageold[i] = image[i];
}
}
/* ----------------------------------------------------------------------
restore state of all atoms to pre-quench state
called after no event detected so can continue
------------------------------------------------------------------------- */
void FixEvent::restore_state()
{
double **x = atom->x;
int *image = atom->image;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
x[i][0] = xold[i][0];
x[i][1] = xold[i][1];
x[i][2] = xold[i][2];
image[i] = imageold[i];
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double FixEvent::memory_usage()
{
double bytes = 6*atom->nmax * sizeof(double);
bytes += atom->nmax*sizeof(int);
return bytes;
}
/* ----------------------------------------------------------------------
allocate atom-based array
------------------------------------------------------------------------- */
void FixEvent::grow_arrays(int nmax)
{
xevent = memory->grow_2d_double_array(xevent,nmax,3,"event:xevent");
xold = memory->grow_2d_double_array(xold,nmax,3,"event:xold");
imageold = (int *)
memory->srealloc(imageold,nmax*sizeof(int),"event:imageold");
// allow compute event to access stored event coords
vector_atom = xevent;
}
/* ----------------------------------------------------------------------
copy values within local atom-based array
------------------------------------------------------------------------- */
void FixEvent::copy_arrays(int i, int j)
{
xevent[j][0] = xevent[i][0];
xevent[j][1] = xevent[i][1];
xevent[j][2] = xevent[i][2];
xold[j][0] = xold[i][0];
xold[j][1] = xold[i][1];
xold[j][2] = xold[i][2];
imageold[j] = imageold[i];
}
/* ----------------------------------------------------------------------
pack values in local atom-based array for exchange with another proc
------------------------------------------------------------------------- */
int FixEvent::pack_exchange(int i, double *buf)
{
buf[0] = xevent[i][0];
buf[1] = xevent[i][1];
buf[2] = xevent[i][2];
buf[3] = xold[i][0];
buf[4] = xold[i][1];
buf[5] = xold[i][2];
buf[6] = imageold[i];
return 7;
}
/* ----------------------------------------------------------------------
unpack values in local atom-based array from exchange with another proc
------------------------------------------------------------------------- */
int FixEvent::unpack_exchange(int nlocal, double *buf)
{
xevent[nlocal][0] = buf[0];
xevent[nlocal][1] = buf[1];
xevent[nlocal][2] = buf[2];
xold[nlocal][0] = buf[3];
xold[nlocal][1] = buf[4];
xold[nlocal][2] = buf[5];
imageold[nlocal] = static_cast<int>(buf[6]);
return 7;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixEvent::write_restart(FILE *fp)
{
int n = 0;
double list[5];
list[n++] = event_number;
list[n++] = event_timestep;
list[n++] = clock;
list[n++] = replica_number;
list[n++] = correlated_event;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
- fwrite(&list,sizeof(double),n,fp);
+ fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixEvent::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
event_number = static_cast<int> (list[n++]);
event_timestep = static_cast<int> (list[n++]);
clock = static_cast<int> (list[n++]);
replica_number = static_cast<int> (list[n++]);
correlated_event = static_cast<int> (list[n++]);
}
diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp
index 0e0aaaf67..8e0dc9ec8 100644
--- a/src/fix_deposit.cpp
+++ b/src/fix_deposit.cpp
@@ -1,441 +1,441 @@
/* ----------------------------------------------------------------------
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 "stdlib.h"
#include "string.h"
#include "fix_deposit.h"
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "domain.h"
#include "lattice.h"
#include "region.h"
#include "random_park.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all("Illegal fix deposit command");
restart_global = 1;
time_depend = 1;
// required args
ninsert = atoi(arg[3]);
ntype = atoi(arg[4]);
nfreq = atoi(arg[5]);
seed = atoi(arg[6]);
if (seed <= 0) error->all("Illegal fix deposit command");
// set defaults
iregion = -1;
globalflag = localflag = 0;
lo = hi = deltasq = 0.0;
nearsq = 0.0;
maxattempt = 10;
rateflag = 0;
vxlo = vxhi = vylo = vyhi = vzlo = vzhi = 0.0;
scaleflag = 1;
// read options from end of input line
options(narg-7,&arg[7]);
// error check on region and its extent being inside simulation box
if (iregion == -1) error->all("Must specify a region in fix deposit");
if (domain->regions[iregion]->interior == 0)
error->all("Must use region with side = in with fix deposit");
xlo = domain->regions[iregion]->extent_xlo;
xhi = domain->regions[iregion]->extent_xhi;
ylo = domain->regions[iregion]->extent_ylo;
yhi = domain->regions[iregion]->extent_yhi;
zlo = domain->regions[iregion]->extent_zlo;
zhi = domain->regions[iregion]->extent_zhi;
if (domain->triclinic == 0) {
if (xlo < domain->boxlo[0] || xhi > domain->boxhi[0] ||
ylo < domain->boxlo[1] || yhi > domain->boxhi[1] ||
zlo < domain->boxlo[2] || zhi > domain->boxhi[2])
error->all("Deposition region extends outside simulation box");
} else {
if (xlo < domain->boxlo_bound[0] || xhi > domain->boxhi_bound[0] ||
ylo < domain->boxlo_bound[1] || yhi > domain->boxhi_bound[1] ||
zlo < domain->boxlo_bound[2] || zhi > domain->boxhi_bound[2])
error->all("Deposition region extends outside simulation box");
}
// setup scaling
if (scaleflag && domain->lattice == NULL)
error->all("Use of fix deposit with undefined lattice");
double xscale,yscale,zscale;
if (scaleflag) {
xscale = domain->lattice->xlattice;
yscale = domain->lattice->ylattice;
zscale = domain->lattice->zlattice;
}
else xscale = yscale = zscale = 1.0;
// apply scaling to all input parameters with dist/vel units
if (domain->dimension == 2) {
lo *= yscale;
hi *= yscale;
rate *= yscale;
} else {
lo *= zscale;
hi *= zscale;
rate *= zscale;
}
deltasq *= xscale*xscale;
nearsq *= xscale*xscale;
vxlo *= xscale;
vxhi *= xscale;
vylo *= yscale;
vyhi *= yscale;
vzlo *= zscale;
vzhi *= zscale;
// random number generator, same for all procs
random = new RanPark(lmp,seed);
// set up reneighboring
force_reneighbor = 1;
next_reneighbor = update->ntimestep + 1;
nfirst = next_reneighbor;
ninserted = 0;
}
/* ---------------------------------------------------------------------- */
FixDeposit::~FixDeposit()
{
delete random;
}
/* ---------------------------------------------------------------------- */
int FixDeposit::setmask()
{
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ----------------------------------------------------------------------
perform particle insertion
------------------------------------------------------------------------- */
void FixDeposit::pre_exchange()
{
int i,j;
int flag,flagall;
double coord[3],lamda[3],delx,dely,delz,rsq;
double *newcoord;
// just return if should not be called on this timestep
if (next_reneighbor != update->ntimestep) return;
// compute current offset = bottom of insertion volume
double offset = 0.0;
if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate;
double *sublo,*subhi;
if (domain->triclinic == 0) {
sublo = domain->sublo;
subhi = domain->subhi;
} else {
sublo = domain->sublo_lamda;
subhi = domain->subhi_lamda;
}
// attempt an insertion until successful
int nfix = modify->nfix;
Fix **fix = modify->fix;
int success = 0;
int attempt = 0;
while (attempt < maxattempt) {
attempt++;
// choose random position for new atom within region
coord[0] = xlo + random->uniform() * (xhi-xlo);
coord[1] = ylo + random->uniform() * (yhi-ylo);
coord[2] = zlo + random->uniform() * (zhi-zlo);
while (domain->regions[iregion]->match(coord[0],coord[1],coord[2]) == 0) {
coord[0] = xlo + random->uniform() * (xhi-xlo);
coord[1] = ylo + random->uniform() * (yhi-ylo);
coord[2] = zlo + random->uniform() * (zhi-zlo);
}
// adjust vertical coord by offset
if (domain->dimension == 2) coord[1] += offset;
else coord[2] += offset;
// if global, reset vertical coord to be lo-hi above highest atom
// if local, reset vertical coord to be lo-hi above highest "nearby" atom
// local computation computes lateral distance between 2 particles w/ PBC
if (globalflag || localflag) {
int dim;
double max,maxall,delx,dely,delz,rsq;
if (domain->dimension == 2) {
dim = 1;
max = domain->boxlo[1];
} else {
dim = 2;
max = domain->boxlo[2];
}
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (localflag) {
delx = coord[0] - x[i][0];
dely = coord[1] - x[i][1];
delz = 0.0;
domain->minimum_image(delx,dely,delz);
if (domain->dimension == 2) rsq = delx*delx;
else rsq = delx*delx + dely*dely;
if (rsq > deltasq) continue;
}
if (x[i][dim] > max) max = x[i][dim];
}
MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world);
if (domain->dimension == 2)
coord[1] = maxall + lo + random->uniform()*(hi-lo);
else
coord[2] = maxall + lo + random->uniform()*(hi-lo);
}
// now have final coord
// if distance to any atom is less than near, try again
double **x = atom->x;
int nlocal = atom->nlocal;
flag = 0;
for (i = 0; i < nlocal; i++) {
delx = coord[0] - x[i][0];
dely = coord[1] - x[i][1];
delz = coord[2] - x[i][2];
domain->minimum_image(delx,dely,delz);
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < nearsq) flag = 1;
}
MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_MAX,world);
if (flagall) continue;
// insertion will proceed
// choose random velocity for new atom
double vxtmp = vxlo + random->uniform() * (vxhi-vxlo);
double vytmp = vylo + random->uniform() * (vyhi-vylo);
double vztmp = vzlo + random->uniform() * (vzhi-vzlo);
// check if new atom is in my sub-box or above it if I'm highest proc
// if so, add to my list via create_atom()
// initialize info about the atoms
// set group mask to "all" plus fix group
if (domain->triclinic) {
domain->x2lamda(coord,lamda);
newcoord = lamda;
} else newcoord = coord;
flag = 0;
if (newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1] &&
newcoord[2] >= sublo[2] && newcoord[2] < subhi[2]) flag = 1;
else if (domain->dimension == 3 && newcoord[2] >= domain->boxhi[2] &&
comm->myloc[2] == comm->procgrid[2]-1 &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0] &&
newcoord[1] >= sublo[1] && newcoord[1] < subhi[1]) flag = 1;
else if (domain->dimension == 2 && newcoord[1] >= domain->boxhi[1] &&
comm->myloc[1] == comm->procgrid[1]-1 &&
newcoord[0] >= sublo[0] && newcoord[0] < subhi[0]) flag = 1;
if (flag) {
atom->avec->create_atom(ntype,coord);
int m = atom->nlocal - 1;
atom->type[m] = ntype;
atom->mask[m] = 1 | groupbit;
atom->v[m][0] = vxtmp;
atom->v[m][1] = vytmp;
atom->v[m][2] = vztmp;
for (j = 0; j < nfix; j++)
if (fix[j]->create_attribute) fix[j]->set_arrays(m);
}
MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world);
break;
}
// warn if not successful b/c too many attempts or no proc owned particle
if (comm->me == 0)
if (success == 0)
error->warning("Particle deposition was unsuccessful");
// set tag # of new particle beyond all previous atoms
// reset global natoms
// if global map exists, reset it now instead of waiting for comm
// since deleting atoms messes up ghosts
if (success && atom->tag_enable) {
atom->tag_extend();
atom->natoms += 1;
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
}
// next timestep to insert
// next_reneighbor = 0 if done
if (success) ninserted++;
if (ninserted < ninsert) next_reneighbor += nfreq;
else next_reneighbor = 0;
}
/* ----------------------------------------------------------------------
parse optional parameters at end of input line
------------------------------------------------------------------------- */
void FixDeposit::options(int narg, char **arg)
{
if (narg < 0) error->all("Illegal fix indent command");
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all("Illegal fix deposit command");
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1) error->all("Fix deposit region ID does not exist");
iarg += 2;
} else if (strcmp(arg[iarg],"global") == 0) {
if (iarg+3 > narg) error->all("Illegal fix deposit command");
globalflag = 1;
localflag = 0;
lo = atof(arg[iarg+1]);
hi = atof(arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"local") == 0) {
if (iarg+4 > narg) error->all("Illegal fix deposit command");
localflag = 1;
globalflag = 0;
lo = atof(arg[iarg+1]);
hi = atof(arg[iarg+2]);
deltasq = atof(arg[iarg+3])*atof(arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"near") == 0) {
if (iarg+2 > narg) error->all("Illegal fix deposit command");
nearsq = atof(arg[iarg+1])*atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"attempt") == 0) {
if (iarg+2 > narg) error->all("Illegal fix deposit command");
maxattempt = atoi(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"rate") == 0) {
if (iarg+2 > narg) error->all("Illegal fix deposit command");
rateflag = 1;
rate = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"vx") == 0) {
if (iarg+3 > narg) error->all("Illegal fix deposit command");
vxlo = atof(arg[iarg+1]);
vxhi = atof(arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"vy") == 0) {
if (iarg+3 > narg) error->all("Illegal fix deposit command");
vylo = atof(arg[iarg+1]);
vyhi = atof(arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"vz") == 0) {
if (iarg+3 > narg) error->all("Illegal fix deposit command");
vzlo = atof(arg[iarg+1]);
vzhi = atof(arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"units") == 0) {
if (iarg+2 > narg) error->all("Illegal fix deposit command");
if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
else error->all("Illegal fix deposit command");
iarg += 2;
} else error->all("Illegal fix deposit command");
}
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixDeposit::write_restart(FILE *fp)
{
int n = 0;
double list[4];
list[n++] = random->state();
list[n++] = ninserted;
list[n++] = nfirst;
list[n++] = next_reneighbor;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
- fwrite(&list,sizeof(double),n,fp);
+ fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixDeposit::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
seed = static_cast<int> (list[n++]);
ninserted = static_cast<int> (list[n++]);
nfirst = static_cast<int> (list[n++]);
next_reneighbor = static_cast<int> (list[n++]);
random->reset(seed);
}
diff --git a/src/fix_nph.cpp b/src/fix_nph.cpp
index 1c91cf516..7692b8e4a 100644
--- a/src/fix_nph.cpp
+++ b/src/fix_nph.cpp
@@ -1,868 +1,868 @@
/* ----------------------------------------------------------------------
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 Stevens (SNL)
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "fix_nph.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "domain.h"
#include "modify.h"
#include "fix_deform.h"
#include "compute.h"
#include "kspace.h"
#include "update.h"
#include "respa.h"
#include "domain.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
enum{XYZ,XY,YZ,XZ,ANISO};
/* ---------------------------------------------------------------------- */
FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 4) error->all("Illegal fix nph command");
restart_global = 1;
box_change = 1;
time_integrate = 1;
scalar_flag = 1;
scalar_vector_freq = 1;
extscalar = 1;
double p_period[3];
if (strcmp(arg[3],"xyz") == 0) {
if (narg < 7) error->all("Illegal fix nph command");
press_couple = XYZ;
p_start[0] = p_start[1] = p_start[2] = atof(arg[4]);
p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[5]);
p_period[0] = p_period[1] = p_period[2] = atof(arg[6]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (domain->dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
}
} else {
if (strcmp(arg[3],"xy") == 0) press_couple = XY;
else if (strcmp(arg[3],"yz") == 0) press_couple = YZ;
else if (strcmp(arg[3],"xz") == 0) press_couple = XZ;
else if (strcmp(arg[3],"aniso") == 0) press_couple = ANISO;
else error->all("Illegal fix nph command");
if (narg < 11) error->all("Illegal fix nph command");
if (domain->dimension == 2 &&
(press_couple == XY || press_couple == YZ || press_couple == XZ))
error->all("Invalid fix nph command for a 2d simulation");
if (strcmp(arg[4],"NULL") == 0) {
p_start[0] = p_stop[0] = p_period[0] = 0.0;
p_flag[0] = 0;
} else {
p_start[0] = atof(arg[4]);
p_stop[0] = atof(arg[5]);
p_flag[0] = 1;
}
if (strcmp(arg[6],"NULL") == 0) {
p_start[1] = p_stop[1] = p_period[1] = 0.0;
p_flag[1] = 0;
} else {
p_start[1] = atof(arg[6]);
p_stop[1] = atof(arg[7]);
p_flag[1] = 1;
}
if (strcmp(arg[8],"NULL") == 0) {
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
} else {
if (domain->dimension == 2)
error->all("Invalid fix nph command for a 2d simulation");
p_start[2] = atof(arg[8]);
p_stop[2] = atof(arg[9]);
p_flag[2] = 1;
}
double period = atof(arg[10]);
if (p_flag[0]) p_period[0] = period;
if (p_flag[1]) p_period[1] = period;
if (p_flag[2]) p_period[2] = period;
}
// process extra keywords
drag = 0.0;
allremap = 1;
int iarg;
if (press_couple == XYZ) iarg = 7;
else iarg = 11;
while (iarg < narg) {
if (strcmp(arg[iarg],"drag") == 0) {
if (iarg+2 > narg) error->all("Illegal fix nph command");
drag = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dilate") == 0) {
if (iarg+2 > narg) error->all("Illegal fix nph command");
if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0;
else error->all("Illegal fix nph command");
iarg += 2;
} else error->all("Illegal fix nph command");
}
// error checks
if (press_couple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
error->all("Invalid fix nph command pressure settings");
if (press_couple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
error->all("Invalid fix nph command pressure settings");
if (press_couple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
error->all("Invalid fix nph command pressure settings");
if (press_couple == XY &&
(p_start[0] != p_start[1] || p_stop[0] != p_stop[1]))
error->all("Invalid fix nph command pressure settings");
if (press_couple == YZ &&
(p_start[1] != p_start[2] || p_stop[1] != p_stop[2]))
error->all("Invalid fix nph command pressure settings");
if (press_couple == XZ &&
(p_start[0] != p_start[2] || p_stop[0] != p_stop[2]))
error->all("Invalid fix nph command pressure settings");
if (p_flag[0] && domain->xperiodic == 0)
error->all("Cannot use fix nph on a non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0)
error->all("Cannot use fix nph on a non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0)
error->all("Cannot use fix nph on a non-periodic dimension");
// convert input periods to frequencies
if ((p_flag[0] && p_period[0] <= 0.0) ||
(p_flag[1] && p_period[1] <= 0.0) || (p_flag[2] && p_period[2] <= 0.0))
error->all("Fix nph periods must be > 0.0");
p_freq[0] = p_freq[1] = p_freq[2] = 0.0;
if (p_flag[0]) p_freq[0] = 1.0 / p_period[0];
if (p_flag[1]) p_freq[1] = 1.0 / p_period[1];
if (p_flag[2]) p_freq[2] = 1.0 / p_period[2];
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
tflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
pflag = 1;
// Nose/Hoover pressure init
omega[0] = omega[1] = omega[2] = 0.0;
omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0;
nrigid = 0;
rfix = NULL;
}
/* ---------------------------------------------------------------------- */
FixNPH::~FixNPH()
{
delete [] rfix;
// delete temperature and pressure if fix created them
if (tflag) modify->delete_compute(id_temp);
if (pflag) modify->delete_compute(id_press);
delete [] id_temp;
delete [] id_press;
}
/* ---------------------------------------------------------------------- */
int FixNPH::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= THERMO_ENERGY;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNPH::init()
{
if (domain->triclinic) error->all("Cannot use fix nph with triclinic box");
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) {
int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
(p_flag[2] && dimflag[2]))
error->all("Cannot use fix npt and fix deform on same dimension");
}
// set temperature and pressure ptrs
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Temperature ID for fix nph does not exist");
temperature = modify->compute[icompute];
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Pressure ID for fix nph does not exist");
pressure = modify->compute[icompute];
// set timesteps and frequencies
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dthalf = 0.5 * update->dt;
double freq = MAX(p_freq[0],p_freq[1]);
freq = MAX(freq,p_freq[2]);
drag_factor = 1.0 - (update->dt * freq * drag);
boltz = force->boltz;
nktv2p = force->nktv2p;
dimension = domain->dimension;
if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
else vol0 = domain->xprd * domain->yprd;
if (force->kspace) kspace_flag = 1;
else kspace_flag = 0;
if (strcmp(update->integrate_style,"respa") == 0) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
step_respa = ((Respa *) update->integrate)->step;
}
// detect if any rigid fixes exist so rigid bodies move when box is remapped
// rfix[] = indices to each fix rigid
delete [] rfix;
nrigid = 0;
rfix = NULL;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) nrigid++;
if (nrigid) {
rfix = new int[nrigid];
nrigid = 0;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
}
}
/* ----------------------------------------------------------------------
compute T,P before integrator starts
------------------------------------------------------------------------- */
void FixNPH::setup(int vflag)
{
// Nkt = initial value for piston mass and energy conservation
// cannot be done in init() b/c temperature cannot be called there
// is b/c Modify::init() inits computes after fixes due to dof dependence
// guesstimate a unit-dependent t_initial if actual T = 0.0
double t_initial = temperature->compute_scalar();
if (t_initial == 0.0) {
if (strcmp(update->unit_style,"lj") == 0) t_initial = 1.0;
else t_initial = 300.0;
}
nkt = atom->natoms * boltz * t_initial;
p_target[0] = p_start[0]; // used by compute_scalar()
p_target[1] = p_start[1];
p_target[2] = p_start[2];
if (press_couple == XYZ) {
double tmp = temperature->compute_scalar();
tmp = pressure->compute_scalar();
} else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
// trigger virial computation on next timestep
pressure->addstep(update->ntimestep+1);
}
/* ----------------------------------------------------------------------
1st half of Verlet update
------------------------------------------------------------------------- */
void FixNPH::initial_integrate(int vflag)
{
int i;
double dtfm;
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = nkt / volume * nktv2p;
for (i = 0; i < 3; i++) {
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
omega[i] += dtv*omega_dot[i];
factor[i] = exp(-dthalf*omega_dot[i]);
dilation[i] = exp(dthalf*omega_dot[i]);
}
// v update only for atoms in NPH group
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (rmass) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
} else {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
}
// remap simulation box and all owned atoms by 1/2 step
remap(0);
// x update by full step only for atoms in NPH group
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
// remap simulation box and all owned atoms by 1/2 step
// redo KSpace coeffs since volume has changed
remap(0);
if (kspace_flag) force->kspace->setup();
}
/* ----------------------------------------------------------------------
2nd half of Verlet update
------------------------------------------------------------------------- */
void FixNPH::final_integrate()
{
int i;
double dtfm;
// v update only for atoms in NPH group
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (rmass) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
}
}
} else {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
}
}
}
// compute new pressure
if (press_couple == XYZ) {
double tmp = temperature->compute_scalar();
tmp = pressure->compute_scalar();
} else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
// trigger virial computation on next timestep
pressure->addstep(update->ntimestep+1);
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = nkt / volume * nktv2p;
for (i = 0; i < 3; i++) {
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
}
}
/* ---------------------------------------------------------------------- */
void FixNPH::initial_integrate_respa(int vflag, int ilevel, int flag)
{
// if flag = 1, then is 2nd call at outermost level from rRESPA
// perform 2nd half of box remap on own + ghost atoms and return
// redo KSpace coeffs since volume has changed
if (flag == 1) {
remap(1);
if (kspace_flag) force->kspace->setup();
return;
}
// set timesteps by level
double dtfm;
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// atom quantities
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// outermost level - update omega_dot, apply to v, remap box
// all other levels - NVE update of v
// x,v updates only performed for atoms in NPH group
if (ilevel == nlevels_respa-1) {
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = nkt / volume * nktv2p;
for (int i = 0; i < 3; i++) {
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
omega[i] += dtv*omega_dot[i];
factor[i] = exp(-dthalf*omega_dot[i]);
dilation[i] = exp(dthalf*omega_dot[i]);
}
// v update only for atoms in NPH group
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
}
// remap simulation box and all owned atoms by 1/2 step
remap(0);
} else {
// v update only for atoms in NPH group
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
}
}
// innermost level - also update x only for atoms in NPH group
if (ilevel == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixNPH::final_integrate_respa(int ilevel)
{
double dtfm;
// set timesteps by level
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// outermost level - update omega_dot, apply to v via final_integrate()
// all other levels - NVE update of v
// v update only performed for atoms in NPH group
if (ilevel == nlevels_respa-1) final_integrate();
else {
// v update only for atoms in NPH group
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixNPH::couple()
{
double *tensor = pressure->vector;
if (press_couple == XYZ)
p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
else if (press_couple == XY) {
double ave = 0.5 * (tensor[0] + tensor[1]);
p_current[0] = p_current[1] = ave;
p_current[2] = tensor[2];
} else if (press_couple == YZ) {
double ave = 0.5 * (tensor[1] + tensor[2]);
p_current[1] = p_current[2] = ave;
p_current[0] = tensor[0];
} else if (press_couple == XZ) {
double ave = 0.5 * (tensor[0] + tensor[2]);
p_current[0] = p_current[2] = ave;
p_current[1] = tensor[1];
} else if (press_couple == ANISO) {
p_current[0] = tensor[0];
p_current[1] = tensor[1];
p_current[2] = tensor[2];
}
}
/* ----------------------------------------------------------------------
change box size
remap owned or owned+ghost atoms depending on flag
remap all atoms or fix group atoms depending on allremap flag
if rigid bodies exist, scale rigid body centers-of-mass
------------------------------------------------------------------------- */
void FixNPH::remap(int flag)
{
int i,n;
double oldlo,oldhi,ctr;
double **x = atom->x;
int *mask = atom->mask;
if (flag) n = atom->nlocal + atom->nghost;
else n = atom->nlocal;
// convert pertinent atoms and rigid bodies to lamda coords
if (allremap) domain->x2lamda(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->x2lamda(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(0);
// reset global and local box to new size/shape
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
oldlo = domain->boxlo[i];
oldhi = domain->boxhi[i];
ctr = 0.5 * (oldlo + oldhi);
domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr;
domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr;
}
}
domain->set_global_box();
domain->set_local_box();
// convert pertinent atoms and rigid bodies back to box coords
if (allremap) domain->lamda2x(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->lamda2x(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(1);
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixNPH::write_restart(FILE *fp)
{
int n = 0;
double list[6];
list[n++] = omega[0];
list[n++] = omega[1];
list[n++] = omega[2];
list[n++] = omega_dot[0];
list[n++] = omega_dot[1];
list[n++] = omega_dot[2];
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
- fwrite(&list,sizeof(double),n,fp);
+ fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixNPH::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
omega[0] = list[n++];
omega[1] = list[n++];
omega[2] = list[n++];
omega_dot[0] = list[n++];
omega_dot[1] = list[n++];
omega_dot[2] = list[n++];
}
/* ---------------------------------------------------------------------- */
int FixNPH::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (tflag) {
modify->delete_compute(id_temp);
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Could not find fix_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all("Fix_modify temperature ID does not compute temperature");
if (temperature->igroup != 0 && comm->me == 0)
error->warning("Temperature for NPH is not for group all");
// reset id_temp of pressure to new temperature ID
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Pressure ID for fix npt does not exist");
modify->compute[icompute]->reset_extra_compute_fix(id_temp);
return 2;
} else if (strcmp(arg[0],"press") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (pflag) {
modify->delete_compute(id_press);
pflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
int icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Could not find fix_modify pressure ID");
pressure = modify->compute[icompute];
if (pressure->pressflag == 0)
error->all("Fix_modify pressure ID does not compute pressure");
return 2;
}
return 0;
}
/* ---------------------------------------------------------------------- */
double FixNPH::compute_scalar()
{
double volume;
if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
else volume = domain->xprd * domain->yprd;
int pdim = p_flag[0] + p_flag[1] + p_flag[2];
double energy = 0.0;
for (int i = 0; i < 3; i++)
if (p_freq[i] > 0.0)
energy += 0.5*nkt*omega_dot[i]*omega_dot[i] /
(p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p);
return energy;
}
diff --git a/src/fix_npt.cpp b/src/fix_npt.cpp
index 3328625a4..a8880e470 100644
--- a/src/fix_npt.cpp
+++ b/src/fix_npt.cpp
@@ -1,1022 +1,1022 @@
/* ----------------------------------------------------------------------
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 Stevens (SNL)
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "fix_npt.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "modify.h"
#include "fix_deform.h"
#include "compute.h"
#include "kspace.h"
#include "update.h"
#include "respa.h"
#include "domain.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
enum{NOBIAS,BIAS};
enum{XYZ,XY,YZ,XZ,ANISO};
/* ---------------------------------------------------------------------- */
FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 7) error->all("Illegal fix npt command");
restart_global = 1;
box_change = 1;
time_integrate = 1;
scalar_flag = 1;
scalar_vector_freq = 1;
extscalar = 1;
t_start = atof(arg[3]);
t_stop = atof(arg[4]);
double t_period = atof(arg[5]);
if (t_start < 0.0 || t_stop <= 0.0)
error->all("Target T for fix npt cannot be 0.0");
double p_period[3];
if (strcmp(arg[6],"xyz") == 0) {
if (narg < 10) error->all("Illegal fix npt command");
press_couple = XYZ;
p_start[0] = p_start[1] = p_start[2] = atof(arg[7]);
p_stop[0] = p_stop[1] = p_stop[2] = atof(arg[8]);
p_period[0] = p_period[1] = p_period[2] = atof(arg[9]);
p_flag[0] = p_flag[1] = p_flag[2] = 1;
if (domain->dimension == 2) {
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
}
} else {
if (strcmp(arg[6],"xy") == 0) press_couple = XY;
else if (strcmp(arg[6],"yz") == 0) press_couple = YZ;
else if (strcmp(arg[6],"xz") == 0) press_couple = XZ;
else if (strcmp(arg[6],"aniso") == 0) press_couple = ANISO;
else error->all("Illegal fix npt command");
if (narg < 14) error->all("Illegal fix npt command");
if (domain->dimension == 2 &&
(press_couple == XY || press_couple == YZ || press_couple == XZ))
error->all("Invalid fix npt command for a 2d simulation");
if (strcmp(arg[7],"NULL") == 0) {
p_start[0] = p_stop[0] = p_period[0] = 0.0;
p_flag[0] = 0;
} else {
p_start[0] = atof(arg[7]);
p_stop[0] = atof(arg[8]);
p_flag[0] = 1;
}
if (strcmp(arg[9],"NULL") == 0) {
p_start[1] = p_stop[1] = p_period[1] = 0.0;
p_flag[1] = 0;
} else {
p_start[1] = atof(arg[9]);
p_stop[1] = atof(arg[10]);
p_flag[1] = 1;
}
if (strcmp(arg[11],"NULL") == 0) {
p_start[2] = p_stop[2] = p_period[2] = 0.0;
p_flag[2] = 0;
} else {
if (domain->dimension == 2)
error->all("Invalid fix npt command for a 2d simulation");
p_start[2] = atof(arg[11]);
p_stop[2] = atof(arg[12]);
p_flag[2] = 1;
}
double period = atof(arg[13]);
if (p_flag[0]) p_period[0] = period;
if (p_flag[1]) p_period[1] = period;
if (p_flag[2]) p_period[2] = period;
}
// process extra keywords
drag = 0.0;
allremap = 1;
int iarg;
if (press_couple == XYZ) iarg = 10;
else iarg = 14;
while (iarg < narg) {
if (strcmp(arg[iarg],"drag") == 0) {
if (iarg+2 > narg) error->all("Illegal fix npt command");
drag = atof(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dilate") == 0) {
if (iarg+2 > narg) error->all("Illegal fix npt command");
if (strcmp(arg[iarg+1],"all") == 0) allremap = 1;
else if (strcmp(arg[iarg+1],"partial") == 0) allremap = 0;
else error->all("Illegal fix npt command");
iarg += 2;
} else error->all("Illegal fix npt command");
}
// error checks
if (press_couple == XY && (p_flag[0] == 0 || p_flag[1] == 0))
error->all("Invalid fix npt command pressure settings");
if (press_couple == YZ && (p_flag[1] == 0 || p_flag[2] == 0))
error->all("Invalid fix npt command pressure settings");
if (press_couple == XZ && (p_flag[0] == 0 || p_flag[2] == 0))
error->all("Invalid fix npt command pressure settings");
if (press_couple == XY &&
(p_start[0] != p_start[1] || p_stop[0] != p_stop[1]))
error->all("Invalid fix npt command pressure settings");
if (press_couple == YZ &&
(p_start[1] != p_start[2] || p_stop[1] != p_stop[2]))
error->all("Invalid fix npt command pressure settings");
if (press_couple == XZ &&
(p_start[0] != p_start[2] || p_stop[0] != p_stop[2]))
error->all("Invalid fix npt command pressure settings");
if (p_flag[0] && domain->xperiodic == 0)
error->all("Cannot use fix npt on a non-periodic dimension");
if (p_flag[1] && domain->yperiodic == 0)
error->all("Cannot use fix npt on a non-periodic dimension");
if (p_flag[2] && domain->zperiodic == 0)
error->all("Cannot use fix npt on a non-periodic dimension");
// convert input periods to frequencies
if (t_period <= 0.0 || (p_flag[0] && p_period[0] <= 0.0) ||
(p_flag[1] && p_period[1] <= 0.0) || (p_flag[2] && p_period[2] <= 0.0))
error->all("Fix npt periods must be > 0.0");
t_freq = 1.0 / t_period;
p_freq[0] = p_freq[1] = p_freq[2] = 0.0;
if (p_flag[0]) p_freq[0] = 1.0 / p_period[0];
if (p_flag[1]) p_freq[1] = 1.0 / p_period[1];
if (p_flag[2]) p_freq[2] = 1.0 / p_period[2];
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
if (strcmp(style,"npt") == 0) newarg[2] = (char *) "temp";
else if (strcmp(style,"npt/asphere") == 0)
newarg[2] = (char *) "temp/asphere";
else if (strcmp(style,"npt/sphere") == 0)
newarg[2] = (char *) "temp/sphere";
modify->add_compute(3,newarg);
delete [] newarg;
tflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
pflag = 1;
// Nose/Hoover temp and pressure init
eta = eta_dot = 0.0;
omega[0] = omega[1] = omega[2] = 0.0;
omega_dot[0] = omega_dot[1] = omega_dot[2] = 0.0;
nrigid = 0;
rfix = NULL;
}
/* ---------------------------------------------------------------------- */
FixNPT::~FixNPT()
{
delete [] rfix;
// delete temperature and pressure if fix created them
if (tflag) modify->delete_compute(id_temp);
if (pflag) modify->delete_compute(id_press);
delete [] id_temp;
delete [] id_press;
}
/* ---------------------------------------------------------------------- */
int FixNPT::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= THERMO_ENERGY;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNPT::init()
{
if (domain->triclinic) error->all("Cannot use fix npt with triclinic box");
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) {
int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
(p_flag[2] && dimflag[2]))
error->all("Cannot use fix npt and fix deform on same dimension");
}
// set temperature and pressure ptrs
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Temperature ID for fix npt does not exist");
temperature = modify->compute[icompute];
if (temperature->tempbias) which = BIAS;
else which = NOBIAS;
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Pressure ID for fix npt does not exist");
pressure = modify->compute[icompute];
// set timesteps and frequencies
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dthalf = 0.5 * update->dt;
double freq = MAX(p_freq[0],p_freq[1]);
freq = MAX(freq,p_freq[2]);
drag_factor = 1.0 - (update->dt * freq * drag);
boltz = force->boltz;
nktv2p = force->nktv2p;
dimension = domain->dimension;
if (dimension == 3) vol0 = domain->xprd * domain->yprd * domain->zprd;
else vol0 = domain->xprd * domain->yprd;
if (force->kspace) kspace_flag = 1;
else kspace_flag = 0;
if (strcmp(update->integrate_style,"respa") == 0) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
step_respa = ((Respa *) update->integrate)->step;
}
// detect if any rigid fixes exist so rigid bodies move when box is remapped
// rfix[] = indices to each fix rigid
delete [] rfix;
nrigid = 0;
rfix = NULL;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) nrigid++;
if (nrigid) {
rfix = new int[nrigid];
nrigid = 0;
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i;
}
}
/* ----------------------------------------------------------------------
compute T,P before integrator starts
------------------------------------------------------------------------- */
void FixNPT::setup(int vflag)
{
t_target = t_start; // used by compute_scalar()
p_target[0] = p_start[0];
p_target[1] = p_start[1];
p_target[2] = p_start[2];
t_current = temperature->compute_scalar();
if (press_couple == XYZ) {
double tmp = pressure->compute_scalar();
} else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
// trigger virial computation on next timestep
pressure->addstep(update->ntimestep+1);
}
/* ----------------------------------------------------------------------
1st half of Verlet update
------------------------------------------------------------------------- */
void FixNPT::initial_integrate(int vflag)
{
int i;
double dtfm;
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
// update eta_dot
t_target = t_start + delta * (t_stop-t_start);
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
eta += dtv*eta_dot;
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = atom->natoms*boltz*t_target / volume * nktv2p;
for (i = 0; i < 3; i++) {
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
omega[i] += dtv*omega_dot[i];
factor[i] = exp(-dthalf*(eta_dot+omega_dot[i]));
dilation[i] = exp(dthalf*omega_dot[i]);
}
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) {
if (which == NOBIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
} else {
if (which == NOBIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
}
// remap simulation box and all owned atoms by 1/2 step
remap(0);
// x update by full step only for atoms in group
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
// remap simulation box and all owned atoms by 1/2 step
// redo KSpace coeffs since volume has changed
remap(0);
if (kspace_flag) force->kspace->setup();
}
/* ----------------------------------------------------------------------
2nd half of Verlet update
------------------------------------------------------------------------- */
void FixNPT::final_integrate()
{
int i;
double dtfm;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) {
if (which == NOBIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
}
}
} else if (which == BIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / rmass[i];
v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
temperature->restore_bias(i,v[i]);
}
}
}
} else {
if (which == NOBIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
}
}
} else if (which == BIAS) {
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = (v[i][0] + dtfm*f[i][0]) * factor[0];
v[i][1] = (v[i][1] + dtfm*f[i][1]) * factor[1];
v[i][2] = (v[i][2] + dtfm*f[i][2]) * factor[2];
temperature->restore_bias(i,v[i]);
}
}
}
}
// compute new T,P
t_current = temperature->compute_scalar();
if (press_couple == XYZ) {
double tmp = pressure->compute_scalar();
} else {
temperature->compute_vector();
pressure->compute_vector();
}
couple();
// trigger virial computation on next timestep
pressure->addstep(update->ntimestep+1);
// update eta_dot
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = atom->natoms*boltz*t_target / volume * nktv2p;
for (i = 0; i < 3; i++) {
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
}
}
/* ---------------------------------------------------------------------- */
void FixNPT::initial_integrate_respa(int vflag, int ilevel, int flag)
{
// if flag = 1, then is 2nd call at outermost level from rRESPA
// perform 2nd half of box remap on own + ghost atoms and return
// redo KSpace coeffs since volume has changed
if (flag == 1) {
remap(1);
if (kspace_flag) force->kspace->setup();
return;
}
// set timesteps by level
double dtfm;
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// atom quantities
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// outermost level - update eta_dot and omega_dot, apply to v, remap box
// all other levels - NVE update of v
// x,v updates only performed for atoms in group
if (ilevel == nlevels_respa-1) {
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
// update eta_dot
t_target = t_start + delta * (t_stop-t_start);
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
eta += dtv*eta_dot;
// update omega_dot
// for non-varying dims, p_freq is 0.0, so omega_dot doesn't change
double f_omega,volume;
if (dimension == 3) volume = domain->xprd*domain->yprd*domain->zprd;
else volume = domain->xprd*domain->yprd;
double denskt = atom->natoms*boltz*t_target / volume * nktv2p;
for (int i = 0; i < 3; i++) {
p_target[i] = p_start[i] + delta * (p_stop[i]-p_start[i]);
f_omega = p_freq[i]*p_freq[i] * (p_current[i]-p_target[i])/denskt;
omega_dot[i] += f_omega*dthalf;
omega_dot[i] *= drag_factor;
omega[i] += dtv*omega_dot[i];
factor[i] = exp(-dthalf*(eta_dot+omega_dot[i]));
dilation[i] = exp(dthalf*omega_dot[i]);
}
// v update only for atoms in group
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor[0] + dtfm*f[i][0];
v[i][1] = v[i][1]*factor[1] + dtfm*f[i][1];
v[i][2] = v[i][2]*factor[2] + dtfm*f[i][2];
}
}
}
// remap simulation box and all owned atoms by 1/2 step
remap(0);
} else {
// v update only for atoms in group
if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / rmass[i];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
}
}
// innermost level - also update x only for atoms in group
if (ilevel == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixNPT::final_integrate_respa(int ilevel)
{
double dtfm;
// set timesteps by level
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// outermost level - update eta_dot and omega_dot,
// apply to v via final_integrate()
// all other levels - NVE update of v
// v update only performed for atoms in group
if (ilevel == nlevels_respa-1) final_integrate();
else {
// v update only for atoms in group
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / rmass[i];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
}
}
}
/* ---------------------------------------------------------------------- */
void FixNPT::couple()
{
double *tensor = pressure->vector;
if (press_couple == XYZ)
p_current[0] = p_current[1] = p_current[2] = pressure->scalar;
else if (press_couple == XY) {
double ave = 0.5 * (tensor[0] + tensor[1]);
p_current[0] = p_current[1] = ave;
p_current[2] = tensor[2];
} else if (press_couple == YZ) {
double ave = 0.5 * (tensor[1] + tensor[2]);
p_current[1] = p_current[2] = ave;
p_current[0] = tensor[0];
} else if (press_couple == XZ) {
double ave = 0.5 * (tensor[0] + tensor[2]);
p_current[0] = p_current[2] = ave;
p_current[1] = tensor[1];
} else if (press_couple == ANISO) {
p_current[0] = tensor[0];
p_current[1] = tensor[1];
p_current[2] = tensor[2];
}
}
/* ----------------------------------------------------------------------
change box size
remap owned or owned+ghost atoms depending on flag
remap all atoms or fix group atoms depending on allremap flag
if rigid bodies exist, scale rigid body centers-of-mass
------------------------------------------------------------------------- */
void FixNPT::remap(int flag)
{
int i,n;
double oldlo,oldhi,ctr;
double **x = atom->x;
int *mask = atom->mask;
if (flag) n = atom->nlocal + atom->nghost;
else n = atom->nlocal;
// convert pertinent atoms and rigid bodies to lamda coords
if (allremap) domain->x2lamda(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->x2lamda(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(0);
// reset global and local box to new size/shape
for (i = 0; i < 3; i++) {
if (p_flag[i]) {
oldlo = domain->boxlo[i];
oldhi = domain->boxhi[i];
ctr = 0.5 * (oldlo + oldhi);
domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr;
domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr;
}
}
domain->set_global_box();
domain->set_local_box();
// convert pertinent atoms and rigid bodies back to box coords
if (allremap) domain->lamda2x(n);
else {
for (i = 0; i < n; i++)
if (mask[i] & groupbit)
domain->lamda2x(x[i],x[i]);
}
if (nrigid)
for (i = 0; i < nrigid; i++)
modify->fix[rfix[i]]->deform(1);
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixNPT::write_restart(FILE *fp)
{
int n = 0;
double list[8];
list[n++] = eta;
list[n++] = eta_dot;
list[n++] = omega[0];
list[n++] = omega[1];
list[n++] = omega[2];
list[n++] = omega_dot[0];
list[n++] = omega_dot[1];
list[n++] = omega_dot[2];
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
- fwrite(&list,sizeof(double),n,fp);
+ fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixNPT::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
eta = list[n++];
eta_dot = list[n++];
omega[0] = list[n++];
omega[1] = list[n++];
omega[2] = list[n++];
omega_dot[0] = list[n++];
omega_dot[1] = list[n++];
omega_dot[2] = list[n++];
}
/* ---------------------------------------------------------------------- */
int FixNPT::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (tflag) {
modify->delete_compute(id_temp);
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all("Could not find fix_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all("Fix_modify temperature ID does not compute temperature");
if (temperature->igroup != 0 && comm->me == 0)
error->warning("Temperature for fix modify is not for group all");
// reset id_temp of pressure to new temperature ID
icompute = modify->find_compute(id_press);
if (icompute < 0) error->all("Pressure ID for fix modify does not exist");
modify->compute[icompute]->reset_extra_compute_fix(id_temp);
return 2;
} else if (strcmp(arg[0],"press") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (pflag) {
modify->delete_compute(id_press);
pflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all("Could not find fix_modify pressure ID");
pressure = modify->compute[icompute];
if (pressure->pressflag == 0)
error->all("Fix_modify pressure ID does not compute pressure");
return 2;
}
return 0;
}
/* ---------------------------------------------------------------------- */
double FixNPT::compute_scalar()
{
double ke = temperature->dof * boltz * t_target;
double keplus = atom->natoms * boltz * t_target;
double volume;
if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
else volume = domain->xprd * domain->yprd;
int pdim = p_flag[0] + p_flag[1] + p_flag[2];
double energy = ke * (eta + 0.5*eta_dot*eta_dot/(t_freq*t_freq));
for (int i = 0; i < 3; i++)
if (p_freq[i] > 0.0)
energy += 0.5*keplus*omega_dot[i]*omega_dot[i] /
(p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p);
return energy;
}
/* ---------------------------------------------------------------------- */
void FixNPT::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dthalf = 0.5 * update->dt;
double freq = MAX(p_freq[0],p_freq[1]);
freq = MAX(freq,p_freq[2]);
drag_factor = 1.0 - (update->dt * freq * drag);
}
diff --git a/src/fix_nvt.cpp b/src/fix_nvt.cpp
index fb41cbbdd..44531aad1 100644
--- a/src/fix_nvt.cpp
+++ b/src/fix_nvt.cpp
@@ -1,689 +1,689 @@
/* ----------------------------------------------------------------------
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: Mark Stevens (SNL)
Andy Ballard (U Maryland) for Nose/Hoover chains
------------------------------------------------------------------------- */
#include "string.h"
#include "stdlib.h"
#include "math.h"
#include "fix_nvt.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "group.h"
#include "update.h"
#include "respa.h"
#include "modify.h"
#include "compute.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{NOBIAS,BIAS};
/* ---------------------------------------------------------------------- */
FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 6) error->all("Illegal fix nvt command");
restart_global = 1;
time_integrate = 1;
scalar_flag = 1;
scalar_vector_freq = 1;
extscalar = 1;
t_start = atof(arg[3]);
t_stop = atof(arg[4]);
double t_period = atof(arg[5]);
drag = 0.0;
chain = 1;
int iarg = 6;
while (iarg < narg) {
if (strcmp(arg[iarg],"drag") == 0) {
if (iarg+2 > narg) error->all("Illegal fix nvt command");
drag = atof(arg[iarg+1]);
if (drag < 0.0) error->all("Illegal fix nvt command");
iarg += 2;
} else if (strcmp(arg[iarg],"chain") == 0) {
if (iarg+2 > narg) error->all("Illegal fix nvt command");
if (strcmp(arg[iarg+1],"yes") == 0) chain = 1;
else if (strcmp(arg[iarg+1],"no") == 0) chain = 0;
else error->all("Illegal fix nvt command");
iarg += 2;
} else error->all("Illegal fix nvt command");
}
// error checks
// convert input period to frequency
if (t_start < 0.0 || t_stop <= 0.0)
error->all("Target T for fix nvt cannot be 0.0");
if (t_period <= 0.0) error->all("Fix nvt period must be > 0.0");
t_freq = 1.0 / t_period;
// create a new compute temp style
// id = fix-ID + temp, compute group = fix group
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
if (strcmp(style,"nvt") == 0)
newarg[2] = (char *) "temp";
else if (strcmp(style,"nvt/asphere") == 0)
newarg[2] = (char *) "temp/asphere";
else if (strcmp(style,"nvt/sllod") == 0)
newarg[2] = (char *) "temp/deform";
else if (strcmp(style,"nvt/sphere") == 0)
newarg[2] = (char *) "temp/sphere";
modify->add_compute(3,newarg);
delete [] newarg;
tflag = 1;
// Nose/Hoover temp init
eta = eta_dot = 0.0;
eta2 = eta2_dot = 0.0;
}
/* ---------------------------------------------------------------------- */
FixNVT::~FixNVT()
{
// delete temperature if fix created it
if (tflag) modify->delete_compute(id_temp);
delete [] id_temp;
}
/* ---------------------------------------------------------------------- */
int FixNVT::setmask()
{
int mask = 0;
mask |= INITIAL_INTEGRATE;
mask |= FINAL_INTEGRATE;
mask |= THERMO_ENERGY;
mask |= INITIAL_INTEGRATE_RESPA;
mask |= FINAL_INTEGRATE_RESPA;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixNVT::init()
{
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Temperature ID for fix nvt does not exist");
temperature = modify->compute[icompute];
if (temperature->tempbias) which = BIAS;
else which = NOBIAS;
// set timesteps and frequencies
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
drag_factor = 1.0 - (update->dt * t_freq * drag);
if (strcmp(update->integrate_style,"respa") == 0) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
step_respa = ((Respa *) update->integrate)->step;
}
}
/* ---------------------------------------------------------------------- */
void FixNVT::setup(int vflag)
{
t_target = t_start; // used by compute_scalar()
t_current = temperature->compute_scalar();
}
/* ---------------------------------------------------------------------- */
void FixNVT::initial_integrate(int vflag)
{
double dtfm;
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
t_target = t_start + delta * (t_stop-t_start);
// update eta, eta_dot, eta2, eta2_dot
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (chain) {
eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4;
factor2 = exp(-dt8*eta2_dot);
eta_dot *= factor2;
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dt4;
eta_dot *= drag_factor;
eta_dot *= factor2;
eta += dthalf*eta_dot;
eta2 += dthalf*eta2_dot;
} else {
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
eta += dtv*eta_dot;
}
factor = exp(-dthalf*eta_dot);
// update v and x of only atoms in group
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
}
if (chain) {
eta_dot *= factor2;
f_eta = t_freq*t_freq * (factor * factor * t_current/t_target - 1.0);
eta_dot += f_eta * dt4;
eta_dot *= factor2;
eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4;
}
}
/* ---------------------------------------------------------------------- */
void FixNVT::final_integrate()
{
double dtfm;
// update v of only atoms in group
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (chain) factor = 1.0;
if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i] * factor;
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / rmass[i] * factor;
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]] * factor;
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]] * factor;
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
}
// compute current T
t_current = temperature->compute_scalar();
// update eta, eta_dot, eta2, eta2_dot
// chains require additional velocity update and recompute of current T
if (chain) {
eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4;
factor2 = exp(-dt8*eta2_dot);
eta_dot *= factor2;
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dt4;
eta_dot *= drag_factor;
eta_dot *= factor2;
eta += dthalf*eta_dot;
eta2 += dthalf*eta2_dot;
factor = exp(-dthalf*eta_dot);
if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0] * factor;
v[i][1] = v[i][1] * factor;
v[i][2] = v[i][2] * factor;
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0] * factor;
v[i][1] = v[i][1] * factor;
v[i][2] = v[i][2] * factor;
temperature->restore_bias(i,v[i]);
}
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0] * factor;
v[i][1] = v[i][1] * factor;
v[i][2] = v[i][2] * factor;
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0] * factor;
v[i][1] = v[i][1] * factor;
v[i][2] = v[i][2] * factor;
temperature->restore_bias(i,v[i]);
}
}
}
}
t_current = temperature->compute_scalar();
eta_dot *= factor2;
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta * dt4;
eta_dot *= factor2;
eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4;
} else {
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
}
}
/* ---------------------------------------------------------------------- */
void FixNVT::initial_integrate_respa(int vflag, int ilevel, int flag)
{
if (flag) return; // only used by NPT,NPH
// set timesteps by level
double dtfm;
dtv = step_respa[ilevel];
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// atom quantities
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// outermost level - update eta_dot and apply to v with factor
// all other levels - NVE update of v (factor = 1)
// innermost level - also update x
if (ilevel == nlevels_respa-1) {
double delta = update->ntimestep - update->beginstep;
delta /= update->endstep - update->beginstep;
t_target = t_start + delta * (t_stop-t_start);
if (chain) {
eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4;
factor2 = exp(-dt8*eta2_dot);
eta_dot *= factor2;
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dt4;
eta_dot *= drag_factor;
eta_dot *= factor2;
eta += dthalf*eta_dot;
eta2 += dthalf*eta2_dot;
} else {
f_eta = t_freq*t_freq * (t_current/t_target - 1.0);
eta_dot += f_eta*dthalf;
eta_dot *= drag_factor;
eta += dtv*eta_dot;
}
factor = exp(-dthalf*eta_dot);
} else factor = 1.0;
if (rmass) {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / rmass[i];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
} else {
if (which == NOBIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
}
}
} else if (which == BIAS) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
temperature->remove_bias(i,v[i]);
dtfm = dtf / mass[type[i]];
v[i][0] = v[i][0]*factor + dtfm*f[i][0];
v[i][1] = v[i][1]*factor + dtfm*f[i][1];
v[i][2] = v[i][2]*factor + dtfm*f[i][2];
temperature->restore_bias(i,v[i]);
}
}
}
}
if (ilevel == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
}
}
}
if (ilevel == nlevels_respa-1 && chain) {
eta_dot *= factor2;
f_eta = t_freq*t_freq * (factor * factor * t_current/t_target - 1.0);
eta_dot += f_eta * dt4;
eta_dot *= factor2;
eta2_dot += (temperature->dof * eta_dot*eta_dot - t_freq*t_freq) * dt4;
}
}
/* ---------------------------------------------------------------------- */
void FixNVT::final_integrate_respa(int ilevel)
{
// set timesteps by level
double dtfm;
dtf = 0.5 * step_respa[ilevel] * force->ftm2v;
dthalf = 0.5 * step_respa[ilevel];
// outermost level - update eta_dot and apply to v via final_integrate()
// all other levels - NVE update of v
if (ilevel == nlevels_respa-1) final_integrate();
else {
double **v = atom->v;
double **f = atom->f;
double *rmass = atom->rmass;
double *mass = atom->mass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
if (rmass) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / mass[type[i]];
v[i][0] += dtfm*f[i][0];
v[i][1] += dtfm*f[i][1];
v[i][2] += dtfm*f[i][2];
}
}
}
}
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixNVT::write_restart(FILE *fp)
{
int n = 0;
double list[4];
list[n++] = eta;
list[n++] = eta_dot;
list[n++] = eta2;
list[n++] = eta2_dot;
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
- fwrite(&list,sizeof(double),n,fp);
+ fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixNVT::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
eta = list[n++];
eta_dot = list[n++];
eta2 = list[n++];
eta2_dot = list[n++];
}
/* ---------------------------------------------------------------------- */
int FixNVT::modify_param(int narg, char **arg)
{
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all("Illegal fix_modify command");
if (tflag) {
modify->delete_compute(id_temp);
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all("Could not find fix_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all("Fix_modify temperature ID does not compute temperature");
if (temperature->igroup != igroup && comm->me == 0)
error->warning("Group for fix_modify temp != fix group");
return 2;
}
return 0;
}
/* ---------------------------------------------------------------------- */
void FixNVT::reset_target(double t_new)
{
t_start = t_stop = t_new;
}
/* ---------------------------------------------------------------------- */
void FixNVT::reset_dt()
{
dtv = update->dt;
dtf = 0.5 * update->dt * force->ftm2v;
dthalf = 0.5 * update->dt;
dt4 = 0.25 * update->dt;
dt8 = 0.125 * update->dt;
drag_factor = 1.0 - (update->dt * t_freq * drag);
}
/* ---------------------------------------------------------------------- */
double FixNVT::compute_scalar()
{
double ke = temperature->dof * force->boltz * t_target;
double energy = ke * (eta + 0.5*eta_dot*eta_dot/(t_freq*t_freq));
return energy;
}
diff --git a/src/fix_ttm.cpp b/src/fix_ttm.cpp
index 3c5c181c7..29d07f41c 100644
--- a/src/fix_ttm.cpp
+++ b/src/fix_ttm.cpp
@@ -1,695 +1,695 @@
/* ----------------------------------------------------------------------
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: Paul Crozier (SNL)
Carolyn Phillips (University of Michigan)
------------------------------------------------------------------------- */
#include "mpi.h"
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "fix_ttm.h"
#include "atom.h"
#include "force.h"
#include "update.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "comm.h"
#include "random_mars.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define MIN(A,B) ((A) < (B)) ? (A) : (B)
#define MAX(A,B) ((A) > (B)) ? (A) : (B)
#define MAXLINE 1024
/* ---------------------------------------------------------------------- */
FixTTM::FixTTM(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 15) error->all("Illegal fix ttm command");
vector_flag = 1;
size_vector = 2;
scalar_vector_freq = 1;
extvector = 1;
nevery = 1;
restart_peratom = 1;
restart_global = 1;
seed = atoi(arg[3]);
electronic_specific_heat = atof(arg[4]);
electronic_density = atof(arg[5]);
electronic_thermal_conductivity = atof(arg[6]);
gamma_p = atof(arg[7]);
gamma_s = atof(arg[8]);
v_0 = atof(arg[9]);
nxnodes = atoi(arg[10]);
nynodes = atoi(arg[11]);
nznodes = atoi(arg[12]);
fpr = fopen(arg[13],"r");
if (fpr == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",arg[13]);
error->one(str);
}
nfileevery = atoi(arg[14]);
if (nfileevery) {
if (narg != 16) error->all("Illegal fix ttm command");
MPI_Comm_rank(world,&me);
if (me == 0) {
fp = fopen(arg[15],"w");
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open fix ttm file %s",arg[15]);
error->one(str);
}
}
}
// error check
if (seed <= 0) error->all("Invalid random number seed in fix ttm command");
if (electronic_specific_heat <= 0.0)
error->all("Fix ttm electronic_specific_heat must be > 0.0");
if (electronic_density <= 0.0)
error->all("Fix ttm electronic_density must be > 0.0");
if (electronic_thermal_conductivity < 0.0)
error->all("Fix ttm electronic_thermal_conductivity must be >= 0.0");
if (gamma_p <= 0.0) error->all("Fix ttm gamma_p must be > 0.0");
if (gamma_s < 0.0) error->all("Fix ttm gamma_s must be >= 0.0");
if (v_0 < 0.0) error->all("Fix ttm v_0 must be >= 0.0");
if (nxnodes <= 0 || nynodes <= 0 || nznodes <= 0)
error->all("Fix ttm number of nodes must be > 0");
v_0_sq = v_0*v_0;
// initialize Marsaglia RNG with processor-unique seed
random = new RanMars(lmp,seed + comm->me);
// allocate per-type arrays for force prefactors
gfactor1 = new double[atom->ntypes+1];
gfactor2 = new double[atom->ntypes+1];
// allocate 3d grid variables
total_nnodes = nxnodes*nynodes*nznodes;
nsum = memory->create_3d_int_array(nxnodes,nynodes,nznodes,"ttm:nsum");
nsum_all = memory->create_3d_int_array(nxnodes,nynodes,nznodes,
"ttm:nsum_all");
T_initial_set = memory->create_3d_int_array(nxnodes,nynodes,nznodes,
"ttm:T_initial_set");
sum_vsq = memory->create_3d_double_array(nxnodes,nynodes,nznodes,
"ttm:sum_vsq");
sum_mass_vsq =
memory->create_3d_double_array(nxnodes,nynodes,nznodes,
"ttm:sum_mass_vsq");
sum_vsq_all = memory->create_3d_double_array(nxnodes,nynodes,nznodes,
"ttm:sum_vsq_all");
sum_mass_vsq_all =
memory->create_3d_double_array(nxnodes,nynodes,nznodes,
"ttm:sum_mass_vsq_all");
T_electron_old =
memory->create_3d_double_array(nxnodes,nynodes,nznodes,
"ttm:T_electron_old");
T_electron =
memory->create_3d_double_array(nxnodes,nynodes,nznodes,"ttm:T_electron");
net_energy_transfer =
memory->create_3d_double_array(nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer");
net_energy_transfer_all =
memory->create_3d_double_array(nxnodes,nynodes,nznodes,
"TTM:net_energy_transfer_all");
flangevin = NULL;
grow_arrays(atom->nmax);
// zero out the flangevin array
for (int i = 0; i < atom->nmax; i++) {
flangevin[i][0] = 0;
flangevin[i][1] = 0;
flangevin[i][2] = 0;
}
atom->add_callback(0);
atom->add_callback(1);
// set initial electron temperatures from user input file
if (me == 0) read_initial_electron_temperatures();
MPI_Bcast(&T_electron[0][0][0],total_nnodes,MPI_DOUBLE,0,world);
}
/* ---------------------------------------------------------------------- */
FixTTM::~FixTTM()
{
if (nfileevery && me == 0) fclose(fp);
delete random;
delete [] gfactor1;
delete [] gfactor2;
memory->destroy_3d_int_array(nsum);
memory->destroy_3d_int_array(nsum_all);
memory->destroy_3d_int_array(T_initial_set);
memory->destroy_3d_double_array(sum_vsq);
memory->destroy_3d_double_array(sum_mass_vsq);
memory->destroy_3d_double_array(sum_vsq_all);
memory->destroy_3d_double_array(sum_mass_vsq_all);
memory->destroy_3d_double_array(T_electron_old);
memory->destroy_3d_double_array(T_electron);
memory->destroy_2d_double_array(flangevin);
memory->destroy_3d_double_array(net_energy_transfer);
memory->destroy_3d_double_array(net_energy_transfer_all);
}
/* ---------------------------------------------------------------------- */
int FixTTM::setmask()
{
int mask = 0;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixTTM::init()
{
if (domain->dimension == 2)
error->all("Cannot use fix ttm with 2d simulation");
if (domain->nonperiodic != 0)
error->all("Cannot use nonperiodic boundares with fix ttm");
if (domain->triclinic)
error->all("Cannot use fix ttm with triclinic box");
// set force prefactors
for (int i = 1; i <= atom->ntypes; i++) {
gfactor1[i] = - gamma_p / force->ftm2v;
gfactor2[i] =
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
net_energy_transfer_all[ixnode][iynode][iznode] = 0;
if (strcmp(update->integrate_style,"respa") == 0)
nlevels_respa = ((Respa *) update->integrate)->nlevels;
}
/* ---------------------------------------------------------------------- */
void FixTTM::setup(int vflag)
{
if (strcmp(update->integrate_style,"verlet") == 0)
post_force_setup(vflag);
else {
((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
post_force_respa_setup(vflag,nlevels_respa-1,0);
((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
}
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force(int vflag)
{
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double gamma1,gamma2;
// apply damping and thermostat to all atoms in fix group
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
if (T_electron[ixnode][iynode][iznode] < 0)
error->all("Electronic temperature dropped below zero");
double tsqrt = sqrt(T_electron[ixnode][iynode][iznode]);
gamma1 = gfactor1[type[i]];
double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
if (vsq > v_0_sq) gamma1 *= (gamma_p + gamma_s)/gamma_p;
gamma2 = gfactor2[type[i]] * tsqrt;
flangevin[i][0] = gamma1*v[i][0] + gamma2*(random->uniform()-0.5);
flangevin[i][1] = gamma1*v[i][1] + gamma2*(random->uniform()-0.5);
flangevin[i][2] = gamma1*v[i][2] + gamma2*(random->uniform()-0.5);
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_setup(int vflag)
{
double **f = atom->f;
int *mask = atom->mask;
int nlocal = atom->nlocal;
// apply langevin forces that have been stored from previous run
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
f[i][0] += flangevin[i][0];
f[i][1] += flangevin[i][1];
f[i][2] += flangevin[i][2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_respa(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixTTM::post_force_respa_setup(int vflag, int ilevel, int iloop)
{
if (ilevel == nlevels_respa-1) post_force_setup(vflag);
}
/* ---------------------------------------------------------------------- */
void FixTTM::reset_dt()
{
for (int i = 1; i <= atom->ntypes; i++)
gfactor2[i] =
sqrt(24.0*force->boltz*gamma_p/update->dt/force->mvv2e) / force->ftm2v;
}
/* ----------------------------------------------------------------------
read in initial electron temperatures from a user-specified file
only called by proc 0
------------------------------------------------------------------------- */
void FixTTM::read_initial_electron_temperatures()
{
char line[MAXLINE];
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
T_initial_set[ixnode][iynode][iznode] = 0;
// read initial electron temperature values from file
int ixnode,iynode,iznode;
double T_tmp;
while (1) {
if (fgets(line,MAXLINE,fpr) == NULL) break;
sscanf(line,"%d %d %d %lg",&ixnode,&iynode,&iznode,&T_tmp);
if (T_tmp < 0.0) error->one("Fix ttm electron temperatures must be > 0.0");
T_electron[ixnode][iynode][iznode] = T_tmp;
T_initial_set[ixnode][iynode][iznode] = 1;
}
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
if (T_initial_set[ixnode][iynode][iznode] == 0)
error->one("Initial temperatures not all set in fix ttm");
// close file
fclose(fpr);
}
/* ---------------------------------------------------------------------- */
void FixTTM::end_of_step()
{
double **x = atom->x;
double **v = atom->v;
double *mass = atom->mass;
double *rmass = atom->rmass;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
net_energy_transfer[ixnode][iynode][iznode] = 0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
net_energy_transfer[ixnode][iynode][iznode] +=
(flangevin[i][0]*v[i][0] + flangevin[i][1]*v[i][1] +
flangevin[i][2]*v[i][2]);
}
MPI_Allreduce(&net_energy_transfer[0][0][0],
&net_energy_transfer_all[0][0][0],
total_nnodes,MPI_DOUBLE,MPI_SUM,world);
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
// num_inner_timesteps = # of inner steps (thermal solves)
// required this MD step to maintain a stable explicit solve
int num_inner_timesteps = 1;
double inner_dt = update->dt;
double stability_criterion = 1.0 -
2.0*inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
if (stability_criterion < 0.0) {
inner_dt = 0.5*(electronic_specific_heat*electronic_density) /
(electronic_thermal_conductivity*(1.0/dx/dx + 1.0/dy/dy + 1.0/dz/dz));
num_inner_timesteps = static_cast<int>(update->dt/inner_dt) + 1;
inner_dt = update->dt/double(num_inner_timesteps);
if (num_inner_timesteps > 1000000)
error->warning("Too many inner timesteps in fix ttm");
}
for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;
ith_inner_timestep++) {
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
T_electron_old[ixnode][iynode][iznode] =
T_electron[ixnode][iynode][iznode];
// compute new electron T profile
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
int right_xnode = ixnode + 1;
int right_ynode = iynode + 1;
int right_znode = iznode + 1;
if (right_xnode == nxnodes) right_xnode = 0;
if (right_ynode == nynodes) right_ynode = 0;
if (right_znode == nznodes) right_znode = 0;
int left_xnode = ixnode - 1;
int left_ynode = iynode - 1;
int left_znode = iznode - 1;
if (left_xnode == -1) left_xnode = nxnodes - 1;
if (left_ynode == -1) left_ynode = nynodes - 1;
if (left_znode == -1) left_znode = nznodes - 1;
T_electron[ixnode][iynode][iznode] =
T_electron_old[ixnode][iynode][iznode] +
inner_dt/(electronic_specific_heat*electronic_density) *
(electronic_thermal_conductivity *
((T_electron_old[right_xnode][iynode][iznode] +
T_electron_old[left_xnode][iynode][iznode] -
2*T_electron_old[ixnode][iynode][iznode])/dx/dx +
(T_electron_old[ixnode][right_ynode][iznode] +
T_electron_old[ixnode][left_ynode][iznode] -
2*T_electron_old[ixnode][iynode][iznode])/dy/dy +
(T_electron_old[ixnode][iynode][right_znode] +
T_electron_old[ixnode][iynode][left_znode] -
2*T_electron_old[ixnode][iynode][iznode])/dz/dz) -
(net_energy_transfer_all[ixnode][iynode][iznode])/del_vol);
}
}
// output nodal temperatures for current timestep
if ((nfileevery) && !(update->ntimestep % nfileevery)) {
// compute atomic Ta for each grid point
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
nsum[ixnode][iynode][iznode] = 0;
nsum_all[ixnode][iynode][iznode] = 0;
sum_vsq[ixnode][iynode][iznode] = 0.0;
sum_mass_vsq[ixnode][iynode][iznode] = 0.0;
sum_vsq_all[ixnode][iynode][iznode] = 0.0;
sum_mass_vsq_all[ixnode][iynode][iznode] = 0.0;
}
double massone;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
double xscale = (x[i][0] - domain->boxlo[0])/domain->xprd;
double yscale = (x[i][1] - domain->boxlo[1])/domain->yprd;
double zscale = (x[i][2] - domain->boxlo[2])/domain->zprd;
int ixnode = static_cast<int>(xscale*nxnodes);
int iynode = static_cast<int>(yscale*nynodes);
int iznode = static_cast<int>(zscale*nznodes);
while (ixnode > nxnodes-1) ixnode -= nxnodes;
while (iynode > nynodes-1) iynode -= nynodes;
while (iznode > nznodes-1) iznode -= nznodes;
while (ixnode < 0) ixnode += nxnodes;
while (iynode < 0) iynode += nynodes;
while (iznode < 0) iznode += nznodes;
double vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
nsum[ixnode][iynode][iznode] += 1;
sum_vsq[ixnode][iynode][iznode] += vsq;
sum_mass_vsq[ixnode][iynode][iznode] += massone*vsq;
}
MPI_Allreduce(&nsum[0][0][0],&nsum_all[0][0][0],total_nnodes,
MPI_INT,MPI_SUM,world);
MPI_Allreduce(&sum_vsq[0][0][0],&sum_vsq_all[0][0][0],total_nnodes,
MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(&sum_mass_vsq[0][0][0],&sum_mass_vsq_all[0][0][0],
total_nnodes,MPI_DOUBLE,MPI_SUM,world);
if (me == 0) {
fprintf(fp,"%d ",update->ntimestep);
double T_a;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
T_a = 0;
if (nsum_all[ixnode][iynode][iznode] > 0)
T_a = sum_mass_vsq_all[ixnode][iynode][iznode]/
(3.0*force->boltz*nsum_all[ixnode][iynode][iznode]/force->mvv2e);
fprintf(fp,"%f ",T_a);
}
fprintf(fp,"\t");
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
fprintf(fp,"%f ",T_electron[ixnode][iynode][iznode]);
fprintf(fp,"\n");
}
}
}
/* ----------------------------------------------------------------------
memory usage of 3d grid
------------------------------------------------------------------------- */
double FixTTM::memory_usage()
{
double bytes = 0.0;
bytes += 5*total_nnodes * sizeof(int);
bytes += 14*total_nnodes * sizeof(double);
return bytes;
}
/* ---------------------------------------------------------------------- */
void FixTTM::grow_arrays(int ngrow)
{
flangevin = memory->grow_2d_double_array(flangevin,ngrow,3,"TTM:flangevin");
}
/* ----------------------------------------------------------------------
return the energy of the electronic subsystem or the net_energy transfer
between the subsystems
------------------------------------------------------------------------- */
double FixTTM::compute_vector(int n)
{
double e_energy = 0.0;
double transfer_energy = 0.0;
double dx = domain->xprd/nxnodes;
double dy = domain->yprd/nynodes;
double dz = domain->zprd/nznodes;
double del_vol = dx*dy*dz;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++) {
e_energy +=
T_electron[ixnode][iynode][iznode]*electronic_specific_heat*
electronic_density*del_vol;
transfer_energy +=
net_energy_transfer_all[ixnode][iynode][iznode]*update->dt;
}
if (n == 0) return e_energy;
if (n == 1) return transfer_energy;
return 0.0;
}
/* ----------------------------------------------------------------------
pack entire state of Fix into one write
------------------------------------------------------------------------- */
void FixTTM::write_restart(FILE *fp)
{
int n = 0;
double list[1 + nxnodes*nynodes*nznodes];
list[n++] = seed;
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
list[n++] = T_electron[ixnode][iynode][iznode];
if (comm->me == 0) {
int size = n * sizeof(double);
fwrite(&size,sizeof(int),1,fp);
- fwrite(&list,sizeof(double),n,fp);
+ fwrite(list,sizeof(double),n,fp);
}
}
/* ----------------------------------------------------------------------
use state info from restart file to restart the Fix
------------------------------------------------------------------------- */
void FixTTM::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
// the seed must be changed from the initial seed
seed = static_cast<int> (2*list[n++]);
for (int ixnode = 0; ixnode < nxnodes; ixnode++)
for (int iynode = 0; iynode < nynodes; iynode++)
for (int iznode = 0; iznode < nznodes; iznode++)
T_electron[ixnode][iynode][iznode] = list[n++];
delete random;
random = new RanMars(lmp,seed+comm->me);
}
/* ----------------------------------------------------------------------
pack values in local atom-based arrays for restart file
------------------------------------------------------------------------- */
int FixTTM::pack_restart(int i, double *buf)
{
buf[0] = 4;
buf[1] = flangevin[i][0];
buf[2] = flangevin[i][1];
buf[3] = flangevin[i][2];
return 4;
}
/* ----------------------------------------------------------------------
unpack values from atom->extra array to restart the fix
------------------------------------------------------------------------- */
void FixTTM::unpack_restart(int nlocal, int nth)
{
double **extra = atom->extra;
// skip to Nth set of extra values
int m = 0;
for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
m++;
flangevin[nlocal][0] = extra[nlocal][m++];
flangevin[nlocal][1] = extra[nlocal][m++];
flangevin[nlocal][2] = extra[nlocal][m++];
}
/* ----------------------------------------------------------------------
maxsize of any atom's restart data
------------------------------------------------------------------------- */
int FixTTM::maxsize_restart()
{
return 4;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
int FixTTM::size_restart(int nlocal)
{
return 4;
}