Compute the rigid body temperature

Dear lammps-user:

I am trying to simulate a methane flow through two montmorillonite (MMT) plates using boundary-driven non-equilibrium molecular dynamics (BD-NEMD). Before that I would like the perform equilibrium MD. As reported in the paper, MMT is treated as rigid body. So I use two separate thermostat to control the system temperature via:

timestep 1.0
fix 1 methane nvt temp 300.0 300.0 100.0
fix 2 boundary rigid/nve group 1 boundary langevin 300.0 300.0 100.0 23234 torque * off off off

I am turning the torque off is because I don’t want the plates to rotate (which has been observed if I don’t put torque off.)

Then, I compute the temperature of two groups via:
compute methane_temp methane temp
compute boundary_temp boundary temp

What I found is a very large temperature fluctuation of boundary temperature.

Step methane_ boundary 2 boundary DOF TEMP PotEng KinEng Press Volume

78900 308.92838 355.39523 355.39523 1.0593657 1.5 355.39523 13494118 1252.5038 2.162838 174636
79050 299.44934 197.80225 197.80225 0.58961093 1.5 197.80225 13494111 1213.6352 -1258.8009 174636
79200 296.72578 902.29643 902.29643 2.6895743 1.5 902.29643 13494096 1204.7022 -1671.1931 174636
79350 298.7758 322.56075 322.56075 0.96149231 1.5 322.56075 13494107 1211.2786 -518.73086 174636
79500 302.69835 65.704396 65.704396 0.19585232 1.5 65.704396 13494121 1226.4029 323.93217 174636
79650 298.29063 106.95878 106.95878 0.3188238 1.5 106.95878 13494109 1208.6706 1175.3383 174636
79800 289.91067 155.88928 155.88928 0.46467635 1.5 155.88928 13494098 1174.8699 -730.99081 174636
79950 305.51408 97.645068 97.645068 0.2910614 1.5 97.645068 13494101 1237.9044 -2336.6775 174636
80100 298.78331 57.395392 57.395392 0.17108476 1.5 57.395392 13494099 1210.5186 480.76507 174636
80250 301.56259 695.2837 695.2837 2.0725086 1.5 695.2837 13494097 1223.6787 -3544.51 174636
80400 318.11962 35.435108 35.435108 0.10562533 1.5 35.435108 13494103 1288.783 -1470.2342 174636
80550 293.26758 92.604261 92.604261 0.27603571 1.5 92.604261 13494137 1188.2798 -711.20646 174636
80700 300.13348 42.105384 42.105384 0.12550815 1.5 42.105384 13494116 1215.9425 -338.41974 174636

I am not sure if my temperature computation for the rigid body is correct. I am looking into the code. There is another available “compute” which is to compute the kinetic energy of rigid body (compute ke/rigid). I then check this value and recalculate the kinetic-energy-based temp. I found that if DOF is 1.5, the computed temp is same as the original compute and fix output. I am not sure if this type of temperature computation is reasonable. Few paper reported the temperature separately for their “rigid” body. Has any one encountered similar problem? Any help is appreciated. Thanks in advance. I will put my input script and data structure as attached if anyone is interested this problem.

Bests,
Shuai

in.large_external_MMTmethane (2.28 KB)

MMT_methane_10_external.data (119 KB)

Dear lammps-user:

I am trying to simulate a methane flow through two montmorillonite (MMT)
plates using boundary-driven non-equilibrium molecular dynamics (BD-NEMD).
Before that I would like the perform equilibrium MD. As reported in the
paper, MMT is treated as rigid body. So I use two separate thermostat to
control the system temperature via:

timestep 1.0
fix 1 methane nvt temp 300.0 300.0 100.0
fix 2 boundary rigid/nve group 1 boundary langevin 300.0 300.0 100.0 23234
torque * off off off

I am turning the torque off is because I don't want the plates to rotate
(which has been observed if I don't put torque off.)

Then, I compute the temperature of two groups via:
compute methane_temp methane temp
compute boundary_temp boundary temp

What I found is a very large temperature fluctuation of boundary
temperature.

​which is not at all surprising. due to making those rigid and removing
rotation, you have only 3 translations left as DOFs. computing a
temperature for such a system makes practically no sense. more importantly,
due to keeping those "boundary objects" rigid, what your methane molecules
are "experiencing"​ is practically indistinguishable from a 0K system. so
why not make them completely immobile and save yourself a lot of trouble?

axel.

Hi Axel,

Thanks for the prompt reply. Why I am doing in this way is because I would like to add an external force field at the boundary to mimic the flow system. Some previous work added the thermostat on the fluid but a lot of works state that the fluid particle should not be perturbed by the thermostat and add the thermostat on “wall” atoms to help remove the heat generated by the external force. If I freeze or tether the wall atoms to its original positions, I do find that the heat cannot be removed and the temperature of methane increases. That’s why I tried “rigid” body. May I ask, if the system only contains 1 large rigid body, does it mean that the thermostating on such rigid body may be meaningless? I am really new to this area, thanks for your help.

Shuai

Hi Axel,

Thanks for the prompt reply. Why I am doing in this way is because I would
like to add an external force field at the boundary to mimic the flow
system. Some previous work added the thermostat on the fluid but a lot of
works state that the fluid particle should not be perturbed by the
thermostat and add the thermostat on "wall" atoms to help remove the heat
generated by the external force. If I freeze or tether the wall atoms to
its original positions, I do find that the heat cannot be removed and the
temperature of methane increases. That's why I tried "rigid" body. May I
ask, if the system only contains 1 large rigid body, does it mean that the
thermostating on such rigid body may be meaningless? I am really new to
this area, thanks for your help.

​yes. just look into a text book on (statistical) thermodynamics​.
temperature is technically only a well defined property in the limit of a
macroscopic system, which means a practically infinite number of particles
(or degrees of freedom). under the assumption of equipartitioning, one can
relate the macroscopic temperature to the microscopic kinetic, but you need
to either average over a large enough number of degrees of freedom or over
a long enough time (which in turn requires that the system is in
equilibrium) to have this assumption to be reasonably well fulfilled.

​one possible approach to immobilize particles and still enable them to
​dissipate kinetic energy, would be to use position restraints (via fix
spring/self) and then apply a dissipative thermostat (fix langevin or fix
csld or fix gle or fix gld) to those particles.

a second option would be to make your walls have multiple layers/zones and
have an immobile zone (OK), on top of that a thermostatted zone, and on top
of that a zone that is not thermostatted. this top zone would than be
interacting with your methane molecules. this is the most complex, but also
the most realistic approach. it is commonly used for slab calculations when
you do have the surface interact with something.

axel.

1 Like