[gmx-users] Free energy calculations (FEP) and soft core potential 1-1-48

[Julian Zachmann]

Dear Gromacs Users,

I am running free energy calculations (FEP) to estimate relative binding
affinities. So far, my results don't match the experimental results (not
even close), so I must be doing something wrong.

My protocol is the following:

- Energy minimisation
- Leap frog minimisation
- 0.1ns  NVT
- 2ns    Simulated Annealing up to 350K
- 0.1ns  NPT
- 1ns     MD production run

During the whole equilibration time I am using umbrella pulling to maintain
an interaction between one residue and the ligand (especially important
during SA) but I remove it for MD because I am afraid it would influence my
results (even though I don't know it for sure).

The simulation time for 1ns is quite short of course but I have 'nstdhdl'
put to 10 to get a lot of output.

The other settings are pasted below. So far I was using 1-1-6 for the soft
core potential but because the results were not good, I was thinking to try
1-1-48. However in this case the simulations always crash. Has anybody run
successfully 1-1-48 soft core potentials in Gromacs and how did you chose
your settings? Like me pasted below?

Any other suggestions about what I am doing wrong?
Could I maintain the pull distance constraint during the MD simulation?

Best regards,
Julian

free_energy              = yes
init_lambda_state        = $i
;                           0      1      2      3      4      5      6
 7      8      9     10     11     12     13     14     15     16     17
  18     19     20     21     22     23     24     25     26     27     28
    29     30     31     32     33     34     35     36     37     38
39     40     41     42     43     44     45     46     47     48     49
  50     51     52     53     54     55     56     57     58     59     60
    61     62     63     64     65     66     67     68     69     70
71     72     73     74     75     76     77     78     79     80     81
  82     83     84     85     86     87     88     89     90     91     92
    93     94     95     96     97     98     99
fep_lambdas             = 0.0000 0.0300 0.0600 0.0900 0.1200 0.1500 0.1800
0.2100 0.2400 0.2700 0.3000 0.3300 0.3600 0.3900 0.4200 0.4500 0.4800
0.5100 0.5400 0.5700 0.6000 0.6300 0.6600 0.6900 0.7200 0.7500 0.7800
0.8100 0.8400 0.8700 0.9000 0.9300 0.9600 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000
vdw_lambdas             = 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0500 0.1000 0.1500
0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500
0.7000 0.7500 0.8000 0.8200 0.8500 0.8700 0.8800 0.8900 0.9000 0.9100
0.9200 0.9300 0.9350 0.9400 0.9450 0.9500 0.9530 0.9570 0.9600 0.9630
0.9670 0.9700 0.9730 0.9770 0.9800 0.9820 0.9840 0.9860 0.9880 0.9900
0.9910 0.9920 0.9930 0.9940 0.9950 0.9955 0.9960 0.9965 0.9970 0.9975
0.9980 0.9983 0.9987 0.9990 0.9992 0.9993 0.9994 0.9995 0.9996 0.9997
0.9998 0.9999 1.0000
; 1-1-6
sc-coul                    = no
sc-alpha                 = 0.5
sc-power                 = 1
sc-sigma                 = 0.3
sc-r-power               = 6
; 1-1-48
;sc-coul                   = no           ;
;sc-alpha                 = 0.0025       ; either 0.5 (sc-r-power = 6) or
0.0025 (sc-r-power = 48)
;sc-power                = 1            ; 1 or 2
;sc-sigma                = 0.3          ; default 0.3
;sc-r-power              = 48           ; possible 6, 48
calc-lambda-neighbors    = -1
nstdhdl                  = 10
dhdl-derivatives         = yes
separate-dhdl-file       = yes