[gmx-users] question of perturbing order in single topology FEP, gromacs
Peng He
hepeng.yanglu at gmail.com
Thu Nov 30 18:17:32 CET 2017
Hi, Asaf,
I totally agreed what you suggested, but I am considering the bonded energy
term changing in the free energy simulation in the single topology, I have
read the paper you suggested and it is using a dual topology method which
also works well in my case, in the single topology method, the bond
equilibrium length is changing according to lambda values. for example, in
ethane -> methanol case, the shrinking(changing the bonded lambdas) of cc
bond to co bond(the second carbon is changing to oxygen gradually and
shares one atom in the topology file) is gonna affect the 1,4 electrostatic
interaction of hydrogens connecting to carbon and carbon/oxygen center.
While in dual topology there is no such problem because the charge is only
disappearing and appearing with no position change. I am wondering if the
single topology method relies on the order of perturbation a lot and which
order is the most correct way.
Best
Peng
2017-11-30 9:27 GMT-05:00 <asaffarhi at post.tau.ac.il>:
> Dear Peng,
>
> There are several ways to do it.
>
> The simplest way is to take each molecule and perturb it separately. You
> have one Ethane molecule and one methanol molecule so you need a file
> containing only Ethane and a file containing only methanol.
>
> The different atoms between the molecules are CH3 in Ethane and OH in
> methanol.
> In each molecule you can remove the VDW and Coul terms of these atoms in
> transformations. That's all. See Fig. 9.
> The transformed systems should be identical in their remaining nonbonded
> terms and bonded terms relating the common atoms so an additional minor
> transformations may be needed.
>
> There may be some issues related to changing the total charge in
> transformations, which can cause artificial changes in the calculated free
> energies.
>
> Best regards,
>
> Asaf
>
>
> Quoting Peng He <hepeng.yanglu at gmail.com>:
>
> 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.
>>>
>>> --
>> 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