[gmx-users] Strange behavior from steepest descents

[Justin A. Lemkul]

I may have discovered the solution to my own problem, so I'll post it here just 
as a follow-up.  It appears that inflation of a POPE bilayer causes extremely 
favorable interactions within the PE headgroup, causing the amine moiety to 
collapse into the phosphate.  This generates the near infinite forces that I'm 
seeing.  I guess somehow the minimization algorithm still finds this 
electrostatic interaction to be somehow favorable, since the potential is 
plummeting.

I suspect that the EM procedure cuts out because the potential is dropping too 
fast for the system to handle, hence why I'm seeing that the system has 
converged to machine precision.  Maybe the numbers are getting to heavy for my 
poor computer's memory :)

Just a few thoughts.

-Justin

Justin A. Lemkul wrote:
> 
> 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

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