3.39 Calculating polarization of solids with FHI-aims

This section describes the relevant keywords connected to the implementation of the Berry-phase formalism within FHI-aims. The theory behind these flags will soon be summarized here.

The polarization direction pol_direction accepts the values 1, 2, or 3, corresponding to the index of the reciprocal space lattice vector along which the polarization is then calculated. The k-point grid defined by the last three fields (n_kpoints_dir1 n_kpoints_dir2 n_kpoints_dir3) should contain integer numbers and defines the k-point density along the reciprocal lattice vectors 1, 2, and 3, respectively. On this mesh, the required eigenstates and wave functions are obtained non-self consistently by Fourier-interpolating the electronic density obtained during the SCF cycle, in close analogy to band-structure or density-of-states calculations.

Note that the polarization is only defined modulo the polarization quantum, therefore, when evaluating polarization differences between different geometries, it is thus necessary to ensure that the polarization changes continuously between different geometries and does not jump by a full quantum. This can be checked by constructing discrete paths between geometries, e.g. between a non-centrosymmetric and centrosymmetric structure. For the latter, the polarization has to vanish. Therefore, relaxed geometries that fulfill all symmetry constraints are necessary to obtain accurate polarization differences.

Multiple polarization calculations can be requested within one FHI-aims run, e.g., to obtain the polarization along all reciprocal lattice vectors:

output polarization 1 20 5 5
output polarization 2 5 20 5
output polarization 3 5 5 20

For each one of these directives, the polarization is evaluated independently. A summary for all directives can be found at the end of the polarization calculations, e.g.,

Summarizing all directives:
Directive 1 in direction of rec. latt. vec. 1 yields the full polarization: 472.792E-03 (C/m²)
Directive 1 in direction of rec. latt. vec. 2 yields the full polarization: 0.0E00 (C/m²)
Directive 1 in direction of rec. latt. vec. 3 yields the full polarization: 0.0E00 (C/m²)

Whenever the polarization is calculated along all three directions within the same run, the code will also output the polarization along the Cartesian axes x, y, z , e.g.,
Cartesian Polarization 403.809855E-06 -684.500436E-06 10.214219E-03

This is particularly useful for non-orthogonal lattice vectors, for which (reciprocal) lattice vectors and Cartesian axes do not coincide.

A quick note on convergence: Generally, a much larger density of k-points is required for polarization calculations than for the SCF cycle. In particular along the reciprocal lattice vector pol$\mathrm{\_}$direction along which the polarization is evaluated, as the example above shows. The reason is that the Berry-connection is evaluated using finite-differences for the wave-function derivatives with respect to this k-path.

The index num_of_planes denotes the following planes in reciprocal space:

1. 1.

   k_x=[-0.5, 0.5] k_y=[0.0, 0.5] k_z= 0.0 || k_x = kparallel / k_y = kperpendicular

2. 2.

   k_x= [-0.5, 0.5] k_y=[0.0, 0.5] k_z= 0.5 || k_x = kparallel / k_y = kperpendicular

3. 3.

   k_x= 0.0 k_y=[-0.5, 0.5] k_z= [0.0, 0.5] || k_y = kparallel / k_z = kperpendicular

4. 4.

   k_x= 0.5 k_y=[-0.5, 0.5] k_z= [0.0, 0.5] || k_y = kparallel / k_z = kperpendicular

5. 5.

   k_x= [0.0, 0.5] k_y=0.0 k_z= [-0.5, 0.5] || k_z = kparallel / k_x = kperpendicular

6. 6.

   k_x= [0.0, 0.5] k_y=0.5 k_z= [0.0, 0.5] || k_z = kparallel / k_x = kperpendicular

The code always iterates over all planes between 1 and num_of_planes and for each plane n the evolution of the WCC can be found in the output file WCCn.dat. For each plane n, the output file WCCn.dat allows to investigate how many times the Berry phase is cycled through the Brillouin zone and thus allows to determine the topological invariant Z2 for this particular plane [162]. In practice, this can be done by plotting the WCC evolution and counting how often the WCC evolution crosses an arbitrary line parallel to the abscissa [324]. An odd number of crossings implies Z2=1 (topological), an even number Z2=0 (trivial insulator). See Fig.3.10 for some examples.

Alternatively, the value of Z2 can be determined automatically using the approach proposed by [283] and implemented in the script get_z2.py that can be found in the utilities folder of the FHI-aims distribution. The main idea behind the algorithm is that topological indexes (in fact winding numbers) can be determined by tracking the biggest gap in phase between the evolution of the individual WCC centers along the WCC path, see [283, 114]. In practice, it is sufficient to copy get_z2.py into the folder containing the WCCn.dat files and to execute it there. Besides plotting the WCC evolution, the determined Z2 indexes for the individual planes n are reported both in the title of the graphes and in the file output_z2.dat.

To distinguish between weak and strong topological insulators [99] Z2 needs to be investigated in multiple planes. For two-dimensional materials with a thick vacuum layer along the Cartesian z axes, it is sufficient to study the WCC evolution in the k_z = 0 plane (num_of_planes=1), since there is no dispersion in z direction. For three-dimensional crystals with three-fold rotational symmetry, it is enough to investigate two planes, i.e., one intersecting the Gamma point and one at the border of the Brillouin zone, by using num_of_planes=2. The topological character of the material is then given by the sum of the two obtained indexes modulo 2.

In the general, three-dimensional case, materials are characterized by a set of 4 indexes v₀;v₁,v₂,v₃ and all six planes (num_of_planes=6 ) in the BZ need to be investigated [99]. The first, strong index v₀ distinguishes between strong and weak topological insulators. The last three indices v₁,v₂,v₃ are called weak topological invariants [100] and denote the value of Z2 for the planes at the BZ boundaries, i.e., those planes in which the index is associated with num_of_planes is 2,4, and 6. For instance, Bi₂Se₃ with $1;0,0,0$ implies that Bi₂Se₃ is a strong topological insulator (there are topologically protected surface states).