Fix NVT with two different temperatures in one crystal

Dear Lammps users,

I performed a simulation with a silicon crystal. I discretized the crystal into 12 slices. In the first slice I fixed NVT at 500K, and for the rest of slices were fixed at 120K. It’s still a perfect crystal. I’d like to see how the temperature changes along the slices. I expected to see the middle slice has temperature of around 120 K instead I found temperature of middle slice decreases from ~120K to 40K after performing 1000000MD steps with timestep of 2fs. I’m wondering if this is the phenomena I should obtain from performing two NVT in one crystal or if I’ve done something wrong with it. Does anyone have any idea about why I got this kind of results?

Here’s my input file:

group atom_h type 1

group atom_c type 2 3 4 5 6 7 8 9 10 11 12

group atom_a type 1 2 3 4 5 6 7 8 9 10 11 12

velocity atom_h create {th} 2344 velocity atom_c create {tc} 5983

fix nvt_h atom_h nvt temp {th} {th} 0.2
fix nvt_c atom_c nvt temp {tc} {tc} 0.2

compute ke_a atom_a ke/atom

variable temp_a atom c_ke_a/1.5/0.000086173303

dump nvt_dump3 atom_a custom 1000 nvt_temp_a.txt id type v_temp_a

Thanks!

Chien

Dear Lammps users,

I performed a simulation with a silicon crystal. I discretized the crystal
into 12 slices. In the first slice I fixed NVT at 500K, and for the rest of
slices were fixed at 120K. It's still a perfect crystal. I'd like to see
how the temperature changes along the slices. I expected to see the middle
slice has temperature of around 120 K instead I found temperature of middle
slice

decreases from ~120K to 40K after performing 1000000MD steps with timestep
of 2fs. I'm wondering if this is the phenomena I should obtain from
performing two NVT in one crystal or if I've done something wrong with it.
Does anyone have any idea about why I got this kind of results?

​if you want to see the distribution of temperature (or rather kinetic
energy) between a "hot​ zone" and a "cold zone" you must not do any
thermostatting in the intermediate regions.

please see examples/KAPPA/in.langevin for how this can be done (in that
case to determine thermal conductivity in a simple liquid). also, i don't
think that a nose-hoover thermostat is a good choice for this kind of
setup. NH thermostats are tuned for equilibrium systems and they do not
behave so well, if they have to continuously add or remove kinetic energy.
a dissipative thermostat is much better suited, as it can transfer kinetic
energy more easily and consistently.

axel.

Hi Axel,

Thanks for your response.

What I’d like to have is a crystal that has a temperature profile with a step change (500K and 120K). I’d like to prepare a crystal with high surface temperature rather than calculating the thermal conductivity using the examples in KAPPA. That’s why I performed two NVT with different temperature on the crystal. Is it impossible to get that temperature profile?

Best,
Chien

Hi Axel,

Thanks for your response.

What I'd like to have is a crystal that has a temperature profile with a
step change (500K and 120K). I'd like to prepare a crystal with high
surface temperature rather than calculating the thermal conductivity using
the examples in KAPPA. That's why I performed two NVT with different
temperature on the crystal. Is it impossible to get that temperature
profile?

​it is unphysical. ​

Hi Axel,

Yes, I understand it’s unphysical. I just like to know whether I can use lammps to create the system like this. It looks like it’s not possible to create a system like this. Thanks!

Chien

Hi Axel,

Yes, I understand it's unphysical. I just like to know whether I can use
lammps to create the system like this. It looks like it's not possible to
create a system like this. Thanks!

​well, it will only work by doing unphysical things, e.g. you would have to
turn off interactions between the two segments​. as soon as the segments
interact, they will couple and thus transfer kinetic energy. or you'll have
to use a so-called massive thermostat, where each degree of freedom is
thermalized independently and with a short time constant (something that
doesn't work well or at all with fix nvt). at any rate, doing something
like this makes any outcome of this questionable, as this would turn a
computer *simulation* into a compute *animation*.

axel.

I see.

Thanks a lot. I’ll look into that and think about whether I should still work on this setup.

Chien