Francesco Aquilante

Jochen Autschbach

Rebecca K Carlson

Liviu F Chibotaru

Mickaël G Delcey

Luca De Vico

Ignacio Fdez

Galván

Nicolas Ferré

Luis Manuel Frutos

Laura Gagliardi

Marco Garavelli

Angelo Giussani

Chad E Hoyer

Giovanni Li Manni

Hans Lischka

Dongxia Ma

Per-Åke Malmqvist

Thomas Müller

Artur Nenov

Massimo Olivucci

Thomas Bondo Pedersen

Daoling Peng

Felix Plasser

Ben Pritchard

Markus Reiher

Ivan Rivalta

Igor Schapiro

Javier Segarra-Martí

Michael Stenrup

Donald G Truhlar

Liviu Ungur

Alessio Valentini

Steven Vancoillie

Valera Veryazov
email: valera.veryazov@teokem.lu.se

Victor P Vysotskiy

Oliver Weingart

Felipe Zapata

Roland Lindh
email: roland.lindh@kemi.uu.se

Molcas 8: New Capabilities for Multiconfigurational Quantum Chemical Calculations across the Periodic Table

In this report we summarize and describe the recent unique updates and additions to the Molcas quantum chemistry program suite as contained in release version 8. These updates include natural and spin orbitals for studies of magnetic properties, local and linear scaling methods for the Douglas-Kroll-Hess transformation, the generalized active space concept in MCSCF methods, a combination of multiconfigurational wave functions with density functional theory in the MC-PDFT method, additional methods for computation of magnetic properties, methods for diabatization, analytical gradients of state average complete active space SCF in association with density fitting, methods for constrained fragment optimization, large-scale parallel multireference configuration interaction including analytic gradients via the interface to the Columbus package, and approximations of the CASPT2 method to be used for computations of large systems. In addition, the report includes the description of a computational machinery for nonlinear optical spectroscopy through an interface to the QM/MM package Cobramm. Further, a module to run molecular dynamics simulations is added and two surface hopping algorithms are included to enable nonadiabatic calculations. Finally, we report on the subject of improvements with respects to alternative file options and parallelization.

I. INTRODUCTION

The Molcas quantum chemistry package was for a long time a programming project confined to Lund University, the group of the late Björn O. Roos, and a few external coworkers (see the Introduction of Ref. 1). During the last seven to eight years, however, this situation has drastically changed, and the project has turned into an international collaboration. Our recent report [START_REF] Aquilante | [END_REF] in 2009 included contributions from nine different research institutes, and we are pleased to note that the current report exceeds this with contributions from some twenty-two different affiliations and thirty-eight contributors. During the last three years the program developers of this project have gathered for annual meetings -Zürich, Switzerland (2013), Alcalá de Henares, Spain (2014), and Siena, Italy (2015)gathering some 40 developers at the last meeting. This reflects the changing nature of program and method development in the field of quantum chemistry. While in the past a small group of developers could harness most of the aspects needed to develop and maintain a quantum chemistry package, this is no longer possible. Today, program development for state-of-the-art quantum chemistry applications, in the context of small groups with limited funding for this task, has to be carried out in collaboration between international partners and with shared responsibility for maintenance, code verification and facilitating a common optimal program development environment. Typically national funding for development of quantum chemical method and program development is too small and short term to support any but the simplest quantum chemistry program packages. While this could be viewed as a criticism of the current granting and financing system of such projects, it has also been an enormous incitement for honest and sometimes even altruistic collaborations in the field. It has fostered a unique understanding among young quantum chemists of the added value of networking, cooperation, and partnership. In this respect the last six years of developments of the Molcas quantum chemistry program package have been spectacular.

In particular, the Molcas development environment has seen a number of long-needed updates and modernizations to support the new form of extensive international collaboration around a single program package. This has enabled and facilitated fast development and a number of new and unique implementations have seen the light of day. With the release of Molcas 8, in 2014, we yet again provided the chemistry community with novel and interesting tools for state-of-the-art computational chemistry. To offer potential users a chance to see all these improvements in a single read, we have gathered a compact summary of these unique features (some already in version 8.0, others to be available in the upcoming release of 8.2) here.

This communication includes reports with respect to alternatives, extensions, and improvements for the treatment of dynamic and nondynamic electron correlation (nondynamic correlation is sometimes called "near-degeneracy correlation", "strong correlation", or "static correlation"), facilities for computation of magnetic properties, methods for scalar rela-tivistic corrections to large molecular systems, and approaches for molecular dynamics for semiclassical or full quantum nuclear dynamics of nonadiabatic processes. Some technical improvements will be given an explicit presentation.

II. ELECTRON CORRELATION METHODS

The multiconfigurational approach is one of the hallmarks of the Molcas quantum chemistry program package. Central to this treatment are the CASSCF and CASPT2 methods. The essential building blocks are the separate treatments of the dynamic and nondynamic electron correlation. In this section we present (a) a brief overview of the scaling properties of a complete and restricted active space, (b) an implementation of a generalization of the CASSCF methods, the GASSCF method, (c) an alternative method for the introduction of dynamic electron correlation already in the CASSCF model, (d) approximations to the CASPT2 method to enable applications to large systems, (e) an implementation of an efficient parallelization of the CASPT2 method, and (f) a job-farm approach for parallel separate MS-RASPT2 states.

A. Scaling properties of the active space

The total number of Slater determinants N SD that can be generated by considering all possible distributions of m electrons in n orbitals with a spin projection M s (which includes the determinants for all S ≥ M s with M s half-integer or integer depending upon odd or even electron cases) is given by the following formula:

N SD (m, n, M s ) = n m 2 + M s n m 2 -M s (1)
When expressed with α and β as the number of α and β electrons respectively, this equation leads to:

N SD (n, α, β) = n α n β ( 2
)
as m = α + β and M s = α-β 2 . The notation N SD (n, α, β) is implied during the remainder of this section. To get the number of configurations with spin S = M s , we only need to eliminate the amount of determinants that originate from higher spins:

N CSF (m, n, S) = N SD (m, n, M s ) -N SD (m, n, M s + 1) = N SD (n, α, β) -N SD (n, α + 1, β -1) = n α n β - n α + 1 n β -1 (3)
Note that we assume here that for the binomial coefficients, n k = 0 when k < 0 or k > n. We've plotted the number of determinants and configurations for a series of increasing sizes of an active space CAS(m,n) (m electrons in n orbitals) for m = n and M s = 0 in Figure 1. The grey area represents the zone where the calculations require significant computational resources. Beyond the grey area, calculations become unfeasible. As can be seen, the number of determinants/configurations grows exponentially with the number of orbitals in the active space. If we rework part of the expression for binomial coefficient, we can set an upper bound:

n k = n! k!(n -k)! = (n -k + 1) • • • n k! ≤ n k k! ( 4
)
For M s = 0, we have α = β, and the number of determinants is thus bounded by:

N SD (n, α, α) = n α n α ≤ n 2α α! 2 (5)
For a constant number of electrons, the number of determinants is bounded by a polynomial. When the number of electrons grows with the number of orbitals, e.g. m = n = 2α as in our example, an upper bound to the number of determinants can be set by using

N SD (n, α, α) ≤ n 2α α! 2 ≈ n 2α 2πα α e 2α = n n πn n 2e n = (2e) n πn ( 6
)
where we used Stirling's approximation for the factorial (n! ≈ √ 2πn( n e ) n ). If instead of an upper bound we use Stirling's approximation exclusively, the number of determinants can be estimated as 2 πn 4 n , and similarly the number of configurations as (1 -( n n+2 ) 2 ) 2 πn 4 n which can even be further approximated as 8 πn 2 4 n for large n, or even 2 n 2 4 n to quickly get the order of magnitude. The approximations are represented by solid lines in Figure 1.

It can be easily shown that the exponential scaling applies to any general situation by considering the total number of possible determinants for any number of electrons in n orbitals. The number of ways one can fill 2n spin orbitals with electrons is equal to 2 2n = 4 n , i.e. each spin orbital can be either occupied or unoccupied, independent of the other spin orbitals.

Sometimes, the scaling of the number of determinants/configurations with the number of active orbitals is called "factorial" scaling (referring to the binomial coefficients), but we do not recommend this term as this would imply an n n growth (i.e. growing faster than exponential).

( n n+2 ) 2 ) πn 4 n ).
For a restricted active space, the total number of determinants will be reduced because of the extra conditions imposed on the distribution of the electrons. Considering again an active space of m electrons in n orbitals, now with n 1 , n 2 , n 3 orbitals in each of the RAS spaces respectively, and a maximum of h holes in RAS1 and p electrons in RAS3. In this case the total number of determinants is given by:

N SD (n, α, β) = n 1 α 1 ,β 1 =n 1 -h α 1 +β 1 ≥2n 1 -h p α 3 ,β 3 =0 α 3 +β 3 ≤p n 1 α 1 n 1 β 1 n 3 α 3 n 3 β 3 n 2 α -α 1 -α 3 n 2 β -β 1 -β 3 (7)
The number of configurations can be obtained in a similar way as for the complete active space, i.e. N CSF = N SD (n, α, β) -N SD (n, α + 1, β -1).

In order to give an idea of the scaling of the number of determinants for certain RAS setups, we plotted the scaling for two reference setups. The first is for a total active space of n electrons in n orbitals, with a fixed number of RAS2 orbitals, and symmetric RAS1 and RAS3 spaces (same number of orbitals and matching holes/excitations) (Figure 2).

With an empty RAS2 space (the yellow circles), we can comfortably use up to 40 active arbitals as long as only up to triple excitations (sdt) are considered. For up to quadruple excitations (sdtq), the active space limit is reached somewhere between 28 and 36 orbitals (without symmetry). For each increase of the inner RAS2 space with 4 orbitals, the number of determinants goes up by 1 to 2 orders of magnitude. With an inner RAS2 space of 12 electrons in 12 orbitals, RAS1/RAS3 excitation levels beyond two electrons do not stretch the limit on the total number of active orbitals any longer.

The second figure is for a fixed RAS2 space of n electrons in n orbitals, with a growing number of empty orbitals in RAS3 (Figure 3). In general, an increase of 4 orbitals in RAS2 causes an increase in the number of determinants of 2 to 3 orders of magnitude. We can see that even with a RAS2 space of 12 electrons in 12 orbitals, the total active space can be stretched up to about 40 orbitals as long as one limits the excitation level to single or double (s, sd). When allowing higher excitation levels (sdt, sdtq), a smaller RAS2 space is needed to keep the calculations feasible for such a large number of active orbitals.

So far, any numbers shown are for non-symmetric systems. It is important to note that any additional symmetry constraints will reduce the number of determinants/configurations too, enabling some calculations that would be unfeasible without symmetry. The computational resources required by the Molcas programs rasscf and caspt2 are discussed in the FIG. 2: Number of determinants for a restricted active space of n electrons in n orbitals (M s = 0), for various sizes of the inner RAS2 and different excitation levels out of/into RAS1/RAS3. The solid line represent the upper bound of the number of determinants.

Supporting Information, in sections S1 and S2 respectively.

B. The Generalized Active Space Self-Consistent Field Method

Strong or nondynamic electron correlation is a key component of any electronic structure calculation that aims for quantitative accuracy in calculated energies. With Molcas it is possible to calculate strong electron correlation by using multiconfigurational self-consistent field methods, MCSCF. Earlier versions of Molcas allowed two possible kinds of MCSCF calculations: the complete-active-space self-consistent-field (CASSCF) [2][3][4] method and the restricted-active-space self-consistent-field (RASSCF) method [5][6][7] . The present version of Molcas allows a more flexible choice, namely the generalized-active-space self-consistentfield (GASSCF) method. 8 This method, like RASSCF, allows restrictions on the active space, but it is more flexible than RASSCF, and it can be applied to larger and more complex systems at affordable cost. If the active space is well chosen and the restrictions are not too severe, MCSCF methods recover most of the nondynamic correlation energy, and part of the dynamic correlation energy. MCSCF methods are generally used when near-degeneracy effects are present, e.g. in many transition metal complexes, excited states, bond breaking and weakening, many structures along reaction paths, and many radicals. In CASSCF the orbital space is divided (by the user) into an inactive space consisting of orbitals that are doubly occupied in all considered configuration state functions, CSFs, an active space consisting of orbitals whose occupation may take any value between 0 and 2, and a virtual space consisting of orbitals that are unoccupied in all considered CSFs. In CASSCF, a number of active electrons are distributed in all possible ways, compatible with the required spin and space symmetry, among the active orbitals, to generate the CI space in which the wave function will be optimized. Both orbital and CI parameters are simultaneously optimized.

The size of the CI expansion increases factorially with the size of the active space. The largest active space currently affordable is formed by 18 electrons in 18 orbital for a singlet spin state. Only a fraction of correlation energy is covered for most of chemical systems of practical interest, while most of the dynamic correlation, which is essential for a quantitative treatment of chemical properties, is not included. Restrictions on the CI space can be applied through the restricted active space, RASSCF, and the generalized active space, GASSCF, methods available in Molcas. The underlying philosophy of these approaches is to partition the active space in subspaces and impose some user-defined restrictions on the electron excitations in such a way that one removes the configurations that contribute only marginally to the total wave function. Without restrictions, both RASSCF and GASSCF methods become equivalent to CASSCF. These restrictions are particularly useful when the cost of using the full CI expansion of the active space is beyond reach. It is noteworthy that truncations of the complete active space, while reducing the size of the CI problem, may result in difficult optimization computations. In RASSCF, the active space is divided into three subspaces, namely the RAS1, RAS2 and RAS3 spaces. The RAS1 orbitals are doubly occupied in the reference wave function. The RAS2 orbitals may have occupation numbers varying from 0 to 2, while the RAS3 orbitals are empty in the reference wave function. The CI expansion is then generated by all the possible excitations in RAS2 plus all those configurations with up to a user-specified number of holes in RAS1 and a user-specified number of particles in RAS3. The GASSCF method is a generalization of the RAS concept. In the GASSCF method, instead of three active spaces, an arbitrary number of active spaces (GAS1, GAS2, . . . ) may in principle be chosen by the user. Instead of a maximum number of holes in RAS1 and particles in RAS3, accumulated minimum and maximum numbers of electrons are specified for GAS1, GAS1+GAS2, GAS1+GAS2+GAS3, etc. in order to define the desired CI expansion. The GAS scheme can be seen as the generalization of the RAS method, and indeed it can reduce to CAS and RAS, respectively when one and three active spaces, respectively, are chosen. The GAS approach has points in common with the occupation-restricted-multiple-active-space (ORMAS) method. 9,10 To illustrate the motivation for and advantages of the GASSCF method, we consider three examples. In the first example, we have an organometallic material or Metal-Organic-Framework (MOF) with separated metal centers such that the orbitals are not delocalized across the metal centers.

Then one can include the near-degenerate orbitals of each center in its own GAS space. This implies that one may choose as many GAS spaces as the number of multiconfigurational centers. A second example is for lanthanide or actinide metal compounds where the f -electrons are near the HOMO-LUMO region but do not participate in bonding. In this case one can put the f-electrons into one or more separated GAS spaces and not allow excitations from and/or to other GAS spaces (see Ref. 8 for details). Thus the wave function would be largely simplified without losing accuracy. A third example would be organic molecules or other molecules where there are localized bonds. Then each bond and its correlating orbital could form a separate GAS space as in GVB. 11 Symmetry considerations may also be invoked for the partitioning of GAS spaces. There is no rigorous scheme to choose a GAS partitioning.

The right GAS strategy is system-specific. This makes the method versatile but at the same time it is not a "black box" method. CASSCF, RASSCF and GASSCF use the Davidson method as CI eigensolver and the direct-CI algorithm by Olsen. 12 New approaches, where very large CI expansions are represented without explicit representation of CI vectors, could be used in MCSCF calculations to replace the deterministic CI eigensolver approaches and are envisioned for the next generations of Molcas. [13][14][15][16] Perturbation Theory Approaches. Due to the limited size of affordable active spaces, most of the dynamic correlation is not recovered by MCSCF methods. The CASPT2 method (perturbation theory through second order on top of a CASSCF wave function) 17,18 and its extension to the Multi-State variant, MS-CASPT2 19 , has been available in Molcas since the nineties. RASPT2 is defined similarly starting with RASSCF wave functions. 7,20 The GASPT2 variant has been implemented locally 21 and it is envisioned to be released in a next version of Molcas. For these perturbative approaches, the preceding MCSCF method provides a well-behaved wave function and should recover most of the strong correlation and part of the dynamic correlation, while PT2 should yield good approximations to the dynamic correlation energy. The CASSCF/CASPT2 method has reached a mature state and has proven successful for atoms and molecules in their ground states as well as excited states.

It is worth mentioning that other post-SCF methods are also available to recover the missing correlation, such as NEVPT2

C. Multi-Configuration Pair-Density Functional Theory

In the present version of Molcas, another method has been recently added for combining MCSCF with a density functional calculation for dynamic correlation energy, the multi-configuration pair-density functional theory (MC-PDFT). 24,25 Dynamic correlation effects not described at the MCSCF level can be recovered by MC-PDFT in a way that is much less expensive than CASPT2. In Kohn-Sham density functional theory, KS-DFT 26 , the electronic energy is expressed as a functional of the spin densities, their gradients, and possibly some orbital-dependent quantities such as Hartree-Fock exchange calculated from occupied orbitals or correlation energy calculated from virtual orbitals. KS-DFT can be used for systems that are too large to be treated with CASSCF/CASPT2, since it uses a single-determinantal description of the electron density. However, with currently available density functionals it does not always give a reliable answer when a single-configuration approximation does not provide a good description of the wave function. In such cases it is often necessary to use ad hoc weighted averages of broken-symmetry wave functions to obtain reasonable energies. [27][28][29][30] As a consequence of using approximate functionals in a spin-

ρ(r) = N |Ψ(x 1 , x 2 , x 3 . . . , x N )| 2 dσ 1 dx 2 . . . dx N r 1 =r (8) Π(r) = N (N -1) 2 |Ψ(x 1 , x 2 , x 3 . . . , x N )| 2 dσ 1 dσ 2 dx 3 . . . dx N r 1 =r 2 =r (9)
The present version of Molcas supports MC-PDFT based on three kinds of MCSCF wave functions: CASSCF, RASSCF and GASSCF. In the MC-PDFT method the energy is calculated as

E = V nn + Ψ| T + Vne |Ψ + V C [ρ] + E ot [ρ, Π] ( 10
)
where |Ψ is the multiconfigurational wave function, V nn is the nuclear repulsion term, T is the kinetic energy operator, Vne is the electron-nuclear interaction operator, and V C [ρ] and E ot [ρ, Π] are respectively the electronic Coulomb energy and the on-top density functional of the total density and the on-top pair density. This energy expression may be written in terms of the one-electron density matrix, D, and the on-top pair density, Π, as

E = V nn + pq h pq D pq + 1 2 pqrs g pqrs D pq D rs + E ot [ρ, Π] (11)
The wave function is thus used to compute the electron-nuclear interaction, the kinetic and their derivatives on a grid and computing the contribution to the energy by numerical integration.

In the current implementation, MC-PDFT can be used also in combination with stateaverage multiconfigurational wave functions, thus allowing one to compute excited states with the same spatial and spin symmetry as the ground state (excited states with a different symmetry than the ground state can be computed by performing separate CASSCF calculations for each symmetry). The method is also very useful for computing ground-state potential energy curves.

D. Multiscale perturbation theory for computational photophysics and photochemistry

Multiconfigurational perturbation theory in the form of CASSCF/CASPT2 is the workhorse of the Molcas package. The method is particularly appreciated for its robustness in tackling problems featuring strong correlations for ground and excited states. It has therefore become a de facto gold standard for computational studies in areas such as photochemistry and photophysics where other electronic structure methods fail to provide comparable uniform accuracy or cannot compete in terms of computational costs.

In fact, thanks to the development of the Cholesky decomposition (CD) approximation, algorithms for CASSCF/CASPT2 calculations allow at present the study of systems comparable in size to those affordable at the DFT level of theory -whenever the active space required is within the current practical limit of 18 electrons in 18 orbitals.

In such cases, the bottleneck in the use of the CD-based CASSCF/CASPT2 method with large basis sets and many electrons resides in the CD-CASPT2 step. Two multi-scale extensions have been explored in order to lower the computational costs associated with the CASPT2 treatment of the dynamic correlation: FNO-CASPT2 31 and FD-CASPT2. 32

The first of these approaches employs the so-called frozen natural orbital (FNO) approximation. In FNO-CASPT2, one truncates the number of virtual orbitals by retaining only those that account for most of the dynamic electron correlation. Such truncation of the virtual space is not possible in the canonical representation, as the resulting loss of accuracy is too severe, but can be effective when the orbital representation is chosen based on different criteria. In particular, one can define a transformation matrix from the (pseudo)canonical orbitals to a set of approximate natural orbitals (NOs) by diagonalizing the following virtualvirtual block of the following MP2-type density matrix:

D(2) ab = ic t ac ii t cb ii , ( 12
)
with the amplitudes defined as:

t ab ii = - (ai|bi) a + b -2 i , ( 13
)
where the index i includes all the inactive orbitals as well as those (pseudo)canonical active orbitals with negative energy. By their nature, even the NOs defined through the diagonalization of an approximate density matrix such as that of Eq. ( 12) allow a hierarchical build-up of dynamic electron correlation, and therefore elimination from the virtual space of NOs with small occupation numbers is possible in this orbital basis. The two-electron integrals included in Eq. ( 13) can be conveniently computed by means of their CD/DF representation, 33 without introducing any fifth-order step that would otherwise make the truncation of the virtual space too costly compared to the full CD-CASPT2 calculation.

The selection of the relevant NOs is performed by means of the following estimator of the correlation energy:

%Tr = 100 ν a η a Tr( D(2) ) , ( 14
)
where the number of retained virtual NOs ν is determined by including more and more virtual NOs -preordered according to the decreasing occupation numbers η's (eigenvalues of D( 2) ) -until a chosen value of %Tr is reached. Typically, even with double-zeta basis sets more than 98% of the trace is determined by only half of these NOs, thus allowing for very significant speed-ups in the calculation. In fact, experience tells that a threshold for %Tr in the range 90-95% is in general sufficiently accurate for vertical excitations whereas 97-99%

provides results that are nearly indistinguishable from those obtained with full CD-CASPT2 even along ground and excited states reaction pathways. An example of performance of FNO-CASPT2 is shown in Fig. 4.

In contrast to the use of NOs, in freeze-and-delete (FD)-CASPT2 the computational costs are lowered by resorting to a localized orbital picture in CD-based CASPT2 calculations. For situations where the active orbitals are localized within a small region of a large molecule, a suitable "active site" (A) can be identified as the collection of atoms where the active orbitals effectively extend. The selection of an atom is based on the measure of the Mulliken charge contribution of the active orbitals to that atom, as compared to a chosen threshold (τ ).

Accordingly, the inactive and secondary orbitals can be separately localized and partitioned between this active site and the remaining atoms of the "environment" (B). The two regions are assumed to be uncoupled, and therefore two separate sets of canonical orbitals can be deduced for the active site and the environment.

Assuming such partitioning of the molecule in active site and environment, we can ap- Dashed lines represent the water donor localised electronic excited state energies whereas full lines refer to the ground state. The geometries of the water dimer at the Franck-Condon region (left-hand side) and at the dissociation limit (right-hand side) have also been given to provide a visual aid for the photo-process.

proximate the first-order correction to the CASSCF wave function Ψ (0) as:

Ψ (1) ≈ [A] pqrs tpr qs â † p â † r âq âs Ψ (0) + [B] aibj tab ij â † a â † b âi âj Ψ (0) . ( 15
)
The hypothesis of decoupled regions allows furthermore to determine the two sets of amplitudes separately. In practice, this means that the standard algorithm for CASPT2 can be used to obtain the (small) set of amplitudes defined in A, whereas the much larger number of amplitudes defined in B correspond to a straightforward MP2 model for which their expression in terms of integrals and orbital energies is known and it is computationally easier than for CASPT2. As shown in Table I, significant speed-ups are therefore possible for large molecules, and with a minor penalty for the accuracy.

E. Parallelization of multiconfigurational methods in Molcas

The rasscf and caspt2 programs have always been among the key components of the implementations, they needed to be adapted to be able to take advantage of a distributed environment. Another indirect consequence is that such a parallel rasscf/caspt2 implementation could also run faster even on a single shared-memory machine by using multiple processors/cores. In this contribution, we focus mainly on the parallelization of the caspt2

program 35 , and at the end include some information about both caspt2 improvements and rasscf parallelization which has already been implemented and will be available soon.

We should first start by saying that multiconfigurational perturbation theory isn't uniquely defined, and several different methods exist, such as complete or restricted active space second-order perturbation theory (CASPT2/RASPT2) 7,[17][18][19][20]36,37 , quasi-degenerate multireference perturbation theory (QDMRPT) [38][39][40] , n-electron valence perturbation theory (NEVPT) [41][42][43] , and occupation-restricted multiple active space perturbation theory (ORMAS-PT) 44 . While methods such as NEVPT have been parallelized 22 , to the best of our knowledge, this is the first fully parallel implementation of the CASPT2 method.

In CASPT2/RASPT2, one solves a set of equations of the Rayleigh-Schrödinger type as a large matrix equation, which eventually reduces to the following expression:

(Λ D + FN -E 0 1) c = -R (16)
where

Ω c = c, R = Ω † R, and Ω = U Λ -1 2 S W such that Λ S = U † SU and Λ D = W † FD W with FD = (U Λ -1 2 S ) † F D U Λ -1 2
S . The present implementation of the CASPT2 method consists of four phases: (i) computing and diagonalizing the S and F D matrices; (ii) constructing the right-hand side (RHS) vector R; (iii) solving the equation system that gives the amplitudes c describing the perturbation; and finally (iv) using these to compute energies and other properties of interest 1 .

Construction and diagonalization of S and F D

In the first phase, the overlap matrices (S) and the diagonal blocks (F D ) of the zerothorder Hamiltonian matrix, that make up the left hand side of the equation, are constructed.

Hereto, several density matrices need to be constructed, of which the largest is the 3-body density matrix Γ. The size of Γ scales roughly with the sixth power of the number of active orbitals, while the amount of floating point operations to compute a single matrix element is proportional to the number of configurations (2 • N CSF ). The total work required to set up Γ thus scales as n 6 a • N CSF . The elements are computed using task-based parallelism, where each task computes a set of elements. After the S and F D matrices have been constructed, they need to be diagonalized, for which the amount of floating point operations scales as n 9 a . The diagonalization is carried out through a call using the scalable linear algebra package (ScaLAPACK) interface, which can then be substituted by any appropriate parallelized implementation of it.

The dominating factor for large complete active spaces (e.g. 16 electrons in 16 orbitals) will be the construction of the density matrix, in particular the amount of work needed per element, while the size of the total matrix remains relatively small. In this case, the diagonalization will not be a significant bottleneck. On the other hand, for a restricted active space reference, the number of active orbitals could easily double (because the number of configurations can be kept small). This could lead to very large density matrices that need to be stored across different machines. In this case, diagonalization will also be a reasonably significant bottleneck.

Construction of the right-hand side

In the second phase the right-hand-side (RHS) vector R is constructed. The latter vector consists of different blocks, each related to a specific excitation case. The size of this vector is dominated by an inactive-virtual block, that is the part of the first-order wave function which is generated by excitations ( XP = Êajcl ) from the inactive orbitals (j, l) to the virtual orbitals (a, c). The size of this block thus scales as n 2 i • n 2 v , the product of the number of inactive (n i ) and virtual (n v ) orbitals squared. For a long time, the RHS by itself was too large to be kept in memory and this resulted in substantial I/O bottleneck. However, the ever increasing size of memory available on single machines has lifted this bottleneck and the RHS can be typically kept fully in memory for medium-sized systems. When going to very large basis sets, we can use the combined resources of multiple machines to distribute the RHS.

For the distribution of the work during construction of the RHS, the parallel implementation of the Cholesky decomposition technique is used 33,45 . The Cholesky vectors are distributed over the available computational units, such that each unit stores a range of columns which are transformed to MO basis at the beginning of the program. Two algorithms are available to compute the RHS. The first algorithm computes for each element of the RHS, a contribution from the local set of Cholesky vectors. The downside of this method is that for a large number of computational units, communicating the entire RHS becomes a communication bottleneck. The second algorithm first collects all the Cholesky vectors on each computational unit and then computes only its local part of the RHS.

Iterative solution of the CASPT2 equation

In the third phase, the CASPT2 equation is solved by means of the preconditioned conjugate gradient (PCG) method 46,47 . This method iteratively solves a system of equations Ax = b, where A is a sparse matrix, which in this case corresponds to the zeroth-order Hamiltonian, x is the solution vector, which corresponds to the first-order wave function, and b corresponds to the RHS. An iterative method is needed because the matrix A is too large to be constructed explicitly and hence cannot be inverted. The basic operations during the PCG iterations consist of one matrix-vector multiplication, vector inner products, scalarvector multiplication, and vector addition. The PCG routine in Molcas stores 6 vectors, for use during the intermediate computations. During the sparse matrix vector multiplication, the different blocks of the RHS vector will need to interact with each other, causing elements from one process to be communicated to another for at least one block. The algorithm is written to avoid communicating any elements from the largest part of the solution vector, i.e. the inactive-virtual block.

Overall performance

For large CAS references, execution time is dominated by the density matrix construction and scaling is good although dependent on the available memory bandwidth. This is clear in Fig. 5b, where the scaling is better for a distributed parallel environment. After the initial speed-up, execution time levels off as the serial part becomes dominant, which shows a serious drop in parallel performance from 16 to 32 processes. The calculation with large RAS reference is completely dominated by memory bandwidth, showing very bad scaling on a single shared-memory machine (Fig. 5c). Execution time is dominated by the diagonalization step, which scales well up to 32 processes in a distributed environment. Finally, for a system with a very small reference but large basis set, total execution time is essentially defined by the RHS step, which is again dependent on the available memory bandwidth and scales badly on a single machine (Fig. 5d). When going to multiple processes, communication overhead becomes too large, and the total execution time will start to increase beyond a certain number of processors. This is the reason to switch to a new method where the Cholesky vectors are collected first, after which no more communication would be needed. This alternative has been implemented but is still in a testing phase and will be available soon.

To conclude this subsection, saturation of memory bandwidth and communication overhead lead in practice to a parallel performance that typically scales well with the number of physical CPUs but not with the number of available cores. Furthermore, I/O bottlenecks on a single machine can be only be further reduced by running on multiple machines, or by providing enough memory and alleviating the I/O bottleneck. Even with these limitations, significant time savings for large calculation can be achieved by increasing the number of processes on a single machine, as long as memory bandwidth allows. Calculations that took more than 3 days on a serial machine, could be performed in less than 5 hours on an InfiniBand cluster, where the individual nodes were not even capable of running the calculation because of memory and I/O requirements. This ensures the continuing study of larger molecular systems by means of CASPT2/RASPT2 through the use of the aggregated computational resources offered by distributed computing systems. Last but not least, the rasscf program which delivers the reference wave functions for the caspt2 program, had been partly parallelized before by making use of the distribution of Cholesky vectors. However, the part that handled the configuration interaction was still serial. Now, this has also been parallelized in the current development version and will appear in a next iteration of Molcas-8. This parallelization is not intended to use aggregate resources but to speed up the configuration interaction part of the method. This greatly benefits studies of systems with very large active spaces, and is the last piece in the parallel execution of the seward/rasscf/caspt2 chain for Cholesky-based calculations.

As can be seen from Figure 5a, the parallel efficiency follows a similar pattern as for the matching CASPT2 calculation (Figure 5b), although performance is slighlty worse. The walltime for the CASSCF part goes from 21% of the total 17h CASSCF+CASPT2 walltime in serial to 27% of the total 1.5h walltime when run in parallel with 32 processes.

F. Parallel computation of separate MS-RASPT2 states using a job-farm

Even considering all the improvements given by the usage of Cholesky decomposition and parallelization of the RASPT2 code, the calculation of RASPT2 energies for large molecules remains a task that may require many hours of CPU time. This is even more true when dealing with MS-RASPT2.

The accurate description of the spectra of many molecules was at the base of the success of the MS-RASPT2 methodology 48,49 . As previously described, a reference multiconfigurational wave function is prepared, as a state average between all states of interest, and subsequent RASPT2 calculations are performed for each root. However, there are cases when a RASPT2 treatment alone is not sufficient. For example, when two (or more) states are energetically very close to each other, as in the case of a conical intersection or an avoided crossing, it is indispensable to employ a multistate treatment.

As multistate RASPT2 (MS-RASPT2) is currently implemented in Molcas, the information about the Hamiltonian effective coupling terms is prepared at each RASPT2 step.

In the final multistate step, this information is collated and symmetrized, in order to produce the MS-RASPT2 results. For this reason, a typical MS-RASPT2 calculation has to be performed all in one go. As an example, if each RASPT2 calculation for a given large system requires one day of CPU time, computing a three roots MS-RASPT2 would require at least 3 days, and the larger/more complex the molecule under scrutiny, the longer would be the necessary time. By increasing the time required to complete a task, the probability that something goes awry with the machine performing the actual calculation also increases.

Unfortunately, in case of computer malfunction, a MS-RASPT2 calculation would stop, and would have to be submitted again from the beginning, since it cannot be restarted from any

given point (at least, not as yet).

However, the implementation of two sets of keywords for the caspt2 module allows to divide a MS-RASPT2 calculation into its RASPT2 components, plus the respective coupling terms. By this simple expedient, an effective coarse parallelization is achieved: following the previous example, the three roots can be submitted in parallel as separate jobs. The coupling terms will then be collected and used to quickly calculate the MS-RASPT2 energies. The so obtained coarse parallelization will take the same amount of time it takes for a single state, as long as the computational resources increase linearly with the requested roots. Each calculation can also still take advantage from the usage of Cholesky decomposition and fine grain parallelization for further speed up, as previously described. The detailed usage of this new possibility is given in a working example in section S5 of the Supporting Information.

The job-farm intended use is for when a multistate treatement is necessary, and each RASPT2 calculation is expected to take a very long time (e.g. even with single-root finegrained parallellization enabled), shortening the overall calculation by days. Possibly, the speed up obtained by splitting a MS-RASPT2 calculation, as described here, together with the fine grain parallelization previously presented, will prompt the usage of this methodology also for larger systems, which were until now considered untreatable with a multistate, multireference perturbation theory approach, because too lengthy to be computed.

III. RELATIVISTIC FEATURES

The non-relativistic electronic Schrödinger Hamiltonian does not provide accurate energies for all elements of the periodic table. In fact, it is only sufficiently accurate for molecules containing atoms with small nuclear charges such as carbon, hydrogen, oxygen, and nitrogen. Qualitatively speaking, the lack of Lorentz invariance of the Schrödinger equation is not significant for molecular systems with small nuclear charges, so relativistic effects on observables may be neglected. A quantum mechanical theory that is valid for the full periodic table of the elements is provided by the (first-quantized) relativistic many-electron theory based on the Dirac one-electron Hamiltonian which accounts for kinematic relativistic effects as well as for spin-orbit coupling. 50

In this section we present (a) a linear scaling implementation of the computation of the DKH correction to one-component wave function model and new capabilities associated with calculations related to magnetic properties, in terms of (b) calculation and visualization of natural orbitals and natural spin orbitals and (c) anisotropic magnetic properties.

A. Local Relativistic Exact-Decoupling for Energies and Properties

The Dirac Hamiltonian is composed of four-dimensional operators. As a consequence, the electronic one-electron wave function must be described by four scalar functions, which form the so-called four-component spinor, or 4-spinor in short. Compared with one-component Schrödinger-based orbitals, these 4-spinors require not only more basis functions, but must also obey a special relation among the different components, called kinetic balance. Although relativistic quantum chemistry has found ways to deal with all difficulties that emerge from such a first-quantized four-component theory, it is desirable to find a representation of the relativistic electronic Hamiltonian, which is sufficiently accurate and easy to interface with the conventional non-relativistic quantum chemical methods. Such a reduced dimensional Hamiltonian, from which the charge-conjugated degrees of freedom have been eliminated

to produce an electrons-only Hamiltonian, features further advantages over the full fourdimensional Dirac Hamiltonian. Most importantly, at most two-component spinors are necessary -only one-component functions if spin-orbit coupling can also be neglectedso that the four-index transformation of electron-electron repulsion integrals, which precedes a standard ab initio electron-correlation calculation, is reduced in computational cost to that of the non-relativistic Schrödinger case.

Our new implementation of relativistic Hamiltonians employed in the Molcas package comprises the one-component exact decoupling X2C 51-54 , Barysz-Sadlej-Snijders (BSS) [55][56][57] and arbitrary-order Douglas-Kroll-Hess (DKH) [58][59][60][61][62][63][64][65] Hamiltonians as well as the atomic mean-field description of spin-orbit coupling interaction treated as perturbation 66 . Moreover, local approaches for calculating and applying the exact-decoupling transformations are employed to reduce the computational cost for polyatomic molecules dramatically. This so-called DLU local decoupling scheme 67 produces a marginal loss of accuracy compared to the full transformation and can be applied to all three Hamiltonians: X2C, BSS, and DKH.

A similar approach was proposed by Seino and Nakai 68,69 for the infinite-order DKH Hamiltonian. For actual calculations, we recommend the X2C Hamiltonian -possibly combined with the DLU scheme for large molecules. As a caveat, we may note that the DLU scheme may require tailored basis sets with not too diffuse functions or the elimination of these very diffuse basis functions from the basis set. 70 This basis-set issue may arise for huge basis sets not optimized for the application in local exact-decoupling calculations.

The central idea of exact decoupling approaches is to block-diagonalize the four-component Dirac Hamiltonian h D by a unitary transformation U The decoupling transformation U in equation ( 17) is not unique. One can easily observe that applying further two independent two-component unitary transformation (one for h + and another one for h -) to h bd still result in a block diagonal structure. In fact, an infinite number of exact decoupling transformations exists. 59 In practice, three principal ways of unitary decoupling have been realized. In DKH theory, the unitary transformation U is decomposed into a sequence of transformations Although the computational cost for relativistic transformations is small compared to that for two-electron integrals, it will increase rapidly if the molecule becomes larger. Since the relativistic effects are mainly localized in atomic centers, the local (atomic) DLU scheme 67

h bd = U h D U † =    h + 0 0 h -    , ( 17
• • • U 2 U 1 U 0 ,
was proposed to reduced the computational cost. It is important to note that the local structure should be exploited at the level of the unitary transformation rather than for the Hamiltonian. 67 The unitary transformation matrix is then decomposed to the following form

U = U AA ⊕ U BB ⊕ U CC ⊕ • • • , ( 18
)
where labels A,B,C denote different atomic blocks. It was found that the DLU approximation introduces very small errors for total energies, which are less than 0.01 milli-hartree for molecules including heavy atoms. The evaluation of relativistic picture change corrected molecular properties had been implemented 63,73,74,76 following the general prescription

ô = ij γ ij ψ i |(U ôU † )|ψ j ( 19
)
where ψ i denotes the relativistic orbitals, ô a one-electron property operator of observable o, and γ ij is a generalized occupation number. In Eq. ( 19), U is directly employed to account for the picture change of molecular property operators. Therefore, the DLU approximated U can be directly substituted to calculate the picture-change-corrected properties with lower cost. The errors for molecular properties due to DLU are also very small. 67 It is important to note that for every physical observable the change of picture must be taken into account. 77

B. Natural orbitals and spin orbitals from SO-RASSI calculations

To facilitate recent studies of magnetic properties of various actinide systems 78-80 , we have implemented the calculation of the u = x, y, and z components of the spin magnetization, and the calculation and visualization of natural orbitals (NOs) and natural 'spin' orbitals (NSOs), within the spin-orbit (SO) restricted active space state-interaction level (RASSI) 66 framework. This allows extracting chemically intuitive information about the electronic structure and magnetic properties directly from SO-RASSI calculations. Illustrative examples are provided below.

In terms of a two-component relativistic SO-RASSI wave function ψ, the electron density ρ(r) and the components m u (r) of the spin magnetization can be defined as

ρ(r) = N ψ † ψ dτ (20a) m u (r) = 2 ψ † Ŝu ψ dτ (20b)
The notation dτ indicates integration over all but one spatial electron coordinates, and the In the SO-RASSI calculations, the electron density ρ and the components m u of the spin magnetization vector field are represented in the atomic orbital (AO) basis set in the form of real symmetric matrices. Diagonalizing these density matrices gives a set of real orthonormal NOs ϕ p with eigenvalues n p (NO occupations) and a set of real orthonormal NSOs ϕ u p with eigenvalues n u p ('spin populations')

notation ψ † • • • ψ indicates 'integration'
ρ(r) = p n p [ϕ p (r)] 2 with p n p = N (21a) m u (r) = p n u p [ϕ u p (r)] 2 with p n u p = 2 S u (21b)
The index p spans the dimension of the molecular orbital basis.

In Equations ( 21), the n p and n u p may be non-integer, and the n u p can be positive as well as negative. For a closed-shell Hartree-Fock (HF) reference, n p = 2 for all occupied orbitals, zero otherwise. With an added unpaired electron, for instance, one gets a doublet. For a restricted open-shell HF (ROHF) scalar relativistic reference (no SO coupling), one would have n p = 1, n z p = ±1 for the unpaired orbital, depending whether M S is ±1/2, and n p = 2, n z p = 0 for the remaining occupied orbitals. A spin-triplet would afford two orbitals with n z p = ±1 for M S = ±1. Dynamic and nondynamic electron correlation as well as SO coupling causes the eigenvalues to deviate from these reference values. SO coupling also causes S z to deviate from M S as the spin projection ceases to be a good quantum number. Zeng et al. 81,82 previously devised two-component 'natural spinors', which are generally complex, and reserved the term 'natural orbital' for an eigenvector of the density matrix of a non-SO calculation. The NOs defined here are determined from the symmetric part of the one-particle density matrix of the SO calculation after integrating over the spin degrees of freedom. This is the density matrix used to calculate expectation values of real spin-free operators in the SO-RASSI step.

Examples are shown in Fig. 6, which are taken from a recent comparative study of the magnetic properties of UO 2 + (5f 1 ), NpO 2 2+ (5f 1 ), and PuO 2 2+ (5f 2 ) and equatorially coordinated complexes of these ions. 79 For reasons of brevity we only discuss m z and the associated NSOs of the actinyl species. In the absence of SO coupling (i.e. in a SR or non-relativistic framework), the 5f φ orbitals of the 5f 1 systems would have n z φ = 1/2 each, due to the orbital degeneracy of the 2 Φ ground states. Under the SO interaction, f φ mixes with f δ of opposite spin, which leads to a reduction of the positive n z φ of the U and Np 5f The action of the crystal field suppresses the effect from SO coupling in this system to a large degree, and correspondingly the n z p are close to the SR limit. Another interesting case, from Ref. 80, is (C 5 (CH 3 ) 4 H) 3 UNO, which is so far the only experimentally characterized U(IV) complex exhibiting temperature-independent paramagnetism up to room temperature. 83 The ground state is non-degenerate and derives from a SR spin-singlet. Selected NOs and their occupation numbers are shown in Fig. 7, comparing SR with SO. The π bonds formed by the uranium 5f π orbitals and the NO ligand are evident. The system has multireference character, as indicated by the relatively large occupations of the antibonding π * combinations, leading to an effective U-NO bond order of only 1.31 instead of 2. At the SO level the bond order is further reduced to 1.18 because SO coupling mixes 5f π with 5f δ , with the latter being non-bonding with respect to the NO ligand. The SO-induced increase of the 5f δ orbitals at the expense of 5f π is clearly visible in Fig. 7.

C. Calculation of magnetic properties in Molcas. The single_aniso module

The multireference CASSCF/(CASPT2)/RASSI-SO 84 spin-orbit calculations proved suitable in providing accurate solutions for arbitrary large complexes, in particular, those containing strongly anisotropic metal sites like lanthanides. This computational approach provides the spin-orbit eigenstates, and the matrix elements of the orbital l and spin ŝ momenta in the basis of multiplet eigenstates of the complex. The accuracy of their estimation depends on the quality of ab initio calculation. Given the strong interest from the molecular magnetism community for accurate calculation of magnetic properties of complexes, we have recently developed and implemented in Molcas computational tools allowing the investigation of anisotropic magnetic properties of complexes using the results of ab initio calculations. These developments for mononuclear complexes have been implemented into the single_aniso module, available in standard version of Molcas package.

The observables calculated by the single_aniso module can be divided into two groups, (i) the magnetic properties and (ii) the parameters of effective (magnetic) Hamiltonians of complexes and fragments (see Fig. 8). Below we review typical predictions obtained with this module. The current version of the single_aniso module is able to compute the following parameters of magnetic Hamiltonians:

Derivation of pseudospin magnetic

• First-rank (linear after pseudospin) Zeeman splitting tensor g αβ , its main values (g X , g Y , g Z ), main magnetic axes (X m , Y m , Z m ) in a given coordinate frame and the sign of the product of the main values g X g Y g Z . 87

• Second-rank (bilinear after pseudospin) zero-field splitting tensor D αβ , its main values (D X , D Y , D Z ) and main anisotropy axes (X a , Y a , Z a ).

• Higher-rank ZFS tensors (D 4 , D 6 , . . ., etc.) and Zeeman splitting tensors (G 3 , G 5 , . . ., etc.) for complexes with S > 3/2 and S > 1, respectively, and weak to moderate spin-orbit coupling effects.

• Angular moments along the main magnetic axes.

An example of derivation of the parameters of magnetic Hamiltonians for a complex with weak spin-orbit coupling effects is shown in Fig. 9. 85 The calculated main values of D and g tensors compare well with experiment. Note that, in general, anisotropy axes and magnetic axes are not parallel, an effect which cannot be caught by ligand field models. 86 An example of calculation of the g tensor for the ground Kramer's doublet of a Ln 3+ complex is shown in Fig. 10.

We should stress high accuracy of calculated magnetic properties of lanthanides by Molcas. This was first proven for individual Dy centers in Dy 3 triangles, for which ab initio calculations predicted the directions of local main magnetic axes 88 very close to the experimental results 89 . Comparable accuracy of ab initio calculated directions of local anisotropy axes have been demonstrated for other mononuclear lanthanide complexes as well by direct comparisons with the results of single-crystal investigations. [90][91][92] Ab initio computation of magnetic properties of complexes. [START_REF] Ungur | Chibotaru Computational Modelling of the Magnetic Properties of Lanthanide Compounds in Lanthanides and Actinides in Molecular Magnetism[END_REF] The code is able to compute all basic static magnetic properties (see Fig. 8): field-and temperature-dependent molar magnetization M mol (B, T ) for powder, magnetization vector for a specified direction of applied magnetic field M( B), Van Vleck magnetic susceptibility tensor χ αβ (T ) and powder susceptibility χ(T ). For the molar magnetization, the Zeeman matrix is built on the basis of several low-lying spin-orbit energy states (user-defined size) and diagonalized. The resulting eigenstates are used to compute molar magnetization where Zeeman states are being populated according to the Boltzmann distribution law. Contribution to the magnetization arising from the states which did not enter explicitly into the Zeeman matrix is considered within the second-order perturbation theory. Molar magnetic susceptibility (tensor and powder) is computed using the zero-field limit of the molar magnetization, where we implemented expressions not requiring Zeeman Hamiltonian to be explicitly computed, enabling much faster computation of these properties compared to the computation of molar magnetization. Intermolecular interactions are accounted for in a mean-field approach, using one single parameter. Fig. 10 shows calculated powder Van Vleck susceptibility and field-dependent magnetization for a DyZn 3 complex and their comparison with experiment.

As a further example, Fig. 11 shows the polar plot of angular dependence of single-crystal susceptibility for two molecules in the unit cell of the crystal. The relative shift of the ab initio and experimental plots is due to the deviation of the calculated main anisotropy axes FIG. 11: a) Orientation of the main magnetic axis in the ground Kramer's doublet of [Cp * ErCOT] ( S = 1/2). 95 Color scheme: Er, purple; C, yellow; H, white. (b) Polar plot of the angular dependence of the single-crystal magnetic susceptibility recorded in the a -b plane, comparison between the results of the ab initio calculation for two molecules (Eq. 22) in the unit cell (left and right plot, respectively) and the contributions extracted from the fit of the experimental data (color code in insets). Computational details: ANO-RCC-VTZP basis for all atoms, active space CAS (11,7), comprising the 4f shell of the Er 3+ , 35 spin quartet states and 112 spin doublet states were optimized in state-average CASSCF calculations and mixed by spin-orbit coupling in rassi. Adapted from Ref. 95. of the two molecules from the experimental ones by a few degrees. 95 These plots are obtained from the calculated susceptibility tensors χ αβ using the expression:

χ( n) = α,β cos θ α χ αβ cos θ β ( 22
)
where n is the direction of the applied field and cos θ α , (α = x, y, z), are the directional cosines.

Ab initio crystal field for lanthanides. [START_REF] Ungur | Chibotaru Computational Modelling of the Magnetic Properties of Lanthanide Compounds in Lanthanides and Actinides in Molecular Magnetism[END_REF]96 The main difficulty encountered by phenomenological crystal-field models is the large number of unknown crystal-field parameters which are to be extracted from available experimental data (27 parameters for lan-thanides). This over-parameterization problem is usually solved by simply neglecting most of the crystal-field parameters, the remaining ones being found by the least-squares fitting procedure. Besides, the problem of determination of the true set of parameters always exists with these models.

Ab initio calculations via RASSCF/RASSI-SO approach are able, in principle, to provide accurate energies and wave functions of low-lying multiplets in lanthanide complexes. Due to the fact that the crystal field description of lanthanide complexes is still of great interest, we have developed a scheme to project the complete ab initio energy matrix (re-written in the basis of J-pseudospin basis) on the crystal-field Hamiltonian, acting on the ground J-multiplet of the lanthanide:

ĤJ CF = n,m (B m n O m n ( Ĵ/J) + C m n Ω m n ( Ĵ/J)), (23)
where O m n and Ω m n are Stevens operators 97 containing components of angular momentum operator Ĵ acting on the ground atomic multiplet, B m n and C m n are crystal-field parameters. The developed approach allows the computation of an unique set of all 27 crystal-field parameters. Another advantage of this procedure is that it permits to project out the crystal field contributions (23) from the total energy matrix and to assess the importance of higher rank (n > 6) non-crystal field contributions. The program also gives the ab initio wave function of the ground J-manifold in the basis of the J-pseudospin eigenstates. [START_REF] Marx | [END_REF][99][100][101] These parameters can be further employed in standard crystal field programs to simulate the low-lying ab initio results. 101 Ab initio blocking barriers for single-molecular magnets. 102 We have also implemented an approach allowing investigation of the structure of the blocking barriers of magnetization in molecular magnets using the results of ab initio calculations. The structure of these barriers is defined as the most efficient relaxation paths from a state with maximal magnetization in the ground exchange doublet to the time-reversed state, corresponding to reversed magnetization. The two basic relaxation mechanisms are due to quantum tunneling of magnetization and the spin-phonon transitions. The rate of the latter is proportional to the matrix elements of the transversal magnetic moment connecting the states with opposite magnetization. The blocking barrier is defined by the shortest path where these quantities are the largest 99,102,103 (see Fig. 12). The method is currently used

IV. MOLECULAR DYNAMICS FEATURES

The ability to simulate dynamical processes is vital to study time-dependent phenomena and can be used in combination with a quantum chemistry code. This capability to run so called "on the fly" molecular dynamics simulations is now present in the new Molcas version. In the following we present some implementations associated with (a) semiclassical adiabatic and nonadiabatic molecular dynamics and tools to generate (b) diabatic states and potential energy surfaces, including (c) an interface to the program Cobramm.

A. Molecular Dynamics Simulations

The new dynamix module provides the capability to run molecular dynamics (MD) simulations. In a MD simulation the nuclei are moved classically exploiting the Born-Oppenheimer (BO) approximation. Their equations of motions are solved numerically using the velocity-Verlet algorithm. 106 The algorithm requires coordinates, velocities and forces as input. The dynamix module can be used with any electronic structure method available in Molcas. In addition the simulations can also account for the surrounding. Solvent effects can be considered implicitly using the reaction field keyword in gateway. Also hybrid QM/MM calculations are feasible in conjunction with the espf program. Temperature can be taken into account through the coupling to a thermostat. 107

In order to introduce the temperature as a constraint in the system and take into account the thermal energy redistribution in condensed media, the Nosé-Hoover chain of thermostats approach has been implemented. It allows one to perform simulations at constant temperature that produce statistics belonging to the canonical ensemble (N V T ) 108 . In this case the initial conditions are prepared by getting random velocities from a Maxwell-Boltzmann distribution.

Nonadiabatic MD simulations

The nonadiabatic molecular dynamics (NAMD) simulation of many-electronic-state problems requires the simultaneous description of different electronic states as well as their corresponding couplings. In these simulations, a trajectory surface hopping (TSH) algorithm allows nonadiabatic transitions between different states.

For this purpose two algorithms are available in Molcas.

The first option takes a deterministic approach to decide about surface hopping based on the approximate evaluation of nonadiabatic coupling as introduced by Robb and coworkers. [109][110][111][112] The details of the derivation can be found in Ref. 113 and in the supporting information of Ref. 114. Briefly, an approximate coupling between two states, k and j, is computed using the numerical approximation

Φ k (t) ∂ ∂t Φ j (t) ≈ Φ k (t) |Φ j (t + ∆t) ∆t ( 24
)
where Φ k is the wave function and ∆t the time step. The coupling is used to estimate the interaction of the two states and to decide if a hop between the states occurs (Fig. 13). If a given threshold is reached a transition between k and j is induced.

The second option is the widely-used Tully's fewest switches algorithm 115 , available in the new module surfacehop. This approach requires the coupled propagation of both the nuclei and the electronic motion. The former is performed by integration of Newton's equations of motion for the nuclei. The step size for integration is on the femtosecond timescale.

FIG. 13: The wave function change of the two states j and k with time along the trajectory, permits the numerical estimation of the nonadiabatic coupling as the scalar product between Φ k (t) and Φ j (t + ∆t).

However, the electron motion is much faster. Hence, it is carried out by integrating the time-dependent Schrödinger equation (TDSE) on the (tens of) attosecond timescale (see eq.

25)

.

ih ċk = j c j V kj -ih Ṙd kj ( 25
)
where c j are the expansion coefficients and Ṙ are the velocities of the nuclei. The indices j and k are labeling different electronic states. The following terms:

d kj = Φ k (r; R)|∇ R Φ j (r; R) (26)
are the nonadiabatic couplings, and

V kj = Φ k (r; R)|H(r; R)Φ j (r; R) (27)
are the electronic Hamiltonian elements for a given geometry R. The wave function in the present implementation, Ψ(r; R, t), is expanded as a linear combination of the different adiabatic electronic states Φ j (r; R) defined by CASSCF or CASPT2 wave functions:

Ψ(r; R, t) = j c j (t)Φ j (r; R).
Within the adiabatic basis currently implemented, the d kj terms are key parameters in determining surface hopping, as they are the only terms promoting the electronic states mixing. Nevertheless, in the integration of equation ( 25) the Ṙd kj terms, rather than the where the Ṙd kj (R) elements are calculated using the right-hand side of equation ( 24). These terms are determined numerically along the trajectory at each step of the nuclear motion integration defined by ∆t coinciding with the exact calculation of the CASSCF or CASPT2 wave function. Nevertheless, the time-step of the electronic integration is typically at least one or two orders of magnitude smaller than ∆t. Consequently, to get accurate results in the integration of the TDSE, an interpolation between consecutive values of Ṙd kj is required.

The implementation chosen for Molcas-8 follows the interpolation-extrapolation scheme of Hammes-Schiffer and Tully 116 for the determination of time-derivative coupling terms used in eq. ( 25). Additionally, in order to preserve the consistency of the population density matrix, we have implemented the decoherence correction proposed by By integrating eq. ( 25) the populations of each electronic state |c i | 2 are obtained. Depending on the variation of these populations, the fewest-switches algorithm 115 stochastically decides whether the trajectory should continue on the same state or should jump onto another one. In case of hop, the next energy gradient for the nuclear motion is calculated using the new state while the nuclear velocities are scaled to preserve the total energy.

Applications

The NAMD implementation in Molcas allows to study problems in photochemistry and photobiology. It is currently being used by several groups for simulating excited state dynamics in many chemical and biological systems. [118][119][120][121][122][123][124][125][126][127] The success stems from the variety of electronic structure methods for excited state calculations in Molcas.

Applications in the gas phase were done to a molecular motor with a chiral hydrogen bond environment (see Fig. 14) 127 and to a chiral model indanylidene-pyrrolinium (NAIP) switch. 125 In these applications CASSCF level of theory was used. However, the first CASPT2 trajectories were performed for a halogenated methyl radical CF 2 I. 122 Recently, also a truncated retinal model with three double bonds was subject to a NAMD investigation. 118 In this study 100 CASSCF and 100 CASPT2 trajectories were computed and analyzed.

In addition, the modular structure of Molcas permits combination with other modules such as espf to do simulations in the QM/MM framework. The first QM/MM application was the calculation of the photoisomerization mechanism in rhodopsin and bathorhodopsin. 124 Further studies were focussed on the retinal in anabaena sensory rhodopsin protein 121 and NAIP molecular switches in methanol solution. 119,126 Most recently, it has been possible to run semiclassical QM/MM trajectories at the CASPT2 level of theory on proteins containing reduced chromophore models. 128 These studies did not only provide a proof-of-principle for the feasibility of these calculations but also were shown to yield mechanistic informations close to the one obtained with the full chromophore (Fig. 15).

B. A New Approach For Generating Diabatic States: The

Dipole-Quadrupole (DQ) Method

Diabatization is a method for transforming from an adiabatic basis to a diabatic basis. In an adiabatic basis, the potential energy surfaces have cuspidal ridges at conical intersection seams, and the nonadiabatic coupling of the adiabatic states is singular at these seams.

A diabatic basis is more convenient for dynamics calculations; diabatic potential energy surfaces are smooth and their couplings are nonsingular. References for various methods for diabatization are assembled and discussed elsewhere [129][130][131][132][133] .

The dipole-quadrupole (DQ) method for direct diabatization was developed recently 134

and has now been added to Molcas. In the DQ scheme, a diabatic state is defined by an adiabatic-to-diabatic transformation where ψ I is an adiabatic state and φ A is a diabatic state; therefore, the diabatic states span the same space as the N chosen adiabatic states. The adiabatic states are eigenfunctions of the Born-Oppenheimer electronic Hamiltonian, and diabatic states have smoothly varying electronic wave functions and energies and negligible momentum couplings. 135 The motiva-tion for working in a basis of diabatic states is the relative ease with which one can obtain an analytical representation and carry out dynamics calculations. In the adiabatic representation, states and surfaces have discontinuous gradients, and the momentum couplings are vectors with singularities. However, diabatic states are smooth by definition and are coupled by the off-diagonal elements of the electronic Hamiltonian, which are smooth, scalar quantities.

φ A = N I=1 ψ I T IA (28)
The DQ method involves the dipole moment with one choice of origin, R 0 , and the quadrupole moment defined with one or more choices of origin; the various origins for the quadrupole terms are denoted as R j with j=1, 2, . . . , N or . The transformation matrix is determined such that the transformation maximizes f DQ , defined by

f DQ = N A   | φ A |R -R 0 | φ A | 2 + Nor j α j φ A |R -R j | 2 φ A 2   , ( 29
)
where α j is a parameter and R specifies the positions of all particles.

This method has been implemented in Molcas as a subroutine in the rassi 66,136,137 module and works with CASSCF, 138 RASSCF, 5 CASPT2, 17,139 or RASPT2 7 wave functions. One first computes dipole and quadrupole matrices in rassi; then these quantities are passed into the DQ subroutine, which finds a transformation matrix that maximizes Equation ( 29). This transformation matrix is T in Equation ( 28), and it yields diabatic states and their energies and couplings at a single point.

The function in Equation ( 29) was inspired by previous work that only used the dipole moment. 140-144 A problem with only using the dipole is the lack of generality when treating reactions because many problems of interest have electronic states with similar dipole moments. By introducing the quadrupole, the DQ method can alleviate this problem and allow diabatic transformations for a larger number of systems or for larger number of states for a given system. This can be illustrated by diabatizing the three lowest-energy

E (eV) R Li-H (Å) V 1 V 2 V 3
FIG. 16: The adiabatic (V 1 , V 2 , and V 3 ) and diabatic (U 1 , U 2 , and U 3 ) potential energy curves of the three lowest-energy 1 Σ + states for LiH. The diabatic states were computed with Boys localized diabatization.

3s Li , and 3p Li . For the DQ method, we used the Li atom as the only choice of origin for the quadrupole integrals (i.e. N or = 1), and we used an α value of 10 a -2 0 . For the dipole integrals, R 0 was the position of the Li atom; however, the results for LiH are independent of R 0 , since it is a neutral molecule. The comparison between using Boys localized diabatization 142 which only uses the dipole (i.e., α = 0), and using the DQ method are shown in Figs. 16 and 17. The adiabatic potential curves (V 1 , V 2 , and V 3 ) are the same in Figs. 16 and17, and the diabatic potential curves (U 1 , U 2 , and U 3 ) are defined by their character near the equilibrium geometry: state 1 is ionic with a 1s 2 H valence configuration, state 2 is covalent with a 2s Li 1s H valence configuration, and state 3 is also covalent with a 2p Li 1s H valence configuration. The dipole moment is able to differentiate between the ionic and covalent states; however, the two covalent states need the quadrupole in addition to the dipole; without the quadrupole one obtains a linear combination of the physical diabats at large inter-nuclear separation, as seen in Fig. 16.

The DQ method is applicable to general photochemical reactions, not just dimers and electron transfer systems. As an example, we also showed that one can also obtain good results for the challenging case of the photodissociation of phenol. 134 This is very encouraging and indicates that the DQ method can be a valuable tool in studying photochemical E (eV)

R Li-H (Å) V 1 V 2 V 3
FIG. 17: The adiabatic (V 1 , V 2 , and V 3 ) and diabatic (U 1 , U 2 , and U 3 ) potential energy curves of the three lowest-energy 1 Σ + states for LiH. The diabatic states were computed with the DQ method with α = 10.0 a -2 0 .

reactions.

C. The Molcas-Cobramm interface

This interface provides easy access to the recently implemented Molcas Cholesky type decomposition-based and FNO-CASPT2 routines in a combined quantum mechanical molecular mechanical framework, with automated MM topology and QM/MM input generation along with a variety of tools to evaluate spectral properties of multichromophoric systems and to process spectral data.

Cobramm is a QM/MM program adopting an ONIOM-like 147 subtractive approach within an electrostatic embedding scheme. It connects the features of various QM codes and the molecular mechanics program Amber (Figure 18). 148 It has been specifically designed to investigate excited state properties and photodynamics of biomolecules. The major aim of

Cobramm is to keep complicated and tedious input and setup procedures at a minimum.

The QM and MM program input is the same as for the specific codes used, and QM and QM/MM optimizations and molecular dynamics can be readily accessed using a minimal number of keywords. Providing a PDB structure and the libraries for non-standard residues (if present) as inputs, the user can easily select the QM and MM regions, the frozen/movable residues and the type of computation to be performed. Molcas is used within the Cobramm module, which is available free of charge upon request. 149

The Molcas-Cobramm interface features usage of the newly developed Cholesky decomposition-based techniques, reaction fields, numerical nonadiabatic couplings and gradients to be used in free and constrained optimizations, numerical frequencies and molecular dynamics. Numerical computations can be performed with a simple parallelization routine distributing the workload among the available CPUs with or without using MPI extensions.

By means of the Molcas-Cobramm interface, the user can take advantage of various helpful features available for the nonadiabatic molecular dynamics. For instance, an energy criterion can be set in order to change the timestep size during simulation (e.g it can be reduced to map with a finer grid strongly coupled regions of the PES), to activate/deactivate the calculation of the derivative couplings (to save computational time), and to switch from state-averaged to state-specific computations (i.e. to avoid possible wave function convergence problems). Moreover, for CASSCF molecular dynamics, CASPT2 single point calculations can be automatically performed on top of the MD snapshots (selected with the desired timestep/frequency). The Molcas-Cobramm interface also offers specific tools allowing automatic analysis of geometrical data and generation of statistical data (including hopping times, distribution of geometry parameters like bond lengths or torsions, etc.).

Finally, the Molcas-Cobramm interface includes a plug-in that makes use of the Spectron code, a program for generating linear and nonlinear response for molecular aggregates and quasi-particles, subject to interaction with ultrashort laser pulses, by solving the Liouville-von Neumann equation for the density matrix. 150 It can be applied to calculate the third-order nonlinear response of (multi)chromophoric systems, allowing for simulations of two-dimensional electronic spectra (2DES), see Figure 18. Also this code is avaliable free of charge.

In the following section, we describe how the recent Molcas developments reported above open the way to accurate first-principles simulations of this emerging nonlinear spectroscopic technique.

Applications: SOS//QM/MM simulations of nonlinear optical spectroscopy

We present here two examples of simulated two-dimensional electronic spectra in the ultraviolet (2DUV) using the aforementioned Molcas-Cobramm interface. 2D electronic FIG. 18: Schematic representation of the Molcas-Cobramm interface and its application to simulations of nonlinear 2D electronic spectra (via the Spectron code plugin).

spectroscopy emerged in the last decade as a nonlinear optical technique based on sequence of ultrashort laser pulses, with high temporal and spectral resolution that provides fundamental insight into coherent excited state dynamics, interstate population transport, solvent reorganization timescales and non-covalent interactions (π-stacking). Experimental 2D electronic spectra contain a wealth of (congested) information associated with time-dependent and state-specific "fingerprint" signals, such as stimulated emission (SE) and excited state absorption (ESA), calling for theoretical interpretations that face several computational challenges. Recently, a computational protocol for simulation of 2DUV spectra of complex systems like peptides or poly-nucleotides has been introduced [151][152][153][154] , based on the sum-overstates (SOS) approach, where an ab initio characterization of the excited state manifolds of multichromophoric systems within a quantum mechanics/molecular mechanics (QM/MM) scheme is coupled with the nonlinear response theory. 155 These electronic structure calculations can be performed by means of the Molcas-Cobramm interface. To this scope, the CASSCF/CASPT2 methodology is the preferred choice, in light of the often multiconfigurational nature of the excited states and the need to describe valence, charge-transfer and doubly excited states at equal footing. The simulations of nonlinear spectra within the CASSCF/CASPT2 framework pose formidable challenges, due to the requirement of large active spaces (necessary to include high-lying excited states, above twice the energy of the first bright state) and the necessity of including up to hundred excited states in the state-average procedure. With the implementation of the state-averaged RASSCF/RASPT2

technique and the parallelization of the CD-RASPT2 method this is no longer an impossible tasks, also for relatively large chromophoric systems. In fact, exploiting Molcas versatile utilities we have recently presented a methodology for performing large active space/large basis set computations for the monomeric building blocks of complex multichromophoric systems (i.e. amino acids and nucleobases) [151][152][153][154] , providing a library of reference values for transition energies and dipole moments that cover the near-IR to far-UV spectral window for application in linear and nonlinear spectroscopy. This set of data has been used to assess the applicability of active space/basis set truncation and semi-empirical correction schemes available in Molcas for reproducing the reference electronic spectra at low cost. A similar strategy has been applied in a recent study by Pulay et al. 156 The low-cost computational recipes have been then employed for simulation of 2D electronic spectra of relatively large chromophoric systems within the SOS//QM/MM scheme. These results point out the potential of UV-active aromatic side chains as local highly specific markers for tracking dynamics and structural rearrangements in proteins by means of 2DUV spectroscopy.

The theoretical study of the model proteic system accounted for the simulation of one of the possible 2DES experiments, where the excited state dynamics is avoided by controlling the time intervals of the laser pulse sequence. However, 2DES can be used also to track ultrafast events along the excited state dynamics of multichromophoric systems. As depicted in Fig. 19b, the nuclear degrees of freedom along the dynamics in the photoactive state(s) can determine electronic fluctuations that can be resolved by 2DES experiments. Using the Molcas-Cobramm interface the fluctuations of higher excited states of pyrene (Fig. 19b), an aromatic chromophore characterized by an ultrafast decay from the bright S 2 (L a ) state to a dark S 1 (L b ) state, have been characterized along a single-trajectory propagated at the SA-CASSCF(2,2) level and nonadiabatic molecular dynamics refined at the SA-60-RAS(4,8;0,0;4,8)/SS-RASPT2 level. 157 Accounting for these fluctuations allowed a proper description of the coherent nuclear motion that strongly affects the time-resolved 2D electronic spectra. The nice agreement between the experimental and theoretical spectra (Fig. 19b) validates the formulated mathematical models and the results obtained with the Molcas-Cobramm interface, including the nonadiabatic molecular dynamics and the stateof-the-art electronic structure computations. The same approach can be employed to study the photochemistry and photophysics of other photoactive molecular systems embedded in complex environments, including biological photoreceptors (e.g. retinal proteins [158][159][160] ) and polynucleotides in DNA/RNA. 153,154

V. MOLECULAR GRADIENTS AND OPTIMIZATIONS

The use of gradients with respect to the nuclei is essential for the determination of molecular structure, to establish reaction paths and as the driving force for molecular dynamics [START_REF] Aquilante | [END_REF]163,164 One challenge remained, however. The CD technique is usually formulated as a numerical procedure, which makes it difficult to express analytical derivatives. 165 As of Molcas-8, analytical derivatives are available for most CD-based methods, as we shall outline below.

CD consists in expressing the ERIs in terms of a sum of outer products of so-called Cholesky vectors:

(µν|λσ) ≈ M J R J µν R J λσ . ( 30
)
The Cholesky vectors are obtained through a recursive numerical procedure that simultaneously determines the number of vectors, M , in accordance with the requested accuracy. 163,164 For typical applications the value of M is 3-5 times the number of atomic orbital (AO) basis functions. Individual ERI errors caused by the decomposition are controlled by a single parameter (the Cholesky threshold), thus offering the user a simple means for accuracy control. With the algorithmic developments 33,166,167 introduced by the Molcas community, CD techniques have witnessed in the last decade a growing interest as means to develop highefficiency and controlled-approximation methods for electronic structure theory calculations.

The computational benefits of using CD-based algorithms vary depending on the type of quantum chemistry method in which CD is introduced. Details on this topic have already been presented in the literature along with extensive benchmarks on CD accuracy. [168][169][170] However, a major breakthrough for CD-based quantum chemistry has been that of realizing the common nature 171 that CD shares with a more popular ERI approximation technique known with the names of resolution of identity (RI) or density fitting (DF). Indeed, it can be shown that the CD representation of Eq. ( 30) is equivalent to RI/DF employing a set of auxiliary product functions (h J , also known as Cholesky basis) defined along the decomposition procedure, namely :

(µν|λσ) ≈ JK C J µν G JK C K λσ , (31)
where

G JK = (h J |h K ) , ( 32
) J C J µν G JK = (µν|h K ) . ( 33
)
One way to fully exploit this connection has been that of using CD to define atomic auxiliary basis sets for standard RI/DF algorithms: the atomic CD (aCD) 172 and later the atomic compact CD (acCD) basis sets. 173 The original CD is now called in Molcas full-CD. The main advantage of the CD-generated auxiliary basis sets compared to the standard predefined auxiliary basis sets is that they are not biased towards any method, and their accuracy can be largely tuned by the choice of the decomposition threshold. [168][169][170] But even more important is the fact that the equivalence with RI/DF allowed to devise a way to compute analytical gradients, which had been seen as a major drawback of CD.

By means of Eq. (31) it is in fact possible to approximate the ERIs' first derivatives (e.g., with respect to a geometrical displacement) as:

(µν|λσ) [START_REF] Aquilante | [END_REF] = J C J µν (λσ|h J ) [START_REF] Aquilante | [END_REF]

+ J (µν|h J ) (1) C J λσ - JK C J µν G (1) JK C K λσ , ( 34
)
where the fitting coefficients can be computed either by inverting 33, or directly from the Cholesky (or RI/DF) vectors, as in Ref. 171. Following this insight, gradients for the Coulomb 171 , exchange 174 and MP2 energies 175 were derived, and efficient algorithms are now present in Molcas-8 that show similar speed-ups as for the energy calculation, while hardly affecting the computed equilibrium structures. In-depth analysis of the above implementations is beyond the scope of the present publication. Here, we will only point out that the two-electron contributions to the CD-energy gradients are computed by building the two-and three-index effective density matrix, P JK and P K µν , respectively, followed by contraction with the integral derivatives: [START_REF] Aquilante | [END_REF] -JK P JK (J|K) [START_REF] Aquilante | [END_REF] .

E (1) ∼ 2 Kµν P K µν (K|µν)
(

) 35
The effective density matrices result from specific tensor contractions involving the fitting coefficients and the usual one-and two-particle density matrices. The gradient code can thus be streamlined in the same way as the conventional code, with first computing the (effective) density matrices and then contracting them on the fly with derivative integrals. A major prefactor in the cost for computing the effective density matrices comes from exchange-type contributions. However, an extension of the LK-screening algorithm 166 has been developed in Molcas-8 that reduces such computational costs for medium to large molecules by resorting to a Cholesky-localized orbital picture. 176

Very recently, single-state CD-CASSCF analytical gradients were for the first time employed in order to analyze the geometry of a ruthenium nitrosyl complex. Both the solution of the linear response equations for the determination of the Lagrange multipliers (mclr) and the gradient formation (alaska) benefit from the CD implementation and overall, the implementation reaches a balance between these two steps and the energy calculation itself as illustrated in Fig. 20. One important feature of this alaska implementation is that the time required for computing the effective density matrices and subsequent contraction with the integral derivatives is nearly independent from the size of the active space, as opposite to the conventional calculation (see Fig. 21).

We are confident that with the current implementation of energy and gradients calculation for CD-based approaches, Molcas-8 could potentially mark a new era in the use of multiconfigurational wave function methods for ab initio molecular dynamics studies in photophysics and photochemistry.

B. Constrained numerical gradients and composite gradients

As pointed out above, one of the essential quantities needed from any quantum chemistry package is the energy gradient. For many electronic structure methods it is possible to obtain the gradient based on an analytical formulation, with a computational cost that is roughly equivalent to that of computing the energy. Also, the range of methods for

which Molcas provides such analytic gradients, and their efficiency, is increasing. However, the more sophisticated methods often lack analytical implementations, and to obtain the gradient one must resort to numerical differentiation, requiring a large number (proportional to the number of atoms in the system) of energy calculations. In spite of this very large computational cost, especially considering that it is required for the highest-level methods, numerical differentiation has the advantage of being easily implemented and almost trivial to parallelize.

Another problem faced by geometry optimization methods is the large dimensionality of extended systems. As the number of atoms increases, so does the number of coordinates and degrees of freedom, and locating an optimum structure becomes increasingly more difficult.

A possible way to overcome this problem is reducing the dimensionality of the system by applying constraints or freezing some degrees of freedom that are considered to be secondary

for the problem at hand. Since several years Molcas includes a powerful and versatile tool to perform optimizations with various types of constraints. 179 If a constrained optimization is combined with a method for which analytic gradients are not available, it should be possible to keep down the number of energy calculations by restricting the differentiation to the degrees of freedom that will be actually optimized, as the gradient along the constrained degrees of freedom is typically not needed. This is now done in Molcas. Thus, whenever numerical gradients are computed with constraints, only the unconstrained coordinates are differentiated. This can represent a very important performance improvement when extensive constraints are used. As an example, in a system formed by two molecules with constraints defined such that the molecules remain rigid and only the intermolecular degrees of freedom are optimized, the effective numerical gradient can be computed with only 13 single-point calculations, regardless of the size of the molecules. This should be compared to the 6N -11 (N , total number of atoms) energy calculations that would be needed if the constraints are not taken into account. FIG. 22: Performance of the present method for constrained optimization with numerical gradients for the S22 set of weakly bound dimers. 180 The mean and median energy deviation (a) is shown together with the mean and median RMS distance (b) with respect to the fully optimized structures. The numbers appearing on the x-axis refer to different values of the energy convergence criterion in units of E h . Starting structures for the constrained optimizations were obtained by randomly modifying the intermolecular coordinates of the fully optimized structures. All computations were performed at the MP2/cc-pVDZ level. Note that for the mean and median energy deviations obtained with internal coordinates, one data point (the formic acid dimer) was removed due to anomalous behavior.

In test calculations on the S22 set of weakly bound dimers 180 , we have shown that when using the partial numerical gradients discussed above, the errors introduced in the final energies and geometries are of the order of the convergence threshold 181 (see Fig. 22). At the same time, the computational cost of the constrained optimization is reduced by up to 95%, depending on the number of constraints defined. We emphasize that this reduction is achieved with no change in the input and virtually no effect on the results. Therefore, when one must use numerical gradients, a case can be made to include as many constraints as is reasonable for the system being studied.

For cases where constraints cannot be used because it is necessary, or at least desirable, to relax all the coordinates in the system, it is still possible to reduce the number of energy evaluations by using what can be seen as hybrid or "composite" gradients. In this case a partial numerical gradient is computed as above, but it is combined with an analytic gradient (computed with some other method for which analytic gradients are available) for the degrees of freedom not included in the numerical differentiation. The separation of the degrees of freedom to be differentiated numerically (with the higher-level method) and analytically (with the lower-level method) should be reasonable, and the combination of the two electronic structure methods must be meaningful.

As an example we can consider the calculation of electronic transition energies in 2,3dimethyl-2-butenal (Fig. 23). 181 The geometries in the ground and excited states can be computed with CASPT2, although since analytic CASPT2 gradients are not available in Molcas, the optimization requires 91 single-point calculations per iteration. Alternatively, the geometries can be optimized at the CASSCF level (with analytic gradients) and the transition energies computed with CASPT2 (CASPT2//CASSCF). Table II shows the effect that the different optimized geometries have on the transition energies. Given that the active space spans only the central π system, it can be assumed that the internal structure of the three methyl groups will have little effect on the transition energies. A composite gradient strategy can therefore be adopted, excluding the internal CH 3 coordinates from the CASPT2 differentiation and using some other method's analytic gradients (like CASSCF or MP2) for them instead. In this way, the number of energy evaluations is reduced to 55 (plus a single gradient calculation at the lower level) per iteration. also shown in Table II, with no significant difference from the full CASPT2 optimization.

It should be noted that when MP2 is used for the CH 3 gradient, this is computed at the ground state, even for the excited state optimization. However, this should not be a problem since the electron density in the methyl groups is not expected to be significantly affected by the electronic state.

The composite gradients and the constrained numerical gradients computed in Molcas can be used in both geometry optimizations, as in the above examples, or, especially the composite gradients, in molecular dynamics simulations (see subsection IV A). In combination with the easy parallelization of the numerical differentiation, we expect that these techniques can help extending the use of higher-level electronic structure methods to larger systems.

C. Geometry optimization based on frozen fragments

A study of interaction between different molecules and/or a molecule and a surface quite often requires geometry optimization with internally frozen fragments. The freezing of coordinates can be used not only for increasing the speed of the calculation, but also can be used to perform a model calculation with very few parameters involved.

In the previous section we introduced a method to take advantage of the frozen coordinates to reduce the computational effort to obtain numerical gradients, but nothing was changed in the optimization algorithm. In this section we refer to a novel optimization algorithm specifically designed for frozen fragments and which bypasses the numerical gradient calculation.

The usage of internal coordinates for geometry optimization is the most preferable way of optimization: it allows to minimize the number of degrees of freedom, and it reduces numerical errors since all displacements are done according to the direction of variation of forces (bond distances, and valence angles). In Ref. 182 we introduced a new procedure for numerical geometry optimization of frozen fragments which is based on displacements of internal coordinates (Molcas module geo). In this approach a partial Z-matrix is constructed and it is used to define Cartesian coordinates which correspond to displacements of internal coordinates. An advantage of the proposed technique is that it can be applied to any computational method which produces a total energy for a chosen geometry.

The method can be used for geometry optimization using CASPT2 or coupled cluster

Hamiltonians. In Ref. 182 we reported optimization of a system containing Li + ion and benzene molecule using the CCSD(T) method. The Li + ion was placed at the benzene plane, and during geometry optimization, the ion moved to the position on top of the benzene ring.

The optimization process took 18 iterations. We should note here that using the standard Molcas tool (slapaf) the convergence was never achieved.

Another group of applications where geometry optimization with frozen fragments is essential is a study of interaction of molecules with a crystalline surface. One relatively simple example of such interaction is the reaction of a single water molecule with an ionic oxide surface, e.g. CaO. The surface of the ionic crystal can be presented as a small cluster embedded into the electrostatic field of point charges. According to the model adsorption, one should control the relaxation of the crystal surface. Different steps can be studied: physical adsorption (no relaxation at all), chemical adsorption (relaxation of the molecule is allowed, while the surface is frozen), chemical reaction (both the molecule and the surface can relax) and possible dissociation of the molecule. Optimization with frozen fragments allows to model all these steps. We performed a study of dissociation of a single water molecule at a CaO surface modelled by a Ç cluster using the CASPT2 Hamiltonian (with active space for CASSCF calculation 12/12). A barrier free dissociation of the water molecule has been observed on the CaO surface, where one Ca atom was allowed to move. However, if all atoms on the CaO surface are frozen the dissociation does not occur.

In this case the possibility to perform optimization with frozen fragments not only allows to use very precise multiconfigurational methods, but also was useful to control different steps of the chemical reaction. The current Molcas interface to external libraries is depicted in Figure 25.

The use of a threaded, external BLAS/LAPACK library allows for a speed-up by using multiple cores on a single shared-memory machine. Alternatively, Molcas can now also be configured to use NVBLAS 185 to run specific BLAS operations on a GPU. Parallel linear algebra operations involving multiple machines are enabled through either GA or DGA (which then use an external BLAS/LAPACK library). Currently, caspt2 is the only code that benefits from using the latter functionality (for very large matrix-matrix multiplication and diagonalization).

C. The Molcas/Columbus link

The purpose of this link is to provide easy access to large scale configuration interaction calculations within the Molcas framework. Both variational and perturbational treatment of spin-orbit (SO) coupling at the MR-CISD level of theory is available. [189][190][191] For practical high-level electron correlation treatmentin particular in the context of variational SO-MRCISD -CSF expansion space dimensions between 10 8 and 10 In n-acenes this splitting is found starting for n=4-5 and for periacenes with 5 armchair benzene rings this value is already reached for a n-zigzag length of n=3. 201 In addition, graphene nanosheets seem to possess strong polyradical character 202 along the zig-zag edges which suggests even richer spin-orbit effects in contact with a copper atom.

Outline of the link

In conclusion, one fruitful possibility for the link between Molcas and Columbus has been demonstrated. The combination of the two program packages is transparent for the user since all necessary restructuring of data files is performed automatically.

VII. SUMMARY

The 8.0 release of the Molcas quantum chemistry package is now available, and version 8.2 will be available soon. In this report we have published the highlights of the new features of these versions together with some details of technical improvements. These can be summarized as follows.

For wave function associated developments we have reported on GASSCF, FD and FNO approximations to the CASPT2 approach. First, for the GASSCF method larger active spaces are accessible at reduced computational cost. Perturbation theory approaches can be applied to GASSCF wave functions to recover dynamic correlation energy, and the GASPT2 extension will be released soon. Second, with the implementation of efficient schemes for compacting the inactive and secondary orbital space (e.g. FD and FNO) the description of the dynamic correlation is no longer the bottleneck for the various Cholesky-based multiconfigurational perturbation theories present in Molcas-8. These developments are comple-mentary to the work aiming at the generalization of the RASSCF/RASPT2 scheme to the GASSCF/GASPT2, as well as to new approaches such as DMRG, which have the potential to overcome the hurdle of the exponential scaling of the number of configurations with the size of the active space.

Also new in Molcas is that wave function theory and density functional theory meet in the multi-configuration pair-density functional theory method, MC-PDFT. This represents a cost-effective and efficient approach to recover the entire correlation energy that can be used as valid alternative to CASPT2 for ground state as well as excited state chemistry.

In the area of new relativistic implementations we report on a linear scaling method for the DKH approach and implementations related to calculations including spin-orbit coupling and magnetization. To be more specific, Molcas supports various scalar-relativistic exactdecoupling schemes. Apart from the DKH Hamiltonians, X2C has become available. In order to facilitate the application of relativistic Hamiltonians to large molecules, the atomic decomposition of the unitary decoupling matrix DLU was developed and implemented in

Molcas for all such Hamiltonians. Moreover, the generation of orbitals in Molcas representing the electron density and spin magnetization components facilitates and simplifies the analysis of electronic structure and bonding at the SO-RASSI level. Additional functionality to represent the orbital magnetization is forthcoming.

For issues related to molecular dynamics simulation we report on a new module for such simulations, interfaces to enable the study of nonlinear optics, and methods for diabatization.

In particular we report that, a new module (dynamix) for molecular dynamics simulation has been added. In addition two surface hopping algorithms were implemented: one based on a deterministic and the other one on the stochastic approach, namely the Tully fewest switches algorithm. This permits the computation of nonadiabatic QM or QM/MM trajectories for different statistical ensembles. In the present implementation, NAMD simulations within the frame of microcanonical and canonical ensembles using multiconfigurational wave functions (CASSCF/RASSCF, CASPT2/RASPT2) are available. Furthermore, a computational machinery for studying nonlinear optical spectroscopy and mainly based on the interaction of different software through the QM/MM package Cobramm is present in Molcas. Therefore, nonlinear optical spectra from accurate experiments can be modeled by coupling CASPT2 calculations with response theory through the sum-over-states (SOS) approach, and also within a quantum mechanics/molecular mechanics scheme for the treat-ment of condensed phase and protein environments. Finally, a new diabatization method has been implemented, DQ, which is a straightforward scheme, applicable with any electronic structure method, not restricted to simple reactions, orbital-free, and not path-dependent.

For the exploration of the molecular potential energy surfaces and computations of stationary points we present improvements in the computation of state-average CASSCF analytical gradients and methods for efficient numerical gradient techniques. By using stateaverage CASSCF in association with density fitting and Cholesky generated auxiliary functions we have demonstrated speed-ups in the range of one order of magnitude reduction of the CPU time. Furthermore, using geometrical constraints Molcas can now greatly improve the performance of optimizations and numerical gradient calculations; together with the composite gradient concept this offers a chance to run geometry optimizations and dynamics at the highest levels of theory and at a reasonable computational cost. Columbus interface when more specialized features are to be required. In this way extended nonrelativistic and spin-orbit coupling MRCI calculations with several billion configurations are accessible using today's standard computer hardware.

In conclusion, we note that the development of the Molcas program package is an ongoing process. Here we have reported the developments that have taken place during the last six years. The software additions and improvements include several unique capabilities that we hope will be useful to the chemistry community for simulations of molecular energies, properties, spectra, and reactivity.

FIG. 1 :

 1 FIG. 1: Number of determinants (blue dots) and configurations (green squares) for a complete active space of m electrons in n orbitals, where m = n and M s = 0. Solid lines represent Stirling's approximation (N SD = 2 πn 4 n , N CSF = (1 -( n n+2 ) 2 ) πn 4 n ).

FIG. 3 :

 3 FIG.3: Number of determinants for a restricted active space of n electrons in n orbitals (M s = 0), for various sizes of the inner RAS2 and different excitation levels into RAS3.The solid line represent the upper bound of the number of determinants.

  unrestricted formalism, energies are obtained that are based on an ambiguous state which results from a solution that does not have correct symmetry properties. MC-PDFT combines the strengths of wave function theory and DFT. It uses spin-and space-adapted wave functions and describes correlation effects inexpensively by functionals of spin free densities. The required density functionals are called on-top density functionals, and the first generation of on-top functionals has been obtained by translation of exchange-correlation functions originally developed for KS-DFT. In MC-PDFT, a multiconfigurational wave function is used to compute total density and on-top pair density of the correct symmetry

  energy, and the classical Coulomb energy, and functionals of the total density and ontop pair-density are used to compute exchange and correlation energy contributions. In the current implementation the on-top pair density functionals are obtained by translation (abbreviated "t") of KS-DFT exchange-correlation functionals, which depend on the total density and spin magnetization density, to functionals that depend on the total density and on-top pair density. Three translated functionals are available in the current version of Molcas: tLSDA, tBLYP and tPBE. The method is still developing but the initial results are quite encouraging. Aside from the cost of the MCSCF wave function there is no extra cost to compute the MC-PDFT energy besides evaluating the density and on-top pair density

FIG. 4 :

 4 FIG. 4: CASPT2(16,12)/ANO-L O[5s4p1d]/H[2s1p] potential energy curves describing the photodissociation of the water dimer with respect to the intramolecular O 1 -H 4 distance. 34Dashed lines represent the water donor localised electronic excited state energies whereas full lines refer to the ground state. The geometries of the water dimer at the Franck-Condon region (left-hand side) and at the dissociation limit (right-hand side) have also been given to provide a visual aid for the photo-process.

FIG. 5 :

 5 FIG. 5: Parallel efficiency for different calculations.

  )to eliminate the small component of the molecular 4-spinors, which would otherwise give rise to unwanted negative-energy (positronic) states. The electronic states are fully described by the two-component operator h + which can be further decomposed into scalar-relativistic and spin-orbit parts. Only large-component basis functions are required for relativistic two-component calculations. There is no need for the small-component basis, but it is implicitly involved in the decoupling transformation step. The decoupling transformation matrix is evaluated for the one-electron Hamiltonian of the Fock operator only and the two-electron contribution is averaged in the atomic mean-field operator. Therefore, the one-component relativistic calculation comprises the same two-electron integral evaluation and self-consistent-field procedure as in the non-relativistic case. The evaluation of the oneelectron relativistic Hamiltonians and any property requires the evaluation of one additional relativistic one-electron integrals matrix each. In addition, relativistic calculations of coreshell properties like contact densities may require more computational effort since large basis sets such as ANO-RCC 71,72 with many steep functions are required to correctly describe the behavior of the relativistic wave function close to the atomic nuclei.73,74

  which results into an order-byorder decoupling in terms of the external potential. It has been shown that the sequence converges to the reference four-component results 61,65 even for systems with very large nuclear charges. The polynomially-scaling arbitrary-order algorithm 65 of the DKH method has been implemented into the Molcas package 75 to replace the old exponentially-scaling one 61 . The expansion technique of the DKH method to approach exact decoupling can be replaced by the iterative technique which yields the BSS method. The so-called X2C method provides another way to achieve exact decoupling, where a preliminary diagonalization of a four-component one-electron matrix is used to obtain the exact decoupling transformation matrix. For a detailed comparison of these approaches, we refer the reader to the numerical examples in Ref. 75, which were obtained with the Molcas package.

  over all spin degrees of freedom. Further, N is the number of electrons and Ŝu = N i=1 Ŝu (i) is the component u of the one-electron spin vector operator. The function m u (r) is a component of the spin magnetization vector field which varies as a function of the electron position r. The reason for the factor of 2 in Equation (20b) is the following: Without SO coupling, and with the usual choice of u = z for the spin quantization axis, the function m z corresponds to the usual spin density ρ α-β , which integrates to 2M S = 2 S z for spin eigenfunctions. For the spin multiplet component with M S = S, two times M S counts the number of excess α spin electrons over β.

FIG. 6 :

 6 FIG. 6: NSOs ϕ z p of Equation (21b) corresponding to the non-bonding 5f orbitals of UO 2 + (5f 1 ), NpO 2 2+ (5f 1 ), and PuO 2 2+ (5f 2 ), and the corresponding contributions n z to 2 S z . Adapted from Ref. 79. A ground state doublet component with S z > 0 was chosen in each case. Isosurface values: ±0.03 atomic units. The isosurfaces shown here are for NpO 2 2+ . Those for the other two systems are visually indistinguishable from the ones shown.

FIG. 7 :

 7 FIG. 7: Selected NOs ϕ p of Equation (21a) of (C 5 (CH 3 ) 4 H) 3 UNO, and the corresponding occupation numbers. Adapted from Ref. 80. Isosurface values: ±0.03 atomic units. The isosurfaces shown here are for the SO calculation. The corresponding SR orbitals have very similar appearances.

FIG. 8 :

 8 FIG.8: The current single_aniso module.

FIG. 9 :

 9 FIG. 9: Main magnetic axes (X m , Y m , Z m , green arrows) and main anisotropy axes (X a , Y a , Z a , purple arrows) of the ground pseudospin S = 1 of the Ni 2+ (3d 8 configuration) with respect to the molecular coordinate system (x, y, z). The results are obtained in a CASSCF/CASPT2/RASSI-SO/single_aniso calculation using the ANO-RCC basis sets for all atoms with VTZP contraction for the Ni 2+ and first coordinated atoms, VDZ contraction for all other atoms, active space CAS(8,10) comprising 3d + 3d' shells of the Ni 2+ , 10 spin triplet and 15 spin singlet states optimized in state-average CASSCF calculations followed by standard Multistate-CASPT2, with subsequent mixing of all spin states by spin-orbit coupling in rassi. Color scheme: turquoise Ni 2+ ; red O; blue N; grey C; white H. Adapted from Ref. 85.

FIG. 10

 10 FIG.10: a) Molecular structure of the macrocyclic DyZn 3 compound.[START_REF] Feltham | [END_REF] Red dashed line shows the main magnetic axis in the ground Kramer's doublet. Color scheme: Dy, violet;Zn, yellow-green; O, red; N, blue; C, grey; H, white. These calculations have been performed using ANO-RCC basis sets of the following contractions: [8s7p5d4f2g1h] for Dy, [5s4p2d1f] for Zn, [4s3p1d] for O and N from the first coordination sphere, [3s2p] for C, and all remaining N and O, and [2s] for H; the active space of the CASSCF method included nine electrons spanning seven 4f orbitals, CAS (9 in 7); b) Ab initio calculated magnetism for 6 H, 6 F and 6 P terms mixed by the spin-orbit coupling. c) the same as (b) for all sextets, 128 quartets and 130 doublets originating from the 4f 9 electronic configuration of the Dy(III). Adapted from Ref. 93.

FIG. 12 :

 12 FIG. 12: Left: Orientation of the main magnetic axes of the ground and first excited doublet states (dashed line) with respect to the molecular structure in [Er(COT) 2 ] -. 99 Right: The magnetization blocking barrier in [Er(COT) 2 ] -anion. 99 The thick black lines represent the components of the low-lying Kramer's doublets as function of their magnetic moment along the axis connecting the centers of cyclooctatetraene (COT) rings. The green dashed lines correspond to diagonal quantum tunneling of magnetization (QTM); the blue dashed lines represent Orbach relaxation processes. The numbers at each arrow stand for the mean absolute value of the corresponding matrix element of transition magnetic moment ((|µ X | + |µ Y | + |µ Z |)/3). The path shown by the red arrows represents the most probable path for magnetic relaxation in this compound. Adapted from Ref. 99.

FIG. 14 :

 14 FIG.14: Energy profiles and electronic-state populations for the ultrafast dynamics of a Z -E photoisomerizing molecular motor involving four electronic singlet states (S 0 -grey-, S 1 -red-, S 2 -blue-, S 3 -green-and the simulation state -violet-).

FIG. 15 :

 15 FIG. 15: Comparison between CASSCF/6-31G*/AMBER and CASPT2/6-31G*/AMBER QM/MM nonadiabatic(Tully) trajectories describing the photoisomerization of the sensory rhodopsin from the cyanobacterium Anabaena PCC 7120. A. Reduced 13-cis-retinal chromophore model used in the calculation. The reduced chromophore model has been inserted in the protein cavity optimized for the full chromophore model (top). The cavity has then been re-defined as comprising all side-chains within 4.0 Å from any atom of the reduced chromophore while the rest of the protein has been kept fixed. B. Evolution of the dihedral angle describing the change in the double bond geometry (red part of the molecule). The set of ten CASSCF-based and CASPT2-based trajectories, selected by a Boltzmann-type equilibrium (to be published) demonstrate a qualitatively consistent behavior of the two levels of theory with a tendency of the CASPT2-based trajectories to decay at a later time and on a wider region.

Figure 19 a

 19 Figure 19 a shows the first example of simulated 2DES spectra obtained with the Molcas-Cobramm interface151,152 . Simulated 2D electronic spectra of a tetrapeptide with two aromatic (phenyl and phenolic) side chains contain enough information to distinguish between configurations with distant and vicinal side chains. SA-CASSCF(14,13)/SS-CASPT2 level of theory has been used to include all valence π-electrons and π-orbitals of both UV chromophores in the active space and 70 states were included in the state-average procedure. The proposed SOS//QM/MM protocol indicated the presence of important signals that are completely neglected using conventional Frenkel exciton Hamiltonian models, including charge-transfer (CT) states that characterize the coupling of UV-chromophores.

FIG. 19 :

 19 FIG. 19: Simulations of 2D electronic spectra obtained with the Molcas-Cobramm interface, including the 2DUV spectra of two different configurations of a tetrapeptide with two aromatic chromophores (panel a). and comparison between experimental and theoretical 2D spectra of pyrene at 1 ps upon photoexcitation (b).

  simulations. Here we present (a) the new state average CASSCF gradients in association with density fitting and subsequently two novel methods useful for (b) more efficient numerical gradient calculations and (c) for optimizations with fixed subfragments of a molecular system. A. Cholesky Decomposition-based Quantum Chemistry: from accurate energetics to structures of large molecules The Molcas package is mostly used for its array of efficient multiconfigurational wave function methods. Before the release of Molcas-7, the existing algorithms implementing these methods suffered from a major shortcoming arising from the calculation, transformation, and storage of electron repulsion integrals (ERIs). Such bottleneck made it difficult to use Molcas for larger systems and molecular dynamics studies, where many energy and gradient calculations are needed. This problem has since Molcas-7 been tackled and effectively solved with the development of techniques based on Cholesky decomposition (CD) of ERIs. The CD approximation to the ERIs had been introduced in quantum chemistry in 1977 by Beebe and Linderberg 161 but almost forgotten due to difficulties in the development of an efficient decomposition algorithm -with N a measure of system size, a naïve implementation would scale as N 6 as compared to the N 4 scaling of simply evaluating the ERIs. Using straightforward Cauchy-Schwartz screening, a cubic-scaling algorithm was presented in 2003 by Koch et al. 162 and further improvements have since been implemented in Molcas-7.

  FIG. 20: Comparison of the time of the different steps in CD-based RI-SA-CASSCF gradients for some Heme models.

  FIG. 23: 2,3-dimethyl-2-butenal

  FIG.25: Overview of the Molcas library interfaces. Optional components are surounded by a dashed line. The GA API must be provided by either the Global Arrays library or the internal DGA implementation. GPU support is possible through the use of NVBLAS, but an additional fully functional BLAS/LAPACK library is still needed.

FIG. 26 :

 26 FIG. 26: Schematically depicted data and program flow for mixed operation of Columbus and Molcas under control of the Molcas-driver. The depicted work flow refers to single-point energy evaluations and structure optimizations. Molcas and Columbus modules as well as the corresponding sections of the input file are depicted in blue and red, respectively.

FIG. 27 : 2 ) 2 )

 2722 FIG. 27: Potential energy curves for the copper-benzene adduct in terms of the perpendicular Cu-Bz distance R. The ground state (blue, Cu( 2 S)⊗Bz( 1 A 1 )) and the lowest excited states derived from Cu( 2 D)⊗Bz( 1 A 1 ) (green), Cu( 2 S)⊗Bz( 3 B 2 ) (black) and Cu( 2 D)⊗Bz( 3 B 2 )(red) are displayed. Solid lines indicate doublets, dashed lines quartets.Energies are given relative to the dissociation products of the ground state. Blue dots indicate the distances for which data have been computed. The inset collects the perturbational spin-orbit MR-CISD results for the three lowest dissociation channels Bz( 1 A 1 ) + Cu( 2 S) (blue), Bz( 1 A 1 )+Cu( 2D52 ) (green) and Bz(1A 1)+Cu( 2 D32 ) (light green).

Finally

  , some technical improvements of Molcas code (including new layers used for I/O and parallelization) open a possibility to run the code at larger varieties of modern hardware architectures. The link between Molcas and the multireference configuration interaction (MRCI) program system Columbus provides the possibility of closely coupled use of the features of both program packages. The construction of this link is quite flexible: a Molcas user can either use the link directly via the Molcas driver utility or via a

TABLE II :

 II The resulting transition energies are main CH 3 ∆E(S 0 → S 1 ) ∆E(S 1 → S 0 ) CASPT2 vertical transition energies of 2,3-dimethyl-2-butenal at different optimized geometries. The first two rows are conventional optimizations, the last two rows use composite gradients, with different methods for the internal methyl coordinates and the rest of the system.

	full CASSCF	3.88	2.50
	full CASPT2	3.84	2.56
	CASPT2 CASSCF	3.84	2.56
	CASPT2 MP2	3.84	2.56

TABLE III :

 III Time needed to generate one-and two-electron integrals (seward), variationally optimize multiconfigurational wave function (rasscf), compute Lagrangian multipliers for the wave function via linear response theory (mclr), and compute dynamic correlation energy by means of the second-order perturbation theory (caspt2). The Molcas performance are measured by using two types of persistent data storage devices, including: 1Tb Western Digital HDD, and 256Gb Intel SSD. The disk usage statistics, and the reference timings correspond to the run when HDD is used for keeping scratch data.The parallel functionality of Molcas has for a long time been exclusively based on a combination of the Global Arrays (GA) library 183 .for memory management and parallellization, and/or allowing the use of an external, threaded BLAS/LAPACK library for linear algebra. GA provides a convenient abstraction layer for operations on arrays that are distributed over different parallel processes. The library supports various one-sided operations for manipulating data in such a global array.The GA library is a large package containing lots of features that are not used by Molcas. Since Molcas fully relied on GA for its parallel functionality, most common problems with running Molcas in parallel were thus related to GA, an external project that Molcas had little to no control over. Considering the relatively small GA feature set used by Molcas, it was decided to develop our own GA-compatible layer implementing a subset of the

	Module	Disk Usage (Gb) HDD a (min) FiM b Time Relative speedup SSD c
	seward	24.3	11.7	1.6×	1.8×
	rasscf	28.5	23.3	1.0×	1.1×
	mclr	52.4	31.1	1.9×	2.2×
	caspt2	58.6	252.0	1.5×	2.7×
	In order to demonstrate the efficiency of FiM, we will provide representative timings
	on a typical Molcas multiconfigurational calculation. The performance of FiM has been
	measured by doing a two-state averaged CASSCF(12,12), mclr, and CASPT2 calculations
	on the 54-atom rhodopsin chromophore molecule (C 21 H 32 N) employing the 6-31G* basis set
	(394 basis functions). As one can see from Table III, the FiM I/O layer performs reasonably

a The WD10EURX-73FH1Y0 model b The FiM timings were obtained by using 28Gb out of 32Gb RAM available for emulating I/O. c The SSDSC2CW240A3 model and I/O access patterns. This I/O interference problem has been mitigated in Molcas-8 by developing a new framework called "Files in Memory" (FiM). The key idea of FiM is to keep a scratch file in RAM entirely instead of using a persistent storage device. Within the FiM layer, the Read/Write operations are just a user-space memory copying with practically instant (nanoseconds) seek time. FiM can be easily activated by modifying certain Molcas configuration files (details can be found in the manual). well compared to SSD. Indeed, within FiM we obtained a speed-up factor of up to 1.9 times in comparison with corresponding HDD disk-based timings. Another feature is that FiM allows users to run Molcas-8 on diskless compute nodes without incurring significant overall performance penalty. B. Parallel framework in Molcas-8 GA API, simply called distributed global arrays (DGA) 184 . While GA supports lots of different communication back-ends, DGA relies solely on MPI for its low-level communication layer.

ACKNOWLEDGMENTS

CEH, LG, and DGT would like to thank the Scientific Discovery through Advanced Computing (SciDAC) program funded by U.S. Department of Energy, Office of Science, Basic Energy Sciences and Advanced Scientific Computing Research under Award No. DE-

Although nowadays a commodity workstation is equipped with a relatively large amount of memory, neither the Molcas program itself nor a buffer cache mechanism of an underlying operating system (OS kernel) can always utilize this memory in an optimal way. On the one hand, the original Molcas code was developed in the beginning of the 90s and it was solely designed to minimize RAM demands at the cost of an extra I/O. On the other hand, the OS kernel caches data from all opened files uniformly, regardless of their sizes

Multi-State Calculations Including Spin-Orbit Interaction on the Benzene-

Copper Complex The motivation for this example is derived from the aim to perform ab initio modeling of interactions of transition metals with graphene nanoribbons. The adsorption process is important for spintronics applications 192 as it transforms the semimetallic graphene to a spin Hall insulator via spin-orbit coupling. 193,194 As starting point for these investigations we choose the adsorption complex of Cu with benzene. Even this relatively small system is quite demanding as the two lowest electronic states of Cu ( 2 S and 2 D), as well as the singlet ground state ( 1 A 1 ) and the first triplet state ( 3 B 2 ) of benzene will be considered. In the scalar relativistic calculations for the two combined subsystems 6 quartets and 12 doublets have been included while in the SO-MRCI approach 24 states were calculated retaining solely one component of each Kramer's pair. The construction of the wave functions is based on a direct product space CAS(11e,6o)⊗CAS(6e,6o) consisting of a CAS for 11 electrons and 6 orbitals referring to the Cu atom which includes the 3d and 4s

shells. The second space relates to the benzene molecule (Bz). More computational details referring to the state-averaging scheme, the construction of MRCI reference functions and the basis set used can be found in the Suppl. Mat. The scalar relativistic calculations have been performed by means of the Columbus/Molcas link using the AO integrals from Molcas integral code seward. For the SO-MRCISD calculations, however, the atomic orbital (AO) integrals (including SO integrals) were calculated using the program Argos developed by

Pitzer 195,196 , since seward does not support spin-orbit relativistic effective core potentials (RECP). Thus, and also due to the direct product type reference space, calculations have been executed under control of the Columbus driver. Otherwise, scalar-relativistic or nonrelativistic as well as variational spin-orbit MR-CISD calculations including appropriate SS/SA-MCSCF calculations may be specified and executed within a Molcas environment (cf. Suppl. Mat. for the actual input).

In