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

Wed Nov 29 19:22:39 CET 2017

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.
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.
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.