A bug is found in pair lj/smooth/linear and problem is cracked

Dear friends,
I found that the total potential energy of my system is dependent on number of processors when I used the potential lj/smooth/linear in my simulation. It suggests there is some problem with this potential because the problem disappears when I use this potential by other force field, for example, lj/smooth, morse. I worked hard to see what was wrong and finally I solved the problem by adding only one line in the source code file pair_lj_smooth_linear.cpp:

if (itype > jtype) cut[itype][jtype]=cut[jtype][itype];

This is because cut[itype][jtype] is endowed a value only when itype is smaller than jtype, and it is essentially 0 when this condition is not satisfied, as I observed during debugging. MD will not be affected other than potential energy evaluation, because this is the only occasion cut[][] array should be explicitly specified.

The original code:

void PairLJSmoothLinear::compute(int eflag, int vflag)
{

.
.
.
for (ii = 0; ii < inum; ii++) {
.
.
.
for (jj = 0; jj < jnum; jj++) {
.
.
.
if (rsq < cutsq[itype][jtype]) {
.
.
.
if (eflag) {
.
.
.
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
evdwl = evdwl - ljcut[itype][jtype]

  • (r-cut[itype][jtype])*dljcut[itype][jtype];
    }
    if (evflag) ev_tally(i,j,nlocal,newton_pair,
    evdwl,0.0,fpair,delx,dely,delz);
    .
    .
    .
    }
    }
    }
    if (vflag_fdotr) virial_fdotr_compute();
    }

The modified code:

void PairLJSmoothLinear::compute(int eflag, int vflag)
{

.
.
.
for (ii = 0; ii < inum; ii++) {
.
.
.
for (jj = 0; jj < jnum; jj++) {
.
.
.
if (rsq < cutsq[itype][jtype]) {
.
.
.
if (eflag) {
.
.
.
if (itype > jtype) cut[itype][jtype]=cut[jtype][itype];

evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
evdwl = evdwl - ljcut[itype][jtype]

  • (r-cut[itype][jtype])*dljcut[itype][jtype];
    }
    if (evflag) ev_tally(i,j,nlocal,newton_pair,
    evdwl,0.0,fpair,delx,dely,delz);
    .
    .
    .
    }
    }
    }
    if (vflag_fdotr) virial_fdotr_compute();
    }

hi steffani,

Dear friends,
       I found that the total potential energy of my system is dependent on
number of processors when I used the potential lj/smooth/linear in my
simulation. It suggests there is some problem with this potential because
the problem disappears when I use this potential by other force field, for
example, lj/smooth, morse. I worked hard to see what was wrong and finally I
solved the problem by adding only one line in the source code file
pair_lj_smooth_linear.cpp:

if (itype > jtype) cut[itype][jtype]=cut[jtype][itype];

this is not a good correct. a better is to not change
the compute() method but rather init_one() method.
e.g. like this:

diff --git a/src/pair_lj_smooth_linear.cpp b/src/pair_lj_smooth_linear.cpp
index 7e8a8ca..d77c8c7 100644
--- a/src/pair_lj_smooth_linear.cpp
+++ b/src/pair_lj_smooth_linear.cpp
@@ -247,6 +247,7 @@ double PairLJSmoothLinear::init_one(int i, int j)
   lj4[j][i] = lj4[i][j];
   ljcut[j][i] = ljcut[i][j];
   dljcut[j][i] = dljcut[i][j];
+ cut[j][i] = cut[i][j];

   return cut[i][j];
}

cheers,
     axel.

Dear Dr. Kohlmeyer,
       thank you for pointing it out. Whatever functions better will be helpful.