[gmx-users] free energy calculations blowing up

Moore, Jonathan (J) JMoore2 at dow.com
Fri Nov 11 22:48:59 CET 2005


My problem is that it's one of someone else's many jobs to maintain the Gromacs on our clusters, and I can't easily just do it myself on those production machines.  Therefore, I try to minimize the frequency of my requests to update to the latest version.

Jonathan 
____________________________ 
Jonathan Moore, Ph.D. 
Research and Engineering Sciences - New Products 
Core R&D 
The Dow Chemical Company 
1702 Building, Office 4E 
Midland, MI 48674  USA 
Phone:  (989) 636-9765 
Fax: (989) 636-4019 
E Mail: jmoore2 at dow.com 

-----Original Message-----
From: David Mobley [mailto:dmobley at gmail.com] 
Sent: Friday, November 11, 2005 4:40 PM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] free energy calculations blowing up


Jonathan,


On 11/11/05, Moore, Jonathan (J) <JMoore2 at dow.com > wrote:
David and David,

Thanks for the suggestions.  I'm not currently using soft cores for fear of
the soft core bug in 3.2.1.
http://www.gromacs.org/pipermail/gmx-users/2004-August/011830.html 

I haven't moved to 3.3 yet on any of my production machines because I'm
holding out for 3.3.1, but maybe I'll have to give up waiting so I can use
soft cores.

You _could_ just change that bit of the source code as described in that post and recompile. I was using a CVS version of 3.2.1 from Nov. 2004 which had that bugfix until very recently and the softcore stuff was working fine for me.

It's worth noting, though, that there are number of other bugfixes implemented in 3.3 which affect free energy calculations; it's not clear to me whether the version I was using of 3.2.1, or the version you are using, had all of these correctly implemented. 

And just by way of advertisement, I'm currently using 3.3 and have been able to compute free energies for turning off charges and vdW for a number of small molecules in water. I've done so with the AMBER force field (from the Pande group) and my results match what I was getting with the AMBER package. This gives me a reasonable amount of confidence that things are working properly. I've also computed the charging free energy of spc water in water, the answer of which is known, and I get the correct answer. 
  
I have, however, included this bugfix http://www.gromacs.org/pipermail/gmx-users/2005-October/017495.html in my 3.3 installation. It's not clear to me that without the bugfix I necessarily would have had problems, but...




David M., I guess this is the paper you were talking about, right?
T. C. Beutler, A. E. Mark, R. C. van Schaik, P. R. Gerber and W. F. van 
Gunsteren. "Avoiding singularities and numerical instabilities in free
energy calculations based on molecular simulations." Chem. Phys. Lett., 222,
529-539 (1994).
         http://dx.doi.org/10.1016/0009-2614(94)00397-1

Yep, that's the one.

David

 


Thanks,
Jonathan
____________________________
Jonathan Moore, Ph.D.
Research and Engineering Sciences - New Products
Core R&D
The Dow Chemical Company 
1702 Building, Office 4E
Midland, MI 48674  USA
Phone:  (989) 636-9765
Fax: (989) 636-4019
E Mail: jmoore2 at dow.com 

-----Original Message-----
From: David Mobley [mailto: dmobley at gmail.com]
Sent: Friday, November 11, 2005 3:27 PM
To: Discussion list for GROMACS users
Subject: Re: [gmx-users] free energy calculations blowing up 


Jonathan, 

To clarify slightly, the singularity for lambda->1 is when the lambda=1
Hamiltonian corresponds to a completely disappeared particle. If you are
changing on vdW radii into another vdW radii things are somewhat more 
complicated and it isn't clear to me that you really have the same sort of
problem. I wouldn't be surprised if you did, though, but it would be at some
intermediate lambda value. I would still try the soft core option, 
especially if your time series of dV/dlambda have anything that looks at all
like shot noise.

David



On 11/11/05, David Mobley < dmobley at gmail.com> wrote: 
Jonathan,

Two comments:
(1) Although this wasn't your question, you might not want to run the
simulations sequentially. I think previous work has shown that you can
introduce some hysteresis this way, although the reference eludes me at the 
moment.
(2) When changing vdW parameters, soft core helps a lot. From what you're
describing it sounds like you are doing that. If you don't use soft core,
dV/dlambda actually has a singularity for particle insertion as lambda->1 
(assuming the exponent of lambda is 1, which it is by default) and is also
very susceptible to shot-type noise, and susceptible to crashes when forces
get too large (which often happens near lambda=1). You can reduce the 
crashes somewhat by using smaller timesteps but this doesn't help with the
shot noise or singularity problems. To understand why this happens you can
look up van Gunsteren's paper from around 1992 introducing the soft core 
potential, or contact me directly for a brief explanation. If you have any
trouble finding the reference I can look it up for you, but I don't have it
in front of me right this instant.

David Mobley
UCSF 





On 11/11/05, Moore, Jonathan (J) <JMoore2 at dow.com > wrote:

I have a question about free energy calculations. 

I'm performing several different simulations where the hydrogen of certain 
hydroxyls on my molecule is replaced with a methyl group.  It's a united
atom model, so the change from hydrogen to methyl doesn't change the number
of sites.  I have a single molecule of interest in SPC water (see more 
details below).  I'm experiencing the problem that on occasion the system
blows up.  I've had no problems like this previously when simulating only
one state or the other (methylated or unmethylated), so I think the problem 
is due to the free energy calculation.  Unless I've made an error in setting
up the free energy calculation, I assume my problem is that the system isn't
stable for the timestep that I'm using (2 fs) and how fast I'm changing 
lambda (I'm incrementing lambda by 0.05 followed by 200 ps of simulation
before increasing it again).

Is it common to either decrease the time step or use lambda increments
smaller than 0.05 when doing free energy calculations this way?  I'm 
planning to try both to see if they fix the problem, but I'm curious if
anyone else has had a problem like this.  I'm not using soft cores, but
maybe I have to?

Also, I'm specifying the atom types like this ( i.e., with the mass being
taken from the .atp file):

[ atoms ]
;   nr    type   resnr  residu    atom    cgnr  charge
     6       H       1    H000     HO3       2   0.410

Therefore, when I specify the B state, I do it like this (i.e. with out
specifying the change in mass, assuming the GROMACS will also get the B
state mass properly from the .atp file):

[ atoms ]
;   nr    type   resnr  residu    atom    cgnr  charge
     6       H       1    H000     HO3       2   0.410     CH3   0.180

Should that work OK?  It wasn't obvious to me where to look in the output to
make sure it is.  I assume I'll be able to answer that question myself when
3.3.1 is released and I can use "gmxcheck -ab" 

More details:
v. 3.2.1, single precision
2 fs timestep
All bond lengths constrained (LINCS)
Nosé-Hoover thermostat with the polymer and solvent coupled separately; 0.4
ps time constant
Berendsen barostat with 0.5 ps time constant and 4.5x10-5 bar -1
compressibility
Triple-range cut-off method

Thanks,
Jonathan

____________________________
Jonathan Moore, Ph.D .
Research and Engineering Sciences - New Products 
Core R&D
The Dow Chemical Company
1702 Building, Office 4E
Midland, MI 48674  USA
Phone:  (989) 636-9765
Fax: (989) 636-4019
E Mail: jmoore2 at dow.com

_______________________________________________ 
gmx-users mailing list
gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-request at gromacs.org .
_______________________________________________ 
gmx-users mailing list
gmx-users at gromacs.org
http://www.gromacs.org/mailman/listinfo/gmx-users
Please don't post (un)subscribe requests to the list. Use the 
www interface or send it to gmx-users-request at gromacs.org.



More information about the gromacs.org_gmx-users mailing list