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

    Magnetic ground state of plutonium dioxide: DFT+U calculations

    2023-03-13 09:19:44YueFeiHou侯躍飛WeiJiang江偉ShuJingLi李淑靜ZhenGuoFu付振國andPingZhang張平
    Chinese Physics B 2023年2期
    關(guān)鍵詞:張平振國

    Yue-Fei Hou(侯躍飛) Wei Jiang(江偉) Shu-Jing Li(李淑靜)Zhen-Guo Fu(付振國) and Ping Zhang(張平)

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

    2Beijing University of Chemical Technology,Beijing 100029,China

    3School of Physics and Physical Engineering,Qufu Normal University,Qufu 273165,China

    Keywords: strongly correlated system,magnetic ground state,noncollinear,MagGene

    1.Introduction

    Some of the actinide oxides (XO2,X= U, Np, Pu,Am,etc.), as typical strongly correlated systems, normally form the cubic fluorite structure with anFmˉ3m(225) space group.[1]These series of compounds are either Mott insulators or charge-transfer insulators and they provide outstanding platforms to study the properties of localized 5f electrons.Now that the 5f orbits are partially filled in these systems,nonzero magnetic moments can be naturally imagined.Meanwhile,with the presence of spin-orbit coupling(SOC)and cubic crystal field effect,the magnetism of these actinide oxides is expected to be even abundant and complex.

    Years before, the magnetic ground states of these actinide oxides had been a fundamental topic causing much discussion.[2]Today experimental evidences of the triple-k(3k)-type magnetic ground states for UO2and NpO2have been reported in literatures.[3-9]However, different from the Curie-Weiss susceptibility found in UO2[10]and NpO2,[11]a temperature-independent susceptibility[12]is observed in PuO2.What is more,the nuclear magnetic resonance(NMR)experiments[13-15]imply no magnetic order occurred in PuO2at low temperature.Surprisingly, a simple crystal field(CF) theory successfully explains the non-magnetic picture of PuO2.[2]In response,inelastic neutron scattering spectrum(INSS)confirms the transition from theΓ1singlet state to theΓ4triplet state (i.e., the CF ground state to the first excited state) of PuO2with a peak at 123 meV.[16]Based on a local density approximation(LDA)calculation to density functional theory(DFT),[17]theΓ1state is 99 meV lower than theΓ4state.Another LDA calculation combined with exact diagonalization(LDA+ED)of the Anderson impurity model[18]suggests a gap of 126 meV between the two states.These theoretical results are considered to be in reasonable agreement with the INSS.Until here, the nonmagnetic ground state of PuO2seems to be clear.However, one abnormal issue is that the experimental susceptibility is much smaller(about 50%)than that of one would expect from the van Vleck coupling of theΓ1with theΓ4.[2,12,16,18]Although the multiconfigurational wave function calculations[19]reproduces the experimental susceptibility of PuO2, more discussions are required to ensure this achievement.Anyway, the magnetism of PuO2needs more investigations from both experimental measurements and theoretical modeling.

    In fact, PuO2is a challenging system to investigate in experiments duo to its nucleonic radioactive decay and the inhomogeneity of the samples.[20-23]As a result, a series of computational works have been carried out to understand the ground state properties of this system in-depth.Before further simulation of its properties, a reasonable magnetic model for the ground state of PuO2is firstly required.For instance, by assuming a simple nonmagnetic (NM), ferromagnetic (FM) or collinear antiferromagnetic (AFM) ground state,many fundamental physical characters of PuO2,including electronic,[24-28]optical,[27,29]thermodynamic,[27,30]lattice dynamic,[24,31]and surface[24,32-34]properties,have been studied through first-principles calculations.Among them,DFT-based calculations are especially popular to study PuO2due to its balance between accuracy and economy.

    It has been revealed that conventional DFT or even DFT+Ucalculations tend to provide rather an FM or a collinear AFM ground state[24-28]than the recognized nonmagnetic state by comparing their free energies.To theoretically obtain a nonmagnetic ground state, the modified LDA calculation[18]seems to succeed while CF parameters obtained from the experiment are applied which makes the computational work not totallyab initio.On the other hand,it may also be meaningful to leave out the energy difference between magnetic or nonmagnetic states and directly impose a nonmagnetic state in the calculation.For instance, by taking into account SOC and HubbardUmodifications,DFT+Ucalculations for the nonmagnetic PuO2are performed[35,36]and useful results are gained.Despite all this, the consistency between the real ground state and the simulated ground state is of great significance to the applicability of DFT.Only with this robust consistency could we then safely predict other similar systems by using the DFT based methodology.Hence,the ground state of PuO2belonging to the DFT framework is required to be determined and examined carefully.Pegget al.[20]calculated the total energy of PuO2with different magnetic structures and found a noncollinear longitudinal 3kAFM ground state by using HSE06 hybrid density functional.Nonmagnetic state matches much higher energy than magnetic states in their calculations.The origin of the inconsistency occurred in a series of DFT-based studies[24-28]needs more investigations.

    In this work,a systematic DFT+Ucalculation for determining the magnetic ground state of PuO2is performed.The SOC effect is included explicitly in our calculations, which has been proved to be an essential effect in PuO2.[35,36]Therefore, the noncollinear magnetic algorithm is employed here.HubbardUmodel is used to describe the Coulomb interactions between correlated electrons.The whole process of this work can also be adopted to study other actinide oxides,such as UO2, NpO2, AmO2,etc.The local density approximation (LDA),[37]the generalized gradient approximation(PBE here),[38]and the strongly constrained and appropriately normed(SCAN)[39]exchange-correlation functionals are applied in our calculations, respectively.The three exchangecorrelation functionals which contain different physical information are belonging to three different levels of the Jacob ladder.[40]A pureab initioDFT+Ucalculation without taking into account any experimental parameters is performed here.

    2.Computational method

    We employ the Viennaab initiosimulation package(VASP)[41]code to achieve DFT+Uscheme.The projected augmented wave (PAW) method[42-44]with the frozen-core approximation is used.The first Brillouin zone of the unit cell is sampled by 6×6×6k-point mesh.The cut-off energy of the plane wave basis for electron wave function is chosen as 520 eV.We have carefully checked the convergence ofk-point mesh and cut-off energy.Atomic structures are fully relaxed until the Feynman-Hellman force on each atom is less than 0.01 eV/°A.The convergence accuracy of the energy for selfconsistent calculation is 10-7eV.

    In the DFT+Ustudy, screened Coulomb parameterUand Hund exchange parameterJare significant.Many efforts have been devoted to determine theUandJin PuO2.By scanning lattice constant,[24,26,27]bulk modulus,[24,26]band gap,[20,27]and magnetic moment[20]in variousUandJ, the best combination of the two parameters can be chosen to match properties observed in experiments.Here, we adopt a scheme from a first-principles calculation[45]which points out the significance of electronic correlation on O-p orbits in actinide oxides.We note that a nonzeroUfor O-p does not change the properties of PuO2qualitatively.The electronic correlation on O-p orbits only slightly modifies the electronic structure (see supporting information S5).Thereby,U=4.7 eV andJ=0.6 eV for Pu-5f orbits andU=3.0 eV andJ=0 eV for O-p orbits are applied.Multiple values ofUfor Pu-f electrons have been used in the literature varying from 4 eV to 6.5 eV.[18,20,27,34-36]As a precaution,a test before further discussions is made and we find that the energy order of the magnetic structures keeps unchanged whenUvaries in the range of 3 eV to 7 eV.An unreasonable metallic ground state which is contrast to the experimental insulating state occurs only whenUis quite small(roughly less than 0.5 eV).

    We make use of the approach originated from Liechtenstein’s[46]DFT+Umethod for the main part of the calculation.The generalized form of the energy functional is given by

    whereρσ(r)is the charge density,andnσis the trace of density matrix for spinσ(↑or↓).The first term on the right-hand side represents standard DFT functional.It can be substituted by LDA, PBE, and SCAN functionals in our work.UandJaffect the total energy independently in the double counting term written as

    wheren=n↑+n↓denotes the total occupation number of both spin-up and spin-down electrons.In addition,Dudarev’s simplified DFT+Uapproach is also employed as a comparison,which will be discussed in Subsection 3.4.

    3.Results

    3.1.Magnetic structures in the cubic fluorite atomic structure

    The Bravais lattice of PuO2is face centered cubic(FCC)with one Pu atom and two O atoms in the primitive cell.In order to produce possible 3kmagnetic structures,the conventional cell (Pu4O8) of PuO2with the simple cubic (SC) lattice is adopted, as shown in the right top of Fig.1.The four first-nearest Pu ions in this unit cell form a regular tetrahedron before relaxing the atomic structure with an imposed magnetic structure.Because of the periodic boundary condition,all second-nearest Pu ions are naturally in FM order.We cautiously do a test and prove that the unit cell adopted here is convincing to describe the magnetic ground state of PuO2(see more details in supporting information S1).The magnetic moments of O ions are ignored.

    Based on former ideas,[2,20]typical magnetic structures are listed in Fig.1.The superiority of our work is that three more collinear magnetic structures as shown in Figs.1(e),1(f),and 1(j) are involved in the calculations.The three magnetic structures are the transformations of transverse 1kAFM or longitudinal 1kAFM structures and they are discussed for the first time.With the existence of strong SOC and remarkable CF effect,[2]the magnetic anisotropy is very likely to play an important role in the formation of the magnetic ground state.As a result,we classify the magnetic structures into 001,011,and 111 structures whose magnetic moments point along〈001〉,〈011〉, and〈111〉, respectively.Details are shown in supporting information S2.

    Fig.1.The atomic structure of cubic fluorite PuO2 (on the right top) and typical collinear and noncollinear magnetic structures [(a)-(l)].Pu atoms represented by the large balls are shown to form magnetic structures.Arrows on the atoms represent magnetic moment vectors.Each color of the arrows corresponds to one specific vector of each magnetic structure.Orientations of the magnetic moments are along(a)-(c)〈001〉,(d)-(h)〈011〉,and(i)-(l)〈111〉.The magnetic structures shown in panels (a), (d), and (i) are FM structures while the others are the so-called AFM structures with the label“AFM”omitted for simplicity.

    3.2.Determination of the magnetic ground state

    The atomic structure of actinide oxides tends to produce Jahn-Teller (JT) distortion when imposing a magnetic structure.For instance, the realistic UO2single crystal lowers the symmetry toPaˉ3 space group at low temperature beneath the N′eel temperature.The distortion on O atoms of about 0.014 °A along〈111〉has been observed within a transverse 3kmagnetic structure.[3]This feature can also be faithfully reflected by our calculations performed on PuO2.The space groups of the distorted crystal with different imposed magnetic structures are listed in Table 1.One can find that the nonmagnetic state and the longitudinal 3kmagnetic state remainFmˉ3mspace group recommended by the experiments.The transverse 3kmagnetic state of PuO2is ofPaˉ3 symmetry.The distortion of O atoms is carefully calculated by employing the three functionals to be around 0.015 °A along〈111〉.This phenomenon is very similar to that in UO2.Meanwhile, imposing this distortion is far not enough to change the intrinsic magnetic ground state of PuO2.Details can be found in supporting information S3.It indicates that the possible atomic distortion in actinide oxides at low temperature is rather a consequence than an motive force of the magnetic structure for the ground state.Dudarevet al.have studied the JT distortion of UO2through LSDA+Ucalculations before.[47]Their conclusions and ours agree well with each other.

    Lattice constant of PuO2calculated by LDA, PBE,SCAN,and HSE06[20]are compared in Fig.2(a).One can see that different imposed magnetic structures only slightly modify the lattice constant.The notable distinction is contributed from the exchange-correlation functional.On one hand,LDA underestimates the lattice constant, while PBE (GGA) overestimates the lattice constant.This is consistent with prior LDA+Uand GGA+Uresults.[26]On the other hand,HSE06 and especially SCAN give out pretty good prediction of the lattice constant.

    The total magnetic moment of a Pu ion is composed of both spin moment and orbital moment originated from 5f electrons.Large magnetic dipole moment of the Pu ion is recommended as shown in Fig.2(b).The magnetic moments of the Pu ions for different magnetic structures are approximately invariant in an acceptable range.This can be thought as a consequence of both strong SOC and on-site Coulomb interaction.These two effects dominate the local electronic configuration of the magnetic ions and tend to keep the magnetic moment unchanged.When rotating the local moments, only the electron-electron interactions between different sites are notably affected.This kind of interactions is relatively weak.As a result, different magnetic structures hardly make huge difference on the local magnetic moments of the Pu ions.

    Table 1.Space group of PuO2 with the magnetic structures imposed.These symmetry properties are only based on the simulated atomic structures.Any atomic distortion less than 10-3 °A is neglected.LDA,PBE,and SCAN as well as HSE06[20] give the consistent results here.

    Fig.2.(a)Lattice constant of PuO2 with different magnetic structures calculated by LDA(red squares), PBE(blue triangles), SCAN(green circles),and HSE06(black crosses).[20]In the abscissa,letters a-l correspond to magnetic structures shown in Figs.1(a)-1(l),and m represents the nonmagnetic situation.Black dashed line denoting the experimental lattice constant from x-ray diffraction.[54] (b)The calculated magnetic moment of a single Pu ion with different magnetic structures imposed.Magnetic structures belonging to 001,011,and 111 types are separated in three regions.(c)Relative energy of PuO2 with different magnetic structures.Nonmagnetic state is not shown as it matches too much higher energy.The total energy of the magnetic ground state for each exchange-correlation functional is set to 0.(d)The lowest energy found in each generation relative to the energy of the magnetic ground state during the MagGene process.

    The total energy of the Pu4O8unit cell varies in a range of about 420 meV among the magnetic structures as shown in Fig.2(c).The energy of the nonmagnetic state in our calculations is almost about 4 eV higher than the other magnetic states no matter which exchange-correlation functional is applied (LDA, PBE, SCAN, as well as HSE06[20]).Thereby,we do not discuss the simulated nonmagnetic states in this work.According to the results for LDA, PBE, and SCAN,the orientation of the local magnetic moments dominate the total energy as we expected.The 111 structures are normally favored in energy and the 011 structures take higher energy than the 111 structures.Finally,the 001 structures cause much higher energy than both 111 structures and 011 structures.Furthermore, the FM configuration(i.e., structures a, d, and i)is also less stable than AFM configuration in a group of structures whose magnetic moments point along the same crystallographic direction family (for〈001〉,〈011〉, and〈111〉, respectively).In the case of HSE06,[20]things are a little different.Referring to Fig.2(c),relative energy of the magnetic structures shows a similar trend to our applied three functionals.However,we notice that the departure in energy between FM structures and AFM structures is larger.From this point of view,HSE06 hybrid functional may take more highly of the interactions from f-p super-exchange and f-f exchange than the other three functionals.Anyway,the longitudinal 3kAFM structure to be the magnetic ground state of PuO2is confirmed by LDA,PBE,SCAN,and HSE06 consistently.

    When SOC is involved in DFT+Ucalculations,total energy is not only the functional of the charge density but also the functional of the magnetization density when no external field is imposed.[48]The magnetic moment of an ion on siteiis the integral of the magnetization density inside the Wigner-Seitz(WS)radius given by

    In our unconstrained DFT+U+SOC calculations realized by VASP,[41]the orientation of the initial local magnetic moment is still needed to be artificially set.Usually, to search the ground state in the subspace of magnetic ordered states,several typical high-symmetry magnetic structures are at least needed to be considered (refer to Fig.1).In this way, the number of considered magnetic structures is limited.To exclude any other possible magnetic structures with lower energy which might be missed,we expand the searching range to a larger space.We firstly let the local moments in the unit cell rotate freely.Random magnetic structures can be then gained at different moments.Meanwhile, to overcome the blindness during the calculations, the MagGene code based on the genetic algorithm is applied to locate the magnetic ground state efficiently.Details can be found in supporting information S4 and also the paper from Zhenget al.[49]At the beginning of the searching process for PuO2, sixteen random magnetic structures are generated.Following that, five structures with the lowest five energies are chosen as parents for the next generation.Based on these parents, six mutation structures and two crossover structures are produced.The last part of the next generation is eight (16-6-2=8) random structures,which help to get rid of the local minima of the energy surface and keep the size of population unchanged.The lowest energy found in each generation is shown in Fig.2(d).When a thirty generations’ searching is finished, the longitudinal 3kAFM ground state is found in the 19th, 22nd,and 23rd generations when LDA,PBE,and SCAN are applied,respectively.

    3.3.Electronic structure of the magnetic ground state

    A large unit cell of PuO2is necessary to describe the longitudinal 3kAFM ground state.The first Brillouin zone of the Pu4O8simple cubic lattice is shown in Fig.3.

    Fig.3.First Brillouin zone and high symmetry path for the simple cubic lattice of the Pu4O8 conventional cell.

    We adoptΓ-X-M-Γ-R-X/M-Rhigh symmetry path[50]to describe the dispersion of the total energy.Density of states(DOS) and electronic band structures calculated by the three functionals are plotted in Fig.4.Our DFT+U+SOC calculations for the longitudinal 3kAFM ground state reasonably suggest PuO2to be a charge-transfer insulator which has been widely identified before.A strong hybridization between Pu-f and O-p orbits is identified.This feature explains why taking into account the correlation on O-p orbits is necessary.Similar to the magnetic ground state calculated by HSE06,[20]a remarkable peak of O-p below the Fermi level in the DOS is located at about-1.1 eV.Another broad peak close to the Fermi level is actually a double peak in our calculation as shown in Figs.4(a)-4(c).One can see from the Figs.4(d)-4(f) that the direct gap of the band structure is 2.46 eV, 2.47 eV, and 2.69 eV corresponding to the LDA,PBE,and SCAN functionals.To our surprise, these values calculated from the longitudinal 3kAFM ground state are even more close to the experimental optical gap which is 2.8 eV[51]than that from the simulated nonmagnetic state(about 1.8 eV).[18,36]Hence, we trust that the longitudinal 3kAFM structure is the best approximation to model the electronic ground state of PuO2so far and it has been adopted to calculate the formation energy.[47]

    Fig.4.Electronic DOS of the longitudinal 3k AFM ground state and the corresponding band structures for LDA[(a)and(d)],PBE[(b)and(e)],and SCAN[(c)and(f)].Blue circles and red circles in panels(d)-(f)represent the contributions from Pu-f and O-p orbitals.The size of the circles shows the weight of each component.

    3.4.A comparison with Dudarev’s DFT+U approach

    Another widely used DFT+Uscheme based on Dudarev’s DFT+Umodel[52]is inspected here.The total energy modified by this model is given by

    whereρσis the density matrix of correlated electrons with spinσ.In this simplified approach,UandJare handled in a spherically averaged form of the screened Coulomb interaction and exchange interaction.We adopt the sameUandJas used in the former Liechtenstein’s model.Total energies of different magnetic structures are calculated again,and the corresponding results are shown in Fig.5(a).We firstly find that the energy range among these structures is about 200 meV,which is only about half of that obtained by Liechtenstein’s approach.For LDA or PBE functional,the〈111〉easy axis of magnetization is reproduced and the longitudinal 3kmagnetic ground state is determined as expected.However, for SCAN functional, structures with atomic moments along〈011〉are favored in energy.

    Fig.5.(a)Relative energy of PuO2 with different magnetic structures obtained from Dudarev’s DFT+U approach.In the abscissa,letters a-l correspond to the magnetic structures shown in Figs.1(a)-1(l).The total energy of the magnetic ground state for each exchange-correlation functional is set to 0.(b)The(1-1 0)section of PuO2 crystal with longitudinal 3k AFM structure.Red arrows represent magnetic moments.(c)The charge density plot on the(1-1 0)section calculated by Liechtenstein’s DFT+U approach and Dudarev’s DFT+U approach.(d)The magnetization densities plot on the(1-1 0)section calculated by Liechtenstein’s DFT+U approach and Dudarev’s DFT+U approach.The direction of the magnetization vector on each grid is not shown for simplicity.

    Moreover, the preference for AFM structures in energy is invalid in this case.As shown in Fig.5(a), the 001-FM structure is the most stable among 001 magnetic structures.In Dudarev’s DFT+Uapproach, it is abnormal for SCAN to get such results which are contrary to LDA and PBE, especially when a consistent conclusion is arrived with Liechtenstein’s DFT+Uapproach for all three exchange-correlation functionals.

    We also would like to compare the ionic states of the longitudinal 3kmagnetic structure calculated by Liechtenstein’s DFT+Uapproach and Dudarev’s DFT+Uapproach.The charge densities of both approaches with SCAN adopted are represented in Fig.5(c).There is only slight difference in the two charge densities, which is originated from Pu ions.However,the magnetization density emphasizes a remarkable distinction.By calculating the integral of the magnetization density inside the Wigner-Seitz radius of Pu, we obtain the magnetic moment to be 3.00μBand 3.72μBfrom Liechtenstein’s approach and Dudarev’s approach.From the distributions of the magnetization densities as shown in Fig.5(d),we notice that two different ionic states are gained from these two DFT+Uapproaches although the same values ofUandJare employed.The corresponding occupation matrices for Pu-f orbits are also shown in supporting information S6 to numerically confirm this character.

    4.Discussion

    We treat PuO2as a standard computational model of actinide oxides with half-filled 5f electrons and expect the methodology here to be expanded to other actinide oxides.Previous DFT calculations have illustrated that a simple collinear magnetic model without SOC is meaningful to describe a series of properties of PuO2as long as a HubbardUmodification is considered.However,this framework of computation is too rough for magnetic properties.With the existence of the on-site CoulombUand strong SOC, the way of occupation for Pu-5f orbits takes on great responsibilities for the total energy.On our calculated NM state, the electronic configuration of the 5f electrons remarkably differs from the magnetic states and that’s why it has much higher energy.Based on the same local occupation for Pu-5f orbits, the total charge density and the magnetization density can be once again modified by a specific magnetic structure.According to our simulations,the total energies slightly fluctuate among different magnetic structures in a range of dozens of meV per PuO2.

    We point out that the typical magnetic structures here combined with their energies can be used to fit an effective pseudo spin Hamiltonian as proposed before.[53]Since the spin-spin interactions in actinide oxides are complex, multiple model parameters have been taken into account in the Hamiltonian.As a result, for the accuracy of the fitting job,the amount of the used magnetic structures should not be too small.Meanwhile, the magnetic structures ought to be with the approximately equal local electronic occupation for Pu-5f orbits.Once any magnetic structures with definitely different ionic states are involved in the fitting job, the fitted parameters would be unphysical and meaningless.At least for PuO2,the magnetic structures simulated by our DFT+Umethod are perfect samples for a high-quality fitting.

    Furthermore,there are at least two possible reasons need to be noticed for understanding the inconsistency of the ground states between DFT calculation and experiments for PuO2.The first reason is that the longitudinal 3kmagnetic state obtained here may be a metastable state which corresponds to a local minima on the energy surface.To validate this, an occupation control procedure for the local electrons should be a promising method to completely search the subspace of the electronic configuration.Another reason is that the physical information of the computational model may be not compatible with the real system.If so, additional modification is required for DFT based calculations to PuO2.Anyway,the origin of this inconsistency needs more intriguing investigations.An universalab initiocomputational framework is needed to reasonably describe the real ground state of actinide oxides based on DFT+Umethod.We will leave this issue to our long-term work.

    5.Conclusions

    We studied the magnetic ordered states of PuO2by using DFT+Umethod with SOC included.The LDA, PBE, and SCAN exchange-correlation functionals are applied, respectively.The typical magnetic structures in the cubic fluorite crystal structure are classified and labeled.The JT distortions are simulated, which can be seen as an appendage of the imposed magnetic structure.The orientations of the local moments in the crystal environment are found to play an important role in the contribution to the total energy.Beyond the usual searching procedure of the magnetic ground state, the MagGene code based on the genetic algorithm is employed to ensure the longitudinal 3kAFM ground state recommended by all the three functionals.After that, the electronic structure of the magnetic ground state is discussed.According to our calculations,PuO2is predicted to be a charge-transfer insulator, whose optical gap agrees well with the experimental value.We believe that the longitudinal 3kAFM structure should be a good approximation to model the magnetic ground state of PuO2.All these results obtained here are in reasonable agreement with that obtained by using a hybrid functional in the prior computational work.Finally, we notice a divergence between Liechtenstein’s DFT+Uapproach and Dudarev’s DFT+Uapproach when both are applied to PuO2.This phenomenon occurs only when SCAN exchangecorrelation functional is adopted.

    Acknowledgements

    Project supported by National Natural Science Foundation of China,(Grant No.12104034).

    We thank professor Fa-Wei Zheng from the School of Physics of Beijing Institute of Technology for his advice to this work.We also thank professor James T.Pegg from Department of Chemistry of University College London for his useful discussion with us.

    猜你喜歡
    張平振國
    Topological properties of Sb(111)surface: A first-principles study
    嘰嘰喳喳的小喜鵲
    愛在拉薩
    這是你爺倆
    金秋(2020年16期)2020-12-09 01:41:50
    我和繼父13年
    我和繼父13年
    文苑·感悟(2019年12期)2019-12-23 07:24:46
    我和繼父13 年
    文苑(2019年23期)2019-12-05 06:50:22
    Enhanced spin-dependent thermopower in a double-quantum-dot sandwiched between two-dimensional electron gases?
    我和繼父的13年
    37°女人(2019年6期)2019-06-10 08:48:11
    張平書法作品選
    菩萨蛮人人尽说江南好唐韦庄 | 精品午夜福利视频在线观看一区| 人人妻人人看人人澡| 亚洲成人av在线免费| 深夜a级毛片| 精品久久久久久久久av| 日韩亚洲欧美综合| 久久精品国产亚洲av天美| 国产伦一二天堂av在线观看| 老熟妇仑乱视频hdxx| 日韩三级伦理在线观看| 国产精华一区二区三区| 色av中文字幕| 色综合色国产| 一区福利在线观看| 综合色丁香网| 又爽又黄无遮挡网站| 可以在线观看毛片的网站| 女同久久另类99精品国产91| 亚洲综合色惰| 亚洲精品久久国产高清桃花| 人妻久久中文字幕网| 看黄色毛片网站| 色哟哟哟哟哟哟| 久久久久久伊人网av| 亚洲国产精品成人综合色| 中文字幕精品亚洲无线码一区| 亚洲人成网站高清观看| 亚洲,欧美,日韩| 成人高潮视频无遮挡免费网站| 国产老妇女一区| 色播亚洲综合网| 少妇人妻精品综合一区二区 | 免费人成在线观看视频色| 亚洲成人av在线免费| 热99在线观看视频| 色播亚洲综合网| 精品日产1卡2卡| 久久久久免费精品人妻一区二区| 日本与韩国留学比较| 女生性感内裤真人,穿戴方法视频| 俺也久久电影网| 欧美高清性xxxxhd video| 国产精品99久久久久久久久| 深夜a级毛片| 深爱激情五月婷婷| 成人二区视频| 99精品在免费线老司机午夜| 亚洲丝袜综合中文字幕| 日韩欧美在线乱码| 两性午夜刺激爽爽歪歪视频在线观看| 国产白丝娇喘喷水9色精品| 国产熟女欧美一区二区| 天天躁夜夜躁狠狠久久av| 九九久久精品国产亚洲av麻豆| 国产精品三级大全| or卡值多少钱| 一区二区三区免费毛片| 免费在线观看成人毛片| 精品久久久久久成人av| 亚洲在线自拍视频| 日本熟妇午夜| av天堂中文字幕网| 人妻夜夜爽99麻豆av| 91午夜精品亚洲一区二区三区| 亚洲精品一区av在线观看| 成人三级黄色视频| 97人妻精品一区二区三区麻豆| 亚洲欧美成人综合另类久久久 | 99久国产av精品| 日本一本二区三区精品| 精品国内亚洲2022精品成人| 伦精品一区二区三区| 人人妻人人澡人人爽人人夜夜 | a级毛片a级免费在线| 我的女老师完整版在线观看| 少妇被粗大猛烈的视频| 色在线成人网| 亚洲精品日韩av片在线观看| 久久亚洲精品不卡| 国产一区二区亚洲精品在线观看| 两个人视频免费观看高清| 日韩亚洲欧美综合| 美女 人体艺术 gogo| 好男人在线观看高清免费视频| 欧美xxxx黑人xx丫x性爽| 春色校园在线视频观看| 天美传媒精品一区二区| 欧美在线一区亚洲| 欧美在线一区亚洲| videossex国产| 69av精品久久久久久| 欧美激情在线99| 日韩高清综合在线| 毛片一级片免费看久久久久| 亚洲七黄色美女视频| 亚洲四区av| 久久人人爽人人片av| 久久人人爽人人片av| 高清毛片免费看| 欧美日本视频| 亚洲中文字幕一区二区三区有码在线看| 欧美日韩精品成人综合77777| 精品乱码久久久久久99久播| 国产精品爽爽va在线观看网站| 少妇的逼好多水| 国产黄片美女视频| 欧美区成人在线视频| 色哟哟哟哟哟哟| 婷婷精品国产亚洲av在线| 日日啪夜夜撸| 国产精品免费一区二区三区在线| 亚洲国产欧洲综合997久久,| av视频在线观看入口| 女生性感内裤真人,穿戴方法视频| 国产成人a∨麻豆精品| 亚洲国产精品成人综合色| 久久精品久久久久久噜噜老黄 | 久久精品国产亚洲av天美| 又黄又爽又刺激的免费视频.| 国产成人freesex在线 | 久久精品国产鲁丝片午夜精品| 欧美+日韩+精品| 嫩草影视91久久| 男女边吃奶边做爰视频| 夜夜夜夜夜久久久久| 日韩欧美在线乱码| 午夜福利18| 狂野欧美激情性xxxx在线观看| 国产麻豆成人av免费视频| 国内精品久久久久精免费| 麻豆av噜噜一区二区三区| 精品久久久噜噜| 身体一侧抽搐| 少妇裸体淫交视频免费看高清| 免费在线观看成人毛片| 亚洲,欧美,日韩| 91av网一区二区| 欧美在线一区亚洲| 国产精品一区二区免费欧美| 日韩av不卡免费在线播放| 噜噜噜噜噜久久久久久91| 麻豆av噜噜一区二区三区| 国产毛片a区久久久久| 成人国产麻豆网| 亚洲高清免费不卡视频| 美女被艹到高潮喷水动态| 此物有八面人人有两片| 99久久无色码亚洲精品果冻| 精品人妻一区二区三区麻豆 | av在线播放精品| 国产91av在线免费观看| 精品一区二区三区视频在线观看免费| 国国产精品蜜臀av免费| 亚洲av不卡在线观看| 人人妻人人澡人人爽人人夜夜 | 亚洲图色成人| 日本免费一区二区三区高清不卡| 两个人视频免费观看高清| 欧美成人免费av一区二区三区| 黄色日韩在线| 男女边吃奶边做爰视频| 国产久久久一区二区三区| 亚洲成人av在线免费| 亚洲18禁久久av| 三级毛片av免费| 熟妇人妻久久中文字幕3abv| 国产成年人精品一区二区| 久99久视频精品免费| 搡老熟女国产l中国老女人| 欧美bdsm另类| 国产真实乱freesex| 欧美日韩国产亚洲二区| 一进一出好大好爽视频| 桃色一区二区三区在线观看| 全区人妻精品视频| 校园人妻丝袜中文字幕| 亚洲欧美日韩高清在线视频| eeuss影院久久| 国产成人影院久久av| 可以在线观看毛片的网站| 国产三级在线视频| 午夜福利视频1000在线观看| 少妇熟女aⅴ在线视频| 一进一出好大好爽视频| 欧美中文日本在线观看视频| 永久网站在线| 在线观看美女被高潮喷水网站| 内地一区二区视频在线| 亚洲激情五月婷婷啪啪| 久久天躁狠狠躁夜夜2o2o| 免费观看精品视频网站| 亚洲美女视频黄频| 韩国av在线不卡| 亚洲av五月六月丁香网| 欧美一区二区国产精品久久精品| 日本三级黄在线观看| 不卡一级毛片| 欧美在线一区亚洲| 91在线观看av| 亚洲精品色激情综合| 少妇熟女欧美另类| 成人特级黄色片久久久久久久| 国产精品日韩av在线免费观看| 人妻久久中文字幕网| 日韩一区二区视频免费看| 亚洲成av人片在线播放无| 成年女人永久免费观看视频| 亚洲国产色片| 大又大粗又爽又黄少妇毛片口| 成人亚洲欧美一区二区av| 日本黄色视频三级网站网址| 亚洲av中文av极速乱| 久久人人精品亚洲av| 午夜影院日韩av| 亚洲内射少妇av| 久99久视频精品免费| 精品久久久噜噜| 白带黄色成豆腐渣| 人人妻,人人澡人人爽秒播| 99热网站在线观看| 有码 亚洲区| 久久国产乱子免费精品| 99国产极品粉嫩在线观看| 高清午夜精品一区二区三区 | 老师上课跳d突然被开到最大视频| 久久久久久久久久久丰满| 国产精品不卡视频一区二区| 日韩人妻高清精品专区| 3wmmmm亚洲av在线观看| 精品一区二区免费观看| 午夜a级毛片| 国产精品免费一区二区三区在线| 91精品国产九色| 中文字幕av在线有码专区| 美女免费视频网站| 亚洲欧美日韩高清专用| 久久精品国产亚洲av涩爱 | 麻豆国产97在线/欧美| 老司机午夜福利在线观看视频| 成人永久免费在线观看视频| 免费不卡的大黄色大毛片视频在线观看 | 一进一出抽搐动态| 欧美高清成人免费视频www| 国产av在哪里看| 成人毛片a级毛片在线播放| 国产欧美日韩精品亚洲av| 精品99又大又爽又粗少妇毛片| 亚洲av五月六月丁香网| 露出奶头的视频| 天堂影院成人在线观看| 老熟妇仑乱视频hdxx| 少妇猛男粗大的猛烈进出视频 | 亚洲人成网站在线播| 午夜福利高清视频| 舔av片在线| 久久久精品大字幕| 97碰自拍视频| 99久久中文字幕三级久久日本| 性欧美人与动物交配| 91久久精品电影网| 乱人视频在线观看| 小蜜桃在线观看免费完整版高清| eeuss影院久久| 热99re8久久精品国产| 国产亚洲精品av在线| 此物有八面人人有两片| 可以在线观看的亚洲视频| 午夜免费男女啪啪视频观看 | 俺也久久电影网| 亚洲在线自拍视频| 日日啪夜夜撸| 亚洲熟妇中文字幕五十中出| 欧美极品一区二区三区四区| 免费观看在线日韩| 亚洲最大成人av| 欧美xxxx性猛交bbbb| 成人特级黄色片久久久久久久| 有码 亚洲区| 久久久久性生活片| 久久久久国产网址| 九九热线精品视视频播放| 久久99热6这里只有精品| 不卡视频在线观看欧美| 免费看美女性在线毛片视频| 亚洲人成网站在线播| 亚洲国产高清在线一区二区三| 美女黄网站色视频| 春色校园在线视频观看| 久久草成人影院| 国产亚洲精品综合一区在线观看| 一级黄色大片毛片| 久久人人精品亚洲av| 日日摸夜夜添夜夜爱| 久久人妻av系列| 成人av一区二区三区在线看| 亚洲av美国av| 校园春色视频在线观看| 欧美一区二区精品小视频在线| 色综合站精品国产| 亚洲丝袜综合中文字幕| 无遮挡黄片免费观看| 欧美xxxx黑人xx丫x性爽| 国产大屁股一区二区在线视频| 97人妻精品一区二区三区麻豆| 毛片一级片免费看久久久久| av专区在线播放| 午夜影院日韩av| 三级毛片av免费| 国产白丝娇喘喷水9色精品| 香蕉av资源在线| 日本欧美国产在线视频| 真人做人爱边吃奶动态| 国产成年人精品一区二区| 午夜爱爱视频在线播放| 日韩国内少妇激情av| 噜噜噜噜噜久久久久久91| 极品教师在线视频| 一区二区三区高清视频在线| 午夜亚洲福利在线播放| 乱码一卡2卡4卡精品| 3wmmmm亚洲av在线观看| 亚洲国产色片| 国内精品一区二区在线观看| 一进一出好大好爽视频| 亚洲图色成人| 久久热精品热| 精品久久久久久久末码| 黄色日韩在线| 97在线视频观看| 欧美成人精品欧美一级黄| 一个人看的www免费观看视频| 欧美日韩乱码在线| 联通29元200g的流量卡| 国产一区二区在线观看日韩| 欧美成人免费av一区二区三区| 高清午夜精品一区二区三区 | 一级毛片久久久久久久久女| 老女人水多毛片| 国产精品av视频在线免费观看| 给我免费播放毛片高清在线观看| 搞女人的毛片| 久久中文看片网| 内地一区二区视频在线| 大型黄色视频在线免费观看| 91在线精品国自产拍蜜月| 成人美女网站在线观看视频| 国产成年人精品一区二区| 国产日本99.免费观看| 国产午夜精品论理片| 免费av毛片视频| 日本欧美国产在线视频| 日本-黄色视频高清免费观看| 亚洲无线在线观看| 欧美国产日韩亚洲一区| 天天躁夜夜躁狠狠久久av| 乱系列少妇在线播放| 国产探花极品一区二区| АⅤ资源中文在线天堂| 色视频www国产| 欧美日本亚洲视频在线播放| 日本一二三区视频观看| 亚洲无线观看免费| 两个人视频免费观看高清| 欧美一级a爱片免费观看看| av中文乱码字幕在线| 亚洲婷婷狠狠爱综合网| 成人特级黄色片久久久久久久| 色综合站精品国产| 亚洲精品乱码久久久v下载方式| 午夜福利在线观看免费完整高清在 | 一区二区三区高清视频在线| 午夜爱爱视频在线播放| 久久欧美精品欧美久久欧美| 在线看三级毛片| 亚洲精品亚洲一区二区| 色视频www国产| 村上凉子中文字幕在线| 男人狂女人下面高潮的视频| 男人舔女人下体高潮全视频| 特级一级黄色大片| 哪里可以看免费的av片| 免费一级毛片在线播放高清视频| 嫩草影院精品99| 精品久久久久久久人妻蜜臀av| 亚洲熟妇熟女久久| 精品久久久久久久人妻蜜臀av| 亚洲无线在线观看| 国产免费男女视频| 女人十人毛片免费观看3o分钟| 色5月婷婷丁香| 综合色丁香网| 国产黄片美女视频| 俺也久久电影网| 亚洲第一区二区三区不卡| 国产精品嫩草影院av在线观看| 不卡一级毛片| 美女 人体艺术 gogo| 亚洲真实伦在线观看| 麻豆精品久久久久久蜜桃| 日产精品乱码卡一卡2卡三| 男女视频在线观看网站免费| 免费av不卡在线播放| 俺也久久电影网| 不卡一级毛片| 午夜福利视频1000在线观看| 一区二区三区高清视频在线| 国产人妻一区二区三区在| 在线a可以看的网站| 亚洲自偷自拍三级| 亚洲人成网站高清观看| 欧美一区二区国产精品久久精品| 国产精品亚洲美女久久久| 欧美日韩乱码在线| 亚洲国产精品久久男人天堂| 亚洲精品成人久久久久久| 高清毛片免费看| 一级a爱片免费观看的视频| 亚洲人成网站在线播| 97碰自拍视频| 欧美丝袜亚洲另类| 国产视频一区二区在线看| 免费观看的影片在线观看| 成人高潮视频无遮挡免费网站| 午夜亚洲福利在线播放| 91狼人影院| а√天堂www在线а√下载| 久久亚洲国产成人精品v| 狠狠狠狠99中文字幕| 搡老熟女国产l中国老女人| 亚洲av第一区精品v没综合| 成人欧美大片| 国产一区二区在线av高清观看| 国产av一区在线观看免费| 一区福利在线观看| 久久人妻av系列| 97人妻精品一区二区三区麻豆| 免费黄网站久久成人精品| 欧美又色又爽又黄视频| 蜜桃亚洲精品一区二区三区| 老熟妇仑乱视频hdxx| 99热这里只有是精品在线观看| 大型黄色视频在线免费观看| 极品教师在线视频| 搡老岳熟女国产| 国产精品久久电影中文字幕| 少妇丰满av| 十八禁国产超污无遮挡网站| 国产精品一及| 成人综合一区亚洲| 亚洲av免费在线观看| 啦啦啦韩国在线观看视频| 国产一区亚洲一区在线观看| 欧美一级a爱片免费观看看| 亚洲精品国产成人久久av| 久久久久久久亚洲中文字幕| 99在线人妻在线中文字幕| 成人综合一区亚洲| 亚洲性夜色夜夜综合| 99久久精品一区二区三区| 黄色视频,在线免费观看| 免费看av在线观看网站| 婷婷六月久久综合丁香| 午夜福利在线观看吧| 变态另类成人亚洲欧美熟女| videossex国产| 日本黄色视频三级网站网址| 国产综合懂色| 美女黄网站色视频| 99久国产av精品| 欧美日韩在线观看h| 极品教师在线视频| 亚洲精品成人久久久久久| 中文在线观看免费www的网站| 久久99热这里只有精品18| 黄色欧美视频在线观看| 亚洲av成人av| 欧美日韩国产亚洲二区| 女人十人毛片免费观看3o分钟| 欧美日韩一区二区视频在线观看视频在线 | 亚洲av中文字字幕乱码综合| 美女xxoo啪啪120秒动态图| 午夜免费男女啪啪视频观看 | 内射极品少妇av片p| 听说在线观看完整版免费高清| 天堂av国产一区二区熟女人妻| 婷婷色综合大香蕉| 国产v大片淫在线免费观看| 亚洲一级一片aⅴ在线观看| 国产av一区在线观看免费| 日产精品乱码卡一卡2卡三| 九九爱精品视频在线观看| 国产精品av视频在线免费观看| 大又大粗又爽又黄少妇毛片口| 免费高清视频大片| 国产成人freesex在线 | 尾随美女入室| 国模一区二区三区四区视频| 22中文网久久字幕| 国产片特级美女逼逼视频| 日韩欧美精品v在线| 干丝袜人妻中文字幕| 色吧在线观看| 中出人妻视频一区二区| 久久亚洲精品不卡| 啦啦啦观看免费观看视频高清| 五月玫瑰六月丁香| 色5月婷婷丁香| 国产精品爽爽va在线观看网站| 网址你懂的国产日韩在线| 天天躁日日操中文字幕| 91午夜精品亚洲一区二区三区| 精品无人区乱码1区二区| 一区福利在线观看| 直男gayav资源| a级毛片免费高清观看在线播放| 久久久欧美国产精品| 久久久久久久午夜电影| 亚洲欧美日韩东京热| 欧美一区二区亚洲| 国产精品久久久久久久电影| 99久久精品国产国产毛片| 国产三级在线视频| 又爽又黄无遮挡网站| 亚洲精品色激情综合| 精品国内亚洲2022精品成人| 国产高清有码在线观看视频| 亚洲精品一卡2卡三卡4卡5卡| 亚洲人与动物交配视频| 国产精品久久久久久亚洲av鲁大| 最近2019中文字幕mv第一页| 99久久精品热视频| 亚洲av成人精品一区久久| 国产精品女同一区二区软件| 久久精品国产亚洲av涩爱 | 高清毛片免费观看视频网站| 日韩一本色道免费dvd| 亚洲人与动物交配视频| 国产女主播在线喷水免费视频网站 | videossex国产| 亚洲精品亚洲一区二区| 99热全是精品| 日本与韩国留学比较| 午夜久久久久精精品| 亚洲一级一片aⅴ在线观看| 国产伦一二天堂av在线观看| 亚洲aⅴ乱码一区二区在线播放| 91久久精品电影网| 麻豆精品久久久久久蜜桃| 久久中文看片网| 精品久久久噜噜| 亚洲欧美日韩卡通动漫| 国产黄片美女视频| 亚洲成人久久爱视频| 色在线成人网| 尤物成人国产欧美一区二区三区| 97超视频在线观看视频| 亚洲av中文字字幕乱码综合| 欧美成人一区二区免费高清观看| 亚洲电影在线观看av| 啦啦啦啦在线视频资源| 亚洲成人中文字幕在线播放| 亚洲国产精品国产精品| 国产精品福利在线免费观看| 国产精品一二三区在线看| 国产黄色视频一区二区在线观看 | 精品国内亚洲2022精品成人| 国产高清视频在线观看网站| 成人综合一区亚洲| 欧美激情在线99| 亚洲精品一卡2卡三卡4卡5卡| 在线看三级毛片| 久久国产乱子免费精品| av专区在线播放| 欧美一区二区国产精品久久精品| 亚洲欧美成人精品一区二区| 国产精华一区二区三区| 1000部很黄的大片| 别揉我奶头 嗯啊视频| 久久中文看片网| 少妇熟女aⅴ在线视频| 一级毛片我不卡| 一级a爱片免费观看的视频| 99久久成人亚洲精品观看| 超碰av人人做人人爽久久| 黄色日韩在线| 国产真实伦视频高清在线观看| 最新在线观看一区二区三区| 国产欧美日韩精品亚洲av| 亚洲经典国产精华液单| 韩国av在线不卡| 又粗又爽又猛毛片免费看| 在线免费观看不下载黄p国产| 日日撸夜夜添| 国产伦精品一区二区三区四那| 97热精品久久久久久| 身体一侧抽搐| 亚洲图色成人| 成人漫画全彩无遮挡| 久久鲁丝午夜福利片| 成人av在线播放网站| 91麻豆精品激情在线观看国产| 两性午夜刺激爽爽歪歪视频在线观看| 99久久无色码亚洲精品果冻| 欧美日本亚洲视频在线播放| 免费人成在线观看视频色| 国产精品久久视频播放| 久久热精品热| 日韩欧美国产在线观看| 波多野结衣高清作品| 日韩,欧美,国产一区二区三区 | 长腿黑丝高跟| 午夜免费男女啪啪视频观看 | 国产伦一二天堂av在线观看| 亚洲国产色片|