Simulations of molecular photodynamics in long timescales.

Nonadiabatic dynamics simulations in the long timescale (much longer than 10 ps) are the next challenge in computational photochemistry. This paper delimits the scope of what we expect from methods to run such simulations: they should work in full nuclear dimensionality, be general enough to tackle any type of molecule and not require unrealistic computational resources. We examine the main methodological challenges we should venture to advance the field, including the computational costs of the electronic structure calculations, stability of the integration methods, accuracy of the nonadiabatic dynamics algorithms and software optimization. Based on simulations designed to shed light on each of these issues, we show how machine learning may be a crucial element for long time-scale dynamics, either as a surrogate for electronic structure calculations or aiding the parameterization of model Hamiltonians. We show that conventional methods for integrating classical equations should be adequate to extended simulations up to 1 ns and that surface hopping agrees semiquantitatively with wave packet propagation in the weak-coupling regime. We also describe our optimization of the Newton-X program to reduce computational overheads in data processing and storage. This article is part of the theme issue 'Chemistry without the Born-Oppenheimer approximation'.

Introduction

Simulation of ultrafast molecular photodynamics is a mature field now [1-3]. Computational chemists have an arsenal of methods and research protocols to predict the molecular evolution within a few picoseconds. Nevertheless, most photoinduced processes occur in much longer timescales, several orders of magnitude over our current research capabilities. A few groups, including ours, are already probing dynamics simulations into these long timescales [4-13].

Although the field of molecular photodynamics in long timescales is still in the early stages, it should soon undergo an explosive rush of publications, primarily fed by the speed-up allowed by machine learning (ML) techniques [14,15]. However, the challenges are immense and go much beyond computational speeding up. From software optimization to methods' adequacy and algorithmic stability, there are issues in every corner. Thus, the maturation of this research area will require tremendous effort in the following years.

Hence, we decided to write a perspective article about photodynamics in long timescales to help frame the scope of the endeavour ahead. Our group was awarded an ERC Advanced Grant in 2019 precisely to work out methods for this field. Our research since then has provided us with a unique perspective on this subject, which we think is worth sharing in this special issue about Chemistry without the Born–Oppenheimer approximation.

Before advancing in the discussion, we would like to delimit the scope of this paper by clearing what we expect from a method to run nonadiabatic dynamics in the long timescale. We loosely define ‘long timescales' as periods much longer than 10 ps. Moreover, we believe that the main advantage of dynamics simulations over statistical approaches (like reaction rate calculations) is to propagate the molecular dynamics in full nuclear dimensionality, without bias towards any set of coordinates. Therefore, our first methodological constraint is to focus on mixed quantum classical (MQC) methods [1], which have full dimensionality as the main feature. (All acronyms used in this paper are defined in the electronic supplementary material, SM-1.) Second, we also primarily aim at general methods, which may be used for any molecular system (as opposed, for instance, to a particular potential energy surface fitted for a specific molecule). Finally, long time-scale simulations should not depend upon unrealistically large computational resources, especially if we consider the high economic and environmental costs of high-performance computing [16]. Thus, our third constraint is that those future methods developed to tackle such simulations should ideally allow us to run dynamics without significantly demanding new computational resources.

Having these problem delimitations in mind, we will survey the following topics impacting simulations of excited-state phenomena occurring over long timescales: (i) electronic structure methods speed-up; (ii) stability of integration methods; (iii) accuracy of nonadiabatic dynamics methods; and (iv) software optimization. All this discussion is supported by a series of simulations of adiabatic, surface hopping and quantum wave packet dynamics using a model Hamiltonian, as well as surface hopping and ML modelling of an actual molecule.

Electronic structure method speed-up

Long time-scale simulations will require a substantial reduction in the cost of quantum mechanical quantities (excited state energies, forces and state couplings) to keep computational resource requirements under control. We envisage three possible strategies to achieve this cost reduction: (i) use of parameterized and approximated electronic structure methods; (ii) use of parameterized full-dimensional Hamiltonian models; and (iii) use of ML as a surrogate model for quantum mechanical predictions. They are discussed in the following subsections.

Parameterized and approximated electronic structure methods

Parameterized electronic structure methods speed up quantum chemical (QC) calculations by replacing computationally intensive steps like integral calculations with precomputed quantities either based on empirical [17-19] or computational data [20,21]. There are several parameterized electronic structure methods currently available for MQC dynamics, allowing low-cost excited state calculations. Some of them are based on configuration interaction procedures, like FOMO-CI [18] and OM2/CI [17]. In turn, TD-DFTB [20] and the semiempirical TD-SCF [3] compute excited states based on linear-response time-dependent theory.

A different strategy to reduce computational costs without parameterization is met in the real-time time-dependent approach based on Kohn–Sham orbitals (TD-KS) [22], which has also been extensively explored for inexpensive MQC dynamics [23]. TD-KS uses KS orbital energy gaps for approximated excitation energies and estimates forces and nonadiabatic couplings from orbital derivatives.

For decades, parameterized electronic structure methods have been pushing the boundaries of computational photochemistry. The first on-the-fly surface hopping simulations were performed with the semiempirical MMVB method [24]. More recently, nonadiabatic dynamics simulations of the order of 10 ps have been reported with FOMO-CI [18], and nonadiabatic dynamics of nanoscopic systems with hundreds of atoms have been treated with real-time TD-KS based on DFTB orbitals [25].

Nevertheless, parameterized and approximated electronic structure methods usually have low accuracy when considering extended regions of the configurational space [26]. Moreover, for the former methods, the transferability of their parameters to simulate different systems is always a significant issue [27].

The final pressing issue with this strategy for cost reduction is the constrain we proposed in §1: long time-scale dynamics should not require significantly more computational resources than conventional ultrafast dynamics. Therefore, we need to compute excited state properties about 1000 times faster than today to satisfy such a condition. Ground state calculations with parameterized electronic structure methods are swift. DFTB, for instance, has been reported to be 1000 times faster than DFT with a hybrid functional [28]. However, the excited state calculations, based, for instance, on CI procedures, may still pose high costs.

Parameterized Hamiltonian models

Another strategy to bypass the high costs of solving the quantum electronic problem is constructing a model Hamiltonian that mimics the physical process of interest accurately and meaningfully. Here again, we reiterate that we wish to explore our system in its full nuclear dimensionality inside MQC. In this subsection, we limit our discussion to two types of model Hamiltonians, which can be used to study multi-state nonadiabatic processes by MQC dynamics in full dimensionality, the spin-boson Hamiltonian (SBH) and the vibronic coupling (VC) models.

SBH is a two-state model initially developed for describing dissipative systems [29,30]. However, due to its inherently generic feature and flexibility regarding the number of degrees of freedom and coupling terms, SBH became one of the favourite toy potentials to model nonadiabatic phenomena in complex macromolecular systems [31]. The SBH model describes a system with two diabatic electronic states that are linearly coupled to a harmonic bath consisting of N bosons with masses M, frequencies, position R and momentum P. Thus, the full Hamiltonian is where

In the isolated system, the energy bias and the electronic interstate coupling between the two diabatic electronic states are and , respectively. are the usual Pauli matrices and is the identity matrix of size 2. The coupling constants between the system and bath are determined by the spectral density, which characterizes the effect of the bath on transitions between the electronic states. Different types of spectral density have been designed, out of which Ohmic and Debye forms are typically used to study nonadiabatic dynamics. More complex and multi-peaked spectral density can also be constructed by fitting the appropriate experimental data [31,32]. After that, to work with this model, discretization of the coupling constants and frequencies can be carried out by standard protocols [33-35]. SBH can be easily transformed into adiabatic representation [36] to obtain the adiabatic energies, gradients and nonadiabatic coupling vectors (see the electronic supplementary material, SM-2).

Despite being extensively explored in MQC and quantum nonadiabatic dynamics [36-43], SBH remains much needed in the future as a model Hamiltonian to study dynamics in the long timescale. Examples are already discussed below in §3a, where we use SBH to propagate classical dynamics for 1 ns. In §4, we simulated quantum wave packet and surface hopping dynamics with SBH for up to 100 ps. Examples of adiabatic potential energy topographies that SBH can describe are shown in electronic supplementary material, SM-2.

One of the main limitations of the SBH model is its restriction to two electronic states. The VC model developed by Koppel et al. [44] (and, for this reason, also known as KDC Hamiltonian) can be an alternative to deal with an arbitrarily large number of states. For electronic states, the diabatic Hamiltonian is decomposed into the kinetic energy of the nuclei and harmonic ground state potential energy around some reference nuclear geometry (Q0), using mass-frequency scaled dimensionless normal-mode coordinates (Q) and a VC matrix (W), which is a low-order Taylor series expansion around the reference nuclear geometry. This prescription renders

The zeroth-order coupling matrix is diagonal and contains the adiabatic excitation energies at the reference point Q0. Truncating the series after the first order leads to the linear VC (LVC) model, which contains the intrastate coupling constant for the electronic state and mode in the diagonal elements of and the interstate coupling constant between states m and n in the off-diagonal elements. That means, and Beyond LVC, higher order VC models, including more terms in the Taylor series expansion (equation (2.3)), have also been constructed to describe the PESs [45,46] better and more accurately.

Usually, the interstate coupling parameters are obtained by fitting ab initio potential energy surfaces. However, parameters have also been successfully determined using excited state Hessian and wave function overlaps [47-49]. Recently, Plasser et al. [50] showed the workability of parameterizing LVC models using only a ground state normal modes calculation and a single-point excited state calculation (including energy gradients) at the equilibrium geometry. Moreover, they proposed a protocol to employ LVC with the standard MQC dynamics methods in adiabatic representation.

It is a common practice to construct the VC model Hamiltonians in reduced dimensionality considering only the important modes for the photophysical process of interest. A 24-dimension quadratic VC model for pyrazine published over 20 years ago still remains an almost solitary example of full-dimensional fitting of a non-trivial molecule [51]. In §1, we stated that we are primarily interested in methods general enough to be applied to any molecule in full nuclear dimensionality. However, the parameterization process is still the main challenge in using VC models for general, full-dimensional dynamics. We include them here because the advances in ML techniques for fitting diabatic surfaces may soon favourably change this scenario, allowing fast model parameterization [45,46].

ML for excited states

In many MQC dynamics applications, the QC properties of the molecular systems are usually calculated for each time step during the simulation. Such a repeatedly large number of calls to a QC program is the main bottleneck to run long timescale in nonadiabatic dynamics simulations. In this scenario, ML has emerged as a promising approach to speed up the prediction of molecular properties required for excited state molecular dynamics while still retaining the same level of accuracy as the underlying QC method used as a reference to train the ML predictor [14,15].

The main advantage of using ML models as a surrogate for the QC method comes from the fact that, once the ML model has been fitted to the reference data to reproduce QC properties within the desired accuracy, the solution of the Schrödinger equation is wholly skipped, thereby resulting in extremely fast predictions for new molecular configurations produced during the dynamics. This type of multi-dimensional regression problem is known as supervised learning in the ML field, where kernel methods and neural networks have become the preferred algorithms for fitting (highly nonlinear) chemical data [52].

However, performing fully ML-driven excited state dynamics is rather challenging [14,15]. One crucial issue is related to the intrinsic multi-dimensional nature of the target properties required for an MQC simulation. The ML model should accurately predict the potential energy of the diverse molecular configurations accessible in each electronic state (which are scalar properties) and interatomic forces and nonadiabatic coupling vectors (vectorial quantities). Also, the ML model should preserve possible correlations between the predicted variables (e.g. interatomic forces are proportional to the energy gradients), so that the ML-driven dynamics still obey certain physical principles such as energy conservation. Finally, the accuracy of the ML model is strictly related to the sampling strategy adopted to generate the training dataset, which should include the most representative molecular configurations accessible during the time evolution of the system on different excited-state surfaces.

Despite these challenges, the highly efficient computational performance of various ML models has been successfully demonstrated for some prototypical photoreactive molecules [9,13,53-55]. In some of these previous studies, it was possible to achieve a nanosecond timescale in the nonadiabatic dynamics, although the molecules investigated undergo decay processes towards the ground state in a much shorter timescale, typically hundreds of femtoseconds. To our knowledge, the potentialities of ML to run nonadiabatic dynamics for more complex systems with intrinsically slow photorelaxation processes such as in fluorescence remains to be explored.

Notably, the most advanced ML models designed for predicting potential energy surfaces already meet the chemical accuracy threshold (1.0 kcal mol−1 = 0.043 eV) [56]. However, this accuracy level depends on the complexity of the molecule and the availability of QC forces to better interpolate the potential energy surfaces in the training process. Moreover, a common problem found in the ML-PES literature is that the models' accuracy concerning energy predictions tends to degrade when going from the ground state to high-order electronic states [53]. This effect can be attributed to the more complex nature and higher diffuseness of the electronic density in excited states, as well as the dense packing of electronic states whose mixing may result in non-smooth PES for the excited states [1].

To illustrate some of the difficulties in learning excited state energies, we computed a dataset of 50 surface hopping trajectories for protonated 7-azaindole at the ADC(2)/cc-pVDZ [57] level (electronic supplementary material, SM-3) to train a kernel-ridge regression model with the RE descriptor, the so-called KREG model [58] (electronic supplementary material, SM-4). The discussion of the physical–chemical properties of this system will be reported elsewhere. Here, we focus exclusively on how well ML can predict its energies.

By plotting the difference between the ML-predicted and the ADC(2) energies as a function of time for each state (figure 1, top panel), one can see that the errors vary approximately in the range of −5 to 5 kcal mol−1 for all states along the first half of the testing trajectory. However, the disagreement between QC energies and ML predictions increases significantly after 600 fs, especially for S2 and S3. This inaccuracy is reflected in the RMSE, which increases by a factor of two when going from the ground state to the highest excited state S3 (see the bottom panels in figure 1). Note that the energy deviation for some specific geometries can be as large as 15 kcal mol−1 when the molecule is in S2 or S3 even though the corresponding RMSE values (averaged over the whole trajectory) are much smaller (2.61 and 3.18 kcal mol−1, respectively). The consequences of such a significant error for a quantitative description of the photodynamics of molecules with ML are a critical issue to be addressed in future studies.

Figure 1

Performance of the KREG model for predicting the ground- and excited-state energies of an entire MD trajectory (test set) of protonated 7-azaindole. (Online version in colour.)

Photochemical reactions in molecular systems are often characterized by multiple decay channels or pathways that may end up with considerably different photoproducts. Within the surface hopping approach, this behaviour can be translated into the proportion of ‘typical’ trajectories representing a particular set of electronic features or molecular geometries in dynamics evolution. Some of these ‘typical’ trajectories may be described as rare events because they represent only a tiny fraction of the total number of simulated trajectories. From the perspective of ML, this discussion raises a central question when developing robust and accurate models for nonadiabatic dynamics: the sampling strategy used to build the training set. Since ML is grounded on statistical principles, the sample of molecular geometries used as a training set should represent as close as possible the distribution of the entire (in principle, unknown) population. However, the unbalanced distribution of ‘typical’ trajectories and molecular configurations often found in nonadiabatic dynamics pose a challenge for constructing high-quality training sets.

The boxplots of figure 2, representing the RMSE distribution for each state of the protonated 7-azaindole dataset, sheds light on this problem. This statistical distribution is obtained using a ‘leave-one-out’ cross-validation scheme in the trajectory space, meaning that one trajectory is selected as the test set for each round of training and evaluation (electronic supplementary material, SM-4). As shown in figure 2, the RMSE values are spread in a wide range, and they increase when moving to high-order excited states (approximately 1.3 kcal mol−1 in S0 to 2.9 kcal mol−1 in S3). A few trajectories appear as outliers in the error distribution (see the black dots in figure 2). These findings indicate that geometries belonging to atypical trajectories are probably underrepresented in the training set. Different sampling strategies [59,60] have been proposed to remedy this issue and improve the quality of the training set, with active learning [61-63] being the most used in nonadiabatic dynamics applications.

Figure 2

Assessment of KREG model uncertainty for predicting ground- and excited-state energies of an entire surface hopping trajectory of protonated 7-azaindole. Each boxplot contains the distribution of RMSE values of energy predictions for one entire trajectory (test set) out of 50 independent MD simulations. The model was trained on 5000 molecular geometries randomly sampled from 49 trajectories after separating one trajectory for testing. This process was repeated until each one of the 50 MD trajectories were evaluated. The white dots in the boxplots represent the mean value of the RMSE distribution, while black dots indicate the outliers. (Online version in colour.)

Stability of integration methods

Even if we manage to substantially bring computational costs down enough to allow the long time-scale propagation, the simulation results must be reliable. Unfortunately, however, long trajectories mean more error accumulation. In the frame of classical nuclear motions, these errors are of two types, first coming from the integration of Newton's equations; second, coming from the zero-point energy (ZPE) spilling. We discuss them in the next subsections.

Errors in the integration methods

Long time-scale simulations put enormous stress on the integration algorithms. Longer trajectories require smaller time steps to reduce error accumulation [64] but at the expense of increased computational costs. This subsection assesses how the velocity Verlet algorithm [65], which is often employed to integrate the classical equation of motions (EOM), behaves in the long timescale.

For our tests, we employed the SBH model discussed in §2b. The model was set with 33 dimensions, equivalent to a molecule with 13 atoms (see details in electronic supplementary material, SM-2). We have run four trajectories of 1 ns each, to test different aspects. All trajectories started with the same initial conditions and were set as microcanonical ensembles. Because our goal is to check the stability of the classical EOM, the simulations were limited to adiabatic dynamics in state 2. The set-up of each trajectory is: The random fluctuations in trajectories 3 and 4 aim at assessing how sensitive the integration is to the precision of the electronic structure calculations. Precision plays a role whenever there is a chance that the electronic structure method may yield a different result when an identical run is repeated. This problem occurs in nonlinear systems converging to different local solutions due to tiny underlying numerical fluctuations [66]. In dynamics, it expresses as small, unphysical discontinuities on the potential energy surface.

Trajectory 1: time step 0.01 fs. This trajectory is taken as the benchmark.

Trajectory 2: time step 0.5 fs. It has the largest time step acceptable for integration of the classical EOMs.

Trajectory 3: time step 0.5 fs plus random fluctuation of Hartree/Bohr added to the SBH adiabatic forces, where r is a uniform random number between 0 and 1.

Trajectory 4: time step 0.5 fs plus random fluctuation of Hartree/Bohr.

Figure 3 shows the total energy conservation level in each trajectory in terms of the absolute total energy variation at time t, . The entire trajectory was sliced in intervals of 50 ps. For each interval, we computed the mean and standard deviation of . Unsurprisingly, with a tiny time step like Δt = 0.01 fs, an excellent level of energy conservation is obtained, with absolute variations of (7 ± 4) × 10−8 Hartree for the entire range up to 1 ns. With Δt = 0.5 fs, the variation increases significantly to (2 ± 1) × 10−4 Hartree. On a positive note, this level of variation remains constant during the entire trajectory, dissipating our initial fear of data degradation at long times. Random fluctuations of about 10−7 Hartree/Bohr do not impact the trajectory. However, random fluctuations 100 times larger, of about 10−5 Hartree/Bohr, degrade the long time-scale energy conservation, with absolute variations reaching a dangerous level of (1.0 ± 0.1) × 10−3 Hartree at 1 ns.

Figure 3

Distribution of total energy deviation as a function of time for different integration set-ups. Each dataset shows the mean (central curve) plus and minus one s.d. (shaded area) for points collected in intervals of 50 ps. Each set-up is characterized by the time step (0.01 or 0.5 fs), and random fluctuation (RF) added to the forces, 10−7 or 10−5 Hartree/Bohr. (Online version in colour.)

ZPE spilling

ZPE is a consequence of the quantum uncertainty principle, and thereby, an inherent error creeps into the propagation of classical trajectories in MQC methods by neglecting it. The classical nuclear motion approximation permits reactions below the quantum threshold of energy leading to unphysical results where products can be formed with less energy than ZPE in the internal vibration degrees of freedom. The problem becomes serious if the ZPE flows out of several modes and pumps into a specific weak bond to make it unphysically hot.

This ZPE spilling phenomenon is well discussed in the literature [67-70]. Nevertheless, it has been systematically neglected in simulations of ultrafast dynamics because of the assumption that it takes longer than the photochemical process to kick in. Whether this hypothesis is justified or not, we obviously cannot lean on that anymore when simulating dynamics in the long timescale.

Different strategies have been designed to reinforce ZPE maintenance during quasi-classical trajectories [71-79]. Most of them [71-75] require calculating instantaneous normal modes and Hessians, which may be cumbersome to be applied in on-the-fly dynamics propagation. In turn, Xie & Bowman proposed [80] a method that, in principle, can be applied to the on-the-fly dynamics of a general polyatomic system. Nevertheless, it has not yet been tested in this context.

Apart from the above discussed active ZPE constraint models, a few passive techniques have been suggested to overcome the ZPE leaking problem in classical trajectory calculations [76-79]. These models discard unphysical trajectories from the statistics that either form products with less than their ZPE or form products where the sum of the product vibrational energies is less than the total product ZPE. The drawback is that passive techniques may bias the results if any of the competing reaction pathways is more prone to undergo ZPE leaking problems than the others.

At this point, we agree with the conclusion of Guo et al. [69] that no satisfactory solution to the infamous ZPE spilling problem in classical dynamics has been found to date. Thus, novel techniques tailored to on-the-fly dynamics must be developed if we aim at long time-scale dynamics [81].

Accuracy of nonadiabatic dynamics

Suppose we managed to reduce computational costs and control the nuclear integration errors in the long time-scale propagation. Are the MQC methods available today tailored to deal with nonadiabatic phenomena in the long timescale? For instance, can surface hopping reasonably predict hops in weakly coupled regions?

This section demonstrates our first attempt to address these questions by investigating the performance of the decoherence-corrected fewest-switches surface hopping (DC-FSSH) dynamics in the long timescale. DC-FSSH results are compared with full-quantum results calculated with multi-configurational time-dependent Hartree (MCTDH) [82,83] and its multi-layer variant ML-MCTDH [84-86]. The tests are based on the SBH model discussed in §2b with 10 harmonic bath modes and Debye spectral density, constructed to describe the nonadiabaticity in the weak coupling limit [87] (see the electronic supplementary material, SM-2).

Methodologies for long timescale, nonadiabatic dynamics in high-nuclear dimensionality are so underdeveloped that we cannot rely on any method as an absolute benchmark. Even the state-of-the-art MCTDH cannot be ensured to perform well in such an unprobed domain. Bearing such limitations in mind, we performed long time-scale dynamics aiming to check whether the results yielded by the different methods are consistent between them. Given that MCTDH is a full-quantum methodology while DC-FSSH is a mixed quantum-classical ad hoc approach [88], we assume without proving that the former should still deliver the most accurate result.

The MCTDH ansatz expands the nuclear wave function into time-dependent single-particle functions (SPFs), where the latter are expressed as a linear combination of time-independent primitive basis. It immediately reduces the number of basis functions required for converged calculation by providing the variationally determined basis for an optimal description of the evolving wave packet. The Dirac–Frenkel variational principle is used to derive coupled nonlinear differential equations as the MCTDH EOM, usually solved by the predictor–corrector integrator scheme.

ML-MCTDH is an extension of the MCTDH approach, in which the multi-dimensional SPFs are propagated themselves within the MCTDH ansatz. The standard propagation method describes a single layer of expansion coefficients for the time-independent basis. By contrast, MCTDH contains a first layer of expansion coefficients for the time-dependent SPF basis and the second layer of time-dependent expansion coefficients that parameterize the time evolution of the SPFs. One can expand the multi-dimensional time-dependent SPFs again using an MCTDH ansatz yielding a three-layer scheme. In other words, the ML-MCTDH method prescribes to apply the MCTDH scheme successively in different layers to propagate the SPFs.

The Heidelberg MCTDH package is used to simulate MCTDH and ML-MCTDH [89,90] dynamics on the 10-dimensional SBH model system. Here, we limit our discussion by presenting only the results keeping the number of SPFs and the grid size Details of the MCTDH and ML-MCTDH set-ups are given in the electronic supplementary material, SM-5.

Trajectory surface hopping [1,91] is a widely known MQC dynamics methodology, which expresses the nuclear wave packet on a particular electronic surface as a swarm of independent classical trajectories obeying Newton's equations, and the force acting on the nuclei is derived as the negative gradient of the potential energy of the corresponding electronic state. At each time step, a stochastic algorithm determines whether the system will propagate on the current electronic state or hop to another one leading to nonadiabatic population transfer between the electronic states. Among various strategies to compute transition probability, Tully's FSSH formulation is probably the most common and extensively reviewed in the literature [92].

We have used Newton-X NS (novel series; see §5) to carry out the DC-FSSH dynamics (electronic supplementary material, SM-6) on the SBH model in the adiabatic representation. Since we performed the DC-FSSH in the adiabatic representation and MCTDH in the diabatic representation, comparing the electronic population obtained from these two methodologies would be unjustified. Thereby, as Landry et al. [39] prescribed, we have calculated the diabatic electronic populations from a swarm of DC-FSSH trajectories using an orthogonal adiabatic to diabatic transformation matrix and mixing angles [93-95]. The diabatic population decay curves are fitted with a bi-exponential function to obtain the time constants associated with the nonadiabatic process, where is the remaining population that does not decay in these timescales. Since the bi-exponential function could not fit the population decay curve obtained from the MCTDH dynamics, the following mono-exponential function is also used to fit only the MCTDH population curve

ML-MCTDH and DC-FSSH yield closely matching time constants, predicting similar kinds of dynamics undergoing with two different timescales (figure 4 and table 1). However, it is surprising to see that, unlike the multi-layer variant, the MCTDH population curve has a single time constant. The long time constant is due to a repopulation of the excited state rather than to a state trapping. We could not yet determine why it shows up in ML-MCTDH but not in MCTDH. We are proceeding with the analysis running new simulations, but discussing these results is out of the scope of this paper. It will require a separate publication exclusively dedicated to this subject. Despite these differences, it can be safely stated that the calculated time constants of the non-radiative relaxation process obtained from quantum dynamics and MQC dynamics are of the same order of magnitude, validating the latter method even in the long timescale.

Figure 4

The diabatic population decay of the excited state for different dynamics methodologies are shown in solid lines, and the corresponding fitted curves are presented in dashed lines. (Online version in colour.)

Table 1

Fitted parameters for the diabatic populations in the long time-scale (100 ps) dynamics of an SBH model. A and are not defined for MCTDH.

	A	τ_1 (ps)	τ_2 (ps)	N0
MCTDH	—	7.85	—	0.003
ML-MCTDH	0.34	7.07	23.71	0.024
DC-FSSH	0.50	6.45	23.31	0.028

Software optimization

Originally, Newton-X, our software program for dynamics simulation on excited states, was developed in the frame of ultrafast processes [96,97]. For this reason, in the choice between computational efficiency and flexibility (modularity to communicate with different electronic structure methods), most of the time, the latter was favoured. Moreover, it is not fit for processing and storing big data. Simulations in the long timescale invert the priorities. Therefore, efficiency and big data become mandatory issues to be addressed.

We have recently started to write a completely new Newton-X version to deal with these issues. The legacy program is now referred to as Newton-X CS (for classical series), while the new program is referred to as Newton-X NS (for novel series). We will present the details of the Newton-X NS implementation in a future publication, but here, we would like to highlight a few points about how it will help move towards long time-scale simulations.

From a technical standpoint, performing nonadiabatic dynamics with FSSH amounts to repeating the following sequence: Going between steps requires conveying information either internally (from (v) to (vi), for instance) or through some external program ((i) to (ii), for example). A straightforward approach is to use a primary driver that will call a set of scripts or programs sequentially, each communicating through text files. Although this approach maximizes modularity, it has clear limits in terms of I/O operations.

Read initial conditions (geometry, velocity, initial electronic state and initial time-dependent wave function coefficients).

Compute electronic energies and energy gradients (usually through an external code).

Obtain time-derivative couplings (either nonadiabatic coupling vectors, through wave function overlap matrices, or approximated methods [98]).

Solve the classical EOM.

Solve the quantum EOM (a locally approximated time-dependent Schrödinger equation).

Decide on whether to hop to another state or not.

Go back to step (i) with the updated set of conditions.

Historically, step (ii) has been the time-limiting step in on-the-fly dynamics simulations, as performing electronic structure computations is much more resource-intensive than the EOM integration. Indeed, most operations in steps (iv) and (v) scale with the number of atoms or the number of states included in the system, while step (ii) scales with a power of the number of electrons, making the I/O operations effectively costless. With the development of fast electronic structure calculations highlighted in §2, this relation does not hold anymore, and the dynamics approach must be adapted to limit the number of I/O operations. This algorithmic restructuring is especially crucial for long time-scale dynamics, where we must tightly optimize the speed of each step. Replacing a set of scripts with a more integrated program can solve this issue by keeping most data in memory during the dynamics and only writing on disk at specific points when unavoidable, like when communicating with an external program.

To deal with this problem, we have entirely restructured Newton-X. While Newton-X CS is a collection of Perl scripts with Fortran programs, Newton-X NS is a complete Fortran program using Perl only for dealing with QM input files. To show the performance of the Newton-X NS, we considered a 1 ps dynamics using the 10D SBH model (electronic supplementary material, SM-2), integrated with a time step of 0.1 fs (resulting in 10 000 simulation steps). In Newton-X CS, such a trajectory takes 20 400 s (5 h 40 min) to complete, while Newton-X NS computes it in only 11 s, astonishingly almost 1850 times faster. Of course, the increase in performance when interfaced with external programs would be less impressive, as the time for the electronic structure computation is likely to be dominant in this case.

The second software-related problem we must deal with when planning for long time-scale dynamics is data storage. Molecular dynamics simulations output a gargantuan amount of data, as the outcome is primarily a long sequence of geometries and energies. This problem is even more pronounced in FSSH, as the algorithm's stochastic nature requires an ensemble of many trajectories. We can make a simple estimate of what these data requirements would represent in the long timescale. Let us assume we performed simulations for a molecule with 20 atoms for a 1 ns trajectory with 0.5 fs time step. The storage of the geometries alone (in double precision) would require around 1 GB of space. Standard FSSH studies routinely need around 100 trajectories to achieve satisfying statistical accuracy, so the geometries alone could take up to 100 GB of space. To this amount, we should still add electronic populations, energies, nonadiabatic couplings and surface hopping-related data (like the transition probabilities). This dataset must be further post-processed to deliver physical–chemical information and made publicly available to comply with open data requirements. It is thus of crucial importance that the data files are uniquely identified and that the set-up for producing the data (software version, operating system, computer used) is also available.

The current state in Newton-X CS is to output these data as a series of text files. Although convenient to write and read, text files have some serious drawbacks. For each new dataset, a new file must be created and documented. This can quickly lead to a multiplication of outputs, making data integrity harder to ensure and data sharing more difficult. Therefore, obtaining the output as a single file containing the entire dataset produced, as well as the condition of production as metadata, would be a desirable feature. This can be achieved with a structured file format. However, the context of MQC dynamics makes the use of structured text formats cumbersome.

In MQC dynamics, we save data as time series. In practice, this means that the effective rank of each dataset is increased by one. Although conventional formats exist for saving two-dimensional data like molecular geometries as a time series (xyz format, for instance), keeping track of data with higher rank, such as energy gradients when several electronic states are involved, quickly becomes cumbersome and decreases the readability associated with text files. Another problem lies in the analysis of trajectories. Each time step will produce heterogeneous data, which must be saved when they are still available in memory, but for convenience of analysis, it is desirable that the data corresponding to the different sets (energies, populations, …) are kept together. Structured text formats exist, like JSON or YAML, but they are not suited for a sequential update of the datasets. One solution would be to save the time steps as a parent field, containing the heterogeneous data as children, post-process the file at the end of the simulation to gather all datasets and produce a new structured file with consistent datasets.

These issues can be solved using the hierarchical data format (HDF) model, which allows the creation of portable, self-described and open binary file formats. The model has clear advantages [99], with complete C, C++ and Fortran APIs, and officially supported bindings with Python, Julia and R. The files created with the underlying library can be easily explored and shared and are self-documented, with the inclusion of metadata directly along with the datasets. Binary files also occupy less disk space than text files, and HDF5 can add compression to the data to lower the footprint even more. We have chosen to adopt the H5MD format [100]—specifically designed to handle molecular dynamics—to Newton-X NS. Using this file format allows us to benefit from all the environments already developed.

For the same simulation as the one showed earlier in this section, set to write the maximum amount of output and using Newton-X NS, the program spends around 3–4 s to write output files, showing no difference between text and H5MD files. However, the files obtained differ in size: the text files together occupy 58 MB of disk space, while the H5MD file occupies only 28 MB. It is also worth noticing that HDF5 files can be read in chunks, thus facilitating the treatment of large datasets, and can be opened only partially just to read a header, creating exciting opportunities for automated archiving and research in databases.

Conclusion

Nonadiabatic dynamics simulations in the long timescale (defined here as much longer than 10 ps) have been an overlooked research area. Several recent works in the field indicate that this situation may be about to change. In this paper, we lay down the elements that should be considered for advancing the research in this field.

First, we delimit the scope of what we expect from methods to run such long time-scale simulations: they should work in full nuclear dimensionality, be general enough to tackle any type of molecular system and be computationally affordable. Bearing these constraints in mind, we examined the main methodological challenges we should face to advance the field: the computational costs of the electronic structure calculations, the stability of the integration methods, the accuracy of the nonadiabatic dynamics algorithms and software optimization. We discuss each of these issues based on simulations of model systems and actual molecules.

We analysed three strategies to speed up electronic structure calculations, using parameterized and approximated methods, employing model Hamiltonians and predicting electronic properties with ML. The first strategy is consistently associated with strong accuracy loss. The second has the main challenge of parameterizing full-dimensional systems, a procedure that has found an ally in ML. For the third strategy, we believe ML may be on the verge of becoming a crucial element for long time-scale dynamics. Nevertheless, it still misses clear protocols for training excited-state energies of large conformational sets and high density of states with the accuracy needed for dynamics.

We showed that conventional methods for integrating the classical EOM simulations should be adequate to extended simulations up to 1 ns, as long as the precision of the electronic structure properties is kept under control. Nevertheless, ZPE leaking must be accounted for. Hessian-free methods for correcting ZPE leaking need to be developed.

We evaluated the performance of decoherence-corrected surface hopping in the long timescale. The results for internal conversion in the weak-coupling regime show that it can make predictions in semiquantitative agreement with wave packet propagation.

Finally, we also described our optimization of the Newton-X program to reduce computational overheads in processing and data storage. In the most favourable case, computational time can be reduced more than 1000 times compared with the costs of simulations done on the current software architecture.