[gmx-users] How to delete bad water molecules

[Chris Neale]

/>>While minimizing,// be sure to include
>>/
//>>/ define=-DFLEXIBLE
/>>
>//>/ in your .mdp file
/
>Just out of curiosity, why do you think minimization should be done
>with flexible water? I've gone the opposite direction and only
>minimize with steepest descents since that's the only minimization
>algorithm in gromacs that will work with constraints on hbonds, which
>is how I run my dynamics, and I want to be minimizing on the same
>potential energy landscape in which I'm running dynamics. That seemed
>important. I'm curious why your strategy is different.

You make a good point. In fact I was not even aware that it was possible to do steep EM without flexible water and
I often do not actually use the l-bfgs step. 

To my mind the importance of minimizing on the same energy landscape as the one that will be used for MD is 
simply as a way to ensure that the MD will not crash. There are two things that I desire from EM: 
1) allows MD not to crash; 2) perturbs the system as little as possible. I believe that flexible waters will allow
the system to energy minimize while perturbing the xray protein structure as little as possible, being less perturbing
than rigid waters. Still, I admit that this is just something that I have decided now since I learned that DFLEXIBLE is not manditory.
 
I have just run the following test and it indicates to me that flexible water can be useful sometimes during EM. 
An EM run with flexible water is ok and without flexible water crashes. The result surprised me a little as I only expected
the flexible water to reduce the energy more than rigid water. The initial spc water box has some contacts that are too close. 
For example SOL292:HW1 is 0.9A from SOL470:HW2.

On the other hand, I don't change my .mdp options in any other way and this means that in l-bfgs I get
(and ignore) the warning message:
WARNING 1 [file em.mdp, line unknown]:
  For efficient BFGS minimization, use switch/shift/pme instead of cut-off.

Chris.

###

Here is the test that I ran:

echo "title">a.gro;echo "0">>a.gro;echo "0 0 0">>a.gro
editconf -f a.gro -o b.gro -box 3 3 3
genbox -cp b.gro -cs ~/gromacs/top/spc216.gro -o c.gro -p c.top
pdb2gmx -f c.gro -o c.pdb -p c.top -ff oplsaa -water spc
grompp_d -f em.mdp -c c.gro -p c.top -o c.tpr
mdrun_d -deffnm d -v

...
Step=  499, Dmax= 1.3e-03 nm, Epot= -4.61644e+04 Fmax= 4.06646e+02, atom= 250
Step=  500, Dmax= 1.6e-03 nm, Epot= -4.61676e+04 Fmax= 1.02330e+03, atom= 250

writing lowest energy coordinates.

Steepest Descents did not converge to Fmax < 10 in 501 steps.
Potential Energy  = -4.61676470800464e+04
Maximum force     =  1.02329833107139e+03 on atom 250
Norm of force     =  2.88377252842351e+03

$cat em.mdp
title               =  EM
cpp                 =  cpp
define              = -DFLEXIBLE
integrator          =  steep
nsteps              =  500
nstlist             =  10
ns_type             =  grid
pbc                 =  xyz
coulombtype         =  PME
rcoulomb            =  0.9
fourierspacing      =  0.12
pme_order           =  4
vdwtype             =  cut-off
rvdw_switch         =  0
rvdw                =  1.4
rlist               =  0.9
DispCorr            =  no
gen_seed            =  9896
constraints         = none

############

Then if I modify em.mdp by removing the define statement:


Steepest Descents:
   Tolerance (Fmax)   =  1.00000e+01
   Number of steps    =          500
Step=    0, Dmax= 1.0e-02 nm, Epot= -2.16756e+04 Fmax= 1.43528e+04, atom= 379
Step=    1, Dmax= 1.0e-02 nm, Epot= -2.47981e+04 Fmax= 5.96853e+03, atom= 2464
Step=    2, Dmax= 1.2e-02 nm, Epot= -2.79251e+04 Fmax= 3.12114e+03, atom= 283
Step=    3, Dmax= 1.4e-02 nm, Epot= -3.04027e+04 Fmax= 1.55279e+03, atom= 2077
Step=    4, Dmax= 1.7e-02 nm, Epot= -3.29729e+04 Fmax= 1.00606e+03, atom= 1409
Step=    5, Dmax= 2.1e-02 nm, Epot= -3.51535e+04 Fmax= 1.04599e+03, atom= 1726
Step=    6, Dmax= 2.5e-02 nm, Epot= -3.64712e+04 Fmax= 2.28944e+03, atom= 1726
Step=    7, Dmax= 3.0e-02 nm, Epot= -3.70197e+04 Fmax= 9.91473e+02, atom= 1726
Step=    8, Dmax= 3.6e-02 nm, Epot= -3.80289e+04 Fmax= 3.64312e+03, atom= 1726
Step=    9, Dmax= 4.3e-02 nm, Epot= -3.83579e+04 Fmax= 9.21816e+02, atom= 283

Back Off! I just backed up step9.pdb to ./#step9.pdb.1#

Back Off! I just backed up step10.pdb to ./#step10.pdb.1#
Wrote pdb files with previous and current coordinates
Segmentation fault