• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2018-10-19 08:00:44XUEBaiCHENTiannanSIEPMANNIlja
    物理化學(xué)學(xué)報(bào) 2018年10期

    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.

    欧美黑人欧美精品刺激| 12—13女人毛片做爰片一| 色尼玛亚洲综合影院| 很黄的视频免费| 亚洲成av人片免费观看| 中亚洲国语对白在线视频| 国产黄色小视频在线观看| 综合色av麻豆| 国产精品九九99| 国产免费av片在线观看野外av| 少妇丰满av| 国产人伦9x9x在线观看| 国产精品永久免费网站| 99久久成人亚洲精品观看| 法律面前人人平等表现在哪些方面| 亚洲男人的天堂狠狠| 性欧美人与动物交配| 性欧美人与动物交配| 老司机午夜福利在线观看视频| 亚洲精品国产精品久久久不卡| 男女下面进入的视频免费午夜| 99re在线观看精品视频| 91av网一区二区| 国产又色又爽无遮挡免费看| 九九久久精品国产亚洲av麻豆 | 人人妻人人澡欧美一区二区| 欧美日本亚洲视频在线播放| 免费看a级黄色片| 美女大奶头视频| 搡老岳熟女国产| 精华霜和精华液先用哪个| 国产不卡一卡二| 神马国产精品三级电影在线观看| 国产高清视频在线播放一区| 夜夜看夜夜爽夜夜摸| 亚洲精品一卡2卡三卡4卡5卡| 国产精品亚洲一级av第二区| 国产精品98久久久久久宅男小说| 婷婷精品国产亚洲av| 欧美日韩精品网址| 白带黄色成豆腐渣| 欧美精品啪啪一区二区三区| 国产成+人综合+亚洲专区| 成人永久免费在线观看视频| 亚洲一区二区三区不卡视频| 国产伦在线观看视频一区| 窝窝影院91人妻| 国产亚洲精品综合一区在线观看| 国产精品野战在线观看| 窝窝影院91人妻| 啦啦啦观看免费观看视频高清| 91老司机精品| 一本精品99久久精品77| 日本一本二区三区精品| 黑人欧美特级aaaaaa片| 亚洲天堂国产精品一区在线| 黑人欧美特级aaaaaa片| 久久久精品大字幕| 窝窝影院91人妻| 国产激情偷乱视频一区二区| 脱女人内裤的视频| 99国产精品99久久久久| 可以在线观看毛片的网站| 夜夜躁狠狠躁天天躁| 久久久精品大字幕| 日韩欧美在线乱码| 久久精品亚洲精品国产色婷小说| 夜夜躁狠狠躁天天躁| 啦啦啦韩国在线观看视频| 免费一级毛片在线播放高清视频| 天堂av国产一区二区熟女人妻| 搞女人的毛片| 9191精品国产免费久久| 午夜日韩欧美国产| 99国产精品一区二区三区| 亚洲va日本ⅴa欧美va伊人久久| 免费无遮挡裸体视频| 三级男女做爰猛烈吃奶摸视频| av福利片在线观看| 亚洲精品一卡2卡三卡4卡5卡| 999久久久精品免费观看国产| ponron亚洲| 亚洲自拍偷在线| 国内揄拍国产精品人妻在线| 精品国产美女av久久久久小说| 欧美黑人欧美精品刺激| 久久久久久大精品| 欧美性猛交╳xxx乱大交人| 国产高清有码在线观看视频| 亚洲国产精品成人综合色| 嫩草影院入口| 欧美黑人巨大hd| 一边摸一边抽搐一进一小说| 丝袜人妻中文字幕| 后天国语完整版免费观看| 女警被强在线播放| 国产单亲对白刺激| www.精华液| 精品久久久久久久久久免费视频| 亚洲人成网站高清观看| 中文字幕av在线有码专区| 亚洲va日本ⅴa欧美va伊人久久| 亚洲自偷自拍图片 自拍| 曰老女人黄片| 国产亚洲精品久久久久久毛片| 精品一区二区三区视频在线观看免费| 黄色 视频免费看| 女人高潮潮喷娇喘18禁视频| 巨乳人妻的诱惑在线观看| 草草在线视频免费看| 黄片小视频在线播放| 波多野结衣高清作品| 久久精品综合一区二区三区| 男女之事视频高清在线观看| 亚洲精华国产精华精| 青草久久国产| 久久这里只有精品19| 99久国产av精品| 国产综合懂色| 日韩大尺度精品在线看网址| 精品国产乱子伦一区二区三区| 国产欧美日韩精品一区二区| 国产伦在线观看视频一区| 女同久久另类99精品国产91| 午夜视频精品福利| 18禁黄网站禁片免费观看直播| 午夜免费观看网址| 久久久久久久午夜电影| 91av网站免费观看| 麻豆av在线久日| 黄色片一级片一级黄色片| 欧美中文综合在线视频| 欧美最黄视频在线播放免费| 免费看光身美女| 日本黄色视频三级网站网址| 国产精品久久久av美女十八| 亚洲美女视频黄频| 黄片小视频在线播放| 美女高潮喷水抽搐中文字幕| 欧美日韩中文字幕国产精品一区二区三区| 日韩欧美免费精品| 中文字幕人妻丝袜一区二区| 国产成年人精品一区二区| 99精品欧美一区二区三区四区| 九九久久精品国产亚洲av麻豆 | 18禁黄网站禁片免费观看直播| 国产三级在线视频| 热99在线观看视频| 亚洲国产欧洲综合997久久,| 国产精品精品国产色婷婷| 日日干狠狠操夜夜爽| 久9热在线精品视频| 欧美3d第一页| 亚洲在线自拍视频| 一个人看视频在线观看www免费 | 欧美日韩精品网址| 日本黄大片高清| 天天躁狠狠躁夜夜躁狠狠躁| 国产精品一及| 俄罗斯特黄特色一大片| 性色avwww在线观看| 99热6这里只有精品| 黄色丝袜av网址大全| 国产伦精品一区二区三区四那| 757午夜福利合集在线观看| 国产午夜福利久久久久久| 亚洲精品一区av在线观看| 久久这里只有精品中国| 久久精品影院6| 亚洲五月天丁香| 久久久久性生活片| 国产高清三级在线| 欧美绝顶高潮抽搐喷水| 国内精品美女久久久久久| 无遮挡黄片免费观看| 草草在线视频免费看| 亚洲精品美女久久久久99蜜臀| 亚洲狠狠婷婷综合久久图片| 国产成年人精品一区二区| 亚洲av电影在线进入| 亚洲aⅴ乱码一区二区在线播放| 国模一区二区三区四区视频 | 免费看光身美女| 在线a可以看的网站| 麻豆久久精品国产亚洲av| 亚洲真实伦在线观看| 精品国产三级普通话版| 长腿黑丝高跟| 国产伦人伦偷精品视频| 男插女下体视频免费在线播放| 免费在线观看成人毛片| 欧美国产日韩亚洲一区| 麻豆一二三区av精品| 99国产精品一区二区三区| 三级男女做爰猛烈吃奶摸视频| 一进一出抽搐gif免费好疼| 精品人妻1区二区| 在线永久观看黄色视频| 久久精品综合一区二区三区| 亚洲av熟女| 狂野欧美白嫩少妇大欣赏| 美女扒开内裤让男人捅视频| 九色国产91popny在线| 国产伦精品一区二区三区四那| 757午夜福利合集在线观看| 国产亚洲精品久久久久久毛片| 国产精品综合久久久久久久免费| 欧美日本亚洲视频在线播放| 国产精品av视频在线免费观看| 国产精品九九99| 99久久无色码亚洲精品果冻| 三级毛片av免费| 国产激情久久老熟女| 中文字幕最新亚洲高清| 国产精品亚洲美女久久久| 999久久久精品免费观看国产| 国产男靠女视频免费网站| 国产蜜桃级精品一区二区三区| 国产成人福利小说| 90打野战视频偷拍视频| 99久久无色码亚洲精品果冻| 可以在线观看的亚洲视频| 中出人妻视频一区二区| 亚洲天堂国产精品一区在线| 亚洲欧美日韩东京热| 国产精品综合久久久久久久免费| avwww免费| 国产真实乱freesex| 亚洲成av人片免费观看| 午夜成年电影在线免费观看| 国产激情欧美一区二区| 国产激情欧美一区二区| 国产精品国产高清国产av| 嫩草影视91久久| 婷婷六月久久综合丁香| 黄色视频,在线免费观看| www.www免费av| 久久久久国产一级毛片高清牌| 国产精品电影一区二区三区| 网址你懂的国产日韩在线| 亚洲精品在线观看二区| 国产高清视频在线观看网站| 一个人观看的视频www高清免费观看 | 香蕉久久夜色| 欧美日韩黄片免| 中文字幕人成人乱码亚洲影| 成在线人永久免费视频| 美女cb高潮喷水在线观看 | 久久久久久人人人人人| 日韩欧美一区二区三区在线观看| 香蕉丝袜av| 男人的好看免费观看在线视频| 成熟少妇高潮喷水视频| 欧美乱码精品一区二区三区| 国产高清视频在线播放一区| 久久天堂一区二区三区四区| 两个人视频免费观看高清| 男人和女人高潮做爰伦理| 日本三级黄在线观看| 特级一级黄色大片| 午夜视频精品福利| 国产主播在线观看一区二区| 极品教师在线免费播放| 香蕉国产在线看| 欧洲精品卡2卡3卡4卡5卡区| 韩国av一区二区三区四区| 亚洲精品在线美女| 美女免费视频网站| 国产成人福利小说| 禁无遮挡网站| 成人三级做爰电影| 99久久99久久久精品蜜桃| 午夜视频精品福利| 亚洲,欧美精品.| 99国产精品一区二区蜜桃av| 99热6这里只有精品| 亚洲在线自拍视频| 麻豆成人av在线观看| 啦啦啦免费观看视频1| 黄色视频,在线免费观看| 午夜福利视频1000在线观看| 午夜成年电影在线免费观看| 又黄又粗又硬又大视频| ponron亚洲| 欧美黄色片欧美黄色片| 国产精品久久久久久人妻精品电影| 性欧美人与动物交配| 看免费av毛片| 欧美日韩黄片免| 亚洲精品美女久久久久99蜜臀| 少妇裸体淫交视频免费看高清| 中国美女看黄片| 欧美极品一区二区三区四区| 久久久久免费精品人妻一区二区| 999久久久国产精品视频| 99精品欧美一区二区三区四区| 搡老妇女老女人老熟妇| 欧美3d第一页| 一级毛片高清免费大全| 悠悠久久av| 亚洲av美国av| 中亚洲国语对白在线视频| 国产视频内射| 久久精品国产综合久久久| 欧美+亚洲+日韩+国产| 亚洲av电影在线进入| 非洲黑人性xxxx精品又粗又长| 精品国内亚洲2022精品成人| 精品不卡国产一区二区三区| 国产精品av视频在线免费观看| 最近最新中文字幕大全免费视频| 少妇的丰满在线观看| 在线国产一区二区在线| 男人舔女人的私密视频| 亚洲av免费在线观看| 这个男人来自地球电影免费观看| 色av中文字幕| 欧美绝顶高潮抽搐喷水| 欧美日本亚洲视频在线播放| 久久久国产精品麻豆| 午夜精品久久久久久毛片777| 欧美一级a爱片免费观看看| 精品国产超薄肉色丝袜足j| 中文字幕人成人乱码亚洲影| 俄罗斯特黄特色一大片| 少妇的丰满在线观看| 91九色精品人成在线观看| 国产乱人视频| 哪里可以看免费的av片| 青草久久国产| 男女那种视频在线观看| 日本成人三级电影网站| 午夜福利在线在线| 日韩欧美一区二区三区在线观看| 69av精品久久久久久| 蜜桃久久精品国产亚洲av| 欧美乱色亚洲激情| 欧美日本亚洲视频在线播放| 免费观看人在逋| 国产97色在线日韩免费| 老司机深夜福利视频在线观看| 国产精品久久久av美女十八| 91在线精品国自产拍蜜月 | 国产真实乱freesex| 国产精品电影一区二区三区| 国内少妇人妻偷人精品xxx网站 | 后天国语完整版免费观看| 久久伊人香网站| 亚洲精品在线美女| 久久精品综合一区二区三区| 欧美中文日本在线观看视频| 9191精品国产免费久久| 日韩欧美国产一区二区入口| 免费观看人在逋| 精品熟女少妇八av免费久了| 一区二区三区国产精品乱码| 久久久国产欧美日韩av| 欧美在线一区亚洲| 国产成人啪精品午夜网站| 在线a可以看的网站| 亚洲熟妇中文字幕五十中出| 日韩有码中文字幕| 最近最新中文字幕大全电影3| 国产一区二区三区视频了| 午夜福利免费观看在线| 99热6这里只有精品| 亚洲最大成人中文| 国产伦精品一区二区三区视频9 | avwww免费| 一边摸一边抽搐一进一小说| 亚洲 欧美一区二区三区| 国产男靠女视频免费网站| 岛国在线免费视频观看| 国产精品 欧美亚洲| 动漫黄色视频在线观看| 一二三四社区在线视频社区8| 很黄的视频免费| 国产熟女xx| 国产真人三级小视频在线观看| 欧美zozozo另类| 国产精品 欧美亚洲| 久久精品91蜜桃| 亚洲精品国产精品久久久不卡| 亚洲人成伊人成综合网2020| 欧美乱码精品一区二区三区| 97人妻精品一区二区三区麻豆| 欧美xxxx黑人xx丫x性爽| 人妻丰满熟妇av一区二区三区| 国产男靠女视频免费网站| 色播亚洲综合网| 淫妇啪啪啪对白视频| 国产三级在线视频| 国产精品 国内视频| bbb黄色大片| 欧美日本视频| 久久久国产欧美日韩av| 欧美zozozo另类| 国产精品永久免费网站| 搞女人的毛片| 丰满人妻一区二区三区视频av | 国产高潮美女av| 亚洲性夜色夜夜综合| 亚洲第一电影网av| 亚洲专区中文字幕在线| 最近视频中文字幕2019在线8| 日本在线视频免费播放| 午夜福利高清视频| 日本a在线网址| 日韩欧美国产在线观看| 久久精品人妻少妇| 亚洲国产精品成人综合色| 国产激情欧美一区二区| 午夜福利欧美成人| 久久九九热精品免费| 搡老熟女国产l中国老女人| 国产私拍福利视频在线观看| av片东京热男人的天堂| 中文亚洲av片在线观看爽| 嫩草影院精品99| 美女高潮的动态| 无限看片的www在线观看| 精品国产美女av久久久久小说| 亚洲欧美激情综合另类| 激情在线观看视频在线高清| 国产一区二区在线av高清观看| 亚洲成a人片在线一区二区| 国产精品1区2区在线观看.| 欧美一区二区国产精品久久精品| 久久久久国产一级毛片高清牌| 黄色 视频免费看| 99久久无色码亚洲精品果冻| 丁香六月欧美| 婷婷六月久久综合丁香| 三级毛片av免费| 欧美激情久久久久久爽电影| 亚洲中文av在线| 亚洲精品乱码久久久v下载方式 | 三级国产精品欧美在线观看 | 欧美日本视频| 国产私拍福利视频在线观看| 欧美日韩中文字幕国产精品一区二区三区| 亚洲色图av天堂| 亚洲美女黄片视频| 亚洲性夜色夜夜综合| 熟女人妻精品中文字幕| 日韩欧美免费精品| 久久久久国内视频| 午夜福利成人在线免费观看| 99久久99久久久精品蜜桃| 日韩精品青青久久久久久| 女生性感内裤真人,穿戴方法视频| 在线观看免费午夜福利视频| 又爽又黄无遮挡网站| 色哟哟哟哟哟哟| 国产欧美日韩精品亚洲av| 青草久久国产| 精品一区二区三区av网在线观看| 三级国产精品欧美在线观看 | 在线观看美女被高潮喷水网站 | 久久久水蜜桃国产精品网| 国模一区二区三区四区视频 | 天堂网av新在线| 亚洲天堂国产精品一区在线| 亚洲无线观看免费| 久久久久久九九精品二区国产| 99在线人妻在线中文字幕| 日韩av在线大香蕉| 岛国在线免费视频观看| 俺也久久电影网| 久久中文字幕人妻熟女| 亚洲精品美女久久久久99蜜臀| 国产真人三级小视频在线观看| 久久天堂一区二区三区四区| 日韩免费av在线播放| 国产乱人视频| 黄片小视频在线播放| 久9热在线精品视频| 脱女人内裤的视频| 中文亚洲av片在线观看爽| 热99在线观看视频| 免费观看的影片在线观看| 国产97色在线日韩免费| 欧美日韩国产亚洲二区| 日本黄色视频三级网站网址| 在线观看一区二区三区| 色播亚洲综合网| 天堂影院成人在线观看| 亚洲成人免费电影在线观看| 97超视频在线观看视频| 亚洲欧美一区二区三区黑人| 人妻丰满熟妇av一区二区三区| 999精品在线视频| 国产av不卡久久| 午夜精品久久久久久毛片777| 成人一区二区视频在线观看| 免费看a级黄色片| 五月玫瑰六月丁香| 婷婷丁香在线五月| 岛国在线观看网站| 草草在线视频免费看| 国产亚洲精品久久久com| 黑人操中国人逼视频| 国产91精品成人一区二区三区| 久久午夜亚洲精品久久| 小蜜桃在线观看免费完整版高清| 动漫黄色视频在线观看| 欧美性猛交黑人性爽| 欧美黄色淫秽网站| 久久久久久人人人人人| 亚洲性夜色夜夜综合| 成人亚洲精品av一区二区| 精品不卡国产一区二区三区| 久久中文看片网| 一卡2卡三卡四卡精品乱码亚洲| 小说图片视频综合网站| 长腿黑丝高跟| 亚洲欧美激情综合另类| 国产乱人伦免费视频| av在线天堂中文字幕| 成熟少妇高潮喷水视频| 在线观看66精品国产| 久久久成人免费电影| 丝袜人妻中文字幕| 天堂网av新在线| 亚洲 国产 在线| 国产亚洲欧美在线一区二区| 嫁个100分男人电影在线观看| 午夜精品一区二区三区免费看| 又黄又爽又免费观看的视频| 国产亚洲欧美在线一区二区| 99久久综合精品五月天人人| 精品不卡国产一区二区三区| 国产欧美日韩精品亚洲av| 在线观看午夜福利视频| 丰满人妻熟妇乱又伦精品不卡| 欧美日韩黄片免| 90打野战视频偷拍视频| 特大巨黑吊av在线直播| 国产99白浆流出| 国产一区二区三区在线臀色熟女| 老司机午夜十八禁免费视频| 又黄又粗又硬又大视频| 女同久久另类99精品国产91| 色播亚洲综合网| 国产精品野战在线观看| 日韩免费av在线播放| 日本a在线网址| 国产视频内射| 老熟妇仑乱视频hdxx| 制服丝袜大香蕉在线| 欧美另类亚洲清纯唯美| 婷婷亚洲欧美| 99在线人妻在线中文字幕| svipshipincom国产片| 欧美一区二区精品小视频在线| 叶爱在线成人免费视频播放| 国产久久久一区二区三区| 曰老女人黄片| 成人午夜高清在线视频| 国产精品一区二区免费欧美| 日本免费一区二区三区高清不卡| 在线观看美女被高潮喷水网站 | 男女做爰动态图高潮gif福利片| 免费看日本二区| 一区二区三区高清视频在线| 我要搜黄色片| 法律面前人人平等表现在哪些方面| 麻豆成人av在线观看| 精品乱码久久久久久99久播| 亚洲成a人片在线一区二区| 又紧又爽又黄一区二区| 一个人免费在线观看的高清视频| 国产 一区 欧美 日韩| 日本免费一区二区三区高清不卡| 18禁美女被吸乳视频| 中文字幕最新亚洲高清| 久久久国产欧美日韩av| 中亚洲国语对白在线视频| 亚洲av五月六月丁香网| 国产亚洲欧美在线一区二区| 欧洲精品卡2卡3卡4卡5卡区| 黄色 视频免费看| 久久久久久久精品吃奶| 国产aⅴ精品一区二区三区波| 日本与韩国留学比较| 亚洲国产欧美一区二区综合| 精品久久久久久久人妻蜜臀av| 国产精品香港三级国产av潘金莲| 后天国语完整版免费观看| 国产成人一区二区三区免费视频网站| 亚洲午夜理论影院| 久久久成人免费电影| 亚洲人与动物交配视频| 久久久久久九九精品二区国产| 夜夜爽天天搞| 一区二区三区国产精品乱码| 香蕉av资源在线| 最近最新中文字幕大全电影3| av中文乱码字幕在线| 午夜精品久久久久久毛片777| 我的老师免费观看完整版| 欧美黑人巨大hd| 国产伦精品一区二区三区四那| 在线永久观看黄色视频| 丝袜人妻中文字幕| 免费无遮挡裸体视频| 欧美一区二区国产精品久久精品| 久久午夜综合久久蜜桃| 久久精品aⅴ一区二区三区四区| 精品无人区乱码1区二区| 变态另类成人亚洲欧美熟女| 色综合亚洲欧美另类图片| 国产又黄又爽又无遮挡在线| 欧美一区二区精品小视频在线| 精品久久久久久久末码|