else virial_thermo = 1;

  // set flags for what arrays to clear in force_clear()

  // need to clear torques if atom_style is dipole

  // need to clear phia if atom_style is granular

  // need to clear torques if array exists

  // don't need to clear f_pair if atom_style is only granular (no virial)

  torqueflag = 0;

  if (atom->check_style("dipole")) torqueflag = 1;

  granflag = 0;

  if (atom->check_style("granular")) granflag = 1;

  if (atom->torque) torqueflag = 1;

  pairflag = 1;

  if (strcmp(atom->atom_style,"granular") == 0) pairflag = 0;

    }

  }

  if (granflag) {

    double **phia = atom->phia;

    for (i = 0; i < nall; i++) {

      phia[i][0] = 0.0;

      phia[i][1] = 0.0;

      phia[i][2] = 0.0;

    }

  }

  // clear f_pair array if using it this timestep to compute virial

  if (vflag == 2 && pairflag) {

	     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

	 these routines CANNOT store atom-based quantities b/c of migration

------------------------------------------------------------------------- */

/* ----------------------------------------------------------------------

/* ----------------------------------------------------------------------

   linemin: use secant approximation to estimate parabola minimum at each step

   should converge more quickly/accurately than "scan", but may be less robust

   initial secant from two points: 0 and sigma0 = mindist

   prevents successvive func evals further apart in x than maxdist

   should converge more quickly/accurately than "scan", but can be less robust

------------------------------------------------------------------------- */

int MinCG::linemin_secant(int n, double *x, double *dir, double eng,

			  double &alpha, int &nfunc)

{

  int i,iter;

  double eta,eta_prev,sigma0,sigmamax,alphadelta,fme,fmax,dsq,e0,tmp;

  double eta,eta_prev,alphamin,alphamax,alphadelta,fme,fmax,dsq,e0,tmp;

  double *f;

  double epssq = SECANT_EPS * SECANT_EPS;

  for (i = 0; i < n; i++) fme += dir[i]*dir[i];

  MPI_Allreduce(&fme,&dsq,1,MPI_DOUBLE,MPI_SUM,world);

  // sigma0 = smallest allowed step of mindist

  // sigmamax = largest allowed step (in single iteration) of maxdist

  // alphamin = smallest allowed step of mindist

  // alphamax = largest allowed step (in single iteration) of maxdist

  fme = 0.0;

  for (i = 0; i < n; i++) fme = MAX(fme,fabs(dir[i]));

  MPI_Allreduce(&fme,&fmax,1,MPI_DOUBLE,MPI_MAX,world);

  if (fmax == 0.0) return 1;

  sigma0 = mindist/fmax;

  sigmamax = maxdist/fmax;

  alphamin = mindist/fmax;

  alphamax = maxdist/fmax;

  // eval func at sigma0

  // test if minstep is already uphill

  // eval func at alphamin

  // exit if minstep is already uphill

  e0 = eng;

  for (i = 0; i < n; i++) x[i] += sigma0*dir[i];

  for (i = 0; i < n; i++) x[i] += alphamin*dir[i];

  eng_force(&n,&x,&dir,&eng);

  nfunc++;

  if (eng >= e0) {

    for (i = 0; i < n; i++) x[i] -= sigma0*dir[i];

    for (i = 0; i < n; i++) x[i] -= alphamin*dir[i];

    eng_force(&n,&x,&dir,&eng);

    nfunc++;

    return 1;

  // secant iterations

  // alphadelta = new increment to move, alpha = accumulated move

  // first step is alpha = 0, first previous step is at mindist

  // prevent func evals for alpha outside mindist to maxdist

  // if happens on 1st iteration, secant approx is likely searching

  //   for a maximum (negative alpha), so reevaluate at alphamin

  f = atom->f[0];

  tmp = 0.0;

  for (i = 0; i < n; i++) tmp -= f[i]*dir[i];

  MPI_Allreduce(&tmp,&eta_prev,1,MPI_DOUBLE,MPI_SUM,world);

  alpha = sigma0;

  alphadelta = -sigma0;

  alpha = alphamin;

  alphadelta = -alphamin;

  for (iter = 0; iter < lineiter; iter++) {

    alpha += alphadelta;

    alphadelta *= eta / (eta_prev - eta);

    eta_prev = eta;

    if (alphadelta*alphadelta*dsq <= epssq) break;

    if (fabs(alphadelta) > sigmamax) {

      if (alphadelta > 0.0) alphadelta = sigmamax;

      else alphadelta = -sigmamax;

    if (alpha+alphadelta < alphamin || alpha+alphadelta > alphamax) {

      if (iter == 0) {

	alpha = alphamin;

	for (i = 0; i < n; i++) x[i] += alphamin*dir[i];

	eng_force(&n,&x,&dir,&eng);

	nfunc++;

      }

      break;

    }

  }

  // if exited loop on first iteration, func eval was at alpha = 0.0

  // else successful line search

  if (iter == 0) return 1;

  return 0;

}

 protected:

  int virial_thermo;          // what vflag should be on thermo steps (1,2)

  int pairflag,torqueflag,granflag;

  int pairflag,torqueflag;

  int neigh_every,neigh_delay,neigh_dist_check;   // copies of reneigh criteria

  int maxpair;                // copies of Update quantities