[gmx-users] Free energy calculations (FEP) and soft core potential 1-1-48
Michael Shirts
mrshirts at gmail.com
Wed Jul 15 16:30:20 CEST 2015
I put in 1-1-48, and although somewhat more efficient (20-30%), it has
a tendency to have floating point errors in single precision. I'll
probably take it out sometime in the next few versions, as the benefit
is not really worth the hassle. If you are off by a lot, it's almost
certainly not going to help.
Whether something matches experiment is independent of whether the
calculation is performed correctly. If it's off by 2-3 kcal/mol, that
could just be the nature of the problem. If it's off by 5-10
kcal/mol, especially for a relative calculation, there is probably
something else wrong. Possibly it has to do with your top file, but
of course hard to say.
100 lambda states is almost certainly more than needed. When I
disappear 60 heavy atoms, I use only 40 lambda points! If you are
changing 1-5 heavy atoms, no more than 20 lambda would be needed,
10-12 is probably enough.
If there is a lot of protein motion, 1 ns is likely not long enough.
> 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).
It probably would affect something during MD, but what the effect
would be depends on many factors.
On Wed, Jul 15, 2015 at 10:20 AM, Julian Zachmann
<FrankJulian.Zachmann at uab.cat> wrote:
> 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
> --
> 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