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

###### Abstract

This article proposes a so-called Quasi Time-Reversible (QTR G-Ext) scheme based on Grassmann extrapolation of density matrices for an accurate calculation of initial guesses in Born-Oppenheimer Molecular Dynamics simulations. The method shows excellent results on four large molecular systems, ranging from 21 to 94 atoms simulated with Kohn-Sham density functional theory surrounded with a classical environment with 6k to 16k atoms. Namely, it clearly reduces the number of self-consistent field iterations, while keeping a similar energy drift as in the extended Lagrangian Born-Oppenheimer method.

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 [29] (SCF). BOMD simulations, that require one to perform tens of thousands of SCF calculations, rely thus heavily on extrapolation techniques [7, 1, 2, 27, 11, 12, 31, 32, 20, 21, 19, 34] that, by using converged solutions from previous iterations, 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 (BOMD) [20].
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 [21, 19].
In this particular case, the time-reversible extrapolation is augmented by the inclusion of a dissipative term, which serves to reduce numerical fluctuations.
In the last few years, the XLBO method has also been proposed in a SCF-less formulation, [18, 22, 34] where the density computed using the XLBO procedure is used directly without further SCF iterations.

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. However, the guess density obtained with XLBO is not exactly idempotent [21], which is in practice a problem that can be either ignored, or easily addressed by, for instance, McWeeny purification [17, 25].

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) [25, 24]. G-Ext has been successfully adopted in the Pisa-group for both ground- and excited-state SCF-based BOMD simulations in a polarizable multiscale framework [16]. While G-Ext is very effective, indeed outperforming XLBO in terms of the average number of SCF iterations required to achieve convergence along a MD trajectory [25], numerical experiments have shown that the extrapolation introduces a bias causing a drift in the total energy for NVE simulations. Such an energy drift is modest (few kcal/mol in 10 ps), but non negligible, at least for not-too-tightly-converged calculations [25]. In this contribution, we address such a limitation by introducing a new strategy to perform the extrapolation, that we name 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 and speed in extrapolating the density matrix during a BOMD simulation, as well as excellent energy conservation properties.

Given a $\mathcal{N}$-dimensional basis, the SCF solves the following nonlinear eigenvalue problem which consists to find a matrix $C$ and a diagonal matrix $E$ such that

$$\{\begin{array}{cc}F(D)C=SCE\hfill & \\ {C}^{T}SC={I}_{N}\hfill & \\ D=C{C}^{T},\hfill & \end{array}$$ |

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 $\stackrel{~}{C}={S}^{1/2}C$. Then the normalized density matrix $\stackrel{~}{D}=\stackrel{~}{C}{\stackrel{~}{C}}^{T}={S}^{1/2}D{S}^{1/2}$ belongs to the manifold

$$\mathcal{G}r(N,\mathcal{N})=\left\{D\in {\mathbb{R}}^{\mathcal{N}\times \mathcal{N}}\right|{D}^{2}=D={D}^{T},\mathrm{Tr}(D)=N\},$$ |

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 $\mathrm{\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

${\mathbb{R}}^{3M}$ | $\u27f6\mathcal{D}\u27f6{\mathcal{T}}_{{D}_{0}}\u27f6\mathcal{G}r$ | (1) | ||

$\bm{R}$ | $\u27fcd\u27fc\mathrm{\Gamma}\phantom{\rule{1em}{0ex}}\u27fcD,$ |

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 [30] $d\in {\mathbb{R}}^{{N}_{\mathrm{QM}}\times {N}_{\mathrm{QM}}}$,

$${(d)}_{ij}=\{\begin{array}{cc}0.5{z}_{i}^{2.4}\hfill & i=j,\hfill \\ \frac{{z}_{i}{z}_{j}}{\Vert \bm{R}({t}_{i})-\bm{R}({t}_{j})\Vert}\hfill & i\ne j,\hfill \end{array}$$ | (2) |

where ${N}_{\mathrm{QM}}$ is the number of atoms treated quantum mechanically and ${z}_{i}$ denotes the nuclear charge of the $i$th atom. Note that other descriptors can also be considered. We will detail the crucial mapping $d\mapsto \mathrm{\Gamma}$ below. The mapping $\mathrm{\Gamma}\mapsto \mathrm{Exp}(\mathrm{\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 \mathrm{Log}(D)=\mathrm{\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 [24, 35, 6]. 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 ${\mathrm{\Gamma}}_{n-i}=\mathrm{Log}({D}_{n-i})$, for $i=1,\mathrm{\dots},q$, the approximation of the density matrix on the tangent space is written as

$${\stackrel{~}{\mathrm{\Gamma}}}_{n}=-{\mathrm{\Gamma}}_{n-q}+\sum _{i=1}^{\stackrel{~}{q}}{\alpha}_{i}\left({\mathrm{\Gamma}}_{n-i}+{\mathrm{\Gamma}}_{n-q+i}\right),$$ | (3) |

where $\stackrel{~}{q}=q/2$ if $q$ is even, while $\stackrel{~}{q}=(q-1)/2$ if $q$ is odd. We remark that if in Eq. (3), the term ${\mathrm{\Gamma}}_{n-q}$ is substituted by ${\stackrel{~}{\mathrm{\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.

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

$$\underset{\bm{\alpha}\in {\mathbb{R}}^{\stackrel{~}{q}}}{\mathrm{min}}\left\{{\Vert {d}_{n}+{d}_{n-q}-\sum _{i=1}^{\stackrel{~}{q}}{\alpha}_{i}\left({d}_{n-i}+{d}_{n-q+i}\right)\Vert}^{2}+{\epsilon}^{2}{\Vert \bm{\alpha}\Vert}^{2}\right\},$$ |

where $\parallel \cdot \parallel $ denotes the ${\mathrm{\ell}}^{2}$-norm and $\epsilon >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

$$\underset{\bm{\alpha}\in {\mathbb{R}}^{\stackrel{~}{q}}}{\mathrm{min}}{\Vert \left[\begin{array}{c}\widehat{d}\\ \mathrm{\U0001d7ce}\end{array}\right]-\left[\begin{array}{c}\widehat{D}\\ \epsilon {I}_{\stackrel{~}{q}}\end{array}\right]\bm{\alpha}\Vert}^{2},$$ |

where the vector $\widehat{d}={d}_{n}+{d}_{n-q}$ is padded with $\stackrel{~}{q}$ zeroes, $\widehat{D}\in {\mathbb{R}}^{{N}_{d}\times \stackrel{~}{q}}$ is the matrix which columns are defined as ${\widehat{D}}_{\cdot ,i}={d}_{n-i}+{d}_{n-q+i}$, and ${I}_{\stackrel{~}{q}}$ is the identity matrix of order $\stackrel{~}{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 \mathrm{\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

${\mathrm{\Gamma}}_{n}$ | $=f({d}_{n})\approx f\left(-{d}_{n-q}+{\displaystyle \sum _{i=1}^{\stackrel{~}{q}}}{\alpha}_{i}\left({d}_{n-i}+{d}_{n-q+i}\right)\right)$ | ||

$=-f\left({d}_{n-q}\right)+{\displaystyle \sum _{i=1}^{\stackrel{~}{q}}}{\alpha}_{i}\left[f\left({d}_{n-i}\right)+f\left({d}_{n-q+i}\right)\right]$ | |||

$=-{\mathrm{\Gamma}}_{n-q}+{\displaystyle \sum _{i=1}^{\stackrel{~}{q}}}{\alpha}_{i}\left({\mathrm{\Gamma}}_{n-i}+{\mathrm{\Gamma}}_{n-q+i}\right)={\stackrel{~}{\mathrm{\Gamma}}}_{n}.$ |

The number $q$ of density matrices taken at previous steps and the value of the regularization parameter $\epsilon $ are chosen in a heuristic manner: we computed the error $\Vert {\mathrm{\Gamma}}_{n}-{\stackrel{~}{\mathrm{\Gamma}}}_{n}\Vert $ for different values of $q$ and $\epsilon $, specifically $q=3,4,\mathrm{\dots},20$ and $\epsilon =0.001,0.002,0.005,0.01,0.02,0.05$, and we selected the combination $(q,\epsilon )$ corresponding to the minimal error. When the SCF convergence threshold is ${10}^{-5}$, we found that good values are $q=5$ and $\epsilon =0.005$, while if it is fixed to ${10}^{-7}$, we found $q=4$ and $\epsilon =0.001,0.002$. Additional details on the selection of $q$ and $\epsilon $ values can be found in Section S1 of the SI.

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, 9]. Some information on the systems is reported in Table 1.

System | ${N}_{\text{QM}}$ | ${N}_{\text{MM}}$ | $\mathcal{N}$ | $N$ | ${N}_{d}$ |
---|---|---|---|---|---|

DMABN | 21 | 6843 | 185 | 39 | 234 |

3HF | 28 | 15046 | 290 | 62 | 409 |

AppA | 31 | 16449 | 309 | 67 | 468 |

OCP | 94 | 6058 | 734 | 154 | 4468 |

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 [10]. This is coupled with a polarizable description of the environment, using the AMOEBA forcefield [26]. 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 [23, 14, 15, 8]. We implemented the QTR G-Ext extrapolation approach in Tinker [28, 13].

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 [33]. 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 [25]

$${\stackrel{~}{\mathrm{\Gamma}}}_{n}=\sum _{i=1}^{q}{\alpha}_{i}{\mathrm{\Gamma}}_{n-i},q=6,$$ |

where the ${\alpha}_{i}$ are computed by solving

$$\underset{\bm{\alpha}\in {\mathbb{R}}^{q}}{\mathrm{min}}\left\{{\Vert {d}_{n}-\sum _{i=1}^{q}{\alpha}_{i}{d}_{n-i}\Vert}^{2}+{\epsilon}^{2}{\Vert \bm{\alpha}\Vert}^{2}\right\},\epsilon =0.01,$$ |

and XLBO without McWeeny purification [21, 19]

$${\stackrel{~}{D}}_{n}=2{\stackrel{~}{D}}_{n-1}-{\stackrel{~}{D}}_{n-2}+\kappa \left({D}_{n-1}-{\stackrel{~}{D}}_{n-1}\right)+c\sum _{i=1}^{8}{\alpha}_{i}{\stackrel{~}{D}}_{n-i},$$ |

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.

DMABN | 3HF | AppA | OCP | |||||
---|---|---|---|---|---|---|---|---|

STF | LTD | STF | LTD | STF | LTD | STF | LTD | |

QTR G-Ext | 0.33 | -0.01 | 0.62 | -0.40 | 0.57 | -0.08 | 0.36 | -0.23 |

G-Ext | 0.35 | -0.43 | 0.61 | -0.94 | 0.56 | -0.93 | 0.38 | -1.38 |

XLBO | 0.32 | 0.01 | 0.57 | -0.42 | 0.59 | 0.14 | 0.39 | -0.28 |

DMABN | 3HF | AppA | OCP | |||||
---|---|---|---|---|---|---|---|---|

STF | LTD | STF | LTD | STF | LTD | STF | LTD | |

QTR G-Ext | 0.37 | 0.01 | 0.59 | -0.30 | 0.53 | 0.18 | 0.38 | -0.16 |

G-Ext | 0.33 | 0.04 | 0.60 | -0.27 | 0.54 | 0.06 | 0.38 | -0.20 |

XLBO | 0.32 | 0.13 | 0.64 | -0.37 | 0.56 | 0.06 | 0.38 | -0.08 |

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 S3 of the SI.

DMABN | 3HF | AppA | OCP | |||||
---|---|---|---|---|---|---|---|---|

$\overline{k}$ | $\sigma $ | $\overline{k}$ | $\sigma $ | $\overline{k}$ | $\sigma $ | $\overline{k}$ | $\sigma $ | |

QTR G-Ext | 3.04 | 0.22 | 2.98 | 0.21 | 3.00 | 0.02 | 2.96 | 0.31 |

G-Ext | 3.55 | 0.85 | 3.16 | 0.69 | 3.03 | 0.54 | 2.91 | 0.41 |

XLBO | 4.00 | 0.05 | 4.00 | 0.00 | 4.00 | 0.07 | 4.00 | 0.01 |

DMABN | 3HF | AppA | OCP | |||||
---|---|---|---|---|---|---|---|---|

$\overline{k}$ | $\sigma $ | $\overline{k}$ | $\sigma $ | $\overline{k}$ | $\sigma $ | $\overline{k}$ | $\sigma $ | |

QTR G-Ext | 5.42 | 0.69 | 5.42 | 0.80 | 5.37 | 0.84 | 4.86 | 0.83 |

G-Ext | 7.33 | 0.63 | 6.96 | 0.79 | 6.56 | 0.75 | 5.83 | 0.87 |

XLBO | 7.51 | 0.65 | 7.45 | 0.65 | 7.43 | 0.80 | 7.21 | 0.75 |

In conclusion, we presented the novel Quasi Time-Reversible Grassmann extrapolation scheme aimed at preserving the energy conservation of Newton’s equations and, at the same time, at keeping low the number of SCF iterations. This scheme is based on the same properties of differential geometry of our previous extrapolation approach, ensuring that our guess density matrices retain all the mathematical properties of a density matrix. The innovation of this contribution lies in the symmetric combination of vectors in the tangent space, which proved to effectively preserve the stability of the total energy during the simulation. To validate its effectiveness, we conducted tests on systems of different sizes, and we obtained excellent results for all of them.

## Acknowledgments

This work was supported by the Italian Ministry of University and Research under grant 2020HTSXMA_002 (PSI-MOVIE) and by the French ‘Investissements d’Avenir’ program, project Agence Nationale de la Recherche (ISITE-BFC) (contract ANR-15-IDEX-0003). ÉP also acknowledges support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 810367–project EMC2) as well as from the Simons Targeted Grant Award No. 896630. BS acknowledges funding by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC 2075 – 390740016. FP is member of the GNCS group of INdAM.

## References

- [1] (1999) Ab initio molecular dynamics, a simple algorithm for charge extrapolation. Comput. Phys. Commun. 118 (1), pp. 31–33. External Links: ISSN 0010-4655, Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [2] (1992) Ab initio molecular-dynamics techniques extended to large-length-scale systems. Phys. Rev. B 45 (4), pp. 1538–1549. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [3] (1993) Density-functional thermochemistry. 3. the role of exact exchange. J. Chem. Phys. 98 (7), pp. 5648–5652. Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [4] (2020) Molecular mechanisms of activation in the orange carotenoid protein revealed by molecular dynamics. J. Am. Chem. Soc. 142 (52), pp. 21829–21841. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [5] (2020) The Multiple Roles of the Protein in the Photoactivation of Orange Carotenoid Protein. Chem 6 (1), pp. 187–203. External Links: Document, ISSN 24519294 Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [6] (1998) The Geometry of Algorithms with Orthogonality Constraints. SIAM J. Matrix Anal. Appl. 20 (2), pp. 303–353. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [7] (2016) On the existence of the optimal order for wavefunction extrapolation in Born-Oppenheimer molecular dynamics. J. Chem. Phys. 144 (24), pp. 244103. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [8] (2020) Gaussian Development Version, Revision J.16. Note: Gaussian, Inc., Wallingford CT, 2020. Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [9] (2021) From crystallographic data to the solution structure of photoreceptors: the case of the AppA BLUF domain. Chem. Sci. 12, pp. 13331–13342. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [10] (1972) Self-Consistent Molecular Orbital Methods. XII. Further Extensions of Gaussian-Type Basis Sets for Use in Molecular Orbital Studies of Organic Molecules.. J. Chem. Phys. 56 (5), pp. 2257–2261. Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [11] (2005) Accelerated, energy-conserving Born–Oppenheimer molecular dynamics via Fock matrix extrapolation. Phys. Chem. Chem. Phys. 7 (18), pp. 3269–3275. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [12] (1994) Exponential transformation of molecular orbitals. J. Chem. Phys. 101 (5), pp. 3862–3865. External Links: ISSN 0021-9606, Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [13] (2018) Tinker-HP: a massively parallel molecular dynamics package for multiscale simulations of large complex systems with advanced point dipole polarizable force fields. Chem. Sci. 9, pp. 956–972. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [14] (2017) Hybrid QM/MM molecular dynamics with AMOEBA polarizable embedding. J. Chem. Theory Comput. 13 (9), pp. 4025–4033. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [15] (2019) Towards large scale hybrid QM/MM dynamics of complex systems with advanced point dipole polarizable embeddings. Chem. Sci. 10 (30), pp. 7200–7211. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [16] (2023) Fast method for excited-state dynamics in complex systems and its application to the photoactivation of a blue light using flavin photoreceptor. J. Phys. Chem. Lett. 14 (5), pp. 1222–1229. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [17] (1960-04) Some recent advances in density matrix theory. Rev. Mod. Phys. 32, pp. 335–369. External Links: Document, Link Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [18] (2017) Next generation extended lagrangian first principles molecular dynamics. J. Chem. Phys. 147 (5), pp. 054103. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [19] (2009) Extended Lagrangian Born–Oppenheimer molecular dynamics with dissipation. J. Chem. Phys. 130 (21), pp. 214109. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics, A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [20] (2006-09) Time-reversible Born-Oppenheimer molecular dynamics. Phys. Rev. Lett. 97, pp. 123001. External Links: Document, Link Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [21] (2008) Extended Born-Oppenheimer Molecular Dynamics. Phys. Rev. Lett. 100 (12), pp. 123004. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics, A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics, A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [22] (2020) Density-matrix based extended lagrangian born–oppenheimer molecular dynamics. J. Chem. Theory Comput. 16 (6), pp. 3628–3640. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [23] (2020) Excited state Born-Oppenheimer molecular dynamics through a coupling between time dependent DFT and AMOEBA. Phys. Chem. Chem. Phys. 22, pp. 19532–19541. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [24] (2020) An approximation strategy to compute accurate initial density matrices for repeated self-consistent field calculations at different geometries. Mol. Phys. 118 (19-20), pp. e1779834. External Links: ISSN 0026-8976, Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics, A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [25] (2021) Grassmann extrapolation of density matrices for Born–Oppenheimer molecular dynamics. J. Chem. Theory Comput. 17 (11), pp. 6965–6973. Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics, A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics, A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [26] (2010) Current status of the AMOEBA polarizable force field. J. Phys. Chem. B 114 (8), pp. 2549–2564. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [27] (2004) Fock matrix dynamics. Chem. Phys. Lett. 386 (4), pp. 272–278. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [28] (2018) Tinker 8: software tools for molecular design. J. Chem. Theory Comput. 14 (10), pp. 5273–5289. Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [29] (1951) New developments in molecular orbital theory. Rev. Mod. Phys. 23, pp. 69–89. Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [30] (2012) Fast and accurate modeling of molecular atomization energies with machine learning. Phys. Rev. Lett. 108, pp. 058301. Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [31] (2003) An efficient orbital transformation method for electronic structure calculations. J. Chem. Phys. 118 (10), pp. 4365–4369. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [32] (2005) Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 167 (2), pp. 103–128. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [33] (1967) Computer “experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys. Rev. 159, pp. 98–103. Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [34] (2017) Performance of extended Lagrangian schemes for molecular dynamics simulations with classical polarizable force fields and density functional theory. J. Chem. Phys. 146 (12), pp. 124115. External Links: Document Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.
- [35] (2019) Manifold interpolation and model reduction. Note: http://arxiv.org/abs/1902.06502 Cited by: A Quasi Time-Reversible scheme based on density matrix extrapolation on the Grassmann manifold for Born-Oppenheimer Molecular Dynamics.

## Determination of optimal $q$ and $\epsilon $ values

The parameters $q$ and $\epsilon $ of the QTR G-Ext method were found by computing the error $\Vert {\mathrm{\Gamma}}_{n}-{\stackrel{~}{\mathrm{\Gamma}}}_{n}\Vert $ and averaging it over the whole simulation for different values of $q$ and $\epsilon $, specifically $q=3,4,\mathrm{\dots},20$ and $\epsilon $ = 0.001, 0.002, 0.005, 0.01, 0.02, 0.05, and we selected the combination $(q,\epsilon )$ corresponding to the minimal error.

DMABN | 3HF | AppA | OCP | |||||

SCF tolerance | ${10}^{-5}$ | ${10}^{-7}$ | ${10}^{-5}$ | ${10}^{-7}$ | ${10}^{-5}$ | ${10}^{-7}$ | ${10}^{-5}$ | ${10}^{-7}$ |

$q$ | 5 | 4 | 5 | 4 | 5 | 4 | 5 | 4 |

$\epsilon $ | 0.005 | 0.001 | 0.005 | 0.001 | 0.005 | 0.002 | 0.005 | 0.002 |

## Supplementary figures: Energy stability and number of SCF iterations

In the following, we report the total energy profile and number of SCF iterations per step for all the simulations performed in this work using the XLBO, G-Ext, QTR, and TR schemes. The time-reversible scheme is obtained by computing the guess density as follows:

$${\stackrel{~}{\mathrm{\Gamma}}}_{n}=-{\stackrel{~}{\mathrm{\Gamma}}}_{n-q}+\sum _{i=1}^{\stackrel{~}{q}}{\alpha}_{i}\left({\mathrm{\Gamma}}_{n-i}+{\mathrm{\Gamma}}_{n-q+i}\right),$$ | (4) |

Eq. 4 is manifestly symmetric, and thus fully time reversible. As mentioned in the main text, the fully TR scheme exhibits excellent stability, but poor performance, as the number of SCF iterations tends to quickly increase along the simulation, to the point that the extrapolation is not anymore beneficial.