Accurate Absorption Energy Calculations in Solution Using the Reference Interaction Site Model Self-Consistent Field Including the Constrained Spatial Electron Density Distribution

The solvation effect is an important factor determining the properties of molecules in solution. The reference interaction site model (RISM) is a powerful method to treat the solvation effect with pair-correlation functions, such as a radial distribution function. This study developed a hybrid method between quantum mechanics and RISM using the spatial electron density distributions on each atomic site (RISM-SCF-cSED). Sophisticated quantum mechanical approaches can be used to consider the solvation effect because the computational cost of RISM-SCF-cSED is reasonable. In this study, the absorption energies of 5-(dimethylamino)-2,4-pentadienal in various solutions were calculated using RISM-SCF-cSED. The experimental data were well reproduced with an average errors of ∼0.06 eV, using multi-reference perturbation theory.


Introduction
The solvation effect is an important factor to determining the properties of molecules in solution, such as the effect of the solvent polarity on the reactivity of chemical reactions in solution [1][2][3], the absorption and emission energies of the imaging molecules [4,5], and the population of isomers [6]. The electronic structure of the solute molecule of interest and the solvation effect needs to be understood to elucidate the reaction mechanism and the energy decay process of these systems.
The statistical approach has great advantages when evaluating the solvation effect [7]. The reference interaction site model (RISM) is one of the methods developed by many scientists [8][9][10]. In the RISM, pair-correlation functions, such as a radial distribution function (RDF), are calculated by solving integral equations. Although these correlation functions can be calculated using molecular dynamics calculations, the trajectories with a long-time simulation are required. However, the trajectories are not required in RISM because it calculates the correlation functions directly by solving integral equations.
The computational efficiency in the RISM became important when it was coupled with quantum mechanical (QM) approaches (RISM-SCF) [11,12]. RISM-SCF has excellent potential in solution chemistry. The electronic structures of a solution in the ground and excited states can be discussed using sophisticated quantum chemical methods. Moreover, the solvation effect can be discussed at the atomic level using RDFs. However, the solvation effect in polar solvent is greatly overestimated in the original RISM-SCF [13,14]; thus, obtaining the converged electronic structure in polar solvents was challenging.
A new method was developed to overcome the challenge in the original RISM-SCF by including the spatial electron density distribution (SED) [13][14][15][16]. In the original RISM-SCF, electrostatic interactions were approximated with point charge interactions. The point charges were fitted to reproduce the electrostatic potential (ESP) that is computed using the QM approaches [11,12]. Although the fitting approach has been employed in J 2021, 4 various atomic charge calculations, unreasonable charge values were assigned for buried atoms [14,17]. Because of the instability in the charge fitting, the polarization of the solute molecules was enhanced in polar solvents. The fitting problem was overcome using the SED, and the SED was introduced into the RISM-SCF framework. As shown in previous studies, the new method (RISM-SCF-cSED) gave reasonable results even for polar solvents, such as ionic liquids [18][19][20], dimethyl sulfoxide (DMSO) [6], and water [5,[21][22][23][24][25][26].
This paper reports the validity of RISM-SCF-cSED by computing the absorption energy of 5-(dimethylamino)-2,4-pentadienal (DAPDA) in solution. This is a good example to show the validity of the method because the absorption energy of DAPDA has been obtained experimentally for a variety of solvents.

Methods
In RISM-SCF-cSED, the electron density of the solute molecule ρ(r) was approximated using the auxiliary basis sets (ABSs) { f i (r)}, as follows: where d are the expansion coefficients and are determined so that the ESP computed withρ(r) reproduces the ESP computed with ρ(r). The electrostatic potential around each atomic site can be defined usingρ(r). The ground state free energy of RISM-SCF-cSED was defined using the following equation [12,15]: where E solu [G] and ∆µ [G] are the solute energy and solvation free energy at the ground state, respectively. The RISM-SCF-cSED was developed by evaluating E solu [G] with various quantum chemical approaches [5,13,15,25,27,28]. When the density functional theory (DFT) is employed, (2) is given by where H core and F are the core Hamiltonian and the Fock matrix defined in the gas phase. The solvated Kohn-Sham equation can be obtained by taking the derivative of (3) with respect to the molecular orbital coefficients C. The free energy gradient was also derived [12,15,28] by taking the derivative of (3) with respect to the atomic coordinates. When calculating the excited state in solution, the dynamics of the solvent molecules in excitation must be considered. For example, in the absorption energy calculations in solution, there is no time for solvent molecules to relax completely around the solute molecules. The excitation process with the RISM was treated by fixing the solvation structure determined at the ground state [5,26,27,29]. The energy in the excited state was defined as where d [κ] is the fitting coefficients in the κ state, and V [κ] is the electrostatic potential on the ith ABS induced by solvent molecules [13,16,30]. ∆µ [G] in (2) was computed using the following equation: where ρ solv s is the number density of solvent at s site; k B is the Boltzmann factor; T is the temperature. {h αs } and {c αs } are the total and direct correlation functions, respectively, and were computed by coupling the following equations, where φ αs (r) is the site-site potential, {ω αγ } is the intramolecular correlation function of the solute molecule, and {χ ts } is the pure solvent site density pair correlation function. The asterisk denotes a convolution integral.

Computational Details
The geometries of molecules in a solution were optimized using the second order Møller-Plesset perturbation theory (MP2) coupled with RISM-SCF-cSED [28]. Figure 1 presents the atomic labels of DAPDA. The intruder state avoidance multireference Møller-Plesset perturbation theory (ISA-MRMP2) [31,32] was used for the absorption energy calculation, and the energy denominator shift was fixed to 0.02. Eight active electrons in seven active orbitals were selected in the ISA-MRMP2 calculations. The usefulness of the ISA-MRMP2 in the absorption energy calculation was determined using time-dependent density functional theory (TD-DFT) with CAM-B3LYP functional [33]. The employed basis set was may-cc-pV(T+d)Z [34] for geometry optimization and energy calculations. In this study, five solvent molecules were selected: tetrachloromethane (TCM), acetone (ACT), acetonitrile (ACN), dimethyl sulfoxide (DMS), and water (WAT). The pure solvent site density pair-correlation functions were calculated using ex-RISM at 300 K [9]. The non-electrostatic site-site potentials were computed using the Lennard-Jones potential, and the parameters were adopted from the OPLS-AA force field [35].
All calculations were performed using the GAMESS program modified by us [36].

Results and Discussion
To discuss the excitation character of DAPDA, natural orbitals at the first excited state were calculated using ISA-MRMP2 and TD-DFT. In Figure 2, the almost singly occupied natural orbitals obtained with ISA-MRMP2 and TD-DFT are shown. Although there are small differences, both results show that the computed excitation had π → π * character.  Figure 3 summarizes the π → π * excitation energies in solvents (TCM, ACT, ACN, DMS, and WAT) computed suing TD-DFT and ISA-MRMP2. Experimental data were used for a comparison [37]. The experimental data showed that the excitation energies decreased with increasing solvent polarity. The solvatochromism in the excitation energy was reproduced even with TD-DFT. However, the TD-DFT overestimated the excitation energies and the difference was ∼0.8 eV. The large error was reduced significantly with ISA-MRMP2. By coupling ISA-MRMP2 with RISM-SCF-cSED, the average error was less than 0.1 eV. The results showed that the multireference perturbation theory coupled with RISM works well for the excitation energy calculations of π-conjugated molecules in solution, as shown in previous studies [5,26]. The origin of solvatochromism can be explained by resonance structures. Scheme 1 presents the resonance structures (neutral and zwitterionic forms). The zwitterionic form has the charge transfer character. The neutral form is dominant in the ground state, whereas the population of the zwitterionic form increases in the excited state because of π → π * excitation. Because the zwitterionic form has a large dipole moment compared to the neutral form, the dipole moment in the excited state is greater than the ground state. Thus, the excitation energy decreases with increasing solvent polarity. The population change in the zwitterionic form can be confirmed using the Wiberg bond index [38]. Figure 4 shows the Wiberg bond index calculated in the excited state for TCM and WAT. In the case of TCM, the bond indices between C 1 − C 2 , C 3 − C 4 , and C 5 − O are greater than those of N − C 1 , C 2 − C 3 , and C 4 − C 5 , which shows the neutral form character. In the case of WAT, however, the bonds between N − C 1 , C 1 − C 2 , C 2 − C 3 , C 3 − C 4 , and C 4 − C 5 have a similar bond index, and the bond index of C 5 − O is small compared to that in TCM. The bond index difference shows that the population of the zwitterionic form in the excited state for WAT is larger than that for TCM.
Using RISM, the solvation structure can be discussed with RDFs. Figure 5 presents the RDFs between the solute oxygen and water hydrogen sites (g OH (r)) and between the solute oxygen and water oxygen sites (g OO (r)). Sharp peaks were observed around r = 2.0 Å in g OH (r) and r = 3.0 Å in g OO (r), respectively. The former corresponds to hydrogen bonding between the carbonyl oxygen and water hydrogen sites. Because the distance between the peak of g OH (r) and g OO (r) is ∼1.0 Å and the bond distance between O and H in water is 1.0 Å, the RDFs shows that the carbonyl oxygen, water hydrogen, and water oxygen sites are in a straight line, as shown in the inset in Figure 5, indicating strong hydrogen bonding between DAPDA and water. The interaction between the carbonyl group and the solvent can explain the solvatochromism in Figure 3.

Conclusions
The excitation energies of DAPDA in five solvent systems obtained using RISM-SCF-cSED were discussed. Using ISA-MRMP2, the average error between the computed and experimental data in the excitation energy was less than 0.1 eV. The solvatochromism in the excitation energies was discussed with the resonance structures, Wiberg bond index, and RDFs. Because of the strong interaction between the carbonyl group of DAPDA and solvent molecules, the zwitterionic form was enhanced in the excited state, leading to a decrease in the excitation in the polar solvent.
This study shows that RISM-SCF-cSED provides reliable results for excitation energy calculations in solution using sophisticated quantum chemical methods, such as ISA-MRMP2. Moreover, the solvation structures can be discussed via correlation functions, such as radial distribution functions, because RISM was derived based on statistical mechanics. We believe that RISM-SCF-cSED is crucial in excited state calculations in solution.
Funding: This research was funded by JSPS KAKENHI Grant Number JP19K05367.

Data Availability Statement:
The authors confirm that the data supporting the findings of this study are available within the article.

Conflicts of Interest:
The authors declare no conflict of interest.