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

src/compute.h

@@ -30,9 +30,9 @@ class Compute : protected Pointers {
 
   int scalar_flag;          // 0/1 if compute_scalar() function exists
   int vector_flag;          // 0/1 if compute_vector() function exists
-  int peratom_flag;         // 0/1 if compute_peratom() function exists
   int size_vector;          // N = size of global vector
-  int size_peratom;         // 0 = just scalar_atom, N = size of vector_atom
+  int peratom_flag;         // 0/1 if compute_peratom() function exists
+  int size_peratom;         // 0 = scalar_atom, N = size of vector_atom
 
   int extensive;      // 0/1 if scalar,vector are intensive/extensive values
   int tempflag;       // 1 if Compute can be used as temperature

src/fix.cpp

@@ -42,7 +42,8 @@ Fix::Fix(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp)
   rigid_flag = 0;
   virial_flag = 0;
   no_change_box = 0;
-  peratom_flag = 0;
+
+  scalar_flag = vector_flag = peratom_flag = 0;
 
   comm_forward = comm_reverse = 0;
 

src/fix.h

@@ -35,6 +35,11 @@ class Fix : protected Pointers {
   int virial_flag;               // 1 if Fix contributes to virial, 0 if not
   int no_change_box;             // 1 if cannot swap ortho <-> triclinic
 
+  int scalar_flag;               // 0/1 if compute_scalar() function exists
+  int vector_flag;               // 0/1 if compute_vector() function exists
+  int size_vector;               // N = size of global vector
+  int scalar_vector_freq;        // frequency compute s/v data is available at
+
   int peratom_flag;              // 0/1 if per-atom data is stored
   int size_peratom;              // 0 = scalar_atom, N = size of vector_atom
   double *scalar_atom;           // computed per-atom scalar
@@ -93,7 +98,8 @@ class Fix : protected Pointers {
   virtual int pack_reverse_comm(int, int, double *) {return 0;}
   virtual void unpack_reverse_comm(int, int *, double *) {}
 
-  virtual double thermo(int) {return 0.0;}
+  virtual double compute_scalar() {return 0.0;}
+  virtual double compute_vector(int) {return 0.0;}
 
   virtual int dof(int) {return 0;}
   virtual void deform(int) {}

src/fix_ave_spatial.cpp

@@ -260,10 +260,9 @@ void FixAveSpatial::init()
     if (ifix < 0) 
       error->all("Fix ID for fix ave/spatial does not exist");
     fix = modify->fix[ifix];
-    if (nevery % modify->fix[ifix]->peratom_freq)
+    if (nevery % fix->peratom_freq)
       error->all("Fix ave/spatial and fix not computed at compatible times");
   }
-
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_ave_time.cpp

@@ -26,29 +26,35 @@
 
 using namespace LAMMPS_NS;
 
+enum{COMPUTE,FIX};
+
 /* ---------------------------------------------------------------------- */
 
 FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
-  if (narg != 9) error->all("Illegal fix ave/time command");
+  if (narg != 10) error->all("Illegal fix ave/time command");
 
   nevery = atoi(arg[3]);
   nrepeat = atoi(arg[4]);
   nfreq = atoi(arg[5]);
 
-  int n = strlen(arg[6]) + 1;
+  if (strcmp(arg[6],"compute") == 0) which = COMPUTE;
-  id_compute = new char[n];
+  else if (strcmp(arg[6],"fix") == 0) which = FIX;
-  strcpy(id_compute,arg[6]);
+  else error->all("Illegal fix ave/time command");
 
-  int flag = atoi(arg[7]);
+  int n = strlen(arg[7]) + 1;
+  id = new char[n];
+  strcpy(id,arg[7]);
+
+  int flag = atoi(arg[8]);
 
   MPI_Comm_rank(world,&me);
   if (me == 0) {
-    fp = fopen(arg[8],"w");
+    fp = fopen(arg[9],"w");
     if (fp == NULL) {
       char str[128];
-      sprintf(str,"Cannot open fix ave/time file %s",arg[8]);
+      sprintf(str,"Cannot open fix ave/time file %s",arg[9]);
       error->one(str);
     }
   }
@@ -59,31 +65,52 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
   if (nfreq < nevery || nfreq % nevery || (nrepeat-1)*nevery >= nfreq)
     error->all("Illegal fix ave/time command");
 
-  int icompute = modify->find_compute(id_compute);
+  int icompute,ifix;
+  if (which == COMPUTE) {
+    icompute = modify->find_compute(id);
     if (icompute < 0) error->all("Compute ID for fix ave/time does not exist");
+  } else {
+    ifix = modify->find_fix(id);
+    if (ifix < 0) error->all("Fix ID for fix ave/time does not exist");
+  }
 
   if (flag < 0 || flag > 2) error->all("Illegal fix ave/time command");
   sflag = vflag = 0;
   if (flag == 0 || flag == 2) sflag = 1;
   if (flag == 1 || flag == 2) vflag = 1;
 
+  if (which == COMPUTE) {
     if (sflag && modify->compute[icompute]->scalar_flag == 0)
       error->all("Fix ave/time compute does not calculate a scalar");
     if (vflag && modify->compute[icompute]->vector_flag == 0)
       error->all("Fix ave/time compute does not calculate a vector");
+  } else {
+    if (sflag && modify->fix[ifix]->scalar_flag == 0)
+      error->all("Fix ave/time fix does not calculate a scalar");
+    if (vflag && modify->fix[ifix]->vector_flag == 0)
+      error->all("Fix ave/time fix does not calculate a vector");
+  }
 
-  if (modify->compute[icompute]->pressflag) pressure_every = nevery;
+  if (which == COMPUTE &&
+      modify->compute[icompute]->pressflag) pressure_every = nevery;
 
   // setup list of computes to call, including pre-computes
 
+  compute = NULL;
+  if (which == COMPUTE) {
     ncompute = 1 + modify->compute[icompute]->npre;
     compute = new Compute*[ncompute];
+  } else ncompute = 0;
 
   // print header into file
 
   if (me == 0) {
+    if (which == COMPUTE)
       fprintf(fp,"Time-averaged data for fix %s, group %s, and compute %s\n",
-	    id,group->names[modify->compute[icompute]->igroup],id_compute);
+	      id,group->names[modify->compute[icompute]->igroup],id);
+    else
+      fprintf(fp,"Time-averaged data for fix %s, group %s, and fix %s\n",
+	      id,group->names[modify->fix[ifix]->igroup],id);
     if (sflag and !vflag)
       fprintf(fp,"TimeStep Value\n");
     else if (!sflag and vflag)
@@ -94,7 +121,8 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
 
   vector = NULL;
   if (vflag) {
-    size_vector = modify->compute[icompute]->size_vector;
+    if (which == COMPUTE) size_vector = modify->compute[icompute]->size_vector;
+    else size_vector = modify->fix[ifix]->size_vector;
     vector = new double[size_vector];
   }
 
@@ -111,7 +139,7 @@ FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
 
 FixAveTime::~FixAveTime()
 {
-  delete [] id_compute;
+  delete [] id;
   if (me == 0) fclose(fp);
   delete [] compute;
   delete [] vector;
@@ -132,7 +160,8 @@ void FixAveTime::init()
 {
   // set ptrs to one or more computes called each end-of-step
 
-  int icompute = modify->find_compute(id_compute);
+  if (which == COMPUTE) {
+    int icompute = modify->find_compute(id);
     if (icompute < 0)
       error->all("Compute ID for fix ave/time does not exist");
   
@@ -146,6 +175,19 @@ void FixAveTime::init()
       }
     
     compute[ncompute++] = modify->compute[icompute];
+  }
+
+  // set ptr to fix ID
+  // check that fix frequency is acceptable
+
+  if (which == FIX) {
+    int ifix = modify->find_fix(id);
+    if (ifix < 0) 
+      error->all("Fix ID for fix ave/time does not exist");
+    fix = modify->fix[ifix];
+    if (nevery % fix->scalar_vector_freq)
+      error->all("Fix ave/time and fix not computed at compatible times");
+  }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -166,8 +208,9 @@ void FixAveTime::end_of_step()
       for (i = 0; i < size_vector; i++) vector[i] = 0.0;
   }
 
-  // accumulate results of compute to local copy
+  // accumulate results of compute or fix to local copy
   
+  if (which == COMPUTE) {
     if (sflag) {
       double value;
       for (i = 0; i < ncompute; i++) value = compute[i]->compute_scalar();
@@ -179,6 +222,13 @@ void FixAveTime::end_of_step()
       for (i = 0; i < size_vector; i++) vector[i] += cvector[i];
     }
 
+  } else {
+    if (sflag) scalar += fix->compute_scalar();
+    if (vflag)
+      for (i = 0; i < size_vector; i++)
+	vector[i] += fix->compute_vector(i);
+  }
+
   irepeat++;
   nvalid += nevery;
 

src/fix_ave_time.h

@@ -29,8 +29,8 @@ class FixAveTime : public Fix {
 
  private:
   int me;
-  int nrepeat,nfreq,nvalid,irepeat;
+  int nrepeat,nfreq,nvalid,irepeat,which,ifix;
-  char *id_compute;
+  char *id;
   FILE *fp;
 
   int sflag,vflag;
@@ -38,6 +38,7 @@ class FixAveTime : public Fix {
   double scalar,*vector;
   int ncompute;
   class Compute **compute;
+  class Fix *fix;
 };
 
 }

src/fix_indent.cpp

@@ -37,6 +37,12 @@ FixIndent::FixIndent(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
   if (narg < 4) error->all("Illegal fix indent command");
+
+  scalar_flag = 1;
+  vector_flag = 1;
+  size_vector = 3;
+  scalar_vector_freq = 1;
+
   k = atof(arg[3]);
 
   // set defaults
@@ -122,7 +128,6 @@ void FixIndent::init()
 
 void FixIndent::setup()
 {
-  eflag_enable = 1;
   if (strcmp(update->integrate_style,"verlet") == 0)
     post_force(1);
   else {
@@ -130,14 +135,12 @@ void FixIndent::setup()
     post_force_respa(1,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
-  eflag_enable = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixIndent::min_setup()
 {
-  eflag_enable = 1;
   post_force(1);
 }
 
@@ -145,10 +148,6 @@ void FixIndent::min_setup()
 
 void FixIndent::post_force(int vflag)
 {
-  bool eflag = false;
-  if (eflag_enable) eflag = true;
-  else if (output->next_thermo == update->ntimestep) eflag = true;
-
   // set current r0
   // for minimization, always set to r0_stop
 
@@ -160,8 +159,10 @@ void FixIndent::post_force(int vflag)
     r0 = r0_start + delta * (r0_stop-r0_start);
   }
 
-  double eng;
+  // indenter values, 0 = energy, 1-3 = force components
-  if (eflag) eng = 0.0;
+
+  indenter_flag = 0;
+  indenter[0] = indenter[1] = indenter[2] = indenter[3] = 0.0;
 
   // spherical indenter
 
@@ -179,7 +180,7 @@ void FixIndent::post_force(int vflag)
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
 
-    double delx,dely,delz,r,dr,fmag;
+    double delx,dely,delz,r,dr,fmag,fx,fy,fz;
 
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
@@ -190,10 +191,16 @@ void FixIndent::post_force(int vflag)
 	dr = r - r0;
 	if (dr >= 0.0) continue;
 	fmag = k*dr*dr;
-	f[i][0] += delx*fmag/r;
+	fx = delx*fmag/r;
-	f[i][1] += dely*fmag/r;
+	fy = dely*fmag/r;
-	f[i][2] += delz*fmag/r;
+	fz = delz*fmag/r;
-	if (eflag) eng -= k3 * dr*dr*dr;
+	f[i][0] += fx;
+	f[i][1] += fy;
+	f[i][2] += fz;
+	indenter[0] -= k3 * dr*dr*dr;
+	indenter[1] -= fx;
+	indenter[2] -= fy;
+	indenter[3] -= fz;
       }
 
   // cylindrical indenter
@@ -220,7 +227,7 @@ void FixIndent::post_force(int vflag)
     int *mask = atom->mask;
     int nlocal = atom->nlocal;
     
-    double delx,dely,delz,r,dr,fmag;
+    double delx,dely,delz,r,dr,fmag,fx,fy,fz;
     
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
@@ -241,14 +248,18 @@ void FixIndent::post_force(int vflag)
 	dr = r - r0;
 	if (dr >= 0.0) continue;
 	fmag = k*dr*dr;
-	f[i][0] += delx*fmag/r;
+	fx = delx*fmag/r;
-	f[i][1] += dely*fmag/r;
+	fy = dely*fmag/r;
-	f[i][2] += delz*fmag/r;
+	fz = delz*fmag/r;
-	if (eflag) eng -= k3 * dr*dr*dr;
+	f[i][0] += fx;
+	f[i][1] += fy;
+	f[i][2] += fz;
+	indenter[0] -= k3 * dr*dr*dr;
+	indenter[1] -= fx;
+	indenter[2] -= fy;
+	indenter[3] -= fz;
       }
   }
-
-  if (eflag) MPI_Allreduce(&eng,&etotal,1,MPI_DOUBLE,MPI_SUM,world);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -265,12 +276,34 @@ void FixIndent::min_post_force(int vflag)
   post_force(vflag);
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   energy of indenter interaction
+------------------------------------------------------------------------- */
 
-double FixIndent::thermo(int n)
+double FixIndent::compute_scalar()
 {
-  if (n == 0) return etotal;
+  // only sum across procs one time
-  else return 0.0;
+
+  if (indenter_flag == 0) {
+    MPI_Allreduce(indenter,indenter_all,4,MPI_DOUBLE,MPI_SUM,world);
+    indenter_flag = 1;
+  }
+  return indenter_all[0];
+}
+
+/* ----------------------------------------------------------------------
+   components of force on indenter
+------------------------------------------------------------------------- */
+
+double FixIndent::compute_vector(int n)
+{
+  // only sum across procs one time
+
+  if (indenter_flag == 0) {
+    MPI_Allreduce(indenter,indenter_all,4,MPI_DOUBLE,MPI_SUM,world);
+    indenter_flag = 1;
+  }
+  return indenter_all[n];
 }
 
 /* ----------------------------------------------------------------------

src/fix_indent.h

@@ -28,12 +28,15 @@ class FixIndent : public Fix {
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);
-  double thermo(int);
+  double compute_scalar();
+  double compute_vector(int);
 
  private:
   int istyle,scaleflag,radflag,thermo_flag,eflag_enable;
-  double k,k3,eng,etotal;
+  double k,k3;
   double x0,y0,z0,r0_stop,r0_start;
+  int indenter_flag;
+  double indenter[4],indenter_all[4];
   int cdim;
   double c1,c2;
   double vx,vy,vz;

src/fix_nph.cpp

@@ -46,6 +46,8 @@ FixNPH::FixNPH(LAMMPS *lmp, int narg, char **arg) :
   restart_global = 1;
   pressure_every = 1;
   box_change = 1;
+  scalar_flag = 1;
+  scalar_vector_freq = 1;
 
   double p_period[3];
   if (strcmp(arg[3],"xyz") == 0) {
@@ -306,7 +308,7 @@ void FixNPH::init()
 
 void FixNPH::setup()
 {
-  p_target[0] = p_start[0];                 // needed by thermo() method
+  p_target[0] = p_start[0];                 // used by compute_scalar()
   p_target[1] = p_start[1];
   p_target[2] = p_start[2];
 
@@ -767,7 +769,7 @@ int FixNPH::modify_param(int narg, char **arg)
 
 /* ---------------------------------------------------------------------- */
 
-double FixNPH::thermo(int n)
+double FixNPH::compute_scalar()
 {
   double volume;
   if (dimension == 3) volume = domain->xprd * domain->yprd * domain->zprd;
@@ -781,6 +783,5 @@ double FixNPH::thermo(int n)
       energy += 0.5*nkt*omega_dot[i]*omega_dot[i] / 
 	(p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p);
 
-  if (n == 0) return energy;
+  return energy;
-  else return 0.0;
 }

src/fix_nph.h

@@ -29,7 +29,7 @@ class FixNPH : public Fix {
   void final_integrate();
   void initial_integrate_respa(int,int);
   void final_integrate_respa(int);
-  double thermo(int);
+  double compute_scalar();
   void write_restart(FILE *);
   void restart(char *);
   int modify_param(int, char **);

src/fix_npt.cpp

@@ -45,6 +45,8 @@ FixNPT::FixNPT(LAMMPS *lmp, int narg, char **arg) :
   restart_global = 1;
   pressure_every = 1;
   box_change = 1;
+  scalar_flag = 1;
+  scalar_vector_freq = 1;
 
   t_start = atof(arg[3]);
   t_stop = atof(arg[4]);
@@ -307,7 +309,7 @@ void FixNPT::init()
 
 void FixNPT::setup()
 {
-  t_target = t_start;                      // used by thermo()
+  t_target = t_start;                      // used by compute_scalar()
   p_target[0] = p_start[0];
   p_target[1] = p_start[1];
   p_target[2] = p_start[2];
@@ -796,7 +798,7 @@ int FixNPT::modify_param(int narg, char **arg)
 
 /* ---------------------------------------------------------------------- */
 
-double FixNPT::thermo(int n)
+double FixNPT::compute_scalar()
 {
   double ke = temperature->dof * boltz * t_target;
   double keplus = atom->natoms * boltz * t_target;
@@ -812,6 +814,5 @@ double FixNPT::thermo(int n)
       energy += 0.5*keplus*omega_dot[i]*omega_dot[i] / 
 	(p_freq[i]*p_freq[i]) + p_target[i]*(volume-vol0) / (pdim*nktv2p);
 
-  if (n == 0) return energy;
+  return energy;
-  else return 0.0;
 }

src/fix_npt.h

@@ -29,7 +29,7 @@ class FixNPT : public Fix {
   virtual void final_integrate();
   void initial_integrate_respa(int, int);
   void final_integrate_respa(int);
-  double thermo(int);
+  double compute_scalar();
   void write_restart(FILE *);
   void restart(char *);
   int modify_param(int, char **);

src/fix_nvt.cpp

@@ -39,6 +39,8 @@ FixNVT::FixNVT(LAMMPS *lmp, int narg, char **arg) :
   if (narg < 6) error->all("Illegal fix nvt command");
 
   restart_global = 1;
+  scalar_flag = 1;
+  scalar_vector_freq = 1;
 
   t_start = atof(arg[3]);
   t_stop = atof(arg[4]);
@@ -129,7 +131,7 @@ void FixNVT::init()
 
 void FixNVT::setup()
 {
-  t_target = t_start;                         // needed by thermo() method
+  t_target = t_start;                         // used by compute_scalar()
   t_current = temperature->compute_scalar();
 }
 
@@ -377,11 +379,9 @@ void FixNVT::reset_target(double t_new)
 
 /* ---------------------------------------------------------------------- */
 
-double FixNVT::thermo(int n)
+double FixNVT::compute_scalar()
 {
   double ke = temperature->dof * force->boltz * t_target;
   double energy = ke * (eta + 0.5*eta_dot*eta_dot/(t_freq*t_freq));
-
+  return energy;
-  if (n == 0) return energy;
-  else return 0.0;
 }

src/fix_nvt.h

@@ -29,7 +29,7 @@ class FixNVT : public Fix {
   virtual void final_integrate();
   virtual void initial_integrate_respa(int,int);
   void final_integrate_respa(int);
-  double thermo(int);
+  double compute_scalar();
   void write_restart(FILE *);
   void restart(char *);
   int modify_param(int, char **);

src/fix_orient_fcc.cpp

@@ -46,6 +46,9 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
 
   if (narg != 11) error->all("Illegal fix orient/fcc command");
 
+  scalar_flag = 1;
+  scalar_vector_freq = 1;
+
   nstats = atoi(arg[3]);
   direction_of_motion = atoi(arg[4]);
   a = atof(arg[5]);
@@ -194,7 +197,6 @@ void FixOrientFCC::init_list(int id, NeighList *ptr)
 
 void FixOrientFCC::setup()
 {
-  eflag_enable = 1;
   if (strcmp(update->integrate_style,"verlet") == 0)
     post_force(1);
   else {
@@ -202,7 +204,6 @@ void FixOrientFCC::setup()
     post_force_respa(1,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
-  eflag_enable = 0;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -211,16 +212,12 @@ void FixOrientFCC::post_force(int vflag)
 {
   int i,j,k,ii,jj,inum,jnum,m,n,nn,nsort,id_self;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  double edelta,added_energy,omega;
+  double edelta,omega;
   double dx,dy,dz,rsq,xismooth,xi_sq,duxi,duxi_other;
   double dxi[3];
   double *dxiptr;
   bool found_myself;
 
-  bool eflag = false;
-  if (eflag_enable) eflag = true;
-  else if (output->next_thermo == update->ntimestep) eflag = true;
-
   // set local ptrs
 
   double **x = atom->x;
@@ -246,7 +243,7 @@ void FixOrientFCC::post_force(int vflag)
   // loop over owned atoms and build Nbr data structure of neighbors
   // use full neighbor list
 
-  if (eflag) added_energy = 0.0;
+  added_energy = 0.0;
   int count = 0;
   int mincount = BIG;
   int maxcount = 0;
@@ -319,16 +316,16 @@ void FixOrientFCC::post_force(int vflag)
 
     if (xi_total < xi0) {
       nbr[i].duxi = 0.0;
-      if (eflag) edelta = 0.0;
+      edelta = 0.0;
     } else if (xi_total > xi1) {
       nbr[i].duxi = 0.0;
-      if (eflag) edelta = Vxi;
+      edelta = Vxi;
     } else {
       omega = (0.5*PI)*(xi_total-xi0) / (xi1-xi0);
       nbr[i].duxi = PI*Vxi*sin(2.0*omega) / (2.0*(xi1-xi0));
-      if (eflag) edelta = Vxi*(1 - cos(2.0*omega)) / 2.0;
+      edelta = Vxi*(1 - cos(2.0*omega)) / 2.0;
     }
-    if (eflag) added_energy += edelta;
+    added_energy += edelta;
   }
   
   if (maxcount) delete [] sort;
@@ -391,11 +388,6 @@ void FixOrientFCC::post_force(int vflag)
     }
   }
 
-  // sum energy across procs
-
-  if (eflag)
-    MPI_Allreduce(&added_energy,&total_added_e,1,MPI_DOUBLE,MPI_SUM,world);
-
   // print statistics every nstats timesteps
 
   if (nstats && update->ntimestep % nstats == 0) {
@@ -433,10 +425,10 @@ void FixOrientFCC::post_force_respa(int vflag, int ilevel, int iloop)
 
 /* ---------------------------------------------------------------------- */
 
-double FixOrientFCC::thermo(int n)
+double FixOrientFCC::compute_scalar()
 {
-  if (n == 0) return total_added_e;
+  MPI_Allreduce(&added_energy,&total_added_e,1,MPI_DOUBLE,MPI_SUM,world);
-  else return 0.0;
+  return total_added_e;
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_orient_fcc.h

@@ -45,7 +45,7 @@ class FixOrientFCC : public Fix {
   void setup();
   void post_force(int);
   void post_force_respa(int, int, int);
-  double thermo(int);
+  double compute_scalar();
   int pack_comm(int, int *, double *, int, int *);
   void unpack_comm(int, int, double *);
   double memory_usage();
@@ -64,9 +64,8 @@ class FixOrientFCC : public Fix {
   char *xifilename, *chifilename;  // file names for 2 crystal orientations
 
   bool use_xismooth;
-  int eflag_enable;
   double Rxi[12][3],Rchi[12][3],half_xi_chi_vec[2][6][3];
-  double xiid,xi0,xi1,xicutoffsq,cutsq,total_added_e;
+  double xiid,xi0,xi1,xicutoffsq,cutsq,added_energy,total_added_e;
   int half_fcc_nn,nmax;
 
   Nbr *nbr;

src/fix_set_force.cpp

@@ -28,6 +28,10 @@ FixSetForce::FixSetForce(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 6) error->all("Illegal fix setforce command");
 
+  vector_flag = 1;
+  size_vector = 3;
+  scalar_vector_freq = 1;
+
   flagx = flagy = flagz = 1;
   if (strcmp(arg[3],"NULL") == 0) flagx = 0;
   else xvalue = atof(arg[3]);
@@ -86,6 +90,7 @@ void FixSetForce::post_force(int vflag)
   int nlocal = atom->nlocal;
 
   foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
+  force_flag = 0;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
@@ -111,6 +116,7 @@ void FixSetForce::post_force_respa(int vflag, int ilevel, int iloop)
     int nlocal = atom->nlocal;
 
     foriginal[0] = foriginal[1] = foriginal[2] = 0.0;
+    force_flag = 0;
     
     for (int i = 0; i < nlocal; i++)
       if (mask[i] & groupbit) {
@@ -132,15 +138,16 @@ void FixSetForce::min_post_force(int vflag)
 }
 
 /* ----------------------------------------------------------------------
-   allow setforce values to be printed with thermo output
+   return components of total force on fix group before reset
-   n = 1,2,3 = components of total force on fix group before reset
 ------------------------------------------------------------------------- */
 
-double FixSetForce::thermo(int n)
+double FixSetForce::compute_vector(int n)
 {
-  if (n >= 1 && n <= 3) {
+  // only sum across procs one time
-    double ftotal;
+
-    MPI_Allreduce(&foriginal[n-1],&ftotal,1,MPI_DOUBLE,MPI_SUM,world);
+  if (force_flag == 0) {
-    return ftotal;
+    MPI_Allreduce(foriginal,foriginal_all,3,MPI_DOUBLE,MPI_SUM,world);
-  } else return 0.0;
+    force_flag = 1;
+  }
+  return foriginal_all[n];
 }

src/fix_set_force.h

@@ -28,12 +28,13 @@ class FixSetForce : public Fix {
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);
-  double thermo(int);
+  double compute_vector(int);
 
  private:
   int flagx,flagy,flagz;
   double xvalue,yvalue,zvalue;
-  double foriginal[3];
+  double foriginal[3],foriginal_all[3];
+  int force_flag;
   int nlevels_respa;
 };
 

src/fix_temp_rescale.cpp

@@ -36,9 +36,13 @@ FixTempRescale::FixTempRescale(LAMMPS *lmp, int narg, char **arg) :
   Fix(lmp, narg, arg)
 {
   if (narg < 8) error->all("Illegal fix temp/rescale command");
+
   nevery = atoi(arg[3]);
   if (nevery <= 0) error->all("Illegal fix temp/rescale command");
 
+  scalar_flag = 1;
+  scalar_vector_freq = nevery;
+
   t_start = atof(arg[4]);
   t_end = atof(arg[5]);
   t_window = atof(arg[6]);
@@ -223,8 +227,7 @@ int FixTempRescale::modify_param(int narg, char **arg)
 
 /* ---------------------------------------------------------------------- */
 
-double FixTempRescale::thermo(int n)
+double FixTempRescale::compute_scalar()
 {
-  if (n == 0) return energy;
+  return energy;
-  else return 0.0;
 }

src/fix_temp_rescale.h

@@ -25,7 +25,7 @@ class FixTempRescale : public Fix {
   int setmask();
   void init();
   void end_of_step();
-  double thermo(int);
+  double compute_scalar();
   int modify_param(int, char **);
 
  private:

src/fix_wall_lj126.cpp

@@ -35,6 +35,11 @@ FixWallLJ126::FixWallLJ126(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 8) error->all("Illegal fix wall/lj126 command");
 
+  scalar_flag = 1;
+  vector_flag = 1;
+  size_vector = 3;
+  scalar_vector_freq = 1;
+
   if (strcmp(arg[3],"xlo") == 0) {
     dim = 0;
     side = -1;
@@ -101,7 +106,6 @@ void FixWallLJ126::init()
 
 void FixWallLJ126::setup()
 {
-  eflag_enable = 1;
   if (strcmp(update->integrate_style,"verlet") == 0)
     post_force(1);
   else {
@@ -109,14 +113,12 @@ void FixWallLJ126::setup()
     post_force_respa(1,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
-  eflag_enable = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallLJ126::min_setup()
 {
-  eflag_enable = 1;
   post_force(1);
 }
 
@@ -124,17 +126,14 @@ void FixWallLJ126::min_setup()
 
 void FixWallLJ126::post_force(int vflag)
 {
-  bool eflag = false;
-  if (eflag_enable) eflag = true;
-  else if (output->next_thermo == update->ntimestep) eflag = true;
-
   double **x = atom->x;
   double **f = atom->f;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  double delta,rinv,r2inv,r6inv,eng;
+  double delta,rinv,r2inv,r6inv,fwall;
-  if (eflag) eng = 0.0;
+  wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
+  wall_flag = 0;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
@@ -145,11 +144,11 @@ void FixWallLJ126::post_force(int vflag)
       rinv = 1.0/delta;
       r2inv = rinv*rinv;
       r6inv = r2inv*r2inv*r2inv;
-      f[i][dim] -= r6inv*(coeff1*r6inv - coeff2) * side;
+      fwall = r6inv*(coeff1*r6inv - coeff2) * side;
-      if (eflag) eng += r6inv*(coeff3*r6inv - coeff4) - offset;
+      f[i][dim] -= fwall;
+      wall[0] += r6inv*(coeff3*r6inv - coeff4) - offset;
+      wall[dim] += fwall;
     }
-
-  if (eflag) MPI_Allreduce(&eng,&etotal,1,MPI_DOUBLE,MPI_SUM,world);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -166,10 +165,32 @@ void FixWallLJ126::min_post_force(int vflag)
   post_force(vflag);
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   energy of wall interaction
+------------------------------------------------------------------------- */
 
-double FixWallLJ126::thermo(int n)
+double FixWallLJ126::compute_scalar()
 {
-  if (n == 0) return etotal;
+  // only sum across procs one time
-  else return 0.0;
+
+  if (wall_flag == 0) {
+    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
+    wall_flag = 1;
+  }
+  return wall_all[0];
+}
+
+/* ----------------------------------------------------------------------
+   components of force on wall
+------------------------------------------------------------------------- */
+
+double FixWallLJ126::compute_vector(int n)
+{
+  // only sum across procs one time
+
+  if (wall_flag == 0) {
+    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
+    wall_flag = 1;
+  }
+  return wall_all[n];
 }

src/fix_wall_lj126.h

@@ -28,13 +28,15 @@ class FixWallLJ126 : public Fix {
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);
-  double thermo(int);
+  double compute_scalar();
+  double compute_vector(int);
 
  private:
   int dim,side,thermo_flag,eflag_enable;
   double coord,epsilon,sigma,cutoff;
   double coeff1,coeff2,coeff3,coeff4,offset;
-  double etotal;
+  double wall[4],wall_all[4];
+  int wall_flag;
   int nlevels_respa;
 };
 

src/fix_wall_lj93.cpp

@@ -31,6 +31,11 @@ FixWallLJ93::FixWallLJ93(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 8) error->all("Illegal fix wall/lj93 command");
 
+  scalar_flag = 1;
+  vector_flag = 1;
+  size_vector = 3;
+  scalar_vector_freq = 1;
+
   if (strcmp(arg[3],"xlo") == 0) {
     dim = 0;
     side = -1;
@@ -98,7 +103,6 @@ void FixWallLJ93::init()
 
 void FixWallLJ93::setup()
 {
-  eflag_enable = 1;
   if (strcmp(update->integrate_style,"verlet") == 0)
     post_force(1);
   else {
@@ -106,14 +110,12 @@ void FixWallLJ93::setup()
     post_force_respa(1,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
   }
-  eflag_enable = 0;
 }
 
 /* ---------------------------------------------------------------------- */
 
 void FixWallLJ93::min_setup()
 {
-  eflag_enable = 1;
   post_force(1);
 }
 
@@ -121,17 +123,14 @@ void FixWallLJ93::min_setup()
 
 void FixWallLJ93::post_force(int vflag)
 {
-  bool eflag = false;
-  if (eflag_enable) eflag = true;
-  else if (output->next_thermo == update->ntimestep) eflag = true;
-
   double **x = atom->x;
   double **f = atom->f;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-  double delta,rinv,r2inv,r4inv,r10inv,eng;
+  double delta,rinv,r2inv,r4inv,r10inv,fwall;
-  if (eflag) eng = 0.0;
+  wall[0] = wall[1] = wall[2] = wall[3] = 0.0;
+  wall_flag = 0;
 
   for (int i = 0; i < nlocal; i++)
     if (mask[i] & groupbit) {
@@ -143,11 +142,11 @@ void FixWallLJ93::post_force(int vflag)
       r2inv = rinv*rinv;
       r4inv = r2inv*r2inv;
       r10inv = r4inv*r4inv*r2inv;
-      f[i][dim] -= (coeff1*r10inv - coeff2*r4inv) * side;
+      fwall = (coeff1*r10inv - coeff2*r4inv) * side;
-      if (eflag) eng += coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset;
+      f[i][dim] -= fwall;
+      wall[0] += coeff3*r4inv*r4inv*rinv - coeff4*r2inv*rinv - offset;
+      wall[dim] += fwall;
     }
-
-  if (eflag) MPI_Allreduce(&eng,&etotal,1,MPI_DOUBLE,MPI_SUM,world);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -164,10 +163,34 @@ void FixWallLJ93::min_post_force(int vflag)
   post_force(vflag);
 }
 
-/* ---------------------------------------------------------------------- */
+/* ----------------------------------------------------------------------
+   energy of wall interaction
+------------------------------------------------------------------------- */
 
-double FixWallLJ93::thermo(int n)
+double FixWallLJ93::compute_scalar()
 {
-  if (n == 0) return etotal;
+  // only sum across procs one time
-  else return 0.0;
+
+  if (wall_flag == 0) {
+    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
+    wall_flag = 1;
+  }
+  return wall_all[0];
 }
+
+/* ----------------------------------------------------------------------
+   components of force on wall
+------------------------------------------------------------------------- */
+
+double FixWallLJ93::compute_vector(int n)
+{
+  // only sum across procs one time
+
+  if (wall_flag == 0) {
+    MPI_Allreduce(wall,wall_all,4,MPI_DOUBLE,MPI_SUM,world);
+    wall_flag = 1;
+  }
+  return wall_all[n];
+}
+
+

src/fix_wall_lj93.h

@@ -29,13 +29,15 @@ class FixWallLJ93 : public Fix {
   void post_force(int);
   void post_force_respa(int, int, int);
   void min_post_force(int);
-  double thermo(int);
+  double compute_scalar();
+  double compute_vector(int);
 
  private:
   int dim,side,thermo_flag,eflag_enable;
   double coord,epsilon,sigma,cutoff;
   double coeff1,coeff2,coeff3,coeff4,offset;
-  double etotal;
+  double wall[4],wall_all[4];
+  int wall_flag;
   int nlevels_respa;
 };
 

src/modify.cpp

@@ -235,14 +235,14 @@ void Modify::end_of_step()
 /* ----------------------------------------------------------------------
    thermo energy call only for relevant fixes
    called by Thermo class
-   arg to Fix thermo() is 0, so fix will return its energy contribution
+   compute_scalar() is fix call to return energy
 ------------------------------------------------------------------------- */
 
 double Modify::thermo_energy()
 {
   double energy = 0.0;
   for (int i = 0; i < n_thermo_energy; i++)
-    energy += fix[list_thermo_energy[i]]->thermo(0);
+    energy += fix[list_thermo_energy[i]]->compute_scalar();
   return energy;
 }
 

src/thermo.cpp

@@ -793,12 +793,20 @@ void Thermo::parse_fields(char *str)
 	  for (int ic = 0; ic < modify->compute[n]->npre; ic++)
 	    int tmp = add_compute(modify->compute[n]->id_pre[ic],
 				  arg_object[nfield]);
+
 	field2object[nfield] = add_compute(id,arg_object[nfield]);
 	addfield(copy,&Thermo::compute_compute,FLOAT);
 
       } else if (word[0] == 'f') {
 	n = modify->find_fix(id);
 	if (n < 0) error->all("Could not find thermo custom fix ID");
+	if (arg_object[nfield] == 0 && modify->fix[n]->scalar_flag == 0)
+	  error->all("Thermo fix ID does not compute scalar info");
+	if (arg_object[nfield] > 0 && modify->fix[n]->vector_flag == 0)
+	  error->all("Thermo fix ID does not compute vector info");
+	if (arg_object[nfield] > 0 && 
+	    arg_object[nfield] > modify->fix[n]->size_vector)
+	  error->all("Thermo fix ID vector is not large enough");
 	field2object[nfield] = add_fix(id);
 	addfield(copy,&Thermo::compute_fix,FLOAT);
 
@@ -1027,7 +1035,8 @@ void Thermo::compute_compute()
 void Thermo::compute_fix()
 {
   int index = field2object[ifield];
-  dvalue = fixes[index]->thermo(arg_object[ifield]);
+  if (arg_object[ifield] == 0) dvalue = fixes[index]->compute_scalar();
+  else dvalue = fixes[index]->compute_vector(arg_object[ifield]-1);
   if (normflag) dvalue /= natoms;
 }
 

src/variable.cpp

@@ -24,6 +24,7 @@
 #include "domain.h"
 #include "modify.h"
 #include "compute.h"
+#include "fix.h"
 #include "output.h"
 #include "thermo.h"
 #include "memory.h"
@@ -407,6 +408,7 @@ void Variable::copy(int narg, char **from, char **to)
      group function = mass(group), xcm(group,x), ...
      atom vector = x[123], y[3], vx[34], ...
      compute vector = c_mytemp[0], c_thermo_press[3], ...
+     fix vector = f_indent[0], f_setforce[3], ...
    numbers start with a digit or "." or "-" (no parens or brackets)
    keywords start with a lowercase letter (no parens or brackets)
    functions contain ()
@@ -414,7 +416,7 @@ void Variable::copy(int narg, char **from, char **to)
    vectors contain []
      single arg must be integer
      for atom vectors, arg is global ID of atom
-     for compute vectors, 0 is the scalar value, 1-N are vector values
+     for compute/fix vectors, 0 is the scalar value, 1-N are vector values
    see lists of valid functions & vectors below
    return answer = value of string
 ------------------------------------------------------------------------- */
@@ -662,6 +664,8 @@ double Variable::evaluate(char *str, Tree *tree)
   // index = everything between [] evaluated as integer
   // if vector name starts with "c_", trailing chars are compute ID
   //   check if compute ID exists, invoke it with index as arg
+  // if vector name starts with "cf_", trailing chars are fix ID
+  //   check if fix ID exists, call its scalar or vector fn with index as arg
   // else is atom vector
   //   find which proc owns atom via atom->map()
   //   grab atom-based value with index as global atom ID
@@ -696,7 +700,7 @@ double Variable::evaluate(char *str, Tree *tree)
       delete [] id;
       modify->compute[icompute]->init();
 
-      // call compute() if index = 0, else compute_vector()
+      // call compute_scalar() if index = 0, else compute_vector()
       // make pre-call to Compute object's pre-compute(s) if defined
 
       int index = atoi(arg);
@@ -713,7 +717,7 @@ double Variable::evaluate(char *str, Tree *tree)
 	answer = modify->compute[icompute]->compute_scalar();
       } else if (index > 0) {
 	if (modify->compute[icompute]->vector_flag == 0)
-	  error->all("Variable compute ID does not compute scalar info");
+	  error->all("Variable compute ID does not compute vector info");
 	if (index > modify->compute[icompute]->size_vector)
 	  error->all("Variable compute ID vector is not large enough");
 	if (modify->compute[icompute]->npre)
@@ -727,6 +731,33 @@ double Variable::evaluate(char *str, Tree *tree)
 	answer = modify->compute[icompute]->vector[index-1];
       } else error->all("Invalid compute ID index in variable");
 
+    } else if (strncmp(vector,"f_",2) == 0) {
+      n = strlen(vector) - 2 + 1;
+      char *id = new char[n];
+      strcpy(id,&vector[2]);
+
+      int ifix;
+      for (ifix = 0; ifix < modify->nfix; ifix++)
+	if (strcmp(id,modify->fix[ifix]->id) == 0) break;
+      if (ifix == modify->nfix)
+	error->all("Invalid fix ID in variable");
+      delete [] id;
+
+      // call compute_scalar() if index = 0, else compute_vector()
+
+      int index = atoi(arg);
+      if (index == 0) {
+	if (modify->fix[ifix]->scalar_flag == 0)
+	  error->all("Variable fix ID does not compute scalar info");
+	answer = modify->fix[ifix]->compute_scalar();
+      } else if (index > 0) {
+	if (modify->fix[ifix]->vector_flag == 0)
+	  error->all("Variable fix ID does not compute vector info");
+	if (index > modify->fix[ifix]->size_vector)
+	  error->all("Variable fix ID vector is not large enough");
+	answer = modify->fix[ifix]->compute_vector(index-1);
+      } else error->all("Invalid fix ID index in variable");
+
     } else if (tree) {
 
       if (strlen(arg)) error->all("Invalid atom vector in variable");