Unexpected geometry change when combing `freezegrps` with `constraints` and running MD

User discussions
1

GROMACS version: 2026.0-dev
GROMACS modification: No
Dear all,
I attempted to combine freezegrps and constraints in an MD simulation of a single GC base pair. Guanine should be frozen, whereas the hydrogen bonds in cytosine should be constrained with LINCS. When I run the simulation, I notice that guanine’s geometry in the first frame of the .trr file differs slightly from the structure used to build the .tpr. Below are the relevant input files.
This is my .mdp file:

title                   = NVT equilibration 
; define                = -DPOSRES -DPOSRES_FC_BB=400.0 -DPOSRES_FC_SC=40.0
; Run parameters
integrator              = md        ; leap-frog integrator
nsteps                  = 500000    ; 1 * 500000 = 500 ps
dt                      = 0.001     ; 1 fs
; Output control
nstxout                 = 500       ; save coordinates every 0.5 ps
nstvout                 = 500       ; save velocities every 0.5 ps
nstenergy               = 500       ; save energies every 0.5 ps
nstlog                  = 500       ; update log file every 0.5 ps
nstxout-compressed      = 500       ; save .xtc file every 0.5 ps
; freeze groups
freezegrps              = DCM
freezedim               = Y Y Y
; Bond parameters
continuation            = no        ; first dynamics run
constraint_algorithm    = lincs     ; holonomic constraints 
constraints             = h-bonds   ; bonds involving H are constrained
lincs_iter              = 1         ; accuracy of LINCS
lincs_order             = 4         ; also related to accuracy
; Nonbonded settings 
cutoff-scheme           = Verlet    ; Buffered neighbor searching
ns_type                 = grid      ; search neighboring grid cells
nstlist                 = 10        ; 20 fs, largely irrelevant with Verlet
rcoulomb                = 1.0       ; short-range electrostatic cutoff (in nm)
rvdw                    = 1.0       ; short-range van der Waals cutoff (in nm)
DispCorr                = EnerPres  ; account for cut-off vdW scheme
; Electrostatics
coulombtype             = PME       ; Particle Mesh Ewald for long-range electrostatics
pme_order               = 4         ; cubic interpolation
fourierspacing          = 0.16      ; grid spacing for FFT
; Temperature coupling is on
tcoupl                  = V-rescale              ; modified Berendsen thermostat
tc-grps                 = DCM      DGM           ; two coupling groups - more accurate
tau_t                   = 0.1      0.1           ; time constant, in ps
ref_t                   = 300      300           ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl                  = no        ; no pressure coupling in NVT
; Periodic boundary conditions
pbc                     = xyz       ; 3-D PBC
; Velocity generation
gen_vel                 = yes       ; assign velocities from Maxwell distribution
gen_temp                = 300       ; temperature for Maxwell distribution
gen_seed                = -1        ; generate a random seed

This is my .ndx file:

[ System ]
   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
  16   17   18   19   20   21   22   23   24   25   26   27   28   29
[ Other ]
   1    2    3    4    5    6    7    8    9   10   11   12   13   14   15
  16   17   18   19   20   21   22   23   24   25   26   27   28   29
[ DCM ]
   1    2    3    4    5    6    7    8    9   10   11   12   13
[ DGM ]
  14   15   16   17   18   19   20   21   22   23   24   25   26   27   28
  29

Below is the .gro used to generate .tpr file:

c_dry_centered, used to generate .tpr
   29
    9DCM     N1    1   1.969   1.801   1.641
    9DCM     C6    2   1.936   1.932   1.637
    9DCM     H6    3   2.019   2.006   1.646
    9DCM     C5    4   1.806   1.969   1.623
    9DCM     H5    5   1.777   2.076   1.620
    9DCM     C4    6   1.709   1.864   1.613
    9DCM     N4    7   1.579   1.892   1.599
    9DCM    H41    8   1.548   1.990   1.595
    9DCM    H42    9   1.509   1.814   1.592
    9DCM     N3   10   1.745   1.735   1.617
    9DCM     C2   11   1.874   1.700   1.631
    9DCM     O2   12   1.912   1.582   1.636
    9DCM    H1'   13   2.068   1.770   1.651
   16DGM    H9'   14   1.275   1.139   1.572
   16DGM     N9   15   1.273   1.243   1.571
   16DGM     C8   16   1.162   1.323   1.559
   16DGM     H8   17   1.060   1.280   1.548
   16DGM     N7   18   1.191   1.450   1.561
   16DGM     C5   19   1.329   1.454   1.575
   16DGM     C6   20   1.419   1.564   1.584
   16DGM     O6   21   1.392   1.685   1.580
   16DGM     N1   22   1.551   1.523   1.598
   16DGM     H1   23   1.625   1.601   1.605
   16DGM     C2   24   1.593   1.392   1.603
   16DGM     N2   25   1.725   1.372   1.617
   16DGM    H21   26   1.758   1.275   1.621
   16DGM    H22   27   1.794   1.452   1.624
   16DGM     N3   28   1.510   1.288   1.595
   16DGM     C4   29   1.381   1.325   1.581
   3.20680   3.20680   3.20680

and below is the .gro structure obtained from the 0-th frame in my .trr file:

t=   0.00000 step= 0
   29
    9DCM     N1    1   1.971   1.801   1.641 -0.0000 -0.0000 -0.0000
    9DCM     C6    2   1.937   1.933   1.637 -0.0000 -0.0000 -0.0000
    9DCM     H6    3   2.018   2.005   1.646 -0.0000 -0.0000 -0.0000
    9DCM     C5    4   1.806   1.970   1.623 -0.0000 -0.0000 -0.0000
    9DCM     H5    5   1.777   2.075   1.620 -0.0000 -0.0000 -0.0000
    9DCM     C4    6   1.709   1.864   1.613 -0.0000 -0.0000 -0.0000
    9DCM     N4    7   1.577   1.892   1.599 -0.0000 -0.0000 -0.0000
    9DCM    H41    8   1.548   1.989   1.595 -0.0000 -0.0000 -0.0000
    9DCM    H42    9   1.511   1.816   1.592 -0.0000 -0.0000 -0.0000
    9DCM     N3   10   1.745   1.735   1.617 -0.0000 -0.0000 -0.0000
    9DCM     C2   11   1.874   1.700   1.631 -0.0000 -0.0000 -0.0000
    9DCM     O2   12   1.912   1.582   1.636 -0.0000 -0.0000 -0.0000
    9DCM    H1'   13   2.066   1.770   1.651 -0.0000 -0.0000 -0.0000
   16DGM    H9'   14   1.275   1.142   1.572  0.9068 -0.4631  2.3615
   16DGM     N9   15   1.273   1.243   1.571 -0.0831 -0.4663 -0.3907
   16DGM     C8   16   1.162   1.323   1.559  0.5987  0.1923  0.3309
   16DGM     H8   17   1.063   1.281   1.548  2.2023 -3.6754 -0.2586
   16DGM     N7   18   1.191   1.450   1.561 -0.1524 -0.3400 -0.4246
   16DGM     C5   19   1.329   1.454   1.575  0.5423  0.6697  0.6809
   16DGM     C6   20   1.419   1.564   1.584 -0.0895  0.1212  0.3290
   16DGM     O6   21   1.392   1.685   1.580  0.1033  0.2589 -0.4021
   16DGM     N1   22   1.551   1.523   1.598 -0.3913 -0.4047 -0.5681
   16DGM     H1   23   1.621   1.596   1.605 -1.2829  0.3268  0.9949
   16DGM     C2   24   1.593   1.392   1.603 -0.3328 -0.5607 -0.0017
   16DGM     N2   25   1.725   1.372   1.617 -0.1477  0.5208  0.3347
   16DGM    H21   26   1.757   1.277   1.621  0.7214  0.7433 -1.0418
   16DGM    H22   27   1.791   1.449   1.624  1.1344 -0.5371  0.0282
   16DGM     N3   28   1.510   1.288   1.595  0.1262  0.1770  0.2288
   16DGM     C4   29   1.381   1.325   1.581 -0.4091  0.1337 -0.0218
   3.20680   3.20680   3.20680

I expect the geometry for guanine remains unchanged, but that’s not the case. The deviations are not due to rigid translation or rotation, and the geometry for guanine remains the same during MD.However, if I disable constraints, then the frozen guanine retains its original geometry. This may suggest a conflict between freezegrps and constraints. Could something be wrong with my .mdp settings? Any advice would be greatly appreciated. Thank you in advance!

2

Yes, there is a conflict. Doesn’t grompp warn about this?

3

Thank you for your explanation! I reran gmx grompp but didn’t seem to find the warning you mentioned:

                    :-) GROMACS - gmx grompp, 2026.0-dev (-:

Executable:   /***/gromacs-main/share/gromacs/bin/gmx
Data prefix:  /***/gromacs-main/share/gromacs
Working dir:  /***/test
Command line:
  gmx grompp -f f_nvt.mdp -c c_dry_centered.gro -r c_dry_centered.gro -p topol.top -o t_nvt.tpr -n index.ndx

Ignoring obsolete mdp entry 'title'
Ignoring obsolete mdp entry 'ns_type'
Setting the LD random seed to 2010320221

Generated 3916 of the 3916 non-bonded parameter combinations
Generating 1-4 interactions: fudge = 0.5

Generated 3916 of the 3916 1-4 parameter combinations

Excluding 3 bonded neighbours molecule type 'Other'

turning H bonds into constraints...

Setting gen_seed to 1809508863

Velocities were taken from a Maxwell distribution at 300 K

NOTE 1 [file f_nvt.mdp]:
  There are 13 atoms that are fully frozen and part of COMM removal
  group(s), removing these atoms from the COMM removal group(s)

Number of degrees of freedom in T-Coupling group DCM is 0.00
Number of degrees of freedom in T-Coupling group DGM is 40.00

The largest distance between excluded atoms is 0.369 nm between atom 18 and 22

Determining Verlet buffer for a tolerance of 0.005 kJ/mol/ps at 300 K

Calculated rlist for 1x1 atom pair-list as 1.000 nm, buffer size 0.000 nm

Set rlist, assuming 4x4 atom pair-list, to 1.000 nm, buffer size 0.000 nm

Note that mdrun will redetermine rlist based on the actual pair-list setup
Calculating fourier grid dimensions for X Y Z
Using a fourier grid of 24x24x24, spacing 0.134 0.134 0.134

Estimate for the relative computational load of the PME mesh part: 0.95

NOTE 2 [file f_nvt.mdp]:
  The optimal PME mesh load for parallel simulations is below 0.5
  and for highly parallel simulations between 0.25 and 0.33,
  for higher performance, increase the cut-off and the PME grid spacing.



This run will generate roughly 5 Mb of data

There were 2 NOTEs

Back Off! I just backed up t_nvt.tpr to ./#t_nvt.tpr.2#

GROMACS reminds you: "I Am Testing Your Grey Matter" (Red Hot Chili Peppers)

Could you please tell me where to find it? And is there a method I may use to avoid this conflict?

4

Are you maybe freezing dimensions along less than three dimensions?

5

Thanks for the suggestion! Unfortunately, I think it doesn’t quite fit my system.I’m simulating a DNA duplex fragment in water and need to freeze one specific base pair in all three dimensions so as to sample several DNA-and-water environments that favor the geometry of that base pair. Restraints may not work here, because I want that base pair geometry to remain completely unchanged.

6

My question was not a suggestion, but for checking why you didn’t get a warning in grompp.

But I realized now that I actually fixed the case of constraints between a frozen and a non-frozen atom. But you have all atoms in a molecule frozen, correct? We should add a grompp warning for that case.

7

That’s correct. In my case, all atoms in one molecule are frozen. Thanks again for help!