# Tutorial: Density of states

## Introduction

In this example, we calculate the phonon (vibrational) density of states (DOS) of graphene at 300 K and zero pressure. The method is based on the velocity auto-correlation (VAC) function. The DOS is calculated as the Fourier transform of the VAC .

## Preparing the Inputs

### The xyz.in file

• The first few lines of the xyz.in file created by the create_xyz.m MATLAB script are:
8640 3 2.1 0 0 0
1 1 0 149.649 155.52 3.35
0 1.24708 0 0 12
0 0 0.72 0 12
0 0 2.16 0 12
0 1.24708 2.88 0 12

• Explanations for the first line:
• The first number tells that the number of particles is 8640.
• The second number 3 in this line is good for graphene described by the Tersoff potential, because no atom can have more than 3 neighbor atoms at room temperature. 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 2.1 means that the initial cutoff distance for the neighbor list construction is 2.1 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 are no grouping methods defined here.
• Explanations for the second line:
• The numbers 1 1 0 mean that the x and y (in-plane) directions are periodic and the z direction is open (free).
• The remaining three numbers are the box lengths in the three directions. The box length in a free direction is chosen based on some convention. This number will only affect the system volume.
• 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 C-12 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/Graphene_Lindsay_2010_modified.txt
velocity    300
# equilibration
ensemble    npt_ber 300 300 0.01 0 0 0 0.0005
time_step   1
dump_thermo 1000
run         1000000
# production 1
ensemble    nve
compute_dos 5 200 400
run         100000
# production 2
ensemble    nve
compute_dos 5 200 400
run         100000
# production 3
ensemble    nve
compute_dos 5 200 400
run         100000
# production 4
ensemble    nve
compute_dos 5 200 400
run         100000
# production 5
ensemble    nve
compute_dos 5 200 400
run         100000

• The second line of the command tells that the velocities will be initialized with a temperature of 300 K.
• There are 6 runs.
• The first run serves as the equilibration stage, where the NPT ensemble (the Berendsen method) is used. The temperature is 300 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 $10^6$ steps (1 ns) for this run and the thermodynamic quantities will be output every 1000 steps.
• The other 5 runs are identical production runs. In each production run, the NVE ensemble is used. The line with compute_dos means that velocities will be recorded every 5 steps (5 fs) and 200 VAC data (the maximum correlation time is then about 1 ps) will be calculated. The last parameter in this line is the maximum angular frequency considered, $\omega_{\rm max} = 2\pi\nu_{\rm max} =400$ THz, which is large enough for graphene. Each production run lasts 100 ps. The major reason for using multiple production runs rather a single one is that computing the VAC requires a lot of memory, which prevents using very long production runs.

## Results and Discussion

• This simulation takes about 6 min when a Tesla K40 is used. The speed of this simulation, being about $4\times 10^7$ atom x step / second, is higher than that of the previous example because the number of neighbors for each atom is smaller here (numbers of atoms are comparable and the potential models are the same).
• Figure to be fixed.
• The figure below shows the calculated VAC and PDOS.
• For 3D isotropic systems, the results along different directions are equivalent and can be averaged, but for 2D materials like graphene, it is natural to consider the in-plane part (the $x$ and$y$ directions in the simulation) and the out-of-plane part (the $z$ direction) separately. It can be seen that the two components behave very differently. We can see that the cutoff frequency for the out-of-plane component (about 40 THz) is smaller than that for the in-plane component (about 52 THz).