[gmx-users] question of perturbing order in single topology FEP, gromacs

asaffarhi at post.tau.ac.il asaffarhi at post.tau.ac.il
Thu Nov 30 15:27:11 CET 2017 

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.

More information about the gromacs.org_gmx-users mailing list