Selected Aspects Concerning the Efficient Calculation of Vibrational Spectra beyond the Harmonic Approximation †

This feature article discusses some selected aspects in the field of vibrational structure calculations based on vibrational self-consistent field, VSCF, and vibrational configuration interaction, VCI, theory. As the quality of such calculations depends strongly on the accuracy of the underlying multidimensional potential energy surface, PES, some techniques will be discussed to establish high-quality PESs in a fully automated manner. As an alternative to VCI theory multiconfiguration self-consistent field, VMCSCF, theory and in particular specific aspects concerning the integral evaluation relevant to both approaches will also be presented. Further aspects concern the efficient calculation of infrared intensities and Franck-Condon factors in vibronic transitions.(doi: 10.5562/cca2149)


INTRODUCTION
Vibrational spectroscopy is one of the very basic spectroscopic methods, which are routinely applied in chemistry and biophysics for the characterization and identification of molecules.Moreover, since it is a fast technique, it can efficiently be used to identify short-living reactive intermediates and can thus be used for the investigation of reaction paths.Students at the undergraduate level learn the techniques of infrared and Raman spectroscopy quite early in order to understand fundamental aspects of the dynamics and structure of molecules.Most vibrational spectra are measured in solution at room temperature with a resolution between 1 and 10 cm −1 .Although these spectra typically suffice to identify the species of interest, they are not accurate enough to allow for an assignment of the most important vibrations of the molecule.In contrast to that, the accurate measurement of vibrational spectra is not a routine application and depends strongly on the physical properties of the probe.For example, in order to measure a gas phase spectrum at low pressures and temperatures the species of interest needs to have a certain vapour pressure, which is rarely fulfilled by larger systems.Matrix isolation techniques are a wonderful tool to trap reactive intermediates and thus to record the spectra of instable systems, 1 but these spectra typically suffer from the interaction between the probe and the matrix, which increases with the polarizability of the rare gas matrix atoms.Consequently, the determination of accurate band positions requires advanced techniques typically applied by some specialists only.Besides this, the correct assignment of the observed peaks to the individual vibrations is a non-trivial task, which has led to a large number of errors in the literature.For that reason, more and more experimentalists augment their measurements by ab initio calculations, which provide the assignment for free.However, as most ab initio methods rely on approximations and thus do not yield exact results, the original problem often is only transferred to the problem of correlating the experimental data with computed ones.Nevertheless, vibrational frequency calculations are usually a very useful tool to guide the assignment of the individual bands.
Most vibrational frequency calculations make use of density functional theory (DFT) within the harmonic approximation.Besides 2 nd order Møller-Plesset perturbation theory (MP2) this essentially is the only approach, which can be applied to systems with more than 100 atoms.However, the accuracy of such calculations is limited, even if subsequent scaling of the frequencies or force constants is applied. 2 Maximum deviations may easily exceed 50 cm −1 , which is of little help once several vibrational transitions can be seen in a small spec-Croat.Chem.Acta 85 (2012) 379.troscopic window.The scaled quantum mechanical force field approach (SQM) of Pulay et al. 2−4 diminishes this error significantly, but relies on a set of empirical scaling factors, which are optimized on the basis of experimental data.For inorganic systems or molecules with unusual bonding patterns as for example in reactive intermediates, usually no reference data are available und thus the applicability of this method remains limited.However, for organic systems with more than 50 atoms it is most likely the method of choice.
In contrast to these empirical approaches, all ab initio methods face the problem that anharmonic contributions to the potential need to be computed explicitly.This constitutes a severe bottleneck and limits these methods to small and medium sized systems only.Independent on the nature of the coordinates used, i.e. rectilinear vs. curvilinear, due to the curse of dimensionality a complete representation of the multidimensional Born-Oppenheimer surface with respect to the number of coordinates is only possible for systems with 3 or 4 atoms.Consequently, any kind of truncated expansion has to be used to allow for the calculation of larger systems.This inevitably leads to deviations from the exact representation and thus the choice of the coordinates.Accordingly, the order of the expansion determines the accuracy of the results.Most programs available use either a Taylor-expansion or a multi-mode expansion of the PES, the latter essentially being a Taylor-expansion with reordered terms.For example, terms of the type 2 ij i j f q q correspond to 2-mode coupling terms in a multimode expansion, but to 3 rd order terms ( iij i i j f q q q ) in a Taylor-expansion. 5,6ypically the multi-mode expansion converges slightly faster and can thus be truncated after the 3 rd order for most systems, while standard Taylor-expansions require 4 th order terms for a proper representation of the potential.These truncations allow for the calculation of multi-dimensional PESs for systems up to 20 atoms at high levels of electronic structure theory.The choice of an appropriate coordinate system to be used within the different expansions of the PES is not as obvious as one may wish.Rectilinear normal coordinates show several advantages: Firstly, the kinetic energy operator in the subsequent vibrational structure calculations becomes fairly simple and easy to implement.Secondly, and more importantly, normal coordinates are unique and well-defined.On the contrary, normal coordinates fail in the description of internal rotations in floppy molecules or molecular clusters.This problem can be solved by the use of curvilinear internal coordinates. 7−9 However, internal coordinates are not unique and due to the fact that most expansions of potentials are not fully converged, this will lead to coordinate dependent results.Consequently, alternative schemes, as for example the use of redundant interatomic distances as proposed by Braams and Bowman, 10,11 are very promising schemes to avoid this problem.Moreover, using internal coordinates the kinetic energy operator becomes very complicated once it shall be represented by analytical expressions. 12,13For that very reason its numerical evaluation became rather popular in the recent past. 14,15nce the potential energy surface has been represented in a proper way, the routes to determine vibrational transitions are manifold.−22 However, the problems of this approach are on hand: with increasing size of the system, the number of (quasi)resonances increases and thus the perturbational approach fails.The identification of resonances and its corresponding treatments 23−26 is a keystep in these approaches but vary from program to program and thus the final results differ significantly for sensitive systems.Moreover, VPT2 theory cannot deal with double-well potentials or potentials with important high-order corrections.Nevertheless, the cost-gain ratio of this method is excellent and its use is absolutely of black-box nature.Alternatively, one may consider variational approaches based on different variants of vibrational selfconsistent field theory, 27−30 VSCF, which are unlimited in accuracy and which are subject of this feature article.In principle these methods are applicable to any kind of PES, but of course require an explicit treatment of vibration correlation effects, which may constitute a severe computational bottleneck for large systems.Strong analogies can be established to electron correlation methods, but a closer look usually reveals that well-established electron correlation methods cannot be transferred to the vibrational problem without trouble.For example, Christiansen showed, that the Møller-Plesset series of methods, VMPn, typically does not converge within the calculation of CH stretchings and thus the lowest order of the series, i.e.VMP2, is already of questionable accuracy for these modes. 31Most of these problems arise from the high density of vibrational states, which is an almost unknown problem in electronic structure theory.Moreover, the VSCF solutions are often poor approximations to the true wave function, which leads to strong correlation effects and very small leading coefficients in the amplitudes of the subsequent correlation treatment.Besides this, in contrast to electronic structure calculations the inclusion of quadruple excitations in the correlation calculations is mandatory in order to obtain reliable results. 32,33This of course is an effect of the high order operators in the Hamiltonian.
The focus of this contribution is on some selected issues concerning the automated generation of potential energy surfaces expanded in terms of normal coordinates and the calculation of vibrational wave functions Croat.Chem.Acta 85 (2012) 379.5][36] Moreover, we will discuss some aspects concerning the efficient calculation of integrals relevant for all vibration correlation methods.Some minor aspects concerning infrared intensities and the calculation of vibronic transitions will also be discussed.We will not consider vibrational perturbation theory and specific approaches for rovibrational calculations of 3or 4-atomic molecules. 37For recent reviews in the field of vibrational structure theory see for example Refs.38 and 39.

AUTOMATED GENERATION OF POTENTIAL ENERGY SURFACES
An accurate multi-dimensional potential energy surface is a mandatory prerequisite for a successful vibrational structure calculation.In the past, such surfaces have been generated in a tedious procedure by hand, i.e. fitting parameters where determined step by step until the surface obeyed all criteria applied.In some cases the ab initio surface was subsequently corrected by experimental data, in order to yield best agreement with spectroscopic measurements. 40Of course this procedure is error-prone, time-consuming and boring.For that reason, an automated generation of highly accurate surfaces by any electronic structure code is desirable.However, this task faces a number of problems, which need to be considered explicitly in the set up of these generators.a. Iterative interpolation: It is well known that certain regions of a PES need to be described more accurately than others.For example, regions around minima need to be determined more precisely than regions of very high energies.As a consequence certain parts need to be represented by more ab initio grid points than others.Therefore, a PES generator, which uses a fixed coarse grained mesh of ab initio points with a subsequent interpolation to a fine grid, will inevitably lead to an inconsistent description of the PES unless the number of ab initio points is vast.Consequently, iterative procedures as the adaptive density-guided approach (ADGA) by Christiansen et al., 39 the adaptive generation of adiabatic PES (AGAPES) of Richter et al. 41 or our iterative interpolation scheme are necessary in order to obtain well-balanced PESs. 42It is the interplay of the surface generator and the fast evaluation of preliminary vibrational wave functions within the iterations, which allow to control important parameters of the surface, e.g. the extension of the surface can be controlled by the density of the vibrational wave function.In this context the exploitation of symmetry relationships is important in both, the potential energy surface and the individual electronic structure calculations, in order to minimize the number of ab initio calculations and to reduce numerical noise. 10,43Equilibration of surfaces: Let us consider a multimode expansion of a potential in grid representation, i.e.   ...
Each of the Δ-surfaces, V i , V ij , ... , is represented by a mesh of grid points.The contributions of higher order are corrections to those of lower order.The Δ-surfaces in grid representation can be transformed to an analytical representation in terms of polynomials and vice versa according to: ...
Within these equations q i denotes the ith coordinate and p the corresponding polynomial coefficients.Typically the vertices of the Δ-surfaces are highest or lowest in energy.Once one of the vertices is significantly higher than the others, i.e. a very unsymmetric potential, it may cause severe trouble, which may lead to a complete crash of the calculation in the worst case or at least to inaccurate results.The worst case typically arises once the underlying Hartree-Fock calculation, HF, does not converge.In such cases, a complete active space self-consistent field calculation, CASSCF, typically still does converge and can be used instead.
Of course the CASSCF and HF energies do not match and thus the HF values cannot simply be replaced by the CASSCF results.Several situations are realistic: (a) if the relative energy is very high, the exact value does not matter too much because the vibrational state energies usually are much lower in energy.(b) Even if the energy is very high and thus irrelevant, the inconsistent value may lead to a problematic fitting and subsequent interpolation.Therefore, the CASSCF wave function may serve as an extremely good initial guess for a new HF calculation, which often indeed does converge.Consequently, this trivial error correction scheme can be used to stabilize the PES generation and to extend the convergence radius of the approach.In order to save computation time this feature will only be active for those electronic structure calculations, in which the HF calculations fail.We have implemented such correction schemes in our PES generator.(c) An equilibration of the Δ-surfaces can be used, in which the reference structure is moved out of the center of the surface.As a result, the troublesome point of high energy can be shifted away and it does not matter if the energy is correct or not.This situation is shown in Figure 1.The shifting of the surface is intrinsically absorbed in the powerful ADGA approach of Christiansen, but of course without the bypass via the CASSCF calculations in problematic situations.In case that a transformation from the grid representation to multivariate polynomials is needed, this transformation has to be corrected for the shifting of the center accordingly.For the case of a shifted 2D Δ-surface, the polynomial representation is given as where   ij rs c denotes the polynomial coefficients and δ i and δ j are the shifts with respect to the coordinates q i and q j .As this form is inconvenient to be used in the integration routines of the vibrational structure calculations, it needs to be transformed to the form shown in equation 3. Consequently the polynomial coefficients with and likewise for all higher terms of the expansion of the potential.As a result the PES generation is significantly stabilized and the energy range within the fitting routines is decreased.

c. Multi-level calculations:
The idea of determining the different orders of the Δ-surfaces, e.g.V i vs. V ij (cf.Eq. 1), by different electronic structure methods is iden-tical to the fundamental ideas of QM/MM calculations or the treatment of distant orbitals pairs in local correlation methods: 44,45 contributions of less importance can be described at a lower level of electronic structure theory.Within the framework of VPT2 theory this idea has been used for a long time.While the harmonic frequencies were computed at high level, e.g.CCSD(T)/cc-pVTZ, the anharmonic corrections were determined at low level, MP2/cc-pVTZ, and then transferred to the high level harmonic frequencies.This leads to significant computational savings and thus allows for the calculation of larger systems.−50 The simplest multi-level scheme one can think of is the use of low quality electronic structure calculations, e.g.MP2/cc-pVTZ, for the high-order terms of the PES expansion.However, for systems with more than 8 atoms, the number of 3D Δ-surfaces becomes very large and even if these terms are treated at very low levels, they will dominate the computational cost of the entire calculation.Actually, most of the time is spent in the unavoidable Hartree-Fock or DFT calculation.For that reason more sophisticated methods have been developed in order to accelerate this step.
In analogy to the modified Shepard interpolation, Matito et al. 51 used low-level Taylor expansions, i.e. gradient and hessian information at the coarse grained mesh of ab initio points, in order to model grid points of higher order.For example, 3D terms can be generated from the gradients and hessians of the 2D surfaces.This approach has the advantage of using the high-level information of the 1D and 2D mode-coupling terms of the potential and thus avoids the choice of a low-level method for the high-order terms, which may be inappropriate for the system of interest.However, additional computational costs arise for the gradient and hessian information.Alternatively we have developed a modeling scheme which uses the high-level information of the low-order terms and which is based on semiempirical MO theory. 52Based on a simple least-squares criterion the parameters within any semiempirical approach can be optimized with respect to the high-level low-order terms of the expansion of the potential.In order to accelerate the non-linear optimization, analytical gradients of the semiempirical energy with respect to the inherent parameters had to be derived.This approach may lead to an acceleration by several orders of magnitude and typically shifts the computational bottleneck back to the calculation of the 1D and 2D terms.However, due to the limitations of standard semiempirical MO theory, elements with d-shell occupations cannot be handled presently.However, an extension to the new generation of semiempirical methods as developed by Thiel et al. can resolve this dilemma. 53,54Within this scheme, semiempirical MO theory is used as a sophisticated set of fit functions, which will be optimized for the individual molecules.Typically this scheme results in mean absolute deviations of 1-2 cm −1 for the final results.
As mentioned above, this modeling scheme shifts the computational bottleneck to the 1D and 2D terms.However, modern electronic structure methods have lowered the demands for these calculations as well.6][57] As a rule of thumb, an explicitly correlated calculation based on a basis set of n-ζ quality corresponds to a standard calculation with a basis set of (n + 2)-ζ quality.Consequently, the basis set limit is reached much earlier.This leads to significant savings in the electron correlation calculation and the underlying Hartree-Fock calculation.Note, that the remaining basis set error in the underlying Hartree-Fock calculation can be diminished by a perturbative correction. 55s large basis sets are mandatory within the generation of accurate potential energy surfaces, explicitly correlated coupled-cluster calculations, e.g.CCSD(T)-F12a, or explicitly correlated multi-reference configuration interaction calculations, i.e.][60] However, for systems with more than 10 atoms even these calculations become prohibitively demanding.However, we have shown that linear scaling local electron correlation methods can be applied to the calculation of local potential energy surfaces without meaningful loss in accuracy. 61The combination of both approaches, i.e. local correlation concepts and explicitly correlated approaches, 62 can be applied to systems with more than 10 atoms without problems.In Table 1 we have shown first benchmark calculations for Li 2 F 2 making use of density fitting explicitly correlated local coupled cluster calculations, 62 DF-LCCSD(T)-F12a, for the 1D and 2D terms and our molecule-specific modeling scheme for the 3D terms. 52The comparison with respect to explicitly correlated coupled-cluster calculations for all three levels, i.e. 1D-3D, shows, that the deviations are very small and cancel each other out to some extent.The corresponding PES calculation for Li 3 F 3 took just 1 day on 4 Intel Xeon cores.This is a speed-up of about 60 with respect to the conventional calculation without the local approximation and the modeling of the 3D terms.These calculations can be further accelerated by exploiting the fact that the construction of PESs can be realized in an embarassingly parallel fashion. 63For example, the calculation for Li 7 F 7 using the same scheme as presented above for Li 2 F 2 took less than one day when using 10.000 cores on a Cray XE6.As a consequence, multi-dimensional potential energy surfaces of coupled-cluster quality can currently be computed for systems of up to 20 atoms.In order to optimize the load balance of such massively parallel calculations, we have interfaced the Molpro suite of ab initio programs 64 to the SEGL grid computing client. 65This combination allows even for the use of distributed architectures.As we believe that computers with several hundreds or even thousands of cores will be the standard in a few years, such calculations will be routine applications in the near future and thus algorithms for the subsequent vibrational structure calculations need to be adapted accordingly.

VIBRATIONAL STATE ENERGIES
0][71][72][73][74][75] Most of them have direct analogues in electronic structure theory, but usually differ in some important aspects.The fastest approach in this context is vibrational Møller-Plesset perturbation theory, VMP, which in contrast to its electronic counterpart is not restricted to double-excitations only. 31,76The use of VMP theory is rather problematic as due to the high state density quasi-degeneracies occur quite frequently in vibrational structure calculations.For this very reason quasi-degenerate perturbational approaches as developed by Yagi et al. 77 or Benoit et al. 78 are more useful.However, independent of this problem and due to the lifted restriction to double excitations, 2 nd order VMP theory, i.e.0][81][82][83] Several implementations have been presented, which make use of this scheme. 69,74,75,81Figure 2 shows the percentage of selected configurations with respect to the initial configuration space needed to obtain converged results for the vibrational ground states of (LiF) n clusters of increasing size.The initial configuration space in these calculations included up to quadruple excitations with a maximum excitation up to the 6 th modal per mode and a maximum overall excitation of 9. Figure 2 shows that only a very small fraction of configurations is important for a proper representation of the vibrational state energies and that the overhead of initially generated configurations becomes huge for larger systems.According to this, configuration selection schemes are very powerful in vibrational structure calculations.However, the development of restrictions to the correlation space prior to the selection is still needed as even the selection may become time consuming for large systems.Ground-state based and state-specific VCI theory and all its variants still are the most widespread and reliable methods for calculating accurate vibrational state energies.However, in the recent past impressive efforts have been undertaken in order to develop vibrational coupled-cluster methods, VCC, which in contrast to VCI theory show the advantage of size-consistency. 33,72,735][86] As several reviews are available on this topic, we will not reconsider all these approaches in detail here, but will focus on some specific aspects of vibrational structure calculations instead. 38,39 Integral evaluation: Let us consider a polynomial multi-mode representation of the potential, i.e.
Moreover, the vibrational wave function Ψ can be expressed in the form

 
with with I  being a configuration represented in the basis of (VSCF) modals i φ .I i n denotes the ith element of the occupation number vector n I belonging to configuration I. Using these definitions integrals of the type 3  3 with occur formally for the 3D terms of the potential function.This equation is relevant to all vibration correlation methods as it is part of the integration module.The integrals IJ ir Q can be kept easily in memory as they depend on just one coordinate q i .As a result one has to consider only the different occupation patterns of the modals for mode i, which in principle can arise for two different configurations rather than storing an array quadratically dependent on the number of configurations.Of course in practice equation 11 will never be computed this way.Different contraction schemes with increasing memory demands can be formulated, i.e.
These integrals can be computed a priori and thus CPU time is traded against memory demands.Note that the memory demands of the Y integrals is about 4 times as large as for the X integrals.This does not lead to a memory bottleneck for the 3D terms of the potential, but may so for the 4D contributions.As a result equation 11 simplifies to , , or , , The speed-up due to these contractions depends strongly on the implementation of the algorithm.Due to the shorter loop lengths for the orders of the polynomials (r, s and t) in contrast to the number of grid points in a grid based implementation, the savings are larger for grid-based versions than polynomial based algorithms.For example, considering a 4D potential for beryllium hydride, 43,87 Be 2 H 4 , the entire polynomial based VCI calculation including an iterative Jacobi-Davidson di-agonalization was accelerated by a factor of 2.8 for a single contraction according to equation 13.A further contraction as shown in equation 14 increased this factor to 3.2.Due to the larger memory demands of the Y integrals, a single contraction appears to be preferable over the twofold contraction.The factor of 2.8 increases to 5.1 for our grid-based implementation.For other test molecules we found even larger accelerations, e.g.CH 4 shows a factor of 8.3.The 2 nd aspect within the calculation of the integrals concerns the product over the Kronecker-δs.In a naive implementation one could calculate this product for each triple ijk and could decide then if the corresponding integral has to be evaluated or not.This of course would result in a very inefficient implementation.Essentially, analogues to the Slater-Condon rules 88 can be established, which determine the elements to be calculated.Considering the 3D terms of the potential, the rules depend on the difference Δn IJ of the two occupation number vectors n I and n J of the two configura-tions, i.e.IJ n denotes the number of different pairs of modals: :  , Δn IJ = 0 describes the contribution to a diagonal element of the VCI matrix and thus is identical to equation 15 without the product over Kronecker-δs as this product will always be 1.Note that in equations 18−20 the differences in the occupation number vectors fix some or all of the mode indices, i.e. i, j and k, respectively.We did not test the performance of the integral evaluation without equation 17 as this check is trivial.However, we tested two implementations with and without unrolling the loops.In dependence on the molecule, the size of the VCI space and other parameters we observed speed-ups of about 1.7 once equations 18−20 have been used for the calculation of the 3D terms and likewise relations for all other terms of the PES expansion.These simple considerations show that nonnegligible speedups can be achieved quite easily.
b. Multi-configuration self-consistent field theory: Most vibrational structure calculations rely on singlereference approaches, i.e. modals will be determined by the VSCF approach and correlation effects will be considered in a subsequent step (VCI, VCC, VMP2).However, the frequent occurrence of near-degeneracies due to the high state density very often leads to very small coefficients of the leading configuration in the correlation treatment.In case of strong (e.g.Fermi) resonances these coefficients can be smaller than 0.7 and thus the optimization of the modal basis to just one configuration leads to inconsistencies once the configuration space is incomplete.This situation corresponds to static electron correlation effects as any kind of resonance cannot be described by the VSCF approach and thus the VSCF wave function is qualitatively incorrect.For that reason Culot and Lievin 89,90 developed the vibrational multiconfiguration self-consistent field approach, VMCSCF, which can be considered the static variant of the multi-configuration time-dependent Hartree method, MCTDH, of Meyer et al. 91−93 Quite recently we presented a novel and efficient implementation of VMCSCF theory 34,35 based on Jacobi 2×2 rotations as originally proposed by Hinze for electronic wave functions. 94nce the Watson operator 95 is written in the form and the wave function is given as the VMCSCF energy can be written as ...
Minimization of the Lagrangian leads to an optimized set of modals r i φ  and VCI- coefficients I C  .The optimized modals can be ex- pressed in the basis of the initial ones via a unitary transformation and the antisymmetrical 2 which contain the rotational angles i rs α in the off- diagonal elements, U is simply given as a 2 × 2 Jacobi rotation matrix.For further details see Ref. 34.Consequently, these equations can be solved by consecutive Jacobi rotations for each mode (microiterations).Once a new set of modals has been obtained, a new VCI matrix will be generated (macroiterations) and the whole procedure will be repeated until self-consistency has been reached.Once a reasonable active space has been cho-sen, the final results typically are very close to those obtained from vibrational full configuration interaction calculations, VFCI.Table 2 shows the results of a VMCSCF calculation for chloroform in comparison to a VCI calculation with the same number of configurations based on a 3D potential of CCSD(T)-F12a/vtz-F12a quality.This example shows that the VMCSCF results are closer to the VFCI calculation, which of course is a result of the modal relaxation and the rather small correlation space.Note that we have chosen a small number of configurations in these calculations by purpose in order to show this effect.An important issue within the framework of (V)MCSCF methods is the convergence behaviour of the different modal relaxation schemes.As our VMCSCF approach outlined above is based on the pioneering working of Hinze for electronic structure methods, 94 it is to some extent surprising that the calculations do converge at all and the question is why?We believe that this question is the same as why statespecific VSCF calculations do converge for highly excited states and do not fall back to the solution for the vibrational ground state.Two differences between the electronic and vibrational structure calculations are rather obvious: (1) due to the simpler shape of the vibrational modals in comparison to the rather complicated structure of molecular valence orbitals, the optimization process in the vibrational structure methods may lead to smaller changes in the modals (relaxation contributions).In other words, the starting guess lies within the convergence radius of the VSCF and VMCSCF modal optimizations.For that reason a state-specific VSCF calculation usually needs only 4−5 iterations, while a corresponding Hartree-Fock calculation needs more

INFRARED INTENSITIES
Once the vibrational wave function has been determined, the calculations of any property is straight forward but requires the evaluation of property surfaces.This of course can be a very demanding task depending on the nature of the property.The calculation of anharmonic infrared intensities requires the determination of dipole surfaces.As the evaluation of molecular electric dipole moments necessitates to solve the coupledperturbed equations in the electron correlation calculation, the computational time roughly doubles for each ab initio calculation.In the worst case this may prohibit a calculation.However, most experiments provide intensities qualitatively rather than quantitatively.Very often they are just provided in terms as low, medium or strong.For that reason harmonic intensities provide a fair estimate for the intensities of the fundamental modes once Fermi resonances are of minor importance.This picture changes completely for vibrational overtones and combination bands.As each electronic structure calculation relies on a Hartree-Fock calculation, which provides the electric dipole moments for free, the determination of RHF dipole surfaces does not take any time except for the fitting.Consequently, estimates of infrared intensities of vibrational overtones and combination bands can be obtained without significant computational effort.The calculation of infrared intensities is based on the integrated absorption coefficients (see for example Refs.96 and 97): i n  and f n  denote the wave functions in the initial and the final vibrational states, respectively.ΔE if is the energy difference between these two states and Δn if is the temperature dependent difference in the fraction of molecules in these states, which is assumed to be one.Note that, once the PES and thus the wave function has been determined by accurate methods, it is only the dipole operator, which is provided at low level, while the wave functions Ψ and the energy difference ΔE if have been obtained at high computational levels.As an example Table 3 shows computed infrared intensities for selected vibrational overtones and combination bands of formaldehyde obtained at different levels of theory. 69,98The intensities obtained from RHF dipole surfaces clearly deviate from those determined at the QCISD(T) level, but may serve as a good guess obtained from negligible additional computational costs.

VIBRONIC TRANSITIONS
The calculation of vibronic transitions is formally a straight forward task once vibrational wave functions have been obtained for the electronic ground and excited states of the molecule of interest.0][101][102][103][104][105] The reason for this arises from the increasing number of basis functions and the Duschinsky rotation, which essentially leads to an exponential scaling of the computational cost for the  Franck-Condon (FC) integrals.Besides this the number of FC integrals becomes huge and thus prescreening schemes are needed to reduce this number to a managable amount.However, established prescreening schemes 99 have focussed primarily on the harmonic approximation.4][105] We focus on the former approach.The FC integral between two electronic states (quantities of the excited state are indicated by a prime, while those belonging to the reference state appear without the prime) essentially can be reduced to a summation over multi-dimensional overlap integrals of the distributed Gaussian basis functions χ μi , i.e.

   d
This integral can be solved as outlined by Doktorov et al., 107,108 but its solution becomes very demanding with increasing size of the system.Moreover, based on VCI wave functions the number of the integrals N I increases for non-linear molecules exponentially with with N States being the number of vibrational states, N Conf being the number of configurations in the VCI calculation and N Basis the number of basis functions and N the number of atoms.The distributed (mode-specific) basis functions are given as with the exponential parameters A μi and the offsets q μi .The normal coordinates of the two different electronic states are connected by the Duschinsky relation, 106 i.e.
   q Sq d (34)   with S being the Duschinsky rotation matrix and d the deplacement vector between the two reference structures.In contrast to other implementations, in our program S includes the Eckart transformation matrix.
and thus allows for screening of the integrals prior to their evaluation.In Table 4 we show the impact of the integral screening due to this threshold on the Franck-Condon factors for the photodetachment process from HS 2 − ( 1 X A'  ) to HS 2 ( 2 X A"  ).Screening of the integrals led to an acceleration of the FCF evaluation by a factor of 2 hardly without loss in accuracy.As can be seen from equation 32, the number of integrals to be evaluated depends strongly on the number of basis functions.As low lying modals and thus low lying vibrational states as well need fewer basis functions for a proper description than high lying ones we introduced a variation of the number of basis functions in order to neglect those, which do not contribute to the individual modals.This has significant impact on equation 32 and reduces the computational effort by factors ranging from 2−5.The impact of this basis reduction on the final Franck-Condon factors is shown in Table 5.As the deviations due to screening and the neglect of unimportant basis functions are negligible and can be controlled by simple thresholds, the evaluation of Franck-Condon integrals can be accelerated by a factor of 4 at least.However, these schemes do not act on the 2 States Conf N N scaling in equation 32 and thus the primary goal of future work within the framework of calculations beyond the harmonic approximations must be the reduction of this factor.

SUMMARY AND CONCLUSIONS
In this feature article we have discussed some specific aspects relevant to vibrational structure calculations based on vibrational self-consistent field and vibrational configuration interaction calculations, which is a strongly evolving field of research.Of course, the topics discussed here are just a small selection of what has been done in this field over the last years, but cover some aspects which have barely been discussed in the literature.We believe that in the near future it will be possible to compute molecules with up to 20 atoms within reasonable time and thus calculations as presented in this study will belong to standard applications.Modern multi-core architectures significantly pushed the limits in that direction and request the adaption of current software.For that reason robustness, black-box nature and extended applicability of the codes will be the primary challenge once computational demands are of less importance.

Figure 2 .
Figure 2. Percentage of selected configurations for the vibrational ground state in VCISDTQ calculations for (LiF) n clusters.

Table 1 .
95pact of local approximations and the modeling of 3D terms on the fundamental modes of Li 2 F 2 .All data were obtained from VCI calculations based on the full Watson Hamiltonian95and are given in cm−1

Table 2 .
Comparison of VCI and VMCSCF results with respect to vibrational full CI data for chloroform.All quantities are given in cm−1

Table 3 .
Comparison of anharmonic infrared intensities (given in km/mol) for selected overtones and combinations bands of formaldehyde obtained from QCISD(T) and RHF dipole surfaces.

Table 4 .
Franck-Condon factors for the photodetachment process from HS 2 Integrals with a norm-estimate being smaller than 10 −8 were neglected.