[gmx-users] Gibbs Energy Calculation, Water/Octanol Partioning
Justin A. Lemkul
jalemkul at vt.edu
Mon Sep 6 03:34:15 CEST 2010
As a follow-up, I found a post that I had remembered seeing before:
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 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
> 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
>> 2 - fill a box randomly with 200 of these minimised octanol
>> 3 - place a single pentane within this octanol box
>> 4 - minimise this pentane / octanol box, L-BFGS then steepest
>> 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
>> 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.
>> 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
ICTAS Doctoral Scholar
Department of Biochemistry
jalemkul[at]vt.edu | (540) 231-9080
More information about the gromacs.org_gmx-users