Modeling the absorption spectrum of the permanganate ion in vacuum and in aqueous solution

TL;DR

This paper demonstrates that combining polarizable embedding with CAS-srDFT and molecular dynamics effectively models the solvent and temperature effects on the MnO₄⁻ absorption spectrum, surpassing traditional implicit solvent approaches.

Contribution

It introduces a novel computational approach integrating PE, CAS-srDFT, and MD simulations to accurately predict solvent and temperature effects on the MnO₄⁻ spectrum.

Findings

Implicit solvent models are insufficient for accurate solvent shift prediction.

The combined PE and CAS-srDFT approach successfully reproduces experimental solvent shifts.

Finite-temperature averaging clarifies the spectral band assignments.

Abstract

The absorption spectrum of the MnO$_{4}$$^{-}$ ion has been a test-bed for quantum-chemical methods over the last decades. Its correct description requires highly-correlated multiconfigurational methods, which are incompatible with the inclusion of finite-temperature and solvent effects due to their high computational demands. Therefore, implicit solvent models are usually employed. Here we show that implicit solvent models are not sufficiently accurate to model the solvent shift of MnO$_{4}$$^{-}$, and we analyze the origins of their failure. We obtain the correct solvent shift for MnO$_{4}$$^{-}$ in aqueous solution by employing the polarizable embedding (PE) model combined with a range-separated complete active space short-range density functional theory method (CAS-srDFT). Finite-temperature effects are taken into account by averaging over structures obtained from ab initio molecular dynamics simulations. The explicit treatment of finite-temperature and solvent effects facilitates the interpretation of the bands in the low-energy region of the MnO$_{4}$$^{-}$ absorption spectrum, whose assignment has been elusive.

Tables

Table 1. Table 1: Vertical excitation energies in eV and oscillator strengths of the lowest intense transitions of the \ce MnO4- ion in vacuum or aqueous solution calculated using geometry-optimized structures

Method	CAS(14,12)-srPBEa				PCMb-CAS(14,12)-srPBEa				RAS(24,17)PT220		experimentc
Geom.	vacuum		PCMb		vacuum		PCMb		vacuum		vacuum21	aq. sol.21	crys.16
State	$E_{\mathrm{vac}}$	$f$	$E_{\mathrm{vac}}$	$f$	$E_{\mathrm{aq}}$	$f$	$E_{\mathrm{aq}}$	$f$	$E_{\mathrm{vac}}$	$f$	$E_{\mathrm{vac}}$	$E_{\mathrm{aq}}$	$E_{\mathrm{xtal}}$
$1^1{T}_2$	2.21	0.095	2.25	0.092	2.16	0.127	2.20	0.123	2.33	0.004	2.09/2.18	2.27/2.35	2.3/2.4
$2^1{T}_2$	3.61	0.010	3.64	0.010	3.61	0.018	3.64	0.018	3.53	0.002	3.55	-	3.6
$3^1{T}_2$	4.01	0.084	4.05	0.082	-	-	-	-	4.20	0.006	3.94	4.00	4.1
$4^1{T}_2$	5.13	0.171	5.19	0.162	-	-	-	-	5.72	0.002	-	-	5.4

Table 2. Table 2: Vertical excitation energies in eV and oscillator strengths of the twelve lowest intense transitions of the \ce MnO4- ion in vacuum or aqueous solution calculated as averages using structures from a molecular dynamics trajectory a

Geom.		QM MDb		QM/MM MDc
Environment		vacuumd		vacuume		PCMf		PEg
Multiplet	State	$E_{\mathrm{vac}}$	$f$	$E_{\mathrm{aq}}$	$f$	$E_{\mathrm{aq}}$	$f$	$E_{\mathrm{aq}}$	$f$
	1	1.93	0.013	1.97	0.008	1.98	0.011	1.91	0.010
11T1	2	2.03	0.005	2.07	0.005	2.07	0.008	2.04	0.003
	3	2.08	0.006	2.13	0.007	2.13	0.008	2.14	0.006
	4	2.19	0.014	2.23	0.017	2.22	0.015	2.30	0.016
11T2	5	2.29	0.019	2.34	0.018	2.32	0.019	2.45	0.014
	6	2.41	0.016	2.46	0.016	2.44	0.019	2.55	0.012
	7	3.36	0.002	-	-	3.42	0.003	3.33	0.003
	8	3.39	0.002	-	-	3.45	0.002	3.40	0.003
21T1/21T2	9	3.42	0.002	-	-	3.48	0.002	3.45	0.003
	10	3.46	0.002	-	-	3.52	0.002	3.50	0.002
	11	3.49	0.002	-	-	3.56	0.002	3.55	0.003
	12	3.54	0.004	-	-	3.61	0.006	3.61	0.004

Table 3. Table 3: Absorption maxima and solvent shifts in eV of the two lowest bands of \ce MnO4- ion obtained from the spectrum in figure  1

Band	$E_{\mathrm{vac}}$	$E_{\mathrm{aq}}^{\mathrm{PCM}}$	$E_{\mathrm{aq}}^{\mathrm{PE}}$	$E_{\mathrm{vac}}^{\exp }$21	$E_{\mathrm{aq}}^{\exp }$21
1st	2.25	2.27	2.40	2.09/2.18	2.27/2.35
2nd	3.44	3.52	3.49	3.55	-
		$\Delta E_{\mathrm{shift}}^{\mathrm{PCM}}$	$\Delta E_{\mathrm{shift}}^{\mathrm{PE}}$		$\Delta E_{\mathrm{shift}}^{\exp }$
1st	-	0.02	0.15	-	0.17/0.18
2nd	-	0.08	0.05	-	-

Full text

University of Southern Denmark] Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark.

Lund University] Department of Chemistry, Lund University, Kemicentrum Sölvegatan 39, Lund, Sweden.

1 Introduction

Inclusion of environment effects in computational chemistry is pivotal to achieve the best possible match between theoretical models and experimental measurements. The presence of an environment influences both energetics and spectroscopic constants. An example of the latter is how the presence of a molecular environment can alter absorption spectra, and many chemical compounds have a characteristic shift of their excitation energies that depends on the specific environment. For solutes, the solvent shift may be defined as $\Delta E_{\mathrm{shift}}=E_{\mathrm{sol}}-E_{\mathrm{vac}}$, where $E_{\mathrm{sol}}$ and $E_{\mathrm{vac}}$ are excitation energies obtained with and without solvent included, respectively. Accurate prediction of excitation energies and solvent shifts poses special demands on the theoretical model. Obviously it is crucial that the theoretical method is able to describe the electronic structure to high accuracy. For a solvated compound, it is also important to accurately model the environment-induced perturbation of the solute’s charge density. Furthermore, it is necessary to include finite-temperature effects which can be done by averaging over several structures from a molecular dynamics (MD) trajectory. Transition metal complexes can generally be expected to display multireference character, and in such cases both dynamical and static correlation must be modeled accurately with multireference methods, such as complete active space second-order perturbation theory (CASPT2).1 A prime example is the \ceMnO4- ion whose absorption spectrum has been a test-bed for quantum-chemical methods2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 for decades. Despite its apparent simplicity, the closed-shell (d0) \ceMnO4- ion has been shown to have sizable non-dynamical correlation.15 Ab initio methods based on a single-reference wave function have resulted in different assignments4, 5, 6, 8, 13 of the experimental spectrum16. The ubiquitous time-dependent density functional theory (TD-DFT) has also been applied, but is noted to significantly overestimate the absolute excitation energies7, 9, 10, although recent developments by Ziegler and coworkers seems to offer some improvement14. An obstacle for the application of multireference methods has been that the required active space has been estimated to be 24 electrons in 17 orbitals, which is rather large.17, 18, 19 A recent study of the electronic absorption spectrum managed to reach the required active space with a restricted active space method (RASPT2).20 Although good agreement with experiment was obtained for the intense transitions, there are still aspects of the \ceMnO4- absorption spectrum that remains puzzling. Transitions in the low-energy region are still not assigned16, 21, and with the measurement of accurate gas-phase data21, experimental estimates of the vacuum to aqueous solvent shift for the first and third electronic bands are available; both are small blue-shifts. Yet previous theoretical studies employing continuum solvation models have predicted either a negligible solvent shift of $0.01-0.02$ eV20, 14 or a somewhat more pronounced red-shift of $-0.13$ eV10, 22, in clear contrast to experiment. It is presently not known whether the discrepancies are caused by inadequate electronic structure methods, solvent models, or both.

In this paper, we address the theoretical description of the \ceMnO4- absorption spectrum in both vacuum and aqueous solvent. To ensure an accurate description it seems that methods that can capture multiconfigurational character are required. Here we use a method that combines Kohn-Sham (KS) density functional theory (DFT) with a multiconfigurational wave function. A number of such schemes are currently in development 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41 and here we employ a range-separation method denoted complete active space short-range DFT (CAS-srDFT).27, 28, 31, 32, 33, 35 The CAS-srDFT method is computationally cheaper than perturbation-based multireference methods and has been shown to be of comparable accuracy35, 36, 34. Further, its simultaneous treatment of dynamical and static correlation enables the use of significantly smaller active spaces without much loss in accuracy. We exploit this here to go beyond what have been done so far, and also consider finite-temperature effects as well as the solvent effects on the electronic absorption spectrum. Thus, we calculate CAS-srDFT excitation energies and oscillator strengths using a series of structures taken from ab initio MD and combined quantum mechanics / molecular mechanics (QM/MM) ab initio MD simulations. From these underlying structures we also show that continuum solvation models describe the electrostatic perturbation by the solvent inadequately. We will instead employ the polarizable embedding (PE) model42, 43 which uses an advanced, classical potential to model the electrostatic solvent effects. With the combined use of QM/MM MD and PE-CAS-srDFT44, we are able to reproduce the experimentally-obtained blue shift of the lowest intense transition. Furthermore, we assign the low-energy parts of the \ceMnO4- absorption spectrum.

2 Results and discussion

First we investigate the performance of the CAS-srDFT method on the \ceMnO4- ion in vacuum and in aqueous solvent described by the polarizable continuum model (PCM). The four lowest intense transitions in Td symmetry (i.e. transitions to the T2 states) are examined. We use symmetric, geometry-optimized structures, i.e., we neglect finite-temperature effects. We performed the calculations based on two geometries: a vacuum geometry (BLYP/6-31G) and a solvent (PCM) geometry (PCM-BLYP/6-31G). The lowest vertical excitation energies, as well as the associated oscillator strengths, were obtained for both structures based on time-dependent CAS(14,12)-srPBE calculations with and without PCM embedding. This allows us to separate the solvent effects on the spectrum into a contribution stemming from the direct electrostatic interactions that polarizes the electron density of the solute, and an indirect contribution which is caused by the solvent-induced change in molecular geometry.

The results for the four lowest intense transitions are shown in table 1. The CAS(14,12)-srPBE calculation based on the vacuum geometry gives the values that can be compared to the accurate results based on RAS(24,17)PT2.20 In terms of excitation energies both methods agree to within 0.1 eV for the first two transitions, and 0.2 eV for the third transition. However, there is a relatively large discrepancy for the fourth transition where they differ by 0.6 eV. We note that Su et al.20 actually report two states with low oscillator strength below $4^{1}\mathrm{T}_{2}$, and the state that we have designated as $4^{1}\mathrm{T}_{2}$ is actually $6^{1}\mathrm{T}_{2}$ in ref. 20. The oscillator strengths are rather different though the same trends are observed except for the fourth state. This is not unprecedented, as we recently showed that CASPT2 and CAS-srDFT can obtain oscillator strengths that differ significantly34, although they usually predict the same trends. Some of the differences have been attributed to the fact that oscillator strengths in CASPT2 (and RASPT2) methods are effectively obtained at the CASSCF level34. For these higher-lying excitations (close to the ionization limit) the inclusion of Rydberg basis sets can be of importance, although it had only little effect in ref. 20. CAS-srDFT has previously been noted to be more sensitive than e.g. CASPT2 to Rydberg-valence mixing35, 36 which could also be a possible explanation for the differences. CAS(14,12)-srPBE compares very well with experimental values for the first three excitations for which there exist experimental values. Although there are no experimental values for the fourth transition in vacuum it has been measured in crystalline phase, and comparing to this value it seems that CAS(14,12)-srPBE slightly underestimates the excitation energy whereas RAS(24,17)PT2 overestimates it by about the same amount.

The most consistent estimates of the solvent shifts are obtained as the difference between PCM-CAS(14,12)-srPBE values using the PCM geometry and the CAS(14,12)-srPBE values using the vacuum geometry. However, previous theoretical studies of \ceMnO4- have occasionally neglected the indirect geometry effect of the solvent.10, 22 The shift of $-0.14$ eV for the first intense transition reported in ref. 10 is based on a vacuum excitation energy that was obtained from a structure optimized in vacuum (BPW91/TZ2P), and the solvated excitation energy was from the work by Seth et al.22 where COSMO-SAOP/TZ2P was used on a structure that was also optimized in vacuum (BP86/TZ2P). Meanwhile, the small shift (0.01 eV) reported by Su et al.20 was obtained with SAOP/TZ2P and COSMO. If both structure and excitation energy are obtained consistently, then our estimate of the solvent shift is only $-0.01$ eV, which is in good agreement with the shift obtained by Su et al.20 but differs somewhat from the red-shift of $-0.14$ eV. The discrepancy can, at least partly, be explained by the neglect of indirect solvent effects since we also obtain a small red-shift of $-0.05$ eV when using structures optimized in vacuum (see table 1). All the same, even with a consistent procedure, implicit solvent models alone do not seem to be able to reproduce the experimental blue-shift which is around 0.1-0.2 eV.21

We now consider \ceMnO4- in an explicit solvent environment. We computed absorption spectra based on the twelve lowest states and in table 2 we report $E_{\mathrm{vac}}$ and $E_{\mathrm{sol}}$ for these twelve states as averages based on structures taken from MD trajectories49. Gaussian convolution of these states leads to the spectrum shown in figure 1. The two intense peaks in this figure correspond to the two lowest intense bands of 1T2 parentage (corresponding to 11T2 and 21T2 for the symmetric structures in table 1). The perturbation induced by the solvent molecules as well as the finite-temperature effects remove the degeneracy of the T1 and T2 states. Thus, the T1 and T2 states are split into six states in total. Further, the transition to the T1 state is dipole-forbidden in Td symmetry, but gains intensity when the symmetry is removed by finite-temperature effects and solvent interactions. The peak-maxima in figure 1 are used to calculate $\Delta E_{\mathrm{shift}}$ which are given in table 3. We start by discussing these, and thereafter we will provide a more detailed analysis of the individual states in table 2.

From the spectrum we obtain an excitation energy of 2.25 eV for the 11T2 state in vacuum. This is 0.04 eV higher than the energy obtained using the geometry-optimized structure (see table 1) which can be attributed to finite-temperature effects in the ground-state. The PCM predicts a very small blue-shift of 0.02 eV (see table 3) and is thus similar to the shift obtained from the geometry-optimized structures. The PE model predicts a more significant blue-shift of 0.15 eV (see table 3) which is in excellent agreement with the experimental shift.

The spectrum in figure 1 can be interpreted on the basis of the individual states given in table 2. We can compare these to the results obtained with a symmetric structure in table 1. For the symmetric structure used in table 1 we obtain an excitation energy of 2.13 eV for the 11T1 multiplet (not shown in table 1), and with the excitation energy of 2.21 eV for the 11T2 state (table 1), the splitting of the 11T1 and 11T2 states in vacuum is only 0.11 eV, similar to results with RASPT220. With such a small energy-gap, the two states are expected to mix when the tetrahedral symmetry is lifted. This is indeed the case, but among the lowest six states, the 11T1 and 11T2 multiplets can still be discerned (see table 2). The three states of 11T1 parentage are still lower in energy than the intense 11T2 band, and can also be discerned in the spectrum in figure 1, as a weak transition at approximately 1.7-1.9 eV (13700-15300 cm*-1*) prior to the large peak from 11T2. Low-intensity peaks have indeed previously been observed in \ceMnO4- spectra around 1.8 eV.16, 21 Analysis have suggested that their origin were either in spin-forbidden 3T2, orbital-forbidden 1T1 states, or both.50, 51, 52. From our calculations here, we cannot exclude the involvement of triplet states, but our results show that the 1T1 states are likely to gain sufficient intensity to be responsible for the observed low-intensity band. This assignment can naturally not be extracted from previous studies22, 10, 20 with symmetric structures, where the intensity of the T1 states are identically zero.

From inspection of table 2 we also observe that the direct electrostatic effect of the solvent is predicted different by the PE model compared to the PCM. In both cases, the indirect solvent effect, coming from the change in geometry of the solvated \ceMnO4- ion, results in a blue-shift of the excitation energies. This can be seen in table 2 by comparing the excitation energies based on QM MD geometries to QM/MM MD geometries both with vacuum environment. For the PCM the direct solvent effects (obtained from column 5 and 7 in table 2) is a (small) red-shift of the excitation energies. Adding this to the indirect solvent shift thus gives the very small blue-shift as discussed above. Meanwhile, the direct solvent effects modeled by the PE model (see table 2) result in a blue-shift of the excitation energy, leading to a more pronounced total blue-shift.

We can also briefly comment on the higher-lying states in table 2 and figure 1. Our CAS-srPBE results from the symmetric structure in table 1 showed that the dipole-allowed 21T2 multiplet is just above the dipole-forbidden 21T1 multiplet: The excitation energy of the former is 3.61 eV (see table 1) whereas the 21T1 multiplet has an excitation energy of 3.48 eV (not shown in table 1). The RAS(24,17)PT2 results from ref. 20 predict a similar situation where the excitation energies are 3.53 eV (21T1) and 3.39 eV (21T1), respectively. These two multiplets are mixed when including finite-temperature and solvent effects and can no longer be discerned (see table 2). From the spectra in figure 1 we obtain an excitation energy of 3.44 eV in vacuum, and blue-shifts of 0.08 and 0.05 for the PCM and the PE model, respectively (see table 3). There is no experimental value for this transition in aqueous solution but a value of 0.09 eV exists for a 0-0 transition in crystalline phase.21, 16 Using this value as reference, both solvent models yield a qualitatively reasonable shift (see table 3), in particular in the light that the shift of excitation energies due to an environment seem to be larger in the crystalline phase than in water.21 Our calculated shifts for water are indeed slightly below the value for the crystalline phase.

It should finally be noted that we do not reproduce the vibrational structure of the spectrum, since our underlying structures only include vibrations from the electronic ground state, while the significant vibrational structure in the \ceMnO4- spectrum is a result of vibrational coupling in the excited states. The correct description of this vibronic coupling was a notable achievement of ref. 10. Although the authors note that the absolute excitation energies are significantly overestimated by the TD-DFT, they also find that the method convincingly reproduces the vibrational structure through an effective Hamiltonian approach. Employing a similar model for CAS-srDFT would be an interesting extension of our current approach but is left for further studies. In any case, the good agreement with experiment obtained here is encouraging for further investigating the PE-CAS-srDFT method for other transition metals, both in solution and in protein systems.

3 Conclusion

In this paper, we present a study of \ceMnO4- in vacuum and in aqueous solution, employing the multiconfigurational CAS-srDFT method combined with implicit and explicit solvation models. In vacuum calculations, we show that CAS-srDFT is at least as accurate as RASPT2, although we can employ significantly smaller active spaces with CAS-srDFT. To address the solvent shift, we employed a combination of CAS-srDFT and the PE model (explicit solvation), as well as the PCM (implicit solvation). Finite-temperature effects are included by using structures obtained from QM and QM/MM MD trajectories. An experimental study has predicted that the lowest intense transition is blue-shifted by about 0.1-0.2 eV compared to vacuum. However, previous theoretical studies that employed continuum solvation models predicted that the solvent shift is essentially zero. We find that the reason lies in an unsatisfactory description of the direct electrostatic interactions between solute and solvent. The solvent-induced perturbation of the solute geometry shifts the excitation energy about 0.05 eV. Adding the direct electrostatic interactions from the solvent via PCM red-shifts the energy by about 0.03 eV, leading to a total blue-shift of only 0.02 eV. Using the PE model instead adds an additional blue-shift of about 0.1 eV, so that the total shift is 0.15 eV. Thus the combined use of the CAS-srDFT and the PE model yield both a solvent shift and absolute excitation energies that are in excellent agreement with experiment.

{acknowledgement}

J. M. H. O. acknowledges financial support from the Danish Council for Independent Research (DFF) through the Sapere Aude research career program (grant id. DFF – 1325-00091 and DFF – 1323-00744) and the Carlsberg Foundation (grant id. CF15-0823). E. D. H. thanks the Carlsberg Foundation (grant id. CF16-0482) and the Villum Foundation for funding.