[gmx-users] Errors in Minimization using grompp

[Justin A. Lemkul]

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

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

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

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