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

[Justin A. Lemkul]

As a follow-up, I found a post that I had remembered seeing before:

http://lists.gromacs.org/pipermail/gmx-users/2003-June/006054.html

Since your system has lots of important LJ interactions between pentane and 
octanol (since pentane is uncharged), you may be seeing effects related to a 
poorly-chosen value of sc_alpha.

-Justin

Justin A. Lemkul wrote:
> 
> 
> 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

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