Water dimer and Free Energy Perturbation (FEP)

Question

I am trying to understand FEP in the context of a water dimer. Let's say I want to calculate $\Delta A$, the free energy change for the following reaction:

$$\ce{H_2O +H_2O\rightarrow H_2O-H_2O},\tag{1}$$

in which the right hand side denotes a dimerized water molecule. Let's call the left side state 1 and right side state 2. In the left side there is no H-bond.

Now as per Zwanzig, the FEP equation is:

$$ \Delta A=kT\ln \left\langle e^{(-[E_2 - E_1]/kT)}\right\rangle_1,\tag{2} $$

in which $E_2$ is the energy of a water dimer, $E_1$ is the energy of two non-interacting molecules, and the sampling is done in two non-interacting water molecule states.

Now the total number of particles is 6 (four H atoms, and 2 O atoms). Now when I am sampling at state 1, the conformational coordinate of say the oxygen atom (in the second water molecule), is also the coordinate of the second O atom in the dimer, because I am passing the same coordinate space to the dimer as well. Here is where I am struggling to understand: The coordinates of 6 atoms are the same in both state 1 and state 2, but the only difference is the H-bond which will come into effect in the dimer. So in the $\lambda$ context, lets say for simplicity we have one H-bond between H atom #3 (of the first water molecule) and O atom #5 (of the second water molecule). Then we can define the Lennard-Jones potential between atom #3 and atom #5 as,

$$\sigma_ {(3,5)}(\lambda)=\lambda \sigma_ {(3,5)},\tag{3}$$ in which $\lambda=0$ denotes state 1, and $\lambda=1$ denotes state 2. This way we will have full blown Lennard-Jones in state 2. Is this how it is modeled in this particular case?

My second question is on electrostatic interaction, say at particular lambda (say $\lambda =0$), we have to run MD on state 1, which at this value of lambda represents non-interacting state. But when we run two such water molecules in MD, the electrostatic interaction will take place nevertheless between these water molecules. I mean I can parametrize away the Lennard-Jones, but it seems I can't kill the inter-molecular electrostatic interaction, as I don't seem to have parameter which I can set at zero at $\lambda=0$. Hence it seems to me I can't satisfy the non-interacting criteria in state 1. Where am I going wrong on electrostatic part?

Answer 1

You are on the right track, but there are a few more details, typically. For one, for a FEP, we calculate the change over many windows. We generally also used soft-core potentials. Soft-core potentials avoid numerical problems at the endpoints when doing perturbations. Also, we don't calculate the energy of a dimer, and subtract the energy of two monomers. You can do this, but it wouldn't be free energy perturbation.

Free energy perturbation, in this case, I believe only needs to sample the differences in potential energy as two monomers go from not interacting, to fully interacting (or vice versa).

Usually we turn off the electrostatic interactions between the two monomers first, typically in 10 windows(small changes provide more accurate results).

NOTE: we keep all interactions alive in each monomer. It is only the interactions between monomers we turn off gradually. This tells us the free energy difference between two monomers fully interacting (a dimer) and two monomers only interacting via LJ potential energy functions.

We then turn off the LJ interaction between the two monomers over (often) about 20 windows. We use soft-core potentials for the LJ interactions between the two monomers to avoid the end-point-catastrophe.

This second part tells us the difference in free energy between two monomers interacting with each other only with LJ potentials, and two monomers not interacting at all with each other.

We can combine the electrostatic decoupling free energy, with the LJ decoupling free energy to get the total free energy difference between a dimer, and two monomers not interacting.

We could, rather than decoupling i.e., going from a dimer to two monomers, start at two monomers, and gradually turn on LJ potential energy interactions between the monomers, and then gradually turn on the electrostatic interactions. All that changes is the sign of the free energy change.

Also, we can also decouple the electrostatic and LJ interactions simultaneously. Infact, there are a lot of games that can be played.

Finally, there are many different Lambda parameters. You do not scale a single one.

See the Gromacs manual and online support (for example) for more details, but here is a screenshot to help illustrate:

We can perturb whatever we want so long as we do it correctly...