[gmx-developers] Optimization of hydrogen positions for the system with constraints

[Leontyev Igor]

Hi,
I need to prepare protein configuration for continuum pKa calculations. A
crystal structure was protonated and hydrogen positions need to be
optimized. Energy minimization with fixed all heavy atoms is carried out for
this purpose. My system consist of the protein 11930 atoms and 4 water
molecules legated to metal centers inside of the protein. I use AMBER
force-field ported to gromacs (TIP3P-ported with amber in ffamber_tip3p.itp
file). The problem is that after the minimization (steep method) positions
of water oxygens are slightly shifted (~0.5 A) from the original
configuration in spite of they are frozen (according to the manual);
whereas, all other protein heavy atoms stay fixed. The shift of 0.5 A for
the legated water oxygens is critical for my continuum calculations. How can
i keep them on original places?
    Excuse me if this question has been already discussed, but I didn't find
directly related questions in the archive. It seems for me that the problem
is related to constraints (settle) used to keep the water rigid. I found in
the posting archive some hints that use of flexible water model may help.
But this does not seem to be the best solution for such trivial problem,
since it was reported in "Revisions for GROMACS version3.3"
http://www.gromacs.org/content/view/18/164/ that "Steepest descents with
constraints now converges properly." Thus, I hope there is some simple
solution.

Bellow is my mdp-file, where the frozen group  "all_heavy_atoms" contains
also those water oxygens.
;====================== Minimization ==========================
integrator = steep    ;algorithm for energy minimization
emtol = 10. ;the convergence criterion for maximum force
emstep = 0.001 ; [nm] - initial step-size
nstcgsteep = 1000 ;frequency of steepest descent step while doing conjugate
gradient
freezegrps = all_heavy_atoms
freezedim = Y Y Y
;==========================================================
dt = 0.002 ; time step
nsteps = 1000 ; number of steps
nstcomm = 1000 ; reset c.o.m. motion
nstxout = 0 ; write coords
nstvout = 0 ; write velocities
nstlog = 1 ; print to logfile
nstenergy = 0 ; print energies
nstlist = 0 ; update pairlist
ns_type = simple ; pairlist method
;================= Cutt off specification =============================
pbc = no ; periodic boundary conditions
coulombtype = Cut-off
rlist = 0 ; cut-off for ns
rvdw = 0 ; cut-off for vdw
rcoulomb = 0 ; cut-off for coulomb
;==============================================================