[gmx-users] Energy min ion problems

[Ken Rotondi]

Hello all,

I'm having a weird and frustrating problem using the l-bfgs minimizer 
with a system containing ions. My system is a protein in a 7.5nm 
dodecahedral box, ~30,000 atoms. It has a total charge of -4, so I used 
genion to add 4 Na+ ions to the system. The resulting system could not 
be E-minimized. Several attempts lead to the "converged to machine 
precision in -1 step", yes, it didn't do a thing. I converted to a .pdb 
and  saw that the added ions were in close proximity to, but not 
overlapping with any of the protein atoms. A subsequent attempt to run 
E-m on this system started, but ran only 68 steps prior to converging 
to machine precision. Final output:

Pot. E.  		-3.12X10**23		
Fmax   		1.42X10**26		HUGE!
norm of force	8.59X10**23		ibid

I went back and ran the minimization on the system to which the 
counter-ions had been added. Same system, just 4 extra HOH's and 4 
fewer Na+'s. It converged to machine precision in 11833 steps:

Pot. E.		-3.23X10**-6
Fmax		717				wanted < 10, but okay
norm of force	38

Has anyone seen similar behavior due to addition of ions? What could 
cause such behavior?  Does anyone have a suggestion as for how to 
proceed? I could just add the ions to the E-mized system and see how 
the resultant system minimizes (which I will, in fact, proceed with 
now), but I think this is a fundamental question worthy of discussion.

I've attached my .mdp file, though I doubt that it's the source of the 
error since it functions fine on the non-ionized (grin) system.

Any help/input would be greatly appreciated.

Many thanks,

Ken

title                    = Chicken
cpp                      = /lib/cpp
define                   = -DFLEXIBLE
constraints              = none
integrator               = l-bfgs
nsteps                   = 10000000
emtol                    = 10           ; max force < (kJ/mol nm)
emstep                   = 0.01         ; initial step size (nm)
nbfgscorr                = 10
rlist                    = 1
nstlist                  = 10
ns_type                  = grid
coulombtype              = switch
rcoulomb                 = 1.2
rcoulomb-switch          = 1
vdwtype                  = switch
rvdw                     = 1.2
rvdw-switch              = 1