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: Sun, 2 Jul 2017 20:05:53 +0000

 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

/* ----------------------------------------------------------------------
   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]);
  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;
	if (force->pair){
		maxEnergy = -999999.0;
		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);			
				
				// SM (1Jul17)
				Compute *compute = modify->compute[icompute];
				compute->compute_peratom();
				double *compute_vector = compute->vector_atom;
				
				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;
	
	double boostAmount = 0.0;	// SM (29Jun17)
	
	 // neighbor->build_one(list->index);	
	
	 // inum = list->inum;
	 // ilist = list->ilist;
	 // numneigh = list->numneigh;
	 // firstneigh = list->firstneigh;
	
	for (ii = 0; ii < nlocal; ii++) {
		i = 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];

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

			  if (mask[j] & jgroupbit & j!=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];
					  
				  // SM: choose alphaB2 
				  if (currentTime < (desorptionTime + intervalTime) ) alphaB2 = 0.0;
				  else alphaB2 = -(alphaB*temp0);
							 
				  if (rsq <= 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*sigma12 - sigma6);		// SM (29Jun17)
						fpair = forcelj*r2inv;
					}

					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*dexp);
						fpair = 2.0 * param1 * param2 * (dexp*dexp - dexp) / r;
					}

					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+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)
				  }
			  }  
			}
		}			
	}
	
	// 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()
{
  // 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[0];
}

/* ----------------------------------------------------------------------
   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)
  
 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.

*/