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

[Dallas Warren]

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.

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
        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.

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.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20100906/7105cfd4/attachment.html>