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