A Python code for calculating the mean-value (Baldereschi’s) point for any crystal structure

In his original paper Baldereschi defined the mean-value point in the Brillouin zone as: “the point such that the value which any given periodic function of wave vector assumes at this point is an excellent approximation to the average value of the same function throughout the Brillouin zone” [1]. The conditions that define the mean-value point (MVP) follow from the expansion of any periodic function defined on the $\mathbf{k}$-space into plane waves:

where $f$ stands for a periodic function of the wave vector $\mathbf{k}$, the $\mathbf{R}$ represent vectors of the direct (Bravais) lattice, which take the role of the wave vectors in the above expansion, and $F(\mathbf{R})$ are the Fourier components of the function $f$. Generally, the function $f$ needs to fulfill Dirichlet conditions for the expansion from eq. (1) to converge to $f$.

From eq. (1) it is easy to show that the integral of $f$ over the entire Brilloiun zone equals the $F(\mathbf{R}=0)$ term in the expansion multiplied by the volume of the Brilloiun zone ($(2\pi)^{3}/V$). Hence, if there were a point $\mathbf{k}=\mathbf{k}_{\mathrm{o}}$ at which the sum of all terms in the expansion above, excluding the $\mathbf{R}=0$ term, vanishes, the value $f(\mathbf{k}_{\mathrm{o}})\,(2\pi)^{3}/V$ would be exactly equal to the integral of $f(\mathbf{k})$ over the entire Brillouin zone. The existence of $\mathbf{k}_{\mathrm{o}}$ is guaranteed by the mean-value theorem of the integral calculus, but its coordinates are dependent on the function $f$. The mean-value point from Baldereschi’s work is a point defined by the symmetry of the crystal that is independent of $f$, at which the value of a function $f$ is an “excellent” approximation of its integral over the Brillouin zone.

Hence, approximating the integral of $f(\mathbf{k})$ over the entire Brillouin zone can be reduced to finding $f(\mathbf{k}_{\mathrm{o}})$, provided that $\mathbf{k}_{\mathrm{o}}$ is known. In the original work of Baldereschi [1], the coordinates of $\mathbf{k}_{\mathrm{o}}$ are evaluated for the three cubic Bravais lattices, face centered cubic (fcc), body centered cubic (bcc), and simple cubic (sc). Subsequently, the mean-value point has been found for some non cubic lattices (see for example [2]). However, to the best of my knowledge, the mean-value point has not been evaluated for all 14 Bravais lattices, and a general, stndalone code for doing so for any input crystal structure is not available. Herein, both the mvp.py code and the resulting MVPs for the 14 Bravais lattice are presented, and various aspects of the procedure are discussed.

The key to evaluating the MVP for a given crystal structure is the requirement that the sum of all terms in eq. (1) for $\mathbf{R}\neq 0$ vanishes at some $\mathbf{k}=\mathbf{k}_{\mathrm{o}}$. These conditions provide a system of equations whose solution(s) define $\mathbf{k}_{\mathrm{o}}$. For the purpose of finding an approximate, function independent $\mathbf{k}_{\mathrm{o}}$ the sum from eq. (1) can be rewritten in a more suitable form using the symmetry of the crystal. As already discussed by Baldereschi in Ref. [1] one can assume, without any loss of generality, that all $F(\mathbf{R})$ are equal for $\mathbf{R}$ vectors that are connected by the point group symmetry operations, that is, that belong to the same star of $\mathbf{R}$. One can then write:

where $s$ goes over all stars of $\mathbf{R}$, $n_{s}$ counts all the lattice vectors within a given star, and $W_{s}(\mathbf{k})$ are the “symmetrized” plane waves belonging to the $\Gamma_{1}$ (fully symmetric) irreducible representation of the point group.

Hence, the $f$-independent conditions that determine $\mathbf{k}_{\mathrm{o}}$ are that all $W_{s}(\mathbf{k}_{\mathrm{o}})=0$ for any $s>0$. This is, of course, too much to ask and this system of equations does not have a solution. However, two simplifications can be made. First, one can assume that $F_{s}$ decrease sufficiently fast with $s$, so that after a relatively small number of terms the non-zero value of $W_{s}$ becomes irrelevant. Second, even if the number of zero terms is smaller than three, which is necessary for all three components of $\mathbf{k}_{\mathrm{o}}$ to be uniquely determined, the values can be chosen so to minimize the first nonzero $W_{s}(\mathbf{k})$. These two assumptions allowed Baldereschi to evaluate the coordinates of the mean-value points for the three cubic systems (sc, bcc, fcc).

The mvp.py code follows the outlined protocol. For any given crystal structure, a set of lattice vectors $\mathbf{R}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2}+n_{3}\mathbf{a}_{3}$ is constructed. Here $\mathbf{a}_{1},\mathbf{a}_{2},\mathbf{a}_{3}$ are the unit cell vectors and the integers $n_{1},n_{2},n_{3}$ are chosen to cover a range of values so that a sufficiently large number of $\mathbf{R}$s is produced; large enough to enclose the first 4 stars of $\mathbf{R}$ ($s=0,1,2,3$). The number of $\mathbf{R}$ vectors is controlled by the R_range input integer variable ($-$ R_range $\leq n_{1},n_{2},n_{3}\leq$ R_range), whose default value is 4. All results presented in this paper are obtained with R_range = 4.

Then the point group symmetry operations are found with the help of spglib library [3]. For this purpose the tolerance factors for the spglib are set to symprec = 1e-02 Å (length) and angle_tolerance = $-1.0$ (angles). The full set of orthogonal transformations (rotations and mirror planes) including those that are parts of screw axes and glide planes is used for the classification of $\mathbf{R}$ vectors into stars. Conveniently, the spglib library provides a separation of the of the space group operations into orthogonal parts (labeled “rotations”) and associated translational parts (fractional translations for the glide planes and screw axes). The mvp.py code takes the entire set of orthogonal parts, and uses that set to classify $\mathbf{R}$ into stars.

In the next step, the $\mathbf{R}$ vectors are classified into classes of symmetry equivalent members forming the stars. This is done by applying the entire set of orthogonal transformations successively on $\mathbf{R}$ vectors and grouping them with those they transform into. The equality of a given $\hat{O}\mathbf{R}$ and some $\mathbf{R}^{\prime}$ is determined using the condition $|\hat{O}\mathbf{R}-\mathbf{R}^{\prime}|=0$ with the tolerance factor gen_tol = 1e-5. Also, applying the symmetry operations may map a given $\mathbf{R}$ to itself. All these duplicate occurrences are removed when constructing the stars. In the next step the functions $W_{i}(\mathbf{k})$ are constructed for $i=1,2,3,4$.

The roots of $W_{i}(\mathbf{k})$ functions are found in the following way. First, a $\mathbf{k}$-point grid is constructed spanning the one unit cell of the reciprocal lattice. The number of $\mathbf{k}$-points is controlled by the nkpts input parameter (default value 2e+7). When constructing the $\mathbf{k}$-point grid the unit vectors of the reciprocal lattice $(\mathbf{b}_{1},\mathbf{b}_{2},\mathbf{b}_{3})$ are divided into $(n_{1},n_{2},n_{3})$ sub-divisions so that the grid spacing is uniform. The grid spacing is computed as $\texttt{step}=[\mathbf{b}_{1}\cdot(\mathbf{b}_{2}\times\mathbf{b}_{3})/\texttt{nkpts}]^{1/3}$ and $n_{i}=|\mathbf{b}_{i}|/\texttt{step}$. Then, the values of $W_{1}(\mathbf{k})$ are evaluated over the entire grid and its zeros are used as starting points near which the simultaneous roots of the $W_{1}(\mathbf{k})$, $W_{2}(\mathbf{k})$ and $W_{3}(\mathbf{k})$ will be evaluated. The situation in which no zeros of $W_{1}(\mathbf{k})$ are found for any grid point is easily remedied by increasing the nkpts.

In the mvp.py code the $\mathbf{k}$-points at which $W_{1}(\mathbf{k})\leq\,$1e-12 (input parameter zero) are extracted and the function fsolve from scipy.optimize package is used to find roots of a vector function $\mathbf{func}_{1}=(W_{1}(\mathbf{k}),W_{2}(\mathbf{k}),W_{3}(\mathbf{k}))$ near every one of those $\mathbf{k}$-points. If the roots of $\mathbf{func}_{1}$ exist then the roots that minimize $W_{4}(\mathbf{k})$ are extracted and those that are possibly outside the reciprocal unit cell are mapped back inside. The list of these $\mathbf{k}$-points then represents the list of mean-value points (symmetry equivalent). User can then extract (on their own) those that belong to the first Brillouin zone, or choose to output only the $\mathbf{k}$-point from this set with the lowest norm (input parameter only_lowest_norm). In case $\mathbf{func}_{1}$ has no roots, the roots of $\mathbf{func}_{2}=(W_{1}(\mathbf{k}),W_{2}(\mathbf{k}))$ are then solved for and those that minimize $W_{3}(\mathbf{k})$ are extracted and mapped back into the reciprocal unit cell. Lastly, our tests indicate that relatively large default value of nkpts is unfortunately required to converge the MVP value. Using different solvers might help converge the MVP faster, but that will be explored in the future.

The mvp.py code is available on github [4] as part of the toolbox package [5] for pylada [6], a collection of python tools for structure prediction (random structure generation), surface cutting, various structure manipulations, etc. All dependencies include: python, numpy, scipy, spglib, pylada, and toolbox. The mvp.py is constructed around the pylada structure object but is easily adaptable to use other structure objects from other codes (e.g ASE [7], aflow [8], pymatgen [9],…), which will remove the pylada dependency.

Finally, the choice to write the mvp.py in Cartesian rather then more elegant crystal coordinates is mainly driven by the differences in how the integers are treated between python2 and python3. In this way the code works with both versions. The price to pay is the introduction of various tolerance factors whose default values are set to best reflect extensive testing. That said, the convergence of the results with respect to various tolerance factors needs to be tested before use.

Table LABEL:tab:results below lists the MVP coordinates for the 14 Bravais lattices. For clarity, the specific unit cells that are used are displayed together with the choice of unit vectors and their coordinates. The corresponding first Brillouin zones are also displayed as well as the location of the MVP points. In the first three rows the original results of Baldereschi are reproduced. The coordinates of the MVPs for the simple cubic, face centered and body centered cubic lattices correspond well to those reported in Ref. [1]. The rest are original results.

Importantly, coordinates of the MVPs in general depend on the choice of the lattice parameters. Table  LABEL:tab:results only lists MVP coordinates for the specific choices provided in the table. For other lattice parameters users should compute the MVPs on their own.

Here I would like to make several remarks, which may be relevant when using the code and/or understanding its structure.

It is first important to note that the systems of equations $\mathbf{func}_{1}=(W_{1}(\mathbf{k}),W_{2}(\mathbf{k}),W_{3}(\mathbf{k}))=(0,0,0)$ or $\mathbf{func}_{2}=(W_{1}(\mathbf{k}),W_{2}(\mathbf{k}))=(0,0)$ are non-linear in $cos(\alpha k_{i})$ with $\alpha$ being an integer multiple of $2\pi$. As there is no general approach of solving such systems of equations (to my knowledge at least) the choice was to use scipy.optimize.fsolve function that finds roots near a certain starting point. It turns out that a relatively large number of those starting points are needed to cover most, if not all, zeros of $W_{1}(\mathbf{k})$, which then allows for the roots of $\mathbf{func}_{1}$ or $\mathbf{func}_{2}$ to be found. Converged values of MVPs in Table LABEL:tab:results are obtained using nkpts between 1e+6 (simple cubic) and 3e+7 (hexagonal lattice). While this might look excessive, it is necessary for the scipy.optimize.fsolve to find correct solutions. The typical time needed for the mvp.py to solve the equations (on a personal computer) is measured in minutes.

Second, by default mvp.py code uses all orthogonal parts (reflections, rotations) belonging to the entire space group rather then the point group operations only. This is done so to account the fact that the entire set of orthogonal transformations may be larger than the point group alone, and that the a group of all pure translations is invariant under the operations of all orthogonal transformations including those belongin to screw axes and glide planes. Hence, when classifying $\mathbf{R}$ vectors into stars all orthogonal transformations need to be used. While this is only relevant when dealing with real crystals (not Bravais lattices), the mvp.py code is written in this way to provide access to all symmetry operations in Cartesian coordinates, which may be useful for other purposes.

Next, the MVPs are not invariant under the supercell transformations in spite of supercells being just different representations of the same crystal/lattice. That is easy to prove just by considering the three cubic lattices. Both fcc and bcc can be represented as simple cubic using their conventional unit cells. However, the $(1/4,1/4,1/4)$ $\mathbf{k}$-point in the crystal coordinates, which is the MPV for simple cubic lattice is not an MVP for neither fcc nor bcc. This is because the first star of $\mathbf{R}$ for both fcc and bcc (set of 12 and 8 lattice vectors, respectively) is not accounted for in the simple cubic lattice and consequently $W_{1}\neq 0$ for both of them at $\mathbf{k}=(1/4,1/4,1/4)$. The second star of $\mathbf{R}$ for fcc and bcc is the same as the first star for simple cubic, which then implies $W_{2}=0$ at $\mathbf{k}=(1/4,1/4,1/4)$ and so on.

Lastly, it is important to keep in mind that the MVP coordinates do in general depend on the choice of the lattice parameters. This is true for Cartesian as well as the crystal coordinates as already discussed in Ref. [2] for tetragonal and rhombohedral lattices. It is precisely this dependence that motivated writing the mvp.py code rather then tabulating all cases for all 14 Bravais lattices.

In conclusion, the mvp.py code is presented for computing the mean-value (Baldereschi’s) point in the Brillouin zone for any crystal structure. The underlying theory is also presented and discussed as are various aspects of the specific implementation in the mvp.py code. The code itself relies on the pylada python package [6] a high-throughput computational physics framework, but is easily adaptable to other similar packages.