add a new potential

You’re probably looping out-of-bounds in one of the arrays, or otherwise reading the wrong memory.

Tim

Thank you. I have checked my code and didn’t find any wrong loops. How can I check the memory reading condition?

double PairWatanabe::coordination(int i)
{
int j,jj,jnum;
double xtmp,ytmp,ztmp,delx,dely,delz,rij,c1,coord;
int ilist,jlist,numneigh,**firstneigh;
int nlocal=atom->nlocal;
// compute coordination number for each atom in group
// use full neighbor list to count atoms less than cutoff
double **x = atom->x;
int nall = atom->nlocal + atom->nghost;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
jlist = list->firstneigh[i];
jnum = list->numneigh[i];
coord = 0;
// printf("coord %d d\n",i,jnum); for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; // if (j >= nall) j = nall;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rij =sqrt( delx
delx + dely
dely + delz
delz);
if(rij<=(R-D)) coord+=1.0;
else if(rij<(R+D))
{ c1=(rij-R+D)/D;
coord+=1.0-0.5c1+sin(PIc1)/(2.0PI); }
else coord+=0.0;
// printf(“coord %d %d %f %f\n”,i,j,rij,coord);
}
return(coord);
}
void PairWatanabe::parag(int itype,double z,double dz,double &g,double &dg)
{
double invz,pinvz,rootz,expinv,fac,expo1,expo2,c1,c2;
if(itype==0) ///Si is type 0 and O is type 1
{
if(z<4.0)
{
rootz=sqrt(z+P2);
fac=P1
rootz-P4;
invz=1.0/(z-4.0);
pinvz=P3invz;
expinv=exp(pinvz);
g=fac
expinv+P4;
dg=expinvdz(0.5P1/rootz-invzpinvzfac);
}
else {g=P4;dg=0.0;}
}
else
{
c1=P8
(z-P9);
expo2=P5exp(c1(z-P9));
c2=exp((P6-z)/P7);
expo1=1.0/(c2+1.0);
g=expo2expo1;
dg=g
dz*(2.0c1+c2expo1/P7);
}
//printf(“g %d %f %f %f %f \n”,itype,z,dz,g,dg);
// g=1.0; dg=0.0;

}
/* ---------------------------------------------------------------------- /
void PairWatanabe::compute(int eflag, int vflag)
{
/////
ytmp = x[i][1];
ztmp = x[i][2];
zi=coordination(i);
printf(“i %d %f\n”,i,zi);
// two-body interactions, skip half of them
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++)
{
j = jlist[jj];
//if(j>=nlocal) j=j%nlocal;
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 = delx
delx + delydely + delzdelz;
r=sqrt(rsq);
ijparam = elem2param[itype][jtype][jtype];
if ((rsq - params[ijparam].cut2sq)>(-1e-6)) continue;
zj=coordination(j);
if(itype==jtype) {gij=1.0;dgij=0.0;}
else
{
parag(itype,zi,dz,gi,dgi);
parag(jtype,zj,dz,gj,dgj);
gij=gigj; dgij=dgigj+gi*dgj;
}

twobody(&params[ijparam],r,gij,dgij,fpair,eflag,evdwl);
// printf(“ij %d %d %d %d %f %f\n”,itype,jtype,i,j,fpair,evdwl);
f[i][0] += delxfpair;
f[i][1] += dely
fpair;
f[i][2] += delzfpair;
f[j][0] -= delx
fpair;
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);
}

Lan

2011/5/9 Timothy Sirk <tsirk@…31…>

Thank you. I have checked my code and didn't find any wrong loops. How can I
check the memory reading condition?

please compile with debug info included and run under the control of a debugger
to get the exact location where your code is segfaulting.

i remember giving detailed instructions on how to do this a couple times
to the mailing list, so please search the mailing list archives.

the fact that keeping printf() in avoids the error is a hint
that you are likely to have an out of bounds memory access.
the printf inhibits the compiler to optimize away some local
variables. the fact that the code "works" with the printf and
even may give correct results, doesn't mean it _is_ correct.

axel.

Axel is right, the debugger is really an essential tool. valgrind’s memcheck will also catch many memory problems, with minimum effort.

good point, however, also valgrind is dependent on compiling with -g
to be able to tell exactly where the memory access violation happened.

axel.

Yes- and memcheck can be noisy by producing alot of unwanted warnings. The debugger is really the best place to start.

Tim

Ok, I found my error. Thank you very much.

Lisa

2011/5/9 Timothy Sirk <tsirk@…31…>

I just installed the newest version of LAMMPS. When I debugged the new function with gdb, there are errors in the below. I even didn’t start the gdb debugger.

gdb ./lmp_linux

GNU gdb 6.6
Copyright © 2006 Free Software Foundation, Inc.
GDB is free software, covered by the GNU General Public License, and you are
welcome to change it and/or distribute copies of it under certain conditions.
Type “show copying” to see the conditions.
There is absolutely no warranty for GDB. Type “show warranty” for details.
This GDB was configured as “x86_64-suse-linux”…
Segmentation fault (core dumped)

When I use gdb for the old version of LAMMPS. It is ok. Is it a debug of new version?

Lisa

2011/5/10 Lisa Lan 5a <5alisalan@…24…>