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

    Accurate Deep Potential model for the Al–Cu–Mg alloy in the full concentration space*

    2021-05-24 02:25:22WanrunJiang姜萬潤YuzhiZhang張與之LinfengZhang張林峰andHanWang王涵
    Chinese Physics B 2021年5期

    Wanrun Jiang(姜萬潤), Yuzhi Zhang(張與之), Linfeng Zhang(張林峰), and Han Wang(王涵)

    1Songshan Lake Materials Laboratory,Dongguan 523808,China

    2Institute of Physics,Chinese Academy of Sciences,Beijing 100190,China

    3Beijing Institute of Big Data Research,Beijing 100871,China

    4Yuanpei College of Peking University,Beijing 100871,China

    5Program in Applied and Computational Mathematics,Princeton University,Princeton,NJ 08544,USA

    6Laboratory of Computational Physics,Institute of Applied Physics and Computational Mathematics,Beijing 100088,China

    Keywords: potential energy surface,deep learning,Al–Cu–Mg alloy,materials simulation

    1. Introduction

    Aluminum, copper, and magnesium (Al–Cu–Mg) based alloys are among the most versatile metallic materials that meet all the requirements of lightweight, high strength, and good fatigue resistance. They are intensively used in automotive, aviation, aerospace, and engineering machinery industries, and expend the scope of activity and the production efficiency of human civilization.[1]The industrial application of Al–Cu–Mg alloys could be traced back to 1906,when age hardening, as the key phenomenon enhancing their specific strength, was first discovered in this material.[2]It led to the well known wrought alloy, duralumin, and soon enabled structural members of Zeppelin airships and the world first passenger airplane without frame or bracing wire to support its wings. After a century, derived commercial aluminum alloy 2000 series are still one of the major material families used for aircraft manufacture.[3]Decades of explorations have linked up the conceptual road-map of age hardening in Al–Cu–Mg alloys, giving a mainline from the super saturated solid solution to a sequence of precipitations that finally impede the movement of dislocations.[4–6]The presence, morphology, and competition of micro-structures were observed by experiments,[6–9]while their chemical formulas,atomic configurations, and relative stability were informed by theoretical calculations.[5,10,11]Phenomenological experiences have been accumulated along the road-map, revealing further complexity of the phenomenon including, but being not limited to,the involvement of distinguishing processes like dislocation–solute interaction and solute clustering in different ageing stages,[12,13]the evolution of micro-structures under various ageing conditions,[14,15]the ageing improvements on desired performance like dimensional stability,[16]and variants of related processes in more complex scenario like Al–Cu–Mg nanocrystals.[17]

    A better understanding of the fundamental principles underlying the experiments needs the atom-resolved tracking of the dynamical processes, which is usually beyond the capability of current experiment equipment. Therefore, molecular simulation techniques like molecular dynamics (MD) are a promising alternative, which is in principle capable to produce atom-scale information of the physical quantities of interest in the investigated processes. MD simulations have provided insights into the precipitation strengthening of Al–Cu alloys, which mechanisms might also be inspiring for understanding related processes in Al–Cu–Mg systems. For example, the dislocations–precipitate interactions, that determine the strengthening effect, are dominated by a looping process for most edge dislocation cases[18]and are significantly influenced by the temperature-dependant cross-slip for screw dislocations.[19]The transition among meta-stable precipitates could be adjusted by the vacancy segregation, while the vacancy densities could further be controlled by heating the surface of the designed nanoscale specimen or exposing bulk counterparts to intense irradiation.[20]Besides,possible common factors for precipitation strengthening in different alloys could also be perceived by MD simulations,such as the importance of local interface curvature in driving the dislocation–precipitate interactions.[21]

    Ideally,ab initio molecular dynamics(AIMD),[22]which evolves atomic systems under the Newton’s law of motion with forces generated on-the-fly by first-principles-based methods, holds the promise to provide an accurate description of Al–Cu–Mg alloy systems. However, due to the high computational cost, first-principles-based methods like density functional theory(DFT)[23]can only be used to simulate materials within the spacial scale of a few hundreds of atoms and the temporal scale of a few dozens of picoseconds, using desktop workstations or ordinary computer clusters. Despite successes in computing basic quantities like structural factors, heat of formation, and cohesive energy for important precipitation,[10,11,24]they are in general not computationally feasible to reveal the atomistic mechanism in large-scale dynamic processes like phase transition, precipitation, and interface migration under working conditions, whose spacialtemporal scales are far beyond the current capability of AIMD.

    In these cases, one has to resort to empirical force fields(EFFs). Representative examples include the Lennard–Jones potential,[25]the Stillinger–Weber potential,[26]the embedded-atom method (EAM) potential,[27]and the modified embedded atom method (MEAM) potential.[28]By approximating the potential energy surface(PES)with relatively simple and analytical functions and optimizing the parameters with a group of target properties,EFF-based MD is orders of magnitude more efficient than AIMD.However,it has been challenging to construct EFF models for multi-component systems like the Al–Cu–Mg alloy. The difficulties are two-folds:(1) The relatively simple function form makes it difficult to fit all target properties in the full concentration space, especially in the case of ternary alloy. (2)The number of tunable parameters in the EFF grows quadratically with respect to the number of components of the system. The number of parameters in an M-component system is M(M+1)/2 times larger than that of a single-component system,so the optimization of the parameters becomes more difficult, if not impossible, for multi-component systems.

    For the Al–Cu–Mg system, most of the existing EFF models are restricted to pure metals[29–36]and binary components.[37–42]Although some of them have been demonstrated to provide decent descriptions of the Al end in the Al–Cu phase diagram[39]and the solid–liquid phase equilibria in the target component range of Al–Mg alloy,[42]the generalization of these single or binary component EFFs to ternary Al–Cu–Mg alloy is not straightforward. One of the few examples of multi-component EFF[43]that covers Al, Cu, and Mg is constructed within the formalism of MEAM potential,[28]and is explicitly optimized for single-element systems of Al,Si, Mg, Cu, and Fe, their binary pairs in a hypothetical NaCl reference structure,as well as some stable binary intermetallic components. The fitting targets include a group of selected properties obtained from ab initio calculations, such as the generalized stacking fault energy for single elements, heat of formation, equilibrium volume, elastic moduli of hypothetical binary pairs,and heat of formation for some stable binary components. The transferability of this MEAM potential to the full concentration space was implicitly suggested,but was found not satisfactory for fundamental properties partially beyond the fitted sets for developing the potential.[44]Therefore,there is a demand for a PES model that is substantially less expensive than ab initio methods,and,in the meanwhile,provides remarkably higher accuracy and wider applicable range than the state-of-the-art EFFs. Ideally, it should describe the static, mechanical, and dynamical properties with a satisfactory accuracy in the full concentration space of the ternary Al–Cu–Mg system.

    Recent advances in the machine learning (ML)techniques[45–50]provide a promising route towards an accurate and efficient PES model for ternary Al–Cu–Mg alloys. Here we consider the smooth edition of Deep Potential(DP),[50]a deep learning based PES representation. It has been demonstrated that, when trained with ab initio data, DP agrees well with an ab initio PES and has an efficiency comparable to that of EFFs. The representation ability of DP stems from the outstanding ability of the deep neural networks to approximate high-dimensional functions,[51,52]and from a group of automatically generated symmetry-preserving descriptors which faithfully capture the information lying in the local environment of atoms. Compared with most ML models, the flexible and trainable descriptors of DP make it more suitable for complex tasks, e.g., multi-component systems like TiO2–H2O interface,[53]strengthening precipitates from aluminum alloys,[20]as well as a (Zr0.2Hf0.2Ti0.2Nb0.2Ta0.2)C high entropy alloy system with six chemical species.[54]To construct the DP model with an optimal set of data,a concurrent learning strategy,the deep potential generator(DP-GEN),[44]is further established. Implemented as a close-loop iterative workflow,DP-GEN could generate the most compact and adequate data set that guarantees the uniform accuracy of DP in the explored configuration space. The effectiveness and efficiency of DPGEN have been validated in the cases of pure Al, pure Mg,and the Al–Mg alloy,[44]as well as pure Cu,[55]where the DP models outperform the state-of-the-art MEAM potentials in almost all examined properties.

    In this paper,we construct the DP model for the Al–Cu–Mg alloys using the DP-GEN scheme for the whole concentration space,i.e.,AlxCuyMgzwith 0 ≤x,y,z ≤1,x+y+z=1, and the configuration space covering a temperature range around 50.0 K to 2579.8 K and a pressure range around 1 bar to 50000 bar. The resulting model gives generally better consistency with DFT results than the state-of-the-art MEAM model[43]when predicting basic energetic and mechanical properties for systems distributed in the full concentration space of Al–Cu–Mg alloy.The tested systems are taken from the Materials Project (MP) database, and most of them are not explicitly covered by the training set. The tested properties include formation energy and volume for equilibrium state,equation of state(EOS),elastic modulus,unrelaxed surface formation energy for both binary and ternary alloys, as well as the relaxed vacancy formation energy and interstitial energy for ternary alloys.

    In the following,we will first give a brief introduction of the DP model and DP-GEN scheme for the Al–Cu–Mg system in Section 2. Then in Section 3,we will focus on the comparison of various properties predicted by DFT,DP,and MEAM.Finally,a conclusion will be given in Section 4.

    2. Method

    The DP-GEN scheme works in an iterative manner. Each iteration is composed of three steps,exploration,labeling,and training. With an initial guess of an ensemble of DP models generated from a relatively simple initial data set, DP-GEN starts to explore the configurational and chemical space by a sampler driven by DP-based simulation. The explored configurations are then classified as accurate, candidate, and failed,according to the model deviation defined as the standard deviation of the atomic force predictions of an ensemble of models.Models in the ensemble use the same neural network architecture and are trained on the same data set. The only difference among them is the random seeds used to initialize parameters of the neural network. Thus,their predictions would converge if the configuration has been well covered by current data set and deviate from each other otherwise. The accurate set consists of the configurations with the model deviation lower than a user-provide lower bound. The failed configurations are those with the model deviation higher than a userprovided upper bound. The failed configurations are typically generated by the DP models with a relatively poor quality in early DP-GEN iterations,and are likely to be unphysical,e.g.,with overlapping atoms. A subset of the candidate configurations with errors between the lower and upper bounds are then sent to the labeling step, at which the labels of a configuration, i.e., the energy, the virial tensor, and the atomic forces,are computed by DFT. Finally, the labeled configurations are added to the existing data set. From the updated data set a new ensemble of DP models are trained,and a new iteration starts from the exploration step. The DP-GEN iterations converge when almost all the explored configurations are classified as accurate. More details of the DP-GEN scheme are referred to previous works.[44,55]

    It should be remarked that(1)the quality of the final DP model is insensitive to the data set that constructs the initial DP model to start the DP-GEN scheme. However,with a poor initial guess of DP,more labeled configurations from later iterations would be needed to improve the model. (2)The conditions under which the DP model is used should be determined by the user beforehand,and the exploration strategy is accordingly designed to ergodically explore the relevant configuration space. In this context,the sampler can be,but is not restricted to, molecular dynamics, Monte–Carlo simulation,structure prediction schemes,or enhanced sampling methods.(3) By classifying and sub-sampling the explored configurations, only the configurations of low accuracy are sent to the labeling step, so the most compact data set, i.e., the minimal data set that helps DP achieve a certain accuracy,can be generated during the DP-GEN iterations.

    The DP-GEN iterations are conducted by the open-source package DP-GEN.[55]More details on initial data set, exploration,labeling,and training steps are provided as follows.

    Initial data set The initial data set used to kick-off the training and exploration can be generated through an automatic workflow integrated in the DP-GEN scheme, which allows the whole scheme to start from little prior knowledge.For alloys, four procedures are sequentially conducted. First,for each user-specific crystalline structure and each possible concentration,the atomic species is randomly assigned on the lattice points. Then, random perturbations are performed on the coordinates by adding values drawn from a uniform distribution in the range[?pa,pa].Perturbations are also performed on the cell vectors by a symmetric deformation matrix that is constructed by adding random noise drawn from a uniform distribution in the range[?pc,pc]to an identity matrix. Next,starting from the perturbed structures, AIMD simulations are conducted to produce labeled data with DFT calculated energy, force, and virial tensor. Finally, the labeled data are used to create the initial data sets. The preparation for singleelement systems follows the same procedure without the need to consider possible concentrations. Further details for preparing initial data sets can be found in previous reports.[44,55]It is worth noting that,users can generate and collect initial data sets through other appropriate manners and tools. The procedure mentioned above just provides a convenient protocol and the quality of the finial model is not sensitive to the initial data sets.

    In this work, for binary and ternary alloys, we consider body-centered-cube (BCC), hexagonal-close-packed (HCP),and face-centered-cube (FCC) structures, within 2×2×2 super-cells. For perturbations, pais set to 0.01 ?A and pcis set to 3%. For single-element cases, training data sets for copper are inherited from previous work,[55]while those for aluminium and magnesium are prepared through individual single-element DP-GEN procedure in similar strategies. Besides, scaling factors compressing the cell are used in the above procedure for binary and ternary alloys. Generated labeled data are added to the final data sets to further cover the high-pressure regime of the configuration space. Specific values are from 0.84 to 0.98 with the interval of 0.02 for binary alloys and from 0.85 to 0.97 with the interval of 0.03 for ternary alloys.

    Exploration In this work, the configuration and chemical space are explored by DP-based molecular dynamics(DeePMD) using the LAMMPS package.[56]Two types of structures, bulk and surface, are considered. For bulk systems,the initial configurations are constructed from randomly perturbed FCC, BCC, and HCP, in 2×2×2 super-cells with lattice positions randomly occupied by Al, Mg, or Cu atoms.It is noted that the 2×2×2 super-cells are large enough to capture the local atomic environments needed to construct the accurate DP model for the Al–Cu–Mg alloys(see the supplementary text, Figs. S1 and S2 for the accuracy of the model when predicting alloy structures with larger super-cells up to 12×12×12). Isothermo-isobaric(NpT)MD simulations are conducted at thermodynamic states ranging from 50.0 K to 2579.8 K (twice the experimental melting point of Cu) and from 1 bar to 50000 bar. In the first 40 iterations (0–39), binary alloys(Al–Cu,Al–Mg,and Cu–Mg)are explored. Then from iterations 40 to 90,configurations of ternary alloys(Al–Cu–Mg)are sampled(see supporting information Table S1 for more details). It is noted that in each iteration, NpT simulations of 2 temperatures and 8 pressures covering the range 1–50000 bar are performed, thus configurations are sampled from 16 different thermodynamic conditions. The temperature increases as the iterations proceed, and the pressures are set to fixed values (1, 10, 100, 1000, 5000, 10000, 20000,50000 bar).The surface initial configurations are generated by displacing two halves of FCC and HCP structures along certain crystallographic directions. The structures are prepared in two ways:(1)Conventional cells are copied along the displacing direction so that the thickness of the surface structures is at least 10 ?A.(2)Conventional cells are replicated by 2×1×1 or 2×2×1 times (see supporting information Table S1 for the assignment of replications to crystals),and then are copied along the displacing direction so that the thickness of surface is at least three layers of atoms. Such sample sizes allow a balance between the rationality in simulating surface environments and the efficiency in successive labeling. We choose{100,110,111}and{001,100,110}as directions of displacing for FCC and HCP, respectively. The atomic positions of these displaced structures are randomly occupied by Al, Mg,or Cu atoms. Simulation cells prepared by(1)are scaled with scaling factors from 0.96 to 1.06 with an interval of 0.02,and those prepared by (2) use scaling factors of 0.96, 1.00, and 1.06. Finally,the atomic positions and the cell shapes are randomly perturbed.Canonical(NVT)simulations are conducted to sample the surface configurations in the temperature range of 50.0–1222.0 K.Binary and ternary alloys are explored from iteration 91 to 94,while Al–Cu and selected ternary alloys are explored from iteration 95 to 103(see supporting information Table S1 for more information). The explored configurations are classified as accurate when the model deviation is smaller than 0.05 eV/?A,as candidate if the model deviation is between 0.05 eV/?A and 0.20 eV/?A,and as failed otherwise.

    Labeling The Vienna ab initio simulation package(VASP)[57,58]is used to conduct the DFT calculation with the generalized gradient approximation(GGA)adopting Perdew–Burke–Ernzerhof (PBE)[59,60]exchange–correlation functional. The projector-augmented-wave (PAW) method[61,62]is used with an energy cut-off of 650 eV for the plane wave basis sets. The K-points in the Brillouin zone are sampled by the Monkhorst–Pack mesh[63]with grid spacing of hk=0.1 ?A?1. The order 1 Methfessel–Paxton smearing method[64]is adopted with σ=0.22 eV.The convergence criterion for the self consistent field (SCF) calculation is 1×10?6eV for the changes of both the total energy and the band structure energy.

    Training The smooth edition of the deep potential model is adopted to construct the PES model.[50]The DeePMD-kit package[65]is used for training. The sizes of the embedding and fitting nets are set to(25,50,100)and(240,240,240),respectively. During the DP-GEN iterations, the cut-off radius is set to 6 ?A. The learning rate starts from 1×10?3and exponentially decays to 3.5×10?8after 1×106training steps.Four models are trained to construct the model ensemble. Using the same architecture and the same training data set,they only differ in the random seeds for initializing model parameters. After the DP-GEN iterations being converged,the production models are trained with the cut-off radius set to 9 ?A and training steps set to 1.6×107. In all training tasks, the Adam stochastic gradient descent method[66]is used with default hyper-parameter settings provided by the TensorFlow package.[67]

    3. Results and discussion

    After the convergence of the DP-GEN process,about 2.73 billion alloy configurations are explored and only a small portion of them(~0.0038%)are selected for labeling. Data sets produced by this work have been uploaded to the online open data repository.[68]

    In the following we systematically benchmark the DP model by testing its accuracy in terms of formation energy,equilibrium volume,equation of state(EOS),elastic modulus,surface formation energy,vacancy formation energy,and interstitial energies. We consider the 58 alloy structures from the Materials Project(MP)database,including 15 Al–Cu,30 Al–Mg, and 7 Cu–Mg binary structures, as well as 6 Al–Cu–Mg ternary structures.The properties are calculated from the alloy structures relaxed at 0 bar.A state-of-the-art MEAM[43]model is used for comparison. DFT calculations of these properties adopt the same numerical setup as those used in the labeling step. When computing the formation energy, equilibrium volume, interstitial energy as well as vacancy defect and surface formation energies,the ionic relaxation stops when all the forces in a configuration are smaller than 1×10?2eV/?A.The elastic moduli are calculated by the strain–stress relation using the finite difference method,so it demands a tighter break criterion for both equilibrium and deformed structures. Here we terminate the ionic step when the energy change in two successive steps is smaller than 1×10?5eV.DP and MEAM ionic relaxations break when either one of the following two criteria is satisfied. (1) The energy change between two successive steps is smaller than one part in 1×1012of the total energy,or(2)the norm of the global force vector(could be approximately considered as the upper bound of any force component on any atom) is smaller than 1×10?6eV/?A. All the computations and analyses of the properties mentioned above are facilitated by the automatic workflow implemented in DPGEN,[55]which invokes VASP[57,58]for DFT calculations and LAMMPS[56]for DP and MEAM calculations. The vacancy and interstitial structures are prepared by exhaustively generating all non-equivalent defect structures from the MP structures by the the python package pymatgen[69](interstital structures are generated by the interstitialcy finding tool(InFit),[70]which is integrated in the pymatgen package).

    3.1. Formation energy and equilibrium volume

    The formation energy of binary alloys,taking Al–Cu for example,is defined as

    and that of the ternary Al–Cu–Mg alloy is defined as

    where E0AlCuand E0AlCuMgdenote the equilibrium energies per atom of the tested Al–Cu and Al–Cu–Mg alloys,respectively.cAland cCudenote the component concentrations of Al and Cu,respectively,in an alloy structure. E0Al, E0Cu,and E0Mgdenote the equilibrium energies per atom of the three pure metals in their stable crystalline structures. Here Al and Cu take the FCC structure,and Mg takes the HCP structure.

    Fig.1.Comparison of the formation energy and the equilibrium volume predicted by DP and MEAM for alloys, with respect to DFT results.(a) Formation energy with dashed lines representing y=x±0.03 and the inserted histogram showing the distribution of absolute errors with respect to the DFT reference. (b)Equilibrium volume with dished lines representing y=x±0.2 and the inserted histogram showing the distribution of errors. 58 structures,including both binary and ternary alloy structures,are adopted from Materials Projects(MP)database.

    The DP predictions of formation energy and equilibrium volume for all tested structures are in satisfactory consistency with the DFT results. In comparison with the MEAM results, the relative errors with respect to DFT references are generally smaller. As shown by Fig. 1(a), the absolute differences between DP and DFT formation energies are below 10 meV·atom?1for 39 MP structures (67% in total 58 MP structures), and no tested case has an error larger than 100 meV·atom?1. As shown in Fig. 1(b), 54 MP structures,or 93% of the tested structures, have absolute errors in equilibrium volume of <0.2 ?A3·atom?1. Compared with MEAM results, the errors in the formation energy (in eV·atom?1) of most of the testing cases are reduced by at least one order of magnitude, and the relative errors in equilibrium volume are generally reduced from several percents to beneath 1%.About 53%of the MEAM results give errors over 100 meV·atom?1.About 17%in total give errors larger than 200 meV·atom?1,so they are beyond the plotted range of the distribution of errors.

    3.2. Equation of state

    Next, we investigate the energy dependence on the volume region around the equilibrium value (80%–120% of the equilibrium volume). The performance of DP on the equation of state (EOS) for both binary and ternary alloys is shown in Fig. 2. The absolute DP root-mean-squared errors(RMSEs) with respect to DFT references are smaller than 10 meV·atom?1for 34 out of 57 tested alloys (MP-1200279 is not included due to the convergence difficulty in the DFT calculation), which is close to the performance on equilibrium energy (39 out of 58 structures with errors smaller than 10 meV·atom?1). This indicates that the DP model generally holds its accuracy for the equilibrium state in this volume range. As a comparison,for more than half of the tested cases, the RMSEs of DP predictions are two orders of magnitudes smaller than MEAM results,thus DP ensures the generally more reliable energy–volume dependence than MEAM for crystals in the whole Al–Cu–Mg concentration space.

    Fig.2. Comparison of DP and MEAM root-mean-squared error(RMSE)of EOS for alloys, with respect to DFT results. Structures are downloaded from the MP database. Volume range for computing RMSE is approximately from 80% to 120% of the DFT equilibrium volume. (a) 5 Al–Cu–Mg structures, 2 of them give MEAM errors larger than the plotting range: (a1) 0.600 eV·atom?1; (a2) 0.663 eV·atom?1. MP-1200279 is not included due to the convergence difficulty in the DFT calculation;(b)7 Cu–Mg structures; (c) 15 Al–Cu structures; (d) 30 Al–Mg structures. The DP and MEAM EOS of MP-1039192 are estimated by single-point energy calculations on converged structures from the DFT EOS computation,because the MP-1039192 structure is not stable in neither DP nor MEAM energy minimization.

    We note that,for all the properties introduced above,despite an overall good agreement with DFT values, DP gives unsatisfactory predictions on a small number of structures. A detailed investigation suggests that this may be attributed to the exploration strategy adopted by the DP-GEN process. The MD simulation for exploring the conformation space starts from structures with FCC,HCP,and BCC lattices,and should include the corresponding crystalline, molten, and surface structures. However, there are local chemical environments that can be hardly sampled using this strategy,so that DP will be in the extrapolation regime when tested on the associated structures.

    Examples include MP-1038818 and MP-1039010 in Fig.2,which are in diamond and simple cubic lattices,respectively. It should be mentioned that, however, if we continue with the DP-GEN process and add related structures into the exploration step, the DP model should be systematically improved. As an evidence, we refer to Ref. [44] for an Al–Mg model generated by a DP-GEN process that takes into considerations the diamond and simple cubic structures. Such a model shows a much better prediction accuracy on MP-1038818 and MP-1039010.

    3.3. Elastic modulus and surface formation energy

    Now we consider the bulk and shear moduli as well as the unrelaxed surface formation energy. The unrelaxed surface formation energy,Esf((lmn)),refers to the energy needed to create a surface with the miller index(lmn)for a given crystal,and is defined as

    Here Es((lmn)) and Nsdenotes the energy and number of atoms of the unrelaxed surface structure with the Miller index (lmn). A denotes the surface area. For each MP case,E0refers to the per-atom potential energy of the relaxed bulk structure (under 0 bar) of this MP crystal. For each surface,energy evaluations of Es((lmn))using three methods are conducted on the same surface geometry,i.e.,the one generated by displacing two halves of the DFT equilibrium bulk structure.This ensures the rationality of the comparison, since a small difference in the equilibrium bulk structure relaxed by different methods and codes could bring distinct surface geometries under the same Miller index(lmn)(surface generation method is similar to that described in Section 2).

    As shown in Fig.3(a),DP generally gives more accurate bulk and shear moduli than MEAM. It presents a more concentrated distribution of data points along the diagonal line denoting the DFT reference. This advantage is more directly revealed by comparing the absolute errors with respect to the DFT reference. Within the leftmost interval of the inserted histogram that allows errors beneath 5 GPa, about 20%more results in total are presented by DP than MEAM.Meanwhile,the tail of the distribution of error is also shorter for DP.Furthermore, both bulk and shear moduli of all 5 tested ternary cases are more accurately predicted by DP.

    DP also outperforms MEAM on predictions of unrelaxed surface formation energy, as shown by Fig. 3(b). Taking DFT results as references and compared with MEAM, there are 50% more data points with absolute errors smaller than 0.05 J/m2when using DP for predictions. However, DP slightly underestimates Esffor some ternary cases and a small portion of binary cases.

    Fig.3. Comparison of the elastic moduli and the unrelaxed surface formation energy predicted by DP and MEAM for alloys, with respect to DFT results. (a) The bulk modulus Bv and the shear modulus Gv of 57 MP structures, with dashed lines representing y=x±10 and the inserted histogram showing the distribution of absolute errors with respect to the DFT reference. MP-1200279 is not included due to the convergence difficulty in the DFT calculation.(b)911 unrelaxed surface formation energies of 56 MP structures,with dashed lines representing y=x±0.1 and the inserted histogram showing distribution of errors.MP-1200279 and MP-1185596 are not included due to the convergence difficulty in the DFT calculation.

    3.4. Vacancy and interstitial defects

    Finally,for ternary alloys,we examine the vacancy defect formation energy and the energy along interstitial relaxation paths. The vacancy formation energy Evfis defined as

    where Ev(Nv) denotes the energy of a relaxed structure with one vacancy defect and Nvatoms. E0is the equilibrium energy per atom of the corresponding pristine crystal.

    For the interstitial defect, due to the generally high energetic instability of the generated interstitial structures, the relaxation may not end up with structures with locally relaxed interstitial point defects. Meanwhile,the relaxation path would be highly sensitive to the practical implementation of the relaxation algorithm.[44]Therefore, instead of comparing the interstitial formation energies of DFT,DP,and MEAM calculated by independent relaxations,we investigate the DP and the MEAM prediction errors along the DFT relaxation path(excluding early high energy geometries).

    Fig.4. Comparison of the vacancy formation energy and the unrelaxed surface formation energy predicted by DP and MEAM for ternary alloys,with respect to DFT results. (a)44 relaxed vacancy formation energies of 6 ternary alloys,with dashed lines representing y=x±0.5 and the inserted panel showing the error distribution.(b)23168 energies per atom along DFT relaxation path of 5 ternary alloys, with dashed lines representing y=x±0.02 and the inserted panel showing the error distribution.MP-1200279 is not included due to the convergence difficulty in the DFT calculation.

    DP is in satisfactory consistency with DFT and shows relative advantages compared with MEAM for both properties.As shown by Fig. 4(a), almost all DP predicted Evfvalues have absolute errors beneath 0.5 eV with respect to DFT values. It can be found that, DP tends to slightly underestimate Evfcompared with the DFT reference. In terms of the energy along DFT interstitial relaxation paths, nearly all DP predictions give absolute errors beneath 1 meV·atom?1. In contrast,MEAM gives about 30% less data points in total within the same error range. Meanwhile,some MEAM predictions have errors of tens of meV·atom?1,as demonstrated by those points far from the DFT reference line in Fig.4(b).

    It is noted that the vacancy and interstitial formation energies calculated by the tested structures may not be well converged due to the finite size effect.However,for the purpose of investigating the accuracy of the DP and the MEAM in terms of the defect formation energy,it is reasonable to compare the results of DP and MEAM with that of the DFT,computed from the defect structures with the same number of atoms.

    4. Conclusion

    Using the deep potential representation and the concurrent-learning scheme, DP-GEN, we develop an accurate PES model for the Al–Cu–Mg alloys in the full concentration range. Systematic benchmarks with DFT results on fundamental energetic, mechanical, and defect properties of 58 MP structures suggest a good accuracy in the full concentration space, and an overall better performance than a stateof-the-art MEAM potential. The training does not rely on any existed first-principles data set. Instead,data is generated on-the-fly through DP-GEN and finally forms a compact set supporting reliable DP model after the scheme is converged.The resulting DP model and data sets are expected to help investigate the physical mechanisms behind complex phenomena of Al–Cu–Mg systems and contribute to the methodology development for materials simulation.

    欧美性猛交╳xxx乱大交人| 亚洲精品乱码久久久v下载方式| 国产激情偷乱视频一区二区| 国产免费一级a男人的天堂| 成人国产麻豆网| 亚洲综合色惰| 色哟哟·www| 最近手机中文字幕大全| 十八禁国产超污无遮挡网站| 欧美最新免费一区二区三区| 日日撸夜夜添| 男人的好看免费观看在线视频| 欧美高清性xxxxhd video| 日韩国内少妇激情av| 亚洲av中文av极速乱| 亚洲熟妇中文字幕五十中出| 1024手机看黄色片| 国产探花在线观看一区二区| 天堂中文最新版在线下载 | 在线观看美女被高潮喷水网站| 成人无遮挡网站| 精品日产1卡2卡| 伦精品一区二区三区| 国产av在哪里看| 亚洲国产欧洲综合997久久,| 国产精品嫩草影院av在线观看| 午夜久久久久精精品| 国语自产精品视频在线第100页| 婷婷色综合大香蕉| 免费看美女性在线毛片视频| 亚洲精品456在线播放app| 亚洲无线观看免费| 亚洲欧美日韩无卡精品| 精品免费久久久久久久清纯| 国产精品久久久久久久久免| 欧美不卡视频在线免费观看| 亚洲美女视频黄频| 高清毛片免费观看视频网站| 国产成人精品一,二区 | 精品日产1卡2卡| 一个人看视频在线观看www免费| 人人妻人人澡欧美一区二区| 综合色av麻豆| 成人二区视频| 国产黄片美女视频| 国产精品av视频在线免费观看| 亚洲七黄色美女视频| 如何舔出高潮| 变态另类成人亚洲欧美熟女| 爱豆传媒免费全集在线观看| 高清毛片免费观看视频网站| 中文字幕熟女人妻在线| 日韩 亚洲 欧美在线| av专区在线播放| 99久久久亚洲精品蜜臀av| 波多野结衣高清作品| 给我免费播放毛片高清在线观看| 丝袜喷水一区| av在线蜜桃| 亚洲18禁久久av| 日本欧美国产在线视频| 国产欧美日韩精品一区二区| 国产极品天堂在线| 国产高清有码在线观看视频| 亚洲av不卡在线观看| 少妇熟女aⅴ在线视频| 成人三级黄色视频| 一个人观看的视频www高清免费观看| 如何舔出高潮| 联通29元200g的流量卡| 精品一区二区三区人妻视频| 亚洲精品色激情综合| 亚洲高清免费不卡视频| а√天堂www在线а√下载| 亚洲aⅴ乱码一区二区在线播放| 乱人视频在线观看| 亚洲国产精品成人久久小说 | 在线观看美女被高潮喷水网站| 亚洲一区高清亚洲精品| 天天一区二区日本电影三级| 婷婷六月久久综合丁香| 国内精品宾馆在线| 国产欧美日韩精品一区二区| 久久国产乱子免费精品| videossex国产| 久久久久久久久久久丰满| 三级男女做爰猛烈吃奶摸视频| 99久久精品热视频| videossex国产| 老师上课跳d突然被开到最大视频| 亚洲成人av在线免费| 99热网站在线观看| 亚洲欧美中文字幕日韩二区| kizo精华| 久久这里只有精品中国| 老司机影院成人| 99热精品在线国产| 久久国产乱子免费精品| 一本久久中文字幕| 亚洲av成人精品一区久久| 亚洲中文字幕一区二区三区有码在线看| 亚洲av免费在线观看| 久久久久网色| 国产精品99久久久久久久久| 丰满的人妻完整版| 能在线免费观看的黄片| 91精品一卡2卡3卡4卡| 亚洲精品日韩在线中文字幕 | 女同久久另类99精品国产91| 哪里可以看免费的av片| 久久久国产成人精品二区| 日本熟妇午夜| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 日本熟妇午夜| 久久精品国产清高在天天线| 午夜激情欧美在线| 婷婷精品国产亚洲av| 97超视频在线观看视频| 精品无人区乱码1区二区| av在线亚洲专区| 哪个播放器可以免费观看大片| 69av精品久久久久久| 久久精品国产亚洲av天美| av天堂中文字幕网| 国国产精品蜜臀av免费| 人妻久久中文字幕网| 亚洲色图av天堂| 免费电影在线观看免费观看| 亚洲国产精品成人久久小说 | 18禁在线无遮挡免费观看视频| 老师上课跳d突然被开到最大视频| 国产免费一级a男人的天堂| a级毛片a级免费在线| 美女高潮的动态| 精品少妇黑人巨大在线播放 | 国产精品久久久久久av不卡| 亚洲av成人精品一区久久| 亚洲精品亚洲一区二区| 亚洲国产精品成人综合色| 中文字幕人妻熟人妻熟丝袜美| 成年女人看的毛片在线观看| 18禁在线播放成人免费| 又黄又爽又刺激的免费视频.| 国产高潮美女av| 综合色av麻豆| 午夜久久久久精精品| 欧美一区二区精品小视频在线| 黄色视频,在线免费观看| 91精品国产九色| 美女脱内裤让男人舔精品视频 | 高清在线视频一区二区三区 | www日本黄色视频网| 嫩草影院入口| 人人妻人人看人人澡| 国产精品无大码| 成人无遮挡网站| 欧美潮喷喷水| 男女下面进入的视频免费午夜| 亚洲va在线va天堂va国产| 久久这里有精品视频免费| 久久99精品国语久久久| 最近手机中文字幕大全| 两性午夜刺激爽爽歪歪视频在线观看| 97超视频在线观看视频| 日本撒尿小便嘘嘘汇集6| 久久久久久久久久久免费av| 久久精品影院6| a级毛色黄片| 老熟妇乱子伦视频在线观看| 综合色av麻豆| 成年av动漫网址| 成年女人看的毛片在线观看| а√天堂www在线а√下载| 国产大屁股一区二区在线视频| 欧美精品一区二区大全| 少妇猛男粗大的猛烈进出视频 | av在线天堂中文字幕| 日本三级黄在线观看| 亚洲中文字幕一区二区三区有码在线看| 国产视频首页在线观看| 老师上课跳d突然被开到最大视频| 久久热精品热| 欧美xxxx黑人xx丫x性爽| 国产高清三级在线| 国产爱豆传媒在线观看| 亚洲一区二区三区色噜噜| av在线天堂中文字幕| 亚洲图色成人| 中文字幕av成人在线电影| 日本av手机在线免费观看| 亚洲国产精品久久男人天堂| 在线观看免费视频日本深夜| 天美传媒精品一区二区| 久久这里有精品视频免费| 亚洲自拍偷在线| 欧美成人a在线观看| 亚洲av成人精品一区久久| av天堂在线播放| 亚洲国产精品成人综合色| 国产日韩欧美在线精品| www日本黄色视频网| 亚洲四区av| 午夜免费男女啪啪视频观看| 亚洲av.av天堂| 午夜爱爱视频在线播放| 在线观看66精品国产| 精品久久久久久久久av| 蜜桃久久精品国产亚洲av| 身体一侧抽搐| 两个人视频免费观看高清| av在线老鸭窝| 日韩欧美精品v在线| 青春草亚洲视频在线观看| 精品一区二区三区视频在线| 日韩,欧美,国产一区二区三区 | 精品免费久久久久久久清纯| 日韩一区二区视频免费看| 亚洲国产精品sss在线观看| 午夜老司机福利剧场| 中文字幕av在线有码专区| 午夜a级毛片| 午夜爱爱视频在线播放| 大香蕉久久网| 亚洲最大成人av| 亚洲不卡免费看| 国内少妇人妻偷人精品xxx网站| 亚洲av一区综合| 精品久久久久久久久久久久久| 又粗又爽又猛毛片免费看| 国产成人a∨麻豆精品| 亚洲av中文字字幕乱码综合| 久久久久久伊人网av| 国产美女午夜福利| 亚洲欧美精品自产自拍| 99久久无色码亚洲精品果冻| 伦精品一区二区三区| 亚洲av免费在线观看| 麻豆av噜噜一区二区三区| 免费观看在线日韩| 欧美xxxx黑人xx丫x性爽| 中国美白少妇内射xxxbb| 亚洲一级一片aⅴ在线观看| 成人鲁丝片一二三区免费| 国产极品精品免费视频能看的| 搡老妇女老女人老熟妇| av专区在线播放| 国产成人91sexporn| 久久精品久久久久久久性| 身体一侧抽搐| 日本一二三区视频观看| 变态另类成人亚洲欧美熟女| 69av精品久久久久久| 久久久久久久久中文| 亚洲综合色惰| 老司机福利观看| 22中文网久久字幕| 久久久色成人| 亚洲国产高清在线一区二区三| 日本爱情动作片www.在线观看| 成人欧美大片| 亚洲精品色激情综合| 亚洲人成网站在线播| 男女啪啪激烈高潮av片| 成人av在线播放网站| 国产高清不卡午夜福利| 99国产精品一区二区蜜桃av| 亚洲在久久综合| 九色成人免费人妻av| 久久精品久久久久久久性| 国产精品1区2区在线观看.| 欧美bdsm另类| 少妇裸体淫交视频免费看高清| 在线观看美女被高潮喷水网站| 亚洲一区高清亚洲精品| 中国美白少妇内射xxxbb| 国内揄拍国产精品人妻在线| 一级毛片我不卡| 精品久久久久久久久久免费视频| 日本黄大片高清| 国产精品爽爽va在线观看网站| a级毛色黄片| 大香蕉久久网| 亚洲av免费在线观看| 久久午夜亚洲精品久久| 99久久精品一区二区三区| 深夜a级毛片| 毛片一级片免费看久久久久| 欧美潮喷喷水| 少妇人妻一区二区三区视频| 国模一区二区三区四区视频| 亚洲av中文字字幕乱码综合| 成人三级黄色视频| 日韩欧美精品免费久久| 99精品在免费线老司机午夜| 亚洲欧美日韩无卡精品| 一级毛片我不卡| 国产激情偷乱视频一区二区| av免费观看日本| 久久久久九九精品影院| 午夜免费男女啪啪视频观看| 一区福利在线观看| 亚洲欧美精品综合久久99| 国产极品天堂在线| 99热全是精品| 男女啪啪激烈高潮av片| 国产麻豆成人av免费视频| 久久精品91蜜桃| 日韩av在线大香蕉| 白带黄色成豆腐渣| 色播亚洲综合网| 免费看光身美女| 国产视频首页在线观看| 久久这里有精品视频免费| 精品午夜福利在线看| 国产精品永久免费网站| 精品国产三级普通话版| 日本在线视频免费播放| 九九久久精品国产亚洲av麻豆| 晚上一个人看的免费电影| 日本与韩国留学比较| 在线免费观看不下载黄p国产| 寂寞人妻少妇视频99o| 亚洲va在线va天堂va国产| 一个人看的www免费观看视频| 日本av手机在线免费观看| 久久亚洲国产成人精品v| 97在线视频观看| 精品国内亚洲2022精品成人| 麻豆国产97在线/欧美| 国产日韩欧美在线精品| 国产成人精品婷婷| av专区在线播放| 久久精品影院6| 精品久久久久久久久久免费视频| 五月玫瑰六月丁香| 精品国内亚洲2022精品成人| 亚洲欧洲日产国产| 人人妻人人看人人澡| 国产一区二区亚洲精品在线观看| 成人一区二区视频在线观看| 欧美一区二区国产精品久久精品| 少妇的逼水好多| 亚洲内射少妇av| 大又大粗又爽又黄少妇毛片口| 国产精品人妻久久久影院| 国产亚洲av片在线观看秒播厂 | 老熟妇乱子伦视频在线观看| 亚洲精品影视一区二区三区av| 国产精品乱码一区二三区的特点| 国产大屁股一区二区在线视频| 菩萨蛮人人尽说江南好唐韦庄 | 欧美潮喷喷水| 久久久精品94久久精品| 国产精品不卡视频一区二区| 综合色丁香网| 人妻制服诱惑在线中文字幕| 欧美激情在线99| 九草在线视频观看| 卡戴珊不雅视频在线播放| 亚洲性久久影院| 偷拍熟女少妇极品色| 91狼人影院| 床上黄色一级片| 成人综合一区亚洲| 久久久久久久午夜电影| kizo精华| 日本爱情动作片www.在线观看| 特级一级黄色大片| 日日干狠狠操夜夜爽| 亚洲最大成人中文| 亚洲高清免费不卡视频| 成人特级黄色片久久久久久久| 亚洲欧美清纯卡通| 国产精品1区2区在线观看.| 欧美成人a在线观看| 黄色日韩在线| 禁无遮挡网站| 欧美区成人在线视频| 最近中文字幕高清免费大全6| 国产久久久一区二区三区| 成人鲁丝片一二三区免费| 欧美高清成人免费视频www| 91aial.com中文字幕在线观看| 国模一区二区三区四区视频| 亚洲,欧美,日韩| 亚洲欧洲日产国产| 在线a可以看的网站| 欧美性猛交黑人性爽| 在线播放国产精品三级| 自拍偷自拍亚洲精品老妇| 午夜爱爱视频在线播放| 亚洲人成网站在线播放欧美日韩| 美女xxoo啪啪120秒动态图| 中文字幕久久专区| 亚洲av男天堂| 午夜福利成人在线免费观看| 2021天堂中文幕一二区在线观| 精品欧美国产一区二区三| 少妇高潮的动态图| 欧美性感艳星| 成人一区二区视频在线观看| 一区福利在线观看| a级毛色黄片| 日韩人妻高清精品专区| 精品欧美国产一区二区三| av福利片在线观看| 国产男人的电影天堂91| 国产乱人偷精品视频| 久久热精品热| 天天躁日日操中文字幕| 国产伦一二天堂av在线观看| 国产男人的电影天堂91| 午夜久久久久精精品| 黄片无遮挡物在线观看| 一个人观看的视频www高清免费观看| 免费人成在线观看视频色| 久久婷婷人人爽人人干人人爱| 51国产日韩欧美| 日本成人三级电影网站| 一级毛片久久久久久久久女| 精品久久久噜噜| 一区福利在线观看| 女人被狂操c到高潮| 神马国产精品三级电影在线观看| 欧美色视频一区免费| 三级经典国产精品| 美女xxoo啪啪120秒动态图| 免费一级毛片在线播放高清视频| 91麻豆精品激情在线观看国产| 99久国产av精品| 人人妻人人澡欧美一区二区| 国产美女午夜福利| 亚洲国产高清在线一区二区三| 69av精品久久久久久| 特级一级黄色大片| 91久久精品国产一区二区三区| 成人综合一区亚洲| 国产精品免费一区二区三区在线| 午夜精品国产一区二区电影 | 亚洲中文字幕日韩| 国产片特级美女逼逼视频| 又粗又爽又猛毛片免费看| 美女高潮的动态| 国内揄拍国产精品人妻在线| 精品久久久噜噜| 淫秽高清视频在线观看| 国产精品国产高清国产av| 能在线免费看毛片的网站| 成人毛片60女人毛片免费| 日韩一区二区三区影片| 一区福利在线观看| 日韩视频在线欧美| 国产精品女同一区二区软件| 国产白丝娇喘喷水9色精品| 禁无遮挡网站| 蜜臀久久99精品久久宅男| 午夜福利成人在线免费观看| 黑人高潮一二区| 波多野结衣巨乳人妻| 搡老妇女老女人老熟妇| 此物有八面人人有两片| 一个人观看的视频www高清免费观看| 夜夜看夜夜爽夜夜摸| 欧美xxxx性猛交bbbb| 黄片无遮挡物在线观看| 精品久久久久久久久久久久久| 天堂√8在线中文| 国产成人午夜福利电影在线观看| 免费人成视频x8x8入口观看| 亚洲一区二区三区色噜噜| 人妻系列 视频| avwww免费| 99热这里只有是精品50| 嫩草影院精品99| 黄色配什么色好看| 一本一本综合久久| 国产精品无大码| av黄色大香蕉| 日韩强制内射视频| 老熟妇乱子伦视频在线观看| 精品一区二区三区视频在线| 日本与韩国留学比较| 99久久人妻综合| 亚洲内射少妇av| 久久精品国产亚洲av涩爱 | 久久久久久九九精品二区国产| 99热这里只有是精品在线观看| 色噜噜av男人的天堂激情| 国产成人a区在线观看| 91久久精品国产一区二区三区| 黄色一级大片看看| 精品一区二区三区视频在线| 18禁黄网站禁片免费观看直播| 国产综合懂色| or卡值多少钱| 伊人久久精品亚洲午夜| 99久国产av精品国产电影| 最好的美女福利视频网| 久久精品国产亚洲av天美| 麻豆国产97在线/欧美| 亚洲最大成人中文| 国产老妇伦熟女老妇高清| 国产精品福利在线免费观看| 国产av一区在线观看免费| 久久精品综合一区二区三区| 久久久久久久久久久免费av| 可以在线观看的亚洲视频| av国产免费在线观看| 国产综合懂色| 久久综合国产亚洲精品| 亚洲成人av在线免费| 三级毛片av免费| 26uuu在线亚洲综合色| 啦啦啦观看免费观看视频高清| 观看美女的网站| 夜夜看夜夜爽夜夜摸| 国产真实乱freesex| 亚洲欧美日韩无卡精品| 国产成人a区在线观看| 久久综合国产亚洲精品| 欧美三级亚洲精品| 婷婷色av中文字幕| 国产男人的电影天堂91| 亚洲最大成人中文| 简卡轻食公司| 亚洲欧美成人综合另类久久久 | 99国产精品一区二区蜜桃av| 亚洲成人久久性| 精品日产1卡2卡| 一边亲一边摸免费视频| 伦理电影大哥的女人| 熟女人妻精品中文字幕| 婷婷六月久久综合丁香| 国产成人午夜福利电影在线观看| 欧美成人一区二区免费高清观看| 国产女主播在线喷水免费视频网站 | 老熟妇乱子伦视频在线观看| 乱系列少妇在线播放| 久久99蜜桃精品久久| 十八禁国产超污无遮挡网站| 国产黄色小视频在线观看| 人体艺术视频欧美日本| 日本撒尿小便嘘嘘汇集6| av免费在线看不卡| 欧美激情国产日韩精品一区| 国产精品嫩草影院av在线观看| av免费观看日本| 99热精品在线国产| 99久久无色码亚洲精品果冻| 亚洲精品亚洲一区二区| 日韩欧美国产在线观看| 我要搜黄色片| 悠悠久久av| 99热这里只有精品一区| 国产精品三级大全| 网址你懂的国产日韩在线| 国产成人aa在线观看| 久久99蜜桃精品久久| 亚洲av电影不卡..在线观看| www.色视频.com| 欧美高清成人免费视频www| 亚洲精品影视一区二区三区av| av.在线天堂| 成人av在线播放网站| av免费观看日本| 99久久精品一区二区三区| 国产精品久久视频播放| 2021天堂中文幕一二区在线观| 国产在线精品亚洲第一网站| 国产黄色小视频在线观看| 边亲边吃奶的免费视频| 色尼玛亚洲综合影院| av在线天堂中文字幕| 亚洲国产高清在线一区二区三| 嘟嘟电影网在线观看| 有码 亚洲区| 三级男女做爰猛烈吃奶摸视频| 久久99热6这里只有精品| 深夜精品福利| 久久鲁丝午夜福利片| 成人性生交大片免费视频hd| 久久精品国产亚洲av香蕉五月| 人人妻人人看人人澡| 最近手机中文字幕大全| 婷婷色综合大香蕉| 日韩一区二区视频免费看| 97热精品久久久久久| 国产真实伦视频高清在线观看| 日韩视频在线欧美| 六月丁香七月| 91aial.com中文字幕在线观看| 成人一区二区视频在线观看| 欧美成人一区二区免费高清观看| 亚洲成人精品中文字幕电影| 啦啦啦韩国在线观看视频| 国产精品久久视频播放| 91aial.com中文字幕在线观看| 成人一区二区视频在线观看| 欧美xxxx黑人xx丫x性爽| 精品人妻视频免费看| 国产探花在线观看一区二区| 久久午夜福利片| 深夜a级毛片| 精品久久久久久久人妻蜜臀av| 男插女下体视频免费在线播放| 日本色播在线视频| 久久久久久久久中文| 99久国产av精品国产电影| 热99在线观看视频| 国产成人一区二区在线| 精品人妻熟女av久视频| 国产精品免费一区二区三区在线| 亚洲国产精品合色在线| 国产日本99.免费观看| 国产高清视频在线观看网站|