# Total Energy and Forces: Ver. 1.2

Taisuke Ozaki, RCIS, JAIST

# Total energy for the collinear case

The OpenMX is based on density functional theories (DFT), the norm-conserving pseudopotentials, and local pseudo-atomic basis functions. The Kohn-Sham (KS) Bloch functions are expanded in a form of linear combination of pseudo-atomic basis functions (LCPAO) centered on site by
 (1)

where and are an expansion coefficient and pseudo-atomic function, a lattice vector, a site index, spin index, an organized orbital index with a multiplicity index , an angular momentum quantum number , and a magnetic quantum number . The charge density operator for the spin index is given by
 (2)

where means the integration over the first Brillouin zone of which volume is , and means the summation over occupied states. The charge density with the spin index is found as
 (3)

with a density matrix defined by
 (4)

Although it is assumed that the electronic temperature is zero in this notes, OpenMX uses the Fermi-Dirac function with a finite temperature in the practical implementation. Therefore, the force on atom becomes inaccurate for metallic systems or very high temperature. The total charge density is the sum of and as follows:
 (5)

Also, it is convenient to define a difference charge density for later discussion as
 (6)

where is an atomic charge density evaluated by a confinement atomic calculations associated with the site . Within the local density approximation (LDA) and generalized gradient approximation (GGA), the total energy of the collinear case is given by the sum of the kinetic energy , the electron-core Coulomb energy , the electron-electron Coulomb energy , the exchange-correlation energy , and the core-core Coulomb energy as
 (7)

The kinetic energy is given by
 (8)

The electron-core Coulomb energy is given by two contributions and related to the local and non-local parts of pseudopotentials:
 (9)

where and are the local and non-local parts of pseudopotential located on a site . Thus, we have
 (10)

 (11)

The electron-electron Coulomb energy is given by
 (12)

where is decomposed into two contributions and coming from the superposition of atomic charge densities and the difference charge density defined by
 (13)

 (14)

Within the LDA, the exchange-correlation energy is given by
 (15)

where is a charge density used for a partial core correction (PCC). The core-core Coulomb energy is given as repulsive Coulomb interactions among effective core charge considered in the generation of pseudopotentials by
 (16)

As discussed in the paper [1], for numerical accuracy and efficiency it is important to reorganize the sum of three terms , , and , as follows:
 (17)

where
 (18)

Therefore, we can reorganize these three terms as follows:
 (19)

 (20) (21) (22)

Following the reorganization of energy terms, the total energy can be given by
 (23)

# Kohn-Sham equation

Considering the orthogonality relation among one-particle wave functions, let us introduce a functional with Lagrange's multipliers :
 (24)

The variation of the functional with respect to the LCPAO coefficients yields the following matrix equation:
 (25)

Moreover, noting that the matrix for Lagrange's multipliers is Hermitian, and introducing diagonalizing the matrix, we can transform above equation as:
 (26)

By renaming by , and defining the diagonal element of to be , we have a well known Kohn-Sham equation in a matrix form as a generalized eigenvalue problem:
 (27)

where the Hamiltonian and overlap matrices are given by
 (28)

 (29)

with a Kohn-Sham Hamiltonian operator
 (30)

and the effective potential
 (31)

# Projector expansion of

The neutral atom potential is spherical and defined within the finite range determined by the cutoff radius of the confinement potential. Therefore, the potential can be expressed by a projector expansion as follows:
 (32)

where a set of radial functions is an orthonormal set defined by a norm for radial functions and , and is calculated by the following Gram-Schmidt orthogonalization [4]:
 (33) (34)

The radial function used is pseudo wave functions for both the ground and excited states under the same confinement potential as used in the calculation of the atomic electron density [2,3]. The most important feature in the projector expansion is that the deep neutral atom potential is expressed by a separable form, and thereby we only have to evaluate the two-center integrals to construct the matrix elements for the neutral atom potential. As discussed later, the two-center integrals can be accurately evaluated in momentum space. For details of the projector expansion method, see also the reference [1].

# Two-center integrals

The overlap integrals, matrix elements for the non-local potentials, for the neutral atom potentials in a separable form discussed in the Sec. 3, and for the kinetic operator consist of two-center integrals. In this section, the evaluation of the two-center integrals is discussed. The pseudo-atomic function in Eq. (1) is given by a product of a real spherical harmonic function and a radial wave function :
 (35)

where means Euler angles, and , for , and radial coordinate. Although the real function is used for the spherical harmonic function in OpenMX instead of the complex function , firstly we consider the case with the complex function as
 (36)

After the evaluation of two-center integrals related to the complex function, they are transformed into two-center integrals for the real function by matrix operations. Then, the pseudo-atomic function given by Eq. (36) can be Fourier-transformed as follows:
 (37)

The back transformation is defined by
 (38)

Noting the Rayleigh expansion
 (39) (40)

and putting Eqs. (36) and (40) into Eq. (37), we have
 (41)

where we defined
 (42)

Using Eqs. (38), (40), and (41), the overlap integral is evaluated as follows:
 (43)

where an integral in the third line was evaluated as
 (44)

and one may notice Gaunt coefficients defined by
 (45)

The matrix element for the kinetic operator can be easily found in the same way as for the overlap integral as follows:
 (46)

# Discretization

The two energy terms and are discretized by a regular mesh in real space [5], while the integrations in , , , and can be reduced to two center integrals which can be evaluated in momentum space. The regular mesh in real space is generated by dividing the unit cell with a same interval which is characterized by the cutoff energies , , and :
 (47)

 (48)

 (49)

 (50)

where , , and are unit cell vectors of the whole system, and , , and are the number of division for , , and -axes, respectively. , , and are determined so that the differences among , , and can be minimized starting from the given cutoff energy. Using the regular mesh , the Hartree energy associated with the difference charge density is discretized as
 (51)

The same regular mesh is also applied to the solution of Poisson's equation to find . Then, the charge density is Fourier transformed by
 (52)

From , we can evaluate by
 (53)

Thus, we see that is also written in a discretized form with the regular mesh as follows:
 (54)

Within the LDA, can be easily discretized as well as by
 (55)

For the GGA, is discretized with the gradient of charge density evaluated with a finite difference scheme in the same way in the LDA. Since the derivative of the charge density with respect to is given by
 (56)

the matrix elements for and in the Kohn-Sham equation Eq. (27) are found by differentiating the energies and with respect to as follows:
 (57)

and
 (58)

where the quantities in the parenthesis correspond to the matrix elements.

# Force on atom

 (59)

The derivative of the kinetic energy with respect to is given by
 (60)

The derivative of the neutral atom potential energy with respect to is given by
 (61)

The derivative of the non-local potential energy with respect to is given by
 (62)

The derivative of the Hartree energy energy for the difference charge density with respect to is given by
 (63)

Considering Eq. (54) and
 (64)

we have
 (65)

and
 (66)

Moreover, the derivative of with respect to is given by
 (67)

The derivative of with respect to is simply given by
 (68)

The derivative of with respect to the atomic coordinate is easily evaluated by
 (69)

The derivative of the screened core-core Coulomb energy with respect to is given by
 (70)

Since the second term is tabulated in a numerical table as a function of distance due to the spherical symmetry of integrands, the derivative can be evaluated analytically by employing an interpolation scheme. The derivatives given by Eqs. (60), (61), (62), (63), and (69) contain the derivative of LCPAO coefficient . The derivative of can be transformed to the derivative of the overlap matrix with respect to as shown below. By summing up all the terms including the derivatives of in Eqs. (60), (61), (62), (63), and (69), we have
 (71)

where is a diagonal matrix consisting of Heaviside step functions. Noting that Eq. (27) can be written by and in a matrix form, the product of two diagonal matrices is commutable, and for any square matrices, Eq. (71) can be written as
 (72)

Moreover, taking account of the derivative of the orthogonality relation with respect to , we have the following relation:
 (73)

Putting Eq. (73) into Eq. (72), we have

where the energy density matrix is given by
 (74)

The terms including the derivative of matrix elements in Eqs. (60), (61), and (62) can be easily evaluated by
 (75)

where
 (76) (77) (78)

The derivatives of these elements are evaluated analytically from the analytic derivatives of Eqs. (43) and (46). The remaining contributions in first terms of Eqs. (63) and (69) are given by
 (79)

The second terms in Eqs. (63) and (69) becomes
 (80)

Thus, we see that the derivative of the total energy with respect to the atomic coordinate is analytically evaluated for any grid fineness as
 (81)

## Bibliography

1
T. Ozaki and H. Kino, Phys. Rev. B 72, 045121 (2005).

2
T. Ozaki, Phys. Rev. B 67, 155108 (2003).

3
T. Ozaki and H. Kino, Phys. Rev. B 69, 195113 (2004).

4
P. E. Blochl, Phys. Rev. B 41, R5414 (1990).

5
J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon, and D. Sanchez-Portal, J. Phys.: Condens. Matter 14, 2745 (2002).

2008-11-23