[gmx-users] Where to stop with TM protein PMF calculations

[Nash, Anthony]

Good morning,

A bit of a long one I am afraid. 

I am simulating a transmembrane dimer, and calculating the association free energy through potential of mean force calculations as a function of interhelical distance. I have got very good umbrella coverage along my reaction coordinate, however, I would like to know where I should stop calculating and normalise to zero? I am comparing two systems of identical composition but with a different conformation. Hence, the normalising step is vital. 

Looking in the literature the cut offs used seem to be around 2nm (two JACS papers come to mind). I've pulled as far as needed to reach a plateau in the PMF curve. This is around 6.5 nm to 8 nm. The problem is if I normalise at the plateau (which seems to be the beliefs of the professors at my institute), the comparative free energy between the two systems is very different to normalising at 2nm (note: the papers in literature do not observe a plateau, there is still an obvious upward trend at their cutoff).

I need justification of where to normalise. According to literature they cut off outside of interaction range. I have decided to test umbrella windows along the reaction coordinate by decomposing the energy of the system. LJ short and Coul short between peptides at a reaction coordinate distance of 8 nm, is 0 as expected. However, I am still getting energy values on Coul. recip even when I decompose (although the gromacs literature says I can't) the system by systematically setting every charge to zero, i.e.,


Peptide A (A), Peptide B (B), SOL, lipid and counter ions have zero charge throughout. 
E_coul_recip_A_B = E_coul_recip_(A_B,A_A,B_B) - (E_coul_recip_A_A + E_coul_recip_B_B)

All charges are zero, with the exception to peptide A in E_coul_recip_A_A and peptide B in E_coul_recip_B_B. and peptide A and B in E_coul_recip_(A_B,A_A,B_B). When I add E_coul_recip_A_B to the LJ short and Coul short I do not get an energy of zero.

So:
1) The coul-recip decomposition, is this a flawed approach? The gromacs manual says it can't be decomposed.
2) Where along the reaction coordinate should I cut off?
3) Why don't authors pull until they reach a plateau?

Many many thanks
Anthony