Box deformation using shear stress in GROMACS depends on forcefield

GROMACS version: 2021.3
GROMACS modification: Yes
Hello,

I have been trying to run a shear simulation on GROMACS (specifically, Couette flow) and I have run into a few issues pertaining to the forcefield that I use. Below is my mdp command for performing the simulations.

;deform
deform = 0 0 0 0.2 0 0

I have performed this under both NVT ( following this reference https://doi.org/10.1016/j.bbagen.2019.06.006) as well as NPT. Here are my thermostat and barostat parameters for the two runs.

; Temperature
Tcoupl = v-rescale
tc_grps = system
tau_t = 0.1
ref_t = 310

; Pressure Pcoupl = Parrinello-Rahman
tau-p = 2.0
ref-p = 0 0 0 0 0 0
pcoupltype = anisotropic
compressibility = 0 0 4.5e-5 0 0 0
refcoord_scaling = com

Initially, I used the CHARMM36 FF for these runs, and I found that for moderate to high shear rates (above 0.2 nm/ps), the simulation failed with this error:


step 575410: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Wrote pdb files with previous and current coordinates

Or


Step 1115555, time 1115.56 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 3862.396973, max 75192.843750 (between atoms 661 and 662) bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
661 662 31.7 0.3379 7496.8257 0.0997
Wrote pdb files with previous and current coordinates
Constraint error in algorithm Lincs at step 1115555
Wrote pdb files with previous and current coordinates

However, the simulation was successful for 0.1 nm/ps.

Then, I tried using the in-built GROMACS forcefields, namely, AMBER99SB-ILDN and OPLS-AA. I ran both under NVT following this reference (https://doi.org/10.1016/j.bbagen.2019.06.006). For these FFs, the simulations were successful irrespective of the shear rate used (I went up to 1 nm/ps).

It is not clear to me why these simulations are FF dependent. For clarification, I have lysozyme in explicit solvent in a large box (1000 nm3). I would prefer to use CHARMM36 so as to remain consistent with another part of my ongoing study. I use GROMACS 2021.3 version. I would appreciate any insight into this issue.

Besides changing force field, have you tried using a more realistic shear rate?