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.