OMNI-P2x universal neural network potential for excited-state simulations

Domain: neuroscience
Published: 2026-04-07
Editorial outcome: minor_revision

Abstract
Photo-active molecular systems play an essential role in modern science and technology, finding applications in solar cells, organic light-emitting diodes, reaction catalysis, photodynamic therapy, and beyond. The rational design of photo-responsive molecules requires an understanding of the photophysical and photochemical processes underlying their operation, which can be gained via the first-principles quantum-mechanical calculations that are prohibitively expensive for high-throughput investigations. To break through this limitation, here we introduce OMNI-P2x: a universal neural network potential for molecular excited and ground electronic states. OMNI-P2x can be used, directly or after fine-tuning, to perform a wide range of photophysical and photochemical simulations. OMNI-P2x approaches the accuracy of time-dependent density functional theory methods at a fraction of the computational cost, while being more accurate and faster than established semiempirical methods for excited-state simulations. Here, we demonstrate its use in UV/Vis absorption spectroscopy, real-time photodynamical simulations, and in the rational design of visible-light-absorbing azobenzene systems.

Introduction
Theoretical photochemistry leverages the ability of quantum-mechanical (QM) methods to accurately predict the properties of electronic excited states of molecules and to provide insights into light-molecule interactions. The QM calculations facilitate better understanding of complex processes such as photoreactivity and photostability [81–84], operation of molecular machines [68, 69], separation of charge carriers in clean energy generation [99, 100], and, in turn, provide rational design principles for photomaterials such as organic light-emitting diodes (OLEDs) [75, 76], molecular solar thermal energy storage (MOST) systems [64, 65], and organic solar cells [78, 79]. Because of the extensive development of computational methods, they can be not only used to explain and verify experimental results, but also to guide novel research in photochemistry.

At the same time, the key obstacle preventing widespread use of theoretical methods for excited-state property prediction is their computational cost, with a schematic representation of their accuracy/computational cost tradeoff presented in Figure 1. While highly-accurate methods, such as equation-of-motion coupled clusters approaches [26] and their variants like ADC(2) [2] exist, they can only be used for medium-sized systems at best, requiring days or weeks of computational time due to their steep computational scaling. On the other hand, time-dependent density functional theory (TD-DFT) [25, 91] provides a more affordable alternative, which typically scales up to systems containing about 100 atoms, however, at a price of decreased reliability. Finally, on the lower end of the accuracy-to-computational-cost spectrum one finds semi-empirical QM methods (SQM), such as a wide family of Hamiltonians based on the neglect of diatomic differential overlap (NDDO) approximation [71] such as AM1 [6], PM3 [80], OM2 [35], and ODM2 [36] used in conjunction with the configuration interaction (CI) treatments, and time-dependent tight-binding DFT methods TD-DFTB [31, 32]. While these methods are certainly more affordable for systems with hundreds and even thousands of atoms, the accuracy of their predictions can often be questionable. [74, 104, 106]

The situation becomes even more complex when real-time simulations of the light-induced processes are required, i.e., for predicting photoswitching quantum yields, two-photon dynamics, or excited-state lifetimes. The most popular and widespread technique for performing such simulations is trajectory surface hopping (TSH) [22] variant of nonadiabatic molecular dynamics (NAMD). In TSH, a swarm of independent, semi-classical trajectories are propagated, with the ability to switch (hop) between adiabatic potential energy surfaces (PESs). The electronic and nuclear degrees of freedom are coupled by so-called interstate couplings, which control the hopping probabilities. Running surface-hopping simulations requires using a QM method that can accurately predict multiple PESs and, ideally, properly describe the multireference character of the wavefunction in the vicinity of conical intersections. This severely limits the choice of theoretical methods and, typically, CASSCF is used, while much slower but higher-accuracy perturbational variants such as CASPT2 [120] are also available. However, as driving dynamics requires calculations of energies and gradients at each time step, with a typical ’minimum’ simulation setup of 100 trajectories, each propagated for 1 ps with a time step of 0.5 fs, requires computing these quantities 1 mln times, which is only possible for simple, often model systems. For simulations of experimentally and technologically relevant compounds, the only available choice are, so far, semi-empirical methods.

In recent years, progress in the field of machine learning (ML) has led to a significant paradigm shift in the area of computational chemistry [15, 23, 47, 116, 119]. Fitting machine learning potentials (MLPs) to reference QM data allows for accurate predictions of molecular properties, such as energies and energy gradients [33, 41, 85, 118, 121], dipole moments [85] or electronic densities [92, 94] at a greatly reduced cost, while retaining the accuracy of the reference method. This allows to break through the limitation of traditional QM methods, making simulations of larger, more complex systems possible at longer timescales [16, 39, 97, 123, 126]. A particularly important, emerging class of potentials is universal MLPs, which we define as models that were trained on vast datasets spanning across chemical space. In turn, they can be used to make out-of-the-box predictions for previously unseen molecules, generalizing from the training data. [8, 27, 28, 67, 77, 128]

Although training and using MLPs for predicting ground-state properties is becoming a rather routine task, using ML methods for electronic excited states is far more challenging, requiring expert human supervision. This is due to the inherent complexity of excited potential energy surfaces, often forming a strongly correlated manifold of states, interacting and mixing with each other in a non-smooth manner. At the same time it should be noted that models capable of providing accurate predictions for excited-state properties, such as SchNarc [96], SpaiNN [102], or PyRAI2MD [132] have been developed and used successfully in the TSH simulations. However, their capabilities are limited to learning one molecule at a time, i.e., they are not transferrable across chemical space.

This has changed with the MS-ANI (Multi-state ANI) model [66], demonstrated to be able to predict nonadiabatic dynamics of three distinct molecules within the same model, paving the way for universal excited-state potentials. More recently, another demonstration of a possibility to create transferable NN for excited-state predictions was shown in Ref. [114]. Transferable nonadiabatic molecular dynamics simulations were also reported by Axelrod et al [18], however, it was limited to the class of diazobenzene photoswitches. Despite all the progress in MLPs for excited states, they have been still lagging behind the ground-state MLPs to such an extent that many experts are questioning whether universal MLPs for excited states are even possible. [49, 109] It is generally agreed that developing such an MLP would be a major breakthrough. [49, 119]

In this work, we present OMNI-P2x: a universal neural network (NN) potential for molecular excited states (Figure 1). We have built OMNI-P2x by combining the ideas of multi-state [66, 110] and all-in-one [29] learning. OMNI-P2x approaches TD-DFT accuracy for excited-state energies and targets CCSD(T) quality for ground-state energies: it can predict the absorption UV/Vis spectra out-of-the-box with higher speed and accuracy than popular semiempirical QM methods, making ML approaches for excited-states a viable alternative to established, cost-efficient QM methods for high-throughput calculations, such as semi-empirical methods.

We leverage the power of OMNI-P2x to perform a high-throughput screening of 500,000 compounds to find azobenzene derivatives featuring red-shifted, visible light absorption. These systems are particularly relevant due to their applications in nonlinear optics [10], photoresponsive polymers [11], and as MOST compounds [12, 13]. The potential of OMNI-P2x in accelerating another type of extensive simulations: nonadiabatic dynamics, is also exploited via fine-tuning of this universal potential within active learning loop. This approach reduces the computational cost and wall-clock time of production-ready ML-TSH simulations, while improving data efficiency and simulation stability.

Results
2.1 Model architecture and training

The architecture of the OMNI-P2x model is based on the idea of multi-state learning[66, 110] and its successful realization in the MS-ANI model designed for accurate and efficient learning of excited-state energies and forces. Multi-state learning can be used to efficiently learn an arbitrary number of electronic states while using the entire neural network to differentiate the molecular electronic states, as opposed to only the last layer used in the multi-output models. This effect is achieved by including an additional state ordering number as an input feature, along the molecular descriptor containing information about the geometrical conformation. Here, we extend this model to learn across multiple QM levels, in the spirit of the all-in-one (AIO) learning paradigm introduced recently [29]. The all-in-one learning has been successfully used to develop the universal OMNI-P1 neural network potential capable of making predictions approaching DFT and CCSD(T) levels, however, it is limited to the ground-state energies and forces. It is worth noting that the ability of learning from multiple levels of theory gives an important advantage of leveraging different data sets of various quantities. Here, we combine both: multi-state learning and all-in-one learning, to obtain OMNI-P2x: the omnipotential-2x, all-in-One Machine-learning Neural-network Interatomic Potential 2 for eXcited states.

OMNI-P2x is based on the established ANI NN potential [7] (Figure 2a). The input layer consists of an atomic environment vector (AEV), calculated using the TorchANI Python package [105], and two additional features: the state ordering number, as in MS-ANI, and a one-hot encoded feature representing the QM level, as in all-in-one learning. The electronic state is encoded as a non-negative integer, consistent with the natural ordering of adiabatic state energies, which we find to be both physically motivated and empirically more accurate than one-hot encoding. In the present work, data at two levels of theory were included: the ground-state-only data with CCSD(T)/CBS accuracy from the ANI-1ccx dataset[7], and TD-DFT B3LYP/6-31+G* excited- and ground-state data from the PubChemQC dataset[88]. As for all ANI-type models, a separate neural network is trained for each chemical element. In the present version, OMNI-P2x supports the following element types: H, C, N, O, F, S, and Cl. The predictions are then made atom-wise and summed up to yield the total energy of the molecule. For the energy part of the predictions, an ensemble of three NN potentials is used, with their predictions averaged to obtain the final predicted energy for all electronic states. A fragment-correction scheme is applied to non-covalently bound structures, to recover the correct, intensive character of excitation energies, for molecular fragments separated by distances greater than the descriptor cutoff radius, with details provided in the methods section, and Suppl. Note. 1. A separate MS-ANI-type NN is trained and used to predict the oscillator strengths of interstate transitions, as shown in 2b). Since only TD-DFT data is used for the excited-state training, this network does not include an additional input feature for the level of theory.

The data for model training has been sourced from two datasets: 1) excited-state data in the PubChemQC project’s dataset [88], and 2) ANI-1ccx, containing about 500,000 ground-state geometries, including out-of-equilibrium conformations, with their total energies computed at the gold-standard CCSD(T)/CBS level of theory. PubChemQC project’s dataset consists of unique, optimized structures of molecules at the TD-B3LYP/6-31+G* level of theory, [1, 20] with the energies and oscillator strengths of the first 10 excited states calculated using TD-DFT. We have used the entirety of the PubChemQC dataset, which meets our filtering criteria (as detailed in the Methods section), yielding 3.1 million molecules with 11 electronic states each. The inclusion of the ANI-1ccx data set provides an important additional information about the effect of out-of-equilibrium geometry changes on energies, which is important for molecular dynamics and other types of simulations. ANI-1ccx data set has been used as a basis for training multiple universal ML-based models [3–5, 8, 9, 107], which are successfully used in molecular dynamics, geometry optimizations, and thermochemistry calculations. In addition, combination of two different types of data set allows to cover larger chemical spaces, i.e., compare the data distribution in both data sets in Figure 3a: the ANI-1ccx data cover small-sized molecules, while the vast majority of the excited-state training data is comprised of medium-sized molecules containing 15 to 45 atoms, with most of their excitation energies located towards the blue end of the visible spectrum, and in the UV range (Figure 3b).

We evaluate the performance of the model using a test set of 300k unseen molecules from the PubChemQC dataset for the excitation energies, as well as a test set of almost 49k ground-state conformations from the ANI dataset. The MAE in energy predictions for the first 11 electronic states with respect to the TD-DFT data shows that the quality of the predictions increases when going to higher excited states, reaching a minimal value for S7 (Figure 4b and c). An outlier here is S0, which has the smallest energy MAE of 1.6 kcal/mol, which we attribute to the improved data-efficiency due to the application of the all-in-one learning to the combined set with the ground-state DFT (from PubChem) and the CCSD(T) (from ANI-1ccx) data, and an overall smoother and lower-variance distribution of ground-state energies that makes them easier for the model to learn. This ground-state MAE can be compared to the similar metrics of other MLIPs: ANI-1 [7] with an RMSE of 1.3 kcal/mol (OMNI-P2x 2.93 kcal/mol), UMA [108] with around 1 to 3 meV/atom (0.023 to 0.069 kcal/(mol · atom), depending on the version (OMNI-P2x 0.11 kcal/(mol · atom) and AIMNet2 [17] reporting a MAE of around 1.25 kcal/mol for neutral compounds. Although those results are not comparable one-to-one, as each of the models was trained and tested on a different dataset, they show that the ground-state error of OMNI-P2x remains within a reasonable range from current state-of-the-art machine-learning derived models, despite not being the main focus of this study. Correlation plots showing the predicted vs true excitation energies for all transitions are available as Suppl. Fig. 4. In all cases, the R2 is above 0.91. For the ground-state energies in the test ANI-1ccx subset, OMNI-P2x has a higher MAE of 4.14 kcal/mol, but a high R2 = 0.9966 for the normalized energies (0.9999 for unnormalized) (Figure 4a). This is anticipated, as the ANI-1ccx training data, which is only used as an auxiliary data set, is much smaller than PubChemQC, and is also featuring a different data distribution (Figure 3a).

Additionally, OMNI-P2x provides an uncertainty quantification (UQ) of all predicted energies, using the standard deviations between quantities predicted by each of the three NN potentials as the error bar. A warning is issued to the user if the UQ of any excited states exceeds a safe threshold, which is determined as the median plus three median absolute deviations of the UQs of the test set, for each electronic state.

Finally, the MAE in OMNI-P2x predictions of oscillator strengths is below 0.04 for all studied electronic states (Figure 4d). One can note, that the highest accuracy in terms of the predicted oscillator strengths is achieved for the excitations to S1 and to higher-lying S8–S10 states.

2.2 UV/Vis absorption spectra

OMNI-P2x is a universal model that can predict both excitation energies and oscillator strengths, i.e., which allows one to obtain the UV/Vis excitation spectra of the QM quality, yet with a much higher speed (see Methods for details).

To test the OMNI-P2x performance for simulating the UV/Vis spectra, we generate the spectra for 50k test molecules from the PubChemQC database and compare them to the reference TD-DFT spectra (see Methods). Let us first examine representative spectra shown in Figure 5. Looking at a spectrum belonging to the first quartile (top 25% spectra with the highest similarity between the OMNI-P2x and TD-DFT based on the similarity score, see Methods and below for details), agreement is observed both: in the intensity and position of the peaks (panel a in Figure 5). For the second quartile, the agreement becomes somewhat worse, which is manifested by the absorption energy red-shift of about 10 nm, however, the intensity of the peak is consistent with the reference method (panel b). This agreement becomes worse in the bottom half of the predicted spectra: while the positions of the predicted peaks are rather accurate, the intensities are different from the reference method. Inspecting the bottom quartile of the predicted spectra, the agreement is still acceptable, with the higher absorption maximum predicted quite accurately, both in terms of the absorption energy and intensity, however, the first absorption peak appears to be blue-shifted due to the incorrect prediction of S2 as a bright state. To provide a benchmark of the predicted spectra that is unaffected by the choice of a broadening factor in the SPC construction, we also rank the spectra by a geometric mean of the excitation energy and oscillator strength MAE, which is presented in Suppl. Fig. 5. Later, we also explore the application of OMNI-P2x for spectra calculations without relying on SPC approximation and the broadening factor.

For the test molecules shown in Figure 5, we tested the effect of finite-temperature distortions on the accuracy of OMNI-P2x and found that higher excited states remain close to equilibrium benchmarks, while the accuracy of the lowest excitations (S1–S3) may decrease more sharply. Violin plots showing the error distribution on non-equlibrium structures are available in Suppl. Figures 6–9. To further probe this behavior, we performed relaxed PES scans of cyclohexanone (C=O stretch) and azobenzene (C–N=N–C torsion), shown in Suppl. Fig. 10. These scans confirm that OMNI-P2x is accurate near equilibrium and along moderate distortions, but errors grow for strongly compressed geometries or when torsion induces a change in excited-state character. Together, these results indicate that OMNI-P2x can extrapolate PES shapes close to equilibrium, but its reliability decreases as the molecular distortion drives the system beyond the domain of the training data.

To put the performance of the universal potential into perspective, we also evaluate two commonly used semiempirical QM methods: ODM2/CIS and TD-DFTB, on the same test of 50k molecules, as shown in Figure 6. The OMNI-P2x spectra are clearly more similar to the reference TD-DFT than the semiempirical ODM2/CIS and TD-DFTB as judged by the distribution of the Spearman correlation coefficients for this test set (Figure 6a). On average, the OMNI-P2x and TD-DFT spectra agree well with each other, with an average correlation coefficient of 0.86, and 80% of the predicted spectra having a correlation coefficient larger than 0.80. In comparison, both the semiempirical QM methods have a substantially worse performance, with an average SCC of 0.53 for TD-DFTB, and 0.50 for ODM2/CIS. Additionally, ODM2’s parameterization is limited to molecules consisting of H, C, N, O, F atoms only (39164 out of 50k structures), i.e., OMNI-P2x has a broader element coverage than ODM2/CIS. Noteworthy, the OMNI-P2x is not just more accurate than semiempirical QM methods, but is also clearly the fastest of the available approaches, with an average calculation time of 225 ms per predicted spectrum, while the ODM2/CIS and TD-DFTB approaches take on average 827 ms and 1152 ms, respectively (Figure 6b). The reference TD-DFT calculations are orders of magnitude slower: more than 357.87 s per predicted spectrum. While the QM methods follow an approximately exponential scaling of cost with respect to accuracy, measured as correlation of the predicted spectra with TD-DFT/B3LYP spectra (note the logarithmic scale on the y axis), OMNI-P2x breaks this unfavorable scaling.

2.3 Nuclear ensemble approach spectra

A more accurate method for predicting UV/Vis absorption spectra is the Nuclear Ensemble Approach (NEA). In this technique, vertical excitation energies and oscillator strengths are computed for an ensemble of hundreds or thousands of molecular geometries [72, 73]. The final spectrum is constructed by summing Lorentzian or Gaussian functions centered at each geometry’s excitation energy, weighted by the corresponding oscillator strength. Compared to single-point convolution, this approach better captures vibronic effects. However, NEA simulations are computationally demanding, as they require quantum mechanical calculations for all structures in the ensemble. To address this cost, ML models offer a natural path for acceleration [52, 60, 61]. Since OMNI-P2x is initially trained only on equilibrium structures, it is not directly applicable to NEA simulations. Nevertheless, fine-tuning enables it to achieve high data efficiency and simulation speed in this context, requiring only 25–50 training points to reproduce reference spectra.

To demonstrate this, we consider two test cases: pyrazine (with 4 excited states) and 9-dicyanomethylene derivative of acridine (9-DCMA, with 30 excited states), the latter of which was previously used to benchmark ML-NEA spectra with kernel ridge regression models [60]. The NEA spectra simulated using reference TD-DFT data, as well as those predicted by the fine-tuned OMNI-P2x model—trained on just 25 data points for pyrazine and 50 for 9-DCMA—are shown in Figure 7. For pyrazine (panel a), the agreement is observed across the entire spectrum. In the more complex case of 9-DCMA, the model accurately captures the intensity and position of the lower-lying excitations. Deviations start emerging below 250 nm, especially in predicted intensities, while the position of the peaks remains rather accurate. We attribute this to the model extrapolating in a regime beyond the number of excited states included in training dataset: OMNI-P2x was pre-trained on 10 excited states, while 30 are considered in this example. Such simulations can be done on commodity hardware (Jupyter notebook using 4 CPUs), in under a minute for pyrazine, and about 40 minutes for 9-DCMA, excluding the labeling time of QM methods.

2.4 Interpertability of ML predictions

To provide interpretability of the ML-learned excitations, we introduce a new approach for analyzing atom-wise contributions to excitation energies. For each atom i, the fractional contribution is defined as |∆Ei|/∆E, where ∆E is the excitation energy for a given transition, and atomic contributions are calculated as the differences between ground- and excited-state atomic energies predicted by the NN ensemble. We show that the OMNI-P2x-derived atom-wise contributions agree well with the reference TD-DFT density difference plots (Fig. 8) Such analysis allows us to obtain the spatial character of the excitation by identifying which atoms are most involved in the transition. This analysis is conceptually related to work by Tkatchenko et al. [98], analyzing atom-wise contributions to ground-state energies.

We demonstrate that this interpretability analysis can work well across chemical space and for different excitations of the same molecule. Panel a) of Fig. 8 shows the S0 → S1 ML ∆E contributions and TD-DFT density differences for several organic molecules. We begin by comparing the examples of benzene and pyrazine, two 6-π electron aromatic compounds. OMNI-P2x correctly localizes the transitions, spreading it evenly across all the benzene carbon atoms, while in pyrazine it is mostly localized on the two nitrogen atoms, reproducing the TD-DFT results. The excitations also exhibit appropriate symmetry: C6 for benzene, and D2h for pyrazine. The next example, azobenzene, shows physically correct localization of the excitation on the N=N double bond, while in p-nitroaniline the partial charge-transfer character of the excited state is well reproduced. Overall, this shows an interpretable character of ML predictions, that is transferable across chemical space.

Atomic contributions to different electronic transitions of the same molecule are shown in Fig. 8b), using azobenzene as an example. OMNI-P2x correctly reproduces the S0 → S1 transition as localized on the central N=N bond, while the S2 transition has larger contributions from the aromatic rings. Finally, a higher excited state (S11) is shown, which is delocalized across the entire conjugated system.

The means of analysis presented above do not exclude the need for dedicated analysis, either manual or using wavefunction-based tools such as TheoDORE [86] or Multiwfn [70], which can provide additional insights beyond spatial localization of the excitations, such as the orbitals’ character. These tools are complementary, as massive-scale simulations enabled by OMNI-P2x can serve as a first step for generating enough statistical data, from where representative samples can be drawn, and post-processed further with the wavefunction-based tools.

2.5 Design and screening of visible-light absorbing azobenzene compounds

The speed and accuracy of OMNI-P2x lend themselves for performing efficient high-throughput screening (HTS) of photo-active molecular systems. To demonstrate an application of OMNI-P2x for such a scenario, we perform the HTS of azobenzene derivatives characterizing with low-energy absorption. Our choice is motivated by the relevance of such compounds for applications in nonlinear optics [10], photoresponsive polymers [11], and as MOST systems [12, 13]. To this end, firstly we have generated a total of half a million azobenzene derivatives (Figure 9) as SMILES strings, using RDKit [87]. In this procedure, we have allowed up to five substitutions with the following functional groups: methyl, ethyl, trifluoromethyl, nitro, hydroxyl, thiol, fluoro, chloro, amine, nitrile, and carboxyl. 3D structures (Cartesian coordinates) have been generated from SMILES using MLatom [62, 130], and optimized using the fast and accurate AIQM3 model [5] (see Methods for details). Subsequently, we have predicted the excitation energies of the optimized structures with OMNI-P2x. To refine the quality of the OMNI-P2x predictions, 1000 random structures have been sampled from the generated set and labeled using the TD-DFT/B3LYP method (see Methods for details on the level of theory). We have then fine-tuned two copies of one of the OMNI-P2x energy models on this set, using different training/validations splits. The procedure of prediction, selection, and optimization has been repeated two more times, at each iteration adding the 1000 most uncertain molecules to the training set, with the uncertainty measured by the deviation between the two trained models. When a correlation coefficient of 0.78 and a MAE of 0.25 eV have been reached as tested on 1000 randomly sampled structures, we have stopped the first phase of the screening procedure, as this accuracy is similar to the accuracy of the target TD-DFT methods. For 1000 lowest-energy absorbing candidates, as predicted by the fine-tuned model, TD-DFT predictions have been made again. Out of these structures, a total of 21 lead candidates with transition energies below 2.0 eV have been determined (see Supporting Material, SM). Due to the potential charge-transfer character of excitations in these molecules, which the B3LYP functional would be unable to capture correctly, we have used at this stage a higher-level electronic structure method, ADC(2) [2], to re-evaluate the excitation energies of the 21 lead candidates.

The three most red-shifted derivatives identified with the described procedure are presented at the bottom of Figure 9, while a database containing all the structures and their computed UV/Vis spectra is available in the SM. A plot comparing the spectra predicted at the three studied levels of theory is available in Suppl. Note 3. The predicted absorption energies of the top three derivatives are lower than absorption energies of reported human-designed red-shifted azobenzene units with the same core structure [45, 54, 112, 113]. At the same time, our top candidates determined via in silico screening show similarities to substitution patterns of the human-designed derivatives, i.e., ortho- and para-placements of electron-donating groups [38, 51].

2.6 Efficient ML-driven nonadiabatic dynamics

As a universal model, OMNI-P2x has the potential to accelerate nonadiabatic dynamics (NAMD) simulations, similarly to how ground-state models are used to efficiently perform ground-state dynamics. However, in this respect it should be noted that NAMD is not just much more computationally demanding than the ground-state dynamics, but also typically requires a system-specific selection of a level of theory, as the common black-box approaches, like TD-DFT, cannot properly describe the key features of excited-state PES, i.e., the first excited to ground-state PES intersection [59]. This often necessitates a suitable choice of more advanced QM methods, such as CASSCF or CASPT2, which require the active spaces specifically selected for a given system. Since our OMNI-P2x is only trained on TD-DFT data, here, instead of directly using it to propagate NAMD, we consider a more realistic scenario, in which we want to fine-tune OMNI-P2x for the specific problem, i.e., the given combination of the desired QM level and system.

To perform the fine-tuning as data-efficient, as possible, we generate data to fine-tune the pre-trained universal potential OMNI-P2x by modifying our previously developed active learning procedure [66], where we use the OMNI-P2x as a starting point in each iteration of active learning loop. Performing direct active learning from scratch, i.e., generating data for a neural network potential with random weights, has been, so far, state of the art in terms of time- and data-efficiency for ML-driven NAMD, but it might still take several days to weeks of wall-clock time to converge, even for relatively simple systems, such as fulvene or azobenzene [66].

As our first test system, we use fulvene, popularly chosen as a prototypical molecule for benchmarking NAMD methods [48]. This choice also allows us to compare our new protocol, based on fine-tuning OMNI-P2x, to the previous state-of-the-art protocol using MS-ANI with random initial weights [66]. As the reference electronic-structure method, we use CASSCF (6,6). The NAMD is performed with the Landau–Zener–Belyaev–Lebedev (LZBL) variant of TSH [19, 56, 57, 115, 136], which does not require the calculation of the nonadiabatic coupling vectors (see Methods).

With active transfer learning employing fine-tuning of the universal potential OMNI-P2x in the active learning loop, we are able to slash the convergence time of fulvene from over three days, to about 24 hours, with a ca. three-fold increase in computational efficiency. Active transfer learning yielded only 3550 conformations, which is about twice fewer than the number of training points required to achieve convergence in active learning procedure performed from scratch, as can be seen in Figure 10a). Moreover, inspecting the time evolution of the electronic state populations, we can see a statistically significant difference in the populations predicted by the fine-tuned model and the previous MS-ANI model. Despite a smaller training set, the TL version of OMNI-P2x better captures the photophysical behavior of the system. A comparison of the dynamics predicted using different ML models and different LZBL propagation settings with reference fewest switches surface hopping dynamics is available in Suppl. Note 4. Further analysis of the NAMD trajectories of fulvene, including comparison of the ML-predicted potential energy surfaces, minimum energy conical intersections, as well as hopping points and deactivation channels is provided in Suppl. Note 4.1. In all cases agreement between the ML model and reference CASSCF calculations is noted.

This improvement is even more substantial when moving to more complex cases of active learning. One such example is the photoisomerization reaction of azobenzene, simulated at the AIQM1/MRCI-SD level of theory (see Methods). Due to a large volume of chemical space covered during such process, learning photoreactions is a particularly difficult task. To assess the improvement in data-efficiency that comes from using a pre-trained universal model, we compare the performance of the OMNI-P2x fine-tuned and MS-ANI trained from scratch on fractions of the full training set required to converge an AL procedure (35071 points, as reported in Ref. 66). Looking at the stability of the dynamics, calculated as the percentage of trajectories finishing without dissociation, shown in Figure 10b, the fine-tuned OMNI-P2x model can propagate stable dynamics given just 0.1% of the dataset, while dynamics propagated using the pure learning approach become consistently stable above 10% of the data. Furthermore, Figure 10c presents the predicted quantum yield of trans to cis photoisomerization, as a function of the training set size. For the TL potentials, all models trained on above 2% of the data predict QYs within the reference range, up to a 95% confidence interval error bar. On the other hand, the performance of the models trained from scratch is much less consistent, with outliers predicting incorrect quantum yields even at high fractions of the training set. An analysis of azobenzene S0 and S1 PES’, optimized MECIs, and deactivation geometries is provided in Suppl. Note 4.2, along with comparison to reference results, which agree with the ML models.

The above findings highlight that pre-training of OMNI-P2x on a vast dataset allows it to gain knowledge on both the chemical compound space and some information about conformational space, greatly reducing the amount of data needed for its fine-tuning. This makes OMNI-P2x a method of choice for propagating ML-accelerated surface hopping dynamics after data-efficient fine-tuning. A similar phenomenon has been reported in Refs. 114, 124, where a neural network potential pre-trained on ground-state data could achieve a lower MAE for excited-state predictions than a potential trained from scratch, in the low-data regime.

Methods
4.1 Model training and data collection

The entire PubChemQC dataset[88] was obtained from the PubChemQC Project(https://nakatamaho.riken.jp/pubchemqc.riken.jp/b3lyp_2017.html). The molecules in the PubChemQC dataset were optimized at the level of B3LYP[20]/6-31G*[1], and excited-state properties were calculated with TD-B3LYP/6-31+G*[1] using GAMESS[40]. For each molecule, first 11 electronic states were calculated, leading to 10 excitation energies and oscillator strengths.

In this study, the molecules were selected from the PubChemQC dataset which satisfied the following conditions: multiplicity of 1, net charge of 0, and TD-DFT data for 11 electronic states, and containing elements C, H, O, N, F, S, and Cl. The resulting data set contained 3.1 million molecules which we compiled with MLatom 3.16.2[62, 130, 136] to MLatom’s native machine-readable database format and saved in JSON files, used for training and testing in this study.

During training, the MSE loss function for energies-only was used:
L = (EML − Eref)^2, (1)
where EML and Eref are the ML-predicted and reference energies respectively.

Analogously, the MSE loss for oscillator strengths was used for a model trained on oscillator strengths.

4.2 UV/Vis spectra prediction

The electronic, UV/Vis absorption spectra are obtained using the standard single-point convolution approximation based on the broadening function as derived in Ref. 21:
L(x) = 0.619/Γ * sum(fi * e^-(x-xi)^2/Γ^2), (2)
where xi are the predicted transition energies, fi are their oscillator strengths, and Γ is a pre-defined full-width at half maximum parameter (FWHM). In this work, a broadening parameter of 0.3 eV was used, as recommended in Ref. 21.

Electronic spectra predicted by different methods were compared using Spearman’s rank correlation coefficient (SCC), r:
r = 1 − 6 * sum(di^2) / (n(n^2−1)), (3)
where di is the difference in rank between the two spectra, and n is the total length of the spectrum. The ODM2/CIS spectra were computed using the MLatom’s interface to the MNDO package [63], and TD-DFTB spectra were obtained using an interface to the DFTB+ program, release 24.1, [46], with the 3ob-3-1 set of parameters [93].

To address non-covalently bound structures separated by distances larger than the molecular descriptor cutoff, we introduce a fragment correction scheme (Supplementary Figure 1). The method checks for non-covalently bound fragments and compares excitation energies of the full system against those of individual fragments to determine whether the state corresponds to localized excitations or an excimer. Validation on a benzene dimer as a function of intermolecular separation shows that the corrected model agrees well with TD-DFT/B3LYP reference values, whereas the uncorrected model diverges when parts of the system fall outside the cutoff radius. Full details are provided in Supplementary Note 1.

For the nuclear ensemble approach spectra, 500 conformations of pyrazine were sampled from 298K Wigner distribution, as implemented in MLatom, and labeled using TD-DFT/B3LYP (6-31+G*) in Gaussian 16.0 (Revision C.01) software package [125]. Data for the 9-DCMA system was extracted from Ref. [60], where it was computed at the CAM-B3LYP [24]/ma-TZVP [127] level of theory.

The thermally distorted structures of molecules presented in Figures S6–S9 were sampled from a 2-ps-long ground-state Born–Oppenheimer molecular dynamics trajectory, with the first 1 ps discarded for thermalization, performed at the GFN2-xTB level of theory [42–44]. All the calculations and plots were prepared with MLatom.

4.3 NAMD simulations

CASSCF calculations used as a reference method for TSH simulations were performed through the interface to the COLUMBUS quantum chemistry package[122], using an active space of 6 electrons in 6 orbitals. TSH simulations were performed using the Landau–Zener–Belyaev–Lebedev formalism [19, 56, 57, 115, 136], with hopping probabilities between states k and j, Pj→k calculated using the following formula:
Pj→k = exp(-π * Zjk^2 / (2 * hbar * Z_dot_jk)), (4)
where Zjk is the energy gap between states j and k, and Z_dot_jk is its second-order time derivative.

The active learning loop uses the standard settings presented in Ref. 66, with initial points sampled from a harmonic approximation Wigner distribution, for a total of 250 points. The maximum propagation time was set to 60 fs with a time step of 0.1 fs. Velocities after hopping were rescaled in the direction of the momentum vector, using a reduced kinetic energy reservoir [134]. In each iteration of the AL procedure, 50 ML-TSH trajectories were run, and 50 gapMD trajectories. 15 additional points were sampled based on hopping probability uncertainty, up to a threshold of 300 points. Testing of the model was performed using 10000 trajectories for the MS-ANI model, as well as the fine-tuned OMNI-P2x model, and 623 CASSCF (6,6) trajectories were used as reference. One model from the energy ensemble of OMNI-P2x was used for fine-tuning in all NAMD applications.

The referenced data for azobenzene dynamics was labeled using the AIQM1/MRCI-SD method [55]. In the semiempirical part of AIQM1, the half-electron restricted open-shell Hartree–Fock formalism [34] was used in the SCF step, with the singly-occupied HOMO and LUMO orbitals. Two additional closed-shell references were added in the MRCI procedure, a HOMO–HOMO configuration and a doubly excited, LUMO–LUMO configuration. The active space consisted of 8 electrons in 10 orbitals (four occupied orbitals, six unoccupied). Single and double excitations within a such-defined active space were allowed. 100 trajectories were ran and analyzed for both the trans→cis and cis→trans photoisomerization reaction. The same number of trajectories was used when testing the model that was trained from scratch.

The simulated quantum yields Φ were calculated as the fraction of the trajectories relaxing to form the photoproduct (the other isomer), Nreactive, to the total number of trajectories Ntraj: Φ = Nreactive/Ntraj.

4.4 Azobenzene screening

All computations were performed with the XACS version of MLatom[89, 129, 136]. We used the MLatom interface to the TorchANI[105] package for all models based on ANI-type networks and their modifications (i.e., MS-ANI, OMNI-P2x, and AIQM3[5]). The AIQM3 calculations also used the MLatom interface to the s-dftd3 program [133] providing the D3(BJ) corrections [30] and the xtb program [135] providing the GFN2-xTB* baseline (GFN2-xTB [42–44] without the D4 dispersion corrections, e.g., as in the AIQM2 method [4]). The azobenzene derivatives were optimized with the AIQM3 version of the UAIQM models [107] targeting the CCSD(T)/CBS accuracy; this model is available as an add-on to MLatom provided via the Aitomic package[117]. For TD-DFT calculations, MLatom’s interface to the Gaussian 16 software package (Revision C.01) [125] was used. ADC(2)-level spectra were computed using the Turbomole [103] software package, with the Dunning correlation-consistent double-ζ basis set (cc-pVDZ) [37].

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript introduces OMNI-P2x as a universal neural network potential for molecular excited and ground state energy prediction. The model is inspired by the multi-state ANI model, based on own prior work by the authors [54], and as a novelty in this manuscript, adds all-in-one learning to simultaneously learn from multiple datasets with different levels of theory.\nThe manuscript is a step forward towards generalizable prediction of excited states energies, but there are a few questions that should be discussed before publication, in order to not only emphasize the strengths but also the limitations of this work.\nPage 6: \"The predictions are then made atom-wise and summed up to yield the total energy of the molecule.\" --> Does that mean that also the excited states energies are predicted by a sum of atom contributions? How are the labels for the excited states pre-processed (relative to the total energy of the ground state or as absolute energies), and how are they normalized?\nA simple test with your code and a duplicated and shifted benzene molecule shows that the excitation energies double when duplicating (and shifting) the input molecular structure. Of course, the model was not trained on such an input system, but the inductive bias given by the summation of single-atom contributions is not suitable for the prediction of excitation energies as non-extensive quantities (in contrast to ground state energies):\nimport mlatom as ml\nomni = ml.models.methods(method='omni-p2x')\nmol = ml.data.molecule()\nmol.read_from_xyz_file(\"Benzene.xyz\")\nomni.predict(molecule=mol, nstates=5)\nmol.excitation_energies\n0.18915822, 0.20142478, 0.2106019, 0.21771565\n\nmol2 = ml.data.molecule()\nmol2.read_from_xyz_file(\"Benzene2.xyz\") # two benzenes, the second is shifted by 10A in +z direction.\nomni.predict(molecule=mol2, nstates=5)\nmol2.excitation_energies:\n0.37831637, 0.40284953, 0.42120379, 0.4354313 # roughly twice the excitation energies of the single benzene molecule.\nFigure 2: How exactly is the state and the level of theory encoded? Figure 2 suggests that the state is an integer number appended to the molecular geometry vector, while the level of theory is represented as multiple one-hot encodings of length two (which seems not intuitive). Why is the state not also represented as a one-hot encoding? Do you have empirical evidence that this is the most informative / data-efficient representation?\nPage 7: \"PubChemQC project’s dataset consists of unique, optimized structures of molecules at the TD-B3LYP/6-31+G* level of theory, with the energies and oscillator strengths of the first 10 excited states calculated using TD-DFT.\" --> How does the model without transfer learning perform on excited states of equilibrium and non-equilibrium structures? Does the model without transfer learning just apply the learned ground state PES also to the excited states, so does it overestimate the correlation between ground and excited states, i.e. parallel ground- and excited states surfaces? Or does it somehow learn to predict excited states PES (at least close to equilibrium structures), even though it was only trained on excited states energies at equilibrium structures?\nFigure 4a shows total energies, but it is unclear if those are ground-state or excited-state total energies.\nFurthermore, calculating the r2 score on non-normalized total energies is misleading. Even just applying the atom-specific normalization to total energies without learning anything will lead to very high r2 scores. The r2 score should be calculated on the predictions on the normalized scale, not on the unnormalized scale.\nFigure 4 shows three different units. It would be better to discuss everything in one unit (e.g. eV), and then additionally convert key quantities to kcal/mol. But showing energies in Hartree in panel (a), MAEs in kcal/mol in (a) and (b), and then MAEs in eV in panel (c) is very confusing. Additionally, the difference between (b) and (c) does not become fully clear from the caption or text.\nExcitation energy predictions (relative to the ground state, separated by state) should be plotted in scatter plots, and r2 values should be calculated. Figure 3b indicates a standard deviation of ~1eV for excitation energies in the dataset, but I assume that the std is smaller if this is separated by state. Figure 4c indicates MAEs of ~0.2 eV, which is significant compared to the variance of the labels. Therefore, I assume that the r2 scores are not 0.9999, as calculated for the non-normalized ground state energies. A more open discussion about the actual accuracy (also in terms of r2 scores) is absolutely necessary to judge the reliability of the results and the impact of this study.\nFigure 6a/b fully focuses on equilibrium structures, which makes sense given the training data of the OMNI-P2x model, but which is not a fully fair comparison. Therefore, this should be explicitly mentioned in the caption and in the text/discussion. If possible, non-equilibrium structures of the same molecules should be sampled at a given temperature to see how fast the OMNI-P2x accuracy decays when going away from equilibrium geometries.\nFigure 6b suggests that ODM2/CIS and TD-DFTB are not on the Pareto front of speed and accuracy. However, it would be interesting to see how a delta learning approach that combines ODM2/CIS or TD-DFTB with OMNI-P2x would perform. The cost is roughly the same as the semi-empirical methods (because OMNI-P2x inference time is negligible), but the overall accuracy might be above the accuracy of OMNI-P2x.\nAll lines in Panel 8a seem more or less identical. It is not clear if the TL approach is statistically significantly better than the No-TL approach. A zoom into a particular part of the plot and the addition of uncertainty intervals is necessary to support the claims. Furthermore, the accuracy of TL and No-TP approaches in the population-time-curve should be quantified using an appropriate metric and plotted as a function of the active learning dataset size to obtain more insight into the advantages of TL here. The examples from 8c/d are not enough.\n(Remarks on code availability)\nI reviewed and tested the code. Everything works, the Jupyter notebook is reproducible. However, the library is not necessarily organised as you would expect from a professional software, which might not be necessary in this context.", "reviewer_2": "In their manuscript, Martyka et al. report the development of two artificial neural networks that, combined, predict approximate TD-DFT excitation energies and oscillator strengths of molecules. This is obtained by combining existing and established architectures by including data at CCSD(T)/CBS and TD-DFT B3LYP/6-31+G* levels of theory. The authors test their model on three applications, namely the convolution of absorption spectra based on predicted energies and oscillator strengths, the rapid screening of azobenzene derivatives to discover new candidates with desired properties, and the use of the model as a starting point to further fine-tune it with active learning to compute energies and gradients for nonadiabatic dynamics. While the manuscript introduces a potentially interesting architecture that contributes to increasing the potential of and reducing the limitations of artificial intelligence in the field of computational photochemistry, this work appears to have some theoretical limits and rather limited impact, breakthrough, and interest to justify publication in Nature Communications. Additionally, the authors claim the universality of the model, which would mislead the wide readership of this journal that might interpret the model as actually universal and flawless, which, in my opinion, it is far from.\nThere are three reasons why I disagree with the universality and the high impact of the model:\n- The model seems definitely transferable (with debatable prediction accuracy) among the chemical space similar to the one used for the training. However, this is limited to a certain type of molecule (only containing H, C, N, O, F, S, and Cl atoms), leaving out large categories of important photoactive molecules. For example, no transition metal complexes are considered, whose inclusion would make the model more universal.\n- There is no information about the character of the states, triplet states energies, nonadiabatic and spin-orbit coupling are not predicted, and there is no information about how to deal with ionic molecules. The model is able to predict vertical singlet excitations with a feature for the state order, but there is still no information about the orbitals involved and the actual character of the excited states. The universality is claimed to be in terms of how many different molecules can be predicted, but not in terms of photophysical properties and accuracy.\n- The training on the TD-DFT data makes for sure possible to predict energies and oscillator strength, but this is still training on wrong results in most cases, as the functional chosen is well known to have very poor reliability for the description of fundamental features of excited states. So a universal model should predict good and physically meaningful results, and not have as a goal to perform better than poorly reliable semi-empirical approaches.\nAdditionally, I found other major issues with the manuscript:\n1. The only task for which I would recommend a trustworthy use of the model is the virtual screening of azobenzene derivatives, where being able to qualitatively and rapidly screen molecules can help identify and select potential candidates with desirable properties. Sadly, the quantitative performance shown in the UV/Vis test case is very poor. Even in the (potentially cherry-picked) example from the top quartile, the prediction is, at least partially, wrong. The oscillator strengths are off, and between 275 and 300 nm there are 4 states predicted vs 1 at TD-DFT/BLYP level. The authors claim an excellent agreement, probably referring to the position of the maxima, but anyway, the convolution of the spectra is highly affected by the artificial broadening of 0.3 eV, which might make the results look more in agreement than they actually are (additionally, the reference to justify 0.3 it is unrelated, as in this manuscript the error is between TD-DFT data and their prediction, and not a discrepancy between computed/computed at different levels or computed/experimental spectra)\n2. The overall idea of the UV/Vis task is to predict the brightness of states in a specific order. This is not enough for this task; it is just a collection of energies with no correlation with their brightness. No information about the character of the states is present, and why a state is bright or not. The adiabatic order of the states (S_1, S_2, etc.) has no correlation with what states they actually are. Two molecules could have S_4 both bright for two completely different and uncorrelated reasons, while two very similar molecules could have S_1 and S_2 with inverted brightness just because, for example, npi* and pipi* states are inverted. So if there is no information on the mapping between excited states, predicting oscillator strength of adiabatic states is useless, if not wrong.\n3. In this logic, I also cannot understand, and it’s not explained, why the highest accuracy is achieved for the S0 excitations to S1 and to higher-lying S8–S10 states. I cannot see any correlation in that, as there is no information on the character of the state in the training. I believe there might just be a high density of states in some energy range that biases the training process.\n4. In the dynamics part, the impact of the new model is limited. There are no particular advances with respect to previous molecule-specific training for nonadiabatic dynamics. The model is not transferable and universal because it needs fine-tuning at a different level of theory, coming back to the problem of single-molecule training. It could proven a computational speed-up, but not a paradigm shift in ML application in NAMD. Additionally, I wonder why the work previously published in Nat. Comms. by Axelrod et al. is not cited, where an example of transferable ML architecture for NAMD is actually reported. (https://www.nature.com/articles/s41467-022-30999-w). Also more recent example of transferable NN for excited states should be mentioned (https://arxiv.org/abs/2502.12870).\n5. The fulvene dynamics look very different from the reference in literature, where a reflection in the S_1 population is present, particularly with trajectory surface hopping (Phys. Chem. Chem. Phys., 2020,22, 15183-15196). This is due to a physical reason, with the molecule encountering two different conical intersections, one of them populating back the first excited states. This is the reason why this model has been chosen to test nonadiabatic methods. Already the reference method used in this manuscript misses this effect. The training with the model actually fully deletes the small reflection at around 15 fs, removing completely the physical behavior expected. I believe this is more due to interpolation reasons than capturing the complexity of the potential energy surface, but in general is an example where the prediction erases the correct physical behavior and potentially leads to misinterpretation. This could be solved using another formulation of TSH, for which nonadiabatic coupling or wavefunction overlap would need to be predicted by the model, but is not available in the architecture presented in this work.\n6. More information about the computational cost of training would be necessary. Yes, the prediction is definitely cheaper than TD-DFTB, but working on a large dataset, what is the computational cost of the training with respect to calculating these excitations with the semi-empirical method?\nMinor points:\n- Is the label of Figure 4a correct? I think the authors show the ground state energies.\n- The statement “While these methods are certainly more affordable for systems with hundreds and even thousands of atoms, the accuracy of their predictions can often be questionable” needs a reference.\n- Some self-glorifying strong claims like “This marks the watershed moment that ML approach for excited-states presents a viable alternative to established QM methods.” should be toned down, as they just represent the author's personal opinions.\n(Remarks on code availability)", "reviewer_3": "This communication presents a novel foundation model / universal model for the excited state potential energy surfaces of organic molecules. The introduction motivates the model well in terms of deficit of the capabilities of existing computational methods for excited states in combination with high-throughput or dynamics-based studies of anything beyond very small molecular systems. The approach taken is to adapt a variant of the well-established ANI neural network potential, recently-introduced by the current authors and other collaborators in [54] npj Computational Materials 11, 132, which introduces a state-index as an extra input to the neural network alongside the geometry and atomic numbers, resulting in a “MS-ANI” approach. The current work builds on the demonstrations in [54] by showing that a foundation model, OMNI-P2X, can be built with this architecture which can be general across a fairly wide range of organic molecules. Overall the model is demonstrated to be quite successful and I think the paper could be quite influential as it represents a first of its kind. However, I feel there is room for improvement in the presentation and demonstrations before publication: in several sections the results themselves are very briefly described and it is hard to discern the level of success of the model.\nTraining data has been sourced from two existing databases, the PubChemQC database of TD-B3LYP/6-31+G* calculations, and the ANI-1ccx dataset of small molecules including off-equilibrium configurations, which has been provided as ground and excited-state DFT and TDDFT data for 10 excitations, and as ground-state-only CCSD(T), allowing the model to predict multiple levels of theory. Results show the S0 ground state has an MAE on energies of 2.63 kcal/mol which makes it comparable to other models trained on large databases of organic molecules, though I would prefer to see an explicit comparison to understand clearly how good the model is for the ground state.\nThe S1 state and higher-lying excitations are substantially worse, at ~5 kcal/mol, which is still useful, but it would be interesting to understand how correlated are the errors between S0 and the other states, ie should we expect errors on S0->Sx transition energies to exceed those of the contributing states, or will there be error cancellation to some degree? Scrutiny of Figs 4b and 4c might provide this information in an averaged sense, but they have different energy units so the comparison is not straightforward (generally, I find the broad selection of energy units (variously Ha, eV, and kcal/mol, sometimes in the same figure) rather troubling.\nFigure 3 shows the distribution of sizes and excitation energies of the training data, indicating that the majority of the PubChemQC data is in the 20-40 atom range whereas that of the ANI-1cxx data is much smaller, peaking at 12 atoms. Figure 4 shows some test results: I note that while figure 4a purports to demonstrate close agreement between predicted and target total energy, the spread of energies is so wide (1000 Hartrees!) that presenting this as a y=x figure is highly redundant – any deviation from y=x visible at this scale would represent an enormous error, and the MAE of 5.78kcal/mol ~ 0.01 Hartree would be invisible on this scale. Perhaps a baseline subtraction of the atomic energies from each case would make this graph more useful?\nRegarding the fact that 10 excited states are calculated: are this many states always going to give meaningful results in these systems with a 6-31+G* basis: in small molecules these might be higher-lying Rydberg states which are not well-represented in a relatively small basis such as 6-31+G*?\nIt is interesting that the Energy MAE on S1 is roughly double that of S0, then further higher-lying excitations have increasingly lower errors. It strikes me that perhaps the increasing Rydberg character of higher-lying states is somehow easier for the model to capture. The authors say “An outlier here is S0, which has the smallest energy MAE of 2.63 kcal/mol, which we attribute to the effect of the added ground-state CCSD(T) data” – I couldn’t see how that could come about – why would the presence of data associated with a second level of theory predictable by the same model improve the accuracy of its predictions for the first?\nSome minor points on figures: “One can note, that the highest accuracy is achieved for the S0 excitations to S1 and to higher-lying S8–S10 states.” Based on Figure 4d this does not appear to be true – do they mean the lowest accuracy? In the stick plot in Figure 5: the baseline of sticks does not align with zero on right hand axis – this needs correcting. Panel 8a is very separate from 8b,8c: they refer to different systems and unrelated calculations.\nI would say that the Statement “OMNI-P2x is approaching the accuracy of time-dependent density functional theory” in the abstract and main text is potentially misleading, despite being carefully written. As Figures 5 and 6 show, while the model does significantly better than semi-empirical methods, at reproducing TDDFT spectra, as measured by the Spearman Correlation Coefficient, the reasons for this in terms of agreement of the underlying stick spectrum with TDDFT are often dubious. In many cases, where TDDFT has a few strong peaks, OMNI-P2X has a forest of intermediate peaks in a similar energy range, with no clear relationship between oscillator strengths. I would want to know what the Mean Relative Error on oscillator strengths was, as well as the MAE, because it seems like in many cases it might be quite large.\nFor a method to be able “approach the accuracy of TDDFT” I would expect individual excitations to be predicted to better than the 5kcal/mol quoted in Figure 4b as the MAE on S1 energies, and for the mean relative error on oscillator strengths not to be too extreme. This current level of energy error corresponds to 0.22eV or 100nm on a S0<-S1 excitation around 1.5eV, which is not really approaching TDDFT accuracy.\nThe section on design and screening of azobenzene compounds provides a proof-of-principle for use in screening large datasets of candidate molecules. In this case azobenzene derivatives are screened to find 21 leads with strongly red-shifted absorption, at least some of which stood up to further scrutiny via calculations at a higher level of theory, ADC(2). I found it frustrating that this section ended without serious investigation of whether OMNI-P2X had really found strong leads: no comparison is made of the spectra from TD-DFT, from ADC(2) or from OMNI-P2X. I would want to see this investigated more deeply.\nThe final section demonstrates what could turn out to be a common use-case for foundation models for the excited state: active-learning-driven fine-tuning / transfer-learning is used to train a OMNI-P2X model at the CASSCF(6,6) level of theory for fulvene. The transfer-learning is shown to substantially accelerate model training, and both the TL and no-TL versions accurate reproduce dynamics of the reference CASSCF calculations. Very brief mention is made of calculations on photoisomerisation of azobenzene, but the results presented are rather cryptic: I was not able to intuit what point was being made about dihedral angles as obtained from the dynamics.\nOne natural question would be how the model responds to large molecular systems and relatively large numbers of states required, and is there consistency of the predictions with respect to the number of states requested. In some testing of the version made available on figshare, I observed that when calling OMNI-PX2 three times, for 2, 5 and 10 states, for a model of a sunscreen candidate molecule with 95 atoms, in each case I got a very different set of excited states predicted (the S1->S0 transition was much lower in the 10 state case than the 2 state case). Is there anything that guarantees consistency with respect to the number of states requested. I note it is relatively stable for the provided example of benzene, though the oscillator strength varies a little.\nOverall, the paper presents a novel methodology and opens up new avenues for investigation, in that it presents the first successful universal excited state model. However, I feel many of the tests and examples need to be refined further and presented more clearly before the community can properly judge the success of the model.\n(Remarks on code availability)\nI downloaded the shared version from figshare, verified that the benzene model reproduced the data provided, and verified the availability of the model's training data. I am not able to review the 25k lines of code associated with MLAtom itself, but the reproducibility was quite straightfoward and MLatom worked out of the box in via torchani.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\n(Remarks on code availability)\nI sat with my supervisor and tested and discussed the code. I did not read the entirety of the codebase."}, "author_response": "Responses to the comments of the\nreviewers\nReviewer comments are in italics, and our responses are in regular, blue font. Our\nchanges to the manuscript are marked in red font.\nReviewer #1\nThe manuscript introduces OMNI-P2x as a universal neural network potential for\nmolecular excited and ground state energy prediction. The model is inspired by the\nmulti-state ANI model, based on own prior work by the authors [54], and as a novelty\nin this manuscript, adds all-in-one learning to simultaneously learn from multiple\ndatasets with different levels of theory.\nThe manuscript is a step forward towards generalizable prediction of excited states\nenergies, but there are a few questions that should be discussed before publication, in\norder to not only emphasize the strengths but also the limitations of this work.\nWe thank you for your positive feedback. In the following revision, we have put\ngreat effort into increasing the clarity of the presented results, showing\napplications that better showcase the power of the OMNI-P2x architecture,\nespecially in NAMD simulations, as well as clearly stating the limitations of the\npresented work. Furthermore, further computational effort allowed us to double\nthe training set size of the ML model and replace the NN potential with an\nensemble of three potentials. This allows not only for more accurate and stable\npredictions of ground and excited state energies, but also to quantify the\nuncertainty of the NN predictions.\nFor the energy part of the predictions, an ensemble of three NN\npotentials is used, with their predictions averaged to obtain the final\npredicted energy for all electronic states.\nWe have used the entirety of the PubChemQC dataset, which meets our\nfiltering criteria (as detailed in the Methods section), yielding 3.1 million\nmolecules with 11 electronic states each.\nAdditionally, OMNI-P2x provides an uncertainty quantification (UQ) of\nall predicted energies, using the standard deviations between quantities\npredicted by each of the three NN potentials as the error bar. A warning\nis issued to the user if the UQ of any excited states exceeds a safe\nthreshold, which is determined as the median plus three median absolute\ndeviations of the UQs of the test set, for each electronic state.\nReviewer 1, Comment 1\nPage 6: \"The predictions are then made atom-wise and summed up to yield the total\nenergy of the molecule.\" --> Does that mean that also the excited states energies are\npredicted by a sum of atom contributions? How are the labels for the excited states pre-\nprocessed (relative to the total energy of the ground state or as absolute energies), and\nhow are they normalized?\nYes, excited state energies are predicted in the same way as ground-state\nenergies, which is a sum of atom-wise contributions. The energy labels of all\nstates are pre-processed together, by means of computing one set of self-atomic\nenergies in the ANI model that is used for all states and then subtracted. This\nensures that the energy gaps between electronic states are conserved in the raw\nANI energies, which would not be the case if self-atomic energies were\ncomputed for each state separately, as tested in our previous work developing\nthe multi-state ANI model.\nReviewer 1, Comment 2\nA simple test with your code and a duplicated and shifted benzene molecule shows that\nthe excitation energies double when duplicating (and shifting) the input molecular\nstructure. Of course, the model was not trained on such an input system, but the\ninductive bias given by the summation of single-atom contributions is not suitable for\nthe prediction of excitation energies as non-extensive quantities (in contrast to ground\nstate energies):\nimport mlatom as ml\nomni = ml.models.methods(method='omni-p2x')\nmol = ml.data.molecule()\nmol.read_from_xyz_file(\"Benzene.xyz\")\nomni.predict(molecule=mol, nstates=5)\nmol.excitation_energies\n0.18915822, 0.20142478, 0.2106019, 0.21771565\n\nmol2 = ml.data.molecule()\nmol2.read_from_xyz_file(\"Benzene2.xyz\") # two benzenes, the second is shifted by 10A\nin +z direction.\nomni.predict(molecule=mol2, nstates=5)\nmol2.excitation_energies:\n0.37831637, 0.40284953, 0.42120379, 0.4354313 # roughly twice the excitation\nenergies of the single benzene molecule.\nWe thank you for pointing out this potential problem and providing a code\nsnippet allowing for its reproduction. OMNI-P2x uses the ANI framework for\npredicting electronic state energies, which generates atom-wise atomic\nenvironment vectors (AEVs) for each atom. These AEVs include information\nabout neighboring atoms up to some cutoff (set to 4.0 A in this work). If parts\nof a system are completely separated by a distance larger than this cutoff, like\nin the benzene dimer example, such pathological behavior can be expected, as\nthe descriptors of each monomer have no information about the rest of the\nsystem. To simulate this behavior, we have computed the S0 → S1 vertical\nexcitation energy for a benzene dimer as a function of the monomer distance,\nR, which is presented below.\nWe can note that while the ML-model initially gives predictions that are\nreasonably close to the DFT results (considering the fact that dimeric systems\nwere not explicitly included in the training set), for distances between 3 and 4 Å\nthe excitation energy begins to rise, as parts of the system move out of the cutoff\nradius of the rest, eventually reaching about double the excitation energy of TD-\nDFT. From this, we conclude that while OMNI-P2x cannot correctly describe\nexcitation energies for systems with large spatial separation due to the intrinsic\nlimitations of the AEV descriptor, it can implicitly learn the correct, intensive\ncharacter of the excitation energy, within the descriptor cutoff.\nIn the limit of infinite separation of chromophores, the OMNI-P2x’s solution\nwill indeed collapse into the single excitations on each of the chromophores.\nThis is unphysical, but it is also an issue with some of the first-principles\ncalculations, i.e., the calculations of a multi-chromophore system in periodic\nboundary conditions, also require deciding which chromophore to treat as the\nexcited one, while others are treated as non-excited chromophores.\nImportantly, the NN architecture itself can learn the correct behavior for\nseparated systems, provided suitable training data. To demonstrate this, we fine-\ntuned OMNI-P2x on a set of 600 benzene dimer conformations sampled via\nmetadynamics, with the results shown below. The agreement between the ML\nprediction and reference is excellent.\nDespite this capability, we acknowledge that the absence of such systems in the\noriginal training data, as well as the inability of the model to infer correct long-\nrange excitation behavior from monomer data alone, results in OMNI-P2x\nproducing unphysical predictions for widely separated chromophores.\nHence, to correct the long-range unphysical behavior, we propose a mitigation\nstrategy, where we automatically decide which chromophore to treat as excited,\nwhile the others are treated in their ground state. Before predicting for a given\nsystem, a check is run to detect non-covalently bound fragments of the system.\nThen, predictions are made separately for each fragment and the entire system.\nIf the excitation energy of the total system is greater than the minimum\nexcitation energy of the subsystems, that means that OMNI-P2x most likely\nincorrectly treated the system as multiply excited, with an excitation localized\non each chromophore. Hence, calculations are performed separately for each\nsub-system and sorted to yield the final spectrum. If, however, the excitation\nenergy of the total system is smaller than the minimal energy for each sub-\nsystem, we identify it as an excimer and keep this result. A visual representation\nof the fragment-correction algorithm is given below:\nWe compare the results obtained with the fragment-corrected approach to\npristine OMNI-P2x and TD-DFT below.\nThe agreement between TD-DFT and the fragment-corrected scheme becomes\nexcellent for distances above the local descriptor cutoff. It should also be noted\nthat DFT is not the most suitable level of theory for describing systems with\nlarge spatial separation, which is also evidenced here by an unphysical bump in\nthe excitation energy between 8 and 9 A.\nThe revised manuscript now reads, with the above discussion recounted in\nsection S1 of the ESI:\nA fragment-correction scheme is applied to non-covalently bound structures, to\nrecover the correct, intensive character of excitation energies, for molecular\nfragments separated by distances greater than the descriptor cutoff radius, with\ndetails provided in Section S1 of the ESI.\nReviewer 1, Comment 3\nFigure 2: How exactly is the state and the level of theory encoded? Figure 2 suggests\nthat the state is an integer number appended to the molecular geometry vector, while\nthe level of theory is represented as multiple one-hot encodings of length two (which\nseems not intuitive). Why is the state not also represented as a one-hot encoding? Do\nyou have empirical evidence that this is the most informative / data-efficient\nrepresentation?\nYour interpretation of Figure 2 is correct, the state ordering number is indeed an\ninteger that is concatenated with the atomic environment vector (AEV), while\nthe level of theory is a one-hot encoded vector. We use non-negative integers to\nencode the state, as the adiabatic state energies are in ascending order, the same\nas the descriptor used. On the other hand, there is no clear ordering between the\nenergies predicted at different levels of theory; hence, we choose a one-hot\nencoding. To have a direct comparison, we have trained two separate NN\nmodels on the fulvene dataset presented in this manuscript: one using integer\nencoding for states, and the other using one-hot encoding. The integer-encoded\nmodel outperformed the one-hot-encoded model, achieving a MAE of 0.0279\nand 0.02590 for S0 and S1 energies, while the one-hot-encoded model had\nMAEs of 0.0330 eV and 0.0326 eV, respectively. We take this benchmark, along\nwith the physically motivated ordering of the descriptor, as evidence to use it in\nour foundational model.\nThe manuscript now reads:\nThe electronic state is encoded as a non-negative integer, consistent with\nthe natural ordering of adiabatic state energies, which we find to be both\nphysically motivated and empirically more accurate than one-hot\nencoding.\nReviewer 1, Comment 4\nPage 7: \"PubChemQC project’s dataset consists of unique, optimized structures of\nmolecules at the TD-B3LYP/6-31+G* level of theory, with the energies and oscillator\nstrengths of the first 10 excited states calculated using TD-DFT.\" --> How does the\nmodel without transfer learning perform on excited states of equilibrium and non-\nequilibrium structures? Does the model without transfer learning just apply the learned\nground state PES also to the excited states, so does it overestimate the correlation\nbetween ground and excited states, i.e. parallel ground- and excited states surfaces?\nOr does it somehow learn to predict excited states PES (at least close to equilibrium\nstructures), even though it was only trained on excited states energies at equilibrium\nstructures?\nThank you for yet another interesting question. To check this, we have\nperformed two relaxed potential energy surface scans, with results presented\nbelow. Panel a) shows the scan along C=O bond stretching in cyclohexanone,\nwhile panel b) shows a torsional scan around the C-N=N-C dihedral angle in\nazobenzene. The geometries were optimized using DFT/B3LYP(6-31+G*),\nand\nresults obtained using TD-DFT/B3LYP and OMNI-P2x (DFT-level predictions)\nare compared for the first three electronic states. For the cyclohexanone\nexample, as expected, in the vicinity of the geometric minimum (1.2 A), the\nmodel is very accurate. This accuracy begins to lower as we move outside of\nthe minimum, although for the stretched C=O bond the accuracy is still rather\nreasonable, without any catastrophic extrapolation failures. Shortening of the\nC=O bond however leads to much larger errors, as such high-energy structures\nare far from the data contained in the training set. Taking a look at the more\ncomplicated example of azobenzene torsion, around the two minima (0 and 180\ndegrees), OMNI-P2x achieves near perfect accuracy for the ground state, and\nan error of about 0.25 eV for the S1 and S2 states. This accuracy is retained up\nto around 25 degrees outside of the minimum, where it starts to rapidly decay\ndue to the change in character of S0 and S1 associated with twisting of the N=N\ndouble bond, which cannot be captured by OMNI-P2x, due to the lack of\ncorresponding structures in the training set. From this test, we conclude that\ndespite training on just equilibrium structures, OMNI-P2x can extrapolate the\nshape of the excited-state PES’s at points close to the geometric minimum;\nhowever, when this geometric change is associated with a substantial change in\ncharacter of the excited state (such as in the azobenzene case), this extrapolation\nfails. The ground-state and excited-state PES are not parallel, although, in line\nwith your suggestion, their correlation is certainly overestimated when moving\nto out-of-equilibrium structures.\nThe main text has been modified to read:\nTo further probe this behavior, we performed relaxed PES scans of\ncyclohexanone (C=O stretch) and azobenzene (C–N=N–C torsion), shown in\nFigure S10. These scans confirm that OMNI-P2x is accurate near equilibrium\nand along moderate distortions, but errors grow for strongly compressed\ngeometries or when torsion induces a change in excited-state character.\nTogether, these results indicate that OMNI-P2x can extrapolate PES shapes\nclose to equilibrium, but its reliability decreases as the molecular distortion\ndrives the system beyond the domain of the training data.\nAt the same time, we show that the capabilities of the model can be easily\nextended to non-equilibrium structures by fine-tuning, even in the ultra-low data\nregime. Another example on which this is demonstrated, aside from NAMD, is\nby simulating nuclear ensemble approach spectra, in a new section “Nuclear\nensemble approach spectra”. By transfer learning with just 25 datapoints, the\nmodel can perfectly recover the dynamically averaged spectrum of pyrazine,\nand with just 50 points can produce accurate results for the much more\ncomplicated example of 9-DCMA. Noteworthy, this example contains 31\nelectronic states, while OMNI-P2x was pre-trained on only 11. Despite that,\nthe\naccuracy of the higher excited states is acceptable (with near-perfect\nreproduction of the lower energy region of the spectrum), having correct\nplacement of the peaks and their relative intensity. This hints that the model can\nalso extrapolate in terms of the number of electronic states requested.\nFurthermore, these results show unprecedented data-efficiency, as previous\napproaches required at least 100 to 200 points to produce acceptable results. The\nmain text now contains a new section, as follows:\nA more accurate method for predicting UV-VIS absorption spectra is the\nNuclear Ensemble Approach (NEA). In this technique, vertical excitation\nenergies and oscillator strengths are computed for an ensemble of\nhundreds or thousands of molecular geometries. The final spectrum is\nconstructed by summing Lorentzian or Gaussian functions centered at\neach geometry’s excitation energy, weighted by the corresponding\noscillator strength. Compared to single-point convolution, this approach\nbetter captures vibronic effects. However, NEA simulations are\ncomputationally demanding, as they require quantum mechanical (QM)\ncalculations for all structures in the ensemble. To address this cost,\nmachine learning (ML) models offer a natural path for acceleration.[78,\n79] Since OMNI-P2x is initially trained only on equilibrium structures, it\nis not directly applicable to NEA simulations. Nevertheless, fine-tuning\nenables it to achieve remarkable data efficiency and simulation speed in\nthis context.\nTo demonstrate this, we consider two test cases: pyrazine (with 4 excited\nstates) and 9-Dicyanomethylene-Acridine (9-DCMA, with 30 excited\nstates), the latter of which was previously used to benchmark ML-NEA\nspectra with kernel ridge regression models.71 The NEA spectra simulated\nusing reference TD-DFT data, as well as those predicted by the fine-tuned\nOMNI-P2x model—trained on just 25 data points for pyrazine and 50 for\n9-DCMA—are shown in Figure 7. For pyrazine (panel a), the agreement\nis excellent across the entire spectrum. In the more complex case of 9-\nDCMA, the model accurately captures the intensity and position of the\nlower-lying excitations. Deviations start emerging below 250 nm,\nespecially in predicted intensities, while the position and relative\nintensities of the peaks remains rather accurate. We attribute this to the\nmodel extrapolating in a regime beyond the number of excited states\nincluded in training dataset: OMNI-P2x was pre-trained on 10 excited\nstates, while 30 are considered in this example. Such simulations can be\ndone on commodity hardware (Jupyter notebook using 4 CPUs), in under\na minute for pyrazine, and about 40 minutes for 9-DCMA, excluding the\nlabeling time of QM methods.\nReviewer 1, Comment 5\nFigure 4a shows total energies, but it is unclear if those are ground-state or excited-\nstate total energies.\nFurthermore, calculating the r2 score on non-normalized total energies is misleading.\nEven just applying the atom-specific normalization to total energies without learning\nanything will lead to very high r2 scores. The r2 score should be calculated on the\npredictions on the normalized scale, not on the unnormalized scale.\nFigure 4 shows three different units. It would be better to discuss everything in one unit\n(e.g. eV), and then additionally convert key quantities to kcal/mol. But showing energies\nin Hartree in panel (a), MAEs in kcal/mol in (a) and (b), and then MAEs in eV in panel\n(c) is very confusing. Additionally, the difference between (b) and (c) does not become\nfully clear from the caption or text.\nExcitation energy predictions (relative to the ground state, separated by state) should\nbe plotted in scatter plots, and r2 values should be calculated. Figure 3b indicates a\nstandard deviation of ~1eV for excitation energies in the dataset, but I assume that the\nstd is smaller if this is separated by state. Figure 4c indicates MAEs of ~0.2 eV, which\nis significant compared to the variance of the labels. Therefore, I assume that the r2\nscores are not 0.9999, as calculated for the non-normalized ground state energies. A\nmore open discussion about the actual accuracy (also in terms of r2 scores) is\nabsolutely necessary to judge the reliability of the results and the impact of this study.\nWe agree with the suggestions concerning the clarity of the results presented in\nFigure 4. In line with this comment, we have made the following changes to\nFigure 4:\n1) Figure 4a has been changed to compare normalized energies (with self-\natomic energies removed), with the R2 calculated.\n2) All units have been converted to kcal/mol.\n3) Scatter plots of predicted vs true excitation energies for each state have\nbeen presented in Figure S4 of the ESI, with calculated R2 scores and\nMAE. In all cases the MAE of the predicted excitation energies is\nsmaller than the reported lower bound of the accuracy of TD-DFT\n(0.2 eV ≈ 4.6 kcal/mol), making the error related to fitting the ML\nmodel lower than the intrinsic accuracy of the method.\nThe main text has been modified as follows:\nWe evaluate the performance of the model using a test set of 300k unseen molecules\nfrom the PubChemQC dataset for the excitation energies, as well as a test set of almost\n49k ground-state conformations from the ANI dataset. The MAE in energy predictions\nfor the first 11 electronic states with respect to the TD-DFT data shows that the quality\nof the predictions increases when going to higher excited states, reaching a minimal\nvalue for S7 (Figure 4b and c). An outlier here is S0, which has the smallest energy\nMAE of 1.6 kcal/mol, which we attribute to the improved data-efficiency due to the\napplication of the all-in-one learning to the combined set with the ground-state DFT\n(from PubChem) and the CCSD(T) (from ANI-1ccx) data, and an overall smoother and\nlower-variance distribution of ground-state energies that makes them easier for the\nmodel to learn. Furthermore, the 8 MAEs for both the predicted total and excitation\nenergies are below 0.2 eV (4.6 kcal/mol), which lies below the lower bound of the\nreported accuracy for the reference TD-DFT method, indicating that the\nfitting error of the ML model is smaller than the intrinsic error of the\nunderlying quantum chemical approach. Correlation plots showing the\npredicted vs true excitation energies for all transitions are available as\nFigure S4 of the ESI. In all cases, the R2 is above 0.91. For\nthe ground-state energies in the test ANI-1ccx subset, OMNI-P2x has a higher MAE of\n4.14 kcal/mol, but still excellent R2 = 0.9957 for the normalized energies (0.9999 for\nunnormalized).\nReviewer 1, Comment 6\nFigure 6a/b fully focuses on equilibrium structures, which makes sense given the\ntraining data of the OMNI-P2x model, but which is not a fully fair comparison.\nTherefore, this should be explicitly mentioned in the caption and in the text/discussion.\nIf possible, non-equilibrium structures of the same molecules should be sampled at a\ngiven temperature to see how fast the OMNI-P2x accuracy decays when going away\nfrom equilibrium geometries.\nWe revise the manuscript by changing the caption of Figure 6 to explicitly\nmention that the evaluation was done on equilibrium structures only:\nThe evaluation of the quality of the UV/Vis absorption spectra predicted\nby OMNI-P2x model on a test set of 50k molecules, at their ground-state\nequilibrium structures.\nSampling a representative number of out-of-equilibrium conformations that include\nfinite-temperature effects for 50k molecules would require on the order of millions of\nadditional reference calculations, which is unfortunately not computationally feasible\nwithin the scope of this work. However, to address this important point, we extend the\nmanuscript with two complementary investigations. Firstly, the newly added “Nuclear\nensemble approach spectra” section shows a workflow where OMNI-P2x can be used\nto predict spectra with explicit inclusion of finite-temperature effects, at a fraction of\nthe cost, which has been discussed in detail in response to your Comment 4.\nFurthermore, the performance of OMNI-P2x when moving away from equilibrium\nstructures has also been quantified by means of PES scans, also described in the\nresponse to Comment 4.\nSecondly, to directly assess the robustness of the base model when moving away\nfrom equilibrium geometries, we sampled 200 thermally distorted\nconformations of the molecules presented in Figure 5 and evaluated excitation\nenergy predictions. The violin plots below (available as Figures S5–S8) show\nthat while errors for higher excited-states remain close to equilibrium-structure\nerrors, lower-lying states (S1–S3), are much more sensitive to thermal\nfluctuations, and thus, their errors can be much higher, with an overall 2 to 4\ntimes lower accuracy of the prediction.\nThose results are referenced in the manuscript as:\nFor the test molecules shown in Figure 5, we tested the effect of finite-\ntemperature distortions on the accuracy of OMNI-P2x and found that\nhigher excited states remain close to equilibrium benchmarks, while the\naccuracy of the lowest excitations (S –S ) may decrease more sharply.\n1 3\nViolin plots showing the error distribution on non-equilibrium structures\nare available in Figures S6–S9 of the ESI.\nAnd in the Methods section:\nThe thermally distorted structures of molecules presented in Figures S6–\nS9 were sampled from a 2-ps-long ground-state Born–Oppenheimer\nmolecular dynamics trajectory, with the first 1 ps discarded for\nthermalization, performed at the GFN2-xTB level of theory.\nReviewer 1, Comment 7\nFigure 6b suggests that ODM2/CIS and TD-DFTB are not on the Pareto front of speed\nand accuracy. However, it would be interesting to see how a delta learning approach\nthat combines ODM2/CIS or TD-DFTB with OMNI-P2x would perform. The cost is\nroughly the same as the semi-empirical methods (because OMNI-P2x inference time is\nnegligible), but the overall accuracy might be above the accuracy of OMNI-P2x.\nIndeed, this is an interesting question. At the moment of preparing this\nmanuscript, we have been working on delta-learning models for excited states\nfor a few months already. We observe that in certain cases delta-learning for\nexcited states can outperform pure ML, but, unfortunately, it requires more\ncareful considerations than delta-learning applied to the ground state problems,\nfor which one of us has extensive expertise. In general, delta-learning for\nexcited states brings about many additional questions and design choices, and\nrequires a dedicated, in-depth study, which we are going to report separately.\nReviewer 1, Comment 8\nAll lines in Panel 8a seem more or less identical. It is not clear if the TL approach is\nstatistically significantly better than the No-TL approach. A zoom into a particular part\nof the plot and the addition of uncertainty intervals is necessary to support the\nclaims. Furthermore, the accuracy of TL and No-TP approaches in the population-time-curve\nshould be quantified using an appropriate metric and plotted as a function of the\nactive learning dataset size to obtain more insight into the advantages of TL here. The\nexamples from 8c/d are not enough.\nWe realize that our original description and figures were not clear enough. We\nre-worked the fulvene example, running 10,000 trajectories with both ML\nmodels (TL and no-TL). The populations obtained using both the models, and\nreference CASSCF results have been plotted in Figure 9, along with 95%\nconfidence intervals shown with shaded areas. An inset zooming in on the\ncritical region between 10 and 20 fs has been added to the plot. As can be seen\nthere, the populations obtained using transfer learning and direct ML models\nexhibit statistically significant differences in the vicinity of the conical\nintersection, and it is the TL model that stays within the 95% confidence\nuncertainty interval of the reference population, while the model trained in an\nAL loop starting from scratch produces results that differ from the reference.\nThe revised manuscript related to the fulvene dynamics now reads:\nBy employing transfer learning from a universal potential, we are able to\nslash the convergence time of fulvene from over three days, to about 24\nhours, with a ca. three-fold increase in computational efficiency.\nAdditionally, the final employed training set contains only 3550\nconformations, which is about twice less as was necessary to achieve\nconvergence in a reference direct active learning procedure, as can be seen\nin Figure 9. Moreover, inspecting the time evolution of the electronic state\npopulations, we can see a statistically significant difference in the\npopulations predicted by the fine-tuned and the previous MS-ANI models.\nDespite a smaller training set, the TL version of OMNI-P2x better captures\nthe photophysical behavior of the system. Additionally, a comparison of\nthe dynamics predicted using different ML models and different LZBL\npropagation settings with reference fewest switches surface hopping\ndynamics is available in Section S4 of the ESI.\nTo increase the clarity of the presented results, numerical experiments\nsimulating the photoisomerisation reaction of azobenzene have been moved to\na separate Figure 10. We trained models using fractions of the full active\nlearning dataset from ref [57] (consisting of ca. 35,000 points), ranging from\n0.1% to 100%. For each fraction, a model was trained either using transfer\nlearning (TL), or from scratch, and used to propagate NAMD starting from trans\nisomer of azobenzene. We then quantify the accuracy of the predicted dynamics,\nby calculating the photoisomerisation quantum yields, as well as the stability of\nthe dynamics (measured as the fraction of the trajectories that remain\nundissociated until the end of the NAMD trajectory). The obtained results show\nthat the transfer learning model can propagate stable dynamics even with as\nlittle as 0.1% of the training points, while the model trained from scratch\nrequires at least 10% of the training set to achieve consistently stable results.\nFurthermore, looking at the quantum yields of the photoreaction predicted by\nthe two types of ML models, all of the TL models trained on a fraction greater\nthan 2% of the training set predict quantum yields that are within a 95%\nconfidence interval of the reference method (with the respective error bar),\nwhile the results obtained using ML models trained from scratch can deviate\nfrom the reference method, depending on the fraction of points use, indicating\nless stable results.\nThe revised manuscript reads:\nTo assess the improvement in data-efficiency that comes from using a pre-\ntrained universal model, we compare the performance of the OMNI-P2x\nfine-tuned and MS-ANI trained from scratch on fractions of the full\ntraining set required to converge an AL procedure (35071 points, as\nreported in Ref. 57). Looking at the stability of the dynamics, calculated\nas the percentage of trajectories finishing without dissociation, shown in\nFigure 10a), the fine-tuned OMNI-P2x model can propagate stable\ndynamics given just 0.1% of the dataset, while dynamics propagated using\nthe pure learning approach become consistently stable above 10% of the\ndata. Furthermore, Figure 10b presents the predicted quantum yield of\ntrans to cis photoisomerization, as a function of the training set size. For\nthe fine-tuned OMNI-P2x potentials, all models trained on above 2% of\nthe data predict QYs within the reference range, up to a 95% confidence\ninterval error bar. On the other hand, the performance of the models\ntrained from scratch is much less consistent, with outliers predicting\nincorrect quantum yields even at high fractions of the training set.\nReviewer #2\nIn their manuscript, Martyka et al. report the development of two artificial neural\nnetworks that, combined, predict approximate TD-DFT excitation energies and\noscillator strengths of molecules. This is obtained by combining existing and\nestablished architectures by including data at CCSD(T)/CBS and TD-DFT B3LYP/6-\n31+G* levels of theory. The authors test their model on three applications, namely the\nconvolution of absorption spectra based on predicted energies and oscillator strengths,\nthe rapid screening of azobenzene derivatives to discover new candidates with desired\nproperties, and the use of the model as a starting point to further fine-tune it with active\nlearning to compute energies and gradients for nonadiabatic dynamics. While the\nmanuscript introduces a potentially interesting architecture that contributes to\nincreasing the potential of and reducing the limitations of artificial intelligence in the\nfield of computational photochemistry, this work appears to have some theoretical\nlimits and rather limited impact, breakthrough, and interest to justify publication in\nNature Communications.\nAdditionally, the authors claim the universality of the model, which would mislead the\nwide readership of this journal that might interpret the model as actually universal and\nflawless, which, in my opinion, it is far from.\nWe respect your opinion and constructive comments below, but with politely\ndisagree with the claims of ‘limited impact, breakthrough’ – our models are\nalready used in many contexts, not just in our groups, but have also attracted\ninterest from other researchers. Reviewer 3 evaluated our work as follows:\n“Overall the model is demonstrated to be quite successful and I think the paper\ncould be quite influential as it represents a first of its kind”. Our models are\ncertainly not ‘flawless’, and we have emphasized that in the revised version of\nthe manuscript. It should also be stressed that hardly any novel tool aiming to\nprovide breakthrough functionality is free from flaws and limitations. Those\nlimitations will likely be gradually removed by future work that aims to improve\nthe current cutting-edge functionality.\nThe term ‘universal potential’ is commonly applied in literature in the context\nof potentials pre-trained on big data sets and generalizing to unseen compounds\nand it is general knowledge that they are not expected to be flawless. In this\nrespect our models are universal as they generalize to unseen compounds, albeit\nnot flawlessly.\nReviewer 2, Comment 1\nThere are three reasons why I disagree with the universality and the high impact\nof the model:\n- The model seems definitely transferable (with debatable prediction accuracy)\namong the chemical space similar to the one used for the training. However,\nthis is limited to a certain type of molecule (only containing H, C, N, O, F, S,\nand Cl atoms), leaving out large categories of important photoactive molecules.\nFor example, no transition metal complexes are considered, whose inclusion\nwould make the model more universal.\nWe agree that inclusion of more elements in the training set would make the\nmodel more universal. However, many (if not most) machine learned\ninteratomic potentials (MLIPs) that are considered universal are limited to a\ncertain set of elements, larger or smaller than OMNI-P2x. Examples include\nANI-1ccx (4 elements), ANI-2 (7 elements), MACE-OFF (10 elements), and\nAIMNet2 (14 elements). Furthermore, the set of elements included in OMNI-\nP2x (H, C, N, O, F, S, Cl) covers the vast majority of organic chemistry,\nencompassing most molecules of practical interest. We are aware that certain\nexperts only consider models universal if they cover most of the periodic table,\nbut it remains a subjective opinion without consensus in the community as no\nterminology rules are officially adopted by any of the regulation bodies (such\nas IUPAC). We and many other groups consider potentials universal if they\ngeneralize to sufficiently large space of chemical compounds to which OMNI-\nP2x certainly belongs.\nWhile excluding transition metals from the training set precludes the usage of\nthe model for an important class of compounds, we would like to point out that\nmost of ground-state MLIPs (all of the listed above) support only non-metallic\nelements. Hence, transition metals complexes remain challenging even for\nground-state ML. This is due to the fact that these systems present fundamental\nchallenges that go beyond adding new elements. Transition metals often exhibit\nmultiple oxidation and spin states, strong multireference character, and\nsignificant relativistic/spin–orbit effects, making their electronic structure much\nmore complicated. For these reasons, universal modeling of transition-metal\ncomplex photochemistry may require dedicated datasets, ML architectures and\ndescriptors, that are beyond the scope of the present work.\nReviewer 2, Comment 2\n- There is no information about the character of the states, triplet states\nenergies, nonadiabatic and spin-orbit coupling are not predicted, and there is\nno information about how to deal with ionic molecules. The model is able to\npredict vertical singlet excitations with a feature for the state order, but there\nis still no information about the orbitals involved and the actual character of\nthe excited states. The universality is claimed to be in terms of how many\ndifferent molecules can be predicted, but not in terms of photophysical\nproperties and accuracy.\nWe acknowledge the limitations of the model. The current version of OMNI-\nP2x supports only molecules within the singlet manifold of states, and does not\npredict nonadiabatic, nor spin-orbit couplings. As almost all universal or non-\nuniversal MLIPs, it does not provide information about the molecular orbitals\nof a given molecule. We added a paragraph clearly stating the limitations of the\nmodel in the “Discussion” section of the manuscript:\nWhile OMNI-P2x is the first universal excited-state neural network\npotential, it naturally has scope for further enhancements. Expanding the\ntraining set could enable support for a wider range of element types and\nout-of-equilibrium structures. Future extensions of the framework will\naim to predict not only for neutral molecules in the singlet manifold, but\nalso for states of higher multiplicity and charged species. Although the\npresent focus on singlet states already enables accurate and chemically\ndiverse predictions for many photophysical applications, the architecture\nis inherently extensible, and these additional capabilities will be the\nsubject of future work.\nReviewer 2, Comment 3\n- The training on the TD-DFT data makes for sure possible to predict energies and\noscillator strength, but this is still training on wrong results in most cases, as the\nfunctional chosen is well known to have very poor reliability for the description of\nfundamental features of excited states. So a universal model should predict good and\nphysically meaningful results, and not have as a goal to perform better than poorly\nreliable semi-empirical approaches.\nThe choice of reference method is obviously important and we fully\nacknowledge that TD-DFT/B3LYP has known limitations in describing certain\nexcited-state properties, particularly charge-transfer and Rydberg states.\nNevertheless, B3LYP remains widely used in photochemical studies and\nprovides a practical balance between accuracy and computational feasibility.\nMoreover, to the best of our knowledge, TD-DFT is the only truly “black-box”\nelectronic-structure approach currently capable of providing excited-state data\nfor a broad range of molecular sizes and compositions across chemical space\nwithout requiring case-by-case adjustments to the calculation protocol. Higher-\nlevel methods such as computationally expensive EOM-CCSD approaches,\nwhile more accurate, have not yet provided data at the scale required for training\na universal model, while multireference methods such as CASSCF or CASPT2\nrequire system-specific settings, making them unsuitable for this task.\nFurthermore, we note that OMNI-P2x is agnostic to the level of theory used and\ncan be re-trained or fine-tuned to any other, potentially more accurate, method,\nwith many examples shown in the manuscript. We also show that such transfer\nlearning, even when starting from TD-DFT/B3LYP data, can greatly enhance\nthe resulting model. This discussion has been added to the main text:\nAs with all machine-learning interatomic potentials, OMNI-P2x can only\nbe as accurate as the reference method on which it is trained. In our case,\nthis is TD-DFT/B3LYP/6-31+G*, which offers a practical balance\nbetween accuracy and computational feasibility across large chemical\nspaces, but is known to struggle with certain excited-state properties—\nmost notably charge-transfer and Rydberg excitations. In small molecules,\nthe use of a modest basis set such as 6-31+G* can further exacerbate this\nissue, as higher-lying excitations may correspond to poorly described\nRydberg states. These limitations are intrinsic to the reference data and\nthus propagate into the trained model. Nevertheless, OMNI-P2x is fully\nagnostic to the level of theory, and can be re-trained or fine-tuned on more\naccurate data when available; in fact, we show in this work that transfer\nlearning from TD-DFT/B3LYP already improves performance when\nadapted to higher-level reference methods.\nReviewer 2, Comment 4\nAdditionally, I found other major issues with the manuscript:\n1. The only task for which I would recommend a trustworthy use of the model is the\nvirtual screening of azobenzene derivatives, where being able to qualitatively and\nrapidly screen molecules can help identify and select potential candidates with\ndesirable properties. Sadly, the quantitative performance shown in the UV/Vis test case\nis very poor. Even in the (potentially cherry-picked) example from the top quartile, the\nprediction is, at least partially, wrong. The oscillator strengths are off, and between\n275 and 300 nm there are 4 states predicted vs 1 at TD-DFT/BLYP level. The authors\nclaim an excellent agreement, probably referring to the position of the maxima, but\nanyway, the convolution of the spectra is highly affected by the artificial broadening of\n0.3 eV, which might make the results look more in agreement than they actually are\n(additionally, the reference to justify 0.3 it is unrelated, as in this manuscript the error\nis between TD-DFT data and their prediction, and not a discrepancy between\ncomputed/computed at different levels or computed/experimental spectra).\nWe respectfully disagree with the remark that only trustworthy use of the model\nis the virtual screening of azobenzene derivatives. Before addressing the\ncomment about spectra, we must point out that the fine-tuning of the model for\nsurface-hopping dynamics is another clear example of the trustworthy use of the\nmodel. As for the UV/Vis spectra, we do agree with you that comparing spectral\nlines may introduce a bias: however, we believe that such a comparison is also\nimportant, as those are the fundamental, observable quantities in experimental\nstudies, as opposed to oscillator strengths and vertical transition energies.\nSingle-point convolution is a common approach and the broadening factor of\n0.3 eV is not arbitrary but was derived in the cited reference by fitting the\nsingle-\npoint convolution (SPC) spectra with those obtained with the nuclear-ensemble\napproach (NEA) on a set of representative molecules. NEA is a better approach\nfor simulating spectra, which we added to the manuscript as another trustworthy\nuse case of our model as described replies to other comments.\nTo provide an alternative ranking of the spectra not based on the criticized\nsimilarity of the spectral shapes, we added an additional benchmark where we\nscore the predictions by the mean absolute error of oscillator strengths and\nvertical transition energies and provide an analogous example in Figure S5 of\nthe ESI, with the main text now reading:\nTo provide a benchmark of the predicted spectra that is unaffected by the\nchoice of a broadening factor in the SPC construction, we also rank the\nspectra by a geometric mean of the excitation energy and oscillator\nstrength MAE, which is presented in Figure S5 of the ESI. Later, we also\nexplore the application of OMNI-P2x for spectra calculations without\nrelying on SPC approximation and the broadening factor.\nReviewer 2, Comment 5\n2. The overall idea of the UV/Vis task is to predict the brightness of states in a specific\norder. This is not enough for this task; it is just a collection of energies with no\ncorrelation with their brightness. No information about the character of the states is\npresent, and why a state is bright or not. The adiabatic order of the states (S_1, S_2,\netc.) has no correlation with what states they actually are. Two molecules could have\nS_4 both bright for two completely different and uncorrelated reasons, while two very\nsimilar molecules could have S_1 and S_2 with inverted brightness just because, for\nexample, npi* and pipi* states are inverted. So if there is no information on the mapping\nbetween excited states, predicting oscillator strength of adiabatic states is useless,\nif not wrong.\nWe thank you for this insightful comment. We agree that OMNI-P2x, like any\nML potential, does not provide explicit information on the orbital character of\nthe ground or excited states (e.g., it cannot distinguish ππ* from nπ* transitions),\nnor can it resolve cases where state ordering differs between molecules. Instead,\nthe model predicts excitation energies and oscillator strengths in the adiabatic\nordering defined by the reference method. While this limits interpretability in\nterms of electronic structure, it is consistent with how UV/Vis spectra are\ntypically constructed from ab initio calculations, where oscillator strengths are\nreported for adiabatic states without explicit orbital labeling. The utility of\nOMNI-P2x is therefore in reproducing excitation spectra (energies and\nintensities) efficiently and at scale, rather than in providing mechanistic orbital\nanalysis.\nReviewer 2, Comment 6\n3. In this logic, I also cannot understand, and it’s not explained, why the highest\naccuracy is achieved for the S0 excitations to S1 and to higher-lying S8–S10 states. I\ncannot see any correlation in that, as there is no information on the character of the\nstate in the training. I believe there might just be a high density of states in some energy\nrange that biases the training process.\nWe note that, contrary to your comment, the S0→S1 excitation energies have\nconsistently exhibited the lowest accuracy in both the original and revised\nmodels, so we are unsure where the claim of “the highest accuracy” originated.\nFor higher-lying excited states, we agree that the higher density of closely\nspaced states may make their energies easier for the ML model to learn,\nresulting in lower prediction errors.\nReviewer 2, Comment 7\n4. In the dynamics part, the impact of the new model is limited. There are no particular\nadvances with respect to previous molecule-specific training for nonadiabatic\ndynamics. The model is not transferable and universal because it needs fine-tuning at\na different level of theory, coming back to the problem of single-molecule training. It\ncould proven a computational speed-up, but not a paradigm shift in ML application in\nNAMD. Additionally, I wonder why the work previously published in Nat. Comms. by\nAxelrod et al. is not cited, where an example of transferable ML architecture for NAMD\nis actually reported. (https://www.nature.com/articles/s41467-022-30999-w). Also\nmore recent example of transferable NN for excited states should be mentioned\n(https://arxiv.org/abs/2502.12870).\nWe acknowledge that our original version might not have clearly demonstrated\nthe advantage of our model for NAMD. We have re-worked the “Efficient ML-\ndriven nonadiabatic dynamics” section to better highlight the advantages\noffered by using a foundational ML model as a basis for training MLPs for\nnonadiabatic molecular dynamics. Our new computational experiment, shown\nin Figure 10, shows that transfer learning allows for obtaining stable and\naccurate TSH simulations with an order of magnitude (2% vs 20%) less data\nthan training from scratch. This is not only a valuable computational speedup,\nbut also a paradigm shift in how MLPs for nonadiabatic molecular dynamics\nare designed, by leveraging the chemical space information contained in a pre-\ntrained model, and fundamentally changing how AL loops are designed for\nNAMD applications. The main text now reads:\nTo assess the improvement in data-efficiency that comes from using a pre-\ntrained universal model, we compare the performance of the OMNI-P2x\nfine-tuned and MS-ANI trained from scratch on fractions of the full\ntraining set required to converge an AL procedure (35071 points, as\nreported in Ref. 57). Looking at the stability of the dynamics, calculated\nas the percentage of trajectories finishing without dissociation, shown in\nFigure 10a, the fine-tuned OMNI-P2x model can propagate stable\ndynamics given just 0.1% of the dataset, while dynamics propagated using\nthe pure learning approach become consistently stable above 10% of the\ndata. Furthermore, Figure 10b presents the predicted quantum yield of\ntrans to cis photoisomerisation, as a function of the training set size. For\nthe TL potentials, all models trained on above 2% of the data predict QYs\nwithin the reference range, up to a 95% confidence interval error bar. On\nthe other hand, the performance of the models trained from scratch is\nmuch less consistent, with outliers predicting incorrect quantum yields\neven at high fractions of the training set.\nThe above findings highlight that pre-training of OMNI-P2x on a vast\ndataset allows it to gain knowledge on both the chemical compound space\nand some information about conformational space, greatly reducing the\namount of data needed for its fine-tuning. This makes OMNI-P2x a\nmethod of choice for propagating ML-accelerated surface hopping\ndynamics after data-efficient fine-tuning. A similar phenomenon has been\nreported in Refs. [58, 99] where a MLP pre-trained on ground-state data\nwas able to achieve higher accuracy for excited-state predictions than an\nMLP trained from scratch, in the low-data regime.\nWe thank you for bringing our attention to above manuscripts. We have included\ncitations to these works where appropriate (one of them in the passage above,\nwith another relevant pre-print from the same group). We would also like to\nmention, that the underlying architecture of OMNI-P2x, the MS-ANI model was\nshown in a proof-of-concept experiment to potentially be able to perform\nNAMD simulations across chemical space, as shown in Ref. 57. The main text\nnow reads:\nThis has changed with the MS-ANI (multi-state ANI) model [57]\ndemonstrated to be able to predict nonadiabatic dynamics of three distinct\nmolecules within the same model, paving the way for universal excited-\nstate potentials. More recently, another demonstration of a possibility to\ncreate transferable NN for excited-state predictions was shown in Ref. 58.\nTransferable nonadiabatic molecular dynamics simulations were also\nreported by Axelrod et al. [59] however, it was limited to the class of\ndiazobenzene photoswitches.\nReviewer 2, Comment 8\n5. The fulvene dynamics look very different from the reference in literature, where a\nreflection in the S_1 population is present, particularly with trajectory surface hopping\n(Phys. Chem. Chem. Phys., 2020,22, 15183-15196). This is due to a physical reason,\nwith the molecule encountering two different conical intersections, one of them\npopulating back the first excited states. This is the reason why this model has been\nchosen to test nonadiabatic methods. Already the reference method used in this\nmanuscript misses this effect. The training with the model actually fully deletes the\nsmall reflection at around 15 fs, removing completely the physical behavior expected.\nI believe this is more due to interpolation reasons than capturing the complexity of the\npotential energy surface, but in general is an example where the prediction erases the\ncorrect physical behavior and potentially leads to misinterpretation. This could be\nsolved using another formulation of TSH, for which nonadiabatic coupling or\nwavefunction overlap would need to be predicted by the model, but is not available in\nthe architecture presented in this work.\nThank you for bringing our attention to this. We would like to point that the\nreflection populating back the excited state is present in the reference dynamics,\nhowever, its magnitude is indeed smaller than in the FSSH dynamics.\nFurthermore, the fine-tuned ML model captures this physical behavior\nsignificantly better than the no-TL model, trained on the same dataset. To\nemphasize this, we have included an inset in Figure 9, zooming in on the 10 fs\nto 20 fs region of the simulated populations. We have added 95% confidence\ninterval error bars showing the statistically significant difference between the\nMLIPs trained from scratch, and using transfer learning, while also increasing\nthe number of trajectories propagated using both ML approaches to 10,000, as\nshown in the new Figure 9.\nThe text accompanying this figure now reads:\nBy employing transfer learning from a universal potential, we are able to\nslash the convergence time of fulvene from over three days, to about 24\nhours, with a ca. three-fold increase in computational efficiency.\nAdditionally, the final employed training set contains only 3550\nconformations, which is about twice less as was necessary to achieve\nconvergence in a reference direct active learning procedure, as can be seen\nin Figure 9. Moreover, inspecting the time evolution of the electronic state\npopulations, we can see a statistically significant difference in the\npopulations predicted by the fine-tuned model and the previous MS-ANI\nmodel. Despite a smaller training set, the TL version of OMNI-P2x better\ncaptures the photophysical behavior of the system.\nThe agreement between the FSSH and LZSH methods of propagating NAMD is greater\nif the reduced kinetic energy reservoir option is disabled. A comparison of populations\nobtained by propagating TSH with three different approaches: ML-LZBL (no reduced\nkinetic energy), ref-LZBL (no reduced kinetic energy) and FSSH is provided as Figure\nS12 in the ESI, showing that the LZBL surface hopping scheme can accurately\nreproduce FSSH. However, as the main goal of designing ML potentials is replicating\nthe underlying reference calculations, we choose to keep the original comparison, as\nthat was the computational setup used in the active learning procedure.\nAdditionally, a comparison of the dynamics predicted using different ML models and\ndifferent LZBL propagation settings with reference fewest switches surface hopping\ndynamics is available in Section S4 of the ESI.\nReviewer 2, Comment 9\n6. More information about the computational cost of training would be necessary. Yes,\nthe prediction is definitely cheaper than TD-DFTB, but working on a large dataset,\nwhat is the computational cost of the training with respect to calculating these\nexcitations with the semi-empirical method?\nAs OMNI-P2x is a foundational model, it can be used to make predictions on\nany dataset without further training; hence, we believe it is fair to compare the\nmodel inference time with the prediction speed of the semi-empirical methods.\nIn the same spirit, semi-empirical methods require time (often years, as in the\ncase of ODM2) to parametrize, which is not included in the total computational\ncost of individual calculations, once the method is already developed and\nparametrized. However, we agree with the referee that including training times\nand hardware information is beneficial to the clarity and fairness of the\nmanuscript; hence, we include an additional “Training times” section in the ESI,\nwith this information provided.\nThe 4 main models of OMNI-P2x (ensemble of three energy-predicting NNs and\noscillator strength NN) have been trained for 100 epochs on 32 NVIDIA GH200 GPUs\nwith 96GB memory each, which totaled to about 48 hours of training time per model.\nFine-tuning of OMNI-P2x for NEA spectra was performed in a Jupyter Lab\nenvironment, with 4 Intel Xeon Gold 6226R CPU cores, taking under a minute for the\npyrazine case, and about 40 minutes for 9-DCMA. The AL loop for NAMD simulations\nof fulvene was run on a computational node with a single NVIDIA 4090 RTX GPU and\nfor NAMD simulations was performed with the same computational resources, taking\nabout 20 minutes per trained model.\nMinor points:\nReviewer 2, Comment 9\nIs the label of Figure 4a correct? I think the authors show the ground state energies.\nThank you for your vigilance. This caption was indeed a typo, and was changed\nto read:\nTest performance of the OMNI-P2x model. Panel (a): correlation plot of\nthe reference~(true) and OMNI-P2x predicted ground-state energies at the\nCCSD(T) level of theory, for the test part of the ANI-1CCx dataset; (b)\nMAE of the electronic state energies; c) MAE of the excitation energies;\nd) MAE of the predicted oscillator strengths.\nReviewer 2, Comment 10\nThe statement “While these methods are certainly more affordable for systems with\nhundreds and even thousands of atoms, the accuracy of their predictions can often be\nquestionable” needs a reference.\nWe thank you for their comment. References have been provided to justify that\nclaim:\nWhile these methods are certainly more affordable for systems with\nhundreds and even thousands of atoms, the accuracy of their predictions\ncan often be questionable. [26–28]\nReviewer 2, Comment 11\nSome self-glorifying strong claims like “This marks the watershed moment that ML\napproach for excited-states presents a viable alternative to established QM methods.”\nshould be toned down, as they just represent the author's personal opinions.\nThis sentence has been re-formulated to read:... making ML approaches for\nexcited-states a viable alternative to established QM methods.\nReviewer #3\nThis communication presents a novel foundation model / universal model for the excited\nstate potential energy surfaces of organic molecules. The introduction motivates the\nmodel well in terms of deficit of the capabilities of existing computational methods for\nexcited states in combination with high-throughput or dynamics-based studies of\nanything beyond very small molecular systems. The approach taken is to adapt a\nvariant of the well-established ANI neural network potential, recently-introduced by\nthe current authors and other collaborators in [54] npj Computational Materials 11,\n132, which introduces a state-index as an extra input to the neural network alongside\nthe geometry and atomic numbers, resulting in a “MS-ANI” approach. The current\nwork builds on the demonstrations in [54] by showing that a foundation model, OMNI-\nP2X, can be built with this architecture which can be general across a fairly wide range\nof organic molecules. Overall the model is demonstrated to be quite successful and I\nthink the paper could be quite influential as it represents a first of its kind. However, I\nfeel there is room for improvement in the presentation and demonstrations before\npublication: in several sections the results themselves are very briefly described and it\nis hard to discern the level of success of the model.\nThank you for the positive evaluation. In the presented revision, we have taken\nnumerous steps to increase the accuracy of the OMNI-P2x model by doubling\nthe training set size, increasing the computational effort used in the training\nprocess, and employing an ensemble of NN models. Furthermore, the\nmanuscript has been revised to present more clearly the advantages and\nlimitations of OMNI-P2x. An additional workflow for fine-tuning the model to\npredict spectra based on non-equilibrium structures using the nuclear ensemble\napproach has been added in a new section, “Nuclear ensemble approach\nspectra”.\nReviewer 3, Comment 1\nTraining data has been sourced from two existing databases, the PubChemQC\ndatabase of TD-B3LYP/6-31+G* calculations, and the ANI-1ccx dataset of small\nmolecules including off-equilibrium configurations, which has been provided as ground\nand excited-state DFT and TDDFT data for 10 excitations, and as ground-state-only\nCCSD(T), allowing the model to predict multiple levels of theory. Results show the S0\nground state has an MAE on energies of 2.63 kcal/mol which makes it comparable to\nother models trained on large databases of organic molecules, though I would prefer\nto see an explicit comparison to understand clearly how good the model is for the\nground state.\nWe agree that such a comparison would be beneficial. However, ground state\nwas not the main focus of this study, and better results can be obtained with\nmore ground-state training data (we used a quite small data set compared to\nother universal ground-state potentials). Comparison to other MLIPs for\nground-state properties has been added to the main text, comparing the ground-\nstate predictions at the DFT level:\nThis ground-state MAE can be put in perspective to the similar metrics of\nother MLIPs: ANI-1 with an RMSE of 1.3 kcal/mol (OMNI-P2x 2.93\nkcal/mol), UMA [76] with around 1 to 3 meV/atom, depending on the\nversion (OMNI-P2x 4.48 meV/atom), and AIMNet2 [77], reporting a\nMAE of around 1.25 kcal/mol for neutral compounds. Although those\nresults are not directly comparable, as each of the models was trained on\na different dataset, they demonstrate that the ground-state error of OMNI-\nP2x exhibits decent performance, even with respect to the current state-\nof-the-art models, despite being trained on much smaller ground-state data,\nas ground-state performance was not the primary focus of this study.\nReviewer 3, Comment 2\nThe S1 state and higher-lying excitations are substantially worse, at ~5 kcal/mol, which\nis still useful, but it would be interesting to understand how correlated are the errors\nbetween S0 and the other states, ie should we expect errors on S0->Sx transition\nenergies to exceed those of the contributing states, or will there be error cancellation\nto some degree? Scrutiny of Figs 4b and 4c might provide this information in an\naveraged sense, but they have different energy units so the comparison is not\nstraightforward (generally, I find the broad selection of energy units (variously Ha,\neV, and kcal/mol, sometimes in the same figure) rather troubling.\nThank you for your insightful questions. The units used in Figure 4 have been\nunified to kcal/mol. Furthermore, the version of the model presented in the\nmanuscript has been revised, replacing the original OMNI-P2x with an\nensemble of three neural networks, newly trained on the entirety of the\nPubChemQC dataset. The ensemble outperforms the previous version of the\nmodel by about 30% in terms of the MAEs. Furthermore, a detailed breakdown\nof the performance of OMNI-P2x predictions of vertical excitation energies has\nbeen provided in the ESI, featuring correlation plots with R2 metrics and MAE\nfor all 10 transitions. The excitation energy errors are similar in magnitude to\nthe target state energy errors and follow the same trend: largest for the S0→S1\ntransition and decreasing for higher-lying transitions.\nReviewer 3, Comment 3\nFigure 3 shows the distribution of sizes and excitation energies of the training data,\nindicating that the majority of the PubChemQC data is in the 20-40 atom range\nwhereas that of the ANI-1cxx data is much smaller, peaking at 12 atoms. Figure 4 shows\nsome test results: I note that while figure 4a purports to demonstrate close agreement\nbetween predicted and target total energy, the spread of energies is so wide (1000\nHartrees!) that presenting this as a y=x figure is highly redundant – any deviation from\ny=x visible at this scale would represent an enormous error, and the MAE of\n5.78kcal/mol ~ 0.01 Hartree would be invisible on this scale. Perhaps a baseline\nsubtraction of the atomic energies from each case would make this graph more useful?\nThank you for this suggestion on improving the presentation of our results.\nPanel a of Figure 4 has been re-worked to show the MAE error on normalized\nenergies (with removed self-atomic energies), and in the main text, we report\nthe MAE, as well as the correlation of the normalized and un-normalized\nenergies:\nFor the ground-state energies in the test ANI-1ccx subset, OMNI-P2x has\na higher MAE of 4.14 kcal/mol, but still excellent R2 = 0.9957 for the\nnormalized energies (0.9999 for unnormalized).\nReviewer 3, Comment 4\nRegarding the fact that 10 excited states are calculated: are this many states always\ngoing to give meaningful results in these systems with a 6-31+G* basis: in small\nmolecules these might be higher-lying Rydberg states which are not well-represented\nin a relatively small basis such as 6-31+G*?\nWe thank you for this comment. We agree that for smaller molecules, TD-DFT\nwith a modest orbital basis set such as 6-31+G* may, in certain cases, yield\nunphysical results for higher-energy transitions. However, selecting a larger\nbasis would lead to prohibitively expensive calculations, considering a training\nset consisting of over 3k medium-sized molecules. It is also unfeasible to define\none excitation-number cutoff limit that would exclude Rydberg states from the\ntraining set, without manual individual, manual inspection. At the end of the day,\nthe applicability of all ML models can never go beyond the applicability of the\nunderlying reference method, and the user should carefully judge if the best-\ncase scenario: exact replication of TD-B3LYP/6-31+G* results is suitable for\ntheir study. The manuscript has been updated with that information:\nAs with all machine-learning interatomic potentials, OMNI-P2x can only\nbe as accurate as the reference method on which it is trained. In our case,\nthis is TD-DFT/B3LYP/6-31+G*, which offers a practical balance\nbetween accuracy and computational feasibility across large chemical\nspaces, but is known to struggle with certain excited-state properties—\nmost notably charge-transfer and Rydberg excitations. In small molecules,\nthe use of a modest basis set such as 6-31+G* can further exacerbate this\nissue, as higher-lying excitations may correspond to poorly described\nRydberg states. These limitations are intrinsic to the reference data and\nthus propagate into the trained model. Nevertheless, OMNI-P2x is fully\nagnostic to the level of theory, and can be re-trained or fine-tuned on more\naccurate data when available; in fact, we show in this work that transfer\nlearning from TD-DFT/B3LYP already improves performance when\nadapted to higher-level reference methods.\nReviewer 3, Comment 5\nIt is interesting that the Energy MAE on S1 is roughly double that of S0, then further\nhigher-lying excitations have increasingly lower errors. It strikes me that perhaps the\nincreasing Rydberg character of higher-lying states is somehow easier for the model to\ncapture. The authors say “An outlier here is S0, which has the smallest energy MAE of\n2.63 kcal/mol, which we attribute to the effect of the added ground-state CCSD(T)\ndata” – I couldn’t see how that could come about – why would the presence of data\nassociated with a second level of theory predictable by the same model improve the\naccuracy of its predictions for the first?\nWe are thankful for this insightful comment. Regarding the prediction accuracy\nof higher excited states, it is indeed possible that their character makes machine\nlearning fitting easier. In particular, the density of states (i.e., the number of\nstates per energy unit) increases with excitation number; if the model captures\nthis trend, it can lead to lower MAEs for higher-lying states.\nConcerning the improved accuracy for S when training on multiple levels of\n0\ntheory, this phenomenon has been discussed in the all-in-one learning (AIO)\nframework (ref. 63, 10.1021/acs.jctc.5c00858.). AIO-trained models often show\nsimilar or even larger improvements compared to transfer learning models pre-\ntrained on DFT data and fine-tuned to coupled-cluster quality. Specifically,\ntransfer learning from ground-state data has been shown to enhance excited-\nstate property predictions (ref. 58, 99), and we expect AIO to yield analogous\nbenefits here. At the same time, the dramatic improvement observed for S ,\n0\nroughly a twofold decrease in MAE, likely also reflects characteristics of the\ntraining data itself, as ground-state energies generally vary more smoothly with\nmolecular geometry, exhibit lower variance, and follow distributions that are\neasier for the model to learn.\nThe manuscript has been updated to read:\nAn outlier here is S , which we attribute to the improved data-efficiency\n0\ndue to the application of the all-in-one learning to the combined set with\nthe ground-state DFT (from PubChem) and the CCSD(T) (from ANI-1ccx)\ndata, and an overall smoother and lower-variance distribution of ground-\nstate energies that makes them easier for the model to learn.\nReviewer 3, Comment 6\nSome minor points on figures: “One can note, that the highest accuracy is achieved for\nthe S0 excitations to S1 and to higher-lying S8–S10 states.” Based on Figure 4d this\ndoes not appear to be true – do they mean the lowest accuracy?\nThank you for catching this mistake. As the entire OMNI-P2x model was\nrevised and re-trained, the sentence was changed to reflect the true accuracy of\nthe predicted oscillator strengths:\nOne can note that the highest accuracy in terms of the predicted oscillator\nstrengths is achieved for the excitations to S and to higher-lying S –S\n1 8 10\nstates.\nReviewer 3, Comment 7\nIn the stick plot in Figure 5: the baseline of sticks does not align with zero on right hand\naxis – this needs correcting.\nThank you for noticing this: the plots have been aligned so that the oscillator\nstrength axis aligns with the bottom of the stick spectrum.\nReviewer 3, Comment 8\nPanel 8a is very separate from 8b,8c: they refer to different systems and unrelated\ncalculations.\nWe agree with you on this matter. We have now separated Figure 8 into Figures\n9 and 10, with a separate figure showing the results of the active learning loop\nfor fulvene, and a separate figure showing the usage of OMNI-P2x in transfer\nlearning models for simulating the photoisomerization dynamics of azobenzene.\nFurther discussion of these plots follows in a response to a question further in\nthis letter.\nReviewer 3, Comment 9\nI would say that the Statement “OMNI-P2x is approaching the accuracy of time-\ndependent density functional theory” in the abstract and main text is potentially\nmisleading, despite being carefully written. As Figures 5 and 6 show, while the model\ndoes significantly better than semi-empirical methods, at reproducing TDDFT spectra,\nas measured by the Spearman Correlation Coefficient, the reasons for this in terms of\nagreement of the underlying stick spectrum with TDDFT are often dubious. In many\ncases, where TDDFT has a few strong peaks, OMNI-P2X has a forest of intermediate\npeaks in a similar energy range, with no clear relationship between oscillator strengths.\nI would want to know what the Mean Relative Error on oscillator strengths was, as well\nas the MAE, because it seems like in many cases it might be quite large.\nFor a method to be able “approach the accuracy of TDDFT” I would expect individual\nexcitations to be predicted to better than the 5kcal/mol quoted in Figure 4b as the\nMAE on S1 energies, and for the mean relative error on oscillator strengths not to be\ntoo extreme. This current level of energy error corresponds to 0.22eV or 100nm on a\nS0<-S1 excitation around 1.5eV, which is not really approaching TDDFT accuracy.\nThank you for pointing this out. We agree that the accuracy of predicting\nspectral lines may not necessarily be correlated with the accuracy of the\npredicted underlying stick spectrum. To quantify this, we reproduce Figure 7\nusing the total error of the absorption energies and oscillator strengths as a\nranking, instead of the Spearman Correlation Coefficient. The results are\navailable in Figure S5 of the ESI.\nWe also agree that the level of error presented in the original model, of 0.22\neV, is rather large, yet it was still comparable to the typical error of TD-DFT\npredictions of 0.2-0.3 eV. However, the revised ensemble of NNs is able to\nachieve a lower error, below 4.0 kcal/mol (0.17 eV) for all states, making the\nerror related to fitting the ML model lower than the inherent error of TD-DFT.\nFurthermore, a more thorough comparison is available in Figure S4 of the ESI,\nshowing correlation plots for excitation energy predictions for all 10 transitions.\nThe main text now reads:\nFurthermore, the MAEs for both the predicted total and excitation energies\nare below 0.2 eV (4.6 kcal/mol), which lies below the lower bound of the\nreported accuracy for the reference TD-DFT method, indicating that the\nfitting error of the ML model is smaller than the intrinsic error of the\nunderlying quantum chemical approach. Correlation plots showing the\npredicted vs true excitation energies for all transitions are available as\nFigure S4 of the ESI. In all cases, the R2 is above 0.91.\nReviewer 3, Comment 10\nThe section on design and screening of azobenzene compounds provides a proof-\nof-principle for use in screening large datasets of candidate molecules. In this case\nazobenzene derivatives are screened to find 21 leads with strongly red-shifted\nabsorption, at least some of which stood up to further scrutiny via calculations at a\nhigher level of theory, ADC(2). I found it frustrating that this section ended without\nserious investigation of whether OMNI-P2X had really found strong leads: no\ncomparison is made of the spectra from TD-DFT, from ADC(2) or from OMNI-P2X.\nI would want to see this investigated more deeply.\nWe agree with you on that. Indeed, a more thorough comparison of the different\nlevels of theory used to screen the azobenzene derivatives would be beneficial.\nFollowing this suggestion, we have now included the SPC spectra for the top 3\ncandidates presented in Figure 7 in the Electronic Supporting Information.\nOverall, we note a good level of agreement between the methods, with S and\n2\nand 0.42 and 0.35 eV between TD-DFT and ADC2 and OMNI-P2x and ADC2,\nrespectively. Furthermore, molecular databases containing molecules with state\nenergies at the three discussed levels of theory are provided, along with a\nJupyter notebook and instructions on how to generate and compare spectra for\nany given compound.\nReviewer 3, Comment 11\nThe final section demonstrates what could turn out to be a common use-case for\nfoundation models for the excited state: active-learning-driven fine-tuning / transfer-\nlearning is used to train a OMNI-P2X model at the CASSCF(6,6) level of theory for\nfulvene. The transfer-learning is shown to substantially accelerate model training, and\nboth the TL and no-TL versions accurate reproduce dynamics of the reference CASSCF\ncalculations. Very brief mention is made of calculations on photoisomerisation of\nazobenzene, but the results presented are rather cryptic: I was not able to intuit what\npoint was being made about dihedral angles as obtained from the dynamics.\nWe agree with your observation and now substantially expand the section\nconcerning the nonadiabatic molecular dynamics, with particular focus set on\nthe example of azobenzene. The results for fulvene have been separated (Figure\n9) from the azobenzene results, presented in Figure 10. We also added a zoom-\nin inset on the crucial state-population region between 10 and 20 fs of the\ndynamics, which clearly shows the difference in populations predicted by the\nTL and no-TL models, with the version derived from OMNI-P2x predicting\nquantitatively accurate dynamics, in agreement with reference populations.\nAdditionally, shaded 95% confidence interval error bars have been added to\nshow the statistical significance of the results. The main text now reads:\nBy employing transfer learning from a universal potential, we are able to\nslash the convergence time of fulvene from over three days, to about 24\nhours, with a ca. three-fold increase in computational efficiency.\nAdditionally, the final employed training set contains only 3550\nconformations, which is about twice less as was necessary to achieve\nconvergence in a reference direct active learning procedure, as can be seen\nin Figure 9. Moreover, inspecting the time evolution of the electronic state\npopulations, we can see a statistically significant difference in the\npopulations predicted by the fine-tuned model and the previous MS-ANI\nmodel. Despite a smaller training set, the TL version of OMNI-P2x better\ncaptures the photophysical behavior of the system. Additionally, a\ncomparison of the dynamics predicted using different ML models and\ndifferent LZBL propagation settings with reference fewest switches\nsurface hopping dynamics is available in Section S4 of the ESI.\nTo assess the improvement in data-efficiency that comes from using a pre-\ntrained universal model, we compare the performance of the OMNI-P2x\nfine-tuned and MS-ANI trained from scratch on fractions of the full\ntraining set required to converge an AL procedure (35071 points, as\nreported in Ref. 57). Looking at the stability of the dynamics, calculated\nas the percentage of trajectories finishing without dissociation, shown in\nFigure 10a, the fine-tuned OMNI-P2x model can propagate stable\ndynamics given just 0.1% of the dataset, while dynamics propagated using\nthe pure learning approach become consistently stable above 10% of the\ndata. Furthermore, Figure 10b presents the predicted quantum yield of\ntrans to cis photoisomerization, as a function of the training set size. For\nthe TL potentials, all models trained on above 2% of the data predict QYs\nwithin the reference range, up to a 95% confidence interval error bar. On\nthe other hand, the performance of the models trained from scratch is\nmuch less consistent, with outliers predicting incorrect quantum yields\neven at high fractions of the training set.\nThe above findings highlight that pre-training of OMNI-P2x on a vast\ndataset allows it to gain knowledge on both the chemical compound space\nand some information about conformational space, greatly reducing the\namount of data needed for its fine-tuning. A similar phenomenon has been\nreported in Refs. [58, 99], where an MLP pre-trained on ground-state data\nwas able to achieve a lower MAE for excited-state predictions than an\nMLP trained from scratch, in the low-data regime.\nFigure 10 now shows a new computational experiment demonstrating the\nimprovements associated with using OMNI-P2x as a pre-trained model for\nconstructing MLIPs for propagating trajectory surface hopping simulations.\nNamely, we measure the stability of the conducted simulations, as well as the\npredicted quantum yields, as a fraction of the full AL dataset from Ref. 57, for\nTL and no-TL models. The results show that OMNI-P2x can perform stable\nsimulations (without any dissociated trajectories) with as little as 0.1% of the\nfull training set and achieves quantitative agreement with the reference\npredicted quantum yields with orders of magnitude less training data than ML\nmodels trained from scratch.\nReviewer 3, Comment 12\nOne natural question would be how the model responds to large molecular systems and\nrelatively large numbers of states required, and is there consistency of the predictions\nwith respect to the number of states requested. In some testing of the version made\navailable on figshare, I observed that when calling OMNI-PX2 three times, for 2, 5 and\n10 states, for a model of a sunscreen candidate molecule with 95 atoms, in each case I\ngot a very different set of excited states predicted (the S1->S0 transition was much\nlower in the 10 state case than the 2 state case). Is there anything that guarantees\nconsistency with respect to the number of states requested. I note it is relatively stable\nfor the provided example of benzene, though the oscillator strength varies a little.\nWe thank you for the question. The OMNI-PX2 model is designed as a multi-\nstate model, in which predictions for each electronic state are made\nconsecutively (i.e., with the neural network called for x = 0, 1, 2 … up to the\nrequested number of states). As a result, the predicted energies are, in principle,\nindependent of the number of states requested.\nWe were able to reproduce this effect using Bemotrizinol (a 95-atom sunscreen\nmolecule). In this case, the predicted energies and oscillator strengths remain\nconstant; however, because the molecule is significantly larger than the\nmolecules in the training set (median ~30 atoms, maximum 65 atoms), the\nordering of higher electronic states may differ from the physically correct order.\nThis can lead to apparent lower-energy peaks in the simulated spectrum, even\nthough the model’s predictions for individual states are consistent.\nThe revised release of OMNI-PX2 addresses such situations through an\nuncertainty quantification system, which issues warnings when predictions are\nhighly uncertain, potentially due to the choice of a system far from the training\nregime.\nThe manuscript has been updated to clearly state:\nFurthermore, the neural network model is expected to produce reliable\npredictions primarily for molecules whose size and chemical composition\nare similar to those represented in the training set; extrapolation to\nsignificantly larger or chemically distinct molecules may lead to\nunphysical or unexpected results. The built-in uncertainty quantification\nmechanism warns the user about predictions that are not covered well by\nthe training set.\nResponses to the comments of the reviewers\nReviewer comments are in italics, and our responses are in regular, blue font. Our changes to the\nmanuscript are marked in red font.\nReviewer #1\nThank you very much for the detailed reply to my comments. After the revision, I think\nthe limitations are more openly communicated and indicated in the manuscript. As the\nmanuscript is certainly a step forward towards generalizable prediction of excited state\nenergies, I support the publication in Nature Communications..\nWe thank you for your positive opinion of the manuscript and your thought-\nprovoking comments, which led to substantial improvements.\nDimers: Thank you for the additional experiments. The fragment-correction scheme is\nvery interesting and a valuable addition to the manuscript.\nLevel of theory encoding and excited states encoding: Thanks for the additional\nablation study.\nOut-of-equilibrium: Thanks for the very interesting additional tests, particularly Figure\nS10 and the fine-tuning study.\nFigure 4: Thanks for the changes.\nFigure 6: It is a very interesting observation that out-of-equilibrium errors are higher\nfor lower excited states. Thanks for adding this information.\nPage 1 of 28\nDelta-Learning: I understand the complexity of this topic, and I am looking forward to\nseeing a more detailed study of this. I agree that it might warrant a separate paper.\nTransfer learning: Thank you for making this clearer and adding the new figure to\nindicate the impact of transfer learning.\nReviewer #2\nI thank the authors for addressing my comments and improving the manuscript, especially thanks to the new analysis\nreported.\nIn particular, the atom-wise energy decomposition and the geometrical analysis of fulvene dynamics make the interpretation\nof the results more reliable and robust. However, for the atom-wise energy composition analysis, a general metric to\nevaluate the accuracy would be more suitable than potentially cherry-picked examples. The comparison between B3LYP\nand ADC(2) is also important, showing that, although correlated, the performances of TD-DFT are systematically worse than\nADC(2).\nI think the manuscript is highly improved through the two revision rounds, and I have no additional comments for the authors\nthat prevent me from recommending publication.\nThank you for your careful review, which led to the manuscript improvement.\nHowever, I believe that my chemistry-related comments have not been sufficiently scientifically addressed. Especially, if the\nmodel has already received high visibility as the authors report, it is even more important to ensure that it delivers physically\ncorrect results.\nBelow, we address the remaining chemistry-related comments and hope that this\nrevision will find your approval.\nPage 2 of 28\nAdditionally, I suggested some potential modifications that might make the impact of the model strong enough to recommend\npublication, like the prediction of the character of the excited states, nonadiabatic coupling, or spin-orbit coupling. However,\nnone of these have been addressed, although in some cases, even acknowledged that architecture is inherently extensible\nfor this. That was not only required to show a justification of the novelty of the model to be published in Nat. Comm., but also\nto avoid potential flaws in the chemical interpretation that can be obtained by using the model as it is.\nWe acknowledge that interpretations are as important as massive simulations,\nwhich our model enables. We agree with you that elucidating the character of\nthe excitations is one of the most crucial tasks in theoretical modeling of\nexcited-state processes. Hence, in the revised manuscript, we show that the\nOMNI-P2x model provides a unique insight into the character of excitations.\nThis insight is based on an analysis of the atom-wise contributions to excitation\nenergies, which are readily obtained from the OMNI-P2x predictions. The\nOMNI-P2x-derived atom-wise contributions to excitation energies agree well\nwith the density difference plots from the reference TD-DFT calculations.\nHence, in the revised manuscript, we have introduced a new interpretable\nmachine learning approach for excited states.\nThe new section of the revised manuscript provides the details of this new\napproach:\nTo provide interpretability of the ML-learned excitations, we introduce\na new approach for analyzing atom-wise contributions to excitation\nenergies. For each atom i, the fractional contribution is defined as\n|ΔEi|/Σ|ΔEi|, where ΔEi is the excitation energy contribution for a given\ntransition from atom i, calculated as the difference between ground- and\nexcited-state atomic energies predicted by the NN ensemble. We show\nPage 3 of 28\nthat the OMNI-P2x-derived atom-wise contributions agree well with the\nreference TD-DFT density difference plots (Fig. 8). Such analysis\nallows us to obtain the spatial character of the excitation by identifying\nwhich atoms are most involved in the transition. Our analysis is\nconceptually related to work by Tkatchenko et al [83], analyzing atom-\nwise contributions to ground-state energies.\nWe demonstrate that this interpretability analysis can work well across\nchemical space and for different excitations of the same molecule. Panel\na) of Fig. 8 shows the S → S ML ΔE contributions and TD-DFT\n0 1\ndensity differences for several organic molecules. We begin by\ncomparing the examples of benzene and pyrazine, two 6-π electron\naromatic compounds. OMNI-P2x correctly localizes the transitions,\nspreading it evenly across all the benzene carbon atoms, while in\npyrazine it is mostly localized on the two nitrogen atoms, reproducing\nthe TD-DFT results. The excitations also exhibit appropriate symmetry:\nC for benzene, and D for pyrazine. The next example, azobenzene,\n6 2h\nshows physically correct localization of the excitation on the N=N\ndouble bond, while in p-nitroaniline the partial charge-transfer character\nPage 4 of 28\nof the excited state is well reproduced. Overall, this shows an\ninterpretable character of ML predictions, that is transferable across\nchemical space.\nAtomic contributions to different electronic transitions of the same\nmolecule are shown in Fig. 8b), using azobenzene as an example.\nOMNI-P2x correctly reproduces the S → S transition as localized on\n0 1\nthe central N=N bond, while the S transition has larger contributions\n2\nfrom the aromatic rings. Finally, a higher excited state (S ) is shown,\n11\nwhich is delocalized across the entire conjugated system.\nI would like to justify my claims about the problems that the authors addressed and that persist in this revised version of the\nmanuscripts.\nWe thank you for acknowledging the improvements in the manuscript, and the\nexpansion of the NAMD results of azobenzene. We would like to clarify, that\nby “extensible”, we meant the extensibility to charged systems, and systems of\nother multiplicity, which were mentioned in your original comment, as stated in\nthe manuscript: “Future extensions of the framework will aim to predict not\nonly for neutral molecules in the singlet manifold, but also for states of higher\nmultiplicity and charged species. Although the present focus on singlet states\nalready enables accurate and chemically diverse predictions for many\nphotophysical applications, the architecture is inherently extensible, and these\nadditional capabilities will be the subject of future work.” As for SOCs and\nNACs, we make no such claims. The problem of predicting couplings that are\nwidely transferable across chemical space remains an open problem. None of\nthe reported models that incorporate non-adiabatic couplings can predict them\nacross chemical space, and the feasibility of such a universal model is unknown.\nWe would like to reiterate that by “universal potential”, we use the common\ndefinition signifying a model that has been pre-trained on a vast region of\nchemical space, and can be used out of the box for unseen compounds, not a\nPage 5 of 28\nmodel that can universally predict every relevant quantity that can be obtained\nwith quantum chemistry. The term universal has been applied to many models\nthat predict only the ground-state energy, and its derivatives (eg. ANI-1,\nMACE, UMA). To avoid misleading anyone, a definition of a universal MLP is included\nin the manuscript:\nA particularly important, emerging class of potentials is universal MLPs,\nwhich we define as models that were trained on vast datasets spanning\nacross chemical space. In turn, they can be used to make out-of-the-box\npredictions for previously unseen molecules, generalizing from the\ntraining data.\nThe authors say:\n“However, many (if not most) machine learned interatomic potentials (MLIPs) that are considered universal are limited to a\ncertain set of elements, larger or smaller than OMNI-P2x.”\nI agree with the authors. However, not clearly going beyond the state-of-the-art in this regard, it shows indeed the limited\nimpact of the model in terms of offering a new paradigm for using MLIP for computational chemistry predictions.\nWe respectfully disagree with the comment that OMNI-P2x does not go beyond\nstate of the art. OMNI-P2x goes beyond the current state of the art for MLPs by\nproviding predictions for excited states, which no universal MLP has done\nbefore. This is a new paradigm in ML-assisted computational chemistry, while\ninclusion of more elements, while valuable, is an incremental change that will\nbe the subject of further work, in the same manner as it is done for ground-state\nmodels.\n.\nPage 6 of 28\nMy concerns about the correctness of the results used for training the model have not been sufficiently addressed or disproved.\nThe authors say: “The choice of reference method is obviously important and we fully acknowledge that TD-DFT/B3LYP has\nknown limitations in describing certain excited-state properties, particularly charge-transfer and Rydberg states. Nevertheless,\nB3LYP remains widely used in photochemical studies and provides a practical balance between accuracy and computational\nfeasibility.“\nFirst of all, I disagree with the statement. Although it could reproduce good energies for some bright states, the danger of\nusing B3LYP in excited states calculations over its benefits is very well established. That’s why several range-separated\nhybrids, or functionals with higher percentages of Hartree-Fock exchange, are strongly recommended to be used for excited\nstates calculations, and they do not add any relevant computational expense. However, the widely used approach in literature\ndoes not justify its correctness, especially because the performances of density functional approximations are highly case-specific.\nIt would have been necessary to prove the accuracy of B3LYP with higher reference calculations or a systematic evaluation of\nthe performance. This would have been very easy to show, for at least a pool of selected molecules, but it has not been done.\nUnfortunately, acknowledging the limitations of B3LYP, but saying that it is highly used, without showing a rigorous evaluation\nof its prediction, does not support the claim supported by scientific evidence and limits the trustworthiness of the accuracy of\nthese calculations that are used to train the model and consequently its reliability to be used for unseen molecules.\nHowever, I thank the author for addressing this problem in the discussion.\nWe agree with the need to provide additional information and benchmarks for\nthe underlying quantum chemistry method of OMNI-P2x. Fortunately, a great\nmany literature benchmarks of B3LYP performance have been established in\nPage 7 of 28\nthe past. Consistently, B3LYP has emerged as one of the better performing\nfunctionals, with MAE errors ranging between 0.2 and 0.4 eV, with respect to\nexperiment (Refs. 101 and 102) or a theoretical best estimate (TBE) (Refs. 103-\n106). We supplement this result with our own benchmark, comparing TD-\nDFT/B3LYP results with energies obtained with ADC(2), which is a higher-\nlevel method, suggested by you in the review, with results provided below:\nOverall, there is a very good correlation between the ADC(2) results and\nB3LYP, with a correlation coefficient of about 0.95 for all electronic states. The\nMAE ranges from about 0.24 eV for S0 -> S1 transitions, up to 0.59 eV for S0-\n> S11. However, it is worth noting the mean signed error of each transition,\nwhich is negative and in almost all cases equal to approximately the MAE value,\nwhich indicates a systematic underestimation compared to ADC(2), rather than\nrandom errors. Full correlation statistics have been presented in Figure S19 and\nTable S3 of the ESI.\nPage 8 of 28\nIn terms of oscillator strengths, benchmark studies by Thiel et al [106] show\ngood correlation between B3LYP and CC2 results, with a MAE of 0.06, with\naccurate predictions of transition dipole moments demonstrated by Robinson et\nal on a smaller test set [109].\nFurthermore, it should be noted that benchmarks comparing ADC(2) and\nB3LYP results show a somewhat better performance of ADC(2) when\npredicting excitation energies (0.1 vs 0.2 eV MAE in Ref. 103, 0.12 vs 0.30 in\nRef. 108, and 0.21 vs 0.45 Ref. 105), the benchmark of Sarkar et al [105] shows\nthat ADC(2) and TD-DFT/B3LYP predict oscillator strengths of about the same\naccuracy, with an error of 15.4% and 17.6% respectively, when compared to\nTBE. Hence, we do not compare oscillator strengths predicted by these two\nmethods, as they are at the same level of error in this regard.\nA discussion on the accuracy of B3LYP and cited benchmarks has been added\nto the manuscript:\nAs such, the expected level of error of TD-DFT/B3LYP, which is usually\nbetween 0.2 and 0.4 eV with respect to experimental\nbenchmarks[103,104], or theoretical best estimates (TBE) [105-108],\nwill be retained. The same is true for the oscillator strengths, which have\nreported MAE of between 0.05 and 0.15, compared to TBE [105,109].\nWe supplement these benchmarks by providing a comparison between\nADC(2) and TD-DFT/B3LYP vertical excitation energies in Fig. S16\nand Table S3 of the ESI. The DFT results correlate well with ADC(2),\nshowing a systematic underestimation of excitation energies, which\nincreases when moving to higher excited states.\nFurthermore, we believe that OMNI-P2x opens the possibility to use higher\nlevels of theory for complicated tasks via efficient fine-tuning, as demonstrated\non the example of nuclear ensemble approach spectra or nonadiabatic molecular\nPage 9 of 28\ndynamics, reducing the number of points that need to be labeled by orders of\nmagnitude.\nAdditionally, the authors claim: “Moreover, to the best of our knowledge, TD-DFT is the only truly “black-box” electronic-structure approach currently capable of providing excited-state data for a broad range of molecular sizes and compositions across chemical space without requiring case-by-case adjustments to the calculation protocol.“. This is not true. ADC(2), available in several software (ORCA, TURBOMOLE, Q-CHEM) is black-box and more accurate. Its error is well-known, it does not require any case-by-case adjustment, and it does not depend on the mathematical form of the functionals chosen. It is yes, more expensive, but I believe affordable for the pool of molecules used in this work.\nYou are right that beyond TD-DFT, there are other ‘black-box’ methods.\nThank you for suggesting the ADC(2) method. We agree with the assessment that\nADC(2) is a black box method that can provide more accurate predictions of\nspectra for organic molecules. However, ADC(2) is prohibitively expensive for\ndata generation at the scale presented in the manuscript with our resources: The\ngeneration of 1000 electronic spectra for the benchmark presented in this\nresponse letter required over 25000 CPU hours (using 12 AMD EPYC 9654\nCPU cores with the TURBOMOLE program, as suggested). By extrapolation,\ncovering the entire training set would require over 76 million CPU hours, almost\ntwo orders of magnitude higher than the cost of TD-DFT, which is\nunfortunately, not feasible for us.\nRegarding the UV/Vis spectra task, the conceptually inaccurate, strong limitation of being just a prediction of a collection of energies and oscillator strengths with poor correlation, with no actual information and correlation with states’ character or wavefunction or orbitals involved, remains after revision.\nWe hope that this revision, introducing elucidation of excitation characters from\nPage 10 of 28\nthe OMNI-P2x predictions, addresses your concerns. We also note that massive\nsimulations enabled by OMNI-P2x can serve for generating enough statistical\ndata, from which representative samples can be drawn for additional manual\nanalysis, or using post-processing tools such as Multwfn and TheoDORE\nmentioned below. Both our analysis tools, based on atom-wise contributions and\nwavefunction-analysis tools, are complementary. We added this discussion to\nthe manuscript.\nThe means of analysis presented above do not exclude the need for\ndedicated analysis, either manual or using wavefunction-based tools such\nas TheoDORE[84] or Multiwfn[85], which can provide additional insights\nbeyond spatial localization of the excitations, such as the orbitals'\ncharacter. These tools are complementary, as massive-scale simulations\nenabled by OMNI-P2x can serve as a first step for generating enough\nstatistical data, from where representative samples can be drawn, and post-\nprocessed further with the wavefunction-based tools.\nThe authors acknowledge this limitation and justify it in their response:\n“While this limits interpretability in terms of electronic structure, it is consistent with how UV/Vis spectra are typically constructed from ab initio calculations, where oscillator strengths are reported for adiabatic states without explicit orbital labeling. The utility of OMNI-P2x is therefore in reproducing excitation spectra (energies and intensities) efficiently and at scale, rather than in providing mechanistic orbital analysis.“\nWe have now amended the manuscript with the newly introduced interpretation\nof the character of excited states with OMNI-P2x.\nPage 11 of 28\nHere, I again disagree with the authors. Even if there is no explicit orbital labeling, an ab initio calculation includes molecular orbital and excited configuration coefficients that allow for determining the character of the excited states. Additionally, well-established post-processing codes like Multiwfn or TheoDORE can efficiently analyze wavefunctions and provide descriptors to train a machine learning model for excited states. This is a fundamental of computational spectroscopy, where the calculations are needed to explain what is behind the shape of the UV/Vis spectra. If this information is missing, the scope of computational chemistry is strongly diminished. That’s why I believe that a ground-breaking architecture should be able to properly reproduce the spectrum because it reproduces correct energy, oscillator strengths, and wavefunction information, not just the shape. If this model deliberately ignores this aspect, it not only reduces the impact of computational chemistry but also its reliability in terms of the physical interpretation of the results.\nWe agree that post-processing tools like Multiwfn and TheoDORE are\ncomplementary to our approaches, particularly to the newly introduced analysis\nof the atom-wise contributions with OMNI-P2x. Analysis with OMNI-P2x of\natom-wise contributions to excitations represents a significant step toward\ninterpretable ML models for excited states, providing chemically meaningful\ninsights about where excitations are localized. This addresses a key use case in\ncomputational spectroscopy: understanding which parts of a molecule are\ninvolved in electronic transitions, while maintaining the outstanding\ncomputational efficiency of OMNI-P2x. As mentioned above, both can be used\ntogether: Multiwfn and TheoDORE can be used for post-processing the single-\npoint calculations done with the electronic structure methods on a fraction of\nstructures generated from the OMNI-P2x simulations. An example would be\ntaking representative snapshots from ML-propagated nonadiabatic dynamics\nfor further analysis.\nPage 12 of 28\nThe authors write: “We note that, contrary to your comment, the S0→S1 excitation energies have consistently exhibited the lowest accuracy in both the original and revised models, so we are unsure where the claim of “the highest accuracy” originated. “\nI just wanted to clarify this point. If the character and order of the excited states are not the target of the prediction, then, in cases reported in Figure 5a,b,c, the prediction is quite acceptable. The shape and the position of the maxima are ok for a qualitatively low-cost prediction, so it is not a problem with the prediction of the lowest states.\nWe thank you for your clarification.\nMy comments about the NAMD example have been insufficiently addressed. I pointed out to the authors the reasons why this molecule was shown, implying a deeper physical analysis to justify the correctness of the simulation predicted. The trajectories are not analyzed to prove that the geometrical evolution of fulvene is well reproduced, passing through the correct conical intersections, and not only the time at which the hop should occur, or the energy gap. This is very important, especially because nonadiabatic couplings or wavefunction overlap are not predicted. I even suggested this point to the authors (since it is done in other implementations). I think this has not been properly addressed in the short paragraph in the discussion that was added.\nWe thank you for this suggestion; it is indeed one of the valuable ways to obtain\nadditional analysis of the OMNI-P2x simulations. We have added the analysis\nof the geometric evolution of fulvene to the revised manuscript as described in\nour reply to the next comment.\nPage 13 of 28\nIndeed, although the shape of the population resembles the CASSCF profile, there is no example of geometrical analysis, no check if the trajectory really passes through the correct conical intersection, or whether it’s just the interpolation that correctly predicts the hops. The populations alone are not a metric of accuracy of a NAMD simulation, as highly sensitive to the simulation parameters, and they are not a measure of correct physical behaviour, as the decay from the excited states could be derived by a wrong deactivation channel. This is particularly important to show, especially because in the previous task on azobenzene, it was shown that out of the equilibrium, the prediction gets significantly worse. Same for the excited state character, there is no analysis to show that there is no problem in this example. Since in NAMD by definition geometries out of the ground state equilibrium are explored, and the excited states’ character will change, this check is fundamental, in my opinion, to show the reliability of the model. It is missing the proof of the correct physical behaviour beyond just interpolation. Additionally, in figure 9, I believe that the y-axis should be the S_1 population..\nWe agree that state-averaged electronic state populations alone are not sufficient\nto judge the quality of ML-NAMD propagation.\nWe start by plotting the S and S PES, as simulated with the transfer learning\n0 1\nML model: panel (a) of Figure S13 of the ESI, and the CASSCF(6,6) reference,\npanel (b), along the two key reaction coordinates of fulvene C=CH2 stretching\nand rotation:\nPage 14 of 28\nThe agreement between both the ML models, and the reference PES is excellent,\nwith the fine-tuned model performing better further along the stretching\ncoordinate.\nAs a next step, we analyze the minimum energy conical intersection (one of the\nmany, at least 5, possible conical intersections of fulvene) optimized using both\nthe ML model and the reference QC method, which we present in Figure S14 of\nthe ESI.\nWe note an excellent agreement between both geometries, with an RMSD of\n0.03 Å.\nPage 15 of 28\nFinally, having confirmed that machine learning reproduces the correct PES\ntopology, and MECI structure for fulvene, we confirm that the deactivation\nproceeds through physically correct deactivation channels. To do that, we\nanalyze the hopping points in the reference CASSCF and ML-driven dynamics.\nA correlation plot between the C=CH2 distance and the mean C=CH2 dihedral\nangle at hopping geometries has been shown in Figure S15 of the ESI.\nThe distribution of the hopping geometries agrees excellently between the\nreference and ML-NAMD dynamics, showing the same pattern of early planar\nhopping, with later twisted-stretched or twisted-shrunk deactivation. The\nfractions of deactivation events proceeding through each channel have been\nquantified in Table S1, with ML-NAMD agreeing with the reference, within a\n95% confidence interval, as shown in Table S2.\nWe also thank the reviewer for noticing the typo in the y-axis of Figure 9 (which\nshould read S population).\nFurther analysis of the NAMD trajectories of fulvene, including\ncomparison of the ML-predicted potential energy surfaces, minimum\nenergy conical intersections, as well as hopping points and deactivation\nchannels, is provided in Section S4.1 of the ESI. In all cases, excellent\nagreement between the ML model and reference CASSCF calculations\nPage 16 of 28\nis noted.\nI appreciate the addition of a second experiment regarding the NAMD. However, the analysis is rather limited, and again, it is not shown, either in the SI, that the propagation that is delivered is correct. More importantly, also considering this example, I once again acknowledge the computational speed-up of the new active learning procedure for NAMD. But this is an improvement on how active learning is designed, but not a paradigm shift for how ML potentials are used in NAMD, since there is still a need for fine-tuning, and the model is not transferable for this task. Additionally, the lack of couplings, available in other architectures, does not push towards a physical improvement in the way MLIPs are used for NAMD.\nWe thank you for acknowledging our extension of the azobenzene transfer\nlearning experiment. We agree that the need for fine-tuning is a computational\nburden that is not removed by the usage of OMNI-P2x, however, the\nconstruction of such a model is also fundamentally difficult from a\ncomputational chemistry perspective. Methods most commonly used for\naccurate and trustworthy NAMD propagation, such as MRCI, CASSCF, or\nperturbation theory methods like CASPT2, require selecting a system-specific\nactive space, and hence cannot be directly replaced by a universal potential that\nhas not been fine-tuned. For this reason, we believe that deriving a procedure\nthat decreases the number of required points for stable and accurate NAMD\npropagation by a factor of 50 (in the case of azobenzene) is a groundbreaking\nmethodological advancement.\nSimilar to the fulvene dynamics, we provide further analysis of the ML-\npredicted PES, as well as the MECI and deactivation pathways. Starting with\nthe analysis of the PES, a scan around the C-N=N-C dihedral angle is presented\nin Figure S16 of the SI, which is the critical reaction coordinate for azobenzene\nphotoisomerization. We note a very good agreement between the ML and\nreference PES, with an MAE of 0.04 eV for S0 and 0.03 eV for S1.\nPage 17 of 28\nThe MECI structures optimized with both the ML model and AIQM1 are\npresented in Figure S17. The structures agree well and represent the same\nphysically correct “twisted-assisted inversion” mechanism.\nPage 18 of 28\nFurthermore, an analysis of hopping geometries, presented as a correlation plot\nbetween the maximum CNN angle and the C-N=N-C dihedral angle (which are\nthe two crucial reaction coordinates in this photoprocess), shows good\nagreement between ML-driven and reference NAMD trajectories, which is\npresented in Figure S18 of the ESI.\nPage 19 of 28\nThe manuscript has been adjusted to reference these results:\nAn analysis of azobenzene S0 and S1 PES', optimized MECIs, and\ndeactivation geometries is provided in Section S4.2, along with a\ncomparison to reference results, which agree excellently with the ML\nmodels.\nIn conclusion, while I agreed with the other reviewers about the advancement in terms of architecture, I also saw similar problems in the results shown. In my first revision, although I recognized that the model works, I underlined some important, potentially wrong interpretations that could be derived from the predictions of the model in computational chemistry applications. These have not been disproven, and the limitations have not been addressed, but remain after revision. Since I find them highly relevant, I still cannot recommend publication in Nature Communications, although I believe this version of the manuscript is improved.\nWe hope that this revision, based on your comments, has addressed your\nremaining reservations about this work.\nPage 20 of 28\nReviewer #3\nThe authors have given responses to each of my comments, and made substantial changes to the manuscript. Many of the\nconcerns I and other referees had are addressed in a fully satisfactory manner, but some of the most important are not, as\ndetailed below. Overall I would say the authors have demonstrated that the model shows some promise for applications in\ntheoretical spectroscopy and non-adiabatic dynamics, and represents an interesting methodological step forward, and the\nML model shows real promise.\nWe thank you for your evaluation that our work represents a methodological\nstep forward and that the ML model shows real promise.\nHowever, all reviewers have expressed a range of different concerns about its accuracy and the overselling of the\ncapabilities of the current model. These concerns add up to an uncertainty that the model exhibits the level of accuracy that\nwould be necessary for an out-of-the-box tool of genuine usefulness to practitioners not expert in analysis of the accuracy of\nML models. Given that it has been presented as a ready-to-use general-purpose tool I am therefore worried that publication\nin Nature Comms would be giving a seal of approval on the accuracy of the model that the current results do not support.\nWe definitely agree that Nature Communications is a journal aiming to present\ntools and methods that push beyond the current state of the art, and at the same\ntime, we believe this is exactly what our model is doing. We acknowledge that\nOMNI-P2x is far from perfect, but we clearly state the limitations in the\nmanuscript, as well as provide uncertainty estimates that warn the user and give\nan idea about its applicability. Furthermore, multiple solutions to mitigate the\ndisadvantages of OMNI-P2x, such as fine-tuning and active learning for\nNAMD, are presented in the manuscript. Hence, we believe the model can be\nuseful to practitioners who are becoming increasingly skilled at utilizing ML\ntools and are aware of the accuracy and applicability of MLIPs. The work also\nPage 21 of 28\nrepresents an important methodological advancement for specialists within the\nfield of ML-assisted (photo)chemistry that will be relevant for the development\nof future models.\nThis is particularly important given the strong statements about the accuracy, some of which the authors have not\nsignificantly toned down despite referees concerns, and which it is very important to avoid being exaggerated. Examples\ninclude “making ML approaches for excited-states a viable alternative to established QM methods” and “The model enables\nout-of-the-box predictions of electronic spectra of organic molecules, with accuracy approaching TD-DFT/B3LYP”. Uncritical\nuse of the model as a substitute for TD-DFT calculations would be extremely risky, given the current results. Therefore I\nwould not support publication of the manuscript in its current form.\nWe have revised the statements about OMNI-P2x accuracy to better reflect the\nactual results, as well as possible applications of the model. The first statement\nhas been revised to read:\n...making ML approaches for excited-states a viable alternative to\nestablished, cost-efficient QM methods for high-throughput calculations,\nsuch as semi-empirical methods.\nWhile the second has been amended to explicitly state the target and accuracy\nof the model:\nThe model enables out-of-the-box predictions of electronic spectra of\norganic molecules, targeting the TD-DFT/B3LYP level of theory, with\nan accuracy of about 0.15 eV. Importantly, OMNI-P2x offers higher-\nquality spectra predictions than commonly used semiempirical QM\nmethods for excited states, while having a substantially lower\ncomputational cost.\nPage 22 of 28\nNote that there is continued mixing of units even when 2 of the referees asked them not to do this as it is confusing for the\nreader. This unit mixing creates the appearance of deliberately obscuring the substantially worse performance of OMNI-P2x\nwith respect to ground-state universal models, which if I am reading it correctly, appears to be a very substantially higher\nerror (albeit with less training data). That is not consistent with statement that “they demonstrate that the ground-state error of\nOMNI-P2x exhibits decent performance, even with respect to the current state-of-the-art models”.\nWe apologize for the usage of multiple units. For other universal, ground-state\nmodels, we supply the error both in the originally reported unit, as well in\nkcal/mol. The exception was UMAs errors, which have been reported in\neV/atom, which has a physically different meaning (quantifying the error in\nprediction per atom of molecule vs error of the total energy). The closest\nconversion available here is kcal/(mol*atom), which has been added to the\nmanuscript. However, the same error metric was provided for OMNI-P2x\nground-state energies as well.\nWe want to assure the reviewer that there was absolutely no intent to obscure\nour model's performance. Indeed, the RMSE of 2.94 kcal/mol (compared to,\ne.g., ANI-1, with an RMSE of 1.3 kcal/mol), and MAE of 1.6 kcal/mol (vs\nAIMNet2 1.25 kcal/mol for neutral, singlet systems) is higher, but we believe it\nis reasonable, given our training set uses equilibrium geometries from\nPubChemQC rather than the extensive conformational sampling in dedicated\nground-state datasets. The performance gap primarily reflects these training\ndata limitations, rather than any fundamental architectural constraints.\nWe have revised the imprecise statement mentioned in the comment to read:\nThis ground-state MAE can be compared to the similar metrics of\nother MLIPs: ANI-1 [68] with an RMSE of 1.3 kcal/mol (OMNI-P2x 2.93\nkcal/mol), UMA [76] with around 1 to 3 meV/atom (0.023 to 0.069\nPage 23 of 28\nkcal/(mol·atom), depending on the version (OMNI-P2x 0.11\nkcal/(mol·atom)) and AIMNet2 [77] reporting a MAE of around 1.25\nkcal/mol for neutral compounds. Although those results are not\ncomparable one-to-one, as each of the models was trained and tested on\na different dataset, they show that the ground-state error of OMNI-P2x\nremains within a reasonable range from current state-of-the-art models,\ndespite not being the main focus of this study.\nI find the decision to replace the model entirely with an ensemble with reduced errors at this stage somewhat surprising – has\nthis substitution been made consistently across all the data in the paper, or just applied selectively to places where the referee\nconcerns were most serious?\nWe thank you for your comment and provide further clarification. The practice\nof creating an ensemble of models is commonly used to increase the accuracy\nand stability of ML predictions, and has been applied in MLIPs such as AIQM1\nor the ANI-1 family of models. Aside from the accuracy increase, this allows us\nto provide an uncertainty quantification, which, as mentioned further in this\nreview, is a significant improvement and would not be possible with just one\nmodel. The new ensemble has been used consistently throughout the manuscript\n(and, as a consequence, e.g., Figure 5 shows different systems and spectra than\nin the original version). For fine-tuning experiments, one model (or two copies\nof the same model in the case of azobenzene screening) was used. Information\nabout this has been added to the manuscript:\nWe have then fine-tuned two copies of one of the OMNI-P2x energy\nmodels on this set, using different training/validation splits...\nand:\nOne model from the energy ensemble of OMNI-P2x was used for fine-\ntuning in all NAMD applications.\nPage 24 of 28\nThe statement “The excitation energy errors are similar in magnitude to the target state energy errors and follow the same\ntrend” doesn’t really answer the original question of “should we expect errors on S0->Sx transition energies to exceed those\nof the contributing state?”\nThe excitation energies errors are expected to be slightly higher than the target\nstate (Sx) energy (about 0.2-0.3 kcal/mol). A correlation plot between the\nexcitation is presented below:\nFigure 4a has been improved and is now more useful in determining the spread of errors.\nWe are glad that our revision was satisfactory.\nPage 25 of 28\nAppropriate cautionary language has been inserted to address the problem of predicting as many as 10 excited states for small\nmolecules, and the fact that increasing success may be due to increasing Rydberg character.\nWe are glad that our revision was satisfactory.\nMinor presentational issues raised here have been acted on appropriately.\nWe are glad that our revision was satisfactory.\nThe main point of this comment was that the statement that the model approaches the accuracy of TDDFT is not justified\nby the data. The reason to use a ML surrogate model is to bypass the computational expense of a given high or moderate\nlevel of theory such as TDDFT or QC methods by replacing it with a model that gives an error compared its ground truth that\nis substantially lower than the error of the ground truth itself: otherwise one is simply compounding the ground truth error.\nThe authors response is simply that the model makes a further error (on top of the TDDFT error) that is comparable to typical\nerrors of TDDFT. I would say this is not acceptable in an ML surrogate model, and significantly undermines the possible uses\nof the model, and I would want to see the statement removed or revised substantially.\nWe agree with the reviewer that this statement might be potentially misleading\nand decided to remove it from the manuscript.\nThis section is more satisfactory now that the SPC spectra for the top 3 candidates is shown. These spectra show there is not\nparticularly good agreement between OMNI-P2x spectra and TDDFT, but the presence of low-energy redshifted absorption in\nboth does at least mean that the screening task of finding such molecules has been successfully achieved.\nWe are glad that our revision was satisfactory.\nPage 26 of 28\nThe new Section S4 of the ESI is very cryptic and does not add to the discussion. Taken at face value S12a simply compares\nLZBL with FSSH and shows they agree for the same PES, which is no great surprise and gives no information about the\nquality of the PES itself. Furthermore, Referee 2’s concerns about the accuracy of the dynamics do not appear to have been\nfully addressed in the rebuttal. However, the demonstrations now included in Figure 10 are fairly convincing evidence that the\nTL based on OMNI-P2x is effective in significantly reducing the amount of training data required to converge the model\npredictions with respect to training set size in such calculations.\nWe thank you for the comments. Indeed, the reason for the inclusion of Section\nS4 into the ESI was to compare the effect of a coupling-free surface hopping\nscheme (LZBL), with FSSH that includes nonadiabatic couplings, as Reviewer\n2 raised concerns about the accuracy of using a method that does not rely on\nnonadiabatic couplings: “This could be solved using another formulation of\nTSH, for which nonadiabatic coupling or wavefunction overlap would need to\nbe predicted by the model...”. We agree that further analysis of the ML-TSH\nresults is necessary, and provide it in section S4 of the manuscript, with a\ndetailed discussion in response to Reviewer’s #2 comment.\nThe fact that the revised model issues warnings when the predictions are highly uncertain (as assessed by committee standard\ndeviation predictions, presumably) is a significant improvement on the original approach.\nWe thank you for acknowledging the positive impact of including an uncertainty\nquantification, which was made possible by employing ensemble averaging.\nIndeed, the committee standard deviation is used, which has been stated in the\nmanuscript:\nAdditionally, OMNI-P2x provides an uncertainty quantification (UQ)\nof all predicted energies, using the standard deviations between\nquantities predicted by each of the three NN potentials as the error bar.\nPage 27 of 28\nA warning is issued to the user if the UQ of any excited states exceeds\na safe threshold, which is determined as the median plus three median\nabsolute deviations of the UQs of the test set, for each electronic state.\nPage 28 of 28\nResponses to the comments of the reviewers\nReviewer comments are in italics, and our responses are in regular, blue font.\nReviewer #2 (Remarks to the Author):\nI thank the authors for addressing my comments and improving the manuscript,\nespecially thanks to the new analysis reported.\nIn particular, the atom-wise energy decomposition and the geometrical analysis of\nfulvene dynamics make the interpretation of the results more reliable and robust.\nHowever, for the atom-wise energy composition analysis, a general metric to\nevaluate the accuracy would be more suitable than potentially cherry-picked examples.\nThe comparison between B3LYP and ADC(2) is also important, showing that,\nalthough correlated, the performances of TD-DFT are systematically worse than\nADC(2).\nI think the manuscript is highly improved through the two revision rounds, and I have\nno additional comments for the authors that prevent me from recommending\npublication.\nWe thank you for your final positive evaluation and previous reviews.\nReviewer #3 (Remarks to the Author):\nI do not have any new comments to add in relation to the latest revision.\nWe thank you for your previous reviews.\nReviewer #4:\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports.\nThis is part of the Nature Communications initiative to facilitate training in peer review\nand to provide appropriate recognition for Early Career Researchers who co-review\nmanuscripts.\nWe thank you for your contribution as a co-reviewer.\nPage 1 of 1"}, {"round": 2, "reviewer_comments": {"reviewer_2": "I thank the authors for addressing my comments and improving the manuscript, especially thanks to the new analysis reported.\nIn particular, the atom-wise energy decomposition and the geometrical analysis of fulvene dynamics make the interpretation of the results more reliable and robust. However, for the atom-wise energy composition analysis, a general metric to evaluate the accuracy would be more suitable than potentially cherry-picked examples. The comparison between B3LYP and ADC(2) is also important, showing that, although correlated, the performances of TD-DFT are systematically worse than ADC(2).\nI think the manuscript is highly improved through the two revision rounds, and I have no additional comments for the authors that prevent me from recommending publication.\n(Remarks on code availability)", "reviewer_3": "I do not have any new comments to add in relation to the latest revision.\n(Remarks on code availability)", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\n(Remarks on code availability)\nReviewed with supervisor."}, "author_response": ""}]

Decision letter