• <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
    張平書法作品選
    成年人黄色毛片网站| 亚洲五月色婷婷综合| 亚洲五月色婷婷综合| 欧美人与性动交α欧美精品济南到| 日本精品一区二区三区蜜桃| 一区二区日韩欧美中文字幕| 精品国产超薄肉色丝袜足j| 色哟哟哟哟哟哟| 亚洲精品成人av观看孕妇| 日本五十路高清| 日韩欧美一区视频在线观看| 人人妻人人爽人人添夜夜欢视频| 亚洲精品粉嫩美女一区| 成年人免费黄色播放视频| 国产精品香港三级国产av潘金莲| 久久久水蜜桃国产精品网| 成人亚洲精品一区在线观看| 99riav亚洲国产免费| 国产成人精品无人区| 99香蕉大伊视频| √禁漫天堂资源中文www| 搡老熟女国产l中国老女人| 国产视频一区二区在线看| av中文乱码字幕在线| 久久欧美精品欧美久久欧美| 999久久久精品免费观看国产| 黄网站色视频无遮挡免费观看| 999精品在线视频| 日韩成人在线观看一区二区三区| 亚洲精品一区av在线观看| 国产不卡一卡二| 女性生殖器流出的白浆| 男女之事视频高清在线观看| 久久中文字幕一级| 免费搜索国产男女视频| 99久久国产精品久久久| 最新美女视频免费是黄的| √禁漫天堂资源中文www| 在线国产一区二区在线| 欧美日韩黄片免| 久久久国产欧美日韩av| 999久久久国产精品视频| 在线观看一区二区三区| 精品国产乱码久久久久久男人| 69精品国产乱码久久久| 国产亚洲欧美精品永久| xxxhd国产人妻xxx| 欧美乱码精品一区二区三区| 国产精品99久久99久久久不卡| 天堂俺去俺来也www色官网| 老汉色∧v一级毛片| 国产精品影院久久| 国产精品香港三级国产av潘金莲| 精品久久久久久久毛片微露脸| 母亲3免费完整高清在线观看| 亚洲精品国产精品久久久不卡| 国产主播在线观看一区二区| 伊人久久大香线蕉亚洲五| 日韩中文字幕欧美一区二区| 两性午夜刺激爽爽歪歪视频在线观看 | 91老司机精品| 午夜久久久在线观看| 在线av久久热| 在线av久久热| 首页视频小说图片口味搜索| 免费不卡黄色视频| 亚洲成a人片在线一区二区| 日本vs欧美在线观看视频| a级毛片在线看网站| 国产主播在线观看一区二区| 国产精品久久久av美女十八| 婷婷精品国产亚洲av在线| 国产野战对白在线观看| 久久久久亚洲av毛片大全| 好男人电影高清在线观看| 亚洲 国产 在线| 久久人人97超碰香蕉20202| 老熟妇仑乱视频hdxx| 悠悠久久av| 悠悠久久av| 日韩免费av在线播放| 一进一出抽搐gif免费好疼 | 国产成人精品在线电影| 欧美日韩乱码在线| 十八禁人妻一区二区| 一进一出抽搐动态| 精品午夜福利视频在线观看一区| 精品欧美一区二区三区在线| 又大又爽又粗| 嫩草影院精品99| 99热只有精品国产| 国产成人系列免费观看| 18禁观看日本| 国产亚洲精品久久久久久毛片| 91精品国产国语对白视频| 最新美女视频免费是黄的| 一级片'在线观看视频| 亚洲人成伊人成综合网2020| 69精品国产乱码久久久| 中文字幕人妻熟女乱码| 免费在线观看黄色视频的| 人人妻人人澡人人看| 午夜老司机福利片| 很黄的视频免费| e午夜精品久久久久久久| 欧美丝袜亚洲另类 | 国产精品久久电影中文字幕| 日本撒尿小便嘘嘘汇集6| 男人操女人黄网站| 国产精品久久视频播放| 性色av乱码一区二区三区2| 欧美日本中文国产一区发布| 国产精品国产av在线观看| 搡老岳熟女国产| 免费看a级黄色片| 成人国产一区最新在线观看| 欧美+亚洲+日韩+国产| 欧美+亚洲+日韩+国产| 国产男靠女视频免费网站| 精品久久久久久,| 免费少妇av软件| 在线观看免费视频网站a站| 国产成人系列免费观看| 纯流量卡能插随身wifi吗| 日韩欧美国产一区二区入口| 日本黄色日本黄色录像| 免费看十八禁软件| 高清在线国产一区| 欧美日韩瑟瑟在线播放| 精品一区二区三区视频在线观看免费 | 在线观看一区二区三区激情| 1024香蕉在线观看| 国产精品九九99| 久久精品国产亚洲av高清一级| 一级片免费观看大全| 精品久久久久久成人av| 国产精品国产高清国产av| 亚洲精品在线美女| 色在线成人网| 日韩成人在线观看一区二区三区| 丁香欧美五月| 国产高清激情床上av| 亚洲片人在线观看| 国产精品国产av在线观看| 国产精品免费一区二区三区在线| 99热只有精品国产| 国产黄a三级三级三级人| 91成年电影在线观看| 久久人妻福利社区极品人妻图片| 男女之事视频高清在线观看| 国产成人一区二区三区免费视频网站| 日韩欧美在线二视频| 亚洲专区字幕在线| 91成人精品电影| 亚洲精品美女久久久久99蜜臀| 亚洲五月天丁香| 欧美日韩亚洲国产一区二区在线观看| xxxhd国产人妻xxx| 黑人巨大精品欧美一区二区mp4| 麻豆成人av在线观看| 一级黄色大片毛片| 免费看a级黄色片| 男女床上黄色一级片免费看| 男女下面插进去视频免费观看| 自拍欧美九色日韩亚洲蝌蚪91| 久久久久久久久免费视频了| 一级黄色大片毛片| 成人18禁在线播放| 丁香六月欧美| 亚洲专区中文字幕在线| 免费av毛片视频| 在线天堂中文资源库| 亚洲一卡2卡3卡4卡5卡精品中文| 黄片播放在线免费| 国产欧美日韩一区二区三区在线| 久久亚洲精品不卡| 免费看十八禁软件| 亚洲在线自拍视频| 亚洲精品中文字幕在线视频| 国产成年人精品一区二区 | 久久精品国产综合久久久| 国产高清视频在线播放一区| 一区福利在线观看| 国产成人精品久久二区二区91| 欧美一级毛片孕妇| 亚洲午夜理论影院| 嫩草影视91久久| videosex国产| 精品一区二区三区四区五区乱码| 国产成人免费无遮挡视频| 搡老熟女国产l中国老女人| 国产亚洲av高清不卡| 一级毛片高清免费大全| 可以免费在线观看a视频的电影网站| 人妻丰满熟妇av一区二区三区| 中文字幕精品免费在线观看视频| 久久香蕉精品热| 午夜a级毛片| av超薄肉色丝袜交足视频| 免费av中文字幕在线| 久久香蕉精品热| 日本黄色日本黄色录像| 国产麻豆69| 久久久久久久久中文| 天天躁狠狠躁夜夜躁狠狠躁| 日日干狠狠操夜夜爽| 国产97色在线日韩免费| 两性午夜刺激爽爽歪歪视频在线观看 | 久久久国产一区二区| 国产成+人综合+亚洲专区| 神马国产精品三级电影在线观看 | 一级作爱视频免费观看| 日日爽夜夜爽网站| 18禁国产床啪视频网站| 欧美成人免费av一区二区三区| 欧美成人午夜精品| 香蕉丝袜av| 久久久久久久久免费视频了| 在线观看一区二区三区激情| 18禁黄网站禁片午夜丰满| 老熟妇乱子伦视频在线观看| 国产av精品麻豆| 最近最新中文字幕大全电影3 | 老司机亚洲免费影院| 长腿黑丝高跟| 欧美大码av| 亚洲美女黄片视频| 国产精品偷伦视频观看了| 成人永久免费在线观看视频| 亚洲中文av在线| 国产成人啪精品午夜网站| 三上悠亚av全集在线观看| 人成视频在线观看免费观看| 国产成人啪精品午夜网站| 嫩草影视91久久| 又大又爽又粗| 欧美日韩av久久| 淫秽高清视频在线观看| 成人18禁高潮啪啪吃奶动态图| 久久人人精品亚洲av| 亚洲精华国产精华精| 精品卡一卡二卡四卡免费| 亚洲视频免费观看视频| 老司机午夜福利在线观看视频| 黄色 视频免费看| 亚洲国产精品一区二区三区在线| netflix在线观看网站| 亚洲成人国产一区在线观看| 成年女人毛片免费观看观看9| 老汉色∧v一级毛片| 交换朋友夫妻互换小说| 看黄色毛片网站| 欧美午夜高清在线| 国产精品久久电影中文字幕| 亚洲欧美日韩无卡精品| 免费在线观看亚洲国产| 纯流量卡能插随身wifi吗| 最近最新免费中文字幕在线| 欧美国产精品va在线观看不卡| 国产午夜精品久久久久久| 国产精品永久免费网站| 国产xxxxx性猛交| 国产精品美女特级片免费视频播放器 | 亚洲欧美一区二区三区黑人| 国产欧美日韩综合在线一区二区| 国产精品一区二区免费欧美| 亚洲精品中文字幕在线视频| 国产午夜精品久久久久久| 国产成人一区二区三区免费视频网站| 99久久99久久久精品蜜桃| 亚洲片人在线观看| 亚洲人成电影免费在线| 亚洲欧美一区二区三区久久| 日韩三级视频一区二区三区| 国产在线精品亚洲第一网站| 日韩欧美一区二区三区在线观看| 亚洲欧美精品综合久久99| 久久久久久免费高清国产稀缺| 满18在线观看网站| 中文字幕精品免费在线观看视频| 俄罗斯特黄特色一大片| 久久国产亚洲av麻豆专区| 婷婷精品国产亚洲av在线| 国产亚洲精品第一综合不卡| 亚洲精品成人av观看孕妇| av中文乱码字幕在线| 久久人妻福利社区极品人妻图片| 国产麻豆69| 老司机在亚洲福利影院| 欧美激情 高清一区二区三区| 久久精品成人免费网站| 亚洲国产中文字幕在线视频| 欧美人与性动交α欧美软件| 在线观看免费视频网站a站| 欧美在线一区亚洲| 两个人看的免费小视频| 久久九九热精品免费| 精品国内亚洲2022精品成人| 国产高清视频在线播放一区| 18禁裸乳无遮挡免费网站照片 | 欧美人与性动交α欧美精品济南到| 变态另类成人亚洲欧美熟女 | 亚洲欧美精品综合久久99| 亚洲av成人不卡在线观看播放网| 国产精品美女特级片免费视频播放器 | 国产精品综合久久久久久久免费 | 操美女的视频在线观看| 国产精品影院久久| 国产麻豆69| 午夜精品国产一区二区电影| 精品熟女少妇八av免费久了| 精品国产亚洲在线| 成人手机av| a级毛片黄视频| 久久人妻福利社区极品人妻图片| 国产熟女xx| 国产三级黄色录像| 免费一级毛片在线播放高清视频 | 欧美人与性动交α欧美精品济南到| 黑人巨大精品欧美一区二区mp4| cao死你这个sao货| 老司机亚洲免费影院| 国产精品av久久久久免费| 国产1区2区3区精品| 国产一卡二卡三卡精品| 正在播放国产对白刺激| 1024香蕉在线观看| 久久影院123| 国产精品久久久人人做人人爽| 久久草成人影院| 免费看a级黄色片| 国产成人av教育| 999久久久国产精品视频| 亚洲色图 男人天堂 中文字幕| 在线观看66精品国产| 在线看a的网站| 亚洲专区字幕在线| 中文字幕人妻丝袜制服| 国产熟女xx| 99在线人妻在线中文字幕| 狂野欧美激情性xxxx| 亚洲精品中文字幕在线视频| 亚洲av片天天在线观看| 大型黄色视频在线免费观看| 日韩三级视频一区二区三区| 国产区一区二久久| 欧美 亚洲 国产 日韩一| 久久狼人影院| 国产国语露脸激情在线看| 亚洲中文字幕日韩| 在线看a的网站| 久久青草综合色| 亚洲视频免费观看视频| 精品一区二区三区四区五区乱码| 两性午夜刺激爽爽歪歪视频在线观看 | 手机成人av网站| 欧美激情极品国产一区二区三区| 99久久综合精品五月天人人| 久久婷婷成人综合色麻豆| 久久精品国产亚洲av高清一级| 免费在线观看日本一区| 国产xxxxx性猛交| 日本vs欧美在线观看视频| 97碰自拍视频| 极品人妻少妇av视频| 伊人久久大香线蕉亚洲五| 国产精品美女特级片免费视频播放器 | 色在线成人网| 亚洲欧美精品综合久久99| 国产精品久久久av美女十八| 高清av免费在线| av在线播放免费不卡| 成熟少妇高潮喷水视频| 黄色女人牲交| www.熟女人妻精品国产| 欧美av亚洲av综合av国产av| av福利片在线| 亚洲av五月六月丁香网| 日韩欧美免费精品| 国产1区2区3区精品| 在线免费观看的www视频| 最新在线观看一区二区三区| 99久久久亚洲精品蜜臀av| av超薄肉色丝袜交足视频| 一本综合久久免费| 国产蜜桃级精品一区二区三区| 亚洲人成电影免费在线| 国产亚洲av高清不卡| 国产伦人伦偷精品视频| 精品国产乱子伦一区二区三区| av在线播放免费不卡| x7x7x7水蜜桃| 一区二区三区激情视频| 亚洲欧美一区二区三区黑人| 一区二区三区国产精品乱码| 国产精品一区二区三区四区久久| 蜜桃久久精品国产亚洲av| 久久精品夜夜夜夜夜久久蜜豆| 欧洲精品卡2卡3卡4卡5卡区| 色吧在线观看| 欧美最黄视频在线播放免费| a级毛片免费高清观看在线播放| 在线国产一区二区在线| 嫩草影院新地址| 成人高潮视频无遮挡免费网站| 亚洲欧美日韩卡通动漫| 亚洲人成网站在线播放欧美日韩| 给我免费播放毛片高清在线观看| 久久久久久九九精品二区国产| 无人区码免费观看不卡| 亚洲欧美日韩无卡精品| 热99re8久久精品国产| 久久精品91蜜桃| 国产精品久久久久久人妻精品电影| www.999成人在线观看| 欧美日韩乱码在线| 国产综合懂色| 午夜免费成人在线视频| 国产精品嫩草影院av在线观看 | 色视频www国产| 欧美一区二区精品小视频在线| 亚洲av免费高清在线观看| 在线天堂最新版资源| 亚洲欧美日韩东京热| 最近视频中文字幕2019在线8| 在线国产一区二区在线| 久久久久久国产a免费观看| 久久久久久久亚洲中文字幕 | www.熟女人妻精品国产| 午夜免费成人在线视频| 久久性视频一级片| 国内精品一区二区在线观看| 亚洲无线在线观看| 99热这里只有是精品50| 草草在线视频免费看| 欧美xxxx性猛交bbbb| 欧美日韩亚洲国产一区二区在线观看| 午夜福利高清视频| 嫁个100分男人电影在线观看| 激情在线观看视频在线高清| 天美传媒精品一区二区| 搡老熟女国产l中国老女人| 深夜精品福利| 久久久久精品国产欧美久久久| 欧美bdsm另类| 欧美中文日本在线观看视频| 中文字幕av成人在线电影| 久久久久国内视频| 国产精品不卡视频一区二区 | 午夜精品一区二区三区免费看| 色综合亚洲欧美另类图片| 麻豆av噜噜一区二区三区| 91午夜精品亚洲一区二区三区 | 特大巨黑吊av在线直播| 91麻豆精品激情在线观看国产| 成人美女网站在线观看视频| 午夜福利免费观看在线| 成人高潮视频无遮挡免费网站| 亚洲精品久久国产高清桃花| 久久久久久九九精品二区国产| 国产视频内射| 99久久精品热视频| 国产av麻豆久久久久久久| 综合色av麻豆| 欧美日韩国产亚洲二区| 精品久久久久久久末码| 成人毛片a级毛片在线播放| 美女高潮的动态| 婷婷六月久久综合丁香| 日日夜夜操网爽| 在线观看美女被高潮喷水网站 | 精品免费久久久久久久清纯| 九九在线视频观看精品| 热99re8久久精品国产| 18禁黄网站禁片免费观看直播| 久久久久精品国产欧美久久久| 亚洲内射少妇av| 久久久久久久久久成人| 成人无遮挡网站| 老熟妇仑乱视频hdxx| 成人鲁丝片一二三区免费| 国产伦精品一区二区三区视频9| 夜夜躁狠狠躁天天躁| 舔av片在线| 日韩精品中文字幕看吧| 午夜精品在线福利| 村上凉子中文字幕在线| 中文在线观看免费www的网站| 国产一区二区激情短视频| 欧美绝顶高潮抽搐喷水| 亚洲五月婷婷丁香| 精品人妻1区二区| 欧美另类亚洲清纯唯美| 亚洲av成人av| 久久精品国产亚洲av涩爱 | av中文乱码字幕在线| 欧美高清性xxxxhd video| 亚洲18禁久久av| 91九色精品人成在线观看| 欧美成人a在线观看| 国产精品综合久久久久久久免费| 国产精品久久久久久久久免 | 精品久久久久久成人av| 十八禁网站免费在线| 91狼人影院| 欧美xxxx黑人xx丫x性爽| 久久热精品热| 亚洲av电影不卡..在线观看| 国产欧美日韩精品亚洲av| 国产精品不卡视频一区二区 | 日本免费一区二区三区高清不卡| 免费看日本二区| 日韩欧美三级三区| 好看av亚洲va欧美ⅴa在| 欧美成人一区二区免费高清观看| 日韩人妻高清精品专区| 欧美高清性xxxxhd video| 有码 亚洲区| 亚洲国产精品合色在线| 国产淫片久久久久久久久 | 国内精品久久久久久久电影| 久久这里只有精品中国| 免费一级毛片在线播放高清视频| 丝袜美腿在线中文| 国产爱豆传媒在线观看| 亚洲国产精品久久男人天堂| 国产精品三级大全| 国语自产精品视频在线第100页| av女优亚洲男人天堂| 五月伊人婷婷丁香| 中文字幕久久专区| 69人妻影院| 我的老师免费观看完整版| 特大巨黑吊av在线直播| 欧美午夜高清在线| 国产熟女xx| 午夜福利视频1000在线观看| 少妇人妻一区二区三区视频| 中文字幕久久专区| 丝袜美腿在线中文| 国产高潮美女av| 色视频www国产| 99热只有精品国产| 91在线观看av| 国产亚洲av嫩草精品影院| 久久国产乱子伦精品免费另类| 国产色爽女视频免费观看| 免费观看人在逋| 国产精品一区二区性色av| 欧美成人免费av一区二区三区| 天堂√8在线中文| 99久久九九国产精品国产免费| 欧美日韩福利视频一区二区| 热99在线观看视频| 成人精品一区二区免费| 日韩高清综合在线| 日韩成人在线观看一区二区三区| 怎么达到女性高潮| 久久久久久久精品吃奶| 亚洲精品亚洲一区二区| 中文在线观看免费www的网站| 免费大片18禁| 国产免费一级a男人的天堂| 亚洲人成网站在线播| 午夜福利免费观看在线| 亚洲成人免费电影在线观看| 亚洲精品456在线播放app | 国产成人av教育| 亚洲综合色惰| 国产免费男女视频| 成人av在线播放网站| 九九久久精品国产亚洲av麻豆| 国产精品不卡视频一区二区 | 国产野战对白在线观看| 夜夜躁狠狠躁天天躁| 日韩av在线大香蕉| 久久久久亚洲av毛片大全| 国产高清三级在线| 久久6这里有精品| 嫩草影视91久久| 久久人妻av系列| 欧美成人性av电影在线观看| 久久热精品热| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 中文字幕人妻熟人妻熟丝袜美| 乱人视频在线观看| 身体一侧抽搐| 久久伊人香网站| 18禁在线播放成人免费| 亚洲成人免费电影在线观看| 午夜激情欧美在线| av欧美777| 99热精品在线国产| 欧美成人性av电影在线观看| 动漫黄色视频在线观看| 老司机午夜十八禁免费视频| av国产免费在线观看| 国产精品一区二区免费欧美| 国产精品精品国产色婷婷| 亚洲精品一卡2卡三卡4卡5卡| 精品一区二区三区视频在线| av中文乱码字幕在线| 欧美另类亚洲清纯唯美| 色吧在线观看| 午夜视频国产福利| 51午夜福利影视在线观看| 一个人看的www免费观看视频| 成年人黄色毛片网站| 亚洲av成人不卡在线观看播放网| 校园春色视频在线观看| 88av欧美| 国产午夜精品论理片| 亚洲精品影视一区二区三区av| 999久久久精品免费观看国产| 精品无人区乱码1区二区|