Modern quantum chemistry with [Open]Molcas.

MOLCAS/OpenMolcas is an ab initio electronic structure program providing a large set of computational methods from Hartree-Fock and density functional theory to various implementations of multiconfigurational theory. This article provides a comprehensive overview of the main features of the code, specifically reviewing the use of the code in previously reported chemical applications as well as more recent applications including the calculation of magnetic properties from optimized density matrix renormalization group wave functions.

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, two surface hopping algorithms are included to enable nonadiabatic calculations, and the DQ method for diabatization is added. Finally, we report on the subject of improvements with respects to alternative file options and parallelization.

First time combination of frozen density embedding theory with the algebraic diagrammatic construction scheme for the polarization propagator of second order.

The combination of Frozen Density Embedding Theory (FDET) and the Algebraic Diagrammatic Construction (ADC) scheme for the polarization propagator for describing environmental effects on electronically excited states is presented. Two different ways of interfacing and expressing the so-called embedding operator are introduced. The resulting excited states are compared with supermolecular calculations of the total system at the ADC(2) level of theory. Molecular test systems were chosen to investigate molecule-environment interactions of varying strength from dispersion interaction up to multiple hydrogen bonds. The overall difference between the supermolecular and the FDE-ADC calculations in excitation energies is lower than 0.09 eV (max) and 0.032 eV in average, which is well below the intrinsic error of the ADC(2) method itself.

Orbital entanglement and CASSCF analysis of the Ru-NO bond in a Ruthenium nitrosyl complex.

Complete active space self-consistent field (CASSCF) wavefunctions and an orbital entanglement analysis obtained from a density-matrix renormalisation group (DMRG) calculation are used to understand the electronic structure, and, in particular, the Ru-NO bond of a Ru nitrosyl complex. Based on the configurations and orbital occupation numbers obtained for the CASSCF wavefunction and on the orbital entropy measurements evaluated for the DMRG wavefunction, we unravel electron correlation effects in the Ru coordination sphere of the complex. It is shown that Ru-NO π bonds show static and dynamic correlation, while other Ru-ligand bonds feature predominantly dynamic correlation. The presence of static correlation requires the use of multiconfigurational methods to describe the Ru-NO bond. Subsequently, the CASSCF wavefunction is analysed in terms of configuration state functions based on localised orbitals. The analysis of the wavefunctions in the electronic singlet ground state and the first triplet state provides a picture of the Ru-NO moiety beyond the standard representation based on formal oxidation states. A distinct description of the Ru and NO fragments is advocated. The electron configuration of Ru is an equally weighted superposition of Ru(II) and Ru(III) configurations, with the Ru(III) configuration originating from charge donation mostly from Cl ligands. However, and contrary to what is typically assumed, the electronic configuration of the NO ligand is best described as electroneutral.

Orthogonality of embedded wave functions for different states in frozen-density embedding theory.

Other than lowest-energy stationary embedded wave functions obtained in Frozen-Density Embedding Theory (FDET) [T. A. Wesolowski, Phys. Rev. A 77, 012504 (2008)] can be associated with electronic excited states but they can be mutually non-orthogonal. Although this does not violate any physical principles--embedded wave functions are only auxiliary objects used to obtain stationary densities--working with orthogonal functions has many practical advantages. In the present work, we show numerically that excitation energies obtained using conventional FDET calculations (allowing for non-orthogonality) can be obtained using embedded wave functions which are strictly orthogonal. The used method preserves the mathematical structure of FDET and self-consistency between energy, embedded wave function, and the embedding potential (they are connected through the Euler-Lagrange equations). The orthogonality is built-in through the linearization in the embedded density of the relevant components of the total energy functional. Moreover, we show formally that the differences between the expectation values of the embedded Hamiltonian are equal to the excitation energies, which is the exact result within linearized FDET. Linearized FDET is shown to be a robust approximation for a large class of reference densities.

Analytical gradients of the state-average complete active space self-consistent field method with density fitting.

An efficient implementation of the state-averaged complete active space self-consistent field (SA-CASSCF) gradients employing density fitting (DF) is presented. The DF allows a reduction both in scaling and prefactors of the different steps involved. The performance of the algorithm is demonstrated on a set of molecules ranging up to an iron-Heme b complex which with its 79 atoms and 811 basis functions is to our knowledge the largest SA-CASSCF gradient computed. For smaller systems where the conventional code could still be used as a reference, both the linear response calculation and the gradient formation showed a clear timing reduction and the overall cost of a geometry optimization is typically reduced by more than one order of magnitude while the accuracy loss is negligible.

A two-scale approach to electron correlation in multiconfigurational perturbation theory.

We present a new approach for the calculation of dynamic electron correlation effects in large molecular systems using multiconfigurational second-order perturbation theory (CASPT2). The method is restricted to cases where partitioning of the molecular system into an active site and an environment is meaningful. Only dynamic correlation effects derived from orbitals extending over the active site are included at the CASPT2 level of theory, whereas the correlation effects of the environment are retrieved at lower computational costs. For sufficiently large systems, the small errors introduced by this approximation are contrasted by the substantial savings in both storage and computational demands compared to the full CASPT2 calculation. Provided that static correlation effects are correctly taken into account for the whole system, the proposed scheme represent a hierarchical approach to the electron correlation problem, where two molecular scales are treated each by means of the most suitable level of theory.

Analytical gradients of complete active space self-consistent field energies using Cholesky decomposition: geometry optimization and spin-state energetics of a ruthenium nitrosyl complex.

We present a formulation of analytical energy gradients at the complete active space self-consistent field (CASSCF) level of theory employing density fitting (DF) techniques to enable efficient geometry optimizations of large systems. As an example, the ground and lowest triplet state geometries of a ruthenium nitrosyl complex are computed at the DF-CASSCF level of theory and compared with structures obtained from density functional theory (DFT) using the B3LYP, BP86, and M06L functionals. The average deviation of all bond lengths compared to the crystal structure is 0.042 Å at the DF-CASSCF level of theory, which is slightly larger but still comparable with the deviations obtained by the tested DFT functionals, e.g., 0.032 Å with M06L. Specifically, the root-mean-square deviation between the DF-CASSCF and best DFT coordinates, delivered by BP86, is only 0.08 Å for S0 and 0.11 Å for T1, indicating that the geometries are very similar. While keeping the mean energy gradient errors below 0.25%, the DF technique results in a 13-fold speedup compared to the conventional CASSCF geometry optimization algorithm. Additionally, we assess the singlet-triplet energy vertical and adiabatic differences with multiconfigurational second-order perturbation theory (CASPT2) using the DF-CASSCF and DFT optimized geometries. It is found that the vertical CASPT2 energies are relatively similar regardless of the geometry employed whereas the adiabatic singlet-triplet gaps are more sensitive to the chosen triplet geometry.

Attractive electron-electron interactions within robust local fitting approximations.

An analysis of Dunlap's robust fitting approach reveals that the resulting two-electron integral matrix is not manifestly positive semidefinite when local fitting domains or non-Coulomb fitting metrics are used. We present a highly local approximate method for evaluating four-center two-electron integrals based on the resolution-of-the-identity (RI) approximation and apply it to the construction of the Coulomb and exchange contributions to the Fock matrix. In this pair-atomic resolution-of-the-identity (PARI) approach, atomic-orbital (AO) products are expanded in auxiliary functions centered on the two atoms associated with each product. Numerical tests indicate that in 1% or less of all Hartree-Fock and Kohn-Sham calculations, the indefinite integral matrix causes nonconvergence in the self-consistent-field iterations. In these cases, the two-electron contribution to the total energy becomes negative, meaning that the electronic interaction is effectively attractive, and the total energy is dramatically lower than that obtained with exact integrals. In the vast majority of our test cases, however, the indefiniteness does not interfere with convergence. The total energy accuracy is comparable to that of the standard Coulomb-metric RI method. The speed-up compared with conventional algorithms is similar to the RI method for Coulomb contributions; exchange contributions are accelerated by a factor of up to eight with a triple-zeta quality basis set. A positive semidefinite integral matrix is recovered within PARI by introducing local auxiliary basis functions spanning the full AO product space, as may be achieved by using Cholesky-decomposition techniques. Local completion, however, slows down the algorithm to a level comparable with or below conventional calculations.

A Multiconfigurational Wave Function Implementation of the Frenkel Exciton Model for Molecular Aggregates.

We present an implementation of the Frenkel exciton model into the OpenMolcas program package enabling calculations of collective electronic excited states of molecular aggregates based on a multiconfigurational wave function description of the individual monomers. The computational protocol avoids using diabatization schemes and, thus, supermolecule calculations. Additionally, the use of the Cholesky decomposition of the two-electron integrals entering pair interactions enhances the efficiency of the computational scheme. The application of the method is exemplified for two test systems, that is, a formaldehyde oxime and a bacteriochlorophyll-like dimer. For the sake of comparison with the dipole approximation, we restrict our considerations to situations where intermonomer exchange can be neglected. The protocol is expected to be beneficial for aggregates composed of molecules with extended π systems, unpaired electrons such as radicals or transition metal centers, where it should outperform widely used methods based on time-dependent density functional theory.

Characterizing Conical Intersections in DNA/RNA Nucleobases with Multiconfigurational Wave Functions of Varying Active Space Size.

We characterize the photochemically relevant conical intersections between the lowest-lying accessible electronic excited states of the different DNA/RNA nucleobases using Cholesky decomposition-based complete active space self-consistent field (CASSCF) algorithms. We benchmark two different basis set contractions and several active spaces for each nucleobase and conical intersection type, measuring for the first time how active space size affects conical intersection topographies in these systems and the potential implications these may have toward their description of photoinduced phenomena. Our results show that conical intersection topographies are highly sensitive to the electron correlation included in the model: by changing the amount (and type) of correlated orbitals, conical intersection topographies vastly change, and the changes observed do not follow any converging pattern toward the topographies obtained with the largest and most correlated active spaces. Comparison across systems shows analogous topographies for almost all intersections mediating population transfer to the dark 1nO/Nπ* states, while no similarities are observed for the "ethylene-like" conical intersection ascribed to mediate the ultrafast decay component to the ground state in all DNA/RNA nucleobases. Basis set size seems to have a minor effect, appearing to be relevant only for purine-based derivatives. We rule out structural changes as a key factor in classifying the different conical intersections, which display almost identical geometries across active space and basis set change, and we highlight instead the importance of correctly describing the electronic states involved at these crossing points. Our work shows that careful active space selection is essential to accurately describe conical intersection topographies and therefore to adequately account for their active role in molecular photochemistry.

The OpenMolcas Web: A Community-Driven Approach to Advancing Computational Chemistry.

The developments of the open-source OpenMolcas chemistry software environment since spring 2020 are described, with a focus on novel functionalities accessible in the stable branch of the package or via interfaces with other packages. These developments span a wide range of topics in computational chemistry and are presented in thematic sections: electronic structure theory, electronic spectroscopy simulations, analytic gradients and molecular structure optimizations, ab initio molecular dynamics, and other new features. This report offers an overview of the chemical phenomena and processes OpenMolcas can address, while showing that OpenMolcas is an attractive platform for state-of-the-art atomistic computer simulations.

Self-consistency in frozen-density embedding theory based calculations.

The bi-functional for the non-electrostatic part of the exact embedding potential of frozen-density embedding theory (FDET) depends on whether the embedded part is described by means of a real interacting many-electron system or the reference system of non-interacting electrons (see [Wesolowski, Phys. Rev. A. 77, 11444 (2008)]). The difference δΔF(MD)[ρ(A)]/δρ(A)(r), where ΔF(MD)[ρ(A)] is the functional bound from below by the correlation functional E(c)[ρ(A)] and from above by zero. Taking into account ΔF(MD)[ρ(A)] in both the embedding potential and in energy is indispensable for assuring that all calculated quantities are self-consistent and that FDET leads to the exact energy and density in the limit of exact functionals. Since not much is known about good approximations for ΔF(MD)[ρ(A)], we examine numerically the adequacy of neglecting ΔF(MD)[ρ(A)] entirely. To this end, we analyze the significance of δΔF(MD)[ρ(A)]/δρ(A)(r) in the case where the magnitude of ΔF(MD)[ρ(A)] is the largest, i.e., for Hartree-Fock wavefunction. In hydrogen bonded model systems, neglecting δΔF(MD)[ρ(A)]/δρ(A)(r) in the embedding potential marginally affects the total energy (less than 5% change in the interaction energy) but results in qualitative changes in the calculated hydrogen-bonding induced shifts of the orbital energies. Based on this estimation, we conclude that neglecting δΔF(MD)[ρ(A)]/δρ(A)(r) may represent a good approximation for multi-reference variational methods using adequate choice for the active space. Doing the same for single-reference perturbative methods is not recommended. Not only it leads to violation of self-consistency but might result in large effect on orbital energies. It is shown also that the errors in total energy due to neglecting δΔF(MD)[ρ(A)]/δρ(A)(r) do not cancel but rather add up to the errors due to approximation for the bi-functional of the non-additive kinetic potential.

Strong correlation treated via effective hamiltonians and perturbation theory.

We propose a new approach to determine a suitable zeroth-order wavefunction for multiconfigurational perturbation theory. The same ansatz as in complete active space (CAS) wavefunction optimization is used but it is split in two parts, a principal space (A) and a much larger extended space (B). Löwdin's partitioning technique is employed to map the initial eigenvalue problem to a dimensionality equal to that of (A) only. Combined with a simplified expression for the (B) portion of the wavefunction, we are able to drastically reduce the storage and computational demands of the wavefunction optimization. This scheme is used to produce reference wavefunctions and energies for subsequent second-order perturbation theory (PT2) corrections. Releasing the constraint of computing the exact CAS energy and wavefunction prior to the PT2 treatment introduces a nonstandard paradigm for multiconfigurational methods. Based on the results of test calculations, we argue that principal parts with only few percents of the total number of CAS configurations could provide final multiconfigurational PT2 energies of the same accuracy as in the standard paradigm. In the future, algorithmic improvements for this scheme will bring into reach active spaces much beyond the present limit of CAS-based methods, therefore allowing for accurate studies of systems featuring strong correlation.

MOLCAS 7: the next generation.

Some of the new unique features of the MOLCAS quantum chemistry package version 7 are presented in this report. In particular, the Cholesky decomposition method applied to some quantum chemical methods is described. This approach is used both in the context of a straight forward approximation of the two-electron integrals and in the generation of so-called auxiliary basis sets. The article describes how the method is implemented for most known wave functions models: self-consistent field, density functional theory, 2nd order perturbation theory, complete-active space self-consistent field multiconfigurational reference 2nd order perturbation theory, and coupled-cluster methods. The report further elaborates on the implementation of a restricted-active space self-consistent field reference function in conjunction with 2nd order perturbation theory. The average atomic natural orbital basis for relativistic calculations, covering the whole periodic table, are described and associated unique properties are demonstrated. Furthermore, the use of the arbitrary order Douglas-Kroll-Hess transformation for one-component relativistic calculations and its implementation are discussed. This section especially focuses on the implementation of the so-called picture-change-free atomic orbital property integrals. Moreover, the ElectroStatic Potential Fitted scheme, a version of a quantum mechanics/molecular mechanics hybrid method implemented in MOLCAS, is described and discussed. Finally, the report discusses the use of the MOLCAS package for advanced studies of photo chemical phenomena and the usefulness of the algorithms for constrained geometry optimization in MOLCAS in association with such studies.

Generating Cu(II)-oxyl/Cu(III)-oxo species from Cu(I)-alpha-ketocarboxylate complexes and O2: in silico studies on ligand effects and C-H-activation reactivity.

A mechanism for the oxygenation of Cu(I) complexes with alpha-ketocarboxylate ligands that is based on a combination of density functional theory and multireference second-order perturbation theory (CASSCF/CASPT2) calculations is elaborated. The reaction proceeds in a manner largely analogous to those of similar Fe(II)-alpha-ketocarboxylate systems, that is, by initial attack of a coordinated oxygen molecule on a ketocarboxylate ligand with concomitant decarboxylation. Subsequently, two reactive intermediates may be generated, a Cu-peracid structure and a [CuO](+) species, both of which are capable of oxidizing a phenyl ring component of the supporting ligand. Hydroxylation by the [CuO](+) species is predicted to proceed with a smaller activation free energy. The effects of electronic and steric variations on the oxygenation mechanisms were studied by introducing substituents at several positions of the ligand backbone and by investigating various N-donor ligands. In general, more electron donation by the N-donor ligand leads to increased stabilization of the more Cu(II)/Cu(III)-like intermediates (oxygen adducts and [CuO](+) species) relative to the more Cu(I)-like peracid intermediate. For all ligands investigated, the [CuO](+) intermediates are best described as Cu(II)-O(*-) species with triplet ground states. The reactivity of these compounds in C-H abstraction reactions decreases with more electron-donating N-donor ligands, which also increase the Cu-O bond strength, although the Cu-O bond is generally predicted to be rather weak (with a bond order of about 0.5). A comparison of several methods to obtain singlet energies for the reaction intermediates indicates that multireference second-order perturbation theory is likely more accurate for the initial oxygen adducts, but not necessarily for subsequent reaction intermediates.

Reaction mechanism of manganese superoxide dismutase studied by combined quantum and molecular mechanical calculations and multiconfigurational methods.

Manganese superoxide dismutases (MnSODs) are enzymes that convert two molecules of the poisonous superoxide radical into molecular oxygen and hydrogen peroxide. During the reaction, the manganese ion cycles between the Mn(2+) and Mn(3+) oxidation states and accomplishes its enzymatic action in two half-cycles (corresponding to the oxidation and reduction of O(2)(-)). Despite many experimental and theoretical studies dealing with SODs, including quantum chemical active-site-model studies of numerous variants of the reaction mechanisms, several details of MnSOD enzymatic action are still unclear. In this study, we have modeled and compared four reaction pathways (one associative, one dissociative, and two second-sphere) in a protein environment using the QM/MM approach (combined quantum and molecular mechanics calculations) at the density functional theory level. The results were complemented by CASSCF/CASPT2/MM single-point energy calculations for the most plausible models to account properly for the multireference character of the various spin multiplets. The results indicate that the oxidation of O(2)(-) to O(2) most likely occurs by an associative mechanism following a two-state (quartet-octet) reaction profile. The barrier height is estimated to be less than 25 kJ.mol(-1). On the other hand, the conversion of O(2)(-) to H(2)O(2) is likely to take place by a second-sphere mechanism, that is, without direct coordination of the superoxide radical to the manganese center. The reaction pathway involves the conical intersection of two quintet states, giving rise to an activation barrier of approximately 60 kJ.mol(-1). The calculations also indicate that the associative mechanism can represent a competitive pathway in the second half-reaction with the overall activation barrier being only slightly higher than the activation barrier in the second-sphere mechanism. The activation barriers along the proposed reaction pathways are in very good agreement with the experimentally observed reaction rates of SODs (k(cat) approximately 10(4)-10(5) s(-1)).

Calculation of protein-ligand interaction energies by a fragmentation approach combining high-level quantum chemistry with classical many-body effects.

We have developed a method to estimate accurate interaction energies between a full protein and a bound ligand. It is based on the recently proposed PMISP (polarizable multipole interaction with supermolecular pairs) method (Soderhjelm, P.; Ryde, U. J. Phys. Chem. A 2009, 113, 617), which treats electrostatic interaction by multipoles up to quadrupoles, induction by anisotropic polarizabilities, and nonclassical interactions by explicit quantum mechanical (QM) calculations, using a fragmentation approach. For a whole protein, electrostatics and induction are treated the same way, but for the nonclassical interactions, a Lennard-Jones term from a standard molecular mechanics (MM) force field (e.g., Amber) is used outside a certain distance from the ligand (4-7 A). This QM/MM variant of the PMISP method is carefully tested by varying this distance. Several approximations related to the classical interactions are also evaluated. It is found that one can speed up the calculation by using density functional theory to compute multipoles and polarizabilities but that a proper treatment of polarization is important. As a demonstration of the method, the interaction energies of two ligands bound to avidin are calculated at the MP2/aug-cc-pVTZ level, with an expected relative error of 1-2%.

Atomic Cholesky decompositions: a route to unbiased auxiliary basis sets for density fitting approximation with tunable accuracy and efficiency.

Cholesky decomposition of the atomic two-electron integral matrix has recently been proposed as a procedure for automated generation of auxiliary basis sets for the density fitting approximation [F. Aquilante et al., J. Chem. Phys. 127, 114107 (2007)]. In order to increase computational performance while maintaining accuracy, we propose here to reduce the number of primitive Gaussian functions of the contracted auxiliary basis functions by means of a second Cholesky decomposition. Test calculations show that this procedure is most beneficial in conjunction with highly contracted atomic orbital basis sets such as atomic natural orbitals, and that the error resulting from the second decomposition is negligible. We also demonstrate theoretically as well as computationally that the locality of the fitting coefficients can be controlled by means of the decomposition threshold even with the long-ranged Coulomb metric. Cholesky decomposition-based auxiliary basis sets are thus ideally suited for local density fitting approximations.