[gmx-users] free energy calculations of methane in water computed by GMX ver 4.5.7 and ver 4.6.2

[Justin Lemkul]

On 6/21/13 11:07 AM, Dwey wrote:
> Hi gmx-users,
>
>       I almost  reproduced  free energy calculations of methane in water on
> Justin's website. First of all, I am able to follow the workflow of
> computing solvation free energy  for several times with Gromacs version
> 4.5.7 and version 4.6.2 installed in two identical Linux boxes.
>
> However.  the output results of GMX ver 4.5.7 and ver 4.6.2 show different
> values of dG
>
> ##################################
> GMX Ver. 4.5.7:
>
> lambda  0.000 -  0.050,   DG  0.05 +/-  0.00
> lambda  0.050 -  0.100,   DG  0.01 +/-  0.00
> lambda  0.100 -  0.150,   DG -0.03 +/-  0.01
> lambda  0.150 -  0.200,   DG -0.08 +/-  0.00
> lambda  0.200 -  0.250,   DG -0.15 +/-  0.00
> lambda  0.250 -  0.300,   DG -0.21 +/-  0.01
> lambda  0.300 -  0.350,   DG -0.28 +/-  0.00
> lambda  0.350 -  0.400,   DG -0.38 +/-  0.00
> lambda  0.400 -  0.450,   DG -0.50 +/-  0.01
> lambda  0.450 -  0.500,   DG -0.66 +/-  0.01
> lambda  0.500 -  0.550,   DG -0.90 +/-  0.01
> lambda  0.550 -  0.600,   DG -1.21 +/-  0.01
> lambda  0.600 -  0.650,   DG -1.37 +/-  0.01
> lambda  0.650 -  0.700,   DG -1.25 +/-  0.01
> lambda  0.700 -  0.750,   DG -0.96 +/-  0.00
> lambda  0.750 -  0.800,   DG -0.62 +/-  0.00
> lambda  0.800 -  0.850,   DG -0.31 +/-  0.00
> lambda  0.850 -  0.900,   DG -0.03 +/-  0.00
> lambda  0.900 -  0.950,   DG  0.20 +/-  0.00
> lambda  0.950 -  1.000,   DG  0.38 +/-  0.00
>
> total   0.000 -  1.000,   DG -8.31 +/-  0.04
>
> ##################################
>
> GMX ver. 4.6.2 ver
>
> point  0.000 -  0.050,   DG  0.00 +/-  0.00
> point  0.050 -  0.100,   DG -0.03 +/-  0.00
> point  0.100 -  0.150,   DG -0.08 +/-  0.00
> point  0.150 -  0.200,   DG -0.14 +/-  0.00
> point  0.200 -  0.250,   DG -0.20 +/-  0.00
> point  0.250 -  0.300,   DG -0.27 +/-  0.00
> point  0.300 -  0.350,   DG -0.34 +/-  0.00
> point  0.350 -  0.400,   DG -0.43 +/-  0.01
> point  0.400 -  0.450,   DG -0.54 +/-  0.01
> point  0.450 -  0.500,   DG -0.71 +/-  0.01
> point  0.500 -  0.550,   DG -0.94 +/-  0.01
> point  0.550 -  0.600,   DG -1.24 +/-  0.02
> point  0.600 -  0.650,   DG -1.39 +/-  0.02
> point  0.650 -  0.700,   DG -1.28 +/-  0.01
> point  0.700 -  0.750,   DG -1.00 +/-  0.00
> point  0.750 -  0.800,   DG -0.67 +/-  0.00
> point  0.800 -  0.850,   DG -0.36 +/-  0.00
> point  0.850 -  0.900,   DG -0.09 +/-  0.00
> point  0.900 -  0.950,   DG  0.14 +/-  0.00
> point  0.950 -  1.000,   DG  0.33 +/-  0.00
>
> total  0.000 -  1.000,   DG -9.23 +/-  0.03
>
> ##################################
>
>
> The value of DG (= -9.23 kJ mol -1)  by GMX ver 4.6.2 is very close to that
> by Justin or Shirts et.al of 2.24 kcal mol-1  (~ -9.36 kJ mol -1), while
> that of GMX ver 4.5.7 ( =-8.31 kJ mol-1) is far away.
>
> I wonder if someone has similar experience to explain the inconsistency
> between the outputs from ver 4.5.7 (~-8.31 kJmol-1)  and ver 4.6.2 (-9.23
> kJmol-1 ), despite  the values of dG computed by others.
>

I can't explain the difference you are seeing, but I know the results have been 
reproduced with versions 4.5.3 and 4.5.5.  It is good to know that the latest 
version (4.6.2) is working as expected, but it might be useful to figure out 
what's up with 4.5.7.  Please let me know if you can reproduce the right result 
with 4.5.5.

>
> Second,  the reason why I almost reproduced dG  but not completely is
>   because I removed the step of  L- BFGS minimization.  I was not able to
> pass this step  for both GMX versions.  Here shows the warning.
>
> ++++++++++++++++++++++++++++
>
> Fatal error:
> The combination of constraints and L-BFGS minimization is not implemented.
> Either do not use constraints, or use another minimizer (e.g. steepest
> descent).
>
> ++++++++++++++++++++++++++++
>
> L-BFGS mdp file can be found  at
> http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/Files/em_l-bfgs.mdp
>

Add "define = -DFLEXIBLE" to the .mdp file.  Water molecules are constrained via 
SETTLE, so that causes issues with L-BFGS.  When I designed the tutorial, grompp 
did not do that check, so it was ignoring those constraints and silently moving 
past it.  I will fix the .mdp file.  Thanks for bringing that to my attention.

-Justin

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================