Unexpected breakdown of MBX-1.3 when doing NPT-PIMD sampling

[Chongbin Wang]

Dear MBX developers,

I have been using MBX-1.3 for some NPT-PIMD simulations of liquid water with our homemade MD engine. However, I occasionally encounter an unexpected breakdown of the simulation after maybe 500-1000 ps of simulation time. Have you ever encountered such issues when using MBX for PIMD simulations? If so, could you please provide some suggestions on how to avoid this issue?

The parameters I used are as follows:
- System: 256 water molecules in a cubic box
  - two body cutoff: 9.0 Å
  - three body cutoff: 4.5 Å
  - dipole iteration method: cg
- Temperature: 298.15 K
- Pressure: 101.325 kPa
- Number of beads: 24
- Time step: 1.0 fs
  - It seems that the time step is too large for PIMD, but it works well for q-spc/fw water model in NVT-PIMD simulations.
  - I have tested smaller time steps (0.5 fs and 0.25 fs), but the breakdown still occurs.
- Integrator: Unified Middle Scheme
  - If you are interested, you can find the implementation details in this paper: DOI: 10.1021/acs.jctc.5c00573
- Thermostat: Langevin
- Barostat: MTTK

I have never encountered such breakdowns when doing NPT-MD (classical MD without path integrals) simulations. I investigated the trajectory carefully and found that the system looks normal before the breakdown, with no abnormal spikes in energy or pressure. The box size also appears reasonable before the breakdown. At the breakdown point, I observe a sudden and significant increase in energy, pressure, and force magnitudes, leading to numerical instability and simulation failure.

I would greatly appreciate any insights or suggestions you may have regarding this issue. Thank you in advance for your help!

Best regards,
Chongbin Wang

[Philip Zhu]

Hello,

We do recognize that with large timesteps, simulations might occasionnaly run into unstable regions of the potential energy surface, particulally with PIMD which allows the structures to explore classically forbidden geometries, and a larger timestep may cause the structure to jump into some "potential energy holes". There are other reports of PIMD crashing but they are generally resolved by using smaller timesteps. Maybe relevant discussions:
https://groups.google.com/g/mbx-users/c/T9SZpiX5Lbc/m/q1gmx3XuCQAJ

https://groups.google.com/g/mbx-users/c/GhcmUDy7pcY/m/gVZIZQXFBAAJ

In these discussions, the system were stable under 0.1fs or 0.05fs. However, if reducing the timestep do not help, you might be seeing different issues. e.g. It might be sensitive to the thermostat or integration method.

We have a recent development that significantly improves the stability issues. (DOI: 10.1021/acs.jctc.5c00407) This model will likely be more tolerant with timestep, integration method, or any other issues, because there shouldn't exist any "holes" to jump into for most systems (unless in extremely high pressure >10GPa etc.).

[Chongbin Wang]

Thank you for your detailed response and suggestions. I will look into the references you provided and consider how to improve the stability of my simulations. Additionally, I will explore the recent development you mentioned in DOI: 10.1021/acs.jctc.5c00407. Will the new technique be available in the upcoming MBX release?

Best regards,
Chongbin Wang

[Philip Zhu]

the model from the above reference is currently available as a separate MBX (zipped) under this github repo https://github.com/paesanilab/Data_Repository/tree/main/MB-pol_PIFS