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();

}