羟基磷灰石表面能的第一性原理计算

羟基磷灰石表面能的第一性原理计算1
王志明，黄远，万怡灶，何芳，王玉林
天津大学材料科学与工程学院，天津（300072）
E-mail：zmwangtju@
摘要：基于密度泛函理论（DFT）框架下的第一性原理，构建层晶（slab）模型，计算了
(100)、(100)、(010)、(001)、(101)晶面的表面能，由此分析了羟基磷灰石（HAP）的
slab模型的原子层厚度对表面能的影响。

采用布拉维法则和唐纳-哈克定律（BFDH）形貌预测方法预测了HAP各晶面的形貌重要性。

计算结果表明：当真空层厚度达12 Å时，slab模型的原子层厚度对表面能影响不大。

表面能计算结果与BFDH方法预测所得结果均表明：(100)面是最稳定晶面，具有最高的形貌重要性，在晶体生长过程中最易成为宏观晶体的表面。

计算结果对于研究HAP作为体内植入物与体液环境的相互作用过程具有重要意义。

关键词：密度泛函理论，表面能，羟基磷灰石
1. 前言
羟基磷灰石的化学组成和晶体结构与人体骨骼和牙齿中主要矿物质相类似，是人体骨骼组织的主要无机成分[1-2]。

同时，由HAP为主要成分构成的生物活性陶瓷具有良好的生物相容性和生物活性，能在骨修复、骨重建等方面起到重要作用，因而受到世界各国科学家的重视而广泛研究[3-5]。

作为重要的植入材料，HAP植入体内后将在体内环境下发生复杂的作用过程，如晶体的降解、晶体离子与体液中离子的交换、新骨组织的形成[6-7]等，HAP的表面特性与这一系列过程都有着直接的关系。

研究HAP的表面特性，将有助于对整个作用过程的理解，而表面能是表征材料表面特性的重要参量，因此，有必要对HAP各晶面的表面能进行研究。

表面能定义为形成单位面积表面时体系能量的增加，为材料表面的基本属性之一。

目前，很难通过实验方法[8-9]精确测定固态材料的表面能，即使通过某些实验方法测得一些简单结构晶体的表面能，误差一般也比较大，而且不能给出究竟是哪个面上的表面能。

对于HAP 这样结构复杂的晶体，目前还没有给出特征晶体表面能的直接测量结果。

因此，通过理论计算的方法来获得HAP的表面能有着重要的意义。

理论上，计算晶体表面能主要应用基于第一性原理的计算方法[10]和半经验分子动力学方法[11]，第一性原理计算方法能够精确到电子层次，计算量较大，但计算结果更为精确。

晶体的宏观形貌主要是由其内部结构决定的，在忽略外部因素影响的情况下，晶体的形态特征服从BFDH法则[12-13]，因此可以通过晶体的内部结构预测其宏观形貌，比较各晶面的形貌重要性。

综上所述，本文采用基于密度泛函理论（DFT）框架下的第一性原理方法，建立 slab 模型，计算了HAP的多个晶面的表面能，分析了slab模型的原子层厚度对表面能的影响。

采用Morphology程序的BFDH方法预测HAP晶体的宏观形貌，以比较HAP各晶面的形貌重要性。

2. 计算方法
本文计算工作采用基于密度泛函(DFT)方法的从头算量子力学程序CASTEP[14] 1本课题得到天津市科技发展计划重点攻关项目(No.043185111-2)的资助。

(Cambridge Serial Total Energy Package)软件包完成。

计算中选用基于广义梯度近似（GGA ）交换关联近似的超软赝势（ultrasoft ）[15]，平面波截断能量E off =300 ev ，迭代过程中的收敛精度为2×10-5 ev /原子，能量计算都在倒易空间中进行。

计算内容分以下几步：
(1) 首先，完成HAP 晶体结构的构建，见图1所表示。

对所构HAP 晶体结构进行晶格常数计算，将计算获得的晶格常数及各原子坐标与实验结果对比，以此验证计算参数的合理性；
(2) 以(100)晶面为例，建立slab 模型。

真空层厚度取12 Å，该厚度可消除slab 模型之间两端的相互作用。

建立不同原子厚度的多个slab 模型，对表面能进行收敛性测试，以确定合适的原子厚度。

(3) 建立其它各个晶面的表面模型，计算相应表面能。

图2所示为所构建的其中三个晶面的表面模型。

各表面模型的真空层厚度均为12 Å。

(4) 采用BFDH 方法预测HAP 宏观形貌，以比较各晶面的形貌重要性。

图1 HAP 的晶体结构模型
(Ca ：绿色；P ：紫色；O ：红色；H ：白色)
上述计算中，表面能采用如下公式[16-17]计算：
()/surf slab bulk E E nE A =− (1)
其中：E surf 是所要计算的表面能，E slab 表示模型能量值，E bulk 表示晶体内部原子数目相同时系统的能量值，A 表示表面积大小，n 为原胞数目。

(100)(100)(100)(010)
图2 HAP各表面层晶模型
(Ca：绿色；P：紫色；O：红色；H：白色)
3. 结果与讨论
实验结果表明：HAP晶体属于P63/m 空间群，六角晶系，化学式为Ca10(PO4)6(OH)2,晶格常数为a=b=9.424Å, c=6.879 Å[18]。

HAP晶格常数计算结果以及晶体中各原子坐标计算值分别见表1和表2。

由表可知，计算结果与实验值较吻合，因此，所选用的计算参数合理，可以用于对HAP系统的计算。

表1 HAP晶体的原子坐标(Å)
表2 HAP晶体的晶格参数(Å)
This work Experiment[18]
a b c a b c
9.433 9.433 6.874 9.424 9.424 6.879
为了保证计算结果的准确性，比较了不同原子厚度的slab模型对表面能的影响，考虑到计算量太大，仅以
(100)面为例进行了表面能的收敛性测试。

图3所示为随原子层厚度的增加而得到的
(100)面的不同表面能值，可见，原子厚度的增加对表面能的影响并不大，因此，为了减少计算量，其他晶面能量的计算均采用一层原子厚度。

This
work Experiment[18]
x y z x y z
O(1) 0.3287 0.48410.25000.32840.48480.2500
O(2) 0.5862 0.46530.25200.58730.46510.2500
O(3) 0.3428 0.25870.06970.34370.25790.0702
P 0.3685
0.36790.25000.39870.36850.2500
Ca(1) 0.3333 0.66670.00230.3333 0.6667 0.0015
Ca(2) 0.2571 0.98780.24920.24680.99340.2500
O(4) 0.0000 0.00000.25430.00000.00000.1950
H 0.0000
0.00000.06730.00000.00000.0608
图3 原子厚度(层数)对HAP
(100)晶面表面能的影响
对表面进行充分弛豫后，计算了各表面模型的总能量，根据式（1）得到了HAP 在(100)、(100)、(010)、(001)、(101)各晶面上的表面能，结果如表3所示：
(h k l) (100)(001)(101) Literature results [19] 0.3322
0.32870.3381- -
由表3可知，计算得到的各表面能的大小与文献结果接近，各晶面表面能大小顺序为： (010) > (100) > ( 001) > (101) >（100）。

其中（100）面表面能最低，是最稳定晶面。

图4 BFDH 方法预测得到的HAP 晶体的宏观形貌
Fig.4 the morphology of HAP crystal derived using BFDH prediction method
所得的HAP 宏观形貌,该图显示（100）、（101）、（001）三个面具有较高的形貌重要性。

计算所得结果按各晶面的形貌重要性顺序排列如下：（100）> (101) > ( 001) > (100)
(101)
(100)>(010)。

由于晶面形貌重要性与其表面能大小成反比，因此，上述表面能计算结果与采用BFDH方法所预测所得的各晶面形貌重要性顺序是一致的。

可见，HAP的（100）面表面能最低，具有最高的形貌重要性，在晶体生长过程中最有可能生长成为宏观晶体的表面。

因此在研究HAP植入物在体液环境下降解过程、以及HAP 与其它离子或小分子等的作用过程时，（100）面是需要考虑的最为重要的一个晶面。

很多实验结果也表明，一些带负电的颗粒物如小分子[20]，聚合体[21]等都以(100）面作为结合位置与HAP发生吸附作用。

4. 结论
本文采用第一性原理总能计算方法、建立slab模型，计算了HAP晶体几个重要晶面的表面能，并分析了所建slab模型的原子厚度与表面能的关系。

采用morphology程序的BFDH 方法对HAP的宏观形貌进行了预测。

结果表明：当slab模型的真空层厚度取12 Å时，原子层厚度对表面能的影响已不大；对HAP各晶面表面能计算结果和用BFDH方法对HAP 晶体宏观形貌的预测结果得到以下结论：各晶面的形貌重要性顺序为（100）> (101) >( 001)> (100)> (010)。

其中（100）表面能最低，具有最高的形貌重要性。

参考文献
[1] R.Z. LeGeros. Calcuim phosphate materials in restorative dentistry: a review.Adv Dent Res,1988,2:164
[2] W. Suchanek, M. Yoshimura. Processing and Properties of Hydroxyapatite-BasedBiomaterials for Use as Hard Tissue Replacement Implants[J]. Mater Res, 1998, 13: 94～117
[3] R.G.T. Geesink. Exprtimental and clinical experience with hydroxyapatitd-coated hip implants. Orthopedics, 1989,12:1239
[4] KG. Nilsson, S. Cajander, J. Karrholm. Early failure of hydroxyapatite coating in total knee arthroplasty[J]. Acta. Orthop Scand. 1994, 65: 212–214.
[5] B. Locardi, UE. Pazzaglia, C. Gabbi, et al. Thermal behavior of hydroxyapatite intended for medical applications. Biomaterials. 1993, 14: 437– 441.
[6] S.K. Lower, P.A. Maurice, S.J. Traina, et al. Simultaneous dissolution of hydroxylapatite and precipitation of hydroxypyromorphite:Direct evidence of homogeneous nucleation [J]. Acta, 1998, 63:1773-1780.
[7] de Groot K, Geesink R, Klein C, Serekian P. Plasma-sprayed coatings of hydroxylapatite [J]. Biomed Mater Res 1987,21:1375-87.
[8] W. Rtyson, W A .Miller. Surface free energies of solid metals: Estimation from liquid surface tension measurements[J].Surf. Sci, 1997,62:67-276.
[9] J. Wan, YL. Fan, DW. Gong eta1.Surface relaxation and stress of fcc metals: Cu, Ag, Au, Ni, Pd, Pt, Al and Pb [J] Mode1.Simu1. Mater. Sci. Eng, 1999, 7:189-206.
[10] M. Methfessel, D. Henning, M. ScheZer. Trends of the surface relaxations, surface energies, and work functions of the 4d transition metals [J].Phys. Rev. B,1992, 46:4816．
[11] J. Kollar, L. Vitos, H. L. Skriver. Surface energy and work function of the light actinides[J]. Phys. Rev. B, 1994, 49:11288-11292.
[12] J.D.H. Donnay, D. Harker. A new of crystal morphology extending the law of bravais, Amer. Mineral, 1937, 22: 446.
[13] P. Hartman, W.G. Perdok, On the relations between structure and morphology of crystals, Acta Cryst., 1955, 8: 525
[14] M. D. Segall, P. J. D. Lindan, M.J. Probert, et al. First-principles simulation: ideas i- ustrations and the CASTEP code [J]. J Phys Cond Matt, 2002, 14(11): 2717-2743.
[15] D. Vanderbilt. Soft self-consistent pseudopotentials in generalized eigenvalue formalism [J]. Phys Rev B,1990, 41: 7892~7895
[16] S. Casassa, C. Pisani. Atomic-hydrogen interaction with metallic lithium: An ab initio embedded - luster study[ J ]. Phys. Rev. B,1995, 51: 7805-7816.
[17] H. M. Polatoglou, MethfesselM, et. al. Vancancy-formation energies at the (111) surface and in bulk Al, Cu, Ag and Rh[J]. Phys. Rev. B, 1993, 48: 1877-1883.
[18] K. Sudarsanan, R.A. Young. Significant Precision in Crystal Structural Details: Holly Springs Hydroxyapatite Acta. Cryst[J]. 1996, B25:1534-1543.
[19] Weihua, Zhu; Ping, Wu. Surface energetics of hydroxyapatite: a DFT study. Chem. Phys. Lett. 2004, 396, 38
[20] J.C. Lieske, R. Leonard, H. Swift. et al. Adhesion of calcium oxalate monohydrate crystals to anionic sites on the surface of renal epithelial cells[J]. AJP. 1996, 270, F192.
[21]E. Bertoni, A. Bigi, G. Falini, et al.: Hydroxyapatite/polyacrylic acid nanocrystals. J Mater Chem 1999, 9:779-782.
First-principles calculations of the energies of
hydroxyapatite surfaces
Wang Zhiming, Huang Yuan, Wan Yizao, He Fang, Wang Yulin School of Materials Science and Engineering, Tianjin University, Tianjin (300072)
Abstract
The surface energies of HAP (100)、HAP (100)
、HAP
(010)
、HAP（001）and HAP（101）surfaces
were theoretically calculated using first-principle approach. The effects of slab thickness on surface energy have been tested. The Bravais-Friedel-Donnay-Harker (BFDH) morphological prediction method was used to predict the morphological importance of the faces of HAP. The results indicated that when the vacuum width is 12Å or more, the influence of slab thickness on surface energy is small. The results also show that the (100) surface is the most stable surface and has the most morphological importance. This conclusion is significant for studying the reactions occurred between HAP implantation and body fluid.
Keywords: DFT; Surface energy; Hydroxyapatite。