6
$\begingroup$

I want to perform NVT simulation to study defect dynamics in sic. I want to study dynamics at 1500K. despite of giving initial and final temperature 1500K (also damping parameter as required part of command) in fix nvt command (input and log file is given below), I can see an exponential increase in temperature. I want to understand the implementation, why lammps is not starting the simulation directly form 1500K.

Input:


units           metal
dimension       3
atom_style      atomic
pair_style      tersoff
boundary        p p p

read_data 3csic555_cvac
#reset_atom_ids sort no

#set atomic attributes

velocity all create 1500.0 348769 rot yes dist gaussian loop geom
#giving atominc interactions
newton on on
pair_style      tersoff
pair_coeff      * * /usr/local/Cellar/lammps/20220623/share/lammps/potentials/SiC.tersoff Si C

#for time intergraion and other constraints

fix             1 all nvt temp 1500.0 1500.0 100.0
#fix             2 all nve
timestep 0.001
dump            mydmp all atom 500 dump.min





thermo 500

run             500000

Logfile:

Per MPI rank memory allocation (min/avg/max) = 3.834 | 3.834 | 3.834 Mbytes
Step Temp E_pair E_mol TotEng Press
       0          500   -6141.8538            0    -6077.353   -76639.045
   10000    276.43509   -6112.3713            0   -6076.7108    -68125.88
   20000    273.32228   -6110.8914            0   -6075.6324   -66875.866
   30000    285.43628   -6111.2779            0   -6074.4562   -66917.301
   40000    278.13412   -6109.1826            0   -6073.3029   -65746.928
   50000    289.22109   -6109.3216            0   -6072.0116   -65597.833
   60000    301.54955   -6109.4815            0   -6070.5811   -64343.259
   70000    306.81328   -6108.5685            0   -6068.9891   -66579.906
   80000    302.67689     -6106.38            0   -6067.3342   -64186.072
   90000    313.68726   -6106.1345            0   -6065.6683   -64413.954
  100000    326.02443   -6105.9705            0   -6063.9128   -63213.492
  110000    328.27444   -6104.3621            0   -6062.0141   -62684.361
  120000    340.34542   -6103.8526            0   -6059.9475   -60997.521
  130000    352.93407   -6103.3545            0   -6057.8254   -62076.109
  140000     357.1562   -6101.3489            0   -6055.2752   -59137.752
  150000    381.54196   -6102.2911            0   -6053.0716   -57748.092
  160000     387.1425   -6100.4918            0   -6050.5498   -58410.882
  170000    380.66263   -6097.0521            0    -6047.946   -54899.128
  180000     403.8425   -6096.9102            0   -6044.8139   -57878.013
  190000    412.23799   -6095.2316            0   -6042.0522   -51955.885
  200000    421.22362   -6093.2239            0   -6038.8854   -51065.847
  210000    433.73796   -6091.7092            0   -6035.7563   -49607.731
  220000    442.23856   -6089.7646            0   -6032.7151   -45590.764
  230000    469.97452   -6089.9588            0   -6029.3313   -45771.141
  240000    498.00504   -6089.9526            0   -6025.7092   -47410.634
  250000    472.52756    -6082.563            0   -6021.6062   -38739.338
  260000    512.67791   -6083.8665            0   -6017.7302   -44930.882
  270000    539.24371   -6083.2532            0   -6013.6898   -37262.625
  280000    547.50321   -6080.2904            0   -6009.6616   -36256.739
  290000    591.67329   -6081.3809            0   -6005.0541   -38453.659
  300000    589.01562    -6076.237            0    -6000.253   -33055.161
  310000    612.20655   -6074.1156            0     -5995.14   -25625.661
  320000     628.9901   -6071.3865            0   -5990.2457   -27308.343
  330000    663.26313   -6070.0644            0   -5984.5024   -21724.991
  340000
Improve this question
$\endgroup$
4
  • $\begingroup$ +1 Welcome to our new community and thank so you so much for contributing your question here! We hope to see much more of you in the future!!! $\endgroup$
    – Nike Dattani
    Oct 26, 2022 at 22:13
  • $\begingroup$ I have used your code on BCC iron of 1000 atom with EAM potential it starts with 1500 and drops to T/2 750K and then slowly heat up. In general if you start with T kelvin initially it will go to T/2. It seems from your output you have given initial velocity at 500K. If possible share your input structure too so I can rerun again $\endgroup$
    – Pranav kumar
    Oct 27, 2022 at 3:29
  • $\begingroup$ Yes it's true I given velocity at 500K, but why it is so ? purpose of giving initial velocities should be just to initialize the system. maintaining the temperature should be linked with NVT ensemble. I was also thinking by default, it uses Noose Hover method to maintain the temperature, does change of method can help ? anyways here is my input structure. $\endgroup$
    – Irslan Ullah Ashraf
    Oct 27, 2022 at 9:06
  • $\begingroup$ units metal dimension 3 atom_style atomic pair_style tersoff boundary p p p read_data 3c_sic_555_si_vac velocity all create 500.0 87287 rot yes dist gaussian loop geom newton on on pair_style tersoff pair_coeff * * SiC.tersoff Si C fix 1 all nvt temp 1500.0 1500.0 90.0 timestep 0.001 dump mydmp all atom 5000 wdout.vac.1500K.dump.min thermo 10000 run 500000 $\endgroup$
    – Irslan Ullah Ashraf
    Oct 27, 2022 at 9:08

2 Answers 2

Reset to default
10
$\begingroup$

The instantaneous temperature will almost always fall to half of its original value when thermal equilibrium is established when we attribute velocities to the atoms in a crystal. The equipartition energy theorem is mainly responsible for this. The Hamiltonian of solids at temperatures much below their melting point mimics a set of coupled harmonic oscillators, hence the average kinetic energy and average potential energy of a harmonic system are equally distributed throughout. In an NVE-ensemble simulation, the total energy is conserved, hence the kinetic energy must fall by half if the potential energy (which starts off at zero). In case of NVT, after each certain time step velocity of the system are rescaled to reach final temperature. In the initial run, when velocity corresponds to 500K was applied, overall system temperature drops to 250K and slowly it start gaining temperature as you run further. NVT can be thought as series of NVE simulation connected by event when velocity is being rescaled.

 fix 1 all nvt temp Tstart Tstop Tdamp

The LAMMPS manual says:

The desired temperature at each timestep is a ramped value during the run from Tstart to Tstop. The Tdamp parameter is specified in time units and determines how rapidly the temperature is relaxed. For example, a value of 10.0 means to relax the temperature in a timespan of (roughly) 10 time units

Improve this answer
$\endgroup$
9
$\begingroup$

As a consequence of these commands:

fix 1 all nvt temp 1500.0 1500.0 100
timestep 0.001

you are running your system with a Nose-Hoover thermostat whose time constant is 100,000 timesteps.

A Nose-Hoover thermostat will correct deviations from the preset temperature within a time scale of a few time constants. Since you begin with a system 1250K(1) away from the preset, it is not surprising that three time constants later your system has not reached its preset.

The LAMMPS documentation has quite a bit of information about what a thermostat actually does and how a proper time constant should be set -- make sure you read it!

(1) You initially set your system to 500K, but most likely the system was exactly at an energy minimum. So within the first few thousand timesteps the thermal kinetic energy was redistributed into potential energy, as the thermo output table indicates.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ thanks for your reply, I am not clear with your statement "it is not surprising that three time constants later your system has not reached its preset.", how you'r saying it "three time constants". Also time constant of 100.00 is the one recommended by Lammps for Nose Hover, as it's mentioned in documentation " A good choice for many models is a Tdamp of around 100 timesteps. ". considering also the start temprature is also 1500K, I think Noose Hoover should correct the temprature early. I am going to try with Noose Hoover chains, as my system(sic diamond) is most harmonic. $\endgroup$
    – Irslan Ullah Ashraf
    Oct 28, 2022 at 9:54
  • 3
    $\begingroup$ (1) the very next sentence in the documentation says: Note that this is NOT the same as 100 time units for most units settings. Pay attention to what units are used in the fix nvt and timestep commands. (2) the fix nvt default is a three-thermostat Nose-Hoover chain, as the Defaults section states (near the bottom of the page). (3) the initial internal thermal equilibrium is not your preset temperature but half that, for reasons stated in both answers on your question so far. $\endgroup$
    – Shern Ren Tee
    Oct 28, 2022 at 10:20
  • 1
    $\begingroup$ Thanks; now I am all clear, and the problem is resolved. your reply is helpful $\endgroup$
    – Irslan Ullah Ashraf
    Oct 31, 2022 at 13:43

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .