[gmx-users] Large differences in calculated and experimental relative hydration free energy.

Mark Abraham mark.j.abraham at gmail.com
Sat Aug 11 01:31:47 CEST 2018


Hi,

Sorry for not having time to dig in detail, but does your protocol work for
a simpler molecule?

Mark

On Fri, Aug 10, 2018, 13:33 Abhishek Acharya <abhi117acharya at gmail.com>
wrote:

> Dear Gromacs Users,
>
> I am trying to calculate the relative hydration free energy of Catechol wrt
> Benzene as a way of validating the parameters of Catechol, which I need for
> further free energy calculations. I am using OPLS-AA for these calculations
> and the parameters for both Catechol and Benzene were generated using the
> LigParGen webserver. Subsequently, I used the alchemical_setup.py script (
> https://github.com/MobleyLab/alchemical-setup) to generate the topology
> for
> transformation of catechol to benzene. This was followed by the necessary
> steps to generate a solvated system.
>
> Relative hydration free energy was calculated in 3 steps:
>
> a) Discharging catechol using the following topology setup done (using 11
> lambda values):
>
> ;   nr       type  resnr    res   atom   cgnr     charge       mass
> typeB    chargeB      massB comments  MCSS = Minimum Common Sub-Structure
>      1   opls_800      1    LIG     C1      1    -0.1361    12.0110
> opls_800     0.0    12.0110 ; MCSS
>      2   opls_801      1    LIG     C2      2    -0.1361    12.0110
> opls_801     0.0    12.0110 ; MCSS
>      3   opls_802      1    LIG     C3      3    -0.1693    12.0110
> opls_802     0.0    12.0110 ; MCSS
>      4   opls_803      1    LIG     C4      4     0.2517    12.0110
> opls_803     0.0    12.0110 ; MCSS
>      5   opls_804      1    LIG     C5      5     0.3382    12.0110
> opls_804     0.0    12.0110 ; MCSS
>      6   opls_805      1    LIG     C6      6    -0.1693    12.0110
> opls_805     0.0    12.0110 ; MCSS
>      7   opls_810      1    LIG    HA1      7     0.1508     1.0080
> opls_808     0.0     1.0080 ; MCSS
>      8 opls_810_dummy      1    LIG    HA2      8    0.00000     1.0080
> opls_810     0.0     1.0080 ; to be appeared
>      9   opls_811      1    LIG    HB1      9     0.1508     1.0080
> opls_811     0.0     1.0080 ; MCSS
>     10   opls_812      1    LIG    HC1     10     0.4461     1.0080
> opls_812_dummy    0.00000     1.0080 ; to be annihilated
>     11   opls_813      1    LIG    HD1     11     0.4422     1.0080
> opls_813_dummy    0.00000     1.0080 ; to be annihilated
>     12   opls_806      1    LIG     O1     12    -0.7156    15.9990
> opls_806_dummy    0.00000    15.9990 ; to be annihilated
>     13   opls_807      1    LIG     O2     13     -0.742    15.9990
> opls_807_dummy    0.00000    15.9990 ; to be annihilated
>     14   opls_808      1    LIG     H1     14     0.1443     1.0080
> opls_808     0.0     1.0080 ; MCSS
>     15   opls_809      1    LIG     H2     15     0.1443     1.0080
> opls_809     0.0     1.0080 ; MCSS
>     16 opls_809_dummy      1    LIG     H3     16    0.00000     1.0080
> opls_809     0.0     1.0080 ; to be appeared
>
> Here only coulomb-lambdas were varied as: 0.00 0.10 0.20 0.30 0.40 0.50
> 0.60 0.70 0.80 0.90 1.00
>
> b) LJ-transformation from catechol to benzene (26 lambda values):
>
> ;   nr       type  resnr    res   atom   cgnr     charge       mass
> typeB    chargeB      massB comments  MCSS = Minimum Common Sub-Structure
>      1   opls_800      1    LIG     C1      1     0.0    12.0110
> opls_800     0.0    12.0110 ; MCSS
>      2   opls_801      1    LIG     C2      2     0.0    12.0110
> opls_801     0.0    12.0110 ; MCSS
>      3   opls_802      1    LIG     C3      3     0.0    12.0110
> opls_802     0.0    12.0110 ; MCSS
>      4   opls_803      1    LIG     C4      4     0.0    12.0110
> opls_803     0.0    12.0110 ; MCSS
>      5   opls_804      1    LIG     C5      5     0.0    12.0110
> opls_804     0.0    12.0110 ; MCSS
>      6   opls_805      1    LIG     C6      6     0.0    12.0110
> opls_805     0.0    12.0110 ; MCSS
>      7   opls_808      1    LIG    HA1      7     0.0     1.0080
> opls_808     0.0     1.0080 ; MCSS
>      8 opls_810_dummy      1    LIG    HA2      8    0.0     1.0080
> opls_810     0.0     1.0080 ; to be appeared
>      9   opls_811      1    LIG    HB1      9     0.0     1.0080
> opls_811     0.0     1.0080 ; MCSS
>     10   opls_812      1    LIG    HC1     10     0.0     1.0080
> opls_812_dummy    0.0     1.0080 ; to be annihilated
>     11   opls_813      1    LIG    HD1     11     0.0     1.0080
> opls_813_dummy    0.0     1.0080 ; to be annihilated
>     12   opls_806      1    LIG     O1     12     0.0    15.9990
> opls_806_dummy    0.0    15.9990 ; to be annihilated
>     13   opls_807      1    LIG     O2     13     0.0    15.9990
> opls_807_dummy    0.0    15.9990 ; to be annihilated
>     14   opls_808      1    LIG     H1     14     0.0     1.0080
> opls_808     0.0     1.0080 ; MCSS
>     15   opls_809      1    LIG     H2     15     0.0     1.0080
> opls_809     0.0     1.0080 ; MCSS
>     16 opls_809_dummy      1    LIG     H3     16    0.0     1.0080
> opls_809     0.0     1.0080 ; to be appeared
>
> Here only vdw-lambdas were varied as:
> vdw_lambdas               = 0.00 0.04 0.08 0.12 0.16 0.20 0.24 0.28 0.32
> 0.36 0.40 0.44 0.48 0.52 0.56 0.60 0.64 0.68 0.72 0.76 0.80 0.84 0.88 0.92
> 0.96 1.00
>
> Meanwhile, coulomb_lambdas were kept 0 throughout, as the charges of State
> A and B have been set to zero in this topology.
>
> c) Recharging to full benzene charges (11 lambda values):
>
> ;   nr       type  resnr    res   atom   cgnr     charge       mass
> typeB    chargeB      massB comments  MCSS = Minimum Common Sub-Structure
>      1   opls_800      1    LIG     C1      1     0.0    12.0110
> opls_800    -0.1474    12.0110 ; MCSS
>      2   opls_801      1    LIG     C2      2     0.0    12.0110
> opls_801    -0.1476    12.0110 ; MCSS
>      3   opls_802      1    LIG     C3      3     0.0    12.0110
> opls_802    -0.1474    12.0110 ; MCSS
>      4   opls_803      1    LIG     C4      4     0.0    12.0110
> opls_803    -0.1472    12.0110 ; MCSS
>      5   opls_804      1    LIG     C5      5     0.0    12.0110
> opls_804    -0.1470    12.0110 ; MCSS
>      6   opls_805      1    LIG     C6      6     0.0    12.0110
> opls_805    -0.1470    12.0110 ; MCSS
>      7   opls_808      1    LIG    HA1      7     0.0     1.0080
> opls_808     0.1474     1.0080 ; MCSS
>      8   opls_810      1    LIG    HA2      8     0.0     1.0080
> opls_810     0.1470     1.0080 ; to be appeared
>      9   opls_811      1    LIG    HB1      9     0.0     1.0080
> opls_811     0.1470     1.0080 ; MCSS
>     10 opls_812_dummy      1    LIG    HC1     10     0.0     1.0080
> opls_812_dummy    0.0     1.0080 ; to be annihilated
>     11 opls_813_dummy      1    LIG    HD1     11     0.0     1.0080
> opls_813_dummy    0.0     1.0080 ; to be annihilated
>     12 opls_806_dummy      1    LIG     O1     12     0.0    15.9990
> opls_806_dummy    0.0    15.9990 ; to be annihilated
>     13 opls_807_dummy      1    LIG     O2     13     0.0    15.9990
> opls_807_dummy    0.0    15.9990 ; to be annihilated
>     14   opls_806      1    LIG     H1     14     0.0     1.0080
> opls_806     0.1474     1.0080 ; MCSS
>     15   opls_807      1    LIG     H2     15     0.0     1.0080
> opls_807     0.1474     1.0080 ; MCSS
>     16   opls_809      1    LIG     H3     16     0.0     1.0080
> opls_809     0.1474     1.0080 ; to be appeared
>
> In this case the topology was changed a bit to have the same LJ parameters
> for both State A and B = Benzene.
> Lambda values used for this transformation were:
> coul_lambdas              = 0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80
> 0.90 1.00
> vdw_lambdas               = 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
> 1.00 1.00 (Although, even a value of 0 throughout should work here)
>
> Finally, the free energy change for the these transformations were
> calculated using the alchemical_analysis script (
> https://github.com/MobleyLab/alchemical-analysis).
> Result for the three steps [MBAR output only (in kcal/mol); statistical
> error in brackets]: Discharging: 27.449 (0.028), LJ-transform: - 0.192
> (0.001) and Recharging: - 0.348 (0.039)
>
> The problem is the the reported experimental value for relative hydration
> free energy of catechol and benzene is around 8.5 kcal/mol which is very
> different from the calculated value (27.6 kcal/mol). The simulations have
> been performed using a similar (cutoff, electrostatics etc.) setup as
> mentioned in the LigParGen paper. The default soft-core parameters were
> used for all three steps.
>
> I also tried a one-step transformation (both vdw and coulomb lambda varied
> simultaneously) and got a value of around 29 kcal/mol.
>
> Since, I am attempting a transmutation of this sort for the first time, I
> wanted to be sure the setup is correct before I contact LigParGen authors.
> Since this seems like a simple transformation, I expected differences
> (between experimental and calculated values) of the order of 1.5-2 kcal/mol
> as seen in case of hydration free energies of test cases used in the
> LigParGen paper. It would be nice to get some pointers as to what might be
> going wrong with my setup.
>
> Thanks in advance.
>
> Abhishek Acharya
> --
> 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