[gmx-users] Free Energy Perturbation plot: symmetry breaking?

[Matteo De Chiara]

Hello everybody,

I 'd like to calculate the Delta G of a Protein-Ligand complex. More
precisely, what I'm trying to evaluate is the Delta-Delta G of a group
of different ligands, in order to investigate their relative binding
affinity.

I performed 100 simulation using a lamba step of 0.01. Then, for each
simulatuion, I run two mdrun using the options -rerun file.trr  and
using in the first run the tpr file derived from previous simulation
and in the second run the tpr file derived from the following
simulation ( i.e. mdrun -rerun 1.trr -s 2.tpr and mdrun -rerun 1.trr
-s 0.tpr). I use g_energy with -fee option to evaluate the Delta G
using the edr file derived from the starting simulation and the edr
file derived from the re-run of the simulation respect the previous
step, as described before. I repeted this evaluation using even the
edr file derived from the re-run of the mdrun using the following
step. ( i.e.  1.edr and 1_rerun2.edr, 1.edr and 1_rerun0.edr).
Finally, I added the first Delta G with the second, then, to the
resulting sum, I added the third delta G, and so on. I plotted these
Delta G-partial-sum against the relative step. In order to validate
the analysis, I plotted two different set of point: for the forward
and backward calculation.

I was expecting to obtain a parabola-like curve, but what I obtained
is a Mexican-hat-like plot (something like the graph below).


      **
    *    *
  *         *
*              *                   ***
                  *              *      *
                      *        *          *
                         * * *              *
                                               *
                                                 *

This should be the typical plot arising from a spontaneous
symmetry-breaking occurs. Is this interpretation correct? If it is,
can someone explain me how and why this happened?
In case, could this be a source of errors or artefacts in my analysis?
And, in case, can someone suggest me how to solve my problem?

Thanks a lot.

Cheers,

--

Matteo De Chiara
Ph.D. Student
Dept. of Molecular Biology, University of Siena
Via Fiorentina 1,
53100 Siena, Italy