Constrained dynamics of oxygen atoms of water molecules

User discussions
1

GROMACS version: 2021.6
GROMACS modification: No

Hi,

I am trying to constrain the positions of certain atoms during an MD simulation (keep cartesian coordinates constant during the whole simulation). The only ways I have come up with so far is by setting the mass very high or by using the freezegrps option on the mdp file, but both fail for constraining the oxygen atoms of water molecules due the error described below.

Whenever I try setting the oxygen mass to 9.999e+8, I get an error message that says:

     Fatal error:
     Can not invert matrix, determinant = 7.313649e+26

Using mdp options freezegrps = OW and freezedim = Y Y Y , I get:

     Fatal error:
     Can not invert matrix, determinant = -inf

Do you have any alternative suggestions on how I can constrain the positions of atoms within gromacs, or how I can make these simulations work.

This is my input file and my system corresponds to a box with 1728 water molecules (TIP3P).

integrator = md
dt = 0.002
nsteps = 50000
nstenergy = 100
nstlog = 100
nstxout-compressed = 100
nstxout = 100
nstvout = 100
nstfout = 100
compressed-x-grps = System
;
cutoff-scheme = Verlet
nstlist = 20
rlist = 1.2
coulombtype = pme
rcoulomb = 1.2
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw_switch = 1.0
rvdw = 1.2
;
tcoupl = v-rescale
tc_grps = System
tau_t = 1.0
ref_t = 303.15
;
pcoupl = no
;
constraints = h-bonds
constraint_algorithm = LINCS
gen_vel = yes
gen_temp = 303.15
;
nstcomm = 100
comm_mode = linear
comm_grps = System
;
refcoord_scaling = com

Thanks!!

2

Have you tried using position restraints? If you use gmx make_ndx to create a group containing only the oxygen atoms, you can then use gmx genrestr to create position restraints for that group, it will place a harmonic force to keep the atoms in place.

I’m not sure about setting the oxygen mass super high, but I’d found from a previous post that freeze groups don’t work for all simulation setups

3

Thanks for the reply!
I saw that gmx genrestr generates an index file that can be used for freezing atoms, by using the -of flag. I tried selecting just the atoms I want to freeze and adding this to my topol.top:

#include “freeze.ndx”

But I get this error when running grompp:

ERROR 1 [file freeze.ndx, line 1]:
Invalid directive freeze

It is not clear to me how I should use this file, how do I give this information to grompp.

And just to clarify, I actually need the atoms to be constrained (fixed coordinates in time). I can’t have them oscillating through the application of an harmonic potential.

Thanks again! I am sharing all the files of a simulation of a protein in water where I am trying to constrain 6 atoms.
freeze.ndx.txt (26 Bytes)
em.gro.txt (116.3 KB)
topol.top (6.8 KB)
md.mdp (965 Bytes)