Simulating a rigid, four-site, polarizable water model with a temperature-grouped thermostat


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 - 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.


Ahhh! Thanks. Documentation – gotta read it.

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

water.dat (1.3 KB)

There are some issues:

  • You need to double-check some labels:
    • 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.

I mean, the abstract of Son’s paper says:

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!