git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6835 f3b2605a-c512-4ea7-a41b-209d697bcdaa

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
 

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 
 

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;
+}

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;
 };
 
 }

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;
+  }
 }
 
 /* ----------------------------------------------------------------------

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

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;
+  }
 }
 
 /* ----------------------------------------------------------------------

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 torqueflag,erforceflag;
-  int erforceflag;
+  int e_flag,rho_flag;
 
   int *newton;                      // newton flag at each level
   class FixRespa *fix_respa;        // Fix to store the force level array

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,6 +399,9 @@ 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
 

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;
+      }
     }
   }
 }

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 torqueflag,erforceflag;
-  int erforceflag;
+  int e_flag,rho_flag;
 
   void force_clear();
 };