Dear LAMMPS developers,

For pair styles which require a kspace solver, such as lj/cut/coul/long, the special neighbors 1-2, 1-3, and 1-4 are included in the neighbor list. The short-range term q[i]*q[j]*erfc(alpha*r)/r is computed for all pairs with r <= rc, including these special neighbors. However, whenever factor_coul < 1, a calculation is done so that the correct interaction is recovered afterwards when kspace includes the long-range term. Thus,

prefactor short-range special neighbor kspace

------^------ ------^------ ---------^------- ------^-----

(q[i]*q[j]/r)*[ erfc(alpha*r) - (1 - factor_coul) + erf(alpha*r) ]

This procedure works because erfc(x) + erf(x) = 1. For instance, interaction will fade out if factor_coul = 0.

What I do not understand is why the same procedure is done for lj/cut/coul/dsf and related pair styles, which do not involve a kspace solver and the term prefactor*erf(alpha*r) will not be added afterwards. In the compute() routine of lj/cut/coul/dsf, the line if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor serves to exclude those interactions which should not have been considered initially. However, they are subtracted without performing any shifting/damping. I’m afraid there is an inconsistency there. Am I right?.

Moreover, I wonder if it is correct to apply any shifting/damping to the special-bond coulombic interactions. To the best of my knowledge, none of the force fields commonly employed in molecular dynamics includes shifting or damping in their parameterization procedures. The 1-4 coulombic interaction is an important part of a dihedral model, so that the fitted parameter values often reflect the presence of such interaction.

Thank you,

Ana Silveira

Dear LAMMPS developers,

For pair styles which require a kspace solver, such as lj/cut/coul/long,

the special neighbors 1-2, 1-3, and 1-4 are included in the neighbor list.

The short-range term q[i]*q[j]*erfc(alpha*r)/r is computed for all pairs

with r <= rc, including these special neighbors. However, whenever

factor_coul < 1, a calculation is done so that the correct interaction is

recovered afterwards when kspace includes the long-range term. Thus,

prefactor short-range special neighbor kspace

------^------ ------^------ ---------^-------

------^-----

(q[i]*q[j]/r)*[ erfc(alpha*r) - (1 - factor_coul) + erf(alpha*r) ]

This procedure works because erfc(x) + erf(x) = 1. For instance,

interaction will fade out if factor_coul = 0.

What I do not understand is why the same procedure is done for

lj/cut/coul/dsf and related pair styles, which do not involve a kspace

solver and the term prefactor*erf(alpha*r) will not be added afterwards.

In the compute() routine of lj/cut/coul/dsf, the line if (factor_coul <

1.0) forcecoul -= (1.0-factor_coul)*prefactor serves to exclude those

interactions which should not have been considered initially. However, they

are subtracted without performing any shifting/damping. I’m afraid there is

an inconsistency there. Am I right?.

no. the kspace term is implicit in coul/dsf or coul/wolf. the purpose of

this special treatment for special neighbors is that in those cases you

must exclude the *full* coulomb interaction to be consistent since you

cannot exclude them correctly in the kspace calculation, regardless of

whether you do that one explicitly or implicitly.

Moreover, I wonder if it is correct to apply any shifting/damping to the

special-bond coulombic interactions. To the best of my knowledge, none of

the force fields commonly employed in molecular dynamics includes shifting

or damping in their parameterization procedures. The 1-4 coulombic

interaction is an important part of a dihedral model, so that the fitted

parameter values often reflect the presence of such interaction.

of course not. that is why there is the special treatment in those pair

styles. with the code how it is currently in LAMMPS, you always exclude (or

scale) the same interaction strength, i.e. the undamped coulomb.

axel.