[lammps-users] code error about a new potential

Hi,
I want to add a new potential to LAMMPS code and I wrote the code on base of stilling-weber potential. But when I compiled it and there is a error: there are 116640 atoms in my sample, but the j of the i list has a number of bigger than 116640. So during the calculation of coordination number of j, the program will break out. I am not sure about the meaning of j.
jlist = firstneigh[i];
jnum = numneigh[i];
printf(“i %d %d %d %d %f %f %f %f %f\n”,inum,jnum,type,i,zi,gi,xtmp,ytmp,ztmp);
for (jj = 0; jj < jnum; jj++)
{
j = jlist[jj];

}
what does j mean? should I use j for the calculation of coordination number and why the j is bigger than the atom number of sample?

Thank you.

Below is part of my code:
double PairWatanabe::coordination(int i)
{
int jjj,jj,jnum,ijparam,itype,jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,rij,rsq,fc,zi,c1;
int *numneigh,**firstneigh,*jlist;
double **x = atom->x;
int *type = atom->type;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itype=map[type[i]];

zi=0;
jlist = firstneigh[i];
jnum = numneigh[i];
printf(“cz %d %d %f %f %f %d\n”,i,itype,xtmp,ytmp,ztmp,jnum);
for (jj = 0; jj < jnum; jj++)
{
jjj = jlist[jj];
jtype = map[type[jjj]];
delx = xtmp - x[jjj][0];
dely = ytmp - x[jjj][1];
delz = ztmp - x[jjj][2];
rsq = delxdelx + delydely + delzdelz;
rij=sqrt(rsq);
ijparam = elem2param[itype][jtype][jtype];
c1=(rij-params[ijparam].bigr+params[ijparam].bigd)/params[ijparam].bigd;
if (rsq > params[ijparam].cutsq) continue;
if(rij<(params[ijparam].bigr-params[ijparam].bigd)) fc=1;
else if(rij<(params[ijparam].bigr+params[ijparam].bigd))
{ fc=1-0.5
c1+sin(PIc1)/(2.0PI);}
else fc=0;
zi+=fc;
printf(“z %d %d %f %f %f %f %f %f\n”,itype,jtype,params[ijparam].bigr,params[ijparam].bigd,rij,c1,fc,zi);
}
return(zi);

}
/* ---------------------------------------------------------------------- */
double PairWatanabe::parag(int itype,double zi,Param parami)
{
double gi;
if(itype==1)
{
if(zi<4) gi=(parami->p1
sqrt(zi+parami->p2)-parami->p4)*exp(parami->p3/(zi-4))+parami->p4;

else gi=parami->p4;
}
else
{gi=parami->p5exp(parami->p8pow(zi-parami->p9,2))/(exp((parami->p6-zi)/parami->p7)+1.0);}
printf(“g %d %f %f %f %f %f %f\n”,itype,gi,parami->p1,parami->p2,parami->p3,parami->p4,parami->p5);
return(gi);
}
/* ---------------------------------------------------------------------- */
void PairWatanabe::compute(int eflag, int vflag)
{
int i,j,k,ii,jj,kk,inum,jnum,jnumm1,itag,jtag,ktag;
int itype,jtype,ktype,ijparam,ikparam,ijkparam,iparam,jparam;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double zi,zj,gi,gj,gij;
double rsq,rsqij,rsqik,rij;
double delrij[3],delrik[3],fj[3],fk[3];
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
int *tag = atom->tag;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over full neighbor list of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
iparam=elem2param[itype][itype][itype];
zi=coordination(i);
gi=parag(itype,zi,&params[iparam]);

// two-body interactions, skip half of them
jlist = firstneigh[i];
jnum = numneigh[i];
printf("i %d %d %d %d %f %f %f %f f\n",inum,jnum,type,i,zi,gi,xtmp,ytmp,ztmp); for (jj = 0; jj < jnum; jj++) { **j = jlist[jj];** jtag = tag[j]; if (itag > jtag) {if ((itag+jtag) 2 == 0) continue;}
else if (itag < jtag) {if ((itag+jtag) % 2 == 1) continue;}
else
{
if (x[j][2] < ztmp) continue;
else if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
else if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp)
continue;
}
jtype = map[type[j]];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delxdelx + delydely + delz*delz;
rij=sqrt(rsq);
ijparam = elem2param[itype][jtype][jtype];
jparam=elem2param[jtype][jtype][jtype];
if (rsq > params[ijparam].cutsq) continue;
printf(“j %d %d %d %f %f %f\n”,jtype,j,jtag,x[j][0],x[j][1],x[j][2]);
zj=coordination(j);
gj=parag(jtype,zj,&params[jparam]);

if(itype!=jtype) gij=gigj;
else gij=1.0;
twobody(&params[ijparam],rsq,gij,fpair,eflag,evdwl);
printf(“ij %d %d %d %d %f %f %f\n”,itype,i,jtype,j,rij,gij,fpair);
f[i][0] += delx
fpair;
f[i][1] += delyfpair;
f[i][2] += delz
fpair;
f[j][0] -= delxfpair;
f[j][1] -= dely
fpair;
f[j][2] -= delz*fpair;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}

Hi,
I want to add a new potential to LAMMPS code and I wrote the code on base of
stilling-weber potential. But when I compiled it and there is a error: there
are 116640 atoms in my sample, but the j of the i list has a number of
bigger than 116640. So during the calculation of coordination number of j,
the program will break out. I am not sure about the meaning of j.

j index values larger than "nall" (=nlocal+nghost) - mind you the total
number of atoms doesn't matter for that, only the total number of atoms
on each MPI task - indicate an "exclusion", i.e. a pair of atoms that
has been flagged as special with a prefactor in the
special_bonds command.

just have a look at the pair_lj_cut.cpp code.

      if (j < nall) factor_lj = 1.0;
      else {
        factor_lj = special_lj[j/nall];
        j %= nall;
      }

so your "real" j is j % nall and the result of j/nall
tells you whether this is a 1-2, 1-3, or 1-4 interaction.

what you do with that depends on how you are
planning to handle bonded interactions.

axel.

Thank you, Axel. In my simulations, there are no special bonds between atoms and the potential consists of two body and three body part. It is very similar with the stilling-weber potential. I am still confusing that why stilling-weber potential is ok when it has no j/nall part while the new potential is not? Both of them don’t have the inner and outer part in the ij potential.

Best Regards,
Lisa

2011/3/18 Axel Kohlmeyer <[email protected]>

Thank you, Axel. In my simulations, there are no special bonds between atoms
and the potential consists of two body and three body part. It is very
similar with the stilling-weber potential. I am still confusing that why
stilling-weber potential is ok when it has no j/nall part while the new
potential is not? Both of them don't have the inner and outer part in the ij
potential.

well, how does your data file look? does it define bonds?

axel.

My data file is just for SiO2 and no special bonds are defined.

2011/3/18 Axel Kohlmeyer <[email protected]>