[gmx-users] Errors in Minimization using grompp

Justin A. Lemkul jalemkul at vt.edu
Wed Jun 13 14:10:22 CEST 2012

On 6/13/12 7:55 AM, reisingere at rostlab.informatik.tu-muenchen.de wrote:
> Hi!
> Thank you for your answer.
>> On 13/06/2012 6:44 PM, reisingere at rostlab.informatik.tu-muenchen.de wrote:
>>> Hi everybody,
>>> I want to do a minimization of the hydrogens of my protein. Only the
>>> hydrogens.
>> Doesn't really matter - nothing will move more than a fraction of an
>> Angstrom unless it's horribly wrong, in which case not moving the heavy
>> atoms won't help you.
>>>    And I want to do this with implicit solvent.
>> Probably not worth the effort if you're just preparing for MD.
> Our goal is not an MD simulation of the protein. We plan a
> Poisson-Boltzmann electrostatics calculation. For that we need the heavy
> atoms as they are in the crystal structure (even 1 Angstroem movement
> would be too much) with good hydrogen atom position.

With restraints on heavy atoms, your positions will not deviate very much at 
all.  A 1 Angstrom movement would be huge in this case; I would expect your 
deviations to be orders of magnitude less.

> We need to minimize, simulate, minimize, simulate the hydrogen atoms.

I don't follow the logic here.  You say you need to do an EM of the H atoms in 
order to do some PBSA calculation with no MD, but here you're doing two 
iterations of EM and MD.

> Implicit solvent is ok for us. (For a later MD run of the complete
> protein, we will use explicit solvent).
> This is the .mdp file:

There are several problems here, most of which I've already stated, but I'll 
recapitulate them again.

> define              =  -DPOSRES
> constraints         =  all-bonds
> integrator          =  steep
> nsteps              =  30000
> vdwtype             =  cutoff
> coulombtype         =  cutoff

As Mark said, these settings should be "cut-off" not "cutoff."

> epsilon_rf = 0
> pbc                 =  xyz

For an implicit system, periodicity should be set to "no" and grompp will warn 
you about using angular COM removal (invoked below) if you leave it set as is.

> nstlist             =  1
> ns_type             =  grid
> rlist               =  1
> rcoulomb            =  1.2
> rvdw                =  1.2

Finite cutoffs with implicit solvent will lead to instability and poor energy 
conservation.  These three lines (rlist, rcoulomb, rvdw) should be set to zero, 
unless you're satisfied with artifacts.

> rvdw_switch     = 0.7

This setting has no effect when using cutoffs.

> comm-mode           =  angular
> comm-grps           =  System
> nstcgsteep = 1000
> emtol               =  5.0

Unless you've compiled with double precision, it is unlikely you will achieve a 
minimal force this low.  The use of restraints is going to preclude most 
movement in the system, so don't be surprised if grompp does not (and cannot) 
achieve this tolerance.

> emstep              =  0.01
> implicit_solvent    =  GBSA
> gb_algorithm        =   HCT
> nstgbradii          =  1
> rgbradii            =  1

This should also be set to zero if the above neighbor searching parameters are zero.

> gb_epsilon_solvent  =  80
> gb_saltconc       =  0
> sa_algorithm        =  Ace-approximation
> sa_surface_tension  = -1
> Additionally I have a question according to the vdwtype and coulombtype.
> Why do I have to set the two parameters to cutoff?

Because settings like PME don't make sense for implicit solvent calculations. 
Plain cutoffs with finite values are almost never adequate for any purpose, however.



Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080


More information about the gromacs.org_gmx-users mailing list