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

[Peng He]

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