modify bond_quartic.cpp to enable useage with angle potential

Dear LAMMPS users:

I am modifying some LAMMPS source codes to allow the quartic style bond to be used with angle potential. Basically I am trying to add lines to the file “bond_quartic.cpp” to do the following: once a bond is broken, modify the type of the angle (associated with the broken bond) to zero. I have a few questions though:

1). I only need to modify the permanent angle type (atom->angle_type[][]), the temporary angle type (neighbor->anglelist[][3]) does not need modifications here. Am I right?
2). Once the angle type is set to zero in the bond routine, will the angle routines (e.g. angle_harmonic.cpp) automatically ignore this type of angles while calculating energy?
3). How can I create a new bond style, e.g. “my_quartic”?

Thanks and Happy NY!

Shaorui

Comments below.

Steve

Dear LAMMPS users:

I am modifying some LAMMPS source codes to allow the quartic style bond to
be used with angle potential. Basically I am trying to add lines to the file
"bond_quartic.cpp" to do the following: once a bond is broken, modify the
type of the angle (associated with the broken bond) to zero. I have a few
questions though:

1). I only need to modify the permanent angle type (atom->angle_type),
the temporary angle type (neighbor->anglelist[3]) does not need
modifications here. Am I right?

Don't think so. The temporary list is used until the next reeighboring.
If you want the angle to be turned off "now", then you'd have to flag
it off. And all angle potentials would have to check that flag. There
should probably be a unique angle_quartic potential to use with this
model, else you have to add code to other angle potentials, which is
not a good idea.

You also need to make sure that all the pairwise interacions are
corrrect both with and without a broken angle. Bond quartic does
this for the broken bond via the "sleight-of-hand" discussed on
its doc page. Something similar needs to be done correctly for
breaking angles, e.g. in an angle_quartic style.

2). Once the angle type is set to zero in the bond routine, will the angle
routines (e.g. angle_harmonic.cpp) automatically ignore this type of angles
while calculating energy?

No, see above.

3). How can I create a new bond style, e.g. "my_quartic"?

The same way any other bond style is added to LAMMPS. You
can derive it from class Bond or from class BondQuartic. See
Section_modify of the manual for more details.

Good luck - this is not an easy task to get right.

Hi Steve:

Thank you very much for your advice. In the bond routine there are a few lines modifying the permanent bond type:

for (m = 0; m < atom->num_bond[i1]; m++)
if (atom->bond_atom[i1][m] == atom->tag[i2])
atom->bond_type[i1][m] = 0;
if (i2 < atom->nlocal)
for (m = 0; m < atom->num_bond[i2]; m++)
if (atom->bond_atom[i2][m] == atom->tag[i1])
atom->bond_type[i2][m] = 0;

My question is why i1 is not checked for being less than atom->nlocal while i2 is checked.

I added my codes changing angle types right after the above lines:

/* --------------ADDED---------------------------------------------*/
for (an = 0; an < nanglelist; an++)
{
if (anglelist[an][3] <= 0) continue;

ai1 = anglelist[an][0];
ai2 = anglelist[an][1];
ai3 = anglelist[an][2];

if (((atom->tag[ai1] == atom->tag[i1])&&(atom->tag[ai2] == atom->tag[i2]))||((atom->tag[ai1] == atom->tag[i2])&&(atom->tag[ai2] == atom->tag[i1])))
{
anglelist[an][3] = 0;

if (ai1 < atom->nlocal)
{
for (am = 0; am < atom->num_angle[ai1]; am++)
{
if (atom->angle_atom2[ai1][am] == atom->tag[ai2])
{
atom->angle_type[ai1][am] = 0;
}
}
}
if (ai2 < atom->nlocal)
{
for (am = 0; am < atom->num_angle[ai2]; am++)
{
if (atom->angle_atom1[ai2][am] == atom->tag[ai1])
{
atom->angle_type[ai2][am] = 0;
}
}
}
if (ai3 < atom->nlocal)
{
for (am = 0; am < atom->num_angle[ai3]; am++)
{
if ((atom->angle_atom1[ai3][am] == atom->tag[ai1])&&(atom->angle_atom2[ai3][am] == atom->tag[ai2]))
{
atom->angle_type[ai3][am] = 0;
}
}
}
}

if (((atom->tag[ai2] == atom->tag[i1])&&(atom->tag[ai3] == atom->tag[i2]))||((atom->tag[ai2] == atom->tag[i2])&&(atom->tag[ai3] == atom->tag[i1])))
{
anglelist[an][3]=0;

if (ai1 < atom->nlocal)
{
for (am = 0; am < atom->num_angle[ai1]; am ++)
{
if ((atom->angle_atom2[ai1][am] == atom->tag[ai2])&&(atom->angle_atom3[ai1][am] == atom->tag[ai3]))
{
atom->angle_type[ai1][am] = 0
}
}
}
if (ai2 < atom->nlocal)
{
for (am = 0; am < atom->num_angle[ai2]; am ++)
{
if (atom->angle_atom3[ai2][am] == atom->tag[ai3])
{
atom->angle_type[ai2][am] = 0;
}
}
}
if (ai3 < atom->nlocal)
{
for (am = 0; am < atom->num_angle[ai3]; am ++)
{
if (atom->angle_atom2[ai3][am] == atom->tag[ai2])
{
atom->angle_type[ai3][am] = 0;
}
}
}
}
}
/* --------------ADDED---------------------------------------------*/

Does this look correct? Anything I missed? Thanks.

Shaorui

i1 is an atom owned by this processor, so it
is always less than nlocal.

Re: the details of what you've done - sorry,
but I don't have time to debug other people's
code. My best suggestion is to understand
the data structures LAMMPS uses to store
angles, both permanently and temporarily,
and update them appropriately, taking care
of parallel issues as well. If you have
a conceptual question, I can try to answer it.

Steve

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] += delxfpair; //delx = x[i1][0]-x[i2][0]
f[i1][1] += dely
fpair; //dely = x[i1][1]-x[i2][1]
f[i1][2] += delzfpair; //delz = x[i1][2]-x[i2][2]
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= delx
fpair;
f[i2][1] -= delyfpair;
f[i2][2] -= delz
fpair;
}

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 = delxdelx + delydely + delzdelz;
r2inv = sigma[type]sigma[type]/rsq;
r6inv = r2inv
r2inv
r2inv;

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] += delxfpair;
f[i1][1] += dely
fpair;
f[i1][2] += delzfpair;
}
if (newton_bond || i3 < nlocal)
{
f[i3][0] -= delx
fpair;
f[i3][1] -= delyfpair;
f[i3][2] -= delz
fpair;
}

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

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\-&gt;single\(\)
// required since special\_bond = 1,1,1
// tally energy/virial in pair, using newton\_bond as newton flag

itype = atom\-&gt;type\[i1\];
jtype = atom\-&gt;type\[i2\];

if \(rsq &lt; cutsq\[itype\]\[jtype\]\) \{
  evdwl = \-force\-&gt;pair\-&gt;single\(i1,i2,itype,jtype,rsq,1\.0,1\.0,fpair\);
  fpair = \-fpair;

  if \(newton\_bond || i1 &lt; nlocal\) \{
    f\[i1\]\[0\] \+= delx\*fpair;                           //delx =

x[i1][0]-x[i2][0]
f[i1][1] += dely*fpair; //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] -= delx*fpair;
f[i2][1] -= dely*fpair;
f[i2][2] -= delz*fpair;
}

  if \(evflag\) force\-&gt;pair\-&gt;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\-&gt;minimum\_image\(delx,dely,delz\);
rsq = delx\*delx \+ dely\*dely \+ delz\*delz;
r2inv = sigma\[type\]\*sigma\[type\]/rsq;
r6inv = r2inv\*r2inv\*r2inv;

if \(rsq &lt; cutsq\[itype\]\[jtype\]\)
\{
    evdwl =

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

this call to pair::single() is bogus. i am surprised you
don't get a segmentation fault. where you have epsilon[type]
and sigma[type] should be 1.0 and 1.0. those are the
arguments for the lj and coul exclusion scaling factors,
which but are required to be to 1.0 in this case. why
did you change these. it is not done like this in bond_quartic.

    fpair = \-fpair;

    if \(newton\_bond || i1 &lt; nlocal\)
    \{
        f\[i1\]\[0\] \+= delx\*fpair;
        f\[i1\]\[1\] \+= dely\*fpair;
        f\[i1\]\[2\] \+= delz\*fpair;
    \}
    if \(newton\_bond || i3 &lt; nlocal\)
    \{
        f\[i3\]\[0\] \-= delx\*fpair;
        f\[i3\]\[1\] \-= dely\*fpair;
        f\[i3\]\[2\] \-= delz\*fpair;
    \}

    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.

this is a case where you should generate a system with a small number
of atoms and print out all forces to exactly debug all properties.

axel.