diff --git a/src/atom.cpp b/src/atom.cpp index b9c9c3e707..a0f65365e6 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -83,6 +83,10 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) eradius = ervel = erforce = NULL; cs = csforce = vforce = ervelforce = NULL; etag = NULL; + rho = drho = NULL; + e = de = NULL; + cv = NULL; + vest = NULL; maxspecial = 1; nspecial = NULL; @@ -106,12 +110,13 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) // customize by adding new flag sphere_flag = ellipsoid_flag = peri_flag = electron_flag = 0; - wavepacket_flag = 0; + wavepacket_flag = sph_flag = 0; molecule_flag = q_flag = mu_flag = 0; rmass_flag = radius_flag = omega_flag = torque_flag = angmom_flag = 0; vfrac_flag = spin_flag = eradius_flag = ervel_flag = erforce_flag = 0; cs_flag = csforce_flag = vforce_flag = ervelforce_flag= etag_flag = 0; + rho_flag = e_flag = cv_flag = vest_flag = 0; // ntype-length arrays diff --git a/src/atom.h b/src/atom.h index 9b56510ca3..98b687ec9a 100644 --- a/src/atom.h +++ b/src/atom.h @@ -56,6 +56,10 @@ class Atom : protected Pointers { double *eradius,*ervel,*erforce,*ervelforce; double *cs,*csforce,*vforce; int *etag; + double *rho, *drho; + double *e, *de; + double **vest; + double *cv; int **nspecial; // 0,1,2 = cummulative # of 1-2,1-3,1-4 neighs int **special; // IDs of 1-2,1-3,1-4 neighs of each atom @@ -81,12 +85,13 @@ class Atom : protected Pointers { // customize by adding new flag int sphere_flag,ellipsoid_flag,peri_flag,electron_flag; - int wavepacket_flag; + int wavepacket_flag,sph_flag; int molecule_flag,q_flag,mu_flag; int rmass_flag,radius_flag,omega_flag,torque_flag,angmom_flag; int vfrac_flag,spin_flag,eradius_flag,ervel_flag,erforce_flag; int cs_flag,csforce_flag,vforce_flag,ervelforce_flag,etag_flag; + int rho_flag,e_flag,cv_flag,vest_flag; // extra peratom info in restart file destined for fix & diag diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp index 7c25edf769..e3d2ae8ead 100644 --- a/src/fix_gravity.cpp +++ b/src/fix_gravity.cpp @@ -34,6 +34,9 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) : if (narg < 5) error->all("Illegal fix gravity command"); time_depend = 1; + scalar_flag = 1; + global_freq = 1; + extscalar = 1; magnitude = atof(arg[3]); @@ -90,6 +93,9 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) : } time_origin = update->ntimestep; + + eflag = 0; + egrav = 0.0; } /* ---------------------------------------------------------------------- */ @@ -98,6 +104,7 @@ int FixGravity::setmask() { int mask = 0; mask |= POST_FORCE; + mask |= THERMO_ENERGY; mask |= POST_FORCE_RESPA; return mask; } @@ -155,6 +162,7 @@ void FixGravity::post_force(int vflag) zacc = magnitude*zgrav; } + double **x = atom->x; double **f = atom->f; double *rmass = atom->rmass; double *mass = atom->mass; @@ -163,6 +171,9 @@ void FixGravity::post_force(int vflag) int nlocal = atom->nlocal; double massone; + eflag = 0; + egrav = 0.0; + if (rmass) { for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { @@ -170,6 +181,7 @@ void FixGravity::post_force(int vflag) f[i][0] += massone*xacc; f[i][1] += massone*yacc; f[i][2] += massone*zacc; + egrav -= massone * (xacc*x[i][0] + yacc*x[i][1] + zacc*x[i][2]); } } else { for (int i = 0; i < nlocal; i++) @@ -178,6 +190,7 @@ void FixGravity::post_force(int vflag) f[i][0] += massone*xacc; f[i][1] += massone*yacc; f[i][2] += massone*zacc; + egrav -= massone * (xacc*x[i][0] + yacc*x[i][1] + zacc*x[i][2]); } } } @@ -188,3 +201,18 @@ void FixGravity::post_force_respa(int vflag, int ilevel, int iloop) { if (ilevel == nlevels_respa-1) post_force(vflag); } + +/* ---------------------------------------------------------------------- + potential energy in gravity field +------------------------------------------------------------------------- */ + +double FixGravity::compute_scalar() +{ + // only sum across procs one time + + if (eflag == 0) { + MPI_Allreduce(&egrav,&egrav_all,1,MPI_DOUBLE,MPI_SUM,world); + eflag = 1; + } + return egrav_all; +} diff --git a/src/fix_gravity.h b/src/fix_gravity.h index 248a9ffa18..3f82d6c7f8 100644 --- a/src/fix_gravity.h +++ b/src/fix_gravity.h @@ -34,6 +34,7 @@ class FixGravity : public Fix { void setup(int); void post_force(int); void post_force_respa(int, int, int); + double compute_scalar(); private: int style; @@ -44,6 +45,8 @@ class FixGravity : public Fix { double degree2rad; int nlevels_respa; int time_origin; + int eflag; + double egrav,egrav_all; }; } diff --git a/src/min.cpp b/src/min.cpp index 7ae69fbb64..55253f11e8 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -137,12 +137,16 @@ void Min::init() ev_setup(); // set flags for what arrays to clear in force_clear() - // need to clear torques,erforce if arrays exists + // need to clear additionals arrays if they exist torqueflag = 0; if (atom->torque_flag) torqueflag = 1; erforceflag = 0; if (atom->erforce_flag) erforceflag = 1; + e_flag = 0; + if (atom->e_flag) e_flag = 1; + rho_flag = 0; + if (atom->rho_flag) rho_flag = 1; // orthogonal vs triclinic simulation box @@ -537,6 +541,16 @@ void Min::force_clear() double *erforce = atom->erforce; for (i = 0; i < nall; i++) erforce[i] = 0.0; } + + if (e_flag) { + double *de = atom->de; + for (i = 0; i < nall; i++) de[i] = 0.0; + } + + if (rho_flag) { + double *drho = atom->drho; + for (i = 0; i < nall; i++) drho[i] = 0.0; + } } /* ---------------------------------------------------------------------- diff --git a/src/min.h b/src/min.h index db9590b1b2..f367fa7b93 100644 --- a/src/min.h +++ b/src/min.h @@ -60,9 +60,10 @@ class Min : protected Pointers { class Compute **vlist_global; class Compute **vlist_atom; + int triclinic; // 0 if domain is orthog, 1 if triclinic int pairflag; int torqueflag,erforceflag; - int triclinic; // 0 if domain is orthog, 1 if triclinic + int e_flag,rho_flag; int narray; // # of arrays stored by fix_minimize class FixMinimize *fix_minimize; // fix that stores auxiliary data diff --git a/src/respa.cpp b/src/respa.cpp index 9786b2875e..83457da738 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -278,12 +278,16 @@ void Respa::init() ev_setup(); // set flags for what arrays to clear in force_clear() - // need to clear torques,erforce if arrays exists + // need to clear additionals arrays if they exist torqueflag = 0; if (atom->torque_flag) torqueflag = 1; erforceflag = 0; if (atom->erforce_flag) erforceflag = 1; + e_flag = 0; + if (atom->e_flag) e_flag = 1; + rho_flag = 0; + if (atom->rho_flag) rho_flag = 1; // step[] = timestep for each level @@ -623,6 +627,16 @@ void Respa::force_clear(int newtonflag) double *erforce = atom->erforce; for (i = 0; i < nall; i++) erforce[i] = 0.0; } + + if (e_flag) { + double *de = atom->de; + for (i = 0; i < nall; i++) de[i] = 0.0; + } + + if (rho_flag) { + double *drho = atom->drho; + for (i = 0; i < nall; i++) drho[i] = 0.0; + } } /* ---------------------------------------------------------------------- diff --git a/src/respa.h b/src/respa.h index e0981c764e..409184c739 100644 --- a/src/respa.h +++ b/src/respa.h @@ -53,8 +53,8 @@ class Respa : public Integrate { private: int triclinic; // 0 if domain is orthog, 1 if triclinic - int torqueflag; // zero out arrays every step - int erforceflag; + int torqueflag,erforceflag; + int e_flag,rho_flag; int *newton; // newton flag at each level class FixRespa *fix_respa; // Fix to store the force level array diff --git a/src/set.cpp b/src/set.cpp index 941a57db75..a541fcca98 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -35,7 +35,8 @@ using namespace LAMMPS_NS; enum{ATOM_SELECT,MOL_SELECT,TYPE_SELECT,GROUP_SELECT,REGION_SELECT}; enum{TYPE,TYPE_FRACTION,MOLECULE,X,Y,Z,CHARGE,MASS,SHAPE, DIPOLE,DIPOLE_RANDOM,QUAT,QUAT_RANDOM, - DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER}; + DIAMETER,DENSITY,VOLUME,IMAGE,BOND,ANGLE,DIHEDRAL,IMPROPER, + MESO_E,MESO_CV,MESO_RHO}; #define BIG INT_MAX @@ -271,6 +272,27 @@ void Set::command(int narg, char **arg) error->all("Invalid value in set command"); topology(IMPROPER); iarg += 2; + } else if (strcmp(arg[iarg],"meso_e") == 0) { + if (iarg+2 > narg) error->all("Illegal set command"); + dvalue = atof(arg[iarg+1]); + if (!atom->e_flag) + error->all("Cannot set this attribute for this atom style"); + set(MESO_E); + iarg += 2; + } else if (strcmp(arg[iarg],"meso_cv") == 0) { + if (iarg+2 > narg) error->all("Illegal set command"); + dvalue = atof(arg[iarg+1]); + if (!atom->cv_flag) + error->all("Cannot set this attribute for this atom style"); + set(MESO_CV); + iarg += 2; + } else if (strcmp(arg[iarg],"meso_rho") == 0) { + if (iarg+2 > narg) error->all("Illegal set command"); + dvalue = atof(arg[iarg+1]); + if (!atom->rho_flag) + error->all("Cannot set meso_rho for this atom style"); + set(MESO_RHO); + iarg += 2; } else error->all("Illegal set command"); // statistics @@ -377,7 +399,10 @@ void Set::set(int keyword) else if (keyword == MASS) atom->rmass[i] = dvalue; else if (keyword == DIAMETER) atom->radius[i] = 0.5 * dvalue; else if (keyword == VOLUME) atom->vfrac[i] = dvalue; - + else if (keyword == MESO_E) atom->e[i] = dvalue; + else if (keyword == MESO_CV) atom->cv[i] = dvalue; + else if (keyword == MESO_RHO) atom->rho[i] = dvalue; + // set shape else if (keyword == SHAPE) diff --git a/src/verlet.cpp b/src/verlet.cpp index 10932c10fb..4fc4024f42 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -63,12 +63,16 @@ void Verlet::init() ev_setup(); // set flags for what arrays to clear in force_clear() - // need to clear torques,erforce if arrays exists + // need to clear additionals arrays if they exist torqueflag = 0; if (atom->torque_flag) torqueflag = 1; erforceflag = 0; if (atom->erforce_flag) erforceflag = 1; + e_flag = 0; + if (atom->e_flag) e_flag = 1; + rho_flag = 0; + if (atom->rho_flag) rho_flag = 1; // orthogonal vs triclinic simulation box @@ -334,6 +338,16 @@ void Verlet::force_clear() for (i = 0; i < nall; i++) erforce[i] = 0.0; } + if (e_flag) { + double *de = atom->de; + for (i = 0; i < nall; i++) de[i] = 0.0; + } + + if (rho_flag) { + double *drho = atom->drho; + for (i = 0; i < nall; i++) drho[i] = 0.0; + } + // neighbor includegroup flag is set // clear force only on initial nfirst particles // if either newton flag is set, also include ghosts @@ -362,6 +376,16 @@ void Verlet::force_clear() for (i = 0; i < nall; i++) erforce[i] = 0.0; } + if (e_flag) { + double *de = atom->de; + for (i = 0; i < nall; i++) de[i] = 0.0; + } + + if (rho_flag) { + double *drho = atom->drho; + for (i = 0; i < nall; i++) drho[i] = 0.0; + } + if (force->newton) { nall = atom->nlocal + atom->nghost; @@ -384,6 +408,16 @@ void Verlet::force_clear() double *erforce = atom->erforce; for (i = atom->nlocal; i < nall; i++) erforce[i] = 0.0; } + + if (e_flag) { + double *de = atom->de; + for (i = 0; i < nall; i++) de[i] = 0.0; + } + + if (rho_flag) { + double *drho = atom->drho; + for (i = 0; i < nall; i++) drho[i] = 0.0; + } } } } diff --git a/src/verlet.h b/src/verlet.h index 8b8c90b7a9..327f72b833 100644 --- a/src/verlet.h +++ b/src/verlet.h @@ -36,8 +36,8 @@ class Verlet : public Integrate { protected: int triclinic; // 0 if domain is orthog, 1 if triclinic - int torqueflag; // zero out arrays every step - int erforceflag; + int torqueflag,erforceflag; + int e_flag,rho_flag; void force_clear(); };