[gmx-users] Re: stepsize too small ... but potential energ negative!

[Anna Marabotti]

Dear Justin,
I would like to say that you were true, but it isn't! :-(
I launched the minimization steps first changing in my em.mdp file the
emstep from 0.1 to 0.01 as you suggested, then using directly your minim.mdp
file from the lysozyme tutorial to create the .tpr:

integrator          =  steep
emtol               =  1000.0
emstep              =  0.01
nsteps              =  50000
nstlist             =  1
ns_type             =  grid
rlist               =  1.0
coulombtype         =  PME
rcoulomb            =  1.0
rvdw                =  1.0
pbc                 =  xyz

In both cases, the calculation stops with the same error as previous. The
log file reports a similar (infinite) value for LJ and Coulomb at the first
step, then nothing more. 
I'm convincing myself that the problem is in the protein structure, not in
the parameters I used for minimization, but I really don't know what's
wrong. Moreover, I don't understand why mdrun goes for 37 steps without
giving me any output values. It doesn't say "I cannot minimize your
structure because it's too wrong for me", it doesn't fail to minimize it
arriving at a very high force or positive potential energy: it apparently
goes for some steps, but actually it seems to do nothing. It's not only the
fact that the structure doesn't minimize that is strange, in my opinion it's
more stranger the fact that I haven't any communication about something
wrong in my structure.

Any other help will be appreciated.
Thank you very much 
Anna
------------------------------

Message: 2
Date: Mon, 24 May 2010 10:08:27 -0400
From: "Justin A. Lemkul" <jalemkul at vt.edu>
Subject: Re: [gmx-users] Re: stepsize too small ... but potential
	energy	negative!
To: Discussion list for GROMACS users <gmx-users at gromacs.org>
Message-ID: <4BFA885B.9040305 at vt.edu>
Content-Type: text/plain; charset=ISO-8859-1; format=flowed

Anna Marabotti wrote:
> Dear Justin, dear Luca,
> here's the answer to your questions:
> - I'm currently using the "classical" forcefield gromos96 43a1 (choice "0"
> in pdb2gmx). After producing the topology, the only warning I see from
> pdb2gmx is this one: WARNING: there were 0 atoms with zero occupancy and
63
> atoms with occupancy unequal to one (out of 1583 atoms). Check your pdb
> file.
> However, atoms with occupancy <1 are present also in "regular" PDB files
(if
> I remember well, also in PDB files I used previously). Is this a problem?
> -  Box: I'm currently using a cubic box and I'm setting 1 nm of distance
> between solute and box with option -d (and also center the box). Looking
at
> the system I can't see any contact between the protein and the box walls.
I
> started by setting 0.8 nm and the problem was the same.
> - .mdp file: here it is:
> cpp                 =
> define              =  -DFLEXIBLE
> constraints         =  none
> integrator          =  steep
> nsteps              =  50000
> emtol               =  1000
> emstep              =  0.1

This step size is far too large!  Try something more like 0.01.

> nstlist             =  5
> ns_type             =  grid
> rlist               =  1.0
> coulombtype         =  PME
> rcoulomb            =  1.0
> rvdw                =  1.2
> fourierspacing      =  0.12
> fourier_nx          =  0
> fourier_ny          =  0
> fourier_nz          =  0
> pme_order           =  4
> ewald_rtol          =  1e-5
> optimize_fft        =  yes
> Tcoupl              =  no
> Pcoupl              =  no
> gen_vel             =  no
> It's very similar to the one Justin suggested in its tutorial.

...except for numerous changes.

> - PR-MD: I'not interested in skipping the minimization and continue with
MD.
> I tried to launch a PR-MD step to see if the error produced in this step
was
> more informative, and try to understand what was the problem on my
> structure. I failed, however...
> 
> I would also add that this is not the first time that I'm using an
homology
> model produced by MODELLER to perform MD. All the checks I made on my
model
> tell me that it is a good model, so I really don't understand what's wrong
> with it.
> 
> Finally, I started from the structure deprived of the first residue in
order
> to see if that residue was the "bad" one, but the problem still persists
and
> the potential energy has the same negative value with the same order of
> magnitude.
> 

Likely due to a flawed .mdp file.  The algorithm is trying to take too large
of 
a step along the potential energy surface, causing bad geometry and infinite

forces.  A bit more finesse should solve this problem.

-Justin