diff --git a/src/REPLICA/neb.cpp b/src/REPLICA/neb.cpp
index b3d08373b..8cb0be746 100644
--- a/src/REPLICA/neb.cpp
+++ b/src/REPLICA/neb.cpp
@@ -1,611 +1,616 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
+// lmptype.h must be first b/c this file uses MAXBIGINT and includes mpi.h
+// due to OpenMPI bug which sets INT64_MAX via its mpi.h
+// before lmptype.h can set flags to insure it is done correctly
+
+#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "neb.h"
#include "universe.h"
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "comm.h"
#include "min.h"
#include "modify.h"
#include "fix.h"
#include "fix_neb.h"
#include "output.h"
#include "thermo.h"
#include "finish.h"
#include "timer.h"
#include "memory.h"
#include "error.h"
#include "force.h"
using namespace LAMMPS_NS;
#define MAXLINE 256
#define CHUNK 1024
#define ATTRIBUTE_PERLINE 4
/* ---------------------------------------------------------------------- */
NEB::NEB(LAMMPS *lmp) : Pointers(lmp) {}
/* ----------------------------------------------------------------------
internal NEB constructor, called from TAD
------------------------------------------------------------------------- */
NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
int n2steps_in, int nevery_in, double *buf_init, double *buf_final)
: Pointers(lmp)
{
double delx,dely,delz;
etol = etol_in;
ftol = ftol_in;
n1steps = n1steps_in;
n2steps = n2steps_in;
nevery = nevery_in;
// replica info
nreplica = universe->nworlds;
ireplica = universe->iworld;
me_universe = universe->me;
uworld = universe->uworld;
MPI_Comm_rank(world,&me);
// generate linear interpolate replica
double fraction = ireplica/(nreplica-1.0);
double **x = atom->x;
int nlocal = atom->nlocal;
int ii = 0;
for (int i = 0; i < nlocal; i++) {
delx = buf_final[ii] - buf_init[ii];
dely = buf_final[ii+1] - buf_init[ii+1];
delz = buf_final[ii+2] - buf_init[ii+2];
domain->minimum_image(delx,dely,delz);
x[i][0] = buf_init[ii] + fraction*delx;
x[i][1] = buf_init[ii+1] + fraction*dely;
x[i][2] = buf_init[ii+2] + fraction*delz;
ii += 3;
}
}
/* ---------------------------------------------------------------------- */
NEB::~NEB()
{
MPI_Comm_free(&roots);
memory->destroy(all);
delete [] rdist;
}
/* ----------------------------------------------------------------------
perform NEB on multiple replicas
------------------------------------------------------------------------- */
void NEB::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"NEB command before simulation box is defined");
if (narg < 6) error->universe_all(FLERR,"Illegal NEB command");
etol = force->numeric(FLERR,arg[0]);
ftol = force->numeric(FLERR,arg[1]);
n1steps = force->inumeric(FLERR,arg[2]);
n2steps = force->inumeric(FLERR,arg[3]);
nevery = force->inumeric(FLERR,arg[4]);
// error checks
if (etol < 0.0) error->all(FLERR,"Illegal NEB command");
if (ftol < 0.0) error->all(FLERR,"Illegal NEB command");
if (nevery == 0) error->universe_all(FLERR,"Illegal NEB command");
if (n1steps % nevery || n2steps % nevery)
error->universe_all(FLERR,"Illegal NEB command");
// replica info
nreplica = universe->nworlds;
ireplica = universe->iworld;
me_universe = universe->me;
uworld = universe->uworld;
MPI_Comm_rank(world,&me);
// error checks
if (nreplica == 1) error->all(FLERR,"Cannot use NEB with a single replica");
if (nreplica != universe->nprocs)
error->all(FLERR,"Can only use NEB with 1-processor replicas");
if (atom->sortfreq > 0)
error->all(FLERR,"Cannot use NEB with atom_modify sort enabled");
if (atom->map_style == 0)
error->all(FLERR,"Cannot use NEB unless atom map exists");
// process file-style setting to setup initial configs for all replicas
if (strcmp(arg[5],"final") == 0) {
if (narg != 7) error->universe_all(FLERR,"Illegal NEB command");
infile = arg[6];
readfile(infile,0);
} else if (strcmp(arg[5],"each") == 0) {
if (narg != 7) error->universe_all(FLERR,"Illegal NEB command");
infile = arg[6];
readfile(infile,1);
} else if (strcmp(arg[5],"none") == 0) {
if (narg != 6) error->universe_all(FLERR,"Illegal NEB command");
} else error->universe_all(FLERR,"Illegal NEB command");
// run the NEB calculation
run();
}
/* ----------------------------------------------------------------------
run NEB on multiple replicas
------------------------------------------------------------------------- */
void NEB::run()
{
// create MPI communicator for root proc from each world
int color;
if (me == 0) color = 0;
else color = 1;
MPI_Comm_split(uworld,color,0,&roots);
int ineb;
for (ineb = 0; ineb < modify->nfix; ineb++)
if (strcmp(modify->fix[ineb]->style,"neb") == 0) break;
if (ineb == modify->nfix) error->all(FLERR,"NEB requires use of fix neb");
fneb = (FixNEB *) modify->fix[ineb];
nall = 4;
memory->create(all,nreplica,nall,"neb:all");
rdist = new double[nreplica];
// initialize LAMMPS
update->whichflag = 2;
update->etol = etol;
update->ftol = ftol;
update->multireplica = 1;
lmp->init();
if (update->minimize->searchflag)
error->all(FLERR,"NEB requires damped dynamics minimizer");
// setup regular NEB minimization
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Setting up regular NEB ...\n");
update->beginstep = update->firststep = update->ntimestep;
update->endstep = update->laststep = update->firststep + n1steps;
update->nsteps = n1steps;
update->max_eval = n1steps;
if (update->laststep < 0 || update->laststep > MAXBIGINT)
error->all(FLERR,"Too many timesteps for NEB");
update->minimize->setup();
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
}
print_status();
// perform regular NEB for n1steps or until replicas converge
// retrieve PE values from fix NEB and print every nevery iterations
// break induced if converged
// damped dynamic min styles insure all replicas converge together
timer->init();
timer->barrier_start();
while (update->minimize->niter < n1steps) {
update->minimize->run(nevery);
print_status();
if (update->minimize->stop_condition) break;
}
timer->barrier_stop();
update->minimize->cleanup();
Finish finish(lmp);
finish.end(1);
// switch fix NEB to climbing mode
// top = replica that becomes hill climber
double vmax = all[0][0];
int top = 0;
for (int m = 1; m < nreplica; m++)
if (vmax < all[m][0]) {
vmax = all[m][0];
top = m;
}
// setup climbing NEB minimization
// must reinitialize minimizer so it re-creates its fix MINIMIZE
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Setting up climbing ...\n");
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,"Climbing replica = %d\n",top+1);
if (universe->ulogfile)
fprintf(universe->ulogfile,"Climbing replica = %d\n",top+1);
}
update->beginstep = update->firststep = update->ntimestep;
update->endstep = update->laststep = update->firststep + n2steps;
update->nsteps = n2steps;
update->max_eval = n2steps;
if (update->laststep < 0 || update->laststep > MAXBIGINT)
error->all(FLERR,"Too many timesteps");
update->minimize->init();
fneb->rclimber = top;
update->minimize->setup();
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,"Step MaxReplicaForce MaxAtomForce "
"GradV0 GradV1 GradVc "
"EBF EBR RDT "
"RD1 PE1 RD2 PE2 ... RDN PEN\n");
}
print_status();
// perform climbing NEB for n2steps or until replicas converge
// retrieve PE values from fix NEB and print every nevery iterations
// break induced if converged
// damped dynamic min styles insure all replicas converge together
timer->init();
timer->barrier_start();
while (update->minimize->niter < n2steps) {
update->minimize->run(nevery);
print_status();
if (update->minimize->stop_condition) break;
}
timer->barrier_stop();
update->minimize->cleanup();
finish.end(1);
update->whichflag = 0;
update->multireplica = 0;
update->firststep = update->laststep = 0;
update->beginstep = update->endstep = 0;
}
/* ----------------------------------------------------------------------
read initial config atom coords from file
flag = 0
only first replica opens file and reads it
first replica bcasts lines to all replicas
final replica stores coords
intermediate replicas interpolate from coords
new coord = replica fraction between current and final state
initial replica does nothing
flag = 1
each replica (except first) opens file and reads it
each replica stores coords
initial replica does nothing
------------------------------------------------------------------------- */
void NEB::readfile(char *file, int flag)
{
int i,j,m,nchunk,eofflag,nlines;
tagint tag;
char *eof,*start,*next,*buf;
char line[MAXLINE];
double xx,yy,zz,delx,dely,delz;
if (me_universe == 0 && screen)
fprintf(screen,"Reading NEB coordinate file(s) ...\n");
// flag = 0, universe root reads header of file, bcast to universe
// flag = 1, each replica's root reads header of file, bcast to world
// but explicitly skip first replica
if (flag == 0) {
if (me_universe == 0) {
open(file);
while (1) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of neb file");
start = &line[strspn(line," \t\n\v\f\r")];
if (*start != '\0' && *start != '#') break;
}
sscanf(line,"%d",&nlines);
}
MPI_Bcast(&nlines,1,MPI_INT,0,uworld);
} else {
if (me == 0) {
if (ireplica) {
open(file);
while (1) {
eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of neb file");
start = &line[strspn(line," \t\n\v\f\r")];
if (*start != '\0' && *start != '#') break;
}
sscanf(line,"%d",&nlines);
} else nlines = 0;
}
MPI_Bcast(&nlines,1,MPI_INT,0,world);
}
char *buffer = new char[CHUNK*MAXLINE];
char **values = new char*[ATTRIBUTE_PERLINE];
double fraction = ireplica/(nreplica-1.0);
double **x = atom->x;
int nlocal = atom->nlocal;
// loop over chunks of lines read from file
// two versions of read_lines_from_file() for world vs universe bcast
// count # of atom coords changed so can check for invalid atom IDs in file
int ncount = 0;
int nread = 0;
while (nread < nlines) {
nchunk = MIN(nlines-nread,CHUNK);
if (flag == 0)
eofflag = comm->read_lines_from_file_universe(fp,nchunk,MAXLINE,buffer);
else
eofflag = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eofflag) error->all(FLERR,"Unexpected end of neb file");
buf = buffer;
next = strchr(buf,'\n');
*next = '\0';
int nwords = atom->count_words(buf);
*next = '\n';
if (nwords != ATTRIBUTE_PERLINE)
error->all(FLERR,"Incorrect atom format in neb file");
// loop over lines of atom coords
// tokenize the line into values
for (i = 0; i < nchunk; i++) {
next = strchr(buf,'\n');
values[0] = strtok(buf," \t\n\r\f");
for (j = 1; j < nwords; j++)
values[j] = strtok(NULL," \t\n\r\f");
// adjust atom coord based on replica fraction
// for flag = 0, interpolate for intermediate and final replicas
// for flag = 1, replace existing coord with new coord
// ignore image flags of final x
// for interpolation:
// new x is displacement from old x via minimum image convention
// if final x is across periodic boundary:
// new x may be outside box
// will be remapped back into box when simulation starts
// its image flags will then be adjusted
tag = ATOTAGINT(values[0]);
m = atom->map(tag);
if (m >= 0 && m < nlocal) {
ncount++;
xx = atof(values[1]);
yy = atof(values[2]);
zz = atof(values[3]);
if (flag == 0) {
delx = xx - x[m][0];
dely = yy - x[m][1];
delz = zz - x[m][2];
domain->minimum_image(delx,dely,delz);
x[m][0] += fraction*delx;
x[m][1] += fraction*dely;
x[m][2] += fraction*delz;
} else {
x[m][0] = xx;
x[m][1] = yy;
x[m][2] = zz;
}
}
buf = next + 1;
}
nread += nchunk;
}
// check that all atom IDs in file were found by a proc
if (flag == 0) {
int ntotal;
MPI_Allreduce(&ncount,&ntotal,1,MPI_INT,MPI_SUM,uworld);
if (ntotal != nreplica*nlines)
error->universe_all(FLERR,"Invalid atom IDs in neb file");
} else {
int ntotal;
MPI_Allreduce(&ncount,&ntotal,1,MPI_INT,MPI_SUM,world);
if (ntotal != nlines)
error->all(FLERR,"Invalid atom IDs in neb file");
}
// clean up
delete [] buffer;
delete [] values;
if (flag == 0) {
if (me_universe == 0) {
if (compressed) pclose(fp);
else fclose(fp);
}
} else {
if (me == 0 && ireplica) {
if (compressed) pclose(fp);
else fclose(fp);
}
}
}
/* ----------------------------------------------------------------------
universe proc 0 opens NEB data file
test if gzipped
------------------------------------------------------------------------- */
void NEB::open(char *file)
{
compressed = 0;
char *suffix = file + strlen(file) - 3;
if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1;
if (!compressed) fp = fopen(file,"r");
else {
#ifdef LAMMPS_GZIP
char gunzip[128];
sprintf(gunzip,"gzip -c -d %s",file);
#ifdef _WIN32
fp = _popen(gunzip,"rb");
#else
fp = popen(gunzip,"r");
#endif
#else
error->one(FLERR,"Cannot open gzipped file");
#endif
}
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
}
/* ----------------------------------------------------------------------
query fix NEB for info on each replica
proc 0 prints current NEB status
------------------------------------------------------------------------- */
void NEB::print_status()
{
double fnorm2 = sqrt(update->minimize->fnorm_sqr());
double fmaxreplica;
MPI_Allreduce(&fnorm2,&fmaxreplica,1,MPI_DOUBLE,MPI_MAX,roots);
double fnorminf = update->minimize->fnorm_inf();
double fmaxatom;
MPI_Allreduce(&fnorminf,&fmaxatom,1,MPI_DOUBLE,MPI_MAX,roots);
double one[4];
one[0] = fneb->veng;
one[1] = fneb->plen;
one[2] = fneb->nlen;
one[nall-1] = fneb->gradvnorm;
if (output->thermo->normflag) one[0] /= atom->natoms;
if (me == 0)
MPI_Allgather(one,nall,MPI_DOUBLE,&all[0][0],nall,MPI_DOUBLE,roots);
rdist[0] = 0.0;
for (int i = 1; i < nreplica; i++)
rdist[i] = rdist[i-1] + all[i][1];
double endpt = rdist[nreplica-1] = rdist[nreplica-2] + all[nreplica-2][2];
for (int i = 1; i < nreplica; i++)
rdist[i] /= endpt;
// look up GradV for the initial, final, and climbing replicas
// these are identical to fnorm2, but to be safe we
// take them straight from fix_neb
double gradvnorm0, gradvnorm1, gradvnormc;
int irep;
irep = 0;
gradvnorm0 = all[irep][3];
irep = nreplica-1;
gradvnorm1 = all[irep][3];
irep = fneb->rclimber;
if (irep > -1) {
gradvnormc = all[irep][3];
ebf = all[irep][0]-all[0][0];
ebr = all[irep][0]-all[nreplica-1][0];
} else {
double vmax = all[0][0];
int top = 0;
for (int m = 1; m < nreplica; m++)
if (vmax < all[m][0]) {
vmax = all[m][0];
top = m;
}
irep = top;
gradvnormc = all[irep][3];
ebf = all[irep][0]-all[0][0];
ebr = all[irep][0]-all[nreplica-1][0];
}
if (me_universe == 0) {
if (universe->uscreen) {
fprintf(universe->uscreen,BIGINT_FORMAT " %12.8g %12.8g ",
update->ntimestep,fmaxreplica,fmaxatom);
fprintf(universe->uscreen,"%12.8g %12.8g %12.8g ",
gradvnorm0,gradvnorm1,gradvnormc);
fprintf(universe->uscreen,"%12.8g %12.8g %12.8g ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++)
fprintf(universe->uscreen,"%12.8g %12.8g ",rdist[i],all[i][0]);
fprintf(universe->uscreen,"\n");
}
if (universe->ulogfile) {
fprintf(universe->ulogfile,BIGINT_FORMAT " %12.8g %12.8g ",
update->ntimestep,fmaxreplica,fmaxatom);
fprintf(universe->ulogfile,"%12.8g %12.8g %12.8g ",
gradvnorm0,gradvnorm1,gradvnormc);
fprintf(universe->ulogfile,"%12.8g %12.8g %12.8g ",ebf,ebr,endpt);
for (int i = 0; i < nreplica; i++)
fprintf(universe->ulogfile,"%12.8g %12.8g ",rdist[i],all[i][0]);
fprintf(universe->ulogfile,"\n");
fflush(universe->ulogfile);
}
}
}
diff --git a/src/REPLICA/prd.cpp b/src/REPLICA/prd.cpp
index 383d07dfd..35fb181a7 100644
--- a/src/REPLICA/prd.cpp
+++ b/src/REPLICA/prd.cpp
@@ -1,903 +1,908 @@
/* ----------------------------------------------------------------------
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)
------------------------------------------------------------------------- */
+// lmptype.h must be first b/c this file uses MAXBIGINT and includes mpi.h
+// due to OpenMPI bug which sets INT64_MAX via its mpi.h
+// before lmptype.h can set flags to insure it is done correctly
+
+#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "prd.h"
#include "universe.h"
#include "update.h"
#include "atom.h"
#include "domain.h"
#include "region.h"
#include "comm.h"
#include "velocity.h"
#include "integrate.h"
#include "min.h"
#include "neighbor.h"
#include "modify.h"
#include "compute.h"
#include "fix.h"
#include "fix_event_prd.h"
#include "force.h"
#include "pair.h"
#include "random_park.h"
#include "random_mars.h"
#include "output.h"
#include "dump.h"
#include "finish.h"
#include "timer.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PRD::PRD(LAMMPS *lmp) : Pointers(lmp) {}
/* ----------------------------------------------------------------------
perform PRD simulation on one or more replicas
------------------------------------------------------------------------- */
void PRD::command(int narg, char **arg)
{
int ireplica;
// error checks
if (domain->box_exist == 0)
error->all(FLERR,"PRD command before simulation box is defined");
if (universe->nworlds != universe->nprocs &&
atom->map_style == 0)
error->all(FLERR,"Cannot use PRD with multi-processor replicas "
"unless atom map exists");
if (universe->nworlds == 1 && comm->me == 0)
error->warning(FLERR,"Running PRD with only one replica");
if (narg < 7) error->universe_all(FLERR,"Illegal prd command");
// read as double so can cast to bigint
int nsteps = force->inumeric(FLERR,arg[0]);
t_event = force->inumeric(FLERR,arg[1]);
n_dephase = force->inumeric(FLERR,arg[2]);
t_dephase = force->inumeric(FLERR,arg[3]);
t_corr = force->inumeric(FLERR,arg[4]);
char *id_compute = new char[strlen(arg[5])+1];
strcpy(id_compute,arg[5]);
int seed = force->inumeric(FLERR,arg[6]);
options(narg-7,&arg[7]);
// total # of timesteps must be multiple of t_event
if (t_event <= 0) error->universe_all(FLERR,"Invalid t_event in prd command");
if (nsteps % t_event)
error->universe_all(FLERR,"PRD nsteps must be multiple of t_event");
if (t_corr % t_event)
error->universe_all(FLERR,"PRD t_corr must be multiple of t_event");
// local storage
int me_universe = universe->me;
int nprocs_universe = universe->nprocs;
int nreplica = universe->nworlds;
int iworld = universe->iworld;
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
// comm_replica = communicator between all proc 0s across replicas
int color = me;
MPI_Comm_split(universe->uworld,color,0,&comm_replica);
// equal_size_replicas = 1 if all replicas have same # of procs
// no longer used
//flag = 0;
//if (nreplica*nprocs == nprocs_universe) flag = 1;
//MPI_Allreduce(&flag,&equal_size_replicas,1,MPI_INT,MPI_MIN,
// universe->uworld);
// workspace for inter-replica communication via gathers
natoms = atom->natoms;
displacements = NULL;
tagall = NULL;
xall = NULL;
imageall = NULL;
if (nreplica != nprocs_universe) {
displacements = new int[nprocs];
memory->create(tagall,natoms,"prd:tagall");
memory->create(xall,natoms,3,"prd:xall");
memory->create(imageall,natoms,"prd:imageall");
}
// random_select = same RNG for each replica, for multiple event selection
// random_clock = same RNG for each replica, for clock updates
// random_dephase = unique RNG for each replica, for dephasing
random_select = new RanPark(lmp,seed);
random_clock = new RanPark(lmp,seed+1000);
random_dephase = new RanMars(lmp,seed+iworld);
// create ComputeTemp class to monitor temperature
char **args = new char*[3];
args[0] = (char *) "prd_temp";
args[1] = (char *) "all";
args[2] = (char *) "temp";
modify->add_compute(3,args);
temperature = modify->compute[modify->ncompute-1];
// create Velocity class for velocity creation in dephasing
// pass it temperature compute, loop_setting, dist_setting settings
atom->check_mass();
velocity = new Velocity(lmp);
velocity->init_external("all");
args[0] = (char *) "temp";
args[1] = (char *) "prd_temp";
velocity->options(2,args);
args[0] = (char *) "loop";
args[1] = (char *) loop_setting;
if (loop_setting) velocity->options(2,args);
args[0] = (char *) "dist";
args[1] = (char *) dist_setting;
if (dist_setting) velocity->options(2,args);
// create FixEventPRD class to store event and pre-quench states
args[0] = (char *) "prd_event";
args[1] = (char *) "all";
args[2] = (char *) "EVENT/PRD";
modify->add_fix(3,args);
fix_event = (FixEventPRD *) modify->fix[modify->nfix-1];
// create Finish for timing output
finish = new Finish(lmp);
// string clean-up
delete [] args;
delete [] loop_setting;
delete [] dist_setting;
// assign FixEventPRD to event-detection compute
// necessary so it will know atom coords at last event
int icompute = modify->find_compute(id_compute);
if (icompute < 0) error->all(FLERR,"Could not find compute ID for PRD");
compute_event = modify->compute[icompute];
compute_event->reset_extra_compute_fix("prd_event");
// reset reneighboring criteria since will perform minimizations
neigh_every = neighbor->every;
neigh_delay = neighbor->delay;
neigh_dist_check = neighbor->dist_check;
if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
if (me == 0)
error->warning(FLERR,"Resetting reneighboring criteria during PRD");
}
neighbor->every = 1;
neighbor->delay = 0;
neighbor->dist_check = 1;
// initialize PRD as if one long dynamics run
update->whichflag = 1;
update->nsteps = nsteps;
update->beginstep = update->firststep = update->ntimestep;
update->endstep = update->laststep = update->firststep + nsteps;
update->restrict_output = 1;
if (update->laststep < 0 || update->laststep > MAXBIGINT)
error->all(FLERR,"Too many timesteps");
lmp->init();
// init minimizer settings and minimizer itself
update->etol = etol;
update->ftol = ftol;
update->max_eval = maxeval;
update->minimize->init();
// cannot use PRD with a changing box
if (domain->box_change)
error->all(FLERR,"Cannot use PRD with a changing box");
// cannot use PRD with time-dependent fixes or regions or atom sorting
for (int i = 0; i < modify->nfix; i++)
if (modify->fix[i]->time_depend)
error->all(FLERR,"Cannot use PRD with a time-dependent fix defined");
for (int i = 0; i < domain->nregion; i++)
if (domain->regions[i]->dynamic_check())
error->all(FLERR,"Cannot use PRD with a time-dependent region defined");
if (atom->sortfreq > 0)
error->all(FLERR,"Cannot use PRD with atom_modify sort enabled");
// perform PRD simulation
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Setting up PRD ...\n");
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,"Step CPU Clock Event "
"Correlated Coincident Replica\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,"Step CPU Clock Event "
"Correlated Coincident Replica\n");
}
// store hot state and quenched event for replica 0
// use share_event() to copy that info to all replicas
// this insures all start from same place
// need this line if quench() does only setup_minimal()
// update->minimize->setup();
fix_event->store_state_quench();
quench();
ncoincident = 0;
share_event(0,0,0);
timer->init();
timer->barrier_start();
time_start = timer->get_wall(Timer::TOTAL);
log_event();
// do full init/setup since are starting all replicas after event
// replica 0 bcasts temp to all replicas if temp_dephase is not set
update->whichflag = 1;
lmp->init();
update->integrate->setup();
if (temp_flag == 0) {
if (universe->iworld == 0) temp_dephase = temperature->compute_scalar();
MPI_Bcast(&temp_dephase,1,MPI_DOUBLE,universe->root_proc[0],
universe->uworld);
}
// main loop: look for events until out of time
// (1) dephase independently on each proc after event
// (2) loop: dynamics, store state, quench, check event, restore state
// (3) share and record event
nbuild = ndanger = 0;
time_dephase = time_dynamics = time_quench = time_comm = time_output = 0.0;
bigint clock = 0;
timer->barrier_start();
time_start = timer->get_wall(Timer::TOTAL);
int istep = 0;
while (istep < nsteps) {
dephase();
if (stepmode == 0) istep = update->ntimestep - update->beginstep;
else istep = clock;
ireplica = -1;
while (istep < nsteps) {
dynamics(t_event,time_dynamics);
fix_event->store_state_quench();
quench();
clock = clock + t_event*universe->nworlds;
ireplica = check_event();
if (ireplica >= 0) break;
fix_event->restore_state_quench();
if (stepmode == 0) istep = update->ntimestep - update->beginstep;
else istep = clock;
}
if (ireplica < 0) break;
// decrement clock by random time at which 1 or more events occurred
int frac_t_event = t_event;
for (int i = 0; i <= fix_event->ncoincident; i++) {
int frac_rand = static_cast<int> (random_clock->uniform() * t_event);
frac_t_event = MIN(frac_t_event,frac_rand);
}
int decrement = frac_t_event*universe->nworlds;
clock -= decrement;
// share event across replicas
// NOTE: would be potentially more efficient for correlated events
// if don't share until correlated check below has completed
// this will complicate the dump (always on replica 0)
share_event(ireplica,1,decrement);
log_event();
int restart_flag = 0;
if (output->restart_flag && universe->iworld == 0) {
if (output->restart_every_single &&
fix_event->event_number % output->restart_every_single == 0)
restart_flag = 1;
if (output->restart_every_double &&
fix_event->event_number % output->restart_every_double == 0)
restart_flag = 1;
}
// correlated event loop
// other procs could be dephasing during this time
int corr_endstep = update->ntimestep + t_corr;
while (update->ntimestep < corr_endstep) {
if (update->ntimestep == update->endstep) {
restart_flag = 0;
break;
}
dynamics(t_event,time_dynamics);
fix_event->store_state_quench();
quench();
clock += t_event;
int corr_event_check = check_event(ireplica);
if (corr_event_check >= 0) {
share_event(ireplica,2,0);
log_event();
corr_endstep = update->ntimestep + t_corr;
} else fix_event->restore_state_quench();
}
// full init/setup since are starting all replicas after event
// event replica bcasts temp to all replicas if temp_dephase is not set
update->whichflag = 1;
lmp->init();
update->integrate->setup();
timer->barrier_start();
if (t_corr > 0) replicate(ireplica);
if (temp_flag == 0) {
if (ireplica == universe->iworld)
temp_dephase = temperature->compute_scalar();
MPI_Bcast(&temp_dephase,1,MPI_DOUBLE,universe->root_proc[ireplica],
universe->uworld);
}
timer->barrier_stop();
time_comm += timer->get_wall(Timer::TOTAL);
// write restart file of hot coords
if (restart_flag) {
timer->barrier_start();
output->write_restart(update->ntimestep);
timer->barrier_stop();
time_output += timer->get_wall(Timer::TOTAL);
}
if (stepmode == 0) istep = update->ntimestep - update->beginstep;
else istep = clock;
}
if (stepmode) nsteps = update->ntimestep - update->beginstep;
// set total timers and counters so Finish() will process them
timer->set_wall(Timer::TOTAL, time_start);
timer->barrier_stop();
timer->set_wall(Timer::DEPHASE, time_dephase);
timer->set_wall(Timer::DYNAMICS, time_dynamics);
timer->set_wall(Timer::QUENCH, time_quench);
timer->set_wall(Timer::REPCOMM, time_comm);
timer->set_wall(Timer::REPOUT, time_output);
neighbor->ncalls = nbuild;
neighbor->ndanger = ndanger;
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,
"Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT
" atoms\n",
timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms);
if (universe->ulogfile)
fprintf(universe->ulogfile,
"Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT
" atoms\n",
timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms);
}
if (me == 0) {
if (screen) fprintf(screen,"\nPRD done\n");
if (logfile) fprintf(logfile,"\nPRD done\n");
}
finish->end(2);
update->whichflag = 0;
update->firststep = update->laststep = 0;
update->beginstep = update->endstep = 0;
update->restrict_output = 0;
// reset reneighboring criteria
neighbor->every = neigh_every;
neighbor->delay = neigh_delay;
neighbor->dist_check = neigh_dist_check;
// clean up
delete [] displacements;
memory->destroy(tagall);
memory->destroy(xall);
memory->destroy(imageall);
delete [] id_compute;
MPI_Comm_free(&comm_replica);
delete random_select;
delete random_clock;
delete random_dephase;
delete velocity;
delete finish;
modify->delete_compute("prd_temp");
modify->delete_fix("prd_event");
compute_event->reset_extra_compute_fix(NULL);
}
/* ----------------------------------------------------------------------
dephasing = one or more short runs with new random velocities
------------------------------------------------------------------------- */
void PRD::dephase()
{
bigint ntimestep_hold = update->ntimestep;
// n_dephase iterations of dephasing, each of t_dephase steps
for (int i = 0; i < n_dephase; i++) {
fix_event->store_state_dephase();
// do not proceed to next iteration until an event-free run occurs
int done = 0;
while (!done) {
int seed = static_cast<int> (random_dephase->uniform() * MAXSMALLINT);
if (seed == 0) seed = 1;
velocity->create(temp_dephase,seed);
dynamics(t_dephase,time_dephase);
fix_event->store_state_quench();
quench();
if (compute_event->compute_scalar() > 0.0) {
fix_event->restore_state_dephase();
update->ntimestep -= t_dephase;
} else {
fix_event->restore_state_quench();
done = 1;
}
if (temp_flag == 0) temp_dephase = temperature->compute_scalar();
}
}
// reset timestep as if dephase did not occur
// clear timestep storage from computes, since now invalid
update->ntimestep = ntimestep_hold;
for (int i = 0; i < modify->ncompute; i++)
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
}
/* ----------------------------------------------------------------------
short dynamics run: for event search, decorrelation, or dephasing
------------------------------------------------------------------------- */
void PRD::dynamics(int nsteps, double &time_category)
{
update->whichflag = 1;
update->nsteps = nsteps;
lmp->init();
update->integrate->setup();
// this may be needed if don't do full init
//modify->addstep_compute_all(update->ntimestep);
bigint ncalls = neighbor->ncalls;
timer->barrier_start();
update->integrate->run(nsteps);
timer->barrier_stop();
time_category += timer->get_wall(Timer::TOTAL);
nbuild += neighbor->ncalls - ncalls;
ndanger += neighbor->ndanger;
update->integrate->cleanup();
finish->end(0);
}
/* ----------------------------------------------------------------------
quench minimization
------------------------------------------------------------------------- */
void PRD::quench()
{
bigint ntimestep_hold = update->ntimestep;
bigint endstep_hold = update->endstep;
// need to change whichflag so that minimize->setup() calling
// modify->setup() will call fix->min_setup()
update->whichflag = 2;
update->nsteps = maxiter;
update->endstep = update->laststep = update->firststep + maxiter;
if (update->laststep < 0 || update->laststep > MAXBIGINT)
error->all(FLERR,"Too many iterations");
// full init works
lmp->init();
update->minimize->setup();
// partial init does not work
//modify->addstep_compute_all(update->ntimestep);
//update->minimize->setup_minimal(1);
int ncalls = neighbor->ncalls;
timer->barrier_start();
update->minimize->run(maxiter);
timer->barrier_stop();
time_quench += timer->get_wall(Timer::TOTAL);
if (neighbor->ncalls == ncalls) quench_reneighbor = 0;
else quench_reneighbor = 1;
update->minimize->cleanup();
finish->end(0);
// reset timestep as if quench did not occur
// clear timestep storage from computes, since now invalid
update->ntimestep = ntimestep_hold;
update->endstep = update->laststep = endstep_hold;
for (int i = 0; i < modify->ncompute; i++)
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
}
/* ----------------------------------------------------------------------
check for an event in any replica
if replica_num is non-negative only check for event on replica_num
if multiple events, choose one at random
return -1 if no event
else return ireplica = world in which event occured
------------------------------------------------------------------------- */
int PRD::check_event(int replica_num)
{
int worldflag,universeflag,scanflag,replicaflag,ireplica;
worldflag = 0;
if (compute_event->compute_scalar() > 0.0) worldflag = 1;
if (replica_num >= 0 && replica_num != universe->iworld) worldflag = 0;
timer->barrier_start();
if (me == 0) MPI_Allreduce(&worldflag,&universeflag,1,
MPI_INT,MPI_SUM,comm_replica);
MPI_Bcast(&universeflag,1,MPI_INT,0,world);
ncoincident = universeflag;
if (!universeflag) ireplica = -1;
else {
// multiple events, choose one at random
// iwhich = random # from 1 to N, N = # of events to choose from
// scanflag = 1 to N on replicas with an event, 0 on non-event replicas
// exit with worldflag = 1 on chosen replica, 0 on all others
// note worldflag is already 0 on replicas that didn't perform event
if (universeflag > 1) {
int iwhich = static_cast<int>
(universeflag*random_select->uniform()) + 1;
if (me == 0)
MPI_Scan(&worldflag,&scanflag,1,MPI_INT,MPI_SUM,comm_replica);
MPI_Bcast(&scanflag,1,MPI_INT,0,world);
if (scanflag != iwhich) worldflag = 0;
}
if (worldflag) replicaflag = universe->iworld;
else replicaflag = 0;
if (me == 0) MPI_Allreduce(&replicaflag,&ireplica,1,
MPI_INT,MPI_SUM,comm_replica);
MPI_Bcast(&ireplica,1,MPI_INT,0,world);
}
timer->barrier_stop();
time_comm += timer->get_wall(Timer::TOTAL);
return ireplica;
}
/* ----------------------------------------------------------------------
share quenched and hot coords owned by ireplica with all replicas
all replicas store event in fix_event
replica 0 dumps event snapshot
flag = 0 = called before PRD run
flag = 1 = called during PRD run = not correlated event
flag = 2 = called during PRD run = correlated event
------------------------------------------------------------------------- */
void PRD::share_event(int ireplica, int flag, int decrement)
{
timer->barrier_start();
// communicate quenched coords to all replicas and store as event
// decrement event counter if flag = 0 since not really an event
replicate(ireplica);
timer->barrier_stop();
time_comm += timer->get_wall(Timer::TOTAL);
// adjust time for last correlated event check (not on first event)
int corr_adjust = t_corr;
if (fix_event->event_number < 1 || flag == 2) corr_adjust = 0;
// delta = time since last correlated event check
int delta = update->ntimestep - fix_event->event_timestep - corr_adjust;
// if this is a correlated event, time elapsed only on one partition
if (flag != 2) delta *= universe->nworlds;
if (delta > 0 && flag != 2) delta -= decrement;
delta += corr_adjust;
// delta passed to store_event_prd() should make its clock update
// be consistent with clock in main PRD loop
// don't change the clock or timestep if this is a restart
if (flag == 0 && fix_event->event_number != 0)
fix_event->store_event_prd(fix_event->event_timestep,0);
else {
fix_event->store_event_prd(update->ntimestep,delta);
fix_event->replica_number = ireplica;
fix_event->correlated_event = 0;
if (flag == 2) fix_event->correlated_event = 1;
fix_event->ncoincident = ncoincident;
}
if (flag == 0) fix_event->event_number--;
// dump snapshot of quenched coords, only on replica 0
// must reneighbor and compute forces before dumping
// since replica 0 possibly has new state from another replica
// addstep_compute_all insures eng/virial are calculated if needed
if (output->ndump && universe->iworld == 0) {
timer->barrier_start();
modify->addstep_compute_all(update->ntimestep);
update->integrate->setup_minimal(1);
output->write_dump(update->ntimestep);
timer->barrier_stop();
time_output += timer->get_wall(Timer::TOTAL);
}
// restore and communicate hot coords to all replicas
fix_event->restore_state_quench();
timer->barrier_start();
replicate(ireplica);
timer->barrier_stop();
time_comm += timer->get_wall(Timer::TOTAL);
}
/* ----------------------------------------------------------------------
universe proc 0 prints event info
------------------------------------------------------------------------- */
void PRD::log_event()
{
timer->set_wall(Timer::TOTAL, time_start);
if (universe->me == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,
BIGINT_FORMAT " %.3f " BIGINT_FORMAT " %d %d %d %d\n",
fix_event->event_timestep,
timer->elapsed(Timer::TOTAL),
fix_event->clock,
fix_event->event_number,fix_event->correlated_event,
fix_event->ncoincident,
fix_event->replica_number);
if (universe->ulogfile)
fprintf(universe->ulogfile,
BIGINT_FORMAT " %.3f " BIGINT_FORMAT " %d %d %d %d\n",
fix_event->event_timestep,
timer->elapsed(Timer::TOTAL),
fix_event->clock,
fix_event->event_number,fix_event->correlated_event,
fix_event->ncoincident,
fix_event->replica_number);
}
}
/* ----------------------------------------------------------------------
communicate atom coords and image flags in ireplica to all other replicas
one proc per replica:
direct overwrite via bcast
else atoms could be stored in different order or on different procs:
collect to root proc of event replica
bcast to roots of other replicas
bcast within each replica
each proc extracts info for atoms it owns using atom IDs
------------------------------------------------------------------------- */
void PRD::replicate(int ireplica)
{
int nreplica = universe->nworlds;
int nprocs_universe = universe->nprocs;
int i,m;
if (nreplica == nprocs_universe) {
MPI_Bcast(atom->image,atom->nlocal,MPI_INT,ireplica,comm_replica);
MPI_Bcast(atom->x[0],3*atom->nlocal,MPI_DOUBLE,ireplica,comm_replica);
} else {
int *counts = new int[nprocs];
if (universe->iworld == ireplica) {
MPI_Gather(&atom->nlocal,1,MPI_INT,counts,1,MPI_INT,0,world);
displacements[0] = 0;
for (i = 0; i < nprocs-1; i++)
displacements[i+1] = displacements[i] + counts[i];
MPI_Gatherv(atom->tag,atom->nlocal,MPI_LMP_TAGINT,
tagall,counts,displacements,MPI_LMP_TAGINT,0,world);
MPI_Gatherv(atom->image,atom->nlocal,MPI_INT,
imageall,counts,displacements,MPI_INT,0,world);
for (i = 0; i < nprocs; i++) counts[i] *= 3;
for (i = 0; i < nprocs-1; i++)
displacements[i+1] = displacements[i] + counts[i];
MPI_Gatherv(atom->x[0],3*atom->nlocal,MPI_DOUBLE,
xall[0],counts,displacements,MPI_DOUBLE,0,world);
}
if (me == 0) {
MPI_Bcast(tagall,natoms,MPI_INT,ireplica,comm_replica);
MPI_Bcast(imageall,natoms,MPI_INT,ireplica,comm_replica);
MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,ireplica,comm_replica);
}
MPI_Bcast(tagall,natoms,MPI_INT,0,world);
MPI_Bcast(imageall,natoms,MPI_INT,0,world);
MPI_Bcast(xall[0],3*natoms,MPI_DOUBLE,0,world);
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < natoms; i++) {
m = atom->map(tagall[i]);
if (m >= 0 && m < nlocal) {
x[m][0] = xall[i][0];
x[m][1] = xall[i][1];
x[m][2] = xall[i][2];
atom->image[m] = imageall[i];
}
}
delete [] counts;
}
}
/* ----------------------------------------------------------------------
parse optional parameters at end of PRD input line
------------------------------------------------------------------------- */
void PRD::options(int narg, char **arg)
{
if (narg < 0) error->all(FLERR,"Illegal prd command");
// set defaults
etol = 0.1;
ftol = 0.1;
maxiter = 40;
maxeval = 50;
temp_flag = 0;
stepmode = 0;
char *str = (char *) "geom";
int n = strlen(str) + 1;
loop_setting = new char[n];
strcpy(loop_setting,str);
str = (char *) "gaussian";
n = strlen(str) + 1;
dist_setting = new char[n];
strcpy(dist_setting,str);
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"min") == 0) {
if (iarg+5 > narg) error->all(FLERR,"Illegal prd command");
etol = force->numeric(FLERR,arg[iarg+1]);
ftol = force->numeric(FLERR,arg[iarg+2]);
maxiter = force->inumeric(FLERR,arg[iarg+3]);
maxeval = force->inumeric(FLERR,arg[iarg+4]);
if (maxiter < 0) error->all(FLERR,"Illegal prd command");
iarg += 5;
} else if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal prd command");
temp_flag = 1;
temp_dephase = force->numeric(FLERR,arg[iarg+1]);
if (temp_dephase <= 0.0) error->all(FLERR,"Illegal prd command");
iarg += 2;
} else if (strcmp(arg[iarg],"vel") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal prd command");
delete [] loop_setting;
delete [] dist_setting;
if (strcmp(arg[iarg+1],"all") == 0) loop_setting = NULL;
else if (strcmp(arg[iarg+1],"local") == 0) loop_setting = NULL;
else if (strcmp(arg[iarg+1],"geom") == 0) loop_setting = NULL;
else error->all(FLERR,"Illegal prd command");
int n = strlen(arg[iarg+1]) + 1;
loop_setting = new char[n];
strcpy(loop_setting,arg[iarg+1]);
if (strcmp(arg[iarg+2],"uniform") == 0) dist_setting = NULL;
else if (strcmp(arg[iarg+2],"gaussian") == 0) dist_setting = NULL;
else error->all(FLERR,"Illegal prd command");
n = strlen(arg[iarg+2]) + 1;
dist_setting = new char[n];
strcpy(dist_setting,arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"time") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal prd command");
if (strcmp(arg[iarg+1],"steps") == 0) stepmode = 0;
else if (strcmp(arg[iarg+1],"clock") == 0) stepmode = 1;
else error->all(FLERR,"Illegal prd command");
iarg += 2;
} else error->all(FLERR,"Illegal prd command");
}
}
diff --git a/src/REPLICA/tad.cpp b/src/REPLICA/tad.cpp
index dfd1ea7e5..80381a7a3 100644
--- a/src/REPLICA/tad.cpp
+++ b/src/REPLICA/tad.cpp
@@ -1,1039 +1,1044 @@
/* ----------------------------------------------------------------------
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: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
+// lmptype.h must be first b/c this file uses MAXBIGINT and includes mpi.h
+// due to OpenMPI bug which sets INT64_MAX via its mpi.h
+// before lmptype.h can set flags to insure it is done correctly
+
+#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "tad.h"
#include "universe.h"
#include "update.h"
#include "atom.h"
#include "domain.h"
#include "region.h"
#include "comm.h"
#include "velocity.h"
#include "integrate.h"
#include "min.h"
#include "neighbor.h"
#include "modify.h"
#include "neb.h"
#include "compute.h"
#include "fix.h"
#include "fix_event_tad.h"
#include "fix_store.h"
#include "force.h"
#include "pair.h"
#include "random_park.h"
#include "random_mars.h"
#include "output.h"
#include "dump.h"
#include "finish.h"
#include "timer.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
TAD::TAD(LAMMPS *lmp) : Pointers(lmp) {}
/* ---------------------------------------------------------------------- */
TAD::~TAD()
{
memory->sfree(fix_event_list);
if (neb_logfilename != NULL) delete [] neb_logfilename;
delete [] min_style;
delete [] min_style_neb;
}
/* ----------------------------------------------------------------------
perform TAD simulation on root proc
other procs only used for NEB calcs
------------------------------------------------------------------------- */
void TAD::command(int narg, char **arg)
{
fix_event_list = NULL;
n_event_list = 0;
nmax_event_list = 0;
nmin_event_list = 10;
// error checks
if (domain->box_exist == 0)
error->all(FLERR,"Tad command before simulation box is defined");
if (universe->nworlds == 1)
error->all(FLERR,"Cannot use TAD with a single replica for NEB");
if (universe->nworlds != universe->nprocs)
error->all(FLERR,"Can only use TAD with 1-processor replicas for NEB");
if (atom->sortfreq > 0)
error->all(FLERR,"Cannot use TAD with atom_modify sort enabled for NEB");
if (atom->map_style == 0)
error->all(FLERR,"Cannot use TAD unless atom map exists for NEB");
if (narg < 7) error->universe_all(FLERR,"Illegal tad command");
nsteps = force->inumeric(FLERR,arg[0]);
t_event = force->inumeric(FLERR,arg[1]);
templo = force->numeric(FLERR,arg[2]);
temphi = force->numeric(FLERR,arg[3]);
delta_conf = force->numeric(FLERR,arg[4]);
tmax = force->numeric(FLERR,arg[5]);
char *id_compute = new char[strlen(arg[6])+1];
strcpy(id_compute,arg[6]);
// quench minimizer is set by min_style command
// NEB minimizer is set by options, default = quickmin
int n = strlen(update->minimize_style) + 1;
min_style = new char[n];
strcpy(min_style,update->minimize_style);
options(narg-7,&arg[7]);
// total # of timesteps must be multiple of t_event
if (t_event <= 0) error->universe_all(FLERR,"Invalid t_event in tad command");
if (nsteps % t_event)
error->universe_all(FLERR,"TAD nsteps must be multiple of t_event");
if (delta_conf <= 0.0 || delta_conf >= 1.0)
error->universe_all(FLERR,"Invalid delta_conf in tad command");
if (tmax <= 0.0)
error->universe_all(FLERR,"Invalid tmax in tad command");
// deltconf = (ln(1/delta))/freq_min (timestep units)
deltconf = -log(delta_conf)*tmax/update->dt;
// local storage
int me_universe = universe->me;
int nprocs_universe = universe->nprocs;
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
delta_beta = (1.0/templo - 1.0/temphi) / force->boltz;
ratio_beta = templo/temphi;
// create FixEventTAD object to store last event
int narg2 = 3;
char **args = new char*[narg2];
args[0] = (char *) "tad_event";
args[1] = (char *) "all";
args[2] = (char *) "EVENT/TAD";
modify->add_fix(narg2,args);
fix_event = (FixEventTAD *) modify->fix[modify->nfix-1];
delete [] args;
// create FixStore object to store revert state
narg2 = 5;
args = new char*[narg2];
args[0] = (char *) "tad_revert";
args[1] = (char *) "all";
args[2] = (char *) "STORE";
args[3] = (char *) "0";
args[4] = (char *) "7";
modify->add_fix(narg2,args);
fix_revert = (FixStore *) modify->fix[modify->nfix-1];
delete [] args;
// create Finish for timing output
finish = new Finish(lmp);
// assign FixEventTAD to event-detection compute
// necessary so it will know atom coords at last event
int icompute = modify->find_compute(id_compute);
if (icompute < 0) error->all(FLERR,"Could not find compute ID for TAD");
compute_event = modify->compute[icompute];
compute_event->reset_extra_compute_fix("tad_event");
// reset reneighboring criteria since will perform minimizations
neigh_every = neighbor->every;
neigh_delay = neighbor->delay;
neigh_dist_check = neighbor->dist_check;
if (neigh_every != 1 || neigh_delay != 0 || neigh_dist_check != 1) {
if (me_universe == 0)
error->warning(FLERR,"Resetting reneighboring criteria during TAD");
}
neighbor->every = 1;
neighbor->delay = 0;
neighbor->dist_check = 1;
// initialize TAD as if one long dynamics run
update->whichflag = 1;
update->nsteps = nsteps;
update->beginstep = update->firststep = update->ntimestep;
update->endstep = update->laststep = update->firststep + nsteps;
update->restrict_output = 1;
if (update->laststep < 0 || update->laststep > MAXBIGINT)
error->all(FLERR,"Too many timesteps");
lmp->init();
// set minimize style for quench
narg2 = 1;
args = new char*[narg2];
args[0] = min_style;
update->create_minimize(narg2,args);
delete [] args;
// init minimizer settings and minimizer itself
update->etol = etol;
update->ftol = ftol;
update->max_eval = maxeval;
update->minimize->init();
// perform TAD simulation
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Setting up TAD ...\n");
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,
"Step CPU N M Status Barrier Margin t_lo delt_lo\n");
if (universe->ulogfile)
fprintf(universe->ulogfile,
"Step CPU N M Status Barrier Margin t_lo delt_lo\n");
}
ulogfile_lammps = universe->ulogfile;
uscreen_lammps = universe->uscreen;
ulogfile_neb = NULL;
uscreen_neb = NULL;
if (me_universe == 0 && neb_logfilename)
ulogfile_neb = fopen(neb_logfilename,"w");
// store hot state and quenched event, only on replica 0
// need this line if quench() does only setup_minimal()
// update->minimize->setup();
// This should work with if uncommented, but does not
// if (universe->iworld == 0) {
fix_event->store_state_quench();
quench();
timer->init();
timer->barrier_start();
time_start = timer->get_wall(Timer::TOTAL);
fix_event->store_event_tad(update->ntimestep);
log_event(0);
fix_event->restore_state_quench();
// do full init/setup
update->whichflag = 1;
lmp->init();
update->integrate->setup();
// main loop: look for events until out of time
// (1) dynamics, store state, quench, check event, restore state
// (2) if event, perform NEB, record in fix_event_list
// (3) if confident, pick earliest event
nbuild = ndanger = 0;
time_neb = time_dynamics = time_quench = time_comm = time_output = 0.0;
timer->barrier_start();
time_start = timer->get_wall(Timer::TOTAL);
int confident_flag, event_flag;
if (universe->iworld == 0) {
while (update->ntimestep < update->endstep) {
// initialize list of possible events
initialize_event_list();
confident_flag = 0;
while (update->ntimestep < update->endstep) {
event_flag = 0;
while (update->ntimestep < update->endstep) {
dynamics();
fix_event->store_state_quench();
quench();
event_flag = check_event();
MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld);
if (event_flag) break;
// restore hot state
fix_event->restore_state_quench();
// store hot state in revert
store_state();
}
if (!event_flag) break;
add_event();
perform_neb(n_event_list-1);
compute_tlo(n_event_list-1);
confident_flag = check_confidence();
MPI_Bcast(&confident_flag,1,MPI_INT,0,universe->uworld);
if (confident_flag) break;
if (universe->iworld == 0) revert_state();
}
if (!confident_flag) break;
perform_event(event_first);
// need to sync timestep with TAD
MPI_Bcast(&(update->ntimestep),1,MPI_INT,0,universe->uworld);
int restart_flag = 0;
if (output->restart_flag && universe->iworld == 0) {
if (output->restart_every_single &&
fix_event->event_number % output->restart_every_single == 0)
restart_flag = 1;
if (output->restart_every_double &&
fix_event->event_number % output->restart_every_double == 0)
restart_flag = 1;
}
// full init/setup since are starting after event
update->whichflag = 1;
lmp->init();
update->integrate->setup();
// write restart file of hot coords
if (restart_flag) {
timer->barrier_start();
output->write_restart(update->ntimestep);
timer->barrier_stop();
time_output += timer->get_wall(Timer::TOTAL);
}
}
} else {
while (update->ntimestep < update->endstep) {
confident_flag = 0;
while (update->ntimestep < update->endstep) {
event_flag = 0;
while (update->ntimestep < update->endstep) {
update->ntimestep += t_event;
MPI_Bcast(&event_flag,1,MPI_INT,0,universe->uworld);
if (event_flag) break;
}
if (!event_flag) break;
perform_neb(-1);
MPI_Bcast(&confident_flag,1,MPI_INT,0,universe->uworld);
if (confident_flag) break;
}
if (!confident_flag) break;
// need to sync timestep with TAD
MPI_Bcast(&(update->ntimestep),1,MPI_INT,0,universe->uworld);
}
}
// set total timers and counters so Finish() will process them
timer->set_wall(Timer::TOTAL, time_start);
timer->barrier_stop();
timer->set_wall(Timer::NEB, time_neb);
timer->set_wall(Timer::DYNAMICS, time_dynamics);
timer->set_wall(Timer::QUENCH, time_quench);
timer->set_wall(Timer::REPCOMM, time_comm);
timer->set_wall(Timer::REPOUT, time_output);
neighbor->ncalls = nbuild;
neighbor->ndanger = ndanger;
if (me_universe == 0) {
if (universe->uscreen)
fprintf(universe->uscreen,
"Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT
" atoms\n",
timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms);
if (universe->ulogfile)
fprintf(universe->ulogfile,
"Loop time of %g on %d procs for %d steps with " BIGINT_FORMAT
" atoms\n",
timer->get_wall(Timer::TOTAL),nprocs_universe,nsteps,atom->natoms);
}
if (me_universe == 0) fclose(ulogfile_neb);
if (me == 0) {
if (screen) fprintf(screen,"\nTAD done\n");
if (logfile) fprintf(logfile,"\nTAD done\n");
}
finish->end(3);
update->whichflag = 0;
update->firststep = update->laststep = 0;
update->beginstep = update->endstep = 0;
update->restrict_output = 0;
// reset reneighboring criteria
neighbor->every = neigh_every;
neighbor->delay = neigh_delay;
neighbor->dist_check = neigh_dist_check;
delete [] id_compute;
delete finish;
modify->delete_fix("tad_event");
modify->delete_fix("tad_revert");
delete_event_list();
compute_event->reset_extra_compute_fix(NULL);
}
/* ----------------------------------------------------------------------
single short dynamics run
------------------------------------------------------------------------- */
void TAD::dynamics()
{
update->whichflag = 1;
update->nsteps = t_event;
lmp->init();
update->integrate->setup();
// this may be needed if don't do full init
//modify->addstep_compute_all(update->ntimestep);
int ncalls = neighbor->ncalls;
timer->barrier_start();
update->integrate->run(t_event);
timer->barrier_stop();
time_dynamics += timer->get_wall(Timer::TOTAL);
nbuild += neighbor->ncalls - ncalls;
ndanger += neighbor->ndanger;
update->integrate->cleanup();
finish->end(0);
}
/* ----------------------------------------------------------------------
quench minimization
------------------------------------------------------------------------- */
void TAD::quench()
{
bigint ntimestep_hold = update->ntimestep;
bigint endstep_hold = update->endstep;
// need to change whichflag so that minimize->setup() calling
// modify->setup() will call fix->min_setup()
update->whichflag = 2;
update->nsteps = maxiter;
update->endstep = update->laststep = update->firststep + maxiter;
if (update->laststep < 0 || update->laststep > MAXBIGINT)
error->all(FLERR,"Too many iterations");
// full init works
lmp->init();
update->minimize->setup();
// partial init does not work
//modify->addstep_compute_all(update->ntimestep);
//update->minimize->setup_minimal(1);
int ncalls = neighbor->ncalls;
timer->barrier_start();
update->minimize->run(maxiter);
timer->barrier_stop();
time_quench += timer->get_wall(Timer::TOTAL);
if (neighbor->ncalls == ncalls) quench_reneighbor = 0;
else quench_reneighbor = 1;
update->minimize->cleanup();
finish->end(1);
// reset timestep as if quench did not occur
// clear timestep storage from computes, since now invalid
update->ntimestep = ntimestep_hold;
update->endstep = update->laststep = endstep_hold;
for (int i = 0; i < modify->ncompute; i++)
if (modify->compute[i]->timeflag) modify->compute[i]->clearstep();
}
/* ----------------------------------------------------------------------
check for an event
return 0 if no event
return 1 if event
------------------------------------------------------------------------- */
int TAD::check_event()
{
int flag;
flag = 0;
if (compute_event->compute_scalar() > 0.0) flag = 1;
return flag;
}
/* ----------------------------------------------------------------------
universe proc 0 prints event info
------------------------------------------------------------------------- */
void TAD::log_event(int ievent)
{
timer->set_wall(Timer::TOTAL, time_start);
if (universe->me == 0) {
double tfrac = 0.0;
if (universe->uscreen)
fprintf(universe->uscreen,
BIGINT_FORMAT " %.3f %d %d %s %.3f %.3f %.3f %.3f\n",
fix_event->event_timestep,
timer->elapsed(Timer::TOTAL),
fix_event->event_number,ievent,
"E ",
fix_event->ebarrier,tfrac,
fix_event->tlo,deltfirst);
if (universe->ulogfile)
fprintf(universe->ulogfile,
BIGINT_FORMAT " %.3f %d %d %s %.3f %.3f %.3f %.3f\n",
fix_event->event_timestep,
timer->elapsed(Timer::TOTAL),
fix_event->event_number,ievent,
"E ",
fix_event->ebarrier,tfrac,
fix_event->tlo,deltfirst);
}
// dump snapshot of quenched coords
// must reneighbor and compute forces before dumping
// addstep_compute_all insures eng/virial are calculated if needed
if (output->ndump && universe->iworld == 0) {
timer->barrier_start();
modify->addstep_compute_all(update->ntimestep);
update->integrate->setup_minimal(1);
output->write_dump(update->ntimestep);
timer->barrier_stop();
time_output += timer->get_wall(Timer::TOTAL);
}
}
/* ----------------------------------------------------------------------
parse optional parameters at end of TAD input line
------------------------------------------------------------------------- */
void TAD::options(int narg, char **arg)
{
if (narg < 0) error->all(FLERR,"Illegal tad command");
// set defaults
etol = 0.1;
ftol = 0.1;
maxiter = 40;
maxeval = 50;
etol_neb = 0.01;
ftol_neb = 0.01;
n1steps_neb = 100;
n2steps_neb = 100;
nevery_neb = 10;
int n = strlen("quickmin") + 1;
min_style_neb = new char[n];
strcpy(min_style_neb,"quickmin");
dt_neb = update->dt;
neb_logfilename = NULL;
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"min") == 0) {
if (iarg+5 > narg) error->all(FLERR,"Illegal tad command");
etol = force->numeric(FLERR,arg[iarg+1]);
ftol = force->numeric(FLERR,arg[iarg+2]);
maxiter = force->inumeric(FLERR,arg[iarg+3]);
maxeval = force->inumeric(FLERR,arg[iarg+4]);
if (maxiter < 0 || maxeval < 0 ||
etol < 0.0 || ftol < 0.0 )
error->all(FLERR,"Illegal tad command");
iarg += 5;
} else if (strcmp(arg[iarg],"neb") == 0) {
if (iarg+6 > narg) error->all(FLERR,"Illegal tad command");
etol_neb = force->numeric(FLERR,arg[iarg+1]);
ftol_neb = force->numeric(FLERR,arg[iarg+2]);
n1steps_neb = force->inumeric(FLERR,arg[iarg+3]);
n2steps_neb = force->inumeric(FLERR,arg[iarg+4]);
nevery_neb = force->inumeric(FLERR,arg[iarg+5]);
if (etol_neb < 0.0 || ftol_neb < 0.0 ||
n1steps_neb < 0 || n2steps_neb < 0 ||
nevery_neb < 0) error->all(FLERR,"Illegal tad command");
iarg += 6;
} else if (strcmp(arg[iarg],"neb_style") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal tad command");
delete [] min_style_neb;
int n = strlen(arg[iarg+1]) + 1;
min_style_neb = new char[n];
strcpy(min_style_neb,arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"neb_step") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal tad command");
dt_neb = force->numeric(FLERR,arg[iarg+1]);
if (dt_neb <= 0.0) error->all(FLERR,"Illegal tad command");
iarg += 2;
} else if (strcmp(arg[iarg],"neb_log") == 0) {
delete [] neb_logfilename;
if (iarg+2 > narg) error->all(FLERR,"Illegal tad command");
if (strcmp(arg[iarg+1],"none") == 0) neb_logfilename = NULL;
else {
int n = strlen(arg[iarg+1]) + 1;
neb_logfilename = new char[n];
strcpy(neb_logfilename,arg[iarg+1]);
}
iarg += 2;
} else error->all(FLERR,"Illegal tad command");
}
}
/* ----------------------------------------------------------------------
perform NEB calculation
------------------------------------------------------------------------- */
void TAD::perform_neb(int ievent)
{
double **x = atom->x;
int nlocal = atom->nlocal;
double *buf_final;
memory->create(buf_final,3*nlocal,"tad:buffinal");
// set system to quenched state of event ievent
if (universe->iworld == 0) {
fix_event_list[ievent]->restore_event();
int ii = 0;
for (int i = 0; i < nlocal; i++) {
buf_final[ii++] = x[i][0];
buf_final[ii++] = x[i][1];
buf_final[ii++] = x[i][2];
}
}
MPI_Bcast(buf_final,3*nlocal,MPI_DOUBLE,universe->root_proc[0],
universe->uworld);
double *buf_init;
memory->create(buf_init,3*nlocal,"tad:bufinit");
// set system to quenched state of fix_event
if (universe->iworld == 0) {
fix_event->restore_event();
int ii = 0;
for (int i = 0; i < nlocal; i++) {
buf_init[ii++] = x[i][0];
buf_init[ii++] = x[i][1];
buf_init[ii++] = x[i][2];
}
}
MPI_Bcast(buf_init,3*nlocal,MPI_DOUBLE,universe->root_proc[0],
universe->uworld);
// create FixNEB object to support NEB
int narg2 = 4;
char **args = new char*[narg2];
args[0] = (char *) "neb";
args[1] = (char *) "all";
args[2] = (char *) "neb";
char str[128];
args[3] = str;
double kspring = 1.0;
sprintf(args[3],"%f",kspring);
modify->add_fix(narg2,args);
fix_neb = (Fix *) modify->fix[modify->nfix-1];
delete [] args;
// switch minimize style to quickmin for NEB
narg2 = 1;
args = new char*[narg2];
args[0] = min_style_neb;
update->create_minimize(narg2,args);
delete [] args;
// create NEB object
neb = new NEB(lmp,etol_neb,ftol_neb,n1steps_neb,
n2steps_neb,nevery_neb,buf_init,buf_final);
// free up temporary arrays
memory->destroy(buf_init);
memory->destroy(buf_final);
// run NEB with NEB timestep
int beginstep_hold = update->beginstep;
int endstep_hold = update->endstep;
int ntimestep_hold = update->ntimestep;
int nsteps_hold = update->nsteps;
if (universe->me == 0) {
universe->ulogfile = ulogfile_neb;
universe->uscreen = uscreen_neb;
}
// had to bypass timer interface
// because timer->array is reset inside neb->run()
// timer->barrier_start();
// neb->run();
// timer->barrier_stop();
// time_neb += timer->get_wall(Timer::TOTAL);
MPI_Barrier(world);
double time_tmp = MPI_Wtime();
double dt_hold = update->dt;
update->dt = dt_neb;
neb->run();
update->dt = dt_hold;
MPI_Barrier(world);
time_neb += MPI_Wtime() - time_tmp;
if (universe->me == 0) {
universe->ulogfile = ulogfile_lammps;
universe->uscreen = uscreen_lammps;
}
// extract barrier energy from NEB
if (universe->iworld == 0)
fix_event_list[ievent]->ebarrier = neb->ebf;
update->beginstep = update->firststep = beginstep_hold;
update->endstep = update->laststep = endstep_hold;
update->ntimestep = ntimestep_hold;
update->nsteps = nsteps_hold;
// switch minimize style back for quench
narg2 = 1;
args = new char*[narg2];
args[0] = min_style;
update->create_minimize(narg2,args);
update->etol = etol;
update->ftol = ftol;
delete [] args;
// clean up
modify->delete_fix("neb");
delete neb;
}
/* ----------------------------------------------------------------------
check if confidence criterion for tstop is satisfied
return 0 if not satisfied
return 1 if satisfied
------------------------------------------------------------------------- */
int TAD::check_confidence()
{
int flag;
// update stopping time
deltstop = deltconf*pow(deltfirst/deltconf, ratio_beta);
flag = 0;
if (deltstop < update->ntimestep - fix_event->event_timestep) flag = 1;
return flag;
}
/* ----------------------------------------------------------------------
store state in fix_revert
------------------------------------------------------------------------- */
void TAD::store_state()
{
double **x = atom->x;
double **v = atom->v;
imageint *image = atom->image;
int nlocal = atom->nlocal;
double **astore = fix_revert->astore;
for (int i = 0; i < nlocal; i++) {
astore[i][0] = x[i][0];
astore[i][1] = x[i][1];
astore[i][2] = x[i][2];
astore[i][3] = v[i][0];
astore[i][4] = v[i][1];
astore[i][5] = v[i][2];
*((imageint *) &astore[i][6]) = image[i];
}
}
/* ----------------------------------------------------------------------
restore state archived in fix_revert
flip sign of velocities to reflect back to starting state
------------------------------------------------------------------------- */
void TAD::revert_state()
{
double **x = atom->x;
double **v = atom->v;
imageint *image = atom->image;
int nlocal = atom->nlocal;
double **astore = fix_revert->astore;
for (int i = 0; i < nlocal; i++) {
x[i][0] = astore[i][0];
x[i][1] = astore[i][1];
x[i][2] = astore[i][2];
v[i][0] = -astore[i][3];
v[i][1] = -astore[i][4];
v[i][2] = -astore[i][5];
image[i] = *((imageint *) &astore[i][6]);
}
}
/* ----------------------------------------------------------------------
Initialize list of possible events
------------------------------------------------------------------------- */
void TAD::initialize_event_list() {
// First delete old events, if any
delete_event_list();
// Create new list
n_event_list = 0;
grow_event_list(nmin_event_list);
}
/* ----------------------------------------------------------------------
Delete list of possible events
------------------------------------------------------------------------- */
void TAD::delete_event_list() {
for (int i = 0; i < n_event_list; i++) {
char str[128];
sprintf(str,"tad_event_%d",i);
modify->delete_fix(str);
}
memory->sfree(fix_event_list);
fix_event_list = NULL;
n_event_list = 0;
nmax_event_list = 0;
}
/* ----------------------------------------------------------------------
add event
------------------------------------------------------------------------- */
void TAD::add_event()
{
// create FixEventTAD object to store possible event
int narg = 3;
char **args = new char*[narg];
char str[128];
sprintf(str,"tad_event_%d",n_event_list);
args[0] = str;
args[1] = (char *) "all";
args[2] = (char *) "EVENT/TAD";
modify->add_fix(narg,args);
if (n_event_list == nmax_event_list)
grow_event_list(nmax_event_list+nmin_event_list);
n_event_list += 1;
int ievent = n_event_list-1;
fix_event_list[ievent] = (FixEventTAD *) modify->fix[modify->nfix-1];
// store quenched state for new event
fix_event_list[ievent]->store_event_tad(update->ntimestep);
// store hot state for new event
fix_event->restore_state_quench();
fix_event_list[ievent]->store_state_quench();
// string clean-up
delete [] args;
}
/* ----------------------------------------------------------------------
compute cold time for event ievent
------------------------------------------------------------------------- */
void TAD::compute_tlo(int ievent)
{
double deltlo,delthi,ebarrier;
ebarrier = fix_event_list[ievent]->ebarrier;
delthi = fix_event_list[ievent]->event_timestep
- fix_event->event_timestep;
deltlo = delthi*exp(ebarrier*delta_beta);
fix_event_list[ievent]->tlo = fix_event->tlo + deltlo;
// update first event
char* statstr = (char *) "D ";
if (ievent == 0) {
deltfirst = deltlo;
event_first = ievent;
statstr = (char *) "DF";
} else if (deltlo < deltfirst) {
deltfirst = deltlo;
event_first = ievent;
statstr = (char *) "DF";
}
// first-replica output about each event
timer->set_wall(Timer::TOTAL, time_start);
if (universe->me == 0) {
double tfrac = 0.0;
if (ievent > 0) tfrac = delthi/deltstop;
if (universe->uscreen)
fprintf(universe->uscreen,
BIGINT_FORMAT " %.3f %d %d %s %.3f %.3f %.3f %.3f\n",
fix_event_list[ievent]->event_timestep,
timer->elapsed(Timer::TOTAL),
fix_event->event_number,
ievent,statstr,ebarrier,tfrac,
fix_event->tlo,deltlo);
if (universe->ulogfile)
fprintf(universe->ulogfile,
BIGINT_FORMAT " %.3f %d %d %s %.3f %.3f %.3f %.3f\n",
fix_event_list[ievent]->event_timestep,
timer->elapsed(Timer::TOTAL),
fix_event->event_number,
ievent,statstr,ebarrier,tfrac,
fix_event->tlo,deltlo);
}
}
/* ----------------------------------------------------------------------
perform event
------------------------------------------------------------------------- */
void TAD::perform_event(int ievent)
{
// reset timestep to that of event
update->ntimestep = fix_event_list[ievent]->event_timestep;
// Copy event to current event
// Should really use copy constructor for this
fix_event->tlo = fix_event_list[ievent]->tlo;
fix_event->ebarrier = fix_event_list[ievent]->ebarrier;
fix_event->event_number++;
fix_event->event_timestep = update->ntimestep;
fix_event_list[ievent]->restore_event();
fix_event->store_event_tad(fix_event_list[ievent]->event_timestep);
// output stats and dump for quenched state
log_event(ievent);
// load and store hot state
fix_event_list[ievent]->restore_state_quench();
fix_event->store_state_quench();
}
/* ----------------------------------------------------------------------
Allocate list of pointers to events
------------------------------------------------------------------------- */
void TAD::grow_event_list(int nmax) {
if (nmax_event_list > nmax) return;
fix_event_list = (FixEventTAD **)
memory->srealloc(fix_event_list,nmax*sizeof(FixEventTAD *),"tad:eventlist");
nmax_event_list = nmax;
}
diff --git a/src/read_data.cpp b/src/read_data.cpp
index 43c290a8b..cf87b2773 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -1,1662 +1,1667 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
+// lmptype.h must be first b/c this file uses MAXBIGINT and includes mpi.h
+// due to OpenMPI bug which sets INT64_MAX via its mpi.h
+// before lmptype.h can set flags to insure it is done correctly
+
+#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "string.h"
#include "stdlib.h"
#include "ctype.h"
#include "read_data.h"
#include "atom.h"
#include "atom_vec.h"
#include "atom_vec_ellipsoid.h"
#include "atom_vec_line.h"
#include "atom_vec_tri.h"
#include "force.h"
#include "molecule.h"
#include "comm.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "special.h"
#include "error.h"
#include "memory.h"
using namespace LAMMPS_NS;
#define MAXLINE 256
#define LB_FACTOR 1.1
#define CHUNK 1024
#define DELTA 4 // must be 2 or larger
#define MAXBODY 20 // max # of lines in one body, also in Atom class
// customize for new sections
#define NSECTIONS 25 // change when add to header::section_keywords
// pair style suffixes to ignore
// when matching Pair Coeffs comment to currently-defined pair style
const char *suffixes[] = {"/cuda","/gpu","/opt","/omp","/kk",
"/coul/cut","/coul/long","/coul/msm",
"/coul/dsf","/coul/debye","/coul/charmm",
NULL};
/* ---------------------------------------------------------------------- */
ReadData::ReadData(LAMMPS *lmp) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
line = new char[MAXLINE];
keyword = new char[MAXLINE];
style = new char[MAXLINE];
buffer = new char[CHUNK*MAXLINE];
narg = maxarg = 0;
arg = NULL;
// customize for new sections
// pointers to atom styles that store extra info
nellipsoids = 0;
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
nlines = 0;
avec_line = (AtomVecLine *) atom->style_match("line");
ntris = 0;
avec_tri = (AtomVecTri *) atom->style_match("tri");
nbodies = 0;
avec_body = (AtomVecBody *) atom->style_match("body");
}
/* ---------------------------------------------------------------------- */
ReadData::~ReadData()
{
delete [] line;
delete [] keyword;
delete [] style;
delete [] buffer;
memory->sfree(arg);
for (int i = 0; i < nfix; i++) {
delete [] fix_header[i];
delete [] fix_section[i];
}
memory->destroy(fix_index);
memory->sfree(fix_header);
memory->sfree(fix_section);
}
/* ---------------------------------------------------------------------- */
void ReadData::command(int narg, char **arg)
{
if (narg < 1) error->all(FLERR,"Illegal read_data command");
// optional args
addflag = mergeflag = 0;
offset[0] = offset[1] = offset[2] = 0.0;
nfix = 0;
fix_index = NULL;
fix_header = NULL;
fix_section = NULL;
int iarg = 1;
while (iarg < narg) {
if (strcmp(arg[iarg],"add") == 0) {
addflag = 1;
iarg++;
} else if (strcmp(arg[iarg],"merge") == 0) {
mergeflag = 1;
iarg++;
} else if (strcmp(arg[iarg],"offset") == 0) {
if (iarg+4 > narg)
error->all(FLERR,"Illegal read_data command");
offset[0] = force->numeric(FLERR,arg[iarg+1]);
offset[1] = force->numeric(FLERR,arg[iarg+2]);
offset[2] = force->numeric(FLERR,arg[iarg+3]);
iarg += 4;
} else if (strcmp(arg[iarg],"fix") == 0) {
if (iarg+4 > narg)
error->all(FLERR,"Illegal read_data command");
memory->grow(fix_index,nfix+1,"read_data:fix_index");
fix_header = (char **)
memory->srealloc(fix_header,(nfix+1)*sizeof(char *),
"read_data:fix_header");
fix_section = (char **)
memory->srealloc(fix_section,(nfix+1)*sizeof(char *),
"read_data:fix_section");
fix_index[nfix] = modify->find_fix(arg[iarg+1]);
if (fix_index[nfix] < 0)
error->all(FLERR,"Fix ID for read_data does not exist");
if (strcmp(arg[iarg+2],"NULL") == 0) fix_header[nfix] = NULL;
else {
int n = strlen(arg[iarg+2]) + 1;
fix_header[nfix] = new char[n];
strcpy(fix_header[nfix],arg[iarg+2]);
}
int n = strlen(arg[iarg+3]) + 1;
fix_section[nfix] = new char[n];
strcpy(fix_section[nfix],arg[iarg+3]);
nfix++;
iarg += 4;
} else error->all(FLERR,"Illegal read_data command");
}
// error checks
if (domain->box_exist && !addflag && !mergeflag)
error->all(FLERR,"Cannot read_data after simulation box is defined");
if (addflag && mergeflag) error->all(FLERR,"Cannot read_data add and merge");
if (domain->dimension == 2 && offset[2] != 0.0)
error->all(FLERR,"Cannot use non-zero z offset in read_data "
"for 2d simulation");
if (domain->dimension == 2 && domain->zperiodic == 0)
error->all(FLERR,"Cannot run 2d simulation with nonperiodic Z dimension");
// perform 1-pass read if no molecular topoogy in file
// perform 2-pass read if molecular topology,
// first pass calculates max topology/atom
int atomflag,topoflag;
int bondflag,angleflag,dihedralflag,improperflag;
int ellipsoidflag,lineflag,triflag,bodyflag;
atomflag = topoflag = 0;
bondflag = angleflag = dihedralflag = improperflag = 0;
ellipsoidflag = lineflag = triflag = bodyflag = 0;
int firstpass = 1;
while (1) {
// open file on proc 0
if (me == 0) {
if (firstpass && screen) fprintf(screen,"Reading data file ...\n");
open(arg[0]);
} else fp = NULL;
// read header info
header();
// problem setup using info from header
// 1st pass only
if (firstpass) {
domain->box_exist = 1;
update->ntimestep = 0;
// apply extra settings before grow(), even if no topology in file
// deallocate() insures new settings are used for topology arrays
// if per-atom topology is in file, another grow() is done below
atom->bond_per_atom = atom->extra_bond_per_atom;
atom->angle_per_atom = atom->extra_angle_per_atom;
atom->dihedral_per_atom = atom->extra_dihedral_per_atom;
atom->improper_per_atom = atom->extra_improper_per_atom;
int n;
if (comm->nprocs == 1) n = static_cast<int> (atom->natoms);
else n = static_cast<int> (LB_FACTOR * atom->natoms / comm->nprocs);
atom->allocate_type_arrays();
atom->deallocate_topology();
atom->avec->grow(n);
domain->print_box(" ");
domain->set_initial_box();
domain->set_global_box();
comm->set_proc_grid();
domain->set_local_box();
}
// customize for new sections
// read rest of file in free format
while (strlen(keyword)) {
// if special fix matches, it processes section
if (nfix) {
int i;
for (i = 0; i < nfix; i++)
if (strcmp(keyword,fix_section[i]) == 0) {
if (firstpass) fix(fix_index[i],keyword);
else skip_lines(modify->fix[fix_index[i]]->
read_data_skip_lines(keyword));
parse_keyword(0);
break;
}
if (i < nfix) continue;
}
if (strcmp(keyword,"Atoms") == 0) {
atomflag = 1;
if (firstpass) {
if (me == 0 && !style_match(style,atom->atom_style))
error->warning(FLERR,"Atom style in data file differs "
"from currently defined atom style");
atoms();
} else skip_lines(atom->natoms);
} else if (strcmp(keyword,"Velocities") == 0) {
if (atomflag == 0)
error->all(FLERR,"Must read Atoms before Velocities");
if (firstpass) velocities();
else skip_lines(atom->natoms);
} else if (strcmp(keyword,"Bonds") == 0) {
topoflag = bondflag = 1;
if (atom->nbonds == 0)
error->all(FLERR,"Invalid data file section: Bonds");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Bonds");
bonds(firstpass);
} else if (strcmp(keyword,"Angles") == 0) {
topoflag = angleflag = 1;
if (atom->nangles == 0)
error->all(FLERR,"Invalid data file section: Angles");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Angles");
angles(firstpass);
} else if (strcmp(keyword,"Dihedrals") == 0) {
topoflag = dihedralflag = 1;
if (atom->ndihedrals == 0)
error->all(FLERR,"Invalid data file section: Dihedrals");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Dihedrals");
dihedrals(firstpass);
} else if (strcmp(keyword,"Impropers") == 0) {
topoflag = improperflag = 1;
if (atom->nimpropers == 0)
error->all(FLERR,"Invalid data file section: Impropers");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Impropers");
impropers(firstpass);
} else if (strcmp(keyword,"Ellipsoids") == 0) {
ellipsoidflag = 1;
if (!avec_ellipsoid)
error->all(FLERR,"Invalid data file section: Ellipsoids");
if (atomflag == 0)
error->all(FLERR,"Must read Atoms before Ellipsoids");
if (firstpass)
bonus(nellipsoids,(AtomVec *) avec_ellipsoid,"ellipsoids");
else skip_lines(nellipsoids);
} else if (strcmp(keyword,"Lines") == 0) {
lineflag = 1;
if (!avec_line)
error->all(FLERR,"Invalid data file section: Lines");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Lines");
if (firstpass) bonus(nlines,(AtomVec *) avec_line,"lines");
else skip_lines(nlines);
} else if (strcmp(keyword,"Triangles") == 0) {
triflag = 1;
if (!avec_tri)
error->all(FLERR,"Invalid data file section: Triangles");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Triangles");
if (firstpass) bonus(ntris,(AtomVec *) avec_tri,"triangles");
else skip_lines(ntris);
} else if (strcmp(keyword,"Bodies") == 0) {
bodyflag = 1;
if (!avec_body)
error->all(FLERR,"Invalid data file section: Bodies");
if (atomflag == 0) error->all(FLERR,"Must read Atoms before Bodies");
bodies(firstpass);
} else if (strcmp(keyword,"Masses") == 0) {
if (firstpass) mass();
else skip_lines(atom->ntypes);
} else if (strcmp(keyword,"Pair Coeffs") == 0) {
if (force->pair == NULL)
error->all(FLERR,"Must define pair_style before Pair Coeffs");
if (firstpass) {
if (me == 0 && !style_match(style,force->pair_style))
error->warning(FLERR,"Pair style in data file differs "
"from currently defined pair style");
paircoeffs();
} else skip_lines(atom->ntypes);
} else if (strcmp(keyword,"PairIJ Coeffs") == 0) {
if (force->pair == NULL)
error->all(FLERR,"Must define pair_style before PairIJ Coeffs");
if (firstpass) {
if (me == 0 && !style_match(style,force->pair_style))
error->warning(FLERR,"Pair style in data file differs "
"from currently defined pair style");
pairIJcoeffs();
} else skip_lines(atom->ntypes*(atom->ntypes+1)/2);
} else if (strcmp(keyword,"Bond Coeffs") == 0) {
if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Invalid data file section: Bond Coeffs");
if (force->bond == NULL)
error->all(FLERR,"Must define bond_style before Bond Coeffs");
if (firstpass) {
if (me == 0 && !style_match(style,force->bond_style))
error->warning(FLERR,"Bond style in data file differs "
"from currently defined bond style");
bondcoeffs();
} else skip_lines(atom->nbondtypes);
} else if (strcmp(keyword,"Angle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Invalid data file section: Angle Coeffs");
if (force->angle == NULL)
error->all(FLERR,"Must define angle_style before Angle Coeffs");
if (firstpass) {
if (me == 0 && !style_match(style,force->angle_style))
error->warning(FLERR,"Angle style in data file differs "
"from currently defined angle style");
anglecoeffs(0);
} else skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"Dihedral Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: Dihedral Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,"Must define dihedral_style before Dihedral Coeffs");
if (firstpass) {
if (me == 0 && !style_match(style,force->dihedral_style))
error->warning(FLERR,"Dihedral style in data file differs "
"from currently defined dihedral style");
dihedralcoeffs(0);
} else skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"Improper Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
error->all(FLERR,"Invalid data file section: Improper Coeffs");
if (force->improper == NULL)
error->all(FLERR,"Must define improper_style before Improper Coeffs");
if (firstpass) {
if (me == 0 && !style_match(style,force->improper_style))
error->warning(FLERR,"Improper style in data file differs "
"from currently defined improper style");
impropercoeffs(0);
} else skip_lines(atom->nimpropertypes);
} else if (strcmp(keyword,"BondBond Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Invalid data file section: BondBond Coeffs");
if (force->angle == NULL)
error->all(FLERR,"Must define angle_style before BondBond Coeffs");
if (firstpass) anglecoeffs(1);
else skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"BondAngle Coeffs") == 0) {
if (atom->avec->angles_allow == 0)
error->all(FLERR,"Invalid data file section: BondAngle Coeffs");
if (force->angle == NULL)
error->all(FLERR,"Must define angle_style before BondAngle Coeffs");
if (firstpass) anglecoeffs(2);
else skip_lines(atom->nangletypes);
} else if (strcmp(keyword,"MiddleBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,
"Invalid data file section: MiddleBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before "
"MiddleBondTorsion Coeffs");
if (firstpass) dihedralcoeffs(1);
else skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"EndBondTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: EndBondTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before EndBondTorsion Coeffs");
if (firstpass) dihedralcoeffs(2);
else skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: AngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before AngleTorsion Coeffs");
if (firstpass) dihedralcoeffs(3);
else skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleAngleTorsion Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,
"Invalid data file section: AngleAngleTorsion Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before "
"AngleAngleTorsion Coeffs");
if (firstpass) dihedralcoeffs(4);
else skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"BondBond13 Coeffs") == 0) {
if (atom->avec->dihedrals_allow == 0)
error->all(FLERR,"Invalid data file section: BondBond13 Coeffs");
if (force->dihedral == NULL)
error->all(FLERR,
"Must define dihedral_style before BondBond13 Coeffs");
if (firstpass) dihedralcoeffs(5);
else skip_lines(atom->ndihedraltypes);
} else if (strcmp(keyword,"AngleAngle Coeffs") == 0) {
if (atom->avec->impropers_allow == 0)
error->all(FLERR,"Invalid data file section: AngleAngle Coeffs");
if (force->improper == NULL)
error->all(FLERR,
"Must define improper_style before AngleAngle Coeffs");
if (firstpass) impropercoeffs(1);
else skip_lines(atom->nimpropertypes);
} else {
char str[128];
sprintf(str,"Unknown identifier in data file: %s",keyword);
error->all(FLERR,str);
}
parse_keyword(0);
}
// error if natoms > 0 yet no atoms were read
if (atom->natoms > 0 && atomflag == 0)
error->all(FLERR,"No atoms in data file");
// close file
if (me == 0) {
if (compressed) pclose(fp);
else fclose(fp);
}
// done if this was 2nd pass
if (!firstpass) break;
// at end of 1st pass, error check for required sections
// customize for new sections
if ((atom->nbonds && !bondflag) || (atom->nangles && !angleflag) ||
(atom->ndihedrals && !dihedralflag) ||
(atom->nimpropers && !improperflag))
error->one(FLERR,"Needed molecular topology not in data file");
if ((nellipsoids && !ellipsoidflag) || (nlines && !lineflag) ||
(ntris && !triflag) || (nbodies && !bodyflag))
error->one(FLERR,"Needed bonus data not in data file");
// break out of loop if no molecular topology in file
// else make 2nd pass
if (!topoflag) break;
firstpass = 0;
// reallocate bond,angle,diehdral,improper arrays via grow()
// will use new bond,angle,dihedral,improper per-atom values from 1st pass
// will also observe extra settings even if bond/etc topology not in file
// leaves other atom arrays unchanged, since already nmax in length
atom->deallocate_topology();
atom->avec->grow(atom->nmax);
}
// create special bond lists for molecular systems
if (atom->molecular == 1) {
Special special(lmp);
special.build();
}
// for atom style template systems, count total bonds,angles,etc
if (atom->molecular == 2) {
Molecule **onemols = atom->avec->onemols;
int *molindex = atom->molindex;
int *molatom = atom->molatom;
int nlocal = atom->nlocal;
int imol,iatom;
bigint nbonds,nangles,ndihedrals,nimpropers;
nbonds = nangles = ndihedrals = nimpropers = 0;
for (int i = 0; i < nlocal; i++) {
imol = molindex[i];
iatom = molatom[i];
nbonds += onemols[imol]->num_bond[iatom];
nangles += onemols[imol]->num_angle[iatom];
ndihedrals += onemols[imol]->num_dihedral[iatom];
nimpropers += onemols[imol]->num_improper[iatom];
}
MPI_Allreduce(&nbonds,&atom->nbonds,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nangles,&atom->nangles,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&ndihedrals,&atom->ndihedrals,1,MPI_LMP_BIGINT,MPI_SUM,world);
MPI_Allreduce(&nimpropers,&atom->nimpropers,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (!force->newton_bond) {
atom->nbonds /= 2;
atom->nangles /= 3;
atom->ndihedrals /= 4;
atom->nimpropers /= 4;
}
if (me == 0) {
if (atom->nbonds) {
if (screen)
fprintf(screen," " BIGINT_FORMAT " template bonds\n",atom->nbonds);
if (logfile)
fprintf(logfile," " BIGINT_FORMAT " template bonds\n",atom->nbonds);
}
if (atom->nangles) {
if (screen)
fprintf(screen," " BIGINT_FORMAT " template angles\n",
atom->nangles);
if (logfile)
fprintf(logfile," " BIGINT_FORMAT " template angles\n",
atom->nangles);
}
if (atom->ndihedrals) {
if (screen)
fprintf(screen," " BIGINT_FORMAT " template dihedrals\n",
atom->nbonds);
if (logfile)
fprintf(logfile," " BIGINT_FORMAT " template bonds\n",
atom->ndihedrals);
}
if (atom->nimpropers) {
if (screen)
fprintf(screen," " BIGINT_FORMAT " template impropers\n",
atom->nimpropers);
if (logfile)
fprintf(logfile," " BIGINT_FORMAT " template impropers\n",
atom->nimpropers);
}
}
}
// for atom style template systems
// insure nbondtypes,etc are still consistent with template molecules,
// in case data file re-defined them
if (atom->molecular == 2) atom->avec->onemols[0]->check_attributes(1);
}
/* ----------------------------------------------------------------------
read free-format header of data file
1st line and blank lines are skipped
non-blank lines are checked for header keywords and leading value is read
header ends with EOF or non-blank line containing no header keyword
if EOF, line is set to blank line
else line has first keyword line for rest of file
------------------------------------------------------------------------- */
void ReadData::header()
{
int n;
char *ptr;
// customize for new sections
const char *section_keywords[NSECTIONS] =
{"Atoms","Velocities","Ellipsoids","Lines","Triangles","Bodies",
"Bonds","Angles","Dihedrals","Impropers",
"Masses","Pair Coeffs","PairIJ Coeffs","Bond Coeffs","Angle Coeffs",
"Dihedral Coeffs","Improper Coeffs",
"BondBond Coeffs","BondAngle Coeffs","MiddleBondTorsion Coeffs",
"EndBondTorsion Coeffs","AngleTorsion Coeffs",
"AngleAngleTorsion Coeffs","BondBond13 Coeffs","AngleAngle Coeffs"};
// skip 1st line of file
if (me == 0) {
char *eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
}
while (1) {
// read a line and bcast length if flag is set
if (me == 0) {
if (fgets(line,MAXLINE,fp) == NULL) n = 0;
else n = strlen(line) + 1;
}
MPI_Bcast(&n,1,MPI_INT,0,world);
// if n = 0 then end-of-file so return with blank line
if (n == 0) {
line[0] = '\0';
return;
}
MPI_Bcast(line,n,MPI_CHAR,0,world);
// trim anything from '#' onward
// if line is blank, continue
if ((ptr = strchr(line,'#'))) *ptr = '\0';
if (strspn(line," \t\n\r") == strlen(line)) continue;
// allow special fixes first chance to match and process the line
// if fix matches, continue to next header line
if (nfix) {
for (n = 0; n < nfix; n++) {
if (!fix_header[n]) continue;
if (strstr(line,fix_header[n])) {
modify->fix[fix_index[n]]->read_data_header(line);
break;
}
}
if (n < nfix) continue;
}
// search line for header keyword and set corresponding variable
// customize for new header lines
if (strstr(line,"atoms")) {
sscanf(line,BIGINT_FORMAT,&atom->natoms);
// check for these first
// otherwise "triangles" will be matched as "angles"
} else if (strstr(line,"ellipsoids")) {
if (!avec_ellipsoid)
error->all(FLERR,"No ellipsoids allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nellipsoids);
} else if (strstr(line,"lines")) {
if (!avec_line)
error->all(FLERR,"No lines allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nlines);
} else if (strstr(line,"triangles")) {
if (!avec_tri)
error->all(FLERR,"No triangles allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&ntris);
} else if (strstr(line,"bodies")) {
if (!avec_body)
error->all(FLERR,"No bodies allowed with this atom style");
sscanf(line,BIGINT_FORMAT,&nbodies);
}
else if (strstr(line,"bonds")) sscanf(line,BIGINT_FORMAT,&atom->nbonds);
else if (strstr(line,"angles")) sscanf(line,BIGINT_FORMAT,&atom->nangles);
else if (strstr(line,"dihedrals")) sscanf(line,BIGINT_FORMAT,
&atom->ndihedrals);
else if (strstr(line,"impropers")) sscanf(line,BIGINT_FORMAT,
&atom->nimpropers);
else if (strstr(line,"atom types")) sscanf(line,"%d",&atom->ntypes);
else if (strstr(line,"bond types")) sscanf(line,"%d",&atom->nbondtypes);
else if (strstr(line,"angle types")) sscanf(line,"%d",&atom->nangletypes);
else if (strstr(line,"dihedral types"))
sscanf(line,"%d",&atom->ndihedraltypes);
else if (strstr(line,"improper types"))
sscanf(line,"%d",&atom->nimpropertypes);
else if (strstr(line,"extra bond per atom"))
sscanf(line,"%d",&atom->extra_bond_per_atom);
else if (strstr(line,"extra angle per atom"))
sscanf(line,"%d",&atom->extra_angle_per_atom);
else if (strstr(line,"extra dihedral per atom"))
sscanf(line,"%d",&atom->extra_dihedral_per_atom);
else if (strstr(line,"extra improper per atom"))
sscanf(line,"%d",&atom->extra_improper_per_atom);
else if (strstr(line,"extra special per atom"))
sscanf(line,"%d",&force->special_extra);
else if (strstr(line,"extra bonds per atom"))
sscanf(line,"%d",&atom->extra_bond_per_atom);
else if (strstr(line,"extra angles per atom"))
sscanf(line,"%d",&atom->extra_angle_per_atom);
else if (strstr(line,"extra dihedrals per atom"))
sscanf(line,"%d",&atom->extra_dihedral_per_atom);
else if (strstr(line,"extra impropers per atom"))
sscanf(line,"%d",&atom->extra_improper_per_atom);
else if (strstr(line,"xlo xhi"))
sscanf(line,"%lg %lg",&domain->boxlo[0],&domain->boxhi[0]);
else if (strstr(line,"ylo yhi"))
sscanf(line,"%lg %lg",&domain->boxlo[1],&domain->boxhi[1]);
else if (strstr(line,"zlo zhi"))
sscanf(line,"%lg %lg",&domain->boxlo[2],&domain->boxhi[2]);
else if (strstr(line,"xy xz yz")) {
domain->triclinic = 1;
sscanf(line,"%lg %lg %lg",&domain->xy,&domain->xz,&domain->yz);
} else break;
}
// error check on total system size
if (atom->natoms < 0 || atom->natoms >= MAXBIGINT ||
atom->nbonds < 0 || atom->nbonds >= MAXBIGINT ||
atom->nangles < 0 || atom->nangles >= MAXBIGINT ||
atom->ndihedrals < 0 || atom->ndihedrals >= MAXBIGINT ||
atom->nimpropers < 0 || atom->nimpropers >= MAXBIGINT)
error->all(FLERR,"System in data file is too big");
// check that exiting string is a valid section keyword
parse_keyword(1);
for (n = 0; n < NSECTIONS; n++)
if (strcmp(keyword,section_keywords[n]) == 0) break;
if (n == NSECTIONS) {
char str[128];
sprintf(str,"Unknown identifier in data file: %s",keyword);
error->all(FLERR,str);
}
// error checks on header values
// must be consistent with atom style and other header values
if ((atom->nbonds || atom->nbondtypes) &&
atom->avec->bonds_allow == 0)
error->all(FLERR,"No bonds allowed with this atom style");
if ((atom->nangles || atom->nangletypes) &&
atom->avec->angles_allow == 0)
error->all(FLERR,"No angles allowed with this atom style");
if ((atom->ndihedrals || atom->ndihedraltypes) &&
atom->avec->dihedrals_allow == 0)
error->all(FLERR,"No dihedrals allowed with this atom style");
if ((atom->nimpropers || atom->nimpropertypes) &&
atom->avec->impropers_allow == 0)
error->all(FLERR,"No impropers allowed with this atom style");
if (atom->nbonds > 0 && atom->nbondtypes <= 0)
error->all(FLERR,"Bonds defined but no bond types");
if (atom->nangles > 0 && atom->nangletypes <= 0)
error->all(FLERR,"Angles defined but no angle types");
if (atom->ndihedrals > 0 && atom->ndihedraltypes <= 0)
error->all(FLERR,"Dihedrals defined but no dihedral types");
if (atom->nimpropers > 0 && atom->nimpropertypes <= 0)
error->all(FLERR,"Impropers defined but no improper types");
if (atom->molecular == 2) {
if (atom->nbonds || atom->nangles || atom->ndihedrals || atom->nimpropers)
error->all(FLERR,"No molecule topology allowed with atom style template");
}
}
/* ----------------------------------------------------------------------
read all atoms
------------------------------------------------------------------------- */
void ReadData::atoms()
{
int nchunk,eof;
if (me == 0) {
if (screen) fprintf(screen," reading atoms ...\n");
if (logfile) fprintf(logfile," reading atoms ...\n");
}
bigint nread = 0;
bigint natoms = atom->natoms;
while (nread < natoms) {
nchunk = MIN(natoms-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
atom->data_atoms(nchunk,buffer);
nread += nchunk;
}
// check that all atoms were assigned correctly
bigint tmp = atom->nlocal;
MPI_Allreduce(&tmp,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " atoms\n",natoms);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " atoms\n",natoms);
}
if (natoms != atom->natoms)
error->all(FLERR,"Did not assign all atoms correctly");
// check that atom IDs are valid
atom->tag_check();
// create global mapping of atoms
if (atom->map_style) {
atom->map_init();
atom->map_set();
}
}
/* ----------------------------------------------------------------------
read all velocities
to find atoms, must build atom map if not a molecular system
------------------------------------------------------------------------- */
void ReadData::velocities()
{
int nchunk,eof;
if (me == 0) {
if (screen) fprintf(screen," reading velocities ...\n");
if (logfile) fprintf(logfile," reading velocities ...\n");
}
int mapflag = 0;
if (atom->map_style == 0) {
mapflag = 1;
atom->map_init();
atom->map_set();
}
bigint nread = 0;
bigint natoms = atom->natoms;
while (nread < natoms) {
nchunk = MIN(natoms-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
atom->data_vels(nchunk,buffer);
nread += nchunk;
}
if (mapflag) {
atom->map_delete();
atom->map_style = 0;
}
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " velocities\n",natoms);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " velocities\n",natoms);
}
}
/* ----------------------------------------------------------------------
scan or read all bonds
------------------------------------------------------------------------- */
void ReadData::bonds(int firstpass)
{
if (me == 0) {
if (firstpass) {
if (screen) fprintf(screen," scanning bonds ...\n");
if (logfile) fprintf(logfile," scanning bonds ...\n");
} else {
if (screen) fprintf(screen," reading bonds ...\n");
if (logfile) fprintf(logfile," reading bonds ...\n");
}
}
// allocate count if firstpass
int nlocal = atom->nlocal;
int *count = NULL;
if (firstpass) {
memory->create(count,nlocal,"read_data:count");
for (int i = 0; i < nlocal; i++) count[i] = 0;
}
// read and process bonds
int nchunk,eof;
bigint nread = 0;
bigint nbonds = atom->nbonds;
while (nread < nbonds) {
nchunk = MIN(nbonds-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
atom->data_bonds(nchunk,buffer,count);
nread += nchunk;
}
// if firstpass: tally max bond/atom and return
if (firstpass) {
int max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
int maxall;
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
maxall += atom->extra_bond_per_atom;
if (me == 0) {
if (screen) fprintf(screen," %d = max bonds/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max bonds/atom\n",maxall);
}
atom->bond_per_atom = maxall;
memory->destroy(count);
return;
}
// if 2nd pass: check that bonds were assigned correctly
bigint n = 0;
for (int i = 0; i < nlocal; i++) n += atom->num_bond[i];
bigint sum;
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
if (!force->newton_bond) factor = 2;
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " bonds\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " bonds\n",sum/factor);
}
if (sum != factor*atom->nbonds)
error->all(FLERR,"Bonds assigned incorrectly");
}
/* ----------------------------------------------------------------------
scan or read all angles
------------------------------------------------------------------------- */
void ReadData::angles(int firstpass)
{
if (me == 0) {
if (firstpass) {
if (screen) fprintf(screen," scanning angles ...\n");
if (logfile) fprintf(logfile," scanning angles ...\n");
} else {
if (screen) fprintf(screen," reading angles ...\n");
if (logfile) fprintf(logfile," reading angles ...\n");
}
}
// allocate count if firstpass
int nlocal = atom->nlocal;
int *count = NULL;
if (firstpass) {
memory->create(count,nlocal,"read_data:count");
for (int i = 0; i < nlocal; i++) count[i] = 0;
}
// read and process angles
int nchunk,eof;
bigint nread = 0;
bigint nangles = atom->nangles;
while (nread < nangles) {
nchunk = MIN(nangles-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
atom->data_angles(nchunk,buffer,count);
nread += nchunk;
}
// if firstpass: tally max angle/atom and return
if (firstpass) {
int max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
int maxall;
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
maxall += atom->extra_angle_per_atom;
if (me == 0) {
if (screen) fprintf(screen," %d = max angles/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max angles/atom\n",maxall);
}
atom->angle_per_atom = maxall;
memory->destroy(count);
return;
}
// if 2nd pass: check that angles were assigned correctly
bigint n = 0;
for (int i = 0; i < nlocal; i++) n += atom->num_angle[i];
bigint sum;
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
if (!force->newton_bond) factor = 3;
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " angles\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " angles\n",sum/factor);
}
if (sum != factor*atom->nangles)
error->all(FLERR,"Angles assigned incorrectly");
}
/* ----------------------------------------------------------------------
scan or read all dihedrals
------------------------------------------------------------------------- */
void ReadData::dihedrals(int firstpass)
{
if (me == 0) {
if (firstpass) {
if (screen) fprintf(screen," scanning dihedrals ...\n");
if (logfile) fprintf(logfile," scanning dihedrals ...\n");
} else {
if (screen) fprintf(screen," reading dihedrals ...\n");
if (logfile) fprintf(logfile," reading dihedrals ...\n");
}
}
// allocate count if firstpass
int nlocal = atom->nlocal;
int *count = NULL;
if (firstpass) {
memory->create(count,nlocal,"read_data:count");
for (int i = 0; i < nlocal; i++) count[i] = 0;
}
// read and process dihedrals
int nchunk,eof;
bigint nread = 0;
bigint ndihedrals = atom->ndihedrals;
while (nread < ndihedrals) {
nchunk = MIN(ndihedrals-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
atom->data_dihedrals(nchunk,buffer,count);
nread += nchunk;
}
// if firstpass: tally max dihedral/atom and return
if (firstpass) {
int max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
int maxall;
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
maxall += atom->extra_dihedral_per_atom;
if (me == 0) {
if (screen) fprintf(screen," %d = max dihedrals/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max dihedrals/atom\n",maxall);
}
atom->dihedral_per_atom = maxall;
memory->destroy(count);
return;
}
// if 2nd pass: check that dihedrals were assigned correctly
bigint n = 0;
for (int i = 0; i < nlocal; i++) n += atom->num_dihedral[i];
bigint sum;
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
if (!force->newton_bond) factor = 4;
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " dihedrals\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " dihedrals\n",sum/factor);
}
if (sum != factor*atom->ndihedrals)
error->all(FLERR,"Dihedrals assigned incorrectly");
}
/* ----------------------------------------------------------------------
scan or read all impropers
------------------------------------------------------------------------- */
void ReadData::impropers(int firstpass)
{
if (me == 0) {
if (firstpass) {
if (screen) fprintf(screen," scanning impropers ...\n");
if (logfile) fprintf(logfile," scanning impropers ...\n");
} else {
if (screen) fprintf(screen," reading impropers ...\n");
if (logfile) fprintf(logfile," reading impropers ...\n");
}
}
// allocate count if firstpass
int nlocal = atom->nlocal;
int *count = NULL;
if (firstpass) {
memory->create(count,nlocal,"read_data:count");
for (int i = 0; i < nlocal; i++) count[i] = 0;
}
// read and process impropers
int nchunk,eof;
bigint nread = 0;
bigint nimpropers = atom->nimpropers;
while (nread < nimpropers) {
nchunk = MIN(nimpropers-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
atom->data_impropers(nchunk,buffer,count);
nread += nchunk;
}
// if firstpass: tally max improper/atom and return
if (firstpass) {
int max = 0;
for (int i = 0; i < nlocal; i++) max = MAX(max,count[i]);
int maxall;
MPI_Allreduce(&max,&maxall,1,MPI_INT,MPI_MAX,world);
maxall += atom->extra_improper_per_atom;
if (me == 0) {
if (screen) fprintf(screen," %d = max impropers/atom\n",maxall);
if (logfile) fprintf(logfile," %d = max impropers/atom\n",maxall);
}
atom->improper_per_atom = maxall;
memory->destroy(count);
return;
}
// if 2nd pass: check that impropers were assigned correctly
bigint n = 0;
for (int i = 0; i < nlocal; i++) n += atom->num_improper[i];
bigint sum;
MPI_Allreduce(&n,&sum,1,MPI_LMP_BIGINT,MPI_SUM,world);
int factor = 1;
if (!force->newton_bond) factor = 4;
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " impropers\n",sum/factor);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " impropers\n",sum/factor);
}
if (sum != factor*atom->nimpropers)
error->all(FLERR,"Impropers assigned incorrectly");
}
/* ----------------------------------------------------------------------
read all bonus data
to find atoms, must build atom map if not a molecular system
------------------------------------------------------------------------- */
void ReadData::bonus(bigint nbonus, AtomVec *ptr, const char *type)
{
int nchunk,eof;
int mapflag = 0;
if (atom->map_style == 0) {
mapflag = 1;
atom->map_init();
atom->map_set();
}
bigint nread = 0;
bigint natoms = nbonus;
while (nread < natoms) {
nchunk = MIN(natoms-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
atom->data_bonus(nchunk,buffer,ptr);
nread += nchunk;
}
if (mapflag) {
atom->map_delete();
atom->map_style = 0;
}
if (me == 0) {
if (screen) fprintf(screen," " BIGINT_FORMAT " %s\n",natoms,type);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " %s\n",natoms,type);
}
}
/* ----------------------------------------------------------------------
read all body data
variable amount of info per body, described by ninteger and ndouble
to find atoms, must build atom map if not a molecular system
if not firstpass, just read but no processing of data
------------------------------------------------------------------------- */
void ReadData::bodies(int firstpass)
{
int i,m,nchunk,nline,nmax,ninteger,ndouble,tmp,onebody;
char *eof;
int mapflag = 0;
if (atom->map_style == 0 && firstpass) {
mapflag = 1;
atom->map_init();
atom->map_set();
}
// nmax = max # of bodies to read in this chunk
// nchunk = actual # read
bigint nread = 0;
bigint natoms = nbodies;
while (nread < natoms) {
if (natoms-nread > CHUNK) nmax = CHUNK;
else nmax = natoms-nread;
if (me == 0) {
nchunk = 0;
nline = 0;
m = 0;
while (nchunk < nmax && nline <= CHUNK-MAXBODY) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
sscanf(&buffer[m],"%d %d %d",&tmp,&ninteger,&ndouble);
m += strlen(&buffer[m]);
onebody = 0;
if (ninteger) onebody += (ninteger-1)/10 + 1;
if (ndouble) onebody += (ndouble-1)/10 + 1;
if (onebody+1 > MAXBODY)
error->one(FLERR,
"Too many lines in one body in data file - boost MAXBODY");
for (i = 0; i < onebody; i++) {
eof = fgets(&buffer[m],MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
m += strlen(&buffer[m]);
}
nchunk++;
nline += onebody+1;
}
if (buffer[m-1] != '\n') strcpy(&buffer[m++],"\n");
m++;
}
MPI_Bcast(&nchunk,1,MPI_INT,0,world);
MPI_Bcast(&m,1,MPI_INT,0,world);
MPI_Bcast(buffer,m,MPI_CHAR,0,world);
if (firstpass) atom->data_bodies(nchunk,buffer,avec_body);
nread += nchunk;
}
if (mapflag && firstpass) {
atom->map_delete();
atom->map_style = 0;
}
if (me == 0 && firstpass) {
if (screen) fprintf(screen," " BIGINT_FORMAT " bodies\n",natoms);
if (logfile) fprintf(logfile," " BIGINT_FORMAT " bodies\n",natoms);
}
}
/* ---------------------------------------------------------------------- */
void ReadData::mass()
{
char *next;
char *buf = new char[atom->ntypes*MAXLINE];
int eof = comm->read_lines_from_file(fp,atom->ntypes,MAXLINE,buf);
if (eof) error->all(FLERR,"Unexpected end of data file");
char *original = buf;
for (int i = 0; i < atom->ntypes; i++) {
next = strchr(buf,'\n');
*next = '\0';
atom->set_mass(buf);
buf = next + 1;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::paircoeffs()
{
char *next;
char *buf = new char[atom->ntypes*MAXLINE];
int eof = comm->read_lines_from_file(fp,atom->ntypes,MAXLINE,buf);
if (eof) error->all(FLERR,"Unexpected end of data file");
char *original = buf;
for (int i = 0; i < atom->ntypes; i++) {
next = strchr(buf,'\n');
*next = '\0';
parse_coeffs(buf,NULL,1);
if (narg == 0) error->all(FLERR,"Unexpected end of PairCoeffs section");
force->pair->coeff(narg,arg);
buf = next + 1;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::pairIJcoeffs()
{
int i,j;
char *next;
int nsq = atom->ntypes* (atom->ntypes+1) / 2;
char *buf = new char[nsq * MAXLINE];
int eof = comm->read_lines_from_file(fp,nsq,MAXLINE,buf);
if (eof) error->all(FLERR,"Unexpected end of data file");
char *original = buf;
for (i = 0; i < atom->ntypes; i++)
for (j = i; j < atom->ntypes; j++) {
next = strchr(buf,'\n');
*next = '\0';
parse_coeffs(buf,NULL,0);
if (narg == 0) error->all(FLERR,"Unexpected end of PairCoeffs section");
force->pair->coeff(narg,arg);
buf = next + 1;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::bondcoeffs()
{
char *next;
char *buf = new char[atom->nbondtypes*MAXLINE];
int eof = comm->read_lines_from_file(fp,atom->nbondtypes,MAXLINE,buf);
if (eof) error->all(FLERR,"Unexpected end of data file");
char *original = buf;
for (int i = 0; i < atom->nbondtypes; i++) {
next = strchr(buf,'\n');
*next = '\0';
parse_coeffs(buf,NULL,0);
if (narg == 0) error->all(FLERR,"Unexpected end of BondCoeffs section");
force->bond->coeff(narg,arg);
buf = next + 1;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::anglecoeffs(int which)
{
char *next;
char *buf = new char[atom->nangletypes*MAXLINE];
int eof = comm->read_lines_from_file(fp,atom->nangletypes,MAXLINE,buf);
if (eof) error->all(FLERR,"Unexpected end of data file");
char *original = buf;
for (int i = 0; i < atom->nangletypes; i++) {
next = strchr(buf,'\n');
*next = '\0';
if (which == 0) parse_coeffs(buf,NULL,0);
else if (which == 1) parse_coeffs(buf,"bb",0);
else if (which == 2) parse_coeffs(buf,"ba",0);
if (narg == 0) error->all(FLERR,"Unexpected end of AngleCoeffs section");
force->angle->coeff(narg,arg);
buf = next + 1;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::dihedralcoeffs(int which)
{
char *next;
char *buf = new char[atom->ndihedraltypes*MAXLINE];
int eof = comm->read_lines_from_file(fp,atom->ndihedraltypes,MAXLINE,buf);
if (eof) error->all(FLERR,"Unexpected end of data file");
char *original = buf;
for (int i = 0; i < atom->ndihedraltypes; i++) {
next = strchr(buf,'\n');
*next = '\0';
if (which == 0) parse_coeffs(buf,NULL,0);
else if (which == 1) parse_coeffs(buf,"mbt",0);
else if (which == 2) parse_coeffs(buf,"ebt",0);
else if (which == 3) parse_coeffs(buf,"at",0);
else if (which == 4) parse_coeffs(buf,"aat",0);
else if (which == 5) parse_coeffs(buf,"bb13",0);
if (narg == 0) error->all(FLERR,"Unexpected end of DihedralCoeffs section");
force->dihedral->coeff(narg,arg);
buf = next + 1;
}
delete [] original;
}
/* ---------------------------------------------------------------------- */
void ReadData::impropercoeffs(int which)
{
char *next;
char *buf = new char[atom->nimpropertypes*MAXLINE];
int eof = comm->read_lines_from_file(fp,atom->nimpropertypes,MAXLINE,buf);
if (eof) error->all(FLERR,"Unexpected end of data file");
char *original = buf;
for (int i = 0; i < atom->nimpropertypes; i++) {
next = strchr(buf,'\n');
*next = '\0';
if (which == 0) parse_coeffs(buf,NULL,0);
else if (which == 1) parse_coeffs(buf,"aa",0);
if (narg == 0) error->all(FLERR,"Unexpected end of ImproperCoeffs section");
force->improper->coeff(narg,arg);
buf = next + 1;
}
delete [] original;
}
/* ----------------------------------------------------------------------
read fix section, pass lines to fix to process
n = index of fix
------------------------------------------------------------------------- */
void ReadData::fix(int ifix, char *keyword)
{
int nchunk,eof;
bigint nline = modify->fix[ifix]->read_data_skip_lines(keyword);
bigint nread = 0;
while (nread < nline) {
nchunk = MIN(nline-nread,CHUNK);
eof = comm->read_lines_from_file(fp,nchunk,MAXLINE,buffer);
if (eof) error->all(FLERR,"Unexpected end of data file");
modify->fix[ifix]->read_data_section(keyword,nchunk,buffer);
nread += nchunk;
}
}
/* ----------------------------------------------------------------------
reallocate the count vector from cmax to amax+1 and return new length
zero new locations
------------------------------------------------------------------------- */
int ReadData::reallocate(int **pcount, int cmax, int amax)
{
int *count = *pcount;
memory->grow(count,amax+1,"read_data:count");
for (int i = cmax; i <= amax; i++) count[i] = 0;
*pcount = count;
return amax+1;
}
/* ----------------------------------------------------------------------
proc 0 opens data file
test if gzipped
------------------------------------------------------------------------- */
void ReadData::open(char *file)
{
compressed = 0;
char *suffix = file + strlen(file) - 3;
if (suffix > file && strcmp(suffix,".gz") == 0) compressed = 1;
if (!compressed) fp = fopen(file,"r");
else {
#ifdef LAMMPS_GZIP
char gunzip[128];
sprintf(gunzip,"gzip -c -d %s",file);
#ifdef _WIN32
fp = _popen(gunzip,"rb");
#else
fp = popen(gunzip,"r");
#endif
#else
error->one(FLERR,"Cannot open gzipped file");
#endif
}
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open file %s",file);
error->one(FLERR,str);
}
}
/* ----------------------------------------------------------------------
grab next keyword
read lines until one is non-blank
keyword is all text on line w/out leading & trailing white space
optional style can be appended after comment char '#'
read one additional line (assumed blank)
if any read hits EOF, set keyword to empty
if first = 1, line variable holds non-blank line that ended header
------------------------------------------------------------------------- */
void ReadData::parse_keyword(int first)
{
int eof = 0;
// proc 0 reads upto non-blank line plus 1 following line
// eof is set to 1 if any read hits end-of-file
if (me == 0) {
if (!first) {
if (fgets(line,MAXLINE,fp) == NULL) eof = 1;
}
while (eof == 0 && strspn(line," \t\n\r") == strlen(line)) {
if (fgets(line,MAXLINE,fp) == NULL) eof = 1;
}
if (fgets(buffer,MAXLINE,fp) == NULL) eof = 1;
}
// if eof, set keyword empty and return
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) {
keyword[0] = '\0';
return;
}
// bcast keyword line to all procs
int n;
if (me == 0) n = strlen(line) + 1;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// store optional "style" following comment char '#' after keyword
char *ptr;
if ((ptr = strchr(line,'#'))) {
*ptr++ = '\0';
while (*ptr == ' ' || *ptr == '\t') ptr++;
int stop = strlen(ptr) - 1;
while (ptr[stop] == ' ' || ptr[stop] == '\t'
|| ptr[stop] == '\n' || ptr[stop] == '\r') stop--;
ptr[stop+1] = '\0';
strcpy(style,ptr);
} else style[0] = '\0';
// copy non-whitespace portion of line into keyword
int start = strspn(line," \t\n\r");
int stop = strlen(line) - 1;
while (line[stop] == ' ' || line[stop] == '\t'
|| line[stop] == '\n' || line[stop] == '\r') stop--;
line[stop+1] = '\0';
strcpy(keyword,&line[start]);
}
/* ----------------------------------------------------------------------
proc 0 reads N lines from file
could be skipping Natoms lines, so use bigints
------------------------------------------------------------------------- */
void ReadData::skip_lines(bigint n)
{
if (me) return;
char *eof;
for (bigint i = 0; i < n; i++) eof = fgets(line,MAXLINE,fp);
if (eof == NULL) error->one(FLERR,"Unexpected end of data file");
}
/* ----------------------------------------------------------------------
parse a line of coeffs into words, storing them in narg,arg
trim anything from '#' onward
word strings remain in line, are not copied
if addstr != NULL, add addstr as extra arg for class2 angle/dihedral/improper
if 2nd word starts with letter, then is hybrid style, add addstr after it
else add addstr before 2nd word
if dupflag, duplicate 1st word, so pair_coeff "2" becomes "2 2"
------------------------------------------------------------------------- */
void ReadData::parse_coeffs(char *line, const char *addstr, int dupflag)
{
char *ptr;
if ((ptr = strchr(line,'#'))) *ptr = '\0';
narg = 0;
char *word = strtok(line," \t\n\r\f");
while (word) {
if (narg == maxarg) {
maxarg += DELTA;
arg = (char **)
memory->srealloc(arg,maxarg*sizeof(char *),"read_data:arg");
}
if (addstr && narg == 1 && !islower(word[0])) arg[narg++] = (char *) addstr;
arg[narg++] = word;
if (addstr && narg == 2 && islower(word[0])) arg[narg++] = (char *) addstr;
if (dupflag && narg == 1) arg[narg++] = word;
word = strtok(NULL," \t\n\r\f");
}
}
/* ----------------------------------------------------------------------
compare two style strings if they both exist
one = comment in data file section, two = currently-defined style
ignore suffixes listed in suffixes array at top of file
------------------------------------------------------------------------- */
int ReadData::style_match(const char *one, const char *two)
{
int i,delta,len,len1,len2;
if ((one == NULL) || (two == NULL)) return 1;
len1 = strlen(one);
len2 = strlen(two);
for (i = 0; suffixes[i] != NULL; i++) {
len = strlen(suffixes[i]);
if ((delta = len1 - len) > 0)
if (strcmp(one+delta,suffixes[i]) == 0) len1 = delta;
if ((delta = len2 - len) > 0)
if (strcmp(two+delta,suffixes[i]) == 0) len2 = delta;
}
if ((len1 == 0) || (len1 == len2) || (strncmp(one,two,len1) == 0)) return 1;
return 0;
}
diff --git a/src/read_dump.cpp b/src/read_dump.cpp
index 14b208e99..45ff1f81b 100644
--- a/src/read_dump.cpp
+++ b/src/read_dump.cpp
@@ -1,960 +1,965 @@
/* ----------------------------------------------------------------------
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: Timothy Sirk (ARL)
------------------------------------------------------------------------- */
+// lmptype.h must be first b/c this file uses MAXBIGINT and includes mpi.h
+// due to OpenMPI bug which sets INT64_MAX via its mpi.h
+// before lmptype.h can set flags to insure it is done correctly
+
+#include "lmptype.h"
#include "mpi.h"
#include "string.h"
#include "stdlib.h"
#include "read_dump.h"
#include "reader.h"
#include "style_reader.h"
#include "atom.h"
#include "atom_vec.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "domain.h"
#include "comm.h"
#include "irregular.h"
#include "error.h"
#include "memory.h"
using namespace LAMMPS_NS;
#define CHUNK 1024
#define EPSILON 1.0e-6
// also in reader_native.cpp
enum{ID,TYPE,X,Y,Z,VX,VY,VZ,Q,IX,IY,IZ};
enum{UNSET,NOSCALE_NOWRAP,NOSCALE_WRAP,SCALE_NOWRAP,SCALE_WRAP};
/* ---------------------------------------------------------------------- */
ReadDump::ReadDump(LAMMPS *lmp) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
dimension = domain->dimension;
triclinic = domain->triclinic;
nfile = 0;
files = NULL;
nfield = 0;
fieldtype = NULL;
fieldlabel = NULL;
fields = NULL;
int n = strlen("native") + 1;
readerstyle = new char[n];
strcpy(readerstyle,"native");
reader = NULL;
fp = NULL;
}
/* ---------------------------------------------------------------------- */
ReadDump::~ReadDump()
{
for (int i = 0; i < nfile; i++) delete [] files[i];
delete [] files;
for (int i = 0; i < nfield; i++) delete [] fieldlabel[i];
delete [] fieldlabel;
delete [] fieldtype;
delete [] readerstyle;
memory->destroy(fields);
delete reader;
}
/* ---------------------------------------------------------------------- */
void ReadDump::command(int narg, char **arg)
{
if (domain->box_exist == 0)
error->all(FLERR,"Read_dump command before simulation box is defined");
if (narg < 2) error->all(FLERR,"Illegal read_dump command");
store_files(1,&arg[0]);
bigint nstep = ATOBIGINT(arg[1]);
int nremain = narg - 2;
if (nremain) nremain = fields_and_keywords(nremain,&arg[narg-nremain]);
else nremain = fields_and_keywords(0,NULL);
if (nremain) setup_reader(nremain,&arg[narg-nremain]);
else setup_reader(0,NULL);
// find the snapshot and read/bcast/process header info
if (me == 0 && screen) fprintf(screen,"Scanning dump file ...\n");
bigint ntimestep = seek(nstep,1);
if (ntimestep < 0)
error->all(FLERR,"Dump file does not contain requested snapshot");
header(1);
// reset timestep to nstep
update->reset_timestep(nstep);
// counters
// read in the snapshot and reset system
if (me == 0 && screen)
fprintf(screen,"Reading snapshot from dump file ...\n");
bigint natoms_prev = atom->natoms;
atoms();
if (me == 0) reader->close_file();
// print out stats
bigint npurge_all,nreplace_all,ntrim_all,nadd_all;
bigint tmp;
tmp = npurge;
MPI_Allreduce(&tmp,&npurge_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
tmp = nreplace;
MPI_Allreduce(&tmp,&nreplace_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
tmp = ntrim;
MPI_Allreduce(&tmp,&ntrim_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
tmp = nadd;
MPI_Allreduce(&tmp,&nadd_all,1,MPI_LMP_BIGINT,MPI_SUM,world);
domain->print_box(" ");
if (me == 0) {
if (screen) {
fprintf(screen," " BIGINT_FORMAT " atoms before read\n",natoms_prev);
fprintf(screen," " BIGINT_FORMAT " atoms in snapshot\n",nsnapatoms);
fprintf(screen," " BIGINT_FORMAT " atoms purged\n",npurge_all);
fprintf(screen," " BIGINT_FORMAT " atoms replaced\n",nreplace_all);
fprintf(screen," " BIGINT_FORMAT " atoms trimmed\n",ntrim_all);
fprintf(screen," " BIGINT_FORMAT " atoms added\n",nadd_all);
fprintf(screen," " BIGINT_FORMAT " atoms after read\n",atom->natoms);
}
if (logfile) {
fprintf(logfile," " BIGINT_FORMAT " atoms before read\n",natoms_prev);
fprintf(logfile," " BIGINT_FORMAT " atoms in snapshot\n",nsnapatoms);
fprintf(logfile," " BIGINT_FORMAT " atoms purged\n",npurge_all);
fprintf(logfile," " BIGINT_FORMAT " atoms replaced\n",nreplace_all);
fprintf(logfile," " BIGINT_FORMAT " atoms trimmed\n",ntrim_all);
fprintf(logfile," " BIGINT_FORMAT " atoms added\n",nadd_all);
fprintf(logfile," " BIGINT_FORMAT " atoms after read\n",atom->natoms);
}
}
}
/* ---------------------------------------------------------------------- */
void ReadDump::store_files(int nstr, char **str)
{
nfile = nstr;
files = new char*[nfile];
for (int i = 0; i < nfile; i++) {
int n = strlen(str[i]) + 1;
files[i] = new char[n];
strcpy(files[i],str[i]);
}
}
/* ---------------------------------------------------------------------- */
void ReadDump::setup_reader(int narg, char **arg)
{
// allocate snapshot field buffer
memory->create(fields,CHUNK,nfield,"read_dump:fields");
// create reader class
// match readerstyle to options in style_reader.h
if (0) return; // dummy line to enable else-if macro expansion
#define READER_CLASS
#define ReaderStyle(key,Class) \
else if (strcmp(readerstyle,#key) == 0) reader = new Class(lmp);
#include "style_reader.h"
#undef READER_CLASS
// unrecognized style
else error->all(FLERR,"Unknown dump reader style");
// pass any arguments to reader
if (narg > 0) reader->settings(narg,arg);
}
/* ----------------------------------------------------------------------
seek Nrequest timestep in one or more dump files
if exact = 1, must find exactly Nrequest
if exact = 0, find first step >= Nrequest
return matching ntimestep or -1 if did not find a match
------------------------------------------------------------------------- */
bigint ReadDump::seek(bigint nrequest, int exact)
{
int ifile,eofflag;
bigint ntimestep;
if (me == 0) {
// exit file loop when dump timestep >= nrequest
// or files exhausted
for (ifile = 0; ifile < nfile; ifile++) {
ntimestep = -1;
reader->open_file(files[ifile]);
while (1) {
eofflag = reader->read_time(ntimestep);
if (eofflag) break;
if (ntimestep >= nrequest) break;
reader->skip();
}
if (ntimestep >= nrequest) break;
reader->close_file();
}
currentfile = ifile;
if (ntimestep < nrequest) ntimestep = -1;
if (exact && ntimestep != nrequest) ntimestep = -1;
if (ntimestep < 0) reader->close_file();
}
MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world);
return ntimestep;
}
/* ----------------------------------------------------------------------
find next matching snapshot in one or more dump files
Ncurrent = current timestep from last snapshot
Nlast = match no timestep bigger than Nlast
Nevery = only match timesteps that are a multiple of Nevery
Nskip = skip every this many timesteps
return matching ntimestep or -1 if did not find a match
------------------------------------------------------------------------- */
bigint ReadDump::next(bigint ncurrent, bigint nlast, int nevery, int nskip)
{
int ifile,eofflag;
bigint ntimestep;
if (me == 0) {
// exit file loop when dump timestep matches all criteria
// or files exhausted
int iskip = 0;
for (ifile = currentfile; ifile < nfile; ifile++) {
ntimestep = -1;
if (ifile != currentfile) reader->open_file(files[ifile]);
while (1) {
eofflag = reader->read_time(ntimestep);
if (iskip == nskip) iskip = 0;
iskip++;
if (eofflag) break;
if (ntimestep <= ncurrent) break;
if (ntimestep > nlast) break;
if (nevery && ntimestep % nevery) reader->skip();
else if (iskip < nskip) reader->skip();
else break;
}
if (eofflag) reader->close_file();
else break;
}
currentfile = ifile;
if (eofflag) ntimestep = -1;
if (ntimestep <= ncurrent) ntimestep = -1;
if (ntimestep > nlast) ntimestep = -1;
if (ntimestep < 0) reader->close_file();
}
MPI_Bcast(&ntimestep,1,MPI_LMP_BIGINT,0,world);
return ntimestep;
}
/* ----------------------------------------------------------------------
read and broadcast and store snapshot header info
set nsnapatoms = # of atoms in snapshot
------------------------------------------------------------------------- */
void ReadDump::header(int fieldinfo)
{
int triclinic_snap;
int fieldflag,xflag,yflag,zflag;
if (me == 0)
nsnapatoms = reader->read_header(box,triclinic_snap,
fieldinfo,nfield,fieldtype,fieldlabel,
scaleflag,wrapflag,fieldflag,
xflag,yflag,zflag);
MPI_Bcast(&nsnapatoms,1,MPI_LMP_BIGINT,0,world);
MPI_Bcast(&triclinic_snap,1,MPI_INT,0,world);
MPI_Bcast(&box[0][0],9,MPI_DOUBLE,0,world);
// local copy of snapshot box parameters
// used in xfield,yfield,zfield when converting dump atom to absolute coords
xlo = box[0][0];
xhi = box[0][1];
ylo = box[1][0];
yhi = box[1][1];
zlo = box[2][0];
zhi = box[2][1];
if (triclinic_snap) {
xy = box[0][2];
xz = box[1][2];
yz = box[2][2];
double xdelta = MIN(0.0,xy);
xdelta = MIN(xdelta,xz);
xdelta = MIN(xdelta,xy+xz);
xlo = xlo - xdelta;
xdelta = MAX(0.0,xy);
xdelta = MAX(xdelta,xz);
xdelta = MAX(xdelta,xy+xz);
xhi = xhi - xdelta;
ylo = ylo - MIN(0.0,yz);
yhi = yhi - MAX(0.0,yz);
}
xprd = xhi - xlo;
yprd = yhi - ylo;
zprd = zhi - zlo;
// done if not checking fields
if (!fieldinfo) return;
MPI_Bcast(&fieldflag,1,MPI_INT,0,world);
MPI_Bcast(&xflag,1,MPI_INT,0,world);
MPI_Bcast(&yflag,1,MPI_INT,0,world);
MPI_Bcast(&zflag,1,MPI_INT,0,world);
// error check on current vs new box and fields
// triclinic_snap < 0 means no box info in file
if (triclinic_snap < 0 && boxflag > 0)
error->all(FLERR,"No box information in dump. You have to use 'box no'");
if (triclinic_snap >= 0) {
if ((triclinic_snap && !triclinic) ||
(!triclinic_snap && triclinic))
error->one(FLERR,"Read_dump triclinic status does not match simulation");
}
// error check on requested fields exisiting in dump file
if (fieldflag < 0)
error->one(FLERR,"Read_dump field not found in dump file");
// all explicitly requested x,y,z must have consistent scaling & wrapping
int value = MAX(xflag,yflag);
value = MAX(zflag,value);
if ((xflag != UNSET && xflag != value) ||
(yflag != UNSET && yflag != value) ||
(zflag != UNSET && zflag != value))
error->one(FLERR,
"Read_dump xyz fields do not have consistent scaling/wrapping");
// set scaled/wrapped based on xyz flags
value = UNSET;
if (xflag != UNSET) value = xflag;
if (yflag != UNSET) value = yflag;
if (zflag != UNSET) value = zflag;
if (value == UNSET) {
scaled = wrapped = 0;
} else if (value == NOSCALE_NOWRAP) {
scaled = wrapped = 0;
} else if (value == NOSCALE_WRAP) {
scaled = 0;
wrapped = 1;
} else if (value == SCALE_NOWRAP) {
scaled = 1;
wrapped = 0;
} else if (value == SCALE_WRAP) {
scaled = wrapped = 1;
}
// scaled, triclinic coords require all 3 x,y,z fields, to perform unscaling
// set yindex,zindex = column index of Y and Z fields in fields array
// needed for unscaling to absolute coords in xfield(), yfield(), zfield()
if (scaled && triclinic == 1) {
int flag = 0;
if (xflag == UNSET) flag = 1;
if (yflag == UNSET) flag = 1;
if (dimension == 3 && zflag == UNSET) flag = 1;
if (flag)
error->one(FLERR,"All read_dump x,y,z fields must be specified for "
"scaled, triclinic coords");
for (int i = 0; i < nfield; i++) {
if (fieldtype[i] == Y) yindex = i;
if (fieldtype[i] == Z) zindex = i;
}
}
}
/* ---------------------------------------------------------------------- */
void ReadDump::atoms()
{
// initialize counters
npurge = nreplace = ntrim = nadd = 0;
// if purgeflag set, delete all current atoms
if (purgeflag) {
if (atom->map_style) atom->map_clear();
npurge = atom->nlocal;
atom->nlocal = atom->nghost = 0;
atom->natoms = 0;
}
// to match existing atoms to dump atoms:
// must build map if not a molecular system
int mapflag = 0;
if (atom->map_style == 0) {
mapflag = 1;
atom->map_init();
atom->map_set();
}
// uflag[i] = 1 for each owned atom appearing in dump
// ucflag = similar flag for each chunk atom, used in process_atoms()
int nlocal = atom->nlocal;
memory->create(uflag,nlocal,"read_dump:uflag");
for (int i = 0; i < nlocal; i++) uflag[i] = 0;
memory->create(ucflag,CHUNK,"read_dump:ucflag");
memory->create(ucflag_all,CHUNK,"read_dump:ucflag");
// read, broadcast, and process atoms from snapshot in chunks
addproc = -1;
int nchunk;
bigint nread = 0;
while (nread < nsnapatoms) {
nchunk = MIN(nsnapatoms-nread,CHUNK);
if (me == 0) reader->read_atoms(nchunk,nfield,fields);
MPI_Bcast(&fields[0][0],nchunk*nfield,MPI_DOUBLE,0,world);
process_atoms(nchunk);
nread += nchunk;
}
// if addflag set, add tags to new atoms if possible
if (addflag) {
bigint nblocal = atom->nlocal;
MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (atom->natoms < 0 || atom->natoms >= MAXBIGINT)
error->all(FLERR,"Too many total atoms");
if (atom->tag_enable) atom->tag_extend();
}
// if trimflag set, delete atoms not replaced by snapshot atoms
if (trimflag) {
delete_atoms();
bigint nblocal = atom->nlocal;
MPI_Allreduce(&nblocal,&atom->natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
}
// can now delete uflag arrays
memory->destroy(uflag);
memory->destroy(ucflag);
memory->destroy(ucflag_all);
// delete atom map if created it above
// else reinitialize map for current atoms
// do this before migrating atoms to new procs via Irregular
if (mapflag) {
atom->map_delete();
atom->map_style = 0;
} else {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
// overwrite simulation box with dump snapshot box if requested
// reallocate processors to box
if (boxflag) {
domain->boxlo[0] = xlo;
domain->boxhi[0] = xhi;
domain->boxlo[1] = ylo;
domain->boxhi[1] = yhi;
if (dimension == 3) {
domain->boxlo[2] = zlo;
domain->boxhi[2] = zhi;
}
if (triclinic) {
domain->xy = xy;
if (dimension == 3) {
domain->xz = xz;
domain->yz = yz;
}
}
domain->set_initial_box();
domain->set_global_box();
comm->set_proc_grid(0);
domain->set_local_box();
}
// move atoms back inside simulation box and to new processors
// use remap() instead of pbc() in case atoms moved a long distance
// adjust image flags of all atoms (old and new) based on current box
// use irregular() in case atoms moved a long distance
double **x = atom->x;
imageint *image = atom->image;
nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) domain->remap(x[i],image[i]);
if (triclinic) domain->x2lamda(atom->nlocal);
domain->reset_box();
Irregular *irregular = new Irregular(lmp);
irregular->migrate_atoms(1);
delete irregular;
if (triclinic) domain->lamda2x(atom->nlocal);
// check that atom IDs are valid
atom->tag_check();
}
/* ----------------------------------------------------------------------
process arg list for dump file fields and optional keywords
------------------------------------------------------------------------- */
int ReadDump::fields_and_keywords(int narg, char **arg)
{
// per-field vectors, leave space for ID and TYPE
fieldtype = new int[narg+2];
fieldlabel = new char*[narg+2];
// add id and type fields as needed
// scan ahead to see if "add yes" keyword/value is used
// requires extra "type" field from from dump file
int iarg;
for (iarg = 0; iarg < narg; iarg++)
if (strcmp(arg[iarg],"add") == 0)
if (iarg < narg-1 && strcmp(arg[iarg+1],"yes") == 0) break;
nfield = 0;
fieldtype[nfield++] = ID;
if (iarg < narg) fieldtype[nfield++] = TYPE;
// parse fields
iarg = 0;
while (iarg < narg) {
int type = whichtype(arg[iarg]);
if (type < 0) break;
if (type == Q && !atom->q_flag)
error->all(FLERR,"Read dump of atom property that isn't allocated");
fieldtype[nfield++] = type;
iarg++;
}
// check for no fields
if (fieldtype[nfield-1] == ID || fieldtype[nfield-1] == TYPE)
error->all(FLERR,"Illegal read_dump command");
if (dimension == 2) {
for (int i = 0; i < nfield; i++)
if (fieldtype[i] == Z || fieldtype[i] == VZ || fieldtype[i] == IZ)
error->all(FLERR,"Illegal read_dump command");
}
for (int i = 0; i < nfield; i++)
for (int j = i+1; j < nfield; j++)
if (fieldtype[i] == fieldtype[j])
error->all(FLERR,"Duplicate fields in read_dump command");
// parse optional args
boxflag = 1;
replaceflag = 1;
purgeflag = 0;
trimflag = 0;
addflag = 0;
for (int i = 0; i < nfield; i++) fieldlabel[i] = NULL;
scaleflag = 0;
wrapflag = 1;
while (iarg < narg) {
if (strcmp(arg[iarg],"box") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
if (strcmp(arg[iarg+1],"yes") == 0) boxflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) boxflag = 0;
else error->all(FLERR,"Illegal read_dump command");
iarg += 2;
} else if (strcmp(arg[iarg],"replace") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
if (strcmp(arg[iarg+1],"yes") == 0) replaceflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) replaceflag = 0;
else error->all(FLERR,"Illegal read_dump command");
iarg += 2;
} else if (strcmp(arg[iarg],"purge") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
if (strcmp(arg[iarg+1],"yes") == 0) purgeflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) purgeflag = 0;
else error->all(FLERR,"Illegal read_dump command");
iarg += 2;
} else if (strcmp(arg[iarg],"trim") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
if (strcmp(arg[iarg+1],"yes") == 0) trimflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) trimflag = 0;
else error->all(FLERR,"Illegal read_dump command");
iarg += 2;
} else if (strcmp(arg[iarg],"add") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
if (strcmp(arg[iarg+1],"yes") == 0) addflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) addflag = 0;
else error->all(FLERR,"Illegal read_dump command");
iarg += 2;
} else if (strcmp(arg[iarg],"label") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal read_dump command");
int type = whichtype(arg[iarg+1]);
int i;
for (i = 0; i < nfield; i++)
if (type == fieldtype[i]) break;
if (i == nfield) error->all(FLERR,"Illegal read_dump command");
int n = strlen(arg[iarg+2]) + 1;
fieldlabel[i] = new char[n];
strcpy(fieldlabel[i],arg[iarg+2]);
iarg += 3;
} else if (strcmp(arg[iarg],"scaled") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
if (strcmp(arg[iarg+1],"yes") == 0) scaleflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) scaleflag = 0;
else error->all(FLERR,"Illegal read_dump command");
iarg += 2;
} else if (strcmp(arg[iarg],"wrapped") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
if (strcmp(arg[iarg+1],"yes") == 0) wrapflag = 1;
else if (strcmp(arg[iarg+1],"no") == 0) wrapflag = 0;
else error->all(FLERR,"Illegal read_dump command");
iarg += 2;
} else if (strcmp(arg[iarg],"format") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal read_dump command");
delete [] readerstyle;
int n = strlen(arg[iarg+1]) + 1;
readerstyle = new char[n];
strcpy(readerstyle,arg[iarg+1]);
iarg += 2;
break;
} else error->all(FLERR,"Illegal read_dump command");
}
if (purgeflag && (replaceflag || trimflag))
error->all(FLERR,"If read_dump purges it cannot replace or trim");
return narg-iarg;
}
/* ----------------------------------------------------------------------
check if str is a field argument
if yes, return index of which
if not, return -1
------------------------------------------------------------------------- */
int ReadDump::whichtype(char *str)
{
int type = -1;
if (strcmp(str,"x") == 0) type = X;
else if (strcmp(str,"y") == 0) type = Y;
else if (strcmp(str,"z") == 0) type = Z;
else if (strcmp(str,"vx") == 0) type = VX;
else if (strcmp(str,"vy") == 0) type = VY;
else if (strcmp(str,"vz") == 0) type = VZ;
else if (strcmp(str,"q") == 0) type = Q;
else if (strcmp(str,"ix") == 0) type = IX;
else if (strcmp(str,"iy") == 0) type = IY;
else if (strcmp(str,"iz") == 0) type = IZ;
return type;
}
/* ----------------------------------------------------------------------
process each of N atoms in chunk read from dump file
if in replace mode and atom ID matches current atom,
overwrite atom info with fields from dump file
if in add mode and atom ID does not match any current atom,
create new atom with dump file field values,
and assign to a proc in round-robin manner
use round-robin method, b/c atom coords may not be inside simulation box
------------------------------------------------------------------------- */
void ReadDump::process_atoms(int n)
{
int i,m,ifield,itype;
int xbox,ybox,zbox;
tagint tag;
double **x = atom->x;
double **v = atom->v;
double *q = atom->q;
imageint *image = atom->image;
int nlocal = atom->nlocal;
tagint map_tag_max = atom->map_tag_max;
for (i = 0; i < n; i++) {
ucflag[i] = 0;
// check if new atom matches one I own
// setting m = -1 forces new atom not to match
// NOTE: atom ID in fields is stored as double, not as ubuf
// so can only cast it to tagint, thus cannot be full 64-bit ID
tag = static_cast<tagint> (fields[i][0]);
if (tag <= map_tag_max) m = atom->map(tag);
else m = -1;
if (m < 0 || m >= nlocal) continue;
ucflag[i] = 1;
uflag[m] = 1;
if (replaceflag) {
nreplace++;
// current image flags
xbox = (image[m] & IMGMASK) - IMGMAX;
ybox = (image[m] >> IMGBITS & IMGMASK) - IMGMAX;
zbox = (image[m] >> IMG2BITS) - IMGMAX;
// overwrite atom attributes with field info
// start from field 1 since 0 = id, 1 will be skipped if type
for (ifield = 1; ifield < nfield; ifield++) {
switch (fieldtype[ifield]) {
case X:
x[m][0] = xfield(i,ifield);
break;
case Y:
x[m][1] = yfield(i,ifield);
break;
case Z:
x[m][2] = zfield(i,ifield);
break;
case VX:
v[m][0] = fields[i][ifield];
break;
case Q:
q[m] = fields[i][ifield];
break;
case VY:
v[m][1] = fields[i][ifield];
break;
case VZ:
v[m][2] = fields[i][ifield];
break;
case IX:
xbox = static_cast<int> (fields[i][ifield]);
break;
case IY:
ybox = static_cast<int> (fields[i][ifield]);
break;
case IZ:
zbox = static_cast<int> (fields[i][ifield]);
break;
}
}
// replace image flag in case changed by ix,iy,iz fields or unwrapping
if (!wrapped) xbox = ybox = zbox = 0;
image[m] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
(((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
(((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
}
}
// create any atoms in chunk that no processor owned
// add atoms in round-robin sequence on processors
// cannot do it geometrically b/c dump coords may not be in simulation box
if (!addflag) return;
MPI_Allreduce(ucflag,ucflag_all,n,MPI_INT,MPI_SUM,world);
int nlocal_previous = atom->nlocal;
double one[3];
for (i = 0; i < n; i++) {
if (ucflag_all[i]) continue;
// each processor adds every Pth atom
addproc++;
if (addproc == nprocs) addproc = 0;
if (addproc != me) continue;
// create type and coord fields from dump file
// coord = 0.0 unless corresponding dump file field was specified
one[0] = one[1] = one[2] = 0.0;
for (ifield = 1; ifield < nfield; ifield++) {
switch (fieldtype[ifield]) {
case TYPE:
itype = static_cast<int> (fields[i][ifield]);
break;
case X:
one[0] = xfield(i,ifield);
break;
case Y:
one[1] = yfield(i,ifield);
break;
case Z:
one[2] = zfield(i,ifield);
break;
}
}
// create the atom on proc that owns it
// reset v,image ptrs in case they are reallocated
m = atom->nlocal;
atom->avec->create_atom(itype,one);
nadd++;
v = atom->v;
q = atom->q;
image = atom->image;
// set atom attributes from other dump file fields
xbox = ybox = zbox = 0;
for (ifield = 1; ifield < nfield; ifield++) {
switch (fieldtype[ifield]) {
case VX:
v[m][0] = fields[i][ifield];
break;
case VY:
v[m][1] = fields[i][ifield];
break;
case VZ:
v[m][2] = fields[i][ifield];
break;
case Q:
q[m] = fields[i][ifield];
break;
case IX:
xbox = static_cast<int> (fields[i][ifield]);
break;
case IY:
ybox = static_cast<int> (fields[i][ifield]);
break;
case IZ:
zbox = static_cast<int> (fields[i][ifield]);
break;
}
// replace image flag in case changed by ix,iy,iz fields
image[m] = ((imageint) (xbox + IMGMAX) & IMGMASK) |
(((imageint) (ybox + IMGMAX) & IMGMASK) << IMGBITS) |
(((imageint) (zbox + IMGMAX) & IMGMASK) << IMG2BITS);
}
}
// invoke set_arrays() for fixes that need initialization of new atoms
// same as in CreateAtoms
nlocal = atom->nlocal;
for (m = 0; m < modify->nfix; m++) {
Fix *fix = modify->fix[m];
if (fix->create_attribute)
for (i = nlocal_previous; i < nlocal; i++)
fix->set_arrays(i);
}
}
/* ----------------------------------------------------------------------
delete atoms not flagged as replaced by dump atoms
------------------------------------------------------------------------- */
void ReadDump::delete_atoms()
{
AtomVec *avec = atom->avec;
int nlocal = atom->nlocal;
int i = 0;
while (i < nlocal) {
if (uflag[i] == 0) {
avec->copy(nlocal-1,i,1);
uflag[i] = uflag[nlocal-1];
nlocal--;
ntrim++;
} else i++;
}
atom->nlocal = nlocal;
}
/* ----------------------------------------------------------------------
convert XYZ fields in dump file into absolute, unscaled coordinates
depends on scaled vs unscaled and triclinic vs orthogonal
does not depend on wrapped vs unwrapped
------------------------------------------------------------------------- */
double ReadDump::xfield(int i, int j)
{
if (!scaled) return fields[i][j];
else if (!triclinic) return fields[i][j]*xprd + xlo;
else if (dimension == 2)
return xprd*fields[i][j] + xy*fields[i][yindex] + xlo;
return xprd*fields[i][j] + xy*fields[i][yindex] + xz*fields[i][zindex] + xlo;
}
double ReadDump::yfield(int i, int j)
{
if (!scaled) return fields[i][j];
else if (!triclinic) return fields[i][j]*yprd + ylo;
else if (dimension == 2) return yprd*fields[i][j] + ylo;
return yprd*fields[i][j] + yz*fields[i][zindex] + ylo;
}
double ReadDump::zfield(int i, int j)
{
if (!scaled) return fields[i][j];
return fields[i][j]*zprd + zlo;
}
diff --git a/src/thermo.cpp b/src/thermo.cpp
index 1a31a1847..e2d131f5b 100644
--- a/src/thermo.cpp
+++ b/src/thermo.cpp
@@ -1,2083 +1,2088 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
+// lmptype.h must be first b/c this file uses MAXBIGINT and includes mpi.h
+// due to OpenMPI bug which sets INT64_MAX via its mpi.h
+// before lmptype.h can set flags to insure it is done correctly
+
+#include "lmptype.h"
#include "mpi.h"
#include "math.h"
#include "stdlib.h"
#include "string.h"
#include "thermo.h"
#include "atom.h"
#include "update.h"
#include "comm.h"
#include "domain.h"
#include "universe.h"
#include "lattice.h"
#include "group.h"
#include "modify.h"
#include "fix.h"
#include "compute.h"
#include "input.h"
#include "variable.h"
#include "neighbor.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "output.h"
#include "timer.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "universe.h"
#include "math_const.h"
using namespace LAMMPS_NS;
using namespace MathConst;
// customize a new keyword by adding to this list:
// step, elapsed, elaplong, dt, time, cpu, tpcpu, spcpu, cpuremain, part
// atoms, temp, press, pe, ke, etotal, enthalpy
// evdwl, ecoul, epair, ebond, eangle, edihed, eimp, emol, elong, etail
// vol, density, lx, ly, lz, xlo, xhi, ylo, yhi, zlo, zhi, xy, xz, yz,
// xlat, ylat, zlat
// bonds, angles, dihedrals, impropers,
// pxx, pyy, pzz, pxy, pxz, pyz
// fmax, fnorm, nbuild, ndanger,
// cella, cellb, cellc, cellalpha, cellbeta, cellgamma
// customize a new thermo style by adding a DEFINE to this list
// also insure allocation of line string is correct in constructor
#define ONE "step temp epair emol etotal press"
#define MULTI "etotal ke temp pe ebond eangle edihed eimp evdwl ecoul elong press"
enum{IGNORE,WARN,ERROR}; // same as several files
enum{ONELINE,MULTILINE};
enum{INT,FLOAT,BIGINT};
enum{SCALAR,VECTOR,ARRAY};
#define INVOKED_SCALAR 1
#define INVOKED_VECTOR 2
#define INVOKED_ARRAY 4
#define DELTA 8
/* ---------------------------------------------------------------------- */
Thermo::Thermo(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
{
MPI_Comm_rank(world,&me);
int n = strlen(arg[0]) + 1;
style = new char[n];
strcpy(style,arg[0]);
// set thermo_modify defaults
modified = 0;
normuserflag = 0;
lineflag = ONELINE;
lostflag = lostbond = ERROR;
lostbefore = 0;
flushflag = 0;
// set style and corresponding lineflag
// custom style builds its own line of keywords
// customize a new thermo style by adding to if statement
// allocate line string used for 3 tasks
// concat of custom style args
// one-time thermo output of header line
// each line of numeric thermo output
// 256 = extra for ONE or MULTI string or multi formatting
// 64 = max per-arg chars in header or numeric output
if (strcmp(style,"one") == 0) {
line = new char[256+6*64];
strcpy(line,ONE);
} else if (strcmp(style,"multi") == 0) {
line = new char[256+12*64];
strcpy(line,MULTI);
lineflag = MULTILINE;
} else if (strcmp(style,"custom") == 0) {
if (narg == 1) error->all(FLERR,"Illegal thermo style custom command");
line = new char[256+narg*64];
line[0] = '\0';
for (int iarg = 1; iarg < narg; iarg++) {
strcat(line,arg[iarg]);
strcat(line," ");
}
line[strlen(line)-1] = '\0';
} else error->all(FLERR,"Illegal thermo style command");
// ptrs, flags, IDs for compute objects thermo may use or create
temperature = NULL;
pressure = NULL;
pe = NULL;
index_temp = index_press_scalar = index_press_vector = index_pe = -1;
id_temp = (char *) "thermo_temp";
id_press = (char *) "thermo_press";
id_pe = (char *) "thermo_pe";
// count fields in line
// allocate per-field memory
// process line of keywords
nfield_initial = atom->count_words(line);
allocate();
parse_fields(line);
// format strings
char *bigint_format = (char *) BIGINT_FORMAT;
char *fformat_multi = (char *) "---------------- Step %%8%s ----- "
"CPU = %%11.4f (sec) ----------------";
sprintf(format_multi,fformat_multi,&bigint_format[1]);
format_float_one_def = (char *) "%12.8g";
format_float_multi_def = (char *) "%14.4f";
format_int_one_def = (char *) "%8d";
format_int_multi_def = (char *) "%14d";
sprintf(format_bigint_one_def,"%%8%s",&bigint_format[1]);
sprintf(format_bigint_multi_def,"%%14%s",&bigint_format[1]);
format_float_user = NULL;
format_int_user = NULL;
format_bigint_user = NULL;
}
/* ---------------------------------------------------------------------- */
Thermo::~Thermo()
{
delete [] style;
delete [] line;
deallocate();
// format strings
delete [] format_float_user;
delete [] format_int_user;
delete [] format_bigint_user;
}
/* ---------------------------------------------------------------------- */
void Thermo::init()
{
int i,n;
// set normvalue to default setting unless user has specified it
if (normuserflag) normvalue = normuser;
else if (strcmp(update->unit_style,"lj") == 0) normvalue = 1;
else normvalue = 0;
// add Volume field if volume changes and not style = custom
// this check must come after domain init, so box_change is set
nfield = nfield_initial;
if (domain->box_change && strcmp(style,"custom") != 0)
addfield("Volume",&Thermo::compute_vol,FLOAT);
// set format string for each field
// include keyword if lineflag = MULTILINE
// add '/n' every 3 values if lineflag = MULTILINE
// add trailing '/n' to last value
char *ptr;
for (i = 0; i < nfield; i++) {
format[i][0] = '\0';
if (lineflag == MULTILINE && i % 3 == 0) strcat(format[i],"\n");
if (format_user[i]) ptr = format_user[i];
else if (vtype[i] == FLOAT) {
if (format_float_user) ptr = format_float_user;
else if (lineflag == ONELINE) ptr = format_float_one_def;
else if (lineflag == MULTILINE) ptr = format_float_multi_def;
} else if (vtype[i] == INT) {
if (format_int_user) ptr = format_int_user;
else if (lineflag == ONELINE) ptr = format_int_one_def;
else if (lineflag == MULTILINE) ptr = format_int_multi_def;
} else if (vtype[i] == BIGINT) {
if (format_bigint_user) ptr = format_bigint_user;
else if (lineflag == ONELINE) ptr = format_bigint_one_def;
else if (lineflag == MULTILINE) ptr = format_bigint_multi_def;
}
n = strlen(format[i]);
if (lineflag == ONELINE) sprintf(&format[i][n],"%s ",ptr);
else sprintf(&format[i][n],"%-8s = %s ",keyword[i],ptr);
if (i == nfield-1) strcat(format[i],"\n");
}
// find current ptr for each Compute ID
// cudable = 0 if any compute used by Thermo is non-CUDA
cudable = 1;
int icompute;
for (i = 0; i < ncompute; i++) {
icompute = modify->find_compute(id_compute[i]);
if (icompute < 0) error->all(FLERR,"Could not find thermo compute ID");
computes[i] = modify->compute[icompute];
cudable = cudable && computes[i]->cudable;
}
// find current ptr for each Fix ID
// check that fix frequency is acceptable with thermo output frequency
int ifix;
for (i = 0; i < nfix; i++) {
ifix = modify->find_fix(id_fix[i]);
if (ifix < 0) error->all(FLERR,"Could not find thermo fix ID");
fixes[i] = modify->fix[ifix];
if (output->thermo_every % fixes[i]->global_freq)
error->all(FLERR,"Thermo and fix not computed at compatible times");
}
// find current ptr for each Variable ID
int ivariable;
for (i = 0; i < nvariable; i++) {
ivariable = input->variable->find(id_variable[i]);
if (ivariable < 0)
error->all(FLERR,"Could not find thermo variable name");
variables[i] = ivariable;
}
// set ptrs to keyword-specific Compute objects
if (index_temp >= 0) temperature = computes[index_temp];
if (index_press_scalar >= 0) pressure = computes[index_press_scalar];
if (index_press_vector >= 0) pressure = computes[index_press_vector];
if (index_pe >= 0) pe = computes[index_pe];
}
/* ---------------------------------------------------------------------- */
void Thermo::header()
{
if (lineflag == MULTILINE) return;
int loc = 0;
for (int i = 0; i < nfield; i++)
loc += sprintf(&line[loc],"%s ",keyword[i]);
sprintf(&line[loc],"\n");
if (me == 0) {
if (screen) fprintf(screen,"%s",line);
if (logfile) fprintf(logfile,"%s",line);
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute(int flag)
{
int i;
firststep = flag;
bigint ntimestep = update->ntimestep;
// check for lost atoms
// turn off normflag if natoms = 0 to avoid divide by 0
natoms = lost_check();
if (natoms == 0) normflag = 0;
else normflag = normvalue;
// invoke Compute methods needed for thermo keywords
for (i = 0; i < ncompute; i++)
if (compute_which[i] == SCALAR) {
if (!(computes[i]->invoked_flag & INVOKED_SCALAR)) {
computes[i]->compute_scalar();
computes[i]->invoked_flag |= INVOKED_SCALAR;
}
} else if (compute_which[i] == VECTOR) {
if (!(computes[i]->invoked_flag & INVOKED_VECTOR)) {
computes[i]->compute_vector();
computes[i]->invoked_flag |= INVOKED_VECTOR;
}
} else if (compute_which[i] == ARRAY) {
if (!(computes[i]->invoked_flag & INVOKED_ARRAY)) {
computes[i]->compute_array();
computes[i]->invoked_flag |= INVOKED_ARRAY;
}
}
// if lineflag = MULTILINE, prepend step/cpu header line
int loc = 0;
if (lineflag == MULTILINE) {
double cpu;
if (flag) cpu = timer->elapsed(Timer::TOTAL);
else cpu = 0.0;
loc = sprintf(&line[loc],format_multi,ntimestep,cpu);
}
// add each thermo value to line with its specific format
for (ifield = 0; ifield < nfield; ifield++) {
(this->*vfunc[ifield])();
if (vtype[ifield] == FLOAT)
loc += sprintf(&line[loc],format[ifield],dvalue);
else if (vtype[ifield] == INT)
loc += sprintf(&line[loc],format[ifield],ivalue);
else if (vtype[ifield] == BIGINT) {
loc += sprintf(&line[loc],format[ifield],bivalue);
}
}
// print line to screen and logfile
if (me == 0) {
if (screen) fprintf(screen,"%s",line);
if (logfile) {
fprintf(logfile,"%s",line);
if (flushflag) fflush(logfile);
}
}
}
/* ----------------------------------------------------------------------
check for lost atoms, return current number of atoms
------------------------------------------------------------------------- */
bigint Thermo::lost_check()
{
// ntotal = current # of atoms
bigint ntotal;
bigint nblocal = atom->nlocal;
MPI_Allreduce(&nblocal,&ntotal,1,MPI_LMP_BIGINT,MPI_SUM,world);
if (ntotal < 0 || ntotal > MAXBIGINT)
error->all(FLERR,"Too many total atoms");
if (ntotal == atom->natoms) return ntotal;
// if not checking or already warned, just return
// reset total atom count
if (lostflag == IGNORE) return ntotal;
if (lostflag == WARN && lostbefore == 1) {
atom->natoms = ntotal;
return ntotal;
}
// error message
if (lostflag == ERROR) {
char str[64];
sprintf(str,
"Lost atoms: original " BIGINT_FORMAT " current " BIGINT_FORMAT,
atom->natoms,ntotal);
error->all(FLERR,str);
}
// warning message
char str[64];
sprintf(str,
"Lost atoms: original " BIGINT_FORMAT " current " BIGINT_FORMAT,
atom->natoms,ntotal);
if (me == 0) error->warning(FLERR,str,0);
// reset total atom count
atom->natoms = ntotal;
lostbefore = 1;
return ntotal;
}
/* ----------------------------------------------------------------------
modify thermo parameters
------------------------------------------------------------------------- */
void Thermo::modify_params(int narg, char **arg)
{
if (narg == 0) error->all(FLERR,"Illegal thermo_modify command");
modified = 1;
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"temp") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
if (index_temp < 0) error->all(FLERR,"Thermo style does not use temp");
delete [] id_compute[index_temp];
int n = strlen(arg[iarg+1]) + 1;
id_compute[index_temp] = new char[n];
strcpy(id_compute[index_temp],arg[iarg+1]);
int icompute = modify->find_compute(arg[iarg+1]);
if (icompute < 0)
error->all(FLERR,"Could not find thermo_modify temperature ID");
temperature = modify->compute[icompute];
if (temperature->tempflag == 0)
error->all(FLERR,"Thermo_modify temperature ID does not "
"compute temperature");
if (temperature->igroup != 0 && comm->me == 0)
error->warning(FLERR,
"Temperature for thermo pressure is not for group all");
// reset id_temp of pressure to new temperature ID
// either pressure currently being used by thermo or "thermo_press"
if (index_press_scalar >= 0) {
icompute = modify->find_compute(id_compute[index_press_scalar]);
if (icompute < 0) error->all(FLERR,
"Pressure ID for thermo does not exist");
} else if (index_press_vector >= 0) {
icompute = modify->find_compute(id_compute[index_press_vector]);
if (icompute < 0) error->all(FLERR,
"Pressure ID for thermo does not exist");
} else icompute = modify->find_compute((char *) "thermo_press");
modify->compute[icompute]->reset_extra_compute_fix(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"press") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
if (index_press_scalar < 0 && index_press_vector < 0)
error->all(FLERR,"Thermo style does not use press");
if (index_press_scalar >= 0) {
delete [] id_compute[index_press_scalar];
int n = strlen(arg[iarg+1]) + 1;
id_compute[index_press_scalar] = new char[n];
strcpy(id_compute[index_press_scalar],arg[iarg+1]);
}
if (index_press_vector >= 0) {
delete [] id_compute[index_press_vector];
int n = strlen(arg[iarg+1]) + 1;
id_compute[index_press_vector] = new char[n];
strcpy(id_compute[index_press_vector],arg[iarg+1]);
}
int icompute = modify->find_compute(arg[iarg+1]);
if (icompute < 0) error->all(FLERR,
"Could not find thermo_modify pressure ID");
pressure = modify->compute[icompute];
if (pressure->pressflag == 0)
error->all(FLERR,"Thermo_modify pressure ID does not compute pressure");
iarg += 2;
} else if (strcmp(arg[iarg],"lost") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
if (strcmp(arg[iarg+1],"ignore") == 0) lostflag = IGNORE;
else if (strcmp(arg[iarg+1],"warn") == 0) lostflag = WARN;
else if (strcmp(arg[iarg+1],"error") == 0) lostflag = ERROR;
else error->all(FLERR,"Illegal thermo_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"lost/bond") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
if (strcmp(arg[iarg+1],"ignore") == 0) lostbond = IGNORE;
else if (strcmp(arg[iarg+1],"warn") == 0) lostbond = WARN;
else if (strcmp(arg[iarg+1],"error") == 0) lostbond = ERROR;
else error->all(FLERR,"Illegal thermo_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"norm") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
normuserflag = 1;
if (strcmp(arg[iarg+1],"no") == 0) normuser = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) normuser = 1;
else error->all(FLERR,"Illegal thermo_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"flush") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
if (strcmp(arg[iarg+1],"no") == 0) flushflag = 0;
else if (strcmp(arg[iarg+1],"yes") == 0) flushflag = 1;
else error->all(FLERR,"Illegal thermo_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"line") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal thermo_modify command");
if (strcmp(arg[iarg+1],"one") == 0) lineflag = ONELINE;
else if (strcmp(arg[iarg+1],"multi") == 0) lineflag = MULTILINE;
else error->all(FLERR,"Illegal thermo_modify command");
iarg += 2;
} else if (strcmp(arg[iarg],"format") == 0) {
if (iarg+3 > narg) error->all(FLERR,"Illegal thermo_modify command");
if (strcmp(arg[iarg+1],"int") == 0) {
if (format_int_user) delete [] format_int_user;
int n = strlen(arg[iarg+2]) + 1;
format_int_user = new char[n];
strcpy(format_int_user,arg[iarg+2]);
if (format_bigint_user) delete [] format_bigint_user;
n = strlen(format_int_user) + 3;
format_bigint_user = new char[n];
char *ptr = strchr(format_int_user,'d');
if (ptr == NULL)
error->all(FLERR,
"Thermo_modify int format does not contain d character");
*ptr = '\0';
sprintf(format_bigint_user,"%s%s%s",format_int_user,
BIGINT_FORMAT,ptr+1);
*ptr = 'd';
} else if (strcmp(arg[iarg+1],"float") == 0) {
if (format_float_user) delete [] format_float_user;
int n = strlen(arg[iarg+2]) + 1;
format_float_user = new char[n];
strcpy(format_float_user,arg[iarg+2]);
} else {
int i = force->inumeric(FLERR,arg[iarg+1]) - 1;
if (i < 0 || i >= nfield_initial)
error->all(FLERR,"Illegal thermo_modify command");
if (format_user[i]) delete [] format_user[i];
int n = strlen(arg[iarg+2]) + 1;
format_user[i] = new char[n];
strcpy(format_user[i],arg[iarg+2]);
}
iarg += 3;
} else error->all(FLERR,"Illegal thermo_modify command");
}
}
/* ----------------------------------------------------------------------
allocate all per-field memory
------------------------------------------------------------------------- */
void Thermo::allocate()
{
// n = specified fields + Volume field (added at run time)
int n = nfield_initial + 1;
keyword = new char*[n];
for (int i = 0; i < n; i++) keyword[i] = new char[32];
vfunc = new FnPtr[n];
vtype = new int[n];
format = new char*[n];
for (int i = 0; i < n; i++) format[i] = new char[32];
format_user = new char*[n];
for (int i = 0; i < n; i++) format_user[i] = NULL;
field2index = new int[n];
argindex1 = new int[n];
argindex2 = new int[n];
// factor of 3 is max number of computes a single field can add
ncompute = 0;
id_compute = new char*[3*n];
compute_which = new int[3*n];
computes = new Compute*[3*n];
nfix = 0;
id_fix = new char*[n];
fixes = new Fix*[n];
nvariable = 0;
id_variable = new char*[n];
variables = new int[n];
}
/* ----------------------------------------------------------------------
deallocate all per-field memory
------------------------------------------------------------------------- */
void Thermo::deallocate()
{
int n = nfield_initial + 1;
for (int i = 0; i < n; i++) delete [] keyword[i];
delete [] keyword;
delete [] vfunc;
delete [] vtype;
for (int i = 0; i < n; i++) delete [] format[i];
delete [] format;
for (int i = 0; i < n; i++) delete [] format_user[i];
delete [] format_user;
delete [] field2index;
delete [] argindex1;
delete [] argindex2;
for (int i = 0; i < ncompute; i++) delete [] id_compute[i];
delete [] id_compute;
delete [] compute_which;
delete [] computes;
for (int i = 0; i < nfix; i++) delete [] id_fix[i];
delete [] id_fix;
delete [] fixes;
for (int i = 0; i < nvariable; i++) delete [] id_variable[i];
delete [] id_variable;
delete [] variables;
}
/* ----------------------------------------------------------------------
parse list of thermo keywords from str
set compute flags (temp, press, pe, etc)
------------------------------------------------------------------------- */
void Thermo::parse_fields(char *str)
{
nfield = 0;
// customize a new keyword by adding to if statement
char *word = strtok(str," \0");
while (word) {
if (strcmp(word,"step") == 0) {
addfield("Step",&Thermo::compute_step,BIGINT);
} else if (strcmp(word,"elapsed") == 0) {
addfield("Elapsed",&Thermo::compute_elapsed,BIGINT);
} else if (strcmp(word,"elaplong") == 0) {
addfield("Elaplong",&Thermo::compute_elapsed_long,BIGINT);
} else if (strcmp(word,"dt") == 0) {
addfield("Dt",&Thermo::compute_dt,FLOAT);
} else if (strcmp(word,"time") == 0) {
addfield("Time",&Thermo::compute_time,FLOAT);
} else if (strcmp(word,"cpu") == 0) {
addfield("CPU",&Thermo::compute_cpu,FLOAT);
} else if (strcmp(word,"tpcpu") == 0) {
addfield("T/CPU",&Thermo::compute_tpcpu,FLOAT);
} else if (strcmp(word,"spcpu") == 0) {
addfield("S/CPU",&Thermo::compute_spcpu,FLOAT);
} else if (strcmp(word,"cpuremain") == 0) {
addfield("CPULeft",&Thermo::compute_cpuremain,FLOAT);
} else if (strcmp(word,"part") == 0) {
addfield("Part",&Thermo::compute_part,INT);
} else if (strcmp(word,"atoms") == 0) {
addfield("Atoms",&Thermo::compute_atoms,BIGINT);
} else if (strcmp(word,"temp") == 0) {
addfield("Temp",&Thermo::compute_temp,FLOAT);
index_temp = add_compute(id_temp,SCALAR);
} else if (strcmp(word,"press") == 0) {
addfield("Press",&Thermo::compute_press,FLOAT);
index_press_scalar = add_compute(id_press,SCALAR);
} else if (strcmp(word,"pe") == 0) {
addfield("PotEng",&Thermo::compute_pe,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"ke") == 0) {
addfield("KinEng",&Thermo::compute_ke,FLOAT);
index_temp = add_compute(id_temp,SCALAR);
} else if (strcmp(word,"etotal") == 0) {
addfield("TotEng",&Thermo::compute_etotal,FLOAT);
index_temp = add_compute(id_temp,SCALAR);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"enthalpy") == 0) {
addfield("Enthalpy",&Thermo::compute_enthalpy,FLOAT);
index_temp = add_compute(id_temp,SCALAR);
index_press_scalar = add_compute(id_press,SCALAR);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"evdwl") == 0) {
addfield("E_vdwl",&Thermo::compute_evdwl,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"ecoul") == 0) {
addfield("E_coul",&Thermo::compute_ecoul,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"epair") == 0) {
addfield("E_pair",&Thermo::compute_epair,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"ebond") == 0) {
addfield("E_bond",&Thermo::compute_ebond,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"eangle") == 0) {
addfield("E_angle",&Thermo::compute_eangle,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"edihed") == 0) {
addfield("E_dihed",&Thermo::compute_edihed,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"eimp") == 0) {
addfield("E_impro",&Thermo::compute_eimp,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"emol") == 0) {
addfield("E_mol",&Thermo::compute_emol,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"elong") == 0) {
addfield("E_long",&Thermo::compute_elong,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"etail") == 0) {
addfield("E_tail",&Thermo::compute_etail,FLOAT);
index_pe = add_compute(id_pe,SCALAR);
} else if (strcmp(word,"vol") == 0) {
addfield("Volume",&Thermo::compute_vol,FLOAT);
} else if (strcmp(word,"density") == 0) {
addfield("Density",&Thermo::compute_density,FLOAT);
} else if (strcmp(word,"lx") == 0) {
addfield("Lx",&Thermo::compute_lx,FLOAT);
} else if (strcmp(word,"ly") == 0) {
addfield("Ly",&Thermo::compute_ly,FLOAT);
} else if (strcmp(word,"lz") == 0) {
addfield("Lz",&Thermo::compute_lz,FLOAT);
} else if (strcmp(word,"xlo") == 0) {
addfield("Xlo",&Thermo::compute_xlo,FLOAT);
} else if (strcmp(word,"xhi") == 0) {
addfield("Xhi",&Thermo::compute_xhi,FLOAT);
} else if (strcmp(word,"ylo") == 0) {
addfield("Ylo",&Thermo::compute_ylo,FLOAT);
} else if (strcmp(word,"yhi") == 0) {
addfield("Yhi",&Thermo::compute_yhi,FLOAT);
} else if (strcmp(word,"zlo") == 0) {
addfield("Zlo",&Thermo::compute_zlo,FLOAT);
} else if (strcmp(word,"zhi") == 0) {
addfield("Zhi",&Thermo::compute_zhi,FLOAT);
} else if (strcmp(word,"xy") == 0) {
addfield("Xy",&Thermo::compute_xy,FLOAT);
} else if (strcmp(word,"xz") == 0) {
addfield("Xz",&Thermo::compute_xz,FLOAT);
} else if (strcmp(word,"yz") == 0) {
addfield("Yz",&Thermo::compute_yz,FLOAT);
} else if (strcmp(word,"xlat") == 0) {
addfield("Xlat",&Thermo::compute_xlat,FLOAT);
} else if (strcmp(word,"ylat") == 0) {
addfield("Ylat",&Thermo::compute_ylat,FLOAT);
} else if (strcmp(word,"zlat") == 0) {
addfield("Zlat",&Thermo::compute_zlat,FLOAT);
} else if (strcmp(word,"bonds") == 0) {
addfield("Bonds",&Thermo::compute_bonds,BIGINT);
} else if (strcmp(word,"angles") == 0) {
addfield("Angles",&Thermo::compute_angles,BIGINT);
} else if (strcmp(word,"dihedrals") == 0) {
addfield("Diheds",&Thermo::compute_dihedrals,BIGINT);
} else if (strcmp(word,"impropers") == 0) {
addfield("Impros",&Thermo::compute_impropers,BIGINT);
} else if (strcmp(word,"pxx") == 0) {
addfield("Pxx",&Thermo::compute_pxx,FLOAT);
index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pyy") == 0) {
addfield("Pyy",&Thermo::compute_pyy,FLOAT);
index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pzz") == 0) {
addfield("Pzz",&Thermo::compute_pzz,FLOAT);
index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pxy") == 0) {
addfield("Pxy",&Thermo::compute_pxy,FLOAT);
index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pxz") == 0) {
addfield("Pxz",&Thermo::compute_pxz,FLOAT);
index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"pyz") == 0) {
addfield("Pyz",&Thermo::compute_pyz,FLOAT);
index_press_vector = add_compute(id_press,VECTOR);
} else if (strcmp(word,"fmax") == 0) {
addfield("Fmax",&Thermo::compute_fmax,FLOAT);
} else if (strcmp(word,"fnorm") == 0) {
addfield("Fnorm",&Thermo::compute_fnorm,FLOAT);
} else if (strcmp(word,"nbuild") == 0) {
addfield("Nbuild",&Thermo::compute_nbuild,BIGINT);
} else if (strcmp(word,"ndanger") == 0) {
addfield("Ndanger",&Thermo::compute_ndanger,BIGINT);
} else if (strcmp(word,"cella") == 0) {
addfield("Cella",&Thermo::compute_cella,FLOAT);
} else if (strcmp(word,"cellb") == 0) {
addfield("Cellb",&Thermo::compute_cellb,FLOAT);
} else if (strcmp(word,"cellc") == 0) {
addfield("Cellc",&Thermo::compute_cellc,FLOAT);
} else if (strcmp(word,"cellalpha") == 0) {
addfield("CellAlpha",&Thermo::compute_cellalpha,FLOAT);
} else if (strcmp(word,"cellbeta") == 0) {
addfield("CellBeta",&Thermo::compute_cellbeta,FLOAT);
} else if (strcmp(word,"cellgamma") == 0) {
addfield("CellGamma",&Thermo::compute_cellgamma,FLOAT);
// compute value = c_ID, fix value = f_ID, variable value = v_ID
// count trailing [] and store int arguments
// copy = at most 8 chars of ID to pass to addfield
} else if ((strncmp(word,"c_",2) == 0) || (strncmp(word,"f_",2) == 0) ||
(strncmp(word,"v_",2) == 0)) {
int n = strlen(word);
char *id = new char[n];
strcpy(id,&word[2]);
char copy[9];
strncpy(copy,id,8);
copy[8] = '\0';
// parse zero or one or two trailing brackets from ID
// argindex1,argindex2 = int inside each bracket pair, 0 if no bracket
char *ptr = strchr(id,'[');
if (ptr == NULL) argindex1[nfield] = 0;
else {
*ptr = '\0';
argindex1[nfield] = input->variable->int_between_brackets(ptr,0);
ptr++;
if (*ptr == '[') {
argindex2[nfield] = input->variable->int_between_brackets(ptr,0);
ptr++;
} else argindex2[nfield] = 0;
}
if (word[0] == 'c') {
n = modify->find_compute(id);
if (n < 0) error->all(FLERR,"Could not find thermo custom compute ID");
if (argindex1[nfield] == 0 && modify->compute[n]->scalar_flag == 0)
error->all(FLERR,"Thermo compute does not compute scalar");
if (argindex1[nfield] > 0 && argindex2[nfield] == 0) {
if (modify->compute[n]->vector_flag == 0)
error->all(FLERR,"Thermo compute does not compute vector");
if (argindex1[nfield] > modify->compute[n]->size_vector)
error->all(FLERR,"Thermo compute vector is accessed out-of-range");
}
if (argindex1[nfield] > 0 && argindex2[nfield] > 0) {
if (modify->compute[n]->array_flag == 0)
error->all(FLERR,"Thermo compute does not compute array");
if (argindex1[nfield] > modify->compute[n]->size_array_rows ||
argindex2[nfield] > modify->compute[n]->size_array_cols)
error->all(FLERR,"Thermo compute array is accessed out-of-range");
}
if (argindex1[nfield] == 0)
field2index[nfield] = add_compute(id,SCALAR);
else if (argindex2[nfield] == 0)
field2index[nfield] = add_compute(id,VECTOR);
else
field2index[nfield] = add_compute(id,ARRAY);
addfield(copy,&Thermo::compute_compute,FLOAT);
} else if (word[0] == 'f') {
n = modify->find_fix(id);
if (n < 0) error->all(FLERR,"Could not find thermo custom fix ID");
if (argindex1[nfield] == 0 && modify->fix[n]->scalar_flag == 0)
error->all(FLERR,"Thermo fix does not compute scalar");
if (argindex1[nfield] > 0 && argindex2[nfield] == 0) {
if (modify->fix[n]->vector_flag == 0)
error->all(FLERR,"Thermo fix does not compute vector");
if (argindex1[nfield] > modify->fix[n]->size_vector)
error->all(FLERR,"Thermo fix vector is accessed out-of-range");
}
if (argindex1[nfield] > 0 && argindex2[nfield] > 0) {
if (modify->fix[n]->array_flag == 0)
error->all(FLERR,"Thermo fix does not compute array");
if (argindex1[nfield] > modify->fix[n]->size_array_rows ||
argindex2[nfield] > modify->fix[n]->size_array_cols)
error->all(FLERR,"Thermo fix array is accessed out-of-range");
}
field2index[nfield] = add_fix(id);
addfield(copy,&Thermo::compute_fix,FLOAT);
} else if (word[0] == 'v') {
n = input->variable->find(id);
if (n < 0)
error->all(FLERR,"Could not find thermo custom variable name");
if (input->variable->equalstyle(n) == 0)
error->all(FLERR,
"Thermo custom variable is not equal-style variable");
if (argindex1[nfield])
error->all(FLERR,"Thermo custom variable cannot be indexed");
field2index[nfield] = add_variable(id);
addfield(copy,&Thermo::compute_variable,FLOAT);
}
delete [] id;
} else error->all(FLERR,"Invalid keyword in thermo_style custom command");
word = strtok(NULL," \0");
}
}
/* ----------------------------------------------------------------------
add field to list of quantities to print
------------------------------------------------------------------------- */
void Thermo::addfield(const char *key, FnPtr func, int typeflag)
{
strcpy(keyword[nfield],key);
vfunc[nfield] = func;
vtype[nfield] = typeflag;
nfield++;
}
/* ----------------------------------------------------------------------
add compute ID to list of Compute objects to call
return location of where this Compute is in list
if already in list with same which, do not add, just return index
------------------------------------------------------------------------- */
int Thermo::add_compute(const char *id, int which)
{
int icompute;
for (icompute = 0; icompute < ncompute; icompute++)
if ((strcmp(id,id_compute[icompute]) == 0) &&
which == compute_which[icompute]) break;
if (icompute < ncompute) return icompute;
int n = strlen(id) + 1;
id_compute[ncompute] = new char[n];
strcpy(id_compute[ncompute],id);
compute_which[ncompute] = which;
ncompute++;
return ncompute-1;
}
/* ----------------------------------------------------------------------
add fix ID to list of Fix objects to call
------------------------------------------------------------------------- */
int Thermo::add_fix(const char *id)
{
int n = strlen(id) + 1;
id_fix[nfix] = new char[n];
strcpy(id_fix[nfix],id);
nfix++;
return nfix-1;
}
/* ----------------------------------------------------------------------
add variable ID to list of Variables to evaluate
------------------------------------------------------------------------- */
int Thermo::add_variable(const char *id)
{
int n = strlen(id) + 1;
id_variable[nvariable] = new char[n];
strcpy(id_variable[nvariable],id);
nvariable++;
return nvariable-1;
}
/* ----------------------------------------------------------------------
compute a single thermodyanmic value, word is any keyword in custom list
called when a variable is evaluated by Variable class
return value as double in answer
return 0 if str is recoginzed keyword, 1 if unrecognized
customize a new keyword by adding to if statement
------------------------------------------------------------------------- */
int Thermo::evaluate_keyword(char *word, double *answer)
{
// turn off normflag if natoms = 0 to avoid divide by 0
// normflag must be set for lo-level thermo routines that may be invoked
natoms = atom->natoms;
if (natoms == 0) normflag = 0;
else normflag = normvalue;
// invoke a lo-level thermo routine to compute the variable value
// if keyword requires a compute, error if thermo doesn't use the compute
// if inbetween runs and needed compute is not current, error
// if in middle of run and needed compute is not current, invoke it
// for keywords that use pe indirectly (evdwl, ebond, etc):
// check if energy was tallied on this timestep and set pe->invoked_flag
// this will trigger next timestep for energy tallying via addstep()
if (strcmp(word,"step") == 0) {
compute_step();
dvalue = bivalue;
} else if (strcmp(word,"elapsed") == 0) {
if (update->whichflag == 0)
error->all(FLERR,
"This variable thermo keyword cannot be used between runs");
compute_elapsed();
dvalue = bivalue;
} else if (strcmp(word,"elaplong") == 0) {
if (update->whichflag == 0)
error->all(FLERR,
"This variable thermo keyword cannot be used between runs");
compute_elapsed_long();
dvalue = bivalue;
} else if (strcmp(word,"dt") == 0) {
compute_dt();
} else if (strcmp(word,"time") == 0) {
compute_time();
} else if (strcmp(word,"cpu") == 0) {
if (update->whichflag == 0)
error->all(FLERR,
"This variable thermo keyword cannot be used between runs");
compute_cpu();
} else if (strcmp(word,"tpcpu") == 0) {
if (update->whichflag == 0)
error->all(FLERR,
"This variable thermo keyword cannot be used between runs");
compute_tpcpu();
} else if (strcmp(word,"spcpu") == 0) {
if (update->whichflag == 0)
error->all(FLERR,
"This variable thermo keyword cannot be used between runs");
compute_spcpu();
} else if (strcmp(word,"cpuremain") == 0) {
if (update->whichflag == 0)
error->all(FLERR,
"This variable thermo keyword cannot be used between runs");
compute_cpuremain();
} else if (strcmp(word,"part") == 0) {
compute_part();
dvalue = ivalue;
} else if (strcmp(word,"atoms") == 0) {
compute_atoms();
dvalue = bivalue;
} else if (strcmp(word,"temp") == 0) {
if (!temperature)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init temp");
if (update->whichflag == 0) {
if (temperature->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(temperature->invoked_flag & INVOKED_SCALAR)) {
temperature->compute_scalar();
temperature->invoked_flag |= INVOKED_SCALAR;
}
compute_temp();
} else if (strcmp(word,"press") == 0) {
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
if (update->whichflag == 0) {
if (pressure->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pressure->invoked_flag & INVOKED_SCALAR)) {
pressure->compute_scalar();
pressure->invoked_flag |= INVOKED_SCALAR;
}
compute_press();
} else if (strcmp(word,"pe") == 0) {
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
if (update->whichflag == 0) {
if (pe->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pe->invoked_flag & INVOKED_SCALAR)) {
pe->compute_scalar();
pe->invoked_flag |= INVOKED_SCALAR;
}
compute_pe();
} else if (strcmp(word,"ke") == 0) {
if (!temperature)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init temp");
if (update->whichflag == 0) {
if (temperature->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(temperature->invoked_flag & INVOKED_SCALAR)) {
temperature->compute_scalar();
temperature->invoked_flag |= INVOKED_SCALAR;
}
compute_ke();
} else if (strcmp(word,"etotal") == 0) {
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
if (update->whichflag == 0) {
if (pe->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pe->invoked_flag & INVOKED_SCALAR)) {
pe->compute_scalar();
pe->invoked_flag |= INVOKED_SCALAR;
}
if (!temperature)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init temp");
if (update->whichflag == 0) {
if (temperature->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(temperature->invoked_flag & INVOKED_SCALAR)) {
temperature->compute_scalar();
temperature->invoked_flag |= INVOKED_SCALAR;
}
compute_etotal();
} else if (strcmp(word,"enthalpy") == 0) {
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
if (update->whichflag == 0) {
if (pe->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pe->invoked_flag & INVOKED_SCALAR)) {
pe->compute_scalar();
pe->invoked_flag |= INVOKED_SCALAR;
}
if (!temperature)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init temp");
if (update->whichflag == 0) {
if (temperature->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(temperature->invoked_flag & INVOKED_SCALAR)) {
temperature->compute_scalar();
temperature->invoked_flag |= INVOKED_SCALAR;
}
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
if (update->whichflag == 0) {
if (pressure->invoked_scalar != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pressure->invoked_flag & INVOKED_SCALAR)) {
pressure->compute_scalar();
pressure->invoked_flag |= INVOKED_SCALAR;
}
compute_enthalpy();
} else if (strcmp(word,"evdwl") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_evdwl();
} else if (strcmp(word,"ecoul") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_ecoul();
} else if (strcmp(word,"epair") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_epair();
} else if (strcmp(word,"ebond") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_ebond();
} else if (strcmp(word,"eangle") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_eangle();
} else if (strcmp(word,"edihed") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_edihed();
} else if (strcmp(word,"eimp") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_eimp();
} else if (strcmp(word,"emol") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_emol();
} else if (strcmp(word,"elong") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_elong();
} else if (strcmp(word,"etail") == 0) {
if (update->eflag_global != update->ntimestep)
error->all(FLERR,"Energy was not tallied on needed timestep");
if (!pe)
error->all(FLERR,
"Thermo keyword in variable requires thermo to use/init pe");
pe->invoked_flag |= INVOKED_SCALAR;
compute_etail();
} else if (strcmp(word,"vol") == 0) compute_vol();
else if (strcmp(word,"density") == 0) compute_density();
else if (strcmp(word,"lx") == 0) compute_lx();
else if (strcmp(word,"ly") == 0) compute_ly();
else if (strcmp(word,"lz") == 0) compute_lz();
else if (strcmp(word,"xlo") == 0) compute_xlo();
else if (strcmp(word,"xhi") == 0) compute_xhi();
else if (strcmp(word,"ylo") == 0) compute_ylo();
else if (strcmp(word,"yhi") == 0) compute_yhi();
else if (strcmp(word,"zlo") == 0) compute_zlo();
else if (strcmp(word,"zhi") == 0) compute_zhi();
else if (strcmp(word,"xy") == 0) compute_xy();
else if (strcmp(word,"xz") == 0) compute_xz();
else if (strcmp(word,"yz") == 0) compute_yz();
else if (strcmp(word,"xlat") == 0) compute_xlat();
else if (strcmp(word,"ylat") == 0) compute_ylat();
else if (strcmp(word,"zlat") == 0) compute_zlat();
else if (strcmp(word,"bonds") == 0) compute_bonds();
else if (strcmp(word,"angles") == 0) compute_angles();
else if (strcmp(word,"dihedrals") == 0) compute_dihedrals();
else if (strcmp(word,"impropers") == 0) compute_impropers();
else if (strcmp(word,"pxx") == 0) {
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
if (update->whichflag == 0) {
if (pressure->invoked_vector != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pressure->invoked_flag & INVOKED_VECTOR)) {
pressure->compute_vector();
pressure->invoked_flag |= INVOKED_VECTOR;
}
compute_pxx();
} else if (strcmp(word,"pyy") == 0) {
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
if (update->whichflag == 0) {
if (pressure->invoked_vector != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pressure->invoked_flag & INVOKED_VECTOR)) {
pressure->compute_vector();
pressure->invoked_flag |= INVOKED_VECTOR;
}
compute_pyy();
} else if (strcmp(word,"pzz") == 0) {
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
if (update->whichflag == 0) {
if (pressure->invoked_vector != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pressure->invoked_flag & INVOKED_VECTOR)) {
pressure->compute_vector();
pressure->invoked_flag |= INVOKED_VECTOR;
}
compute_pzz();
} else if (strcmp(word,"pxy") == 0) {
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
if (update->whichflag == 0) {
if (pressure->invoked_vector != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pressure->invoked_flag & INVOKED_VECTOR)) {
pressure->compute_vector();
pressure->invoked_flag |= INVOKED_VECTOR;
}
compute_pxy();
} else if (strcmp(word,"pxz") == 0) {
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
if (update->whichflag == 0) {
if (pressure->invoked_vector != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pressure->invoked_flag & INVOKED_VECTOR)) {
pressure->compute_vector();
pressure->invoked_flag |= INVOKED_VECTOR;
}
compute_pxz();
} else if (strcmp(word,"pyz") == 0) {
if (!pressure)
error->all(FLERR,"Thermo keyword in variable requires "
"thermo to use/init press");
if (update->whichflag == 0) {
if (pressure->invoked_vector != update->ntimestep)
error->all(FLERR,"Compute used in variable thermo keyword between runs "
"is not current");
} else if (!(pressure->invoked_flag & INVOKED_VECTOR)) {
pressure->compute_vector();
pressure->invoked_flag |= INVOKED_VECTOR;
}
compute_pyz();
}
else if (strcmp(word,"fmax") == 0) compute_fmax();
else if (strcmp(word,"fnorm") == 0) compute_fnorm();
else if (strcmp(word,"nbuild") == 0) {
compute_nbuild();
dvalue = bivalue;
} else if (strcmp(word,"ndanger") == 0) {
compute_ndanger();
dvalue = bivalue;
}
else if (strcmp(word,"cella") == 0) compute_cella();
else if (strcmp(word,"cellb") == 0) compute_cellb();
else if (strcmp(word,"cellc") == 0) compute_cellc();
else if (strcmp(word,"cellalpha") == 0) compute_cellalpha();
else if (strcmp(word,"cellbeta") == 0) compute_cellbeta();
else if (strcmp(word,"cellgamma") == 0) compute_cellgamma();
else return 1;
*answer = dvalue;
return 0;
}
/* ----------------------------------------------------------------------
extraction of Compute, Fix, Variable results
compute/fix are normalized by atoms if returning extensive value
variable value is not normalized (formula should normalize if desired)
------------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
void Thermo::compute_compute()
{
int m = field2index[ifield];
Compute *compute = computes[m];
if (compute_which[m] == SCALAR) {
dvalue = compute->scalar;
if (normflag && compute->extscalar) dvalue /= natoms;
} else if (compute_which[m] == VECTOR) {
dvalue = compute->vector[argindex1[ifield]-1];
if (normflag) {
if (compute->extvector == 0) return;
else if (compute->extvector == 1) dvalue /= natoms;
else if (compute->extlist[argindex1[ifield]-1]) dvalue /= natoms;
}
} else {
dvalue = compute->array[argindex1[ifield]-1][argindex2[ifield]-1];
if (normflag && compute->extarray) dvalue /= natoms;
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_fix()
{
int m = field2index[ifield];
Fix *fix = fixes[m];
if (argindex1[ifield] == 0) {
dvalue = fix->compute_scalar();
if (normflag && fix->extscalar) dvalue /= natoms;
} else if (argindex2[ifield] == 0) {
dvalue = fix->compute_vector(argindex1[ifield]-1);
if (normflag) {
if (fix->extvector == 0) return;
else if (fix->extvector == 1) dvalue /= natoms;
else if (fix->extlist[argindex1[ifield]-1]) dvalue /= natoms;
}
} else {
dvalue = fix->compute_array(argindex1[ifield]-1,argindex2[ifield]-1);
if (normflag && fix->extarray) dvalue /= natoms;
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_variable()
{
dvalue = input->variable->compute_equal(variables[field2index[ifield]]);
}
/* ----------------------------------------------------------------------
one method for every keyword thermo can output
called by compute() or evaluate_keyword()
compute will have already been called
set ivalue/dvalue/bivalue if value is int/double/bigint
customize a new keyword by adding a method
------------------------------------------------------------------------- */
void Thermo::compute_step()
{
bivalue = update->ntimestep;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_elapsed()
{
bivalue = update->ntimestep - update->firststep;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_elapsed_long()
{
bivalue = update->ntimestep - update->beginstep;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_dt()
{
dvalue = update->dt;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_time()
{
dvalue = update->atime + (update->ntimestep-update->atimestep)*update->dt;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cpu()
{
if (firststep == 0) dvalue = 0.0;
else dvalue = timer->elapsed(Timer::TOTAL);
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_tpcpu()
{
double new_cpu;
double new_time = update->ntimestep * update->dt;
if (firststep == 0) {
new_cpu = 0.0;
dvalue = 0.0;
} else {
new_cpu = timer->elapsed(Timer::TOTAL);
double cpu_diff = new_cpu - last_tpcpu;
double time_diff = new_time - last_time;
if (time_diff > 0.0 && cpu_diff > 0.0) dvalue = time_diff/cpu_diff;
else dvalue = 0.0;
}
last_time = new_time;
last_tpcpu = new_cpu;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_spcpu()
{
double new_cpu;
int new_step = update->ntimestep;
if (firststep == 0) {
new_cpu = 0.0;
dvalue = 0.0;
} else {
new_cpu = timer->elapsed(Timer::TOTAL);
double cpu_diff = new_cpu - last_spcpu;
int step_diff = new_step - last_step;
if (cpu_diff > 0.0) dvalue = step_diff/cpu_diff;
else dvalue = 0.0;
}
last_step = new_step;
last_spcpu = new_cpu;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cpuremain()
{
if (firststep == 0) dvalue = 0.0;
else dvalue = timer->elapsed(Timer::TOTAL) *
(update->laststep - update->ntimestep) /
(update->ntimestep - update->firststep);
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_part()
{
ivalue = universe->iworld;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_atoms()
{
bivalue = atom->natoms;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_temp()
{
dvalue = temperature->scalar;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_press()
{
dvalue = pressure->scalar;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pe()
{
dvalue = pe->scalar;
if (normflag) dvalue /= natoms;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_ke()
{
dvalue = temperature->scalar;
dvalue *= 0.5 * temperature->dof * force->boltz;
if (normflag) dvalue /= natoms;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_etotal()
{
compute_pe();
double ke = temperature->scalar;
ke *= 0.5 * temperature->dof * force->boltz;
if (normflag) ke /= natoms;
dvalue += ke;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_enthalpy()
{
compute_etotal();
double etmp = dvalue;
compute_vol();
double vtmp = dvalue;
if (normflag) vtmp /= natoms;
compute_press();
double ptmp = dvalue;
dvalue = etmp + ptmp*vtmp/(force->nktv2p);
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_evdwl()
{
double tmp = 0.0;
if (force->pair) tmp += force->pair->eng_vdwl;
MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
if (force->pair && force->pair->tail_flag) {
double volume = domain->xprd * domain->yprd * domain->zprd;
dvalue += force->pair->etail / volume;
}
if (normflag) dvalue /= natoms;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_ecoul()
{
double tmp = 0.0;
if (force->pair) tmp += force->pair->eng_coul;
MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
if (normflag) dvalue /= natoms;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_epair()
{
double tmp = 0.0;
if (force->pair) tmp += force->pair->eng_vdwl + force->pair->eng_coul;
MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
if (force->kspace) dvalue += force->kspace->energy;
if (force->pair && force->pair->tail_flag) {
double volume = domain->xprd * domain->yprd * domain->zprd;
dvalue += force->pair->etail / volume;
}
if (normflag) dvalue /= natoms;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_ebond()
{
if (force->bond) {
double tmp = force->bond->energy;
MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
if (normflag) dvalue /= natoms;
} else dvalue = 0.0;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_eangle()
{
if (force->angle) {
double tmp = force->angle->energy;
MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
if (normflag) dvalue /= natoms;
} else dvalue = 0.0;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_edihed()
{
if (force->dihedral) {
double tmp = force->dihedral->energy;
MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
if (normflag) dvalue /= natoms;
} else dvalue = 0.0;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_eimp()
{
if (force->improper) {
double tmp = force->improper->energy;
MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
if (normflag) dvalue /= natoms;
} else dvalue = 0.0;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_emol()
{
double tmp = 0.0;
if (atom->molecular) {
if (force->bond) tmp += force->bond->energy;
if (force->angle) tmp += force->angle->energy;
if (force->dihedral) tmp += force->dihedral->energy;
if (force->improper) tmp += force->improper->energy;
MPI_Allreduce(&tmp,&dvalue,1,MPI_DOUBLE,MPI_SUM,world);
if (normflag) dvalue /= natoms;
} else dvalue = 0.0;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_elong()
{
if (force->kspace) {
dvalue = force->kspace->energy;
if (normflag) dvalue /= natoms;
} else dvalue = 0.0;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_etail()
{
if (force->pair && force->pair->tail_flag) {
double volume = domain->xprd * domain->yprd * domain->zprd;
dvalue = force->pair->etail / volume;
if (normflag) dvalue /= natoms;
} else dvalue = 0.0;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_vol()
{
if (domain->dimension == 3)
dvalue = domain->xprd * domain->yprd * domain->zprd;
else
dvalue = domain->xprd * domain->yprd;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_density()
{
double mass = group->mass(0);
compute_vol();
dvalue = force->mv2d * mass/dvalue;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_lx()
{
dvalue = domain->xprd;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_ly()
{
dvalue = domain->yprd;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_lz()
{
dvalue = domain->zprd;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_xlo()
{
dvalue = domain->boxlo[0];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_xhi()
{
dvalue = domain->boxhi[0];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_ylo()
{
dvalue = domain->boxlo[1];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_yhi()
{
dvalue = domain->boxhi[1];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_zlo()
{
dvalue = domain->boxlo[2];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_zhi()
{
dvalue = domain->boxhi[2];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_xy()
{
dvalue = domain->xy;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_xz()
{
dvalue = domain->xz;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_yz()
{
dvalue = domain->yz;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_xlat()
{
dvalue = domain->lattice->xlattice;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_ylat()
{
dvalue = domain->lattice->ylattice;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_zlat()
{
dvalue = domain->lattice->zlattice;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_bonds()
{
bivalue = atom->nbonds;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_angles()
{
bivalue = atom->nangles;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_dihedrals()
{
bivalue = atom->ndihedrals;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_impropers()
{
bivalue = atom->nimpropers;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pxx()
{
dvalue = pressure->vector[0];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pyy()
{
dvalue = pressure->vector[1];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pzz()
{
dvalue = pressure->vector[2];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pxy()
{
dvalue = pressure->vector[3];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pxz()
{
dvalue = pressure->vector[4];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_pyz()
{
dvalue = pressure->vector[5];
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_fmax()
{
double **f = atom->f;
int nlocal = atom->nlocal;
double max = 0.0;
for (int i = 0; i < nlocal; i++) {
max = MAX(max,fabs(f[i][0]));
max = MAX(max,fabs(f[i][1]));
max = MAX(max,fabs(f[i][2]));
}
double maxall;
MPI_Allreduce(&max,&maxall,1,MPI_DOUBLE,MPI_MAX,world);
dvalue = maxall;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_fnorm()
{
double **f = atom->f;
int nlocal = atom->nlocal;
double dot = 0.0;
for (int i = 0; i < nlocal; i++)
dot += f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
double dotall;
MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
dvalue = sqrt(dotall);
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_nbuild()
{
bivalue = neighbor->ncalls;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_ndanger()
{
bivalue = neighbor->ndanger;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cella()
{
dvalue = domain->xprd;
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellb()
{
if (!domain->triclinic)
dvalue = domain->yprd;
else {
double* h = domain->h;
dvalue = sqrt(h[1]*h[1]+h[5]*h[5]);
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellc()
{
if (!domain->triclinic)
dvalue = domain->zprd;
else {
double* h = domain->h;
dvalue = sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]);
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellalpha()
{
if (!domain->triclinic)
dvalue = 90.0;
else {
// Cos(alpha) = (xy.xz + ly.yz)/(b.c)
double* h = domain->h;
double cosalpha = (h[5]*h[4]+h[1]*h[3])/
sqrt((h[1]*h[1]+h[5]*h[5])*(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]));
dvalue = acos(cosalpha)*180.0/MY_PI;
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellbeta()
{
if (!domain->triclinic)
dvalue = 90.0;
else {
// Cos(beta) = xz/c
double* h = domain->h;
double cosbeta = h[4]/sqrt(h[2]*h[2]+h[3]*h[3]+h[4]*h[4]);
dvalue = acos(cosbeta)*180.0/MY_PI;
}
}
/* ---------------------------------------------------------------------- */
void Thermo::compute_cellgamma()
{
if (!domain->triclinic)
dvalue = 90.0;
else {
// Cos(gamma) = xy/b
double* h = domain->h;
double cosgamma = h[5]/sqrt(h[1]*h[1]+h[5]*h[5]);
dvalue = acos(cosgamma)*180.0/MY_PI;
}
}