NPT Step Issue with Protein 1QS4 (4 Waters and Mg²⁺ in Active Site)

I am currently working on a MD simulation involving the protein 1QS4 with GROMACS, which contains four crystallographic water molecules and a Mg²⁺ ion in its active site. I have performed QM/MM calculations using Gaussian to generate the topology and GRO file for my system.

To stabilize the active site, I defined a position restraint file (posre.itp) for the protein, the four waters, and the Mg²⁺ ion, excluding hydrogen atoms, with a restraint force constant of 100,000. Using this setup, I could completed the energy minimization (EM) and NVT equilibration steps.

However, during the NPT equilibration step, I encountered LINC warnings, and the simulation stops. I have taken the following steps to resolve the issue, but without success:

  1. Changed the barostat to Berendsen.
  2. Increased the pressure coupling time constant (tau_p) to 5 ps.

Despite these changes, the issue persists. I am using GROMACS version 2018 for my simulations.
I suspect the issue might be related to the restrained water molecules or Mg²⁺ ion in the active site during NPT equilibration, but I am uncertain how to proceed.
I will be grateful some one can guide me for this issue

Going from a restraint force constant of 100,000 to 0 is quite a large step. If you needed such a high force constant for NVT, I would recommend running two more simulations of NVT equilibration with restraint force constants of 10,000 and 1,000, respectively. Apart from that, how long was your NVT equilibration?

Thank you very much for your guidance, i did NVT step with nsteps: 500000 and dt: 0.002

I also had to add a force of 100,000 during the EM step aslo before starting NVT step; otherwise, the simulation crashes. When I reduced the force to 10,000, this step also crashed.

Assuming that your force field parameters are correct, it seems like the system is unstable or in a very unfavourable state. If the simulation crashes when reducing the restraint force to 10,000 there are two options I’d try:

  1. Running another step of EM (after the first) with a lower restraint force constant, e.g. 10,000 or 50,000.
  2. See if it works running NVT with a restraint force of 50,000 (after running with 100,000).
    or
  3. Extend the NVT simulation to, e.g., 10,000,000 steps (20 ns equilibration is still not extremely much by any means) with the same restraint force constant (100,000).