Spatial localization in nuclear spin-induced circular dichroism – a quadratic response function analysis †

The paper reports a new approach for calculating nuclear spin-induced circular dichroism (NSCD), an optical effect arising when nuclear spins are anisotropically oriented along the light beam. The results show that the strength of the NSCD signal of a particular nucleus at a given excitation wavelength reveals whether the nucleus belongs to an excited chromophore. NSCD could thus be used as an experimental probe for localization of excited states with spatial resolution of individual atoms. Nuclear magneto-optic (NMO) eﬀects are recently described phenomena originating from the interaction of light with local magnetic fields produced by nuclear spins. The phenomena border nuclear magnetic resonance and optical spectroscopy and are expected to provide rather unique spectroscopic features, borrowing from both localized response of the atomic nuclei as well as more global excitation properties of the whole molecule or its chromophore moieties. A number of quantum-chemical computational studies have been carried out, oﬀering a reasonable agreement with nuclear magneto-optics experiments performed so far. However, the detailed structure-spectra relation is still poorly understood. In this report we address the question of locality of one of the NMO eﬀects, namely nuclear spin-induced circular dichroism (NSCD). We implement an alternative computational approach for calculation of the NSCD intensities, based on residues of quadratic response functions, and use it to investigate the NSCD response of diﬀerent nuclei in a model molecular system with well-defined separate chromophores. The results show that significant NSCD at a given energy only occurs at the nuclei which are located in the chromophore that is excited. We rationalize these findings using analysis via diﬀerence densities, and approximate sum-over-states calculations. This behaviour of NSCD opens a way to experimental studies of localization of excited states in molecules, potentially with resolution down to the order of bond-length.


Introduction
Nuclear magneto-optic spectroscopy (NMOS) is an umbrella term for a group of spectroscopic effects that have been experimentally and theoretically investigated in recent years. In general, NMOS effects manifest themselves as changes in polarization state of the light upon interaction with a sample possessing a macroscopic nuclear magnetization. The induced optical effect depends on the mutual orientation of the light beam and sample magnetization vector and whether the molecule is in resonance with the light beam or the effect takes place in the dispersive energy region.
The first described, and so far the only experimentally observed, NMOS effect is nuclear spin-induced optical rotation (NSOR). [1][2][3][4][5] NSOR is the rotation of the plane of polarization of linearly polarized light induced by the nuclear spin magnetization oriented along the light beam. It can be thought of as a nuclear magneto-optic analogue to Faraday rotation, a classical magnetooptic effect caused by the external magnetic field parallel with the direction of the light beam. Following the experimental observation, a number of studies on the topic of NSOR have been published, 3,6-10 successfully rationalizing the measured NSOR values.
In addition to NSOR, several other NMOS effects have been theoretically predicted. These include three Cotton-Mouton-like effects: nuclear spin-induced Cotton-Mouton effect (NSCM), 11,12 nuclear spin-induced Cotton-Mouton effect in external magnetic field (NSCM-B), 13,14 and nuclear quadrupole-induced Cotton-Mouton effect (NQCM). 15,16 Most recently, nuclear spin-induced circular dichroism (NSCD) has been theoretically investigated. [17][18][19][20] NSCD is intimately connected to NSOR, stemming from the same physical origins. However, in contrast to NSOR, which is a dispersive effect, NSCD occurs at energies of electronic excitations, i.e., where the molecule has optical absorption bands. The relation between NSOR and NSCD can thus be compared to, e.g., that of optical rotation dispersion and circular dichroism in chiral compounds. Similarly, much like NSOR can be viewed as a nuclear magneto-optic analogue of Faraday rotation, NSCD can be seen as an NMO variation of magnetic circular dichroism.
It can be reasonably expected that NMOS in general should exhibit nucleus-specific features, owing to the presence of local interactions, such as the paramagnetic spin-orbit operator ĥ PSO , which appears in the non-relativistic description of NSOR and NSCD. 17 Such local properties have been briefly discussed before in the case of NSOR, which has been shown to provide a strong enhancement of signal near optical resonances, in a nucleusspecific manner. 9 However, available theoretical analyses have dealt with this topic only in a passing manner and the experimental measurements so far have not been performed near electronic resonances in order to confirm this prediction. Moreover, it has also been shown that the overall shape of the electron density and NSCD spectra seem to qualitatively correlate, 20 leading to the hypothesis that the analysis of the excited state densities may hold a key to new, nucleus-specific information hidden in NSCD.
To date, NSCD spectra have been obtained via the complex polarization propagator (CPP) approach 17,18,20 also known as damped response theory. [21][22][23] The CPP method has the advantage of taking into account contributions from all excited states at any calculated energy point, but it comes with the requirement of selecting beforehand a peak broadening parameter g, and the need to calculate many discrete points on a frequency grid to obtain a sufficiently smooth spectral curve. Moreover, in order to extract the NSCD contribution from a particular excited state, one needs to resort to numerical fitting of the CPP curve. 20 The results of such fitting are prone to errors due to the fact that even low energy regions of the CPP curve contain contributions from ''tailing'' peaks of excited states lying at higher energies.
The purpose of this study is twofold: First, we present an alternative method to CPP for calculating NSCD, based on residues of quadratic response functions, [24][25][26] to tackle the issue of obtaining analytical NSCD intensities for individual excited states. This approach allows us to calculate the NSCD intensity in a manner similar to calculations of the B term of magnetic circular dichroism. [25][26][27] Second, we use this alternative response function formalism as well as a sum-over-states approach, to investigate the NSCD response with respect to the location of the nuclei within the excited chromophore. We note in passing that even if both approaches to large extent build upon existing response function residue computational machinery, the specific formulas for NSCD were derived and coded for the first time for this work. The results show that the strongest NSCD signals come from nuclei in the regions of significant electron density changes during excitation. The NSCD can thus in practice offer an experimental tool to probe, nucleus by nucleus, the spatial extent of electronic excitations.

NSCD intensities from residues of quadratic response functions
The NSCD ellipticity (per unit sample length, spin polarization, and nuclear concentration) has earlier been shown 17 to be given by where o is the angular frequency of the incident polarized light, m 0 is the vacuum permeability, N A the Avogadro number, I K the nuclear spin quantum number and e abz is the Levi-Civita antisymmetric tensor, and implicit summation over repeated indices in implied (here and throughout). The quantity hhm a ;m b ,ĥ hf K,z ii g o,0 indicates the damped complex quadratic response function of the (Cartesian components a, b and z of) electric dipolel ¼ Àe P i r i and hyperfine interaction ĥ hf K operators. The latter corresponds, in non-relativistic theory and for a closed shell system, to the paramagnetic spin-orbit operatorĥ PSO ¼ e hm 0 4pm e g K P i l K;i r Ki 3 , with l K,i as the angular momentum of electron i about the position R K of nucleus K.
The expression in eqn (1) was derived by analogy with the one for the Magnetic Circular Dichroism (MCD) ellipticity (per unit sample length, concentration and magnetic field strength) 26,28 where the PSO operator is replaced with the orbital Zeeman (OZ) interaction operator component ĥ OZ z (responsible for Zeeman splitting). For a closed-shell system without degenerate excited states, the latter term in the MCD ellipticity corresponds to 29 where B(0f ) is known as ''the B term of MCD'', 27,30,31 for a transition from (ground) state 0 to excited state f, and a f (o,o f ) is a (Lorentzian) line-shape function with g as half-width at half-maximum (HWHM). A sum-over-state expression for B(0f ) was derived long ago by Stephens and Buckingham, 30-32 and it has previously been shown to be attainable from the single residues of (standard) quadratic response functions [24][25][26]29 Bð0 which can also be recast as the derivative of a (magneticfield perturbed) one-photon dipole transition strength, 29,33 The expression in eqn (5) allows one to straightforwardly decompose and rationalize the MCD spectral bands in terms of (signed) intensity sticks associated with specific electronic transitions 0f. To take advantage of a similar possibility for the NSCD signals, we introduce here the NSCD B term, B K (0f ), and compute it via residues of the regular quadratic response functions hhl;l,ĥ PSO K ii o,0 implemented in DALTON 34,35 as a generalization of the MCD computational protocols presented earlier. 25,36 The results were checked for consistency by comparing the computed (total) ellipticity spectrum of (active) nucleus K (per unit path length and unit of spin polarization) obtained as with the one obtained using the CPP formalism. 17 In other words, NSCD spectra were generated by Lorentzian broadening the computed stick spectra {o f , B K (0f )} using the same value of the (inverse) life-time parameter g used in the CPP calculations. The prefactor L in eqn (8) is a product of physical constants, see ESI. †

Computational details
All molecular structures were optimized using Turbomole 7.0.2 at B3LYP/def-TZVP/vacuum level with m5 grid 37,38 using available symmetry point groups (D 2h for ethene, C 2v for water and pyridine, D 6h and benzene, C s for the model system described below). The local minimum geometry was confirmed by calculations of the harmonic vibrational frequencies. The choice of method for calculation of excitation properties and NSCD depended on the purpose of the calculation, differing between benchmark tests and production calculations for investigations of the NSCD localization.

Theoretical models
2.3.1 Benchmark. To assess the consistency of the results with those obtained using the CPP protocols, the implementation was first tested on simple model molecules: water, ethene, pyridine and benzene. The latter was slightly distorted to avoid the occurence of degenerate states. The basis sets employed varied from 6-31G to aug-cc-pVTZ and the computational methods included Hartree-Fock and two density functionals -BHandHLYP and B3LYP. [39][40][41] A complete set of tests is given in Table 1.
The excitation energies and NSCD intensities for investigations of localization of the NSCD response in our model molecule (vide infra) were calculated using the BHandHLYP functional with the def2-SVPDD-0 basis set, 20 with a development version of DALTON. 34,35 2.3.2 The model molecule -PPT. The chosen model system for investigating the locality of the NSCD effect was 2-(3 0 -phenylpropyl)-1,3,5-triazine, labelled PPT throughout the following. The initial calculation of the excited states showed that the two isolated chromophoric units offer a number of excitations well-localized on either ring, suitable for testing the hypothesis of locality of NSCD. To facilitate the discussion of NSCD signals, the atoms were numbered as shown in Fig. 1.
The strength of the calculated NSCD signal of the individual nuclei has been assessed with respect to the local electronic density changes r diff, f upon excitation, defined as where the electron densities r f and r 0 of the excited state f and the ground state 0, respectively, were obtained in Turbomole 7.0.2 at the same level of theory as the corresponding NSCD.
The default Turbomole spatial grid for the electron densities was used -that is, a rectangular grid with spacing between points of approximately 0.1 Â 0.1 Â 0.1 Å. Molecular visualizations were prepared using the UCSF Chimera package. 42

Benchmarking the B-term formalism
As anticipated in the Methods section, the new implementation was tested for consistency via comparison of NSCD computed via eqn (7) and (8) with that obtained using the CPP protocol. The NSCD terms B K (0f ) for excited states f calculated via the B-term-like formalism have been broadened by the Lorentzian band of the same half-width g = 0.00456E H = 1000 cm À1 as used in the CPP calculations, according to eqn (8).
The results for ethene and pyridine at the HF/6-31G level are shown in Fig. 2. Clearly, the new B-term formalism allows to faithfully reproduce the CPP curve. Moreover, in this case of rather well separated excited states, the calculation is much more efficient than with CPP, resulting in significant reduction of computational time. It should be noted, however, that this advantage may not hold for systems with many closely spaced excited states or when a high-energy region of the spectrum is desired, since in that case many lower lying excited states would need to be evaluated, potentially making the B-term calculation too demanding. Moreover, the CPP approach is robust even in cases of degenerate states, where the residue approach may become numerically unstable. This is for instance the case for benzene in a fully symmetric D 6h geometry. For molecules belonging to the D 6h point group, degenerate electronic excited states occur, leading to divergent NSCD B-term values, on Table 1 Molecules, methods and basis sets used for benchmarking of the B-term approach to NSCD. Graphical plots are shown in Fig. 2, 3 (11) below). At the same time, degenerate excited states give rise to A-term type of signals, which can indeed be regarded as special occurrences of the B terms, stemming from the excitation to a degenerate excited state. 29 The signals due to the A terms can be pragmatically obtained as pseudo-A terms by computing the B-term contributions at a slightly distorted geometry, 43 see ESI. † Additional benchmark tests with similar level of agreement are shown in ESI. † As a second example, we have performed ad hoc comparisons of NSCD results for two carbon atoms (C1 and C7) of the PPT molecule (Fig. 3). It is apparent how inclusion of the B K terms from more excited states into the spectral ellipticity curve improves the convergence of the results towards the CPP result. It can be seen in the case of C7 that even excited states 13 to 17 contribute to the spectral curve in the lower energy region, showing the effect of ''tailing''.
Based on these results, we conclude that the B-term computational protocol provides correct NSCD intensities and signs and can be reliably used for the analysis presented in the following sections.

Localization of NSCD in molecules
The model for our investigation of locality of NSCD was the PPT molecule. The properties of its first nine excited states are summarized in Table 2 and their difference densities are shown in Fig. 4. As can be clearly seen, the first four excitations are localized in the triazine moiety, while the following three are in the phenyl ring. Higher excitations are delocalized over the  whole molecule. The numbers in Fig. 4 indicate the calculated strengths of the NSCD B terms, multiplied by 1000 for convenience.
A trend can be observed in Fig. 4. In the case of the first excited state, the carbon atoms of triazine show two orders of magnitude larger NSCD B terms than the ones in the phenyl group. In addition, the magnitude of the NSCD for the carbon atoms forming the bridge between the two chromophores decreases going from the triazine side towards the phenyl. A somewhat similar situation is observed in the case of the second excited state, with the remarkable exception of the carbon nucleus C9 in triazine, that shows negligible NSCD.
The third excited state also shows stronger signals in the excited chromophore, although the absolute numbers are smaller than in the previous cases, likely because of the small oscillator strength, which modulates the total NSCD intensity of the transition (see eqn (11) below). Likewise, the fourth state behaves in comparable manner, with the NSCD signal also significantly extending towards the adjacent bridge atom. On the other hand, excited states 5, 6 and 7, which are localized on the phenyl ring, manifest strong NSCD signals for the phenyl ring carbon nuclei and much smaller ones for the nuclei located in triazine. Excited states 8 and 9 are spread around the whole molecule and likewise is the overall NSCD intensity.
From these observations, a tentative conclusion can be made: the excited state location, determined through electron density changes, influences the magnitude of the NSCD response for the individual nuclei. In general, large NSCD signals appear in the regions of significant changes in electron density upon excitation. This behaviour is consistent with our previous observation, where atoms located outside of the chromophore showed significantly lower NSCD compared to those located within. 20 It should be noted, however, that large change in the local electron density is apparently a necessary, but not sufficient, condition for large NSCD, since in some individual cases, small NSCD could be observed even in such regions.

Rationalization of the NSCD relation to the electron density
It is worth discussing the underlying origin of the correlation between the local electron density and the NSCD magnitude.
The B K (0f ) has the SOS expression where E k indicates the energy of state k within a manifold of excited states, and all other quantities have been defined earlier. The SOS expression above was derived by analogy with the SOS expression for the MCD B term, and it can be considered at the basis of the quadratic response residue algorithm presented here. We have shown before 20 that it qualitatively reproduces the spectral shapes of exact CPP calculations. It is thus justified to use this expression to analyze and rationalize the relation between NSCD and the electron density. Starting from eqn (11), we can isolate the term k = 0 in the first summation, and the term k = f in the second one The last term connects the NSCD intensity for the excited state f to the difference in expectation value of the electric dipole moment operator (component) between the ground state 0 and the excited state f. Each expectation value can be written as contraction of the dipole moment integrals with the one particle density matrix of the given state, whereas the electron density for each state in eqn (9) is a contraction of the density matrix of that state with two orbitals. Thus the two quantities, density differences in eqn (9) and dipole moment differences in eqn (12), are related via the density matrix differences. Values of the last term in eqn (12) for the first nine excited states and all carbon atoms of PPT are plotted in Fig. 5.
As can be seen, there is a correspondence between the appearance of appreciable NSCD signals and the size of this ''density difference term''. However, it should be noted that the differences in dipole moments do not refer to any specific nucleus by themselves and thus can not be responsible for the differences between different atoms. That is modulated by the element h0|ĥ PSO K | f i, as we will demonstrate below. 3.3.1 Role of different matrix elements for total NSCD. Let us now elucidate in more detail how the local density difference plays such significant role in determining the NSCD intensity for specific nuclei, as seen in Fig. 4. Since ĥ PSO K is a local operator with r À3 dependence on the distance from the nucleus K, we will discuss the excitations in terms of changes happening ''near'' the nucleus. For this purpose, let us divide the molecule into two chromophores -triazine (T), and phenyl (P) and write the excited states as being composed of two separate regions of excitation where f T and f P can be either f or 0, depending on whether they are excited or not, respectively. For example, the first excited state can be expressed as |1i = |1 T ,0 P i since the phenyl moiety can be thought of as remaining in the ground state. Due to the local nature of the operator ĥ PSO K acting on nucleus K, we can then approximate the corresponding matrix elements taking into account only the electronic state of the chromophore to which the specific nucleus belongs. Thus, for example, integral h0|ĥ PSO C1 |1i for nucleus C1, which is in the phenyl chromophore, reduces to 0ĥ PSO where ''X'' indicates that the state of this chromophore is not taken into account by the operator, because of its remote position. Considering that hi|ĥ PSO K |ii = 0 (imaginary operator), and on account of the long distance to the excited chromophore, the ĥ PSO C1 is ''blind'' to the state of the more distant chromophore, and the matrix element tends to be relatively small in this case.
Based on this argument, integrals of the form hi|ĥ PSO K | ji can only be significant if the states i and j differ locally in the vicinity of nucleus K. From Fig. 4, we can see that, for example, states 0, 1, 2, 3 and 4 can not mutually produce significant integrals for nuclei in the phenyl group since they all share locally the same ground state density. The same argument applies to states 5, 6 and 7 for the case of triazine ring.
In short, this effectively means that if the nucleus K resides on a non-excited chromophore, the element h0|ĥ PSO K | f i is small and the state is ''locked out'' of interaction with several energetically close excited states, that all share locally unexcited ground-state wavefunction. On the other hand, the excited chromophore does not face such restrictions as excited state densities are, in general, different from both ground state, as well as from one another.

Relative significance of individual operators and their integrals
Expanding on the previous section, we can present yet another qualitative picture of the dependence of NSCD on the local properties in terms of the relative magnitudes of different molecular integrals appearing in the expression for the NSCD.
For this purpose, we interpret eqn (11) in literal sum-overstates terms. First, we need to assess the level of agreement between the B-term (QR-residue) and sum-over-states calculations. Fig. 6 shows the correlation of NSCD calculated using the two approaches for the first ten excited states and all twelve carbon atoms (120 points in total, 100 excited states used in the summation). As can be seen, the results are not completely converged even for 100 states, but they are steadily improving with the number of states (see ESI †). Nevertheless, there is a significant semi-quantitative agreement justifying our following discussion.
Þ (the term on the last line of eqn (12)) for the first nine excited states for all carbon atoms of the PPT molecule. Numbers in panes of excited states 5, 6, and 9 indicate multiplication factor for those values. We can expand and partition the NSCD expression in eqn (11) in the following way |fflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflffl ffl} |fflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflffl ffl} The first point to note is that the reciprocal energy prefactor on each line governs the overall weight of the individual contribution to NSCD. The energy differences between two excited states appearing on the first line can be arbitrarily small, and closely spaced excited states will thus in general contribute more than ones lying farther apart. On the other hand, the expression on the second line is weighted by the reciprocal value of the excited states energies. The latter are commonly sizable values, and their reciprocal value will then be (relatively) small. For example, in the current case of PPT, the lowest excitation energy is 5.11 eV, while the energy separation of the individual excited states is on the order of 0.1 eV. Thus, the terms on the first line have much larger overall weight.
Since a detailed analysis of the interplay of the molecular integrals based on full expansion of such expression becomes quickly rather cumbersome for more than a few states, it might be prudent to construct some auxiliary measure that would capture qualitative differences of various contributions. For this purpose, sums of magnitudes of four matrix elements appearing in eqn (12), i.e. linear and quadratic response function residues for both operators, have been calculated. The sums are summarized in Table 3 along with some of their basic properties that will be discussed. The table indicates the spatial nature of the matrix element, where elements with the ĥ PSO K operator are deemed ''local'' due to their r À3 dependence on the distance from their respective nuclei, and elements with electric dipole moment operatorl are described as ''global'' to reflect their independence on the specific nucleus.
The choice of the sums being composed of these four types of molecular integrals is motivated by the following: (a) they are four naturally occurring components in the expression for NSCD (eqn (11)); (b) they have different spatial extent (''local'' vs. ''global''); (c) they differ in the extend to which their values change when the final excited state f changes fg, since they stem from either linear or quadratic response functions. Summing together their absolute magnitudes thus gives some qualitative insight into the relative significance of that particular type of integral in determining of the NSCD values for different excited states and nuclei.
We stress again that these sums are only a qualitative tool that does not correlate directly with the NSCD intensity, as the latter contains the mixed products of such matrix elements. Nevertheless, the sums can be used as auxiliary measures of the variability of the average total magnitude of individual elements appearing in the expression for NSCD intensity, with respect to different excited states f and nuclei K. The properties of some of these sums can be inferred at a glance, whereas for others they are more involved and we thus present them graphically in Fig. 7. Let us inspect them.
The elements h0|l|ki form a sum that is obviously independent of the nucleus and is only slightly dependent on the final excited state f. The dependence arises because sums are composed of (n À 1) identical terms, except for a single element h0|l| f i, which is excluded in each case. In comparison, elements of the form hk|l| f i sum to a quantity that is, predictably, also nucleus-independent, but depends strongly on the final state f. This is due to the fact that all (n À 1) elements hk|l| f i in the summation are different for different states f, with the exception of the pair hk|l| f i = h f |l|ki which naturally appears also when roles of states k and f in the summation are exchanged.
Concerning the ĥ PSO K operator, it can be seen that the is dependent on the final excited state f only through a single element (as discussed above for the dipole moment operator), but predictably depends considerably on the nucleus. This is expected as the operator is nucleusspecific. Finally, the sums of elements of the form hk|ĥ PSO K | f i exhibit the largest variation in magnitude both between excited states as well as between individual nuclei. This is a consequence of ĥ PSO K being a nucleus-specific operator and every single element in the sum changing for different final state f.
Remarkably, it can also be seen that the variation in intensity of this sum somewhat correlates with the intensity of NSCD for a particular combination of nucleus and excited state. It should be pointed out again, that since the elements hk|ĥ PSO K | f i are combined with other contributions, the sum does not correlate with the actual NSCD intensity, but merely gives the qualitative information about the relative importance of the excited states at the site of the probed nucleus. Of particular interest here is excited state number 2 for nucleus C9. As it was pointed out in the discussion of Fig. 4, this nucleus resides inside the volume of large density change, but it exhibits very weak NSCD. As it can be seen here, the magnitudes of the h2|ĥ PSO C2 |ki elements are significant. However, due to their combination with other smaller integrals in the sum and variation in signs, the net result is accidentally almost vanishing in this particular case. The results suggest that the integrals hk|ĥ PSO K | f i play a significant role in determining the total NSCD intensity due to their sensitivity to both to excited states as well as different nuclei. This effect appears to be a necessary, if not sufficient, condition for appreciable NSCD and its impact is potentially further enhanced by the large reciprocal energy weighting prefactor.

Conclusions
We have presented a new computational protocol for calculating excited-state specific, B-term-like, contributions to the NSCD ellipticity. The approach is based on single residues of quadratic response functions involving two electric dipole and one PSO operators and it allows to analytically evaluate the NSCD for a given electronic transition. Benchmark tests have shown that the values thus obtained yield NSCD spectra identical to the ones calculated via a previously established complex polarization propagator approach.
Using the new method, we have investigated the question of local nature of the excitation in NSCD and its relation to the signal from the nuclei residing in the chromophore being excited. The results show that the strong NSCD signals form only for nuclei residing in the excited chromophores, defined as a place of significant change in electronic density. The findings were rationalized based on an analysis of the individual terms contributing to the NSCD intensities, and found to originate from the ĥ PSO K operator due to its localized nature. In particular, second order perturbations involving two excited states appear to be important for correct description of this phenomenon.
The presented results show that NSCD has the ability to distinguish between nuclei based on their location in different chromophores. This opens up interesting possibilities for a two-dimensional nuclear-resonance/optical-resonance spectroscopy for probing localization of excitation states in molecules. The results also establish a basis for a rational search for suitable molecular systems to target in first observations of this exciting phenomenon.

Conflicts of interest
There are no conflicts to declare.

Acknowledgements
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement NMOSPEC, No. 654967 (P. Š.) and ETN COSINE, No. 765739 (S. C.). We also thank the Magnus Ehrnrooth Foundation for financial support (P. Š.). The authors acknowledge financial support from the Kvantum institute (University of Oulu) and from Academy of Finland (Grant 316180) (P. Š.). We acknowledge grants of computer capacity from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533). S. C. acknowledges support from DTU Chemistry and from the Independent Research Fund Denmark -DFF-Forskningsprojekt2 grant no. 7014-00258B.