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

    Unveiling the Underlying Mechanism of Transition Metal Atoms Anchored Square Tetracyanoquinodimethane Monolayers as Electrocatalysts for N2 Fixation

    2022-07-04 09:12:46ShengyaoLvChunxiangHuangGuoliangLiandLimingYang
    Energy & Environmental Materials 2022年2期

    Shengyao Lv , Chunxiang Huang, Guoliang Li*, and Liming Yang*

    1. Introduction

    Conversion of dinitrogen(N2)to ammonia(NH3)is one of the important processes from the perspectives of agriculture, industry, and human development. NH3is the critical ingredient for the fertilizer production and the synthesis of nucleotides, amino acids, and some other main organic compounds.[1–3]NH3is also a potential energy carrier and transportation fuel.[4,5]In the past century, NH3has always been synthesized by N2and H2using the industrial Haber-Bosch(HB)approach.[6]However, N2is difficult to be utilized due to its extreme stability and large amount of energy will be required to destroy its inert triple N≡N bond. Traditional Haber-Bosch process operates at high temperature and pressure and consumes the world energy usage of approximately 2%.[7,8]Furthermore, the HB process emits a large amount of CO2and causes greenhouse gas effects.[9]Thus, it is highly desirable to design and develop new routes for N2reduction to NH3at mild conditions.

    Electrocatalytic N2fixation is an environmentfriendly process and has received widespread concern.Green,efficient,and sustainable electrocatalytic N2reduction reaction(NRR)provides a bright prospect for N2fixation and NH3synthesis,because it can reduce the energy for breaking the N≡N bond to a certain extent.[10–12]In the electrochemical NRR process, the NH3production efficiency can be improved by choosing effective catalyst,controlling the operating potential,pH,and electrolyte.[13,14]There is no doubt that highly efficient catalysts, which can increase the reaction rate,enhance the selectivity,and reduce energy consumption,are the core of electrocatalytic NRR process.[15–26]

    2D materials are an important platform to design novel NRR catalysts.[21,24–28]As a unique ligand, 7,7,8,8-tetracyanoquinodimethane(TCNQ) has a high electron affinity of 3.383 eV,[29]which makes it a strong organic electron acceptor and an ideal component for the construction and synthesis of metal-organic semiconductor charge-transfer materials.Indeed,a large number of TCNQ-based complexes have been synthesized. Recent experimentally fabricated 2D square Ni- and Mn-TCNQ materials and their intriguing physicochemical properties on Au(111)and Ag(100)surfaces attracted our attention.[30,31]Based on our experience, 2D square TM-TCNQ monolayers (here labeled as TMsTCNQ) might be potential single-atom catalysts (SACs) for electrochemical NRR. Indeed, some TM-sTCNQ materials have been reported for electrochemical oxygen evolution reaction (OER).[32–34]Inspired by these intriguing properties and promising progress and applications of the TM-TCNQ materials, we carried out a comprehensive investigation on the electrochemical mechanism of NRR on 3d, 4d, 5d TM atoms embedded in square TCNQ networks by high-throughput screenings and first-principals calculations. Considering the computational costs and the purpose of finding a promising NRR catalyst, our research follows the following steps. First, the TM-sTCNQ monolayer was subjected to preliminary high-throughput screening, and the TMsTCNQ SAC with good electrocatalytic NRR performance was selected.Then,comprehensive studies were carried out on the whole electrocatalytic NRR process,the geometries,stabilities,electronic properties,and magnetic moments for the selected TM-sTCNQ monolayers.

    2. Computational Details

    All the computations were performed with Vienna Ab Initio Simulation Package (VASP)[35]and the spin-polarized density functional theory (DFT) method was used. The projector-augmented wave(PAW)[36]pseudopotentials were to employed describe the ionelectron interactions. The exchange-correlation functional of Perdew?Burke?Ernzerhof (PBE) with generalized gradient approximation (GGA) was utilized to describe the electron-electron interactions.[37]Considering the influence of possible dispersion interactions on Gibbs free energy, the dispersion corrected DFT-D3 method was employed.[38]In order to ensure the computational accuracy, the kinetic cutoff energy of 450 eV was set for the plane wave function, and the convergence threshold for the energy and force in structural optimization process was 10?5eV and 0.02 eV ?A?1, respectively. To prevent the interaction between adjacent periodic atoms, a vacuum height of over 15 ?A was used on the z axis. For structural optimizations and electronic structure calculations, the Brillouin area used Γ-centered Monkhorst?Pack (MP)k-point grids of 3 × 3 × 1 and 7 × 7 × 1, respectively. To further evaluate the thermodynamic stability of TM-sTCNQ, ab initio molecular dynamics (AIMD) simulations were performed at a temperature of 500 K for a 10 ps time period with the step being 2 fs. A large number of relevant experiences have shown that the effect of solvation on NRR is ~0.1 eV.[39]Due to its large computational costs, the solvation effects have not been considered in the present research.[24,40,41]

    A single TM atom was anchored to the crystal surface of the sTCNQ film to simulate SAC, and the binding energy (Eb) between TM and substrate can be obtained by Equation (1):

    with ETM-sTCNQ, ETM, and EsTCNQbeing the energies of the TMsTCNQ unit, the TM atom, and the sTCNQ unit, respectively.

    The cohesive energy (Ec) of the TM atoms was calculated using Equation (2):

    with ETM(bulk), ETM(single), and n being the energy of a bulk metal,a TM atom in vacuum, and the number of TM atoms in the bulk section, respectively.

    The Gibbs free energy change(ΔG)at each elementary step was calculated by Equation (3) based on the computational hydrogen electrode(CHE)model:[42]

    with ΔE being the energy difference between the product and reactant;∫CpdT being the enthalpy correction (enthalpy correction improves the accuracy of computational results by 0.03–0.11 eV); EZPEand S being the zeropoint energies and entropies,respectively,given by vibrational frequency computations;T = 298.15 K herein.When calculating the frequencies, we fixed the TM-sTCNQ catalytic material, and only the adsorbed chemical species were allowed to vibrate.The EZPE,Cp,and S values for small gas molecules (H2, N2, NH3)come from the NIST database.[43]U is the required applied potential. ΔGpHis the correction of H+free energy with concentration,ΔGpH= 2.303 × kB× T × pH,in which kBis the Boltzmann constant, and pH is set to 0 in the present research.

    The overpotential (η) equals to the difference between Uequilibrium(the equilibrium potential) and Uonset(the onset potential),η = Uequilibrium?Uonset, where Uonset= ?ΔGmax/e and Uequilibrium= 0.17 V for NRR.

    Figure 1. a) The structure of the TM-sTCNQ monolayer, with the unit cell of TM-sTCNQ represented by the blue dotted square. N2 is adsorbed on TM-sTCNQ by end-on and side-on forms, as shown by the pink dashed cycles above. Green, brown, light blue, and light pink balls indicate the TM, C, N, and H atoms,respectively.Schematic diagram of adsorbed N2 bonding to the TM atom on TM-sTCNQ by b)end-on and c)side-on forms.

    3. Results and Discussion

    Before starting the electrocatalyst screening and mechanism investigation, we first analyze the topology motif and structures of 2D square TM atoms embedded TCNQ networks (i.e., TM-sTCNQ). Following the structures reported experimentally,[30,31]we construct the TM-sTCNQ monolayer via Materials Studio. The unit cell of 2D TM-sTCNQ monolayer (the blue dash square in the left part of Figure 1a) includes two TM atoms and two TCNQ molecules with the chemical composition as(C12H4N4TM)2and the nominal TM:TCNQ stoichiometry is 1:1.The central TM atom is coordinated by four TCNQ in a “windmill” form and extended to the 2D sheet with the same coordination pattern. All TM atoms are equivalent and all TCNQ molecules are also equivalent in 2D TM-sTCNQ monolayers. For ease of discussion and analysis, the TM atoms in the unit cell can be “formally” distinguished by the unit cell frame as the central TM atom and the four TM atoms on the frame border even though they are equivalent. The central TM atom is solely belonged to the unit cell,whereas the four TM atoms on the frame border are equally shared by four adjacent unit cells, thus contribute only one(4 × 1/4 = 1) TM atom to the unit cell. Each TCNQ molecule is shared by two adjacent unit cells, thus, totally, two (4 × 1/2 = 2)TCNQ molecules belong to the unit cell.

    Now, we turn to discuss the N2adsorption on the TM-sTCNQ monolayer.Figure 1a shows the two adsorption modes,end-on and side-on, and Figure 1b,c give the orbital interactions between adsorbed N2and the active center TM atom of electrocatalysts TMsTCNQ. The empty TM-d orbitals can accept electrons from N2. At the same time, the filled TM-d orbital can provide electrons into the π*antibonding orbital of N2.Affected by the metal-ion type,coordination number, and steric effect, when the TM-sTCNQ monolayer adsorbs N2in side-on form, there may also be interactions between the π-bonding orbital of the N2molecule and the TM-d orbital(Figure 1c).This“acceptance-donation”interplay leads to the activation of adsorbed N2molecules.

    All the optimized structures of the TM-sTCNQ monolayers are planar, except that the Au atoms in Au-sTCNQ are slightly out of the plane,as shown in Table S1. The Au atom and the TCNQ molecule in Au-sTCNQ have a two-coordinate structure, in agreement with related experimental results.[44]

    3.1. Initial Selection of the TM-sTCNQ Electrocatalysts

    Electrocatalytic NRR usually refers to a six proton-coupled electron transfer process (N2+ 6H++ 6e?→2NH3). The N2adsorption, the first protonation (*N2+ H++ e?→*N2H), and the last protonation(*NH2+ H++ e?→*NH3)steps are the three key steps for electrocatalytic NRR.The N2adsorption step determines whether NRR can launch spontaneously.

    Figure 2. Gibbs free energy change ΔG*N2 and N≡N bond distance RN–N for N2 adsorbed on various TM-sTCNQ monolayers in a) end-on and b) side-on forms. The black line represents ΔG*N2, and the blue line represents the RN–N. Calculated values of ΔG*N2 and the RN–N in N2-TM-sTCNQ are listed in Table S3.

    Figure 3. Gibbs free energy change ΔG*N2→*N2H and ΔG*NH2→*NH3 for N2 adsorbed on different TMsTCNQ monolayers in a) end-on and b) side-on forms. The green columns represent ΔG*N2→*N2H, and the red columns represent ΔG*NH2→*NH3. The ΔG*N2→*N2H and ΔG*NH2→*NH3 values are listed in Table S4.

    3.1.1. N2Adsorption

    The prerequisite for the occurrence of NRR on the TM-sTCNQ surface is that the TM-sTCNQ monolayers can effectively capture and activate N2molecules, N2+ * →*N2(* means the TM-sTCNQ monolayer). The free energy change for N2adsorption (ΔG*N2) and the extension of the N≡N bond distance (RN–N)are the basic standards to measure the catalytic ability of a material to the NRR process.In this article, we consider two common N2adsorption modes (Figure 1a and Figure S1): 1)end-on mode, in which N2is placed upright relative to the TM-sTCNQ monolayer with only one TM–N bond formed; 2) side-on mode,in which N2is parallelly placed on the TM-sTCNQ monolayer with both N atoms bonded to the TM active center. Considering 3d, 4d, 5d elements, 30 kinds of TM atoms,and two N2adsorption configurations, we have constructed 60 kinds of N2-TM-sTCNQ structures. Figure 2 shows the N2adsorption energies and the RN–Nchanges after N2is adsorbed on different TM-sTCNQ monolayers.Some attempted N2-TM-sTCNQ structures desorb N2during the structural optimizations(Table S2),indicating that they cannot adsorb N2spontaneously, and will not be discussed further. As shown in Figure 2, when TM =Sc,Ti,V,Cu,Y,Zr,Nb, Mo, Tc,Lu,Hf,Ta,W,and Re,and N2adsorbed in end-on form,ΔG*N2< 0; when TM = Sc, Y, Zr, Nb, Lu,Hf, Ta, and W, and N2adsorbed in side-on form, ΔG*N2< 0. ΔG*N2< 0 indicates the TM-sTCNQ monolayer can adsorb N2molecules spontaneously. Compared with the side-on modes, we find that TM-sTCNQ is more inclined to adsorb N2molecules with end-on modes because of its smaller steric hindrance.Compared with the bond distance of 1.115 ?A in free N2molecule,the RN–Nin adsorbed N2increases,indicating that N2has been activated effectively, which is beneficial to subsequent N≡N cleavage and further NRR process.

    3.1.2. The First and Last Protonation Steps of NRR

    The first and last protonation steps are the two steps which may have the largest Gibbs free energy change(ΔGmax)among the NRR pathway and thus be the potential determining step (PDS), thus become the focus of this study.Based on the results above,we continue to explore the TM-sTCNQ monolayers with ΔG*N2< 0, as shown in Figure 2.Although there are many TM-sTCNQ monolayers that can capture and activate N2molecules,most of them cannot form*N2H or make*NH2stable. Using ΔG < 0.8 eV as a criterion,only the Nb-,Mo-,Tc-,W-,and Re-sTCNQ monolayers (N2adsorption in end-on form) and the Sc-, Nb-, and W-sTCNQ monolayers (N2adsorption in side-on form)show good NRR catalytic activity(Figure 3).

    3.1.3. Hydrogen Evolution Reaction

    The hydrogen evolution reaction (HER) is the main competitive reaction of electrocatalytic NRR.[39]The selectivity of the catalyst directly affects their catalytic performance for NRR. If protons preferentially occupy the TM active site of the TM-sTCNQ monolayers during the NRR process,the TM-sTCNQ monolayers cannot effectively capture N2molecules,and the NRR catalytic efficiency will be greatly reduced.So,good selectivity is very important for evaluating the performance of a catalyst. The competition between NRR and HER on the TM-sTCNQ monolayers can be compared by calculating ΔG*N2and ΔG*H, with ΔG*N2< ΔG*H(ΔG*N2is more negative) indicating that N2is more easily adsorbed on the TM-sTCNQ surface. As shown in Figure 4, for Nb-and Mo-sTCNQ,ΔG*N2is more negative than ΔG*H,so they have good selectivity for NRR. Thus, the Nb- and Mo-sTCNQ monolayers may be promising electrocatalysts for NRR.For Tc-sTCNQ,ΔG*N2and ΔG*Hare comparable though ΔG*His slightly more negative, and we will continue to study its electrocatalytic performance in depth.

    Figure 4. Gibbs free energy change ΔG*N2 and ΔG*H on different TMsTCNQ monolayers. According to the high-throughput screening results above, the N2 adsorption in end-on configuration is considered on the Nb-,Mo-, Tc-, W-, and Re-sTCNQ monolayers, and that in side-on configuration is considered on the Sc-, Nb-, and W-sTCNQ monolayers. The red, fuchsia,and blue colors represent ΔG*N2 (side-on), ΔG*N2 (end-on), and ΔG*H,respectively. The corresponding ΔG values are listed in Table S5.

    Figure 5. A simplified diagram for the possible reaction paths of nitrogen reduction reaction (NRR).The distal, alternating, consecutive, enzymatic, and mixed pathways are shown in green, blue, dark purple, red, and light purple, respectively. For simplicity, the NRR path only presents the bonding relationship between adsorbed species and the active center TM atom.

    3.2. Catalytic Mechanism for NRR on the TM-sTCNQ Electrocatalysts

    During NRR, the N2molecules absorbed on electrocatalysts will be hydrogenated gradually by proton-coupled electron transfers, forming NRR intermediates adsorbed on the TM site,until each N2molecule is converted into two NH3molecules. Generally, the NRR process has four types of reaction pathways (Figure 5).When N2is adsorbed in end-on configuration,it has a distal or alternating pathway; when N2is adsorbed in side-on form, it has consecutive or enzymatic pathway. The distal and consecutive pathways are similar, in which the protonelectron pairs first attack one of the N atoms continuously,releasing one NH3,then attack the other N atom to generate the second NH3. On the contrary, in the alternating and enzymatic pathways, the proton-electron pairs alternately attack two N atoms of N2,and successively generate two NH3on the last two protonation steps.In addition to four commonly known pure reaction pathways above, we also consider the mixed pathways so as to comprehensively looking for the low onset potential paths. As shown in Figure 2, N2can be adsorbed spontaneously(ΔG*N2< 0) on Nb-sTCNQ in both side-on and end-on forms, while on Mo-sTCNQ and Tc-sTCNQ only in end-on form.Thus,in the following computations,only the NRR pathways following these N2-TM-sTCNQ intermediates are explored.

    Figure 6 shows Gibbs free energy diagrams for NRR on the NbsTCNQ monolayer through various possible pathways.Through the distal pathway, for the N2capture and the first protonation steps, the ΔG*N2and ΔG*N2Hare ?0.41 and 0.45 eV,respectively,and the N–N bond lengths are extended to 1.144 and 1.244 ?A, respectively. The*N2H formation needs to overcome a small barrier.The next four protonation steps of NRR are all exothermic, with ΔG*N-NH2, ΔG*N,ΔG*NH, and ΔG*NH2being ?0.54, ?0.47, ?0.96, and ?0.05 eV,respectively. The last protonation step(*NH2+ H++ e?→*NH3) is the PDS, with the ΔGmaxbeing 0.48 eV. The release of NH3caused a slight upward slope(0.65 eV)of free energy.Thus,the NH3generated can be removed in a flowing N2environment,and the NRR recycle can be easily realized.

    Figure 6. Gibbs free energy diagrams for nitrogen reduction reaction on Nb-sTCNQ through a) distal, b)alternating, c) distal-alternating, d) consecutive, and e) enzymatic pathways, at both zero and onset potentials. The black indicates the zero potential diagram, while the red shows the onset potential diagram. The two values involved in the PDS are highlighted in blue.

    As shown in Figure 6b, in the alternating pathway, the N2adsorption and *N2H formation are the same as those in the distal pathway.The Gibbs free energy of the second protonation step still rises(ΔG*NH–NH= 0.50 eV), and the N–N bond length is elongated to 1.302 ?A. For the *NH–NH + H++ e?→*NH–NH2and *NH–NH2+ H++ e?→*NH2–NH2steps, the free energy changes are?0.73 and 0.57 eV, respectively, and the N–N bond lengths are extended to 1.425 and 1.471 ?A, respectively. The remaining steps are the same as the distal pathway. The PDS of the alternating pathway is*NH–NH2+ H++ e?→*NH2–NH2,with ΔGmaxbeing 0.57 eV.

    Various possible distal-alternating mixed paths have also been considered. Figure 6c shows the optimal mixed path, where PDS is the*NH2+ H++ e?→*NH3step with onset potential equal to?0.48 V.In general,the Nb-sTCNQ still prefers to follow a distal path for NRR,when N2is adsorbed in end-on configuration.

    In the consecutive and enzymatic pathways, N2is adsorbed on the Nb-sTCNQ monolayer in side-on form. In these two pathways, only the*N2→*N–NH*and *NH2→*NH3steps are endothermic,with ΔG*N–NH*and ΔG*NH3being 0.02 and 0.48 eV,respectively.The PDSs for these two pathways are both *NH2+ H++ e?→*NH3, with ΔGmax= 0.48 eV.Although the possibility of capturing N2on the NbsTCNQ monolayer in side-on form is lower (ΔG*N–N*= ?0.07 eV)than that in end-on form,the consecutive and enzymatic pathways are very smooth for the NRR process. Thus, the Nb-sTCNQ monolayer is promising for the electrocatalytic NRR process.

    The Mo-sTCNQ monolayer can also spontaneously adsorb N2(in end-on form), with ΔG*N2=?0.43 eV and RN–N=1.141 ?A. As shown in Figure S2, the formation of *N–NH needs to overcome a barrier of 0.50 eV,and RN–Nfurther extends to 1.231 ?A.In the distal pathway,the*N–NH →*N–NH2→*N →*NH →*NH2steps are all exothermic, with the free energy changes being ?0.28, ?0.84, ?0.40,and ?0.06 eV, respectively. The PDS is the*N–N + H++ e?→*N–NH step, with ΔGmax= 0.50 eV. Similarly, we explored the distal-alternating pathway on the MosTCNQ monolayer. The onset potential of the distal-alternating pathway is reduced to?0.50 V relative to ?0.55 V of the alternating pathway. On the Mo-sTCNQ monolayer, the NH3desorption energy in the last step is only 0.38 eV, so NH3can be released quickly from the catalyst surface under onset potential.

    For NRR on Tc-sTCNQ(Figure S3),the PDS is the *N2+ H++ e?→*N–NH step with Uonset= ?0.77 V through distal, alternating and mixed pathways. The reaction process is similar to that on the Mo-sTCNQ monolayer.

    We also analyzed the N–N distance changes along the NRR process on the Nb-,Mo-,and TcsTCNQ monolayers through different reaction pathways,with the results shown in Figures S4–S6.As NRR progresses,N–N gradually extends,vividly reflecting the process of gradual activation and reduction of N2.

    The PDS determines the minimum applied potential (onset potential) required by NRR,Uonset= ?ΔGmax/e. The smaller the negative onset potential, the smaller the overpotential and the barrier that needs to be overcome inthe reaction, which is more conducive to NRR. In this regard, after systematically exploring whole NRR processes on the Nb-, Mo-, and Tc-sTCNQ monolayers, the NRR PDSs through different reaction pathways on different catalysts are summarized and compared. As shown in Table 1, the NRR PDSs on the Nb-, Mo-, and Tc-sTCNQ monolayers are basically the first protonation step and the last protonation step, similar to previous research results. Thus, the first and last protonation steps of NRR are very important.

    Table 1. The PDS and onset potential of NRR on the Nb-, Mo-, and TcsTCNQ monolayers through different reaction paths. The first protonation step (*N2 + H+ + e?→*N–NH) is marked in red, the last protonation step (*NH2 + H+ + e?→*NH3) is marked in cyan, and the remaining reaction steps are shown in black.

    3.3. Origin of Catalytic Performance of TM-sTCNQ

    To understand the relationship between the NRR catalytic performance of the TM-sTCNQ monolayers and their magnetic moments,the change of the NRR onset potentials with relation to the magnetic moments of TM-sTCNQ (TM = Nb, Mo, Tc) was shown in Figure S7. Similar trend was found between the onset potentials and the TM-sTCNQ magnetic moments. Strong interactions between the TM atoms and the sTCNQ monolayers make the TM-sTCNQ(TM = Nb, Mo, Tc) materials have strong magnetic moments,which promotes the charge transfer between the NRR intermediates and the substrate, thereby reducing the onset potential of NRR.We also found that there is a similar changing trend between the onset potentials and the d band centers (εd) of the TM active center of the TM-sTCNQ catalysts, which is consistent with the influence of the magnetic moment of the catalyst on NRR. The TMsTCNQ catalyst can activate N2and reduce the barrier of the NRR process. As the εdof the Tc, Mo, and Nb active center increases,being ?0.61, 0.09, and 0.81 eV, respectively, the energy of the anti-bond orbital produced by the interaction between the active center and the adsorbed molecule will also increase.

    Figure 7. Charge variations of the three moieties along the a) distal, b) alternating, c) distal-alternating, d) consecutive, and e) enzymatic pathways on the Nb-sTCNQ catalyst. f) Definition of the three moieties of the NxHy-TM-sTCNQ intermediates: moiety 1 is the TM-sTCNQ substrate with its TMN4 excluded,moiety 2 is TMN4, and moiety 3 is the adsorbed NxHy species. Bader charge analysis diagrams for NRR on Mo-sTCNQ and Tc-sTCNQ are given in Figure S8,and corresponding values are shown in Tables S7–S9.

    Figure 8. The spin densities (in top and side view) for N2 adsorbed on, a, a′) Nb-sTCNQ (end-on), b, b′) Nb-sTCNQ (side-on), c, c′) Mo-sTCNQ (end-on)and d, d′) Tc-sTCNQ (end-on). The charge density differences (in top and side view) for N2 adsorbed on e, e′) Nb-sTCNQ (end-on), f, f′) Nb-sTCNQ (sideon), g, g′) Mo-sTCNQ (end-on) and h, h′) Tc-sTCNQ (end-on). The electron density accumulation area and depletion are shown in yellow and blue,respectively.

    The charge transfers of N2on TM-sTCNQ (TM = Nb, Mo, Tc)in the NRR process were explored by Bader charge analysis. Each reaction intermediate is divided into three parts,[24]the TM-sTCNQ substrate without TMN4(TM = Nb, Mo, Tc; N4= four surrounding N atoms) (moiety 1), TMN4(moiety 2), and adsorbed NxHyspecies (moiety 3), as presented in Figure 7f. Taking Nb-sTCNQ as an example, when adsorbed on Nb-sTCNQ, N2obtains 0.37(end-on) or 0.64 (side-on) electrons, and these electrons are all provided by NbN4(Figure 7). The charge density differences diagram of N2-TM-sTCNQ (Figure 8e–h′) more intuitively shows the charge transfer between N2and TM-sTCNQ. In subsequent hydrogenation steps along NRR, the charges of the three parts fluctuate significantly. For example, among the distal path on Nb-sTCNQ(Figure 7a), N2H obtains 0.15 electrons from Nb-sTCNQ in the first protonation step, while 0.19 electrons are transferred to NbsTCNQ in the fourth protonation step. During the whole NRR process on TM-sTCNQ (TM = Nb, Mo, Tc), TMN4is the major part involved in the charge transfer process, the TM-sTCNQ substrate without TMN4also plays some roles, and the cooperation between TM and sTCNQ promotes the reduction of N2to NH3.

    Table 2. The total magnetic moment of N2-TM-sTCNQ and its distribution.TM1 represents the center TM of the crystal lattice (with nitrogen); TM2 represents the corner TM atom of the lattice (without nitrogen); N or C(TM1-CN) means the N or C atom of the four CN groups connected to TM1, N or C (TM2-CN) means the N or C atom of the four CN groups connected to TM2, and N (N2) is the N atom of adsorbed N2 molecule.

    Figure 9. The band structure of the a) Nb-sTCNQ, b) Mo-sTCNQ and c) Tc-sTCNQ catalysts. Γ (0.0, 0.0, 0.0), X (0.5, 0.0, 0.0) and M (0.5, 0.5, 0.0) refer to high symmetry points in the first Brillouin zone in reciprocal space. The projected density of states (PDOS) of d) Nb-sTCNQ, e) N2-Nb-sTCNQ (end-on), f)N2-Nb-sTCNQ (side-on), g) Mo-sTCNQ, h) N2-Mo-sTCNQ (end-on), j)Tc-sTCNQ, k) N2-Tc-sTCNQ (end-on), and i) free N2. The PDOSs of TM-d and N-2p are plotted in blue and red, respectively.

    To further illustrate the electron acceptance/donation between the Nb-,Mo-,and Tc-sTCNQ catalysts and the N2molecules,the spin densities and charge density differences for N2-TM-sTCNQ (TM = Nb,Mo, Tc) are shown in Figure 8. For N2-Nb-sTCNQ (end-on) and N2-Mo-sTCNQ (end-on), the spin densities are mainly concentrated around the central TM atom (Figure 8a,a′,c,c′), due to the interactions between central TM atom and N2.Moreover,the N atom at the end of N2also contributes part of the spin densities, indicating there are unpaired electrons on the N2molecule, and N2is efficiently activated.Compared with N2-Nb-sTCNQ (end-on) and N2-Mo-sTCNQ (endon), the spin densities around center TM atoms and the N2molecules in the N2-Nb-sTCNQ(side-on)and N2-Tc-sTCNQ(end-on)structures are significantly smaller (Figure 8b,b′,d,d′), which is consistent with the ΔG*N2values. Table 2 lists the main magnetic moment contributions and total magnetic moment values for the N2-TM-sTCNQ structures.The magnetic moment difference between TM1and TM2can also reflect the interaction between N2and the active center TM atom. It shows that the N2-Nb-sTCNQ and N2-Mo-sTCNQ catalysts have stronger N2adsorption and activation ability, and are more inclined to adsorb N2in end-on configuration.As shown in the charge density differences diagram,the TM–N bond accumulates charges,while the N–N bond consumes charges,which also shows strong interactions between TM and N2and the activation of the N≡N bond during the adsorption process.

    To get a further understanding into the interactions between N2and TM-sTCNQ, the band structures and the projected density of states(PDOS)were investigated.As presented in Figure 9a–c,the Nb-,Mo-,and Tc-sTCNQ monolayers have negligible energy gaps(at least in one spin channel) with some bands cross the Fermi level, indicating threematerials are metallic. In addition, we use hybrid functional HSE06 to further verify the accuracy and reliability of the bandgap on the TMsTCNQ catalysts. The results show that the band gaps increase as expected. See Figure S9 and Table S10 for details. The valence electron configuration of free N2molecule is (1σg)2(1σu)2(1πu)4(2σg)2. The strong adsorption of N2on the TM-sTCNQ(TM = Nb,Mo,Tc)monolayer may be due to the "acceptance-donation" relationship between the central TM atom and adsorbed N2molecule. The PDOSs of Nb-,Mo-, and Tc-sTCNQ show that there are orbital interactions between TM-d and N-p on TCNQ (Figure 9d,g,j), and the interactions become much stronger after the catalyst adsorbs N2(Figure 9e,f,h,k).Thus,the TM-sTCNQ (TM = Nb, Mo, Tc) catalysts can effectively interact with N2, which is beneficial to the activation of N≡N, the charge transfer between TM-sTCNQ and N2,and the NRR process.

    Table 3. The binding energy (Eb, in eV) and the cohesive energy (Ec, in eV)of the TM-sTCNQ (TM = Nb, Mo, Tc) monolayers.

    3.4. The Stability of the TM-sTCNQ Monolayers

    Structural stability is the basic prerequisite of electrocatalytic materials.Especially,the active center TM atom can be embedded stably in the substrate and not reunite easily,which is closely related to the experimental feasibility and the long-term catalytic activity of SAC. Comparing the binding energies Ebof the Nb,Mo,and Tc atoms anchored on the sTCNQ monolayers with the cohesive energies Ecof the Nb,Mo,and Tc atoms,the stabilities of TM-sTCNQ(TM = Nb,Mo,Tc)are analyzed.As shown in Table 3,the Ebvalues of the TM-sTCNQ(TM = Nb,Mo,Tc)monolayers are ?10.58,?9.72,and ?9.77 eV,respectively,and the Ecvalues of Nb,Mo,and Tc are 7.44,6.85,and 7.38 eV,respectively.For all the three catalysts,Eb+ Ec< 0,indicating that the sTCNQ substrate can stabilize the Nb,Mo,and Tc atoms well,and the TM-sTCNQ(TM = Nb,Mo,Tc)monolayers have synthetic prospects.

    For practical application of the TM-sTCNQ (TM = Nb, Mo, Tc)monolayers under actual electrocatalytic reaction conditions, the thermal stability at room temperature or higher is very important. Thus,we performed AIMD simulations at 500 K to evaluate the thermal stability of TM-sTCNQ (TM = Nb, Mo, Tc) monolayers. When simulating (no symmetry constraints), a 2 × 2 × 1 supercell was employed to explore possible structural reconstruction and to minimize the effects of periodic boundary conditions.Our AIMD simulations of 10 ps time period at 500 K are given in Figures S10–S12.It is clear to see that the skeleton and integrity of the structures maintain well at 500 K although there are some wrinkles and undulations. The high thermal stability makes TM-sTCNQ (TM = Nb, Mo, Tc) fabricable and applicable as electrocatalysts under real reaction conditions.

    4. Conclusions

    In this work, spin-polarized DFT studies have been carried out on the structures and N2adsorption properties of the TM-sTCNQ monolayers, and high-throughput screenings of highly active NRR electrocatalysts have been conducted among them. Two promising NRR electrocatalysts have been obtained through the screenings, i.e.,the Nb-sTCNQ and Mo-sTCNQ monolayers. Especially, Nb-sTCNQ exhibits nice performance with a lower onset potential of only?0.48 V. Bader charge, charge differential, band structure, and PDOS analyses prove that TM-sTCNQ (TM = Nb, Mo, Tc) have good NRR catalytic activity. Binding energy analyses show that TMsTCNQ (TM = Nb, Mo, Tc) are stable materials and should be able to be synthesized experimentally. Among them, the NRR selectivity of Tc-sTCNQ material is poor. We hope that our research can promote further experimental and theoretical exploration of more 2D materials as efficient NRR electrocatalysts.

    Acknowledgements

    S.L., C.H., and L.Y. gratefully acknowledge support from the National Natural Science Foundation of China (22073033, 21873032, 21673087, 21903032), startup fund (2006013118 and 3004013105) from Huazhong University of Science and Technology, the Fundamental Research Funds for the Central Universities(2019kfyRCPY116), and the Innovation and Talent Recruitment Base of New Energy Chemistry and Device(B21003).S.-Y.L.,C.-X.H.,and G.L.gratefully acknowledge support from Guangdong Basic and Applied Basic Research Foundation(2021A1515010382) and the computational resources from the computing cluster at the Key Laboratory of Theoretical Chemistry of Environment, Ministry of Education&School of Chemistry,South China Normal University.The work was carried out at the LvLiang Cloud Computing Center of China, and the calculations were performed on TianHe-2. The computing work in this paper is supported by the Public Service Platform of High Performance Computing by Network and Computing Center of HUST.

    Conflict of Interest

    The authors declare no conflict of interest.

    Supporting Information

    Supporting Informationis available from the Wiley Online Library or from the author.

    Keywords

    2D TM-sTCNQ monolayers, density functional theory method,electrochemical nitrogen reduction reaction, high-throughput screening,single-atom catalysts

    Received: July 7, 2021

    Revised: September 4, 2021

    Published online: September 7, 2021

    [1] V. Smil, Nature 1999, 400, 415.

    [2] J. W. Erisman, M. A. Sutton, J. Galloway, Z. Klimont, W. Winiwarter,Nat. Geosci. 2008, 1, 636.

    [3] R. Schlogl, Angew. Chem. Int. Ed. 2003, 42, 2004.

    [4] S.L.Foster,S.I.P.Bakovic,R.D.Duda,S.Maheshwari,R.D.Milton,S.D.Minteer,M.J.Janik,J.N.Renner,L.F.Greenlee,Nat.Catal.2018,1,490.

    [5] S. Giddey, S. P. S. Badwal, C. Munnings, M. Dolan, ACS Sustainable Chem. Eng. 2017, 5, 10231.

    [6] A. J. Mart′?n, T. Shinagawa, J. P′erez-Ram′?rez, Chem 2019, 5, 263.

    [7] J. G. Chen, R. M. Crooks, L. C. Seefeldt, K. L. Bren, R. M. Bullock, M. Y.Darensbourg, P. L. Holland, B. Hoffman, M. J. Janik, A. K. Jones, M. G.Kanatzidis, P. King, K. M. Lancaster, S. V. Lymar, P. Pfromm, W. F. Schneider, R. R. Schrock, Science 2018, 360, eaar6611.

    [8] L. Li, C. Tang, X. Cui, Y. Zheng, X. Wang, H. Xu, S. Zhang, T. Shao, K.Davey, S. Z. Qiao, Angew. Chem. Int. Ed. 2021, 60, 14131.

    [9] J.N.Galloway,A.R.Townsend,J.W.Erisman,M.Bekunda,Z.C.Cai,J.R.Freney,L.A.Martinelli,S.P.Seitzinger,M.A.Sutton,Science 2008,320,889.

    [10] B. H. R. Suryanto, H.-L. Du, D. Wang, J. Chen, A. N. Simonov, D. R.MacFarlane, Nat. Catal. 2019, 2, 290.

    [11] H. Xu, K. Ithisuphalap, Y. Li, S. Mukherjee, J. Lattimer, G. Soloveichik, G.Wu, Nano Energy 2020, 69, 104469.

    [12] D. Z. Yao, C. Tang, L. Q. Li, B. Q. Xia, A. Vasileff, H. Y. Jin, Y. Z. Zhang,S. Z. Qiao, Adv. Energy Mater. 2020, 10, 2001289.

    [13] M. A. Shipman, M. D. Symes, Catal. Today 2017, 286, 57.

    [14] X. W. Lv, C. C. Weng, Z. Y. Yuan, ChemSusChem 2020, 13, 3061.

    [15] Y. Qiu, X. Peng, F. Lu, Y. Mi, L. Zhuo, J. Ren, X. Liu, J. Luo, Chem. Asian J. 2019, 14, 2770.

    [16] X.Liu,Y.Jiao,Y.Zheng,M.Jaroniec,S.Z.Qiao,J.Am.Chem.Soc.2019,141,9664.

    [17] X. F. Li, Q. K. Li, J. Cheng, L. Liu, Q. Yan, Y. Wu, X. H. Zhang, Z. Y.Wang, Q. Qiu, Y. Luo, J. Am. Chem. Soc. 2016, 138, 8706.

    [18] H.Y.Zhou,J.C.Li,Z.Wen,Q.Jiang,Phys.Chem.Chem.Phys.2019,21,14583.

    [19] X. Yao, Z. Chen, Y. Wang, X. Lang, W. Gao, Y. Zhu, Q. J. Jiang, Mater.Chem. A 2019, 7, 25961.

    [20] Q. Tang, D. E. Jiang, ChemPhysChem 2019, 20, 3141.

    [21] L. Xu, L. M. Yang, E. Ganz, ACS Appl. Mater. Interfaces 2021, 13, 14091.

    [22] X. Liu, Y. Jiao, Y. Zheng, S.-Z. Qiao, ACS Catal. 2020, 10, 1847.

    [23] X. Guo, J. Gu, S. Lin, S. Zhang, Z. Chen, S. Huang, J. Am. Chem. Soc.2020, 142, 5709.

    [24] C. X. Huang, G. Li, L. M. Yang, E. Ganz, ACS Appl. Mater. Interfaces 2021, 13, 608.

    [25] S. Y. Lv, C. X. Huang, G. Li, L. M. Yang, ACS Appl. Mater. Interfaces 2021, 13, 29641.

    [26] M.R.Zhao,B.Song,L.M.Yang,ACS Appl.Mater.Interfaces 2021,13,26109.

    [27] L. M. Yang, V. Bacic, I. A. Popov, A. I. Boldyrev, T. Heine, T. Frauenheim,E. Ganz, J. Am. Chem. Soc. 2015, 137, 2757.

    [28] B. Y. Song, Y. Zhou, H. M. Yang, J. H. Liao, L. M. Yang, X. B. Yang, E.Ganz, J. Am. Chem. Soc. 2019, 141, 3630.

    [29] G. Z. Zhu, L. S. Wang, J. Chem. Phys. 2015, 143, 221102.

    [30] T.-C. Tseng, N. Abdurakhmanova, S. Stepanow, K. Kern, J. Phys. Chem. C 2011, 115, 10211.

    [31] N. Abdurakhmanova, T. C. Tseng, A. Langner, C. S. Kley, V. Sessi, S. Stepanow, K. Kern, Phys. Rev. Lett. 2013, 110, 027202.

    [32] X. Ren, X. Q. Ji, Y. C. Wei, D. Wu, Y. Zhang, M. Ma, Z. Liu, A. M. Asiri,Q. Wei, X. P. Sun, Chem. Commun. 2018, 54, 1425.

    [33] Y. C. Wei, X. Ren, H. M. Ma, X. Sun, Y. Zhang, X. Kuang, T. Yan, D. Wu,Q. Wei, Chem. Eur. J. 2018, 24, 2075.

    [34] M. W. Xie, X. L. Xiong, L. Yang, X. F. Shi, A. M. Asiri, X. P. Sun, Chem.Commun. 2018, 54, 2300.

    [35] G. Kresse, J. Hafner, Phys. Rev. B: Condens. Matter. 1993, 47, 558.

    [36] G. Kresse, D. Joubert, Phys. Rev. B: Condens. Matter. 1999, 59, 1758.

    [37] J. P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 1996, 77, 3865.

    [38] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 2010, 132,154104.

    [39] J. H. Montoya, C. Tsai, A. Vojvodic, J. K. Norskov, ChemSusChem 2015,8, 2180.

    [40] D. Ma, Z. Zeng, L. Liu, X. Huang, Y. Jia, J. Phys. Chem. C 2019, 123,19066.

    [41] C. Ling, X. Bai, Y. Ouyang, A. Du, J. Wang, J. Phys. Chem. C 2018, 122,16842.

    [42] J. K. Norskov, J. Rossmeisl, A. Logadottir, L. Lindqvist, J. R. Kitchin, T. Bligaard, H. Jonsson, J. Phys. Chem. B 2004, 108, 17886.

    [43] Computational Chemistry Comparison and Benchmark Database,http://cccbdb.nist.gov/ (accessed: June 2020).

    [44] M. N. Faraggi, N. Jiang, N. Gonzalez-Lakunza, A. Langner, S. Stepanow,K. Kern, A. Arnau, J. Phys. Chem. C 2012, 116, 24558.

    我的亚洲天堂| 亚洲精品在线美女| www.自偷自拍.com| 久久精品影院6| 成在线人永久免费视频| 国产精品九九99| 波多野结衣高清无吗| 每晚都被弄得嗷嗷叫到高潮| 久久国产精品男人的天堂亚洲| 久久久久久国产a免费观看| 十八禁人妻一区二区| 国产精品亚洲av一区麻豆| 日韩大尺度精品在线看网址 | av有码第一页| 日韩有码中文字幕| 日日夜夜操网爽| 亚洲va日本ⅴa欧美va伊人久久| 久久欧美精品欧美久久欧美| 俄罗斯特黄特色一大片| 成人国语在线视频| 一本大道久久a久久精品| 大型黄色视频在线免费观看| 亚洲av日韩精品久久久久久密| 亚洲色图 男人天堂 中文字幕| 欧美另类亚洲清纯唯美| 亚洲av五月六月丁香网| 久久婷婷人人爽人人干人人爱 | 久久国产精品影院| 国产精品久久久久久精品电影 | 欧美丝袜亚洲另类 | 日本精品一区二区三区蜜桃| 黄片大片在线免费观看| xxx96com| 亚洲电影在线观看av| 黄色女人牲交| 丝袜人妻中文字幕| 夜夜躁狠狠躁天天躁| 午夜两性在线视频| 亚洲久久久国产精品| 午夜免费激情av| 色尼玛亚洲综合影院| 亚洲一区中文字幕在线| 日韩欧美一区二区三区在线观看| av网站免费在线观看视频| 午夜精品在线福利| 国产av又大| 三级毛片av免费| 国产aⅴ精品一区二区三区波| 精品福利观看| 一进一出好大好爽视频| a级毛片在线看网站| 校园春色视频在线观看| 亚洲三区欧美一区| 亚洲激情在线av| 最近最新免费中文字幕在线| av天堂久久9| 精品不卡国产一区二区三区| 男女床上黄色一级片免费看| 欧美色视频一区免费| 黄片播放在线免费| 黄色 视频免费看| 岛国视频午夜一区免费看| www.www免费av| 在线播放国产精品三级| 亚洲一区二区三区色噜噜| 午夜福利,免费看| 午夜福利视频1000在线观看 | videosex国产| 欧美丝袜亚洲另类 | 午夜福利在线观看吧| 日本三级黄在线观看| 国产乱人伦免费视频| 黄色女人牲交| 国产91精品成人一区二区三区| 色在线成人网| 黄色女人牲交| 99香蕉大伊视频| 国产精品电影一区二区三区| 男女午夜视频在线观看| 日韩av在线大香蕉| 国产激情久久老熟女| 无限看片的www在线观看| 色在线成人网| 国产精品久久久久久人妻精品电影| 熟女少妇亚洲综合色aaa.| 一边摸一边抽搐一进一出视频| 国产xxxxx性猛交| 国产成人精品无人区| 午夜福利成人在线免费观看| 国产亚洲av高清不卡| 宅男免费午夜| bbb黄色大片| 看片在线看免费视频| 18禁观看日本| 国产精品av久久久久免费| 女人被狂操c到高潮| 老司机深夜福利视频在线观看| 日韩成人在线观看一区二区三区| 一级,二级,三级黄色视频| 久久久久精品国产欧美久久久| 亚洲人成网站在线播放欧美日韩| 变态另类丝袜制服| 女人被狂操c到高潮| 日韩精品中文字幕看吧| 欧美老熟妇乱子伦牲交| 无人区码免费观看不卡| 欧美精品啪啪一区二区三区| 欧美黑人欧美精品刺激| 性少妇av在线| 欧美国产精品va在线观看不卡| 久久狼人影院| 国产av又大| 中文字幕高清在线视频| 涩涩av久久男人的天堂| 色精品久久人妻99蜜桃| 日日爽夜夜爽网站| 波多野结衣高清无吗| 久久人妻熟女aⅴ| 自拍欧美九色日韩亚洲蝌蚪91| 一边摸一边抽搐一进一出视频| 免费av毛片视频| 一区二区三区精品91| 高潮久久久久久久久久久不卡| 亚洲国产欧美日韩在线播放| 国产精品av久久久久免费| 精品日产1卡2卡| 亚洲精品一卡2卡三卡4卡5卡| 美女午夜性视频免费| 长腿黑丝高跟| 他把我摸到了高潮在线观看| 欧美黄色片欧美黄色片| 黑丝袜美女国产一区| 69av精品久久久久久| 热99re8久久精品国产| 国产aⅴ精品一区二区三区波| 日韩国内少妇激情av| 亚洲全国av大片| 欧美黄色片欧美黄色片| 国产精品久久久久久人妻精品电影| 两人在一起打扑克的视频| 精品国内亚洲2022精品成人| 老司机在亚洲福利影院| 99国产极品粉嫩在线观看| 99在线人妻在线中文字幕| 国产精品秋霞免费鲁丝片| 午夜福利影视在线免费观看| 久久精品aⅴ一区二区三区四区| 桃红色精品国产亚洲av| 国产免费男女视频| 99精品欧美一区二区三区四区| 亚洲人成电影免费在线| 老司机在亚洲福利影院| 一个人免费在线观看的高清视频| 日本在线视频免费播放| 午夜免费成人在线视频| 亚洲七黄色美女视频| 久久午夜综合久久蜜桃| 涩涩av久久男人的天堂| tocl精华| 在线国产一区二区在线| 日韩欧美一区二区三区在线观看| 99久久精品国产亚洲精品| 亚洲国产精品合色在线| 久久精品国产清高在天天线| 亚洲欧美激情综合另类| 97人妻天天添夜夜摸| 长腿黑丝高跟| 亚洲国产精品久久男人天堂| 亚洲精品在线观看二区| 国产精品久久视频播放| 亚洲午夜理论影院| 国产成年人精品一区二区| 人人妻人人澡人人看| 无限看片的www在线观看| 一进一出好大好爽视频| 首页视频小说图片口味搜索| 妹子高潮喷水视频| 免费女性裸体啪啪无遮挡网站| 国产精品 欧美亚洲| 欧美日韩黄片免| 天天添夜夜摸| 欧美久久黑人一区二区| 亚洲av成人不卡在线观看播放网| 18美女黄网站色大片免费观看| 啦啦啦韩国在线观看视频| 欧美日本亚洲视频在线播放| 真人一进一出gif抽搐免费| 91在线观看av| 亚洲av成人av| 中出人妻视频一区二区| 成人国产一区最新在线观看| 一个人观看的视频www高清免费观看 | 麻豆一二三区av精品| 国产精品香港三级国产av潘金莲| bbb黄色大片| 成年版毛片免费区| 久久久精品欧美日韩精品| 国产成人啪精品午夜网站| 女人被躁到高潮嗷嗷叫费观| cao死你这个sao货| 国产成人精品久久二区二区免费| 久久亚洲精品不卡| 日韩一卡2卡3卡4卡2021年| 又黄又爽又免费观看的视频| 97人妻精品一区二区三区麻豆 | 欧美黄色片欧美黄色片| 久久久久九九精品影院| 免费久久久久久久精品成人欧美视频| 欧美日本视频| 国产一区二区三区视频了| 一a级毛片在线观看| 性色av乱码一区二区三区2| 国产国语露脸激情在线看| 久久午夜综合久久蜜桃| 精品国产国语对白av| 国产在线观看jvid| 黑丝袜美女国产一区| 精品国产一区二区久久| 亚洲av成人一区二区三| 久久国产乱子伦精品免费另类| 国产三级黄色录像| 91国产中文字幕| 亚洲色图综合在线观看| 亚洲精品一区av在线观看| 大型av网站在线播放| 成年女人毛片免费观看观看9| av电影中文网址| 国产野战对白在线观看| 国产成人精品久久二区二区91| 欧美日韩亚洲国产一区二区在线观看| 两性午夜刺激爽爽歪歪视频在线观看 | 亚洲成av人片免费观看| 中文字幕久久专区| 国产亚洲精品第一综合不卡| 免费搜索国产男女视频| 高潮久久久久久久久久久不卡| 午夜福利18| 极品教师在线免费播放| 波多野结衣一区麻豆| 亚洲国产精品久久男人天堂| 搞女人的毛片| 男人舔女人下体高潮全视频| 亚洲av第一区精品v没综合| 狂野欧美激情性xxxx| 日韩精品中文字幕看吧| 精品人妻在线不人妻| 久久人妻熟女aⅴ| 一卡2卡三卡四卡精品乱码亚洲| 国产1区2区3区精品| 日本五十路高清| 久久精品国产亚洲av高清一级| 自拍欧美九色日韩亚洲蝌蚪91| 国产乱人伦免费视频| 18禁国产床啪视频网站| 午夜免费激情av| 怎么达到女性高潮| 欧美日本亚洲视频在线播放| 男女之事视频高清在线观看| 91在线观看av| 亚洲中文av在线| 色尼玛亚洲综合影院| 午夜影院日韩av| 久久香蕉激情| 亚洲成人免费电影在线观看| 亚洲国产欧美一区二区综合| 美国免费a级毛片| 人人妻人人爽人人添夜夜欢视频| 国产激情久久老熟女| 亚洲少妇的诱惑av| 国产成人一区二区三区免费视频网站| 少妇粗大呻吟视频| 真人一进一出gif抽搐免费| 亚洲精品一卡2卡三卡4卡5卡| 国内毛片毛片毛片毛片毛片| 看黄色毛片网站| 亚洲成人国产一区在线观看| 欧洲精品卡2卡3卡4卡5卡区| av中文乱码字幕在线| 给我免费播放毛片高清在线观看| 亚洲专区国产一区二区| 亚洲aⅴ乱码一区二区在线播放 | 九色国产91popny在线| 日韩 欧美 亚洲 中文字幕| 久久久精品国产亚洲av高清涩受| 国产欧美日韩一区二区三区在线| 欧美成狂野欧美在线观看| 国语自产精品视频在线第100页| 亚洲狠狠婷婷综合久久图片| 热re99久久国产66热| 黄片小视频在线播放| 色综合欧美亚洲国产小说| 伊人久久大香线蕉亚洲五| 一级,二级,三级黄色视频| 亚洲色图综合在线观看| 色哟哟哟哟哟哟| 国产激情欧美一区二区| 国产精品99久久99久久久不卡| 成人三级黄色视频| 国产精品久久久久久人妻精品电影| www.精华液| 亚洲熟女毛片儿| 国产私拍福利视频在线观看| 欧美色视频一区免费| 国产亚洲精品av在线| 搡老岳熟女国产| 中国美女看黄片| 午夜免费鲁丝| 国产在线观看jvid| 黄片大片在线免费观看| 精品高清国产在线一区| 国产成人系列免费观看| 午夜福利在线观看吧| 国产激情欧美一区二区| 在线视频色国产色| 久久天躁狠狠躁夜夜2o2o| 亚洲色图综合在线观看| 免费少妇av软件| 给我免费播放毛片高清在线观看| 夜夜看夜夜爽夜夜摸| 中亚洲国语对白在线视频| 欧美日韩瑟瑟在线播放| 18禁美女被吸乳视频| 精品久久久精品久久久| 亚洲狠狠婷婷综合久久图片| 国产亚洲精品av在线| 色综合站精品国产| 久久 成人 亚洲| 中文字幕人妻丝袜一区二区| 我的亚洲天堂| 变态另类丝袜制服| 看片在线看免费视频| 99久久99久久久精品蜜桃| 午夜两性在线视频| 成人手机av| 精品福利观看| 1024视频免费在线观看| 中文亚洲av片在线观看爽| 大型av网站在线播放| 日韩国内少妇激情av| 美女高潮喷水抽搐中文字幕| 亚洲色图 男人天堂 中文字幕| 搡老妇女老女人老熟妇| 国产熟女午夜一区二区三区| 亚洲中文日韩欧美视频| 美女免费视频网站| 精品欧美国产一区二区三| 亚洲精品av麻豆狂野| 免费在线观看影片大全网站| 国产精品av久久久久免费| 少妇的丰满在线观看| 一二三四在线观看免费中文在| 美女高潮喷水抽搐中文字幕| 日本撒尿小便嘘嘘汇集6| 9热在线视频观看99| 一个人免费在线观看的高清视频| 午夜久久久在线观看| 18禁观看日本| 日韩av在线大香蕉| 日日夜夜操网爽| 久久狼人影院| bbb黄色大片| 女同久久另类99精品国产91| 亚洲性夜色夜夜综合| 91大片在线观看| 亚洲精品一区av在线观看| 这个男人来自地球电影免费观看| 99国产精品99久久久久| 免费av毛片视频| 精品国产一区二区久久| 99精品在免费线老司机午夜| 国产av一区在线观看免费| 久久久久亚洲av毛片大全| 免费人成视频x8x8入口观看| 中文字幕av电影在线播放| 午夜福利欧美成人| 青草久久国产| 韩国av一区二区三区四区| 国产真人三级小视频在线观看| tocl精华| 国产99久久九九免费精品| 波多野结衣av一区二区av| 桃红色精品国产亚洲av| 天堂影院成人在线观看| 大型黄色视频在线免费观看| 欧美日韩亚洲综合一区二区三区_| 中文字幕最新亚洲高清| 国产男靠女视频免费网站| 18美女黄网站色大片免费观看| 精品高清国产在线一区| 国产高清激情床上av| 亚洲人成电影免费在线| 中文字幕精品免费在线观看视频| 88av欧美| aaaaa片日本免费| 欧美日韩瑟瑟在线播放| 亚洲精品粉嫩美女一区| 大码成人一级视频| 亚洲一区二区三区色噜噜| 免费看美女性在线毛片视频| 九色国产91popny在线| av中文乱码字幕在线| 这个男人来自地球电影免费观看| 国产激情欧美一区二区| 在线观看免费视频网站a站| 精品熟女少妇八av免费久了| 国语自产精品视频在线第100页| 91成年电影在线观看| 亚洲avbb在线观看| 老司机靠b影院| 亚洲av熟女| 啦啦啦韩国在线观看视频| 搡老岳熟女国产| 男人操女人黄网站| 欧美精品亚洲一区二区| 国产黄a三级三级三级人| 人妻久久中文字幕网| 久久婷婷成人综合色麻豆| а√天堂www在线а√下载| 脱女人内裤的视频| 亚洲国产看品久久| 成人手机av| 美国免费a级毛片| 一二三四社区在线视频社区8| 中出人妻视频一区二区| 免费在线观看影片大全网站| 亚洲欧美一区二区三区黑人| 欧美一级a爱片免费观看看 | 午夜精品在线福利| 天天一区二区日本电影三级 | 日韩大尺度精品在线看网址 | 少妇的丰满在线观看| 亚洲三区欧美一区| 51午夜福利影视在线观看| 亚洲精品av麻豆狂野| 欧美黑人精品巨大| 国产97色在线日韩免费| 成人特级黄色片久久久久久久| 久久人妻熟女aⅴ| 窝窝影院91人妻| 久久狼人影院| 久久这里只有精品19| 欧美精品亚洲一区二区| 国产av精品麻豆| 人人澡人人妻人| 午夜福利欧美成人| e午夜精品久久久久久久| 久久久久久免费高清国产稀缺| 久久人妻熟女aⅴ| 中文字幕av电影在线播放| 色综合亚洲欧美另类图片| 成人国语在线视频| 丝袜在线中文字幕| 亚洲欧美一区二区三区黑人| 成人精品一区二区免费| 免费av毛片视频| 国内精品久久久久久久电影| 天天一区二区日本电影三级 | 夜夜躁狠狠躁天天躁| 成人三级黄色视频| 午夜精品在线福利| 国产精品98久久久久久宅男小说| 免费高清在线观看日韩| 淫妇啪啪啪对白视频| av欧美777| 操美女的视频在线观看| 变态另类成人亚洲欧美熟女 | 可以在线观看毛片的网站| 国产精品久久久久久精品电影 | 亚洲免费av在线视频| 美女大奶头视频| 色在线成人网| 亚洲人成电影免费在线| 不卡av一区二区三区| 欧美日本中文国产一区发布| 精品日产1卡2卡| 亚洲av成人不卡在线观看播放网| 给我免费播放毛片高清在线观看| 国产精品乱码一区二三区的特点 | 高潮久久久久久久久久久不卡| 久久久国产欧美日韩av| 成年人黄色毛片网站| 欧美不卡视频在线免费观看 | 中文字幕高清在线视频| www.999成人在线观看| 多毛熟女@视频| 久久欧美精品欧美久久欧美| 亚洲色图 男人天堂 中文字幕| 久久国产精品影院| 这个男人来自地球电影免费观看| 日日夜夜操网爽| 国产91精品成人一区二区三区| 大码成人一级视频| 男人舔女人的私密视频| 久久精品成人免费网站| 熟妇人妻久久中文字幕3abv| 久久久久久大精品| 欧美一级a爱片免费观看看 | 国产单亲对白刺激| 一进一出抽搐gif免费好疼| 久久中文字幕一级| av天堂久久9| 亚洲片人在线观看| 国产精品久久电影中文字幕| 亚洲片人在线观看| 国产av一区在线观看免费| 91麻豆精品激情在线观看国产| aaaaa片日本免费| 侵犯人妻中文字幕一二三四区| 国产精品精品国产色婷婷| 看免费av毛片| 国产欧美日韩一区二区三| 亚洲第一欧美日韩一区二区三区| 黄色 视频免费看| 黄片大片在线免费观看| 亚洲色图av天堂| а√天堂www在线а√下载| 一级毛片精品| 一二三四社区在线视频社区8| 精品第一国产精品| av视频免费观看在线观看| 午夜久久久在线观看| 在线视频色国产色| 亚洲成a人片在线一区二区| 嫩草影视91久久| 欧美成人一区二区免费高清观看 | 女人爽到高潮嗷嗷叫在线视频| 国产黄a三级三级三级人| www.熟女人妻精品国产| 色在线成人网| 色尼玛亚洲综合影院| 亚洲av成人不卡在线观看播放网| 久久久国产成人免费| 一边摸一边做爽爽视频免费| 高清毛片免费观看视频网站| 久久久久国产精品人妻aⅴ院| 大陆偷拍与自拍| 日日摸夜夜添夜夜添小说| 亚洲精品美女久久av网站| 亚洲在线自拍视频| 美女午夜性视频免费| 精品久久蜜臀av无| 欧美+亚洲+日韩+国产| 可以在线观看的亚洲视频| 色哟哟哟哟哟哟| 麻豆国产av国片精品| 国产乱人伦免费视频| 亚洲aⅴ乱码一区二区在线播放 | 午夜a级毛片| 色综合欧美亚洲国产小说| 午夜影院日韩av| 国产精品秋霞免费鲁丝片| 此物有八面人人有两片| 99香蕉大伊视频| 精品国产超薄肉色丝袜足j| 99riav亚洲国产免费| 亚洲va日本ⅴa欧美va伊人久久| 久久久久久久午夜电影| 怎么达到女性高潮| 久久精品国产清高在天天线| 亚洲国产中文字幕在线视频| 91成人精品电影| 夜夜夜夜夜久久久久| 久久久国产欧美日韩av| 两人在一起打扑克的视频| 村上凉子中文字幕在线| 亚洲色图 男人天堂 中文字幕| 激情视频va一区二区三区| 色播亚洲综合网| 老司机在亚洲福利影院| 免费人成视频x8x8入口观看| 色婷婷久久久亚洲欧美| 国产精品精品国产色婷婷| 精品国内亚洲2022精品成人| 欧美老熟妇乱子伦牲交| www.熟女人妻精品国产| 国产免费av片在线观看野外av| 国产97色在线日韩免费| 91成人精品电影| 亚洲欧美激情综合另类| 亚洲国产欧美日韩在线播放| 在线观看日韩欧美| 神马国产精品三级电影在线观看 | 黑人巨大精品欧美一区二区蜜桃| 给我免费播放毛片高清在线观看| 不卡一级毛片| 中文字幕人妻熟女乱码| 一进一出抽搐gif免费好疼| 国产欧美日韩综合在线一区二区| 嫩草影视91久久| 精品熟女少妇八av免费久了| 国产欧美日韩一区二区三区在线| 桃红色精品国产亚洲av| 国语自产精品视频在线第100页| 亚洲精品中文字幕一二三四区| 在线十欧美十亚洲十日本专区| 色精品久久人妻99蜜桃| 国产99白浆流出| 麻豆国产av国片精品| 自线自在国产av| 欧美黄色淫秽网站| 国产又爽黄色视频| av天堂久久9| 亚洲激情在线av| 我的亚洲天堂| 女性生殖器流出的白浆| 好看av亚洲va欧美ⅴa在| 黄频高清免费视频| 九色亚洲精品在线播放| 18美女黄网站色大片免费观看| 成人特级黄色片久久久久久久| 波多野结衣一区麻豆| 亚洲中文字幕一区二区三区有码在线看 | 一本大道久久a久久精品| 国产精品精品国产色婷婷| 日韩国内少妇激情av|