[gmx-users] question of perturbing order in single topology FEP, gromacs
asaffarhi at post.tau.ac.il
asaffarhi at post.tau.ac.il
Sun Dec 3 08:49:50 CET 2017
Hi Peng,
You can compare the relative free energy with the other topology to
see which order probably gives you the correct results. I don't know
what is your project and for free energies associated with bonds where
the atom (Carbon or Oxygen) is without nonbonded interactions there
are analytic results
http://iopscience.iop.org/article/10.1088/1367-2630/18/2/023039/meta
Best,
Asaf
Quoting Peng He <hepeng.yanglu at gmail.com>:
> 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.
>>
> --
> 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