4D17周期境界条件(PBC)計算に対するエネルギー密度解析(EDA)の開発とその応用（早大理工） ○倉林 佑二、中井 浩巳【緒言】
近年固体物理で注目される現象のなかで、構造欠陥やドーピング、複雑な形をした固体表面などは、特異な物性を示すことが実験・理論の両面から数多く報告されている。このような現象のバリエーションは数知れず、無限の可能性を秘めているという点で大変興味深い。そのメカニズムを解明するためには、構造欠陥が生じた箇所やドーピングされた箇所とともに、その周辺に及ぼす影響を局所的に調べる必要がある。しかし、このような効果は、実験により見積もることは難しいので、理論的に固体結晶の局所的な効果を評価できる新たな解析手法の開発が望まれる。
固体モデルに対して理論的な取り扱いをする際には、クラスターモデルと周期モデルが用いられる。クラスターモデルでは、系を再現するためには大規模なサイズが必要となる。しかし、計算コストとの兼ね合いから結晶の一部を切り出した構造を用いることになる。その際に、クラスターの適切なサイズ及び形を決定することは困難であり、クラスターの端の効果により電子状態の記述が悪くなる。一方、周期モデルでは、固体結晶は単位構造の繰り返しであるという周期性を考慮し、単位構造に対して周期的に寄与を取り入れて計算する。よって、無限系を表現することができ、クラスターでの端による悪さも生じない。周期モデルは固体結晶に対して有効な計算方法であり、近年その適用例は増えてきている。
化学現象の解明には局所的な相互作用を見積もる必要があるが、そのための適当な解析手法がないのが現状である。そこで、このような問題を解決すべく、最近、当研究室ではクラスターモデルに対するエネルギー密度解析(EDA)[1]が開発された。EDA は、電子状態計算によって得られた系の全エネルギーを、定量性を失うことなくその構成原子に分割するという解析手法である。
そこで本研究では、固体結晶の記述がよい周期境界条件(PBC)[2]計算に対して、局所的な解析手法である EDA を適用した。本手法により、固体モデルにおいて局所的な相互作用をエネルギーとして見積もることが可能となり、固体モデルに対する現象を解明するための強力なツールとなり得る。
今回の発表では、複雑な表面を持つ MgO(n10) (n=1, 2, 3, 4)面に対して適用し、面の局所的な安定性について考察する。

【理論】
本研究での PBC 計算は、従来の分子軌道法の考え方を周期系に拡張したものである。したがって、Gauss 関数により基底関数展開を行うが、周期境界条件が課されるため、Bloch の定理に従う。
±∞

ϕ µ(k ) (r ) = ∑ e ik⋅g χ µ (r − g )g ; lattice vector, k ; wave number vector, χ ; Gaussian functiong

結晶軌道は Bloch Gaussian の線形結合をとる。
NBasis ±∞NBasisψ n(k ) (r ) = ∑ Cµ(kn )ϕ µ(k ) (r ) = ∑ ∑ Cµ(kn )eik⋅g χ µ (r − g ) C ; crystal orbital coefficientµµg

EDA はいくつかの分割スキームを用いている。電子部分に関しては Mulliken の電子密度解析[3]の手法に倣い、エネルギーを分割する。たとえば、運動エネルギー項の分割は、次の通りである。
±∞ NBasis NBasis

AEkin=∑ ∑gµ∈ APνµ( )Tµν( )∑νgg⎛ 1 ⎞Tµν(g ) = ∫ χ µ (r )⎜ − ∇ r2 ⎟ χν (r − g )dr⎝ 2 ⎠また、核電子引力エネルギーは、核と電子の 2 体間の相互作用であるから、相互作用を半分にして電子と核に分割する。
AE NucAtt=

1 NAtoms ±∞ NBasis NBasis (g ) ±∞ (g ,g′ ) 1 ±∞ NBasis NBasis (g ) ±∞ (g ,g′ )V A, µν + ∑ ∑ ∑ Pνµ ∑ VA, µν∑ ∑ ∑ ∑ Pνµ ∑2 A g µ∈B ν
2 g µ νg′g′,g′ )VA(g, µν= ∫ χ µ (r )− ZAχν (r − g )drr − R A − g′

DFT 計算における交換相関エネルギーはグリッドを用いた数値積分により求められる。
UnitCellAE XC=∫

A(r )dr ≅ ∑∑ω g p A,g (rg )ε XC (rg )ρ (r )ε XC
grid ±∞
ggw (r )p A,g (r ) = NAtoms ±A∞,g∑ ∑ wB,m (r )Bmω g はグリッドの範囲の重み、 p A,g (r ) は原子 A の partition function、 wA,g (r ) は原子 A の重み関数、
ε XC (r ) は交換相関エネルギー密度である。
【結果・考察】MgO 結晶のさまざまな表面について、2 次元 PBC 計算を行った。計算方法は BLYP、基底関数はMg には 6-21G、O には 6-31G を用いた。Fig. 1 に MgO(310)面の断面図を示す。(310)面は階段状になっており、表面の各サイトにはそれぞれ名前がついており、キンク(凸)とキンク(凹)の間には、2つのテラスサイトがある。Table 1 に、MgO(n10) (n=1, 2, 3, 4)、(100)表面での、各サイトにおける原子化エネルギーを示す。それぞれ 5 種類の面において、キンク（凸）からテラスを経て、キンク（凹）へと順番に数値を示す。
まず、キンク（凸）に関しては、(n10)面は n が大きくなると、表面が(100)面をとるテラスの露出が多くなる。よって、(n10)面の原子化エネルギーは(110)面と(100)面との間の値をとっている。テラスに関しては、表面が(100)面の構造となるので、(n10)の原子化エネルギーは(100)面に近い値をとっている。
キンク（凹）に関しても、(100)面のテラ

キンク（角状の箇所）凸キンク（角状の箇所）凹

スとステップに囲まれるため、(100)面に近い値をとる。(210)面に関しては、

ステップ（階段の段差）

MgO (310)面

この規則には当てはまらず、特殊な物

テラス（階段の踊り場）
Fig. 1. MgO(310) surface site

性を示す可能性があると考えられる。
以上のことから、PBC 計算に対して
EDA を適用したことにより、複雑な構造の表面における局所的なサイトの安定性を見積もることができた。

【参考文献】

Table 1. MgO atomization energy (kcal/mol) of each surface site.
surface
キンク（凸）
テラス

キンク（凹）

(110)
213
―――
232

[1] H. Nakai, Chem. Phys. Lett., 363, 73 (2002).
[2] R. Dovesi et al., Phys. Stat. Sol., 217, 63 (2000).
[3] R. S. Mulliken, J. Chem. Phys., 23, 1833 (1955).

(210)
228
180―
―
179(310)
220
238
235―
240(410)
220
237
245
234
238

(100)
―
246
―
――