Skip to content

Exercise 7: The Role of the Atomic Motion

Estimated total CPU time: 70 min

Warning

In the following exercises, computational settings including the reciprocal space grid (tag k_grid), the basis set, and supercell's size, have been chosen to allow for a rapid computation of the exercises in the limited time and within the CPU resources available during the tutorial session. Without loss of generality, these settings allow to demonstrate trends of the lattice dynamics of materials. In the production calculation, all computational parameters should be converged.

In this exercise, you will:

  • Learn how to create thermal displacements in supercells.
  • Investigate how the electronic band structure changes due to the atomic motion.

The lattice expansion is not the only effect that can alter the band gap: as a matter of fact, also the atomic motion leads to (instantaneous) changes in the electronic structure. In an experiment, which is usually performed on a time scale that is orders of magnitude larger than the period of the typical vibration in a solid, we thus only measure the thermodynamic average of the electronic structure. To investigate this aspect, we could perform Molecular Dynamics simulations at different temperatures, as you have learned it in the Tutorial IV. However, we will not be able to perform ab initio MD calculations due to the limited time and computational resources available during our tutorial session. Furthermore, we are not interested in the dynamical properties of the moving atoms, but merely in representative snapshots of the atomic configurations at each time step when the system is in thermal equilibrium.

In general, reasonable displacement which mimic thermal motion could be constructed by utilizing harmonic force constants 4. Superposing obtained normal modes with randomized phase factors allows one to produce structures representative of different temperatures without introducing additional parameters. Therefore, in this exercise to generate configurations without any MD simulation, we will instead use the force constants obtained in previous calculations and generate configuration according to the following scheme 4:

In the harmonic approximation, where the force on an atom is fully determined by the force constants and the displacement of the other atoms,

\[ \begin{aligned} \mathbf{F}_I = - \sum\limits_{J} \Phi_{IJ} \Delta\mathbf{R}_{J}~, \qquad (11) \end{aligned} \]

the equation of motion for the displacement \(\Delta \mathbf{R}_i\) of atom \(i\) in a supercell with \(N\) atoms in total are given by

\[ \begin{aligned} \begin{pmatrix} m_1 \, \Delta \ddot{\mathbf{R}}_1 \\ m_2 \, \Delta \ddot{\mathbf{R}}_2 \\ \vdots \\ m_N \, \Delta \ddot{\mathbf{R}}_N \end{pmatrix} = \begin{pmatrix} \Phi_{11} & \Phi_{12} & \cdots & \Phi_{1N} \\ \Phi_{21} & \Phi_{22} & \cdots & \Phi_{2N} \\ \vdots & \vdots & \ddots & \vdots \\ \Phi_{N1} & \Phi_{N2} & \cdots & \Phi_{NN} \end{pmatrix} \begin{pmatrix} \Delta {\mathbf{R}}_1 \\ \Delta {\mathbf{R}}_2 \\ \vdots \\ \Delta {\mathbf{R}}_N \end{pmatrix} ~. \qquad (12) \end{aligned} \]

By diagonalizing the hessian, one obtains \(3N\) eigenvalues \(\omega_s^2\) and eigenvectors \(\boldsymbol{e}_s\). This is equivalent to the solution of Eq. (4) if only the \(\Gamma\) point (\(\mathbf{q} = 0\)) is considered and the supercell is treated as one big unit cell.

The solution of the equations of motion given by Eq. (12) is then

\[ \Delta \mathbf{R}_i (t) = \frac{1}{\sqrt{m_i}} \sum_{s=1}^{N} \boldsymbol{e}_{is}\,A_{s}~\sin (\omega_s t + \varphi_s) \qquad (13) \]
\[ \Delta \dot{\mathbf{R}}_i (t) = \frac{1}{\sqrt{m_i}} \sum_{s=1}^{N} \boldsymbol{e}_{is}\,A_{s}\,\omega_s\,\cos (\omega_s t + \varphi_s) ,\qquad (14) \]

where \(A_{s}\) denotes the amplitude of phonon mode \(s\), and \(\phi_s\) is a random phase factor that can only be determined by boundary conditions. Since we are interested in representative samples, we could choose the displacements obtained from Eq. (13) at arbitrary times \(t\), so the factor \(\sin (\omega_s t + \varphi_s)\) will always just contribute a random number in the range \([-1, 1]\). 1 On the other hand, we can assign values to the amplitudes \(A_{s}\) by noting that, in thermal equilibrium, the equipartition theorem holds which states that each mode contributes an equal amount to the energy. As you know, the total energy of a set of harmonically coupled quantum oscillators is given by

\[ \begin{aligned} E = \sum_s \hbar \omega_s (n_\text{B} (\omega_s, T) + \frac{1}{2})~. \nonumber \end{aligned} \]

By this argument one finds that 5

\[ \begin{aligned} \left\langle A^2_{s} \right\rangle = \frac{\hbar}{\omega_{s}} \left( n_\text{B} (\omega, T) + \frac{1}{2} \right)~, \qquad (15) \end{aligned} \]

where \(n_\text{B}\) denotes the Bose-Einstein distribution function,

\[ \begin{aligned} n_\text{B} (\omega, T) = \frac{1}{\text{e}^{\frac{\hbar \omega}{k_\text{B} T}}-1}~, %~~\stackrel{T \rightarrow \infty}{\longrightarrow} %~~\frac{k_\t{B} T}{\hbar \omega} \end{aligned} \]

since phonons are bosons. With the mean amplitude given by Eq. (15), we can generate a Gaussian distribution of random amplitudes that give the correct average. By this means, we can generate thermodynamically correct displacements within the harmonic approximation from two random numbers. 2

!!! important Remark: Besides creating the samples for free, we can also mimic the effect of zero point motion at vanishing temperature: Since the amplitude was determined by a consideration rooted in quantum mechanics, it does not vanish when the temperature goes to zero:

\[ \begin{aligned} \langle A_s \rangle ~\stackrel{T \rightarrow 0}{\longrightarrow} ~\frac{\hbar}{2 \omega_s}~, \label{eq:A0K} \end{aligned} \]

whereas it would go to zero in the classical case. Verify by taking the \({T \rightarrow 0}\) limit in Eq. (15).

In order to investigate the influence of atomic motion on the band gap, we proceed in the following way: First we need to generate an ensemble of supercells with displaced atoms according to Eq. (13) for each temperature of interest, and compute the bandgap of each sample 3. Afterwards we can collect the band gaps and form the "thermodynamic average".

Create Samples

Our samples are created from force constants (harmonic) computed with Phonopy, which were already obtained before, during phonon spectrum calculations. Please, copy the results of the seconds exercise into current directory.

In order to produce files that contain harmonic force constants, we utilize interface between FHI-vibes and Phonopy. Go to the folder Tutorial/electron-phonon-coupling/7_bandgap_temperature/input. There you would find folder phonopy with the phonons calculatiosn of the Si, which we did before. Enter phonopy folder and type:

vibes output phonopy trajectory.son --full

Then, the file output/FORCE_CONSTANT is generated. Note that the file FORCE_CONSTANT uses the Phonopy file format for force constants, which saves the force constants in a reduced form. We need to remap them to 3\(N\)\(\times\)3\(N\) matrix, where \(N\) is the number of atoms in the supercell. This can be achieved by a FHI-vibes utility executed within the folder output:

vibes utils fc remap FORCE_CONSTANTS

which produces a file named FORCE_CONSTANTS_remapped. This can now be used to set up samples according to the trick explained above. To achieve this, we use the FHI-vibes tool create-samples. The basic command to create a supercell with displacements from force constants is

vibes utils create-samples geometry.in.supercell -fc FORCE_CONSTANTS_remapped -T 100 -n 10

This will create ten samples corresponding to the thermal displacements at 100 K.

Furthermore, one can request to use the Bose-Einstein distribution function to estimate the amplitudes by adding the flag --quantum:

vibes utils create-samples geometry.in.supercell -fc FORCE_CONSTANTS_remapped -T 100 -n 10 --quantum

Compute Band structures

We now want to investigate the impact of nuclear motion on the electronic band gap of silicon at different temperatures, for the classical case. We should create at least 10 samples (number of samples should be converged in production calculations) for each temperature between 0 K and 800 K in steps of 100 K, using the classical distribution. Store the samples in dedicated folders samples_000K, samples_100K ... samples_900K.

vibes utils create-samples geometry.in.supercell -fc FORCE_CONSTANTS_remapped -T 100 -n 10
mkdir samples_100K
mv geometry.in.supercell.0100K.* samples_100K
...

Once all the samples are generated for each temperature from the range of interest, we could start our calculations. There is a script run.py, which could help to do the work. The script is self-explanatory, please inspect it and run preferable in a separate folder, but note that folders with the samples should be presented there. For example, in the directory of the exercise do:

mv phonopy/output/samples* .
python run.py

Depending on your computer, the calculation of one set (10 samples) for one temperature might already take several minutes (you could also utilize data from solutions folder). After the calculations are finished, use the script postprocess.py to read in the band gaps and compute mean plus standard deviation for each temperatures. 

python postprocess.py

It will read in the band gaps, do the math and make a plot.

The script is again self-explanatory enough and helps you with extracting the band gaps and plot temperature dependence of the band gap. Which trends do you observe? You could generate samples with the Bose-Einstein distribution function and check what would be the differences between the classical and quantum case?

The plot shows the temperature dependence of the band gap in the classical limit. In the plot, the dots are the averaged band gaps among the ensemble for each temperature and the bars represent the corresponding standard deviation for each temperature.

As a final question: How do you think will the band gaps and error bars evolve when you go towards convergence, i. e., you increase the supercell size and the number of samples?

  1. Note that the displacements obtained from Eq. 13 are conceptually different from the single displaced atoms we encountered in the finite differences approach in exercises 1-3. ↩

  2. We use pseudo-random numbers. In real world applications of statistical sampling techniques, so called quasi-random sequences typically ensure better convergence behavior. Please refer to https://en.wikipedia.org/wiki/Low-discrepancy sequence for a general introduction to the topic and to the appendix of 6 for a concise discussion in a closely related context. As a matter of fact, you can use quasi random numbers with vibes by adding the flag --sobol. But with this few samples the differences should be negligible. ↩

  3. Remark: Since we destroy the local symmetry of the crystal by displacing the atoms, the minimal HOMO-LUMO gap will in general not be located on one of the paths connecting high symmetry points of the Brillouin zone, as we have calculated it in exercise 4. In this task, we will instead estimate the bandgap by simply calculating the HOMO-LUMO gap on a denser mesh of k-points. ↩

  4. D. West, and S. K. Estreicher, Phys. Rev. Lett. 96, 115504 (2006). ↩↩

  5. M. Dove, Introduction to Lattice Dynamics, Cambridge University Press (1993). ↩

  6. S. E. Brown, I. Georgescu, and V. A. Mandelshtam, J. Chem. Phys. 138, 044317 (2013) ↩