A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics

Ab-initio Born-Oppenheimer Molecular Dynamics (BOMD) is a very powerful and versatile tool to simulate molecular processes where the quantum nature of the system is not negligible. Unfortunately, this comes at a high computational price, which stems from the necessity of solving the quantum mechanical (QM) equations, typically Kohn-Sham Density Functional Theory (KS-DFT) equations, to compute the energy and forces at every time step. Such equations are nonlinear and are solved using a fixed-point iterative method known as Self-Consistent Field [33] (SCF). Achieving SCF convergence typically requires, in a standard single-point run, up to 20 iterations, making the MD simulation very expensive, as in turn the SCF has to be performed tens of thousands of times. Two main families of methods have been developed to address such a limitation. In extended Lagrangian methods, such as Car-Parrinello Molecular Dynamics (CPMD)[6] or Atom-centered Density matrix propagation (ADMP)[35], the electronic degrees of freedom are propagated, avoiding thus the need of solving the SCF problem. This requires one to endow the electronic degrees of freedom with a fictitious mass, that needs to be small enough to keep the trajectory close to its Born-Oppenheimer counterpart. As a consequence, rather short time steps need to be used. A different strategy relies on developing extrapolation techniques[8, 1, 2, 31, 12, 13, 36, 37, 22, 23, 21, 39] for BOMD that allow one to converge the SCF in a limited number of iterations. In this work, we choose the second strategy, which is particularly effective for calculations using localized basis sets, e.g., Gaussian-type orbitals. The extrapolation techniques used in BOMD use converged solutions from previous MD steps to compute an accurate guess for the SCF, limiting thus the number of iterations required to achieve convergence. A significant contribution to this field was given by Niklasson and co-workers in 2006, with their work on the time-reversible extrapolation for Born-Oppenheimer Molecular Dynamics [22]. The core concept involves generating a guess density matrix by combining the density matrices from previous steps in a symmetric and time-reversible manner. However, numerical applications showed that enforcing an exact time-reversibility can lead to errors accumulating in long-time simulations, spoiling thus the convergence properties of the algorithm in the long run. This led to the development of the Extended Lagrangian Born-Oppenheimer approach (XLBO) in 2008 [23, 21, 20, 24, 39]. In this particular case, the time-reversible extrapolation is augmented by the inclusion of a dissipative term, which serves to reduce numerical fluctuations. XLBO can be seen as an intermediate strategy between Car-Parrinello like approaches and extrapolation techniques for BOMD, as it indeed propagates an auxiliary density matrix that can either be used directly in a CPMD spirit[25, 14], possibly after refining the density using an approximate SCF solver, or used as a guess for the SCF[23]. Here, we focus on the latter approach.

In Niklasson’s XLBO scheme, the guess density is propagated in time subject to a potential that forces it to be close to the converged density. The result is a guess density that is accurate enough to achieve reasonable SCF convergence (e.g., $10^{-5}$ RMS norm of the density matrix change) in as little as four iterations: Niklasson’s pioneering work has therefore been crucial in extending the applicability of BOMD. However, the XLBO method suffers of a few shortcomings. First, the guess density obtained with XLBO is not exactly idempotent [23], unless it is postprocessed using, e.g., McWeeny purification[19, 29]. Second, its performance degrades if a tightly converged SCF solution is required, as it is the case when a post Hartree-Fock method is used to compute the energy and forces (e.g., in a time-dependent DFT excited state simulation).

Recently, we proposed a different strategy to compute a guess density using linear extrapolation. This is non-trivial, because in general a linear combination of density matrices does not preserve idempotency or, in other words, density matrices belong to a differentiable manifold called Grassmann manifold and not to a vector space. Our approach uses tools from differential geometry to map the Grassmann manifold onto its tangent space, which is a vector space. It then performs a linear extrapolation on the tangent space, and then maps back the extrapolated density to the manifold. We named such a method Grassmann extrapolation (G-Ext) [29, 28]. G-Ext is an accurate and efficient strategy for ab-initio MD simulations, that has been shown to outperform XLBO, especially if a tight SCF convergence is required[29]. It has been successfully adopted in the Pisa-group for both ground- and excited-state SCF-based BOMD simulations in a polarizable multiscale framework [4, 10, 18, 26]. Unfortunately, G-Ext suffers from a serious shortcoming. Numerical experiments have shown that the extrapolation introduces a bias causing a drift in the total energy for NVE simulations[29]. Such an energy drift is modest in absolute terms (few kcal/mol in 10 ps, to be compared with total energies of hundreds of thousands kcal/mol), but large if compared with the energetics of typical chemical processes. While using a tight convergence criterion for the SCF solves the problem[29], this is not an option for expensive, production simulations, limiting thus the gains introduced with the overall approach.

In this contribution, we not only address such a limitation by introducing a new strategy to perform the extrapolation, but further improve the performance of the method. We name the new strategy Quasi Time-Reversible Grassmann extrapolation method (QTR G-Ext). This approach leverages the principles of differential geometry, similarly to the previous method, but offers enhanced accuracy, improved performances, and excellent energy conservation properties. Given a $\mathcal{N}$-dimensional atomic orbitals (AO) basis, the SCF solves the following nonlinear eigenvalue problem which consists to find a matrix $C$ and a diagonal matrix $E$ such that

where $C\in\mathbb{R}^{\mathcal{N}\times N}$ contains the $\mathcal{N}$ coefficients of the $N$ occupied molecular orbitals, $D\in\mathbb{R}^{\mathcal{N}\times\mathcal{N}}$ is the density matrix, $E\in\mathbb{R}^{N\times N}$ is a diagonal matrix which entries are the energy levels, $F$ denotes the DFT operator, $S\in\mathbb{R}^{\mathcal{N}\times\mathcal{N}}$ is the overlap matrix, and $I_{N}$ denotes the identity matrix of order $N$.

We assume that the density matrix is orthogonal. In any case, it can be transformed into such matrix by considering the Löwdin factorization of the overlap matrix $S$ and consequently the modified coefficient matrix $\widetilde{C}=S^{1/2}C$. Then the normalized density matrix $\widetilde{D}=\widetilde{C}\widetilde{C}^{T}=S^{1/2}DS^{1/2}$ belongs to the manifold

which is isomorphic to the so-called “Grassmann manifold”, therefore we identify $\mathcal{G}r$ by this name. From now on, we assume that the density matrix has been orthonormalized and we denote it by $D$.

Since $\mathcal{G}r$ is a differential manifold, given a point $D_{0}\in\mathcal{G}r$, there exists a tangent space $\mathcal{T}_{D_{0}}\subset\mathbb{R}^{\mathcal{N}\times N}$, such that tangent vectors $\Gamma(D)\in\mathcal{T}_{D_{0}}$ can be associated to nearby points $D\in\mathcal{G}r$.

In MD, $t\to\bm{R}(t)$ represents the trajectory of the nuclei. The transformation of the electronic structure can be interpreted as a trajectory denoted by $t\to D_{\bm{R}(t)}$ on the manifold. In order not to burden the notation, we simply indicate $D$ in place of $D_{\bm{R}(t)}$. The objective is to determine a suitable approximation for the density matrix at the next step of the molecular dynamics trajectory by extrapolating the densities from previous steps. Since the tangent space $\mathcal{T}_{D_{0}}$ is a vector space, we approximate the density matrix on $\mathcal{T}_{D_{0}}$. In order to solve the extrapolation problem, we decompose the mapping $\bm{R}\to D$ as a composition of several maps

where the first function $\bm{R}\mapsto d$ is a map from atomic positions to molecular descriptors. Here, as a descriptor, we use the Coulomb matrix [34] $d\in\mathbb{R}^{N_{\rm QM}\times N_{\rm QM}}$,

where $N_{\rm QM}$ is the number of atoms treated quantum mechanically and $z_{k}$ and $\bm{R}_{k}$ denotes the nuclear charge and the position of the $k$th atom, respectively. Note that other descriptors can also be considered. We will detail the crucial mapping $d\mapsto\Gamma$ below. The mapping $\Gamma\mapsto\operatorname{Exp}(\Gamma)=D$ is the so-called Grassmann exponential which maps tangent vectors on $\mathcal{T}_{D_{0}}$ to $\mathcal{G}r$, and it is a locally bijective function in a neighborhood of $D_{0}$. Its inverse $D\mapsto\operatorname{Log}(D)=\Gamma(D)$ is the Grassmann logarithm. These mappings are computed by means of the singular value decomposition (SVD). For mathematical details, the interested reader is referred to [28, 40, 7]. In our method, during the MD, we use a fixed reference point $D_{0}$ to construct the tangent space $\mathcal{T}_{D_{0}}$.

Let $n$ be the current time step of the MD. Given previous $q$ snapshots $\Gamma_{n-i}=\operatorname{Log}(D_{n-i})$, for $i=1,\ldots,q$, the approximation of the density matrix representative on the tangent space is written as

where $\widetilde{q}=q/2$ if $q$ is even, while $\widetilde{q}=(q-1)/2$ if $q$ is odd. We remark that if in Eq. (3), the term $\Gamma_{n-q}$ is substituted by $\widetilde{\Gamma}_{n-q}$, a “fully” time-reversible approach (instead of quasi time-reversible) is obtained. Numerical experiments with the fully time-reversible approach, that are reported in the Supporting Information (SI), showed good behavior for total energy conservation, but unfortunately a strong increase in the number of performed SCF iterations. This is consistent with what has been observed by Niklasson and coworkers[21], who remark that exact time-reversibility under noisy conditions (e.g., not fully converged SCF) can lead to error accumulations and significantly worse SCF convergence.

The descriptors are involved in the computation of the coefficients $\bm{\alpha}=[\alpha_{1},\ldots,\alpha_{\widetilde{q}}]^{T}$ appearing in Eq. (3). Indeed, they are computed by solving the least-squares problem with Tikhonov regularization

where $\|\cdot\|$ denotes the $\ell^{2}$-norm and $\varepsilon>0$ is the regularization parameter. Since the Coulomb matrix (2) is symmetric, in the above formula $d_{j}$ represents the vectorized Coulomb matrix considering the lower triangle. In matrix form, it corresponds to solving the following least-squares problem

where the vector $\bm{b}=d_{n}+d_{n-q}$ is padded with $\widetilde{q}$ zeroes, $A\in\mathbb{R}^{N_{d}\times\widetilde{q}}$ is the matrix which columns are defined as $A_{\cdot,i}=d_{n-i}+d_{n-q+i}$, and $I_{\widetilde{q}}$ is the identity matrix of order $\widetilde{q}$. Then the initial guess for the density matrix is obtained as the composition of the three maps in (1), where the second map $d\mapsto\Gamma$ is given by (3). Note that if this second map denoted by $f$ was linear, then the guess would be close to exact, namely

After computing the coefficients $\alpha_{i}$ by solving (4) and the tangent vector $\widetilde{\Gamma}_{n}$ by Eq. (3), we obtain the sought guess density matrix for the SCF iterative method as $\widetilde{D}_{n}=\operatorname{Exp}(\widetilde{\Gamma}_{n})$.

The number $q$ of density matrices taken at previous steps and the value of the regularization parameter $\varepsilon$ are chosen in a heuristic manner: we computed the error $\|\Gamma_{n}-\widetilde{\Gamma}_{n}\|$ for different values of $q$ and $\varepsilon$, specifically $q=3,4,\ldots,20$ and $\varepsilon=0.001,0.002,0.005,0.01,0.02,0.05$, and we selected the combination $(q,\varepsilon)$ corresponding to the minimal error. When the SCF convergence threshold is $10^{-5}$, we found that good values are $q=5$ and $\varepsilon=0.005$, while if it is fixed to $10^{-7}$, we found $q=4$ and $\varepsilon=0.001,0.002$. Additional details on the selection of $q$ and $\varepsilon$ values can be found in Section S1 of the SI. The computational cost to compute the extrapolation coefficients $\bm{\alpha}$ is negligible compared to the time for a single MD step. Thanks to the symmetric property of the coefficients, the size of the system (5) is $(N_{d}+\tilde{q})\times\tilde{q}$, and $\tilde{q}$ is a small number (in our simulations $\tilde{q}=2$, as $q=4$ or $q=5$).

The QTR G-Ext approach is tested on four different systems. The first system is dimethylaminobenzonitrile (DMABN) in methanol. The second system is 3-hydroxyflavone (3HF) in acetonitrile. The last two systems (OCP and AppA) are chromophores embedded in a biological matrix-namely, a carotenoid in the orange carotenoid protein (OCP) and a flavin in the AppA Blue-Light Using Flavin photoreceptor [5, 4, 10]. Some information on the systems is reported in Table 1.

KS-DFT has been adopted to describe the QM subsystem, with the B3LYP hybrid functional [3] and the 6-31G(d) Pople’s basis set [11]. This is coupled with a polarizable description of the environment, using the AMOEBA forcefield [30]. For each system, we performed a QM/AMOEBA geometry optimization until a root-mean-square norm on the forces of 4 kcal/mol/Å is found and finally a 2 ps QM/AMOEBA NVT equilibration to obtain the starting point of the simulations presented in this work.

All simulations have been performed using the Gaussian-Tinker interface [27, 16, 17, 9]. We implemented the QTR G-Ext extrapolation approach in Tinker [32, 15].

To assess the quality of the guess density obtained by the QTR G-Ext extrapolation, we performed 10 ps BOMD simulations, with 0.5 fs time step, in the NVE ensemble, using the velocity Verlet integrator [38]. All systems were tested with an SCF convergence threshold fixed to $10^{-5}$ and $10^{-7}$ with respect to the RMS variation of density. We compare our approach in terms of energy stability and number of iterations required to reach convergence with other two extrapolation schemes, which are the G-Ext scheme [29]

where the $\alpha_{i}$ are computed by solving

where $\varepsilon=0.01$, and XLBO[23, 21]

with fixed parameters $\kappa=1.86$, $c=0.0016$, and $\bm{\alpha}=(-36,99,-88,11,32,-25,8,-1)$.

Figure 1 provides the plot of the total energy along the DMABN simulation, with a $10^{-5}$ SCF convergence threshold. Despite the non-fully time-reversible formulation of our newly implemented approach, we observe a great improvement with respect to the G-Ext scheme. In particular, the QTR G-Ext method resembles the fully time-reversible scheme XLBO. The same behaviour is almost imperceptible when the SCF convergence is set to $10^{-7}$ (Figure 2), since the accumulation of errors that generates the energy drift when G-Ext is used is lower, so we can appreciate the same trend with all the extrapolation schemes. Analogous figures are reported in Section S2 of the SI for all tested systems. To better evaluate the energy stability, we consider the average short-time fluctuation (STF) of the energy, which is computed by getting the RMS of the energy fluctuation every 50 fs and averaging over the trajectory, and the long-time drift (LTD) for a long-time analysis, that is the slope of the linear regression line of the energy. Tables 2 and 3 disclose STF and LTD for convergence thresholds $10^{-5}$ and $10^{-7}$, respectively. QTR G-Ext, G-Ext, and XLBO show comparable STF, which is specific for the system and is related to the time step for the integration. On the other hand, the absolute value of LTD is in general higher for $10^{-5}$ simulations, in particular for G-Ext. We can state that the QTR G-Ext method solves the energy-drift issue of G-Ext, showing an LTD that is always similar to the XLBO one, suggesting again a good time-reversible behaviour.

The gain of our new methodology is not only in terms of accuracy (energy stability), but also in terms of the computational time of the simulation. Tables 4 and 5 report the average number of SCF iterations required to achieve convergence, as well as the standard deviation for $10^{-5}$ and $10^{-7}$ SCF thresholds, respectively. We remark that each strategy requires $q$ previous density matrices, before having them available a standard SCF is performed. Therefore, for the computation of average and standard deviation, we discard the first $q$ points. The two tables show that for all the tested systems, the QTR G-Ext method requires the lowest number of SCF iterations, for both convergence thresholds. Moving averages of SCF iteration numbers during the simulations for all systems and with both SCF convergence thresholds are reported in Section S2 of the SI.

The performances of the QTR G-Ext guess are maintained also for larger and smaller time steps. We compared QTR G-Ext and XLBO for time steps of 0.1, 0.25, 0.75, and 1 fs by running MD simulations for the DMABN system with an SCF converge threshold of $10^{-5}$. All the results can be found in Section S3 of the SI. Both methods show excellent energy conservation for the smaller time steps, and afford reasonably stable simulations even for the larger ones, which is remarkable, as such simulations employ a time step that is too large to accurately sample molecular vibrations involving protons and are in general very noisy. For all time steps, QTR G-Ext requires a smaller average number of SCF iterations than XLBO. Finally, we tested the method for a looser SCF convergence of $10^{-4}$, again, a value that should not be used for production applications, as the error on the SCF solution transfers to the forces, affecting thus the quality of the dynamics. The results are reported in Section S4 of the SI. Again, good energy conservation is shown for both methods, and QTR G-Ext outperforms XLBO in terms of average number of SCF iterations required.

In conclusion, we presented the Quasi Time-Reversible Grassmann Extrapolation scheme, a new extrapolation method for ab-initio molecular dynamics that not only allows for energy-conserving simulations, but exhibits overall excellent performances. Our numerical tests, performed on large, complex systems described with a polarizable multiscale strategy and taken from real-life production applications, show that QTR G-Etx is able to provide a guess density to BOMD simulations that allows the convergence of the SCF procedure in about 3 iterations on average for convergence thresholds that are typical of ground-state production runs, which is a 25% gain with respect to the state-of-the-art XLBO method. Tighter convergences, that are required for, e.g., time-dependent DFT excited state simulations, can also be achieved in as little as 5-6 iterations. Furthermore, our numerical tests show that the new method does not introduce any significant bias in the guess density and thus exhibits very good energy conservation properties. This can be clearly seen by comparing the long-term drift observed in simulations for the two different SCF convergence thresholds used in our tests: while the previous G-Ext method shows a sharp increase in the drift going from $10^{-7}$ to $10^{-5}$ SCF convergence threshold, this is not the case for the QTR G-Ext method. We stress here that, due to the cost of BOMD simulations, every gain in performances is important, as it can easily translate in hundreds or thousands of saved CPU hours. The QTR G-Ext method is easy to implement and does not introduce any significant computational overhead, and represents therefore an effective strategy to extend the applicability of BOMD simulations to larger and more complex systems.