[gmx-users] Gibbs Energy Calculation, Water/Octanol Partioning

[Justin A. Lemkul]

Dallas Warren wrote:
> I am working with Emanuel on this, so I will try here to pull all the 
> information together, clarify some things and address all the questions 
> people have asked re further information.
> 
> We are attempting to calculate the partition coefficient (logP between 
> water and 1-octanol) of small molecules.  As a starting point we are 
> attempting to reproduce the results published by Garrido et al. 
> (http://dx.doi.org/10.1021/ct900214y) where they generate the Gibbs 
> energy of solvation in water and 1-octanol, for the alkanes from methane 
> to octane, then this is used to calculate the partition coefficient.
> 

Have you tried using the sd integrator?  That's what the authors of this paper 
used.  I don't know if there are stability issues with md vs. sd, but it's a 
starting point.

> The forcefield we are using is ffG53a6 (which Garrido used as well).  We 
> have been successfully using this same forcefield and Gibbs energy 
> calculation procedure to determine the Gibbs energy of hydration 
> (solvation in water) of small molecules such as methanol, toluene and 
> mono-ethylene glycol.
> 
> The simulation procedure is:
>         1 - minimise a single octanol molecule, L-BFGS then steepest decent
>         2 - fill a box randomly with 200 of these minimised octanol 
> molecules
>         3 - place a single pentane within this octanol box
>         4 - minimise this pentane / octanol box, L-BFGS then steepest decent
>         5 - turn on temperature coupling (v-rescale, T=298K, tau_t = 
> 0.5), constant volume for 50,000 steps, 2fs
>         6 - turn on pressure coupling (Berendsen tau_p = 2.0, P = 1.0), 
> temperature coupling same, for 50,000 steps, 2fs
>         7 - change pressure coupling to Parrinello-Rahman (tau_p = 0.5), 
> temperature coupling same, for 50,000 setps, 2fs
>         8 - production run, 5ns, constant volume / no pressure coupling, 
> temperature coupling as per before, 2,500,000 steps
> 
> A State is the "normal" pentane molecule, B State is the pentane 
> molecule made up of dummy atoms, with zero interactions and charges, but 
> the same mass.  Can see that in the pentane topology file below.
> 
> Follow these links to see the topology files / parameter files:
>         octanol http://hydra.pharm.monash.edu.au/md_project/octanol.txt

Just FYI, this URL is not working for some reason, so I haven't been able to 
look at it.

>         pentane http://hydra.pharm.monash.edu.au/md_project/pentane.txt
>         production run parameter file 
> http://hydra.pharm.monash.edu.au/md_project/production.txt
> 
> We have started with pentane, so Emanuel has run pentane in vacuum for 5 
> ns, and it is stable and the molecular confirmations and movements 
> appear sane.  Moving to pentane in octanol, it is being simulated with 
> 16 λ (lambda) values (reproducing Garrido).  For λ=0. 0.05, 0.1, 0.2, 
> 0.3, 0.4 and 0.5 the 5 ns simulations of pentane in octanol are stable.  
> Once lambda reaches 0.6 and above, the simulation ends with a LINCS error:
> 
> "       Fatal error:
> Too many LINCS warnings (1001)
> If you know what you are doing you can adjust the lincs warning 
> threshold in your mdp file or set the environment variable 
> GMX_MAXCONSTRWARN to -1, but normally it is better to fix the problem."
> 
> The atoms involved in this error are all the single pentane molecule and 
> it appears to be essentially folding up on itself, appearing to be due 
> to the repulsion of the neighbouring alkane atoms of the octanol.  At 
> lambda 0.6 it fails at 1.7ns, and the time taken to fail decrease with 
> increasing lambda to 0.48 ns when lambda is 1.0
> 
> Our initial thoughts was that it may be something to do with the soft 
> core parameter, causing some discontinuity in the interactions between 
> the atoms, though that does not appear to be the case.  In order to may 
> be provide some insight into what is the cause of the issue, the lambda 
> 0.6 simulations have been repeated changing the time step and value of 
> the soft core parameter (sc_alpha).
>         1 - 2fs and sc_alpha = 1.51, fails at 1.7 ns, original situation 
> mentioned above
>         2 - 1fs and sc_alpha = 1.51, completes fine
>         3 - 2fs and sc_alpha = 0.50, fails at 0.214 ns
>         4 - 1fs and sc_alpha = 0.50, completes fine
> 
> It did appear that using the time step allows it to complete, but the 
> thought was that this was simply covering up the issue.  So, repeating 
> all 16 lambda using 1fs and sc_alpha 1.51, now it is fine up to lambda 
> 0.75 and fails now for 0.80, 0.85, and 0.90.  When lambda is 0.95 and 
> 1.00, it completes fine.
> 

As noted in the Garrido et al. paper (citations to Michael Shirts' and David 
Mobley's work), as well as the free energy tutorial that I believe David Mobley 
wrote, combining sc_alpha != 0.5 and sc_power == 1 is not stable.  A combination of:

integrator = sd
...
sc_power = 1
sc_alpha = 0.5

should be the most stable.  I am unsure of the underlying theoretical basis for 
these assessments, but they seem to be quoted frequently.  Maybe someone else 
can comment.

-Justin

> Repeating these Gibbs energy calculation of pentane in water (versus 
> octanol above) runs fine for all lambda values without any issues.  A 
> box of pentane only is also stable without any problems.  From this it 
> appears that it is not an issue with the pentane topology, and it 
> shouldn't be since that is such a simple molecule and are using the same 
> alkane parameters published many times by others.
> 
> Our thinking is that it may be due to the fact that the softcore is 
> having difficulty with the alkane chains of the octanol as the atoms are 
> appearing / disappearing.  Or is it a bug with the software?
> 
> Catch ya,
> 
> Dr. Dallas Warren
> Medicinal Chemistry and Drug Action
> Monash Institute of Pharmaceutical Sciences, Monash University
> 381 Royal Parade, Parkville VIC 3010
> dallas.warren at monash.edu
> +61 3 9909 9304
> ---------------------------------
> When the only tool you own is a hammer, every problem begins to resemble 
> a nail.
> 

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

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

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