Constraints Setup for Small Molecule Solvent o-Dichlorobenzene in Martini3 Force Field

[zhangtrzhangtr]

Hello, when I completed the Martini3 model of my organic polymer and mixed it with a large amount of the small molecule solvent o-dichlorobenzene, I used the o-dichlorobenzene model you previously published. During the molecular dynamics simulation process with a timestep of 30fs, when the simulation reached between 40 to 50 nanoseconds, it's possible that due to the interaction between the organic polymer and the small molecule o-dichlorobenzene, the constraints in the o-dichlorobenzene's top file could not be satisfied. This led to the termination of the simulation due to a large number of LINCS WARNINGS. I would like to ask if it is reasonable to replace the constraints in the o-dichlorobenzene top file with bond parameters and at the same time set a larger force constant to make the bond lengths within the o-dichlorobenzene molecule more rigid.
;;;;;; 1,2-DICHLOROBENZENE or 1,3-DICHLOROBENZENE

[moleculetype]
; molname nrexcl
DCLBZ 1

[atoms]
; id type resnr residu atom cgnr charge
1 SX3 1 DCLBZ R1 1 0
2 SX3 1 DCLBZ R2 2 0
3 TC5 1 DCLBZ R3 3 0

[constraints] replace [bonds]
; i j funct length
1 2 1 0.313 ; cog average of 1,2- and 1,3-dichlorobenzene 1 2 1 0.313 8000
1 3 1 0.305 ; cog average of 1,2- and 1,3-dichlorobenzene 1 3 1 0.305 8000
2 3 1 0.305 ; cog average of 1,2- and 1,3-dichlorobenzene 2 3 1 0.305 8000

[ricalessandri]

@zhangtrzhangtr we don't recommend running systems that contain a lot of tiny beads with a time step > 20 fs. The 30 fs timestep is the parameter that may cause issues. Since you closed this Discussion, let me know if this problem still persists (with dt = 20 fs) and I can try to help you further.

[zhangtrzhangtr]

If I adjust the timestep to 20fs and then conduct an NPT molecular dynamics simulation mixing my created organic molecule model with 3000 o-dichlorobenzene molecules, LINCS warnings still occur when the simulation reaches 45 nanoseconds. Since I have defined define=-DFLEXIBLE in the mdp file, my organic molecule will have its constraints converted to bond length parameters during the molecular dynamics simulation, so this will not cause LINCS warnings. However, GROMACS indicates that the constraints between the beads constituting o-dichlorobenzene cannot be satisfied, resulting in numerous LINCS warnings and causing the simulation to fail. I aim to achieve as high accuracy as possible for my created molecular model, so I have used tiny beads throughout, with the entire molecule containing about 40 tiny beads and 10 virtual sites. I am considering whether I can also replace the constraints in the o-dichlorobenzene molecules with bond lengths. Additionally, I have conducted a 30ns test, and the bond length distribution of my created molecular model shows good consistency with the bond length distribution obtained through all-atom molecular dynamics simulations.

[zhangtrzhangtr]

I believe that converting constraints into bond length parameters can make the model more flexible and avoid the instability caused by molecules being too rigid, which results in the inability to satisfy the constraints between them. If I replace the constraints in the o-dichlorobenzene molecules with bond lengths, and the bond length distribution obtained from the tests shows good consistency with the bond length distribution obtained from all-atom molecular dynamics , would this approach be reasonable? Thank you for your guidance.

[zhangtrzhangtr]

It is also possible that there are issues with the molecular model I designed because I found that even after converting the constraints of the o-dichlorobenzene molecules to bond lengths, the simulation still fails after about 50 nanoseconds at a 20 fs timestep. Perhaps the design of the bond lengths and angles between the structures is not reasonable.

[zhangtrzhangtr]

Oh , I've finally found the reason. It was a problem with the GROMACS version. I was using GROMACS 2022, and I conducted coarse-grained molecular dynamics using my own Martini 3 molecular model, your previously published caffeine model, and various solvent mixtures. After about 40 ns, the system would crash. When I reinstalled GROMACS 2018, my model was able to run long-term dynamics stably and normally. It seems that the latest version of GROMACS is temporarily unable to properly utilize the Martini 3 force field.

[ricalessandri]

Thanks for the updates, @zhangtrzhangtr. Indeed some Gromacs version have been showing some issues with virtual sites and/or constraints. As far as I (and some colleagues) have tested, also Gromacs 2023.3 should be good.

If I replace the constraints in the o-dichlorobenzene molecules with bond lengths, and the bond length distribution obtained from the tests shows good consistency with the bond length distribution obtained from all-atom molecular dynamics , would this approach be reasonable? Thank you for your guidance.

Regarding this, such an approach would indeed be reasonable.