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

src/compute_variable_atom.cpp

@@ -0,0 +1,96 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include "string.h"
+#include "compute_variable_atom.h"
+#include "atom.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+ComputeVariableAtom::ComputeVariableAtom(LAMMPS *lmp, int narg, char **arg) :
+  Compute(lmp, narg, arg)
+{
+  if (narg != 4) error->all("Illegal compute variable/atom command");
+
+  // store variable name
+
+  int n = strlen(arg[3]) + 1;
+  varname = new char[n];
+  strcpy(varname,arg[3]);
+
+  peratom_flag = 1;
+  size_peratom = 0;
+
+  nmax = 0;
+  result = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+ComputeVariableAtom::~ComputeVariableAtom()
+{
+  delete [] varname;
+  memory->sfree(result);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeVariableAtom::init()
+{
+  // set ivariable used by this compute
+
+  ivariable = input->variable->find(varname);
+  if (ivariable < 0)
+    error->all("Could not find compute variable name");
+
+  // test if variable of correct ATOM type
+}
+
+/* ---------------------------------------------------------------------- */
+
+void ComputeVariableAtom::compute_peratom()
+{
+  // grow result array if necessary
+
+  if (atom->nlocal > nmax) {
+    memory->sfree(result);
+    nmax = atom->nmax;
+    result = (double *) memory->smalloc(nmax*sizeof(double),
+					"compute/variable/atom:result");
+    scalar_atom = result;
+  }
+
+  // parse variable once to create parse tree
+  // evaluate tree for all atoms, will be zero for atoms not in group
+  // free parse tree memory stored by Variable
+
+  input->variable->build_parse_tree(ivariable);
+  input->variable->evaluate_parse_tree(igroup,result);
+  input->variable->free_parse_tree();
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based array
+------------------------------------------------------------------------- */
+
+int ComputeVariableAtom::memory_usage()
+{
+  int bytes = nmax * sizeof(double);
+  return bytes;
+}

src/compute_variable_atom.h

@@ -0,0 +1,37 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef COMPUTE_VARIABLE_ATOM_H
+#define COMPUTE_VARIABLE_ATOM_H
+
+#include "compute.h"
+
+namespace LAMMPS_NS {
+
+class ComputeVariableAtom : public Compute {
+ public:
+  ComputeVariableAtom(class LAMMPS *, int, char **);
+  ~ComputeVariableAtom();
+  void init();
+  void compute_peratom();
+  int memory_usage();
+
+ private:
+  int nmax,ivariable;
+  char *varname;
+  double *result;
+};
+
+}
+
+#endif

src/fix_ave_spatial.cpp

@@ -0,0 +1,439 @@
+/* ----------------------------------------------------------------------
+   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: Pieter in't Veld (SNL)
+------------------------------------------------------------------------- */
+
+#include "stdlib.h"
+#include "string.h"
+#include "fix_ave_spatial.h"
+#include "atom.h"
+#include "update.h"
+#include "domain.h"
+#include "lattice.h"
+#include "modify.h"
+#include "compute.h"
+#include "group.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+enum{LOWER,CENTER,UPPER,COORD};
+enum{DENSITY_MASS,DENSITY_NUM,VX,VY,VZ,FX,FY,FZ,COMPUTE};
+enum{SAMPLE,ALL};
+
+/* ---------------------------------------------------------------------- */
+
+FixAveSpatial::FixAveSpatial(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg < 11) error->all("Illegal fix ave/spatial command");
+
+  nevery = atoi(arg[3]);
+  nfreq = atoi(arg[4]);
+
+  if (strcmp(arg[5],"x") == 0) dim = 0;
+  else if (strcmp(arg[5],"y") == 0) dim = 1;
+  else if (strcmp(arg[5],"z") == 0) dim = 2;
+  else error->all("Illegal fix ave/spatial command");
+
+  if (strcmp(arg[6],"lower") == 0) originflag = LOWER;
+  if (strcmp(arg[6],"center") == 0) originflag = CENTER;
+  if (strcmp(arg[6],"upper") == 0) originflag = UPPER;
+  else originflag = COORD;
+  if (originflag == COORD) origin = atof(arg[6]);
+
+  delta = atof(arg[7]);
+
+  MPI_Comm_rank(world,&me);
+  if (me == 0) {
+    fp = fopen(arg[8],"w");
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open fix ave/spatial file %s",arg[8]);
+      error->one(str);
+    }
+  }
+
+  if (strcmp(arg[9],"density") == 0) {
+    if (strcmp(arg[10],"mass") == 0) which = DENSITY_MASS;
+    else if (strcmp(arg[10],"number") == 0) which = DENSITY_NUM;
+    else error->all("Illegal fix ave/spatial command");
+
+  } else if (strcmp(arg[9],"atom") == 0) {
+    if (strcmp(arg[10],"vx") == 0) which = VX;
+    else if (strcmp(arg[10],"vy") == 0) which = VY;
+    else if (strcmp(arg[10],"vz") == 0) which = VZ;
+    else if (strcmp(arg[10],"fx") == 0) which = FX;
+    else if (strcmp(arg[10],"fy") == 0) which = FY;
+    else if (strcmp(arg[10],"fz") == 0) which = FZ;
+    else error->all("Illegal fix ave/spatial command");
+
+  } else if (strcmp(arg[9],"compute") == 0) {
+    which = COMPUTE;
+    int n = strlen(arg[10]) + 1;
+    id_compute = new char[n];
+    strcpy(id_compute,arg[10]);
+
+  } else error->all("Illegal fix ave/spatial command");
+
+  // parse optional args
+
+  normflag = ALL;
+  int scaleflag = 1;
+
+  int iarg = 11;
+  while (iarg < narg) {
+    if (strcmp(arg[iarg],"norm") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
+      if (strcmp(arg[iarg+1],"all") == 0) normflag = ALL;
+      else if (strcmp(arg[iarg+1],"sample") == 0) normflag = SAMPLE;
+      else error->all("Illegal fix ave/spatial command");
+      iarg += 2;
+    } else if (strcmp(arg[iarg],"units") == 0) {
+      if (iarg+2 > narg) error->all("Illegal fix ave/spatial command");
+      if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
+      else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
+      else error->all("Illegal fix ave/spatial command");
+      iarg += 2;
+    } else error->all("Illegal fix ave/spatial command");
+  }
+
+  // if density, no normalization by atom count should be done
+  // thus ALL and SAMPLE should give same answer, but code does normalize
+  // thus only ALL is computed correctly, so force norm to be ALL
+
+  if (which == DENSITY_MASS || which == DENSITY_NUM) normflag = ALL;
+
+  // setup scaling
+
+  if (scaleflag && domain->lattice == NULL)
+    error->all("Use of fix ave/spatial with undefined lattice");
+
+  if (scaleflag) {
+    xscale = domain->lattice->xlattice;
+    yscale = domain->lattice->ylattice;
+    zscale = domain->lattice->zlattice;
+  }
+  else xscale = yscale = zscale = 1.0;
+
+  // apply scaling factors
+
+  double scale;
+  if (dim == 0) scale = xscale;
+  if (dim == 1) scale = yscale;
+  if (dim == 2) scale = zscale;
+  delta *= scale;
+  if (originflag == COORD) origin *= scale;
+
+  // setup and error check
+
+  if (nevery <= 0) error->all("Illegal fix ave/spatial command");
+  if (nfreq < nevery || nfreq % nevery)
+    error->all("Illegal fix ave/spatial command");
+
+  if (delta <= 0.0) error->all("Illegal fix ave/spatial command");
+  invdelta = 1.0/delta;
+
+  nvalues = 1;
+  if (which == COMPUTE) {
+    int icompute = modify->find_compute(id_compute);
+    if (icompute < 0)
+      error->all("Compute ID for fix ave/spatial does not exist");
+    if (modify->compute[icompute]->peratom_flag == 0)
+      error->all("Fix ave/spatial compute does not calculate per-atom info");
+    nvalues = compute_size_peratom = modify->compute[icompute]->size_peratom;
+    if (nvalues == 0) nvalues = 1;
+  }
+
+  if (me == 0) {
+    fprintf(fp,"Spatial-averaged data for fix %s, group %s, and %s %s\n",
+	    id,group->names[igroup],arg[9],arg[10]);
+    fprintf(fp,"TimeStep Number-of-layers (one per snapshot)\n");
+    fprintf(fp,"Layer Coord Atoms Value(s) (one per layer)\n");
+  }
+  
+  nsum = nlayers = maxlayer = 0;
+  coord = NULL;
+  count_one = count_many = count_total = NULL;
+  values_one = values_many = values_total = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixAveSpatial::~FixAveSpatial()
+{
+  if (which == COMPUTE) delete [] id_compute;
+  if (me == 0) fclose(fp);
+
+  memory->sfree(coord);
+  memory->sfree(count_one);
+  memory->sfree(count_many);
+  memory->sfree(count_total);
+  memory->destroy_2d_double_array(values_one);
+  memory->destroy_2d_double_array(values_many);
+  memory->destroy_2d_double_array(values_total);
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixAveSpatial::setmask()
+{
+  int mask = 0;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveSpatial::init()
+{
+  // set ptrs to current compute and precompute
+
+  if (which == COMPUTE) {
+    int icompute = modify->find_compute(id_compute);
+    if (icompute < 0) 
+      error->all("Compute ID for fix ave/spatial does not exist");
+    compute = modify->compute[icompute];
+    
+    if (compute->id_pre) {
+      icompute = modify->find_compute(compute->id_pre);
+      if (icompute < 0)
+	error->all("Precompute ID for fix ave/spatial does not exist");
+      precompute = modify->compute[icompute];
+    } else precompute = NULL;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveSpatial::end_of_step()
+{
+  int i,j,m,ilayer,nstide;
+  double lo,hi;
+
+  // if computing the first sample, setup layers
+  // compute current origin = boundary for some layer
+  // lo = layer boundary immediately below boxlo
+  // hi = layer boundary immediately above boxhi
+  // allocate and initialize arrays based on new layer count
+
+  if (nsum == 0) {
+    double *boxlo = domain->boxlo;
+    double *boxhi = domain->boxhi;
+    if (originflag == LOWER) origin = boxlo[dim];
+    else if (originflag == UPPER) origin = boxhi[dim];
+    else if (originflag == CENTER)
+      origin = 0.5 * (boxlo[dim] + boxhi[dim]);
+
+    if (origin < domain->boxlo[dim]) {
+      m = static_cast<int> ((domain->boxlo[dim] - origin) * invdelta);
+      lo = origin + m*delta;
+    } else {
+      m = static_cast<int> ((origin - domain->boxlo[dim]) * invdelta);
+      lo = origin - m*delta;
+      if (lo > domain->boxlo[dim]) lo -= delta;
+    }
+    if (origin < domain->boxhi[dim]) {
+      m = static_cast<int> ((domain->boxhi[dim] - origin) * invdelta);
+      hi = origin + m*delta;
+      if (hi < boxhi[dim]) hi += delta;
+    } else {
+      m = static_cast<int> ((origin - domain->boxhi[dim]) * invdelta);
+      hi = origin - m*delta;
+    }
+
+    offset = lo;
+    nlayers = static_cast<int> ((hi-lo) * invdelta + 0.5);
+    double volume = domain->xprd * domain->yprd * domain->zprd;
+    layer_volume = delta * volume/domain->prd[dim];
+
+    if (nlayers > maxlayer) {
+      maxlayer = nlayers;
+      coord = (double *) memory->srealloc(coord,nlayers*sizeof(double),
+					  "ave/spatial:coord");
+      count_one = (double *) 
+	memory->srealloc(count_one,nlayers*sizeof(double),
+			 "ave/spatial:count_one");
+      count_many = (double *) 
+	memory->srealloc(count_many,nlayers*sizeof(double),
+			 "ave/spatial:count_many");
+      count_total = (double *) 
+	memory->srealloc(count_total,nlayers*sizeof(double),
+			 "ave/spatial:count_total");
+      values_one = memory->grow_2d_double_array(values_one,nlayers,nvalues,
+						"ave/spatial:values_one");
+      values_many = memory->grow_2d_double_array(values_many,nlayers,nvalues,
+						 "ave/spatial:values_many");
+      values_total = memory->grow_2d_double_array(values_total,nlayers,nvalues,
+						"ave/spatial:values_total");
+    }
+
+    for (m = 0; m < nlayers; m++) {
+      coord[m] = offset + (m+0.5)*delta;
+      count_many[m] = count_total[m] = 0.0;
+      for (i = 0; i < nvalues; i++) values_many[m][i] = 0.0;
+    }
+  }
+  
+  // zero out arrays for one sample
+
+  nsum++;
+  for (m = 0; m < nlayers; m++) {
+    count_one[m] = 0.0;
+    for (i = 0; i < nvalues; i++) values_one[m][i] = 0.0;
+  }
+
+  // perform the computation for one sample
+  // sum within each layer, only include atoms in fix group
+  // DENSITY_MASS adds mass
+  // DENSITY_NUM adds 1 to values
+  // ATOM adds atom vector to values
+  // COMPUTE adds its vector to values
+
+  double **x = atom->x;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+
+  if (which == DENSITY_MASS) {
+    int *type = atom->type;
+    double *mass = atom->mass;
+    double *rmass = atom->rmass;
+
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	ilayer = static_cast<int> ((x[dim][i] - offset) * invdelta);
+	if (ilayer < 0) ilayer = 0;
+	if (ilayer >= nlayers) ilayer = nlayers-1;
+	count_one[ilayer] += 1.0;
+	if (mass) values_one[ilayer][0] += mass[type[i]];
+	else values_one[ilayer][0] += rmass[i];
+      }
+    }
+
+  } else if (which == DENSITY_NUM) {
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	ilayer = static_cast<int> ((x[dim][i] - offset) * invdelta);
+	if (ilayer < 0) ilayer = 0;
+	if (ilayer >= nlayers) ilayer = nlayers-1;
+	count_one[ilayer] += 1.0;
+	values_one[ilayer][0] += 1.0;
+      }
+    }
+
+  } else if (which != COMPUTE) {
+    double *vector;
+    int nstride = 3;
+    if (which == VX) vector = &atom->v[0][0];
+    else if (which == VY) vector = &atom->v[0][1];
+    else if (which == VZ) vector = &atom->v[0][2];
+    else if (which == FX) vector = &atom->f[0][0];
+    else if (which == FY) vector = &atom->f[0][1];
+    else if (which == FZ) vector = &atom->f[0][2];
+
+    m = 0;
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	ilayer = static_cast<int> ((x[dim][i] - offset) * invdelta);
+	if (ilayer < 0) ilayer = 0;
+	if (ilayer >= nlayers) ilayer = nlayers-1;
+	count_one[ilayer] += 1.0;
+	values_one[ilayer][0] += vector[m];
+      }
+      m += nstride;
+    }
+
+  } else {
+    if (precompute) precompute->compute_peratom();
+    compute->compute_peratom();
+    double *scalar = compute->scalar_atom;
+    double **vector = compute->vector_atom;
+
+    m = 0;
+    for (i = 0; i < nlocal; i++) {
+      if (mask[i] & groupbit) {
+	ilayer = static_cast<int> ((x[dim][i] - offset) * invdelta);
+	if (ilayer < 0) ilayer = 0;
+	if (ilayer >= nlayers) ilayer = nlayers-1;
+	count_one[ilayer] += 1.0;
+	if (compute_size_peratom == 0) values_one[ilayer][0] += scalar[i];
+	else
+	  for (j = 0; j < nvalues; j++)
+	    values_one[ilayer][j] += vector[i][j];
+      }
+    }
+  }
+
+  // average a single sample
+
+  if (normflag == ALL) {
+    for (m = 0; m < nlayers; m++) {
+      count_many[m] += count_one[m];
+      for (j = 0; j < nvalues; j++)
+	values_many[m][j] += values_one[m][j];
+    }
+  } else {
+    MPI_Allreduce(count_one,count_many,nlayers,MPI_DOUBLE,MPI_SUM,world);
+    for (m = 0; m < nlayers; m++) {
+      if (count_many[m] > 0.0)
+	for (j = 0; j < nvalues; j++)
+	  values_many[m][j] += values_one[m][j]/count_many[m];
+      count_total[m] += count_many[m];
+    }
+  }
+
+  // output the results
+  // time average across samples
+  // if density, also normalize by volume
+
+  if (update->ntimestep % nfreq == 0) {
+    if (normflag == ALL) {
+      MPI_Allreduce(count_many,count_total,nlayers,MPI_DOUBLE,MPI_SUM,world);
+      MPI_Allreduce(&values_many[0][0],&values_total[0][0],nlayers*nvalues,
+		    MPI_DOUBLE,MPI_SUM,world);
+      for (m = 0; m < nlayers; m++) {
+	if (count_total[m] > 0.0)
+	  for (j = 0; j < nvalues; j++)
+	    values_total[m][j] /= count_total[m];
+	count_total[m] /= nsum;
+      }
+    } else {
+      MPI_Allreduce(&values_many[0][0],&values_total[0][0],nlayers*nvalues,
+		    MPI_DOUBLE,MPI_SUM,world);
+      for (m = 0; m < nlayers; m++) {
+	for (j = 0; j < nvalues; j++)
+	  values_total[m][j] /= nsum;
+	count_total[m] /= nsum;
+      }
+    }
+
+    if (which == DENSITY_MASS || which == DENSITY_NUM) {
+      for (m = 0; m < nlayers; m++)
+	values_total[m][0] *= count_total[m] / layer_volume;
+    }
+
+    if (me == 0) {
+      fprintf(fp,"%d %d\n",update->ntimestep,nlayers);
+      for (m = 0; m < nlayers; m++) {
+	fprintf(fp,"  %d %g %g",m+1,coord[m],count_total[m]);
+	for (i = 0; i < nvalues; i++) fprintf(fp," %g",values_total[m][i]);
+	fprintf(fp,"\n");
+      }
+    }
+
+    nsum = 0;
+ }
+}

src/fix_ave_spatial.h

@@ -0,0 +1,52 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_AVE_SPATIAL_H
+#define FIX_AVE_SPATIAL_H
+
+#include "stdio.h"
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixAveSpatial : public Fix {
+ public:
+  FixAveSpatial(class LAMMPS *, int, char **);
+  ~FixAveSpatial();
+  int setmask();
+  void init();
+  void end_of_step();
+
+ private:
+  int me;
+  int nfreq;
+  int dim,originflag,which,normflag;
+  double origin,delta;
+  char *id_compute;
+  FILE *fp;
+
+  int nlayers,nvalues,nsum,maxlayer;
+  int compute_size_peratom;
+  double xscale,yscale,zscale;
+  double layer_volume;
+  double *coord;
+  double *count_one,*count_many,*count_total;
+  double **values_one,**values_many,**values_total;
+  double offset,invdelta; 
+  class Compute *compute;
+  class Compute *precompute;
+};
+
+}
+
+#endif

src/fix_ave_time.cpp

@@ -0,0 +1,166 @@
+/* ----------------------------------------------------------------------
+   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: Pieter in't Veld (SNL)
+------------------------------------------------------------------------- */
+
+#include "stdlib.h"
+#include "string.h"
+#include "fix_ave_time.h"
+#include "update.h"
+#include "modify.h"
+#include "compute.h"
+#include "group.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+/* ---------------------------------------------------------------------- */
+
+FixAveTime::FixAveTime(LAMMPS *lmp, int narg, char **arg) :
+  Fix(lmp, narg, arg)
+{
+  if (narg != 8) error->all("Illegal fix ave/time command");
+
+  nevery = atoi(arg[3]);
+  nfreq = atoi(arg[4]);
+
+  int n = strlen(arg[5]) + 1;
+  id_compute = new char[n];
+  strcpy(id_compute,arg[5]);
+
+  int flag = atoi(arg[6]);
+
+  MPI_Comm_rank(world,&me);
+  if (me == 0) {
+    fp = fopen(arg[7],"w");
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open fix ave/time file %s",arg[7]);
+      error->one(str);
+    }
+  }
+
+  // setup and error check
+
+  if (nevery <= 0) error->all("Illegal fix ave/time command");
+  if (nfreq < nevery || nfreq % nevery)
+    error->all("Illegal fix ave/time command");
+
+  int icompute = modify->find_compute(id_compute);
+  if (icompute < 0) error->all("Compute 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 (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");
+
+  if (modify->compute[icompute]->pressflag) pressure_every = nevery;
+
+  if (me == 0) {
+    fprintf(fp,"Time-averaged data for fix %s, group %s, and compute %s\n",
+	    id,group->names[modify->compute[icompute]->igroup],id_compute);
+    if (sflag and !vflag)
+      fprintf(fp,"TimeStep Value\n");
+    else if (!sflag and vflag)
+      fprintf(fp,"TimeStep Vector-values\n");
+    else if (!sflag and vflag)
+      fprintf(fp,"TimeStep Scalar-value Vector-values\n");
+  }
+
+  nsum = 0;
+  scalar = 0.0;
+  vector = NULL;
+  if (vflag) {
+    size_vector = modify->compute[icompute]->size_vector;
+    vector = new double[size_vector];
+    for (int i = 0; i < size_vector; i++) vector[i] = 0.0;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixAveTime::~FixAveTime()
+{
+  delete [] id_compute;
+  if (me == 0) fclose(fp);
+  delete [] vector;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixAveTime::setmask()
+{
+  int mask = 0;
+  mask |= END_OF_STEP;
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveTime::init()
+{
+  // set ptrs to current compute and precompute
+
+  int icompute = modify->find_compute(id_compute);
+  if (icompute < 0) error->all("Compute ID for fix ave/time does not exist");
+  compute = modify->compute[icompute];
+
+  if (compute->id_pre) {
+    icompute = modify->find_compute(compute->id_pre);
+    if (icompute < 0)
+      error->all("Precompute ID for fix ave/time does not exist");
+    precompute = modify->compute[icompute];
+  } else precompute = NULL;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixAveTime::end_of_step()
+{
+  int i;
+
+  if (precompute) {
+    if (sflag) double tmp = precompute->compute_scalar();
+    if (vflag) precompute->compute_vector();
+  }
+
+  nsum++;
+  if (sflag) scalar += compute->compute_scalar();
+  if (vflag) {
+    compute->compute_vector();
+    double *cvector = compute->vector;
+    for (i = 0; i < size_vector; i++) vector[i] += cvector[i];
+  }
+
+  if (update->ntimestep % nfreq == 0) {
+    if (me == 0) {
+      fprintf(fp,"%d",update->ntimestep);
+      if (sflag) fprintf(fp," %g",scalar/nsum);
+      if (vflag)
+	for (i = 0; i < size_vector; i++) fprintf(fp," %g",vector[i]/nsum);
+      fprintf(fp,"\n");
+    }
+
+    nsum = 0;
+    scalar = 0.0;
+    if (vflag) 
+      for (i = 0; i < size_vector; i++) vector[i] = 0.0;
+ }
+}

src/fix_ave_time.h

@@ -0,0 +1,46 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef FIX_AVE_TIME_H
+#define FIX_AVE_TIME_H
+
+#include "stdio.h"
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixAveTime : public Fix {
+ public:
+  FixAveTime(class LAMMPS *, int, char **);
+  ~FixAveTime();
+  int setmask();
+  void init();
+  void end_of_step();
+
+ private:
+  int me;
+
+  int nfreq;
+  char *id_compute;
+  FILE *fp;
+
+  int sflag,vflag;
+  int size_vector,nsum;
+  double scalar,*vector;
+  class Compute *compute;
+  class Compute *precompute;
+};
+
+}
+
+#endif

src/min_cg.cpp

@@ -72,12 +72,12 @@ void MinCG::init()
   delete [] fixarg;
   fix_minimize = (FixMinimize *) modify->fix[modify->nfix-1];
 
-  // zero gradient vectors before first atom exchange
+  // zero local vectors before first atom exchange
 
-  setup_vectors();
+  set_local_vectors();
   for (int i = 0; i < ndof; i++) h[i] = g[i] = 0.0;
 
-  // virial_thermo is how virial should be computed on thermo timesteps
+  // virial_thermo = how virial computed on thermo timesteps
   // 1 = computed explicity by pair, 2 = computed implicitly by pair
 
   if (force->newton_pair) virial_thermo = 2;
@@ -127,28 +127,31 @@ void MinCG::init()
 
 void MinCG::run()
 {
-  double tmp,*f;
+  int i;
+  double tmp;
 
   // set initial force & energy
 
   setup();
-  setup_vectors();
   output->thermo->compute_pe();
-  ecurrent = output->thermo->potential_energy;
+  energy = output->thermo->potential_energy;
+  energy += energy_extra;
 
   // stats for Finish to print
 	
-  einitial = ecurrent;
+  einitial = energy;
 
-  f = atom->f[0];
   tmp = 0.0;
-  for (int i = 0; i < ndof; i++) tmp += f[i]*f[i];
+  for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
   MPI_Allreduce(&tmp,&gnorm2_init,1,MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < ndof_extra; i++) gnorm2_init += fextra[i]*fextra[i];
   gnorm2_init = sqrt(gnorm2_init);
 
   tmp = 0.0;
-  for (int i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
+  for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
   MPI_Allreduce(&tmp,&gnorminf_init,1,MPI_DOUBLE,MPI_MAX,world);
+  for (i = 0; i < ndof_extra; i++)
+    gnorminf_init = MAX(gnorminf_init,fabs(fextra[i]));
 
   // minimizer iterations
 
@@ -169,16 +172,19 @@ void MinCG::run()
     output->next_dump_any = update->ntimestep;
     if (output->restart_every) output->next_restart = update->ntimestep;
     output->next_thermo = update->ntimestep;
-    int ntmp;
+    eng_force();
-    double *xtmp,*htmp,etmp;
-    eng_force(&ntmp,&xtmp,&htmp,&etmp);
     output->write(update->ntimestep);
   }
   timer->barrier_stop(TIME_LOOP);
 
   // delete fix at end of run, so its atom arrays won't persist
+  // delete extra arrays
 
   modify->delete_fix("MINIMIZE");
+  delete [] xextra;
+  delete [] fextra;
+  delete [] gextra;
+  delete [] hextra;
 
   // reset reneighboring criteria
 
@@ -188,17 +194,20 @@ void MinCG::run()
 
   // stats for Finish to print
 	
-  efinal = ecurrent;
+  efinal = energy;
 
   f = atom->f[0];
   tmp = 0.0;
-  for (int i = 0; i < ndof; i++) tmp += f[i]*f[i];
+  for (i = 0; i < ndof; i++) tmp += f[i]*f[i];
   MPI_Allreduce(&tmp,&gnorm2_final,1,MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < ndof_extra; i++) gnorm2_final += fextra[i]*fextra[i];
   gnorm2_final = sqrt(gnorm2_final);
 
   tmp = 0.0;
-  for (int i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
+  for (i = 0; i < ndof; i++) tmp = MAX(fabs(f[i]),tmp);
   MPI_Allreduce(&tmp,&gnorminf_final,1,MPI_DOUBLE,MPI_MAX,world);
+  for (i = 0; i < ndof_extra; i++)
+    gnorminf_final = MAX(gnorminf_final,fabs(fextra[i]));
 }
 
 /* ----------------------------------------------------------------------
@@ -207,6 +216,20 @@ void MinCG::run()
 
 void MinCG::setup()
 {
+  // allocate extra arrays
+  // couldn't do in init(), b/c modify and fixes weren't yet init()
+  // set initial xextra values via fixes
+
+  ndof_extra = modify->min_dof();
+  xextra = new double[ndof_extra];
+  fextra = new double[ndof_extra];
+  gextra = new double[ndof_extra];
+  hextra = new double[ndof_extra];
+
+  modify->min_xinitial(xextra);
+
+  // perform usual setup
+
   if (comm->me == 0 && screen) fprintf(screen,"Setting up minimization ...\n");
 
   // setup domain, communication and neighboring
@@ -222,7 +245,7 @@ void MinCG::setup()
   comm->borders();
   neighbor->build();
   neighbor->ncalls = 0;
-  setup_vectors();
+  set_local_vectors();
 
   // compute all forces
 
@@ -247,6 +270,7 @@ void MinCG::setup()
   if (force->newton) comm->reverse_communicate();
 
   modify->setup();
+  energy_extra = modify->min_energy(xextra,fextra);
   output->setup(1);
 }
 
@@ -259,14 +283,14 @@ void MinCG::iterate(int n)
 {
   int i,gradsearch,fail;
   double alpha,beta,gg,dot[2],dotall[2];
-  double *f;
 
-  f = atom->f[0];
+  for (i = 0; i < ndof; i++) h[i] = g[i] = f[i];
-  for (int i = 0; i < ndof; i++) h[i] = g[i] = f[i];
+  for (i = 0; i < ndof_extra; i++) hextra[i] = gextra[i] = fextra[i];
 
   dot[0] = 0.0;
   for (i = 0; i < ndof; i++) dot[0] += f[i]*f[i];
   MPI_Allreduce(dot,&gg,1,MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < ndof_extra; i++) gg += fextra[i]*fextra[i];
 
   neval = 0;
   gradsearch = 1;
@@ -277,8 +301,8 @@ void MinCG::iterate(int n)
 
     // line minimization along direction h from current atom->x
 
-    eprevious = ecurrent;
+    eprevious = energy;
-    fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval);
+    fail = (this->*linemin)(neval);
 
     // if max_eval exceeded, all done
     // if linemin failed or energy did not decrease sufficiently:
@@ -287,8 +311,8 @@ void MinCG::iterate(int n)
 
     if (neval >= update->max_eval) break;
 
-    if (fail || fabs(ecurrent-eprevious) <= 
+    if (fail || fabs(energy-eprevious) <= 
-    	update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) {
+    	update->tolerance * 0.5*(fabs(energy) + fabs(eprevious) + EPS)) {
       if (gradsearch == 1) break;
       gradsearch = -1;
     }
@@ -299,13 +323,16 @@ void MinCG::iterate(int n)
     // force new search dir to be grad dir if need to restart CG
     // set gradsearch to 1 if will search in grad dir on next iteration
 
-    f = atom->f[0];
     dot[0] = dot[1] = 0.0;
     for (i = 0; i < ndof; i++) {
       dot[0] += f[i]*f[i];
       dot[1] += f[i]*g[i];
     }
     MPI_Allreduce(dot,dotall,2,MPI_DOUBLE,MPI_SUM,world);
+    for (i = 0; i < ndof_extra; i++) {
+      dotall[0] += fextra[i]*fextra[i];
+      dotall[1] += fextra[i]*gextra[i];
+    }
     
     beta = MAX(0.0,(dotall[0] - dotall[1])/gg);
     gg = dotall[0];
@@ -319,6 +346,10 @@ void MinCG::iterate(int n)
       g[i] = f[i];
       h[i] = g[i] + beta*h[i];
     }
+    for (i = 0; i < ndof_extra; i++) {
+      gextra[i] = fextra[i];
+      hextra[i] = gextra[i] + beta*hextra[i];
+    }
 
     // output for thermo, dump, restart files
 
@@ -334,9 +365,11 @@ void MinCG::iterate(int n)
    set ndof and vector pointers after atoms have migrated
 ------------------------------------------------------------------------- */
 
-void MinCG::setup_vectors()
+void MinCG::set_local_vectors()
 {
   ndof = 3 * atom->nlocal;
+  x = atom->x[0];
+  f = atom->f[0];
   if (ndof) g = fix_minimize->gradient[0];
   else g = NULL;
   if (ndof) h = fix_minimize->searchdir[0];
@@ -346,11 +379,11 @@ void MinCG::setup_vectors()
 /* ----------------------------------------------------------------------
    evaluate potential energy and forces
    may migrate atoms
-   new energy stored in ecurrent and returned (in case caller not in class)
+   energy = new objective function energy = poteng of atoms + eng_extra
-   negative gradient will be stored in atom->f
+   atom->f, fextra = negative gradient of objective function
 ------------------------------------------------------------------------- */
 
-void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng)
+void MinCG::eng_force()
 {
   // check for reneighboring
   // always communicate since minimizer moved atoms
@@ -375,7 +408,7 @@ void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng)
     timer->stamp(TIME_COMM);
     neighbor->build();
     timer->stamp(TIME_NEIGHBOR);
-    setup_vectors();
+    set_local_vectors();
   }
 
   // eflag is always set, since minimizer needs potential energy
@@ -409,21 +442,17 @@ void MinCG::eng_force(int *pndof, double **px, double **ph, double *peng)
     timer->stamp(TIME_COMM);
   }
 
-  // fixes that affect minimization
+  // min_post_force = forces on atoms that affect minimization
+  // min_energy = energy, forces on extra degrees of freedom
 
   if (modify->n_min_post_force) modify->min_post_force(vflag);
+  if (modify->n_min_energy) energy_extra = modify->min_energy(xextra,fextra);
 
   // compute potential energy of system via Thermo
 
   output->thermo->compute_pe();
-  ecurrent = output->thermo->potential_energy;
+  energy = output->thermo->potential_energy;
-
+  energy += energy_extra;
-  // return updated ptrs to caller since atoms may have migrated
-
-  *pndof = ndof;
-  *px = atom->x[0];
-  *ph = h;
-  *peng = ecurrent;
 }
 
 /* ----------------------------------------------------------------------
@@ -487,25 +516,15 @@ void MinCG::force_clear(int vflag)
 
 /* ----------------------------------------------------------------------
    line minimization methods
-   find minimum-energy starting at x along dir direction
+   find minimum-energy starting at x along h direction
-   input: n = # of degrees of freedom on this proc
+   update atom->x by alpha, call eng_force() for result
-          x = ptr to atom->x[0] as vector
+   alpha = distance moved along h to set x to minimun-energy configuration
-	  dir = search direction as vector
+   return 0 if successful move, 1 if failed (no move)
-	  eng = current energy at initial x
+   insure last call to eng_force() is consistent with return
-	  min/max dist = min/max distance to move any atom coord
-   output: return 0 if successful move, set alpha
-           return 1 if failed, no move, no need to set alpha
-           alpha = distance moved along dir to set x to min-eng config
-           caller has several quantities set via last call to eng_force()
-	     INSURE last call to eng_force() is consistent with returns
      if fail, eng_force() of original x
-	       if succeed, eng_force() at x + alpha*dir
+     if succeed, eng_force() at x + alpha*h
-             atom->x = coords at new configuration
+   eng_force() may migrate atoms due to neighbor list build
-	     atom->f = force (-Grad) is evaulated at new configuration
+     therefore linemin routines CANNOT store atom-based quantities
-	     ecurrent = energy of new configuration
-   NOTE: when call eng_force: n,x,dir,eng may change due to atom migration
-	 updated values are returned by eng_force()
-	 this routine CANNOT store atom-based quantities b/c of migration
 ------------------------------------------------------------------------- */
 
 /* ----------------------------------------------------------------------
@@ -515,52 +534,55 @@ void MinCG::force_clear(int vflag)
    quit as soon as energy starts to rise
 ------------------------------------------------------------------------- */
 
-int MinCG::linemin_scan(int n, double *x, double *dir, double eng,
+int MinCG::linemin_scan(int &nfunc)
-			double mindist, double maxdist,
-			double &alpha, int &nfunc)
 {
   int i;
-  double fmax,fme,elowest,alphamin,alphamax,alphalast;
+  double fmax,fme,elowest,alpha,alphamin,alphamax,alphalast;
 
   // alphamin = step that moves some atom coord by mindist
   // alphamax = step that moves some atom coord by maxdist
 
   fme = 0.0;
-  for (i = 0; i < n; i++) fme = MAX(fme,fabs(dir[i]));
+  for (i = 0; i < ndof; i++) fme = MAX(fme,fabs(h[i]));
   MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world);
-  if (fmax == 0.0) return 1;
+  for (i = 0; i < ndof_extra; i++) fmax = MAX(fmax,fabs(hextra[i]));
 
-  alphamin = mindist/fmax;
+  if (fmax == 0.0) return 1;
-  alphamax = maxdist/fmax;
+  alphamin = dmin/fmax;
+  alphamax = dmax/fmax;
 
   // if minstep is already uphill, fail
   // if eng increases, stop and return previous alpha
   // if alphamax, stop and return alphamax
 
-  elowest = eng;
+  elowest = energy;
   alpha = alphamin;
 
   while (1) {
-    for (i = 0; i < n; i++) x[i] += alpha*dir[i];
+    for (i = 0; i < ndof; i++) x[i] += alpha*h[i];
-    eng_force(&n,&x,&dir,&eng);
+    for (i = 0; i < ndof_extra; i++) xextra[i] += alpha*hextra[i];
+    eng_force();
     nfunc++;
 
-    if (alpha == alphamin && eng >= elowest) {
+    if (alpha == alphamin && energy >= elowest) {
-      for (i = 0; i < n; i++) x[i] -= alpha*dir[i];
+      for (i = 0; i < ndof; i++) x[i] -= alpha*h[i];
-      eng_force(&n,&x,&dir,&eng);
+      for (i = 0; i < ndof_extra; i++) xextra[i] -= alpha*hextra[i];
+      eng_force();
       nfunc++;
       return 1;
     }
-    if (eng > elowest) {
+    if (energy > elowest) {
-      for (i = 0; i < n; i++) x[i] += (alphalast-alpha)*dir[i];
+      for (i = 0; i < ndof; i++) x[i] += (alphalast-alpha)*h[i];
-      eng_force(&n,&x,&dir,&eng);
+      for (i = 0; i < ndof_extra; i++)
+	xextra[i] += (alphalast-alpha)*hextra[i];
+      eng_force();
       nfunc++;
       alpha = alphalast;
       return 0;
     }      
     if (alpha == alphamax) return 0;
 
-    elowest = eng;
+    elowest = energy;
     alphalast = alpha;
     alpha *= SCAN_FACTOR;
     if (alpha > alphamax) alpha = alphamax;
@@ -574,43 +596,44 @@ int MinCG::linemin_scan(int n, double *x, double *dir, double eng,
    prevents successvive func evals further apart in x than maxdist
 ------------------------------------------------------------------------- */
 
-int MinCG::linemin_secant(int n, double *x, double *dir, double eng,
+int MinCG::linemin_secant(int &nfunc)
-			  double mindist, double maxdist,
-			  double &alpha, int &nfunc)
 {
   int i,iter;
-  double eta,eta_prev,sigma0,sigmamax,alphadelta,fme,fmax,dsq,e0,tmp;
+  double eta,eta_prev,sigma0,sigmamax,alpha,alphadelta,fme,fmax,dsq,e0,tmp;
-  double *f;
   double epssq = SECANT_EPS * SECANT_EPS;
 
   // stopping criterion for secant iterations
 
   fme = 0.0;
-  for (i = 0; i < n; i++) fme += dir[i]*dir[i];
+  for (i = 0; i < ndof; i++) fme += h[i]*h[i];
   MPI_Allreduce(&fme,&dsq,1,MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < ndof_extra; i++) dsq += hextra[i]*hextra[i];
 
   // sigma0 = smallest allowed step of mindist
   // sigmamax = largest allowed step (in single iteration) of maxdist
 
   fme = 0.0;
-  for (i = 0; i < n; i++) fme = MAX(fme,fabs(dir[i]));
+  for (i = 0; i < ndof; i++) fme = MAX(fme,fabs(h[i]));
   MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world);
-  if (fmax == 0.0) return 1;
+  for (i = 0; i < ndof_extra; i++) fmax = MAX(fmax,fabs(hextra[i]));
 
-  sigma0 = mindist/fmax;
+  if (fmax == 0.0) return 1;
-  sigmamax = maxdist/fmax;
+  sigma0 = dmin/fmax;
+  sigmamax = dmax/fmax;
 
   // eval func at sigma0
   // test if minstep is already uphill
 
-  e0 = eng;
+  e0 = energy;
-  for (i = 0; i < n; i++) x[i] += sigma0*dir[i];
+  for (i = 0; i < ndof; i++) x[i] += sigma0*h[i];
-  eng_force(&n,&x,&dir,&eng);
+  for (i = 0; i < ndof_extra; i++) xextra[i] += sigma0*hextra[i];
+  eng_force();
   nfunc++;
 
-  if (eng >= e0) {
+  if (energy >= e0) {
-    for (i = 0; i < n; i++) x[i] -= sigma0*dir[i];
+    for (i = 0; i < ndof; i++) x[i] -= sigma0*h[i];
-    eng_force(&n,&x,&dir,&eng);
+    for (i = 0; i < ndof_extra; i++) xextra[i] -= sigma0*hextra[i];
+    eng_force();
     nfunc++;
     return 1;
   }
@@ -618,24 +641,25 @@ int MinCG::linemin_secant(int n, double *x, double *dir, double eng,
   // secant iterations
   // alphadelta = new increment to move, alpha = accumulated move
 
-  f = atom->f[0];
   tmp = 0.0;
-  for (i = 0; i < n; i++) tmp -= f[i]*dir[i];
+  for (i = 0; i < ndof; i++) tmp -= f[i]*h[i];
   MPI_Allreduce(&tmp,&eta_prev,1,MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < ndof_extra; i++) eta_prev -= fextra[i]*hextra[i];
 
   alpha = sigma0;
   alphadelta = -sigma0;
 
   for (iter = 0; iter < lineiter; iter++) {
     alpha += alphadelta;
-    for (i = 0; i < n; i++) x[i] += alphadelta*dir[i];
+    for (i = 0; i < ndof; i++) x[i] += alphadelta*h[i];
-    eng_force(&n,&x,&dir,&eng);
+    for (i = 0; i < ndof_extra; i++) xextra[i] += alphadelta*hextra[i];
+    eng_force();
     nfunc++;
 
-    f = atom->f[0];
     tmp = 0.0;
-    for (i = 0; i < n; i++) tmp -= f[i]*dir[i];
+    for (i = 0; i < ndof; i++) tmp -= f[i]*h[i];
     MPI_Allreduce(&tmp,&eta,1,MPI_DOUBLE,MPI_SUM,world);
+    for (i = 0; i < ndof_extra; i++) eta -= fextra[i]*hextra[i];
 
     alphadelta *= eta / (eta_prev - eta);
     eta_prev = eta;

src/min_cg.h

@@ -28,31 +28,38 @@ class MinCG : public Min {
 
  protected:
   int virial_thermo;         // what vflag should be on thermo steps (1,2)
-  int pairflag,torqueflag,granflag;
+  int pairflag,torqueflag,granflag;               // force clear flags
   int neigh_every,neigh_delay,neigh_dist_check;   // copies of reneigh criteria
 
   int maxpair;               // copies of Update quantities
   double **f_pair;
 
   class FixMinimize *fix_minimize;  // fix that stores gradient vecs
-  double ecurrent;            // current potential energy
   double mindist,maxdist;    // min/max dist for coord delta in line search
 
-  int ndof;               // # of degrees-of-freedom on this proc
+  int ndof;                  // 3N degrees-of-freedom on this proc
-  double *g,*h;           // local portion of gradient, searchdir vectors
+  double *x;                 // vec of 3N coords, ptr to atom->x[0]
+  double *f;                 // vec of 3N forces, ptr to atom->f[0]
+  double *g;                 // vec of 3N old forces, ptr to fix_minimize::g
+  double *h;                 // vec of 3N search dir, ptr to fix_minimize::h
 
-  typedef int (MinCG::*FnPtr)(int, double *, double *, double,
+  int ndof_extra;            // extra degrees of freedom to include in min
-			      double, double, double &, int &);
+  double energy_extra;       // extra potential energy
+  double *xextra;            // extra vecs associated with ndof_extra
+  double *fextra;
+  double *gextra;
+  double *hextra;
+
+  double energy;             // potential energy of atoms and extra dof
+
+  typedef int (MinCG::*FnPtr)(int &);
   FnPtr linemin;                        // ptr to linemin functions
-
+  int linemin_scan(int &);
-  int linemin_scan(int, double *, double *, double,
+  int linemin_secant(int &);
-		   double, double, double &, int &);
-  int linemin_secant(int, double *, double *, double,
-		     double, double, double &, int &);
 
   void setup();
-  void setup_vectors();
+  void set_local_vectors();
-  void eng_force(int *, double **, double **, double *);
+  void eng_force();
   void force_clear(int);
 };
 

src/min_cg_fr.cpp

@@ -39,14 +39,14 @@ void MinCGFR::iterate(int n)
 {
   int i,gradsearch,fail;
   double alpha,beta,gg,dot,dotall;
-  double *f;
 
-  f = atom->f[0];
   for (int i = 0; i < ndof; i++) h[i] = g[i] = f[i];
+  for (i = 0; i < ndof_extra; i++) hextra[i] = gextra[i] = fextra[i];
 
   dot = 0.0;
   for (i = 0; i < ndof; i++) dot += f[i]*f[i];
   MPI_Allreduce(&dot,&gg,1,MPI_DOUBLE,MPI_SUM,world);
+  for (i = 0; i < ndof_extra; i++) gg += fextra[i]*fextra[i];
 
   neval = 0;
   gradsearch = 1;
@@ -57,8 +57,8 @@ void MinCGFR::iterate(int n)
 
     // line minimization along direction h from current atom->x
 
-    eprevious = ecurrent;
+    eprevious = energy;
-    fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval);
+    fail = (this->*linemin)(neval);
 
     // if max_eval exceeded, all done
     // if linemin failed or energy did not decrease sufficiently:
@@ -67,8 +67,8 @@ void MinCGFR::iterate(int n)
 
     if (neval >= update->max_eval) break;
 
-    if (fail || fabs(ecurrent-eprevious) <= 
+    if (fail || fabs(energy-eprevious) <= 
-    	update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS)) {
+    	update->tolerance * 0.5*(fabs(energy) + fabs(eprevious) + EPS)) {
       if (gradsearch == 1) break;
       gradsearch = -1;
     }
@@ -79,10 +79,10 @@ void MinCGFR::iterate(int n)
     // force new search dir to be grad dir if need to restart CG
     // set gradsesarch to 1 if will search in grad dir on next iteration
 
-    f = atom->f[0];
     dot = 0.0;
     for (i = 0; i < ndof; i++) dot += f[i]*f[i];
     MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
+    for (i = 0; i < ndof_extra; i++) dotall += fextra[i]*fextra[i];
 
     beta = dotall/gg;
     gg = dotall;
@@ -96,6 +96,10 @@ void MinCGFR::iterate(int n)
       g[i] = f[i];
       h[i] = g[i] + beta*h[i];
     }
+    for (i = 0; i < ndof_extra; i++) {
+      gextra[i] = fextra[i];
+      hextra[i] = gextra[i] + beta*hextra[i];
+    }
 
     // output for thermo, dump, restart files
 

src/min_sd.cpp

@@ -35,10 +35,9 @@ void MinSD::iterate(int n)
 {
   int i,fail;
   double alpha,dot,dotall;
-  double *f;
 
-  f = atom->f[0];
   for (int i = 0; i < ndof; i++) h[i] = f[i];
+  for (i = 0; i < ndof_extra; i++) hextra[i] = fextra[i];
 
   neval = 0;
 
@@ -48,28 +47,29 @@ void MinSD::iterate(int n)
 
     // line minimization along direction h from current atom->x
 
-    eprevious = ecurrent;
+    eprevious = energy;
-    fail = (this->*linemin)(ndof,atom->x[0],h,ecurrent,dmin,dmax,alpha,neval);
+    fail = (this->*linemin)(neval);
 
     // if max_eval exceeded, all done
     // if linemin failed or energy did not decrease sufficiently, all done
 
     if (neval >= update->max_eval) break;
 
-    if (fail || fabs(ecurrent-eprevious) <= 
+    if (fail || fabs(energy-eprevious) <= 
-    	update->tolerance * 0.5*(fabs(ecurrent) + fabs(eprevious) + EPS))
+    	update->tolerance * 0.5*(fabs(energy) + fabs(eprevious) + EPS))
       break;
 
     // set h to new f = -Grad(x)
     // done if size sq of grad vector < EPS
 
-    f = atom->f[0];
     dot = 0.0;
     for (i = 0; i < ndof; i++) dot += f[i]*f[i];
     MPI_Allreduce(&dot,&dotall,1,MPI_DOUBLE,MPI_SUM,world);
+    for (i = 0; i < ndof_extra; i++) dotall += fextra[i]*fextra[i];
     if (dotall < EPS) break;
 
     for (i = 0; i < ndof; i++) h[i] = f[i];
+    for (i = 0; i < ndof_extra; i++) hextra[i] = fextra[i];
 
     // output for thermo, dump, restart files
 

src/style_class2.h

@@ -1,56 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifdef AngleInclude
-#include "angle_class2.h"
-#endif
-
-#ifdef AngleClass
-AngleStyle(class2,AngleClass2)
-#endif
-
-#ifdef BondInclude
-#include "bond_class2.h"
-#endif
-
-#ifdef BondClass
-BondStyle(class2,BondClass2)
-#endif
-
-#ifdef DihedralInclude
-#include "dihedral_class2.h"
-#endif
-
-#ifdef DihedralClass
-DihedralStyle(class2,DihedralClass2)
-#endif
-
-#ifdef ImproperInclude
-#include "improper_class2.h"
-#endif
-
-#ifdef ImproperClass
-ImproperStyle(class2,ImproperClass2)
-#endif
-
-#ifdef PairInclude
-#include "pair_lj_class2.h"
-#include "pair_lj_class2_coul_cut.h"
-#include "pair_lj_class2_coul_long.h"
-#endif
-
-#ifdef PairClass
-PairStyle(lj/class2,PairLJClass2)
-PairStyle(lj/class2/coul/cut,PairLJClass2CoulCut)
-PairStyle(lj/class2/coul/long,PairLJClass2CoulLong)
-#endif

src/style_dpd.h

@@ -1,28 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifdef AtomInclude
-#include "atom_vec_dpd.h"
-#endif
-
-#ifdef AtomClass
-AtomStyle(dpd,AtomVecDPD)
-#endif
-
-#ifdef PairInclude
-#include "pair_dpd.h"
-#endif
-
-#ifdef PairClass
-PairStyle(dpd,PairDPD)
-#endif

src/style_granular.h

@@ -1,50 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifdef AtomInclude
-#include "atom_vec_granular.h"
-#endif
-
-#ifdef AtomClass
-AtomStyle(granular,AtomVecGranular)
-# endif
-
-#ifdef FixInclude
-#include "fix_freeze.h"
-#include "fix_gran_diag.h"
-#include "fix_nve_gran.h"
-#include "fix_pour.h"
-#include "fix_shear_history.h"
-#include "fix_wall_gran.h"
-#endif
-
-#ifdef FixClass
-FixStyle(freeze,FixFreeze)
-FixStyle(gran/diag,FixGranDiag)
-FixStyle(nve/gran,FixNVEGran)
-FixStyle(pour,FixPour)
-FixStyle(SHEAR_HISTORY,FixShearHistory)
-FixStyle(wall/gran,FixWallGran)
-#endif
-
-#ifdef PairInclude
-#include "pair_gran_hertzian.h"
-#include "pair_gran_history.h"
-#include "pair_gran_no_history.h"
-#endif
-
-#ifdef PairClass
-PairStyle(gran/hertzian,PairGranHertzian)
-PairStyle(gran/history,PairGranHistory)
-PairStyle(gran/no_history,PairGranNoHistory)
-#endif

src/style_kspace.h

@@ -1,38 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifdef KSpaceInclude
-#include "ewald.h"
-#include "pppm.h"
-#include "pppm_tip4p.h"
-#endif
-
-#ifdef KSpaceClass
-KSpaceStyle(ewald,Ewald)
-KSpaceStyle(pppm,PPPM)
-KSpaceStyle(pppm/tip4p,PPPMTIP4P)
-#endif
-
-#ifdef PairInclude
-#include "pair_buck_coul_long.h"
-#include "pair_lj_cut_coul_long.h"
-#include "pair_lj_cut_coul_long_tip4p.h"
-#include "pair_lj_charmm_coul_long.h"
-#endif
-
-#ifdef PairClass
-PairStyle(buck/coul/long,PairBuckCoulLong)
-PairStyle(lj/cut/coul/long,PairLJCutCoulLong)
-PairStyle(lj/cut/coul/long/tip4p,PairLJCutCoulLongTIP4P)
-PairStyle(lj/charmm/coul/long,PairLJCharmmCoulLong)
-#endif

src/style_manybody.h

@@ -1,28 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifdef PairInclude
-#include "pair_eam.h"
-#include "pair_eam_alloy.h"
-#include "pair_eam_fs.h"
-#include "pair_sw.h"
-#include "pair_tersoff.h"
-#endif
-
-#ifdef PairClass
-PairStyle(eam,PairEAM)
-PairStyle(eam/alloy,PairEAMAlloy)
-PairStyle(eam/fs,PairEAMFS)
-PairStyle(sw,PairSW)
-PairStyle(tersoff,PairTersoff)
-#endif

src/style_molecule.h

@@ -1,116 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifdef AngleInclude
-#include "angle_charmm.h"
-#include "angle_cosine.h"
-#include "angle_cosine_squared.h"
-#include "angle_harmonic.h"
-#include "angle_hybrid.h"
-#endif
-
-#ifdef AngleClass
-AngleStyle(charmm,AngleCharmm)
-AngleStyle(cosine,AngleCosine)
-AngleStyle(cosine/squared,AngleCosineSquared)
-AngleStyle(harmonic,AngleHarmonic)
-AngleStyle(hybrid,AngleHybrid)
-#endif
-
-#ifdef AtomInclude
-#include "atom_vec_angle.h"
-#include "atom_vec_bond.h"
-#include "atom_vec_full.h"
-#include "atom_vec_molecular.h"
-#endif
-
-#ifdef AtomClass
-AtomStyle(angle,AtomVecAngle)
-AtomStyle(bond,AtomVecBond)
-AtomStyle(full,AtomVecFull)
-AtomStyle(molecular,AtomVecMolecular)
-#endif
-
-#ifdef BondInclude
-#include "bond_fene.h"
-#include "bond_fene_expand.h"
-#include "bond_harmonic.h"
-#include "bond_hybrid.h"
-#include "bond_morse.h"
-#include "bond_nonlinear.h"
-#include "bond_quartic.h"
-#endif
-
-#ifdef BondClass
-BondStyle(fene,BondFENE)
-BondStyle(fene/expand,BondFENEExpand)
-BondStyle(harmonic,BondHarmonic)
-BondStyle(hybrid,BondHybrid)
-BondStyle(morse,BondMorse)
-BondStyle(nonlinear,BondNonlinear)
-BondStyle(quartic,BondQuartic)
-#endif
-
-#ifdef DihedralInclude
-#include "dihedral_charmm.h"
-#include "dihedral_harmonic.h"
-#include "dihedral_helix.h"
-#include "dihedral_hybrid.h"
-#include "dihedral_multi_harmonic.h"
-#include "dihedral_opls.h"
-#endif
-
-#ifdef DihedralClass
-DihedralStyle(charmm,DihedralCharmm)
-DihedralStyle(harmonic,DihedralHarmonic)
-DihedralStyle(helix,DihedralHelix)
-DihedralStyle(hybrid,DihedralHybrid)
-DihedralStyle(multi/harmonic,DihedralMultiHarmonic)
-DihedralStyle(opls,DihedralOPLS)
-#endif
-
-#ifdef DumpInclude
-#include "dump_bond.h"
-#endif
-
-#ifdef DumpClass
-DumpStyle(bond,DumpBond)
-#endif
-
-#ifdef FixInclude
-#endif
-
-#ifdef FixClass
-#endif
-
-#ifdef ImproperInclude
-#include "improper_cvff.h"
-#include "improper_harmonic.h"
-#include "improper_hybrid.h"
-#endif
-
-#ifdef ImproperClass
-ImproperStyle(cvff,ImproperCvff)
-ImproperStyle(harmonic,ImproperHarmonic)
-ImproperStyle(hybrid,ImproperHybrid)
-#endif
-
-#ifdef PairInclude
-#include "pair_lj_charmm_coul_charmm.h"
-#include "pair_lj_charmm_coul_charmm_implicit.h"
-#endif
-
-#ifdef PairClass
-PairStyle(lj/charmm/coul/charmm,PairLJCharmmCoulCharmm)
-PairStyle(lj/charmm/coul/charmm/implicit,PairLJCharmmCoulCharmmImplicit)
-#endif

src/style_opt.h

@@ -1,30 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifdef PairInclude
-#include "pair_eam_opt.h"
-#include "pair_eam_alloy_opt.h"
-#include "pair_eam_fs_opt.h"
-#include "pair_lj_charmm_coul_long_opt.h"
-#include "pair_lj_cut_opt.h"
-#include "pair_morse_opt.h"
-#endif
-
-#ifdef PairClass
-PairStyle(eam/opt,PairEAMOpt)
-PairStyle(eam/alloy/opt,PairEAMAlloyOpt)
-PairStyle(eam/fs/opt,PairEAMFSOpt)
-PairStyle(lj/cut/opt,PairLJCutOpt)
-PairStyle(lj/charmm/coul/long/opt,PairLJCharmmCoulLongOpt)
-PairStyle(morse/opt,PairMorseOpt)
-#endif

src/style_xtc.h

@@ -1,20 +0,0 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifdef DumpInclude
-#include "dump_xtc.h"
-#endif
-
-#ifdef DumpClass
-DumpStyle(xtc,DumpXTC)
-#endif