Investigation of the reactions of U, U⁺ and U²⁺ with ammonia: mechanisms and topological analysis

Abstract

Reactions of uranium atom and ions (U, U⁺ and U²⁺) with NH₃ have been studied in the gas phase using density functional calculations. The detailed description of reaction mechanisms and diverse topological analysis provide insights into the chemistry of uranium species. The reaction pathways leading to different products are closely examined. In the case of U and U⁺ reacting with NH₃, the isomerization channel (the products are H₂UNH and H₂UNH⁺, respectively) and H₂ elimination channel (decomposed to UNH + H₂ and UNH⁺ + H₂) were investigated. In the U²⁺ + NH₃ reaction, UNH²⁺ was found to be the only product generated by the dehydrogenation of the HUNH₂²⁺ intermediate. The obtained results are compared to recent experimental studies. Bonding evolution along with the reactions was investigated using diverse topological analyses, including electron localization function, atoms in molecules, and natural bond orbital. A comparison with other results reported in the literature for similar systems, such as the reactions of uranium species with NH₃, H₂O and CH₄, is made to underline the similarities and differences of the reaction mechanisms.

---

1. Introduction

During the two past decades, the reactions of actinide with small molecules in the gas phase have attracted considerable attention both in the theoretical^1–7 and experimental investigations.^8–12 These studies provide insights into the novel electronic structures and chemical properties of the actinide species. Ammonia is widely used as a nitrogen source in the production of metal nitrides and amides,^13–16 hence the reaction of ammonia with actinide can provide insights into the formation mechanisms of actinide imines. It is found that the lone pair electrons of nitrogen plays a vital role in the bonding and structure of these metal imine molecules.¹⁷ In addition, the ammonia N–H bond activation by actinide atom and cations is also an important topic for the fundamental understanding of the activation of prototypical chemical bonds, and is conducive for exploring the influence of the electronic state of actinide species on the reactivity and reaction channel.

Laser-ablated uranium atoms react with ammonia in the gas phase to form H₂N–UH and HN=UH₂,¹⁷ which were detected and characterized by infrared spectroscopic method. The possible reaction channel was proposed as following equation:

U + NH₃ → U–NH₃ → H₂NUH → HN=UH₂ (1)

This experiment also indicated that H₂N–UH and HN=UH₂ have comparable energies but visible light irradiation decreases the energy of H₂N–UH and increases that of HN=UH₂.¹⁷ In addition to uranium, thorium atoms can activate N–H bonds in ammonia¹⁸ to form thorimine (HN=ThH₂). Further analyses indicate that there are similarities between the reactions of Th and U with NH₃.

Experimental investigations mentioned above provide original information on the gas phase reactivity of U atoms. However, limitation exists in this reaction. To obtain detailed information on intermediate mechanisms as well as the nature of the bonding of the complexes, theoretical computation is of very important. On the other hand, quantum chemical and density functional theory (DFT) methods have been proven to be a successful method for the study of actinide compounds.¹⁹ Therefore, theoretical studies can be particularly valuable in providing additional detailed information regarding reaction mechanisms.²⁰ Andrews and co workers studied the bond character of HN=ThH₂ and HN=UH₂ using DFT.^17,18 They revealed a new type of multiple bond from actinide atoms, an additional π bonding between the actinide and nitrogen due to an important contribution from the 5f and 6d orbitals.

To the best of our knowledge, experimental data on the products of the reaction of U cations with NH₃ are scarce, and this is the first mechanistic study of these reactions. A systematic investigation of the reactivity of U atom and cations toward NH₃ should provide a complete mechanistic scheme of this chemistry.

The main focus of the present work is to perform a thorough investigation of the reactions between U atom and cations (U⁺ and U²⁺) with NH₃. The complete reaction mechanism for both high- and low-spin states and molecular structures and energetics properties of all the intermediates are reported. Bonding evolution during the reaction pathways were investigated in terms of diverse analyses, including electron localization function (ELF), atoms in molecules (AIM) and natural bond orbital (NBO). In addition, we compare the similarities and differences of the energy and mechanism of the reaction between uranium atom and cations and NH₃, H₂O and CH₄ to provide insights into the chemistry of the uranium species.

2. Computational details

On the basis of the performance observed in previous DFT studies of the reaction between uranium atom and cations and H₂O, N₂O and CH₄,^{1,3,5,12,21–24} different levels of theory were applied in our study to optimize the geometries of the minima and transition states on the potential energy surface of all the possible spin states. The B3LYP,^25–27 PW91PW91,^28,29 TPSSTPSS,³⁰ PBE0,³¹ PBEPBE³² and BMK³³ functionals were used along with the Stuttgart relativistic effective core potential (RECP)³⁴ for the U atom. This small-core RECP treats 60 electrons as the “core”, named as SDD, leaving the n ≥ 5 shell (5s, 5p, 5d, 5f, 6s, 6p, 6d, and 7s) as the valence electrons. The SDD RECP constitutes a good balance between the relativistic corrections introduced via the RECP and valence electrons.¹ For N and H atoms, the 6-311++G(d,p) basis sets by Pople and co-workers were employed.³⁵ These calculations were carried out with the Gaussian 03 package.³⁶ We ensured that the transition state obtained on each potential energy surface has only one imaginary frequency, and intrinsic reaction coordinate (IRC) calculation proved that each imaginary frequency correctly connects the reactants and products. The vibrational zero point energy (VZPE) corrections were included in all the relative energies. We investigated the expectation value of S² to avoid spin contamination impact on the results, and ensured that the calculated values have low spin contamination.

The investigation of bonding characteristic along the pathway is based on the analysis of the electron localization function (ELF),^37,38 which were performed with the Multiwfn³⁹ package. The topological properties of the (3, −1) bond critical points (bcp) in the gradient field of the electron density were analyzed using the atoms-in-molecules (AIM) theory,⁴⁰ as implemented in the Multiwfn code. The bond order analysis was implemented within the natural bond orbital (NBO) scheme^41,42 to gather detailed information of the evolution of chemical bonds during the reaction.

3. Results and discussion

To evaluate the accuracy of current computational schemes, first, we calculated the relative energies of atomic energy levels as well as the first and second adiabatic ionization potentials (IP) for the U atom. The results are listed in Table 1 along with the corresponding experimental values. As can be seen, results obtained by the approaches used in this study are in reasonable agreement with available experimental data. In addition to the PBE/SDD method, other methods correctly predict the ordering of electronic states. The PW91/SDD and BMK/SDD calculations are in good agreement with the experimental values of the relative energy levels. In the case of ionization potentials, these calculations also agree well with experimental results, especially the B3LYP/SDD method. Although, note that the BMK/SDD underestimates the IP1.

Table 1 Relative energies (in kcal mol⁻¹) of energy levels and adiabatic ionization potentials (in eV) for the U atom

	B3LYP/SDD	PW91/SDD	TPSS/SDD	PBE0/SDD	PBE/SDD	BMK/SDD	Exptl
a Statistically averaged spin orbit energy levels obtained from http://web2.lac.u-psud.fr/lac/Database/Contents.html.b Ref. 44.
Triplet	34.727	10.423	39.786	24.738	0.000	13.531	16.475^a
Quintet	0.000	0.000	0.000	0.000	30.794	0.000	0.000^a
IP1	6.24	5.96	6.14	6.25	6.11	3.45	6.19^b
IP2	11.58	11.21	10.93	11.16	11.17	12.52	11.59^b

---

In addition, we compared the calculated vibrational frequencies of the H₂NUH and HNUH₂ complexes with available experimental data. The results are collated in Table 2, along with the computational values reported by Andrews et al.¹⁷ Our results indicate that the adopted computational methods perform well in frequency calculation. In particular, the B3LYP and PW91 approaches give rise to a better agreement with the experimental values for both the complexes.

Table 2 Comparison of observed values and calculated harmonic frequencies for H₂NUH and HNUH₂^a

	B3LYP	PW91	TPSS	PBE0	PBE	BMK	Other^b	Other^c	Exptl^d
a All calculations used the SDD for U and 6-311++G(d,p) for H and N atoms. All frequencies are in cm^−1.b Ref. 14. Calculation at the B3LYP/TZP/SDD level of theory.c Ref. 14. Calculation at the PW91/TZP/SDD level of theory.d Ref. 14. Observed argon matrix frequencies.
H2NUH	3552.2	3477.6	3481.7	3601.8	3474.0	3639.5	3357.8	3486.1	—
	3466.0	3381.3	3392.9	3508.2	3377.9	3549.4	3468.2	3387.6	—
	1537.8	1481.8	1528.5	1536.5	1477.9	1544.3	1540.9	1484.1	1488.6
	1339.2	1316.9	1350.6	1375.7	1315.1	1382.4	1382.2	1387.8	1349.8
	527.1	528.9	534.7	538.2	528.1	505.9	525.9	532.9	—
	505.9	482.5	523.9	515.1	480.2	480.4	489.4	465.7	508.5
	429.0	416.5	422.6	435.7	416.0	402.5	409.5	392.9	—
	231.2	222.8	327.1	241.3	228.9	235.1	252.0	237.2	—
	194.1	134.4	203.1	200.3	147.2	114.5	159.2	169.1	—
HNUH2	3560.7	3463.5	3470.2	3602.6	3440.7	3620.6	3558.8	3470.4	—
	1441.2	1435.9	1449.3	1473.7	1442.4	1461.7	1465.2	1446.4	1436.3
	1381.8	1394.9	1374.7	1417.5	1412.6	1404.4	1420.4	1418.8	1403.0
	827.3	805.9	796.9	855.6	817.6	879.8	837.9	821.0	819.9
	572.0	523.3	548.6	657.6	532.7	576.9	520.9	514.8	495.7
	516.0	492.5	506.4	531.2	489.9	545.6	483.8	486.1	495.5
	470.3	470.0	478.5	473.9	480.0	480.8	471.6	473.0	465.1
	328.0	268.9	264.8	340.2	261.7	384.3	316.1	286.0	—
	248.3	239.7	198.3	266.1	254.2	378.0	258.1	244.6	—

---

3.1 Reaction mechanisms

3.1.1 Neutral uranium reaction. Intermediate complexes, transition states and product molecules were optimized at different levels of theory. Geometric parameters corresponding to the lowest energy spin states are compared in Fig. 1, and show relatively little dependence on the computational method. The potential energy profiles of the U + NH₃ reaction are presented in Fig. 2. Activation barriers for all the TSs are reported in Table 3.

Fig. 1 Structures and selected geometric parameters of stationary points on the U–NH₃ potential energy surface optimized at the B3LYP/SDD, PW91/SDD, TPSS/SDD, PBE0/SDD, PBE/SDD and BMK/SDD levels of theory (from top to bottom rows, respectively). Bond distances are in Å and angles are in degrees.

Fig. 2 Potential energy profile for the reaction of U + NH₃ computed at the B3LYP/SDD and PW91/SDD (in parentheses) levels of theory.

Table 3 Activation barriers (kcal mol⁻¹) of the transition states (see Fig. 2, 3 and 5)^a

Reaction	TS	B3LYP	PW91	TPSS	PBE0	PBE	BMK
a SDD for U, U^+ and U^2+ and 6-311++G(d,p) for H and O atoms.b Calculated as the energy difference between the first transition state (TS1) and the first complex I.c Calculated as the energy difference between the second transition state (TS2a) and the intermediate II.d Calculated as the energy difference between the third transition state (TS2b) and the intermediate II.
U + NH3	TS1^b	4.10	10.51	8.09	12.03	9.50	20.40
	TS2a^c	30.93	21.87	25.09	31.15	21.01	32.92
	TS2b^d	22.72	15.57	18.85	21.32	13.95	27.37
U^+ + NH3	TS1^b	28.20	20.63	18.34	25.60	20.73	26.67
	TS2a^c	57.57	44.40	43.87	43.20	43.18	64.80
	TS2b^d	24.70	13.10	12.77	6.43	13.34	25.33
U^2+ + NH3	TS1^b	61.26	46.49	42.30	49.71	46.22	25.95
	TS2^c	13.97	6.68	7.42	10.22	6.41	15.23

---

As can be seen in Fig. 2, the first step of the reaction involves the exothermic formation of an initial complex U–NH₃. This complex has a ⁵A′ ground state and lies at −7.79 kcal mol⁻¹ below the ground state reactants. Then, the U atom inserts into one of the N–H bond to form the intermediate (II) HUNH₂. In this process, the system needs to overcome the first transition state TS1 (⁵A) with an activation barrier of 4.10 kcal mol⁻¹ at the B3LYP/SDD level (10.51 kcal mol⁻¹ at PW91/SDD method). The imaginary frequency of −1290.55 cm⁻¹ at B3LYP/SDD level (−1208.52 cm⁻¹ at PW91/SDD level) correspond primarily to the migration of an H atom from the nitrogen atom to the uranium. This intermediate HUNH₂ is the global minimum of the potential energy profile and lies at −50.65 kcal mol⁻¹ below the asymptote of the reactants.

From the HUNH₂ intermediate, the system can occur along two reaction pathways. The first one is the isomerization. This reaction proceeds to form the uranimine HNUH₂ after the system surpasses the second transition state, TS2a, which is labeled by an imaginary frequency of −1276.30 cm⁻¹ at B3LYP/SDD level (−1197.24 cm⁻¹ at PW91/SDD level), corresponding largely to the sequential movement of H atom to the U atom. Our calculations on HNUH₂ predicted that the molecule should be non-planar and have a triplet (³A) ground state. It is worth noting that there is an intersystem crossing between the quintet and triplet spin surfaces. After this crossing point, the rest of the isomerization reaction evolves along the triplet spin surface, and the quintet PES marches at higher energy with respect to the triplet PES. This feature indicates that this isomerization is a two-state reaction (TSR). Further, we performed ab initio molecular dynamics (MD) simulation to confirm the thermal stability of this product H₂UNH. Simulations were carried out for almost 5 ps using the Born–Oppenheimer molecular dynamics (BOMD)⁴³ method in Gaussian03. Our MD runs indicate that the H₂UNH structure maintains its identity at both room temperature (300 K) and high temperature (500 K). The details of this analysis are included in the ESI.†

Another reaction channel is toward the formation of dehydrogenation products. After the formation of the intermediate II, a second H atom migrating from nitrogen to uranium yields the molecular hydrogen complex, (H₂)UNH (5A′′, complex IV in Fig. 1), overcoming a four-center transition state TS2b. From (H₂)UNH, the system can proceed without any barrier toward UNH + H₂. This transition state has an imaginary frequency of −1256.99 cm⁻¹ at B3LYP/SDD level, corresponding to the expected motion that puts the two hydrogen atoms together. The TS2b lies at −20.59 kcal mol⁻¹ below the reactant asymptote. In the H₂–UNH structure, the H–H distance is 0.750 Å and 0.786 Å at B3LYP/SDD and PW91/SDD methods, respectively. In this case, the hydrogen molecule has already formed and electrostatically interacted with the UNH. The overall dehydrogenation of U + NH₃ was calculated to be exothermic at −38.06 kcal mol⁻¹ at B3LYP level of theory. Our results indicated that the entire lowest energy reaction pathway on this PES is lower in energy than the reactant asymptote.

The possible formation of UNH₂ (quartet state) was also considered. This process may occur from the H–UNH₂ intermediate by the cleavage of the H–U bond. However, we could not locate the transition state in this case, and the scan calculations performed by varying the H–U bond length indicated that this hydrogen–uranium bond breaking process is a barrierless intrinsic transition.

3.1.2 Uranium ions reaction. The potential energy profiles of the reaction of U⁺ and U²⁺ with NH₃ are shown in Fig. 3 and 5, respectively. Fig. 4 and 6 illustrate the geometrical parameters for all the lowest spin-state species.

Fig. 3 Potential energy profile for the reaction of U⁺ + NH₃ computed at the B3LYP/SDD and PW91/SDD (in parentheses) levels of theory.

Fig. 4 Structures and selected geometric parameters of stationary points on the U⁺–NH₃ potential energy surface optimized at the B3LYP/SDD, PW91/SDD, TPSS/SDD, PBE0/SDD, PBE/SDD and BMK/SDD levels of theory (from top to bottom rows, respectively). Bond distances are in Å and angles are in degrees.

Fig. 5 Potential energy profile for the reaction of U²⁺ + NH₃ computed at the B3LYP/SDD and PW91/SDD (in parentheses) levels of theory.

Fig. 6 Structures and selected geometric parameters of stationary points on the U²⁺–NH₃ potential energy surface optimized at the B3LYP/SDD, PW91/SDD, TPSS/SDD, PBE0/SDD, PBE/SDD and BMK/SDD levels of theory (from top to bottom rows, respectively). Bond distances are in Å and angles are in degrees.

Our results indicated that the lowest-energy reactions of U⁺ and U²⁺ with NH₃ evolve along a similar dehydrogenation channel as the U atom. Moreover, it is worth pointing out that these dehydrogenation channels occur under a mechanism similar to the insertion mechanism for the reaction of the bare actinide cations with H₂O and CH₄.^1–4 The formation of U²⁺–NH₃ is exothermic by as much as −72.65 kcal mol⁻¹, which is larger than that computed for the U⁺–NH₃ complex (−42.12 kcal mol⁻¹) and significantly larger than that for the U–NH₃ complex (−7.97 kcal mol⁻¹). This may attributed to the extra charge on U²⁺, which provides significant electrostatic interaction for the U²⁺–NH₃ complex. In addition, the possible formation of UNH₂⁺ and UNH₂²⁺ (quintet and quartet state, respectively) from the H–UNH₂ intermediate by the cleavage of the H–U bond were also considered. Our scan calculation results indicate that this hydrogen–uranium bond breaking process is a barrierless intrinsic transition. It is similar to the formation process of UNH₂ during the reaction of U + NH₃.

Despite these qualitative similarities in the dehydrogenation channels, some differences exist in the mechanism. First, in the case of U⁺ + NH₃, the quartet state remains as the entire lowest-energy reaction path. After the HUNH₂⁺ intermediate, the system does not occur along the two reaction pathways at lowest-energy spin surface such as those found for the U + NH₃ case. The U⁺ reaction leads to the dehydrogenation on the entire lowest-energy reaction path. The isomerization channel competes only at the higher energy doublet spin surface. Second, our study indicated that U²⁺ reacts with NH₃ to produce only UNH²⁺, implying that dehydrogenation is the only channel in this reaction. Unfortunately, despite careful reviews, all the attempts to localize TS1 along the quintet spin surface were unsuccessful. Although this deficiency prevents us from characterizing the quintet spin surface, the TS1 structure is definitely very high in energy, and it is certain that this reaction requires a spin crossing from the quintet spin surface to the triplet spin surface of the products. The change in spin state probably occurs in the vicinity of the first transition state before the formation of the intermediate complex II. The rest of the dehydrogenation reaction evolves along the triplet spin surface.

A comparison of the activation barrier heights (in Table 3) indicates that the reaction of U²⁺ and NH₃ has to surpass the highest first activation barrier (around 57.16 and 33.06 kcal mol⁻¹ higher than U atom and U⁺ at B3LYP/SDD method, respectively). On the other hand, the exothermic energy of the dehydrogenation reaction decreases on going from U to U⁺ and finally U²⁺. The reactivity of uranium species toward ammonia N–H bond activation should be related to the electrons in 5f orbitals because evidence is emerging that the 5f electrons increase the size of activation barrier in the early actinide species and the direct participation of 5f orbitals makes the insertion process unfavorable.⁵

To conclude this section, our results reveal that the uranimine HN=UH₂ is the lowest energy final product in the reaction of the U atom and NH₃, whereas the dehydrogenation products have the lowest energies in the case of U⁺ and U²⁺. The reaction mechanism of the formation of primary reaction products are given below along with reaction energies computed by the relativistic B3LYP/SDD methods.

U + NH₃ → U–NH₃ → HUNH₂ → HNUH₂, ΔE = −47.72 kcal mol⁻¹ (2)

U + NH₃ → U–NH₃ → HUNH₂ → UNH + H₂, ΔE = −38.06 kcal mol⁻¹ (3)

U⁺ + NH₃ → U⁺–NH₃ → HUNH₂⁺ → HNUH₂⁺, ΔE = −6.07 kcal mol⁻¹ (4)

U⁺ + NH₃ → U⁺–NH₃ → HUNH₂⁺ → UNH⁺ + H₂, ΔE = −33.17 kcal mol⁻¹ (5)

U²⁺ + NH₃ → U²⁺–NH₃ → HUNH₂²⁺ → UNH²⁺ + H₂, ΔE = −25.13 kcal mol⁻¹ (6)

3.2 Bonding analysis

To evaluate the bonding evolution of each species involved in the studied reaction pathways, different topological methods such as ELF, AIM and NBO analysis were performed. The ELF calculation did not provide a clear core–valence separation on the U center when the effective core potential (ECP) was used. However, evidence is emerging that electron density obtained from a small-core ECP calculation is sufficient to reproduce not only the correct topology of the uranium–ligand bonds but also the values of the local properties^3,21,45 because the deeper uranium core has little impact on the topological properties of the valence. Therefore, the small-core RECP used in this work is probably responsible for these observations. ELF and AIM analysis were performed using wave functions computed at the PW91/SDD method along with the 60 electron SDD small-core RECP.

The ELF localization domains (η = 0.70) of the lowest-energy minima and transition states, corresponding to the U + NH₃ reaction pathway are collated in Fig. 7, whereas those corresponding to U⁺⁽²⁺⁾ + NH₃ are reported in the Supporting Information (Fig. S1 and S2†). We have also listed the AIM parameters of these species in Tables 4–6. It is worth mentioning that these AIM parameters of bond critical points (bcp) are powerful tools in investigating the characterization and evolution of chemical bonds.^39,46–48

Fig. 7 ELF localization domains (η = 0.70) of the lowest energy minima and transition states corresponding to the U + NH₃ reaction pathway.

Table 4 Topological properties of the charge density calculated at the (3, −1) bond critical points for all the species involved in the U + NH₃ reaction pathway^a

Bond	I (^5A′)		TS1 (^5A)		II (^5A)		TS2a (^3A)		III (^3A)		TS2b (^5A)		IV (^5A′′)
	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)
a ρ(bcp) and ▽^2ρ(bcp) in au. For H2 ρ(bcp) = 0.259 au; ▽^2ρ(bcp) = −1.021 au. For the N–H bond in NH3 ρ(bcp) = 0.329 au; ▽^2ρ(bcp) = −1.398 au.
U–N	0.049	+0.141	0.089	+0.234	0.119	+0.264	0.190	+0.341	0.209	+0.426	0.165	+0.363	0.189	+0.423
N–H1	0.318	−1.382	0.147	−0.136
N–H2	0.318	−1.382	0.319	−1.377	0.316	−1.290	0.128	−0.050			0.132	−0.070
N–H3	0.318	−1.382	0.320	−1.391	0.319	−1.313	0.322	−1.423	0.317	−1.365	0.320	−1.365	0.317	−1.331
U–H1					0.084	+0.039	0.089	+0.016	0.093	+0.005	0.060	+0.136	0.033	+0.101
U–H2									0.093	+0.005
H1–H2											0.132	−0.244	0.237	−0.868

---

Table 5 Topological properties of the charge density calculated at the (3, −1) bond critical points for all the species involved in the U⁺ + NH₃ reaction pathway^a

Bond	I (^4A′)		TS1 (^4A)		II (^4A)		TS2a (^2A)		III (^2A)		TS2b (^4A′)		IV (^4A′)
	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)
a ρ(bcp) and ▽^2ρ(bcp) in au. For H2 ρ(bcp) = 0.259 au; ▽^2ρ(bcp) = −1.021 au. For the N–H bond in NH3 ρ(bcp) = 0.329 au; ▽^2ρ(bcp) = −1.398 au.
U–N	0.063	+0.148	0.119	+0.226	0.140	+0.267	0.222	+0.386	0.229	+0.446	0.184	+0.383	0.213	+0.431
N–H1	0.318	−1.438	0.129	−0.080
N–H2	0.318	−1.437	0.317	−1.434	0.315	−1.365	0.102	+0.046			0.135	−0.079
N–H3	0.318	−1.437	0.317	−1.434	0.316	−1.365	0.317	−1.512	0.308	−1.444	0.318	−1.457	0.316	−1.417
U–H1					0.097	+0.012	0.114	−0.053	0.111	−0.066	0.072	+0.119
U–H2							0.089	+0.069	0.111	−0.066			0.039	+0.138
H1–H2											0.127	−0.219	0.234	−0.853

---

Table 6 Topological properties of the charge density calculated at the (3, −1) bond critical points for all the species involved in the U²⁺ + NH₃ reaction pathway^a

Bond	I (^5A′)		TS1 (^3A)		II (^3A)		TS2 (^3A′′)		IV (^3A′′)		UNH (^3A′′)
	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)	ρ(bcp)	▽^2ρ(bcp)
a ρ(bcp) and ▽^2ρ(bcp) in au. For H2 ρ(bcp) = 0.259 au; ▽^2ρ(bcp) = −1.021 au. For the N–H bond in NH3 ρ(bcp) = 0.329 au; ▽^2ρ(bcp) = −1.398 au.
U–N	0.061	+0.134	0.130	+0.235	0.173	+0.241	0.203	+0.377	0.246	+0.463	0.252	+0.467
N–H1	0.317	−1.462	0.104	−0.001
N–H2	0.317	−1.462	0.309	−1.457	0.299	−1.372	0.180	−0.293
N–H3	0.317	−1.462	0.310	−1.468	0.304	−1.448	0.304	−1.556	0.302	−1.509	0.301	−1.528
U–H1					0.125	−0.129	0.095	+0.012
U–H2									0.038	+0.098
H1–H2							0.097	−0.081	0.237	−0.897

---

3.2.1 U–NH₃ and U⁺⁽²⁺⁾–NH₃. The ELF analysis of these complexes indicates the absence of a disynaptic valence basin between the U and ammonia nitrogen atom, proving that there is no covalent bond formation between the fragments. The interaction between these fragments can be considered to be physically corresponding to an electrostatic interaction. This fact is also supported by the AIM analysis, which shows a very low corresponding charge density (ρ(bcp) = 0.049 au in the case of U–NH₃, 0.063 au in the case of U⁺–NH₃ and 0.061 au in the case of U²⁺–NH₃). The ▽²ρ(bcp) values are small and positive (0.141 au in the case of U–NH₃, 0.148 au in the case of U⁺–NH₃ and 0.134 au in the case of U²⁺–NH₃).

3.2.2 TS1 species. The ELF analysis shows that the first N–H bond breaking occurs at this stage in all the three cases, as evidenced by the lack of the V(N, H1) basin, which is replaced by a trisynaptic V(U, H1, N) basin. According to the AIM analysis, the ρ(bcp) of the N–H1 bond is strongly reduced, although the ▽²ρ(bcp) is still negative but relatively smaller than that of the rest of N–H bonds. These analyses also indicate that the U–N bond is not formed yet but have an important impact on the increase in the electronic density of the U–N bcp. We found that the ρ(bcp) of U–N continues to grow along the reaction paths.

3.2.3 H–U–NH₂ and H–U–NH₂⁺⁽²⁺⁾. The formation of this complex II involves a further isomerized structure and remarkable topological changes. Our ELF analysis reveals that the trisynaptic V(U, H1, N) basin has been replaced by a disynaptic V(H1, U) basin, as indicated by the formation of the U–H1 bond. AIM analysis also illustrates this fact, which shows the presence of a (3, −1) bcp between the An and H1 atoms, even with a low ρ(bcp) and a positive ▽²ρ(bcp) in the cases of U and U⁺. In the case of H–U²⁺–NH₂ complex, the electronic density on this bcp is slightly higher (0.125 au) than that of U (0.084 au) and U⁺ (0.097 au), and the ▽²ρ(bcp) is negative.

3.2.4 TS2a species. The ELF topological analysis shows that the second N–H bond breaks during this process. This fact is proven by the lack of the V(N, H2) basin, which gives its place to another trisynaptic V(U, H2, N) basin. The AIM analysis shows that N–H2 (ρ(bcp) = 0.128 au (U) and 0.102 (U⁺)) is significantly lower than the N–H3 basin, which indicates the breaking of the N–H2 bond. It is also noticed that the ρ(bcp) and ▽²ρ(bcp) of U–N bond increase in all three cases with the ▽²ρ(bcp) values becoming increasingly positive. The higher positive values of ▽²ρ(bcp) suggest a more ionic than covalent nature of these bonds.

3.2.5 HNUH₂ and HNUH₂⁺. These species are the final products of the isomerization channels of the reaction between U and U⁺ and NH₃. The ELF analysis indicates that the N–H2 bond is completely broken during the formation step of HNUH₂, as evidenced by the disappearance of the V(U, H2, N), leading to the formation of a disynaptic V(U, H2) basin. According to AIM, the ▽²ρ(bcp) of U–H1 is equal to the value of U–H2 in both the U and U⁺ cases, which is in line with the axisymmetric nature of the two hydrogen atoms about the U–N bond axis. The ρ(bcp) of U–H increases up to 0.093 au (U) and 0.111 au (U⁺). The difference is that ▽²ρ(bcp) of U–H in HNUH₂⁺ are negative, indicating a covalent nature of this bond.

3.2.6 TS2b species. The noteworthy topological characteristic of these species is the presence of a trisynaptic basin, V(H1, H2, N). In this stage, two hydrogen atoms gradually approach each other. Accordingly, there is an increase in the ρ(bcp) at the H1–H2 in all the three cases (see Tables 4–6). The AIM analysis shows that the N–H2 bond is not completely broken but has an important impact on the decrease in the electronic density of the corresponding bcp.

3.2.7 H₂UNH and H₂UNH⁺⁽²⁺⁾. The U–H1 and N–H2 covalent bonds are completely broken at this step and replaced by the H1–H2 bond. This is evidenced by the lack of two disynaptic valence basins V(U, H1) and V(N, H2), and the formation of a disynaptic V(H1, H2) basin, as shown in the ELF analysis. We also found that the type of basins in these species is similar to those corresponding to the UNH and U⁺⁽²⁺⁾NH combined with a disynaptic V(H1, H2) basin. These facts are supported by the AIM analysis, which indicates that the ρ(bcp) and ▽²ρ(bcp) values for H1–H2 bcp are comparable to that of the free H₂ molecule. In addition, our computation indicates that there is an increase in the ρ(bcp) at the U(or U⁺ and U²⁺)–N bcp, and these values increase steadily from the initial complex I to the imine molecules and dehydrogenation products. The ▽²ρ(bcp) for U–N bcps are still positive and around 0.4 au in all the three cases, which shows that those U–N bond have a dominant ionic character.

3.2.8 NBO analysis. Fig. 8 collates the intrinsic reaction coordinate energy and Wiberg bond order⁴⁹ along with the intrinsic reaction coordinate s of U + NH₃ reaction, which were calculated at the B3LYP/SDD and PW91/SDD levels, while those corresponding to U⁺ and U²⁺ + NH₃ are reported in Fig. S2 and S3 in the ESI.† There has been evidence that for the actinide molecules the hybrid DFT predicts their structures to be more “ionic” and less “covalent” than does a pure GGA.¹⁹ Taking into account that the B3LYP and PW91 methods performed well in the geometry optimization and frequency calculations, as discussed above, we chose these two methods to complement each other to provide a more credible bond analysis.

Fig. 8 IRC energy and Wiberg bond order along the reaction coordinate s calculated at the B3LYP/SDD and PW91/SDD levels corresponding to the U + NH₃. Black solid curves are the energy and colored curves are the Wiberg bond order.

Fig. 8 clearly shows that as the reaction proceeds, the bond order of breaking bonds gradually decreases to zero, while that of the new bonds gradually increases. Furthermore, the position of the transition state is close to the crossing region of these two sets of bond order (the forming bond and the breaking bond). As can be seen in the TS2b process of Fig. 8, the U–N bond order curve has a small valley in the vicinity of the transition state. This may due to the fact that the trisynaptic basin V(H1, H2, N) gathers more electrons of system compared with other basins. Wiberg bond order analysis suggests that the “U–N” bond in those products are in fact showing multiple bond character, with a high bond order (around 3.0) in HNMH₂ (M = U and U⁺) and H₂AnNH (An = U, U⁺ and U²⁺).

To analyze the relative contribution of 5f electrons/orbitals in the title reaction, we performed the NBO analysis for final complexes (III and IV), and the results are summarized in Table 7. As seen in this table, a direct participation of 5f orbitals and a df hybridization with a large contribution of 6d orbitals is observed in these final complexes, as observed during the reaction of early actinide ions with CH₄.⁵ Significant contribution of 5f orbitals is found in the chemical bond of the U²⁺ complex. Therefore, these features suggest that the 5f-orbital participation is important in the reactivity of uranium species.

Table 7 Natural bond orbital analysis for the final complex (III and IV)

Reaction	Complex	Atom	Charge	Electron configuration	Bond character
Isomerization	HNUH2 (III)	U	−0.359	5f^2.31 6d^0.75	f[22.83%] d[62.70%]
				7s^0.33 7p^0.03	s[14.24%] p[0.22%]
		N	−0.487	2s^0.79 2p^2.20	s[56.40%] p[ 43.59%]
	HNUH2^+ (III)	U	1.938	5f^2.67 6d^1.18	f[60.68%] d[38.72%]
				7s^0.27 7p^0.07	s[0.29%] p[0.30%]
		N	−0.900	2s^0.79 2p^2.10	s[2.74%] p[97.25%]
Dehydrogenation	H2–UNH (IV)	U	0.928	5f^3.12 6d^1.10	f[29.66%] d[61.89%]
				7s^0.90 7p^0.03	s[8.18%] p[0.27%]
		N	−1.283	2s^1.62 2p^4.65	s[53.20%] p[46.80%]
	H2–UNH^+ (IV)	U	1.784	5f^3.27 6d^0.91	f[42.37%] d[57.54%]
				7s^0.11 7p^0.02	s[0.00%] p[0.09%]
		N	−1.213	2s^1.63 2p^4.57	s[0.00%] p[99.99%]
	H2–UNH^2+ (IV)	U	2.395	5f^2.90 6d^0.74	f[68.10%] d[31.83%]
				7s^0.05 7p^0.01	s[0.00%] p[0.07%]
		N	−0.945	2s^1.63 2p^4.30	s[0.00%] p[99.98%]

---

3.3 A comparison between the reaction of uranium species with H₂O, CH₄ and NH₃

Water, methane and ammonia are isoelectronic molecules.^14,50 The mechanisms of the gas phase reaction between the U atom with H₂O has been studied by Andrews and co workers¹² by reacting laser-ablated U atoms with H₂O during condensation under excess argon. The isomerization channel (H₂UO was the product) and the H₂ elimination channel (decomposed to UO + H₂) were observed. These two similar pathways have also been addressed in our present study in the U + NH₃ reaction. To the best of our knowledge, the gas-phase study of U atom + CH₄ is scarce, and little is known about their intermediates and mechanisms.

The reaction of U cations (U⁺ and U²⁺) with H₂O or CH₄ have been well-studied by several earlier investigations.^1,5,21,24,51 Comparisons of the obtained results show that the reaction mechanisms are similar but the energetics are significantly different.

As seen in Fig. 3 and 5 in the present work and the PES of U⁺⁽²⁺⁾ reacting with H₂O and CH₄ in the literature mentioned above, the behavior of the mechanisms is similar in all the systems considered, especially for the dehydrogenation channel. It can be summarized that in the first stage the reactants form a stable ion–dipole complex (I), U–MH_x⁺⁽²⁺⁾ (M = O, N and C; x = 2, 3 or 4, respectively). Then, the U⁺⁽²⁺⁾ inserts into the M − H bond, followed by the first M − H bond breaking, which is realized by a first transition state TS1. After that, the system forms the inserted intermediate (II) HUMH_y⁺⁽²⁺⁾ (y = 1, 2 or 3), which lies lower in energy than the reactants. After passing through a four-center transition state TS2b, the reaction proceeds to yield a molecular hydrogen complex, (H₂)UMH_z⁺⁽²⁺⁾ (z = 1 or 2). Finally, UMH_z⁺⁽²⁺⁾ + H₂ products are formed directly from this hydrogen complex.

Another similarity is that the entire dehydrogenation path occurs along the ground spin surface, and the intermediate II is the global minimum in U⁺ + H₂O, CH₄ and NH₃. Whereas in the case of U²⁺, the dehydrogenation channel evolves an intersystem, crossing between the quintet and triplet state, and the global minimum is the complex I.

In addition to the similarities between mechanisms, there are differences in the energetics. The relative energy of complex I and dehydrogenation products with respect to the ground state reactants as well as the activation barriers of the transition states for the reaction of U⁺⁽²⁺⁾ + H₂O, CH₄ and NH₃ are listed in Table 8. To guarantee the reliability, we selected the reported values that were computed by the same method B3LYP/SDD for the comparison. It can be seen that the relative energy of complex I decreases on going from NH₃ to H₂O and finally CH₄. The reaction with NH₃ has the highest activation barrier in TS1, whereas the highest activation barrier of TS2b is found in the reaction with H₂O. Both NH₃ and H₂O cases are exothermic, and the most exothermic is the reaction with H₂O. On the contrary, U⁺⁽²⁺⁾ + CH₄ were computed to be endothermic, particularly the U²⁺ should not activate CH₄, which is in accordance with the experimental observations.²⁴

Table 8 The relative energy and activation barriers (in kcal mol⁻¹) in the dehydrogenation of U⁺⁽²⁺⁾ + H₂O, CH₄ and NH₃. Reaction computed at the B3LYP/SDD level of theory

Reaction	ΔEI^a	TS1^b	TS2b^c	ΔEp^d
a Energy difference between the complex I and the ground state reactants.b Energy difference between the transition state TS1 and the second complex I.c Energy difference between the transition state TS2b and the third complex II.d Energy difference between the dehydrogenation products and the ground state reactants.e in the present work.f Ref. 1.g Ref. 21.h Ref. 24.
U^+ + NH3^e	−42.12	28.20	24.69	−33.17
U^2+ + NH3^e	−72.65	61.26	13.96	−25.13
U^+ + H2O^f	−30.91	14.52	31.56	−49.65
U^2+ + H2O^f	−52.26	34.27	29.79	−30.28
U^+ + CH4^g	−21.6	18.88	27.48	16.94
U^2+ + CH4^h	−24.61	39.67	9.08	24.13

---

4. Conclusions

The theoretical calculations have been performed to investigate the mechanism of the gas phase reaction of U and U⁺⁽²⁺⁾ with NH₃. The present results contribute to a better understanding of the two types of reaction channels, isomerization and dehydrogenation, for the title reactions. The bonding nature was analyzed on the basis of ELF, AIM and NBO calculations. The main conclusions can be summarized as follows.

(i) For the U + NH₃ reaction, there is a spin surface crossing in the isomerization pathway. It is different from that of the U⁺ + NH₃ and U²⁺ + NH₃, in which no spin crossing was found on the lowest energy surface. Our calculations indicate that the lowest energy pathway of both the channels is exothermic, which does not involve reaction barriers over the reactant asymptotes. The formation of HNUH₂ is the reaction with the higher exothermicity, calculated at −47.72 kcal mol⁻¹ at B3LYP/SDD level of theory.

(ii) In the case of U⁺ + NH₃, the minimum energy pathway was computed as the dehydrogenation. The dehydrogenation pathway evolves along the quartet state from the initial complex to the final products. The isomerization channel competes only at the higher energy doublet spin surface.

(iii) For the U²⁺ + NH₃ reaction, UNH²⁺ was found to be the only product generated from the dehydrogenation. This reaction has to surpass the highest barrier for the first transition state compared to that of the U and U⁺ cases. On the contrary, the second barrier height of U²⁺ and NH₃ is the lowest.

(iv) ELF analysis indicated that the U–N bond in all species is characterized by high ionic character. This fact is supported by the AIM results, which provide important insights into the chemical bond evolution of the reactions. The results indicate that the ▽²ρ(bcp) for U–N bcps becomes more positive, indicating that the ionic character of the U–N chemical bonds becomes stronger.

Acknowledgements

We are very grateful to Dr Sobereva for many helpful discussions and who provided us the Multiwfn package. Computer time made available by the Center of High Performance Computing at Physics discipline of Sichuan University is gratefully acknowledged.

Notes and references

1. M. D. C. Michelini, N. Russo and E. Sicilia, Angew. Chem., Int. Ed., 2006, 45, 1095–1099.
2. M. D. C. Michelini, N. Russo and E. Sicilia, J. Am. Chem. Soc., 2007, 129, 4229–4239.
3. M. E. Alikhani, M. D. C. Michelini, N. Russo and E. Sicilia, J. Phys. Chem. A, 2008, 112, 12966–12974.
4. G. Mazzone, M. D. C. Michelini, N. Russo and E. Sicilia, Inorg. Chem., 2008, 47, 2083–2088.
5. K. J. de Almeida and H. A. Duarte, Organometallics, 2009, 28, 3203–3211.
6. K. J. de Almeida and H. A. Duarte, Organometallics, 2010, 29, 3735–3745.
7. J. Zhou and H. B. Schlegel, J. Phys. Chem. A, 2010, 114, 8613–8617.
8. M. Santos, J. Marcüalo, A. P. de Matos, J. K. Gibson and R. G. Haire, J. Phys. Chem. A, 2002, 106, 7190–7194.
9. G. P. Jackson, F. L. King, D. E. Goeringer and D. C. Duckworth, J. Phys. Chem. A, 2002, 106, 7788–7794.
10. B. Y. Liang, L. Andrews, J. Li and B. E. Bursten, J. Am. Chem. Soc., 2002, 124, 6723–6733.
11. J. K. Gibson, R. G. Haire, M. Santos, J. Marcüalo and A. P. de Matos, J. Phys. Chem. A, 2005, 109, 2768–2781.
12. B. Y. Liang, R. D. Hunt, G. P. Kushto, L. Andrews, J. Li and B. E. Bursten, Inorg. Chem., 2005, 44, 2159–2168.
13. M. H. Chen, H. Lu, J. Dong, L. Miao and M. F. Zhou, J. Phys. Chem. A, 2002, 106, 11456–11464.
14. N. Russo and E. Sicilia, J. Am. Chem. Soc., 2001, 122, 2588–2596.
15. M. D. C. Michelini, N. Russo and E. Sicilia, Inorg. Chem., 2004, 43, 4944–4952.
16. M. F. Zhou, M. H. Chen, L. N. Zhang and H. Lu, J. Phys. Chem. A, 2002, 106, 9017–9023.
17. X. F. Wang, L. Andrews and C. J. Marsden, Chem.–Eur. J., 2008, 14, 9192–9201.
18. X. F. Wang, L. Andrews and C. J. Marsden, Chem.–Eur. J., 2007, 13, 5601–5606.
19. G. A. Shamov, G. Schreckenbach and T. N. Vo, Chem.–Eur. J., 2007, 13, 4932–4947.
20. W. Koch and M. C. Holthausen, A Chemist's Guide to Density Functional Theory, Wiley-VCH, New York, 2000.
21. E. Di Santo, M. D. C. Michelini and N. Russo, Organometallics, 2009, 28, 3716–3726.
22. E. Di Santo, M. D. C. Michelini and N. Russo, J. Phys. Chem. A, 2009, 113, 14699–14705.
23. P. Li, W. X. Niu, X. F. Tian, T. Gao and H. Y. Wang, Int. J. Quantum Chem., 2004, 114, 760–768.
24. E. Di Santo, M. Santos, M. D. C. Michelini, J. Marçalo, N. Russo and J. K. Gibson, J. Am. Chem. Soc., 2011, 133, 1955–1970.
25. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377.
26. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
27. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
28. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671–6682.
29. J. P. Perdew, K. Burke and Y. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 16533–16539.
30. J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
31. C. Adamo and V. Barone, J. Chem. Phys., 1999, 110, 6158–6169.
32. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
33. A. D. Boese and M. L. M. Jan, J. Chem. Phys., 2004, 121, 3405–3416.
34. W. Kuchle, M. Dolg, H. Stoll and H. Preuss, J. Chem. Phys., 1994, 100, 7535–7542.
35. (a) R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 1980, 72, 650–654; (b) J. P. Blaudeau, M. P. McGrath, L. A. Curtiss and L. Radom, J. Chem. Phys., 1997, 107, 5016–5021.
36. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, J. A. Montgomery, Jr, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzalez and J. A. Pople, Gaussian 03, Revision E.01, Gaussian Inc., Wallingford CT, 2004.
37. A. D. Becke and K. E. Edgecombe, J. Chem. Phys., 1990, 92, 5397–5403.
38. A. Savin, R. Nesper, S. Wengert and T. R. Fassler, Angew. Chem., Int. Ed. Engl., 1997, 36, 1808–1832.
39. T. Lu and F. Chen, J. Comput. Chem., 2012, 33, 580–592.
40. R. F. W. Bader, Atoms in Molecules. A quantum theory, Clarendon, Oxford, 1990.
41. A. E. Reed and F. Weinhold, J. Chem. Phys., 1985, 83, 1736–1740.
42. A. E. Reed, L. A. Curtiss and F. Weinhold, Chem. Rev., 1988, 88, 899–926.
43. V. Bakken, J. M. Millam and H. B. Schlegel, J. Chem. Phys., 1999, 111, 8773–8777.
44. J. Blaise and J.-F. Wyart, International Tables of Selected Constants, Energy Levels and Atomic Spectra of Actinides, Vol. 20, Tables of Constants and Numerical Data, Paris, 1992.
45. S. F. Vyboishchikov, A. Sierraalta and G. Frenking, J. Comput. Chem., 1996, 18, 416–429.
46. C. H. Suresh and N. Koga, Inorg. Chem., 2002, 41, 1573–1578.
47. Y. Zhao and D. G. Truhlar, J. Phys. Chem. A, 2004, 108, 6908–6918.
48. C. H. Suresh and P. K. Sajith, Inorg. Chem., 2011, 50, 8085–8093.
49. K. B. Wiberg, Tetrahedron, 1968, 24, 1083–1096.
50. H. A. Bent, J. Chem. Educ., 1996, 43, 170–186.
51. P. Li, W. X. Niu, X. F. Tian, T. Gao and H. Y. Wang, J. Phys. Chem. A, 2013, 117, 3761–3770.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4ra03525a

---