[lammps-users] Possible error in pair_dipole_cut.cpp

\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

hello gerald,

good catch. i agree with your findings.
whether the error will show up or not
will depend a lot on the compiler used
and particularly whether somebody compiles
with high optimization or not. but this is
definitely a bug and your bugfix seems
to be the correct solution.

cheers,
    axel.

yes - that's a bug - I'll post a patch.

The coulomb/dipole cutoff is typically >= the
LJ cutoff, in which case it wouldn't matter,
but it is a bug.

Thanks,
Steve