I am trying to apply a directional thermostat (partial temperature control) to a group of rigid molecules. My goal is to thermostat the system only in the y and z directions, leaving the x direction unthermostatted.
I attempted to achieve this using compute temp/partial in combination with fix rigid/nve/small and fix langevin. My commands are as follows:
compute 4 h2 temp/partial 0 1 1
fix 5 h2 rigid/nve/small molecule
fix 6 h2 langevin 323 323 100.0 48279
fix_modify 6 temp 4
I also tried a similar setup using fix rigid/nvt directly (modifying its temperature compute). In both cases, the simulations run smoothly without throwing any warnings or errors.
However, because rigid body integration handles translational and rotational degrees of freedom differently, I am not entirely sure if the partial thermostat is truly taking effect as intended on the rigid bodies.
Does LAMMPS physically support direction-dependent (partial) temperature control for rigid bodies?Is my current approach (using fix_modify to pass temp/partial to fix langevin or fix rigid/nvt) the correct way to implement this?
Any insights or suggestions would be greatly appreciated. Thank you for your time!
I can say for certain that a combination of a rigid fix and a per-atom thermostat like fix langevin cannot work. The rigid fix maintains an internal state of the rigid bodies and that only uses the per-atom forces as input. The per-atom coordinates and per-atom velocities are in every step regenerated from the rigid body positions, orientations, center of mass velocities, and orientational velocities. Only when the initial rigid body information is assembled, will the per-atom coordinates and velocities be used. Please see the documentation for confirmation.
As for applying fix_modify, I suggest you make a careful test. I recommend to make a test simulation where the effect is more visible. For example, have an initial system with all velocities set to zero and then use compute temp/region on a subset of the system, so you should have part of your system immobile and part moving.
Please note that the plain rigid fixes have a langevin keyword to apply a langevin thermostat to translation and rotation of rigid bodies similar to how fix langevin operates on atoms.
Thank you so much for the detailed explanation and the definitive answer regarding the conflict between fix rigid and fix langevin.
Following your advice about making a careful test, I will set up an extreme test case with zero (or near-zero) initial velocity to see the physical effect. I plan to test if the partial temperature can be correctly passed to the rigid NVT integrator using this corrected syntax:
After verifying this, I will also implement and test your recommended approach using the built-in langevin keyword within the plain fix rigid command to compare the trajectories and thermodynamic behavior.
I also have a related follow-up question about my specific H2 model. I am using a 3-site H2 model with a virtual center of mass (M) for electrostatics, so the H-M-H angle is exactly 180 degrees. Is it true that such a linear 3-site model can only be constrained using fix rigid? I understand that fix shake might encounter mathematical singularities or become numerically unstable when trying to constrain an angle of exactly 180 degrees. Could you confirm if fix rigid (or fix rigid/small) is my only robust option for maintaining the geometry of this specific linear model?
Thank you again for your time and invaluable insights!
fix rigid/nvt combined with fix_modify temp/partial 0 1 1.
The built-in langevin keyword within fix rigid/small.
In both cases, the unthermostatted X-direction rapidly heated up to the target temperature. It is completely clear to me now that the rigid body integrator inherently couples rotational and translational degrees of freedom, making it physically impossible to decouple them using fix_modify.
This leads to my ultimate dilemma. I am currently simulating fluid mixtures (including this 3-site linear H2 model) flowing through a nano-slit. In such Non-Equilibrium MD (NEMD) flow simulations, it is standard practice to thermostat only the directions perpendicular to the flow (e.g., Y and Z) to avoid biasing or damping the natural flow velocity profile in the X-direction.
Given the two hard constraints:
fix rigid inherently does not support directional thermostatting.
fix shake mathematically diverges for a strictly linear (180.0 degree) 3-site model.
What is the recommended best practice in LAMMPS to achieve directional temperature control for such rigid linear molecules in a flow simulation? Would introducing a tiny, artificial bend (e.g., setting the angle to 179.9 degrees) to safely bypass the SHAKE singularity—thus allowing me to use fix shake coupled with a standard fix langevin + fix_modify—be considered a acceptable compromise ? Or is there another advanced technique I am missing?
The more I think about this, the more I feel that this is not such a good idea for two reasons.
as you have already noticed, you cannot decouple translation and rotation and thus there will always be an energy transfer. This would not be restricted to rigid bodies, but basically any molecule and thus using a compute temp/partial bias only makes sense to me for atomic systems.
I am not convinced that when simulating a flow through a pore, there is a justification for any kind of thermostat. What would be the physical reasoning where the kinetic energy that is exchanged comes from. For a Nose-Hoover thermostat, the justification is that you have only a small bulk(!) system in your simulation and you want it to sample the phase space of a much larger system and thus you couple weakly(!) to a reservoir. For a Langevin thermostat, the usual reasoning is that you use it to have the effect of an implicit solvent. In many non-equilibrium simulations it is also used to just efficiently remove excess energy that it added to the system.
Thus my conclusion from that is that with a directional thermostat, you would be taking a convenient, but not necessarily realistic, shortcut. I think a better approach would be to quantify the kinetic energy added from the induced flow and then just use a Nose-Hoover thermostat to thermostat to a suitably adjusted target temperature. Nose-Hoover is a global thermostat and thus should not result in a significant bias on the system, if the reservoirs should be large enough. At the very least I would expect the effect not to be worse than that of applying a bias from compute temp/partial due to the coupling between directions from the bonds/bodies.
I don’t recall that this topic has ever come up. But keep in mind that my focus is on code development and my knowledge of the published literature and thus the overall consensus on what should be done for what application is limited to cases that I am exposed to for one reason or another. There are essentially two kind of issues here:
what is the best way to thermostat such a system from the statistical mechanical point of view?
what can be done in the code to address its current limitations (e.g. not being able to use fix shake on linear molecules)?
For point 1. you need to talk to somebody else. For point 2. I could see a workaround in internally projecting any linear molecule to a diatomic where the two atoms are chosen to closely represent the moments of inertia of the lineal molecule, then apply the SHAKE constraint and then project the resulting constraint forces/velocities back to the constituent atoms, very similar to how fix rigid computes positions and velocities of its constituent atoms. This is a non-trivial programming project and also requires some careful math beforehand to make sure that my guesstimates are actually true.
While this may be standard practice, partial thermostats introduce inaccuracies when used on rigid molecules with long-range ordering (such as from the density fluctuations in a slit), since the remaining degrees of freedom after rigidification are not equally partitioned between x, y, and z axes: https://pubs.acs.org/doi/10.1021/acs.jctc.4c00957
There are complicated temperature computes available such as compute temp/profile command — LAMMPS documentation to calculate temperature commensurate with NEMD flow profiles. But honestly, if you have walls, the most defensible choice is to simply thermostat the walls (with something like fix spring/self to stop them from drifting) while omitting any thermostatting whatsoever on the gas particles.
As the NEMD simulation progresses it should reach a steady state in which the internal energy of the flowing gas is constant on average and the average power input from the NEMD force is equal to the average thermal power dissipation through the wall thermostat, and at that point the power source, the gas, and the wall must be at the same “temperature” (for whatever definition of temperature makes sense in a flagrantly non-equilibrium situation!).