# Configurational forces in electronic structure calculations using Kohn-Sham density functional theory

###### Abstract

We derive the expressions for configurational forces in Kohn-Sham density functional theory, which correspond to the generalized variational force computed as the derivative of the Kohn-Sham energy functional with respect to the position of a material point x. These configurational forces that result from the inner variations of the Kohn-Sham energy functional provide a unified framework to compute atomic forces as well as stress tensor for geometry optimization. Importantly, owing to the variational nature of the formulation, these configurational forces inherently account for the Pulay corrections. The formulation presented in this work treats both pseudopotential and all-electron calculations in single framework, and employs a local variational real-space formulation of Kohn-Sham DFT expressed in terms of the non-orthogonal wavefunctions that is amenable to reduced-order scaling techniques. We demonstrate the accuracy and performance of the proposed configurational force approach on benchmark all-electron and pseudopotential calculations conducted using higher-order finite-element discretization. To this end, we examine the rates of convergence of the finite-element discretization in the computed forces and stresses for various materials systems, and, further, verify the accuracy from finite-differencing the energy. Wherever applicable, we also compare the forces and stresses with those obtained from Kohn-Sham DFT calculations employing plane-wave basis (pseudopotential calculations) and Gaussian basis (all-electron calculations). Finally, we verify the accuracy of the forces on large materials systems involving a metallic aluminum nanocluster containing 666 atoms and an alkane chain containing 902 atoms, where the Kohn-Sham electronic ground state is computed using a reduced-order scaling subspace projection technique [40].

^{†}

^{†}preprint:

## I Introduction

Electronic structure calculations based on density functional theory (DFT) have played an important role in providing significant insights into a wide variety of materials properties. The Kohn-Sham approach [1; 2] to DFT provides an efficient formulation to compute the ground-state properties of materials systems by reducing the many-body Schrödinger problem of interacting electrons into an equivalent problem of noninteracting electrons in an effective mean field that is governed by the electron-density. While this effective single-electron formulation has no approximations, and is exact, in principle, the quantum mechanical interactions between electrons manifest in the form of an unknown exchange-correlation functional, which is modeled in practice. While improving the accuracy of the exchange-correlation functional is still an active area of research, the widely used models for the exchange-correlation functional [3] have shown to predict a range of materials properties across various materials systems with good accuracy.

An important aspect of electronic structure calculations using DFT is the determination of accurate quantum mechanical forces on the atomic nuclei and stress tensor in periodic systems. The forces and stresses are used in many aspects of electronic structure calculations, including geometry optimization, calculation of the dynamical matrix and phonon spectra, determination of elastic properties and ab-initio molecular dynamics calculations. In DFT, the atomic force on a nuclei can be computed by taking recourse to the Hellmann-Feynman theorem [4; 5], which relates the force acting on a nucleus to the electrostatic field at the nucleus due to the electron-density and other nuclei [6]. However, in practical DFT calculations, the variational force can differ from the Hellmann-Feynman force if the numerical basis sets employed in the DFT calculations explicitly depend on atomic positions [7; 8; 9] (Pulay force correction due to incomplete basis-set), and, moreover, the Hellmann-Feynman theorem does not provide elastic stresses on periodic systems. Stress is computed as the derivative of the energy with respect to strain, and the practical numerical implementation of these derivatives relies very much on the particular basis set type employed in the DFT calculation and has to account for possible strain dependence on the basis sets employed, which manifest as Pulay stresses. To this end, there have been significant efforts in the past few decades to implement variationally consistent forces [10; 11; 12; 13; 14; 16; 15; 17; 18; 19; 20; 21; 22; 23; 24; 25] and stresses [26; 27; 28; 29; 30; 31; 32; 33; 16] in ab-initio codes employing various basis sets.

We note that the general approach to computing atomic forces or stresses in the existing numerical implementations rely on the outer variations of the Kohn-Sham energy functional with respect to the change in the position of atoms (for computing forces) or the lattice vectors (for computing stresses). Thus, most of the approaches in the literature rely on separate formulations for computing forces and stresses. In this work, we propose a configurational force approach to compute both atomic forces and stresses in periodic systems under a single framework using Kohn-Sham density functional theory. The idea of configurational force dates back to Eshelby’s formulation [34], and is widely used in continuum physics [35] to study problems involving material inhomogenities like coherent phase interfaces, evolving surfaces, cracks etc., in materials. In contrast to the existing methods employed in electronic structure calculations, the proposed approach is based on the inner variations of the Kohn-Sham energy functional resulting in configurational forces. These forces correspond to the generalized variational derivative of the Kohn-Sham energy functional with respect to the position of a material point x. Hence, the proposed approach provides a generalized force with respect to internal positions of atoms as well as the periodic cell, in the case of periodic calculations, thus providing a unified framework to compute atomic forces as well as stresses for geometry optimization. Moreover, the proposed formulation allows us to treat both pseudopotential and all-electron calculations in a single framework. Furthermore, we note that the configurational forces inherently account for Pulay corrections owing to the variational nature of the formulation, and no separate treatment is required to account for these corrections.

The proposed configurational force approach provides a natural way to evaluate atomic forces and stresses in the context of finite-element discretization. The finite-element basis [36], a piecewise polynomial basis, offers some unique advantages for Kohn-Sham DFT calculations. Finite-element basis can handle arbitrary boundary conditions, which enables the consideration of isolated, semi-periodic, as well as, periodic systems under a single framework. Further, the finite-element basis is amenable to adaptive spatial resolution, which can be exploited for handling all-electron DFT calculations [37; 38; 39; 40; 41] as well as the development of coarse-graining techniques [42; 43; 44]. Moreover, the locality of these basis functions provides good scalability on parallel computing platforms. Though significant efforts have been focused in the recent past to develop finite element based methodologies for conducting Kohn-Sham DFT calculations [45; 46; 47; 48; 49; 50; 51; 52; 53; 54; 55; 56; 37; 38; 57; 58; 59; 60; 61; 62; 39; 40; 41; 63], most of these efforts have been confined to computing the electronic ground-state. To the best of our knowledge, expressions for computing both stresses and forces that are variationally consistent have not been developed for the finite-element discretization. The proposed work fills this important gap to enable large-scale electronic structure calculations using finite-element basis. Further, the developed configurational forces also provide the generalized forces corresponding to the nodal positions of the underlying finite-element triangulation, which can in turn be effectively used to determine the optimal positions of the finite element nodes— -adaptivity, an a posteriori mesh adaption technique [64].

In this work, the proposed configurational force approach is presented in the context of spectral finite-element discretization. However, the method presented is general, and applicable to any real-space basis employed in the solution of the Kohn-Sham DFT problem. Keeping in mind the reduced-order scaling methods [65; 40], the formulation has been developed using non-orthogonal electronic wavefunctions. We begin by considering the formulation of the Kohn-Sham density functional theory at finite temperature in terms of non-orthogonal wavefunctions, where the electronic ground state, for a given position of atoms, is governed by a variational problem involving the minimization of the Kohn-Sham free energy functional with respect to the wavefunctions and the density matrix subject to the constraint on the total number of electrons in the system. Subsequently, we discuss the local variational real-space formulation of Kohn-Sham DFT, where the computation of the Kohn-Sham electronic ground-state for a given position of atoms can be formulated as a local variational problem involving wavefunctions, density matrix and electrostatic potentials. We then derive the expressions for the configurational forces corresponding to geometry optimization. Using the approach of inner variations, we evaluate the generalized forces from the locally reformulated energy functional corresponding to perturbations of the underlying space. We note that the final form of the configurational force expression involves integrals of Eshelby tensors contracted with the gradient of the generator () associated with the perturbation of the underlying space. The derived configurational forces provide a unified framework to compute both forces on atoms and stresses in periodic systems. To elaborate, the force on any given atom is given by the configurational force computed using a generator whose compact support includes only the atom of interest, while the stress tensor is computed by using a generator corresponding to affine deformations of the underlying space. Subsequently, we present the details involved in the numerical implementation of the configurational force expressions within the framework of finite-element discretization, and thereby discuss a computationally efficient strategy to conduct geometry optimization using finite-element basis.

Finally, we present the accuracy of the proposed configurational force approach on sample benchmark problems. To this end, we consider both non-periodic and periodic material systems, and conduct both pseudopotential and all-electron calculations. We begin by examining the rates of convergence of the finite-element approximation in the computed forces and stresses. We observe a close to convergence in all benchmark problems. Further, in order to asses the variational nature of the computed forces and stresses, we compare these with those obtained from finite differencing energies, which are in excellent agreement. As a further test, we compute the forces and stresses as a function of bond length or lattice parameter and compare these with the derivative of polynomial fits to the energy dependence on these parameters. These extensive tests have ascertained the variational nature of the configurational forces in our numerical implementation. We also compare the accuracies of forces and stresses obtained from the proposed configurational force approach with those obtained using external DFT packages, and find very good agreement. Finally, we consider large materials systems comprising of an aluminum nanocluster containing unit-cells ( atoms) and an alkane chain containing atoms, where the finite-element discretized Kohn-Sham problem is solved using the reduced order scaling subspace projection method [40] involving non-orthogonal localized wavefunctions. We find that the configurational forces computed using the reduced-order scaling method are in very good agreement with those computed without the additional approximations used in the reduced-order scaling method, with differences being well within the desired chemical accuracy.

The remainder of the paper is organized as follows: Section II describes the formulation of the Kohn-Sham DFT problem in the context of non-orthogonal wavefunctions, followed by the local real-space reformulation of electrostatic interactions in Section III. The configurational force expressions for both non-periodic and periodic cases are presented in Section IV, with the details of the derivation presented in the Appendix. Section V describes the numerical implementation of forces and stresses in the context of finite-element discretization. Section VI presents the results on various systems demonstrating the accuracy and variational nature of the computed forces and stresses. Finally, we conclude with a summary and outlook in Section VII.

## Ii Kohn-Sham density functional theory

Consider a materials system consisting of electrons with nuclei whose position vectors are denoted by . Neglecting spin, the free energy of the system in Kohn-Sham density functional theory(DFT) [1; 2] at finite temperature [66] can be written as

(1) |

where () denotes the set of electronic wavefunctions. In all generality, we assume that these wavefunctions are non-orthogonal in order to realize a reduced-order scaling numerical implementation of DFT [65; 40]. Here denotes the matrix corresponding to the single particle density operator () expressed in the non-orthogonal basis , i.e. , where are the matrix elements of the inverse of the overlap matrix given by

(2) |

We note that the superscript * in the above, and all the equations subsequently, denotes the complex conjugate. The electron-density in equation (1) can be expressed in terms of the density matrix and the non-orthogonal wavefunctions as

(3) |

Further, , which denotes the kinetic energy of non-interacting electrons, is given by

(4) |

The exchange-correlation energy, which incorporates all the quantum-mechanical interactions, is denoted by . While the explicit form of remains elusive, various approximations have been developed over the past decades, with the local density approximation (LDA) [6] and generalized gradient approximation (GGA) [67; 68] being adopted across a range of materials systems. In the present work, LDA exchange-correlation energy is adopted which has the following functional form:

(5) |

where . In particular, we employ the Slater exchange and Perdew-Zunger [69; 70] form of the correlation functional.

The electrostatic interaction energy, , represents the interactions between electrons and nuclei, which can be further decomposed as

(6) |

where the first term is the Hartree energy representing the electrostatic interaction energy between electrons, and denotes the repulsive energy between nuclei. These are given by

(7) |

with denoting the charge on the nucleus. in equation (6) denotes the classical interaction energy between electrons and nuclei, and is given by

(8) |

As chemical bonding in many material systems is not influenced by the tightly bound core electrons close to the nucleus of an atom, these core electrons may not play a significant role in governing many materials properties. Hence, the pseudopotential approach is commonly adopted, where only wavefunctions for the valence electrons are computed in an effective potential of the nucleus and core electrons given by a pseudopotential. The pseudopotential is often defined by the operator , where is the local part of the pseudopotential operator and is the nonlocal part of the operator. In pseudopotential Kohn-Sham DFT, is given by

(9) |

The norm-conserving Troullier-Martins pseudopotential [71] in the Kleinman-Bylander form [72] is employed in this work, where the action of the operator on the wavefunction is given by

(10) | |||

(11) |

with

In the above, is the local potential of atom , denotes the pseudopotential component of atom corresponding to the azimuthal quantum number , and is the corresponding single-atom pseudo-wavefunction with azimuthal quantum number and magnetic quantum number .

In a non-periodic setting, representing an isolated atomic system, all integrals in equations (4)-(11) are over and the summations include all atoms in the system. In the case of an infinite periodic crystal, all integrals involving x in equations (4)-(11) are over the unit-cell, whereas the integrals involving y are over . Further, the summation over is on atoms in the unit-cell, and the summation over extends over all lattice sites.

The electronic entropy contribution, in equation (1), is given by

(12) |

where with denoting the Boltzmann constant and denoting the electronic temperature.

Finally, the ground state in Kohn-Sham DFT is governed by the variational problem

(13) |

where with denoting the Lagrange multiplier (Fermi-energy) enforcing the constraint on number of electrons, i.e., . Equation (13) indicates that the electronic ground-state needs to be computed for every configuration of the nuclei encountered during the minimization procedure over the nuclear positions.

## Iii Real-space formulation

In this section, we present the local variational real-space formulation of Kohn-Sham DFT, which is subsequently used to derive the expressions for configurational forces associated with internal atomic relaxations and cell relaxation. We note that the various components of the electrostatic interaction energy (cf. equation (6)) are non-local in real-space. These extended interactions are typically computed in Fourier space, using Fourier transforms. However, Fourier space formulations employing plane-waves provide only uniform spatial resolution and restrict simulation domains to periodic geometries and boundary conditions. Thus, they are not well suited for material systems involving molecules, nano-clusters, or, systems containing defects. Further, the non-locality of the plane-wave basis in real-space limits the scalability of computations on parallel computing platforms.

The real-space formulation discussed here is devoid of the aforementioned limitations of a Fourier space formulation. The proposed formulation extends the recent efforts in the local reformulation of Kohn-Sham DFT [51; 57; 39], and differs from prior efforts in the way the extended electrostatic interactions have been treated. In particular, the proposed formulation is crucial to developing a unified framework to compute the configurational forces associated with both atomic relaxations and cell relaxations, discussed subsequently. We note that the treatment of extended electrostatic interactions is similar to Das et al. [73], which was proposed in the context of orbital-free DFT formulation. Here, we provide the relevant details in the context of Kohn-Sham DFT.

Let denote a regularized Dirac distribution located at . Thus, the charge distribution of the nuclear charge is given by . We define the nuclear charge distribution and to reformulate the repulsive energy as

(14) |

where denotes the self energy of the nuclear charges which depends only on the nuclear charge distribution and can be expressed as

(15) |

with denoting the electrostatic potential corresponding to the nuclear charge and is given by

(16) |

Since the kernel corresponding to extended electrostatic interactions is the Green’s function of the Laplace operator, the evaluation of self energy as well as the potential can be reformulated as the following local variational problem:

(17) |

with being the minimizer of the variational problem in equation (17) associated with the nucleus. In the preceding equation, denotes the Hilbert space of functions such that the functions and their first-order derivatives are square integrable on . Further, the local part of in equation (9) can be rewritten as

(18) |

We note that in (and in ) denote the valence charges of the (and ) nucleus in the case of pseudopotential calculations. While, for all-electron calculations, they denote the atomic number with . The kernel corresponding to the extended electrostatic interactions in the expressions for , and is the Green’s function of the Laplace operator. Thus, can be rewritten as the following local variational problem:

(19) |

where denotes the trial function for the total electrostatic potential due to the electron-density and the nuclear charge distribution, and is a suitable function space corresponding to the boundary conditions of the problem, discussed subsequently. Using the local reformulation in equations (17) and (III), the electrostatic interaction energy in DFT can now be expressed as the following variational problem:

(20) |

where

with

(21) |

(22) |

with denoting the vector containing the trial electrostatic potentials corresponding to all nuclear charges in the simulation domain. Further, the minimization over in equation (20) refers to a simultaneous minimization over all these electrostatic potentials in . We further note that in the case of all-electron Kohn-Sham DFT calculations. Thus, this local reformulation also provides a unified framework for both pseudopotential as well as all-electron DFT calculations.

Finally, using the local reformulation of the extended electrostatic interaction energy, the computation of the electronic ground-state for a given position of atoms in equation (13) can be formulated as the following local variational problem in wavefunctions, density matrix and electrostatic potentials:

(23) |

(24) | ||||

(25) |

In the above, denotes a suitable function space that guarantees the existence of minimizers. Further, we remark that numerical computations involve the use of bounded domains, which in non-periodic calculations correspond to a large enough domain containing the compact support of the wavefunctions, and, in periodic calculations, correspond to the supercell. Denoting such an appropriate bounded domain by subsequently, in the case of non-periodic problems, and in the case of periodic problems. We note that, in practice, the solution to the variational problem in equation (23) is computed by taking recourse to the Kohn-Sham equations, which constitutes a non-linear eigenvalue problem solved using self-consistent field iteration on the electron density (cf. [57; 39]).

## Iv Configurational forces

We now derive the expressions for the configurational forces corresponding to geometry optimization. To this end, we employ the approach of inner variations, where we evaluate the generalized forces corresponding to perturbations of underlying space, which provides a unified expression for the generalized force corresponding to the geometry of the simulation cell—internal atomic positions, as well as, the simulation domain boundary. Further, owing to the real-space formulation in Section III, the derived expressions are applicable to both pseudopotential and all-electron calculations.

In order to derive the expressions for the configurational forces, we first define a bijective mapping which represents the infinitesimal perturbation of the underlying space, mapping a material point x to a new point such that . Further, we define the generator of this mapping as . We note that the mapping should be constrained to rigid body deformations in the compact support of the regularized nuclear charge distribution in order to preserve the integral constraint . To this end, in the compact support of for , where is unitary, and , are independent of x denoting rotation and translation operations, respectively.

### iv.1 Non-periodic DFT calculations

We first discuss the case of non-periodic problems with bounded domains , where is large enough that the values of wavefunctions are negligible outside of , i.e., we assume the wavefunctions have a compact support on . We will subsequently discuss the case of periodic calculations with periodic boundary conditions. In view of the finite-element basis employed to discretize the DFT functional subsequently, we recast in equation (24) as the following:

In the above, we employ the divergence theorem and make use of the fact that for non-periodic problems.

Let us consider the perturbation of the underlying space given by , which maps a material point x in to in (which denotes the image of under the perturbation). The ground-state energy in the perturbed space is given by

(26) |

where , , and are the solutions of the saddle point problem (23) on the perturbed space. We now evaluate the configurational force by computing the Gâteaux derivative of , i.e., , to arrive at (cf. Appendix for details of the derivation):

(27) |

where and denote Eshelby tensors whose expressions in terms of the solutions of the saddle point problem (23) on the original space (, , , ) are provided below. We note that in the above expressions, and subsequently, the outer product between two vectors is denoted by “”, the dot product between two vectors by “.” and dot product between two tensors by “:”. Dropping the superscript on the electronic fields for notational convenience, the Eshelby tensors are given by:

where

Further, the other terms in equation (27) are given by

where

We note that the terms and do not appear in the expressions for and , respectively, owing to the restriction that corresponds to rigid body deformations in the compact support of the nuclear charge distribution . Hence, in the compact support of , we have . We note that the second term in equation (27) involves an integral over which appears intractable. However, we can split this integral on a bounded domain , containing the compact support of , and its complement which can in turn be computed as a surface integral. Thus,

which can be written as the following equation:

(28) |

where denotes the outward normal to the surface . The last equality follows from the fact that on . Though the expression for configurational force in equation (27) is derived in the case of pseudopotential calculations, the expression for all-electron calculations is obtained by using , and . Finally, we note that the force on any given atom is computed by choosing such that its compact support only contains the atom of interest. We refer to the Appendix for the detailed derivation of equation (27).

The computation of electronic ground-state in equation (23) is equivalent to finding the occupied eigen subspace spanned by the eigenfunctions (canonical wavefunctions) corresponding to the smallest eigenvalues of the Kohn-Sham self-consistent Hamiltonian. Consequently, if represent the orthonormal canonical eigenfunctions spanning the occupied subspace of the Kohn-Sham Hamiltonian, the expressions for the Eshelby tensor E and in equation (27) reduce to the following form in terms of , (orbital occupancy functions), (Kohn-Sham eigenvalues) and (electrostatic potential):

(29) |

where

Further,

where