Dear lammps user,
First of all I wanted to thanks in advance.
I’m trying to write my own pair style.
In particular it should be the combination of a Lennard jones 5-10 and a yukawa pair style.
I wrote the .h and .cpp file and I put the in the EXTRA-PAIR directory of my lammps distribution directory.
The installation of the EXTRA-PAIR package and the re-compilation of lammps via the command make.“machines” works well.
When I try to use this new pair style in my script this leads to an error in the definition of coefficients.
the error is : segmentation fault
With some test my guess that the failure is due to the neighbour list or the coefficients but I don’t understand where.
I’m using the distribution of 29-SEPT-2021.
I’m trying this only implementig the basic stuff, as write in lammps dopcumentation, such as allocate, setting, coefficient and compute.
I’ll attach the script below.
Many thanks in advance
Daniele
file .h pair_daniele.h
#ifdef PAIR_CLASS
// clang-format off
PairStyle(daniele,PairDaniele);
// clang-format on
#else
#ifndef LMP_PAIR_DANIELE_H
#define LMP_PAIR_DANIELE_H
#include “pair.h”
namespace LAMMPS_NS {
class PairDaniele : public Pair {
public:
PairDaniele(class LAMMPS *);
virtual ~PairDaniele();
virtual void compute(int, int);
virtual void settings(int, char **);
virtual void coeff(int, char **);
protected:
double kappa;
double cut_global;
double **cut;
double **a, **epsilon, **sigma;
double **lj1,**lj2,**lj3,**offset;
void allocate();
};
} // namespace LAMMPS_NS
#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: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair cutoff < Respa interior cutoff
One or more pairwise cutoffs are too short to use with the specified
rRESPA cutoffs.
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: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
*/
######### file .cpp pair_daniele.cpp ##############
// clang-format off
#include “pair_daniele.h”
#include
#include
#include “atom.h”
#include “comm.h”
#include “force.h”
#include “neighbor.h”
#include “neigh_list.h”
#include “neigh_request.h”
#include “update.h”
#include “math_const.h”
#include “memory.h”
#include “error.h”
using namespace LAMMPS_NS;
using namespace MathConst;
/* ---------------------------------------------------------------------- */
PairDaniele::PairDaniele(LAMMPS lmp) : Pair(lmp)
{
writedata = 0;
}
/ ---------------------------------------------------------------------- */
PairDaniele::~PairDaniele()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(cut);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(lj1);
memory->destroy(lj2);
memory->destroy(lj3);
memory->destroy(offset);
memory->destroy(a);
}
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairDaniele::allocate()
{
allocated = 1;
int n = atom->ntypes+1;
memory->create(setflag,n,n,“pair:setflag”);
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(cutsq,n,n,“pair:cutsq”);
memory->create(cut,n,n,“pair:cut”);
memory->create(epsilon,n,n,“pair:epsilon”);
memory->create(sigma,n,n,“pair:sigma”);
memory->create(lj1,n,n,“pair:lj1”);
memory->create(lj2,n,n,“pair:lj2”);
memory->create(lj3,n,n,“pair:lj3”);
memory->create(a,n,n,“pair:a”);
memory->create(offset,n,n,“pair:offset”);
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairDaniele::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,“Illegal pair_style command”);
kappa = utils::numeric(FLERR,arg[0],false,lmp);
cut_global = utils::numeric(FLERR,arg[1],false,lmp);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairDaniele::coeff(int narg, char **arg)
{
if (narg < 4 || narg > 5)
error->all(FLERR,“Incorrect args for pair coefficients”);
if (!allocated) allocate();
int ilo,ihi,jlo,jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp);
double sigma_one = utils::numeric(FLERR,arg[3],false,lmp);
double a_one = utils::numeric(FLERR,arg[4],false,lmp);
double cut_one = cut_global;
if (narg == 8) cut_one = utils::numeric(FLERR,arg[5],false,lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
a[i][j] = a_one;
cut[i][j] = cut_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR,“Incorrect args for pair coefficients”);
}
/* force computing */
void PairDaniele::compute(int eflag, int vflag)
{
int i,j,ii,jj,inum,jnum,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double r,rinv,screening,rsq,r2inv,rgamR,rgamA,forcedaniele,factor;
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
ev_init(eflag,vflag);
double **x = atom->x;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
double *special_lj = force->special_lj;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
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 < jnum; jj++) {
j = jlist[jj];
factor = special_lj[sbmask(j)];
j &= NEIGHMASK;
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];
if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq;
r = sqrt(rsq);
rinv = 1.0/r;
screening = exp(-kappa*r);
lj1[itype][jtype]=pow((sigma[itype][jtype]/r),10);
lj2[itype][jtype]=pow((sigma[itype][jtype]/r),5);
lj3[itype][jtype]=4*epsilon[itype][jtype];
forcedaniele =(a[itype][jtype]*screening*(r2inv+(rinv*kappa)))+(lj3[itype][jtype]*rinv*((10*lj1[itype][jtype])-(5*lj2[itype][jtype])));
fpair = factor*forcedaniele;
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
}
if (eflag) {
evdwl = ((a[itype][jtype] * screening * rinv)+(lj3[itype][jtype] * (lj1[itype][jtype]-lj2[itype][jtype])) -
offset[itype][jtype]);
evdwl *= factor;
}
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}