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-requ...@gromacs.org.

Reply via email to