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

XUE Bai , CHEN Tiannan ,?, SIEPMANN J. Ilja ,2,＊

1 Department of Chemistry and Chemical Theory Center, University of Minnesota, 207 Pleasant Street SE, Minneapolis, MN

55455-0240, USA.

2 Department of Chemical Engineering and Materials Science, University of Minnesota, 421 Washington Avenue SE, Minneapolis,MN 55455-0132, USA.

? Current address: Department of Computer Science and Engineering, University of Minnesota, 200 Union Street SE, Minneapolis,MN 55455-0154, USA.

Abstract: Dye-sensitized solar cells (DSSCs) are one of the most promising renewable energy technologies. Charge transfer and charge transport are pivotal processes in DSSCs, which govern solar energy capture and conversion. These processes can be probed using modern electronic structure methods. Because of the heterogeneity and complexity of the local environment of a chromophore in DSSCs (such as solvatochromism and chromophore aggregation), a part of the solvation environment should be treated explicitly during the calculation. However, because of the high computational cost and unfavorable scaling with the number of electrons of high-level quantum mechanical methods, approaches to explicitly treat the local environment need careful consideration. Two problems must be tackled to reduce computational cost. First,the number of configurations representing the solvent distribution should be limited as much as possible. Second, the size of the explicit region should be kept relatively small. The purpose of this study is to develop efficient computational approaches to select representative configurations and to limit the explicit solvent region to reduce the computational cost for later (higher-level) quantum mechanical calculations. For this purpose, an ensemble of solvent configurations around a 1-methyl-8-oxyquinolinium betaine (QB) dye molecule was generated using Monte Carlo simulations and molecular mechanics force fields. Then, a fitness function was developed using data from inexpensive electronic structure calculations to reduce the number of configurations. Specific solvent molecules were also selected for explicit treatment based on a distance criterion, and those not selected were treated as background charges. The configurations and solvent molecules selected proved to be good representatives of the entire ensemble; thus, expensive electronic structure calculations need to be performed only on this subset of the system, which significantly reduces the computational cost.

Key Words: Monte Carlo simulation; Chromophore; Spectra; Solution

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 requirement1. First generation silicon solar cells are based on solid-state junctions and have a conversion efficiency of about 10% to 28%2. However, these cells are expensive, energy-consuming to manufacture, and their rigid shape limits them to specific applications3. In 1991, Graetzel’s research group reported the first dye-sensitized solar cell(DSSC)4. 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 (TiO2) electrode5. It has been shown that the formation of dye-aggregates on the TiO2surface will affect the solar cell efficiency6–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 aggregation8.

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 Seoud9, 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 extent11–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.18. 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 solvent19. 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 solvation20. 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 molecules21,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.

Fig. 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 field27,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 SM829in Gaussian 0930. Then the CM5 charges were obtained by using CM5PAC31,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 field33,34, and the TIP4P model35was used for water. The Lorentz-Berthelot combining rules26were 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 corrections23,24were used for the Lennard-Jones interactions,and the Ewald summation method24,25was 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 method36,37using Gaussian 09.The absorption spectra were generated as follows: first, every excitation energy obtained from the ZINDO calculation was broadened with a Lorentz distribution38,

where E0is 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 Rad ial 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.

Fig. 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).

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

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 symmetric39. 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.2 Develop ment 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 shift40, 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, qi,GSand qjare the ground-state partial charges of QB and the solvent’s partial charges, respectively, and rijis 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 qi,ESis 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 Epredis the predicted value of the ZINDO excitation energy, V refers to either VCoul,GSor VCoul,ES, and α1and α2are the regression coefficients. These two coefficients were obtained by the following procedure: first, 400 configurations were used as a training set. Their VCoul,GSand VCoul,ESwere calculated and then correlated with the corresponding ZINDO excitation energies by using the linear regression method. The correlations between predicted ZINDO excitation energies Epredand real ZINDO excitation energies EZINDOare shown in Table 1 and Fig. 4. For both of the solvents, the fitness functions based on VCoul,GShave reasonably high correlation coefficients,that are both larger than 0.7. On the other hand, the fitness functions based on VCoul,EShave much smaller correlation coefficients compared with VCoul,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 VCoul,GSmay yield improved performance over a random selection of configurations, while VCoul,ESalone is a poor predictor.

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

where β1, β2, and β3are three coefficients. If there are no solvent molecules around the chromophore, which corresponds to VCoul,GS= 0 and VCoul,ES= 0, it is obvious that β3is equivalent to the gas-phase ZINDO excitation energy of the chromophore,so the equation can be rewritten as

where Eg,ZINDOis the gas phase ZINDO excitation energy of QB. The other two coefficients, β1and β2were 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 Epredis 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, β1/β2, 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 β1and β2are 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 β1and β2with the average values obtained here for water and acetonitrile. Note that after this replacement, Epredis 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 β1= ?10.3 and β2=9.90, it has little detrimental effect on the rank of the configurations. However, it is obvious that the absolute value of Epredis changed by a certain amount, but this is of little interest to us. We demonstrate that having β1and β2set to their averages can still yield satisfactory results and it sheds light on the development of a universal fitness function for all solvents.

Table 1 Correlation coefficients (R) between E pred and E ZINDO when using V Coul,GS and V Coul,ES

Fig. 4 Correlation between E pred and E ZINDO using fitness functions based on V Coul,GS and V Coul,ES.Only the training sets are shown here, and the red lines denote the linear regression lines.

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

Fig. 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.(top) solvation in water (left) and acetonitrile (right) with β1 and β2 obtained from the regression method for a specific solvent; (bottom) solvation in water (left) and acetonitrile with universal β1 and β2 values.

3.3 Selection of rep resentative 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 Epred.Second, they were divided into Nconfigsections containing each an equal number of configurations. Nconfigis 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 Epredvalue closest to the mean Epredvalues 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 Epredvalues 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 Fn(x) and Gn(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 Fn(x)and Gn(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 EZINDOas 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 Dbiased/Drandis also shown in Fig. 6. The first thing to notice is that the Dbiased/Drandratios 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 Dbiased/Drandratios 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.

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.

Fig. 6 D and D* for two biased selection schemes based on predicted and pre-known excitation energies, and D biased /D rand.Uncertainties are only shown when they are larger than the symbol size.

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

Fig. 7 Average ZINDO excitation energies of selected configurations and uncertainties.top: water; bottom: acetonitrile. Dashed lines denote the averages for the full 1600 configurations.

3.4 Selection of solvent molecules

Fig. 8 Correlation between E pred and E TD-DFT for water and acetonitrile using Eqs. (6) and (8).

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 CH3group, 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 Epredwas correlated to ETD-DFTaccording to Eq. (6) with one revision of changing Eg,ZINDOto Eg,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,

Epred= EZINDO? Eg,ZINDO+ Eg,TD-DFT(8)

According to Table 3 and Fig. 8, Epredcalculated from Eq. (8)has a high correlation with ETD-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 Epredand EZINDOindicates 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.