I am working with a system that includes a flexible MOF framework and rigid water molecules, using the Drude oscillator approach. The atoms are divided into three groups:
“ATOMS_mof”: includes non-polarizable atoms and Drude cores (DC) in the MOF,
“ATOMS_w”: includes non-polarizable atoms (Hw, Ow, M) in SWM4-NDP water molecules,
“DRUDES”: includes all Drude particles (DP) for both MOF and water (i.e. Ow_DP).
To simulate the system under NPT conditions, I applied the following fixes:
fix LANG all langevin/drude 300. 100. 12435 1. 20 13977
fix RIGID ATOMS_w rigid/nve/small molecule
fix NVE DRUDES nve
fix NPH ATOMS_mof nph iso 1. 1. 500
However, I noticed that this setup applies an NVT-like integration to the rigid water (rigid/nve/small) and an NPT barostat to the MOF atoms. This seems to cause a sudden temperature increase at the start of the simulation.
From what I understand, it is not possible to apply two separate NPT fixes to interacting subsystems (e.g., MOF and water), as this could lead to conflicting box scaling and instabilities.
So I am posting this question to ask what is the recommended way to handle such a Drude system that includes both a rigid water model and a flexible MOF under NPT conditions?
It is not correct to apply an external thermostat to atoms within a rigid fix. Those require thermostats for the rigid objects, not individual atoms, so you would probably want to use fix rigid/nvt/small
I would first run for a bit without a barostat. (i.e. use fix nve instead of fix nph) until the system has equilibrated and stabilized somewhat. Whether you see a change in the temperature and/or kinetic energy of the run depends on how and how well the system was equilibrated before and how much the potential energy was changed in between the two runs.
As a superset of @akohlmey 's advice, if you don’t know for sure that you need a barostat, then you should try running without one. Monitor the average pressure of your system and, if it stays close to your desired pressure, you can present that in the SI and chant “thermodynamic equivalence of ensembles, thermodynamic equivalence of ensembles” until your reviewers go away.
(The bit about the reviewers is is meant to be a joke, but the rest of it is sincere advice.)
I think I need to explain that I am using the SWM4-NDP water model, as the manual states:
NVT ensemble using Langevin thermostat to rigid bodies:
comm_modify vel yes
fix LANG all langevin/drude 300. 100 12435 1. 20 13977
fix RIGID ATOMS rigid/nve/small molecule
fix NVE DRUDES nve
I did the same; the non-polarizable H and Drude core for Ow were grouped into ATOMS_w, and I have kept it rigid here, moving as a rigid object, while allowing the Drude particle(DP) for Ow to be flexible and be thermalized to 1 K.
I tried with this advice, the temperatures started being abnormal from the 1800 step, and finally the bond atoms missing error occurred. I have attached the input and log here. nvt.in (1.9 KB) 1.data (641.3 KB) 1.log (13.8 KB)
This does not make sense when I look at the implementation of the rigid fixes. All rigid fix styles maintain and update, center of mass position and velocity of each rigid body as well as orientation and rotational velocity. The per-atom velocities are recreated from those in every step, so any velocity updates from a langevin fix will be either ignored or inconsistent with the rigid body motion.
Fix rigid has a langevin option to apply a langevin thermostat to those rigid body DOFs.