Hi all,
I’m just starting to run MD simulations with lammps and my thermostats (Berendsen and NH) aren’t thermostatting. I’m working with a classic Si tersoff system (216 atoms). I’ll give it an initial temp with the velocity command and use a fix nvt command set to 300. When I let the system run the temp will drop to half of what I initially give it with the velocity command and ignore the fix nvt command. I think that the temperature always drops to half of what I initially set it because half the KE is transitioning to PE.
So… bottom line, I’m looking for clues as to why my thermostats aren’t working? My input and log file are below.
Thanks,
Zack
Input
#Exploration of different temp fixes
#initialization
boundary p p p
units metal
echo screen
newton on
log log.si_temp
#definitions
lattice diamond 5.43
region all prism 0.0 3 0.0 3 0.0 3 0.0 0.0 0.0
create_box 1 all
create_atoms 1 region all
pair_style tersoff
pair_coeff * * Si.tersoff Si
neighbor 2.0 bin
neigh_modify every 1
mass 1 28.0
dump 1 all xyz 100 position.xyz
thermo 100
#thermostatting
velocity all create 1000.0 429349 dist gaussian
fix 1 all nvt 300 300 10
timestep .001
run 10000
Output
Step Temp E_pair E_mol TotEng Press
0 1000 -999.9909 0 -972.19997 7952.3189
100 613.55782 -989.23983 0 -972.18849 8398.2308
200 532.03165 -986.97415 0 -972.1885 8989.3852
300 573.17419 -988.12181 0 -972.19276 8817.7695
400 430.86479 -984.16833 0 -972.19419 9158.1246
500 478.74264 -985.50663 0 -972.20193 11194.764
600 491.4608 -985.8719 0 -972.21374 9397.6173
700 498.28864 -986.07553 0 -972.22762 8836.9486
800 495.56364 -986.01114 0 -972.23896 9942.3414
900 496.32827 -986.04765 0 -972.25423 10558.343
1000 482.62876 -985.68497 0 -972.27227 9364.1523
1100 518.41167 -986.69762 0 -972.29048 9883.3242
1200 482.57125 -985.72277 0 -972.31167 9772.3385
1300 582.36152 -988.52095 0 -972.33658 9087.0875
1400 550.32941 -987.65295 0 -972.35878 10031.153
1500 447.90941 -984.83049 0 -972.38267 9708.5834
1600 541.28672 -987.45493 0 -972.41207 9617.1049
1700 482.00202 -985.83551 0 -972.44022 9030.7912
1800 454.06124 -985.08882 0 -972.47004 10499.057
1900 433.06498 -984.53659 0 -972.50131 9687.6716
2000 468.1601 -985.54799 0 -972.53738 8314.7335
2100 501.24598 -986.50325 0 -972.57316 10302.656
2200 468.85276 -985.63712 0 -972.60726 10278.153
2300 535.6572 -987.53453 0 -972.64812 10713.986
2400 479.53007 -986.01454 0 -972.68795 10002.678
2500 372.97779 -983.08997 0 -972.72457 10019.359
2600 477.04398 -986.02905 0 -972.77155 9618.5371
2700 507.30769 -986.91157 0 -972.81301 10203.965
2800 484.81736 -986.32957 0 -972.85604 10240.417
Hi all,
I’m just starting to run MD simulations with lammps (newest version) and my thermostats (Berendsen and NH) aren’t thermostatting. I’m working with a classic Si tersoff system (216 atoms). I’ll give it an initial temp with the velocity command and use a fix nvt command set to 300. When I let the system run the temp will drop to half of what I initially give it with the velocity command and ignore the fix nvt command. I think that the temperature always drops to half of what I initially set it because half the KE is transitioning to PE.
So… bottom line, I’m looking for clues as to why my thermostats aren’t working? My input and log file are below.
Thanks,
Zack
Input
#Exploration of different temp fixes
#initialization
boundary p p p
units metal
echo screen
newton on
log log.si_temp
#definitions
lattice diamond 5.43
region all prism 0.0 3 0.0 3 0.0 3 0.0 0.0 0.0
create_box 1 all
create_atoms 1 region all
pair_style tersoff
pair_coeff * * Si.tersoff Si
neighbor 2.0 bin
neigh_modify every 1
mass 1 28.0
dump 1 all xyz 100 position.xyz
thermo 100
#thermostatting
velocity all create 1000.0 429349 dist gaussian
fix 1 all nvt 300 300 10
timestep .001
run 10000
Output
Step Temp E_pair E_mol TotEng Press
0 1000 -999.9909 0 -972.19997 7952.3189
100 613.55782 -989.23983 0 -972.18849 8398.2308
200 532.03165 -986.97415 0 -972.1885 8989.3852
300 573.17419 -988.12181 0 -972.19276 8817.7695
400 430.86479 -984.16833 0 -972.19419 9158.1246
500 478.74264 -985.50663 0 -972.20193 11194.764
600 491.4608 -985.8719 0 -972.21374 9397.6173
700 498.28864 -986.07553 0 -972.22762 8836.9486
800 495.56364 -986.01114 0 -972.23896 9942.3414
900 496.32827 -986.04765 0 -972.25423 10558.343
1000 482.62876 -985.68497 0 -972.27227 9364.1523
1100 518.41167 -986.69762 0 -972.29048 9883.3242
1200 482.57125 -985.72277 0 -972.31167 9772.3385
1300 582.36152 -988.52095 0 -972.33658 9087.0875
1400 550.32941 -987.65295 0 -972.35878 10031.153
1500 447.90941 -984.83049 0 -972.38267 9708.5834
1600 541.28672 -987.45493 0 -972.41207 9617.1049
1700 482.00202 -985.83551 0 -972.44022 9030.7912
1800 454.06124 -985.08882 0 -972.47004 10499.057
1900 433.06498 -984.53659 0 -972.50131 9687.6716
2000 468.1601 -985.54799 0 -972.53738 8314.7335
2100 501.24598 -986.50325 0 -972.57316 10302.656
2200 468.85276 -985.63712 0 -972.60726 10278.153
2300 535.6572 -987.53453 0 -972.64812 10713.986
2400 479.53007 -986.01454 0 -972.68795 10002.678
2500 372.97779 -983.08997 0 -972.72457 10019.359
2600 477.04398 -986.02905 0 -972.77155 9618.5371
2700 507.30769 -986.91157 0 -972.81301 10203.965
2800 484.81736 -986.32957 0 -972.85604 10240.417
Hi all,
I'm just starting to run MD simulations with lammps (newest version) and my
thermostats (Berendsen and NH) aren't thermostatting. I'm working with a
classic Si tersoff system (216 atoms). I'll give it an initial temp with the
velocity command and use a fix nvt command set to 300. When I let the
system run the temp will drop to half of what I initially give it with the
velocity command and ignore the fix nvt command. I think that the
temperature always drops to half of what I initially set it because half the
KE is transitioning to PE.
So... bottom line, I'm looking for clues as to why my thermostats aren't
working? My input and log file are below.
zack,
your system will have to equilibrate. that takes time.
have you checked with nve integration, if your system
conserves energy?
a nose-hoover thermostat does not exchange energy
that quickly. if you want to have a faster energy exchange,
you might want to try out the langevin thermostat or
use temp/rescale. you *do* start from 1000K and that
will stir up your system a bit. i assume that is intended,
but then you also have to give it time to settle.
cheers,
axel.
Your thermostat time constant is 10, which is 10,000
timesteps in metal units. So you probably haven't run
long enough to see equilibration. You can probably use
a time constant of 0.1 or 100 steps.
Steve
Thanks Axel and Steve for replying. I ran my simulation for 1,000,000 iterations (1 nanosecond altogether) and the thermostat worked. I’d still like to have a better understanding of which thermostats to use where and how to estimate the time to equilibration. Axel, the temp/rescale and langevin seem completely unphysical but in the end do they accomplish the same thing as the NH thermostat? I understand that even the NH thermostat is unphysical but at least it mirrors the canonical ensemble distribution. If I’m just prepping a system for NVE time integration, why would I want to use the NH thermostat if temp/rescale works faster?
Steve, where did I set the thermostat time constant? Is that the TDamp argument in the nvt fix? What does the time constant mean (is it the time to reduce temperature by some amount)? Like I mentioned earlier, I’m new to the computer simulation game so I apologize in advance for asking trivial questions.
Best,
Zack
I would say temp/rescale is unphysical, but Langevin is physical
in the sense it mimics coupling to a background solvent. The time
constant is a param in the fix nvt command - see the doc page
for what the "time" associated with it means. If you start
far from your target temp, then it can take a few time constants
to get there. Usually a time constant of about 100 timesteps is
a good value to use. With the NH chains option (on by default)
it seems to take longer to move to a target temp (from an initially
different temp), so you can try turning that off it well. However chains
does damp unwanted oscillations.
Steve