[gmx-users] question of perturbing order in single topology FEP, gromacs
Peng He
hepeng.yanglu at gmail.com
Wed Nov 29 20:56:41 CET 2017
Hi, Asaf,
Thank you very much for your suggestions, I am not sure I understand it
correctly.
2017-11-29 13:22 GMT-05:00 <asaffarhi at post.tau.ac.il>:
> Dear Peng,
>
> You can transform each of the molecules seperately.
You can remove the nonbonded interactions of the atoms that are different
> between the compared molecules.
>
Does it mean when going from A-B to A-C, decouple the B's nonbonded than
coupling the C's nonbonded interaction?
> Then, if you have only one dihedral relating the different atoms to the
> common atoms in each molecule it does not need to be removed. Also, the
> bond angle and covalent bond terms do not need to be removed. This is based
> on NVT but may be also accurate for NPT.
>
In the Ethane -> methanol case, the C-C bond length decreases from 1.5 to
1.4 when perturbing to C-O bond, and the kappa for the bond also changes, I
am not sure but I think we can use bonded lambdas to perturb the bond, I
did not remove any bond or angle in this case.
I think what you suggest is dual topology method which contains a dummy CH3
group in Methanol and a dummy OH group in Ethane, in that case I don't need
to change the bond length and angle, and it would not create such different
free energy values. I really test it in dual topology method and there is
little effect from the order of perturbations.
>
> We also saw that the order of transforming matters but maybe you can try
> without removing the bonded terms.
> For more details you can read:
> http://www.sciencedirect.com/science/article/pii/S0010465516303411
>
> Best,
> Asaf
>
>
> Quoting Peng He <hepeng.yanglu at gmail.com>:
>
> HI, all,
>>
>> I am recently doing some test FEP calculations to get familiar, I am
>> preparing a system of ethane -> methanol using Amber gaff force field and
>> single topology method.
>>
>> " 4 [ atomtypes ]
>> 5 ;name bond_type mass charge ptype sigma epsilon
>> Amb
>> 6 c3 c3 0.00000 0.00000 A 3.39967e-01
>> 4.57730e-01 ; 1.91 0.1094
>> 7 hc hc 0.00000 0.00000 A 2.64953e-01
>> 6.56888e-02 ; 1.49 0.0157
>> 8 oh oh 0.00000 0.00000 A 3.06647e-01
>> 8.80314e-01 ; 1.72 0.2104
>> 9 h1 h1 0.00000 0.00000 A 2.47135e-01
>> 6.56888e-02 ; 1.39 0.0157
>> 10 dh hc 0.00000 0.00000 A 0.00000 0.00000
>> 11 ho ho 0.00000 0.00000 A 0.00000e+00
>> 0.00000e+00
>> 12 [ moleculetype ]
>> 13 ;name nrexcl
>> 14 UNK 3
>> 15
>> 16 [ atoms ]
>> 17 ; nr type resi res atom cgnr charge mass ; qtot
>> bond_type
>> 18 1 c3 1 UNK C1 1 -0.095100 12.01000 c3
>> 0.116700 12.01000
>> 19 2 c3 1 UNK C2 2 -0.095100 12.01000 oh
>> -0.598801 12.01000
>> 20 3 hc 1 UNK H1 3 0.031700 1.00800 h1
>> 0.028700 1.00800
>> 21 4 hc 1 UNK H2 4 0.031700 1.00800 h1
>> 0.028700 1.00800
>> 22 5 hc 1 UNK H3 5 0.031700 1.00800 h1
>> 0.028700 1.00800
>> 23 6 hc 1 UNK H4 6 0.031700 1.00800 ho
>> 0.396000 1.00800
>> 24 7 hc 1 UNK H5 7 0.031700 1.00800 dh
>> 0.000000 1.00800
>> 25 8 hc 1 UNK H6 8 0.031700 1.00800 dh
>> 0.000000 1.00800
>> 26
>> 27 [ bonds ]
>> 28 ; ai aj funct r k
>> 29 1 2 1 1.5375e-01 2.5179e+05 1.4233e-01
>> 2.6501e+05
>> 30 1 3 1 1.0969e-01 2.7665e+05 1.0969e-01
>> 2.7665e+05
>> 31 1 4 1 1.0969e-01 2.7665e+05 1.0969e-01
>> 2.7665e+05
>> 32 1 5 1 1.0969e-01 2.7665e+05 1.0969e-01
>> 2.7665e+05
>> 33 2 6 1 1.0969e-01 2.7665e+05 9.7300e-02
>> 3.1079e+05
>> 34 2 7 1 1.0969e-01 2.7665e+05 1.0969e-01
>> 2.7665e+05
>> 35 2 8 1 1.0969e-01 2.7665e+05 1.0969e-01
>> 2.7665e+05
>> ...
>> ..."
>>
>> I know to perturbing from large molecule(ethane) to small(methanol), I
>> need
>> to perturb the electrostatic interaction first than vdw, But I don't know
>> where to perturb the bonded lambdas. I have tested different orders of
>> Bonded lambdas, and I have got a quite different answer which seems odd to
>> me. I have a table for that
>> stage1 stage2 stage3 dG(kj) dG(kcal)
>> single topology solv q/bonded/vdw -14.36 -3.43
>> bonded q vdw -13.79 -3.30
>> bonded/q vdw -14.91 -3.56
>> q bonded vdw -17.13 -4.09
>> q bonded/vdw -20.22 -4.83
>> q vdw bonded -25.2 -6.02
>> single topology vac q/bonded/vdw 10.32 2.47
>> bonded q vdw 9.53 2.28
>> bonded/q vdw 9.1 2.17
>> q bonded vdw 8.71 2.08
>> q bonded/vdw 9.64 2.30
>> q vdw bonded 9.84 2.35
>> The desolvation free energy of ethane is -2.48 and methanol 3.48kcal.
>> The best result which close the cyclo for ethane,vac -> ethane, solv ->
>> methanol, solv -> methanol,vac -> ethane,vac is perturbing all three
>> together(q/bonded/vdw) followed by perturbing bonded/q together then vdw.
>> The worst is to decouple bonded lambdas the last. When I try dual topology
>> method, there is no such difference where to place the bonded lambdas...
>>
>> I am confused that which one should used and why others are not OK?
>>
>> Thank you
>> Best
>> Peng
>> --
>> 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