[gmx-users] free energy calculations blowing up

[David]

On Fri, 2005-11-11 at 08:59 -0800, Moore, Jonathan (J) 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
You may want longer time constant for the barostat.

Check the soft-core settings, this may help, there was a change in one
of the default settings in 3.3 that improved stability.

20 lambda points looks good.

> 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.
-- 
David.
________________________________________________________________________
David van der Spoel, PhD, Assoc. Prof., Molecular Biophysics group,
Dept. of Cell and Molecular Biology, Uppsala University.
Husargatan 3, Box 596,          75124 Uppsala, Sweden
phone:  46 18 471 4205          fax: 46 18 511 755
spoel at xray.bmc.uu.se    spoel at gromacs.org   http://xray.bmc.uu.se/~spoel
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++