# Tutorial: Thermal expansion

Jump to navigation
Jump to search
### The

### The

## Contents

## Introduction

- A given crystal should have a well defined average lattice constant at a given pressure and temperature. Here we use silicon as an example to show how to calculate lattice constants using
`GPUMD`

. We use a cubic system (of diamond structure) consisting of [math]10^3\times 8 = 8000[/math] silicon atoms and use the minimal Tersoff potential [1].

- All the input and output files can be found here.

## Preparing the Inputs

### The `xyz.in`

file

- The first few lines of the
`xyz.in`

file are:

8000 4 3 0 0 0 1 1 1 54.3 54.3 54.3 0 0 0 0 28 0 0 2.715 2.715 28 0 2.715 0 2.715 28 0 2.715 2.715 0 28

- Explanations for the first line:
- The first number tells that the number of particles is 8000.
- The second number 4 in this line is good for silicon crystal described by the Tersoff potential, because no atom can have more than 4 neighbor atoms in the temperature range studied. One can make this number larger, which only results in using more memory. If this number is not large enough,
`GPUMD`

will give an error message and exit. - The next number 3 means that the initial cutoff distance for the neighbor list construction is 3 angstrom. The point here is that we only need to consider the first nearest neighbors. Therefore, any number larger than the first nearest neighbor distance and smaller than the second nearest neighbor distance is OK here. Note that we will also not update the neighbor list. There is no such need in this problem.
- The remaining three zeros in the first line mean:
- the box is orthogonal;
- the initial velocities are not contained in this file;
- there is no grouping method defined here.

- Explanations for the second line:
- The first three ones mean that all three directions are periodic.
- The remaining three numbers are the box lengths in the three directions. It can be seen that we have used an initial lattice constant of 5.43 angstrom to build the model.

- Starting from the third line, the numbers in the first column are all 0 here, which means that all the atoms are of type 0 (single atom-type system). The next three columns are the initial coordinates of the atoms. The last column gives the masses of the atoms. Here, we consider isotopically pure Si-28 crystal. In some applications, one can consider mass disorder in a flexible way.

### The `run.in`

file

- The
`run.in`

input file is given below:

potential potentials/tersoff/Si_Fan_2019.txt 0 velocity 100 ensemble npt_ber 100 100 0.01 0 0 0 0.0005 time_step 1 dump_thermo 10 run 20000 ensemble npt_ber 200 200 0.01 0 0 0 0.0005 dump_thermo 10 run 20000 ensemble npt_ber 300 300 0.01 0 0 0 0.0005 dump_thermo 10 run 20000 ensemble npt_ber 400 400 0.01 0 0 0 0.0005 dump_thermo 10 run 20000 ensemble npt_ber 500 500 0.01 0 0 0 0.0005 dump_thermo 10 run 20000 ensemble npt_ber 600 600 0.01 0 0 0 0.0005 dump_thermo 10 run 20000 ensemble npt_ber 700 700 0.01 0 0 0 0.0005 dump_thermo 10 run 20000 ensemble npt_ber 800 800 0.01 0 0 0 0.0005 dump_thermo 10 run 20000 ensemble npt_ber 900 900 0.01 0 0 0 0.0005 dump_thermo 10 run 20000 ensemble npt_ber 1000 1000 0.01 0 0 0 0.0005 dump_thermo 10 run 20000

- The first line of command tells that the potential to be used is specified in the file Si_Fan_2019.txt.

- The second line of the command tells that the velocities will be initialized with a temperature of 100 K.

- Then, the next 4 lines tell how to do the first run. This run will be in the NPT ensemble, using the Berendsen method. The temperature is 100 K and the pressures are zero in all the directions. The coupling constants are 0.01 (dimensionless) and 0.0005 (in the natural unit system adopted by GPUMD) for the thermostat and the barostat, respectively. The time step for integration is 1 fs. There are [math]2 \times 10^4[/math] steps for this run and the thermodynamic quantities will be output every 10 steps.

- After this run, there are 5 other runs with the same parameters but the target temperature. Note that the time step only needs to be set once if one wants to use the same time step in the whole simulation. In contrast, one has to use the dump_thermo keyword for each run in order to get outputs for each run. That is, we can say that the time_step keyword is propagating and the dump_thermo keyword is non-propagating.

## Results and Discussion

- It takes less than 1 min to run this example when a Tesla K40 card is used. The speed of the run is about [math]3 \times 10^7[/math] atom x step / second. Using a Tesla P100, the speed is close to [math]10^8[/math] atom x step / second.

- The output file
`thermo.out`

contains many useful data. The results are shown in the following figure:

- (a): The temperature for each run quickly reaches the target temperature (with fluctuations).

- (b): The pressure (averaged over the three directions) for each run quickly reaches the target pressure zero (with fluctuations).

- (c): The lattice constant (averaged over the three directions) for each run reaches a plateau (with fluctuations) after some steps.

- (d): We calculate the average lattice constant at each temperature by averaging the second half of the data for each run. The average lattice constants at different temperatures can be well fit by a linear function, with the thermal expansion coefficient being estimated to be [math] \alpha \approx 7.2 \times 10^{-6} [/math] K
^{-1}.

## Related pages

- Tutorial: Density of states
- Tutorial: Thermal conductivity from EMD
- Tutorial: Thermal transport from NEMD and HNEMD

## References

[1] Zheyong Fan, Yanzhou Wang, Xiaokun Gu, Ping Qian, Yanjing Su, and Tapio Ala-Nissila, A minimal Tersoff potential for diamond silicon with improved descriptions of elastic and phonon transport properties, accepted to JPCM. https://doi.org/10.1088/1361-648X/ab5c5f