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