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

    Microscopic Investigation of Ethylene Carbonate Interface:A Molecular Dynamics and Vibrational Spectroscopic Study

    2018-10-19 08:00:46WANGLinXINLiangISHIYAMATatsuyaPENGQilingYEShenMORITAAkihiro
    物理化學學報 2018年10期

    WANG Lin ,XIN Liang ,ISHIYAMA Tatsuya ,PENG Qiling ,YE Shen ,MORITA Akihiro ,?

    1 Department of Chemistry,Graduate School of Science,Tohoku University,Aoba-ku,Sendai980-8578,Japan.

    2 Elements Strategy Initiative for Catalystsand Batteries(ESICB),Kyoto University,Kyoto 615-8520,Japan.

    3 Department of Applied Chemistry,Graduate School of Scienceand Engineering,University of Toyama,Toyama 930-8555,Japan.4 Institute for Catalysis,Hokkaido University,Kita-ku,Sapporo 001-0021,Japan.

    Abstract:Ethylene carbonate(EC)liquid and its vapor-liquid interface were investigated using a combination of molecular dynamics(MD)simulation and vibrational IR,Raman and sum frequency generation(SFG)spectroscopies.The MD simulation was performed with a fl exible and polarizable model of the EC molecule newly developed for the computation of vibrational spectra.The internal vibration of the model was described on the basis of the harmonic couplings of vibrational modes,including the anharmonicity and Fermi resonance coupling of C=O stretching.The polarizable model was represented by the charge response kernel(CRK),which is based on ab initio molecular orbital calculations and can be readily applied to other systems.The fl exible and polarizable model can also accurately reproduce the structural and thermodynamic properties of EC liquid.Meanwhile,a comprehensive set of vibrational spectra of EC liquid,including the IR and Raman spectra of the bulk liquid as well as the SFG spectra of the liquid interface,were experimentally measured and reported.The set of experimental vibrational spectra provided valuable information for validating the model,and the MD simulation using the model comprehensively elucidates the observed vibrational IR,Raman,and SFG spectra of EC liquid.Further MD analysis of the interface region revealed that EC molecules tend to orientate themselves with the C=O bond parallel to the interface.The MD simulation explains the positive Im[χ(2)](ssp)band of the C=O stretching region in the SFG spectrum in terms of the preferential orientation of EC molecules at the interface.This work also elucidates the distinct lineshapes of the C=O stretching band in the IR,Raman,and SFG spectra.The lineshapes of the C=O band are split by the Fermi resonance of the C=O fundamental and the overtone of skeletal stretching.The Fermiresonance of C=O stretching was fully analyzed using the empirical potential parameter shift analysis(EPSA)method.The apparently different lineshapes of the C=O stretching band in the IR,Raman,and SFG spectra were attributed to the frequency shift of the C=O fundamental in different solvation environments in the bulk liquid and at the interface.This work proposes a systematic procedure for investigating the interface structure and SFG spectra,including general modeling procedure based on ab initio calculations,validation of the model using available experimental data,and simultaneous analysis of molecular orientation and SFG spectra through MD trajectories.The proposed procedure provides microscopic information on the EC interface in this study,and can be further applied to investigate other interface systems,such as liquid-liquid and solid-liquid interfaces.

    Key Words:Ethylene carbonate;SFG;Fermiresonance;EPSA

    1 Introduction

    Organic carbonates,such as ethylene carbonate(EC),propylenecarbonate(PC),dimethyl carbonate(DMC)and their mixtures,have been widely used as solvents in modern batteriesdueto anumber of their suitableproperties,including high dielectric constants,wide electrochemical window and stability1.It is known that the electrode/electrolyte solution interfaces are quite important in the performance of batteries,as the electrochemical reactions together with ion interactions occur at the interface2–5.Therefore,many studies have investigated the electrode/electrolyte solution interfaces by various spectroscopic techniques,such as X-ray6,7,infrared(IR)8,9,Raman10,11,and sum frequency generation(SFG)spectroscopies12–16.

    Among thosespectroscopic techniques,SFGisaparticularly powerful tool to detect the molecules at interfaces due to its high interface sensitivity and selectivity17–19.However,it is often diきcultto precisely interpretthe SFGspectra,and reliable theoretical analysis is strongly required to extract microscopic information from observed spectra.In recent years,theoretical calculationmethodsof SFGspectraby molecular dynamic(MD)simulation have been developed19–30and successfully applied to investigateaqueousinterfaces19–22,25–27,29,30aswell asorganic interfaces31–34.This work presents our analysis of ECliquid in combination of MDsimulation and experimental measurements of vibrational spectra.We report experimental IR,Raman and SFG spectra of EC liquid,and these experimental spectra are analyzed by the MD simulation.We are particularly interested in the interface and SFG spectra,while the experimental IR and Raman spectra are utilized to examine the reliability of the molecular model and MDsimulation.

    In the MD calculation of SFG spectra,molecular modeling is of key importance.In this work,we propose a fl exible and polarizable model for the EC molecule.The fl exible model includes the anharmonicity and Fermi resonance coupling of C=O stretching to properly reproduce its vibrational spectra,while the polarizable model is based on the Charge Response Kernel(CRK)theory35,36.The CRK theory provides non-empirical description of intramolecular polarization on the basis of site representation.It is fully based on ab initio molecular orbital or the density functional theory(DFT),and thus it is easily applicable to general molecules.A remarkable advantage of the CRK theory in the SFG analysis is that the CRK model is designed to properly describe instantaneous dipole moment and polarizability of molecules during MD simulation,which are requisite to the SFG analysis based on the time correlation function formula21,24.So far we have applied the CRK model to various molecules,including water37,sulfate ion38,methanol39,benzene32,ethanol40,PC and DMC33.Besides calculating vibrational spectra,CRK model was also used to investigate the diあusion of radicals in solvent35,41.

    The vibrational modeling of EC molecule in the present paper is performed after that of PC and DMC in our previous work33.However,the modeling of EC is more complicated than PC or DMC,as the C=O stretching band of EC is obviously split by the Fermi resonance.Therefore,correct description of the Fermi resonance is necessary to the vibrational modeling of EC.The mechanism of Fermi resonance is further analyzed with the empirical potential parameter shift analysis(EPSA)method31,in which one systematically modifi es some parameters of intramolecular force fi eld and observes the response of calculated vibrational spectra.This analysis is useful for interpreting complicated vibrational spectra,especially when the Fermi resonance is involved.

    The remainder of this paper is constructed as follows.Section 2 presents the details of MD simulation,including the force field and MD conditions.Section 3 presents the experimental conditions of spectroscopic measurements.Then we describe the analysis methods of vibrational spectra in Section 4.Section 5 is devoted to the results and discussion about the IR,Raman,SFG spectra and surface orientation analysis of EC.Brief summary and conclusion follow in Section 6.

    2 Molecular dynamics simulation

    2.1 Force fi eld

    The total potential of ECliquid,Utotal,is represented by

    where Uintrais the intramolecular potential,UCis the electrostatic interaction and ULJis Lennard-Jones(LJ)interaction.The function forms are explained in detail in the following.

    2.1.1 Intramolecular potential

    The intramolecular potential Uintrain Eq.(1)is represented asthesum of thoseof constituent molecules,

    Table1 Natural internal coordinates S1–S24 of EC.

    Fig.1 Molecular structureof EC.Thered vectorsillustratethedefi nition of molecule-fi xed coordinates(ξ,η,ζ).

    where uintra(i)is the intramolecular potential of molecule i.uintraof each molecule is represented with a set of the natural internal coordinates42(NICs)to describeall theintramolecular degrees of freedom.The structure of ECmolecule is illustrated in Fig.1.A EC molecule consisting of Natom=10 atoms has 3Natom? 6=24 NICs,and the defi nition of NICs,S1–S24,are listed in Table 1.The functional form of uintrais given as follows:In theright sideof Eq.(3),thefi rst term(I)representsharmonic potential,whereΔSi=Si?denotes the displacement of Sifrom its equilibrium value,and ki,jis a symmetric matrix of the harmonic force constants.The second term(II)describestheanharmonic potential along C=Ostretching(S1),with k3and k4being higher-order force constants.The third term(III)representsthe Fermiresonance(FR)couplingbetween the C=Ostretching(S1)and theovertoneof skeletal stretching(S4).Thesecond(II)and third(III)termsareenergetically minor yet important to represent the vibrational spectra of EC in the C=O stretching region.

    The parameters in Eq.(3)are determined with quantum mechanical(QM)calculations in the following procedure.are obtained from geometry optimization of EC molecule and ki,jare determined from the Hessian at the optimized geometry.The optimized geometry of EC molecule in the molecule-fi xed coordinates(ξ,η,ζ,see Fig.1)is shown in Supporting Information(SI).k3,k4andκ144are calculated from corresponding higher-order derivatives of energy with respect to the internal coordinates.The parameters related to the C=O stretching,k11,k44,k3,k4andκ144,are further tuned so as to be consistent with the experimental IR,Raman and SFG spectra in the C=O stretching region.The QM calculations are performed at B3LYP43,44/aug-cc-pVTZ45level using Gaussian 09 package46.The parameters thus obtained aresummarized in SI.

    In the above modeling,we included an empirical fitting procedure to tune several parameters(k11,k44,k3,k4and κ144)according to experimental spectra.In previous PC and DMC study33,the modeling procedure did not include any empirical fi ttings.And the obtained model can successfully reproduce experimental vibrational spectra.In this work,however,the C=O stretching of EC molecular is clearly split by Fermiresonance,which is diあerent from PCand DMC.The splitting of Fermiresonance isquitesensitive to thevibrational frequency of C=O stretching and corresponding overtone.In order to obtain reasonable Fermi splitting in C=O stretching,we slightly tuned the parameters that are related to the vibrational frequency of S1and overtone of S4.It should be noted that this treatment will not infl uence the objective analysisof interfacestructureand vibrational spectra.

    2.1.2 Electrostatic interaction

    The electrostatic interaction UCin Eq.(1)isgiven by

    where Qaiand Vaiare respectively the partial charge and electrostatic potential at site a of molecule i.Kabiis the CRK of molecule i,which is defi ned as Kabi=(?Qai/?Vbi).The fi rst term of the right-hand side of Eq.(4)indicates the intermolecular site-site Coulombic interactions,and the second term represents the electronic reorganization energy.The electrostatic potential Vaiiswritten by

    where rai,bjis the distance between site a of molecule i and site b of molecule j.fai,bjisthedamping factor of short-range Coulomb interaction39.In thepresent ECliquid fai,bj=1(no damping)is employed,which leads to stable MD trajectories.The partial charge Qaichanges with the electrostatic potential by

    where we omit the suきx of molecule i.in Eqs.(7)and(8)are the partial charge and CRK at the equilibrium conformation{ΔSt=0}.The present model can express the fl uctuating charge distribution in response to the vibration and electrostatic potential.

    Table2 ESPchargesand LJ parameters,σ (?)and ε (J·mol?1),of EC.

    All the ingredient parameters in Eqs.(7,8),Qeqa,Keqab,(?Qa/?St) and (?Kab/?St), are calculated with B3LYP/aug-cc-pVTZ using Gaussian 09 package46with our extension to incorporate the CRK calculation.The partial charges are determined via the electrostatic potential(ESP)fi tting using ChelpG method47,with employing the damping treatment48to avoid ill-defi ned behavior of ESP charges at buried sites.The damping parameter is set to beλ=0.25.The fi rst order derivatives(?Qa/?St)and (?Kab/?St)are calculated by 4-point numerical diあerentiation method.The calculated ESPcharges are listed in Table 2.Other parameters obtained assuch aregiven in SI.

    The above procedure of fl exible and polarizable modeling can be readily transferred to other organic molecules.The parametersof target moleculescan beobtained through aset of QM calculations.It isalso known that QM calculation for large molecules is quite time consuming.Currently we are working on thetransferable CRK model for largeorganic molecules.

    2.1.3 LJ interaction

    The LJinteraction ULJin Eq.(1)isrepresented by

    where rai,bjis the distance between the site a of molecule i and thesite b of molecule j.The LJparametersof ECaretaken from the OPLS force fi eld49.The LJ parameters for distinct atom species are giv■e by the Lorentz-Berthelot rule:σab=(σa+ σb)/2,εab=.All the LJparameters are given in Table 2.

    2.2 MD conditions

    The present MD simulation examines the properties of bulk ECliquid and interface,and accordingly wecarried out two sets of simulationtoinvestigatebulk liquid and interfacerespectively.The MD simulation of bulk liquid deals with IR and Raman spectra besides some molecular properties in the liquid,while the vapor-liquid simulation treats SFG spectra and molecular orientation at theinterface.

    In the simulation of bulk liquid,325 EC molecules were initially placed in a cubic cell of 33?×33?×33?(1?=0.1 nm).Then MD simulation under isothermal-isobaric(NPT)ensemble was carried out for 0.3 ns of equilibration and 1.8 ns of production.In vapor-liquid interface,325 EC molecules were initially placed in the same cubic region of a rectangular cell of Lx×Ly×Lz=33?×33?×165?.Therefore,the liquid forms a slab with two vapor-liquid interfaces normal to the z axis.For the interface simulation,48 independent initial confi gurations were prepared and treated in parallel.MD simulations under canonical(NVT)ensemble were carried out for 0.3 ns of equilibration and for 0.6 ns of production for each confi guration.Therefore,the total production time amounts to 0.6 ns×48=28.8 ns for the interfacesimulation.

    In either simulation,three-dimensional periodic boundary conditions were employed.The initial confi gurations were generated using the Packmol package50.The temperature was 313.15 K(40°C)using the Nosè-Hoover thermostat51,52with a coupling constant of 0.4 ps.The temperature was set to be slightly above the melting point 39°C53.The pressure in the bulk simulation was set to 1 bar using the Hoover barostat54with a coupling constant of 2 ps.The equations of motion were integrated using the velocity Verlet algorithm with a time step of 0.61 fs.The Ewald summation method55,56with a separation parameter of 0.242 ??1was used for long range electrostatic interactions.The LJ interactions were evaluated using 14? cutoあ,with adopting tail corrections in energy and pressure calculation55.

    3 Experimental conditions

    3.1 Materials

    Special grade regent of ethylene carbonate(EC)was purchased from Kanto Chemical Co.,Inc.(Tokyo,Japan).The chemical wasused without further purifi cation.

    3.2 Infrared measurements

    A Fourier transform infrared(FTIR)spectrometer,JASCO 6100 with a DTGS detector was used in the study.The IR spectrum for each sample was obtained using an attenuated total refl ection(ATR)accessary with a ZnSe crystal(JASCO).IR characterizations were carried out at 40°C to melt EC into liquid phase.Typically,8 interferograms were accumulated for one spectrum with a resolution of 4 cm?1.All IR spectra shown in the study were after the ATRcorrection.

    3.3 Raman measurements

    The Raman spectra were recorded using a Renishaw inVia microscopewithatemperature-controllablethin-layer liquid cell at 40°C to melt EC into liquid phase.The excitation light(Nd:YAG laser at 532 nm,5%of intensity)was focused on the liquid surface through a 20×objective and the acquisition time was extended to 15 s.The spectral resolution of the Raman spectra in the study was ca.1.0 cm?1.Details were described elsewhere57.

    3.4 SFG measurements

    The details about our femtosecond broadband SFG system were described elsewhere33,58–61.Briefl y,a Ti:sapphire oscillator(MaiTai,Spectra Physics)provides an 800 nm seed beam with pulse duration of ca.100 fs.The seed beam is amplifi ed by a 1 kHz,2 W Ti:sapphire regenerative amplifi er(Spitfi re Pro,Spectra Physics)pumped by a Nd:YLF laser(EMPower,Spectra Physics)to produces 120 fs pulses at 800 nm.Half of the amplifier’s output is directed into an optical parametric amplifi er(OPA)followed by diあerence frequency generation(DFG)to generate mid-IR pulses(TOPAS,Light Conversion Inc.).The IRpulses are tunable from 2.5 to 10μm with a spectral width of~200 cm?1.The other half of the amplifi er’s output is directed into a homemade pulse shaper to generate picosecond pulse at 800 nm (an FWHM of~10 cm?1)as the visible beam.The SFG signal was recorded by a CCD detector(DU420-BV,Andor Technology)attached to a spectrograph(MS3504,Solar-TII).A specially designed cell with a heating mechanism was employed for the SFG measurement on the EC/air interface at 40°C.The cell was covered by a thin CaF2window plate.The SFG spectra were recorded by spatially and temporally overlapping the IR and visible pulses on the surface of acetonitrile/water liquid under the ssp(i.e.,s-SFG,s-visible,p-IR)and sps polarization combinations.The angles of incidence of IRand visible beams were 50°and 70°,respectively.All the SFG spectra of the liquid surface were normalized by a SFG spectrum of a GaAs.It took 20 min to integrate one SFG spectrum on the liquid/air interface.

    4 Analysis of vibrational spectra

    The detailed analysis method of vibrational IR,Raman and SFG spectra calculated in MD simulation can be found in our previous work31,33,34.Here we briefl y outline the analysis method.

    4.1 IR and Raman spectra

    The vibrational IR,isotropic Raman and anisotropic Raman spectraarecalculated fromthefollowinglineshapefunctions62:

    where〈 〉denotes the ensemble average.M is the dipole moment vector of the system,ˉA=(1/3)Tr A is the isotropic part of polarizability tensor A,and B=A?ˉA I is the traceless anisotropic part of the polarizability tensor.Raman spectrum is given by the istropic and anisotropic line shape functions62,

    where n is the refractive index,c is the speed of light,kBis the Boltzmann constant,and T is the temperature.ω0is the frequency of incidentradiation.Inthefollowingcalculations,we determine it to be consistent with the experimental wavelength of 532 nm.

    4.2 χ(2)and SFG spectra

    Theintensitiesof SFGspectraaregiven by the squareof the second-order susceptibilityχ(2)19:

    where the ssp and sps denote the polarization combinations of SFG,visible and IR lights,respectively.(For example,ssp meansthecombination of s-polarized SFG,s-polarized visible,and p-polarized IR.)The diあerent polarization combinations correspond to the diあerent tensor elements of,where the suきxes p,q,r denote the space fi xed coordinates x–z19.Note that Eq.(14)omits the prefactor of|χ(2)|2,but it is useful to analyze the line shape(frequency dependence)of the SFG spectra.

    Since Eq.(16)includes the time correlation function of A and M,itisdecomposed into theself and crosscorrelationparts of constituent molecules33,

    whereαpq,jandμr,jare the polarizability and dipole moment of molecule j under the space-fi xed coordinates,respectively.Eq.(18)representsas the sum of vibrational resonant component of molecular hyperpolarizability(j),and accordingly allow us to decomposeinto the contributions of molecules.

    Each molecular hypderpolarizability tensorin Eq.(18)isexpressed for the C=Ostretching band by24,33

    where the suきx j is omitted.qCOis the normal mode coordinate of the C=O stretching.mCOis the reduced mass,ωCOis the harmonic frequency,andΓCOis the damping factor for the mode.(?αpq/?qCO)and(?μr/?qCO)correspond to the transition dipole moment and the transition polarizability for the C=O stretching mode in the space-fi xed coordinates,respectively,and they are related to the transition properties in themolecule-fi xed coordinateswith the rotation matrix D by63

    where p,q,r denote the suきxes of space-fi xed coordinates(x,y,z),and p′,q′,r′denote those of molecule-fi xed coordinates(ξ,η,ζ).The superscript “mol” in Eq.(21)indicates tensor elements defi ned in the molecule-fi xed coordinates.

    Eq.(20)indicates that the sign of Im[]is governed by the sign of(?αpq/?qCO)(?μr/?qCO)in the space-fi xed coordinates,which changes with molecular orientation.By analyzing the relation between the factor(?αpq/?qCO)(?μr/?qCO)and molecular orientation at interface,we can interpret the imaginary part of SFG spectra and the molecular orientation at interface.

    4.3 Empirical parameter shift analysis(EPSA)

    The Fermi resonance in the C=O stretching band of EC is analyzed withtheEPSA procedure31,34.Asshownby Eq.(3),the Fermi resonance is described by the term(III),κ144ΔS1Δwhich considers the anharmonic coupling between the C=O stretchingmode(S1)and theovertoneof skeletalstretching(S4).Therefore,we gradually change the harmonic frequency of S4and trace the response of the band shapes of the IR,Raman and SFG spectra.The modifi cation of S4frequency is done by multiplying the harmonic force constant k4,4by a scaling factor(1.0–1.12).A scaling factor of 1.0 retains the original spectra,whileafactor of 1.12 shiftstheovertoneof S4by about 1800×?1.0)≈105cm?1.Thisscalinghassignifi cant infl uence on the mixing mechanism of the Fermi resonance in the C=O stretching,and the analysis of the spectral response helps understanding the mechanism of Fermi splitting in the C=Ostretchingof ECinbulk liquid and atinterface.Theresults of EPSA arediscussed in Section 5.5.

    Table3 Calculated and experimental resultsof density(ρ),heat of vaporization(ΔH v),and surfacetension(γ)of ECliquid.

    5 Results and discussion

    Here we discuss the computational and experimental results of liquid EC and its interface.Section 5.1 presents the computational results of some thermodynamic properties of liquid EC to validate the accuracy of the MD simulation.Section 5.2 compares the experimental and simulated IR and Raman spectra for the bulk liquid EC.Then we treat the interfaceof ECliquid and the SFGspectrain Secs.5.3 and 5.4.The Fermi resonance mechanism in the C=O band of EC in the vibrational spectra is elucidated in Section 5.5.

    5.1 Thermodynamic properties

    In order to validate the performance of our polarizable force fi eld,we fi rst calculated thermodynamic properties of EC liquid.The calculated results,together with available experimental data64–66,are summarized in Table 3.

    The heat of vaporizationΔHvwas calculated from the bulk simulation by thefollowing equation:49

    where UCand ULJare defi ned in Eqs.(4)and(9),respectively,and N is the number of molecules.The surface tensionγwas calculated from the interface MD simulation by the diあerence of normal and tangential pressure tensors67:

    where Pxx,Pyy,and Pzzarethe xx,yy,and zz componentsof the pressure tensor,respectively.

    As shown in Table 3,the calculated density of ECis in good agreement with experimental value.However,the calculated results of heat of vaporization and surface tension are slightly overestimated,which suggests the intermolecular interactions between EC molecules tend to be overestimated in our MD simulation.Such behaviors were also found in the previous PC and DMC simulation33.The overestimation of three organic carbonates is arguably due to the LJ parameters we used.We have employed the general OPLSparameters,which have been optimized to non-polarizable force fi eld,and thus they may tend to overestimate intermolecular interactions in combination of our polarizable model.We nevertheless adopted the OPLS parameters to facilitate transferability of the present model,though there may be some room for optimization.The current model is accurate enough to analyze the experimental vibrational spectra and the structure of interface,as we discuss below.

    Fig.2 Calculated(a,c)and experimental(b,d)IR(a,b)and Raman(c,d)spectra of liquid EC.

    5.2 Infrared and Raman spectra

    Fig.2 displays the calculated and experimental IR and Raman spectra of EC liquid.The overall features of IR and Raman spectra are reproduced by the present calculation results.The C=O stretching band(1700–1900 cm?1)shows stronger intensity than C–H stretching band(2800–3100 cm?1)in the IR spectrum,while vice versa in the Raman.These behaviors are well represented by the calculated results.It should be noted that the experimental C–H stretching band is somewhat red shifted than the calculated one.This deviation is because the current model did not incorporate the anharmonic couplings in the C–H stretching modes.It is also found that larger discrepancies occur in the lower frequency region(1000–1500 cm?1),especially in the Raman spectra.These discrepanciesarecaused by thefact that thecurrent model does not consider conformer dependence of vibrational modes.In previous study of ethanol40,we have confi rmed that the neglect of conformer dependence causes the overestimation of transition polarizability and Raman intensity in the lower frequency region.The analysis of the C–H stretching band and lower frequency bands are not the present purpose,and these departures do not infl uence the following analysis of C=O stretching and the related molecular orientation at the interface.

    Then we focus on the C=O stretching band,and compare the calculated and experimental IRand Raman band in detail in Fig.3.The experimental IR and Raman spectra(panels(b,d))show two peaks in the C=O band,and the intensity of the lower frequency peak is slightly stronger than the higher frequency one.The calculated spectra (panels (a,c))qualitatively reproduce the two-band structure in C=O stretching region.The peak positions and relative intensity of the two peaks are properly reproduced,which further validate the present model in representing vibrational spectra of EC liquid.The two-band structure of the C=O stretching band is due to the Fermiresonance,as will be explained in Section 5.5.

    5.3 Molecular orientation at interface

    Fig.3 Calculated(a,c)and experimental(b,d)IR(a,b)and Raman(c,d)spectra of liquid EC in the C=O stretching region.Thepeak frequencies arelabeled in the panels.

    Fig.4 (a)Calculated density profileof EC interfaceasa function of?z.Bluelinedenotesthefitting result by Eq.(24).(b)Definition of theorientation angleθC O of an EC molecule.(c)Product ofdensityρand cosC O by Eq.(25)asa function of?z.(d)Probability density distributions P(cosC O,{z 1<?z<z 2})in thebulk and interfaceby Eq.(26).

    Before analyzing the SFG spectra of the EC interface,we analyze the interfacial structure through MD simulation trajectories.Fig.4a shows the calculated density profi le at the vapor-liquid interface of liquid EC.Here we define the depth coordinate?z,where?z=0 denotes the position of Gibbs dividing surface,and?z>0 refersto thevapor side while?z<0 to the liquid side.The?z coordinates of molecules are defi ned with their center of mass positions.Then density profile along the depth coordinate?z ispresented by the following formula68:

    whereρ0indicatesthebulk density of theliquid,and z10-90isthe 10-90 thicknessof thesurface.Theoptimized parametersin Eq.(24)by the least square fi tting are ρ0=1.330 g·cm?3,z10-90=3.41?for the EC interface.The fi tting curve is also shown in Fig.4a.

    Next we analyzed the molecular orientation of EC at interface.The orientation angle of an EC molecule is defi ned by thetilt angleofθCObetween the C=Obond and thesurface normal z,as illustrated in Fig.4b.Then the ensemble average ofρ〈cosθCO〉as a function of?z was calculated through

    whereρ(cosθCO,?z)is the density distribution as a function of cosθCOand?z.Theresult ofρ〈cosθCO〉(?z)isshown in Fig.4c.Itisclear thatρ〈cosθCO〉(?z)convergestozeroin?z>5?,asthe density is drastically reduced in the vapor region.On the other hand,ρ〈cosθCO〉(?z)→0 is also found in the suきciently deep liquid region(?z<?10?),because of the random orientation in the bulk liquid.In the interface region(?10?<?z<5?),ρ〈cosθCO〉shows an oscillating behavior.At the topmost layer(?z>?5?)it is found thatρ〈cosθCO〉(?z)<0,implying the C=Obond slightly pointing to thebulk phase.

    Thedistributionof theorientationalangleθCOattheinterface is further investigated.The normalized probability density of θCOin asegmented surfacelayer z1<?z<z2isdefined by

    Eq.(26)is calculated in the bulk({?z<?10?})and at the interface({?10?<?z<?5?})of liquid EC,and the results of distribution are shown in Fig.4d.P(cosθCO,{?z<?10 ?})in the bulk is fl at over cosθCO,which indicates the random distribution in the bulk phase.At the interface region({?10?<?z<?5?}),however,the random distribution clearly breaks and cosθCOshows a preferable distribution at cosθCO≈ 0.This result indicates that at interface region EC molecules tend to orient themselves with their C=O bonds parallel to the interface,as illustrated in Fig.4b.The P(cosθCO,{?10?<?z<?5?})profi le at the surface shows slightly larger probability for cosθCO< 0 than for cosθCO> 0,which is consistent with the negative ρ〈cosθCO〉(?z)at the topmost layer in Fig.4c.These results of molecular orientation will be used to discuss the phase information of the SFGspectra in the following.

    5.4 SFG spectra

    Fig.5 showsthecalculated and experimental SFGspectrafor the ECvapor-liquid interfaceunder the ssp and sps polarization combinations.Itcanbeseenthatboth ssp and sps polarized SFG spectra in the C=O stretching region show two-peak structure,which is well reproduced by the calculated results.Besides,the calculated spectra also well reproduce the relative intensity between ssp and sps spectra.While the two-peak structure of the C=Oband iscommonly seen in the IR,Raman(Fig.3)and SFG(Fig.5),it is interesting to fi nd in Fig.5 that the higher frequency peak isstronger than thelower frequency onein both the ssp and sps SFGspectra,which isin contrast to the IRand Raman spectrain Fig.3.Themechanism of thediあerencein the bulk and interfacespectrawill bediscussed in Section 5.5.

    Fig.5 Calculated(a)and experimental(b)SFG intensity spectra of EC vapor-liquid interfaceunder ssp(red)and sp s(blue)polarizations.Themain peak frequencies are labeled.

    Fig.6 (a)Calculated Im[]spectra of EC interface.Diあerent colorsdenotethecalculated resultsfrom diあerent interfaceregions.(b)Self part spectrum Im[]by Eq.(18).

    Regarding the peak position,one fi nds a noticeable discrepancy in the higher frequency peak of ssp spectrum between simulation(a,1819 cm?1)and experiment(b,1790 cm?1).The reason of this discrepancy is not clear at present.Further investigation including the bulk quadrupole eあects on the SFGmay be required to solve this discrepancy.

    Table4 Derivativesof dipolemoment(?μr/?q CO)mol and polarizability(?αp q/?q CO)mol with respect to the C=O stretching normal mode q CO of EC molecule.

    The prefactor is analyzed by calculating(?μr/?qCO)moland(?αpq/?qCO)molof an EC molecule under the molecule-fi xed coordinates(ξ,η,ζ)in Fig.1.The ζaxis is along the carbonyl C1=O1bond,the η axis is normal to the carbonate O2–C1–O3plane,and the remaining ξaxis is given to be orthogonal to the other two.Table 4 summarizes thecalculated(?αpq/?qCO)moland(?μr/?qCO)molof the EC molecule with the B3LYP/aug-cc-pVTZ.The transition dipole is dominated by theζcomponent,and the transition polarizability is dominated by the diagonal,positive components.Therefore,the sign of Im[] ~ (?αyy/?qCO)(?μz/?qCO)in the space-fixed coordinates is mainly determined by the direction of theζaxis or the tilt angleθCO(see Fig.4b)33.Since the(?μζ/?qCO)molcomponent is negative in Table 4,the positive Im[]indicates negative cosθCO(downward C=O bond).This is fully consistent to the orientation analysis in Fig.4c,indicating that EC molecules exhibit slightly C=O downward orientation at the topmost layer of interface.

    5.5 Fermiresonance analysis

    The Fermi splitting of the C=O stretching band of EC molecule is analyzed with EPSA procedure in Section 4.3,where the harmonic force constant k4,4is changed to sk4,4with a scaling factor s(=1.0–1.12),and the responses of the line shape of C=O stretching band of IR,Raman and SFG are displayed in Fig.7(a–c),respectively.s=1.0(black lines)denotes the original spectra,while s=1.12(green lines)denotes the end point of the scaling where the overtone of skeletal stretching is completely shifted out of the C=O stretching region.It is clear that in the beginning of EPSA(black lines),all spectra show two band structures in C=O stretching region.As we shifted the overtone mode to higher frequency,one of the bands increases its intensity while the satellite band decreases.This is a clear characteristic of the merge of two Fermi split bands in the EPSA analysis,as illustrated in Fig.7d31.At the end of EPSA(green lines)one peak remains,which should beattributed to the fundamental of the C=Ostretching.between the two Fermi splitt components in IR,Raman and SFGspectra,as mentioned in Section 5.4.In the IRand Raman bands in Fig.3 the lower-frequency component is stronger,whereas vice versa in the SFG bands in Fig.5.Fig.7d illustrates this diあerence as the consequence of the Fermi splitting.The relative level of the fundamental of C=O stretching νC=Oand the overtone 2νskelchanges between the bulk and interface,due to the diあerent amount of the red shift of thefundamental.Thechangeof therelativelevelscausesthe reversal of relative intensities of the Fermi splitting components in the IR and Raman spectra in the bulk liqud and in the SFGspectraat theinterface.

    Fig.7 Calculated IR(a),Raman(b),and Im[](c)spectraof liquid ECin the EPSA procedure.Diあerent colorsdenote diあerent scaling factors s on theforceconstant k 4,4.(d)Fermi resonancemechanism of thethe C=O stretching fundamental(νC=O)and theovertoneof skeletal stretching(2νskel).Theupper and lower panels indicatethecasesof s=1.0 and s=1.12,respectively.Left panels illustrate theenergy diagramsof split levels,while the right panelsillustratetheconsequent lineshapes.Theblack and bluelines correspond to thecasesof thebulk liquid(IR,Raman)and interface(SFG),respectively.

    One notices in Fig.7 that the remaining bands at s=1.12 have slightly diあerent peak frequencies among the IR,Raman and SFGspectra.In particular,thepeak frequency of SFGband islocated atat~1815 cm?1,slightly abovethepeak frequencies of IR and Raman.This diあerence is caused by the diあerent solvation environment between the bulk phase and interface.The solvation at the interface is not as strong as that in the bulk phase,sincethemoleculestend tobeexposed to thevapor phase at the interface.The weak solvation at the interface results in a relatively small red shift of the C=O stretching peak than that in theinterior of bulk liquid.

    The diあerent solvation environment between bulk and interface also elucidates the reversal of relative intensities

    6 Conclusions

    In this work,the vapor liquid interface of EC is fully investigated by MD simulation and SFG measurement.A fl exible and polarizable model of EC molecule is proposed on the basis of CRK theory.The anharmonicity and Fermi resonance coupling of C=O stretching mode are properly included in our model.The experimental IR,Raman and SFG(ssp and sps)vibrational spectra of EC are also reported,which provide useful test of the performance of the proposed models.The model is validated by calculating thermodynamic properties,IR,Raman and SFG spectra,and the calculated resultsshow reasonableagreement with experimental data.The C=O stretching bands of the vibrational spectra are properly described by the present model,including the feature of Fermi splitting.

    MD simulations revealed that EC molecules show oscillating orientation asafunction of depth from theinterface.In the topmost layer of EC interface,EC molecules tend to orientate themselves with C=O bond parallel to the interface with a slightly downward C=O orientation.This orientation causes the positive Im[χ(2)]band in C=Ostretching region.

    The Fermi splitting of the C=O stretching band is fully analyzed by the EPSA procedure.The results demonstrate that the C=O stretching fundamental exhibits a larger amount of red shift in the bulk phase than at the interface,due to diあerent solvation environments.This diあerence explains the reversal of relative intensity of Fermi splitting components in the IR,Raman and SFGspectra.

    Acknowledgement:The numerical computations were performed with the supercomputers in Research Center for Computational Science,Okazaki,Japan.

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

    精品一区二区免费观看| 天天躁夜夜躁狠狠久久av| 最近2019中文字幕mv第一页| 小说图片视频综合网站| 久久午夜亚洲精品久久| av在线老鸭窝| 欧美成人一区二区免费高清观看| 国产 一区精品| 久久精品久久久久久噜噜老黄 | 一区福利在线观看| 激情 狠狠 欧美| 男人狂女人下面高潮的视频| 两个人视频免费观看高清| АⅤ资源中文在线天堂| 中国国产av一级| 性插视频无遮挡在线免费观看| 最近手机中文字幕大全| 国产69精品久久久久777片| 99热全是精品| 真实男女啪啪啪动态图| 97超碰精品成人国产| 久久鲁丝午夜福利片| 人妻夜夜爽99麻豆av| 国产免费一级a男人的天堂| 精品一区二区三区av网在线观看| 成人性生交大片免费视频hd| av在线亚洲专区| 97人妻精品一区二区三区麻豆| 天天一区二区日本电影三级| 中国美白少妇内射xxxbb| 国产探花极品一区二区| 久久人人爽人人片av| 精品久久久久久成人av| 色吧在线观看| 亚洲天堂国产精品一区在线| 久久欧美精品欧美久久欧美| 国产视频一区二区在线看| 成人美女网站在线观看视频| 免费av观看视频| 日韩,欧美,国产一区二区三区 | 最近中文字幕高清免费大全6| 久久人人爽人人爽人人片va| 男人舔奶头视频| 日韩强制内射视频| 国产免费一级a男人的天堂| 日日撸夜夜添| 国产精品久久久久久久电影| 亚洲av五月六月丁香网| 精品久久久久久久久av| 网址你懂的国产日韩在线| 国产激情偷乱视频一区二区| 色吧在线观看| 免费看a级黄色片| 国内精品久久久久精免费| 尾随美女入室| 晚上一个人看的免费电影| 熟女人妻精品中文字幕| 欧美区成人在线视频| 国产亚洲av嫩草精品影院| 午夜视频国产福利| 乱系列少妇在线播放| 亚洲美女黄片视频| 中文资源天堂在线| 国产色婷婷99| 99热这里只有是精品50| 日本欧美国产在线视频| 久久久久九九精品影院| 久久久久久大精品| 日本免费a在线| 免费搜索国产男女视频| 一本久久中文字幕| 人人妻人人澡欧美一区二区| 观看免费一级毛片| 久久久精品94久久精品| 色吧在线观看| 97人妻精品一区二区三区麻豆| 久久人人精品亚洲av| 亚洲美女黄片视频| 18禁黄网站禁片免费观看直播| 午夜视频国产福利| 午夜免费男女啪啪视频观看 | 老司机福利观看| 国产精品永久免费网站| 中国美白少妇内射xxxbb| 91久久精品电影网| 国产女主播在线喷水免费视频网站 | 91久久精品国产一区二区成人| 日本a在线网址| 哪里可以看免费的av片| 亚洲精品国产成人久久av| 观看美女的网站| 欧美成人a在线观看| 在线免费观看的www视频| 久久国产乱子免费精品| 久久精品夜夜夜夜夜久久蜜豆| 国产男人的电影天堂91| 欧美一区二区国产精品久久精品| а√天堂www在线а√下载| 精品免费久久久久久久清纯| 国内精品宾馆在线| 亚洲欧美精品综合久久99| 精品国内亚洲2022精品成人| 免费看光身美女| a级毛片免费高清观看在线播放| 乱系列少妇在线播放| 国内久久婷婷六月综合欲色啪| 久久精品国产亚洲av涩爱 | 嫩草影院入口| 99热这里只有精品一区| 久久久久久久久中文| 插逼视频在线观看| 天天一区二区日本电影三级| 欧美极品一区二区三区四区| av在线天堂中文字幕| 午夜影院日韩av| 亚洲欧美日韩无卡精品| 亚洲性久久影院| 亚洲欧美精品综合久久99| 久久精品国产亚洲av天美| 亚洲av第一区精品v没综合| 色吧在线观看| 久久人人精品亚洲av| 日韩中字成人| 国产成人精品久久久久久| 看十八女毛片水多多多| 插逼视频在线观看| 精品人妻视频免费看| 午夜影院日韩av| 日本黄大片高清| 欧美高清成人免费视频www| 欧美另类亚洲清纯唯美| 女的被弄到高潮叫床怎么办| 国产久久久一区二区三区| 欧美极品一区二区三区四区| 午夜久久久久精精品| 国产色爽女视频免费观看| 欧美绝顶高潮抽搐喷水| 国产极品精品免费视频能看的| 日韩精品青青久久久久久| 久久久久久伊人网av| 欧美一区二区亚洲| 国产高清视频在线播放一区| 国产亚洲精品久久久久久毛片| 又粗又爽又猛毛片免费看| 网址你懂的国产日韩在线| 久久久久性生活片| 日本爱情动作片www.在线观看 | 欧美成人精品欧美一级黄| 一个人观看的视频www高清免费观看| 免费看美女性在线毛片视频| 久久久久久久午夜电影| 激情 狠狠 欧美| 免费看美女性在线毛片视频| 日本爱情动作片www.在线观看 | 国产精品不卡视频一区二区| 成熟少妇高潮喷水视频| 亚洲四区av| 久久精品夜夜夜夜夜久久蜜豆| 成人永久免费在线观看视频| 极品教师在线视频| 久久午夜亚洲精品久久| 国产美女午夜福利| 亚洲av不卡在线观看| 久久久久性生活片| 国产又黄又爽又无遮挡在线| 99热这里只有精品一区| 淫妇啪啪啪对白视频| 搡老熟女国产l中国老女人| h日本视频在线播放| 欧美日韩乱码在线| 免费看a级黄色片| 亚洲专区国产一区二区| 悠悠久久av| 日韩中字成人| 天堂动漫精品| 亚洲精品日韩av片在线观看| 免费看光身美女| 日本免费一区二区三区高清不卡| 91麻豆精品激情在线观看国产| 嫩草影院新地址| 给我免费播放毛片高清在线观看| 午夜福利成人在线免费观看| 中文字幕av在线有码专区| 大型黄色视频在线免费观看| 亚洲欧美成人综合另类久久久 | 少妇裸体淫交视频免费看高清| 久久久成人免费电影| 国产真实乱freesex| 成人鲁丝片一二三区免费| av专区在线播放| 日韩在线高清观看一区二区三区| 日日摸夜夜添夜夜添小说| 在线观看免费视频日本深夜| av在线观看视频网站免费| 亚洲精品亚洲一区二区| 欧美精品国产亚洲| 午夜老司机福利剧场| 亚洲国产精品久久男人天堂| 不卡视频在线观看欧美| 天堂√8在线中文| 久99久视频精品免费| 国产又黄又爽又无遮挡在线| 国产精品亚洲美女久久久| 村上凉子中文字幕在线| 亚洲av熟女| 国产精品久久电影中文字幕| 国产精品99久久久久久久久| 久久久久久九九精品二区国产| 国产精品爽爽va在线观看网站| 女人十人毛片免费观看3o分钟| 99热全是精品| 亚洲精品乱码久久久v下载方式| 小说图片视频综合网站| 国产男人的电影天堂91| 欧美zozozo另类| 中文字幕精品亚洲无线码一区| 真人做人爱边吃奶动态| .国产精品久久| 久久久久久伊人网av| 亚洲人与动物交配视频| 国产色爽女视频免费观看| 成年版毛片免费区| 亚洲精品一区av在线观看| 亚洲av免费在线观看| АⅤ资源中文在线天堂| 欧美xxxx黑人xx丫x性爽| 国产真实乱freesex| 啦啦啦啦在线视频资源| 国产成人a区在线观看| 国产黄片美女视频| 美女免费视频网站| 欧美一区二区国产精品久久精品| 91久久精品国产一区二区成人| 成人av在线播放网站| 91精品国产九色| 日韩大尺度精品在线看网址| 伊人久久精品亚洲午夜| 亚洲av中文av极速乱| 哪里可以看免费的av片| 美女 人体艺术 gogo| 亚洲精华国产精华液的使用体验 | 久久久久久久久中文| eeuss影院久久| 精品国产三级普通话版| 色播亚洲综合网| 别揉我奶头~嗯~啊~动态视频| 欧美高清性xxxxhd video| 日本黄色片子视频| 午夜日韩欧美国产| 色在线成人网| 国产精品女同一区二区软件| 女人被狂操c到高潮| 亚洲成a人片在线一区二区| 国产亚洲精品av在线| 少妇被粗大猛烈的视频| 免费大片18禁| 老司机午夜福利在线观看视频| av黄色大香蕉| 成人毛片a级毛片在线播放| 身体一侧抽搐| 精品久久久久久久久亚洲| 国产成人影院久久av| 日本在线视频免费播放| 麻豆一二三区av精品| 亚洲最大成人中文| 欧美日本视频| 中文资源天堂在线| 国产熟女欧美一区二区| 国产真实伦视频高清在线观看| 男女视频在线观看网站免费| 晚上一个人看的免费电影| 国产精品亚洲一级av第二区| 97超视频在线观看视频| 99热这里只有是精品在线观看| 白带黄色成豆腐渣| 亚洲av美国av| 成年免费大片在线观看| 国产精品一区二区三区四区免费观看 | 国产黄a三级三级三级人| 亚洲高清免费不卡视频| 国产日本99.免费观看| 国产在线男女| 日韩国内少妇激情av| 亚洲美女黄片视频| 国产精品精品国产色婷婷| 啦啦啦啦在线视频资源| 看片在线看免费视频| av在线观看视频网站免费| 国产精品日韩av在线免费观看| 日韩高清综合在线| 亚洲国产精品久久男人天堂| 人妻制服诱惑在线中文字幕| 日本成人三级电影网站| 一区福利在线观看| 欧美一区二区亚洲| 狂野欧美白嫩少妇大欣赏| 成人一区二区视频在线观看| 插逼视频在线观看| 国产不卡一卡二| 国产精品精品国产色婷婷| 一区二区三区免费毛片| 欧美bdsm另类| 舔av片在线| 最近在线观看免费完整版| 成人美女网站在线观看视频| 久久国产乱子免费精品| 久久韩国三级中文字幕| 老女人水多毛片| 青春草视频在线免费观看| 波野结衣二区三区在线| 日韩三级伦理在线观看| 亚洲欧美中文字幕日韩二区| 免费看光身美女| 精品国产三级普通话版| 男人的好看免费观看在线视频| 成人欧美大片| 村上凉子中文字幕在线| 久久热精品热| 日本色播在线视频| 少妇的逼好多水| 日本免费a在线| aaaaa片日本免费| 蜜桃久久精品国产亚洲av| 全区人妻精品视频| 99国产极品粉嫩在线观看| or卡值多少钱| 午夜老司机福利剧场| 99热这里只有是精品在线观看| 欧洲精品卡2卡3卡4卡5卡区| 日韩一区二区视频免费看| 无遮挡黄片免费观看| 丰满乱子伦码专区| 久久精品久久久久久噜噜老黄 | 成年av动漫网址| 国产老妇女一区| 一区福利在线观看| 亚洲av免费高清在线观看| av黄色大香蕉| 精品不卡国产一区二区三区| 久久精品夜夜夜夜夜久久蜜豆| 俺也久久电影网| 国产av麻豆久久久久久久| 国产熟女欧美一区二区| 久久久久久久久久黄片| 99在线视频只有这里精品首页| 啦啦啦韩国在线观看视频| 亚洲av成人精品一区久久| 三级经典国产精品| 色综合亚洲欧美另类图片| 禁无遮挡网站| 国产伦一二天堂av在线观看| 色哟哟哟哟哟哟| 香蕉av资源在线| 中国美女看黄片| 啦啦啦韩国在线观看视频| av.在线天堂| 香蕉av资源在线| 美女xxoo啪啪120秒动态图| 亚洲成a人片在线一区二区| 高清毛片免费看| 欧美日韩国产亚洲二区| av在线亚洲专区| 日韩欧美国产在线观看| 精品无人区乱码1区二区| 免费观看的影片在线观看| 亚洲人成网站在线播| 国产精品电影一区二区三区| 波多野结衣高清无吗| 一级毛片电影观看 | 色视频www国产| 日本成人三级电影网站| 波多野结衣高清作品| 51国产日韩欧美| 大香蕉久久网| 国产男人的电影天堂91| 亚洲成a人片在线一区二区| 黄色一级大片看看| 精品久久国产蜜桃| 日韩 亚洲 欧美在线| 99久国产av精品| 亚洲最大成人手机在线| 成年女人毛片免费观看观看9| 亚洲国产精品国产精品| 波野结衣二区三区在线| 精品乱码久久久久久99久播| 最近2019中文字幕mv第一页| 久久久欧美国产精品| 少妇丰满av| 欧美zozozo另类| 久久精品国产亚洲av天美| 国产极品精品免费视频能看的| 欧美性猛交黑人性爽| 亚洲国产精品久久男人天堂| 午夜福利高清视频| 成人无遮挡网站| 欧美潮喷喷水| 亚洲国产精品sss在线观看| 国产亚洲av嫩草精品影院| 日韩成人伦理影院| 色吧在线观看| 国产视频内射| 国产视频一区二区在线看| 亚洲精品粉嫩美女一区| 午夜激情福利司机影院| 不卡视频在线观看欧美| 国产成人福利小说| 精品99又大又爽又粗少妇毛片| 你懂的网址亚洲精品在线观看 | 久久久久久久亚洲中文字幕| 国产三级在线视频| 亚洲高清免费不卡视频| 亚洲国产日韩欧美精品在线观看| 国内揄拍国产精品人妻在线| 最新在线观看一区二区三区| 日韩在线高清观看一区二区三区| 精品人妻偷拍中文字幕| 中国国产av一级| 久久精品夜色国产| 亚洲专区国产一区二区| 精品一区二区免费观看| 91久久精品国产一区二区成人| 免费av观看视频| 两个人的视频大全免费| 女生性感内裤真人,穿戴方法视频| 亚洲欧美精品综合久久99| 亚洲,欧美,日韩| 99久久中文字幕三级久久日本| 性色avwww在线观看| 99久久成人亚洲精品观看| 22中文网久久字幕| 最新中文字幕久久久久| 一级毛片久久久久久久久女| 国产久久久一区二区三区| 日本与韩国留学比较| 亚洲国产精品成人久久小说 | 女同久久另类99精品国产91| 嫩草影院精品99| 亚洲在线观看片| 国产av在哪里看| 午夜爱爱视频在线播放| 国产亚洲精品久久久久久毛片| 狂野欧美激情性xxxx在线观看| 国产在线精品亚洲第一网站| 国产高清激情床上av| 亚洲精品456在线播放app| 日韩av不卡免费在线播放| 国产精华一区二区三区| 黄色配什么色好看| 亚洲精品粉嫩美女一区| 级片在线观看| 婷婷六月久久综合丁香| 麻豆精品久久久久久蜜桃| 在线观看66精品国产| 亚洲七黄色美女视频| 成人一区二区视频在线观看| 97超视频在线观看视频| 亚洲人成网站在线播| 亚洲国产精品sss在线观看| 日韩亚洲欧美综合| 欧美绝顶高潮抽搐喷水| 欧美性猛交黑人性爽| 欧美最黄视频在线播放免费| 亚洲欧美日韩无卡精品| 一级毛片我不卡| 国产精品一区二区三区四区免费观看 | 尾随美女入室| 国产成人91sexporn| 中文字幕精品亚洲无线码一区| 国产高潮美女av| 免费看美女性在线毛片视频| 最近最新中文字幕大全电影3| 色综合亚洲欧美另类图片| 国产91av在线免费观看| 午夜免费男女啪啪视频观看 | 婷婷六月久久综合丁香| 一级黄片播放器| 午夜福利在线观看免费完整高清在 | 亚洲av五月六月丁香网| 日本精品一区二区三区蜜桃| 韩国av在线不卡| 在线天堂最新版资源| 我要搜黄色片| 精品久久久久久久久久免费视频| av在线观看视频网站免费| 春色校园在线视频观看| 91狼人影院| 天堂动漫精品| 晚上一个人看的免费电影| 成人av一区二区三区在线看| 嫩草影院新地址| 亚洲中文日韩欧美视频| 高清日韩中文字幕在线| 成人二区视频| 亚洲欧美精品自产自拍| 性插视频无遮挡在线免费观看| 欧美潮喷喷水| 久久精品国产清高在天天线| 99热这里只有是精品在线观看| 亚洲国产精品合色在线| 精品久久久久久久久av| 最新在线观看一区二区三区| 不卡一级毛片| 18禁裸乳无遮挡免费网站照片| 久久精品国产亚洲av涩爱 | 精品一区二区免费观看| 少妇人妻精品综合一区二区 | 亚洲av免费高清在线观看| 久久久久国产网址| 中文字幕熟女人妻在线| 国产探花在线观看一区二区| 色吧在线观看| 丝袜喷水一区| 在线观看一区二区三区| 免费看日本二区| 自拍偷自拍亚洲精品老妇| 麻豆成人午夜福利视频| 日本-黄色视频高清免费观看| 久久综合国产亚洲精品| 国产aⅴ精品一区二区三区波| 欧美日韩精品成人综合77777| 最近2019中文字幕mv第一页| 亚洲电影在线观看av| 日日摸夜夜添夜夜爱| 亚洲,欧美,日韩| 舔av片在线| 麻豆一二三区av精品| 亚洲精品粉嫩美女一区| 国内精品久久久久精免费| 丰满的人妻完整版| 九九爱精品视频在线观看| 日日啪夜夜撸| 麻豆国产av国片精品| 一级黄色大片毛片| 午夜福利在线观看吧| 国产精品国产高清国产av| 亚洲欧美日韩卡通动漫| 国产毛片a区久久久久| 国产黄色小视频在线观看| 国产伦一二天堂av在线观看| 国产精品一区二区三区四区久久| 观看免费一级毛片| 精品午夜福利在线看| 精品久久久噜噜| 长腿黑丝高跟| 成人欧美大片| 亚洲最大成人av| 真人做人爱边吃奶动态| 舔av片在线| 91在线观看av| 最近中文字幕高清免费大全6| 搞女人的毛片| av视频在线观看入口| 国产成人freesex在线 | 又爽又黄a免费视频| 色综合亚洲欧美另类图片| 精品人妻一区二区三区麻豆 | 日韩欧美免费精品| 成人高潮视频无遮挡免费网站| 美女内射精品一级片tv| 国产午夜福利久久久久久| 99久久精品国产国产毛片| 全区人妻精品视频| 日本欧美国产在线视频| 亚洲中文字幕日韩| 国产欧美日韩一区二区精品| 22中文网久久字幕| 午夜激情欧美在线| 亚洲精品日韩在线中文字幕 | 91精品国产九色| eeuss影院久久| 99久久成人亚洲精品观看| 色吧在线观看| 国产精品电影一区二区三区| 中国国产av一级| 国产成人freesex在线 | 欧美激情国产日韩精品一区| av在线播放精品| 网址你懂的国产日韩在线| 欧美xxxx黑人xx丫x性爽| 精品一区二区三区av网在线观看| 精品久久久久久久久亚洲| 一进一出抽搐动态| 亚洲中文字幕日韩| 国产真实乱freesex| 国产精品久久久久久久电影| 少妇猛男粗大的猛烈进出视频 | 两个人视频免费观看高清| 白带黄色成豆腐渣| videossex国产| 精品福利观看| 身体一侧抽搐| 一进一出抽搐动态| 色在线成人网| 中文字幕免费在线视频6| 在线观看一区二区三区| 性欧美人与动物交配| 日本免费一区二区三区高清不卡| 99热6这里只有精品| 久久这里只有精品中国| 国产激情偷乱视频一区二区| 嫩草影院新地址| 夜夜看夜夜爽夜夜摸| 午夜精品在线福利| 波野结衣二区三区在线| videossex国产| 国产精品永久免费网站| h日本视频在线播放| 欧美xxxx黑人xx丫x性爽| 亚洲色图av天堂| 欧美xxxx黑人xx丫x性爽| 在线观看av片永久免费下载| 国产亚洲精品久久久com| 国产精品综合久久久久久久免费| 国产真实乱freesex| 久久久欧美国产精品|