Theoretical Analysis on Temperature- and Pressure-Dependences of NO-CO-O₂ Reaction on Rh(111) Surface

Abstract

三元触媒による自動車排ガス浄化反応において，NO_xの還元反応の活性や選択性は反応条件に敏感なことが知られている．この反応は，低温あるいは高温条件下，または高O₂分圧下で活性が低下する．本研究では，Rh(111)表面におけるNO-CO-O₂モデル反応系に対して，密度汎関数理論に基づいて反応速度を解析した．特に，吸着子の表面被覆率を考慮することで，固気平衡に対する温度および圧力の依存性をあらわに取り込んだ．NO還元反応の転化率が温度に対して最大値をとること，温度の上昇に伴ってNO還元反応のメカニズムがN + NO再結合からN + N再結合に変化することを見出した．吸着子の表面被覆率に対する詳細な解析から，N原子の被覆率がNO還元反応の活性および選択性に支配的な因子であることを見出した．

1 はじめに

三元触媒は，ガソリンエンジン車の排気ガス中に含まれるCO, NO_x, 燃え残った燃料(HC; hydrocarbon)を，同時に酸化還元することで無害なCO₂, H₂O, N₂に変換する．三元触媒の構成要素のひとつである金属Rh触媒はN-O結合解離能を有し，NO_x還元反応に高い活性を示す．このような背景から，金属Rh表面でのNO分子の振る舞いが，触媒化学のみならず表面化学の観点からも詳細に検討されてきた [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]．

三元触媒の構成要素であるPtおよびPdとの比較では，Pt(111)表面およびPd(111)表面ではNO解離が不活性である一方で，Rh(111)表面では極めて効率的に表面NO分子がN原子とO原子に解離することが見出された [1,2,3]．NO_x還元反応全体に対しては，実験的研究と理論的研究の双方でNO-CO反応，あるいは共存酸素を考慮したNO-CO-O₂反応がモデル反応系として頻繁に取り扱われてきた [4,5,6,7,8,9,10]．

NO-CO (-O₂)反応機構に対してこれまでに以下の知見が得られている．

(1) NO分子およびCO分子がRh表面上の活性サイトに吸着し，NO*, CO*となる．

(2) NO*がN*, O*に解離する．

(3) N*が他のN*あるいはNO*と再結合しN₂*またはN₂O*を生成する．

(4) N₂O*がN₂*とO*に解離する．

(5) O*がCO*と再結合しCO₂*を生成する．

(6) N₂*とCO₂*が表面から脱離し，外部に放出される．

ここで，アスタリスク(*)付きの化学式は吸着子を意味する．

三元触媒反応の効率は，反応ガスの温度と組成に依存する [11,12,13,14,15,16,17]．NO_x還元反応では，低温条件およびO₂分圧が高い条件において活性が大きく低下する．これは，ガソリンエンジン車の燃費低減の妨げとなっている．空燃比(燃焼ガス中の燃料質量に対する空気質量の割合)を高めたリーンバーンエンジンでは，排気ガスが酸素過剰となる．また，電動モーターを組み合わせたハイブリッド車では排気ガスが低温である場合が多い．加えて，NO_x還元反応には最適温度が存在し，低温条件のみならず極度の高温でも反応活性が低下する．NO_x還元反応における温度および圧力に対する依存性は，三元触媒反応の効率化，さらにガソリンエンジンの燃費低減に重要な課題である．また，N*が直接他のN*と再結合してN₂*を生成するか(N + N機構)，NO*と再結合してN₂O*を生成したのちに解離反応によってN₂*を生成するか(N + NO機構)が，温度に依存して変化することが知られており，このメカニズムも基礎的な観点から興味深い．

密度汎関数理論(DFT)は，表面反応の微視的解析に有効な手段である．最近では，DFT計算と反応速度論シミュレーション，化学工学シミュレーションを組み合わせることで，反応器中の反応ガスの濃度変化に対する微視的検討なども行われている [17]．

既往の研究では，DFT計算において清浄な表面を仮定することが多い．すなわち，(擬) 2次元周期的境界条件を課された金属スラブ上における単分子(あるいは反応に関与する最低限の反応分子)のみが考慮される．結果として，吸着子の被覆率はスラブモデルのサイズから自動的に決定される．しかしながら，実際には，反応ガスの温度や圧力の効果が吸着平衡を通じて表面被覆率に直接影響し，これを反映して反応エネルギーも変化する．特に，吸着(脱離)エネルギーは表面被覆率に敏感で，昇温脱離(TPD; temperature programmed desorption)スペクトルのシミュレーションでは，清浄表面を仮定すると定性的にも誤った結果を与える．この問題は，吸着エネルギーの被覆率依存性を考慮することで解決される [18]．

このような背景のもと，本研究ではDFTに基づいて表面被覆率をあらわに考慮してRh(111)面上におけるNO-CO-O₂反応を検討した．反応速度解析によって，温度および圧力に依存した反応転化率変化，また種々の表面吸着子の被覆率変化を明らかにした．

本論文の構成は，次の通りである．第2章では，理論的背景と計算方法について述べる．特に，被覆率を考慮した反応速度解析について詳細に述べる．第3章では，計算結果を示し，本反応における温度と圧力の効果を議論する．特に，反応機構の温度依存性とO₂分圧依存性について考察する．第4章では，本論文の結論を述べる．

2 理論的背景・計算方法

2.1 被覆率依存の反応速度解析

本研究では，既報 [17]を参考にNO-CO-O₂反応に不可欠な10の素反応を検討した．ここで，気相分子として反応物，生成物，および共存酸素を検討した．すなわち，NO, CO, O₂, N₂, CO₂が吸着平衡にあるとした．また，酸素過剰条件ではRhが酸化され部分的にRh₂O₃となることが，operandoのX線吸収分光法によって示されているが [16]，本研究ではO₂分圧によらず金属Rh上での反応を検討した．    NO + * → NO*    (1)       NO* + * → N* + O*    (2)       N* + NO* → N₂O* + *    (3)       2N* → N₂* + *    (4)       N₂O* + * → N₂* + O*    (5)       N₂* → N₂ + *    (6)       CO + * → CO*    (7)       CO* + O* → CO₂* + *    (8)       CO₂* → CO₂ + *    (9)       O₂ + 2* → 2O*    (10)    ここで，アスタリスク(*)は表面上の空きサイトを意味する．N*種の吸着エネルギーが大きいことを反映して，式(3), (4)の反応を律速段階と仮定した．

反応速度解析に際して，定常状態近似を適用した．すなわち，律速段階である式(3), (4)以外の素反応は，すべて平衡に至っていると仮定した．各素反応における平衡定数Kは，各吸着子の被覆率θ, 分圧p, 反応エネルギーEによって，次のように表される．ただし，θ_*は表面上の空きサイトの割合を意味する．下付きの数字は，上記の素反応の式番号に対応する．   

$K_{\text{1}}=\frac{{\theta }_{\text{NO}}}{p_{\text{NO}}{\theta }_{*}}=\exp \left(\frac{-E_1\left(1-{\theta }_{*}\right)}{RT}\right)$(11)

  

$K_{\text{2}}=\frac{{\theta }_{\text{N}}{\theta }_{\text{O}}}{{\theta }_{\text{NO}}{\theta }_{*}}=\exp \left(\frac{-E_2\left(1-{\theta }_{*}\right)}{RT}\right)$(12)

  

$K_{\text{5}}=\frac{{\theta }_{{\text{N}}_2}{\theta }_{\text{O}}}{{\theta }_{{\text{N}}_2\text{O}}{\theta }_{*}}=\exp \left(\frac{-E_5\left(1-{\theta }_{*}\right)}{RT}\right)$(13)

  

$K_{\text{6}}=\frac{p_{{\text{N}}_2}{\theta }_{\text{*}}}{{\theta }_{{\text{N}}_2}}=\exp \left(\frac{-E_6\left(1-{\theta }_{*}\right)}{RT}\right)$(14)

  

$K_{\text{7}}=\frac{{\theta }_{\text{CO}}}{p_{\text{CO}}{\theta }_{\text{*}}}=\exp \left(\frac{-E_7\left(1-{\theta }_{*}\right)}{RT}\right)$(15)

  

$K_{\text{8}}=\frac{{\theta }_{{\text{CO}}_2}{\theta }_{\text{*}}}{{\theta }_{\text{CO}}{\theta }_{\text{O}}}=\exp \left(\frac{-E_8\left(1-{\theta }_{*}\right)}{RT}\right)$(16)

  

$K_{\text{9}}=\frac{p_{{\text{CO}}_2}{\theta }_{\text{*}}}{{\theta }_{{\text{CO}}_2}}=\exp \left(\frac{-E_9\left(1-{\theta }_{*}\right)}{RT}\right)$(17)

  

$K_{\text{10}}=\frac{{\theta }_{\text{O}}{}^2}{p_{{\text{O}}_2}{\theta }_{\text{*}}{}^2}=\exp \left(\frac{-E_{10}\left(1-{\theta }_{*}\right)}{RT}\right)$(18)

式(11)-(18)より，以下の解析解が得られる．   

$\begin{array}{c}{\theta }_{\text{*}}=1/(1+K_1{p}_{\text{NO}}+\frac{p_{\text{NO}}{K}_1{K}_2}{\sqrt{p_{{\text{O}}_2}{K}_{10}}}+\frac{p_{{\text{N}}_2}\sqrt{p_{{\text{O}}_2}{K}_{10}}}{K_5{K}_6}\\ +\frac{p_{{\text{N}}_2}}{K_6}+p_{\text{CO}}{K}_7+\sqrt{p_{{\text{O}}_2}{K}_{10}})\end{array}$(19)

  

${\theta }_{\text{NO}}=K_1{p}_{\text{NO}}{\theta }_{\text{*}}$(20)

  

${\theta }_{\text{N}}=\frac{p_{\text{NO}}{K}_1{K}_2}{\sqrt{p_{{\text{O}}_2}{K}_{10}}}{\theta }_{*}$(21)

  

${\theta }_{{\text{N}}_{\text{2}}\text{O}}=\frac{p_{{\text{N}}_{\text{2}}}\sqrt{p_{{\text{O}}_{\text{2}}}{K}_{10}}}{K_5{K}_6}{\theta }_{*}$(22)

  

${\theta }_{{\text{N}}_2}=\frac{p_{{\text{N}}_{\text{2}}}}{K_6}{\theta }_{*}$(23)

  

${\theta }_{\text{CO}}=p_{\text{CO}}{K}_7{\theta }_{*}$(24)

  

${\theta }_{{\text{CO}}_2}=\frac{p_{{\text{CO}}_{\text{2}}}}{K_9}{\theta }_{*}$(25)

  

${\theta }_{\text{O}}=\sqrt{p_{{\text{O}}_2}{K}_{10}}{\theta }_{*}$(26)

また，律速段階における表面サイトあたりの反応速度(TOF; turn-over frequency) r₃, r₄は以下の式から計算できる [19]．    r₃ = k₃θ_Nθ_NO    (27)       r₄ = k₄θ_N²    (28)    式(19)より空きサイトの割合θ_*は平衡定数K₁, K₂, K₅, K₆, K₇, K₁₀に依存する．一方，それぞれの平衡定数は，式(11)-(15), (18)より，空きサイトの割合θ_*に依存する反応エネルギーE(1−θ_*)から決定される．結果として，式(19)は両辺にθ_*を含む非線形方程式となる．このため，繰り返し手続きによって自己無撞着な解を得る必要がある．本研究では，Python言語によって本手続きのプログラムを実装した．

2.2 計算条件

式(1)-(10)に示す素反応の反応エネルギーを，吸着子の被覆率をあらわに考慮してDFT計算から算出した．表面被覆率は，各素反応に関与しないNO分子をスラブ上に配置することで考慮した．

DFT計算は，Vienna ab initio Simulation Package (VASP) [20,21,22,23]を用いて行った．交換相関汎関数にはrevPBE [24]，基底関数にはカットオフエネルギー400 eVの平面波をそれぞれ用いた．内殻電子はprojector augmented wave (PAW)法によって記述した．Rh(111)面に対して，(4×4)ユニットセルおよび(4×2)ユニットセル4層からなるスラブモデルを構築し，k点サンプリングはそれぞれ5×5×1, 4×8×1とした．真空層は20 Åとした．構造最適化の際には，スラブモデルの下部を固定した．気相分子に対しては，1辺20 Åのスーパーセルを用いた．ステップ間のエネルギー差が10⁻⁵ eV以下であることを条件として，構造最適化の収束を判定した．

3 結果と考察

3.1 NOの吸着構造および反応エネルギーの被覆率依存性

Rh(111)面に対するNOの吸着構造を，NOの被覆率(θ_NO)が0.06 (1/16)から0.75 (12/16)までの12通りに対して網羅的に検討した．すなわち，(2×2)-3NO構造を満足する配置のうち，被覆率0.5 ML以下でatop位置に対する吸着を除外したすべての配置，合計38種類を検討した．それぞれの被覆率において最安定な吸着構造をFigure 1に示す．また，対応する吸着エネルギーの被覆率依存性をFigure 2に示す．

Figure 1.

 Stable conformations of NO adsorption on Rh(111) surface at each surface coverage of 0.06 (a) − 0.75 (l).

Figure 2.

 Relationship between the surface coverage of NO molecules on Rh(111) surface and NO adsorption energies (eV).

NOの被覆率が最も小さい条件(θ_NO = 0.06)では，吸着エネルギーは–2.12 eVで，発熱的に吸着した．被覆率の増加に伴って吸着による安定化が単調に減少し，最も被覆率の高い条件(θ_NO = 0.75)では，吸着エネルギーは–1.65 eVであった．これは，吸着子間の斥力的な相互作用による吸着の不安定化を意味する．吸着エネルギーと被覆率との間には，次式の良好な線形の相関が見られた．    E₁(1 − θ_*) = 0.695(1 − θ_*) − 2.192    (29)   

NO吸着反応(式(1))と同様に，式(2)-(10)の反応に対しても線形の相関を得た．Table 1に，それぞれの傾きおよび切片をまとめる．いずれの場合においても，斥力的な吸着子間相互作用を示唆する結果を得た．すなわち，空きサイトを消費する反応(式(1), (2), (5), (7), (10))に対しては傾きが正であった．また逆に，空きサイトが生成される反応(式(3), (4), (6), (8), (9))は傾きが負となった．

Table 1.  Gradients and intercepts for the relationships between surface coverage and reaction energies (eV).

	Gradient	Intercept
(1)	0.695	−2.192
(2)	1.914	−0.783
(3)	−0.980	1.556
(4)	−3.327	−0.344
(5)	0.090	−2.524
(6)	−0.649	0.400
(7)	0.688	−1.839
(8)	−1.223	1.125
(9)	−0.151	−0.471
(10)	10.383	−5.243

3.2 NO-CO-O₂反応の反応速度解析

Rh(111)表面上のNO-CO-O₂反応において，反応速度と表面吸着種の被覆率が温度および酸素分圧にどのように依存するかを検討した．ここで，三元触媒反応では反応ガス(自動車排ガス)が定常的に供給されていることを反映して，気相にある化学種の分圧は一定であると仮定した．化学量論を考慮して，COおよびNOの分圧はそれぞれ1.0 × 10⁻⁴, 5.0 × 10⁻⁵ atm, CO₂の分圧は0.17 atmとした．O₂分圧は，10⁻², 10⁻¹, 1 atmを検討した．また，式(10)の素反応を除いたモデルに対しても検討した．これは気相にO₂分子が存在しない条件に相当する．

Figure 3に各酸素分圧に対するTOFの温度依存性を示す．ここで，青色はr₃で表されるN + NO反応の速度を，オレンジ色はr₄で表されるN + N反応の速度を，緑色はr₃ + r₄で表される全反応速度をあらわす．また，Table 2にそれぞれの反応速度が最大となるピーク位置の温度とTOFを，Figure 4に吸着子の被覆率について示した．

Figure 3.

 Temperature dependences of TOFs for (a) no gaseous oxygen condition, and the conditions with O₂ partial pressures of (b) 10⁻², (c) 10⁻¹, and (d) 1 atm, respectively. TOFs for N + N and N + NO reactions are shown in orange and blue lines, respectively, while that for total is shown in the green line.

Table 2.  Peak temperatures (K) of TOFs (site⁻¹ s⁻¹) for N +N, N + NO, and total reactions.

pO2 (atm)	N + N			N + NO			total
	Temp.	TOF		Temp.	TOF		Temp.	TOF
No O2 gas	780	78.4		715	74.1		750	121.3
10−2	745	120.4		705	91.4		735	186.2
10−1	710	4.3		710	43.8		710	48.1
1	705	0.0		700	5.1		700	5.1

Figure 4.

 Temperature dependences of coverages of surface species for (a) no gaseous oxygen condition, and the conditions with O₂ partial pressures of (b) 10⁻², (c) 10⁻¹, and (d) 1 atm, respectively. NO, N, and O adsorbates and vacancy site are shown in orange, gray, green, and blue lines, respectively.

気相のO₂分子を考慮しない条件において，反応全体のTOFは750 Kで最大値121.3 site⁻¹ s⁻¹となった．TOFが温度に対して最大となることは，低温条件および極度に高温の条件のいずれにおいてもNO還元反応が不活性となるという実験事実と対応する．N + NO反応とN + N反応からの寄与に注目すると，それぞれ715, 780 Kで最大値74.1, 78.4 site⁻¹ s⁻¹となった．また，それぞれの反応速度の大小関係は740 Kで入れ替わった．すなわち，740 K以下ではN + NO機構が優勢であるのに対し，それ以上の温度ではN + N機構が優勢であることがわかった．

TOFの振る舞いは，表面吸着種の被覆率によって理解できる．NOの被覆率は温度上昇に伴って急速に減少した(Figure 4 (a))．これは，NO*の脱離(式(1)の逆反応)やNO解離反応(式(2))，N₂O生成反応(式(3))によってNO*が消費されたことを意味する．一方，N*はNO解離反応によって生成され，およそ670 Kから増加し760 Kで最大となる．N + NO機構の反応速度はθ_Nとθ_NOの積に比例するため(式(27))，最大値を持つ．さらに高温(> ca. 800 K)では，θ_Nも急速に減少する．結果として，N + NO機構とN + N機構のどちらも反応速度が低下する．

類似の傾向が，O₂分圧が10⁻² atmの条件でも見られた．反応全体のTOFは735 Kで最大値186.2 site⁻¹ s⁻¹となった．また，720 KでN + NO機構とN + N機構が入れ替わった．一方で，高温(> ca. 900 K)におけるO*種の被覆率において顕著な差が見られた．すなわち，気相のO₂を考慮しない条件(Figure 4 (a))では単調に増加する一方で，O₂分圧が10⁻² atmの条件(Figure 4 (b))では緩やかに減少した．また，空きサイトの割合は，O*の被覆率を補償する振る舞いを示した．これは，前者のモデルにおいてO₂の脱離反応を考慮しなかったことに起因するアーティファクトである．

実験的には，酸素過剰条件においてNO還元活性が急激に失われることが知られている．Rh(111)表面 [14]およびRh/Al₂O₃系 [15]において，NOに対するO₂の濃度がおよそNO : O₂ = 1 : 2以上では，N₂生成が観測されなくなる．これは，RhがRh₂O₃に酸化され [16]，NO解離反応およびN₂生成反応が困難になるためである，と説明される．本研究で検討したO₂分圧(10⁻² − 1 atm)は，一般的な実験の条件よりも極めて高いにも関わらず，一定程度NO還元反応が進行する結果が得られた．これはO₂分圧によらず金属Rh上での反応を検討したことに起因するアーティファクトであると考えられる．

反応ガス組成依存性に関しては，O₂分圧の増加に従ってθ_Nが顕著に減少した．すなわち，式(10)の反応によって供給されたO*により，式(2)で表されるNO*の解離平衡が反応物(NO)側に傾いたために，θ_Nが減少した．これを反映して，全反応速度もO₂分圧が10⁻² atmにおける186.2 site⁻¹ s⁻¹から1 atmの5.1 site⁻¹ s⁻¹に単調に減少した．結果として，仮にRh表面が金属状態に保たれていた場合でも，O₂分圧の増加により，N*の再結合(式(3), (4))が抑制されるためにNO還元反応が不活性となることが見出された．

4 結論

Rh表面上でのNO-CO-O₂反応は，三元触媒中でのNO_x還元反応のモデル反応として，実験的にも理論的にも重要である．本研究では，DFT計算および反応速度解析に基づいて本反応の温度および圧力に対する依存性を理論的に検討した．特に，反応エネルギーの表面被覆率依存性をあらわに取り込むことで，気相分子の吸着平衡を適切に考慮した．

反応速度解析から，反応温度に対してTOFが最大となる温度が存在すること，全反応速度に対するN + N機構およびN + NO機構の寄与がおよそ720 − 740 Kで切り替わることを確かめた．吸着種の表面被覆率に対する解析から， N*およびNO*の被覆率がNO-CO-O₂反応の温度・圧力依存性に重要な因子であることを見出した．

謝辞

本研究は，文部科学省・元素戦略研究拠点「実験と理論計算科学のインタープレイによる触媒・電池の元素戦略研究拠点(ESICB)」，早稲田大学・特定課題研究助成費の支援を受けて実施された．また，本研究の計算の一部は，自然科学研究機構(NINS)・計算科学研究センター(RCCS)および東京大学物性研究所(ISSP)・スーパーコンピュータシステムの計算機を利用して行った．

参考文献

• [1]   R. J. Gorte, L. D. Schmidt, J. L. Gland, Surf. Sci., 109, 367 (1981). doi:10.1016/0039-6028(81)90494-5
• [2]   T. W. Root, L. D. Schmidt, G. B. Fisher, Surf. Sci., 134, 30 (1983). doi:10.1016/0039-6028(83)90310-2
• [3]   D. T. Wickham, B. A. Banse, B. E. Koel, Surf. Sci., 243, 83 (1991). doi:10.1016/0039-6028(91)90347-U
• [4]   C. T. Campbell, J. M. White, Appl. Surf. Sci., 1, 347 (1978). doi:10.1016/0378-5963(78)90037-5
• [5]   L.H. Dubois, P. K. Hansma, G. A. Somorjai, J. Catal., 65, 318 (1980). doi:10.1016/0021-9517(80)90309-7
• [6]   S. B. Schwartz, G. B. Fisher, L. D. Schmidt, J. Phys. Chem., 92, 389 (1988). doi:10.1021/j100313a030
• [7]   C. S. Gopinath, F. Zaera, J. Catal., 186, 387 (1999). doi:10.1006/jcat.1999.2561
• [8]   C. S. Gopinath, F. Zaera, J. Phys. Chem. B, 104, 3194 (2000). doi:10.1021/jp993103g
• [9]   K. C. Taylor, Catal. Rev. - Sci. Eng., 35, 457 (1993). doi:10.1080/01614949308013915
• [10]   M. Shelef, G. W. Graham, Catal. Rev. - Sci. Eng., 36, 433 (1994). doi:10.1080/01614949408009468
• [11]   M. Shelef, R. W. McCabe, Catal. Today, 62, 35 (2000). doi:10.1016/S0920-5861(00)00407-7
• [12]   J. Kašpar, P. Fornasiero, N. Hickey, Catal. Today, 77, 419 (2003). doi:10.1016/S0920-5861(02)00384-X
• [13]   M. Okada, N. Koike, H. Suzuki, Chemosphere, Glob. Chang. Sci., 2, 413 (2000). doi:10.1016/S1465-9972(00)00042-8
• [14]   C. S. Gopinath, F. Zaera, J. Catal., 200, 270 (2001). doi:10.1006/jcat.2001.3215
• [15]   T. Yamamoto, T. Tanaka, R. Kuma, S. Suzuki, F. Amano, Y. Shimooka, Y. Kohno, T. Funabiki, S. Yoshida, Phys. Chem. Chem. Phys., 4, 2449 (2002). doi:10.1039/b201120b
• [16]   H. Asakura, S. Hosokawa, T. Ina, K. Kato, K. Nitta, K. Uera, T. Uruga, H. Miura, T. Shishido, J. Ohyama, A. Satsuma, K. Sato, A. Yamamoto, S. Hinokuma, H. Yoshida, M. Machida, S. Yamazoe, T. Tsukuda, K. Teramura, T. Tanaka, J. Am. Chem. Soc., 140, 176 (2018). doi:10.1021/jacs.7b0711429224338
• [17]   A. Ishikawa, Y. Tateyama, J. Phys. Chem. C, 122, 17378 (2018). doi:10.1021/acs.jpcc.8b05906
• [18]   C.-C. Wang, J.-Y. Wu, J.-C. Jiang, J. Phys. Chem. C, 117, 6136 (2013). doi:10.1021/jp309394p
• [19]   Fundamental Concepts in Heterogeneous Catalysis, J.K. Nørskov, F. Studt, F. Abild-Pedersen, T. Bligaard (Eds.), Wiley, New Jersey (2014).
• [20]   G. Kresse, J. Hafner, Phys. Rev. B, 47, 558 (1993). doi:10.1103/PhysRevB.47.558
• [21]   G. Kresse, J. Hafner, Phys. Rev. B, 49, 14251 (1994). doi:10.1103/PhysRevB.49.14251
• [22]   G. Kresse, J. Furthmüller, Comput. Mater. Sci., 6, 15 (1996). doi:10.1016/0927-0256(96)00008-0
• [23]   G. Kresse, J. Furthmüller, Phys. Rev. B, 54, 11169 (1996). doi:10.1103/PhysRevB.54.11169
• [24]   Y. Zhang, W. Yang, Phys. Rev. Lett., 80, 890 (1998). doi:10.1103/PhysRevLett.80.890