[gmx-users] Strange behavior from steepest descents

[Justin A. Lemkul]

Hi all,

I'm preparing a few membrane protein systems using the InflateGRO script from 
Tieleman's site.  The bilayer consists of 340 POPE lipids, using structures and 
topologies from the same site.  I've prepared several systems (with other 
lipids) using the methodology described, but am running into some very strange 
behavior with the POPE system.

I scaled the coordinates for the lipids by a factor of 4 and attempted a 
steepest descents minimization, my usual procedure.  This step is essentially an 
in vacuo minimization of scaled lipids in the presence of a protein.  The 
potential energy decreases nicely, and levels off.  Suddenly, the potential 
plummets again, after which mdrun cuts out, reporting that it has converged to 
machine precision (although it has not reached my desired Fmax < 1000).  The 
forces were incredibly large (greater than 10^11), so there is definitely a 
problem with the system that I have yet to figure out.

What troubles me is that mdrun believes the steepest descents algorithm has 
converged, when clearly it has not.  I have posted an image here to illustrate 
the potential energy:

http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/Images/Gromacs/weird_energy.jpg

Conjugate gradients does not give the same result, instead converging to a 
potential energy of roughly -4 x 10^5, a reasonable magnitude, given the size of 
the system.  The forces, however, remain large (10^5 or so), indicating I 
probably still have some clashes to resolve.

Still I am troubled by the steepest descents result.  Is this indicative of some 
problem?  I usually use steepest descents, and have never had a problem with it 
before.  I see the same result in serial or parallel (on 2 CPU's of my dual core 
laptop), using single-precision GMX 3.3.3.

Thanks for reading yet another long email from me.  The .mdp file I used for 
steepest descents (that gives the weird behavior) is below.  The behavior 
persists even if I lower emtol and emstep, it just takes a bit longer to arrive 
at the same conclusion :)

-Justin

=======em.mdp========

; LINES STARTING WITH ';' ARE COMMENTS
title           = Minimization  ; Title of run

; Run control
cpp             = /usr/bin/cpp  ; Preprocessor
define          = -DSTRONG_PR
include         = -I../         ; Directories to include in the topology format

; Parameters describing what to do, when to stop and what to save
integrator      = steep
emtol           = 1000.0        ; Stop minimization when Fmax < 1000
emstep          = 0.01          ; Energy step size
nsteps          = 50000         ; Maximum number of steps
energygrps      = system        ; Which energy group(s) to write to disk

; Neighbor searching and electrostatics
nstlist         = 1
ns_type         = simple
rlist           = 1.0
coulombtype     = cut-off
rcoulomb        = 1.0
rvdw            = 1.0
pbc             = xyz

=======================

-- 
========================================

Justin A. Lemkul
Graduate Research Assistant
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================