general, accurate, and computationally efficient method to evaluate NAC is desirable.

In this Letter, we describe a concentric approximation (CA) for the evaluation of NAC within the PAW method and demonstrate that it reproduces the AE NAC to within 0.1 meV at the computational expense of the PS NAC. We demonstrate that small displacements in the atom center of partial waves in the core region can be ignored in the evaluation of NAC over an MD time step dt that is as large as 10 fs. The present work merges the efficiency of the PS NAC and the accuracy of the exact AE NAC. The method is validated with two prototypical periodic condensed matter systems containing transition metals, bulk TiO 2 and Cu, and a molecule/TiO 2 interface. The PS NAC exhibits 1 meV errors, and the CA method reduces the error by an order of magnitude. We also demonstrate a modification of the CA method, in which the atomic displacement is monitored, the CA is used for atoms with smaller displacements, and the full AE calculation is performed for atoms with larger displacements. We show that the CA NAC accuracy is robust with respect to variations of the MD time step, such that the time step is limited by other considerations.

Next, we briefly review NAC in NAMD. For simplicity we do not use many-body states here, since NACs between manybody states can be reduced to NACs between KS orbitals. 9 The many-body states can be represented by Slater determinants of KS orbitals, or by their linear superpositions, as in linear response TDDFT and time-dependent GW-BSE theory. More details can be found in ref 9. TDSE describes evolution of the electronic subsystem in NAMD:

Implemented within the TDKS theory, the TD wavefunction Ψ(r,R,t) can be expanded using a set of instantaneous adiabatic KS orbitals ψ i (r;R(t)), which are obtained by solving the stationary KS equations at atomic configuration R(t):

With eq 2, eq 1 yields

Here,

is the NAC matrix element. NAC calculation plays a key role in NAMD. However, an analytic solution for NAC is not available in many cases. Even if the

term can be obtained analytically, its calculation may be expensive, as it is a highly multi-dimensional vector. Hence, the scalar NAC of eq 4 is routinely evaluated numerically using the Hammes-Schiffer-Tully formulation: 25

(5)

In eq 5, ψ j can be either the PS wavefunction or the AE wavefunction. Because PS wavefunctions are not orthonormal, a post-orthonormalization, e.g., Loẅdin transformation, is required, 26 or two additional terms should be introduced to correct the non-orthogonality, 22 as was done early on with both PAW and Vanderbilt pseudopotentials. 27,28 The AE wavefunction ψ k is constructed in the PAW method from the PS wavefunction ψ̃k using the following linear transformation:

The AE partial waves ϕ i a are typically chosen to be the solutions of the radial Schrodinger equation for atom a for a specific angular momentum, labeled by i. The index i is a shorthand for angular momentum numbers (l,m) and reference energy (ϵ). The details of index i depend on construction of a specific pseudopotential. The transformation from the PS wavefunction to the AE wavefunction is linear, and both PS and AE wavefunctions can be expanded around each atom site a into corresponding partial waves using a set of projector functions, |pĩ a ⟩. The sum

a k (7) represents the AE wavefunction inside the augmentation sphere of the site a, and the sum

represent the PS wavefunction inside the augmentation sphere of the site a.

The projector functions are dual to the partial waves:

Under the PAW formalism, for each atom site a, the AE partial waves |ϕ i a ⟩ and their PS counterparts are pairwise identical beyond the augmentation sphere with the cutoff radius r c :

The PS wavefunction ψ̃k can be expanded locally inside the augmentation sphere of atom a into the PS partial waves:

For any operator A, 29 eqs 6-8 lead to

The Journal of Physical Chemistry Letters pubs.acs.org/JPCL Letter ∑ ∑

If A is a local operator, then only term 13a needs to be explicitly evaluated. This is because, according to eq 10, χ j a -χj a in term 13b vanishes outside the augmentation sphere (the | r -R 0 | > r c region in Figure 1a, or the gray region in Figure 1b), while at the same time, according to eq 11, ψ̃j -χj a vanishes inside the augmentation sphere (the |r -R 0 | < r c region in Figure 1a, or the blue and orange region in Figure 1b). Therefore, these terms are never nonzero simultaneously in the same region in the space.

Similarly, in term 13c, there is no region for which χ j a -χj a and χ k a ′ -χk a ′ are both nonzero when a and a′ are different atom sites, unless the augmentation spheres (blue and orange regions in Figure 1b) of different sites overlap. However, such situations are excluded from the PAW approach, which fails in such a situation.

Thus, for a local operator, eq 13 can be reduced to

If an operator is nonlocal, such as the time-derivative operator

, the CA introduces an error, which is small if the displacement of nuclei within time dt along a trajectory are small compared with the radius of the augmentation sphere. In such case the bra and ket wavefunctions in term 13b still vanish in the majority region (|r -R 0 | < r c and |r -R 0 ′| > r c regions in Figure 1c). Only in the r c -|R 0 -R 0 ′| < |r -R 0 | < r c region, the bra and ket wavefunctions are nonzero simultaneously. Therefore, the term 13b is negligible. Term 13c is also small, unless the augmentation sphere at t + dt overlaps with the augmentation sphere of a different site at t, Figure 1d. However, this situation does not happen if dt is small. Note that the augmentation spheres can overlap if the material is studied at a very high pressure. In such cases the PAW method breaks down in general.

In the implementation, ψ̃j and ψ̃k are expanded in a planewave basis set, while ϕ i a , ϕ ̃i a , and p̃i a are atom-centered functions defined for each atom. The exact evaluation of the AE NAC with eq 13 needs integration between different grids, which is computationally complex and expensive. The computation cost of the exact AE NAC can be a substantial limitation during NAMD, which can require evaluation of millions of NAC matrix elements. Note that NAC evaluation requires more accurate convergence of the electronic wavefunction than force evaluation during adiabatic MD, further adding to the computational cost.

The numerical evaluation of NAC is performed over a small atomic time step, on the order of 1 fs, during which atoms move much less than the size of the augmentation sphere (core radius). A small time step is needed for an accurate evaluation of NAC with any electronic structure method, regardless of the PAW approach. Similarly, when the nuclear derivative vector, ⟨ψ i |∇ R |ψ j ⟩ in eq 4 is needed explicitly, it can be computed numerically over atomic displacements dR that are much smaller than the core radius. Therefore, it is reasonable to treat b), the nucleus is stationary. The blue and red solid lines represent the pseudo (PS) wavefunctions and all-electron (AE) partial waves, respectively. The corresponding dashed lines represent the PS partial waves. R 0 is the atomic center for a specific atom site a. r c is the augmentation sphere cutoff for this atom site. In (c) and (d), the nucleus undergoes a small displacement after dt. The corresponding PS wavefunction and partial waves are indicated for times t and t + dt, respectively. R 0 and R 0 ′ are the atomic centers for times t and t + dt. In (b) and (d), the circles indicate different atomic sites and augmentation spheres. The AE and PS wavefunctions are identical outside the augmentation sphere. Inside the sphere, the AE and PS wavefunctions are expanded in atom-centered partial waves. Since atomic displacement is negligible for a small dt, we consider the two atomic-centered AE wavefunctions concentric inside the core region, to achieve easier and efficient integration. Importantly, these wavefunctions are only used to correct the compensation part within the augmentation sphere (r c ), and the main contribution to NAC arising from the PS wavefunction is unchanged.

The Journal of Physical Chemistry Letters the center of the atomic partial waves (χ j a and χj a ) as stationary over the time step dt.

With the CA approximation, the exact overlap encountered in eq 5 is simplified to

in which we defined

The NAC, which is routinely evaluated numerically using the Hammes-Schiffer-Tully formulation, is reduced to

Approximating eq 13 with eq 15 constitutes the CA. Eq 15 has a form similar to the overlap of the eigen-wavefunctions derived in the original work by Blochl. 17 However, ψ j and ψ k of eq 15 are not from the same structure, nor they are orthonormal, and eq 16 changes correspondingly. The PS wavefunction contribution to the overlap, ⟨ψ̃j|ψ̃k⟩, is not modified, and the impact of atomic movements on the electronic properties in the interstitial region is adequately accounted for.

In the implementation of the CA method, we assume that the center of the radial functions (ϕ i a and ϕ ̃i a ) has not moved during dt (Figure 1b). Therefore, the integration of ⟨ϕ i a |ϕ i′ b ⟩ -⟨ϕ ̃i a |ϕ ̃i′ b ⟩ comes to a constant for a specific atom species and can be precalculated. The coefficients ⟨pĩ a |ψ̃k⟩ are already calculated within every ab initio calculation employing PAW, and therefore, they are readily available without further cost. Thus, our method does not introduce additional computation cost relative to the PS NAC calculation. In the current implementation, the CA NAC computation cost is about 3-4 orders of magnitude lower than calculation of the exact AE NAC and is determined by the disk I/O operations.

In addition to the accurate and fast evaluation of the NAC itself, the accuracy improvement provided by the CA NAC over the PS NAC is important for other advanced techniques developed for NAMD. For example, adiabatic KS orbitals are defined up to a phase, and phase-inconsistency of the adiabatic orbitals computed along MD trajectories results in an erroneous evaluation of the NAC, subsequently causing problems in the numerical integration of the NAMD equations and making NAMD results not reproducible. 30 Further, the trivial 31 or "unavoided" 32 crossing problem arises when adiabatic wavefunctions are localized in different parts of a system and their corresponding eigen-energies become quasidegenerate, e.g., with energy difference less than 5 meV. Both NAC phase-consistency and proper ordering of quasidegenerate states are required for correct NAC evaluation using eq 5 and rely on accurate calculation of the overlap ⟨ψ j |ψ k ⟩. 7,[30][31][32][33] The error introduced by the PS wavefunction can potentially amplify these problems, and the CA method allows for a more robust phase-consistency and state-tracking.

The developed approach is general and can be implemented with any PAW software package. Here, we use the Vienna Ab Initio Simulation Package. [34][35][36][37] More simulation details are provided in the Supporting Information (SI). To evaluate the performance of the CA method, we investigate three prototypical systems containing transition metals, namely, bulk Cu and rutile TiO 2 , and a CH 3 OH/TiO 2 interface (see SI). The optimized atomic structures and their corresponding densities of states (DOS) at 0 K are shown in Figure 2 and SI, Figure S1. In the previous work, 22 we found larger PS NAC errors in these types of systems containing transition metals. The deviations between the PS and AE NACs are notable for electronic transitions associated with sub-shell d electrons and can reach 2 meV. 22 The error is smaller if the system does not include a transition metal. Therefore, we illustrate the method with systems containing transition metals. At the same time, the efficient and accurate evaluation of NAC is important for other systems as well, in particular in situations when one encounters trivial crossings, and when NAC acquires a random sign/phase that has to be kept track of along the trajectory, as discussed above.

In the following, we calculate three types of NAC, i.e., the PS NAC, the CA NAC, and the exact AE NAC. The overlap between the exact AE wavefunctions employed for evaluation of the AE NAC is evaluated with PAWPYSEED. 38,39 We consider NAC between all pairs of orbitals in the range from VBM -10 to VBM + 10 in both TiO 2 and Cu at the Γ-point of the simulation supercells (Figure 3). VBM stands for valence band maximum, coinciding with the Fermi level. Similarly to what we found in the previous work, 22 PS NACs show notable error for transitions involving sub-shell d electrons. These errors can exceed 1 meV for both Cu and TiO 2 , raising a concern for employing PS NACs for systems with transition metals. NACs evaluated with the CA method show significant improvements. The deviations from the exact NAC are negligible, regardless whether or not sub-shell d electrons are involved.

In order to get a more quantitative picture on the improvement of the CA method over the PS approximation, we present average absolute and relative deviation of the NACs obtained with the two methods (Figure 4). The average NAC deviations obtained with the CA method remains within 0.05-0.1 meV for Cu and 0.02-0.04 meV for TiO 2 , while the corresponding deviations for the PS NACs reach 0.3-0.7 and 0.6-1.0 meV, respectively. Because the absolute NAC errors are generally independent of the NAC value itself, the relative error can be large for very small NAC. However, this is not a matter of concern, because transitions between states coupled by small NACs do not contribute to the overall dynamics, which is governed by pathways with large NACs. Typical NACs for interband and intraband relaxation range within 1-5 and 10-100 meV, respectively. [40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55] As demonstrated above, the CA method shows negligible difference from the exact AE NAC calculation, because atomic movements within dt are sufficiently small compared to the size of the augmentation sphere. We anticipate situations in which some atomic motions may be large, for example, with longer time steps, unstable systems undergoing large-scale motions, or light and fast hydrogen atoms. Note that in many cases hydrogen orbitals do not contribute to investigated processes. If some atoms move over relatively large distances, it may be desirable to perform the exact AE evaluation for these atoms and apply the CA for the other atoms. While this idea requires further testing with appropriately chosen systems for which the CA method leads to significant errors, we state it here and demonstrate it with the current systems (Table 1). We also increased dt to 10 fs, which is too large to be used in practical NAMD simulations, and found that, remarkably, the CA method still proved accurate (Figure 5). The average atom displacements in the current systems are 0.005 Å when a 1 fs time step is used. The largest displacement within the 1 fs time step is 0.012 Å. The results reported in the "partial CA" line in Table 1 were obtained by using the exact AE evaluation for atoms moving more than 0.005 Å within the 1 fs time step and the approximate CA evaluation for atoms moving less than 0.005 Å. The data obtained with the partial and full CA methods are comparable. The computational expense of the partial CA method scales linearly with the number of atoms requiring the AE treatment.

Finally, we examine the dependence of the deviation on the time step used in the discretization for evaluating NAC (Figure 5). All MD trajectories were obtained with the 1 fs time step, and we only varied the dt used in eq 17 for NAC, from the same trajectory data; e.g., we used every 10th point from the MD trajectories to obtain NACs for dt = 10 fs. Surprisingly, the deviations barely change when dt increases from 1 to 5 fs, although the mean ion displacement increases more than 5fold. A notably increased deviation is observed with the 10 fs time step for NACs larger than 10 meV. However, the deviation remains within 0.1 meV, and because the absolute NAC values are larger than 10 meV, the relative error remains under 1%. Therefore, we conclude that the CA method is robust to time step variations and that the choice of the MD time step is determined by other factors.

The NAC is the scalar product of the NAC vector (NACV),

, with the classical nuclear velocity, eq 4.

Some NAMD techniques, such as the Ehrenfest approach and the velocity rescaling step in surface hopping, require knowledge of NACV. 56,57 Although the CA can be used for numerical evaluation of NACV using finite displacements of atoms, an analytic evaluation of NACV is more efficient. Since NACV is a multidimensional vector, including hundreds of components, while NAC is a scalar, NAC calculation is preferable if possible. In surface hopping, one can compute that the scalar NAC is the vast overwhelming majority of the cases and obtain NACV only to rescale velocity. In such cases, the computational efficiency can be greatly increased by computing NAC using the CA method. In many nanoscale systems, atomic motions are dominated by thermal fluctuations, and the amplitudes of thermal atomic fluctuations are larger than the differences in system geometries between ground and excited states. This situation allows for the classical path approximation (CPA), which provides huge computational savings. 9 Instead of having to calculate hundreds or thousands of excited-state trajectories, the CPA allows one to precompute a single ground-state trajectory. Under the CPA, one can perform ab initio NAMD for systems involving hundreds of atoms for many picoseconds. Calculation of NACVs is not required under the CPA, while efficient evaluation of NAC is very important and can be achieved using the CA method in the PAW approach. The CPA methodology developed in our group 9 has demonstrated great successes in applications to many important classes of systems and processes. The calculation of NAC constitutes the main additional cost on top of the standard DFT calculation, and the developed CA method provides an important computational advance.

To summarize, the CA method demonstrated here allows efficient and accurate evaluation of NAC with the PAW method that is widely used in modeling electronic structure of condensed matter, nanoscale, and molecular systems. The CA NAC shows an order of magnitude improvement over the widely used PS NAC, making the accuracy comparable to the exact AE NAC. At the same time, the computational cost of the CA NAC evaluation is similar to that of the PS NAC and is 3-4 orders of magnitude lower than the evaluation of the AE NAC. While the PS NAC ignores the core region, the CA NAC explicitly includes evaluation of the wavefunction overlap in the core. The computational simplification is achieved by noting that atomic displacements over typical time steps used in MD simulations are much smaller than the core region size. Therefore, one can assume that the properties of the core region change little over the MD time step. As a result, the computationally expensive terms in the AE NAC can be either dropped or precomputed and used for the whole MD simulation. Additionally, we put forth an intermediate approximation, in which NAC contributions from atoms that move significantly are treated at the exact AE level, while the CA is applied to the remaining atoms. We demonstrate the accuracy of the CA method with two typical condensed matter systems containing sub-shell d electron levels that contribute to NA transitions. As demonstrated previously, such systems exhibit the largest NAC errors. The CA method allows consistent evaluation of NAC with 0.1 meV accuracy that is robust with changing MD time step. Remarkably, even for TiO 2 containing light O atoms, the CA NAC reproduces the The Journal of Physical Chemistry Letters AE NAC with better than 2% accuracy even for a 10 fs time step. Besides the NAC value itself, the accurate NAC evaluation can be important for NAMD algorithms that address the phase-consistency and trivial crossing issues in the NAC computation. The developed approach constitutes an important advance for successful application of NAMD to modern problems in physical chemistry and materials science.