[gmx-users] Bennet error in FEP-calculations for charged ligands.
Michael Shirts
mrshirts at gmail.com
Wed Jan 23 17:00:15 CET 2019
As free energies get larger, then the error gets less accurate. So if it is
reporting 18.51 +/- 2.95 for one of the intervals, that likely suggest
there is very little overlap in that area.
A total free energy difference of in -6.12 +/- 3.19 for the
transformation indicates that the result is not very certain; you're within
two standard deviations, and again, the Bennett error formula is inaccurate
for larger error. I would suggest adding 1-2 more intermediates between
those points.
https://github.com/MobleyLab/alchemical-analysis provides some of the tools
to check overlap.
I would also suggest doing a multi stage transformation to increase overlap
- for atoms that are disappearing or appearing that are charged, turn off
charges off, change vdw, then turn charges back on. In some cases,
especially when the atoms involved have different sizes, sc-coul can lead
to some semi-pathological cases when vdw softcore and coul softcore result
in some very low/high potentials at intermediates.
On Wed, Jan 23, 2019 at 8:21 AM Artem Shekhovtsov <job.shekhovtsov at gmail.com>
wrote:
> Hi all!
> I encountered an error in my relative free energy calculations and do not
> know how to fix it.
> Molecules for which I want to carry out calculations contain a carboxyl
> group.
> To validate the protocol, I tried to run the fep-calculations of symmetric
> molecules for which the change in free energy will be zero.
> During validation, I was faced with the fact that the convergence error for
> charged ligands greatly exceeds that for neutral ones.
> For example, for the case of m-methyl benzoic acid [O-]C(= O)c1cc(C)ccc1:
> Mutation (-CH3 + H) for one meta position (-H + CH3) for another meta
> position relative to the carboxyl group.
>
> Ionized acid (solvent leg, 5 ns):
> point 0 - 1, DG 28.74 +/- 0.02
> point 1 - 2, DG 53.79 +/- 0.41
> point 2 - 3, DG 27.76 +/- 2.45
> point 3 - 4, DG 18.51 +/- 2.95
> point 4 - 5, DG 8.80 +/- 0.79
> point 5 - 6, DG 4.04 +/- 0.07
> point 6 - 7, DG -0.06 +/- 0.04
> point 7 - 8, DG -2.62 +/- 0.33
> point 8 - 9, DG -8.95 +/- 0.27
> point 9 - 10, DG -23.89 +/- 1.61
> point 10 - 11, DG -29.32 +/- 1.29
> point 11 - 12, DG -54.15 +/- 0.08
> point 12 - 13, DG -28.79 +/- 0.03
>
> total 0 - 13, DG -6.12 +/- 3.19
>
> Unionized acid (solvent leg, 5 ns):
> point 0 - 1, DG 0.08 +/- 0.01
> point 1 - 2, DG -6.88 +/- 0.10
> point 2 - 3, DG -14.69 +/- 0.05
> point 3 - 4, DG -21.07 +/- 0.03
> point 4 - 5, DG -13.21 +/- 0.02
> point 5 - 6, DG -7.46 +/- 0.02
> point 6 - 7, DG -0.09 +/- 0.03
> point 7 - 8, DG 7.37 +/- 0.02
> point 8 - 9, DG 13.24 +/- 0.03
> point 9 - 10, DG 21.15 +/- 0.06
> point 10 - 11, DG 14.72 +/- 0.03
> point 11 - 12, DG 6.98 +/- 0.07
> point 12 - 13, DG -0.05 +/- 0.01
>
> total 0 - 13, DG 0.09 +/- 0.22
>
> For a neutral molecule containing a charged carboxyl and amino groups
> ([O-]C(=O)c1cc(C)c([NH3+])cc1) the error is small:
> point 0 - 1, DG -4.65 +/- 0.01
> point 1 - 2, DG -19.95 +/- 0.04
> point 2 - 3, DG -28.67 +/- 0.13
> point 3 - 4, DG -57.84 +/- 0.19
> point 4 - 5, DG -44.18 +/- 0.08
> point 5 - 6, DG -26.84 +/- 0.05
> point 6 - 7, DG 0.10 +/- 0.20
> point 7 - 8, DG 26.92 +/- 0.06
> point 8 - 9, DG 44.28 +/- 0.13
> point 9 - 10, DG 57.79 +/- 0.13
> point 10 - 11, DG 28.68 +/- 0.09
> point 11 - 12, DG 19.98 +/- 0.12
> point 12 - 13, DG 4.66 +/- 0.03
>
> total 0 - 13, DG 0.28 +/- 0.21
>
> Adding Na and Cl ions to ([O-]C(=O)c1cc(C)c([NH3+])cc1) does not cause an
> increase in error.
> point 0 - 1, DG -4.61 +/- 0.01
> point 1 - 2, DG -19.96 +/- 0.07
> point 2 - 3, DG -28.68 +/- 0.04
> point 3 - 4, DG -57.67 +/- 0.04
> point 4 - 5, DG -44.24 +/- 0.01
> point 5 - 6, DG -26.91 +/- 0.04
> point 6 - 7, DG 0.03 +/- 0.03
> point 7 - 8, DG 26.89 +/- 0.03
> point 8 - 9, DG 44.20 +/- 0.08
> point 9 - 10, DG 57.75 +/- 0.09
> point 10 - 11, DG 28.65 +/- 0.06
> point 11 - 12, DG 20.02 +/- 0.05
> point 12 - 13, DG 4.67 +/- 0.01
>
> total 0 - 13, DG 0.12 +/- 0.18
>
> FEP-part of *.mdp:
> free-energy = yes
> sc-power = 1
> sc-alpha = 0.5
> sc-coul = yes
> fep-lambdas = 0.0 0.01 0.05 0.1 0.2 0.3 0.4 0.6 0.7
> 0.8 0.9 0.95 0.99 1.0
> coul-lambdas = 0.0 0.01 0.05 0.1 0.2 0.3 0.4 0.6 0.7
> 0.8 0.9 0.95 0.99 1.0
> vdw-lambdas = 0.0 0.01 0.05 0.1 0.2 0.3 0.4 0.6 0.7
> 0.8 0.9 0.95 0.99 1.0
>
> What is the reason of such error and how I can fix it?
>
> mdp, itp, gro, xvg files by link -
> https://drive.google.com/open?id=1MiepOQb2QAZ9rclpY13owCnl8QWcTHnf
> Ready to provide any additional information.
>
> Thank you,
> Shekhovtsov Artem
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list