Towards the origin of effective An(III)/Ln(III) separation by tridentate N-donor ligands: a theoretical study on atomic charges and polarisabilities for Cm(III)/Gd(III) separation

ABSTRACT Soft N-donor ligand have shown to separate An(III) from Ln(III). The origin of the selectivity has not been entirely identified, and similar ligands show very different separation qualities. In this study we present a theoretical investigation of several relevant N-donor ligands in terms of atomic charges and polarisabilities obtained from an atoms in molecules approach. These allow new insights into the bonds between the ligands’ nitrogen atoms and the metal cation and explain a major part of the selectivity towards actinide ions. We deduct the superiority of 2,6-bis(1,2,4-triazine-3-yl)pyridines in separation quality compared to similar ligands for the Cm(III)/Gd(III) separation. Furthermore, improvements of existing ligands are developed that allow not only a direct experimental confirmation but also a systematic experimental study of the interactions and their influence on the selectivity.


Introduction
In spent nuclear fuel, the long-term radiotoxicity is dominated by plutonium and the minor actinides. Selective separation of these elements is the first step of the Partitioning and Transmutation (P & T) scheme, a strategy to reduce the long-term radiotoxicity. Highly selective extracting ligands are required for the separation of the trivalent actinides (An(III)) from the fission lanthanides (Ln(III)) due to their high similarity in both chemical properties and ionic radii. It has been shown that high selectivity is achievable using soft Sor N-donor molecules [1][2][3][4]. Numerous SANEX (selective actinide extraction) process options based on Ndonor ligands were tested in order to achieve an efficient An(III)/Ln(III) separation. A variety of such extracting agents has been developed that favour the formation of an An(III) over the Gd(III) complex in solution compared to the respective aquo ions. Among these heteroaromatic nitrogen donor ligands 2,6-bis(1,2,4-triazine-3yl)pyridines (BTPs) were the first extractants to achieve separation factors (SFs) of more than 100 for Am(III) and Eu(III) from acidic solutions without the necessity of adding a synergistic ligand [5,6]. BTP-type ligands form highly symmetric 1:3 complexes [M(BTP) 3 ] 3 + with An(III) and Ln(III) ions (see ref. [4] and refs therein). In addition, some of them (i.e. CA-BTP [7]) meet additional important requirements for an efficient extractant (e.g. reasonably fast extraction kinetics, stability towards high radiation, acidity and high solubility). The reason for the selectivity of the BTP ligands from both theoretical and experimental side still has not been understood entirely. For a further development and optimisation of new ligands, however, it is of utmost importance to gain a deeper understanding on the molecular level. Different spectroscopic methods like extended X-ray absorption fine structure, time-resolved laser fluorescence spectroscopy, X-ray crystallography and nuclear magnetic resonance have been used to study the An(III)/Ln(III) complexes of several ligands [4,[8][9][10][11][12][13][14][15][16][17][18]. Theoretical approaches on this topic have tried to highlight differences in covalency, structure and pH-dependence [19][20][21][22][23][24][25][26][27][28][29][30][31]. The size of the ligands and their metal complexes restrict the possibility to apply wavefunction-based quantum-chemical (QC) methods beyond the MP2 level. On the other hand, molecular dynamics (MD) simulations allowing the investigation of large system sizes are based on classical energy terms. Those often lack the complexity to describe the demanding interactions apparent between the heavy actinides and their surrounding and hence have been restricted to Ln(III) ions [19,20]. In this paper we will study the effect of different side rings (diazene, triazine and tetrazine) as well as the influence of different side chains and their positions on the atomic charges and polarisabilities in the molecule. We use the Hirshfeld method [32] an atoms in molecules (AIM) method to compute accurate atomic data which will explain differences in structure and complex stability.
The main goal of this article is to establish a theory for gas-phase separation that will be the last step towards accurate force-field adjustment for An(III)/Gd(III) ions interacting with N-donor ligands. In future MD simulations based on these force fields this strategy will be transferred to liquid and organic phases.

Molecular structures
For all investigated ligands, the structure was optimised using density functional theory (DFT) with the B3-LYP functional [33,34] as implemented in the TURBOMOLE software package [35]. Basis sets of triple zeta quality have been used for all atoms during the geometry optimisation. To study torsional and side-chain effects separately, all ligands (if possible) have been optimised in C 2 as well as C 2v symmetry. Basis set convergence as well as torsional dependence of charges and polarisabilities were studied on a training set consisting of four ligands: r 2,6-bis(1,2-diazin-3-yl)pyridine (BDP) r 2,6-bis(1,2,4-triazin-3-yl)pyridine (BTP) For these four ligands we considered increasing basisset sizes aug-cc-pVTZ, aug-cc-pVQZ and aug-cc-pV5Z on all atoms, respectively [36]. Using the aug-cc-pVTZ basis set and the optimised C 2v structure of the ligands both N 1 CCN 2 -torsion angles were varied from 0°to 10°a nd 20° (Figure 1). In a second step, we considered several side chains substituting a hydrogen atom at the R 4 , R 5 and  R 6 = ethyl r 2,6-bis(5,6-diethyl-1,2-diazin-3-yl)pyridine (EtBDP56) R 5 = R 6 = ethyl r 2,6-bis(5,6-dipropyl-1,2,4-triazin-3-yl)pyridine (nPrBTP) R 5 = R 6 = n − propyl r 2,6-bis(5,6-diisopropyl-1,2,4-triazin-3-yl)pyridine (iPrBTP) R 5 = R 6 = i − propyl r 2,6-bis(5-methyl-1,2,4-triazin-3-yl)pyridine (mBTP5) R 5 = methyl r 2,6-bis(6-methyl-1,2,4-triazin-3-yl)pyridine (mBTP6) R 6 = methyl r 2,6-bis(5,6-dimethyl-1,2,4-triazin-3-yl)pyridine (mBTP56) R 5 = R 6 = methyl r 2,6-bis(1,2,4-triazin-3-yl)4-methyl-pyridine (mBTPp4) P 4 = methyl r 2,6-bis(3-propyl-1,2,4,5-tetrazin-3-yl)pyridine (nPrBQP) R 6 = ethyl For the methyl-substituted BTP-ligands mBTP5, mBTP6, mBTP56 and mBTPp4 as well as the ethylsubstituted BDP-ligands EtBDP5, EtBDP6 and EtBDP56, we also optimised the corresponding 1:3 complex with Cm(III) and Gd(III). For both ions, small-core relativistic pseudopotentials with corresponding basissets (Cm/ECP60MWB and Gd/ECP28MWB [37]) were employed. Due to the poor convergence behaviour of the B3-LYP functional for complexes including heavy ions, we chose the BH-LYP functional [38] for the geometry optimisations. On the resulting structures, a single-point MP2 calculation followed determining the gas-phase interaction energies E g corrected for basis-set superposition errors by the counterpoise method [39]: Based on the difference in interaction energies E g = E(CmL 3 ) − E(GdL 3 ), a first approximative SF for the exchange reaction (3) can be calculated according to Note, that in [40] we computed the difference of the corresponding aquo ions as 17.7 kcal/mol leading to a positive E in the exchange reaction (3) This approximation has to be understood in terms of an energy decomposition. Petit et al. [29] suggest a division of the total bonding energy E bonding into parts according to As we will focus on the latter two parts in the gas phase, we will not consider any solvent molecules nor study temperature dependence or entropy effects. These effects should be treated by MD techniques in a fully solvated system. Furthermore, as already pointed out by Bryantsev and Hay [24], zero-point and thermal corrections only lead to small changes in G. We have determined the corrections for some model metal/ligand exchange reaction to be smaller than 0.2 kcal/mol at ambient conditions, which is smaller than the accuracy those values can be calculated with in gas phase. Also it should be noted, that the motivation to study Cm/Gd separation is mostly based on the ions' well-defined 8 S 7/2 ground state that allows the usage of single-reference methods. Although Cm/Gd separation is relevant, most experimental data focusses on Am/Eu and Cm/Eu separation. Panak et al. [4] showed, however, that differences in G for Am/Eu and Cm/Eu separation compared to Cm/Gd separation are rather small for N-donor ligands (< 0.2 kcal/mol). Hence the obtained SFs between Gd(III) and Cm(III) can easily be compared to the experimental results for Am/Eu and Cm/Eu systems. Although there have been theoretical studies concerning these systems, multi-configurational methods would be necessary to describe them properly.

Hirshfeld method
Many approaches have been pursued in order to divide molecular electronic properties into atomic contributions, especially for accurate force-field development. In general, all AIM techniques origin in a division of the space into atomic regions, i.e. the Hirshfeld [32], LoProp [41], Mulliken [42] or Bader basin approaches [43], to name a few. In this study, will we employ the Hirshfeld method to derive charges and dipole polarisabilities for a set of relevant N-donor ligands. This method has already been used successfully to study polarisabilities in aromatic systems [44].
The Hirshfeld charge q i for atom i is defined as using the atomic weighting functions and the difference of molecular density ρ and promolecular density ρ 0 Here, ρ 0 i (r) is the spherical ground-state density of the free atom i. Lu et al. [45] point out that for the carbon atom, the s 2 p 2 ground-state density is not spherical and hence the s 1 p 3 state is used throughout their investigation. Although the charges are not effected much by this approximation, the polarisabilities strongly depend on the choice made here. Hence, we use fractional occupation numbers as implemented in TURBOMOLE to create an s 2 p 2/3 x p 2/3 y p 2/3 z state which has a spherical density distribution and is closer to the actual ground state of the atom. Having defined the Hirshfeld charges, an atomic dipole is computed.
where R γ i is the γ component of the positional vector. From this, the diagonal elements of the atoms' static dipole polarisability tensor are calculated as a limit of the electric field F γ Here, μ γ i (0) represents the dipole moment of atom i without any electric field, e.g. F = 0. The isotropic polarisability α i is obtained by averaging over x-, y-and zcomponents. Later, we will investigate α i as well as its' zcomponent α zi where the z-axis is chosen to be the ligand's C 2 symmetry axis.

Induced dipoles
The perturbation of the electron cloud experiencing an electric field can be described considering a set of induced dipoles {μ i } located on the atomic centres. Given the electric field, they are calculated according to Here, α i are the static dipole polarisabilities corresponding to atom i, E i is the electric field generated by the charges q i and T i is the dipole tensor. Introducing a Thole damping term [46] the electric field and the dipole tensor become: Assuming the polarisabilities to be isotropic and independent of the electric field the total polarisation interaction can be described by the sum This energy is the part of equation (4) that is investigated here. In order to estimate the part of the separation originating from the differences in induced dipoles from both ions and ligands we carried out two short gas-phase MD runs of 10 ps using the POLARIS(MD) software package [47]. Here, we use the same force field for both ions so the only differences are their polarisabilities and masses. For the ligand/ligand interaction, the standard CHARMM parameters were employed [48]. Remaining ion/ligand parameters were adjusted to reproduce the gas-phase metal-nitrogen distances and will only be used for this qualitative analysis. The MD runs were performed

Results and discussion
The nPrBTP and iPrBTP ligands are well-known extracting ligands [4], whereas apart from nPrBQP, EtBDP6 [12] and mBTP56 [49] the other ligands are up to our knowledge purely theoretical model systems. Based on the results of this section, provided similar solubilities of the systems, a series of experiments can be designed to verify changes in complexation strengths directly.

Charges
It has been shown that one of the advantages of the Hirshfeld method is the fast convergence with the size of the basis set [45,50]. Accordingly we found the charges for the ligands to be converged using the aug-cc-pVTZ basis set. As Hirshfeld charges are known to be too small, we will also investigate the atomic dipole moment corrected Hirshfeld charges (ADCH) [45].
The nitrogen charges q(N 1 ) and q(N 2 ) depicted in Figure 2 show a rather weak torsion-angular dependence. With increasing torsional angle from 0°to 20°, q(N 1 ) Table . Atomic Hirshfeld charges (e) (ADCH charges) and isotropic polarisabilities (Å  ) for both N  and N  atoms in the C v optimised structures.  Table . Atomic Hirshfeld charges (e) (ADCH charges) and isotropic polarisabilities (Å  ) for both N  and N  atoms, ring-ring torsion angle θ (rad) and energy difference to the C v structures (kcal/mol) in the C  optimised structures. increases by 5% for all ligands, whereas the q(N 2 ) charges decrease by about 10%. The bond between heavy ions and N-donor ligands is known to be mainly ionic [51]; hence, the charges will determine the coordinating structure and dominate the binding energy. In accordance, we have computed one of the highest binding energies for the Terpy 1:3 complex with Cm(III) and Gd(III) among the investigated ligands in our recent study [40]. As in the case of the Terpy ligand, however, this is not sufficient to yield effective Ln(III)/An(III) separation. The introduction of electron-pushing side chains increases the electron density over the whole aromatic system. Accordingly, the partial charges of the coordinating nitrogen atoms change depending of the substituent's position ( Table 1). Differences of the effect on charges are small for the three possible side-chain sites R 4 , R 5 and R 6 .
We notice a more pronounced change in the ADCH charges. Especially for the C 2 structures the dipole correction sometimes seems to overshoot ( Table 2). This probably origins from the different dipoles of the molecules due to their different torsional angles. Because of this inconvenience we will continue the analysis based on the unscaled Hirshfeld charges.

Polarisabilites
Along the Ln(III) and An(III) series, despite having the same charge, the polarisabilities of the ions change significantly. Induced by the diffuse 5f-orbitals of the An(III) ions, their polarisability is in general larger than the corresponding Ln(III) ion. For example, we found α Cm = 1.165Å 3 [52] and α Gd = 0.86Å 3 .
Polarisability studies generally focus on isotropic values only, but with regard of covalent bonds more focus should be put on the part of the polarisability tensor pointing towards the virtual metal-ion center, here the zaxis (Figure 3). Like Hirshfeld charges, the isotropic polarisabilities converge very fast with increasing basis sets. Computed changes from aug-cc-pVTZ to aug-cc-pV5Z basis sets are less than 2% for the BDP, BTP and BQP ligands. The z-components α z , however, show a slightly less stable behaviour and changes up to 10% were found. Within these uncertainties we will continue the discussion on polarisabilities using the aug-cc-pVTZ basis sets. Figures 4 and 5 highlight the change of the polarisabilities and their z-component with changing torsional angle. First we notice, that α(N 1 ) is significantly smaller compared to α(N 2 ). Whereas the first appear to be angularindependent, the latter change non-monotonically up to 30%. The z-components show a similar behaviour. In particular for the N 1 nitrogen atoms, the z-axis points exactly in the ideal bond-direction; hence, the potential of the BTP pyridine-nitrogen to form a covalent bond is supposed to be slightly higher compared to the other ligands as the higher N 1 z-polarisability allows an easier deformation of the electron density. For the N 2 nitrogen atoms, the bond direction is not as well defined due to different bond lengths and angles in different complexes; hence, an evaluation of the isotropic values is preferred. Again, the BTP  nitrogen atoms exhibit the highest polarisabilities, especially for small torsional angles. It is important to notice that the isotropic molecular polarisabilities decrease with the amount of nitrogen atoms, i.e. for Terpy, BDP, BTP and BQP, we obtain values of 31.6, 29.9, 28.1 and 26.7 A 3 , respectively, from a B3-LYP/aug-cc-pVTZ calculation. Hence, the computation of atomic polarisabilities is mandatory for a detailed investigation of the metalligand interaction.
As already shown for the charges, we find a slight effect on the polarisabilities depending on the position of the substituent. In this case, however, the effect of the R 5 position seems to be overall stronger. Both N 1 and N 2 polarisabilities increase, which is partially compensated when a second side chain is introduced. An alkyl substituent on the R 6 position, although increasing the nitrogen electron density has a negative effect on both α(N 1 ) and α(N 2 ). It would hence be very interesting to investigate a mono-substituted species at the R 5 position experimentally. Such systems have already been synthesised for bis-triazine bis-pyridine (BTBP) based ligands but no comparison to the alternative substitution at R 6 has been drawn [53]. The R 4 position, up to our knowledge, has not been target to substitutions, yet. Although an effect on α(N 1 ) similar to the R 5 substitution is observed in C 2v it vanishes when the ligand is relaxed in C 2 symmetry. This proves that the effect different to R 5 originated by a deformation of the pyridine ring by the close proximity of the side chain due to the symmetry restrictions. For all substituents, we have determined the main effect on the polarisabilities to be on the meta-positions in the aromatic ring.

Separation
Within the 10 ps time-frame of the MD runs, the 1:3 complex was stable although the geometry, in particular the bis-angle, varies with changing polarisabilities. It will be interesting to see how solvent molecules influence the dynamics in future studies, after a proper adjustment of all relevant force fields. Considering classic chargedipole and dipole-dipole interactions via the formula (15) a separation between the two ions is achieved due to the larger interaction of the Cm(III) ion with the ligand. It should hence in principle be possible to assess and improve the separation quality by tuning the polarisabilities of the coordinating nitrogen atoms. The optimal separation area based on the ADCH charges computed for BTP is located around α(N 1 ) = 0.35 and α(N 2 ) = 0.95 ( Figure 6). Using the generally smaller unscaled Hirshfeld charges instead the profile flattens significantly and the minimum is shifted towards α(N 1 ) = 0.31 and α(N 2 ) = 0.65 (see supplemental information Figure 1). Within these uncertainties we continue the discussion based on the ADCH charges. In Figure 6 we also plot the position of some ligands on the contour map. Although solvation  effects are not included it can already be seen that nPrBTP is clearly superior to most other ligands considered. The mBTP5 species is, behind nPrBTP, the next closest to the determined minimum in the contour plot. Its qualities also show in the corresponding 1:3 complexes. The ligand-series mBTP5, mBTP6 and mBTP56 allows a very constructive analysis of the roles that charges and polarisabilities play for the structure and binding energy. Table 3 summarises the optimised distances d 1 and d 2 of the metal ion to the coordinating N 1 and N 2 nitrogen atoms as well as the interaction energies E g for both Cm and Gd complexes.
As predicted by the polarisability analysis, the mBTP5 ligand has better separation qualities in gas phase with a E g of 15.84 kcal/mol compared to mBTP6 (16.58 kcal/mol). Table 3 also highlights the differences in binding energy gained by adding side chains. Again, the R 5 position induces the positive effect in both energy and separation quality. The larger energy gain might be an indicator for covalency induced by the higher polarisability.
The contour map (Figure 6) also suggests a possible improvement on the BDP ligand by substituting at the R 5 rather than at the R 6 position as investigated by Beele et al. [12] Optimised gas-phase complexes confirm an increase of the SF by a factor of 2 from 44 (EtBDP6) to 80 (EtBDP5). It is further interesting, that the substitution at both R 5 and R 6 positions does not lead to a further improvement as in the BTP case as polarisabilities drop out of the optimal separation zone ( Table 4).
Regarding the polarisabilities of nPrBTP, nPrBQP and EtBDP6, we see increasing distances to the optimal separation zone in Figure 6. These correlated well to the experimental results obtained for these ligands [4,12] with SFs of SF(nPr-BTP) = 133-193, SF(nPrBQP) = 9 and SF(EtBDP6) = 5. A direct comparison of the gasphase data for ligands from different families (BDP, BTP, etc.), however, is unadvisable as their interaction with the solvent differs and hence demand a study in solution. Additionally, the SFs obtained by Beele et al. on nPrBQP and EtBDP6 were obtained at very low ligand concentrations (<10 mM), whereas gas-phase calculations correspond to highly saturated ligand solutions.

Conclusions
We have shown that charges and polarisabilities determined by the Hirshfeld method can be connected to separation quality in gas phase. As already proven for Hirshfeld charges, also the polarisabilities converge very fast with increasing basis set and we obtained stable values using the aug-cc-pVTZ basis sets. A rather weak torsionangular dependency of both charges and polarisabilities was determined allowing the usage of those values in MD simulations, where the ligand's structure will change rather quickly due to thermodynamic effects.
This study is an essential step towards the creation of accurate force fields. Corresponding MD simulations will include dynamic and solvent effects in the investigation and help understand the extraction processes. The optimal area of polarisabilities in gas phase was determined by short MD runs and will be recomputed in solution in future studies to include temperature and entropy changes.
We have also shown, that the inclusion of side chains, especially their position and length is of high importance and should not be omitted when performing studies on separation. For example, we have deducted that the performance of the EtBDP ligand can be improved significantly by changing the position of the ethyl residue. Based on the presented results, experiments have been designed to confirm the findings.