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

[Justin A. Lemkul]

Anna Marabotti wrote:
> 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.
> 

This situation is really bizarre.  Can you send me your coordinate file and 
topology (off-list) so I can have a look?  It will be much faster if I can see 
for myself what's going on.

-Justin

> 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
> 
> 

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

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

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