how to obtain right region temperature for the system constrained by fix shake

[111190]

Dear all,

I’m trapped in the calculation of region temperature. The version of 20160216 lammps is employed. Two methods are used in my simulation about inorganic aqueous solution. One is by the command of compute temp/region and the other is by the command of compute ke/atom. The specific procedures are as follows:

Both methods give the consistent temperatures in hot or cold region. But it is very confusing that the calculated temperature is wrong. For example, if the system is kept at 300K by Nose-Hoover thermostat, the region temperature will be about 212K, which is quite unreasonable.

In order to find out what should account for the wrong region temperature, I simulated pure SPC/E water maintained at 300K by ose-Hoover thermostat under NVT ensemble. Two conditions are simulated. A part of the log file is showed :

Results of the first condition using special_bonds lj/coul 0 0 0.5

Results of the second condition using default special_bonds lj/coul 0 0 0

From results above, mainly two problems exist.

Firstly, from http://lammps.sandia.gov/doc/compute_ke.html, thermo_ke corresponds to thermo_temp and is different from compute_ke. But in the first condition, thermo_ke doesn’t agree with thermo_temp. For both conditions, thermo_ke equals compute_ke completely.

Secondly, from http://lammps.sandia.gov/doc/compute_temp_region.html, compute temp/region differs from other compute temp because it does not subtract out degrees-of-freedom due to fixes that constrain motion, such as fix shake. In the first condition, temp/region and ke give wrong overall temperature and ebond as well as eangle equals 0. In the second condition, temp/region and ke give right overall temperature but ebond and eangle are not 0. So I doubted whether both ebond and eangle should be 0 when fix shake is used to constrain bonds and angles. The only difference of inputfile in two conditions above is the special_bonds settings. In the first condition, I think the auto-generation special neigbors in lammps is wrong. From auto-generation in logfile, when using lj/coul 0 0 0.5, 1-4 neigbors are 1 for water but it is 0 actually.

So, now, I have no idea how to get the right region temperature under the condition when fix shake is used, ebond and eangle are always 0 regardless of any special_bonds settings，which happens in my aqueous solution simultion.

Can anyone tell me how to solve the problem? Thank you very much for attention.

Best regards,

Tingting Chen

[akohlmey]

Dear all,

I’m trapped in the calculation of region temperature. The version of
20160216 lammps is employed. Two methods are used in my simulation about
inorganic aqueous solution. One is by the command of compute temp/region
and the other is by the command of compute ke/atom. The specific
procedures are as follows:

Both methods give the consistent temperatures in hot or cold region. But
it is very confusing that the calculated temperature is wrong. For example,
if the system is kept at 300K by Nose-Hoover thermostat, the region
temperature will be about 212K, which is quite unreasonable.

In order to find out what should account for the wrong region temperature,
I simulated *pure SPC/E water maintained at 300K *by ose-Hoover
thermostat under NVT ensemble.

Two conditions are simulated. A part of the log file is showed :

Results of the first condition using special_bonds lj/coul 0 0 0.5

Results of the second condition using default special_bonds lj/coul 0 0 0

From results above, mainly two problems exist.

Firstly, from http://lammps.sandia.gov/doc/compute_ke.html, thermo_ke
corresponds to thermo_temp and is different from compute_ke. But in the
first condition, thermo_ke doesn’t agree with thermo_temp. For both
conditions, thermo_ke equals compute_ke completely.

Secondly, from http://lammps.sandia.gov/doc/compute_temp_region.html,
compute temp/region differs from other compute temp because it does not
subtract out degrees-of-freedom due to fixes that constrain motion, such as
fix shake. In the first condition, temp/region and ke give wrong overall
temperature and ebond as well as eangle equals 0.

​your conclusion is incorrect. only compute temp/region is not producing
the desired output, and that is for the reason that you already mentioned,
i.e. that it doesn't account for removed DOFs due to fix shake.​

In the second condition, temp/region and ke give right overall temperature
but ebond and eangle are not 0. So I doubted whether both ebond and
eangle should be 0 when fix shake is used to constrain bonds and angles.
The only difference of inputfile in two conditions above is the
special_bonds settings. In the first condition, I think the auto-generation
special neigbors in lammps is wrong. From auto-generation in logfile,
when using lj/coul 0 0 0.5, 1-4 neigbors are 1 for water but it is 0
actually.

So, now, I have no idea how to get the right region temperature under the
condition when fix shake is used, ebond and eangle are always 0 regardless
of any special_bonds settings，which happens in my aqueous solution
simultion.

Can anyone tell me how to solve the problem? Thank you very much for
attention.

​you are not providing complete input decks, so one can only speculate, but
your second output shows, that fix shake does not constrain any bonds or
angles. hence having a non-zero bond and angle energy is correct and for
the same reason compute temp/region​ is able to provide correct output,
since no DOFs are removed. i can only speculate, that there is an
inconsistency in your data file, e.g. that the bond and atom types are not
the ones you provide to fix shake.

BTW: your concerns about special_bonds are invalid. for technical reasons
the maximal number of excluded neighbors has to be at least 1.

​axel.​

[111190]

Dear Professor,

Thanks for your prompt reply and useful suggestions.

In the previous simulation, I created water box in lammps by combining lattice command, create box command and water molecule template. Considering that you think it may be due to the data file, I use the data file generated by write_data in lammps at this time. It is fortunate that the same results are obtained regardless of the special_bonds settings. For convenience, the data file, input file and log file have been attached.

But the problem still exists that the system temperature calculated by temp/region is 200K when the thermo_temp is 300K. I think that although compute temp/region doesn’t subtract out DOF due to fix shake, there should be only a subtle deviation between thermo_temp and the temperature calculated by temp/region. But in my simulation a deviation of 100K occurs so this indicates that the right region temperature can’t be achieved for the constrained system.

In addition, thermo_ke always equals compute_ke and both ke agree with 200K when the termo_temp is 300K. This also indicates that the right region temperature can’t be obtained by compute ke/atom.

So, how can I get the right region temperature? This is vital for temperature profile of the system.

Best regards,

Tingting Chen

temptest.in (1.24 KB)

water.data (191 KB)

log.lammps (4.61 KB)

[akohlmey]

Dear Professor,

Thanks for your prompt reply and useful suggestions.

In the previous simulation, I created water box in lammps by combining
lattice command, create box command and water molecule template.
Considering that you think it may be due to the data file, I use the data
file generated by write_data in lammps at this time. It is fortunate that
the same results are obtained regardless of the special_bonds settings. For
convenience, the data file, input file and log file have been attached.

But the problem still exists that the system temperature calculated by
temp/region is 200K when the thermo_temp is 300K. I think that although
compute temp/region doesn’t subtract out DOF due to fix shake, there
should be only a subtle deviation between thermo_temp and the temperature
calculated by temp/region. But in my

​that is nonsense. please look up in a text book how temperature is
computed. ​the ratio of 200K to 300K is *exactly* what you should get, if
you do not remove the constrained DOFs. you have 9 DOFs (3x3) per water
molecule and you have 3 constraints (2x bond 1x angle), that makes a total
of 6 DOFs.

simulation a deviation of 100K occurs so this indicates that the right
region temperature can’t be achieved for the constrained system.

​that conclusion is incorrect.
​

In addition, thermo_ke always equals compute_ke and both ke agree with
200K when the termo_temp is 300K. This also indicates that the right
region temperature can't be obtained by compute ke/atom.

​that conclusion is incorrect as well. the documentation of compute
temp/region *​explicitly* says that it does not handle constrainted DOFs
*at all* and that is what you see.

So, how can I get the right region temperature? This is vital for
temperature profile of the system.

​first off. the main issue here is due to your not understanding the
problem correctly and completely.​ LAMMPS does exactly what its
documentation says it will do and thus is working correctly. to have this
work automatically and correctly in regions for all kinds of systems
(including correct handling of extended particles) would require a
significant programming effort and add a lot more complexity to something
that is already quite complex. hence the decision to - until somebody finds
a better way - just document the situation and let people correct for it
after the fact, where possible.

thus for the case of a pure water system with a rigid water potential using
fix shake, you *do* know the amount of DOFs removed per atom, you can
correct your temperature profile after the fact. if you set your thermostat
to 200K and use compute temp/region for biasing, you should get the biasing
corresponding to 300K. i would also consider not removing the DOFs due to
center of mass translation (default is on) with compute_modify <compute ID>
extra 0.

in the specific case of trying to match compute/temp region with compute
temp for the entire system, you have to consider that compute temp will
remove 3 extra DOFs, but due to having to rescale the temperature for
constraints, the equivalent amount for compute temp/region would be 4.5.

thus to set up compute temp/region to output exactly the same property that
compute temp computes for the entire system, you need to change your input
as follows:

compute t all temp/region box #calculate overall temperature of
the system
compute_modify t extra 4.5
variable t equal 1.5*c_t

from the thermo output, you can see that the value of v_t exactly matches
the value for the entire box including the COM removal.

Step Temp c_t v_t KinEng c_ke E_bond E_angle
       0 300 200 300 893.34877 893.34877
         0 0
    1000 290.38226 193.58817 290.38226 864.70877 864.70877
         0 0

however, regional temperature is usually not to be considered translational
invariant, so the correct code for (smaller) regions would be

compute t all temp/region <region ID>
compute_modify t extra 0
variable t equal 1.5*c_t

axel.

p.s.: please note that using fractional numbers with compute_modify extra
will require an update to a newer version of LAMMPS.

[111190]

Dear Professor,

I am truly grateful for your explanations in detail which benefit me a lot.

Maybe now I have understood the problem to some degree. For example, I have a system containing 100 water molecules, 10 Li⁺ and 10 Br^-. In simulation, only water molecules are constrained by fix shake. Then I will get a ratio of (6100+310+310)/( 9100+310+310) between temperature from compute temp/region and thermo_temp. Is that right?

Moreover, if I use fix langevin to maintain a region to 290K and another region to310K after the pure water system was equilibrated at 300K, I can get the right temperature profile by times a ratio. I will also get both negative values of f_hot and f_cold directly as follows.

compute t1 all temp/region box1

compute t2 all temp/region box2

compute_modify t1 extra 0

compute_modify t2 extra 0

variable t1 equal c_t1*1.5

variable t2 equal c_t2*1.5

fix 2 all nve

fix hot all langevin 310 310 100 59804 tally yes

fix cold all langevin 290 290 100 287859 tally yes

fix_modify hot temp t1

fix_modify cold temp t2

thermo_style custom step temp v_t1 v_t2 f_hot f_cold

run 200000

PPPM initialization …

G vector (1/distance) = 0.277089

grid = 15 15 15

stencil order = 5

estimated absolute RMS force accuracy = 0.0155483

estimated relative force accuracy = 4.68231e-005

using double precision FFTs

3d grid and FFT values/proc = 10648 3375

Memory usage per processor = 8.7137 Mbytes

Step Temp t1 t2 hot cold

100000 308.31886 287.50315 311.05773 0 0

101000 306.644 310.69684 300.46558 -1875.86 -1725.275

102000 299.91158 303.73445 292.76817 -3736.2948 -3501.3728

103000 305.25255 323.19538 279.59118 -5553.8168 -5214.8888

104000 309.45164 328.9651 314.42733 -7438.607 -6909.2763

105000 308.4474 305.87449 305.87099 -9271.7098 -8656.4075

But in fact, f_hot should be negative and f_cold should be positive if fix shake isn’t applied. Can I get the actual accumulated energy corresponding to 300K by processing some thermodynamic parameters? This makes me confused.

Best regards,

Tingting Chen

langevin.in (1.83 KB)

water.data (191 KB)

[akohlmey]

Dear Professor,

I am truly grateful for your explanations in detail which benefit me a
lot.

Maybe now I have understood the problem to some degree. For example, I
have a system containing 100 water molecules, 10 Li+ and 10 Br-. In
simulation, only water molecules are constrained by fix shake. Then I will
get a ratio of (6*100+3*10+3*10)/( 9*100+3*10+3*10) between temperature
from compute temp/region and thermo_temp. Is that right?

Moreover, if I use fix langevin to maintain a region to 290K and another
region to310K after the pure water system was equilibrated at 300K, I can
get the right temperature profile by times a ratio. I will also get both
negative values of f_hot and f_cold directly as follows.

​there are two problems here:

1) your thermostats should be using a target temperature that is divided by
1.5
2) langevin is not a good choice in this case here, due to its use of
random velocity changes.​ on average about 1/3rd of those will be along the
constrained coordinates and thus ineffectual and thus explain the imbalance
in the net thermostat energy

the following modified input segment addresses both issues:
- it computes the required target temperature with equal style variables
(thot and tcold)
- it uses temp/csvr, a velocity rescaling thermostat, that is compatible
with shake.

now the tallied thermostat energies match in magnitude and have opposite
signs and the kinetic energy profile along the z axis is consistent with
the hot-cold zones and the energy exchange between them.

axel.

compute t1 all temp/region box1
compute t2 all temp/region box2

compute_modify t1 extra 0
compute_modify t2 extra 0

variable t1 equal c_t1*1.5
variable t2 equal c_t2*1.5
variable thot equal 310/1.5
variable tcold equal 290/1.5
compute ke all ke

# equilibration and thermalization

fix 1 all shake 0.0001 20 0 b 1 a 1
fix 2 all nvt temp 300 300 100.0
velocity all create 300 102486 mom yes rot yes dist gaussian

timestep 1
thermo 1000
thermo_style custom step temp c_t1 v_t1 c_t2 v_t2 ke c_ke
run 100000
unfix 2

# additional short equilibration with hot/cold region

fix 2 all nve
fix hot all temp/csvr \{thot\} {thot} 100 151465
fix cold all temp/csvr \{tcold\} {tcold} 100 786452
fix_modify hot temp t1
fix_modify cold temp t2
run 10000

# reset accumulated thermostat energy statistics and timestep for
production run
unfix hot
unfix cold
run 0
reset_timestep 0

fix hot all temp/csvr \{thot\} {thot} 100 151465
fix cold all temp/csvr \{tcold\} {tcold} 100 786452
fix_modify hot temp t1
fix_modify cold temp t2

# output kinetic energy profile
compute ke_pa all ke/atom
compute bin all chunk/atom bin/1d z 0.0 1.0 units lattice
nchunk once
fix dist all ave/chunk 10 1000 10000 bin c_ke_pa &
                                       file 1d-spatial-distribution.dat

thermo_style custom step temp c_t1 v_t1 c_t2 v_t2 f_hot f_cold
run 500000

[111190]

Thank you very much. I have used this method to simulate the conductivity. The results agree well with that from literatures.