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

    First-principles calculations of K-shell x-ray absorption spectra for warm dense ammonia*

    2021-05-24 02:23:36ZiLi李孜WeiJieLi李偉節(jié)CongWang王聰DafangLi李大芳WeiKang康煒XianTuHe賀賢土andPingZhang張平
    Chinese Physics B 2021年5期
    關(guān)鍵詞:張平

    Zi Li(李孜), Wei-Jie Li(李偉節(jié)), Cong Wang(王聰),?, Dafang Li(李大芳), Wei Kang(康煒),Xian-Tu He(賀賢土),, and Ping Zhang(張平),,?

    1Institute of Applied Physics and Computational Mathematics,Beijing 100088,China

    2Center for Applied Physics and Technology,Peking University,Beijing 100871,China

    Keywords: first-principles,warm dense,x-ray absorption spectrum,ammonia

    1. Introduction

    Ammonia is one of the major constituents of the giant planets in our solar system. Specifically, in Uranus and Neptune, it is supposed that there is a thick layer of “hot ices”,predominantly made up of water,hydrocarbons,and ammonia in solar proportions (molar fractions: 56% H2O, 36% CH4,and 8% NH3).[1]From theoretical estimations, the pressure in this layer could reach 300 GPa, and the temperature is between 2000 K and 5000 K as evinced by the warm adiabatic model,[2–4]or between 7000 K and 14000 K as evinced by a model including a thermal layer.[4,5]The knowledge of the physico–chemical properties of ammonia in these extreme conditions is critical to the understanding of the observable properties (gravitational moments, atmospheric composition,and magnetic field)of the planets.[2,6–8]

    At low temperature and ambient pressure,ammonia crystallizes in a proton-ordered cubic phase(phase I),[9,10]where each molecule forms six equivalent weak hydrogen bonds,bent through an angle of 160°. Within 20 GPa and 400 K,a total of five solid phases of ammonia have been discovered,where phase II[11]and phase III[12]are rotationally disordered(plastic) hexagonal close-packed (hcp) and face-centered cubic (fcc) structures, respectively, and phase IV[13,14]and phase V[15–18]are proton-ordered orthorhombic structures. At higher temperature and pressure, a new phase with superionic characteristics above 700 K and 57 GPa[19]was found experimentally, which coincides with the prediction made in a previous theoretical study;[20]however, another experiment could not identify this superionic phase at pressures reaching 60 GPa, and it is speculated that this new phase may only be realized at a higher pressure, due to the chemical dissociation of ammonia.[21]In most of these static experiments on the exploring of the phase diagram, diamond(or tungsten carbide) anvil cells are usually employed to achieve such high pressures,[12–19,21]and the techniques of neutron diffraction(for deuteroammonia ND3),[9,11,13,18]x-ray diffraction[12,16,19,21]and Raman spectroscopy[14,15,17–19,21]are applied to detect the atomic structure and to confirm the phase transitions.

    Shock-wave experiments were also conducted on ammonia and mixtures containing ammonia, to establish the equation of states and salient conductivity data under warm dense conditions.[8,22–27]High-explosive techniques were first used to obtain the Hugoniot curve of ammonia up to 39 GPa,and no transition is observed therein within this pressure region.[22]A two-stage light-gas gun was used to reach a pressure of 64 GPa.[23]The first shock temperature measurement on ammonia was performed at pressures of up to 61 GPa to diagnose the behavior under these extreme conditions.[24]The conductivities of the ammonia and the “synthetic Uranus” (a mixture of water, isopropanol, and ammonia)were measured under shock-compressed conditions,[8,25,26]and the low-pressure conductivity was found to be dominated by conduction via hydrogen ions from the dissociation of the molecules.[26]A recent experiment employing laser-driven shock compression could achieve extreme conditions up to 260 GPa and 2.5×104K for a mixture containing ammonia.[27]In these dynamic experiments,the pressure and density of the thermodynamic states are usually determined by measuring shock velocities and fluid velocities, and the temperature is identified by measuring the radiation of shock-heated samples.

    On the theoretical side,first-principles molecular dynamics methods have been broadly utilized to examine the phase diagram,equation of states,and chemical reactions of ammonia under extreme conditions. A superionic phase of ammonia, characterized by highly-mobile hydrogen ions diffusing through a nitrogen lattice with hcp structure, was predicted above 1200 K and 60 GPa,[20]a decade before the experimental observation,[19]confirming the predictive power of theoretical methods. The equation of states and the phase diagram of ammonia were then determined up to 330 GPa and 1.0×104K.[28]We investigated the thermophysical and transport properties of shocked liquid ammonia with the pressure and temperature reaching 1.3 TPa and 1.2×105K,and found that liquid ammonia undergoes three structural changes of regime along the principal Hugoniot curve: molecular fluid, mixture,and plasma.[29,30]A combined experimental and theoretical study has examined the physical properties of “synthetic Uranus”, indicating the formation of large clusters of carbon and nitrogen under the thermal conditions deep within the interior of Uranus.[26]At zero temperature,an ionic crystal,characterized by alternate layers of NH+4and NH2ions, was suggested to exist at pressures greater than 90 GPa;[31]other possible zero temperature structures were also discussed.[32]

    In order to detect warm dense matter in static and dynamic experiments,the x-ray absorption spectrum(XAS)provides high temporal and spatial resolution.[33–36]The XAS directly relates to the electron transition from the atomic inner shell to unoccupied states, which provides important information on electronic structures and atomic environments. Previous theoretical studies have shown that, under extreme conditions, the changes in the XAS could signal a phase transition,[37,38]nonmetal–metal transition,[39,40]magnetic ordering,[41]or temperature evolution.[42–44]Focusing on nitrogen, the x-ray spectrum has been employed to determine the chemical bonding properties of nitrogen atoms in organic molecules[45,46]and nano-materials;[47]however, there is little information about the XAS of ammonia under extreme conditions,which hinders the application of this powerful tool in relevant experiments. Here, to reveal the importance of XAS in matter detection and provide information for future experimental studies, we perform first-principles calculations to determine the nitrogen K-shell XAS for the ammonia under thermodynamic conditions along the principal Hugoniot curve. A correlation between the shape of XAS curves and the warm dense state is established,and a detailed analysis of atomic structures is undertaken to explain the XAS changes on these Hugoniot states.

    2. Computational method and models

    The simulations of first-principles molecular dynamics(MD)and the nitrogen K-shell XAS are performed on a series of 10 thermodynamic conditions along the principal Hugoniot curve of ammonia, which has been theoretically determined in our previous work.[30]The density of the conditions increases from 0.8 g/cc in condition 1, to 2.6 g/cc in condition 10, in increments of 0.2 g/cc for intermediate conditions; the corresponding temperature ranges from 250 K to 1.12663×105K, as labelled in Fig. 1(a), with external pressures increasing from 1.2 GPa to 1.3 TPa.[30]The VASP code[48,49]is employed for the MD simulations with a canonical (or NVT) ensemble to achieve the corresponding temperature and density, with the projector augmented wave(PAW)pseudopotentials[50]and the Perdew–Burke–Ernzerhof exchange–correlation functional.[51]The calculations are performed at the gamma point only,and the energy cutoff of the plane-wave basis is 550 eV. Each MD simulation evolves for 8 ps,with a timestep of 0.5 fs–1 fs. A cubic supercell containing 64 ammonia molecules is considered in these calculations,except under conditions 9 and 10,where 32 and 27 molecules are considered,respectively,due to the large computational effort involved.A total of five randomly selected snapshots from each MD trajectory under thermal equilibrium are selected for the analysis of the atomic structure and the determination of XAS.

    The K-shell XAS of the nitrogen atom in warm dense ammonia is determined by using the XSPECTRA package[52,53]supplied by the Quantum Espresso distribution.[54]The XAS corresponds to the electron transition from a core 1s orbital on a specified nitrogen atom (referred to as an absorbing atom) to the valance states of the whole system with a core hole on the absorbing atom. For the determination of the XAS, the Troullier–Martins form of norm-conserving pseudopotential[55]is employed and the 1s electron of the nitrogen atom is set as the frozen core. The gamma point is sampled in the Brillouin zone and the energy cutoff of the plane-wave basis is 120 Ry (1 Ry=13.605692312 eV). The calculation is conducted repeatedly by setting each nitrogen atom in the supercell as the absorbing atom at one time, and the final XAS curve(referred to as the total XAS)is an average over both atoms and snapshots. The better to describe the change of XAS under warm dense conditions, an additional all-electron calculation with a larger energy cutoff of 220 Ry is performed on each MD snapshot to determine the absolute energy position of the XAS curve (the eigenvalue of the nitrogen 1s orbital),and a pre-factor of Fermi–Dirac function is added to the XAS.Also,a rigid blue shift of 20 eV is made for all the XAS curves the better to reproduce the experimental results. All these procedures are the same as those described in our previous study on warm dense nitrogen.[56]

    3. Results and discussion

    The total XAS curves at the ten typical thermodynamic conditions are shown in Fig. 1(a). The main characteristic of the XAS curve is a peak structure located at around 400 eV, followed by a plateau above 420 eV, similar to the experimental results for nitrogen compounds under ambient conditions.[46,57]The position and shape of the peak change significantly with the conditions, making the XAS a good indicator of the warm dense states. According to the change of the XAS,the Hugoniot curve could be classified into three stages. The XAS peak is very sharp and high in conditions 1–4. While in conditions 5–7, the peak becomes broader and lower, and a shoulder structure could be observed at around 405 eV;in particular in condition 7,due to the peak broadening and the shoulder structure, the XAS curve is almost flat between 400 eV and 410 eV. Finally, under conditions 8–10,the peak is even lower and almost undistinguishable from the plateau in the higher energy region.

    In Fig.1(b),the relative peak intensity(peak height)and the peak position for the XAS curves are illustrated. The peak intensity decreases monotonically from conditions 1–10, and the intensity of condition 10 is only one-third of that of condition 1. For the peak positions, from conditions 1–7, the peak of the XAS curve shifts slightly to higher energy (i.e., a blue shift), from 399.2 eV in condition 1 to 399.9 in condition 4,and 401.8 eV in condition 7.This shift relates to the shift of the nitrogen 1s energy level,due to the stronger binding between the nucleus and the inner 1s electron at higher temperatures,as pointed out in our previous work.[58]The larger position shift in conditions 5–7 is also related to the presence of the shoulder structure at around 405 eV.While in conditions 8–10,the peak position undergoes a red-shift. This is partially attributed to the uncertainty of the peak position due to the low XAS intensity, as well as the formation of new products under extreme conditions. The width of the XAS peak is shown in Fig.1(c).Here, the width is defined as the energy difference from the peak maximum to the position where the XAS intensity decreases to 1/e(with e the base of the natural logarithm)times the maximum; both the widths from the left and the right of the peak maximum are determined. The width increases from conditions 1–7,due to the increase of structural disorder. The fluctuation of the left-hand width in conditions 8–10 is related to the blue shift of the peak position; the right-hand width in conditions 9–10 could not be determined due to the low peak intensity.

    Fig.1. The XAS in conditions 1–10 along the Hugoniot curve. (b)The peak position (black) and peak intensity (blue) of the XAS curves. (c) The XAS peak width from the left(black)and the right(blue)of the peak maximum.

    We then examine the contribution of the total XAS curve by each nitrogen atom (referred to as the atomic XAS). For each atomic XAS, we determine its peak position and right width as done for the total XAS curve and demonstrate their distributions for typical conditions in Figs.2(a)and 2(b). Under conditions 1,3,and 5,the peak position is mainly located at around 400 eV and is slightly delocalized with increasing temperature and density. In condition 7,the positions are distributed at almost equal intervals from 400 eV to 410 eV.The peak width is distributed around 1.5 eV for the thermodynamic conditions considered, and the percentage of the peak widths exceeding 2 eV slightly increases from conditions 1–7.

    The variety of peak positions and widths of the atomic XAS are attributed to different chemical environments of the nitrogen atoms. Here, we identify the local atomic structures at the Hugoniot conditions to understand the changes in the total XAS curves and deduce connections among XAS data,the atomic structure,and the thermodynamic conditions.Firstly,we examine the local atomic structures under Hugoniot conditions. In Fig. 3, the pair correlation functions (PCFs)of nitrogen–hydrogen (N–H), nitrogen–nitrogen (N–N), and hydrogen–hydrogen (H–H) in typical conditions are plotted,and changes in the nature of the bonding along the Hugoniot curve could be observed.In condition 3,clear peaks are shown in the PCF curves. The PCF of N–H peaks at 1.0 ?A, corresponding to the N–H bond;the maximum of N–N PCF around 3.0 ?A reflects the typical separation of NH3molecules;the H–H PCF shows a peak at 1.6 ?A,related to the H–H separation in one NH3molecule. In conditions 5 and 7,an extra peak of N–N PCF at ~1.2 ?A indicates the formation of N–N bond, and a shoulder character in H–H PCF at 0.75 ?A indicates a small amount of H2molecules. In condition 10,no peak appears in the PCFs,indicating a disordered local structure at such a high temperature(105K).

    Fig. 2. The distributions of peak position (a) and peak width (b) for conditions 1,3,5,and 7.

    The products under the thermodynamic conditions are assessed: the possible products are classified into three types:the pristine ammonia molecule NH3, a fragment NHx(with x=0, 1, or 2) from dissociation, and a nitrogen cluster (Ncluster), where two or more nitrogen atoms are bonded together (hydrogen atoms may also be present therein). Here,two atoms are recognized as being bound if the interatomic distance is smaller than a specified cutoff radius, which is set to 15% greater than the typical bond length. For the N–H bond, the cutoff radius is 1.15 ?A, corresponding to typical bond lengths of 1.02 ?A.For the N–N bond,two cutoff values of 1.7 ?A and 1.4 ?A are considered,corresponding to N–N single bond(1.46 ?A)and double bond(1.2 ?A).The value of 1.7 ?A is used to justify if a nitrogen cluster is formed,and in a nitrogen cluster,if one of the N–N bond lengths is less than 1.4 ?A,it is called a type-I nitrogen cluster(or N-cluster I),otherwise,it is called a type-II nitrogen cluster(or N-cluster II).

    Fig. 3. The PCFs in (a) condition 3, (b) condition 5, (c) condition 7, and(d)condition 10. The black,red,and blue curves correspond to the PCFs of N–H,N–N,and H–H,respectively.

    In Fig. 4(a), the percentage of nitrogen atoms in different types of products is demonstrated. The percentage of NH3is 100% in condition 1 and decreases to 67% in condition 4(while remaining dominant), and to 11% in condition 7; finally,in conditions 8–10,there is only a trace amount of NH3molecules (less than 3%). The fragment NHxappears from condition 2,reaching a maximum of 57%in condition 7,and then decreases from condition 8(42%)to condition 10(23%).The formation of N-clusters begins in condition 5, and the amount increases along the Hugoniot curve (39% for type-I and 36% for type-II in condition 10), with similar ratios for type-I and type-II N-clusters in each condition. With increasing temperature and density along the Hugoniot curve,firstly,the N–H bond breaks and the NHxfragment appears(in condition 2);then,the NHxfragments begin to form N-clusters in condition 5,leading to the proportional amount of NHxchanging slightly from conditions 5–7; finally, in conditions 8–10,the proportional amount of NHxfragments decreases due to the formation of more N-clusters. During the three stages representing change to the XAS curves, the dominant product is the molecule NH3(conditions 1–4), the fragment NHx(conditions 5–7)and the N-cluster(conditions 8–10), respectively. The three stages correspond to the three structural regimes pointed out in our previous work,[29]i.e., molecular fluid regime (dominated by NH3molecules), mixture regime(presence of NH3, NHx, and a small amount of N-clusters),and a plasma regime(characterized by an extremely high temperature and short-lived chemical bonds).

    Fig. 4. The analysis of the atomic structure of the MD snapshots in thermodynamic conditions 1–10. (a) The percentages of NH3 molecules, NHx fragments, and N-clusters; (b) the distributions of the distance dNH in each ammonia molecule at typical conditions; (c) the percentages of isolated N atoms, NH, and NH2; (d) the percentages of N-clusters with different numbers of nitrogen atoms.

    We then examine the structural deformation of NH3molecules under different thermodynamic conditions. Figure 4(b) shows the distribution of the distance dNHbetween the nitrogen atom and the hydrogen plane(the plane where the three hydrogen atoms of the molecule located), which is labelled as dNHherein. In conditions 1 and 2, dNHis mainly distributed around 0.4 ?A, close to its value in the gaseous phase (0.38 ?A, corresponding to the H–N–H angle of 107°).While from condition 3, the distance varies from 0.0 ?A to 0.6 ?A, indicating severe deformation of the NH3structure.Figure 4(c) shows the ratios of various types of fragments NHxfor x=0, 1, and 2. The NH2is first formed in condition 2 (2%) and reaches a maximum of 41% in condition 5;the maximum percentages of the NH fragment and isolated N atom are achieved at a higher temperature and density (conditions 7 and 8, respectively), implying that the NH and N atom should be formed from the decomposition of the NH2.In conditions 9 and 10, the NHxfragment-content decreases to 5%–10%, due to the formation of nitrogen clusters. Figure 4(d) shows the distribution of N-clusters with two to five nitrogen atoms, where the solid part and the empty part of the bar indicate the type-I and type-II clusters, respectively.In conditions 5–8,the N-clusters mainly contain two nitrogen atoms, while in conditions 9 and 10, larger N-clusters with three or more nitrogen atoms dominate, due to the formation of more N–N bonds at higher temperature and density;type-I and type-II clusters show similar distributions on the cluster size in each thermodynamic condition. Although the chemical bonds are short-lived in conditions 9 and 10 as pointed out in our previous work,[29]there should be still a significant number of N-clusters at each MD snapshot due to the extremely high temperature and density,which could affect the shape of the XAS curve.

    Next, we construct the connection between the atomic structure and the XAS curve.In Figs.5(a)–5(d),the total XAS under typical conditions is decomposed into the contributions from the products of NH3, NHx, and N-clusters (referred to as the decomposed XAS), which is determined by the summation of the atomic XAS of the corresponding product type.In all these conditions,the decomposed XAS curves from the same type of product are similar;specifically,the XAS curves contributed to by NH3molecules show a sharp peak at around 400 eV; the peak widths of XAS curves from the NHxfragment are much larger than those of NH3;for the type-I cluster,the XAS peak is replaced by a broad plateau ranging from 400 eV to 410 eV;the XAS curves of type-II clusters are similar to those of NH3molecules (a sharp peak), but slightly shifted to a lower energy.We also examine the electron density of states(DOS)to show the correlation between the XAS and the electronic structure. In Figs. 5(e) and 5(f), we show the local DOS curves contributed by the products of NH3, NHx,and N-clusters, which are obtained by the summation of the corresponding atomic local DOS.Here,for a specified type of product,the characteristics of the DOS curves above the Fermi energy are very similar to those of the XAS curves: compared to the local DOS of the NH3molecule, the peak of the DOS from the NHxis broader and the peak of the DOS from the type-II cluster shifts to a lower energy;the plateau structure of XAS curves from the type-I cluster is also observed in its DOS curves.

    The shape of the total XAS could be understood from the decomposed XAS curves. In conditions 1–4, NH3molecule is the dominant product,giving the total XAS a sharp peak at around 400 eV;the increase(decrease)of the peak width(intensity)is attributed to the deformation of the molecule,as illustrated in Fig.4(b).The broad peak in conditions 5–7 and the shoulder structure slightly above 400 eV is caused by the dominant contribution from the NHx, as well as from the plateau structure of type-I clusters in condition 7. In conditions 8–10,the large number of type-I clusters and the extremely high temperature smear out the peak of the XAS curves,and the presence of the type-II clusters slightly shifts the XAS to the lower energy region. The effects of N-clusters in warm dense ammonia are different from that in warm dense nitrogen, where the formation of the N-clusters induces multiple small XAS peaks at energies between 405 eV and 420 eV, as examined in our previous theoretical study.[56]This difference may be attributed to the presence of hydrogen atoms in warm dense ammonia,that could saturate the dangling bonds in N-clusters and thus eliminate the possible XAS peaks.

    Fig. 5. The contributions to the total XAS curve from ammonia molecules NH3 (black), fragments NHx (red), type-I (blue), and type-II (green) Nclusters in (a) condition 3, (b) condition 6, (c) condition 7, and (d) condition 9;the total XAS is plotted as a black dotted curve in each figure,(e)and(f):the electron local DOS from various kinds of products in(e)condition 6 and(f)condition 7;the Fermi energy is set to zero.

    Furthermore, the effects of the characteristic parameters(peak position, width, and intensity) of atomic XAS on the shape of the decomposed XAS are explored. Here,we reconstruct the decomposed XAS by replacing the atomic XAS with a model Gaussian function(GF).In model I,the center of the GF is located at the peak position of the atomic XAS,and the width and height are set to an averaged value of all the atomic XAS at each thermodynamic condition. In model II, all the characteristic parameters of GF are the same as those of the atomic XAS curve. The reconstructed XAS curves predicted by the two models are shown in Figs.6(a)–6(d)in conditions 7 and 9. Compared with the first-principles results in Figs.5(c)and 5(d), important characteristics of the XAS curves could be reproduced by model I, such as the broad peak of NHx,the plateau structure of the type-I cluster, and the red shift of the XAS peak for the type-II cluster. When including all of the characteristic parameters in model II, the resulting XAS curves are similar to those derived from first-principles calculations, except for the plateau structure in the high energy region above 420 eV, since the GF used in the models does not consider the long tail of the XAS curves in this region.Since the total XAS curve is the sum of decomposed XAS curves, it is concluded that the distribution of the peak position is the dominant factor affecting the shape of the total XAS curves,and inclusion of all parameters pertaining to peak position,width,and height could reproduce the total XAS curves with accuracy. The peak position is related to the energy level of valance states around the specified nitrogen atom. For the NHxfragment and type-I N-cluster, the presence of the unpaired electron or the nitrogen-nitrogen double bond renders the energy level very sensitive to the atomic structure, thus a broad energy distribution and a large peak width could be expected in the XAS curve.

    Fig. 6. The reconstructed XAS contributions from NH3 (black), NHx (red),N cluster-I(blue),and N cluster-II(green)based on modelling methods. (a)model I in condition 7;(b)model II in condition 7;(c)model I in condition 9;(d)model II in condition 9.

    4. Conclusions

    In summary, we have performed first-principles molecular dynamics and XAS simulations for warm dense ammonia along the Hugoniot curve. The calculated nitrogen Kshell XAS curves show a main peak at around 400 eV, and with increasing temperature and density,the peak has a lower intensity and it broadens. According to the changes in the XAS curves,the Hugoniot curve could be classified into three stages. The analysis of the atomic structures indicates that,along the Hugoniot curve,the NH3molecules could firstly dissociate to NHxfragments then form nitrogen clusters (type-I and type-II). At the highest temperature and density considered(condition 10),the preserved NH3molecules account for less than 3%of the total,and the amounts of type-I and type-II nitrogen clusters are 39% and 36% of the total, respectively.For the three stages of the Hugoniot curve,the dominant products (in descending order) are: NH3molecules, NHxfragments, and nitrogen clusters, respectively. The contributions of the three types of products to the total XAS are quite different, inducing the change of the XAS along the Hugoniot curve. By reconstructing the XAS curves by using Gaussian functions, the distribution of the peak position of the atomic XAS is found to be the dominant factor affecting and determining the shape of the total XAS curve.

    猜你喜歡
    張平
    Magnetic ground state of plutonium dioxide: DFT+U calculations
    Topological properties of Sb(111)surface: A first-principles study
    嘰嘰喳喳的小喜鵲
    Magnetic phase diagram of single-layer CrBr3?
    High-pressure elastic anisotropy and superconductivity of hafnium:A first-principles calculation*
    這是你爺倆
    金秋(2020年16期)2020-12-09 01:41:50
    高中生英語(yǔ)學(xué)習(xí)焦慮現(xiàn)狀研究
    Enhanced spin-dependent thermopower in a double-quantum-dot sandwiched between two-dimensional electron gases?
    千年煙云大西安
    金秋(2019年2期)2019-04-16 09:14:42
    張平書法作品選
    日日爽夜夜爽网站| 亚洲国产精品sss在线观看| 人人妻人人澡欧美一区二区| 身体一侧抽搐| 成人午夜高清在线视频 | 国产成人欧美在线观看| 丝袜在线中文字幕| 午夜久久久久精精品| 色尼玛亚洲综合影院| 99国产精品一区二区蜜桃av| 在线天堂中文资源库| 久久婷婷成人综合色麻豆| 国产极品粉嫩免费观看在线| bbb黄色大片| 精品久久久久久久末码| 免费电影在线观看免费观看| 啦啦啦 在线观看视频| cao死你这个sao货| 久9热在线精品视频| 很黄的视频免费| 操出白浆在线播放| 51午夜福利影视在线观看| 桃红色精品国产亚洲av| 99热6这里只有精品| 欧美成狂野欧美在线观看| 真人一进一出gif抽搐免费| 99久久无色码亚洲精品果冻| 久久国产亚洲av麻豆专区| 欧美黑人精品巨大| 777久久人妻少妇嫩草av网站| 日韩国内少妇激情av| 50天的宝宝边吃奶边哭怎么回事| 日本在线视频免费播放| 热99re8久久精品国产| 一个人观看的视频www高清免费观看 | 窝窝影院91人妻| 亚洲国产欧洲综合997久久, | 久9热在线精品视频| 可以免费在线观看a视频的电影网站| 日韩欧美国产在线观看| 国产亚洲精品综合一区在线观看 | 久久久久久久午夜电影| 深夜精品福利| 国产真人三级小视频在线观看| 欧美成人免费av一区二区三区| 99国产综合亚洲精品| 久久精品aⅴ一区二区三区四区| 韩国精品一区二区三区| 少妇 在线观看| 国产精品久久久av美女十八| 国产亚洲精品久久久久久毛片| 午夜老司机福利片| 国产成人精品久久二区二区免费| 动漫黄色视频在线观看| 一级作爱视频免费观看| 此物有八面人人有两片| 午夜精品在线福利| 一进一出抽搐动态| 精品午夜福利视频在线观看一区| 日本撒尿小便嘘嘘汇集6| 日韩一卡2卡3卡4卡2021年| 欧美色视频一区免费| 亚洲一卡2卡3卡4卡5卡精品中文| 两个人视频免费观看高清| 老司机在亚洲福利影院| 国产激情久久老熟女| 亚洲成人国产一区在线观看| 久久欧美精品欧美久久欧美| 日韩精品青青久久久久久| 久久国产精品影院| 免费看日本二区| 一级作爱视频免费观看| 99国产精品99久久久久| 欧美又色又爽又黄视频| 一夜夜www| 国产精品亚洲av一区麻豆| 黑人操中国人逼视频| 嫩草影视91久久| av有码第一页| 欧美黄色淫秽网站| 757午夜福利合集在线观看| 午夜福利18| 免费高清视频大片| 男女之事视频高清在线观看| 99国产综合亚洲精品| 视频在线观看一区二区三区| 一本久久中文字幕| 精品国产乱子伦一区二区三区| x7x7x7水蜜桃| 欧美在线一区亚洲| 满18在线观看网站| av超薄肉色丝袜交足视频| 高清毛片免费观看视频网站| 很黄的视频免费| 久久人妻av系列| av电影中文网址| 天天躁夜夜躁狠狠躁躁| 国产精品影院久久| 亚洲美女黄片视频| 欧美国产精品va在线观看不卡| 9191精品国产免费久久| 欧美日韩瑟瑟在线播放| 欧美性长视频在线观看| 国产伦一二天堂av在线观看| 一边摸一边做爽爽视频免费| 国产黄片美女视频| 波多野结衣高清无吗| 在线观看www视频免费| 成人国产一区最新在线观看| 亚洲自拍偷在线| 亚洲中文字幕一区二区三区有码在线看 | 波多野结衣av一区二区av| 97碰自拍视频| 制服丝袜大香蕉在线| 国产亚洲欧美98| 可以在线观看毛片的网站| 国产高清视频在线播放一区| 亚洲国产看品久久| 久久久久国内视频| 久久人妻av系列| 欧美绝顶高潮抽搐喷水| 免费在线观看完整版高清| 久久久久九九精品影院| 国产久久久一区二区三区| 亚洲国产欧美网| 在线观看66精品国产| 精品国产亚洲在线| 亚洲国产欧洲综合997久久, | 亚洲色图 男人天堂 中文字幕| 国产日本99.免费观看| 给我免费播放毛片高清在线观看| 欧美中文日本在线观看视频| 一级毛片高清免费大全| 久久久久久免费高清国产稀缺| 中文字幕人成人乱码亚洲影| 国产野战对白在线观看| 亚洲中文字幕日韩| 国产黄片美女视频| 亚洲国产看品久久| 久久久久久国产a免费观看| 制服丝袜大香蕉在线| 亚洲成人国产一区在线观看| 大型黄色视频在线免费观看| 国产一区二区三区视频了| 精品国产亚洲在线| 欧美激情 高清一区二区三区| 一边摸一边抽搐一进一小说| 久久青草综合色| 韩国精品一区二区三区| 日本 欧美在线| 国产精品一区二区免费欧美| 亚洲人成网站在线播放欧美日韩| 村上凉子中文字幕在线| 成人免费观看视频高清| 国产熟女午夜一区二区三区| 一本精品99久久精品77| 手机成人av网站| 亚洲第一av免费看| 亚洲精品中文字幕一二三四区| 99久久国产精品久久久| 国内精品久久久久精免费| 搞女人的毛片| 国产成人欧美| 亚洲,欧美精品.| 18美女黄网站色大片免费观看| 三级毛片av免费| 精品国产乱子伦一区二区三区| 黄网站色视频无遮挡免费观看| www日本黄色视频网| 黄色视频,在线免费观看| 最近在线观看免费完整版| 欧美中文日本在线观看视频| 日本三级黄在线观看| 在线观看舔阴道视频| 国产99白浆流出| 男人舔女人下体高潮全视频| 色播亚洲综合网| 观看免费一级毛片| 波多野结衣av一区二区av| 夜夜爽天天搞| 中文字幕高清在线视频| 一a级毛片在线观看| 嫁个100分男人电影在线观看| 国产激情偷乱视频一区二区| 男女下面进入的视频免费午夜 | 久久久久免费精品人妻一区二区 | 国产区一区二久久| 一级黄色大片毛片| 免费高清在线观看日韩| 亚洲天堂国产精品一区在线| 自线自在国产av| 国产成人精品久久二区二区91| 天堂影院成人在线观看| 18禁观看日本| 精品久久久久久久末码| 精品国产国语对白av| 一区福利在线观看| 在线观看免费日韩欧美大片| 中文字幕精品亚洲无线码一区 | 伊人久久大香线蕉亚洲五| 国产精品九九99| 久久九九热精品免费| 欧美国产日韩亚洲一区| 亚洲精品国产区一区二| 亚洲人成77777在线视频| 成熟少妇高潮喷水视频| 欧美绝顶高潮抽搐喷水| 国产亚洲精品第一综合不卡| 免费高清视频大片| 国产成人影院久久av| 国产成人欧美在线观看| 一区二区三区国产精品乱码| 亚洲精品一卡2卡三卡4卡5卡| 国产精品一区二区三区四区久久 | 中文亚洲av片在线观看爽| 亚洲精品美女久久av网站| 成年免费大片在线观看| 日本黄色视频三级网站网址| 51午夜福利影视在线观看| 精品国产乱子伦一区二区三区| 美女大奶头视频| 又大又爽又粗| 老熟妇乱子伦视频在线观看| 99国产综合亚洲精品| 两个人看的免费小视频| 波多野结衣av一区二区av| 成人三级黄色视频| 午夜免费激情av| 99国产精品99久久久久| 国产私拍福利视频在线观看| 色av中文字幕| 免费看a级黄色片| a在线观看视频网站| 91九色精品人成在线观看| 国产精品免费一区二区三区在线| 成人亚洲精品av一区二区| 黄色 视频免费看| 免费女性裸体啪啪无遮挡网站| 国产精品一区二区免费欧美| 国产真人三级小视频在线观看| 久久久国产成人免费| 老熟妇乱子伦视频在线观看| 91成年电影在线观看| 丰满人妻熟妇乱又伦精品不卡| 神马国产精品三级电影在线观看 | 欧美久久黑人一区二区| 久久精品国产亚洲av香蕉五月| 日韩一卡2卡3卡4卡2021年| 熟女电影av网| 99久久无色码亚洲精品果冻| 日日夜夜操网爽| 黄色女人牲交| 特大巨黑吊av在线直播 | 十分钟在线观看高清视频www| 亚洲精品在线美女| 久久精品成人免费网站| 美女国产高潮福利片在线看| 少妇裸体淫交视频免费看高清 | 久久精品亚洲精品国产色婷小说| 午夜福利免费观看在线| 免费在线观看日本一区| 精品国产一区二区三区四区第35| 久久精品成人免费网站| 香蕉久久夜色| 欧美日韩中文字幕国产精品一区二区三区| 99精品久久久久人妻精品| 一本综合久久免费| 淫秽高清视频在线观看| 中文在线观看免费www的网站 | www日本在线高清视频| 国产一级毛片七仙女欲春2 | 国产成人av教育| 精品国产国语对白av| 2021天堂中文幕一二区在线观 | 国产精品 欧美亚洲| 色哟哟哟哟哟哟| 亚洲中文av在线| 免费看日本二区| 两性午夜刺激爽爽歪歪视频在线观看 | 精品第一国产精品| 欧美三级亚洲精品| 母亲3免费完整高清在线观看| 观看免费一级毛片| 免费无遮挡裸体视频| 亚洲性夜色夜夜综合| 国产97色在线日韩免费| 欧美三级亚洲精品| 日韩欧美免费精品| 9191精品国产免费久久| 国产真人三级小视频在线观看| 1024视频免费在线观看| 精品国内亚洲2022精品成人| 欧美又色又爽又黄视频| 亚洲精华国产精华精| 99国产精品一区二区三区| 黄色女人牲交| 黄色片一级片一级黄色片| 熟女电影av网| 欧美日韩中文字幕国产精品一区二区三区| 亚洲国产中文字幕在线视频| 亚洲第一av免费看| 别揉我奶头~嗯~啊~动态视频| 亚洲成国产人片在线观看| 黄网站色视频无遮挡免费观看| 日韩av在线大香蕉| e午夜精品久久久久久久| 动漫黄色视频在线观看| 十八禁人妻一区二区| 国产成人精品久久二区二区免费| 视频区欧美日本亚洲| 热99re8久久精品国产| 中文资源天堂在线| 黄色 视频免费看| 啪啪无遮挡十八禁网站| 国产亚洲精品av在线| 国产蜜桃级精品一区二区三区| 精品无人区乱码1区二区| 亚洲精品一卡2卡三卡4卡5卡| 男人舔奶头视频| 午夜精品久久久久久毛片777| 国产激情偷乱视频一区二区| 成人免费观看视频高清| 狂野欧美激情性xxxx| 好男人电影高清在线观看| 亚洲欧美一区二区三区黑人| 黄色毛片三级朝国网站| 精品国产一区二区三区四区第35| 麻豆成人午夜福利视频| 欧美国产日韩亚洲一区| 成人av一区二区三区在线看| 亚洲人成77777在线视频| 国产av一区在线观看免费| 欧美在线黄色| 91麻豆av在线| 欧美亚洲日本最大视频资源| 极品教师在线免费播放| 精品久久久久久久久久免费视频| 宅男免费午夜| 久久久久国产一级毛片高清牌| 此物有八面人人有两片| 欧美成人性av电影在线观看| 中文字幕精品免费在线观看视频| 欧美黑人欧美精品刺激| 伦理电影免费视频| 色在线成人网| 中出人妻视频一区二区| 久久婷婷成人综合色麻豆| 禁无遮挡网站| 欧美国产精品va在线观看不卡| 欧美日韩亚洲综合一区二区三区_| 亚洲天堂国产精品一区在线| 国产免费av片在线观看野外av| 欧美日韩一级在线毛片| 欧美精品亚洲一区二区| 99国产综合亚洲精品| 一本精品99久久精品77| 男女午夜视频在线观看| 后天国语完整版免费观看| 好看av亚洲va欧美ⅴa在| 天天躁夜夜躁狠狠躁躁| 中国美女看黄片| 精品国内亚洲2022精品成人| 女性生殖器流出的白浆| 最近最新中文字幕大全电影3 | 制服人妻中文乱码| 国产av一区二区精品久久| 满18在线观看网站| 国产熟女xx| 精品乱码久久久久久99久播| 黄片播放在线免费| 99re在线观看精品视频| 男人操女人黄网站| 国产一区二区三区视频了| 亚洲精品中文字幕一二三四区| 亚洲免费av在线视频| 久久人人精品亚洲av| 中国美女看黄片| 欧美乱色亚洲激情| 国产精品九九99| 99久久99久久久精品蜜桃| 美女免费视频网站| 麻豆成人午夜福利视频| 国产高清视频在线播放一区| 日本 欧美在线| 真人做人爱边吃奶动态| 亚洲av电影在线进入| 日本在线视频免费播放| 最好的美女福利视频网| 欧美中文日本在线观看视频| 18禁观看日本| 日本一区二区免费在线视频| 男人舔女人的私密视频| 日韩视频一区二区在线观看| www.精华液| 久久精品亚洲精品国产色婷小说| 亚洲精品一区av在线观看| 天天添夜夜摸| 人人澡人人妻人| 99久久久亚洲精品蜜臀av| 亚洲av熟女| 日本在线视频免费播放| 国产亚洲精品久久久久久毛片| 欧美日韩亚洲综合一区二区三区_| 欧美大码av| 午夜福利在线在线| 欧美色视频一区免费| 国产成人欧美| 久久性视频一级片| 亚洲熟女毛片儿| 性色av乱码一区二区三区2| 女性被躁到高潮视频| 97人妻精品一区二区三区麻豆 | 色av中文字幕| 欧美激情久久久久久爽电影| 国产99白浆流出| 久久久国产精品麻豆| 一进一出抽搐动态| 一级毛片女人18水好多| 久久国产乱子伦精品免费另类| 香蕉丝袜av| 少妇的丰满在线观看| 国产成人欧美| 99在线人妻在线中文字幕| 97超级碰碰碰精品色视频在线观看| 波多野结衣av一区二区av| 亚洲成人久久爱视频| av中文乱码字幕在线| 国产精华一区二区三区| 每晚都被弄得嗷嗷叫到高潮| 亚洲精品久久国产高清桃花| 午夜视频精品福利| 18禁裸乳无遮挡免费网站照片 | 日韩精品中文字幕看吧| 中国美女看黄片| 欧美性长视频在线观看| 亚洲精华国产精华精| 香蕉久久夜色| a级毛片在线看网站| 国产伦在线观看视频一区| 日韩欧美三级三区| 国产久久久一区二区三区| 无遮挡黄片免费观看| 欧美人与性动交α欧美精品济南到| 999久久久国产精品视频| 十八禁人妻一区二区| 亚洲狠狠婷婷综合久久图片| 可以免费在线观看a视频的电影网站| 伊人久久大香线蕉亚洲五| 又大又爽又粗| 午夜视频精品福利| 国产精华一区二区三区| 国产精品亚洲一级av第二区| 中文字幕人妻熟女乱码| 又紧又爽又黄一区二区| 欧美性猛交╳xxx乱大交人| 国产精品国产高清国产av| 日韩视频一区二区在线观看| 搡老熟女国产l中国老女人| 久久久精品欧美日韩精品| xxx96com| 丰满的人妻完整版| 国产久久久一区二区三区| 亚洲成av片中文字幕在线观看| 后天国语完整版免费观看| 熟女电影av网| 国产精品 国内视频| 欧美三级亚洲精品| 又大又爽又粗| 欧美精品啪啪一区二区三区| 日韩大码丰满熟妇| 在线观看66精品国产| 香蕉国产在线看| 国产乱人伦免费视频| 露出奶头的视频| 精品高清国产在线一区| 免费女性裸体啪啪无遮挡网站| netflix在线观看网站| 69av精品久久久久久| 国产伦在线观看视频一区| 最近最新免费中文字幕在线| 又黄又爽又免费观看的视频| 制服诱惑二区| 一级黄色大片毛片| 99久久国产精品久久久| 国产一区在线观看成人免费| 国产v大片淫在线免费观看| 国产免费av片在线观看野外av| 精品欧美一区二区三区在线| 国产私拍福利视频在线观看| 亚洲人成电影免费在线| 99re在线观看精品视频| 中文在线观看免费www的网站 | 久久精品影院6| 91在线观看av| 亚洲精品国产精品久久久不卡| 国产亚洲精品综合一区在线观看 | 88av欧美| 免费一级毛片在线播放高清视频| 最近最新中文字幕大全免费视频| 搞女人的毛片| 国产色视频综合| 在线观看66精品国产| 国产极品粉嫩免费观看在线| 欧美成人午夜精品| 99在线视频只有这里精品首页| 精品久久久久久,| 亚洲一区高清亚洲精品| 欧美激情 高清一区二区三区| 熟女少妇亚洲综合色aaa.| 久久精品91蜜桃| 免费在线观看亚洲国产| 亚洲真实伦在线观看| 日韩欧美在线二视频| 97人妻精品一区二区三区麻豆 | 午夜两性在线视频| 免费看十八禁软件| 精品人妻1区二区| 热99re8久久精品国产| 久久精品夜夜夜夜夜久久蜜豆 | 在线观看日韩欧美| 欧美色视频一区免费| 妹子高潮喷水视频| 日韩有码中文字幕| 伊人久久大香线蕉亚洲五| 国产激情久久老熟女| av中文乱码字幕在线| 亚洲第一欧美日韩一区二区三区| 男女做爰动态图高潮gif福利片| 夜夜躁狠狠躁天天躁| 亚洲精品久久成人aⅴ小说| 日日爽夜夜爽网站| 久久亚洲真实| 国产精品亚洲美女久久久| 桃红色精品国产亚洲av| 村上凉子中文字幕在线| 久久久久国产一级毛片高清牌| 国产精品一区二区三区四区久久 | 很黄的视频免费| www.自偷自拍.com| 伦理电影免费视频| 久久久国产精品麻豆| 真人做人爱边吃奶动态| 成在线人永久免费视频| 精品久久久久久久久久久久久 | 男女做爰动态图高潮gif福利片| 女警被强在线播放| 不卡av一区二区三区| 亚洲熟妇中文字幕五十中出| 最新美女视频免费是黄的| 亚洲国产看品久久| 亚洲avbb在线观看| 人人妻人人澡欧美一区二区| 两人在一起打扑克的视频| 日本黄色视频三级网站网址| 亚洲成av片中文字幕在线观看| 日韩大尺度精品在线看网址| 男人舔女人下体高潮全视频| 欧洲精品卡2卡3卡4卡5卡区| 村上凉子中文字幕在线| 亚洲第一欧美日韩一区二区三区| 久久精品人妻少妇| 桃色一区二区三区在线观看| 国产精品综合久久久久久久免费| 99久久国产精品久久久| 变态另类成人亚洲欧美熟女| 国产欧美日韩精品亚洲av| 中文字幕久久专区| avwww免费| 国产区一区二久久| 欧美日本亚洲视频在线播放| 欧美日韩一级在线毛片| 国产精品98久久久久久宅男小说| 人人妻,人人澡人人爽秒播| 日本免费a在线| 中文字幕最新亚洲高清| 香蕉av资源在线| 成人国语在线视频| 国产成人系列免费观看| 精品一区二区三区四区五区乱码| 婷婷精品国产亚洲av在线| 午夜精品在线福利| netflix在线观看网站| 91av网站免费观看| 精品第一国产精品| 女生性感内裤真人,穿戴方法视频| 成年免费大片在线观看| 巨乳人妻的诱惑在线观看| 在线观看日韩欧美| 久久久久亚洲av毛片大全| 嫁个100分男人电影在线观看| 久久婷婷成人综合色麻豆| 俺也久久电影网| 伦理电影免费视频| 嫩草影视91久久| 成熟少妇高潮喷水视频| 亚洲激情在线av| 国内久久婷婷六月综合欲色啪| 亚洲av美国av| 日韩欧美国产在线观看| 啦啦啦 在线观看视频| 观看免费一级毛片| 男女之事视频高清在线观看| 色老头精品视频在线观看| x7x7x7水蜜桃| 国产真人三级小视频在线观看| 中文字幕高清在线视频| 欧美日韩亚洲综合一区二区三区_| 免费在线观看影片大全网站| 午夜精品久久久久久毛片777| 亚洲五月色婷婷综合| 久9热在线精品视频| 国产乱人伦免费视频| 欧美色欧美亚洲另类二区| 国产成人系列免费观看| 久久青草综合色|