Improved Segmented All-Electron Relativistically Contracted Basis

Feb 3, 2016 - CERES: An ab initio code dedicated to the calculation of the electronic structure and magnetic properties of lanthanide complexes. Simon...
0 downloads 0 Views 971KB Size
Article pubs.acs.org/JCTC

Improved Segmented All-Electron Relativistically Contracted Basis Sets for the Lanthanides Daniel Aravena,†,‡ Frank Neese,*,† and Dimitrios A. Pantazis*,† †

Max Planck Institut für Chemische Energiekonversion, Stifstr. 34-36, 45470 Mülheim an der Ruhr, Germany Facultad de Química y Biología, Universidad de Santiago de Chile (USACH), Casilla 40, Correo 33, Santiago, Chile



S Supporting Information *

ABSTRACT: Improved versions of the segmented all-electron relativistically contracted (SARC) basis sets for the lanthanides are presented. The second-generation SARC2 basis sets maintain efficient construction of their predecessors and their individual adaptation to the DKH2 and ZORA Hamiltonians, but feature exponents optimized with a completely new orbital shape fitting procedure and a slightly expanded f space that results in sizable improvement in CASSCF energies and in significantly more accurate prediction of spin−orbit coupling parameters. Additionally, an extended set of polarization/correlation functions is constructed that is appropriate for multireference correlated calculations and new auxiliary basis sets for use in resolution-ofidentity (density-fitting) approximations in combination with both DFT and wave function based treatments. Thus, the SARC2 basis sets extend the applicability of the first-generation DFT-oriented basis sets to routine all-electron wave function-based treatments of lanthanide complexes. The new basis sets are benchmarked with respect to excitation energies, radial distribution functions, optimized geometries, orbital eigenvalues, ionization potentials, and spin−orbit coupling parameters of lanthanide systems and are shown to be suitable for the description of magnetic and spectroscopic properties using both DFT and multireference wave function-based methods.

1. INTRODUCTION All-electron basis sets for the lanthanides have become more common in the past few years, in connection with their use in computational treatments that incorporate relativistic effects explicitly within the Hamiltonian.1−6 To address the need for compact all-electron basis sets that are adapted to the popular scalar relativistic second-order Douglas−Kroll−Hess (DKH2) and zero-order regular approximation (ZORA) Hamiltonians,7−15 we had proposed a family of segmented all-electron relativistically contracted (SARC) basis sets.16 Their design followed simple rules and an even-tempered construction in order to maintain simplicity, while their overall size was kept small to facilitate routine application in molecular chemistry. They were adapted separately in their contraction coefficients for the DKH2 and ZORA Hamiltonians and were targeted toward density functional theory (DFT) calculations, which are known to converge relatively quickly for molecular properties with respect to basis set size. For this reason, although the basis sets were shown to be flexible enough for DFT calculations, © 2016 American Chemical Society

they did not include higher angular momentum functions that would be necessary to sufficiently recover correlation energy in post Hartree−Fock methods. Since the original report, the SARC basis sets for the lanthanides have been used successfully in a variety of DFT applications by many research groups, see for examples, refs 17−29. As the field of applications became wider, some limitations of the original construction were also noticed. Specifically, the growing interest in the computational modeling of magnetic and spectroscopic properties of lanthanide complexes by multiconfigurational methods, such as the complete active space self-consistent field method (CASSCF),30−32 requires a type of basis set construction that is better optimized for this type of calculation than the original DFT-oriented SARC basis. An important factor in this respect is the f shell, which should be flexible enough to simultaneously Received: November 4, 2015 Published: February 3, 2016 1148

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156

Article

Journal of Chemical Theory and Computation

procedure, we adopted an even-tempered sequence for the exponent generation for each angular momentum, similar to the original SARC basis set, through the following series:

describe electronic excitations and also the state mixing due to spin−orbit coupling, which critically depends on the shape of the f orbitals. The relative inflexibility of the f shell in the original SARC formulation was also noticed in their suboptimal performance for the fourth ionization energy.33 Furthermore, the inclusion of dynamical correlation by means of perturbative methods (i.e., N-electron valence state perturbation theory, NEVPT2, 34−37 or complete active space second-order perturbation theory, CASPT238−41) poses additional requirements in the form of high angular momentum functions for the description of orbital excitations. An example of all-electron basis sets typically used in wave function treatments of lanthanide systems are the high quality relativistic atomic natural orbital (ANO-RCC) basis sets.6 Being a generally contracted basis, the ANO basis allows for the accurate description of molecular systems with a small number of functions that are in turn composed by a large number of primitives, but this feature requires dedicated integral handling routines to be efficiently implemented. The above considerations lead us to propose here a secondgeneration version of the SARC basis sets, SARC2. Our aim is to enable routine but reliable wave function calculations of lanthanide complexes and to improve accuracy for DFT calculations, while retaining the compactness and efficiency of the first-generation basis sets as much as possible. It is our stated goal to maintain the construction principles originally employed, i.e., to build an even-tempered basis with a relatively small number of functions. We also maintain the previously used contraction pattern to achieve a relatively decontracted basis set, where only the tightest functions are contracted in order to not compromise the description of the chemically relevant valence part. On the other hand, the parameter optimization approach was completely modified to an orbital shape-fitting procedure and a polarization set is proposed based on the maximization the dynamical correlation energy recovered by NEVPT2. Finally, auxiliary basis sets for resolution of the identity (RI) methods for the approximation of two-electron integrals are presented to work for NEVPT2 and CASSCF calculations.

ζl , i = αlx−i

i = 1, ..., Nl

(1)

where ζl,i is the i-th exponent for angular momentum l, αl represents the tightest exponent, x is related to the ratio between consecutive exponents for a given angular momentum, and Nl is the number of basis functions for a given angular momentum. Consequently, for a fixed number of functions, the basis set generation procedure involves the optimization of only two parameters per angular momentum. The parameters αl and x were optimized to fit the radial part of the occupied orbitals obtained from spin-averaged Hartree− Fock (SAHF) calculations for lanthanide atoms and ions in different electronic configurations by employing a large, totally decontracted basis set. For this purpose, tighter and looser functions were successively added to the universal Gaussian basis set (UGBS) of de Castro and Jorge42 until convergence of the total energies to below 1 mEh was obtained for all configurations (this extended reference basis will be further referenced as “UGBS+8”). We considered three atomic configurations (4fn+16s2, 4fn5d16s2, and 4fn6s26p1) and also one ionic (4fn) configuration in the generation of the reference orbitals, as the +3 oxidation state is the most common in molecular compounds. All calculations were performed by employing the DKH2 Hamiltonian for the consideration of scalar relativistic effects. ZORA versions of the SARC2 basis sets were also created and benchmarked. The ZORA-SARC2 basis sets resulted in very similar performance as the DKHSARC2 basis sets, and these results are presented in the Supporting Information. Once the reference orbitals were generated for all configurations, the norm of the projection of each orbital on the target basis set was calculated, and the error function was expressed as log(1−c), where c is the product of the projection norms for each orbital for a given angular momentum. The rationale about this choice of error function can be understood as follows: in practice, the tested basis sets must be able to represent the reference functions as close as possible to yield orbital energies close to the reference values. The overlap between reference and trial functions is always close to 1, while apparently small differences in overlap can be important for basis set performance, for instance, between overlaps of 0.999 and 0.9999. In order to set an optimization function that is sensitive to such small differences in value, we decided to use log(1−c) as the error function, as it will yield an increasingly small number as the overlap approaches 1 over several orders of magnitude. When several functions had to be optimized, c was taken as the product of the overlaps of individual functions, as 1 would still represent a perfect match, which will decrease if the description of any of the functions deteriorates. By scanning the value of this error function in the phase space defined by the higher and lower exponents for a definite basis set size, we were able to derive a consistent set of parameters that maximized the overlap with respect to the reference functions while presenting a regular progression along the lanthanide series (Table 1). Not unexpectedly, we observe that the addition of an extra f function allows for a tighter highest exponent and a looser smallest exponent in comparison with the original SARC basis. An example of the error function constructed under this approach is presented in Figure 1.

2. CONSTRUCTION OF THE BASIS SETS Following the general features of the original SARC basis set, we began a reconstruction of the basis sets by first keeping the same basis size (23s16p12d6f) and contraction scheme of the previously proposed basis to a final [18s12p9d3f] pattern. As one of the stated purposes is to avoid any significant increase in computational cost, we carefully explored the benefits in chemically relevant properties by minimal additions to the original construction. One of the outcomes of this search was that a sizable improvement in CASSCF energies and spin−orbit coupling (SOC) properties, as well as a uniformly better performance for higher ionization energies, is obtained if the number of f functions is increased. Therefore, we present here a new version of the basis sets, called SARC2, which is based on a (23s16p12d7f)/[18s12p9d4f] contraction (SARC2-QZV). Compared to the original SARC basis sets, this version contains an additional uncontracted f function and reoptimized exponents for all primitives. Results from the SARC2 basis set were considered to be close enough to the basis set limit in terms of relativistic and nonrelativistic excitation energies since larger basis sets resulted in much smaller corrections for these properties relative to the original SARC bases. In order to keep a low number of fitting parameters for the optimization 1149

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156

Article

Journal of Chemical Theory and Computation

SARC2 basis and found that a correction of around 5% can result for the largest polarization sets (i.e., 3g2h and 3g2h1i) (black line in Figure 2). As expected, given their similar

Table 1. Highest Exponents per Angular Momentum for SARC2 Basis Sets (DKH2 Hamiltonian) La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu

ζs,i

ζp,i

ζd,i

ζf,i

3908799.087 4164893.933 4449443.762 4733993.591 4863007.097 5612784.160 5701463.010 6038159.408 6758316.703 6818227.038 7263979.969 8067496.753 8211204.948 8712676.996 9214149.044

9595.61858 10250.85729 10853.28165 11485.00946 12181.13671 12843.87165 13576.79857 14322.67353 15102.28152 15963.66658 16777.00761 17747.42876 18686.84276 19664.07866 20638.81581

435.686930 471.324765 494.626426 513.816029 533.005633 550.824550 567.272782 580.979641 594.686501 604.281302 612.505418 616.617476 619.358848 619.358848 627.582964

36.296483 56.040299 59.487632 63.561753 67.322479 70.769812 75.157327 78.918054 81.425205 82.678781 87.066295 89.886840 96.781506 95.841325 99.915445

Figure 2. Differences in excitation energies of the lanthanide trihalides (%) between the SARC2-QZV basis and the same basis with sets of polarization functions of different size (black squares). Differences in excitation energies (%) between two consecutive basis sets (red circles). We considered the excitation energies of all possible highest multiplicity states stemming from CASSCF(n,7) calculations, where the active space encompasses the 4f shell. All excitations were equally weighted.

correction with respect to SARC2, the 3g2h and 3g2h1i polarization sets yield similar excitation energies, with an average deviation of 0.16% between them (last point in red line of Figure 2). The maximum deviation between these two orbital sets in terms of excitation energy is 31 cm−1 in an excitation of 35,000 cm−1, while the largest percentage deviation corresponds to 4% (10 cm−1 in an excitation of 230 cm−1). Although these deviations are small, they can still be considered as outliers given that the vast majority of excitation energies show deviations well below 1%. Given the similarity between 3g2h and 3g2h1i, we consider the former to be already converged and choose the smaller set (3g2h) as the polarization functions most appropriate for dynamically correlated calculations with the SARC2 basis set (SARC2-QZVP). The performance of the polarization set was compared with results from the ANO-RCC basis.6 We considered the f−f and f−d excitations of PrIIICl3 as a test case. The ANO basis for ErIII is composed by 301 primitives and was completely decontracted to a size of 25s22p15d11f4g2h. CASSCF results are similar for both calculations, with a mean average deviation of 140.8 cm−1 in excitation energies that reduces to 18.2 cm−1 if we consider only excitations below 1000 cm−1, as the deviations are fairly proportional to the excitation energies, which cover a range from a few tens of cm−1 to 60,000 cm−1 (Table S13). Interestingly, NEVPT2 corrections for d−f excitations are larger in SARC2-QZVP (30%) than in ANO-RCC (23%). It must be noted that the exponents for the polarization functions of SARC2-QZVP (g: 19.1, 4.75, 1.18 and h: 10.1, 3.8) are significantly larger than the corresponding ANO values (g: 2.24, 1.13, 0.45, 0.18 and h: 1.36, 0.54) but closer to the exponents in the ECP-based def2-TZVP basis set (g: from 19.87 to 0.24 and h: 1.4).5 The above observations are

Figure 1. Exponent scan of the error function log(1−c). The presented example corresponds to the Erbium DKH2 basis, where the f function set is optimized considering seven primitives with respect to their minimum and maximum exponents.

After the optimization of the exponents, contraction coefficients were obtained from SAHF calculations of the trivalent ions. We chose this electronic configuration as the most relevant for chemical applications involving lanthanide complexes and to avoid any bias that could be induced by the lack of regularity in electronic configurations of other oxidation states (0, + 1, + 2) along the lanthanide series. 2.1. Polarization Functions for Correlated Methods. For the development of the polarization functions sets, we considered the LnX3 trihalides as adequate systems for exponent optimization. Using geometry optimized structures for the trichlorides, a set of polarization functions for the SARC2 basis set was developed for the inclusion of dynamical correlation treatments, such as NEVPT2. For this purpose, several polarization function sets of increasing size were fitted to maximize the NEVPT2 correction to state energies on top of CASSCF(n,7) calculations, with the 4f shell as the active space. We considered six different polarization sets in a decontracted form (1g, 2g, 2g1h, 3g1h, 3g2h, and 3g2h1i). Once the function coefficients were optimized, we compared the excitation energies of the lanthanide trihalides with respect to the bare 1150

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156

Article

Journal of Chemical Theory and Computation

comparison to light elements. We also tried to build smaller basis sets by removing the tighter exponents for each angular momentum from l = 0 to l = 6. Although acceptable excitation energies were still possible to attain with smaller auxiliary basis sets, we decided to stick to the larger set because differences in total energies were significantly higher and convergence to the desired state was delicate in some cases when employing the smaller basis. In the case of the ORCA program,52 various levels of RI type approximations are available for CASSCF and NEVPT2 calculations. First, it is possible to not fit the Fock operator but only the integrals necessary for the CASSCF gradient and also for NEVPT2. A further speedup can be obtained by fitting both the Coulomb (J) and exchange (K) part of the Fock operator (RI-JK) or evaluating the exchange term numerically by the chain of spheres approximation (RIJCOSX). The latter options are also available for Hartree−Fock and hybrid functionals, while pure DFT functionals can take advantage from Coulomb fitting (RI-J). The SARC2-QZVP/JK basis (in conjunction with SARC2-QZVP) was tested in terms of excitation energies of the lanthanide trichlorides under the RIJK and RIJCOSX approximations (Tables S13 and S14). In all cases, the resolution of the identity approximation was applied both for the Fock operator and the gradients. Deviations in CASSCF excitation energies under the RI-JK approximation remained below 6 cm−1, with an average of 0.29 cm−1. Results for RIJCOSX were also satisfactory, although slightly less accurate, with an average deviation of 1.8 cm−1 and a maximum of 10.7 cm−1. The SARC2-QZV/JK was also tested for CASSCF excitation energies and performed similarly well. Considering the impact of polarization functions for the description of NEVPT2 excitation energies (Section 2.1), it is important to assess the degree of accuracy that can be achieved using the SARC2-QZVP//SARC2-QZVP/JK basis. In comparison with CASSCF, deviations on NEVPT2 energies were higher, with an average of 29 cm−1 and a maximum of 110 cm−1 for RI-JK (RIJCOSX average and maximum deviations were 30 and 113 cm−1, respectively). The errors were proportional to the magnitude of the excitation energies, which range from some tens to several thousand wavenumbers. To complement the error calculation based on the absolute value of the excitation energy, we calculated the percentage deviation for every energy difference. For RI-JK, the average deviation is 0.18%, which is reasonable for general purposes, with a maximum of 0.39%. In the case of the RIJCOSX approximation, the average deviation was 0.51% and the maximum 7.5%; due to the presence of 5−10 cm−1 deviations on top of very low energy excitations ranging from 50 to 250 cm−1, higher excitation energies presented deviations well below 1%. This degree of accuracy should be adequate for studies of excitation energies between terms (i.e., optical spectroscopy), although it could be not sufficient for situations where very accurate estimations (i.e., on the order of 1 cm−1) for the first excitation energies are required. The performance of RI approximations in the context of DFT calculations was also tested through bond energy calculations of the LnCl3 models (Table S16). We employed the hybrid PBE0 functional to test the RI-JK approximations. A similar trend as that observed for CASSCF calculations was encountered in this case, as the RI-JK method yielded very small deviations (average 0.14 kcal/mol) with respect to the non-RI calculations.

presumably related to the fact that the SARC2 exponents were optimized for NEVPT2 correlation energies without imposing any frozen core approximation; hence, they are likely recovering a part of (semi)core correlation. To test both sets of polarization exponents, additional NEVPT2 calculations were performed using the SARC2-QZV basis and defining a polarization set consisting on the g and h primitives of both the ANO-RCC or the SARC2-QZVP. NEVPT2 excitation energies for this basis resemble more the SARC2-QZVP results (MAE = 148.8 cm−1) than the decontracted ANO-RCC (3068 cm−1). It is interesting to note that ANO-RCC and SARC2-QZVP perform similarly with respect to the basis set including either polarization functions for excitations below 1000 cm−1 (MAE = 18.2 and 19.4 cm−1 for SARC2-QVZP and ANO-RCC, respectively). Larger excitations, related with f−d transitions, concentrate most of the deviation observed for ANO-RCC. 2.2. Auxiliary Basis Sets. Resolution of the identity (RI), or density-fitting, techniques are widely employed as an effective approach to speed up electronic structure calculations without compromising accuracy.43−51 This approximation is exploited to simplify the evaluation of the Coulomb part (RI-J), exchange term (RI-K), and also post-Hartree−Fock methods (RI-C). Depending on the exact field of application and method, different accuracies for the auxiliary basis are required. For instance, pure density functionals can be approximated by fitting the Coulomb part by a relatively small basis (/J type), while in hybrid functionals or Hartree−Fock methods, both the exchange and Coulomb term can be fitted, requiring a comparatively larger auxiliary basis (usually denoted/JK). Finally, dedicated auxiliary basis sets for correlated methods can be developed (/C type). As our goal is to provide a fitting basis set that performs adequately for routine applications of wave function methods, we generated two different auxiliary basis sets (SARC2-QZV/JK and SARC2-QZVP/JK) intended to be employed in conjunction with SARC2-QZV and SARC2QZVP, respectively. Although not specifically tailored for correlation methods, both auxiliary basis perform accurately with the corresponding versions of the SARC2 basis for NEVPT2 calculations. For the development of the auxiliary basis set, we first generated the full product basis for all possible angular momentum values up to l = 7 and then grouped the functions according to their exponent values. Blocks of primitives were identified as sets of functions in which the ratio between consecutive exponents is below a given threshold (in this case, 2.0). Then, function blocks were replaced with an eventempered sequence with a ratio between consecutive functions equal to the aforementioned threshold. This procedure led to an auxiliary basis of around 850 primitives for SARC2-QZV/JK. The same procedure resulted in an auxiliary basis close to 1000 primitives for SARC2-QZVP/JK (patterns for the auxiliary basis are presented in the Supporting Information). The main difference between SARC2-QZV/JK and SARC2-QZVP/JK is the presence of l = 7 functions in the latter, which are important for accurate excitation energies for RI-NEVPT2. In this way, the auxiliary basis were 4.7 (for SARC2-QZVP/JK) and 5.6 (for SARC2-QZV/JK) times larger than the size of the main basis, which is adequate for an RI-JK auxiliary basis set (presented in the Supporting Information for both ZORA and DKH2 Hamiltonians). It has to be considered that all-electron basis sets for heavy atoms span a very broad exponent range and present higher angular momentum functions, extending the function space that has to be covered by the product basis in 1151

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156

Article

Journal of Chemical Theory and Computation

Table 2. Radial Expectation Values of Innermost Orbitals for Each Angular Momentum from Spin-Average DKH2-ROHF Calculations for +3 Oxidation State (Bohr); Comparison of Basis Set Limit (UGBS+8) and SARC2 Results UGBS+8 La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu MAE

SARC2

< rs>

< rp>

< rd>

< r f>

< rs>

< rp>

< rd>

< rf>

0.024394 0.023889 0.023399 0.022925 0.022465 0.022018 0.021585 0.021163 0.020753 0.020355 0.019967 0.019589 0.019221 0.018862 0.018512 −

0.093930 0.092036 0.090209 0.088445 0.086741 0.085094 0.083501 0.081959 0.080466 0.079020 0.077617 0.076257 0.074937 0.073656 0.072411 −

0.259119 0.252828 0.246843 0.241142 0.235704 0.230509 0.225542 0.220788 0.216228 0.211854 0.207654 0.203617 0.199734 0.195996 0.192394 −

− 1.018704 0.972202 0.932320 0.897360 0.866211 0.838130 0.812563 0.795743 0.778491 0.761302 0.744471 0.728159 0.712422 0.697308 −

0.024395 0.023890 0.023400 0.022926 0.022466 0.022019 0.021586 0.021164 0.020754 0.020356 0.019968 0.019590 0.019222 0.018863 0.018513 0.000001

0.093945 0.092051 0.090224 0.088461 0.086757 0.085110 0.083517 0.081975 0.080482 0.079036 0.077634 0.076274 0.074954 0.073673 0.072428 0.000016

0.259178 0.252882 0.246897 0.241197 0.235761 0.230569 0.225605 0.220853 0.216296 0.211925 0.207728 0.203695 0.199815 0.196079 0.192476 0.000067

− 1.018708 0.972133 0.932163 0.897153 0.865969 0.837844 0.812248 0.795424 0.778180 0.760961 0.744135 0.727759 0.712083 0.696964 0.000262

2.3. Computational Details. All electronic structure calculations were performed using the ORCA 3.0.3 program.52 For the UGBS+8 reference orbitals, spin average restricted open shell Hartree−Fock53 calculations were converged for four different configurations, as indicated in the text. Once the reference orbitals were obtained, no more calculations were required for the fitting of the SARC2-QZV exponents, as the optimal exponents are obtained to maximize the overlap between the target basis and the reference orbitals. These calculations are very fast and allowed for an initial parameter scan that comprised several orders of magnitude in the exponent values. In the first step, we considered deliberately too large and small boundary values, taking as reference the exponents of the reference UGBS+8 basis. Once the relevant regions of the parameter space were located, more accurate scans were performed to restrict the location of the minimum, while the final step comprised a multivariate optimization. Ionization potentials were obtained from CASSCF calculations considering the same electronic configurations as ref 33. The relativistically recontracted versions of the def2-TZVP basis set were employed for all LnX3 calculations,54 while tests related to the performance of RI approximations were carried out with the auxiliary basis sets def2-TZVP/JK, as def2-TZVP/C led to larger and more variable deviations in total energies, as expected from the better description of core functions in the former. Results for RIJCOSX calculations were dependent on the choice of the numerical grid for the exchange contribution. To assess the performance of the auxiliary basis set without any additional error associated with numerical integration, we set very accurate radial and angular grids and disabled the overlap fitting method (keywords IntAccX 10.0, Grid_x 7, and UseSFitting false in the %method block of ORCA). In order to obtain a single parameter to account for the strength of the SOC interaction, the ζ parameter was fitted to reproduce the matrix elements of the QDPT matrix in an analogous procedure to the ab initio Ligand Field approach for the calculation of ligand field parameters.55 Dipole moment calculations of a set of small lanthanide containing molecules (BP86 functional)56,57 were performed with SARC2-QZV, SARC2-QZVP, decontracted ANO-RCC, and the reference basis employed by Gulde et al.5 For consistency with the

original paper, the Cl def2-QZVPP basis was employed. Integration grids for all DFT calculations were set tighter to avoid any dependence in numerical settings (IntAcc 10.0 and Grid7 in ORCA).

3. ASSESSMENT OF PROPERTIES 3.1. Radial Expectation Values. The correct reproduction of the spatial extent of atomic orbitals compared with reference, basis-set-limit values, is a strong indication of the construction quality for a basis set of reduced size. Table 2 presents the radial expectation values for the 1s, 2p, 3d, and 4f orbitals from spinaveraged Hartree−Fock calculations for the chemically relevant +3 oxidation state of all lanthanides. It is clearly shown that the SARC2 basis sets yield similar values in comparison with our reference UGBS+8 basis for the 1s, 2p, and 3d orbitals for all elements, with slightly larger deviation for the 4f orbitals. However, this deviation represents an enormous improvement compared to the previous SARC basis sets, which display a mean absolute deviation of 0.017245 Bohr compared to the UGBS+8 reference (Table S3, Supporting Information, for details). This improvement from 1.96% to 0.03% in the average mean absolute deviation with respect to UGBS+8 results for the 4f shell confirms that the procedure followed in the present work results in a much more accurate description of the important f space than what was achieved previously. 3.2. Total Energies. As basis set optimization with respect to total energies tends to produce a better representation of the core electrons due to the higher (in absolute number) energies of inner orbitals, there is not an a priori bias toward core or valence orbitals on our approach, as it relies in orbital shape rather than energy. We found that equal weighting resulted in total energies, which were several Hartrees above the basis set limit due to the optimization of too small tighter exponents. Therefore, we decided to first optimize the tightest exponent to maximize the overlap of the innermost function of each angular momentum with respect to their corresponding reference shapes, while keeping a smaller exponent close to the UGBS limit. In this way, the tightest exponent is optimized to represent the tightest function, and then the cutoff exponent is optimized with respect to all occupied orbitals. This procedure led to total energy differences with respect to the basis set limit 1152

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156

Article

Journal of Chemical Theory and Computation

energy from the basis set limit related with the use of SARC2 mostly affect the chemically inert core electrons. If we compare the new results with the original SARC basis set, we observe that deviations from reference energies decrease also in the valence part, indicating that the SARC2 basis provides a better description of both core and valence electrons than the former basis set. The eigenvalues of the outermost shells for every angular momentum are significantly better described in the SARC2 case, in particular for the 4d (with a deviation of 2.2 mEh in SARC2 in comparison to 28.4 mEh for SARC) and 4f (1.6 mEh for SARC2 and 5.3 mEh in SARC). 3.4. Ionization Energies. The original SARC basis sets already performed very well for the first two ionization potentials, with the MAD increasing significantly for the third and becoming disconcertingly large for the fourth. The SARC2 basis sets yield significantly improved results that almost completely remove the deterioration observed previously on progressing from IP1 to IP4 (Table 5). The most impressive improvements, by an order of magnitude, are obtained for IP3 and IP4. These results are equivalent or better in terms of MAD with the corresponding values obtained with the Cologne all-

that were around 1 Eh higher than the basis set limit, significantly smaller than the original SARC energy differences (Table 3). We consider these results strongly improved, but we Table 3. Spin-Average DKH2-ROHF Energies (Hartree) for Neutral Ground State of Lanthanides Using Different Basis Sets UGBS+8 La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu MAD

SARC

SARC2

E

E

Error

E

Error

−8486.5228 −8853.3056 −9229.7463 −9616.0150 −10012.1851 −10418.4433 −10834.9461 −11261.6990 −11698.6279 −12146.3735 −12604.8661 −13074.2524 −13554.7321 −14046.4966 −14549.6201

−8485.4450 −8852.1031 −9228.3766 −9614.4994 −10010.5150 −10416.6038 −10832.9266 −11259.5487 −11696.1826 −12143.6954 −12601.9350 −13071.0477 −13551.2398 −14042.6904 −14545.6044

1.0779 1.2025 1.3698 1.5156 1.6701 1.8395 2.0196 2.1502 2.4453 2.6780 2.9311 3.2047 3.4923 3.8062 4.0158 2.3612

−8485.9666 −8852.7212 −9229.1208 −9615.3427 −10011.4612 −10417.6702 −10834.1093 −11260.7943 −11697.6570 −12145.3178 −12603.7228 −13073.0176 −13553.3888 −14045.0277 −14548.0263

0.5562 0.5845 0.6255 0.6722 0.7239 0.7731 0.8368 0.9047 0.9708 1.0557 1.1433 1.2348 1.3433 1.4689 1.5938 0.9658

Table 5. First (IP1) to Fourth (IP4) Ionization Potentials (eV) Calculated at CASSCF Level for Different Basis Sets (DKH2 Hamiltonian) IP1

would also like to point out that total energies are not a safe indicator for the performance of a basis set in chemically relevant applications. For example, since we are interested among others in an accurate basis set for the representation of excitation energies of lanthanide complexes with multireference treatments, relative energies would be a more crucial parameter than total energies in order to evaluate the performance of the basis set. 3.3. Orbital Eigenvalues. In order to evaluate how these differences in total energy would affect calculated properties, we compared the orbital energies obtained from SARC2 and SARC with respect to the large UGBS+8 basis (Table 4). It can be observed that higher deviations in orbital energies are associated with core orbitals, smoothly decreasing through valence shells, indicating that most of the differences in total

La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu MAD

Table 4. Mean Absolute Deviations for Orbital Energies (mEh) with Respect to UGBS+8 Basis Set for the Neutral Lanthanide Atoms from DKH2 CASSCF(n,7) Calculations, where Active Space Comprises the 4f Shell

1s 2s 2p 3s 3p 4s 3d 4p 5s 4d 5p 6s 4f

SARC

SARC2

288.7 81.6 98.0 104.3 94.9 52.2 93.4 38.3 8.9 28.4 3.5 1.5 5.3

16.1 23.8 9.3 10.9 3.4 2.9 3.6 7.4 2.6 2.2 2.2 1.3 1.6

La Ce Pr Nd Pm Sm Eu Gd Tb Dy Ho Er Tm Yb Lu MAD 1153

IP2

UGBS+8

SARC

SARC2

UGBS+8

SARC

SARC2

4.38 4.56 4.43 4.47 4.52 4.56 4.59 4.79 4.76 4.83 4.91 4.99 5.06 5.14 4.35 −

4.40 4.57 4.44 4.48 4.52 4.56 4.60 4.80 4.76 4.84 4.92 5.00 5.07 5.15 4.34 0.009 IP3

4.37 4.56 4.43 4.47 4.51 4.55 4.59 4.79 4.75 4.83 4.90 4.98 5.05 5.13 4.35 0.006

10.36 11.15 9.89 10.04 10.19 10.34 10.48 11.35 10.71 10.82 10.93 11.04 11.16 11.26 12.59 −

10.34 11.47 9.91 10.06 10.21 10.36 10.50 11.37 10.73 10.84 10.95 11.07 11.18 11.28 12.61 0.039 IP4

10.36 11.14 9.88 10.03 10.18 10.33 10.47 11.35 10.70 10.81 10.92 11.03 11.14 11.25 12.59 0.009

UGBS+8

SARC

SARC2

UGBS+8

SARC

SARC2

18.10 18.24 19.68 19.80 19.83 21.05 22.72 19.53 17.62 18.84 18.50 18.10 19.22 20.92 19.94 −

18.07 17.97 19.39 19.67 19.80 21.07 22.78 19.54 17.84 19.07 18.74 18.35 19.46 21.16 19.97 0.152

18.10 18.25 19.69 19.82 19.85 21.07 22.75 19.54 17.64 18.87 18.52 18.13 19.25 20.95 19.94 0.020

49.89 34.90 37.23 38.93 39.16 39.31 40.81 42.83 35.25 37.47 38.95 38.73 38.46 39.84 41.84 −

49.89 34.30 36.76 38.82 39.40 39.79 41.42 42.93 36.45 38.62 40.10 39.92 39.66 41.01 42.18 0.667

49.89 34.92 37.26 38.96 39.21 39.36 40.86 42.89 35.31 37.53 39.01 38.80 38.54 39.91 41.93 0.052

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156

Article

Journal of Chemical Theory and Computation

calculations, specifically CASSCF and NEVPT2. Despite the lack of a ligand field, CASSCF calculations on free ions provide a useful benchmark for the evaluation of spin−orbit coupling and multiplet splitting, which are representative of more complex systems given the low sensitivity of lanthanide systems to the ligand field. The consistency of these results will be further corroborated in the following with respect to LnCl3 calculations. Table 6 presents the excitation energies for all high-spin multiplets emerging from 4f configurations. It can be observed

electron basis sets (Table S6, Supporting Information). The significantly improved description of the fourth ionization potential of the SARC2 basis in comparison to the original SARC is evidently related with the greatly improved description of the f space by inclusion of an extra f function to some extent and the new optimized f exponents. It can be observed that the tighter f functions were necessary to give flexibility to the description of f electrons for different chemical environments. 3.5. Optimized Geometries and Dipole Moments. Although the development of the SARC2 basis set was mainly motivated by the current interest for multireference calculations on lanthanide systems, it should be confirmed that the basis set can also perform accurately for typical DFT calculations. Therefore, we evaluated the new SARC2 basis in terms of optimized geometric parameters (bond distances and angles) and dipole moments. Optimized geometries were obtained using the PBE0 functional58 and were compared with geometries obtained with the Cologne basis sets for the lanthanides,33 which were previously reported to compare favorably with reference data.59,60 We performed optimizations with both lanthanide basis sets using the exact same settings to avoid any effects arising from different basis sets on the halides or from other potential technical differences between electronic structure packages. It is worth noting that we did not impose any symmetry constraints in the geometry optimizations, and hence, significant angle deviations from C3v symmetry were often observed. We adopt this approach as a more practical and realistic test for geometry optimizations, being more representative of actual DFT practice than constrained optimizations. Bond distances and angles for LnX3, (X = F, Cl, Br and I) are presented in the Supporting Information and indicate a mean average error on bond distances lower than 0.0011 Å and angles that were in agreement by less than 1° in most cases, with an average of 0.46° (Table S9). Therefore, it is clear that SARC2 compares well with other basis sets with respect to reference values of LnX3 systems for the optimization of molecular structures at the DFT level. A good benchmark set for the assessment of the performance of the basis set with respect to DFT-derived molecular properties is the list of dipole moments for small lanthanide molecules presented by Gulde et al.5 as part of the evaluation of the def2-nZVP/PP basis sets, which are designed to work in conjunction with small-core effective core potentials (ECPs). The small molecule set contains lanthanide molecules with coordination numbers 1 to 3, which are mainly hydrides, oxides, and halides. In this way, we can test our basis with respect to a property which is fairly sensitive to the calculated charge distribution and check if the SARC2 basis compares well with current ECP basis sets. As the basis presented here is intended to be employed in conjunction with the DKH2 Hamiltonian, we performed again the reference calculations with the larger basis presented by Gulde et al. to have comparable reference values. The mean absolute deviation in the dipole moments was 0.050 and 0.055 D for SARC2-QZVP and SARC2-QZV, respectively, reasonably comparing to the numbers obtained for the def2 bases (0.035 D for def2-QZVP and 0.071 for def2-TZVP). The complete list of calculated dipole moments is presented Table S10 of the Supporting Information. 3.6. CASSCF Excitation Energies. After accounting for the main limitations of the original SARC basis and its performance toward DFT calculations, we turn our attention to the performance of the SARC2 basis with respect to wave function

Table 6. CASSCF Excitation Energies (cm−1) for Highest Multiplicity Terms of Lanthanide Trivalent Free Ions Pr

Nd

Pm

Sm

Dy

Ho

Er

Tm

MAD

3

H 3 F 3 P 4 I 4 G 4 F, 4S 4 D 5 I 5 G 5 F, 5S 5 D 6 H 6 F 6 P 6 H 6 F 6 P 5 I 5 G 5 F, 5S 5 D 4 I 4 G 4 F, 4S 4 D 3 H 3 F 3 P

UGBS+8

SARC-TZV

SARC2-QZV

0.0 5981.1 27911.9 0.0 14577.2 22907.0 37484.1 0.0 15164.7 23830.3 38995.1 0.0 6739.8 31452.6 0.0 7499.1 34995.6 0.0 17887.5 28109.0 45996.5 0.0 18289.8 28741.0 47030.8 0.0 8013.7 37397.4

0.0 5719.2 26689.8 0.0 14051.0 22080.1 36131.1 0.0 14700.8 23101.3 37802.2 0.0 6561.6 30620.8 0.0 7369.2 34389.8 0.0 17615.1 27680.9 45296.0 0.0 18041.7 28351.3 46393.1 0.0 7916.4 36943.1 577.5

0.0 5982.4 27917.7 0.0 14581.6 22913.9 37495.5 0.0 15170.3 23839.1 39009.5 0.0 6742.7 31465.9 0.0 7503.4 35016.0 0.0 17898.7 28126.5 46025.2 0.0 18301.1 28758.9 47060.1 0.0 8019.4 37423.7 12.4

that the new SARC2 basis yields a large improvement in mean absolute deviations with respect to the reference UGBS+8 basis results to a value of 12.4 cm−1, while the deviation for SARCTZV was 577.5 cm−1. In percentage, SARC2-QZV MAD corresponds to 0.05% with respect to UGBS+8 and SARC-TZV raises to 2.43%. Since excitation energies of several thousands of wavenumbers are involved in this test, the agreement for the SARC2 basis is considered satisfactory in every respect. 3.7. Spin−Orbit Coupling Constants. One of the distinctive aspects of the modeling of systems including lanthanide ions is the need for inclusion of spin−orbit coupling (SOC) effects due to the high atomic number of the lanthanides and the presence of unquenched angular momentum in the ground state in most cases. Bond distances and vibrational frequencies can be reasonably modeled by the consideration of scalar relativistic effects by the ZORA or DKH2 Hamiltonians. Scalar relativistic effects require basis sets that can properly describe the changes in orbital shape under 1154

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156

Article

Journal of Chemical Theory and Computation

SARC basis significantly underestimates this effect due to insufficient flexibility for the description of f orbitals, while SARC2 is close to the basis set limit. An even closer value requires the inclusion of extra f functions, achieving a slow convergence toward UGBS+8 results. For instance, the difference in the ζ parameter between UGBS+8 and SARC2 for the case of Tm3+ is 22.7 cm−1. The inclusion of an eighth f function reduces the error to 10.9 cm−1, and the deviation with a ninth f set is 2.7 cm−1. To put these numbers in context, it must be noted that the error of the original SARC is 287 cm−1, while a completely reoptimized basis with six f functions yields a deviation of 73.5 cm−1. The SARC2 average deviation for all elements in Table 7 is 16.6 cm−1, combining a very good description of SOC with a small basis set size.

these Hamiltonians, affecting mostly the tighter basis functions accounting for the innermost orbitals. On the other hand, magnetic and spectroscopic properties require the explicit description of state mixing due to the spin−orbit coupling operator. It is then expected that the strength of the spin−orbit coupling will require an accurate account of the shape of the orbitals responsible for the spin−orbit coupling mixing. It is important to stress that state interaction due to spin−orbit coupling has a larger effect than ligand field in the energy splitting of electronic levels. In this way, an accurate description of SOC is crucial for a reasonable representation of spectroscopic properties. Furthermore, the magnetic properties of most lanthanide complexes cannot be described without considering SOC effects, even at a qualitative level. Although variational approaches for the inclusion of spin−orbit coupling in multireference calculations have been proposed and tested,61 most current methodologies based on wave function methods include SOC effects by means of quasi-degenerate perturbation theory (QDPT) in a two-step procedure. First, a regular CASSCF calculation is performed to obtain wave functions and state energies that serve as the basis and diagonal matrix elements of the state interaction matrix, respectively. Then, the one-electron SOC operator is calculated between all pairs of CASSCF states. In the QDPT matrix, nonrelativistic states are split in their Ms components according to their multiplicity, and the matrix elements connecting different microstates corresponding to a pair of CASSCF states only differ by their corresponding Clebsch−Gordan coefficients. In this way, the matrix elements of the SOC operator can be fitted to a model matrix in a procedure analogous to the ab initio Ligand Field (AILF) method.55 A natural choice for the SOC operator in the model matrix is

Ĥ = ζ ∑ li·si

4. CONCLUSIONS An entirely redesigned version of the SARC basis sets for lanthanides is presented, with exponents optimized according to a new orbital shape fitting procedure and a significantly improved f shell. The new SARC2 basis sets improve upon the performance of the first-generation basis sets for DFT calculations and safely extend the applicability of the basis sets to wave function-based approaches such as CASSCF and CASPT2 or NEVPT2. According to the presented benchmark tests for both CASSCF and DFT calculations, the new SARC2QZV and the polarized version SARC2-QZVP perform well in the modeling of lanthanide systems for a wide array of metrics. Importantly, wave function-based calculations including scalar and angular spin−orbit effects can be accurately carried out to explore spectroscopic and magnetic properties of lanthanide systems using the SARC2 basis sets. Resolution of the identity (RI) approximations can be accurately employed for CASSCF and NEVPT2 using the reported auxiliary basis sets (SARC2QZV/JK and SARC2-QZVP/JK). In terms of the basis building strategy, exponent optimization based on the radial distribution function shape was successful in producing a balanced basis set in the description of both core and valence regions, avoiding the complications associated with total energy or target property optimized basis sets. The SARC2 basis sets thus represent a reliable choice for multiconfigurational calculations of magnetic and spectroscopic properties of lanthanide systems, as they are able to accurately describe excitation energies, spin− orbit coupling-induced state mixing, and energy corrections associated with dynamical correlation (NEVPT2) in a consistent and reliable way.

(2)

i

where l and s are the one-electron orbital and spin angular moments, and ζ is the spin−orbit coupling parameter. Further details about the representation of the SOC operator go beyond the scope of this paper, and we recommend refs 62 and 63 for the interested reader. Table 7 presents the value of the ζ parameter, which represents the magnitude of the SOC interaction, calculated by the fitting of the matrix elements of the SOC operator in the QDPT matrix, as described earlier. We observe that the original



Table 7. Spin−Orbit Coupling (SOC) ζ Parameter (cm−1) for Trivalent Lanthanide Ions at Different Basis Sets

Ce Pr Nd Pm Sm Eu Tb Dy Ho Er Tm Yb MAD

UGBS+8

SARC-TZV

SARC2-QZV

686.7 809.9 942.9 1086.6 1241.7 1409.1 1768.8 1964.8 2177.5 2406.8 2653.2 2917.3 −

558.0 645.8 773.2 909.8 1056.6 1214.6 1541.9 1724.3 1922.1 2135.8 2365.8 2612.4 217.1

678.2 800.3 932.6 1075.3 1229.3 1395.9 1752.3 1944.8 2156.4 2383.1 2630.5 2887.6 16.6

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b01048. Detailed information about the benchmark tests (PDF) DKH2 and ZORA optimized SARC2-QZV and SARC2QZVP basis for elements from La to Lu, auxiliary SARC2-QZV/JK, and SARC2-QZVP/JK basis sets (TXT)



AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (F.N.). *E-mail: [email protected] (D.A.P.). Notes

The authors declare no competing financial interest. 1155

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156

Article

Journal of Chemical Theory and Computation



(31) Roos, B. O.; Malmqvist, P.-A. Phys. Chem. Chem. Phys. 2004, 6, 2919−2927. (32) Roos, B. O.; Fülscher, M.; Malmqvist, P.-Å.; Merchán, M.; Serrano-Andrés, L. In Quantum Mechanical Electronic Structure Calculations with Chemical Accuracy; Langhoff, S., Ed.; Understanding Chemical Reactivity; Springer: The Netherlands, 1995; Vol. 13, pp 357−438. (33) Dolg, M. J. Chem. Theory Comput. 2011, 7, 3131−3142. (34) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. Chem. Phys. Lett. 2001, 350, 297−305. (35) Angeli, C.; Cimiraglia, R.; Malrieu, J.-P. J. Chem. Phys. 2002, 117, 9138. (36) Angeli, C.; Bories, B.; Cavallini, A.; Cimiraglia, R. J. Chem. Phys. 2006, 124, 054108. (37) Angeli, C.; Cimiraglia, R.; Evangelisti, S.; Leininger, T.; Malrieu, J.-P. J. Chem. Phys. 2001, 114, 10252. (38) Andersson, K.; Malmqvist, P.; Roos, B. O. J. Chem. Phys. 1992, 96, 1218. (39) Roos, B. O.; Andersson, K. Chem. Phys. Lett. 1995, 245, 215− 223. (40) Finley, J.; Malmqvist, P.-Å.; Roos, B. O.; Serrano-Andrés, L. Chem. Phys. Lett. 1998, 288, 299−306. (41) Ghigo, G.; Roos, B. O.; Malmqvist, P.-Å. Chem. Phys. Lett. 2004, 396, 142−149. (42) de Castro, E. V. R.; Jorge, F. E. J. Chem. Phys. 1998, 108, 5225− 5229. (43) Baerends, E. J.; Ellis, D. E.; Ros, P. Chem. Phys. 1973, 2, 41−51. (44) Dunlap, B. I.; Connolly, J. W. D.; Sabin, J. R. J. Chem. Phys. 1979, 71, 3396. (45) Van Alsenoy, C. J. Comput. Chem. 1988, 9, 620−626. (46) Kendall, R. A; Fruchtl, H. A. Theor. Chem. Acc. 1997, 97, 158− 163. (47) Eichkorn, K.; Treutler, O.; Ö hm, H.; Häser, M.; Ahlrichs, R. Chem. Phys. Lett. 1995, 240, 283−290. (48) Eichkorn, K.; Weigend, F.; Treutler, O.; Ahlrichs, R. Theor. Chem. Acc. 1997, 97, 119−124. (49) Whitten, J. L. J. Chem. Phys. 1973, 58, 4496. (50) Neese, F. J. Comput. Chem. 2003, 24, 1740−1747. (51) Laikov, D. N. Chem. Phys. Lett. 1997, 281, 151−156. (52) Neese, F. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 73− 78. (53) Stavrev, K. K.; Zerner, M. C. Int. J. Quantum Chem. 1997, 65 (5), 877−884. (54) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297−3305. (55) Atanasov, M.; Ganyushin, D.; Sivalingam, K.; Neese, F. In Molecular Electronic Structures of Transition Metal Complexes II; Mingos, D. M. P., Day, P., Dahl, J. P., Eds.; Structure and Bonding; Springer: Berlin Heidelberg, 2012; pp 149−220. (56) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (57) Perdew, J. P. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 33, 8822−8824. (58) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. (59) Kovács, A.; Konings, R. J. M. J. Phys. Chem. Ref. Data 2004, 33, 377−404. (60) Hargittai, M. Coord. Chem. Rev. 1988, 91, 35−88. (61) Ganyushin, D.; Neese, F. J. Chem. Phys. 2013, 138, 104113. (62) Atanasov, M.; Aravena, D.; Suturina, E.; Bill, E.; Maganas, D.; Neese, F. Coord. Chem. Rev. 2015, 289−290, 177−214. (63) Neese, F. J. Chem. Phys. 2005, 122, 034107.

ACKNOWLEDGMENTS The authors thank the Max Planck Society for financial support of this work. D.A. thanks CONICYT+PAI “Concurso nacional de apoyo al retorno de investigadores/as desde el extranjero, convocatoria 2014 82140014”. Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02).



REFERENCES

(1) Pantazis, D. A.; Neese, F. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2014, 4, 363−374. (2) Jensen, F. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2013, 3, 273− 295. (3) Peterson, K. A.; Dyall, K. G. In Computational Methods in Lanthanide and Actinide Chemistry; Dolg, M., Ed.; John Wiley & Sons, Ltd., 2015; pp 195−216. (4) Weigend, F. In Computational Methods in Lanthanide and Actinide Chemistry; Dolg, M., Ed.; John Wiley & Sons, Ltd., 2015; pp 181−194. (5) Gulde, R.; Pollak, P.; Weigend, F. J. Chem. Theory Comput. 2012, 8, 4062−4068. (6) Roos, B. O.; Lindh, R.; Malmqvist, P. Å.; Veryazov, V.; Widmark, P. O.; Borin, A. C. J. Phys. Chem. A 2008, 112, 11431−11435. (7) Rösch, N.; Matveev, A.; Nasluzov, V. A.; Neyman, K. M.; Moskaleva, L.; Krüger, S. In Relativistic Electronic Structure Theory Part 2. Applications; Schwerdtfeger, P., Ed.; Elsevier, 2004; Vol. 14, pp 656− 722. (8) Nakajima, T.; Hirao, K. Chem. Rev. 2012, 112, 385−402. (9) Rösch, N.; Krüger, S.; Mayer, M.; Nasluzov, V. A. In Recent Developments and Applications of Modern Density Functional Theory; Seminario, J., Ed.; Elsevier, 1996; Vol. 4, pp 497−566. (10) Reiher, M. Theor. Chem. Acc. 2006, 116, 241−252. (11) Chang, C.; Pelissier, M.; Durand, P. Phys. Scr. 1986, 34, 394− 404. (12) Reiher, M. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2012, 2, 139−149. (13) Saue, T. ChemPhysChem 2011, 12, 3077−3094. (14) van Lenthe, E.; Baerends, E. J.; Snijders, J. G. J. Chem. Phys. 1993, 99, 4597. (15) van Lenthe, E.; Baerends, E. J.; Snijders, J. G. J. Chem. Phys. 1994, 101, 9783. (16) Pantazis, D. A.; Neese, F. J. Chem. Theory Comput. 2009, 5, 2229−2238. (17) DeBeer George, S.; Neese, F. Inorg. Chem. 2010, 49, 1849− 1853. (18) Regueiro-Figueroa, M.; Platas-Iglesias, C. J. Phys. Chem. A 2015, 119, 6436−6445. (19) Kaneko, M.; Miyashita, S.; Nakashima, S. Inorg. Chem. 2015, 54, 7103−7109. (20) Rodríguez-Rodríguez, A.; Esteban-Gómez, D.; De Blas, A.; Rodríguez-Blas, T.; Botta, M.; Tripier, R.; Platas-Iglesias, C. Inorg. Chem. 2012, 51, 13419−13429. (21) Zhang, Y.; Yang, Y.; Jiang, H. J. Phys. Chem. A 2013, 117, 13194−13204. (22) Jones, M. B.; Gaunt, A. J.; Gordon, J. C.; Kaltsoyannis, N.; Neu, M. P.; Scott, B. L. Chem. Sci. 2013, 4, 1189. (23) Singh, S. K.; Rajaraman, G. Dalton Trans. 2013, 42, 3623−3630. (24) Rios, C.; Salcedo, R. Polyhedron 2013, 57, 47−51. (25) Behrle, A. C.; Barnes, C. L.; Kaltsoyannis, N.; Walensky, J. R. Inorg. Chem. 2013, 52, 10623−10631. (26) Deng, Q.; Popov, A. a. Chem. Commun. 2015, 51, 5637−5640. (27) Cimpoesu, F.; Dragoe, N.; Ramanantoanina, H.; Urland, W.; Daul, C. Phys. Chem. Chem. Phys. 2014, 16, 11337−11348. (28) Kaneko, M.; Miyashita, S.; Nakashima, S. Dalton Trans. 2015, 44, 8080−8088. (29) Zhang, Y.; Krylov, D.; Rosenkranz, M.; Schiemenz, S.; Popov, A. A. Chem. Sci. 2015, 6, 2328−2341. (30) Malmqvist, P.-Å.; Roos, B. O. Chem. Phys. Lett. 1989, 155, 189− 194. 1156

DOI: 10.1021/acs.jctc.5b01048 J. Chem. Theory Comput. 2016, 12, 1148−1156