Coupling two ions to lambda=0.5 cuts their charge-charge interaction force by 1/2 not 1/4

User discussions
1

GROMACS version: 2023
GROMACS modification: No

Hello,

I have a question about charge coupling to lambda. In a system with two ions of opposite charge that are both coupled to lambda, we have seen that these two ions seek one another out and rapidly form direct interactions at intermediate lambda values (though not at lambda=0 or lambda=1), in a box of water with other non-lambda-couple-ions. As part of this investigation, I realized that two ions whose charges are both half-decoupled do not interact with one another in the same way that two non-lambda-coupled ions with half-charges do. I am unable to logically connect what I have found with these force calculations to to the observed formations of ion pairs, but at the moment this is my only handle for what appears to be pathological behavior in the longer simulations where decoupled ions rapidly associate.

My two questions are:
(a) is this somehow expected?
(b) is there a workaround to get two lambda-coupled ions to feel each other’s charges as they would feel non-lambda-coupled ions with the given effective charge?

Worked example:
Take two atoms, LJ epsilon = 0 and put them 0.5 nm apart in a cube with edge length 10 nm. Put the separation vector on Cartesian X to easily get the force vector.

Assign atom 1 a charge of +1 and assign atom 2 a charge of -1 and don’t use the free energy code.
Forces (kJ/mol/nm) are:
Atom 1: 558.373
Atom 2: -558.373

Make the charges +/- 0.5.
Forces (kJ/mol/nm) are:
Atom 1: 139.593
Atom 2: -139.593

Make the charge definitions +/- 1.0 e but make both B-state charges 0.0 and set lambda = 0.5.
Forces (kJ/mol/nm) are:
Atom 1: 279.186
Atom 2: -279.186

In fact, the forces don’t fall to +/-139.593 until lambda = 0.75

Now combine: Make the anion a static -0.5 e charge not dependent on lambda and make the cation +1.0 e in the A state and 0 e in the B state and set lambda = 0.5.
Forces (kJ/mol/nm) are:
Atom 1: 140.456
Atom 2: -139.593

Or the combination above when the anion has a static -1 e charge.
Forces (kJ/mol/nm) are:
Atom 1: 279.186
Atom 2: -280.913

I am guessing that the issue is that with integer state-A charges and 0 state-B charges and lambda=0.5, (1-lambda)q_iq_j is not equal to (1-lambda)q_i(1-lambda)*q_j . Stated another way, I am guessing that each lambda-scaled charge is not aware when interacting partners also have a lambda-scaled charge.

We can reproduce these two ions coming together in a box of water within 0.1 ns whereas other non-lambda-coupled ions tend to remain separated on that timescale (gromacs 2023, 2019.6, 2016.6, 5.1.6), but for force/energy simplicity what I worked up above is for ions in constant volume simulations without water.

Runtime options:

integrator = md
dt = 0.002
nsteps = 0
comm-mode = Linear
nstcomm = 100
nstxout = 1
nstvout = 1
nstfout = 1
nstlog = 1
nstenergy = 1
nstxout-compressed = 1
cutoff-scheme = verlet
nstlist = 10
pbc = xyz
rlist = 0.8
coulombtype = PME
rcoulomb = 0.8
vdw-type = cut-off
rvdw = 0.8
fourierspacing = 0.12
pme_order = 4
tcoupl = no
gen-vel = no
continuation = no
constraints = none

free-energy = yes
init-lambda = 0.5
delta-lambda = 0
sc-function = gapsys
; Gapsys Soft Core Parameters
sc-gapsys-scale-linpoint-lj = 0.85
sc-gapsys-scale-linpoint-q = 0.1
sc-gapsys-sigma-lj = 0.3

Thank you for your advice,
Chris.

2

PS: seems possible that the effect we’re observing (seemingly unphysical attraction between lambda-coupled Na and Cl in water at lambda=0.5) is an indirect effect based on (a) the fact that the lambda-scaled real-space Coulombic term is related to (1-lambda) * q_i * q_j not (1-lambda) * q_i * (1-lambda) * q_j such that it doesn’t matter whether one or both charges are coupled to lambda, and (b) the instantaneous dielectric between two ions is expected to be reduced by lambda coupling of either (with even more reduction if both are coupled) because lambda-scaled groups have less charge density so organize their local environments less so have a smaller local, instantaneous dielectric around them. This would lead to 3 cases: (a) two uncoupled ions with high instantaneous dielectric and full Coulombic attraction, (b) only one of the ions is coupled to lambda so the Coulombic attraction is cut in half and the instantaneous dielectric is reduced around one of the two ions (call it a medium-strength dielectric), and (c) both ions coupled to lambda so the Coulombic attraction is still only cut in half but now the instantaneous dielectric is reduced around both ions. The leap of faith here is that this case C leads to the bizarre attractive motion we’re seeing in our simulations.