Perspectives from Quantum Chemistry for molecules and nanograins with astrophysical Interest

Quantum Chemistry methods offer powerful tools for the computation of various properties of molecules and nanoparticles or clusters such as structure, thermodynamics, infrared and electronic spectroscopy, in some cases on a routine basis. However, neither the accuracy nor the applicability of the various methods are uniform, and may depend strongly on the nature and the size of the compounds and the desired properties. The ab initio methods for approaching electronic structure can be classified into two classes. The first class includes wavefunction-based methods, namely post-Hartree Fock schemes in the framework of Configuration Interaction or Coupled Cluster schemes which can now be used for molecules containing up to a few ten atoms (for single geometry calculations) and are likely to provide accurate results whenever applicable. The second class of methods is that of density-based methods which cover systems between a few tens up to a few hundreds of atoms. It is often of broad applicability and satisfactorily accurate in many cases. For more extended systems or extensive calculations, tight-binding type approximations or evolved force-fields (based on bond-orders for instance) are also quite promising. In this overview, we will briefly remind the basics of standard Quantum Chemistry methods and provide a survey of present trends which should be of interest for astrochemistry simulations at the atomic scale.


INTRODUCTION
More than a hundred molecules have now been clearly assigned in the interstellar medium (ISM).They range from small diatomics and triatomics (neutrals, cations, anions and radicals) to larger polyatomic molecules, such as hydrocarbons and even small aminoacids.It has been put forward that grains play a considerable role in the chemistry of the ISM, presumably responsible for the ratio of molecular/atomic hydrogen for instance.Focusing on the small size range, very small grains (VSG) are presumed to include silica, carbonaceous systems, possibly ices.As an example of the challenge, the PAH (polycyclic aromatic hydrocarbons) hypothesis has been invoked to explain the presence of Aromatic Infra-Red bands (AIB) in the mid-IR range.The size of the molecules, the arrangement of the cycles, their hydogenenation and their ionization states offer a variety of parameters to investigate.However until now, no single PAH molecule has been assigned in the ISM.Assumptions of quite large PAH molecules, with a few hundreds to a thousand carbon atoms, have been proposed, as well as clusters of such molecules and nanograins.Other assumptions are proposed for the composition (complexation with heteroatoms such as Si or Fe) or morphology of such grains.A variety of possibilities are opened due to the versatile nano-organization of carbon (molecules, fullerenes, diamond nanoclusters, graphene sheets, amorphous clusters), including grains with clean or polluted surfaces.The filiation between small and larger complexes, via for instance carbon clusters or hydrocarbon chains, is still in question.The knowledge of IR spectroscopy features is of course essential to fingerprint and assign various a e-mail: fernand.spiegelman@irsamc.ups-tlse.frThis is an Open Access article distributed under the terms of the Creative Commons Attribution-Noncommercial License 3.0, which permits unrestricted use, distribution, and reproduction in any noncommercial medium, provided the original work is properly cited.molecules and complexes.This clearly relies on laboratory astrophysics, namely molecular spectroscopy characterization, and also on molecular simulation.Besides spectroscopy, reactivity is essential to understand formation and filiation.However even for molecules, reactivity, formation mechanisms and rates, fragmentation and recombination processes, interaction with the UV radiation are not as well known as structural and/or spectroscopic features, except for the smallest systems.Fewer results are available on complexes and even small clusters.

EPJ Web of Conferences
Quantum Chemistry offers precious tools to investigate the properties of a variety of species ranging from molecules to nanograins.One of the challenges is to bridge the gap between small molecular systems with around a dozen of atoms and nanoparticules or nanograins with a few hundreds or thousand atoms, described at the microscopic level (see for instance the system in Figure 1).The understanding of the gas-grain interaction, namely the possibility to simulate reactivity at or with surfaces is also needed.Accuracy and size seem to be contradictory goals in the description of electronic structure.The situation may also vary depending on the purpose: single point calculation, local or global geometry optimization, normal mode analysis, Monte Carlo or Molecular Dynamics simulations, Quantum molecular Dynamics.For instance if one refers to dynamical properties, the concept of precalculating potential energy surfaces breaks down beyond 10-20 degrees of freedom, and one has to consider onthe fly methods where the electronic structure and forces are computed at each step of the dynamics.Moreover in order to achieve representative statistical sampling, for instance via the initial conditions at a given temperature, or to understand the thermal behaviour, one needs to repeat energy calculations several 10 6 -10 8 times or more.The most advanced ab initio methods for electronic structure based on the determination of the wavefunction, although in principle exact, show drastic size dependance.It is a challenge to approach even medium size systems (10-100 atoms, depending on the nature of the atoms and on the number of the electrons explicitely considered).Their use in extensive molecular dynamics is highly computer-time demanding.Alternatively, methods based on the Density Functional Theory (DFT) are relatively more computer-efficient and routinely used, but despite recent progresses, show some drawbacks which prevent their direct use to study some specific, however important, situations.This is the case for instance of weak intermolecular interactions (such as Van der Waals forces), or charge resonance interactions which allow the interplay between electron delocalization, polarization and electron-electron correlation in charged molecular clusters.Multi-open shell systems and the correct account of spin multiplicity is also a difficulty.
The present survey aims at a brief description of the standard methods used to investigate the electronic structure of molecules, grains and surfaces and an introduction to some current trends and advances in the field of Quantum Chemistry, underlining some relevant theoretical challenges that have to be taken up.The overview is mainly focused on gas phase systems, but references to methods relevant for periodic systems are also given.This survey intends to provide a brief insight in the following approaches -explicitely correlated methods, -ab initio wavefunction-based (WF) methods developed in fragment-localized orbitals, which may be expected to achieve computational time showing linear scaling as a function of the size of the systems, -progresses in Density Functional Theory (DFT), and in particular the definition of new functionals for the description of dispersion forces,

06003-p.2
Chemistry in Astrophysical Media, AstrOHP 2010 -hybrid methods combining wavefunction-and density-based methods, also called short-range/longrange (sr/lr) methods, which open a way to correct the treatment of dispersion forces, resonance and open-shell systems beyond the usual Kohn-Sham single determinant-based formalism, -tight-binding type approximations of DFT and their developments, -reactive potentials, which tend to incorporate many-body quantum effects under the form of classical force fields.

STANDARD METHODS: WAVEFUNCTION VERSUS DENSITY FUNCTIONAL APPROACH
The electronic problem in its time-independent Born-Oppenheimer formulation is controlled by the electronic Schrödinger equation R stands for the coordinates of the nuclei (classical parameters in the BO approximation), and r describes the quantum positions of the electrons.E(R) is the potential energy surface (energy landscape) describing the forces between the nuclei, while (R, r) is the electronic wavefunction.
There are essentially two ways to approach the solution of this equation, at least for the electronic ground state.1.The traditional approach (see the article of Dabia Talbi in this volume) is to undertake the direct determination of the eigenvalue(s) and multi-electronic wavefunction(s) (WF) of the molecular system under study [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15].Even in the case of the ground state, this is a difficult task which can only be reached approximatively.The usual way is to develop the multielectronic wavefunction as a configuration expansion, where k is either a Slater determinant or a spin configuration (minimal combinations of determinants, eigenfunction of the total spin S 2 ).This expansion defines the general Configuration Interaction (CI) framework [16].Using such expansion, one can alternatively express the eigenvalue problem as a minimization problem searching the wavefunction expansion which minimizes the expectation value with a normalization constraint imposed upon .If this combination is limited to a single determinant 0 , one recovers the well-known Hartree-Fock (HF) solution [17], which reduces the initial multi-electronic problem to the determination of singleelectron spinorbitals i defining 0 This approximation is a mean-field approach, as the spinorbital attributed to a given electron is determined in the field of the mean Coulomb and exchange (v H F x ) contributions of the other electrons in addition to the one-electron kinetic − 2 and nuclear attraction contributions v ne of the Hamiltonian.This equation is solved iteratively, which at convergence yields the self-consistent field (SCF) Hartree-Fock solution.Today, the HF-SCF solution essentially serves as a conceptual background or an initial step for more elaborate methods to approach electronic structure with WF methods.One should nevertheless mention that the HF theory is based on the concept of one-electron molecular orbitals either occupied or empty (virtual) in a reference determinant optimized within a SCF procedure, concept that also underlies the Kohn-Sham approach to Density Functional Theory (see below).For finite systems, various post Hartree-Fock methods are now routinely available which allow to treat electron-electron interaction beyond the mean-field Hartree-Fock method.Most of them imply truncations/partitions of the CI space.In many treatments, CI truncations are defined according to selections based on classification of the molecular orbitals.Moreover, in a finite CI, the molecular orbitals may be simultaneously optimized.According to their participation in the excitation processes, the molecular orbitals are then labelled as frozen (not involved, not optimized), inactive (not involved, optimized), active (involved and optimized), external (possibly involved in larger CI or MRCI, without further orbital optimization).Full CI within a finite subspace of molecular orbitals defines a Complete Active Space (CAS) calculation.Such subspace can be further reduced as in the Restricted Active Space (RAS) version which only involve some types of excitations in the active space.Truncation may actually also be defined according to the excitation degree reached from single, double, triple, quadruple... excitations operator T kl... ij ... with respect to the reference Hartree Fock determinant 0 , expressing the above the CI expansion as N is a normalization constant.The CI space reduction can also be guided according to quantitative criterions such as determining perturbative estimates of the weights of the determinants and retaining only the largest ones in the final CI expansion.Finally, hierarchical partitions can also be defined with respect to the type of CI treatment (for instance variational vs perturbative).It is not the scope of the present review to explain in details all these methods.Let us only mention some of the most popular ones.One class is that of perturbative post Hartree-Fock methods, such as the second order or fourth-order Møller-Plesset Perturbation Theory [18][19][20], respectively MP2 or MP4 which are efficient for non-degenerate situations.Another important scheme, variational, is provided by the multiconfigurational versions of the self-consistent field approach [21] (MCSCF), in which a finite CI is combined with SCF optimization of the orbitals.This general scheme includes various versions such as the Complete Active Space Self-Consistent Field (CASSCF) scheme or the Restricted Active Space Self-Consistent Field (RASSCF) scheme.One can still improve such calculations by subsequent Multi-Reference CI (MRCI) in which all single and double excitations from a multireference space are included and so is the part of correlation brought by the excitations towards the external orbitals.Finally one may cite the multi-reference perturbative methods CIPSI [10], CAS-PT2 and RAS-PT2 which combine perturbative approach of the external space with variational treatment of the main one (CAS or RAS for instance [2,13,14]).Those methods allow for the determination of excited states.

06003-p.4
Chemistry in Astrophysical Media, AstrOHP 2010 An alternative attractive way is offered by Coupled Clusters (CC) theories [3,4,[22][23][24][25], which use an exponential expansion of the wavefunction.For a single reference, the CC expansion is with CC methods often provide highly accurate description of the ground state.They are particularly suited to describe weak interactions because of nice size-consistency and additivity properties consecutive to the exponential form of the wavefunction.2. The second approach [26][27][28][29][30] relies on the calculation of the properties, not from the multi-electronic wavefunction, but from the total electronic density (r) which depends on the three coordinates of R 3 only, and hence is a much simpler function.Density Functional Theory (DFT) has strongly developed since its formulation by Hohenberg and Kohn [26] and the mean-field type formulation provided by Kohn and Sham [27].The theory relies on two theorems initially derived by Hohenberg and Kohn [26] and generalized by Levy and Lieb [31,32].The first theorem states that the energy E 0 of the ground state of any electronic system is a functional of its density 0 .
= min = F [ 0 ] + v ne (r) 0 (r)dr (10) F [ ] is an universal functional for electronic systems and v ne is the one-electron potential exerted by the nuclei.The second theorem is a variational principle with respect to the density, enabling the determination of the energy E 0 of the ground state.
While the two theorems provide an existence framework and a variational principle to find the density, a crucial problem stems from the fact that the functional E[ ] (namely F [ ]) is not known.Kohn and Sham [27] made a decisive advance by establishing the mathematical similarity of the equations for the actual system with that of a fictitious system of non-interacting electrons in the same nuclei field.The wavefunction of the latter is known to be a single Slater determinant with electrons in individual orbitals i , solutions of the one-electron Schrödinger equation and associated with a total density (r) The density of the actual interacting electrons can be derived in a similar way, provided that v ne is complemented by an effective potential v KS incorporating the Coulomb, exchange and electronelectron correlation contributions.v KS is thus usually decomposed into a Coulomb term and the 06003-p.5 This defines the Kohn-Sham equations, which are exact, provided the exchange-correlation potential v xc (or the functional E xc [ ]) is known.It is seen that Density-Functional Theory also provides a mean-field type SCF equation in the sense that each electron is determined in the Kohn-Sham field.Correlation (electron fluctuations) is in principle fully accounted for in DFT, while Hartree-Fock only accounts for exchange correlation between electrons with identical spins.Unfortunately, the exchange-correlation functional is only approximated.As a consequence, an important drawback is that the self-interactions of the electron which appears in the Coulomb term of the KSequation (since integration is performed over total densities), is very unperfectly cancelled by the exchange contribution (due to the approximation in the exchange functional).There have been numerous proposals, such as for instance hybrid functionals (mixing exchange functional contribution with exact HF exchange), or explicit corrections [33][34][35][36][37][38][39][40][41] to self-interaction errors (SIE).
Let us just now briefly mention that in practical calculations, most methods (HF, post HF and DFT) are based on expansion of the molecular orbitals i in terms of single electron basis functions.In the case of Quantum Chemistry codes, the Linear Combination of Atomic Orbitals scheme (LCAO) is generally used, spanning the orbitals on basis functions l localized on the atoms a, where the basis functions l can be either Gaussian Type Orbitals (GTO) l (x, y, z)r l exp(− r 2 ) most frequently used or sometimes Slater Type Orbitals (STO) l (x, y, z)r n exp(− r).The l functions provide the angular dependance (for instance cartesian combinations indexed by of spherical harmonics Y lm ).In the case of extended materials such as periodic solids, 2D or 1D nanostructures, the symmetry is conveniently described in the reciprocal space.Efficient implemention may involve plane waves (PW) basis functions or Bloch type function basis functions where the l are local basis functions (possibly Gaussian or Slater functions) centered on the atoms of the elementary crystal cell, and T is a periodic translation vector in the direct space.In the DFT framework, one may also mention schemes expressed in real space 3D grids [42][43][44][45], combinations of plane waves with localized function (Augmented Plane Waves or APW scheme), or wavelets [46] which are promising for spatial multi-scale representations.In the WF approach the one-electron basis set determines the number of molecular virtual orbitals and the spatial range of possible excitations.The dimension of the total CI Hilbert space increases with the number of electrons and the number of accessible orbitals.The quality of the results clearly depends on the extension (flexibility) of the one-particle basis set, of the size and the content of the configuration interaction space that is actually reached in a given method, and on the accuracy of the method itself (variational vs perturbative for instance) as schematized in Figure 3.  Let us mention some possible applications of the various methods.The WF methods such as variational CI or CC are generally limited to finite systems.CC methods for instance provide in general very accurate descriptions of the ground state potential energy surfaces, especially when combined with explicite correlation factors (see below) and have been used to determine highly accurate ro-vibrational spectroscopy or reaction rates [48][49][50][51].To address gas-grain interaction, one may replace a presumably infinite grain surface with a more or less extended but finite cluster (cluster model for a surface), as shown for instance in figure 4 where a finite 2D cluster is used to model the interaction of beryllium nanochains on infinite graphene or graphite surface.
DFT has become an efficient and routinely used method to investigate structural, spectroscopic and electronic properties of molecules and materials.It has been extensively applied to investigate the structural and IR spectra (in the harmonic approximation) of large PAHs molecules up to few hundreds of atoms [52,53].From higher order energy derivatives, anharmonic effects in IR spectra are also available.In the case of periodic systems, properties such as chemisorption, physisorption, reactions of molecules or adsorbates (clusters) on or with surfaces can also be obtained, including surface and undersurface relaxation properties.Pathways for the adsorption of atomic and molecular hydrogen on various surfaces have been determined with either LDA (local density approximation) or GGA (generalized gradient approximation) functionals, allowing the determinations of various properties.
Finally, let us note that DFT is also used in on-the fly molecular dynamics contexts, such as for instance in the powerful Car-Parinello classical molecular dynamics method (CPMD) [54].Nonharmonic IR spectra can be actually be derived MD [55,56].Quantum versions of CPMD such as the Path Integral scheme [57] also incorporate quantum vibrational dynamics.

CORRELATION AND WF SCHEMES WITH LOCAL ORBITALS
Generally, the orbital space of a system can be separated into valence orbitals qualitatively directly involved in the chemical bond, and orbitals external to the valence shell (resulting from the extension of the basis set).The external orbitals, although not qualitatively essential for the description of the bonding of the valence states are necessary to achieve quantitative accuracy in the description of electronelectron correlation.They are also needed to describe non-valence excited states.This hence determines a partition of the CI space spanned by the valence configurations (or determinants) possibly chosen as the active space (whenever not too large), and the external subspace (generally much larger).The (until now) standard and most widely used post-Hartree-Fock schemes are achievied using symmetry molecular orbitals possibly delocalizing over the whole molecule, originating from HF-SCF or MC-SCF procedure.However Configuration Interaction can be conducted either using delocalized MO's or any other (orthonormalized) MO set.Obviously, the partition between valence and external subspace remains valid in either case, provided that a unitary transformation is achieved within each subspace independantly.This is the essence of the early dual interpretation of molecular bonding, either in terms of molecular orbitals with paired electron, or alternatively following the valence-bond formulation of Heitler and London [58,59] in terms of covalent vs electron-transfer configuration.One can furthermore define bond orbitals in polyatomic systems, or fragment-localized orbitals which for instance may be suited for molecules in molecular clusters or functional groups.
This dual description will be illustrated via the example of the simplest molecule with two active electrons namely H 2 (Fig. 5).While very basic, this example provides some quite general insight.For sake of simplicity, the question of the overlaps will be avoided in this essentially qualitative analysis.The atomic valence orbitals for the H 2 molecule are 1s a and 1s b (for the present purpose, one can refer to the 1s orbitals of atomic hydrogen, but more flexible localized descriptions are possible, involving optimization and relaxation in the molecular field) [60].This generates two molecular orbitals, a bonding one ( g , in-phase gerade combination) and an antibonding one ( u , out-of phase ungerade combination).Hence, for singlet states and total spin projection M S = 0, the complete active valence space (CAS) one can either be generated from localized orbital description or alternatively from the set of molecular delocalized orbitals (equivalent while smaller thanks to symmetry benefits) , label the electron spin projections.The former description allows to analyze electron localization, discriminating configurations either having electrons on different sites (covalent forms), or two electrons on the same site (H + H − electron-transfer configurations or zwitterionic forms).The valence CI governs the interplay between the various localized configurations and thus fluctuation of the two electrons over the molecule, defining the so-called non dynamical correlation (Fig. 6).The same interplay can be obtained with either scheme, using valence localized or symmetry orbitals, the analysis in the latter case being however less straightforward.
Non-dynamical correlation (totally or partially accounted for in CASSCF schemes, depending on the choice of the active orbitals), is obviously important to describe the covalent/charge transfer weights at equilibrium.Furthermore, it turns out to be crucial to properly describe bond breaking and formation which occur along reaction paths and in dissociative process: the breaking of a single bond requires a minimal CI involving two bond orbitals (a bonding one and an antibonding one).For a triple bond, one needs the determinantal set defined by six orbitals the , * , x , * x , y and * y .The complete explosion of a cluster may require the inclusion of full valence CAS-CI or CAS-SCF.Thus, in general, the nondynamical correlation refers to important combinations between quasi-degenerate configurations, which deserve to receive variational treatment, at least in some geometrical regions of the potential energy surface.In H 2 , the non-dynamical correlation involves fluctuations between electronic distributions on different atoms and is therefore a long-range contributions.Note that in the case of incomplete many-fold atomic valence shells (transition metals), short range intra-valence contributions may also appear.The interplay is governed by the relative energy ordering of the various configurations and the magnitudes of their couplings.The energies depend on the electron-electron Coulomb penalty for electron pairs (large for on-site pairs, smaller for inter-site pairs) on the resulting net electrostatic balance brought by the fluctuation with respect to the neutral situation, and on the exchange contributions between electron pairs with same spins (the most strongly localized, the largest).As for the couplings, the largest ones originate from kinetic energy two-center integrals which favor delocalization.
After analysing the electron-electron correlation in the valence shell, let us examine fluctuations of the electron distributions of valence states out of the valence shell.This fluctuation determines the dynamical correlation (Fig. 6).Note than in the case of H 2 , we may start by considering the interaction of determinants such as |  excited states.The coupling matrix element with the ground state corresponds to the interaction between two transition dipole moments located on the two sites, and is therefore a 1/R 3 interaction.Since this coupling occurs between non-degenerate states, the effect of the excited configuration may be accounted for perturbatively and provides a 1/R 6 behaviour contribution to the ground state energy which is characteristic of Van der Waals dispersion forces.Dispersion forces are caused by the interaction between fluctuating electrons at different sites: it is a case of long range electron-electron correlation.One may moreover also consider the influence of double excitations on the same site towards higher shells.This lowers in particular the energy of configurations with two electrons on the same site (H + H − ).determine their on-site radial correlation.This part of the dynamical correlation corresponds to intraatomic (short range) correlation.Its effect may be to change the relative positioning of the valence configurations, for example the gap betwen covalent and charge-transfer configurations, and hence affect the configuration mixing in the valence space.In general, non-dynamical correlation (possibly large), can be accounted for perturbatively since it results from numerous excitations implying however significant energy differences.This provides a justification of CAS (valence)-PT2 (external) schemes.
Let us note however that the situation can be more delicate in some cases due to the ambiguity of the definition of the valence shell.For instance intra-atomic short range correlation can also be observed in the non-dynamical contribution, for instance in Beryllium compounds where both 2s and 2p shells must be considered as defining the valence active space.2s 2 a → 2p 2 a excitations then belong to the non-dynamical valence CI space.In this case, their contributions appear as non-dynamical longrange contributions, and actually their account must more safely be taken into account variationally, inducing larger molecular stability than usual van der Waals interactions.The obligation to consider several atomic shells as defining valence is also true for metals or mixed valence compounds.Thus, the concepts dynamical-non-dynamical vs short range -long-range are actually qualitatively fruitful but non-univoquely defined.
Not only achieving CI within localized scheme is relevant from a conceptual and interpretative point of view, but is also of significant interest to develop new numerical schemes [61,62].One of them consists in working with atom-or fragment-localized orbitals to perform post-Hartree Fock calculations (Fig. 7).The valence CAS-CI itself very quickly overpasses the present limits.For instance, the valence space of the C 3 molecule (considering the 12 molecular orbitals generated by the 12 electrons originating from 2s and 2p, and M S = 0) has a size of 853776 (not considering symmetry).In a delocalized framework, it is essentially possible to select the active orbitals for the main CI (MCSCF or CASSCF) according to their energetics.The interactions spread over the whole system, so that the two-electron integrals to be calculated may involve sites far away from each other.Correlatively, the CI size increases drastically with the the density of the one-electron levels.This provides computer-time scaling with sometimes quite high powers of the number of electrons, typically N 3−6 depending on the method.One may also sometimes face convergence problems in the MCSCF codes due to the difficulty to isolate an active subspace, due to degeneracy with higher configurations (so-called intruder-states).In order to

06003-p.10
Chemistry in Astrophysical Media, AstrOHP 2010 extend the WF methods towards large molecules, clusters and nanograins, one may take the advantage of working within a localized framework.This brings several advantages : 1. a better definition of the active space; chemical bonding or bond breaking are often local processes, 2. a possible classification of the orbitals according to spatial localization, and not necessarily energy, 3. only treating explicitly the region of interest, keep frozen the rest of the system not involved in the system, 4. avoiding calculation of two-electron integrals over the whole system.Thus, if each regional part receives an independent treatment of correlation, one may achieve linear scaling computational behaviour [61][62][63][64].This opens the use of CI to larger systems, as illustrated in Figure 8, where the interaction of Na with a local orbital of C 60 is shown [65].Such calculations can be performed at various post-Hartree Fock levels, namely MP2, CI, and are currently being developed with MCSCF and MRCI methods.
Finally, let us mention that while the ab initio standard treatment of periodic systems have for a long time essentially be confined to HF or LDA, combinations of DFT with the N-body problem are currently proposed [192][193][194].Closer to the quantum-chemical framework, recent developments now also incorporate the MP2 methods, based on the use of Wannier functions [66], which provide an orthonormal set of direct space cell-localized functions [67][68][69][70].

EXPLICITELY CORRELATED WF METHODS
An important class of methods to be mentionned in the WF category is that of explicitely correlated methods.This class of methods stems from the fact that the Coulomb hole which describes the behaviour of two electrons avoiding each-other via the Coulomb interaction cannot be analytically correctly reproduced via any determinantal or multi-determinantal wavefunction.Actually, the quantitative convergence of CI for a two-electron atom for instance is very slow, in particular as a function of the l angular quantum number of the basis set used to generate the determinants.The origin of this slow convergency is the local behaviour of the two-electron wavefunction in the limit where the interelectronic distance tends to zero, which, following Kato [71], reads meaning that the wavefunction has a linear dependance in the vicinity of r 12 → 0, namely exhibits a cusp.No determinantal neither multi-determinantal wavefunction can exhibit such analytical behaviour (showing always a zero first derivative at this limit), as shown in Figure 9. Thus, using the concepts of the previous section, it is seen that the explicitely correlated methods aim at curing the very slow convergence of the short-range correlation with high l atomic functions.In order to overpass this difficulty, following Hylleras [72] and Slater [73], Kutzelnigg and coworkers [74] proposed to introduce explicit pair correlation factors complementing the usual expansions of correlated wavefunctions, namely where Q 12 is a projector necessary to avoid redundancy with excitations already accounted for in the standard part of the expansion.The pair correlation factor f (r 12 ) initially introduced by Kutzelnigg [74] was restricted to the first order term r 12 .More recently, an exponential correlation factor exp(− r 12 ) was proposed [75,76] for an better account of the cusp.Those explicitely correlated expansions have been implemented within different approximations of the CI expansion, for instance in the perturbative schemes [77] such as MP2-R12 , MP2-F12, or in various approximations of the Coupled Cluster scheme [78][79][80][81] CC-R12 or CC-F12.Those methods have been used for small polyatomic molecules of astrophysical interest, such as the H 3 O + molecular ion [82], the calculation of the 9-dimensional potential energy surface of the H 2 O-H 2 system (used to determine the inelastic collisionnal transfer) [83,84], that of the CO-H 2 colliding system [85], accurate investigation of the ammonia molecule [86].When extensive basis sets and proper extrapolation procedures are used, such methods may reach extremely high accuracy in the description on the potential energy surfaces of some systems (claimed accuracy of the order of 1 cm −1 ).

COMBINING CONFIGURATION INTERACTION WITH DFT
Hybrid WF/DFT methods were first developed by the groups of Stoll and Savin [87][88][89][90].They rely on the combination of density functional theory methods to account for the short range correlation with

06003-p.12
Chemistry in Astrophysical Media, AstrOHP 2010 wavefunction methods to account for the long range correlation.First, the electron-electron interaction is separated into two parts.
V ee = V SR ee + V LR ee (22) both terms being defined as follows and erf is the standard error function and a coupling parameter between DFT and WF treatments.V SR ee will be treated via a Kohn-Sham type treatment while V LR ee will receive a WF type solution.The softness of the separation is controlled via the parameter ( → 0 yields a pure DFT scheme, while → ∞ represents a pure WF scheme) Emphasizing the wavefunction description for the long-range part and the density description for the short-range part the variational principle can be rewritten as In the last contribution appears a functional with respect to the short-range electron-electron interaction, which is partitionned into a Coulomb term 1 − erf ( r 12 ) r 12 (r 1 ) (r 2 )dr 1 dr 2 (26) and a short-range exchange-correlation trem E SR xc .The variational principle can thus be expressed as ) which describes a wavefunction variational problem with long-range coulombic interaction only, in the field of the short-range interaction taken as a functional of the density.Minimizing E 0 is equivalent to the search of an eigenfunction LR for the following hamiltonian ) where the two last terms are one electron-potentials derived from the short-range functional contributions, namely The CI solution and the density are to be derived consistently.remains an arbitrary parameter which is generally taken close to 0.4.Various WF/DFT combinations can be explored using different CI schemes and different functionals, the scope being to achieve more efficient convergence towards the correct results without using the most sophisticated (expensive) DFT functional/CI method combinations.

TREATMENT OF LONG-RANGE INTERACTIONS WITHIN DFT
The basic problem within DFT is that the exact functional is unknown.Various functionals [91][92][93][94][95][96][97][98][99][100] have been derived beyond the local approximation [101], either involving the functional gradient or hybrid functionals including an amount of exact exchange.The search for new functionals including non-covalent interactions is very active [40,[102][103][104][105][106][107][108].Recently proposed functionals have shown to provide more satisfactory description of dispersion forces for a wide range of van der Waals 06003-p.13systems including benzene and polycyclic aromatic hydrocarbon PAH dimers.For instance, the M06-2X functional [102] is a highly non local functional with an enhanced contribution of non-local exchange completed with a meta-GGA correlation functional and was shown to provide reasonable results for noncovalent systems.In a more empirical scheme, the DFT energy can be corrected by an additional atom pair contribution to account for London dispersion through a −f (R)C 6 /R 6 term involving a short range cut-off f (R) [109][110][111][112][113][114][115][116].Table 1 compares the binding energies for selected isomers of coronene dimer obtained with Symmetry-Adapted Perturbation Theory-DFT (see also below), which can be considered as an accurate reference, with empirically corrected DFT (BLYP-D [115], PBE-D [115], B97-D [117], PBE1-D [118]) and with the M06-2X functional [119].The binding energies are in qualitatively and sometimes quantitatively good agreement with those of the reference SAPT(DFT) calculation [120].

EPJ Web of Conferences
The geometrical parameters (although not shown here) also agree with the SAPT results.

TIGHT-BINDING APPROXIMATION OF THE DFT THEORY
Tight-binding approximations have been rationalized by many authors [122][123][124][125].One of the most systematic formalisation is that of Density-Functional based Tight-Binding (DFTB) and self-consistent DFTB derivation [126][127][128][129][130], which relies on the following approximations 1. the development of the fluctuation of electronic molecular density about the atomic densities 0 and the neglect of terms beyond second order 2. the neglect of three-body terms 3. the expression of the Kohn-Sham equation in a minimal valence basis set 4. the discretization of the molecular density on the atomic sites, through the use of atomic charges fluctuations q a (usually Mulliken charges [131]) with respect to neutral atoms The expression of the energy becomes a,b∈ atoms ab q a q b (31) V R accounts for the Coulomb+Pauli repulsion of the core electrons and possibly other scalar corrections (dispersion terms for instance, see results c-DFTB-D in table 1 H 0 ( 0 ) is the one-electron operator at the reference density and ab represent the on-site and inter-site Coulomb interactions incorporating formally second order exchange-correlation contributions.

06003-p.14
Chemistry in Astrophysical Media, AstrOHP 2010 5. the precalculation of the two-center integrals as functions of interatomic distances (until now using DFT/LDA calculations) and their tabulation via analytical functions.All those approximations make the calculation of the energy and its gradient particularly efficient while keeping a quantum description.Also, the presence of the second-order term makes the method able to account for charge fluctuations and even various ionization states (for instance see a variant application to C 2+ n clusters [132]).As DFT, DFTB is a self-consistent method (SCC-DFTB [129]), however with respect to the atomic charges fluctuations instead of the density.In concern with the repulsive part, several parametrizations do exist with different reference data (bulk properties, biological systems, etc..).As in many parametric methods, partial reparametrization may sometimes be needed.For instance, the derivation of accurate IP differences is crucial in order to describe the interaction of a silicium atom with a PAH cation, since the Si-PAH complexes are characterized by the significant influence of the charge transfer situation at dissociation (Si + + PAH) only lies with 3 eV above the (Si + PAH + ) dissociation limit.The balance between covalent and charge transfer character depends on the former.Figure 11 illustrates the transferability of a simple correction (shifting atomic energies) to the PAH-Si difference for increasing size PAH's [56].
Fe-PAH complexes have also been investigated [133] since they appear as other tentative candidate systems for some of the near-IR AIBs. Figure 12 shows the comparison between Density Functional pathway for Fe at the surface of a coronene molecule and that determined with DFTB after some parameter adjustment.This shows that the DFTB PES is realistic.The inclusion of DFTB in classical on-the-fly molecular dynamics simulations allows to investigate the temperature behavior of the IR spectrum via the calculation of the Fourier transform of the dipole autocorrelation function.This technique is used with DFTB as it has been used with DFT for instance in the context of the Car-Parrinello method [54,55].The advantage of DFTB stems from the computer efficiency and the higher sampling that can be achieved.Figure 14 shows the temperature dependence of the IR spectra for CH modes, comparing the evolution for bare neutral and cationic coronene and [Fe-coronene] + .From such simulations, it is possible to derive the anharmonic band shifting and broadening.Figure 13 illustrates the evolution of the CH frequencies for the same systems as above.As another example, one can also follow the evolution of the linear anharmonicity coefficient in various PAHs as a function of size as shown in Figure 15

VBCI EXTENSIONS OF DFT AND DFTB
Most of the traditionally used functionals undergo troubles whenever the wavefunction is essentially multiconfigurational as it is the case for the dissociation of ionized molecular cluster.In particular, the  well known deficiencies of many functionals related to the self-interaction error lead to an unphysical repulsive behavior of the potential energy curve in the long range dissociation region (see Figure 16).Alternatively Valence Bond (VB) [58,59,[134][135][136][137][138][139] models do not meet such drawback and are well suited to deal with charge resonance.In such models, the wavefunction of the system is spanned over a reduced basis { A } of configurations, and the expansion coefficients are obtained by diagonalizing a hamiltonian matrix built from the charge-localized configurations energies and charge-transfer integrals.Each A corresponds to a configuration where the charge is localized on a given fragment A of the system.Such configurations do not present charge delocalization and can be obtained from constrained DFT calculations.This is the basis of the CDFT-CI method [140][141][142][143][144] which has been recently adapted for the DFTB framework (DFTB-VBCI see [145]).As can be seen from Figure 16, a correct behavior of the dissociation energy curve is recovered and the binding energy (corresponding to a so-called parallel displaced isomer) is in agreement with the experimental results.

BOND ORDER AND REACTIVE POTENTIALS
Bond-order potentials (BOP) were introduced to take into account the neighboring chemical environment of a given bond and describe several hybridization states in various geometrical configurations with different coordinations, still retaining a classical force-field formulation.Its derivation can be justified from a quantum framework and aims at incorporating in the force field nonadditive interactions with essentially quantum origin, showing in this respect connexion with the manybody embedded atom potentials [149,150] developed in the context of bulk materials.BOP were first 06003-p.17[145] approaches, compared with experimental binding energies/enthalpy taken from (a) (Ref [146]), (b) (Ref [147]) and (c) (Ref [148]).defined by Abell [151] under the form Only the attractive bond energy is described by a bond-order expression, close to the pair-embedding form, while pair-additivity is assumed for the repulsive short range potential between atoms i and j .The B ij are the many-body bond order functions, to be taken as functions of local coordination.Tersoff [152] popularized BOPs by parametrizing forms for bulk systems including C, Si and Ge.Pettifor [153][154][155] derived a systematic way to determine approximations for the potentials directly from the tight-binding formalism.The tight-binding cohesive energy is where ij , H ij are the elements of the single particle-density matrix, and tight-binding hamiltonian matrix respectively, i the atomic one-electron energies.Apart from the short range potential, the delocalization (or band) energy can thus be decomposed into a bond inter-site contribution and a on-site promotion contribution.In the following, we focus on the former.

06003-p.18
Chemistry in Astrophysical Media, AstrOHP 2010 The individual contribution e ij of the bond linking sites i and j can be expressed with the singleparticle density matrix in the atomic basis: with the bond order defined as i ,jν = 2 i ,jν .
The dependence of bond-order potentials on geometry is made explicit, through the parametrization of the bond orders i jν .In the case of s, p systems, it is necessary to rotate the atomic functions, so that the p functions are aligned along the bond, while the p x and p y lie perpendicular to it.The bond energy can then be separated as Potentials and bond orders must in principle be derived for both the and bonds in the pair, taking hybridization into account, see for instance Pettifor and coworkers [154][155][156][157]. Description of formation and breaking of chemical bonds is a further challenge.Reactive potentials are a recent extension of bond-order potentials that attempt to provide a general framework for simulating complex materials involving many atoms, long time scales, and chemical reactions.
The reactive empirical bond-order potential (REBO) developed by the Brenner group [158,159] is based on the Abell-Tersoff bond-order formalism [151,152], and is devoted to hydrocarbon chemistry.The interactions are written as for the repulsive and attractive components, respectively.The function f c (r) is 1 for the nearestneighbors, zero otherwise.The bond-order function B ij includes triplets (i, j , k), involving an angle dependence G(cos ij k ) used to detect the level of hybridization of each separate atom.While this approach allowed a good account of the diamond and graphite properties, small ring hydrocarbons and undercoordinated carbon atoms appear to be poorly described, and a more complex form was proposed [159].Brenner and coworkers also modeled the effects of conjugation by measuring the number N conj ij of neighbors connected to all pairs of bonded atoms.A conjugation energy associated with bond i-j is then defined to depend on the total coordination of atoms i and j , as well as the conjugation number N conj ij .Finally, the contribution of double bonds between carbon atoms is included, in order to reproduce the barrier for rotation around these bounds.The corresponding energy depends on the dihedral angles ij k and interpolations between integer values of the coordination and conjugation numbers [159].
In order to be used in molecular simulation, the coordination number N i is extrapolated to noninteger values using a sum of continuous distance-dependent switching functions, possibly restricted to specific types of atoms (C or H): The REBO potential is most efficient close to its training set, that is for crystals and hydrocarbon molecules.Its main drawback is that it lacks long-range interactions, which are important for surface or liquid state properties.Additional Lennard-Jones terms, with extra suitable switching functions have been proposed in the literature, such as the adaptive intermolecular reactive bond-order potential (AIREBO) of Stuart and coworkers [160], or the version derived by Los and Fasolino [161,162].

EPJ Web of Conferences
The most sophisticated approach is that developed in ReaxFF approach [163], where the energy of the hydrocarbonaceous system is partitioned into various contributions explicitly denoted as follows: This rather exhaustive decomposition takes into account : 1. the bond energy using a standard bond-order expression 2. an overcoordination penalty 3. an under-coordination penalty 4. a valence energy E val , which contains the main angular dependence related to (i, j , k) triplets.The effects of a possible lone electron pair are accounted as well. 5. an additional penalty energy E pen , designed to reproduced the stability of systems with two double bonds sharing an atom in a valency angle 6. a contribution E tors for torsion angles, which depends on the bond orders of the four atoms involved in the dihedral ij k , and includes a correction for overcoordination between two central sp 3 carbon atoms.7. a conjugation energy term E conj , which is maximum when successive bonds have bond-order values of 1.5 (similarly to the (AI)REBO potentials) 8. a long-range van der Waals interaction potential.9. an electrostatic energy contribution E Coulomb calculated using fluctuating charges obtained with a shielded Coulomb form.Even in its earliest version, the ReaxFF force field has a large number of parameters (more than 30), which explains its versatile character.Beyond hydrocarbons, ReaxFF has been extended and applied to pure and hydrogenated silica [164], metals and transition metals [165][166][167], metal hydrides [168,169], explosive materials [170][171][172], boron nitrides [173,174], and silicon fracture [175].Even though it has been used essentially for condensed materials, it incorporates a lot of molecular properties in its training sets.Reactive force fields along the lines of ReaxFF could thus become a method of choice also for nanoscale systems in the forthcoming years.
In the domain which is interesting to molecular astrophysics, one can mention the work of Bachellerie et al. [176][177][178], who modeled the recombination of gaseous hydrogen into molecular hydrogen at the surface of a model nanograins described by a graphene sheet using molecular dynamics.The bond-order model of Brenner was used, modified and reparametrized from DFT calculations to achieve agreement with DFT calculations following selected pathways for the H-H-coronene system.Bachellerie et al. were able to investigate molecular dynamics corresponding to the Eley Rideal or Langmuir -Hinshelwood candidate mechanisms for model systems with a graphene sheet of 200 carbon atoms, considering surface thermalization.Molecular dynamics showed the importance of considering surface relaxation in the simulation.Examples can also be found in other domains, for instance MD simulation of the formation and thermal cracking of hydrocarbons [179].

OUTLOOK AND FURTHER DEVELOPMENTS
In the present review, we have essentially considered the ability of Quantum Chemsitry methods to address the electronic ground state.This include both ab initio methods such as WF and DFT methods, as well as non ab initio methods such as tight-binding or non-additive force field methods.
As far as ab initio WF type methods are concerned, the present survey must be completed by the mention of Symmetry-Adapted Perturbation Theory [180] which is a scheme where the intermolecular energy is computed via Perturbation Theory defined from the fragment states, quite relevant to determine intermolecular interactions, the order of which is several orders of magnitude less than the total energies.In the other limit, one must also mention Quantum Monte Carlo approach [181][182][183], the aim of which

06003-p.20
Chemistry in Astrophysical Media, AstrOHP 2010 is to compute the total energy via stochastic sampling of the electronic wavefunction and integrals (properties) as mean values.The non ab-initio methods, much faster from the computer-time point of view, may offer proper tools (still within quantum framework for electrons) in intensive simulations of the dynamics of complex systems with numerous atomic degrees of freedom.Their reliability has been strongly enhanced by the use of parametrizations originating from ab initio calculations, usually pair matrix elements, which in DFTB for instance have been parametrized from LDA calculations.Parametrization should certainly benefit from fitting data obtained with more accurate functionals.Clearly, reactive potentials provide another option for the simulation of complex systems, within a classical framework but at the expense of loosing the explicit quantum description and manipulating analytical expressions with numerous parameters.Their determination also relies on accurate ab initio knowledge.Finally sophisticated fitting procedures are also being considered, such as on-the fly fitting in MD or MC simulations, dynamically rescaling the parameters with respect to intermediate ab initio level calculations along the simulations.This is also interesting in the context of multi-scale / multi-methods investigations.
Beyond IR spectroscopy, UV/visible spectroscopy assignment is certainly also of great interset and implies the theoretical determination of excited states.Let us just mention here a few schemes.As already mentionned, excited states can clearly be obtained in the WF approach, such as for instance in the CAS-PT2 scheme, where the CAS-CI excited roots are treated variationnally within the CAS and further perturbed up to second order Perturbation Theory.There is some active work in the framework of multi-reference Coupled Cluster Theory to address excited states [23][24][25][184][185][186], however such formulation is not yet quite standard.In the context of Coupled Cluster theories, an alternative approach is to use the equation of motion (EOM) and response theory following a CC determination of the ground state [187,188].As for the ground state, WF schemes are in principle exact but become rapidly extremely computer-time consuming, the difficulty being still increased by degeneracy and the intrinsic open-shell and multi-configurational character of most excited states.Some of the treatments exposed above for the ground state might be useful for excited states.Others may be difficult to transpose.For instance local orbital schemes can be used for valence type excitations, and even sometimes excitons, hardly for Rydberg states.
In the density formalism approach, the most widely used approach to excited states is certainly provided by the Time-dependent Density Functional Theory (TDDFT).The general TDDFT formalism was introduced by Runge and Gross [189] 1. an existence theorem stating that the action of en electronic system is a unique functional of the time-dependant density [189][190][191] density (r, t) 2. a variational theorem stating that the time-dependant density can be obtained by minimizing the action integral.Let us mention that the use of the standard leads to unphysical behavior.A solution was provided by van Leeurwen [198] who introduced the Keldish action.TDDFT allows to define the equation for the time-evolution of the density.In the Kohn-Sham formulation and incorporating the interaction with an external field w(r, t), the electron dynamics is given by i * (r, t) *t = [− 2 + v ne (r + w(r, t) The functional defining the action is non-local in both space and time, and so is its functional derivative.Often, the adiabatic approximation and the usual time-independant functionals are used without modifications, namely Recovering the steady states (excited eigenstates electronic potential energy surfaces) via the timedependant formalism is generally achieved within the response theory framework and in its usual

06003-p.21
EPJ Web of Conferences formulation via a linearization (small perturbing fields) and a Fourier transform in the frequency space.Steady states TDDFT calculations are now routine in most Quantum Chemistry packages, and have been for instance successfully used to help to analyze the electronic spectra of PAH molecules [195][196][197].Again, while TDDFT is much more tractable than WF function schemes, it suffers from the ignorance of a correct functional (which should be non-local both in time and in space).The defects already present in the time-independant formulation, in particular the SIC error, naturally affects TDDFT.The SIC error for instance induces an incorrect potential for excited electrons, departing from the expected long-range Coulomb behaviour, which must be cured.Charge transfer excited states also are described quite poorly.Nevertheless, TDDFT methods offer a promising tool to molecular response properties and also to non adiabatic time-evolution of electrons in excited states, involving highly non-linear physics in intense fields.Conceptual, methodological and algorithmic works dedicated to the improvement of the functionals, numerical implementation and the coupling with the nuclei dynamics, are in progress in several groups.

Figure 3 .
Figure 3. Key elements for a CI calculation.

Figure 4 .
Figure 4. HOMO and LUMO orbitals for small beryllium nanochains deposited on a graphene surface.Figure taken from Ref. [47].

Figure 6 .
Figure 6.H 2 molecule : configurations describing non-dynamical (two electrons in the s orbitals) vs dynamical correlation (two electrons out of the valence shell, here in the p orbitals).

2p a 2p b
| with covalent ones | 1s a 1s b |.This fluctuation process involves simultaneous local single excitations of two electrons from their 1s grounds states to 2p
The two electrons of H − initially the 1s-like orbital undergo fluctuation towards configurations of the 2p shell like | 2p b 2p b | which determine their angular correlation, while excitations towards | 2s b 2s b |

Figure 11 .
Figure 11.Experimental and calculated values (in eV) of the ionization potential differences (I P ) = I P (P AH ) − I P (Si) with various methods (shift : with parameter corrections) from Ref. [56].

Figure 12 .
Figure 12.Potential energy surface along a selected pathway for Fe at the surface of a coronene molecule.Comparison between DFT and DFTB (taken from Ref.[133]).

Figure 13 .
Figure 13.Frequency dependence of the CH band as a function of temperature (from MD/DFTB calculations) (taken from Ref. [133]).

Figure 15 .
Figure 15.Evolution of anharmonicity factors characterizing the CH and ν CC modes as a function of PAH size (from MD/DFTB calculations; full line : PAH; dotted line : PAH + ).

EPJ
Web of Conferences