Constraining terminal basepairs in DNA

GROMACS version: 2019.6
GROMACS modification: No
I am simulating a short (10bp) DNA molecule using the CHARMM36 force field. The DNA is melting, as is known to happen in this force field. Previous work has gotten around this by constraining the hydrogen bonds in the terminal bases. I have attempted this by first defining groups in an index that represent each atom in the hydrogen bonds and then using the “pull” code below:

pull_ncoords            = 4         ; four hydrogen bonds
pull_ngroups            = 8         ; eight groups defining four reaction coordinate 
pull_group1_name        = H62 
pull_group2_name        = O4
pull_group3_name        = N1 
pull_group4_name        = H3
pull_group5_name        = O4b 
pull_group6_name        = H61b
pull_group7_name        = H3b 
pull_group8_name        = N1b
pull_coord1_type        = constraint  ; hard constraint
pull_coord1_geometry    = distance 
pull_coord1_dim         = Y Y Y
pull_coord1_groups      = 1 2
pull_coord1_start       = yes       ; define initial COM distance > 0
pull_coord1_rate        = 0      ; Not actually pulling

(I repeat the last parts 3 more times for each hydrogen bond.)

However, when I run this, I get an error similar to
"3 particles communicated to PME rank 9 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension y."
and the simulation fails. This does not happen when I do not use the constraints. I have tried constraining the heavy molecules (exchanging the hydrogens for the heavy atoms they are bound to) and the residues themselves and have tried using the “harmonic” constraint rather than the “constraint”, however none of them work. One thing to note is that I did not equilibrate with these constraints, but since the heavy atoms are position restrained, I assumed this was not necessary.
Any insight would be greatly appreciated.

I would only try to restrain distances between heavy atoms. If you apply biases to H atoms, then the integrator is simultaneously trying to constrain the associated covalent bond while applying the biasing force. Restrain the distance between the donor and acceptor heavy atoms instead and it should be more stable.

Thank you. I thought this might be the case. However, I tried this (using the corresponding heavy atoms) and still received the same error. It is strange to me that these rigid constraints are causing these errors when the unconstrained system works well.

The issue has become much stranger. As I simply am trying to keep the terminal basepairs from melting, I changed my pulling type to “flat-bottom” with a 1.0nm initial distance (the reference distances are ~0.3nm). GROMACS still terminated with the error above. However, if I check the “pullx.xvg” and “pullf.xvg” files, it looks as if the constraint was never even triggered. (All x values are less than the constraint and all f values are 0.)