LAMMPS WWW Site - LAMMPS Documentation - LAMMPS Mailing List Archives
Re: [lammps-users] Fw: Reading Energy in Custom Pair Function
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [lammps-users] Fw: Reading Energy in Custom Pair Function


From: S M <smatrix_2012@...8...>
Date: Tue, 11 Jul 2017 02:50:16 +0000

Hi, last week I wrote a fix to apply a boost force to a selected atom in a group. Now, I am trying to create a neighbor list in the same fix. I have used the following syntax (lines 346-353 in attached code):


class NeighList *list;

neighbor->build_one(list->index);

inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;

But I get "segmentation fault (core dumped)." I am guessing that the neighbor array is not allocated properly. I have looked at other fixes but have not found any example to follow. Can anybody help me correct my code?


The .cpp file is attached.


P.S. Also, for some reason, this fix functions about three times slower than a similar fix (e.g. fix addforce, which was modified into this fix). I am not sure what makes it so slow, and I doubt that the absence of a neighbor list is the reason (I tested by bypassing the loop). If anybody can spare the time to go through my code, please help me figure out how to make it more efficient.


Thanks,

S




From: Steve Plimpton <sjplimp@...24...>
Sent: Monday, July 3, 2017 10:36 AM
To: S M
Cc: Axel Kohlmeyer; lammps
Subject: Re: [lammps-users] Fw: Reading Energy in Custom Pair Function
 
Any fix can produce output which other fixes, computes, output access.
See fix indent as a simple example.

Steve

On Sun, Jul 2, 2017 at 2:05 PM, S M <smatrix_2012@...8...> wrote:

 Thanks, Axel, I have followed the example in "fix_ave_atom.cpp" and implemented necessary computes in my fix. It works as expected now.


*Anybody in the mailing list interested in the BOND-BOOST METHOD is welcome to use this fix*


On a less important note:

This fix writes the amount of repulsive energy applied ("boostAmount") to a text file, which is in turn read by other fixes (e.g. fix_nve.cpp). Instead, I am thinking of whether boostAmount can be stored as a variable and compute, which can be accessed by other fixes without the need of a separate text file.


Can anybody help me out here?

Thanks,
S



From: Axel Kohlmeyer <akohlmey@...24...>
Sent: Saturday, July 1, 2017 6:38 PM
To: S M
Cc: lammps
Subject: Re: [lammps-users] Fw: Reading Energy in Custom Pair Function
 
On Sat, Jul 1, 2017 at 5:47 PM, S M <smatrix_2012@...8...> wrote:
>
> Hi, here is my update:
>
> I am writing a fix (fix_addboost) that identifies the particle in a group
> possessing the highest pair potential energy, and then applies a repulsive
> potential to that particle only (this principle is used in the bond-boost
> method, especially of Miron-Fichthorn).
>
> My code compiles, but I get a "segmentation fault (core dumped)" error upon
> running. I have located the problem in the following lines (lines 277-278 in
> the full code attached) that try to read "eatoms[i]".
>
> fprintf(screen,"eatom[i] = %f \n", eatom[i]);
> currentEnergy = eatom[i];
>
> I am guessing that eatoms[i] have not been called properly. I have followed
> the syntax in "compute_pe_atom.cpp" and used the following in my code:
>
> if (force->pair){
> double *eatom = force->pair->eatom;
>         ....
> }
>
> I have also included pair.h and force.h at the beginning. I have gone
> through the LAMMPS developer manual also, but I still cannot figure out what
> is going wrong.

with your current setup, LAMMPS doesn't know, that you want to use
per-atom energies, so the array is not allocated and the energy not
tallied into that array.
the clean way to deal with this would be:
- create an instance of compute pe/atom with the suitable arguments
from within your fix
- read the data from the compute's vector_atom member after having
called its compute_peratom() method
- call modify->clearstep_compute() /
modify->addstep_compute(nextstep), where nextstep is the next
timestep, where you need access to the data from the compute.

see other fixes that maintain their own compute instances for more details.

axel.

>
> Can anybody help me solve this problem?
>
> Thanks,
> S
>
>
> ________________________________
> From: S M
> Sent: Sunday, June 18, 2017 2:16 PM
> To: lammps
>
> Subject: Re: [lammps-users] Reading Energy in Custom Pair Function
>
>
>
> Ok, here is the full scenario:
>
> I have a group of atoms at low temperature that interact with each other
> only via a LJ potential. At every timestep, I want to read the potential
> energy (from the LJ interaction) of every atom. Then, I want to identify the
> atom with the maximum potential energy, and weaken the forces acting on it,
> so that it is able to rapidly escape from the vicinity of the rest of the
> atoms.
>
> - When writing my previous email, I was thinking of applying a repulsive
> pair potential to the atom having the maximum potential energy, that would
> effectively weaken the original LJ interaction. That is why I thought about
> writing a custom pair potential that applies only to that specific atom in
> the group.
>
> - Now, I am thinking that the same can be achieved by applying an additional
> force on that same atom (and subsequent reaction forces on the other atoms),
> similar to the 'fix addforce' function. In this case, I would have to write
> a custom fix.
>
> In either of the above, I have to identify the atom with the maximum
> potential energy at every timestep.
>
> Can anybody provide some advice about how to proceed?
>
> Thanks,
> S


------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
lammps-users mailing list
lammps-users@...6297....sourceforge.net
https://lists.sourceforge.net/lists/listinfo/lammps-users


/* ----------------------------------------------------------------------
   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 "stdlib.h"
#include "fix_addboost.h"
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "domain.h"
#include "region.h"
#include "respa.h"
#include "input.h"
#include "variable.h"
#include "memory.h"
#include "error.h"
#include "force.h"

//SM (27Jun17)
#include "group.h"	
#include "timer.h"
#include "run.h"
#include <string>
#include "math.h"

// SM (29Jun17)
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "comm.h"
#include "pair.h" 
#include "compute.h"

using namespace LAMMPS_NS;
using namespace FixConst;

enum{NONE,CONSTANT,EQUAL,ATOM};

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

FixAddBoost::FixAddBoost(LAMMPS *lmp, int narg, char **arg) :
  Fix(lmp, narg, arg)
{
  if (narg < 6) error->all(FLERR,"Illegal fix addboost command");

  dynamic_group_allow = 1;
  scalar_flag = 1;
  vector_flag = 1;
  size_vector = 3;
  global_freq = 1;
  extscalar = 1;
  extvector = 1;

 
  // -----------------
  // SM (29Jul17): 
  // -----------------
  
  xstr = ystr = zstr = NULL;
  
  // SM: set second group for interaction
  int n = strlen(arg[3]) + 1;
  xstr = new char[n];
  strcpy(xstr,arg[3]);

  jgroup = group->find(xstr);
  if (jgroup == -1)
    error->all(FLERR,"FixAddboost second group ID does not exist");
  jgroupbit = group->bitmask[jgroup];
  
		// fprintf(screen,"Second group = %d \n", jgroup);
  
  // SM (1Jul17): Read per-atom energy  
  if (strncmp(arg[4],"c_",2) == 0) {
	  int n = strlen(arg[4]);
      char *suffix = new char[n];
      strcpy(suffix,&arg[4][2]);

      char *ptr = strchr(suffix,'[');
      if (ptr) {
        if (suffix[strlen(suffix)-1] != ']')
          error->all(FLERR,"Illegal fix addboost command");
        *ptr = '\0';
      }

      n = strlen(suffix) + 1;
      ystr = new char[n];
      strcpy(ystr,suffix);
      delete [] suffix;
	  
      icompute = modify->find_compute(ystr);
      if (icompute < 0)
        error->all(FLERR,"Compute ID for fix addboost does not exist");
      if (modify->compute[icompute]->peratom_flag == 0)
        error->all(FLERR,
                   "Fix addboost compute does not calculate per-atom values");
  }
  
  // SM (29Jun17): set type of boost potential (LJ, Morse, vdW) and parameters and cutoff
  if (strstr(arg[5],"LJ") == arg[5]) potentialType = 1;
  else if (strstr(arg[5],"morse") == arg[5]) potentialType = 2;
  else if (strstr(arg[5],"vdW") == arg[5]) potentialType = 3;  
  else error->all(FLERR,"Potential Type for fix addboost must be LJ, morse or vdW");

		// fprintf(screen,"Potential Type = %d \n", potentialType);
  
  param1 = force->numeric(FLERR,arg[6]);
  param2 = force->numeric(FLERR,arg[7]);
  param3 = force->numeric(FLERR,arg[8]);

  cutoff = force->numeric(FLERR,arg[9]);
  
		// fprintf(screen,"cutoff = %f \n", cutoff);
  
  alphaB = force->numeric(FLERR,arg[10]);
  intervalTime = force->inumeric(FLERR,arg[11]);
  
  //--------------------------------------
  
  // optional args

  iregion = -1;
  idregion = NULL;
  estr = NULL;

  int iarg = 12;
  while (iarg < narg) {
    if (strcmp(arg[iarg],"region") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix addboost command (insufficient arguments for region)");
      iregion = domain->find_region(arg[iarg+1]);
      if (iregion == -1)
        error->all(FLERR,"Region ID for fix addboost does not exist");
        int n = strlen(arg[iarg+1]) + 1;
        idregion = new char[n];
        strcpy(idregion,arg[iarg+1]);
      iarg += 2;
	  
		// fprintf(screen,"Region = %d \n", iregion);	  
	  
    } else if (strcmp(arg[iarg],"energy") == 0) {
      if (iarg+2 > narg) error->all(FLERR,"Illegal fix addboost command");
      if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
        int n = strlen(&arg[iarg+1][2]) + 1;
        estr = new char[n];
        strcpy(estr,&arg[iarg+1][2]);
      } else error->all(FLERR,"Illegal fix addboost command");
      iarg += 2;
    } else error->all(FLERR,"Illegal fix addboost command");
  }

  force_flag = 0;
  foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;

  maxatom = atom->nmax;
  memory->create(sforce,maxatom,4,"addboost:sforce");
}

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

FixAddBoost::~FixAddBoost()
{
  delete [] xstr;
  delete [] ystr;
  delete [] zstr;
  delete [] estr;
  delete [] idregion;
  memory->destroy(sforce);
}

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

int FixAddBoost::setmask()
{
  int mask = 0;
  mask |= POST_FORCE;
  mask |= THERMO_ENERGY;
  mask |= POST_FORCE_RESPA;
  mask |= MIN_POST_FORCE;
  return mask;
}

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

void FixAddBoost::init()
{

  // SM (29Jun17): deleted variable check

  // set index and check validity of region
	if (iregion >= 0) {
		iregion = domain->find_region(idregion);
		if (iregion == -1)
		  error->all(FLERR,"Region ID for fix addboost does not exist");
	}
	
	// SM (1Jul17): check validity of compute
	int icompute = modify->find_compute(ystr);
	if (icompute < 0)
		error->all(FLERR,"Compute ID for fix addboost does not exist");

	
  /*
  if (xstyle == ATOM || ystyle == ATOM || zstyle == ATOM)
    varflag = ATOM;
  else if (xstyle == EQUAL || ystyle == EQUAL || zstyle == EQUAL)
    varflag = EQUAL;
  else varflag = CONSTANT;

  if (varflag == CONSTANT && estyle != NONE)
    error->all(FLERR,"Cannot use variable energy with "
               "constant force in fix addboost");
  if ((varflag == EQUAL || varflag == ATOM) &&
      update->whichflag == 2 && estyle == NONE)
    error->all(FLERR,"Must use variable energy with fix addboost");

  if (strstr(update->integrate_style,"respa"))
    nlevels_respa = ((Respa *) update->integrate)->nlevels;
  */

}

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

void FixAddBoost::setup(int vflag)
{
  if (strstr(update->integrate_style,"verlet"))
    post_force(vflag);
  else {
    ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
    post_force_respa(vflag,nlevels_respa-1,0);
    ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
  }
}

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

void FixAddBoost::min_setup(int vflag)
{
  post_force(vflag);
}

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

void FixAddBoost::post_force(int vflag)
{
  double **x = atom->x;
  double **f = atom->f;
  int *mask = atom->mask;
  imageint *image = atom->image;
  int nlocal = atom->nlocal;

 
  //-------------------
  // SM (26Jun17)
  //-------------------
	int i,j;
	int npair = nlocal;
	double maxEnergy, currentEnergy;	
	tagint *molecule = atom->molecule;	  
	int *type = atom->type;	
	int desorptionTime;			  			
	int currentTime = update->ntimestep;	  	
	double alphaB2;			
	
	modify->clearstep_compute();  				// SM (1Jul17)
	// modify->addstep_compute_all(currentTime);	// SM (1Jul17)
	
	// SM: read initial temperature
	float temp0;
	FILE *fileTemp;
	fileTemp = fopen("current_temp.txt" , "r");
	if (fileTemp==NULL) error->all(FLERR,"Error opening current_temp.txt"); 
	else 
	while (! feof(fileTemp) ) fscanf(fileTemp,"%f\n",&temp0); 
	fclose(fileTemp);

	// SM: read latest desorption time
	FILE *myFile;
	myFile = fopen("desorption_time.txt" , "r");
	if (myFile==NULL) error->all(FLERR,"Error opening desorption_time.txt"); 
	else 
	while (! feof(myFile) ) fscanf(myFile,"%d\n",&desorptionTime); 
	fclose(myFile);

	// SM (26Jun17): Identify the particle which is located inside given region and has highest pair energy
	int molID;
	int maxE = 0;

	// SM (5Jul17): choose alphaB2 and set flag for bypass
	bool bypass;
	if (currentTime < (desorptionTime + intervalTime) ) {
		alphaB2 = 0.0;
		bypass = true;
	}
	else {
		alphaB2 = -(alphaB*temp0);	
		bypass  = false;
	}
		
	maxEnergy = -999999.0;
	if (bypass == false) {		// SM (5Jul17): bypass loop if boost = 0
		// SM (1Jul17) [see fix_ave_atom.cpp]
		Compute *compute = modify->compute[icompute];
		compute->compute_peratom();
		double *compute_vector = compute->vector_atom;

		for (i = 0; i < nlocal; i++) {
		  if (mask[i] & groupbit) {
			if ( domain->regions[iregion]->match(x[i][0],x[i][1],x[i][2]) ) {			
					// fprintf(screen,"Region match: particle %d \n", i);			
				
				currentEnergy = compute_vector[i];
					// fprintf(screen,"Current energy = %f \n", currentEnergy);
				
				if (currentEnergy > maxEnergy) {
					maxE = i;
					maxEnergy = currentEnergy;
				}
			}
		  } 
		}  
	}
	molID = molecule[maxE];
	
				// fprintf(screen,"Boost Molecule = %d \n", maxE);
	
    // SM (29Jun17): apply boost pair potential
	int ii,jj,inum,jnum,itype,jtype;
	double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
	double rsq,r2inv,r6inv,forcelj,factor_lj,energylj;
	int *ilist,*jlist,*numneigh,**firstneigh;
	
	boostAmount = 0.0;	// SM (29Jun17)

	/*
	class NeighList *list;	

	neighbor->build_one(list->index);	

	inum = list->inum;
	ilist = list->ilist;
	numneigh = list->numneigh;
	firstneigh = list->firstneigh;
	*/
	
	/*
	for (ii = 0; ii < nlocal; ii++) {
		i = ilist[ii];
		if ((mask[i] & groupbit) && molecule[i]==molID) {
			xtmp = x[i][0];
			ytmp = x[i][1];
			ztmp = x[i][2];
			itype = type[i];
			// jlist = firstneigh[i];
			// jnum = numneigh[i];
			
					// fprintf(screen,"molID[i] matched = %d \n", i);
	*/
		
	i = maxE;
	
	xtmp = x[i][0];
	ytmp = x[i][1];
	ztmp = x[i][2];

	if (bypass == false) {		// SM (5Jul17): bypass loop if boost = 0
			// fprintf(screen,"   bypass = FALSE \n");	

		for (jj = 0; jj < nlocal; jj++) {
		  j = jj;
		  // j &= NEIGHMASK;

		  if ((mask[j] & jgroupbit) && molecule[j]!=molecule[i]) {
			  delx = xtmp - x[j][0];
			  dely = ytmp - x[j][1];
			  delz = ztmp - x[j][2];
			  rsq = delx*delx + dely*dely + delz*delz;
			  jtype = type[j];

					// fprintf(screen,"   rsq = %f \n", rsq);	
			  
			  if (rsq <= cutoff*cutoff) {

				if (potentialType == 1) {
					// LJ: param1 = epsilon; param2 = sigma
					double sigma6, sigma12;		// SM (29Jun17)
					r2inv = 1.0/rsq;
					r6inv = r2inv*r2inv*r2inv;
					sigma6 = param2*param2*param2*param2*param2*param2;
					sigma12 = sigma6*sigma6;
					energylj = 4.0 * param1 * r6inv * (sigma12*r6inv - sigma6);		// SM (29Jun17)
					forcelj = 24.0 * param1 * r6inv * (2.0*sigma12*r6inv - sigma6);	// SM (29Jun17)
					fpair = forcelj*r2inv;
						// fprintf(screen,"   particle[j] = %d \n", j);	
				}

				else if (potentialType == 2) {
					// MORSE: param1 = De; param2 = alpha; param3 = R0
					double r = sqrt(rsq);
					double dr = r - param3;
					double dexp = exp(-param2 * dr);
					energylj = param1 * (dexp*dexp - 2.0*dexp);
					fpair = 2.0 * param1 * param2 * (dexp*dexp - dexp) / r;
				}

				else if (potentialType == 3) {
					// GRIMME'S vdW: param1 = c6; param2 = Rr; param3 = s6
					double dscale = 20.0;
					double r = sqrt(rsq);
					double ddexp = exp(-dscale * (r/param2 - 1.0));
					double invexp = 1.0/(1.0+ddexp);
					r2inv = 1.0/rsq;
					r6inv = r2inv*r2inv*r2inv;
					energylj = -param3*param1*r6inv*invexp;					
					fpair = param3 * param1 * (invexp*invexp*ddexp*(-dscale/param2)*r6inv) / r;
					fpair -= param3 * param1 * (6.0 * invexp * r6inv * r2inv);
				}
				
				fpair *= alphaB2;	// SM (22Apr17)

				f[i][0] += delx*fpair;
				f[i][1] += dely*fpair;
				f[i][2] += delz*fpair;
			
				f[j][0] -= delx*fpair;
				f[j][1] -= dely*fpair;
				f[j][2] -= delz*fpair;
				
				boostAmount += energylj*alphaB2;	// SM (29Jun17)
				
					// fprintf(screen,"Boost Amount = %f \n", boostAmount);
			  }
		  }  
		}
	}
		
	// SM (29Jun17): write new boost in file
	FILE *FP;
	FP = fopen("current_boost.txt" , "w");
	fprintf(FP, "%f", boostAmount);
	fclose(FP);

	// SM (29Jun17): append boost every 1 timestep in file
	FILE *FP2;
	FP2 = fopen("total_boost.dat", "a");
	fprintf(FP2, "%f \n", boostAmount);
	fclose(FP2);
	
    modify->addstep_compute(currentTime+1);		// SM (1Jul17)
    modify->addstep_compute_all(currentTime+1);	// SM (1Jul17)	 

	//-------------------	 
	
 /*
  // reallocate sforce array if necessary

  if ((varflag == ATOM || estyle == ATOM) && nlocal > maxatom) {
    maxatom = atom->nmax;
    memory->destroy(sforce);
    memory->create(sforce,maxatom,4,"addboost:sforce");
  }

  // foriginal[0] = "potential energy" for added force
  // foriginal[123] = force on atoms before extra force added

  foriginal[0] = foriginal[1] = foriginal[2] = foriginal[3] = 0.0;
  force_flag = 0;

  // constant force
  // potential energy = - x dot f in unwrapped coords

  if (varflag == CONSTANT) {
    double unwrap[3];
    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit & molecule[i]==molID) {	// SM (26Jun17)
        if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
        domain->unmap(x[i],image[i],unwrap);
        foriginal[0] -= xvalue*unwrap[0] + yvalue*unwrap[1] + zvalue*unwrap[2];
        foriginal[1] += f[i][0];
        foriginal[2] += f[i][1];
        foriginal[3] += f[i][2];
        f[i][0] += xvalue;
        f[i][1] += yvalue;
        f[i][2] += zvalue;
	  }

  // variable force, wrap with clear/add
  // potential energy = evar if defined, else 0.0
  // wrap with clear/add

  } else {

    modify->clearstep_compute();

    if (xstyle == EQUAL) xvalue = input->variable->compute_equal(xvar);
    else if (xstyle == ATOM)
      input->variable->compute_atom(xvar,igroup,&sforce[0][0],4,0);
    if (ystyle == EQUAL) yvalue = input->variable->compute_equal(yvar);
    else if (ystyle == ATOM)
      input->variable->compute_atom(yvar,igroup,&sforce[0][1],4,0);
    if (zstyle == EQUAL) zvalue = input->variable->compute_equal(zvar);
    else if (zstyle == ATOM)
      input->variable->compute_atom(zvar,igroup,&sforce[0][2],4,0);
    if (estyle == ATOM)
      input->variable->compute_atom(evar,igroup,&sforce[0][3],4,0);

    modify->addstep_compute(update->ntimestep + 1);

    for (int i = 0; i < nlocal; i++)
      if (mask[i] & groupbit & molecule[i]==molID) { // SM (26Jun17)
        if (region && !region->match(x[i][0],x[i][1],x[i][2])) continue;
        if (estyle == ATOM) foriginal[0] += sforce[i][3];
        foriginal[1] += f[i][0];
        foriginal[2] += f[i][1];
        foriginal[3] += f[i][2];
        if (xstyle == ATOM) f[i][0] += sforce[i][0];
        else if (xstyle) f[i][0] += xvalue;
        if (ystyle == ATOM) f[i][1] += sforce[i][1];
        else if (ystyle) f[i][1] += yvalue;
        if (zstyle == ATOM) f[i][2] += sforce[i][2];
        else if (zstyle) f[i][2] += zvalue;
      }
  }
  */
}

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

void FixAddBoost::post_force_respa(int vflag, int ilevel, int iloop)
{
  if (ilevel == nlevels_respa-1) post_force(vflag);
}

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

void FixAddBoost::min_post_force(int vflag)
{
  post_force(vflag);
}

/* ----------------------------------------------------------------------
   potential energy of added force
------------------------------------------------------------------------- */

double FixAddBoost::compute_scalar()
{
  // SM (5Jul17):
  return boostAmount;
}

/* ----------------------------------------------------------------------
   return components of total force on fix group before force was changed
------------------------------------------------------------------------- */

double FixAddBoost::compute_vector(int n)
{
  // only sum across procs one time

  if (force_flag == 0) {
    MPI_Allreduce(foriginal,foriginal_all,4,MPI_DOUBLE,MPI_SUM,world);
    force_flag = 1;
  }
  return foriginal_all[n+1];
}

/* ----------------------------------------------------------------------
   memory usage of local atom-based array
------------------------------------------------------------------------- */

double FixAddBoost::memory_usage()
{
  double bytes = 0.0;
  if (varflag == ATOM) bytes = atom->nmax*4 * sizeof(double);
  return bytes;
}
/* ----------------------------------------------------------------------
   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 FIX_CLASS

FixStyle(addboost,FixAddBoost)

#else

#ifndef LMP_FIX_ADDBOOST_H
#define LMP_FIX_ADDBOOST_H

#include "fix.h"

namespace LAMMPS_NS {

class FixAddBoost : public Fix {
 public:
  FixAddBoost(class LAMMPS *, int, char **);
  ~FixAddBoost();
  int setmask();
  void init();
  void setup(int);
  void min_setup(int);
  void post_force(int);
  void post_force_respa(int, int, int);
  void min_post_force(int);
  double compute_scalar();
  double compute_vector(int);
  double memory_usage();
  
  double alphaB, param1, param2, param3, cutoff;	// SM (28Jun17)
  int intervalTime;									// SM (29Jun17)	
  // char *group2;										// SM (29Jun17)		
  int jgroup,jgroupbit;								// SM (29Jun17)	
  int potentialType;								// SM (29Jun17)	
  // class NeighList *list;							// SM (29Jun17)
  int icompute;										// SM (1Jul17)
  double boostAmount;								// SM (5Jul17)
  
 private:
  double xvalue,yvalue,zvalue;
  int varflag,iregion;
  char *xstr,*ystr,*zstr,*estr;
  char *idregion;
  int xvar,yvar,zvar,evar,xstyle,ystyle,zstyle,estyle;
  double foriginal[4],foriginal_all[4];
  int force_flag;
  int nlevels_respa;

  int maxatom;
  double **sforce;
};

}

#endif
#endif

/* ERROR/WARNING messages:

E: Illegal ... command

Self-explanatory.  Check the input script syntax and compare to the
documentation for the command.  You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.

E: Region ID for fix addforce does not exist

Self-explanatory.

E: Variable name for fix addforce does not exist

Self-explanatory.

E: Variable for fix addforce is invalid style

Self-explanatory.

E: Cannot use variable energy with constant force in fix addforce

This is because for constant force, LAMMPS can compute the change
in energy directly.

E: Must use variable energy with fix addforce

Must define an energy vartiable when applyting a dynamic
force during minimization.

*/