\lammps\lammps-3Oct10\src\DIPOLE\pair_dipole_cut.cpp
Hello,
I have probably found an error within the PairDipoleCut::compute procedure:
If the condition if (rsq < cut_coulsq[itype][jtype]) fails, the forcecoulx , tixcoul are not set to zero, and therefore th old values are added to the forces and torques
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
torque[j][0] += fq*tjxcoul;
torque[j][1] += fq*tjycoul;
torque[j][2] += fq*tjzcoul;
I think the error could be removed by putting the assignments
forcecoulx = forcecouly = forcecoulz = 0.0;
tixcoul = tiycoul = tizcoul = 0.0;
tjxcoul = tjycoul = tjzcoul = 0.0;
outside (above) the if (rsq < cut_coulsq[itype][jtype]) condition.
Here the code with the possible error
void PairDipoleCut::compute(int eflag, int vflag)
{
...
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
...
for (jj = 0; jj < jnum; jj++) {
...
if (rsq < cutsq[itype][jtype]) {
...
// atom can have both a charge and dipole
// i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole
if (rsq < cut_coulsq[itype][jtype]) { // HERE
forcecoulx = forcecouly = forcecoulz = 0.0;
tixcoul = tiycoul = tizcoul = 0.0;
tjxcoul = tjycoul = tjzcoul = 0.0;
... // Calculation of forcecoulx , tixcoul and so on
}
// LJ interaction
if (rsq < cut_ljsq[itype][jtype]) {
...
} else forcelj = 0.0;
// total force
fq = factor_coul*qqrd2e;
fx = fq*forcecoulx + delx*forcelj; // here
fy = fq*forcecouly + dely*forcelj; // here
fz = fq*forcecoulz + delz*forcelj; // here
// force & torque accumulation // here
f[i][0] += fx;
f[i][1] += fy;
f[i][2] += fz;
torque[i][0] += fq*tixcoul;
torque[i][1] += fq*tiycoul;
torque[i][2] += fq*tizcoul;
if (newton_pair || j < nlocal) {
f[j][0] -= fx;
f[j][1] -= fy;
f[j][2] -= fz;
torque[j][0] += fq*tjxcoul;
torque[j][1] += fq*tjycoul;
torque[j][2] += fq*tjzcoul;
}
if (eflag) {
if (rsq < cut_coulsq[itype][jtype]) {
ecoul = qtmp*q[j]*rinv;
if (dipole[itype] > 0.0 && dipole[jtype] > 0.0)
ecoul += r3inv*pdotp - 3.0*r5inv*pidotr*pjdotr;
if (dipole[itype] > 0.0 && q[j] != 0.0)
ecoul += -q[j]*r3inv*pidotr;
if (dipole[jtype] > 0.0 && qtmp != 0.0)
ecoul += qtmp*r3inv*pjdotr;
ecoul *= factor_coul*qqrd2e;
} else ecoul = 0.0;
if (rsq < cut_ljsq[itype][jtype]) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) -
offset[itype][jtype];
evdwl *= factor_lj;
} else evdwl = 0.0;
}
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,ecoul,fx,fy,fz,delx,dely,delz);
}
}
}
if (vflag_fdotr) virial_compute();
}
Here, what should be done, if I am right:
void PairDipoleCut::compute(int eflag, int vflag)
{
...
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
...
for (jj = 0; jj < jnum; jj++) {
...
if (rsq < cutsq[itype][jtype]) {
...
// atom can have both a charge and dipole
// i,j = charge-charge, dipole-dipole, dipole-charge, or charge-dipole
// new position of the assignments to zero
forcecoulx = forcecouly = forcecoulz = 0.0;
tixcoul = tiycoul = tizcoul = 0.0;
tjxcoul = tjycoul = tjzcoul = 0.0;
if (rsq < cut_coulsq[itype][jtype]) {
// old position of the assignments to zero
... // Calculation of forcecoulx , tixcoul and so on
}
...
Please give me feedback, because I am going to modify the dipole-dipole interaction for my own purpose.
Best wishes
Gerald Rosenthal