Computational Details

In this study, two models have been used to understand the arrangement of the environment of Bk(IV) in carbonate and carbonate-hydroxide solutions. All calculations herein, except ligand-field DFT (LFDFT), were performed using the ORCA 4.1 package.29 For DFT optimizations of the molecular structures, we utilized the scalar relativistic all-electron Douglas-Kroll-Hess Hamiltonian (DKH), the PBE0 hybrid functional, and a SARC-DKH-TZVP basis set for the Bk(IV) atom and def2-TZVP for all the other atoms. No symmetry constraints were used. TD-DFT was used to model the UV-vis spectrum of the two optimized models. Kohn-Sham density functional theory with approximate functionals has proven to be remarkably successful given its relative simplicity and efficiency.  However, it is well known that DFT produces an important delocalization error responsible for overestimation of covalency, which has an impact in the correct prediction of the ground state and spectroscopic properties. These problems can be solved in part using detailed parameterizations principally in the interelectronic terms. A good example of these functionals is the Coulomb attenuating method, CAM-B3LYP functional,30 which was developed to minimize deviations in charge-transfer excitation energies and used in this work to reproduce the absorption spectra of both carbonate compounds. It is important to note the recent study of Yang et al.31, which is the closer DFT study to heavy actinides using this class of functionals. In this context, this work represents one of the first studies where the performance of CAM-B3LYP functional is applied to heavy actinides.
The Complete Active Space Self-Consistent Field (CASSCF) approach combined with NEVPT2 for the inclusion of both, static and dynamic correlation, allows the construction of the configuration of all possible multiplicities within a determined active orbital space. Additionally, state mixing due to spin-orbit (SO) coupling can be accounted for via quasi-degenerate perturbation theory approach (QDPT). The effect of the dynamic correlation is included in the QDPT step as a diagonal correction of the nonrelativistic state energies. This combination of tools has been shown to be an effective methodology in the description of actinide compounds. The CASSCF/NEVPT2 approach provides a reliable treatment of the multiconfigurational character present in open-shell actinide systems and can be used with all-electron relativistic Hamiltonians. The CASSCF calculations were carried out using, initially, a minimal active space, which corresponds to the seven 5f electrons of the Bk(IV) ion in seven 5f orbitals. This active space was then augmented to a CAS(11,9) and finally to a CAS(13,10), including two and three symmetry adapted ligand orbitals respectively, which were selected based on its Bk–O bonding or antibonding character and 5f contribution. The 6d shell of orbitals was also considered. However, the occupation numbers of these orbitals were less than 0.02, so they were omitted in the results and discussion section.
Ligand-field theory was also considered for this study from two perspectives, ab initio (AILFT) and DFT (LFDFT)32. The former was calculated from the CASSCF matrix within ORCA, whereas the latter was performed in ADF201933,34. LFDFT was performed based on Kohn-Sham orbitals calculated at ZORA/PBE0/TZP level of theory.
To obtain a better understanding of the nature of the chemical bond, Natural Bond Orbital (NBO) analyses were performed on correlated wavefunctions to localize orbitals. The NBO analysis is helpful to understand the Lewis picture of bonding by performing a transformation of the wavefunction into a localized form. NBOs are conceived as a set of localized bonds and lone pairs as basic units, where delocalization effects are missing from the molecular orbital (MO) perspective. However, there is an intermediate approach that recovers the MO perspective in the natural localized framework called the natural localized molecular orbitals (NLMOs).35 These calculations were carried out using the Weinhold’s package NBO636 attached to ORCA.
Another tool especially developed to gain a deeper comprehension on the nature of bonding in molecular systems is the Bader’s quantum theory of atoms in molecules (QTAIM).37 QTAIM parameters at the bond critical points (BCPs) such as the electron density ρBCP(r), Laplacian of the electron density ∇ρBCP(r), delocalization δBCP(r) and localization l indices, and total energy densities HBCP(r) are commonly used to characterize the nature of bonds. It is well-known that covalent bonds or open-shell interactions, in terms of QTAIM metrics, present ρBCP(r) > 0.2 a.u. with ∇ρBCP(r) < 0 and HBCP(r) < 0, while ρBCP (r) < 0.1 a.u. with ∇ρ (r) > 0 and HBCP(r) > 0 represent ionic, van der Waals or hydrogen bonds (closed-shell interactions). However, as Kerridge pointed out, there is no a priori  reason to assume this convention for f-element chemistry.13 On the other hand, delocalization and localization indices are useful to quantify the bond order and calculate the partial oxidation states, respectively.38
Additionally, the Interacting Quantum Atom (IQA)27energy partition analysis have been carried out. This complementary tool to understand the chemical bond under the framework of nonoverlapping quantum atoms decomposes the total energy of interaction between two topological atoms into Coulomb and exchange-correlation contributions. These calculations are obtained by 6-dimensional integration over the whole electron density between the two topological atoms, which uses the 2-electron density matrix (2EDM).
QTAIM parameters and IQA energies were calculated using the AIMAll software developed by Todd Keith.39 All parameters reported herein were based on multiconfigurational wavefunctions. For post-HF wavefunctions, AIMAll uses the Muller approximation40 to obtain the 2EDM due to the high computational costs involved in obtaining the exact matrix.