Martini3-IDP and lincs errors?

[jessniblo]

Hello!

I've been running simulations of synthetic IDPs with no problems, but have been running into issues using biological IDPs with the Martini3-IDP extension. Currently I am running simulations with mpirun -np 12 and 4 threads per task, and my workflow is:

1. Generate initial conformation via a short Campari simulation
2. Martinize the all atom conformation using the developer version on martinize, run an energy minimization on just the protein, use insane to solvate the box
3. Run a minimization (steep, emtol=10, emstep=0.02)
4. Run relaxation (dt=0.002, nsteps = 5000000 (10ns), v-rescale on the system with tau_t=1.0, ref_t=298, parrinello-rahman with tau_p=12)
5. Run Relaxation 2 (dt=0.01, nsteps = 5000000 (50ns), v-rescale on the system with tau_t=1.0, ref_t=298, parrinello-rahman with tau_p=12)

I'm currently seeing that my simulations fail late within the second relaxation (step 4724457/5000000) with errors like:

Step 4724457, time 47244.6 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.055367, max 0.271242 (between atoms 38 and 39)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
38 39 90.0 0.3050 0.3877 0.3050

Step 4724458, time 47244.6 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 1.550811, max 7.597391 (between atoms 38 and 39)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
38 39 90.0 0.3877 2.6222 0.3050
Wrote pdb files with previous and current coordinates

Step 4724459, time 47244.6 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.104112, max 0.510044 (between atoms 38 and 39)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
32 33 37.6 0.2700 0.2700 0.2700
38 39 90.0 2.6222 0.4606 0.3050
Wrote pdb files with previous and current coordinates

Eventually this LINCS error transforms to:

step 4724500:
New coordinates: 35910.852 -24566.178 1.637
Old cell boundaries in direction Y: 3.900 7.800
New cell boundaries in direction Y: 3.903 7.805
Atom 10622 moved more than the distance allowed by the domain decomposition (3.902529) in direction Y
Atom 9966 moved more than the distance allowed by the domain decomposition (3.902529) in direction Y
distance out of cell -14.462538
New coordinates: 57.769 -14.463 10.702
Old cell boundaries in direction Y: 0.000 3.900
distance out of cell 15.198914

The lincs errors match my IDP, but the domain decomposition error corresponds to water beads.

I've been trying to work this out to no avail for days now. Does anyone have any tips on how to get Martini3-IDP to work without the lincs error?

Thank you!
Jess

[WangLiguo-kyrie]

My suggestions would be: 

1. make sure position restraints when "run an energy minimization on just the protein";
2. add one more NVT eq between EM and NPT stage for better relaxation;
3. change PR pressure coupling to Berendsen during NPT eq, because PR would cause stronger fluctuation when pressure changing;
4. Usually  we don't need very careful eq relaxation with gradually increasing step size, so the NPT eq with stepsize=0.002 wouldn't be necessary;
5. Because CG models use much longer MD stepsize, usually 20fs in Martini, so sometimes we would meet Lincs warnings. If these problems occur infrequently, I  think it's still okay, you can just continue your simulation based on the updating checkpoint file. If you want to avoid Lincs warnings, adjust 'lincs-warnangle' in mdp file would help, like adjust the default warning angle from 30 to 90, but please make sure it doesn't harm your simulation results.

Hope this helps!

[WangLiguo-kyrie]

I have run plenty of proteins consisting of Trp, they were all okay, so TRP is not the reason.
The virtual site in TRP has specific bead type, so it's meaningless to define VS 0.00 0.000 V 0.0 0.0 in atomtypes

[jessniblo]

I do think this is an issue with my pipeline, not the parameters, but would really appreciate some help identifying what component of my pipeline is causing the issue. I have simulated ~150 IDPs with both the original Martini3 and new Martini3-IDP parameters, and while sequences without a virtual bead work beautifully, I have not been successful in getting a sequence with a virtual bead to complete a production run with either set of parameters. Additionally, I've had issues using carbohydrates included in Martini3 [https://github.com/orgs/Martini-Force-Field-Initiative/discussions/46], and these carbohydrates all include VS too. Mixing a sequence that completed a production run that doesn't have a VS with sucrose that has a VS causes a domain decomposition error, so I do think my pipeline is having issues handing the virtual sites properly.

I've attached a zip file (jkniblo_issueRun.zip) including my entire workflow and the error files that I am seeing. I'm using Gromacs 2023.2 on my university's HPC system, which uses a slurm manager. I submit MS_2task.slurm to the queue, so this file will give you all the commands that I am using, and SystemPrep.sh will give all system preparation. The production run here crashed immediately so the errors are all seg faults (run_prod.out) and in run.err at line 82497 you can see the prior warning about the oscillation period.

Thank you for possibly looking through this and any advice you can give on how to get my pipeline to work with all sequences!

[WangLiguo-kyrie]

I just had a testing using your protein itp, and I can have a stable simulation for 100ns NPT eq and 100ns NPT production. Since you mentioned you have tried one more NVT stage and Berendsen pressure coupling wiht 4ps constant according to my suggestions, I think our pipeline should be close in principle. But I still attached my mdp files for reference.

Another point is that, when at first I used your IDR.gro file, which generated by polyply, it indeed crashed quickly. And then I regenerate the initial protein conformation by polyply again to avoid protein splitted by pbc box. So the initial conformation maybe the issue, making VS calculation unstable when box size fluctuation? Maybe you could check that further.
test.zip

[jessniblo]

Just wanted to send you a quick note thanking you so much for all of your help! My NVT pipeline was indeed very similar to your workflow (just longer run times), but I did notice that you were running everything threaded and without MPI. I gave that a try with Gromacs 2023.2 that was built on my university's HPC system, and had a few runs work but most still failed with the same issues. I built a copy of Gromacs 2021.5, and that was able to get everything working with no problem. I'm not sure if this is a specific issue with my HPC system's 2023 version, but it does seem that this was a Gromacs version issue.

[WangLiguo-kyrie]

glad to know you find the setting allowing stable simulations. It would be valuable to figure out if there are some Virtual site related issues in latest Gromacs versions.