Special bonds factor_coul in coul/long

Hello,

I’ve come across a counter-intuitive behavior in the coul/long and lj/cut/coul/long pairstyles, and I’m trying to work out whether it’s a bug or a feature I don’t understand. If it’s a feature, I would appreciate some help understanding exactly how it works.

The pairstyles in the KSPACE package based on coul/long seem to handle bonded interactions differently from the corresponding truncated pairstyles in the main src directory.
Compare the following:

lj/cut/coul/cut (line numbers taken from the github commit on 24 May 2021).
(117) if (rsq < cut_coulsq[itype][jtype])
(118) forcecoul = qqrd2e * qtmp*q[j]*sqrt(r2inv);
(119) else forcecoul = 0.0;

(126) fpair = (factor_coul*forcecoul + factor_lj*forcelj) * r2inv;

with lj/cut/coul/long:
(147) prefactor = qqrd2e * qtmp*q[j]/r;
(148) forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
(149) if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;

(171) fpair = (forcecoul + factor_lj*forcelj) * r2inv;

The description of special_bonds in the manual suggests that it should apply equally to all pairstyles, but these two examples are clearly different. In the truncated potential, factor_coul multiplies the entire realspace force but in coul/long it subtracts off only a part of the force. In particular, when factor_coul=0, which is the default for bonded interactions, the coulomb force will be 0 in coul/cut but in lj/cut/coul/long it will be prefactor*(erfc + EWALD_F*grij*expm2 - 1).

I couldn’t find a reason for this difference in the documentation nor is there any explanation in the code comments. I found a reference to it in a non-official guide from 2010 suggesting that factor_coul is used to cancel out that atom pair’s contribution to the kspace long-range force. But I assume that value would depend on the specific method used to evaluate the kspace part (e.g. conventional Ewald vs PPPM), which this doesn’t seem to.

Does anyone know whether the unusual behavior for KSPACE pairstyles is intended? If it is, could you provide any information on what the additional force term is supposed to represent? An equation reference in a paper or link to a manual page would be perfect, but I haven’t been able to find either.

Thanks,
Sam

This is the intended behavior.
The reason for this is the calculation of the coulomb interaction for the /long and equivalent pair styles is split into two parts: the short range part in the pair style and the long-range part in the kspace style. The long range part includes all interactions of all periodic replica of atom i with all periodic replica of atom j, except for the damped short range part, which needs to be computed only for all pairs within the cutoff radius.

Now, if you have an exclusion (i.e. need to apply factor_coul), this applies only to the specific bond/angle/dihedral (or special 1-2, 1-3, or 1-4 pair) and for that you need to exclude the full (not damped) coulomb interaction, and this is what is done in the pair styles. If the code would be structured just like the cutoff code, the part of the coulomb interaction that is “damped away” would not be removed and that would be incorrect.

Hi Axel,

Thanks for your quick response, that makes sense.

I think my confusion arose because it isn’t clear that prefactor is actually the exact coulomb force computed in the minimum-image convention, so I didn’t understand how the quantity being subtracted off related to the original force. Would it be helpful for me to make a pull request to add some comments explaining these lines, or do you prefer to keep the code free of user comments like that?

Thanks again,

Sam

There are too many places where this code construct is used. This would be more something for a chapter in the programmer guide in the manual together with explaining the Special class and related issues.