Efficient Calculation of Absorption Spectra in Solution: Approaches for Selecting Representative Solvent Configurations and for Reducing the Number of Explicit Solvent Molecules

1.   Introduction

The increasing demand for energy, the depletion of fossil fuels, and the harmful impacts of fossil fuel combustion have created an urgent need to find alternative energy supplies, which are environmentally friendly and renewable. Among all renewable energy technologies, solar energy appears most promising, since the amount of the energy emitted from the Sun onto the Earth per year is several thousands times higher than the global energy requirement ¹. First generation silicon solar cells are based on solid-state junctions and have a conversion efficiency of about 10% to 28% ². However, these cells are expensive, energy-consuming to manufacture, and their rigid shape limits them to specific applications ³. In 1991, Graetzel's research group reported the first dye-sensitized solar cell (DSSC) ⁴. Since then, DSSCs have attracted significant attentions because of their low cost and environmentally friendly properties, and have the potential to replace traditional silicon solar cells. In a DSSC, dyes are often adsorbed onto a non-crystalline titanium dioxide (TiO₂) electrode ⁵. It has been shown that the formation of dye-aggregates on the TiO₂ surface will affect the solar cell efficiency ^6-8. In general, dye-aggregation is an undesirable phenomenon in DSSCs. However, in some limited cases, photocurrent generation can be enhanced by careful control of aggregation ⁸.

Measuring the electronic absorption spectra is crucial in designing new dyes, since it can provide guidance to find the improved sensitizers. A chromophore is the part of the dye that can transit, reflect, or absorb certain wavelengths of visible light. Due to its simplicity in structure, many dye studies begin with a single chromophore. In the gas phase (at finite temperature), the absorption spectrum of the chromophore usually exhibits several discrete peaks. However, when measured in solution, these peaks often broaden considerably and even merge into a single continuous band due to the solvent effect. Meanwhile, the center of the peak may also be shifted, and this phenomenon is called solvatochromism. In addition, it has been shown that when a chromophore is dissolved into a multi-component solvent mixture, some property changes of the chromophore, such as the solvatochromic shift, may not have a linear dependence on the composition of the mixture, and this phenomenon is called "preferential solvation" ^{9, 10}. This presents a challenge to the study of chromophore solvation, since it may indicate that the local environment of the chromophore is not uniformly distributed. In the work of El Seoud ⁹, three possible explanations for preferential solvation are proposed. The first one is the "dielectric enrichment", which means the local environment of the chromophore is enriched by the solvent with a higher dielectric constant. The second is the presence of specific interactions between chromophore and solvent. The third explanation is solvent microheterogeneity. El Seoud derived an empirical equation to correlate the solvatochromic shift with the binary solvent composition. However, this study was not based on a molecular level description, so the underlying reason of preferential solvation was not explained and no universal approach to predict solvatochromic shift was suggested.

Computer simulations provide a promising alternative to empirical methods in studying solvation problems. The existing methods of treating solvents can be divided into two main categories. The first is to treat solvent molecules explicitly. An example is the use of full quantum mechanical (QM) methods. These methods explicitly model the chromophore and surrounding solvent molecules. One often needs to be careful in selecting the number of solvent molecules, since the amount necessary to mimic a real solution may be different from case to case. Meanwhile, a Monte Carlo (MC) or a molecular dynamics (MD) simulation is often carried out to sample the phase space. If such QM calculations are affordable, they should provide the highest accuracy. However, the large system size required to reproduce a real solution and large number of configurations necessary to cover the entire phase space often make full QM methods implausible. Recent advances in full QM methods try to partition the large solution system into molecular fragments, and the computational cost is reduced to some extent ^11-14. However, the problem of configurational sampling still exists. An alternative way to deal with the solvent distribution is to combine the advantages of molecular mechanics (MM) and quantum mechanics. One such method is the sequential MM/QM method, which was introduced by Coutinho et al. ^15-17. In their approach, a long simulation using MM force fields for a system containing a chromophore solvated in explicit solvent molecules is performed to generate uncorrelated configurations. Then, about 100 configurations of the chromophore and only the solvent molecules in the first solvation shell are selected, and electronic structure calculations are carried out on each of these configurations. A similar approach suggested by Aidas et al. ¹⁸. increased the number of configurations and solvent molecules by a factor of 10, but the following calculations of the absorption spectrum had to be limited to a QM/MM approach. This illustrates a trade-off between choosing either a large system size and relatively inexpensive QM calculations or only including the local solvent molecules around the chromophore and higher level QM calculations.

The second category is to treat the solvents implicitly. These methods avoid two issues of full QM methods by replacing the solvent molecules with a continuum medium and the solute being considered to be in a cavity of the solvent ¹⁹. Despite their advantages in computational cost, implicit solvation methods do not consider specific solvent-solute interactions, thus failing to explain fully some observations such as preferential solvation ²⁰. In some recent studies, a new scheme that is able to combine both explicit methods and implicit methods was developed. In this scheme, the chromophore molecule is surrounded by only a few explicit solvent molecules in the QM region and then embedding the whole QM region in a continuous medium to account for the contributions from more distant solvent molecules ^{21, 22}. Results show that employing these "mix-discrete" approaches can yield good agreement with experiments. Overall, these studies have demonstrated that using pure implicit solvation methods cannot achieve satisfactory results. At least some of the local solvent molecules need to be treated explicitly, which generates two challenges; determining which solvent molecules should be selected and how to sample phase space efficiently.

In this work, we develop a method that explicitly treats the local solvation environment of a chromophore while still efficiently handling the problem raised by inefficient sampling of the phase space. In previous studies, too many solvent configurations (100–1000) were selected so that high-level QM methods would be prohibitively expensive in the later calculation of the absorption properties. Our goal in this work is to develop a fitness function, based on which we can select a subset of solvent configurations from MC simulations and still reproduce the excitation energy distribution of the entire ensemble. This article is organized as follows: the next section describes the computational details. The third section reports the results obtained from our new method. The last section provides conclusions of this work.

2.   Computational details

MC simulations were carried out with the "Monte Carlo for Complex Chemical Systems-Minnesota" (MCCCS-MN, Version 16.1) software suite developed in house by the Siepmann research group. Simulations for systems containing one chromophore molecule solvated in 500 solvent molecules were carried out in the isothermal-isobaric (NpT) ensemble at T = 298 K and p = 0.1 MPa. The chromophore under investigation in this study is 1-methyl-8-oxyquinolinium betaine (QB), and its structure is shown in Fig. 1. The solvent consisted of either water or acetonitrile (ACN). Since there was only one QB molecule, it was fixed at the center of the simulation box, and translational and rotational moves were only carried out on solvent molecules. The probability for volume moves (involving scaling of the solvent center-of-mass coordinates) was set to yield approximately one accepted volume move per MC cycle (MCC, one MCC consists of N Monte Carlo moves, where N = 501 is the number of molecules in the system). The remaining moves were evenly divided between rigid-body translations and rotations of the solvent molecules. The maximum displacements were adjusted during the equilibration period to ensure a 50% acceptance rate for the latter two move types.

图 1

图 1  (left) Structure of QB and atom numbering where cyan, white, red, and blue spheres indicate carbon, hydrogen, oxygen, and nitrogen atoms, respectively; (right) chemical structure of QB and ground state (GS) and excited state (ES) partial charges of oxygen and nitrogen atoms calculated at the TD-CAM-B3LYP/6-311+G(d, p) level of theory.

Figure 1.  (left) Structure of QB and atom numbering where cyan, white, red, and blue spheres indicate carbon, hydrogen, oxygen, and nitrogen atoms, respectively; (right) chemical structure of QB and ground state (GS) and excited state (ES) partial charges of oxygen and nitrogen atoms calculated at the TD-CAM-B3LYP/6-311+G(d, p) level of theory.

下载: 全尺寸图片 幻灯片

The atoms on the aromatic ring of the QB molecule were modeled by the TraPPE-EH force field ^{27, 28}. Because the TraPPE-EH model does not specify transferable partial charges for aromatic compounds, the following procedure was used to obtain the partial charges: first, the QB molecule was optimized at the M06-2X/6-31+G(d, p) level of theory in implicit 1-octanol using SM8 ²⁹ in Gaussian 09 ³⁰. Then the CM5 charges were obtained by using CM5PAC ^{31, 32}. The complete set of force field parameters for QB is provided in the Supporting Information. The methyl group in QB and in the acetonitrile molecule was described by the TraPPE-UA force field ^{33, 34}, and the TIP4P model ³⁵ was used for water. The Lorentz-Berthelot combining rules ²⁶ were used for the interactions between unlike non-bonded atoms. Periodic boundary conditions were used and only the interactions between the nearest images were calculated (minimum image convention). A spherical truncation at 1.4 nm and analytical tail corrections ^{23, 24} were used for the Lennard-Jones interactions, and the Ewald summation method ^{24, 25} was used for Coulomb interactions. Since the methyl group is represented by only one bead in the TraPPE-UA model, but the QM calculations require explicit hydrogen atoms, three virtual hydrogen atoms (consistent geometry but without Lennard-Jones site or partial charges) were added to each methyl group.

16 independent simulations were carried out, with each one consisting of 15000 MCCs for the equilibration period, followed by 100000 MCCs for the production period. Configurations were recorded every 1000 production cycles, generating 1600 uncorrelated configurations for the subsequent analysis.

The excitation energy for each configuration (QB and 500 solvent molecules without periodic boundary condition, but with the location of each complete solvent molecule determined by its center-of-mass position in the supercell) was calculated through the ZINDO method ^{36, 37} using Gaussian 09. The absorption spectra were generated as follows: first, every excitation energy obtained from the ZINDO calculation was broadened with a Lorentz distribution ³⁸,

where E₀ is the single-point excitation energy obtained from the ZINDO calculation, and γ is the scale parameter that is set to 0.005 in the present work. Then, all individual Lorentzian peaks were added together and normalized by the number of configurations and the total area under the curve to obtain the full spectrum of the solvated chromophore molecule at the ZINDO level.

3.   Results and discussion

3.1   Radial distribution function and ZINDO excitation energy spectra

As we have mentioned, an important aspect of the sequential MM/QM method is the validity of the molecular mechanics simulation. Before calculation of the absorption spectra, we first analyzed the simulation trajectories. The MC simulations were carried out in the isobaric-isothermal ensemble and, hence, the box lengths fluctuated throughout the trajectories. The average length of simulation boxes were 2.4814 ± 0.0001 nm and 3.5052 ± 0.0001 nm for water and acetonitrile, respectively, and the fluctuations in lengths were less than 2%. The most polar atom in QB is the oxygen atom. For the two solvents, the O(QB)-O(water) radial distribution function (RDF) in water and the O(QB)-Me(ACN) RDF in acetonitrile and their corresponding number integrals were analyzed. As can be seen from Fig. 2, the RDF with water exhibits two peaks at 0.29 and 0.46 nm, and the number of water molecules in the first solvation shell of O(QB) is about 3. The RDF with acetonitrile also shows two peaks, but at considerably large distances of 0.34 and 0.76 nm, and the number of acetonitrile molecules in the first solvation shell is about 5. Both RDFs indicate specific interactions of the solvent molecules with the QB molecule.

图 2

图 2  (top) O(QB)-O(water) radial distribution function (RDF, left) and corresponding number integral (nint, right); (bottom) O(QB)-Me(ACN) radial distribution function (RDF, left) and corresponding number integral (nint, right) (1 = 0.1 nm).

Figure 2.  (top) O(QB)-O(water) radial distribution function (RDF, left) and corresponding number integral (nint, right); (bottom) O(QB)-Me(ACN) radial distribution function (RDF, left) and corresponding number integral (nint, right) (1 = 0.1 nm).

下载: 全尺寸图片 幻灯片

The ZINDO absorption spectra obtained from 1600 configurations for each solvent are shown in Fig. 3. For both water and acetonitrile, the distributions are continuous bands, which range from 1.90 eV to 2.38 eV. The average ZINDO excitation energy for QB in acetonitrile is 2.12 eV, while that for solvation in water is 2.14 eV. Both distributions show a regular peak shape; i.e., an indication that the 1600 configurations collected from the MM simulations is good representatives of the entire ensemble. Since it is not feasible to afford high-level quantum calculations on all 1600 configurations, it is necessary to select a subset of configurations that can still reproduce the distributions shown in Fig. 3. In previous work by Canuto et al. ^15-17, the authors only calculated the average value of excitation energies from selected configurations. However, a single property such as the average is often not sufficient enough to depict the entire absorption spectrum. Sometimes, due to specific solvent-solute interactions, the absorption spectra may not be symmetric ³⁹. In our work, we are not only interested in the average values of excitation energies, but also reproducing the complete spectra shown in Fig. 3.

图 3

图 3  ZINDO excitation energy spectra of the QB chromophore in water and acetonitrile based on 1600 configurations (with a Lorentzian broadening of γ = 0.005).

Figure 3.  ZINDO excitation energy spectra of the QB chromophore in water and acetonitrile based on 1600 configurations (with a Lorentzian broadening of γ = 0.005).

下载: 全尺寸图片 幻灯片

3.2   Development of the fitness function

In the next step, a fitness function was developed to correlate some molecular mechanics properties with the ZINDO excitation energies. The vertical excitation energy is related to the energy difference between the ground state and the excited state of the molecule, while holding chromophore and solvent configuration fixed. Since the solute-solvent Coulomb interactions are responsible for up to 95% of the total solvatochromic shift ⁴⁰, a direct way to develop a fitness function is to investigate the first-order electrostatic interactions between the ground state partial charges of QB and the partial charges of the solvent molecules, and similarly for the excited state partial charges of QB. The ground state and excited state partial charges of QB were obtained from the CM5 calculation at the TD-CAM-B3LYP/6-311+G(d, p) level of theory in the gas phase. In Fig. 1, we highlighted the partial charges of oxygen and nitrogen atoms in QB. (For a complete set of partial charges, please refer to the Supporting Information.) Although the bonding arrangement in QB supports a zwitterionic character, the N atom is found to carry a negative partial charge regardless of QB being in the ground or excited state. Similarly, for both ground and excited states, the partial charge of the O atom is considerably smaller in magnitude than the expected formal charge of -1|e|.

The first-order electrostatic interaction energy between the ground state partial charges of QB and the solvent partial charges is calculated from

where n is the number of beads in each solvent molecule, N is the number of solvent molecules in the simulation box, q_{i, GS} and q_j are the ground-state partial charges of QB and the solvent's partial charges, respectively, and r_ij is the distance between the corresponding solute site and the solvent site. The first summation is over all interaction sites of QB, and the second summation is over all interaction sites of all solvent molecules. The first-order electrostatic interaction energy between QB and the solvent molecules is calculated in an analogous way,

where q_{i, ES} is the excited-state solute charge. These two quantities were investigated to see if any of them alone could be treated as a good fitness function. To this extent, a simple least-squares linear regression method was utilized to correlate these two quantities with the ZINDO excitation energy.

where E_pred is the predicted value of the ZINDO excitation energy, V refers to either V_{Coul, GS} or V_{Coul, ES}, and α₁ and α₂ are the regression coefficients. These two coefficients were obtained by the following procedure: first, 400 configurations were used as a training set. Their V_{Coul, GS} and V_{Coul, ES} were calculated and then correlated with the corresponding ZINDO excitation energies by using the linear regression method. The correlations between predicted ZINDO excitation energies E_pred and real ZINDO excitation energies E_ZINDO are shown in Table 1 and Fig. 4. For both of the solvents, the fitness functions based on V_{Coul, GS} have reasonably high correlation coefficients, that are both larger than 0.7. On the other hand, the fitness functions based on V_{Coul, ES} have much smaller correlation coefficients compared with V_{Coul, GS}. The R value of 0.346 for solvation in water is already quite small, but the R value of 0.103 for acetonitrile is close to zero, i.e., the two properties are almost uncorrelated (see also Fig. 4). As a result, a fitness function using only V_{Coul, GS} may yield improved performance over a random selection of configurations, while V_{Coul, ES} alone is a poor predictor.

表 1

表 1  Correlation coefficients (R) between E_pred and E_ZINDO when using V_{Coul, GS} and V_{Coul, ES}

Table 1.  Correlation coefficients (R) between E_pred and E_ZINDO when using V_{Coul, GS} and V_{Coul, ES}

下载: 导出CSV

Solvent	VCoul, GS	VCoul, ES
water	0.739	0.346
ACN	0.715	0.103

图 4

图 4  Correlation between E_pred and E_ZINDO using fitness functions based on V_{Coul, GS} and V_{Coul, ES}.

Figure 4.  Correlation between E_pred and E_ZINDO using fitness functions based on V_{Coul, GS} and V_{Coul, ES}.

下载: 全尺寸图片 幻灯片

Next, these two quantities were combined into one fitness function with the following form,

where β₁, β₂, and β₃ are three coefficients. If there are no solvent molecules around the chromophore, which corresponds to V_{Coul, GS} = 0 and V_{Coul, ES} = 0, it is obvious that β₃ is equivalent to the gas-phase ZINDO excitation energy of the chromophore, so the equation can be rewritten as

where E_{g, ZINDO} is the gas phase ZINDO excitation energy of QB. The other two coefficients, β₁ and β₂ were calculated by taking 400 configurations out of the 1600 generated ones as the training set and then applying a multidimensional linear regression method, that uses ordinary least squares method to solve for the two coefficients. The numerical results are provided in Table 2, and the correlation is illustrated in Fig. 5. As one can see, after including both quantities in the fitness function, there is a significant improvement in correlation coefficients for both water and acetonitrile. The R value for acetonitrile increases from 0.715 to 0.959. The data points for acetonitrile have a narrow distribution around the linear fit line, and the difference in scatter between the training set and the other 1200 configurations is negligible. These indicate that E_pred is an excellent fitness function. Meanwhile, for water, the increase and absolute value of the correlation coefficient are not as large as for acetonitrile, with the R value improving from 0.739 to 0.875. It can be seen that, for solvation in both water and acetonitrile, the ratio between the first two coefficients of the fitness functions, β₁/β₂, is close to -1. This indicates that our fitness function may have a clear physical meaning—it tells us that the solvatochromic shift mainly arises from the difference in first-order electrostatic interaction energies between solvent with ground state chromophore and solvent with excited state chromophore. It should also be noted that the first two coefficients β₁ and β₂ are about 40% smaller in magnitude for solvation in water than for acetonitrile. If we want this method to be universal for all solvents and also for solvent mixtures, then only one set of universal coefficients can be used. Since our interest is to select representative configurations, the prediction of the relative order of the excitation energy is more important than the prediction of the absolute values of the excitation energy. A feasible approach to make coefficients universal is to replace β₁ and β₂ with the average values obtained here for water and acetonitrile. Note that after this replacement, E_pred is no longer directly a prediction of the ZINDO excitation energy, but it can still work as a "rank" to distinguish the relative order of ZINDO excitation energies among different configurations. As can be seen from Table 2 and Fig. 5, the correlation coefficients for both solvent water and acetonitrile only decrease by less than 0.01. This indicates that even if we set β₁ = -10.3 and β₂ = 9.90, it has little detrimental effect on the rank of the configurations. However, it is obvious that the absolute value of E_pred is changed by a certain amount, but this is of little interest to us. We demonstrate that having β₁ and β₂ set to their averages can still yield satisfactory results and it sheds light on the development of a universal fitness function for all solvents.

表 2

表 2  Correlation and coefficients of fitness function based on both V_{Coul, GS} and V_{Coul, ES}.

Table 2.  Correlation and coefficients of fitness function based on both V_{Coul, GS} and V_{Coul, ES}.

下载: 导出CSV

Solvent	R	β1	β2	Eg, ZINDO	β1/β2
water	0.875	−8.52	8.20	1.77	−1.04
ACN	0.959	−12.1	11.6	1.77	−1.05
water	0.875	−10.3	9.90	1.77	−1.04
ACN	0.953	−10.3	9.90	1.77	−1.04

图 5

图 5  Correlation between E_pred and E_ZINDO using the fitness function combining V_{Coul, GS} and V_{Coul, ES} for solvent water and acetonitrile.

Figure 5.  Correlation between E_pred and E_ZINDO using the fitness function combining V_{Coul, GS} and V_{Coul, ES} for solvent water and acetonitrile.

下载: 全尺寸图片 幻灯片

3.3   Selection of representative configurations

After the rank of all configurations was obtained, a subset of configurations from the 1600 for solvent water and acetonitrile were selected and assessed for their ability to reproduce the complete spectra shown in Fig. 3. Three biased selection schemes were developed, and their results were compared with a random selection. The first two steps of all three biased selection schemes are the same. First, the 1600 configurations obtained from the MC simulations were sorted based on E_pred. Second, they were divided into N_config sections containing each an equal number of configurations. N_config is the number of configurations that will be passed on to the subsequent high-level quantum mechanical calculations. The last step differs between the three biased schemes. For the scheme denoted as S-mean, the configurations with the E_pred value closest to the mean E_pred values of each specific section were selected as representative configurations. For the scheme denoted as S-median, and the configurations at the mid-point of each section were selected. If this section has an even number of configurations, we randomly selected the configuration just before or after the mid-point. For the third biased scheme denoted as S-random, one configuration was randomly selected from each section. The spectra of the subsets of configurations obtained by these three biased selection schemes were compared with the one of an entirely random selection, and their similarity indices with the spectra of the complete set of 1600 configurations were calculated. It should be noted that the selection based on E_pred values of specific configurations only indirectly accounts for the volume fluctuations encountered in the MC trajectories, i.e., no attempt is made to bias the selection in a manner so that the average box length of the selected configurations matches the ensemble average.

In our opinion, assessing the similarity of two distributions by comparing only their mean values and standard deviations would be insufficient, since these two properties cannot tell the subtle differences between two distributions. Hence, a new similarity index, D, that is based on Kolmogorov-Smirnov test, was introduced,

where F_n(x) and G_n(x) are cumulative distribution functions, and sup is the supremum of this set. In the case of discrete variables, it becomes the largest value in this set. The mathematical meaning of D is the largest difference between cumulative probabilities of two distributions. A smaller D value means that two distributions are more similar. In order to have a meaningful D value, the two cumulative distributions F_n(x) and G_n(x) should have equal number of data points and the same range in x.

We tested the selection of 5, 10, 20, 40, 80, 100, 160 and 200 representative configurations and used the same procedure as when treating 1600 configurations to plot the spectra. This ensures that the distribution of the subset and the distribution of 1600 configurations are comparable if D was utilized. For selecting a certain number of representative configurations for each solvent, the three biased selection schemes and the random selection were carried out, and their D values were calculated. Note that since S-median, S-random, and complete random do not yield a unique selection, the D values of these schemes were obtained by averaging over 100 test cases. The error bars for these schemes reflect the standard deviation of the 100 test cases. A control experiment was also constructed, which was denoted D^* here: even if the fitness function is perfect, there would be some information lost when the number of configurations is reduced from 1600 to 200, 160, 100, 80, 40, 20, 10, and 5. Consequently, the lower bound of D would not be exactly 0. The values of D^* were calculated by selecting configurations directly based on E_ZINDO as if we know the excitation energies of all the configurations beforehand, and this gives us the base level of the D value. The previously discussed selection methods were performed on the ranked configurations, and the results are shown in Fig. 6. Note that S-mean scheme performs in many cases quite similar to the S-median schemes, but in a few cases performs poorly, so only the S-median scheme is discussed here. (Results for the S-mean scheme are provided in the Supporting Information.) First, one should notice that both D and D^* increase as one decreases the number of configurations selected, indicating that more information is lost when only a small fraction of the total number of configurations was selected. As indicated by Fig. 6, in most cases, the two biased selection schemes yield significantly smaller D values, indicating that they clearly outperform the random selection and are able to select more representative configurations. In addition, the statistical uncertainty in D for the biased selections, especially when the number of configurations selected is small, is less than that of the random selection. This implies that the smart selection schemes can consistently select more representative configurations over the 100 test cases. In order to quantitatively compare the performance of the two biased selection schemes, the ratio D_biased/D_rand is also shown in Fig. 6. The first thing to notice is that the D_biased/D_rand ratios are consistently smaller than 1 for both biased selection schemes, i.e., they indeed have better performance than the random selection. If the results of the two biased selections are compared, one should notice that the S-median scheme performs either equally well or, in most cases, better than the S-random scheme. The average D_biased/D_rand ratios for water and acetonitrile for the S-median scheme are 0.63 and 0.53, respectively. For the S-random scheme, the corresponding ratios are 0.71 and 0.59, respectively. The smaller values observed for solvation in acetonitrile versus water are an outcome of the better performance of the prediction scheme (see comparison of R values in Table 3). The smaller values found for S-median (and also S-mean) versus S-random indicate that ranking the configurations in each subset provides some benefit over just sorting into subsets.

图 6

图 6  D and D^* for two biased selection schemes based on predicted and pre-known excitation energies, and D_biased/D_rand.

Figure 6.  D and D^* for two biased selection schemes based on predicted and pre-known excitation energies, and D_biased/D_rand.

下载: 全尺寸图片 幻灯片

表 3

表 3  Correlation and coefficients of fitness function based on selected solvent molecules.

Table 3.  Correlation and coefficients of fitness function based on selected solvent molecules.

下载: 导出CSV

Solvent	R	β1	β2	Eg, TD-DFT	β1/β2
water	0.918	−10.3	9.90	2.05	−1.04
ACN	0.887	−10.3	9.90	2.05	−1.04
water	0.926		use Eq. (8)
ACN	0.898		use Eq. (8)

A comparison of the D^* value for each selection scheme shows that the S-median scheme yields a significantly smaller D^* value when the number of configurations selected is relatively small. This also demonstrates that S-median is a better selection scheme than S-random. In addition, we want to emphasize that the D^* values for selections based on pre-known excitation energies are much smaller than the D values based on the predicted excitation energies. Thus, additional gains in the performance of the selection schemes should arise from improving the fitness function.

Finally, the average ZINDO excitation energies of selected configurations are shown in Fig. 7. As one can see, for both solvents and the three selection schemes, the differences between average ZINDO excitation energies are less than 0.01 eV. Thus, one can hardly distinguish performance by only looking at the average values. However, it should be noted that the uncertainty in average excitation energies for the S-median scheme is the smallest among the three, meaning it can provide the most consistent results over different test cases.

图 7

图 7  Average ZINDO excitation energies of selected configurations and uncertainties.

Figure 7.  Average ZINDO excitation energies of selected configurations and uncertainties.

下载: 全尺寸图片 幻灯片

3.4   Selection of solvent molecules

A simple criterion to reduce the number of solvent molecules was investigated, that is the distance to the mid point of the C(2)―C(3) bond in QB. For solvation in acetonitrile, 50 solvent molecules that have the smallest distances to the mid point were selected for explicit inclusion, and 160 water molecules were also chosen using the same method to represent the local solvation in water. The remainder of the solvent molecules was included as background using the partial charges from the MM representation, i.e., the acetonitrile molecule was represented by three partial charges at the positions of the CH₃ group, C and N atoms, and the water molecule was replaced by three partial charges at the positions of two H atoms and the M site. The QM calculations on these smaller systems were performed at the TD-CAM-B3LYP/3-21G level of theory, and E_pred was correlated to E_TD-DFT according to Eq. (6) with one revision of changing E_{g, ZINDO} to E_{g, TD-DFT} (calculated at the same level of theory). The results are provided in Table 3 and Fig. 8. It should be noticed that although the fitness function was developed based on ZINDO excitation energies, it performs remarkably well when applied to the TD-DFT excitation energies. It indicates that the fitness function may be transferable regardless of the level of theory used. Although the explicit solvent region for water is reduced from 500 molecules to 160, the correlation of the fitness function is raised from 0.870 to 0.918, and that for acetonitrile drops from 0.953 to 0.887. These indicate that reducing solvent molecules does not undermine the overall behavior of the fitness function. Since the time spent on calculating the perturbation by the background charges is negligible compared to that spent on explicit molecules with optimizable charge distributions, one can save a large amount of computational time by only explicitly representing a small region of solvent molecules and then immersing the entire system into the long-range background charges. Finally, it should be emphasized that only force field charges were used here as the background in our test calculations. The next step for improvement may be replacing them with ab initio charges. The correlation between ZINDO excitation energies and TD-DFT excitation energies were also investigated. Their correlation was described by the following equation,

图 8

图 8  Correlation between E_pred and E_TD-DFT for water and acetonitrile using Eqs. (6) and (8).

Figure 8.  Correlation between E_pred and E_TD-DFT for water and acetonitrile using Eqs. (6) and (8).

下载: 全尺寸图片 幻灯片

According to Table 3 and Fig. 8, E_pred calculated from Eq. (8) has a high correlation with E_TD-DFT, which means that the excitation energies calculated from ZINDO and TD-DFT method are highly correlated. Thus, ZINDO excitation energy itself can serve as a predictor for TD-DFT (or even higher level methods).

4.   Conclusions

A new scheme that can increase the accuracy of selecting representative solvent configurations has been developed in this work. Unlike previous studies, which only select configurations at regular intervals from the MM trajectory, our new scheme specifically targets those configurations that are the representative of the distribution found in the entire ensemble. Since the number of configurations used in the later QM calculation is reduced, we expect a great improvement in the calculation efficiency. In this work, we investigated the solvatochromism of chromophore QB in water and acetonitrile. In the first part, the fitness functions were developed separately for solvation in water and acetonitrile. Although the fitness functions only includes terms that account for the Coulomb interactions between the solute and solvent molecules, the high correlation between E_pred and E_ZINDO indicates that the fitness functions already capture the most important contributions to solvatochromism. The success of developing a common fitness function for solvation in both water and acetonitrile implies that it may be possible to use a universal fitness function for all kinds of solvents and solvent mixtures. We also introduced three biased selection schemes that pick representative configurations based on the newly developed fitness function, and compared them with a random selection. Our results show that sorting all configurations into ranked subsets and randomly selecting one configuration from each subset does already offer some improvement over a completely random selection from the entire set of configurations. However, selecting the configuration closes to the median or mean of each subset provides a far more reliable selection or representative configurations. The benefit of a biased selection increases with a decrease in the target number of selected configurations. Based on the D and D^* similarity indices, we recommend that the S-median scheme should be used. Meanwhile, by applying a distance-based criterion to select explicit solvent molecule for the local environment and replacing the remainder with background partial charges, allows for significant reduction in the computational cost.

Acknowledgment: This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program, under Award No. DE-SC0008666. The authors acknowledge the Minnesota Supercomputing Institute (MSI) at the University of Minnesota for providing computational resources that contributed to this work.

Supporting Information: available free of charge via the internet at http://www.whxb.pku.edu.cn.

References

1. [1]

   Ladomenou, K.; Kitsopoulos, T. N.; Sharma, G. D.; Coutsolelos, A. G. RSC Adv. 2014, 4, 21379. doi: 10.1039/c4ra00985a  doi: 10.1039/c4ra00985a
   

2. [2]

   Green, M. A.; Emery, K.; Hishikawa, Y.; Warta, W.; Dunlop, E. D. Prog. Photovolt Res. Appl. 2015, 23, 1. doi: 10.1002/pip.2637  doi: 10.1002/pip.2637
   

3. [3]

   Graetzel, M. Nature 2001, 414, 338. doi: 10.1038/35104607.  doi: 10.1038/35104607
   

4. [4]

   O'Regan, B.; Graetzel, M. Nature 1991, 353, 737. doi: 10.1038/353737a0  doi: 10.1038/353737a0
   

5. [5]

   Gong, J.; Liang, J.; Sumathy, K. Renew. Sust. Energ. Rev. 2012, 16, 5848. doi: 10.1016/j.rser.2012.04.044  doi: 10.1016/j.rser.2012.04.044
   

6. [6]

   Deing, K. C.; Mayerh ffer, U.; Würthner, F.; Meerholz, K. Phys. Chem. Chem. Phys. 2012, 14, 8328. doi: 10.1039/c2cp40789b  doi: 10.1039/c2cp40789b
   

7. [7]

   Luo, L.; Lin, C. -J.; Tsai, C. -Y.; Wu, H. -P.; Li, L. -L.; Lo, C. -F.; Lin, C. -Y.; Diau, E. W. -G. Phys. Chem. Chem. Phys. 2010, 12, 1064. doi: 10.1039/b919962d  doi: 10.1039/b919962d
   

8. [8]

   Pastore, M.; De Angelis, F. ACS Nano 2010, 4, 556. doi: 10.1021/nn901518s  doi: 10.1021/nn901518s
   

9. [9]

   El Seoud, O. A. Pure Appl. Chem. 2007, 79, 1135. doi: 10.1351/pac200779061135  doi: 10.1351/pac200779061135
   

10. [10]

    Tada, E. B.; Novaki, L. P.; El Seoud, O. A. J. Phys. Org. Chem. 2000, 13, 679. doi: 10.1002/1099-1395(200011)13:11<679::AID-POC299>3.0.CO;2-R  doi: 10.1002/1099-1395(200011)13:11<679::AID-POC299>3.0.CO;2-R
    

11. [11]

    Gao, J.; Zhang, J. Z. H.; Houk, K. N. Accounts Chem. Res. 2014, 47, 2711. doi: 10.1021/ar500293u  doi: 10.1021/ar500293u
    

12. [12]

    Li, S.; Li, W.; Ma, J. Accounts Chem. Res. 2014, 47, 2712. doi: 10.1021/ar500038z  doi: 10.1021/ar500038z
    

13. [13]

    Wang, B.; Yang, K. R.; Xu, X.; Isegawa, M.; Leverentz, H. R.; Truhlar, D. G. Accounts Chem. Res. 2014, 47, 2731. doi: 10.1021/ar500068a  doi: 10.1021/ar500068a
    

14. [14]

    He, X.; Zhu, T.; Wang, X.; Liu, J.; Zhang, J. Z. H. Accounts Chem. Res. 2014, 47, 2748. doi: 10.1021/ar500077t  doi: 10.1021/ar500077t
    

15. [15]

    Coutinho, K.; De Oliveira, M. J.; Canuto, S.Int. J. Quantum Chem. 1998, 66, 249. doi: 10.1002/(SICI)1097-461X(1998)66:3<249::AID-QUA6>3.0.CO;2-V  doi: 10.1002/(SICI)1097-461X(1998)66:3<249::AID-QUA6>3.0.CO;2-V
    

16. [16]

    Jaramillo, P.; Pérez, P.; Fuentealba, P.; Canuto, S.; Coutinho, K. J. Phys. Chem. B 2009, 113, 4314. doi: 10.1021/jp808210y  doi: 10.1021/jp808210y
    

17. [17]

    Barreto, R. C.; Coutinho, K.; Georg, H. C.; Canuto, S. Phys. Chem. Chem. Phys. 2009, 11, 1388. doi: 10.1039/b816912h  doi: 10.1039/b816912h
    

18. [18]

    Aidas, K.; Kongsted, J.; Osted, A.; Mikkelsen, K. V.; Christiansen, O. J. Phys. Chem. A 2005, 109, 8001. doi: 10.1021/jp0527094  doi: 10.1021/jp0527094
    

19. [19]

    Christopher, C. Essentials of Computational Chemistry: Theories and Models; John Wiley & Sons: Chichester, UK, 2013.

20. [20]

    Masunov, A.; Tretiak, S.; Hong, J. W.; Liu, B.; Bazan, G. C. J. Chem. Phys. 2005, 122, 224505. doi: 10.1063/1.1878732  doi: 10.1063/1.1878732
    

21. [21]

    Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2015, 119, 958. doi: 10.1021/jp506293w  doi: 10.1021/jp506293w
    

22. [22]

    Murugan, N. A. J. Phys. Chem. B 2011, 115, 1056. doi: 10.1021/jp1049342  doi: 10.1021/jp1049342
    

23. [23]

    Wood, W. W.; Parker, F. R. J. Chem. Phys.1957, 27, 720. doi: 10.1063/1.1743822  doi: 10.1063/1.1743822
    

24. [24]

    Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: Oxford, UK, 1987.

25. [25]

    Ewald, P. Ann. Phys. 1921, 64, 253. doi: 10.1002/andp.19213690304  doi: 10.1002/andp.19213690304
    

26. [26]

    Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham, A. Intermolecular Forces: Their Origin and Determination; Pergamon Press: Oxford, UK, 1987.

27. [27]

    Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2013, 117, 273. doi: 10.1021/jp307328x  doi: 10.1021/jp307328x
    

28. [28]

    Zhang, L.; Siepmann, J. I. Theor. Chem. Acc. 2006, 115, 391. doi: 10.1007/s00214-005-0073-1  doi: 10.1007/s00214-005-0073-1
    

29. [29]

    Marenich, A. V; Olson, R. M.; Kelly, C. P.; Cramer, C. J.; Truhlar, D. G. J. Chem. Theory Comput. 2007, 3, 2011. doi: 10.1021/ct7001418  doi: 10.1021/ct7001418
    

30. [30]

    Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, Revision D.01;Gaussian Inc.: Wallingford, CT, USA, 2013.

31. [31]

    Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. CM5PAC; Uniersity of Minnesota: Minneapolis, MN, USA, 2011.

32. [32]

    Marenich, A. V.; Jerome, S. V.; Cramer, C. J.; Truhlar, D. G. J. Chem. Theor. Comput. 2012, 8, 527. doi: 10.1021/ct200866d  doi: 10.1021/ct200866d
    

33. [33]

    Martin, M. G.; Siepmann, J. I. J. Phys. Chem. B 1998, 102, 2569. doi: 10.1021/jp972543+  doi: 10.1021/jp972543+
    

34. [34]

    Wick, C. D.; Stubbs, J. M.; Rai, N.; Siepmann, J. I. J. Phys. Chem. B 2005, 109, 18974. doi: 10.1021/jp0504827  doi: 10.1021/jp0504827
    

35. [35]

    Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. doi: 10.1063/1.445869  doi: 10.1063/1.445869
    

36. [36]

    Thompson, M. A.; Zerner, M. C. J. Am. Chem. Soc. 1991, 113, 8210. doi: 10.1021/ja00022a003  doi: 10.1021/ja00022a003
    

37. [37]

    Zerner, M. C.; Loew, G. H.; Kirchner, R. F.; Mueller-Westerhoff, U. T. J. Am. Chem. Soc. 1980, 102, 589. doi: 10.1021/ja00522a025  doi: 10.1021/ja00522a025
    

38. [38]

    Hollas, J. M. Modern Spectroscopy; John Wiley & Sons: Chichester, UK, 2004.

39. [39]

    Voityuk, A. A.; Kummer, A. D.; Michel-Beyerle, M. -E.; Rösch, N. Chem. Phys. 2001, 269, 83. doi: 10.1016/S0301-0104(01)00334-2  doi: 10.1016/S0301-0104(01)00334-2
    

40. [40]

    Lin, Y. L.; Gao, J. J. Chem. Theory Comput. 2007, 3, 1484. doi: 10.1021/ct700058c  doi: 10.1021/ct700058c