Hi Steve:

It’s still about the task of enabling quartic bond being used with angle potential. In the bond_quartic.cpp, there are a few lines of codes doing the job of subtracting pairwise potential of two bonded atoms:

// subtract out pairwise contribution from 2 atoms via pair->single()

// required since special_bond = 1,1,1

// tally energy/virial in pair, using newton_bond as newton flag

itype = atom->type[i1];

jtype = atom->type[i2];

if (rsq < cutsq[itype][jtype]) {

evdwl = -force->pair->single(i1,i2,itype,jtype,rsq,1.0,1.0,fpair);

fpair = -fpair;

if (newton_bond || i1 < nlocal) {

f[i1][0] += delx*fpair; //delx = x[i1][0]-x[i2][0]*

f[i1][1] += delyfpair; //dely = x[i1][1]-x[i2][1]

f[i1][2] += delz*fpair; //delz = x[i1][2]-x[i2][2]*

}

if (newton_bond || i2 < nlocal) {

f[i2][0] -= delxfpair;

f[i2][1] -= dely*fpair;*

f[i2][2] -= delzfpair;

}

if (evflag) force->pair->ev_tally(i1,i2,nlocal,newton_bond,

evdwl,0.0,fpair,delx,dely,delz);

}

In my own harmonic angle potential source code (“angle_my_harmonic.cpp”), I did the similar to subtracting pairwise potential of atom 1, 3 of an angle:

itype = atom->type[i1];

jtype = atom->type[i3];

delx = x[i1][0] - x[i3][0];

dely = x[i1][1] - x[i3][1];

delz = x[i1][2] - x[i3][2];

domain->minimum_image(delx,dely,delz);

rsq = delx*delx + dely*dely + delz*delz;*

r2inv = sigma[type]*sigma[type]/rsq;*

r6inv = r2invr2invr2inv;

if (rsq < cutsq[itype][jtype])

{

evdwl = -force->pair->single(i1,i3,itype,jtype,rsq,epsilon[type],sigma[type],fpair);

fpair = -fpair;

if (newton_bond || i1 < nlocal)

{

f[i1][0] += delx*fpair;*

f[i1][1] += delyfpair;

f[i1][2] += delz*fpair;*

}

if (newton_bond || i3 < nlocal)

{

f[i3][0] -= delxfpair;

f[i3][1] -= dely*fpair;*

f[i3][2] -= delzfpair;

}

if (evflag)

{

force->pair->ev_tally(i1,i3,nlocal,newton_bond,evdwl,0.0,fpair,delx,dely,delz);

}

}

Using the newly created LAMMPS executable, my system got blown up once MD starts. Then I found that if the pairwise forces subtracted on atom 1,3 were negated (i.e., should be f[i1][*] -= and f[i3][*] +=), the system stabled. However, after some theoretical derivation, the original expression (f[i1][*] += and f[i3][*] -=) seems correct. So I am confused here.

Sorry to bother you again with these specific codes, but this problem has being puzzling me for a while. I really appreciate your time of looking into it.

Best

Shaorui