[gmx-users] free energy calculations blowing up

[Moore, Jonathan (J)]

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