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

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


From: S M <smatrix_2012@...8...>
Date: Sat, 1 Jul 2017 21:47:46 +0000


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.

Can anybody help me solve this problem?

Thanks,



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



From: jewett.aij@...24... <jewett.aij@...24...> on behalf of Andrew Jewett <jewett@...1937...>
Sent: Saturday, June 17, 2017 1:55 AM
To: S M
Cc: lammps
Subject: Re: [lammps-users] Reading Energy in Custom Pair Function
 


On Jun 16, 2017 5:45 PM, "S M" <smatrix_2012@...8...> wrote:
Hi,

I am writing a custom pair potential, where I first need to read the energy of every particle

Can you clarify what you mean here?
Is the pair style you are writing not already calculating the (pair contribution to) the energy per particle?  (... that is...excluding many-body forces or long range electrostatics.... Are you using these?) 
What contributions to the energy do you mean?
Contributions to the energy from bonds, angles, dihedrals... containing that atom?  (What about fixes like "fix spring" or "fix wall"?)


the custom pair potential is applied, and then apply the custom pair potential to certain particles based on their energy readings.

I am guessing that I need to use the "pair->vdW" command.

?
This is not an input script command (and this string does not appear anywhere in the source code)


Can anybody suggest a definite strategy in terms of writing necessary codes?

My guess would be what you are trying to do is not easy to do in a general way.  But perhaps if you clarify what you want, someone can suggest something.

Cheers

Andrew



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" 

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): 
  // -----------------
  
  // SM: set second group for interaction
  
  xstr = ystr = zstr = NULL;
  
  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 (29Jun17): set type of boost potential (LJ, Morse, vdW) and parameters and cutoff
  
  if (strstr(arg[4],"LJ") == arg[4]) potentialType = 1;
  else if (strstr(arg[4],"morse") == arg[4]) potentialType = 2;
  else if (strstr(arg[4],"vdW") == arg[4]) 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[5]);
  param2 = force->numeric(FLERR,arg[6]);
  param3 = force->numeric(FLERR,arg[7]);

  cutoff = force->numeric(FLERR,arg[8]);
  alphaB = force->numeric(FLERR,arg[9]);
  intervalTime = force->inumeric(FLERR,arg[10]);
  
  //--------------------------------------
  
  // optional args

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

  int iarg = 11;
  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;
  // delete [] group2;	// SM (29Jun17)
  // delete [] list;	// SM (29Jun17)
  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");
  }
  
  /*
  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;
	if (force->newton) npair += atom->nghost;
	double maxEnergy, currentEnergy;	
	tagint *molecule = atom->molecule;	  
	int *type = atom->type;	
	int desorptionTime;			  			
	int currentTime = update->ntimestep;	  	
	double alphaB2;			

	// 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){
		double *eatom = force->pair->eatom;
		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);			

				// for (j = 0; j < npair; j++) {
					
				fprintf(screen,"eatom[i]: \n");
				fprintf(screen,"eatom[i] = %f \n", eatom[i]);
				currentEnergy = eatom[i];
				
				// }
				if (i==0) maxEnergy = 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);
	
	//-------------------
  
    
/*
  // 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)
  
 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.

*/