7.45. More on the Excited State Dynamics module¶

ORCA has now a module designed to calculate properties related to excited states named ORCA_ESD. It can be used to predict absorption/emission spectra, transition rates, resonant Raman, and MCD spectra, based on a path integral approach to the dynamic process [199]. It has some of the functionalities of ORCA_ASA and even more, as it will be discussed. What we do here is NOT a conventional dynamics with trajectories along time points, we rather solve the equation for the transition rates or intensities depending on the different cases considered.

This formulation works because there is an analytic solution to the path integral of the Multidimensional Harmonic Oscillator and the assumption of Harmonic nuclear movement is critical. In many cases that approximation does hold and the results are in very good agreement with the experiment. The general usage of the ORCA_ESD module and some examples are already presented on Sec. Excited State Dynamics and it is recommended to read that before going into the details here. We now will discuss the specifics and keywords related to of each part of the module. A complete keyword list can be found at the end of this section.

7.45.1. Absorption and Emission Rates and Spectrum¶

7.45.1.1. General Aspects of the Theory¶

The idea behind calculating the absorption or emission rates starts with the equation from the quantization of the electromagnetic field for the transition rates between and initial and a final state:

(7.275)¶\[ k(\omega)_{if} = \frac{4\omega^3 n^2}{3\hbar c^3} | \langle \Psi_i | \widehat{\mu} | \Psi_f \rangle |^2 \delta(E_i - E_f \pm \hbar\omega) \]

with \(\hbar\omega\) being the energy of the photon, \(\widehat\mu\) the dipole operator and \(n\) the refractive index of the solvent, as suggested by Strickler and Berg [830].

One way to obtain \(k(\omega)\) is to compute it in the frequency domain, by calculating the Franck-Condon Factors between all initial and final states that satisfy the Dirac delta in Eq. (7.275), considering the thermally accessible initial states with the appropriate weight,

(7.276)¶\[ k_{}^{obs} = \int k(\omega) d\omega, \qquad k(\omega) = \sum_{if}P_i(T)k_{if}(\omega), \]

where \(P_i(T) = e^{-\frac{\epsilon_i}{k_BT} }/Z\) is the Boltzmann population of a given initial state at temperature \(T\), \(\epsilon_i\) is the total vibrational energy of state \(i\) and \(Z\) is the vibrational partition function. However, this can lead to a very large number of states to be included, particularly if there are low frequency modes. In this work we will take the a different approach and switch to the time domain, by using the Fourier Transform representation of the Dirac delta,

(7.277)¶\[ \delta(\omega) = \frac{1}{2\pi} \int_{-\infty}^{+\infty} e^{i\omega t} dt, \]

so that the equation to solve, in atomic units, is:

(7.278)¶\[ k(\omega) = \frac{2\omega^3}{3\pi c^3Z} \sum_{if}e^{-\frac{\epsilon_i}{k_BT} } \langle\Theta_i|\vec\mu^e |\bar{\Theta}_f \rangle \langle\bar{\Theta}_f|\vec\mu^e |\Theta_i \rangle \nonumber \int e^{i(E_i-E_f - \omega)t} dt, \]

with \(\vec\mu^e\) being the “electronic transition dipole” and \(|\Theta \rangle\) the vibrational wavefunction of the initial or final state.

After some extra steps, redefinition of the time variable and insertion of a resolution of identity, it can be shown that this equation is ultimately simplified to a Discrete Fourier Transform (DFT) of a correlation function \(\chi(t)\) with a timestep \(\Delta t\), multiplied by a prefactor \(\alpha\) [199]:

(7.279)¶\[\begin{split}\begin{aligned} k(\omega) &= \alpha \int Tr(\vec\mu^e e^{-i\widehat{H}\tau}\vec{\mu^e} e^{-i\widehat{\overline{H} }\bar{\tau} })e^{i \Delta E t}e^{-i\omega t} dt\nonumber \\ &= \alpha \int \chi(t)e^{-i \omega t} dt\nonumber\\ &= 2\alpha \hspace{5pt} \mathfrak{Re} \int_0^{\infty} \chi(t)e^{-i\omega t} dt\nonumber\\ &\simeq 2\alpha \Delta t \hspace{5pt} \mathfrak{Re} \hspace{5pt} DFT\{\chi(t)\},\end{aligned} \end{split}\]

and this correlation function is then calculated using path integrals analytically at each time point \(t\).

If one considers that the electronic part of the transition dipole varies with nuclear displacements and we allow for it to depend on the normal coordinates (\(\mathbf{Q}\)), such as:

(7.280)¶\[ \vec\mu^e(\mathbf{{Q} }) = \vec\mu_0^e +\sum_{i} \left.\frac{\partial{\vec\mu^e} }{\partial{Q_i} }\right|_{\mathbf{Q}=0} Q_i + \ldots, \]

we can even include vibronic coupling or the so-called Herzberg-Teller (HT) effect. The Frank-Condon (FC) approximation keeps only the coordinate-independent term. The correlation function for the HT cases can then be derived recursively from the FC one and the calculation is done quite efficiently. It is important to say the one must choose ONE set of coordinates in order to expand the transition dipole. In our formulation, it is always that of the FINAL state and that has some implications discussed below.

Other important aspect of this theory is that, in order to solve the path integrals, one has to work in one set of coordinates, the initial (\(\mathbf{Q}\)) or the final state ones (\(\mathbf{\bar{Q} }\), with the bar indicating final coordinates). As we have a transition matrix element, one set of coordinates have to be transformed into the other and it is easy to show that they are related by

(7.281)¶\[ \mathbf{Q} = \mathbf{J\bar{Q} } + \mathbf{K}. \]

That was first proposed by Duschinsky in the late 1930s[228] with the Duschinsky rotation matrix \(\mathbf{J}\) and the displacement vector \(\mathbf{K}\) defined as

\[\mathbf{J} = \mathbf{L}_x^{\text{T} }\mathbf{\bar{L} }_x, \qquad \mathbf{K}=\mathbf{L}_x^{\text{T} }(\mathbf{\bar{q}_0}-\mathbf{q_0}),\]

with \(\mathbf{L}_x\) being the matrix containing the normal modes, here described in Cartesian coordinates (\(x\)), and \(\mathbf{q_0}\) begin mass weighted coordinates (\(q_i = \sqrt{m_i}x_i\)).

The program runs by first reading and obtaining the initial and final state geometries and Hessians, then computes the Duschinsky rotation matrix and displacement vector, calculates the derivatives for the transition dipoles and computes the correlation function. After that, the DFT is done and the rates are obtained and printed when necessary. As the intensities observed experimentally are proportional to the rates, the spectrum is also calculated and printed on a BASENAME.spectrum file. If PRINTLEVEL HIGH is requested under %ESD, the correlation function is also printed on a BASENAME.corrfunc file.

OBS: The units for the Emission spectra are rather arbitrary, but for Absorption they are the experimental “molar absorptivity (\(\varepsilon\))” in L mol\(^{-1}\)cm\(^{-1}\) [199]. Be aware that these are dependent on the line width of the curves.

7.45.1.2. Approximations to the excited state PES¶

As already mentioned in Sec. Excited State Dynamics, in order to predict the rates we need at least a ground state (GS) and an excited state (ES) geometry and Hessian. In ORCA, we have seven different ways to approximate this ES PES: AHAS, VH, VG, HHBS, HHAS, UFBS and UFAS. Those can can be choosen by setting the HESSFLAG under %ESD. If you actually optimize the ES geometry and input the Hessian, that will be called an Adiabatic Hessian (AH) method an no keyword must be given on the input.

OBS: In the present version, these approximations are only available for Absorption, Fluorescence and resonant Raman. They cannot be directly used for phosphorescence and ISC rate calculations. However, one can do the latter calculations by generating approximate ES Hessians from “fake” fluorescence or absorption runs, and then using the ES Hessians in AH calculations (vide infra).

The idea behind these approximations is to do a geometry update step (\(\mathbf{\Delta S = -gH^{-1} }\) for Quasi-Newton and \(\mathbf{\Delta S = -g(H+S)^{-1} }\) for Augmented Hessian) to obtain the ES structure and somehow approximate the ES Hessian. The gradient (\(\mathbf{g}\)) and Hessian (\(\mathbf{H}\)) used on the step are on column Step of Table Table 7.26 below, with a description of the final ES Hessian:

Table 7.26 Methods used to estimate the ES PES¶
Method	Step	ES Hessian
AHAS	ES grad + GS Hessian	calculated on the ES geometry
VH	ES grad + ES Hessian at GS geometry	calculated on the GS geometry
VG (default)	ES grad + GS Hessian	equal to GS Hessian
VGFC	ES grad + GS Hessian	equal to GS Hessian (+ APPROXADEN TRUE)
HHBS	ES grad + Hybrid ES Hessian on GS geometry	Hybrid Hessian on GS geometry
HHAS	ES grad + GS Hessian	Hybrid Hessian on ES geometry
UFBS	ES grad + Updated frequencies ES Hessian on GS geometry	Updated frequencies ES Hessian on GS geometry
UFAS	ES grad + GS Hessian	Updated frequencies ES Hessian on ES geometry

OBS: Always use the GS geometry on the input file, equal to the one in the GSHESSIAN! If you asked for OPT FREQ at the input, a .xyz file is generated with the same geometry found on the .hess. If you picked the geometry from the .hess file, remember that it is in atomic units, so you have to use BOHRS on the main input.

After the calculation of the ES PES, a file named BASENAME.ES.hess is printed and can be used in future calculations. For example, one can use it in a separate AH calculation, as if the file contains the optimized ES geometry and the exact ES Hessian (of course, in reality both are approximate). This allows one to perform e.g. phosphorescence and ISC rate calculations with approximate ES PESs even though ESD does not support these calculations directly. For example, the phosphorescence rate from a triplet state to the \(S_0\) state, with the VG approximation for the triplet PES and the ES gradient computed at the TDDFT level, can be computed as follows:

Run an ESD fluorescence calculation using the VG approximation, but where “IROOTMULT TRIPLET” is added to the %TDDFT block. The resulting rate and spectrum are not meaningful, but the BASENAME.ES.hess file generated from this calculation is the correct VG approximation to the triplet PES.
Run an ESD phosphorescence calculation using the AH method, with the aforementioned BASENAME.ES.hess file as the TSHESSIAN.

This way, the first ESD calculation only uses the information of the scalar triplet state in computing the VG step, while the second ESD calculation only uses the SOC-corrected states. The same ES Hessian file can therefore be used in the calculations of all three spin sublevels of the triplet state.

In addition to the BASENAME.ES.hess file, if there was any updates on the GS Hessian, like transition dipole derivatives, a BASENAME.GS.hess is also printed.

The step can be controlled with the GEOMSTEP flag, with QN or AUGHESS options.
Currently all ES Hessians are calculated numerically, if you want to change the parameters related to the frequency calculations, you can do that under %FREQ (Sec. Vibrational Frequencies). The numerical gradient settings are under %NUMGRAD (Sec. Numerical Gradients).
The ES hybrid Hessian is calculated in the same way as described in Sec. Hybrid Hessian, except that the “model” Hessian is the GS one.
The ES Hessian with updated frequencies is recalculated from the same GS normal modes, after an update on the frequencies, as \(\mathbf{H}_{up} = \mathbf{L\omega}^2_{up}\mathbf{L}^T\). With \(\mathbf{L}\) being the normal modes and \(\omega_{up}\) the updated frequencies, with negative sign being kept after the square.
The frequencies are updated depending on a calculation of the energy after a given step. If the ES modes are equal to the GS, then a step over a coordinate \(\delta q_i\) that would result in an energy difference \(\delta E\) is given by \(\delta q_i = (-\frac{g_i}{\sqrt{\omega_i} } + \frac{\sqrt{g_i^2} }{\omega_i} + 2\delta E \omega_i)/\omega_i\). The default \(\delta E\) used is \(10^{-4}\) Eh, in general above the error of the methods. If the error in energy after the step is larger than a threshold given by the UPDATEFREQERR flag (default 0.20 or 20%), the gradients are calculated and the frequency recomputed. If not not, that mode and frequency are assumed to be the same.
The Updated Frequencies method can greatly accelerate the calculation of the Hessian, for much fewer gradient calculation are necessary, although you do not correct the modes. Also, the derivatives over the modes are already computed simultaneously.
The expected energy error \(\delta E\) can be changed using the UF_DELE flag.
The default method is the VG, but the AHAS is more trustworthy for unknown systems, although a lot heavier (Sec. A better model, Adiabatic Hessian After a Step (AHAS) and [199]).
Always check the sum of \(K_i^2\) printed on the output. If that number is too high (above 8 or so), it means that the geometries are too displaced and the theory might not work on these cases (check for different coordinate systems then, Sec. Normal modes coordinate systems).
Also check for RMSD between the geometries after a step. If the difference is too big, there might be problems with the step.

7.45.1.3. Mixing methods¶

In principle, it is possible to use different methods to compute different parts needed for the ORCA_ESD module. You could, for instance use (TD)DFT analytical gradients for the ground/excited state geometries and Hessians and a more elaborate method such as STEOM-CCSD to get the energies and transition dipoles. If you want to do that, just input the Hessians and use the DELE flag for the energy difference between the states - at their own geometry! - and TDIP x,y,z to input the transition dipole. If there is SOC and the transition dipole is complex, use TDIP x.real, x.imag, y.real, y.imag, z.real, z.imag. The program will automatically detect each case. If you don’t input these, they will be obtained by the module during the run, so you can set the excited method you want and let the program figure out DELE and TDIP.

OBS: It is not advisable to mix different levels of theory during a geometry step though. If you obtained a GS Hessian from DFT, doing a step based on a CASSCF ES numerical gradient might not lead to reasonable results. The same would be problematic even for different DFT functionals.

7.45.1.4. Removal of frequencies¶

If, after calculating an ES Hessian you end up with negative frequencies, the calculation of the correlation function might run into trouble. The default for the module is to turn all negative frequencies positive, printing a warning if any of them was lower than -300 cm\(^{-1}\). If that is the case, you are probably on a saddle point and not even a minimum, so results should be taken with care,

You can also choose to completely remove the negative frequencies (and the corresponding from the GS), by setting IFREQFLAG REMOVE or leave them as negative with IFREQFLAG LEAVE under %ESD.

Sometimes, low frequencies have displacements that are just too large (check on the \(\mathbf{K}\) vector), or the experimental low frequency modes are too anharmonic and you might want to remove them. It is also possible to do that setting the TCUTFREQ flag (in cm\(^{-1}\)), and all frequencies below the given threshold will be removed.

7.45.1.5. Normal modes coordinate systems¶

When working with systems with weak bonds, such as hydrogen bonds and \(\pi\) stacking, or with biphenyls and similar planar molecules it is common that there will be low frequency-high amplitude modes related to the angular bending. It has been shown that, in some cases, the normal modes transformed from Cartesian coordinates might be not sufficient to describe systems with these large amplitude motion [149]. In those, the definition of normal modes in terms of some curvilinear set of coordinates such as the internal ones are more suitable.

The first order transformation from Cartesian (\(\mathbf{x}\)) to internal (\(\mathbf{s}\)) coordinates is given by Wilson’s \(\mathbf{B}\) matrix[426] as

(7.282)¶\[ \mathbf{s} = \mathbf{B(x-x_0) }, \]

and here we use Baker’s[68] delocalized internal coordinates as default. First, a redudant set is build and an singular value decomposition of the \(\mathbf{G} = \mathbf{BB}^T\) matrix is performed to obtain the non-redundant set. The latter can be generated by \(\mathbf{B' = U}^T\mathbf{B}\), where \(\mathbf{U}\) are the eigenvectors correponding to non-zero eigenvalues of \(\mathbf{G}\). Then an orthogonal set is contructed from \(\mathbf{B'' = G'}^{-1/2}\mathbf{U}^T\mathbf{B}\). As the eigenvectors are not conitnuous functions of the coordinates, in order to avoid using an arbitrary selection, we will follow the work of Reimers[714] and set \(\mathbf{G}^{-1/2}\mathbf{U}^T = \mathbf{\bar{G} }^{-1/2}\mathbf{\bar{U} }^T\), or use the same transformation for the initial an final coordinates. Please note that this may lead to numbers larger than 1 on the Duschinsky rotation matrix, for it is an approximation and the calculated rates may vary a little. The normal modes in internal coordinates (\(\mathbf{L}_s\)) are then obtained from those in Cartesian ones (\(\mathbf{L}_x\)) as

(7.283)¶\[ \mathbf{L}_s = \mathbf{B''M}^{1/2}\mathbf{L}_x, \]

and the Duschinsky relation ((7.281)) still holds[149], with the displacement vector being simply

(7.284)¶\[ \mathbf{K}_s = \mathbf{L}_s^T(\mathbf{\bar{s} } - \mathbf{s}). \]

The options available for coordinate systems can be set under COORDSYS, and can be CARTESIAN, INTERNAL (for Baker delocalized - default), WINT (for weighted internals following Swart and Bickelhaupt [837]), FCWL (force constant weighted following Lindh [526]) and FCWS (same as before, but using Swart’s model Hessian).

OBS: Calculating in internal coordinates is usually better but not necessarily. If something goes wrong, you may also want to try the Cartesian option.

7.45.1.6. Geometry rotation and Duschinsky matrices¶

The electronic transition should not to take place simultaneously with translations and rotations[746] of the molecular structure. Before further calculations take place, the initial and final state structures are superimposed to satisfy Eckart’s conditions by obtaining a rotation matrix \(\mathbb{B}\) that ensures \(\sum m_i(\mathbb{B}\mathbf{R_i} \times \mathbf{\bar{R}_i})=0\) [239], with \(\mathbf{R}\) being Cartesian coordinates. As the initial geometry is rotated, so must be the corresponding normal modes \(\mathbf{L}_x\) also. This can be turned of by setting the flag USEB FALSE.

By the default the Duschinsky rotation matrix is set to Identity, to take advantage of our much faster algorithm. To turn that on, just set USEJ TRUE. You can check the “diagonality” of this matrix by looking at the Diagonality Index (DI), here defined as \(\sqrt{\sum_i \mathbf{J}_{ii}^2 / \sum_{ij}\mathbf{J}_{ij}^2}\). A DI=1 would be a perfectly diagonal matrix. The amount of mixing between modes is rpresented by the Mixing Index, with is here is given by \(\langle |J_{max}| \rangle\), or the average value of the maximum \(J_i\) of each line. If DI=1, it means each normal coordinate from the initial state is equal to a mode of the final state. When USEJ=TRUE, the largest components of the \(\mathbf{J}\) matrix are printed along with the \(\mathbf{K}\) vector, so you can have a better idea of how the mixing is occuring.

7.45.1.7. Derivatives of the transition dipole¶

The derivatives of the transition dipoles with respect to the normal coordinates of the final state can be obtained directly from the derivatives with respect to the Cartesian coordinates as

(7.285)¶\[ \mathbf{U}(\mathbf{\bar{Q} }) = \mathbf{\bar{L} }_x^T \mathbf{M}^{-1/2} \mathbf{U(\bar{R} }), \]

\(\mathbf{U}\) being the matrix of the x,y and z components of the derivative, \(\mathbf{M}\) a \(3N\times 3N\) matrix with the atomic masses along the diagonal. Also, in case one already has the derivatives with respect to the initial state , those can be transformed into the derivatives with respect to the final state by using the Duschinsky relation, assuming that \(\vec\mu^e_0(\mathbf{\bar{Q} }) + \sum_i \frac{\partial \vec\mu^e}{\partial \bar{Q}_i}\bar{Q}_i = \vec\mu^e_0(\mathbf{Q}) + \sum_i \frac{\partial \vec\mu^e}{\partial Q_i}Q_i\), so that

(7.286)¶\[ \frac{\partial \vec\mu^e}{ \partial \bar{Q}_k} = \sum_j \mathbf J_{jk}\frac{\partial \vec\mu^e}{\partial Q_j}. \]

By default, this transformation is NOT done, since Eq. (7.286) is an approximation. If you want to turn it on, set CONVDER TRUE under %ESD.

OBS: Remember that, if you already have the Cartesian derivatives over the final state, like if you use AHAS for an absorption spectrum, the conversion should be exact (although there might be numerical issues, always use a larger GRAD for frequencies!).

Alternatively, these can be calculated numerically from displacements over each normal mode. In this case, it is convenient to use the dimensionless normal coordinates \(q_i = \omega_i^{1/2}Q_i\) which represent proportional displacements on the PES [677]. We use \(\Delta q = 0.01\) by default, but this can be changed setting the DER_DELQ flag.

Again, DO NOT MIX different coordinates systems. If the derivatives were calculated over one coordinate set and you decide to change it, it has to be recalculated. You can manually delete them from the BASENAME.ES.hess file.
For hybrid functionals, you can choose to use DFT for the gradient, energy and transition dipole, and the fast simplified TDA (Sec. Simplified TDA and TD-DFT) only for the derivatives by seting STDA TRUE under %ESD.
A simple trick can be used to accelerate the computation of derivatives. If the first displacement gives a transition dipole that is too close to the reference, then the derivative can be assumed to be small and just the plus displacement may be taken to compute the derivative (with an usually small error). If it is large enough, then the minus displacement is also done and central differences is used. This is the default method and can be turned off by setting FASTDER to FALSE.
The central differences option can be altogether turned off by setting CENTRALDIFF FALSE under %ESD.
If you are having problems, set a larger PRINTLEVEL to check the individual calculation of the derivatives, you might be having some kind of root flipping during the displacement, or some other issue.

7.45.1.8. The Fourier Transform step¶

After the calculation of the correlation function, it is necessary to do a Fourier Transform (FT) step. To do that numerically, it is needed to correctly choose the grid in wich the time points will be computed, for that affects how the results will be obtained in the frequency domain. We have developed a method to generate an optimal set of parameters, depending on the final spectral resolution desired [199] and it will be used by default. Even so, you can choose your own grid by setting the NPOINTS and MAXTIME (in atomic units!) flags under %ESD. There are a few comments related to that:

Because we solve the FT using a very efficient Cooley-Tukey algorithm, the NPOINTS should be always multiple of two. You can put any number on the input, but the next larger multiple of two will be calculated and set.
The MAXTIME should be enough so that the correlation function goes to zero. If anything goes wrong, please check the BASENAME.corrfunc file for that.
The finer the spectral resolution chosen with SPECRES, the largest MAXTIME must be.
If you have a larger MAXTIME, you also must increase NPOINTS, otherwise the grid will be too sparse and many oscillations will be skipped.

7.45.1.9. Spectrum options¶

The final spectrum is saved in a BASENAME.spectrum file, with the total spectrum, the FC and HT parts discriminated, as explained in Sec. Excited State Dynamics. He are some detailed about it:

The range for which the spectrum is saved is given by default, but it can be set using SPECRANGE flag under %ESD, as SPECRANGE 10000,70000.
All of the INPUT units should always be in CM-1, but you can choose the OUTPUT units by setting the UNIT flag to CM-1, NM or EV.
In order to better converge the correlation function and approximate experimental spectra, a lineshape function can be used instead of the delta. The default is to use a LOREZENTIAN lineshape, but LINES can be set to DELTA, LORENTZ, GAUSS or VOIGT.
The DELTA lineshape might lead to a correlation function that oscillates forever, so please take care with that option.
The default line widths are LINEW 50 and INLINEW 250.
If you use a VOIGT lineshape, the Gaussian width can be controlled separately using the INLINEW flag. By default, it will be proportial to the Lorezntian to reach the same FWHM.
The LINEW and INLINEW are NOT the full width half maximum (\(FWHM\)) of these curves. However they are related to them by: \(FWHM_{lorentz} = 2\times LINEW\) and \(FWHM_{gauss} = 2.355 \times INLINEW\). For the VOIGT curve, it is a little more complicated but in terms of the other FWHMs, it can be aproximated as \(FWHM_{voigt} = 0.5346\times FHWM_{lorentz} + \sqrt{(0.2166\times FWHM_{lorentz}^2 + FWHM_{gauss}^2) }\).
The resolution of the spectrum can be modified with the SPECRES flag. By default it is a fraction of the LINEW. Please be aware that higher resolution (smaller SPECRES), mean a larger grid for the correlation function and more time points to calculate on.

7.45.1.10. General¶

The temperature is accounted for exactly on Absorption and Emission [199] and can be set using the TEMP flag.
PRINTLEVEL can be set to HIGH in order to print more details, including Huang-Rhys factors which are useful for rationalizing the contribution of different vibrational modes to the rate/spectrum.
The frequencies read from the Hessian files can be scaled by any number by setting the SCALING flag under %ESD. The default is 1.0.
If necessary, the transition dipole can also be scaled by setting the TDIPSCALING flag.
If you just want to compute an ES PES and stop, set WRITEHESS to TRUE and the correlation function will be skipped.
In order to make use of the fastest algorithm, set SAMEFREQ to TRUE and the DO method will be used, assuming equal Hessians between initial and final states and maximizing the efficiency when calculating the correlation function.
If you want to calculate phosphorescence rates, you MUST input the adiabatic energy difference DELE manually (the energy difference between each state at its own geometry). And, of course, don’t forget to set the SOC module to true.

7.45.2. Intersystem crossing rates¶

7.45.2.1. General Aspects of the Theory¶

Intersystem crossing (ISC) rates between a given initial state \(i\) and a final state \(f\) can be calculated from Fermi’s Golden rule:

(7.287)¶\[ k(\omega)_{if} = \frac{2 \pi}{\hbar} | \langle \Psi_f | \widehat{H}_{SO} | \Psi_i \rangle |^2 \delta(E_i - E_f), \]

which is quite similar to the Eq. (7.275) for Fluorescence, except for the frequency term. The same trick used there can be applied here to swtich to the time domain. Then, we are left with a simple time integration, which is not anymore difficult to solve than the equations above.

One can use all of the infrastructure already presented to compute these ISC rates, including Duschisnky rotation, vibronic coupling effects, use of different coordinate systems and so on. Right now, its use is optimized for CIS/TDDFT, as explained in Section ISC, TD-DFT and the HT effect, but it can be applied in general by combining simpler methods to obtain the geometries and Hessians with more advanced methods to compute the SOC matrix elements, when needed.

7.45.2.2. Tips and Tricks¶

The DELE must be given when using ESD(ISC), it is not automatically computed.That is the energy of the initial state minus the energy of the final state, each at its own geometry.
A SOC matrix element calculated from any method can be given on the input using the SOCME Re, Im flag, where these are the real and imaginary parts of that number.
The SOCMEs from TD-DFT are not bad, maybe except for those between the ground singlet and the triplets. In that case, a multireference calculation might be the preferred option.
If the final state is higher in energy than then initial state, the DELE is a negative number. In that case, there is barrier to go up when doing the ISC and the rates becomes more sensitive to the temperature.
In contrast to Fluorescence, the ISC rates depend strongly on the inclusion of Duschisnky rotations, please take care when using USEJ FALSE.
The default LINES is GAUSS, and the default INLINEW of 250 \(cm^{-1}\) might be too large. One should always investigate it by varying the width a bit. Other LINES can increase the error, please take care when changing it.

7.45.3. Resonant Raman Spectrum¶

7.45.3.1. General Aspects of the Theory¶

Raman intensities can be obtained in many different ways, depending on the experimental set up. As discussed extensively by D. A. Long [197, 533], the part of it that is “set up independent” is the Scattering Factor (or Raman activity), given by:

(7.288)¶\[ S = 45a^2 + 7\gamma \hspace{10 pt} \left( \frac{C^2 m^2}{amu V^2} \right), \]

where the \(a\) is related to the isotropic part of the “transition polarizability” between an initial state and a final state with a different vibrational quantum number \(\langle \Psi_f| \hat{\alpha} | \Psi_i \rangle = \alpha_{if}\):

\[a = \frac{1}{3}|(\alpha_{xx}+\alpha_{yy}+\alpha_{zz})|\]

and \(\gamma\) is also related to its off-diagonal elements:

\[\gamma = \frac{1}{2} [|(\alpha_{xx}-\alpha_{yy})|^2 + |(\alpha_{yy}-\alpha_{zz})|^2 + |(\alpha_{zz}-\alpha_{xx})|^2 +6(|\alpha_{xy}|^2 + |\alpha_{yz}|^2 + |\alpha_{xz}|^2)].\]

This transition polarizability itself can be computed using Kramers, Heisenberg and Dirac (KHD) formalism and it can be shown that each of its Cartesian components can be calculated as a sum-over-states:

(7.289)¶\[ \alpha_{\rho\sigma}^{if} = \frac{1}{\hbar} \sum_{n \neq i,f} \left( \frac{\langle \Psi_f | \mu_\rho | \Psi_n \rangle \langle \Psi_n | \mu_\sigma | \Psi_i \rangle }{\Delta E_{ni} - \omega_{laser} + i\gamma_{lt} } + \frac{\langle \Psi_f | \mu_\rho | \Psi_n \rangle \langle \Psi_n | \mu_\sigma | \Psi_i \rangle }{\Delta E_{ni} + \omega_{laser} + i\gamma_{lt} } \right) \]

In (7.289), the sum is over any number of electronic excited states \(n\), \(\Delta E_{in}\) is the energy difference between the initial state and the excited, \(\omega_{laser}\) is the laser energy and \(\gamma_{lt}\) is the lifetime broadening to avoid a zero on the denominator. If we work with a laser for which the frequency is close to the excited state energy difference, the first term is much larger than the second and can approximate alpha by

(7.290)¶\[ \alpha_{\rho\sigma}^{if} \simeq \frac{1}{\hbar} \sum_{n \neq i,f} \left( \frac{\langle \Psi_f | \mu_\rho | \Psi_n \rangle \langle \Psi_n | \mu_\sigma | \Psi_i \rangle }{\Delta E_{ni} - \omega_{laser} + i\gamma_{lt} } \right). \]

This equation can be solved using a path integral approach by switching to the integral form of \(1/x\):

(7.291)¶\[ \frac{1}{x} = \frac{i}{\hbar} \int_0^{\infty}e^{-ixt/\hbar}dt \]

So that the components of \(\alpha_{if}\) can be given by:

(7.292)¶\[ \alpha_{\rho\sigma}^{if} \simeq \sum_{n \neq i,f} \frac{i}{\hbar^2}\int_{0}^{\infty}\langle \Psi_f | \mu_\rho | \Psi_n \rangle \langle \Psi_n | \mu_\sigma | \Psi_i \rangle e^{-it(\Delta E_{in} - \omega_{laser} - i\gamma_{lt}) }dt \]

From here on, it is possible to show that the \(\alpha_{\rho\sigma}^{if}\) can be calculated as an integral of a correlation function [197], which is similar to the one previously discussed. In order to calculate that, a path integral scheme is also used and a geometry and Hessian for the ES are needed. The ORCA_ESD module predicts the ES PES (if not inputed), computes the \(\alpha_{if}\) and then prints the Scattering factor on a spectrum named BASENAME.spectrum.LASERE.

OBS.: The actual Raman Intensity collected with any polarization at 90 degrees, the I(\(\pi/2; \parallel^{s}+\perp^{s},\perp^{i}\) [533]), can be obtained by setting RRINTES to TRUE under %ESD.

OBS2.: In the current implementation, if a multistate calculation is requested, Eq. (7.290) is solved for each state and all spectra are summed up afterwards.

7.45.3.2. Specific Keywords and Details¶

In order to solve Eq. (7.292), the same information as for Absorption/Emission is needed and to compute the ES PES all of the above approximations are also valid. The main difference here is that a laser energy, given by the LASERE flag, should be given. If it is not given, the default is to set it to the 0-0 energy difference between the ground and the excited state. As mentioned before, more information can be found on Sec. Resonant Raman Spectrum.

You can choose more than one laser energy by setting LASERE 1,2,3,4 and etc. If so, each spectrum will be saved on a different BASENAME.spectrum.LASERE file.
You can also choose more than one excited state to be accounted for with the flag STATES 1,2,3, etc. and the final spectrum will be the sum of all of Scattering Factors given by the \(\alpha_{if}\)s. You can NOT choose several states and laser energies currently.
The automatic selection for the integral grid is also done based on the same ideas as mentioned before.
The default lineshape for resonant Raman is VOIGT.
The lineshape of the RR spectra will be taken from the RRSLINEW flag. In this case, LINEW and INLINEW are used only in the calculation of the correlation function.
Currently the only temperature for which this model works is at 0K.
In terms of which vibrationally excited states to be considered, currently it goes up to Raman Order 2, which means fundamentals, first overtones and combination bands (up to a total quantum number of 2). You can reduced that using RORDER flag.
It is also possible to include HT effect here for weak transitions, but be aware the calculation is much more costly. Due to technical reasons, the data is saved only on memory so, if you plan to go being 300 modes and do HT, there should be A LOT of memory available, about \(8 \times NMODES^4\). Also, you should expect a VERY long time for the computation of the correlation function. We are currently working on ways to accelerate this particular case.
As it is explained on the reference paper [197], the calculations using both Duschinsky rotation and HT effect can be greatly accelerated setting cutoffs for the derivatives and J matrix elements. The RRTCUTDER is a ratio with respect to the transition dipole moment below which the derivatives will be ignored and RRTCUTJ is a cutoff for J matrix elements. As the paper suggests, RRTCUTDER = 0.001 and RRTCUTJ = 0.01 are in general good numbers. We do recommend using these, but please be aware of the specific needs of your system.

7.45.4. Circular Polarized Spectroscopies¶

7.45.4.1. General Aspects of the Theory¶

Starting from ORCA 6, the ESD module has been expanded to include calculations of CP rates and spectra for chemical system. This enables the computation of the vibronic effect on ECD, CPF, and CPP spectra. The methodology is generally similar to the interaction of unpolarized light with matter, including absorption and photoluminescence. However, when a CP photon interacts with an optically active system, the electric field of the photon induces a linear displacement of the charge (transition electric dipole), while the magnetic field induces a circulation of the charge (transition magnetic dipole). These combined interactions cause an electron to be excited in a helical motion, involving both translation and rotation, along with their associated operators.

In a commonly use practice by defining a laboratory frame in which the z-axis defines the direction of the light trajectory, CP light interactions can be generated with the use of the complex vectors \(\mathcal{E} \pm= \frac{1}{\sqrt{2} } (\hat{x} \pm i \hat{y})\)

In this framework the FFMIO operator transforms as:

(7.293)¶\[ {\mathrm{T}}_{\mathrm{IF}}^\mathrm{\pm} =1/2\sum_{j=1}^{N} \left\langle\mathrm{I}\middle|\mathrm{e}^{{\mathrm{-ikr}}_\mathrm{j}} \left(\mathrm{\varepsilon\bullet}{\hat{\mathrm{p}}}_\mathrm{x}\right)\middle| \mathrm{F}\right\rangle\mathrm{\pm}\left\langle\mathrm{I} \middle|\mathrm{e}^{{\mathrm{-ikr}}_\mathrm{j}}\left(\mathrm{\varepsilon\bullet}{\hat{\mathrm{p}}}_\mathrm{y}\right)\middle|\mathrm{F}\right\rangle \]

In both ECD and CPL (CPF or CPP) spectroscopies the measured intensities are related to the difference of absorption or luminescence of the left and right polarized transition moments given by:

(7.294)¶\[\Delta_{IF}^{L\pm R}(k,\epsilon)={{|T}_{IF}^-|}^2\pm{{|T}_{IF}^+|}^2 \]

which leads to the following expressions for the sum and the difference of the square moduli \(|{T_{IF}^\pm|}^2\):

(7.295)¶\[\Delta_{IF}^{L+R}(k,\epsilon)=1/2\left\langle I\middle|\sum_{j=1}^{N}{e^{{-ikr}_j}\left(\varepsilon\bullet{\hat{p}}_x\right)}\middle| F\right\rangle\left\langle I\middle|\sum_{j=1}^{N}{e^{{-ikr}_j}\left(\varepsilon\bullet{\hat{p}}_y\right)}\middle| F\right\rangle \]

(7.296)¶\[\Delta_{IF}^{L-R}(k,\epsilon)=-\mathbf{Im}\left(\left\langle I\middle|\sum_{j=1}^{N}{e^{{-ikr}_j}\left(\varepsilon\bullet{\hat{p}}_x\right)}\middle| F\right\rangle\left\langle I\middle|\sum_{j=1}^{N}{e^{{-ikr}_j}\left(\varepsilon\bullet{\hat{p}}_y\right)}\middle| F\right\rangle\right) \]

Hence within the ED approximation, ECD and CPL radiative transition rates can be calculated through the orientational average of Equation (4), employing the Fermi’s Golden rule:

(7.297)¶\[ k_{ECD}\left(\omega\right)=\frac{{16\pi}^2\omega_{ECD}}{3}\ \sum_{F}\mathbf{Im}\left(\left|\left\langle\Psi_I\left|\hat{\mu}\right|\left.\Psi_F\right\rangle\right.\left\langle\Psi_F\left|\hat{m}\right|\left.\Psi_I\right\rangle\right.\right|\right)\delta\left({E_{FI}\pm\omega}_{ECD}\right) \]

(7.298)¶\[ k_{CPL}\left(\omega\right)=\frac{{16\omega_{CPL}}^3n^2}{3\hbar c^3}\ \sum_{F}\mathbf{Im}\left(\left|\left\langle\Psi_I\left|\hat{\mu}\right|\left.\Psi_F\right\rangle\right.\left\langle\Psi_F\left|\hat{m}\right|\left.\Psi_I\right\rangle\right.\right|\right)\delta\left({E_{FI}\pm\omega}_{CPL}\right) \]

where \(\omega_{ECD/CPL}\) are the excitation and emission CP photon energies, respectively while \(\omega_{FI}\) are the energies between the initial and final states reached in the absorption or the photoluminescent processes. Similarly \(E_{FI}\) is the transition energy and \(\delta\) refers to the line-broadening mechanism arising from the lifetimes of the relevant final states and c is the speed of light. In the above expressions \(\hat{\mu}\) defines electric dipole operator while \(\hat{m}\) is the respective magnetic dipole operator \(\hat{m}=\frac{1}{2m_ec}\sum_{i}{r_i\times\widehat{p_i}}\) and \(m_e\) is the electron mass. In the above expressions \(\mathbf{Im}\left(\left|\left\langle\Psi_I|\hat{\mu}|\left.\Psi_F\right\rangle\right.\left\langle\Psi_F|\hat{m}|\left.\Psi_I\right\rangle\right.\right|\right)\) represents the rotatory strength (\(R_{IF}\)). As discussed above the transition rates including the vibronic coupling, on the Frank-Condon (FC) and Herzenberg-Teller (HT) limits, can be efficiently proceed through the path integral approach[199] in which it is possible to calculate \(k_{ECD/CPL}^{obs}\) from the Fourier Transform (FT) of the respective correlation function \(\chi_{(t)}\) that is computed from the path integral of the multidimensional harmonic oscillator according to[784]:

(7.299)¶\[ k_{ECD/CPL}^{obs}(\omega)=2\alpha\mathcal{R}e\int_{0}^{\infty}\chi_{\left(t\right)}e^{\pm i\omega t}dt\]

with \(\alpha\) being a collection of constants and for CP transition one-photon rates (ECD, CPL) considering electric dipole and magnetic dipole interactions in the expression of the rotatory strengths it takes the form:[784]

(7.300)¶\[ \chi_{\left(t\right)}=e^{\pm i\omega t}[ \mathbf{Im}\left[\mu_em_e^\ast\right]\rho^{FC}\left(t\right) +\sum_{k}\mathbf{Im}\left[\mu_e\frac{\partial m_e^\ast}{\partial{\bar{Q}}_k}\right]\rho_k^{\frac{HT}{FC}}\left(t\right) +\sum_{k}\mathbf{Im}\left[\frac{\partial\mu_e}{\partial{\bar{Q}}_k}m_e^\ast\right]\rho_k^{\frac{HT}{FC}}\left(t\right) +\sum_{kl}\mathbf{Im}\left[\frac{\partial\mu_e}{\partial{\bar{Q}}_k}\frac{\partial m_e^\ast}{\partial{\bar{Q}}_l}\right]\rho_k^{HT}\left(t\right)] \]

where \(\mu_e\) and \(m_e\) represent the respective transition dipole \({\langle\Psi}_I|\hat{\mu}|\left.\Psi_F\right\rangle\) and magnetic dipole \(\left\langle\Psi_F|\hat{m}|\left.\Psi_I\right\rangle\right.\) moment integrals between initial and final states \(I, F\) while:

(7.301)¶\[ \rho^{FC}=Tr(e^{-i\hat{H}\tau}e^{-i\hat{H}\tau}) \]

(7.302)¶\[ \rho_k^{HT/FC}=Tr({{\bar{Q}}_ke}^{-i\hat{H}\tau}e^{-i\hat{H}\tau}) \]

(7.303)¶\[ \rho_{kl}^{HT}=Tr({{\bar{Q}}_ke}^{-i\hat{H}\tau}{\bar{Q}}_le^{-i\hat{H}\tau}) \]

where, these traces can be evaluated following the approach discussed in Ref[199].

Finally, it is quite commonly that the ECD and CPL spectral intensities are represented against normalized absorption and photoluminescent intensities defining, similar expressions for, the dissymmetry factors \(g_{abs}\) and \(g_{lum}\):

(7.304)¶\[ g_{abs/lum}\ =2\frac{I_{LCP}\ -\ I_{RCP}}{I_{LCP}\ +\ I_{RCP}}\ \ ~\ \left.\frac{4R}{D}\right|_{GS/ES\ Structure};\ \ \ -2<g_{abs/lum}<2 \]

Implying that the associated dissymmetry spectra can also be calculated, where D and R the square of the transition dipole and the rotatory strength, respectively.

7.45.5. Magnetic Circular Dichroism¶

7.45.5.1. General Aspects of the Theory¶

The formulation presented for the calculation of the absorption spectrum may be extended to the absorption of circularly polarized light (CPL) of a system under the effect of an external magnetic field in order to compute the MCD spectrum.[268] By assuming an electric dipolar approximation under a length formulation for the light-matter interaction, the molar absorptivity contribution by the transition between an initial state i and a final state f of one fixed oriented molecule may be expressed as:

(7.305)¶\[ \epsilon(\omega)_{if} = \alpha \omega |\mathcal{E} \cdot \langle \Psi_i| \hat{\mu} | \Psi_f \rangle|^2 \delta(E_i - E_f \pm \hslash \omega) \]

where \(\mathcal{E}\) is the polarization vector of the incident light, and \(\alpha\) is a positive collection of constants.

By considering the case in which the light is propagating in the laboratory fixed \(\hat{z}''\) direction, the circularly polarized light is described by \(\mathcal{E} = \frac{1}{\sqrt{2} } (\hat{x}'' \pm i \hat{y}'')\) where the “\(-\)” sign corresponds to the left circularly polarized light and the “\(+\)” to the right circularly polarized light.

Similarly, as the MCD calculations presented in section Simulation of (Magnetic) Circular Dichroism and Absorption Spectra, the total absorptivity may be obtained by summing over all possible transitions and averaging the molecular orientations by using 3 rotation angles \(\theta\), \(\phi\), and \(\chi\).

(7.306)¶\[ \epsilon(\omega) = \alpha \omega \int_0^{2\pi}\int_0^{2\pi}\int_0^{\pi} \sum_{if} \frac{e^{- \frac{\epsilon_i}{k_B T} }}{Z} |\mathcal{E} \cdot \langle \Psi_i| \hat{\mu} | \Psi_f \rangle|^2 sin(\theta) \delta(E_i - E_f \pm \hslash \omega) d\theta d\phi d\chi \]

Considering \(\chi\) the angle which rotates the molecule in the plane perpendicular to the magnetic field perturbation (which is colinear with the incident light propagation direction), the integrals of eq. (7.306) may be computed conveniently in a seminumerical scheme by using an intermediate \(r'\) frame.

(7.307)¶\[ \epsilon(\omega) = \alpha \omega \sum_p w_p \sum_{if} \frac{e^{- \frac{\epsilon_i}{k_B T} }}{Z} \sum_{\beta \gamma = x',y'} \overline{\overline{\mathcal{E} }}_{\beta \gamma} \langle \Psi_i| \hat{\mu}_\beta | \Psi_f \rangle \langle \Psi_f| \hat{\mu}_\gamma | \Psi_i \rangle \delta(E_i - E_f \pm \hslash \omega) \]

where \(\theta\) and \(\phi\) values are defined by a point p in a numerical grid, \(x'\), and \(y'\) are determined by \(\theta\) and \(\phi\) values according to eq. (7.308), and \(\chi\) has been integrated analytically in the \(\overline{\overline{\mathcal{E} }}_{\beta \gamma}\) values.

(7.308)¶\[\begin{split} \begin{split} \hat{\mu}_{x'} &= cos(\theta) cos(\phi) \mu_{x''} + cos(\theta) sin(\phi) \mu_{y''} - sin(\theta) \mu_{z''} \\ \hat{\mu}_{y'} &= -sin(\phi) \mu_{x''} + cos(\phi) \mu_{y''} \end{split} \end{split}\]

(7.309)¶\[ \overline{\overline{\mathcal{E} }}_{\beta \gamma}^\pm = \int_0^{2\pi} \mathcal{E}_\beta^\pm \mathcal{E}_\gamma^{\pm*} d\chi \]

Finally, by applying the Bohr-Oppenheimer approximation and writing the Dirac delta function in the time domain, we get:

(7.310)¶\[ \epsilon(\omega) = \alpha \omega \sum_p w_p \sum_{if} \frac{e^{- \frac{\epsilon_i}{k_B T} }}{Z} \sum_{\beta \gamma = x',y'} \overline{\overline{\mathcal{E} }}_{\beta \gamma} \langle \Theta_i| \hat{\mu}_\beta^e | \Theta_f \rangle \langle \Theta_f| \hat{\mu}_\gamma^e | \Theta_i \rangle \int e^{i (E_i-E_f-\omega)t} dt \]

and by taking the difference of absorbance between left and right CPL produce we reach an expression of the MCD intensities:

(7.311)¶\[ \Delta \epsilon(\omega) = \epsilon^-(\omega) - \epsilon^+(\omega) = -2 \frac{\alpha}{Z} \omega \sum_{if} \sum_p w_p Im \Big[\int_{-\infty}^\infty \tilde{\chi} (t,x',y') e^{-i \hslash \omega t} dt \Big] \]

where \(\tilde{\chi}\) under a first-order approximation of the transitions moments with respect to the nuclear displacement is:

(7.312)¶\[ \tilde{\chi} (t,\beta,\gamma)= e^{i \Delta E t} \Big[\hat{\mu}_{0\beta}\hat{\mu}_{0\gamma}^* \rho^{FC}(t) + \sum_k \frac{\partial \hat{\mu}_{\beta} }{\partial Q_k}\hat{\mu}_{0\gamma}^* \rho^{FC/HT}(t) + \sum_k \hat{\mu}_{0\beta} \frac{\partial \hat{\mu}_{\gamma}^*}{\partial Q_k} \rho^{FC/HT}(t) + \sum_{kl} \frac{\partial\hat{\mu}_{\beta} }{\partial Q_k} \frac{\partial\hat{\mu}_{\gamma}^*}{\partial Q_l} \rho^{HT}(t) \Big] \]

Similarly, as section Simulation of (Magnetic) Circular Dichroism and Absorption Spectra, the transition moments under the effect of an external magnetic field perturbation may be estimated by using a QDPT (eq. (7.274)), and the derivatives are approximated numerically in a similar way as the ESD-Absorption case.

7.45.5.2. Specific Keywords and recommendations¶

Once selected the ESD(MCD) calculation two variables need to be defined in the %esd block:

The intensity of the magnetic field in Gauss “B”
The grid to make the molecular orientational average “LEBEDEVINTEGRATIONPOINTS”

The MCD signals use to be much more sensitive than the corresponding absorption intensities. As result, the MCD calculations are much more sensitive to errors in the electronic structure. So it is highly recommended to first fully understand the electronic structure problem and verify there are no important problems. In this direction we suggest the following recommendations:

Be sure the obtained electronic states are in the correct order by assigning point group symmetry labels and comparing them with better electronic structure methods. Due to the MCD intensities emerging as a result of the interaction of states by the magnetic field perturbation, a wrong-located state in energy may affect the MCD intensities, even if it is not in the energy range you are computing. We do not recommend using the method as a black box.
It is recommended to do first an ESD-absorption calculation on the exact same level of theory, verify the intensities and also solve any problem related to the PES representation.
The lineshape for the best agreement of the MCD intensities compared against the experimental measurements may differ for the best value for absorption intensities.

7.45.6. Complete Keyword List for the ESD Module¶

%ESD                             #The booleans are the defaults
  ESDFLAG    ABS   (default)     #Which calculation to make?
             ECD
             FLUOR (default)
             CPF/CPFLUOR
             PHOSP
             CPP/CPPHOSP
             MCD                 #Only available with TDDFT
             RR 
             ISC
             IC

  GSHESSIAN  "BASENAME.hess"     #The ground state Hessian
  ESHESSIAN  "BASENAME_S1.hess"  #The excited state Hessian
  TSHESSIAN  "BASENAME_T.hess"   #The triplet state Hessian

  HESSFLAG   AHAS                #How to obtained the ES PES?
             VH
             VG (default)
             VGFC                # VG + APPROXADEN TRUE
             HHBS
             HHAS
             UFBS
             UFAS

  STATES     1,2,3,4             #IROOTS to compute

  MODELIST   4,5,6               # only include the 4th, 5th and 6th vibrational modes
                                 # in the ESD calculation
  SINGLEMODE 1,5,10              # run three ESD calculations, each considering only
                                 # one of the 1st, 5th and 10th vibrational modes

  DOHT       FALSE               #Do HT effect?
  FASTDER    TRUE                #Use the fast derivatives algorithm?
  DELQ       0.01 (default)      #Dimensionless displacemente for derivatives
  STDA       FALSE               #Use sTDA during derivatives?

  APPROXADEN FALSE               #Compute the DELE by extrapolating info from
                                 #the Hessians, avoiding second single-point 
                                 #at the ES geometry                

  USEJ       FALSE               #Consider Duschinsky rotations?
  USEB       TRUE                #Rotate the initial state?
  SAMEFREQ   TRUE                #Use DO method and J=1.

  DELE       12000               #Custom energy difference
  TDIP       x,y,z               #Custom transition electric dipole
             x.re,x.im,y.re,y.im,z.re,z.im
  TMDIP      x,y,z               #Custom transition magnetic dipole
             x.re,x.im,y.re,y.im,z.re,z.im
  SOCME      x.re,x.im           #Custom spin-orbit coupling matrix elements, in a.u.
  LASERE     34000               #The laser energy for RR

  B          3000.0              #Magnetic field in Gauss for MCD
  LEBEDEVINTEGRATIONPOINTS 14    #Lebedev Grid for MCD molecular
                                 #orientational average

  GEOMSTEP   AUGHESS (default)   #Geometry step?
             QN
  STEPSCALING  1.0               #A number for scaling the steps
  STEPCONSTR   0,2,5             #A list of atoms that will not be moved
  COORDSYS   CART                #Coordinate system for the normal modes?
             INT (default)
             WINT
             FCWL
             FCWS

  TCUTFREQ   100                 #Cutoff for frequencies
  IFREQFLAG  POSITIVE (default)  #What to do with negative frequencies?
             LEAVE
             REMOVE

  UF_DELE    1E-4                #Energy difference for updated freq.
  UFFREQERR  0.2                 #Tolerated percentual error

  TEMP       298.15 (default)    #Temperature to consider
  UNIT       WN                  #wavenumbers  (Output units - input still in cm-1!)
             NM                  #nanometers
             EV                  #electron volts
  CENTRALDIFF  TRUE              #Central differences?
  CONVDER    FALSE               #Convert derivatives between state
  SCALING    1.0 (default)       #Scaling for frequencies
  TDIPSCALING  1.0 (default)     #Scaling for the transition dipole

  LINES      DELTA               #The lineshape function
             LORENTZ (default)
             GAUSS
             VOIGT (default for RR)
  LINEW      50 (default)
  INLINEW    250 (default)
  RRSLINEW   10 (default)
  
  RORDER     2 (default)         #The Raman Order for RR
  RRINTENS   false               #Calculate the intensities instead
  RRTCUTDER  0 (default)         #A cutoff for derivatives
  RRTCUTJ    0 (default)         #A cutoff for J matrix elements

  WRITEHESS  FALSE               #Make ES PES and leave

  MAXTIME    12000               #Max time (atomic units!) for the FT
  NPOINTS    131072              #Total number of points

  SPECRANGE  10000,50000         #Spectral range
  SPECRES    1.0                 #Spectral resolution
    
  PRINTVIB   FALSE               #Save a .xyz file per each vibational mode.
                                 #(Requires DOHT true)
  PRINTLEVEL 0                   #If set to 1 (or 2, 3, ...), prints additional
                                 #details of the calculation
END