O2,CO2,and H2O Chemisorption on UN(001)Surface:Density Functional Theory Study

Ru-song Li,Bin He,Fei Wng,Xu Peng,Hu Wng

a.Xi’an Research Institute of Hi-Tech,Hongqing Town,Xi’an 710025,China

b.Type Office of China Academy of Engineering Physics,Mianyang 621900,China

O2,CO2,and H2O Chemisorption on UN(001)Surface:Density Functional Theory Study

Ru-song Lia?,Bin Hea,Fei Wanga,Xu Penga,Hua Wangb

a.Xi’an Research Institute of Hi-Tech,Hongqing Town,Xi’an 710025,China

b.Type Office of China Academy of Engineering Physics,Mianyang 621900,China

We performed density functional theory calculations of O2,CO2,and H2O chemisorption on the UN(001)surface using the generalized gradient approximation and PW91 exchangecorrelation functional at non-spin polarized level with the periodic slab model.Chemisorption energies vs.molecular distance from UN(001)surface were optimized for four symmetrical chemisorption sites.The results showed that the bridge parallel,hollow parallel and bridge hydrogen-up adsorption sites were the most stable site for O2,CO2,and H2O molecular with chemisorption energies of 14.48,4.492,and 5.85 kJ/mol,respectively.From the point of adsorbent(the UN(001)surface),interaction of O2with the UN(001)surface was of the maximum magnitude,then CO2and H2O,indicating that these interactions were associated with structures of the adsorbate.O2chemisorption caused N atoms on the surface to migrate into the bulk,however CO2and H2O had a moderate and negligible effect on the surface,respectively.Calculated electronic density of states demonstrated the electronic charge transfer between s,p orbital in chemisorption molecular and U6d,U5f orbital.

Chemisorption,Density functional theory,Geometric relaxation,Electronic density of state

I.INTRODUCTION

Uranium mononitride(UN)is considered as a promising fuel for the fast nuclear generation IV reactors[1, 2].Compared to many uranium and plutonium oxide nuclear fuel,UN has several advantages[2-7],such as higher melting point and smaller lattice constant,which indicate that UN has stronger anti-corrosion capability.Many research groups have identif i ed the important role of UN in anti-corrosion application since the 1960s [8-11],and performed several experimental and theoretical study.The bulk properties of actinide nitrides have been investigated,especially the elastic and magnetic properties.However,contrary to a large number of available experimental data[12-19],there were few reports related to the f i rst principle calculations on uranium nitrides,particularly chemisorption behaviors of molecular or atoms on the surface of uranium nitrides. Only recently,some groups have simulated the reactivity of molecular/atoms with the surface of uranium nitrides[20-22].These reports indicated that molecule O2would spontaneously dissociate after chemisorption on the(001)surface of UN,then the produced O atoms exhibited a strong chemisorption behavior.To deeply understand molecular O2,CO2and H2O chemisorption on the(001)surface of UN,we apply the density functional theory(DFT)method to perform electronic structure calculations in this work.

II.METHODOLOGY

Uranium mononitride(UN)has Nacl-type structure (fcc)with lattice constant a=0.4889 nm,and interatomic dist√ance√ between U and N atoms is 0.248 nm. A periodicthree-layer slab supercell model and single-sided chemisorption mode(a molecule was placed on one side of the slab model,namely 0.5 ML adsorptivity)were used to study molecular chemisorption on the UN(001)surface in all calculations.A vacuum layer of 2.0 nm was added to the unit cell of the layers.

In this work,we used the generalized gradient approximation(GGA)within the framework of density functional theory(DFT)and PW91 exchange-correlation functional[23-25]with the periodic slab model to simulate the chemisorption for gaseous molecular on the UN(001)surface.The outer fourteen electrons (6s26p65f66d17s2)of U were treated as valence electrons and the remaining seventy-eight electrons were treated as core.DFT semi-core pseudopotentials(DSPP)and a double numerical basis set with polarization functions (DNP)have been used to treat core electrons and valence electrons,respectively.All electron basis sets were used for C,H and O atoms.3×3×1 Monkhorst-Pack kpoint meshes were applied in the Brillouin zone(BZ).Aplane-wave cutof fenergy Ecut=500 eV was chosen.The convergence of self-consistent f i eld(SCF)was less than 1.0×10-5eV/atom.Nomagnetic configuration was appropriate for the light actinide U element from the point of total energy,so U5f electrons were in the delocalized 5f3electronic configuration in this work.

Single molecule,one per unit cell,was allowed to approach the UN(001)surface along four different symmetrical positions,namely(i)directly on top of a U atom(U-top position);(ii)directly on top of a N atom (N-top position);(iii)on the middle of two nearest neighbor U atoms(bridge position);(iv)adsorption molecule saw a U atom located on the layer directly below the surface hollow site(hollow position).The chemisorption energy ECwas optimized with respect to the height R of the chemisorbed molecule above the surface,and was given by[26]:

where E(M)is the total energy of the bare UN(001) slab,E(X)is the total energy of the isolated molecule, and E(M+X)is the total energy of the molecule chemisorbed on the surface.

The relative change for bond length was used to describe the change for U-N bond after molecule chemisorption,and was given by

where Riis the bond length for the central N atom on the UN(001)surface and the ith U atom,R0is the U-N bond length before chemisorption,and?denoted the relative change for Ri.

III.SURFACE CONFIGURATION FOR THE UN(001) SURFACE

To test validity of the computational parameters,we fi rst cleaved the UN(001)surface,then checked convergence of the total energy of the UN(001)surface with di ff erent vacuum thicknesses.We considered the vacuum thickness test to be convergent as long as the total energy was less than 10 meV,and the calculation result is shown in Fig.1.The result indicated that the total energy of the system was convergent when the vacuum thickness was larger than 1.8 nm.Therefore,we added a vacuum layer of 2.0 nm onto the unit cell of the slabs in order to reduce the in fl uence of boundary condition on the calculation,and a model for the UN(001)surface is plotted in Fig.2.

The interactions of the surface atoms with UN matrix atoms would be unstable due to the absence of adjacent atom,which was contrary to the matrix atoms.Meanwhile,these non-equilibrium interactions might cause the surface atoms to relax,reconstruct, fi nd a new equilibrium site after cleaving,and fi nally lower the total energy of the surface system.Moreover,this relaxation behavior might change U-N bond length.A configuration model for relaxation calculation is shown in Fig.3, and the results are listed in Table I.

FIG.1 Convergence test of the vacuum thickness for the UN(001)surface.

FIG.2 A calculation model for the UN(001)surface.

In this work,relaxation was def i ned as the relative change of U-N bond length.We f i xed the lowest-layer atoms during relaxation calculation for the UN(001) surface,and presented in Table I.The total energy for this system reduced by about 0.739 eV.However,the relative relaxation was ignorable(the maximum relative 1.109%),so we f i xed the atoms in two low-lying layers to calculate in the following section,otherwise particular declaration.

IV.CHEMISORPTION BEHAVIOR OF GASEOUS MOLECULAR ON THE UN(001)SURFACE

The energy minimum principle demonstrates that the higher the system symmetry,the lower the system energy,and the more stable the system.Therefore,we preferred the gaseous molecular to be chemisorbed onto the high symmetrical position in the crystal surface. UN crystal has face-centered cubic(fcc)structure,andseveral high symmetrical chemisorption positions exist on the surface,we considered four representative positions in this work,namely bridge(B),hollow(H),U-top(U)and N-top(N),as shown in Fig.4.The central atom or the geometrical centre of the chemisorbed molecule was directly placed on the top of individual position to study the chemisorption behavior of molecule on these positions.Chemisorption parameters included not only position,but also chemisorption orientation and height of the chemisorbed molecule to the surface. The chemisorption orientation was associated with the molecule structure,while the chemisorption height depended on stability of the system configuration.

TABLEI Calculated results for configuration relaxation. Ri(i=1,2,in nm)and dij(i=1,2,j=1,2,in nm)represent U-N bond length in the intralayer and interlayer,respectively.

FIG.3 A configuration model for relaxation calculation.Ri(i=1,2)and dij(i=1,2,j=1,2)represent U-N bond length in the intralayer and interlayer,respectively.

A.Chemisorption behavior of molecule O2on the UN(001)surface

As we all know,the chemisorption behavior of molecule O2on the UN(001)surface crystal would be crucial to understand the anti-corrosion mechanisms for UN compound.Experiment result showed that O2is a linear molecule,and O-O bond length is 0.1209 nm. We considered two chemisorption modes for molecule O2approaching the UN(001)surface to investigate O2chemisorption behavior:parallel(P)and vertical(V). The P configuration calculation for U-top chemisorption position was not fulfilled in one single cell.Therefore,we only considered the V configuration for U-top chemisorption position to save the computational resource.

FIG.4 Diagram of four symmetrical chemisorption positions for the UN(001)surface:bridge,hollow,U-top,and N-top.

FIG.5 System configuration of O2molecule on UN(001) surface(a)before and(b)after chemisorption.

For different chemisorption height from the UN(001) surface,the total energy of O2-UN(001)system would change.Chemisorption height was def i ned as the nearest distance for molecule O2from the UN(001)surface, and expressed by the fractional coordinates.The results for the total energy of O2-UN(001)system are listed in Table II.We optimized the system configuration with the minimum energies,the results for O2chemisorption energies,O-O bond lengths,and Mulliken charges of individual configuration are presented in Table III.We could see that the most stable chemisorption position for O2was BP.The total energy of this system would decrease after optimization.As depicted in Fig.5,compared to the non-optimized configurations,optimization induced molecule O2to migrate towards U atom,while N atom moved into the bulk.We labeled the atoms in the calculation configuration to describe the relaxation behavior,and the atom sequence was marked in Fig.5(b).The distance of O from the ith O was used to describe the relative displacement for the atoms in the calculated configuration.The result showed that O-O bond length increased,and O2obviously moved outwards.However,N on the UN(001)surface migratedinto the bulk,as shown in Table IV.To further understand the interaction of O2with the UN(001)surface, we analyzed the projected density of states(PDOSs) before and after O2chemisorption in term of the electronic structure calculations,and shown in Fig.6.

TABLE II Total energies of O2-UN(001)system.h denotes the chemisorption height for O2,and the BP,BV,HP,HV, NP,NV and UV configurations represent approaching the bridge position in parallel manner,the bridge position in vertical manner,the hollow position in parallel manner,the hollow position in vertical manner,the N-top position in parallel manner, the N-top position in vertical manner,the U-top position in vertical manner,respectively.

TABLE III Chemisorption configurations,chemisorption energies Echemisorption(in kJ/mol),O-O bond lengths dO-Oand Mulliken charges Q for O2.

TABLE IV Relative displacements for the atoms in the optimized configuration.

O2s and O2p orbital PDOSs shifted towards lower energy band after O2chemisorption,indicating that the red shift effect occurred in O2frequency,as shown in Fig.6(a).We did not consider other s,p orbitals in U atom because U6d and U5f orbitals dominate the electronic properties of U atom.PDOSs of U6d and U5f orbitals before and after chemisorption are depicted in Fig.6(b).After chemisorption,peak of U6d PDOS apparently widened,which showed that U6d orbital had a very strong interaction with the substrate.Peak value of U5f PDOS decreased from 247.36 electron/eV to 221.25 electron/eV,the peak position shifted from 0.35 eV to 0.546 eV,and the peak area of U5f PDOS diminished,showing that 5f orbital lost electrons.A new f state appeared in the energy range of-1.5 eV to -0.5 eV,and a new d state emerged in the energy range from-5 eV to-4 eV,which implied that O2s or O2p electrons contributed to U6d and U5f orbital,as shown in Fig.6(b).Figure 6(c)and(d)depicted U6d,U5f,O2s orbital and U6d,U5f,O2p orbital PDOSs for O atom being the nearest neighbor of U atom,respectively.As shown in Fig.6(c),O2s orbital hybridized with U6d orbital,and produced a very small peak.While O2s orbital PDOS did not obviously overlapped with U5f orbital.However,O2p orbital and U6d,U5f orbital formed an unambiguous hybridization peak,as shown in Fig.6(d),which was in sharp contrast with Fig.6(c). The larger the overlapping area of PDOS,the higher the hybridized bonding.Therefore,the electronic charge of U6d orbital transferred to O2s and O2p orbital(mainly 2p orbital),while U5f orbital transferred to O2p orbital.

B.Chemisorption behavior of CO2on the UN(001) surface

FIG.6 Projected density of states(PDOSs)of O2s,O2p orbital(a),U6d,U5f orbital(b),O2s,U6d,U5f orbital(c),and O2p,U6d,U5f orbital(d)before and after O2chemisorption on the BP position of the UN(001)surface.The Fermi energy stands at 0 eV.

TABLE V Calculated total energies of CO2-UN(001)system.h denoted the chemisorption height for CO2.

As we all know,CO2is a linear molecule,and the experimental C-O bond length is 0.116 nm.We considered the same approaching mode for CO2as O2,namely P and V modes.The central atom or the geometrical centre of the chemisorbed CO2was directly placed on the top of individual chemisorption positions.We fi rst calculated the total energy of CO2-UN(001)system through seven di ff erent con fi gurations,then chose the chemisorption height h for the lowest total energy as the optimal height for the chemisorption con fi guration,and the calculated results of the total energy are presented in Table V, fi nally,we optimized these con fi gurations to search the optimal chemisorption position.

The geometrical parameters,chemisorption energies,bond lengths,and Mulliken charges for CO2chemisorbed on different positions on the UN(001)surface are listed in Table VI.According to the minimum energy principle,the configuration with the maximum chemisorption energy was the most stable chemisorption configuration.Therefore,the most stable configuration for CO2chemisorbed on the UN(001)surface was HP,as shown in Table VI.To further understand the interaction of CO2with the UN(001)surface,we an-alyzed PDOSs before and after CO2chemisorption in term of the electronic structure calculations,and shown in Fig.7.

FIG.7 PDOSs of(a)C/O 2s,2p orbital,(b)U6d,U5f orbital,(c)O2s,U6d,U5f orbital,and(d)O2p,U6d,U5f orbital before and after CO2chemisorption on the HP position of the UN(001)surface.The Fermi energy stands at 0 eV.

TABLE VIChemisorption configurations,chemisorption energy Echemisorption(in kJ/mol),O-O bond lengths dO-Oand Mulliken charge Q for molecule CO2.

C/O 2s,2p orbital PDOSs shifted towards lower energy band after O2chemisorption,implying that the red shift effect occurred in CO2frequency,as shown in Fig.7(a).We also did not consider other s,p orbital in U because U6d and U5f orbital dominated the electronic properties of U atom.PDOSs of U6d and U5f orbital before and after chemisorption were depicted in Fig.7(b).After chemisorption,the peak value of U6d PDOS in the energy range of 0-1 eV decreased from 17.4 electron/eV to 13.9 electron/eV,while the peak value of PDOS in the energy range of-5 eV to-4 eV decreased from 17.0 electrons/eV to 16.0 electrons/eV, showing that U6d orbital lost electrons,as shown in Fig.7(b).Figure 7(c)and(d)depicted U6d,U5f,O2s orbital and U6d,U5f,O2p orbital PDOSs for O atom being the nearest neighbor of U atom,respectively.As shown in Fig.7(c),O2s orbital and U6d orbital formed a very small hybridization peak,and O2s orbital PDOS did not obviously overlap with U5f orbital.However, O2p orbital and U6d,U5f orbital produced an unambiguous hybridization peak,as shown in Fig.7(d),which was in sharp contrast with Fig.7(c).The larger the overlapping area of PDOS,the higher the hybridized bonding.Therefore,the electronic charge of U6d orbital transferred to O2s and O2p orbital(mainly O2p orbital),while U5f orbital transferred to O2p orbital, these results were similar to O2chemisorption.

C.Chemisorption behavior of molecule H2O on the UN(001)surface

Previous report showed that H2O chemisorbed on U metal would result in formation of UO2passivation film on the metal surface,and this film could prevent U metal from further oxidization[7].Therefore,investigation of H2O chemisorption on the UN surface would be helpful for deeply understanding the anti-corrosion mechanisms of U metal because UN passivation film wasalso a corrosion-resistant material.

TABLE VII Calculated total energies of H2O-UN(001)system.h denotes the chemisorption height for CO2.BU,BD,HU, HD,NU and ND configurations represent bridge H-up,bridge H-down,hollow H-up,hollow H-down,N-top H-up and N-top H-down,respectively.

FIG.8 System configuration before(a)and after(b) chemisorption.

As we all know,H2O is a planar V-type molecule,and the experimental H-O bond length is 0.10 nm,∠HOH bond angle is 104.5?.We considered three symmetrical positions for studying H2O chemisorption on the UN(001)surface:H-up,H-down and H-par,where H-up and H-down denoted that H2O plane was vertical with the UN(001)surface from the upper side and the lower side,respectively,while H-par indicated that H2O plane was parallel with the UN(001)surface.For different chemisorption height from the UN(001)surface, the total energy of H2O-UN(001)system would change. According to the minimum energy principle,the total energy of a system would be the minimum value for H2O chemisorption on the optimal position,and the system would be the most stable one.

A chemisorption height h was def i ned as the nearest distance for the central O atom in H2O from the UN(001)surface.The results for the total energy of H2O-UN(001)system with different h are listed in Table VII.We optimized the systems with the minimum energies,and calculated results for H2O chemisorption energies,H-O bond lengths,∠HOH bond angles and Mulliken charges of individual configuration are listed in Table VIII.From Table VIII we could see that the most stable chemisorption position for H2O was BU(maximum chemisorption energy 5.85 kJ/mol),so H in H2O was vertically chemisorbed on the UN(001)surface from the upper side,as shown in Fig.8(b).Z coordinate of U atom increased by about 0.1%,indicating that H2O chemisorption had a negligible effect on the UN(001) surface.However,H2O not only migrated,but also rotated around O atom.∠HOH bond angle increased by about 5?,while U-O bond length decreased from 0.2932 nm to 0.2634 nm.We labeled every atom in H2O to observe the position changes of H2O before and after optimization,as shown in Fig.8.Calculated results of the relative change for X,Y,and Z coordinates were presented in Table IX,where?X,?Y,and?Z denoted the relative change for X,Y,and Z coordinates, respectively.

To further understand the interaction of H2O with the UN(001)surface,we analyzed PDOSs before and after H2O chemisorption in term of the electronic structure calculations,and shown in Fig.9.

FIG.9 Projected density of states(PDOSs)of(a)H1s,O2s,O2p orbital,(b)U6d and U5f orbital,(c)O2s,U6d,and U5f orbital and(d)O2p,U6d,U5f orbital before and after H2O chemisorption on the BU position of the UN(001)surface.The Fermi energy stands at 0 eV.

TABLE VIII Chemisorption configurations,chemisorption energy Echemisorption,O-O bond lengths dO-O,bond angles, and Mulliken charges Q for H2O.

TABLE IX Relative changes for X,Y and Z coordinates of H and O atoms in the optimized H2O

Figure 9(a)depicted H1s,O2s,and O2p PDOSs before and after chemisorption.H1s,O2s,O2p PDOSs shifted towards lower energy band after H2O chemisorption,indicating that the red shift effect occurred in H2O frequency,as shown in Fig.9(a).PDOSs of U6d and U5f orbital before and after chemisorption are plotted in Fig.9(b).Chemisorption had a negligible effect on the shapes of PDOSs of U6d and U5f orbital,however,the peak area of U6d PDOS in the energy range of-1 and 1 eV obviously diminished,showing that U6d orbital lost electrons,as shown in Fig.9(b).The peak position of U5f orbital shifted from 0.35 eV to 0.468 eV,showing that U5f orbital also lost electrons.O2s PDOS did not obviously overlap with U6d and U5f orbital,as shown in Fig.9(c).While in the energy range of-7 eV to-6 eV, -5 eV to-4 eV,O 2p orbital and U6d orbital formed a very small hybridization peak,as shown in Fig.9(d). Therefore,O atom in H2O and U atom on the UN(001) surface formed a covalent bond,and a small amount of U6d orbital electrons transferred to O2p orbital.Thesmaller the overlapping area,the weaker the chemical bonding,and the lower the transferred charges,which was in consistent with the Mulliken charge analysis.

V.CONCLUSION

We performed the density functional theory calculations of O2,CO2,and H2O chemisorption on the UN(001)surface using the generalized gradient approximation(GGA)and PW91 exchange-correlation functional at non-spin polarized level with the periodic slab model,the results showed that(i)bridge parallel(BP), hollow parallel(HP)and bridge hydrogen-up(BU)adsorption sites were the most stable site for O2,CO2and H2O,respectively.(ii)O2s and O2p orbital PDOSs shifted towards lower energy band after O2chemisorption,and the red shift effect occurred in O2frequency. O2p orbital and U6d,U5f orbital formed an unambiguous hybrid peak,and the electronic charge of U6d and 5f orbital mainly transferred to O2p.(iii)C/O 2s,2p orbital PDOSs shifted towards lower energy band after O2chemisorption,and the red shift effect also occurred in CO2frequency.O2s orbital and U6d orbital formed a very small hybrid peak,and O2s orbital PDOS did not obviously overlapped with U5f orbital.However,O2p orbital and U6d,U5f orbital formed an unambiguous hybrid peak,the electronic charge of U6d and U5f orbital mainly transferred to O2p orbital.(iv)H1s,O2s, O2p PDOSs shifted towards the lower energy band after H2O chemisorption,indicating that the red shift effect also occurred in H2O frequency.U6d and U5f orbital both lost electrons.O2s PDOS did not obviously overlapped with U6d and U5f orbital,that is to say,O atom in H2O and U atom on the UN(001)surface formed a covalent bond,and a small amount of U6d orbital electrons transferred to O2p orbital.In the future,we plan to investigate other molecular and atoms chemisorption on the surface of actinide compounds(especially the more promising nuclear fuel,such as Pu and U compounds),provide the corresponding anti-corrosion techniques,and improve the operational life and efficiency.

VI.ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China(No.51271198)and Self-Topics Fund of Xi’an Research Institute of High Technology(No.YX2012cxpy06).Ru-song Li would like to thank Wen Li from Xi’an Research Institute of Hi-Tech for useful discussions and studentship support.

[1]P.D.Wilson,The Nuclear Fuel Cycle,Oxford:Oxford University Press,(1996).

[2]D.Bocharov,D.Gryaznov,Y.F.Zhukovskii,and E.A. Kotomin,Surf.Sci.605,396(2011).

[3]M.B.Shuai,H.R.Hu,X.Wang,P.J.Zhao,and A.M. Tian,J.Mol.Struct.:THEOCHEM 536,269(2001)

[4]P.P.Dholabhai and A.K.Ray,J.Alloys Compd.444-445,356(2007)

[5]C.J.Burns,Science 309,1823(2005).

[6]E.N.Hodkin and M.G.Nicholas,J.Nucl.Mater.47, 23(1973).

[7]E.N.Hodkin,J.Nucl.Mater.67,171(1977).

[8]H.Shibata,T.Tsuru,M.Hirata,and Y.Kaji,J.Nucl. Mater.401,113(2010).

[9]P.F.Weck,E.Kim,N.Balakrishnan,F.Poineau,C. B.Yeamans,and K.R.Czerwinski,Chem.Phys.Lett. 443,82(2007).

[10]R.A.Evarestov,A.V.Bandura,M.Losev,E.A.Kotomin,Y.F.Zhukovskii,and D.Bocharov,J.Comput. Chem.29,2079(2008).

[11]L.Petit,A.Svane,Z.Szotek,W.M.Temmerman,and G.M.Stocks,Phys.Rev.B 80,045124(2009).

[12]S.Sunder and N.H.Miller,J.Alloys Compd.271,568 (1998).

[13]M.S.Brooks and D.Glotzel,Physica B 102,51(1980).

[14]M.S.Brooks,J.Phys.F 14,639(1984).

[15]D.Sedmidubsky,R.J.Konings,and P.Novak,J.Nucl. Mater.344,40(2005).

[16]E.A.Kotomin,R.W.Grimes,Y.Mastrikov,and N.J. Ashley,J.Phys.:Condens.Matter.19,106208(2007).

[17]A.F.Raymond and A.K.Ray,Phys.Rev.B 76,115101 (2007).

[18]A.G.Ritchie,J.Nucl.Mater.102,170(1981).

[19]K.Winer,C.A.Colmenares,and R.L.Smith,Surf. Sci.183,67(1987).

[20]R.A.Evarestov,A.V.Bandura,M.V.Losev,E.A.Kotomin,Yu.F.Zhukovskii,and D.Bocharov,J.Comput. Chem.29,2079(2008).

[21]Yu.F.Zhukovskii,D.Bocharov,E.A.Kotomin,R. A.Evarestov,and A.V.Bandura,Surf.Sci.603,50 (2009).

[22]Yu.F.Zhukovskii,D.Bocharov,and E.A.Kotomin,J. Nucl.Mater.393,504(2009).

[23]B.Delley,Comput.Mate.Sci.117,122(2000).

[24]B.Delley,J.Chem.Phys.113,7756(2000).

[25]B.Delley,Phys.Rev.B 66,155125(2002).

[26]A.F.Raymond and A.K.Ray,Phys.Rev.B 75,195112 (2007).

ceived on January 30,2013;Accepted on December 5,2013)

?Author to whom correspondence should be addressed.E-mail：rusong231@126.com