Temperature profile problems


I am trying to calculate the temperature profile across my simulation box, which consists of liquid confined within two walls on the top and bottom. The left and right boundaries are periodic. The two walls have different temperatures fixed to them so one would expect the temperature profile to be a straight line (after a sufficient amount of time).

To calculate the temperature profile I use the following commands.

compute ke flow ke/atom
variable temp_profile atom c_ke/1.5
fix c9n flow ave/spatial 1 200000 200000 y lower 0.025 density/number vx v_temp_profile units reduced file temp_profile.dat.

However, the temperature of the liquid is not only far from a straight line, but it is also not within the temperatures I fixed for the walls (deviating by a big margin). I also went as far as to fix the temperature of the liquid itself but that also doesn’t produce reasonable results. Has anyone stumbled on similar problems? Is there something obviously wrong with the way I am calculating the temperature?

Thanks in advance for any help

This is only true if you’re in the linear response regime of the liquid. Could we have a look at your input files and an example temperature profile? Also, what do you mean that fixing the temperature of the liquid doesn’t give reasonable results? How are you fixing it and why isn’t it working?

Dear Niall,

Thank you very much for your response. Unfortunately I am on a different computer at the moment which does not have the temperature profile. However, I am attaching the input file (it might be slightly out of date but more or less the same) and I should also be able to send you the temperature profile later today.

Regarding the fix, i used the command (which is commented out in the input file I am attaching)

fix 6 flow temp/rescale 1 1.0 2.0 0.02 1.0 # Re-adjust the temperature of the fluid every time-step

to keep it within the required temperatures. However, the temperature profiles show a temperature of over 3.0.

Thanks again for your help

flow.in (2.45 KB)

Dear Michael,

Your ave/spatial command is averaging from the very start of the simulation (there's no equilibration period). I assume you want the temperature profile at the non-equilibrium steady state, so you've got to allow the system some time to get there before you start averaging the temperature. Equilibrate your system first, by running for maybe 100000 timesteps (to start with, at least). You can look at things like energy conservation and the radial distribution function of the fluid component to see when the fluid has reached a steady state. You can then create a restart file and run a new simulation from that point to calculate your temperature profile.

I also don't think you need pair_style hybrid. Because all of the particles interact via the Lennard Jones potential, you can just use:
- pair_style lj/cut 2.2
- pair_coeff 1 1 1.0 1.0
- pair_coeff 1 2 1.0 1.0
- pair_coeff 2 2 0.0 0.0

You can also use "neigh_modify exclude type 2 2" to make the neighborlist ignore the 2-2 interactions. However, I'm quite new to LAMMPS, so it's probably best if somebody else can corroborate that!


Dear Niall,

Since I am running it for 1200000 timesteps (which I assume is enough), shouldn’t the profiles produced by ave/spatial during the last timesteps be correct anyway? Do I have to start taking the averages after it has reached steady state? I am not expecting the temperature profiles during the first timesteps to be any good!

fix 6 flow temp/rescale 1 1.0 2.0 0.02 1.0 # Re-adjust the temperature
of the fluid every time-step

If you're doing this, it will mess up the linear temperature profile
you are hoping
to establish.


Hi Steve and thanks for your help.

I know it messes the temperature profiles. That fix is actually commented out. I only used it at some point hoping that it would shed some light on my problem (which it didn’t).

Note that the ave/spatial is calculating average kinetic energy not
temperature. Only in the unit system where kboltzman=1 they will be
identical. One of your problems could be units related. Try running a
bulk simulation of the liquid and compute the temperature with the fix
as well as with the default thermo command. That should show you if
there is any problems with the units while help you realizing the unit
factor connecting energy and temperature.
Related to creating the temp gradient, why bothered using two walls?
Just divide the box in that direction into bins and place two Langevin
thermostats at different temperatures separated by the appropriate
distance. Yet, as Nail said make sure your gradient is small enough so
that your liquid will be in the linear regime.

Sorry, that's my mistake- I thought ave/spatial defaulted to the "running average" mode. I would try Carlos's suggestion to see if that helps, and also possibly make the ave/spatial bins a little larger (you'll lose some spatial resolution, but gain more particles contributing to the average in each one).

Compute temp/profile will calculate the actual temperature profile
on a per-bin basis if you give it the out bin option. You can
then output this via fix ave/time (w/out averaging if desired),
and plot it.