I am working on simulating a system that includes hexadecane molecules sandwiched between two iron oxide surfaces, with the latest lammps version (2 Apr 2025). The simulation box is triclinic with no initial tilt.
I am thermostatting the system by applying a langevin thermostat on the middle region inside both iron oxide surfaces and I want to deform the simulation box in the xz plane to simulate shearing. During these shearing runs I want to monitor and output certain parameters that I need for my study, namely the mass density profile of the hexadecane molecule film, the temperature profile as well as the velocity profile.
For the density profile I used the following commands:
variable xz equal xz
region hexadecane_prism prism 0 54 0 50.5 23 40 0 v_xz 0 units box
compute density_hexadecane_bins hexa_molecules chunk/atom bin/1d z lower 0.00842 units reduced region hexadecane_prism
fix density_profile_hexa hexa_molecules ave/chunk 100 10 1000 density_hexadecane_bins density/mass file density_profile.txt
My aim with this is to have several bins with a thickness of 1 Angstrom in the z direction (as the simulation box is roughly around 119 Angstrom in the z direction) to monitor how the distribution of atoms is like inside the film as the box is deforming in the xz plane.
My question is as follows:
To monitor the temperature that is independent of the box deformation/shearing that the system is undergoing, my understanding is that I need to use compute temp/profile as is suggested in the documentation. However, how can I implement this command with the specific binning style that I used in the fix density/mass command to end up with a temperature profile similar to the mass density profile of the hexadecane film? As I want to plot the temperature distribution inside the hexadecane film and see how it varies.
The phrasing of your question is a bit confusing to me, but I think you can simplify the input to get what you want. Using a region is a bit of a stretch in your case. If you’re interested in your hexadecane film, simply using your hexa_molecules group will only store contributions from your hexadecane atoms. Spanning the whole box is not a problem then. You can use the z lower 1. units box with no region to get exactly the bin size you want.
Then you can simply make something like variable nx equal "ceil(lz)" that will give you the upper integer value just above your lz length so you’ll end up with the approximate number of bins. Applying it on the hexadecane group will compute the temperature of the group only.
As I understand what you define as “the temperature that is independent of the box deformation/shearing that the system is undergoing”, this is the biased temperature from temp/profile, right? You can write it on a regular basis using a fix ave/time command directly.
Hello @Germain, Thank you very much for your response.
Concerning this, you are right using region is not necessary, I only added it before as an added layer of “security” as I wanted to make sure that the bins are evolving dynamically with the continuous tilting of the simulation box. But after re-reading the documentation thanks to you bringing it up, I now know that the region command is not needed.
As a side note, concerning the units, if I understood correctly according to the documentation for this command, for triclinic boxes only the reduced units may be used. In this regard, the reduced units would only apply to the box dimensions for defining the size of the bins and not to other quantities correct? As in the mass density would still have the units specified originally by the units command, in my case it’s real units, or would it then be in reduced units as well?
Yes exactly, that is what I am referring to, as I want to remove the contribution from the shear induced velocity to have an accurate temperature profile for the molecules.
You are right. And this makes things even simpler as you can get the bin size by dividing the length by the number of bins you get from the above formula.
Yes, this option only concerns the size of the bins. Bins are only a label attributed to atoms so there is no reasons other quantities are affected.