# Water dimer and Free Energy Perturbation (FEP)

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?

• It's generally discouraged to cross post a question in multiple places: chemistry.stackexchange.com/questions/166924/…. You should either close one of the questions or include links to any alternate versions in each post.
– Tyberius
Aug 7, 2022 at 0:39

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...

• can you plz elaborate on the inter-molecular electrostatic interaction part? So to keep all forces active in each monomer, we can't change the atomic charges on each monomer. And electrostatic-energy is $q_i q_j / r_{ij}$. Now $r_{ij}$ comes from geometry and can't be managed through lambda. Only way I think we can kill inter-molecular electrostatic interaction, if the program knows which $q_i$ belongs to which molecule. If they belong to different molecule then interaction is $\lambda_{inter-molec-coulomb} q_i q_j / r_{ij}$. Aug 7, 2022 at 5:14
• The program most definitely does know which charge belongs on each atom. Generally each atom has a type, and the program looks up the LJ sigma and epsilon for that type. For charges, I personally just store a vector of charges, where charge 1 is for atom 1 etc., however, it doesn't matter how this is stored, it is most certainly stored. The program knows everything. The user says which molecule is special, and the program scales its inter-interactions with other molecules Aug 8, 2022 at 12:45
• @user3001408 In short - your comment is correct. Just be aware that most real simulations don't calculate electrostatics based on q1q2/r12, we use some form of ewalds summation to account for periodic boundaries. The mathematical form gets more complicated, but the concept it exactly the same as you have said. Aug 8, 2022 at 12:47