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.)

In a similar case, I was able to successfully prevent the double helix from unwinding by constraining the distances between the heavy atom centroids of the corresponding terminal bases to their starting distance with a small margin to allow for normal thermal oscillations, along the following lines.

in index.ndx

[ baseA_1 ]
11 12 14 16 17 18 19 20 
[ baseB_28 ]
1753 1754 1755 1758 1759 1760 1762 1763 1764 1765 1766 

in production.mdp

pull_group6_name = baseA_1
pull_group9_name = baseB_28

pull_coord4_type = flat-bottom
pull_coord4_geometry = distance
pull_coord4_groups = 6 9
pull_coord4_dim = Y Y Y
pull-coord4-init = 0.2
pull_coord4-start = yes
pull_coord4_k = 15000 

Thank you @genie. This got me a lot further than most attempts. However, the simulation keeps stopping halfway because it says that the distance between two pull groups is 3 nm, which is larger than 0.49 times the size of the box. This is strange to me because if I look at the pullx.xvg file, it never shows the bases separated by more than about 0.4nm. It’s all pretty frustrating. Thanks for the help with this, though!

I see you are using GROMACS 2019.6. I would suggest updating to the latest version. There was a bug in pull code that was causing it to misreport the groups with too large distance. Basically, for more complex geometries having several defining vectors (such as direction-relative, angle, or dihedral), each vector is checked but the code always reported groups 1 and 2 in the error message. This was fixed sometime in 2020.
Once you know the actual offending distance, you could adjust the group definitions (to make them closer) or set additional restraints to prevent them from drifting apart.

Thank you, genie. I’ll see if I can update. I’m running on a colleagues cluster, so I’ll have to do some bribing! :) I’ll repost here if I can get it working.

I just wanted to give a final update on this.

I never did get the “pull” function to work. As @genie mentioned, it may have been due to my Gromacs install, but as I am running on someone else’s cluster, I was not able to modify that. I tried multiple different definitions for my groups and nothing seemed to work.

However, I recently found out about the “merge” function in pdb2gmx. I used this function to combine my two DNA strands into one topology. (Previously, they were two separate topologies that were included as .itp files in the top-level topology file.)

By merging the two DNA strands into one topology, I was able to create a bonded (restraint) interaction between the first and last basepairs (using the bond “10” function). Despite the fact that this should be identical (at least as far as I understand the documentation), this approach worked perfectly where the “pull” function had failed. The DNA did not melt and I did not receive any errors in the simulation.

Just in case anyone is trying to do the same thing, this is the code I added to the bottom of the [bonds] section of my topology (topol.top):

;	atom_i atom_j	func	r0	r1	r2 k(kJ mol^-1 nm^-2)
   22	618		10	0.0	0.5	0.6	10000
   16	615		10	0.0	0.5	0.6	10000
  300	331		10	0.0	0.5	0.6	10000
  303 	337		10	0.0	0.5	0.6	10000

The atoms indicated are the heavy atoms (nitrogen and oxygens) that form the hydrogen bonds in my terminal A-T bases.

Hope this helps somebody!

Hello Kandresen

I am using the same code you used, I added this code to the bond section in the topology file.
I have one question regarding the values 0.0 , 0.5 , 0.6 ? Are selecting these values critical?

Thank you so much!

Amnah,
From what I understand, the restraint that is being used is a “U”-shaped potential. See the “Distance Restraints” section of the link below:
(Restraints — GROMACS 2022 documentation)

I basically didn’t want any forces imposed unless the ends started to fray. From what I understand, setting the “r0” to zero should get rid of the “U” shape and simply make it a quadratic force that turns on at “r1”. I checked the distances between the atoms I wanted to impose the restraint on and they were all around 0.3nm. Therefore, I (relatively arbitrarily, as I remember) chose the quadratic force to “turn on” at 0.5nm and become linear at 0.6nm. In my case, I was largely interested in inter-DNA (rather than intra-DNA) interactions, and so I really just needed to keep the DNA from melting/fraying during the simulation. The literature I was following seems similarly ambivalent to the details of the restraints, so I feel at least semi-justified in this approach.
So, I guess this is a long way of saying that for me the values of r1 and r2 didn’t seem matter too much as long as the restraints only did what they were supposed to do which was to prevent the DNA from fraying when it exhibited those tendencies.
Full disclosure: I am relatively new to the simulation game, so others may have better/more nuanced answers.