The Computation of Local Stress in ab initio Molecular Simulations

TL;DR

This paper develops formulas for calculating local stress in ab initio molecular dynamics models, specifically for Born-Oppenheimer and Ehrenfest dynamics, and compares these with full quantum simulations.

Contribution

It introduces new formulas for local stress computation in ab initio models, applicable to tight-binding and real-space methods, enhancing analysis accuracy.

Findings

Formulas for local stress are derived for different ab initio models.

Comparisons show the formulas closely match full ab initio simulations.

The work improves understanding of stress distribution in quantum-based molecular dynamics.

Abstract

Motivated by the increasingly more important role of ab initio molecular dynamics models in material simulations, this work focuses on the definition of local stress, when the forces are determined from quantum-mechanical descriptions. Two types of ab initio models, including the Born-Oppenheimer and Ehrenfest dynamics, are considered. In addition, formulas are derived for both tight-binding and real-space methods for the approximations of the quantum-mechanical models. The formulas are examined via comparisons with full ab initio molecular simulations.

Figures

Tables

Table 1. Table 1: This table shows the results in the verification of the properties ( 10 ) and ( 11 ).

Force Contributions	${𝒇}^{\mathrm{loc}}$	${𝒇}^{\mathrm{nloc}}$	${𝒇}^{i-i}$
$\frac{\sqrt{{\sum }_I{\left|{\sum }_J{𝒇}_{IJ}-{𝒇}_I\right|}^2}}{\sqrt{{\sum }_I{𝒇}_I^2}}$	0.0040	0.0043	1.2380$\times {10}^{-9}$
$\frac{\sqrt{{\sum }_{I,J}{\left|{𝒇}_{IJ}+{𝒇}_{JI}\right|}^2}}{\sqrt{{\sum }_{I,J}{\left|{𝒇}_{IJ}\right|}^2}}$	0.0284	0.0734	5.1974$\times {10}^{-23}$

Equations

$$\left\{\begin{array}[]{rl}\dot{R}_{I}&=P_{I}/m_{I},\\ \dot{P}_{I}&=-\frac{\partial}{\partial R_{I}}E_{\mathrm{tot}},\quad I=1,2,\cdots,N_{\mathrm{ion}}.\end{array}\right.$$

$$\hat{H}{\psi }_{\ell }={\epsilon }_{\ell }{\psi }_{\ell },\ell =1,2,\cdots ,N_{\mathrm{ele}}.$$

$$\hat{H}=-\frac{{\nabla }^2}{2}+{\hat{V}}_H[n]+{\hat{V}}_{\mathrm{xc}}[n]+{\hat{V}}_{ext}.$$

$${\hat{V}}_{KS}={\hat{V}}_H[n]+{\hat{V}}_{\mathrm{xc}}[n].$$

$$E_{\mathrm{tot}}=\sum _{\ell =1}^{N_{\mathrm{ele}}}{n}_{\ell }{\epsilon }_{\ell },{\mathbf{f}}_I=-\frac{\partial E_{\mathrm{tot}}}{\partial R_I},I=1,2,\cdots ,N_{\mathrm{ion}}.$$

$$i{\partial }_t{\psi }_{\ell }=\hat{H}{\psi }_{\ell },$$

$$\boldsymbol{j}(\boldsymbol{x},t)=\sum_{I=1}^{N_{\mathrm{ion}}}P_{I}(t)\varphi\big{(}\boldsymbol{x}-R_{I}(t)\big{)}.$$

$$\begin{array}[]{rcl}\partial_{t}\boldsymbol{j}&=&\displaystyle-\nabla_{\boldsymbol{x}}\cdot\Big{[}\sum_{I=1}^{N_{\mathrm{ion}}}P_{I}(t)\otimes P_{I}(t)\varphi\big{(}\boldsymbol{x}-R_{I}(t)\big{)}\Big{]}\\ &&+\displaystyle\sum_{I=1}^{N_{\mathrm{ion}}}\boldsymbol{f}_{I}(t)\varphi\big{(}\boldsymbol{x}-R_{I}(t)\big{)}.\\ \end{array}$$

$${\partial }_t\mathbf{j}={\nabla }_{\mathbf{x}}\cdot \mathbf{\sigma },$$

$${\mathbf{f}}_I=\sum _{J\ne I}{\mathbf{f}}_{IJ}$$

$${\mathbf{f}}_{IJ}=-{\mathbf{f}}_{JI}.$$

$$\begin{array}[]{rcl}&&\displaystyle\sum_{I=1}^{N_{\mathrm{ion}}}\sum_{J=1}^{N_{\mathrm{ion}}}\boldsymbol{f}_{IJ}\varphi\big{(}\boldsymbol{x}-R_{I}(t)\big{)},\\ &=&\displaystyle\frac{1}{2}\sum_{I=1}^{N_{\mathrm{ion}}}\sum_{J=1}^{N_{\mathrm{ion}}}\Big{[}\boldsymbol{f}_{IJ}\varphi\big{(}\boldsymbol{x}-R_{I}(t)\big{)}+\boldsymbol{f}_{JI}\varphi\big{(}\boldsymbol{x}-R_{J}(t)\big{)}\Big{]}\\ &=&\displaystyle\frac{1}{2}\sum_{I=1}^{N_{\mathrm{ion}}}\sum_{J=1}^{N_{\mathrm{ion}}}\boldsymbol{f}_{IJ}\Big{[}\varphi\big{(}\boldsymbol{x}-R_{I}(t)\big{)}-\varphi\big{(}\boldsymbol{x}-R_{J}(t)\big{)}\Big{]}\\ &=&\displaystyle-\nabla\cdot\frac{1}{2}\sum_{I=1}^{N_{\mathrm{ion}}}\sum_{J=1}^{N_{\mathrm{ion}}}\boldsymbol{f}_{IJ}\otimes R_{IJ}B_{IJ}.\end{array}$$

$$B_{IJ}(\boldsymbol{x})=\int_{0}^{1}\varphi\big{(}\boldsymbol{x}-R_{J}-\lambda R_{IJ}\big{)}d\lambda.$$

$$\boldsymbol{\sigma}=-\sum_{I=1}^{N_{\mathrm{ion}}}P_{I}(t)\otimes P_{I}(t)\varphi\big{(}\boldsymbol{x}-R_{I}(t)\big{)}-\frac{1}{2}\sum_{I=1}^{N_{\mathrm{ion}}}\sum_{J=1}^{N_{\mathrm{ion}}}\boldsymbol{f}_{IJ}\otimes R_{IJ}B_{IJ}(\boldsymbol{x}).$$

$${\psi }_{\ell }\approx \sum _I\sum _{\alpha \in occ(R_I)}{c}_{\ell ,I\alpha }{\varphi }_{\alpha }(\cdot ,R_I),\forall \ell =1,2,\cdots ,N_e,$$

$$\begin{array}[]{rl}S_{I\alpha,J\beta}=&\langle\phi_{\alpha}(\cdot,R_{I})|\phi_{\beta}(\cdot,R_{J})\rangle,\\ H_{I\alpha,J\beta}=&\langle\phi_{\alpha}(\cdot,R_{I})|\hat{H}|\phi_{\beta}(\cdot,R_{J})\rangle.\end{array}$$

$$M_{\alpha ,\beta }(R_{IJ})=M_{\beta ,\alpha }(R_{JI}).$$

$$H{\mathbf{c}}_{\ell }={\epsilon }_{\ell }S{\mathbf{c}}_{\ell }.$$

$$HC=SC\Lambda .$$

$$E_{\mathrm{tot}}=\sum _{\ell }\sum _{I\alpha }\sum _{j\beta }{H}_{I\alpha ,J\beta }{c}_{\ell ,I\alpha }{c}_{\ell ,J\beta }{n}_{\ell }.$$

$$\begin{array}[]{rl}\boldsymbol{f}_{I}=&-2\displaystyle\sum_{J\neq I}\sum_{\alpha\in occ(R_{I})}\sum_{\beta\in occ{R_{J}}}\frac{\partial}{\partial R_{I}}H_{I\alpha,J\beta}n_{\ell}c_{\ell,I\alpha}c_{\ell,J\beta}^{*}\\ &-2\displaystyle\sum_{J\neq I}\sum_{\alpha\in occ(R_{I})}\sum_{\beta\in occ{R_{J}}}\frac{\partial}{\partial R_{I}}S_{I\alpha,J\beta}n_{\ell}\varepsilon_{\ell}c_{\ell,I\alpha}c_{\ell,J\beta}^{*}.\end{array}$$

$$\begin{array}[]{lr}\boldsymbol{f}_{IJ}=&\displaystyle-2\sum_{\alpha\in occ(R_{I})}\sum_{\beta\in occ{R_{J}}}\frac{\partial}{\partial R_{I}}H_{I\alpha,J\beta}n_{\ell}c_{\ell,I\alpha}c_{\ell,J\beta}^{*}\\ &-2\displaystyle\sum_{\alpha\in occ(R_{I})}\sum_{\beta\in occ{R_{J}}}\frac{\partial}{\partial R_{I}}S_{I\alpha,J\beta}n_{\ell}\varepsilon_{\ell}c_{\ell,I\alpha}c_{\ell,J\beta}^{*}.\end{array}$$

$$\begin{array}[]{rcl}\rho_{I,J}&=&\displaystyle\sum_{\alpha\in occ(R_{I})}\sum_{\beta\in occ{R_{J}}}n_{\ell}c_{\ell,I\alpha}c_{\ell,J\beta}^{*},\\ \Gamma_{I,J}&=&\displaystyle\sum_{\alpha\in occ(R_{I})}\sum_{\beta\in occ{R_{J}}}n_{\ell}\varepsilon_{\ell}c_{\ell,I\alpha}c_{\ell,J\beta}^{*}.\end{array}$$

$$\boldsymbol{f}_{IJ}=-2\mathrm{tr}\big{(}\rho_{I,J}\frac{\partial}{\partial{R_{I}}}H_{I,J}\big{)}-2\mathrm{tr}\big{(}\Gamma_{I,J}\frac{\partial}{\partial{R_{I}}}M_{I,J}\big{)}.$$

$${\partial }_t{\mathbf{j}}_k+{\mathbf{\tau }}_{k+1/2}-{\mathbf{\tau }}_{k-1/2}=0.$$

$$E_{\mathrm{tot}}(R_1,R_2,\cdots ,R_N)=\sum _{\ell }{n}_{\ell }{\epsilon }_{\ell },$$

$$\frac{\partial }{\partial \lambda }{E}_{\mathrm{tot}}=\sum _{\ell }{n}_{\ell }⟨{\varphi }_{\ell }|\frac{\partial \hat{H}}{\partial \lambda }|{\varphi }_{\ell }⟩.$$

$$\begin{array}[]{rl}V_{\mathrm{ext}}(\boldsymbol{r})=&\displaystyle\sum_{I}\int w(\boldsymbol{r}-R_{I})n(\boldsymbol{r},t)d\boldsymbol{r}\\ +&\displaystyle\sum_{I}\sum_{\ell}\int\int\psi_{\ell}(\boldsymbol{r},t)^{*}U(\boldsymbol{r}-R_{I},\boldsymbol{r}^{\prime}-R_{I})\psi_{\ell}(\boldsymbol{r}^{\prime},t)d\boldsymbol{r}d\boldsymbol{r}^{\prime}\\ +&\displaystyle\frac{1}{2}\sum_{I=1}^{N}\sum_{J\neq I}\frac{Z_{I}Z_{J}}{|R_{I}-R_{J}|}.\end{array}$$

$$n(\mathbf{r},t)=\sum _{\ell =1}^{N_e}{n}_{\ell }|{\psi }_{\ell }(\mathbf{r},t){|}^2.$$

$${\mathbf{f}}_{IJ}^{\mathrm{i}-\mathrm{i}}=\frac{Z_I{Z}_J(R_I-R_J)}{|R_I-R_J{|}^3}.$$

Full text

Keywords: Ab initio Molecular Dynamics, Local stress.

1 Introduction

The last few decades have witnessed considerable interest in ab initio molecular dynamics (AIMD) as a simulation tool in chemistry and material science [32, 33, 8]. Compared to classical molecular dynamics (MD), where the interatomic potential is constructed in advance, AIMD determines on-the-fly the forces on the ions from quantum-mechanical (QM) degrees of freedom, which remain active throughout the computation. AIMD is particularly suitable for problems where classical MD models are difficult to construct, e.g., alloys with multiple species and materials interfaces, and problems where there are significant changes in the chemical environment, e.g., bonding patterns. There have been great recent progress in the development of efficient AIMD models, see the reviews [32, 43] for more details. Various software packages have also been developed and widely used [5, 30, 17].

Due to the advent of growing computing powers, it is conceivable that AIMD will soon be applied to large-scale mechanical systems, at which point, studying the elastic field induced by lattice defects, boundary condition, or electric field, would be of interest. While computing the total stress in QM models is relatively straightforward, see the monograph [31], to our knowledge, the computation of local stress has not investigated or implemented.

Meanwhile, for classical MD models, there have been consistent mathematical framework for the definition and computational methods to define the local stress, mostly motivated by the Irvine-Kirkwood formalism [23]. A popular approach is by Hardy [21, 22], where the local density and momentum distributions are represented by certain smooth kernel functions. Following the Newton’s equations of motion in the MD model, fundamental conservation laws can be derived, and the local stress can be identified from the momentum balance. There has been tremendous recent progress toward the definition and computation of such local stress, see [1, 4, 3, 10, 11, 12, 14, 13, 16, 15, 29, 37, 36, 40, 42, 46, 45, 48, 49, 50, 52, 51] and the references therein for recent development.

The primary goal of this paper is to establish the notion of local stress in AIMD models, where the electronic degrees of freedom are described by the density-functional theory [25] and its extensions [19, 44]. Due to the wide variety of AIMD models, we will focus on two popular methods, including the Born-Oppenheimer MD (BOMD) and the Ehrenfest dynamics. The BOMD approach involves solving the eigenvalue problem at each time step, while in Ehrenfest dynamics, the wave functions are determined by solving the time-dependent Schrödinger equation (TDSE). In both cases, the formulas for computing the total force on an ion are similar. Of particular importance for the definition of the local stress is the distribution of the force among each pair of ions. We show that the force decomposition critically depends on how the quantum-mechanical models are discretized in space. We discuss two specific implementations, the tight-binding (TB) approach, where the wave functions are approximation by linear combinations of atomic orbitals (LCAO), and the real-space approximations, e.g., finite-difference or finite element approximations. We will show that the force decomposition is rather straightforward in the former approach, since the interatomic distance is naturally involved in the Hamiltonian matrix. However, for real space methods, since the grid points are not attached to ion positions, such force decomposition is rather non-trivial.

The rest of the paper is organized as follows. In section 2, we review the AIMD models, as well as the definition of stress in classical MD models. The tight-binding approximation, along with the stress definition, will be presented in section 3. In section 4, we discuss the real space approximation methods, and we present some numerical results.

2 Ab initio Molecular Models and Fundamental Conservation Laws

Throughout the paper, we will denote the ion positions by $R_{I}$, and the coordinate for the electronic wave functions by $\boldsymbol{r}.$ The elastic field will be defined in terms of the spatial variable $\boldsymbol{x}.$ We assume that the system under consideration consists of $N_{\mathrm{ion}}$ ions, and $N_{\mathrm{ele}}$ electrons.

2.1 Two AIMD models

Let us start with the Newton’s equations of motion for the ions,

[TABLE]

Here $m_{I}$ is the mass; $R_{I}$ and $P_{I}$ represent the coordinate and momentum of the nuclei, respectively. The main departure of AIMD from conventional MD models is that the potential energy $E_{\mathrm{tot}}$ are determined from an underlying quantum-mechanical description, instead of using a pre-determined interatomic potential. How the electronic degrees of freedom are introduced, on the other hand, leads to different formalisms of AIMD. We will consider two models that have been widely implemented. Extensions to other models, e.g., [8, 38, 2, 26, 27, 28], especially the Car-Perrinello approach [8], are likely to be straightforward.

Born-Oppenheimer MD. First, we consider the Born-Oppenheimer MD (BOMD) model, where at each step, one solves the eigenvalue problem,

[TABLE]

Again we denote the number of electrons by $N_{\mathrm{ele}}$.

For the hamiltonian operator $\hat{H}$, we will consider the density-functional theory [25], which uses a noninteracting system with an exchange-correlation function that reproduces the electronic density of the full, interacting system. The Hamiltonian operator $\hat{H}$ (in reduced units) consists of the kinetic energy, Hartree and exchange correlation, and the external potential,

[TABLE]

We will denote the second and third terms, which are functionals of the electron density, as Kohn-Sham potential,

[TABLE]

The Hamiltonian operator $\hat{H}$ depends on the ions’ positions, which contributes to the electronic structure as an external potential $\hat{V}_{ext}$. This will be explained in more details in the next section. As a result, the eigenvalue problems need to be solved repeatedly, since the position of the ions is changing continuously. After each eigenvalue problem is solved, the total energy and force on each atom are computed as follows,

[TABLE]

The coefficients $n_{\ell}$ are the occupation numbers.

The coupled system (2) and (1) are the basis of BOMD models, which can be interpreted as differential-algebraic equations (DAE). The main assumptions is that the electronic states are instantaneously relaxed to ground states during the movement of the nuclei. Typically, the explicit formulas for computing the forces (5) are derived based on the Hellmann-Feymann theorem. We will discuss these formulas in the context of tight-binding and real-space approximations.

Ehrenfest Dynamics. Another type of AIMD models is the Ehrenfest dynamics, e.g., [27], where wave functions are determined from the time-dependent Schrödinger equations,

[TABLE]

in conjunction with the Newton’s equations of motion for the ions (1). The coupled dynamics describe an ion-electron coupling that occurs continuously in time. For the time-dependent model, we consider the time-dependent density-functional theory (TDDFT) [19, 39, 44].

In this case, the formula for computing the forces in (5) will be identically the same as in the static case, but for the reason that the wave functions and ion positions are independent variables .

2.2 Definition of local stress based on momentum balance and force decomposition

Since the Newton’s equation (1) remains in all AIMD models, the fundamental conservation law of momentum still holds. More specifically, given the momenta and coordinates of the nuclei, we define the local momentum,

[TABLE]

This is the starting point in Hardy’s derivation of stress [21], which we will follow here. The function $\varphi$ is introduced as an approximation to the dirac delta function in the original Irvine-Kirkwood formalism [23], and it is assumed to be non-negative with integral being 1. More specifically, by taking time derivative, one gets,

[TABLE]

The conservation of momentum at the continuum level asserts that

[TABLE]

with $\boldsymbol{\sigma}$ being the total stress tensor. So the goal is to express the right hand side of the equation (8) into a divergence form.

This will be true if the force $\boldsymbol{f}_{I}$ exhibits the following patterns

The force can be decomposed as:

[TABLE] 2. 2.

Each component is skew-symmetric,

[TABLE] 3. 3.

$\boldsymbol{f}_{IJ}$ is local: It should decay quickly as a function of the interatomic distance.

The first two conditions are critical (and sufficient) to ensure the conservation of momentum. It is worthwhile to point out, although it is trivial to many, that the fact that the force can be decomposed as in (10) does not by any means imply that the interaction is pairwise. In fact, the force component $\boldsymbol{f}_{IJ}$ may depend on atoms other than the two atoms $I$ and $J$. The last property is not mandatory. However, when the forces are highly nonlocal, the computation would be much more expensive, since many particles have to be visited.

To see why such force decomposition is necessary, we can substitute (10) into the second term on the right hand side of (8), which yields,

[TABLE]

Here in the last step, we used the fundamental theorem of calculus, and converted the difference between the two terms into a line integral along $R_{IJ}=R_{I}-R_{J}$,

[TABLE]

This derivation reveals a divergence term, from which one can identify the local stress,

[TABLE]

It is clear that such definition is always up to a divergence-free term, which makes it non-unique [1]. The momentum in the first term can be further decomposed into a mean ($\boldsymbol{j}(\boldsymbol{x},t)$) and relative momentum and they lead to the convection term, typically appears on the left hand side of the momentum balance equation, and a kinetic stress term, e.g., see the derivations in [36]. But our focus will be on the potential part of the stress (the last term).

In the next two sections, we will discuss how the force decomposition can be obtained when the force is determined from quantum-mechanical models, including the algebraic description (2) and the dynamics model (6).

3 Force Decomposition and Local Stress in Tight-binding models

3.1 Tight-binding methods

Tight-binding (TB) methods represent an important class of approximations that are constructed using basis functions centered at atoms, here denoted by $\phi_{\alpha}(\boldsymbol{r}-R)$. We also choose $\alpha\in occ(R_{I})$ with $occ(R_{I})$ indicating the set of local orbitals. The wave function $\psi_{\ell}$ will be approximated by a linear combination of the atomic orbitals (LCAO),

[TABLE]

The first step in implementing TB is assembling the overlap and Hamiltonian matrices, defined as follows,

[TABLE]

The matrix elements only depend on the relative position of the atoms; $R_{IJ}=R_{I}-R_{J}$. Namely, we can write $S_{I\alpha,J\beta}:=M_{\alpha,\beta}(R_{IJ})$. It satisfies the symmetry condition,

[TABLE]

In practice, they are precomputed in advance and represented in parametric forms.

For BOMD, the next step is usually the eigenvalue problem. Using a projection into the subspace spanned by the local orbitals, one can reduce the eigenvalue problem into a finite-dimensional generalized eigenvalue problem,

[TABLE]

Here $c_{\ell}$ is an eigenvector, with elements denoted by $c_{\ell,I\alpha}$

Let $C$ be the matrix that contains the eigenvectors as columns, then the diagonalization can be written simply as,

[TABLE]

The diagonal matrix $\Lambda$ contains the eigenvalues $\varepsilon_{\ell}$’s.

3.2 The force decomposition

In order to define the force on an ion, we first observe that the total energy, $E_{\mathrm{tot}}=\sum_{\ell}\varepsilon_{\ell}n_{\ell},$ can be written in terms of the Hamiltonian matrix and the eigenvectors,

[TABLE]

By taking the derivative with respect to the ion position $R_{I}$, and using the orthogonality condition in (18), one can derive the Hellmann-Feymann formula, which in this case, is given by,

[TABLE]

The derivative of the eigenvectors with respect to the ion positions does not appear due to the orthogonality conditions (18), but not because they are not dependent of them.

This provides a natural decomposition,

[TABLE]

This is exactly how the forces are evaluated in the DFTB+ model [5]. The formulas can be simplified, by introducing the density-matrix $\rho$, and the energy density-matrix $\Gamma$,

[TABLE]

Then the force can be written as,

[TABLE]

The fact that $\boldsymbol{f}_{IJ}$ is skew-symmetric is also evident.

3.3 A numerical test

To verify the desired property of the force decomposition, we consider a silicon nanowire with 512 atoms. The dimension of this quasi one-dimensional system is $96.96\AA\times 10.87\AA\times 10.87\AA$. The nanowire is divided into 16 block along the longitudinal direction. By integrating the conservation law (9) in a block, we have,

[TABLE]

Here $\boldsymbol{j}_{k}$ is the total momentum in the $k$th block. $\boldsymbol{\tau}$ is the projection of the stress in the longitudinal direction, i.e., $\boldsymbol{\tau}=\sigma\cdot(1,0,0)$, and it represents the traction between two adjacent blocks.

We follow the atomic units used in the software package, and modified the code to generate the force components $\boldsymbol{f}_{IJ}.$ We ran a BOMD simulation in DFTB+ for 100 steps with step size $\Delta t=2.067$ (0.05 femto-second). The initial velocity is randomly chosen. Fig. 1 shows the total momentum in a block in the longitudinal direction $j^{(1)}$, and a traverse direction $j^{(2)}(t)$. To examine the calculation of $\boldsymbol{f}_{IJ}$, we computed the traction $\boldsymbol{\tau}$ from the AIMD simulation, and then using the balance equation (24), we reconstructed the velocity at the same time steps. Excellent agreement has been found.

4 Stress Calculation in Real Space Models

In real space methods, the Hamiltonian operator will be discretized using finite difference [7] or finite element methods [35]. The Hamiltonian operator is then represented at the grid points as a matrix, here denoted by $H.$ The numerical procedure in the quantum mechanical model would involve either solving an eigenvalue problem (2) or the time-dependent Schrödinger equations (6). The main difficulty is that the electronic degrees of freedom are defined at grid point, and they are not specifically tied to the position of the ions.

4.1 Force decomposition and stress calculations

After the eigenvalues are computed, the total energy can be expressed as,

[TABLE]

with $n_{\ell}$’s being the occupation numbers.

The Hellmann-Feymann formula states that the change of the total energy with respect to the change of any parameter $\lambda$ is given by [6, 18],

[TABLE]

In real-space methods, the Hamiltonian operator is approximated at grid points. The only part of the Hamiltonian that explicitly depends on the nuclei position is the external potential,

[TABLE]

Here $n(\boldsymbol{r},t)$ is the electron density, given by,

[TABLE]

Remark: It is worthwhile to point out that in principle, for ground state calculations, the electron density has an implicit dependence on the ion position. However, the Hellmann-Feymann theorem asserts that the derivatives only need to taken with respect to the $R_{I}$ in the exterma; potential.

We will denote the three terms on the right hand side of (27) as $V^{\mathrm{loc}}$, $V^{\mathrm{nloc}}$, and $V^{\mathrm{i-i}}$, respectively. The term $V^{\mathrm{i-i}}$ is a classical Coulomb interaction and the corresponding force decomposition is trivial, since it has a classical pairwise form and no explicit involvement of the wave functions:

[TABLE]

Meanwhile, the first two terms stem from the pseudopotential approximation, e.g., norm-preserving potentials [20], which represents the influence of inner-shell electrons and reduce the problem to valence electrons. It represents the ion-electron interactions via a local term, here denoted by $w(\boldsymbol{r})$, which usually only depends on the distance from an ion, and a nonlocal term, which consists of projections to local orbitals. Here we have represented the summation of the projectors by the function $U$. Interested readers are referred to [20] for specific forms of these terms. The details of the force calculations in this case has been explained in the monograph [24].

It suffices to consider the derivatives of the first two terms with respect to the ion positions. The Hellmann-Feymann formula suggests to define,

[TABLE]

However, this formula does not reveal the interactions among the neighboring atoms, i.e., the decomposition $\boldsymbol{f}_{IJ}.$ It is tempting to separate the integrals by dividing the entire domain into non-intersecting subdomains, each of which contains a nucleus. However, our numerical tests suggest that the symmetry condition is not fulfilled. More importantly, in the nonlocal term, the projector $U$ only depends on points very close to the atom $R_{I}$. As a result, separating the integral would not yield a term that explicitly depends on a neighboring atom $R_{J}.$ In fact, when the second formula in (30) is implemented, the force becomes an on-site force, for which a force decomposition is not possible due to (10).

This difficulty has prompted us to examine the formula that determines the force and derive alternative expressions. In principle, the wave functions depend on the ion position as well, due to the external potential. Such dependence is implicit in the Feymann-Heymann theorem, since it is cancelled when the equation is left-multiplied by an eigen-function. To gain more insight, we observe the explicit dependence of the external potential on the relative positions to the atoms, $\{\boldsymbol{r}-R_{I},I=1,2,\cdots,\},$ and we will write the wave function as follows,

[TABLE]

Similarly we write the electron density as,

[TABLE]

It is, however, not necessary to know the explicit forms of the functions $\phi_{\ell}$ and $d.$ They are introduced here only to show explicitly the dependence on the ion positions. Now, we turn to the local force (30). Using integration by parts, we obtain,

[TABLE]

The last step is due to the form (31) and (32). This motivates us to define the force decomposition as follows,

[TABLE]

Similarly, we can define the force decomposition for the nonlocal part,

[TABLE]

The remaining issue is how to determine the change of the wave functions. This is done by using the density-functional theory perturbation (DFTP) approach [6]. Let us denote the change of the wave functions, due to an instantaneous change of the position of an ion in one direction, by $\delta\psi$. For the eigenvalue problem that needs to be solved in BOMD, the perturbation yields,

[TABLE]

for $\quad\ell=1,2,\cdots,N_{\mathrm{ele}}.$ Here $n_{0}(\boldsymbol{r})$ is the ground state electron density, and

[TABLE]

will correspond to the term $\frac{\partial}{\partial R_{J}}n(\boldsymbol{r},t)$ in (34). $\hat{V}_{KS}$ is the part of the Hamiltonian that depends on the electron density (4). The term $\delta\hat{V}_{ext}$ refers to the derivative of the external potential with respect to the ionic position. The second condition is to ensure the orthogonality of the wave functions.

The equations form a linear system for the wave functions, known as the Sternheimer equation, which can be solved via iterative methods. The coefficients of the linear system only depend on the ground states. So one does not need to solve the eigenvalue problem again. Such an approach has been an important route to determine phonon spectrum, polarizability, dielectric constants, etc [47].

At this point, it is clear that the formulas (34) and (35) would satisfy the property (10). But the skew-symmetric property (11) does not seem to be a direct consequence of our derivation. Here we offer a heuristic argument: suppose that one of the ions, $R_{J}$, undergoes an infinitesimally small displacement. In the Sternheimer equation (36), $\hat{V}_{ext}$ would be a function of $\boldsymbol{r}-R_{J}.$ As a result, we assume that the perturbed electron density will also be centered around $R_{J}$. We assume that it is asymmetric with respect to $\boldsymbol{r}-R_{J},$ e.g., $\delta n(R_{J}-\boldsymbol{r})=-\delta n(\boldsymbol{r}-R_{J})$. Therefore, we can write (34) as

[TABLE]

The local potential $w$ only depends on the relative distance. So by changing variables, $\boldsymbol{r}-R_{I}=R_{J}-\boldsymbol{r}^{\prime}$, we find that,

[TABLE]

The same argument can be made toward the nonlocal part of the pseudopotential (35) using a similar assumption on the density matrix.

For Ehrenfest dynamics models, one can again take the derivative of the wave function with respect to the ion position. This procedure yields,

[TABLE]

In practical implementations, one can solve this linear Schrödinger equation to determine $\delta\psi_{\ell},$ and implement the formulas (29), (34) and (35).

4.2 A numerical test

We consider a two-dimensional system – a single layer graphene sheet with 32 atoms. The computation is done by using OCTOPUS [30], a real space implementation of the ground state and time-dependent density-functional theory. Again, we use the standard atomic units. The length unit is Bohr radius and the energy unit is Hartree. All the results will be given in terms of these two units. For the computation, the grid size in the finite-difference approximation is chosen to be $0.2066,$ and the simulation domain consists of three-dimensional $90\times 78\times 40$ grid with center of the rectangular domain shifted to the origin. The dimension of the computational domain is $18.595\times 16.116\times 8.265$ (in Bohr radius). Periodic boundary conditions are applied in the first two space dimensions. We chose the standard pseudo-potential set in OCTOPUS.

To create a non-homogeneous stress, we set the initial displacement as follows,

[TABLE]

The atoms will be displaced according to this field in the first two dimensions with $\boldsymbol{x}$ as being the reference coordinates. We set the velocity to zero, so the stress is only due to the elastic field. We modified the part of the OCTOPUS code that computes the phonon spectrum. In particular, we followed the linear response calculations and computed $\delta\psi_{\ell}.$ The results are then used in to compute the force decomposition in (34) and (35).

In the first row of the table 1, we examine the force decomposition property (10). More specifically, we add up $\boldsymbol{f}_{IJ}$ and compare the sum with $\boldsymbol{f}_{I}$, for each of the three contributions, (29), (34) and (35). In the second row, we verify the skew-symmetric property (11), and again for each contribution. More specifically, we compute the Frobenius norm of the matrix $\boldsymbol{f}_{IJ}+\boldsymbol{f}_{JI}.$ Even though the system has non-homogeneous deformation, one can see that the two properties (10) and (11) are satisfied up to some small error. The error from the first row is clearly numerical error, since from the derivations of (34) and (35), the condition (10) should hold exactly. One may also choose to maintain the skew-symmetric property (11) exactly by defining $\tilde{\boldsymbol{f}}_{IJ}=(\boldsymbol{f}_{IJ}-\boldsymbol{f}_{JI})/2.$ In this case, the first condition (29) will hold approximately. Fig. 2 displays plots of the matrices $\boldsymbol{f}^{\mathrm{loc}}_{IJ}$ and $\boldsymbol{f}^{\mathrm{nloc}}_{IJ}$ to provide a more direct view of the skew-symmetric structure.

We now show in Fig. 3 how the magnitude of the force components $\boldsymbol{f}_{IJ}$ changes with respect to the distance between the two atoms ($r_{IJ}=|R_{IJ}|$). Interestingly, $\boldsymbol{f}_{IJ}$ exhibits a clear decay pattern, but it only becomes significantly smaller beyond the fifth neighbor, and negligible when the distance is around 10 Bohr (7th nearest neighbors). This is a lot larger than the cut-off radius of the Tersoff potential [41], also shown in the same figure. The force decomposition for the Tersoff potential can be found in [46]. In fact, there have been observations that the phonon spectrum of graphene will not be captured accurately, unless fifth neighbors included in the first-principle calculation [34].

Finally, we computed the local Hardy stress using the following two-dimensional kernel function, with a cut-off radius $r_{\mathrm{cut}}=4$ Bohr,

[TABLE]

We compared our results to that from the empirical Tersoff potential, and they exhibit very good agreement, as shown in Fig. 4. In principle, this should not be a means to validate the computational results. But in this particular case, since the displacement is quite smooth, we expect that the two results should be similar.

5 Summary and discussions

We have investigated in this paper the appropriate force decomposition in ab initio molecular dynamics models. The goal is to be able to define the local stress that can either be used to analyze the elastic field, or be integrated to a continuum description to build a multiscale method. For tight-binding models, it turns out that such decomposition is rather straightforward. For real-space methods, however, this becomes a non-trivial task. We argued that one must estimate the change in the wave functions with respect to the ion displacement in order to obtain such a force decomposition. This perturbation can be computed within the linear response framework by solving the Sternheimer equation.

6 Acknowledgments

This research was supported by NSF under grant DMS-1522617 and DMS-1819011. We would also like to thank the developers of DFTB+ [5] and OCTOPUS [9] to make the software packages available to our group.