[gmx-users] Bennet error in FEP-calculations for charged ligands.
Artem Shekhovtsov
job.shekhovtsov at gmail.com
Thu Jan 24 13:07:10 CET 2019
Michael, Thank you so much for such a quick and detailed response.
I tried to increase the number of intermediate lambda windows, but
unfortunately, it was not possible to reduce the error to an acceptable
level.
coul-lambdas = 0.0 0.02 0.05 0.08 0.12 0.16 0.2 0.3 0.4 0.6 0.7 0.8 0.84
0.88 0.92 0.95 0.98 1.0
point 0 - 1, DG 47.68 +/- 0.13
point 1 - 2, DG 34.31 +/- 0.56
point 2 - 3, DG 19.42 +/- 0.97
point 3 - 4, DG 15.56 +/- 1.07
point 4 - 5, DG 8.45 +/- 0.03
point 5 - 6, DG 3.53 +/- 0.02
point 6 - 7, DG 4.99 +/- 0.02
point 7 - 8, DG 2.30 +/- 0.02
point 8 - 9, DG -0.02 +/- 0.03
point 9 - 10, DG -2.79 +/- 0.31
point 10 - 11, DG -5.74 +/- 0.72
point 11 - 12, DG -6.58 +/- 0.45
point 12 - 13, DG -8.02 +/- 0.03
point 13 - 14, DG -8.60 +/- 0.03
point 14 - 15, DG -15.92 +/- 0.04
point 15 - 16, DG -33.08 +/- 1.38
point 16 - 17, DG -47.34 +/- 1.22
total 0 - 17, DG 8.15 +/- 1.69
As for the three-level transformation, I agree that it should help, but it
seems I don’t quite understand how to implement this behavior in GROMACS.
I don't know by which options of the .mdp file it should be implemented and
how to build topology files for such transformations.
If you provided me with a small example of an mdp and itp file for such a
transformation, I would be very grateful to you.
Thank you,
Shekhovtsov Artem
On Wed, Jan 23, 2019 at 7:00 PM Michael Shirts <mrshirts at gmail.com> wrote:
> 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.
> >
> --
> 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