Development and Application of a Complete Active Space Spin-Orbit Configuration Interaction Program Designed for Single Molecule Magnets.

We present a spin-orbit configuration interaction program which has been tailored for the description of the magnetic properties of polynuclear metal complexes with partially filled d- and f-shells. The spin-orbit operators are directly included in the configuration interaction program based on Slater-determinants. The lowest states are obtained by a Block-Davidson-type diagonalisation. The usage of localised active orbitals enables the construction of start vectors from tensor products of single-center wave functions that already include spin-orbit interaction. This allows for an analysis of the role and the interplay of the different metal centres. Furthermore, in case of weak coupling of the metal centres these tensor products are already close to the final wave functions ensuring fast convergence. In combination with a two-layer hybrid parallelisation, this makes the program highly efficient. Based on the spin-orbit coupled wave functions, magnetic D-tensors, g-tensors and temperature-dependent susceptibilities can be calculated. The applicability and performance of the program is shown exemplarily on a trinuclear transition metal (CoII VII CoII ) complex.

Introduction

Aggregates of several exchange‐coupled open‐shell transition metal spin centres show a large variety of magnetic properties, that depend on the attributes of the individual spin centres and how they interact. Depending on their magnetic behaviour, such aggregates can be used for magnetic cooling, information storage, or quantum computing. Quantum chemical calculations can be very useful to understand and improve the property one wants to achieve. In the field of information storage (single molecule magnets), for example, it is important that the molecule keeps its magnetisation as long as possible.

The magnetic properties of these systems are usually aroused by metal centres with partially filled d‐ or f‐shells which show spin and orbital degeneracies. Often, several magnetic centres are combined in one molecule and the unpaired electrons of these centres are coupled to almost degenerate spin states (exchange coupling). Because of the complexity of the electronic structure, wave function based multi‐configuration methods are required which allow to consider ligand‐field splittings, spin‐orbit coupling, and exchange coupling with the same accuracy. This is especially true if there are spatial near‐degeneracies within the individual spin centres. Several approaches to describe single molecule magnets are based on complete active space self consistent field (CASSCF) reference wave functions. The complete active space (CAS) wave functions, which include static correlation, can qualitatively describe the ligand‐field splitting and the quasi degenerate spin states.

CASSCF has severe deficiencies in describing exchange couplings in multinuclear transition metal complexes. Antiferromagnetic couplings in particular, commonly denoted as super‐exchange, are severely underestimated in many cases, often by a factor of 2–3.[ , , ] This is especially true if only a minimal set of metal‐centred active orbitals (5 active orbitals per d‐block centre, and 7 active orbitals for each f‐block centre) is used. Extending such an active space is not an option for multinuclear complexes since then one rather quickly hits computational limits. To improve upon CASSCF, it is usually considered mandatory to include dynamic correlation, either by perturbation theory or approaches of configuration interaction (CI) type. However, including dynamic correlation on top of a large active space CASSCF calculation is generally very costly. For example in Ref. [9], for the magnetic exchange coupling in the well‐studied copper acetate monohydrate molecule, the computational demanding difference dedicated CI (DDCI3) method had to be used to get reasonable agreement to experiment, while multireference methods based on second order perturbation theory failed to describe the ground‐state magnetic coupling. Moreover, exchange coupling constants with the density matrix renormalization group method can only be obtained at a high computational effort. Thus, such methods are not applicable if one aims for the coupling of more spin centers. Therefore we follow a pragmatic and cost‐effective strategy to improve CASSCF exchange couplings: It has been demonstrated that one of the main reasons for the underestimation of super‐exchange in CASSCF‐type calculations is the inflexibility of the CAS ground state wave functions to describe charge‐transfer (CT) configurations. By correcting the weight of the CT configurations in the wave functions, the modified complete active space configuration interaction methods (MCASCI) (cf. Section 2.4) can be used to describe large exchange coupled systems which are not accessible with traditional multi reference methods.

The inclusion of spin‐orbit coupling is compulsory in many applications. For example, g‐matrix anisotropy or zero‐field splitting would not exist (or be very small) without spin‐orbit coupling. In systems where the ground state is quasi‐degenerate in the absence of spin‐orbit coupling (that is, there are several spin multiplets with an energy spacing that is not much larger than the strength of the spin‐orbit interaction) these close‐lying states are strongly mixed when spin‐orbit coupling is switched on, such that one has to consider the whole quasi‐degenerate manifold in order to get a qualitatively correct description of the electronic and magnetic properties.

However, if one includes spin‐orbit coupling in the Hamiltonian, its symmetry is reduced and the calculations get much more involved. For this reason, in the common spin‐orbit configuration interaction ansatz (SOCI), one starts the calculation without spin‐orbit coupling and includes this in a second stage. To this end, one splits the Hamiltonian into a non‐relativistic (or scalar‐relativistic) part and a spin‐orbit part , and starts from CI‐type wave functions that diagonalise :

This method is very well suited to treat quasi‐degeneracies provided that the wave functions Ψ encompass at least the whole quasi‐degenerate set of eigenstates of the scalar Hamiltonian. This is the main drawback of the method because this number can become very large for oligonuclear transition metal complexes (it scales exponentially with the number of spin centres).

As an illustrative example, the spin‐orbit coupling in mono‐ and polynuclear high‐spin CoII compounds will be discussed. Since the ligand‐field states are most important for magnetic properties, is usually diagonalised with an active space containing seven electrons in the five 3d‐orbitals (7‐in‐5), i. e., with 120 Slater determinants if the blocking of the different m states is not taken into account. Depending on the ligand field, the ground state of for a single high‐spin CoII centre can be spatially degenerate (octahedral field) or non‐degenerate (tetrahedral field). In the octahedral case, there is a configuration with a ground state and the main part of the spin‐orbit coupling can be accessed in a SOCI where only the lowest three quartet states of are considered. The SOCI ground state is then an isolated Kramers doublet. Note, however, that this treatment lacks higher‐order spin‐orbit effects from the other ligand‐field states.

In the tetrahedral case we have a configuration with a ground state that is well separated from all other ligand‐field states. Spin‐orbit coupling affects ground‐state properties only in second order and it is necessary to include all ligand‐field states (10 quartets and 40 doublets) in the SOCI if one wants to calculate the zero‐field splitting of the quartet state. Now in case of a compound containing three of these CoII centres, an active space including the electrons and d‐orbitals of all three centres (21‐in‐15) must be used in the scalar relativistic calculation. This leads to large CI dimensions but this can be handled by modern full CI codes. The problem with SOCI is different: In the octahedral case, the three quartets at each CoII site couple to 27 spin ladders with one decet, two octet, three sextet, four quartet, and two doublet states each. Spin‐orbit coupling of all these 324 scalar‐relativistic states (encompassing 1728 micro‐states) yields four Kramers doublets at low energies. It is quite an effort (but still possible) to generate the CI vectors of all these states that are then used to calculate the SOCI matrix elements. In the tetrahedral case, however, one must include all ligand‐field states. In the trinuclear case, there are 564.000 such states (including 202.000 doublets, 188.000 quartets and 75.000 sextets), and it is no longer possible to calculate and store the CI vectors of all these states which is a prerequisite to perform the SOCI calculation.

An alternative approach, followed in our work and denoted CASOCI below, is to include directly in the Hamiltonian, , diagonalising both operators at the same time and the same space as proposed by Ganyushin and Neese in Ref. [14]. The CI dimension is larger than in the scalar case because spin‐orbit coupling mixes several m components and the CI coefficients are no longer real, but this does not impose a significant burden. What one gains is that one has to construct only a relatively small number of CI vectors if one is only interested in the low‐lying manifold, because the spin‐orbit effects are automatically included and one never has to construct high‐lying ligand‐field states as in the SOCI case. In our trinuclear CoII example with tetrahedral coordination, one needs to converge 64 roots (micro‐states) arising from the exchange coupling of the three local states, and in the octahedral case one only needs 8 micro‐states assuming that spin‐orbit coupling in the local states is stronger than exchange coupling. We use a determinant‐based CI approach and calculate the required roots efficiently via a Block‐Davidson‐type diagonalisation in a spin‐string basis using an effective one‐electron spin‐orbit operator. This approach has been tailored in order to calculate the magnetic properties of molecules comprising up to three exchange‐coupled open‐shell metal ions. It can be combined with the MCASCI approach in order to treat exchange coupled systems. If the exchange coupling is relatively weak (which is usually the case), rather good trial wave‐functions (Davidson start vectors) for the full system can be obtained from tensor products of single‐centre SOCI wave functions that are quite easily calculated through a full diagonalisation of the single‐centre matrix representations of the (spin‐orbit containing) Hamiltonian. In the next section, details of the method and the implementation are described. This includes also the extraction of zero‐field splitting tensors and magnetic properties from the CASOCI wave functions. Finally, a demonstrative application of the new program on a trinuclear CoIIVIICoII complex is presented.

Theory

Complete Active Space Spin‐Orbit Configuration Interaction for Magnetic Properties

In order to obtain thermodynamic data such as (temperature‐dependent) magnetic susceptibilities, and get information about possible electronic relaxation pathways, several of the lowest eigenstates of the molecular system of interest have to be calculated. In this work, the wave functions and energies of these electronic states are obtained in a direct CI calculation by applying the well‐established (multi‐root) Davidson‐Liu method.[ , ] Spin‐orbit coupling (SOC) is of paramount importance for many magnetic properties and has to be included in the Hamiltonian. For the calculation of magnetic susceptibilities one also needs to include the Zeeman interaction with an external magnetic field. We use a valence‐space Hamiltonian in the space of the active orbitals, where the effect of the doubly occupied inactive electrons are folded in through the core energy and Fock operators:

where denotes the usual scalar electronic Born‐Oppenheimer Hamiltonian, which reads

where refer to active orbital indices. includes the interaction between the core and valence electrons, and contains, besides the nuclear repulsion energy , the interaction between the core electrons:

(h and g are the matrix elements of the one‐particle Hamiltonian and the electron repulsion integrals). The indices i,j run over the inactive (doubly occupied in all Slater determinants) orbitals. The spin‐orbit and Zeeman‐terms require to go beyond the spin‐free single‐replacement operators and are thus defined entirely through the spin‐dependent single‐replacement operators (where refer either to α or β spin):

with the external (homogeneous) magnetic field vector. Note that we only consider (effective) single‐particle spin‐orbit contributions through the vector‐type matrix elements , which however contain the effect of the two‐electron spin‐orbit interaction in an averaged manner (cf. Section 2.5. for details). The Zeeman operator in its most general form is determined (for each cartesian component of the magnetic field) by the scalar ( ) and vector‐type ( ) matrix elements. Note that both can have the values x,y,z, but for λ these are cartesian components in spin space while for κ these are in the real space. In the nonrelativistic limit the Zeeman part takes a simpler form, since spin and orbital dependent parts separate

where denotes Bohr's magneton and is the free‐electron g‐factor. As it stands, the orbital part of the Zeeman operator depends on the origin chosen for the angular momentum operator. This so‐called gauge origin problem is a major obstacle if the response of a wave function to the magnetic field has to be calculated, but less so when calculating first‐order magnetic properties from wave functions obtained at zero field: Expectation values of the Zeeman operator with unperturbed wave functions, or matrix elements involved in (quasi‐) degenerate first‐order perturbation approaches, show a rather weak gauge dependence. For example, the magnetisation curve of the CoIIVIICoII complex (see Section 3.1) has been calculated using the central V atom or either of the two Co atoms as the gauge origin, and the plots of all three curves superimpose.

The last contribution to the Hamiltonian, , implements a shift of diagonal elements which corrects the energy of charge transfer states, as detailed in Section 2.4.

Direct Spin‐Orbit CI in Spin String Formulation

The key step in the iterative Davidson‐Liu algorithm is the calculation of the sigma vectors[ , ] which are necessary to update the CI vectors:

where denotes the CI matrix based on determinants and belongs to the set of expansion vectors constructed in the previous iteration. In order to make efficient use of pipelining and vector processing of modern computers, an integral‐driven (direct) algorithm has been implemented to construct the sigma vectors . The number of vectors to be constructed in each iteration is equal or similar to the number of roots to be determined, it can be smaller close to convergence and larger when restarting the Krylov space expansion (see below).

For a scalar Hamiltonian, the seminal work of Siegbahn as well as Knowles and Handy provides us with a well‐established and efficient full CI algorithm. The (string based) Knowles/Handy algorithm is probably the best full CI algorithm if the number of electrons and active orbitals is of about the same size, which is the case for typical valence‐type CASSCF calculations. Sometimes the use of spin‐adapted configuration state functions (CSFs) is advocated since this reduces the CI dimension. This however holds to a much lesser degree if spin‐orbit interaction is included, which mixes different‐spin CSFs.

Since the Knowles/Handy algorithm is well established, we will only dwell on its concept as much as required to demonstrate how to extend this to the spin‐orbit and Zeeman operators. Central to the algorithm is to write Slater determinants as direct products of so‐called α and β strings and that determine which of the α and β spin‐orbitals are occupied in the determinant

so that we can label the coefficients of a CI vector with a compound index denoting which of the α and β strings have to be combined to form the corresponding Slater determinant. Let us first consider the contribution from the operators

so the sum over disappears. Furthermore, for Slater determinants the coupling coefficients are very sparse, for a given pair p, q and α‐string there is at most one with a non‐zero coupling coefficient . It is easy to construct a table containing this information in compressed form, that is a table of integers . These are defined such that

This table is easily constructed and can be held in memory. In the worst case it contains about entries for N act active orbitals, which is a very small number compared to the CI dimension. The outermost loops go over to make the updates on the σ vector as memory‐local as possible, and the inner loops go over p,q. There is no loop over which (together with a parity factor) is simply . We note in passing that for different , different regions of the σ vector are updated. This allows us to distribute the work load associated with different on threads running in parallel without the need to synchronise their memory updates. A more detailed account of the parallelisation of the program is given in the Supporting Information.

So far, everything is fairly well known, but it sets the scene for calculating matrix elements of operators that change the number of α spin electrons. Consider the contribution

Again, memory locality considerations suggest that the outermost loops are over and the inner loops are over p,q. There is no loop over and since e. g. for p,q and given, there is no or only one for which a non‐zero value and that is . The only thing one needs are two additional tables

which are even easier to construct and store since they have only one active orbital index. Note that R + and R – connect strings where the number of electrons differ by so this contribution to the σ vector leads to a mixing of Slater determinants with different m values.

The lowest eigenvalues and eigenvectors of the CI matrix are obtained with the iterative multi‐root Davidson‐Liu method[ , ] which is a Lanczos Krylov‐space expansion method using a diagonal preconditioner. We have implemented both the original Davidson preconditioner as well as a modification due to Olsen . Besides the handling of complex vectors, the special feature in the target applications is that one may need to calculate a fairly large number of low‐energy roots (N roots) to treat multinuclear exchange‐coupled complexes. Since in each iteration up to N roots new vectors are added to the Krylov space, the expansion vectors can produce quite some memory demand, and one has to restrict this number. If adding the new expansion vectors would expand the Krylov space beyond this limit, one has to “collapse” the space and continue with the N roots best linear combinations of expansion vectors, augmented by the new expansion vectors. Our observation is that the Davidson convergence is slowed down if the Krylov space is restricted to less then expansion vectors, and our standard setup is to limit the Krylov space to expansion vectors. Note that in most iterations, the number of σ vectors that have to be constructed simultaneously is up to N roots, but in an iteration following a “collapse” it can be up to .

Start Vectors for the Direct CI Diagonalisation

Being an iterative method, the convergence of the Davidson procedure benefits from well‐chosen starting vectors, that is, trial vectors that have a significant overlap with the converged solutions. Apart from that, a bad starting guess may lead to convergence to excited states that are not desired in the calculation. Common procedures include first solving the eigenvalues problem in a subspace of important configurations, usually chosen by absolute values of diagonal elements of the Hamilton matrix, and then projecting them onto the full space in order to obtain a set of starting vectors. However, this approach does not account for strong coupling between different determinants. We found that these approaches have very limited success when treating oligonuclear complexes, and that it is vital here that the starting vectors to a fair degree span the space of the low‐energy solutions sought.

If there was no couplings between the metal centres, the total wave function would be the tensor product of the single‐centre multiplets, at least in the absence of internal redox reactions where one electron jumps from one centre to the other. These tensor products span, to a good approximation, the space of the low‐energy wave functions of the oligonuclear complex with weakly interacting metal centres. To be able to construct the single‐centre multiplets, the orbitals must be localised on the metal centres. For a single centre, the full CI dimension is rather small (d  5 case: dimension 252, f  7 case: dimension 3432) therefore the fragment CI wave functions can be calculated by conventional CI matrix diagonalisation. For the CoII d  7 case with tetrahedral coordination outlined above, the CI dimension is 120. There will be four low‐energy micro‐states (the ) well separated from all the others, and the tensor product of these four micro‐states on all three centres yields 64 wave functions which already have a large overlap with the final 64 low‐energy solutions.

We want to stress that using these tensor product start vectors is a key factor to achieve fast and smooth convergence in the direct CI iterations. Another view to look at these start vectors is to note that their number matches the number of desired roots, and that they do not contain contributions from wave functions with local ligand‐field excitations.

Modified CASCI for an Efficient Treatment of Exchange Coupling

When calculating magnetic properties of oligonuclear transition metal complexes a decent description of the exchange coupling is quite important, since errors made here can lead to a calculated ground state multiplet with the wrong total spin. While the magnetic properties of isolated centres are often reasonably well described at the CAS level, this is not at all the case for exchange couplings. The role of different contributions to the magnetic exchange coupling has been summarised in a recent review of Malrieu et al. With CAS‐type wave functions that miss dynamical correlation, one usually gets much too small (often by a factor of 3) antiferromagnetic contributions to the exchange couplings. However, sophisticated calculations including dynamic correlation effects on multinuclear exchange coupled metal complexes are rather scarce because of the massive computational requirements. One of the reasons for this is that dynamical correlation is usually taken into account before one includes spin‐orbit coupling, which leads to the spin‐orbit CI procedure whose problems for oligonuclear complexes have been outlined above.

An alternative approach to a reasonable description of the magnetic exchange coupling is the modified CASCI (MCASCI) method. This pragmatic and cost‐effective method is rooted in the observation that the coupling of the metal centres mixes highly excited charge‐transfer states into the low‐energy manifold, and that this mixing makes an important contribution to the exchange interaction. Since the orbital relaxation that would result if one electron jumps from one metal centre to the other cannot happen in the CAS‐type wave function, the charge‐transfer states are much too high in energy, therefore their mixing with the low‐energy manifold (and hence the exchange interaction) is underestimated. Using valence orbitals that are localised on the metal centres, one can easily identify the charge‐transfer Slater determinants and shift the corresponding diagonal element of the CI matrix by an independently calculated single‐centre relaxation energy Δ. While not being perfect, the MCASCI method often greatly improves exchange couplings at virtually no extra cost in the CASOCI calculation itself.

The non‐trivial part of the MCASCI approach is the calculation of the relaxation energies Δ for all possible “jumps” of an electron from one centre to another. In practice, these values are obtained in a series of open‐shell Hartree‐Fock or state averaged CASSCF calculations either on mono‐nuclear model complexes or on systems where all but one spin centre is replaced by a suitable diamagnetic ion, and to calculate the amount of orbital relaxation energy when adding or removing one electron from the open‐shell centre. The relaxation energies do not need to be very accurate, they just need to be in the right energy range. To our experinece, the approach works better if the orbitals for the full complex are obtained for the high‐spin states. Problems can occur when the charge transfer states are very low in energy and in cases where the high‐spin states themselves can be stabilised by charge transfer.

Relativistic Corrections and Spin‐Orbit Hamiltonian

Without spin‐orbit interaction, all g values would be isotropic, and there would be no zero field splittings of multiplets. Therefore, taking spin‐orbit coupling into account is of prime importance when calculating magnetic properties, even if metals are involved for which scalar relativistic effects are of lesser importance. For those systems, it is usually sufficient to use the non‐relativistic scalar Hamiltonian together with the leading‐order Breit‐Pauli spin‐orbit operator, which implies that nonrelativistic orbitals are used. To avoid explicit treatment of two‐electron spin‐orbit integrals in the CI calculation, an effective one‐particle spin‐orbit mean‐field (SOMF) operator[ , ] is used, which contracts the two‐electron spin‐orbit integrals with a spin‐free density matrix from a reference configuration. Note that without the two‐electron part of the SOMF operator, the spin‐orbit operator would not contain any shielding of the nuclear potential by the electrons, therefore it cannot be neglected. A systematic comparison of the SOMF approach, as opposed to using the full microscopic two‐electron spin‐orbit operator, has been presented recently by Netz et al. For mononuclear 3d complexes, the deviation for g matrices was negligible, and for the D tensors it was (with one exception) significantly below 10 %.

In compounds with elements heavier than metals, scalar‐relativistic effects have to be included already in the orbital optimisation step. Here our default option is the Douglas‐Kroll Hamiltonian, but it is also possible to use quasirelativistic effective core potentials. The (scalar) Hamiltonian used in the orbital optimisation step, together with the corresponding spin‐orbit part, is then used in the CASOCI calculation.

Implementation

The CASOCI program itself is a stand‐alone program which needs one‐ and two‐electron matrix elements in the basis of the localised active orbitals. These integrals are calculated with the integral code from TURBOMOLE[ , ] but could also be calculated with any other integral code. The molecular orbitals coefficients and the specification of active and inactive orbital sets have to be provided to the integral code, and then all necessary quantities such as the core Fock‐operator , the core energy , the two electron integrals g, and the spin orbit and Zeeman matrix elements , and can be calculated and stored. The CASOCI program only needs these matrix elements, whose number is very small since only active orbital indices are involved. So far the program has been used with CASSCF orbitals obtained either with MOLPRO or with the program of Meier and Staemmler , where we have programmed an interface to store orbitals and basis set data in TURBOMOLE format.

Magnetic Properties

In recent years, methods based on ab initio wave functions have been established as standard tools in single molecule magnetism. In this respect, molecular orbital and generally wave function analyses together with occupations or energies can be very helpful to interpret and finally understand magnetic data. However, to directly compare to and finally understand the experiment, it is necessary to be able to compute magnetic properties from ab initio data. Within single molecule magnetism, these properties generally divide into i. thermodynamic ensemble properties and ii. effective Hamiltonian parameters.

Magnetisation and Susceptibility

In the CASOCI calculation, there is a single molecule. Measurements of magnetisation and magnetic susceptibilities are made on powder samples which are in thermal equilibrium with the environment, so we have to consider the free energy of an ensemble of randomly oriented molecules. Instead of orienting the molecule in all possible directions, we can consider a molecule in a fixed orientation and average over all possible directions of the magnetic field. Boltzmann and rotational averaging then leads to a free energy expression as function of the magnitude B of the magnetic field and the temperature T:

with k the Boltzmann constant and a unit vector pointing to points on the unit sphere:

In principle one has to sum over all eigenstates of the molecule so the requirement is that the CASOCI spectrum encompasses all thermally accessible states. Note that the free energy only depends on the size of the magnetic field, while the energy levels of the molecule in a fixed orientation depend both on the size and direction. The (macroscopic) magnetisation M is now simply obtained as

where we have introduced the relative Boltzmann populations p which depend on both magnetic field and temperature. In an isotropic sample, the magnetisation (vector) is necessarily parallel to the applied magnetic field so we get a scalar quantity here, and the sign of M denotes whether the magnetisation is parallel or anti‐parallel to the applied field.

The magnetic susceptibility is defined as the second derivative and we get

The last term involving the second derivatives of the CASOCI energies w.r.t. the magnetic field contributes diamagnetism and temperature‐independent paramagnetism. This term is not included in our calculations as we cannot calculate the response of the orbitals to the applied magnetic field. Technically, we construct the matrix elements of the field‐dependent Hamiltonian in the space of our CASOCI wave functions and diagonalise the matrix representation to obtain the field‐dependent E, while their derivatives are obtained as expectation values of the resulting (field‐dependent) wave functions with the Zeeman operator. The integral performing the rotational averaging is evaluated by numerical quadrature on a unit sphere. We use Lebedev integration grids which provide a quite regular triangularisation of the sphere. A 110‐point grid is sufficient in all cases studied so far. Note that the integrand is symmetric w.r.t. inversion, so for any pair of grid points related by inversion, only one needs to be treated explicitly (this reduces the effective number of grid points to 55).

Experimental measurements are usually performed in magnetic fields of about 0.1 to 1 Tesla, in this region the magnetisation is nearly linear in the magnetic field such that and χ are numerically very close. The ratio is measured in experiments with a constant magnetic field where one measures the ratio of the magnetisation of an unknown and a reference sample, while χ (the derivative of M) is measured when one modulates the applied magnetic field and measures the modulation of the magnetisation. This subtle difference is only visible at very low temperatures when saturation sets in.

Spin Hamiltonian Parameters

For the analysis of magnetic properties effective Hamiltonian theory is usually employed, e. g. in order to characterise, compare and assess the spectroscopic signatures from electron paramagnetic resonance (EPR) experiments. In this approach, the measured data is modelled by means of the properties of a parametrised model Hamiltonian acting solely on pseudo‐spin functions. E.g., for single ion magnets or in the strong exchange limit

Here, the single‐ion D‐tensor describes the splitting of a pseudo‐spin multiplet at zero‐field, the g‐matrix parametrises the coupling to an applied magnetic field . When no high resolution high‐field EPR data is available, the parameters in Eq. (23) are usually obtained from fitting the (rather featureless) temperature dependence of the magnetic susceptibility, and from symmetry arguments. In these cases it is very helpful if such parameters can be extracted from ab initio wave functions. However, it has to be kept in mind that such a mapping between ab initio and spin Hamiltonian wave functions requires that one starts with a selection of ab initio micro‐states forming a multiplet, and this is only possible if there is a selection of ab initio micro‐states Ψ with energies E that are close in energy and well separated from all other states. What is needed then to extract spin Hamiltonian parameters are these energies as well as three matrices containing the matrix representation of the cartesian components of the Zeeman operator, namely

Magnetic axes and g‐values can be computed from the Abragam‐Bleaney tensor . Based on the ansatz from Gerloch and McMeeking , a straightforward extension to arbitrary pseudo‐spins has been given by Chibotaru:

Note that is invariant to an arbitrary unitary transformation of the multiplet functions Ψ . Diagonalising yields the magnetic axes and squares of g‐factors.

In order to extract the D‐tensor and the g‐matrix from ab initio data, one has to start with a mapping of unitary linear combinations of the ab initio multiplet functions onto eigenfunctions of the pseudo‐spin operator , denoted as , conforming to the Condon‐Shortley phase convention

with a unitary matrix . One can then construct the zero‐field and Zeeman Hamiltonian matrix elements in the basis of the and determine the spin Hamiltonian parameters such that the ab initio and spin Hamiltonian matrix elements match as good as possible . In the general case, an exact match can only be obtained if one extends the spin Hamiltonian by higher‐order (Stevens) operators. The parts of the spin Hamiltonian involving and are the lowest‐order terms of such an expansion. Determining the parameters of these lowest‐order terms is equivalent to obtaining a least‐squares match of the non‐zero matrix elements of the spin Hamiltonian with the ab initio matrix elements. Expressions how to calculate and from the zero‐field and Zeeman matrix elements are given in Ref. [36].

The results of this procedure depend on the mapping chosen (matrix ). If one wants to treat a Kramers doublet well separated from all other states using a pseudo‐spin , it is not clear from the outset which linear combination of the two degenerate ab initio micro‐states is “spin up”. For cases with very weak spin‐orbit splitting a mapping may be chosen which produces approximate eigenfunctions of the true (microscopic) spin operator , for lanthanides the atomic total angular momentum operator can be used instead.

A universal approach for obtaining an unique mapping has been given by Chibotaru. In this ansatz, one first defines a cartesian reference frame (often the magnetic axes of the G‐tensor are used) and diagonalises the Zeeman operator for a magnetic field in the z‐direction of the reference frame, and sorts the resulting eigenfunctions by their Zeeman interaction. The rationale for such a mapping is that it produces a g‐matrix with and . At the end, each mapping leads to different that are equally “true” if it comes to reproducing the physics, but what one wants are spin Hamiltonian parameters that one can also interpret, and the behaviour of a pseudo‐spin with very large off‐diagonal g‐matrix elements is against our intuition which tells us that a magnetic field should align a spin in its direction.

The diagonalisation has not yet specified the (relative) phase of the functions . To fix the phases, one considers for the matrix elements

for the Zeeman operators for a magnetic field in the x and y direction of the reference frame. One then tries to choose the relative phases such that the matrix elements X are on the positive real and Y are on the positive imaginary axis in the complex plane. X real leads to and implies , while Y purely imaginary yields and the imaginary part being positive makes (see Eq. (33) in Ref. [36]). Of course, with a single phase factor one cannot fulfil all four conditions, so one can, for example, choose the phase as to minimise and fix the sign to make or positive, whichever is larger in absolute value. In normal cases (with positive g values), this mapping will produce a g‐matrix where the diagonal elements are positive, g and g are zero, and g and g are small.

Results

Zero‐Field Splitting in a Sulfur‐Bridged Co‐V‐Co Complex

To test the functionality of the new CASOCI implementation, we investigated the zero‐field splitting in a trinuclear CoIIVIICoII complex [Co2V(SPh)8]2− shown in Figure 1, the coordinates are given in Section S2.1 in the Supporting Information. This quasi‐centrosymmetric system was studied before by means of mono‐ and binuclear model complexes using the standard two‐step SOCI approach, for more details see Ref. [32].

Figure 1

Trinuclear CoIIVIICoII system used to test the functionalities of the CASOCI program. For visualisation, VMD has been used. The active space contains 17 electrons in 13 orbitals. Co: magenta, V: blue, S: yellow, C: black, H: white.

Without spin‐orbit coupling, all three ions have 4 A ground states arising from the d  7 configuration in tetrahedral coordination of CoII and the d  3 configuration in octahedral coordination of VII, respectively. The next states follow at 2600 cm−1 for CoII and at 17000 cm−1 for VII. The corresponding active space is constructed from 17 electrons in 13 orbitals (five 3d orbitals at each CoII and three 3d orbitals at VII). Because of the energetic separation from ligand‐field excited states, magnetic anisotropy in this complex is generated by second order spin‐orbit coupling. While the contribution from quasi‐octahedral VII is negligible, at least the seven quartet states corresponding to the 4 F state of the free ion have to be considered for each tetrahedral CoII, yielding micro‐states that need to be considered in a SOCI. At the time of our former study, such a calculation was not possible and computations in smaller subspaces, e. g., the first spin‐ladder (64 micro‐states) were futile. Consequently, we deduced the electronic structure of the trinuclear complex from the spin‐Hamiltonian

(with indices 1,3 referring to CoII and 2 to VII) whose parameters were obtained from mono‐ and binuclear model complexes. The result was a well‐separated pseudo‐quartet state which is split by zero‐field splitting into two Kramers doublets in agreement with experimental results. Thus the magnetic properties at low temperatures of this compound are dominated by 4 spin‐orbit coupled states only.

With the new CASOCI implementation, such a calculation is possible to carry out for the whole trinuclear complex and has thus been performed using def2‐TZVP basis sets for Co, V and S and def2‐SVP for C and H. Orbitals were obtained from a state‐averaged CASSCF calculation for decet states, originating from coupling the 4 F ground states of two CoII ions with 4 A of VII. The energies of the CASSCF states are given in the Supporting Information, Section S2.2, the corresponding active orbitals in Figure S2. The CASOCI calculations were performed with and without the MCASCI approach. In the latter, a relaxation energy of 0.4 Hartree for charge transfer between Co and V was used. All energies obtained with the MCASCI approach are denoted as'+ Shift’ in the following. The CASOCI subspace included micro‐states from the ground manifold, start vectors were obtained from the tensor product states, information about the zero‐field splitting of the individual centres was thus obtained simultaneously in this computation and is listed, together with the (very similar) data from the calculations on mono‐nuclear model complexes in Table 1. The relative energies of the 64 lowest states for the different approaches are summarized in the Supporting Information, Section S2.3. Results of the subspace CI are given in Section S2.4 of the Supporting Information together with results on complexes were two metal centers were substituted either by the diamagnetic MgII or ZnII ions. The CASOCI+Shift energies are plotted in Figure 2 against pseudo‐spin multiplicities. The lowest states reflect the spin‐ladder spectrum generated by the spin Hamiltonian Eq. 29 fairly well, confirming the second‐order nature of spin‐orbit coupling in this complex. From around 2500 cm−1 and above, such a mapping is more difficult due to mixing with the lowest ligand‐field states (indicated by the dashed line). The convergence behaviour of the CASOCI calculation again demonstrates the advantage of using tensor product states. Before the Davidson procedure, the energy splitting of these 64 states is caused by spin‐orbit coupling only and is in the range of 700 cm−1. However, already after the first iteration, the states are reordered due to exchange coupling and show the typical spin‐ladder spectrum with a span of about  cm−1 with and 1300 cm−1 without shift, respectively. Note that in this complex, we have a rather special case with an exceptionally strong exchange coupling (  cm−1, Ref. [32]) between adjacent centres and at the same time very low‐lying (  cm−1) single‐centre ligand‐field excited states, therefore the latters can intrude into the CASOCI states in the MCASCI approach, causing convergence problems for the highest roots. The 4 lowest micro‐states, however, which dominate the magnetic properties at low temperatures, converge in only 27 Davidson iterations which is quite rapid considering the strong exchange coupling. The zero‐field splitting parameters computed from the ab initio wave functions are  cm−1,  cm−1 with MCASCI shifts and  cm−1,  cm−1 without. These results are close to the previously reported values obtained by fitting the temperature dependence of the magnetic moment obtained experimentally (  cm−1, E was not considered in Ref. ).

Table 1

Zero‐field splittings ( ) of the trinuclear CoIIVIICoII system. D is obtained as half the energy difference between the Kramers doublets. In Ref. [31], a scaled nucleus approach with a scaling factor of 0.61 was used for the spin‐orbit integrals.

D_Co	D V	D_Co_2V	Reference
−31.6	0.05	-29.9^a	[31]
-27.1^b	0.01^b	−26.8	this work (CASSCF)

From projection of three‐center onto single‐ion spin Hamiltonian. From subspace SOCI (see Table S4). The sign of D cannot be obtained from the subspace SOCI.

Figure 2

CASOCI energies (with shift) of the trinuclear CoIIVIICoII system, mapped to pseudo‐spin multiplicities. Above  cm−1, such a mapping is difficult due to mixing with ligand‐field excited states.

The magnetic powder data was computed from the CASOCI approach by equation 22 using the lowest 18 states and is depicted in Figure 3. As can be seen in Figure 2 and Section S2.3 in the Supporting Information, there is an energy gap after 18 states and the higher states should not contribute to the magnetic susceptibility up to room temperature. The calculated magnetic susceptibility is larger than the experimental one at all temperatures. For further analysis of the influence of different parameters, a one‐centre and a three‐centre spin Hamiltonian were used respectively. Details are summarized in the Supporting Information, Section S2.5 and Figure S3. In the one centre Hamiltonian, g 1 =2.31 (2.32), g 2 =2.34 (2.40), g 3 =2.87 (2.91) was obtained for the lowest quartet state. The numbers in parentheses refer to CASOCI without MCASCI shifts. For the three‐centre model, the values for CoII are taken from a mono‐nuclear model where VII and one of the CoII were diamagnetically substituted by MgII.

Figure 3

Magnetic susceptibility based on 18 CASOCI states without shift (blue line) and with shift (black line) of the trinuclear CoIIVIICoII system. Experimental values (black crosses) taken from Ref. [31], 3‐centre Hamiltonian with (magenta line).

Compared to other systems, the usage of the MCASCI approach did not improve the results. It seems that in this case with the very strong coupling along the metal‐metal axis, the antiferromagnetic coupling is overestimated by MCASCI. The slight increase of χT at higher distances is caused by higher lying states and is better described by the calculation without shift. Using the spin Hamiltonian parameters obtained in the way described in Section 2.7, the susceptibility is similar to the calculated values. However, a global reduction of the three components of g for the quartet ground state in the one centre model or for the ground state of CoII in the three‐centre Hamiltonian with a factor of 0.9 significantly improves the agreement with the experimental data.

Note that in Ref. [32], the optimal fit of the temperature dependence of the effective magnetic moment was obtained with a very small isotropic g‐values of 1.83 for VII and 2.12 for CoII, respectively. In particular the g‐factor for VII is in disagreement with the calculations where a g‐factor very close to g is obtained as expected for a d 3 system in quasi‐octahedral environment.

Within the giant spin approximation (strong exchange limit), one obtains

from the spin projection coefficients for the lowest state. Therefore, setting and leads to almost the same molecular g value, and adjusting  cm−1 then yields almost the same temperature dependent susceptibility as the original fit of the experimental data (see Figure 3). With these parameters also the slope of the curve in the low temperature regime is well described, while the curve obtained from the ab initio data has a slightly different slope in this regime which seems to be connected to the anisotropy of the calculated g‐factor.

However, the calculated anisotropy of g in the mono‐nuclear CoII model is in the same range as the values obtained in Ref. [39] where the magneto‐structural correlation of [Co is investigated and CASSCF and n‐electron valence state perturbation theory (NEVPT2) results are compared, indicating only a small decrease of the zero‐field splitting and of the g‐factor anisotropy. Taken into account that a minimal active space was used in a purely ab initio approach, and also uncertainties in the experimental measurements cannot be excluded, the agreement between experiment and theory is reasonably well.

Conclusion

We present a new SOCI program, termed CASOCI, which is designed for the description of magnetic properties of oligonuclear 3d‐ and 4 f‐metal compounds. In the CASOCI program, spin‐orbit interaction is directly included into the diagonalisation, coupling all micro‐states in the CAS. This approach bypasses the necessity of the prevalent two‐step procedure to compute a whole set of scalar quasi‐degenerate states before the SOCI. In oligonuclear complexes, this number grows exponentially, whereas there are typically only a few low‐lying spin‐orbit coupled states, thus favouring the one‐step approach. Our approach is accompanied by the MCASCI ansatz for treating strongly exchange coupled systems efficiently.

In contrast to Ganyushin and Neese in Ref. [14], our program is based on Slater determinants, which has the advantage that the corresponding algebra is ideally suited for modern computers. The efficiency of the program for oligonuclear complexes is further increased by the use of spin‐orbit coupled tensor product start vectors which solve the magnetically uncoupled problem and are easily generated at the beginning of the calculation, with the additional benefit of giving insight into the zero‐field splitting of the individual centres.

To arrive at a program valuable for application calculations, the functionalities of the program were extended to compute magnetic properties directly from the ab initio wave functions. To demonstrate this, a trinuclear CoIIVIICoII complex was chosen as example. Using the conventional two‐step procedure would require to couple more than 3000 scalar ligand‐field states, while the magnetic properties are well described by a much smaller number of spin‐orbit coupled states, a typical example where the the one‐step approach is well‐suited. The D‐tensor and g‐matrix of the lowest pseudo‐quartet, as well as magnetic powder data for a temperature up to 300 K was obtained based on ab initio data using the new program.

Though we have shown only results for d‐element compounds, the program was successfully used on systems where one of the magnetic centers is an organic radical ligand, and on compounds with centres. For ‐compounds, configuration‐averaged methods like restricted open shell Hartree Fock (ROHF) are an excellent and cost‐effective alternative to state‐averaged CASSCF for the optimisation of the molecular orbitals. This opens the possibility for treating extended systems containing weakly coupled lanthanide ions, that have recently gained popularity in the single‐molecule magnetism community due to their potentially highly axial magnetic anisotropy.

Supporting Information

Supporting Information is available free of charge (PDF document). Supporting information contains a detailed description of the hybrid OpenMP and MPI parallelisation of the CASOCI program including a performance measurement, as well as computational details about molecular structures, active orbitals, CASSCF and CASOCI total energies, and simulation of magnetic properties of the CoIIVIICoII complex.

Conflict of interest

The authors declare no conflict of interest.

As a service to our authors and readers, this journal provides supporting information supplied by the authors. Such materials are peer reviewed and may be re‐organized for online delivery, but are not copy‐edited or typeset. Technical support issues arising from supporting information (other than missing files) should be addressed to the authors.

Supporting Information

Click here for additional data file.