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.5c1+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->p1sqrt(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,¶ms[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,¶ms[jparam]);
if(itype!=jtype) gij=gigj;
else gij=1.0;
twobody(¶ms[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] += delxfpair;
f[i][1] += delyfpair;
f[i][2] += delzfpair;
f[j][0] -= delxfpair;
f[j][1] -= delyfpair;
f[j][2] -= delz*fpair;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}