[gmx-users] cutoff questions

[Ken Rotondi]

Hello all,

I'm working on minimizing and annealing to room temp. a system of 
protein, ions and solvent. I run the minimizations in three steps with 
the protein frozen, restrained, and free. A sample .mdp file is shown:

title            = Chicken
cpp              = /lib/cpp
define           = -DFLEXIBLE -DPOSRES 	; needed for minimization & 
protein restraint
constraints      = none		;
integrator       = l-bfgs	;
nsteps           = 10000	;
emtol            = 10		; max force < (kJ/mol nm)
emstep           = 0.01		; initial step size (nm)
nbfgscorr		 = 10
rlist			 = 1.4
nstlist			 = 10
ns_type			 = grid
coulombtype		 = switch
rcoulomb		 = 1.0
rcoulomb-switch	 = 0.9
vdwtype			 = switch
rvdw			 = 1.0
rvdw-switch		 = 0.9

The only difference amongst the min.mpd files is the presence of freeze 
statements and omission of the '-DPOSRES' flag as needed for the 
'frozen' and free minimizations. The minimizations proceed quite well, 
converging below the desired Fmax each time:

Maximum force     =  9.4426432e+00 on atom 11099
Maximum force     =  8.8269501e+00 on atom 3824
Maximum force     =  6.9491901e+00 on atom 6710

Problems arise in the first 2.5 ps run at 2.5K, at step 48 where I get:

t = 0.096 ps: Water molecule starting at atom 10406 can not be settled.
Check for bad contacts and/or reduce the timestep.Wrote pdb files with 
previous and current coordinates

This message occurs for a different water for the next 4 steps followed 
by these LINKS warnings:

Step 53, time 0.106 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 112.400597 (between atoms 786 and 787) rms 6.366230

Followed by:

Step 54, time 0.108 (ps)  LINCS WARNING
relative constraint deviation after LINCS:
max 7585225216.000000 (between atoms 786 and 787) rms 677410432.000000

Suggesting that they've moved apart 7.5 METERS in two steps... aka, 
blows-up!

My real confusion results from identical runs (a script is used) made 
with the three minimization.mdp files altered only in rvdw, rcoulomb, 
their respective switch functions, nsteps and emtol:  eg.,

title            = Chicken
cpp              = /lib/cpp
define           = -DFLEXIBLE  	; needed for minimization
constraints      = none		;
integrator       = l-bfgs	;
nsteps           = 100000	;
emtol            = 400		; max force < (kJ/mol nm)
emstep           = 0.01		; initial step size (nm)
nbfgscorr		 = 10
rlist			 = 1.4
nstlist			 = 10
ns_type			 = grid
coulombtype		 = switch
rcoulomb		 = 1.4
rcoulomb-switch	 = 1.2
vdwtype			 = switch
rvdw			 = 1.4
rvdw-switch		 = 1.2
freezegrps		 = PROTEIN
freezedim		 = Y Y Y

In this case, the minimization converges to below machine precision 
each time without reaching the desired Fmax:

Maximum force     =  7.3445282e+02 on atom 11387
Maximum force     =  9.4466376e+02 on atom 6077
Maximum force     =  9.8665173e+02 on atom 4949

However, this run anneals smoothly to room temperature producing a 
structure with very low RMSD to the crystal structure.

I don't understand why structures that are better minimized blow up. I 
don't understand why it takes 49 steps for problems to crop up. Why 
such a difference? I appear to be missing something fundamental here.

Thanks, and apologies for the long post,

Ken

-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: text/enriched
Size: 3747 bytes
Desc: not available
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20050817/71f054ea/attachment.bin>