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,
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
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
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
By this argument one finds that ^{5}
where \(n_\text{B}\) denotes the BoseEinstein distribution function,
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:
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 FHIvibes and Phonopy. Go to the folder
Tutorial/electronphononcoupling/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 FHIvibes 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 FHIvibes tool createsamples
. The basic command to create a supercell
with displacements from force constants is
vibes utils createsamples 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 BoseEinstein distribution
function to estimate the amplitudes by adding the flag
quantum
:
vibes utils createsamples 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 createsamples 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
selfexplanatory, 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 selfexplanatory 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 BoseEinstein 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?

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 13. ↩

We use pseudorandom numbers. In real world applications of statistical sampling techniques, so called quasirandom sequences typically ensure better convergence behavior. Please refer to https://en.wikipedia.org/wiki/Lowdiscrepancy 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. ↩ 
Remark: Since we destroy the local symmetry of the crystal by displacing the atoms, the minimal HOMOLUMO 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 HOMOLUMO gap on a denser mesh of kpoints. ↩

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

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

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