[gmx-users] Expanded Ensemble and Gromacs 4.6

[Joakim Jämbeck]

Hi,

I am trying to use the Expanded Ensemble (EE) method to compute the free energy of solvation of a small organic molecule.

Basically I am playing around but I cannot get the simulations to run.

Here are my EE and free energy settings:

%-------------------------------------------------------
free-energy 		= expanded	; Expanded ensemble
couple-moltype 		= C1X		; Molecule to introduce
couple-lambda0 		= none		; Go from no interactions with solvent...
couple-lambda1 		= vdw-q		; ...to full interactions
couple-intramol 	= no		; Do not decouple internal interactions
init_lambda_state 	= 0		; Start from the first column of the lambda vectors
delta-lambda 		= 0		; No increments in lambda 
nstdhdl 		= 100		; Frequency for writing dH/dlambda
coul-lambdas 		= 0.0 0.25 0.5 0.75 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
vdw-lambdas 		= 0.0 0.0 0.0 0.0 0.0 0.2 0.4 0.6 0.7 0.8 0.9 1.0
dhdl-print-energy 	= yes		; Include total energy in the dhdl file	
sc-alpha 		= 0.5		; Soft core alpha parameter
sc-power 		= 1		; Power of lambda in soft core function 
sc-r-power 		= 6		; Power of radial term in the soft core function
sc-sigma 		= 0.3		; Soft core sigma
sc-coul 		= no		; No soft core for Coulomb
separate-dhdl-file 	= yes		; Seperate dhdl files
dhdl-derivatives 	= yes		; Print derivatives of the Hamiltonian

; expanded ensemble variables
nstexpanded 		= 100		; Number of steps between attempts to change the Hamiltonian
lmc-stats 		= wang-landau	; WL algorithm to explore state space
lmc-move 		= gibbs		; Decides which state to move to
lmc-weights-equil 	= number-steps	; EE weight updating stops after... 
weight-equil-number-steps = 5000000	; ...10ns of simulation

; Seed for Monte Carlo in lambda space
lmc-seed 		= 1993
lmc-repeats 		= 1
lmc-gibbsdelta 		= -1
lmc-forced-nstart 	= 0
symmetrized-transition-matrix = no
nst-transition-matrix 	= -1
mininum-var-min 	= 100
wl-scale 		= 0.6
wl-ratio 		= 0.8
init-wl-delta 		= 1
wl-oneovert 		= yes
%-------------------------------------------------------

Besides this I use LINCS to constrain all bonds, sd-integrator and a "normal" cut-off scheme with a 1 fs time step. 

However, once I try to run the files on the cluster I always end up with LINCS warnings and after a few seconds the program crashes due to too many LINCS warnings. If I increase the time step to 2 fs I run into problems SETTLE for the water. I always start from a nice structure that is taken as the final snap shot from a simple 1 ns MD simulation of my system.

What could be the problem?

If I use free_energy = yes instead of EE things work fine. 

Did I perhaps mess up the EE settings or something?

Thanks in advance.

Best,
Joakim