I am attempting to simulate a water model with a TIP4P-style geometry, and this model has Drude oscillators attached to the hydrogens and M-site (reference for model - A systematic development of a polarizable potential of water | The Journal of Chemical Physics | AIP Publishing). My simulations using the Extended Lagrangian approach in LAMMPS are showing signs of violating equipartition as described by Son et. al. for Dual Langevin or NH thermostats, probably due to the complexity of this model compared to SWM4 water (reference for Son et. al. investigation of Drude thermostatting - http://dx.doi.org/10.1021/acs.jpclett.9b02983). It seems as though the temperature-grouped thermostat will solve these issues, but I am unsure of how to apply it while keeping my water geometry rigid. Due to the water having a 4-site geometry, it cannot be constrained by SHAKE to my understanding. Also, fix rigid/small seems to apply NVE integration and thus layering the temperature-grouped thermostat on top would integrate atomic positions twice. Is there a way using the current features of LAMMPS to implement this? Thank you so much!
Son’s temperature-grouped thermostat (if correctly implemented – but I would assume it is unless proven otherwise) can’t really work for your system. As the documentation says, it works by grouping the thermal kinetic energies of the system into the molecular COMs, intramolecular atomic COMs, and Drude oscillator COMs.
But a rigid molecule has no kinetic energy in its intramolecular atomic COMs, so I don’t see how a temperature-grouped thermostat for your system would work differently from the suggestion in the Drude Tutorial:
fix DIRECT all drude/transform/direct
fix RIGID ATOMS rigid/nvt/small molecule temp 300. 300. 100
fix NVT DRUDES nvt temp 1. 1. 20
fix INVERSE all drude/transform/inverse
Also – in case you need further help, please give further description of your simulations, especially how you have quantitatively measured “equipartition failure”. Equipartition in rigid systems can be very subtle. You may want to post your script here, “fencing” it with triple backticks
```
like this
```
to enable verbatim formatting – but if it is complex and not well annotated (which happens to all of us!) please include comments as needed to further explain what’s going on.
By the way – fix rigid/small only modifies the equations of motion to rigidify molecules. The integration is done separately by fix nvt (or similar). Only the integrator-included versions, such as fix rigid/small/nvt, perform time integration as well as modifying the equations of motion.
See below. I remembered wrongly – both fix rigid/small and fix rigid/nvt/small do time integration; the latter just uses a different integrator scheme.
I now see that perhaps the temperature-grouped thermostat is not a promising solution due to the rigid water frame, compared to the problematic models for acetone and ionic liquids examined by Son where only C-H bonds were rigid. However, I was looking into this because my system was displaying the same signs of ‘equipartition violation’ as those reported in Fig. 2 of Son et. al. Meaning that applying a Langevin dual thermostat on my system overestimated the bulk density, achieved the target core atom temperature ~300K, but never achieved the intended Drude temperature (hovered around ~1.5K rather than 1K). The dual Nose-Hoover thermostat application leads to an underestimation of the bulk density, doesn’t achieve core temp (~320K rather than ~300K), but achieves the intended Drude temperature of 1K. This is exactly what Son et. al. observes for equipartition failures in molecules with high-frequency intramolecular modes, which my model doesn’t have. I would then like to modify my question to if you have any intuition as to why I am getting similar failures in my different scenario. It is also worth noting that in order to achieve a better separation of time scales for the light hydrogen-drude pairs, I am equalizing my hydrogen and oxygen masses (I am only interested in thermodynamic properties, not transport).
I apologize for not including my input script, I have pasted in my input script (including commands for
both thermostatting techniques) and attached my single molecule data file. Thank you so much for taking the time to respond, and I apologize for my initial lack of clarity -
# BK3 Water Model
units real
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
pair_style buck/gauss/long 10.0
special_bonds lj/coul 0.0 0.0 0.0
read_data water.dat
replicate 8 8 8
reset_timestep 0
# FF parameters for custom buckingham + gaussian charge potential
# 1 = Oxygen, 2 = M-site, 3 = Hydrogen, 4 = M-site Drude, 5 = H Drude
pair_coeff 1 1 76959.847 0.2808988764 793.499 10.0 0.0 # A, B, C, rcut, alp_ij
pair_coeff 1 2 0.0 0.01 0 10.0 0.0
pair_coeff 1 3 0.0 0.01 0 10.0 0.0
pair_coeff 1 4 0.0 0.01 0 10.0 0.0
pair_coeff 1 5 0.0 0.01 0 10.0 0.0
pair_coeff 2 2 0.0 0.01 0 10.0 0.0
pair_coeff 2 3 0.0 0.01 0 10.0 0.0
pair_coeff 2 4 0.0 0.01 0 10.0 0.0
pair_coeff 2 5 0.0 0.01 0 10.0 0.0
pair_coeff 3 3 0.0 0.01 0 10.0 0.0
pair_coeff 3 4 0.0 0.01 0 10.0 0.0
pair_coeff 3 5 0.0 0.01 0 10.0 0.0
pair_coeff 4 4 0.0 0.01 0 10.0 0.70422
pair_coeff 4 5 0.0 0.01 0 10.0 0.86769
pair_coeff 5 5 0.0 0.01 0 10.0 1.25
kspace_style ewald 1.0e-6
pair_modify tail yes
group ATOMS type 1:3
group CORES type 2 3
group DRUDES type 4 5
variable dt equal 0.5
variable N_step equal 1000000
variable TK equal 298.15 # core temperature
variable TDK equal 1.0 # drude temperature
variable PBAR equal 1.0 # system pressure
variable TD equal ${dt}*100 # thermostat damping parameter for cores
variable TDD equal ${dt}*10 # thermostat damping parameter for drudes
variable PD equal ${dt}*1000 # barostat damping parameter
neighbor 1 bin
neigh_modify every 1 delay 1
timestep ${dt}
fix DRUDE all drude N C C D D
velocity ATOMS create ${TK} 12345
velocity DRUDES create ${TDK} 12345
delete_bonds ATOMS multi
comm_modify vel yes
compute TATOM ATOMS temp
compute TEMP all temp/drude
dump 1 all custom 10 dump.coord id type x y z
dump_modify 1 sort id format line "%5d %3d %13.8f %13.8f %13.8f" flush yes
thermo_style custom step c_TEMP[1] c_TEMP[2] pe density vol
thermo 100
# OPTION 1: Dual-Langevin Thermostat, Nose-Hoover Barostat
fix 1 all langevin/drude ${TK} ${TD} 123459 ${TDK} ${TDD} 123456 zero yes
fix 2 ATOMS rigid/nph/small molecule iso ${PBAR} ${PBAR} ${PD#}
fix 3 DRUDES nve
# OPTION 2: Dual-Nose-Hoover Thermostat, Nose-Hoover Barostat
fix DIRECT all drude/transform/direct
fix RIGID ATOMS rigid/npt/small molecule temp ${TK} ${TK} ${TD} iso ${PBAR} ${PBAR} ${PD}
fix_modify RIGID temp TATOM press thermo_press
fix NVT DRUDES nvt temp ${TDK} ${TDK} ${TDD}
fix INVERSE all drude/transform/inverse
run ${N_step}
write_data out.lammps
The data file has atom types 2 and 3 as main atoms with significant mass but your comments and pair_coeff have type 1 as oxygen
The data file Bond Coeffs section labels the 0.9750 angstrom bond as O-M and 0.2661 angstrom bond as O-H, they’re probably swapped.
I’m not sure I recognize the pair style buck/gauss/long. If that’s been self-constructed, you need to verify that its code is correct. In particular, incorrect virial calculations together with npt integrators will jeopardize your simulation correctness.
You’ve got the M-site as a massive site and negligible mass at the oxygen’s position? I’m really not sure you can have such a large reorganization and still keep the model’s structural properties. At the very minimum, you should go back and work with mass-16 oxygens and mass-1 hydrogens (and a smaller time step if needed) and verify that the mass redistribution is sensible. I know protein people sometimes use mass-4 hydrogens – but that’s when they’re oversimulating massive proteins, not pure water boxes.
You are starting from an extremely highly-ordered arrangement with a density of about 500 kg/m^3, hoping to get to a final density of 1000 kg/m^3. That’s a big jump of equilibration, isn’t it? You could try the create_atoms mol option to create an initial random configuration much closer to your final density – maybe even create it at the correct density and use nvt integrators, to make sure it’s not a barostatting problem.
Heat flow from the real degrees of freedom to the Drude degrees of freedom leads to a steady temperature gradient and puts the system at an incorrect effective temperature.
so I could totally handwave-handwave and say: you’re starting with a highly-ordered system whose initial volume is 8 times its equilibrium volume with high cohesive energies at equilibrium (given both water’s fixed dipole moment and the non-additivity of polarizable interactions). That’s a lot of heat and energy to disperse through baro- and thermostatting, so that excess heat is accumulating in the Drude DOFs before it can be thermostatted away.
But I really don’t know. I hope the points above give you some food for thought though. All the best!