Four-Reference State-Specific Brillouin-Wigner Coupled- Cluster Method: Study of the IBr Molecule

We implemented the state-specific Brillouin–Wigner coupled-cluster method for the complete model space spanned by four reference configurations generated by two electrons in two active orbitals. We applied the method (together with the previously suggested a posteriori size-extensivity correction) to the calculation of spectroscopic constants of the IBr molecule, using averaged relativistic effective core potential.


Introduction
This paper is a continuation of our previous studies [1][2][3][4][5][6][7][8][9] on the development of a multi-reference coupled-cluster (MRCC) method that would be free of the problem of intruder states and that would be amenable to treatment of systems requiring more than two reference configurations.Avoidance of intruder states was achieved by using the Brillouin-Wigner (BW) resolvent, and simplicity of the method and feasibility of calculations were achieved by subjecting the BWCCSD method to a state-specific form [2][3][4].Although the method so developed was general in respect to the number of reference configurations, its implementation [4] in the ACES II program [10] was limited to two-reference cases, where HOMO and LUMO orbitals have different spatial symmetry and only two closed-shell configurations can contribute to the wave function.In order to be able to treat molecules with quasidegenerate HOMO and LUMO orbitals with the same spatial symmetry, we extended our original implementation of the BWCCSD method for open shell reference configurations.For merits of the statespecific BW theory we paid the price that the method was not size extensive any longer.We developed therefore an a posteriori correction [5,6] for eliminating the size-inconsistent terms in the amplitudes.Again, the correction in its original implementation [6] was confined to two-reference cases.In this paper it has been generalized and calculated for four reference configurations.
The IBr molecule was selected for testing.Its HOMO and LUMO are of the same symmetry and the two singly excited configurations (determinants) HOMO → LUMO have to be included in the reference space.One of us (J.P.) was a coauthor of a paper [11] dealing with single-reference calculations on IBr.Results of these calculations are used in this paper for judging the performance of BWCCSD.Finally, IBr was selected for that reason that this molecule is also interesting experimentally and the MR BWCCSD calculations could be helpful for the theoretical interpretation.
The IBr has been studied for several decades.The early measurements of spectra of IBr in visible [12] and ultraviolet [13] have found bands corresponding to transitions to higher electronic states.Further studies confirmed these findings [14] and gave refined values of rotational and vibrational constants from high resolution far infrared absorption spectra [15].The value of dissociation energy was obtained from isotopically selective velocity mapping measurements [16].The predissociation of IBr has been studied by resonance Raman spectroscopy [17], pump-probe experiments [18] and time dependent wave packet dynamics [19].Recently ion imaging methods were used to investigate the photolysis of IBr [20].
Early theoretical calculations of electronic structure on IBr [21] were performed at Hartree-Fock level.Subsequent ab initio studies with inclusion of relativistic effects provided values of spectroscopic [11,22,23] and electric [22,24,25] properties.

Theory
Since the derivation of the state-specific Brillouin-Wigner coupled-cluster theory has been presented for the general multireference case previously [3,5], and details of the implementation for closed-shell molecules have also been published [4], we give here only a brief review of the method and the relevant working equations.
We use the notation that i, j and a, b indices represent occupied and unoccupied spin-orbitals, respectively, while p, q are generic spin-orbital indices.Moreover, we reserve the letters k and c for denoting the internal excitations (thus these indices are fixed for a given pair of reference configurations µ, ν), while i, j, a, b are used as general or summation indices.For many systems, including the IBr molecule, k is the HOMO and c is the LUMO orbital and the four (spin-unrestricted) reference configurations are generated by all possible occupations of these two orbitals by two electrons.
The model function for the ground state is assumed in the form where M is the number of reference configurations.For the exact ground-state wave function Ψ 0 and exact energy 0 (2) and 0 0 0 The "effective" Hamiltonian H eff is defined as where P is the projection operator onto the model space spanned by M configurations Assuming the CCSD cluster operator T(µ) = T 1 (µ) + T 2 (µ) with respect to configuration µ as Fermi vacuum, the wave operator Ω 0 in the Hilbert space ansatz [26] is subject to the state-specific Brillouin-Wigner analogue of the Bloch equation where B 0 is the Brillouin-Wigner resolvent 0 0 Here 0 E % denotes the exact energy (lowest eigenvalue of H eff ), while E q are the unperturbed energies corresponding to Φ q similarly as in the Rayleigh-Schrödinger resolvent, and the notation q >M in Eq. (8) means that internal excitations (the ones relating Φ µ and Φ ν , µ,ν ≤ M) are excluded.
The diagonal matrix elements of the effective Hamiltonian read where H µµ is Hamiltonian expectation value for reference configuration Φ µ and H N (µ) is the Hamiltonian normally ordered with respect to Fermi vacuum Φ µ .They can be computed by a modified routine for energy in single-reference CCSD [4].Coupling between the reference configurations is furnished by the off-diagonal elements These have to be treated separately for cases when Φ ν and Φ µ are single-or double-excitations with respect to each other.In the former case we can write either and the matrix elements are obtained as r.h.s. of the T 1 amplitude update equations.In the latter and the matrix element is obtained from the T 2 amplitude equations as described in [4].Presently, our implementation does not support reference configurations which would be mutually more than biexcited.
The final equations for T 1 amplitudes read { } where the r.h.s. is the same as in the case of single reference CCSD, except that the amplitudes of internal excitations are set to zero.In the T 2 case we have where P ij is antisymmetrization operator acting on the i, j indices.The r.h.s.term in curly brackets can be identified with the r.h.s. of the single reference CCSD equations computed with internal amplitudes set to zero.As described in [6], the amplitudes obtained by this method need to be corrected for retaining sizeextensivity.The idea is to make BWCC a posteriori close to its Rayleigh-Schrödinger (RS) version.
Since our BW variant of the Bloch equation ( 7) contains only a denominator of the type 0 ( ) where The first two terms on the r.h.s. of Eq. ( 13) are size-extensive because they may be viewed as the analog of the usual RS form of the Bloch equation.However, the "additional" third term on the r.h.s. of Eq. ( 13) gives, on iterating the BWCCSD equations, size-inextensive terms.Computationally, it is simplest to identify the size-inextensive terms in the amplitude update equations and eliminate them in an additional a posteriori iteration.For T 1 amplitudes the only term responsible for size-inextensivity is on the l.h.s of Eq. (11).With the T 2 amplitudes it is the 0 ( )( ) term on the l.h.s of Eq. ( 12) and the last term on the r.h.s. of Eq. (12).
From the corrected amplitudes we construct a new H eff matrix, and by its diagonalization we obtain the final energy.Note that the size-inextensive terms cannot be eliminated during iterations, since the only coupling between amplitudes of different reference configurations, provided by 0 E % , would be lost.
The BWCCSD calculation proceeds iteratively.After a standard initial guess of amplitudes, the H eff matrix elements are calculated and H eff is diagonalized.The lowest eigenvalue is then used as 0 E % and new T 1 and T 2 amplitudes for all reference configurations are computed according to Eqs. ( 11) and ( 12).This procedure is repeated, employing the DIIS convergence acceleration, until amplitudes converge.Subsequently the a posteriori size-extensivity correction is performed in an additional iteration.
Note that the absence of explicit coupling of amplitudes from different reference sets makes the method simple and computationally feasible for larger molecules.The computational demands per iteration for the M-reference BWCCSD method are thus only marginally higher than M times the demands for single-reference CCSD procedure.

Computational
In molecules containing heavy atoms, such as iodine or bromine, relativistic effects, and spin-orbit (LS) interactions in particular, cannot be neglected and in some cases even first-order perturbative treatment is not accurate enough.Moreover, in order to obtain potential energy (PE) curve accurate enough in the whole range of internuclear distances, inclusion of both statical and dynamical electron correlation in a sophisticated way is absolutely essential.The ground state of the IBr molecule is 1 Σ+.
Since it has zero angular and spin momentum, it is not subject to spin-orbit splitting, however, scalar relativistic effects play an important role.Therefore in our calculations the inner-shell electrons have been treated using recent averaged relativistic effective core potentials (AREP) [27].We employed two valence basis sets: the basis set supplied together with the AREP of the size (6s6p1d/3s3p1d) for both I and Br labeled A, and basis set augmented by diffuse functions of size (8s8p3d/5s5p3d) for both atoms, labeled B, which is described in [11] and serves mainly for comparison with that paper.For each internuclear geometry, we performed first a CASSCF calculation with two electrons in the HOMO and LUMO active orbitals, i.e. the same active space as in BWCCSD.The canonical CASSCF orbitals were used in the subsequent BWCCSD treatment based on the four spin-unrestricted reference configurations possible in this active space.Since the BWCCSD calculation is performed in a spin-unrestricted form, three eigenvalues of H eff correspond to singlet states and one to a triplet.It has been checked that no spin contamination occurs in numerical procedure.The CASSCF was chosen for the description of the reference state because it is size extensive and correctly describes the static correlation.
For all coupled cluster calculations reported in this paper, we used our implementation of the BW CCSD method [4] and of the a posteriori correction [6] within the ACES II program [10].The MR-CISD results reported in this work were obtained using the GAMESS-UK program [19].
The spectroscopic constants were obtained from the polynomial expansion of the energy by means of the Dunham analysis.D e values were obtained as difference between the minimum and the dissociation limit of the potential energy curves.Rovibrational levels (cf.Table 2) have been calculated by the program LEVEL [17].

Results and Discussion
The MR-BWCCSD potential energy curves of IBr computed in basis sets A and B are shown in Figs. 1 and 2, respectively.CCSD curves are shown for comparison and it can be seen that the single reference method fails to meet the correct dissociation limit indicated by the dashed line.This clearly shows the need for a multireference description for larger interatomic distances, while near the equilibrium the single-reference description is quite reliable.
The calculated spectroscopic constants are compared with the sample results of previous works [11,22,23] in Table 1.We also computed a set of rovibrational levels, given in Table 2.
First we tested the size extensivity, and we found that the size-consistency error was small.It was only 0.10 kcal/mol with the basis set A and 0.11 kcal/mol with the basis set B. This is because the CASSCF orbitals are localized for large interatomic distances, which is beneficial for size consistency, making thus BWCCSD near size consistent even without the a posteriori correction.
On comparing the calculated spectroscopic constants with those obtained from single-reference calculations [11], using both the same basis set B and the pseudopotential, it can be seen that the

Dissociation curve of IBr
Basis Set A Energy (kcal/mol) Interatomic distance (a.u.)  dissociation energy D e has decreased from 47.7 kcal/mol (CCSD) to 38.0 kcal/mol (MR BWCCSD), due to the multireference description at large internuclear separation.The equilibrium bond length r e has increased by approximately 0.01 Å, while the vibrational frequency ω e has decreased slightly and the anharmonicity ω e x e remained almost same.They are also very similar to values obtained by the MR-CISD method (cf.Table.1).The lack of improvement of the multi-reference method for r e , ω e and ω e x e with respect to single reference CCSD can be explained by a low contribution of the excited references near the equilibrium geometry and therefore effective reduction of BWCCSD to a single reference CCSD, since the values of r e , ω e , and ω e x e are determined by the shape of the PES curve near equilibrium.Additional inclusion of dynamical correlation by means of connected triples, when implemented, should result in an improvement of the obtained spectroscopic constants.As our experience with MR-BWCCSD shows, our previous results are very close to those of stateuniversal MR-CCSD [4,8,31,32].We believe that in case of IBr this will be similar.Unfortunately we did not find such a result from state-universal MR-CCSD method for comparison.
The dependence of the calculated values of r e , ω e and ω e x e show a very small difference between the two basis sets used.That means that adding extra diffuse functions to the basis set does not have a large effect on the ground state wave function in the equilibrium interatomic distance range.However the D e value, which is slightly more sensitive, has changed by approximately 1.5 kcal mol -1 towards the experimental value.

Conclusions
The BWCCSD method developed previously has been shown in this paper to be computationally feasible for a four-reference problem.The calculations performed for the IBr molecule yielded a smooth potential curve in the whole range of interatomic distances assumed (4 to 80 au).The potential curve has the correct dissociation limit, and the spectroscopic constants compare well with the experimental and best available theoretical data.The results obtained for IBr also provide another numerical evidence that the size-extensivity error of the MR BWCCSD method is small.

Figure 1 .Figure 2 .
Figure 1.Potential energy curves of the IBr molecule calculated by the CCSD and MR-BWCCSD methods in the basis set A (see text) shown by solid lines.Energy scale is relative to the sum of energies of isolated atoms, which is indicated by the dashed line.Notice the incorrect dissociation limit of the single-reference CCSD.

Table 1 .
Spectroscopic constants of IBr b Using a relativistic effective core potential.c The same reference configurations as in MR-BWCCSD; without size-extensivity correction.d The same reference configurations as in MR-BWCCSD; with Davidson correction.e Based on RHF molecular orbitals.f Based on CASSCF molecular orbitals.

Table 2 .
Rovibrational energy levels of IBr obtained from potential curve calculated by MR BWCCSD in basis B (cm -1 )