[lammps-users] help with temp/region and extra option

Hi,

I am trying to assign different temperatures for different regions in a 3D simulation box. The temperature calculations should not include any drifting velocity. I am having some troubles to make it work. Would you please help?

Here is what I did,

#assign region:

region 25o cylinder z 0 0 25 -40 40 side out units box
region 26i cylinder z 0 0 26 -40 40 side in units box
region 26o cylinder z 0 0 26 -40 40 side out units box
region 26 intersect 2 26i 25o units box

region 27i cylinder z 0 0 27 -40 40 side in units box
region 27o cylinder z 0 0 27 -40 40 side out units box
region 27 intersect 2 27i 26o units box

region 28i cylinder z 0 0 28 -40 40 side in units box
region 28o cylinder z 0 0 28 -40 40 side out units box
region 28 intersect 2 28i 27o units box

region 29i cylinder z 0 0 29 -40 40 side in units box
region 29 intersect 2 29i 28o units box

#temperature calculations:

fix d26 drop1 temp/rescale 1 1.3 1.3 0.01 1 region 26 # assign temperatures
compute_modify d26_temp extra 3 #take off drifting velocities
compute cr26 drop1 temp/region 26 #checking temperature with compute command
compute_modify cr26 extra 3 #take off drifting velocities

fix d27 drop1 temp/rescale 1 1.4 1.4 0.01 1 region 27
compute_modify d27_temp extra 3
compute cr27 drop1 temp/region 27
compute_modify cr27 extra 3

fix d28 drop1 temp/rescale 1 1.5 1.5 0.01 1 region 28
compute_modify d28_temp extra 3
compute cr28 drop1 temp/region 28
compute_modify cr28 extra 3

fix d29 drop1 temp/rescale 1 1.6 1.6 0.01 1 region 29
compute_modify d29_temp extra 3
compute cr29 drop1 temp/region 29
compute_modify cr29 extra 3

#output results

thermo_style custom step c_d26_temp c_d27_temp c_d28_temp c_d29_temp c_cr26 c_cr27 c_cr28 c_cr29

AND here are the results,
1495000 1.6 1.6 1.6 1.6 1.3736374 1.3389688 1.3441566 1.3349029
1500001 1.6 1.6 1.6 1.6 1.3743383 1.3437465 1.3394622 1.3441895

I thought,
1) _d26_temp c_d27_temp c_d28_temp c_d29 should be about 1.3, 1.4, 1.5 and 1.6
2)c_d2*_temp should have the same value as c_cr2* if they are in the same region.

What did I do wrong? Any help is highly appreciated.

Thanks,
Shihai

So here is the issue. When a temp is used by a thermostat
it is evaluated in the middle of the timestep. Then the velocities
are adjusted by the thermostat, but the temp is not re-evaluated.

At the end of the step, the thermo output evaluates temp computes.
However the thermostat one knows it has already been evaluated
on this timestep, so it isn't re-evaluated. Instead the old value
is used. Thus another temp for the same region or atoms will
show a different (current) temperature.

I'm of mixed opinion whether this is a bug or a feature. We could
reset the timestep flag for the compute temp so it would be
re-evaluated by thermo. But then there would be no way to output
what the temp was before the thermostat changed the velocities.
You can always define another compute temp to use with thermo
output which is what you did. If we leave it this way, we should
document that it stores the old temperature, so it won't be confusing
as it was for you (and me).

Steve

Steve,

Thanks for you help . Just to help to make myself clear about the temperatures out puts:

For example, in my case,

"c_d26_temp" is the temperature output for region 26 evaluated by "fix d26 drop1 temp/rescale 1 1.3 1.3 0.01 1 region 26 & compute_modify d26_temp extra 3 ", which is BEFORE the velocity adjustment.

"c_cr26" is the temperature output for region 26 computed by "compute cr26 drop1 temp/region 26 & compute_modify cr26 extra 3", which is the AFTER the velocity adjustment.

Am I right? I agree this explains why the temperature outputs from "c_d26_temp" and "c_cr26" are different.

However I am still confused by fact that I didn't see a temperature difference in different regions. Meaning I used,

"
fix d26 drop1 temp/rescale 1 1.3 1.3 0.01 1 region 26
  fix d27 drop1 temp/rescale 1 1.4 1.4 0.01 1 region 27
fix d28 drop1 temp/rescale 1 1.5 1.5 0.01 1 region 28
ix d29 drop1 temp/rescale 1 1.6 1.6 0.01 1 region 29
"
to set different temperatures for different regions, and the output showed the same temperature regardless what region it is computed by thermostat after many many time steps (1500000),

command = "c_d26_temp c_d27_temp c_d28_temp c_d29_temp " value = "1.6 1.6 1.6 1.6"

and from computing command "c_cr26 c_cr27 c_cr28 c_cr29" value="1.3743383 1.3437465 1.3394622 1.3441895"

It appears to me there is some kind of bugs in the " fix .... temp/rescale.... region.." command. The LAST" fix .... temp/rescale.... region.." command overrides the all previous " fix .... temp/rescale.... region.." commands and all the temperatures in different regions are rescaled to the same temperature setted in the last command which is 1.6 in my case based on the output from "c_d26_temp c_d27_temp c_d28_temp c_d29_temp " .

I am also confused by the values shown from "c_cr26 c_cr27 c_cr28 c_cr29" calculation, it seems the difference between the setting temperature (1.3, 1.4, 1.5 1.6) and actually values (1.3743383 1.3437465 1.3394622 1.3441895) is too big. I was using 0.01 for criterial.

On a separate question about the extra option in compute_modify command, for example if the number of atoms is 100 initially and changes during simulation in region 26, should I use "compute_modify d26_temp extra 3" or "compute_modify d26_temp dynamic yes extra 300"?

Thanks again for you help. Any commends are highly appreciated.

Best regards,
Shihai

I looked at this some more. Two points.
Your input script has commands that look like they
are for an older version of LAMMPS. So try everything with
the current version. Second, what I said in my previous email
about the temperature computed in the fix being out-of-sync
with what is printed at the end of step, was wrong. I assumed
it from looking at your output, but didn't look at the code.

The code appears to do the right thing. This example script
(much simpler than yours, modified version of bench/in.lj)
works just as I would expect for the temp of the 3 regions.

# 3d Lennard-Jones melt

units lj
atom_style atomic

lattice fcc 0.8442
region box block 0 20 0 20 0 20
create_box 1 box
create_atoms 1 box
mass 1 1.0

velocity all create 1.44 87287 loop geom

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5

neighbor 0.3 bin
neigh_modify delay 0 every 20 check no

fix 1 all nve

region box1 block 0 5 0 5 0 5
region box2 block 10 15 10 15 10 15
region box3 block 15 20 15 20 15 20

compute temp1 all temp/region box1
compute temp2 all temp/region box2
compute temp3 all temp/region box3

fix scale1 all temp/rescale 10 2.0 2.0 0.5 1.0
fix scale2 all temp/rescale 10 2.5 2.5 0.5 1.0
fix scale3 all temp/rescale 10 3.0 3.0 0.5 1.0

fix_modify scale1 temp temp1
fix_modify scale2 temp temp2
fix_modify scale3 temp temp3

compute temp1_final all temp/region box1
compute temp2_final all temp/region box2
compute temp3_final all temp/region box3

thermo_style custom step temp pe press &
    c_temp1 c_temp1_final &
    c_temp2 c_temp2_final &
    c_temp3 c_temp3_final

thermo 10

run 100

So unless you can convince me otherwise with the current
version, I think there is no bug.

Steve