[gmx-users] Large differences in calculated and experimental relative hydration free energy.
Abhishek Acharya
abhi117acharya at gmail.com
Fri Aug 10 12:32:35 CEST 2018
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
More information about the gromacs.org_gmx-users
mailing list