[gmx-users] Gibbs Energy Calculation, Water/Octanol Partioning
Justin A. Lemkul
jalemkul at vt.edu
Mon Sep 6 03:15:07 CEST 2010
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
> 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
> 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
> 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
> 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