A Mathematical and Computational Proposal for the Development of Tight-Binding Order-N Density Matrix Methodology

ABSTRACT

The aim of this paper is to present, in a broad context, a proposal of mathematical foundation that serves as framework for the development of the electronic structure calculation methodology called “Density Matrix Tight Binding method (DMTB)”, which was originally presented in the literature by the group led by physicist David Vanderbilt in 1993⁹, as well as its relationship to the computational modeling that can be chosen for its implementation. The approach adopted makes it clear that the final mathematical formulation of this methodology is completely dependent on the computational strategies chosen for its effective implementation. Thus, we put the DMTB as a mathematical-computational model of variable final formulation. Finally, we propose an implementation based on the nonlinear conjugate gradient method. The final model obtained is slightly different from the DMTB that was originally presented in 1993, in agreement with the version presented by Millam and Scuseria in 1997¹¹. The approach used develops the mathematical aspects, aiming at the effective computational implementation of the methodology.

Keywords:
tight-binding density matrix; mathematical and computational modeling; Lagrange multipliers

1 INTRODUCTION

In general, poly-electronic systems have a rather intricate mathematical structure and quantum mechanical equations can only be solved using simplifying hypotheses such as the Born-Oppenheimer approximation and the independent electron approximations. Currently the theoretical investigations of electronic properties are made through computational simulations of mathematical models based on the basic methodologies of the quantum condensed matter theory: the semi-empirical and first principles methodologies. Many proposals for solving mathematical models that have quantum description in their formulations, run into the high computational complexity of the algorithms involved. This limits the application of these models to systems with many atoms, preventing more realistic simulations. Several mathematical-computational strategies with order complexity O(N) have been proposed, where N is the number of atoms involved. One is the semi-empirical methodology called the Tight-Binding Density Matrix Method (DMTB). Density matrix based methods for obtaining the ground state energy were proposed simultaneously in 1993 by Li, Nunes and Vanderbilt⁹ and Daw³ (based on different arguments and techniques), where the crystalline electronic states can be explicitly described in terms of atomic orbitals and the crystalline Hamiltonian can be constructed from semi-empirical parameters and thus only implicitly dependent on atomic orbitals. This methodology presents the possibility of computational implementation of linear complexity in relation to the number of electrons of the modeled system. DMTB has been widely used successfully for the study of electronic structure in the last two decades with methodological formulations different from those presented in the original article. An example can be seen in¹¹, for example.

This paper presents, in a broad context, a proposal of mathematical foundation for the DMTB. The approach used develops the mathematical aspects aiming the effective computational implementation of the methodology. In this construction, the final mathematical formulation of this method is completely dependent on the computational strategies chosen for its implementation. Thus, here, the DMTB is established as a mathematical-computational model of final variable formulation.

2 METHODOLOGY

Let’s consider the representation of crystal structures by supercells. We can then formulate the problem we are interested in as follows: Given a structure generated by repeating a supercell with N atoms and M electronic orbitals per atom, we obtain the quantum mechanical energy of the ground state of this structure.

The quantum mechanics described by the formalism of quantum mechanical states is based directly on the Schrödinger equation. However, there is a more general formulation, using a technique called the Density Operator (or Density Matrix). This is a deep topic that we cannot address here. For an initial reading, we recommend chapter 9 of¹. We will simply say that, in that approach, the concept of quantum state is generalized to elements of a new vector space over C, where states can be seen as particular cases, and the ρ density matrix is an projection operator onto subspace of the Occupied States of the H.

In this formalism, the Electronic Density Operator is given by

$\rho =\sum _n{f}_n\left|{\psi }_n><{\psi }_n\right|$

where n is the Hamiltonian discrete spectrum index, f _n is the occupancy number and the symbolism introduced by Dirac is used, in which the symbol | >, called ket, represents a vector of the vector space of quantum states, and the symbol < | , called bra, represents a vector of the dual topological of this space. For this symbolism, we recommend⁴, on page 18, and⁷, on page 32.

Writing the quantum states as a linear combination of the finite base of the Tight-Binding method gives the operator matrix ρ, called the electron density matrix ρ_e . Assuming that the orbital base forms an orthonormal set, we can write the number of electrons - N _e - and the total electronic energy of the system E as a function of the density matrix ρ_e , using the matrix trace operator

$N_e=tr\left({\rho }_e\right)=\sum _i{{\rho }_e}_{ii};E=tr\left({\rho }_{e}H\right)=\sum _{i,j}{{\rho }_e}_{ij}{H}_{ji}.$

Also, since the density operator is a projector, the density matrix is idempotent

${\rho }_e^2={\rho }_e.$

Recalling that in this type of methodology (Tight-Binding + Density Matrix) the total energy of a system with N _e electrons is usually given by:

$E_{tot}=E+E_{rep}+E_{0}N,$ (2.1)

where E _rep represents the repulsive potential, N is the number of system atoms and E ₀ is an energy constant per atom. For a presentation of the Tight-Binding method using the Electronic Density Operator formalism, as well as the demonstrations of the above identities, we recommend¹². In the approach of this reference, as in ours, the results were obtained assuming that the orbital base forms an orthonormal set.

3 PROBLEM FORMULATION

Since the quantum ground state is the minimum energy state, and since the ground state energy is a function of electron density (E = tr(ρ_e H)), we can now reshape the problem we are attacking: given a structure generated by repeating a supercell with N atoms and M electronic orbitals per atom with a Hamiltonian matrix H (order square NM), get the ρ_e electron density matrix that minimizes E(ρ) = tr(ρH). Note that we are considering energy as a real matrix variable function. From now on, we will be identifying the square matrix space M _n (ℝ) (n = NM ²) with the space ${\mathrm{ℝ}}^{n^2}$ and thus the energy E(ρ) will be considered a scalar field in ${\mathrm{ℝ}}^{n^2}$.

3.1 The Unrestricted Minimization Model

The minimization indicated in the formulation above cannot be done unrestrictedly. As we said, the quantum mechanics used in constructing the problem imposes that the number of electrons present in our structure must be reproduced by the ground state density matrix. In addition, the ground state density matrix must be idempotent. Finally, the electron density matrix must commute with the tight-binding Hamiltonian matrix. Given the computational strategies of interest, our goal is to obtain an unrestricted minimization problem. The electron number bond was introduced in the model through the Lagrange Multiplier Theorem:

Theorem 3.1. (Lagrange’s Theorem) Suppose that f : U → ℝ belongs to Class function C ^k (k ≥ 1) in the open set U ⊂ ℝ ⁿ⁺¹ , and M = g ⁻¹ (0) a hyper face contained in U, inverse image of the regular value 0 by a function g : U → ℝ, also class C ^k . Then p ∈ M is a critical point of f | _M if and only if there is a real number µ (called a Lagrange multiplier) such that

$\nabla f\left(p\right)+\mu \nabla g\left(p\right)=0.$

Considering the identification of M _n (ℝ) (n = NM ²) with the space ${\mathrm{ℝ}}^{n^2}$, the previous theorem guarantees us the following: ρ₀ is a critical point of the function E(ρ) = tr(ρH), surface restricted g ⁻¹(0) = {ρ ∈ M _n (ℝ); tr(ρ) − N _e = 0}, if, and only if, there is a real number µ that makes the gradient of the function null in ρ₀

$\Omega \left(\rho \right)=tr\left(\rho H\right)+\mu \left(tr\right(\rho )-N_e),$ (3.1)

as long as the g link gradient is not null at ρ. In addition, the g gradient can be easily obtained:

$\nabla \left(g\right(\rho \left)\right)=\nabla \left(\right(tr\left(\rho \right)-N_e\left)\right)=I^t=I,$

showing that the use of Lagrange Multipliers is allowed.

It is worth noting here that Vanderbilt’s group originally conceives scalar µ as a potential chemical in⁹. Thus, some procedure is necessary to update the chemical potential value during the search for the density matrix. In contrast, in our approach, the interpretation of µ as a Lagrange multiplier, which has the role of selecting arrays that reproduce the correct number of electrons, allows us to use the combination of computational implementation based on nonlinear conjugate gradient method, which was reported in⁵, requiring physical consistency for the number of electrons, to explicitly determine the µ multiplier as a function of ρ in each iteration step. Therefore, as we will see, in the above model, µ will not be an independent unknown. This same interpretation for µ is also reported in¹¹.

The constraint of idempotence was roughly and implicitly imposed by the so-called “McWeeny purification transformations”¹⁰, defined in the square matrix space, as follows. Consider the mapping

$f\left(x\right)=3x^2-2x^3.$

defined in square matrix space. It is evident that every idempotent matrix is a fixed point of this mapping, and it is easy to see that if λ is an eigenvalue of an array X, then f(λ) is eigenvalue of F(X), where f is the polynomial f(x) = 3x ² − 2x ³.

As we know, the eigenvalues of an idempotent matrix belong to the set {0, 1}. Let B be a ”approximately idempotent” matrix with an absolute maximum error of the order of O(ε), in the sense that that for a certain ε > 0, its eigenvalues are given by h or 1 + h, with h ∈ and 0 < |h| < ε. With that, eigenvalues of F(B) are f(h) or f(1 + h). From Taylor’s development of order 1 to f, we see that the matrix F(B) will be approximately idempotent with absolute maximum error of the order of O(ε²). Therefore, using F(ρ) instead of ρ in the equation (3.1) presents a proposal for the implicit and approximate introduction of the idempotency link in our model.

The mapping F is known as the MacWneey Purification Transformation¹⁰ and was incorporated into the Vanderbilt group model. Analyzing the purification polynomial f, we see that with, the eigenvalues of a matrix in the range [−0.5, 1.5], the eigenvalues of the purified matrix will remain in the range [0, 1].

Figure 1
The purifying polynomial f(x) = 3x ² − 2x ³.

With this, the function that will be minimized, unrestrictedly, is given by

$\Omega \left(\rho \right)=tr\left(\right[3{\rho }^2-2{\rho }^3\left]H\right)+\mu \left(tr\right(\rho )-N_e),$ (3.2)

where the multiplier µ will be defined next.

Still considering the matrix space identification and considering Ω also as a scalar field, we get that

$\nabla \Omega =3(H\rho +\rho H{)}^t-2(H{\rho }^2+\rho H\rho +{\rho }^{2}H{)}^t+\mu I,$ (3.3)

which is fundamental to the computational methods of unrestricted optimization.

3.2 Computational Strategies

In the basic algorithms of unrestricted scalar field optimization Ω(ρ) (steepest descent, Newton, Quasi-Newton, etc.) line search routines produce an iterative sequence that converges to the wanted stationary point

${\rho }_{k+1}={\rho }_k+\alpha D_k,$

where α is the step taken in the search direction D _k . In our context, we must require that in each iteration the number of electrons system is kept constant, that is, for every k, we have to tr(ρ_k+1 ) = tr(ρ_k ). So let’s impose, for every k, that

$tr\left(D_k\right)=0.$

Without going into too much detail, a nonlinear conjugate gradient method (NLCG) solves the problem of unrestricted minimization of a differentiable function Ω : ℝ^d → ℝ building a vector sequence using line search

${\rho }_{k+1}={\rho }_k+{\alpha }_k{D}_k,$ (3.4)

with search direction defined by

$D_k=\left\{\begin{array}{l}-\nabla \Omega \left({\overrightarrow{x}}_k\right),k=1\\ -\nabla \Omega \left({\overrightarrow{x}}_k\right)+{\beta }_k{\overrightarrow{d}}_{k-1},k\ge 2,\end{array}\right.$ (3.5)

where β_k are real numbers. There are several (not equivalent) proposals for the β_k scalar reported in the Numerical Optimization literature. Each proposal generates a different NLCG method. There are several proposals for choosing the NLCG step reported in the literature and this is still the subject of scientific discussions (see, for example,^{14), (2), (16} e¹⁷). Linear searches have an influence on theorems of global convergence for NLCG algorithms. The global convergence study for the NLCG still receives many current contributions. It is a β-dependent study. In addition, the β parameter is critical for updating direction d _k . There are several choices for the beta parameter, being the¹³ and⁶ the most classic ones. Fixing an NLCG algorithm, using the (descent) search direction definition given in (3.5) and our physical requirement that tr(D _k ) = 0, along with the gradient given in (3.3) we find that we must define

${\mu }_k=\frac{6}{N_e}\left[tr\right({\rho }_k^{2}H)-tr({\rho }_{k}H\left)\right].$ (3.6)

NLCG-based methods have been applied in Mathematical and Computational Modeling of Condensed Matter Physics phenomena since at least the 1980s (see¹⁵, for example). The conjugate gradient method (nonlinear) has some implementation possibilities with linear complexity. The canonical isomorphism between the matrix space and the Euclidean space ℝ^d has computational complexity of the order d ². This complexity can be reduced to linear complexity by using a cutting radius for the density matrix , as suggested in⁹, as follows. Assuming a ℝ^d implementation for NLCG containing d steps (asymptotically optimal algorithm), taking into account our identification between arrays and vectors, the DMTB-adapted version would use d ² steps, since, for the calculation of the system energy, in principle it would be necessary to compute all elements of the density matrix. However, it is known that

$\underset{R_{ij}\to \infty }{\lim }{\rho }_{ij}=0,$

where R _ij is the distance between the i and j orbitals (⁹). This convergence to zero has order of polynomial convergence in metals and exponential in insulators. Given this fact, Vanderbilt’s group proposed introducing a R _c > 0 cutting radius so that if R _ij > R _c then ρ_ij = 0. Thus, fixed an atom in the lattice, let L be the number of lattice atoms contained in the radius sphere R _c , centered on the atom that was previously fixed (due to lattice symmetry, L does not depend on the fixed atom then being a constant for each choice of cutting radius). With that, to form an element ρ_ij eventually not null, i may occupy any of the NM positions available, but j must range from LM positions. Hence, if we consider only these eventually non-null elements of the density matrix, our matrix vectoring process produces a vector with NLM ² coordinates. Having this vector as input, the CG algorithm would use NLM ² steps, thus presenting complexity O(N). Of course, R _c must be chosen to obtain density matrices close to true densities. In the implementation of⁵M = 4 (one s orbital and three p orbitals) was used.

4 FINAL CONSIDERATIONS

For this µ (equation (3.6)) construction to have the effect we want, it is critical that the initial ρ₀ matrix already has the correct number of electrons in the supercell. Thus, the expression

$\mu \left(tr\right(\rho )-N_e),$

that was already null in ρ₀, remains null throughout the process. So the expression of the functional Ω becomes just

$\Omega \left(\rho \right)=tr\left(\right(3{\rho }^2-2{\rho }^3\left)H\right).$

With this approach, the Lagrange multiplier has no effect on the functional objective Ω. Its effect, however, manifests itself in the gradient ∇Ω with each new iteration.

Therefore, the final problem is the unrestricted minimization (via nonlinear conjugate gradient algorithms) of

$\Omega \left(\rho \right)=tr\left(\right(3{\rho }^2-2{\rho }^3\left)H\right),$

using the equations (3.3) and (3.6).

Finally, Quantum Mechanics still requires us to have the commutativity of the product ρH, when ρ = ρ_e is the true density matrix. This property directly affects the choice of stopping criteria for our NLCG-based algorithm. In fact, assuming that the built arrays are becoming approximately idempotent enough that we can consider µ = 0, we see that the gradient becomes

$\nabla \Omega =(H\rho +\rho H-2\rho H\rho {)}^t.$

Of course, if ρ and H switch, then ∇Ω = 0. On the other hand, let’s assume that ∇Ω = 0. Then,

$H\rho +\rho H=2\rho H\rho .$

Multiplying this equation by ρ on the left, we see that

$\rho H\rho +{\rho }^{2}H=2{\rho }^{2}H\rho .$

Considering ρ₂ = ρ, we get ρH = ρHρ. Similarly, doing the right multiplication now gives Hρ = ρHρ. And so ρH = Hρ. With this, we see that ρ will be a critical point of Ω if, and only if, it switches with the Hamiltonian Tight-Binding matrix. This leads us to conclude that we should not use a maximum number of iterations in our algorithm. The stopping criterion should only be linked to the value tolerance ||∇Ω||.

Below is the pseudo-code based on the NLCG methodology for minimizing the functional Ω(ρ) used in⁵.

Algorithm 1
NLCG for energy functional (Polak-Ribière-type)

It is therefore an iterative process that builds the density matrix of the ground state. ρ_e .

These calculations are performed with parametrized tight-binding Hamiltonian for specific atomic species. It should be clear that, although based on NLCG, this computational approach does not fit the general definition of conjugate gradient methods, since the gradient of the field in question is updated with each iteration of the process. Therefore, the known global convergence theorems for this category of algorithms cannot be applied. Thus, the correctness of the algorithm (in the computers science sense) based on this proposal is a mathematical challenge to be solved. In addition to the issue of proposed iterative process for the minimization, the algorithm depends on the parametrization of the Hamiltonian that is used and, as we know, the accuracy and applicability of Hamiltonian parametrizations is quite relative. Therefore, the numerical methodology proposed here it lacks Numerical Analysis and we believe that the formalism developed is adequate to receive such treatment. However, the convergence and validation of the results of a particular case can be appreciated in⁵, where the implementation with silicon crystals was made rigorously along the lines of the present work, using the Kwon’s Hamiltonian⁸. The computational implementation was done in C++, under the object-oriented paradigm.

For silicon, many measurements, obtained experimentally and obtained by ab initio methods, are in agreement with the predictions obtained by DMTB. Without going into too much detail, the validation type proposed in the reference⁵ was as follows: we can take the experimental lattice parameter of Si crystals. It is known that this lattice parameter is worth approximately 5.43 angstroms. Therefore, the DMTB algorithm applied to silicon (ie, with a Si parametrized Tight-Binding Hamiltonian), when fed with this grid parameter value, must have full energy (E _tot ) lower than the other values for this parameter. More precisely, the total energy given by the equation (2.1), seen now as a function of the lattice parameter, should be a minimum of 5.43 Å. Since both the crystalline Hamiltonian matrix and E _rep are given in a parametrized form, this type of validation is not only linked to the DMTB, but also to the chosen Hamiltonian. Below is the graph obtained using Kwon’s Hamiltonian⁸.

Figure 2
Experimental Lattice Parameter Calculation.

For this graph, 10 values for energy were obtained from DMTB powered with lattice parameter values ranging from 5.20 to 5.65, from 0.5 Å to 0.5 Å. Then we adjusted the points by least squares to get a parabolic function. The minimum of the parabola occurs at 5.44 AA. Thus, we see that the relative error is approximately 0.2 %. For further details on this result, see⁵.

Still on the issue of computational complexity, we see that, even with the use of the cutting radius, the above algorithm has no linear complexity because, in line 13, we have the product trace of two matrices. However, there are well-established parallel programming proposals for this operation. In addition, manipulation methods for sparse arrays are required. Therefore, the final analysis of algorithmic complexity must be done under all these considerations.

REFERENCES

• ¹
  B. Amaral, T.A. Baraviera & M.O.T. Cunha. “Mecânica Quântica para Matemáticos em Formação”. IMPA (2011).
• ²
  P. Armand. Modification of the Wolfe Line Search Rules to Satisfy the Descent Condition in the Polak-Ribiere-Polyak Conjugate Gradient Method. Journal of Optimization Theory and Applications, 2 (2007), 132. doi: 10.1007/s10957-006-9123-7.
  » https://doi.org/10.1007/s10957-006-9123-7
• ³
  M.S. Daw. Model for energetics of solids based on the density matrix. Phys. Rev. B, 47(16) (1993), 10895-10898.
• ⁴
  P. Dirac. “The Principles of Quantum Mechanics”. Oxford University Press (1947).
• ⁵
  F.L.S. Ferreira, J.M.B. de Barros Filho & M.M. de Araújo. Gradiente Conjugado não Linear sobre um Ambiente “Matriz densidade Tight Binding ordem N” para o Estudo de Semi Condutores. In “Anais do ERMAC - Encontro Regional de Matemática Aplicada e Computacional / II Simpósio de Matemática Aplicada e Computacional da UFRRJ”, volume 3. SBMAC (2017), p. 1-7.
• ⁶
  R. Fletcher & C.M. Reeves. Function minimization by conjugate gradients. The Computer Journal, 7(2) (1964), 149-154. doi: 10.1093/comjnl/7.2.149. URL http://dx.doi.org/10.1093/comjnl/7.2.149
  » https://doi.org/10.1093/comjnl/7.2.149» http://dx.doi.org/10.1093/comjnl/7.2.149
• ⁷
  J.M. Jauch. “Foundation of Quantum Mechanics”. Addison-Wesley Publisching Company (1968).
• ⁸
  I. Kwon, R. Biswas, C.Z. Wang, K.M. Ho & C.M. Soukoulis. Transferable tight-binding models for silicon. Phys. Rev. B, 49(11) (1994), 7242-7250.
• ⁹
  X.P. Li, R.W. Nunes & D. Vanderbilt. Density-matrix electronic-structure method with linear system-size scaling. Phys. Rev. B, 47(16) (1993), 10891-10894.
• ¹⁰
  R. McWeeny. Some Recent Advances in Density Matrix Theory. Rev. Mod. Phys., 32(2) (1960), 335-369.
• ¹¹
  J.M. Millam & G.E. Scuseria. Liner scaling density matrix. J. Chem. Phys., 106(13) (1997), 5569-5577.
• ¹²
  A.T. Paxton. “An Introduction to the Tigth Binding aproximation-implementation by diagonalisation”. John von Neumann Institute for Computing (2009).
• ¹³
  E. Polak & G. Ribiere. Note sur la convergence de méthodes de directions conjuguées. ESAIM: Mathematical Modelling and Numerical Analysis - Modelisation Mathematique et Analyse Numerique, 3(R1) (1969), 35-43. URL http://eudml.org/doc/193115
  » http://eudml.org/doc/193115
• ¹⁴
  M. Sun & J. Liu. Three modified Polak-Ribiere-Polyak conjugate gradient methods with sufficient descent property. Journal of Inequalities and Applications, 2015(1) (2015), 125. doi: 10.1186/s13660-015-0649-9. URL http://dx.doi.org/10.1186/s13660-015-0649-9
  » https://doi.org/10.1186/s13660-015-0649-9» http://dx.doi.org/10.1186/s13660-015-0649-9
• ¹⁵
  M.P. Terter. Iterative minimization techniques for ab initio total-energy calculations: molecular dynamics and conjugate gradients. Phys. Rev. B, 40 (1989), 1225-1263.
• ¹⁶
  G. Yu, L. Guan & Z. Wei. Globally Convergent Polak-Ribière-Polyak Conjugate Gradient Methods Under a Modified Wolfe Line Search. Appl. Math. Comput., 215(8) (2009), 3082-3090. doi: 10.1016/ j.amc.2009.09.063. URL http://dx.doi.org/10.1016/j.amc.2009.09.063
  » https://doi.org/10.1016/ j.amc.2009.09.063» http://dx.doi.org/10.1016/j.amc.2009.09.063
• ¹⁷
  Y. Zhang, H. Zheng & C. Zhang. Global Convergence of a Modified PRP Conjugate Gradient Method. Procedia Engineering, 31 (2012), 986-995. doi: http://dx.doi.org/10.1016/j.proeng.2012.01.1131. URL http://www.sciencedirect.com/science/article/pii/S1877705812011551
  » https://doi.org/10.1016/j.proeng.2012.01.1131» http://www.sciencedirect.com/science/article/pii/S1877705812011551

Publication Dates

• Publication in this collection
  06 Sept 2021
• Date of issue
  2021

History

• Received
  15 Jan 2020
• Accepted
  22 Feb 2021