Methionine in CHARMM36 crashes with HMR

User discussions
1

GROMACS version: 2024.2
GROMACS modification: No

Hi, our lab members have been running some protein simulations in CHARMM36 with hydrogen mass repartitioning + 4 fs timestep and there has been a consistent problem with the SD-CE-HE[123] bonds in methionines - they throw LINCS warnings until (often after 100-200 ns) one of them turns into an error and crashes the simulation. Many different particular atoms in the system show this issue, but always involving the same atoms from methionine. We tried both implicit (through .mdp) and explicit (in .top) versions of HMR and the issue happened in both.

I could reproduce this error quickly with an NVT simulation of ~25 capped methionines in vacuum - it seems that the mass imbalance between the sulfur and the repartitioned methyl group makes the methyl C-H oscillations unstable. I tested an ad-hoc solution where another 4 atomic mass units are transferred from the sulfur to the methyl carbon (now also added as an option in Gromologist), and it removed the instabilities.

What surprises me though is that this has not been reported anywhere, so I’m thinking if anyone else saw that behavior, or perhaps there’s something wrong with our setup? I used very “vanilla” .mdp settings, played around with a few of them, and it wasn’t too sensitive to any obvious parameters.

Attaching a minimal structure + topology (with HMR, without the Met modification) if anyone wants to reproduce it independently.
merged_hmass.top (55.9 KB)
25met.pdb.txt (56.1 KB)

Best,
Miłosz

2

AFAIK, using HMR with CHARMM is largely untested. It seems to be widely used in the AMBER community. I wonder what their force constant is for the SD-CE bond compared to CHARMM. The fact that you can alleviate the problem by increasing the CE mass suggests that its oscillatory frequency goes berserk when repartitioning the mass (which is not entirely surprising).

3

Thanks Justin!, good to know it’s mostly untested - I thought there was something wrong with the setup.

I’ve just checked and the difference between Amber and CHARMM is surprisingly small, around 5%:

  S   CT3     1   0.18160   200832.0  ; CHARMM36m
  CT  S       1   0.18100   189953.6  ; Amber99SB-ILDN

but given the LINCS error doesn’t happen a lot, perhaps that’s enough to trigger the instability?

EDIT: I did some simple tests and changing the force constant from 200k to 190k didn’t help that much, at 180k there were just 3 LINCS warnings over 1500 ns*molecules (10 ns, 6 replicas, 25 molecules per replica, all in vacuum) so perhaps somewhere slightly below that value would yield acceptable stability.

That being said, I don’t know why CHARMM should be more sensitive to such instabilities - the S-C-H angle and CG-SD force constants are actually lower in CHARMM36 than in Amber (with no extra U-B terms for angles), so there might be some non-trivial oscillatory modes at play.