第19回CMSI神戸ハンズオン
OpenMXチュートリアル
２０１４年10月10日（金）CMSI神戸拠点講師
東京大学物性研究所
東京大学物性研究所
東京大学物性研究所

特任教授
尾崎泰助
拠点研究員 Truong Vinh Truong Duy
拠点研究員 坂下達哉

プログラム
13:00-13:35

13:35-13:50
13:50-14:40
14:40-14:50
14:50-15:50
15:50-17:00
17:00-17:30

OpenMXの概要
実習内容(1)の説明
インストール方法
各自実習(1)
休憩／コーヒーブレーク
計算の実際と実習内容(2)の説明
各自実習(2)
質疑・応答

チュートリアルコースの目的
密度汎関数理論に基礎付けられた第一原理電子状態計算手法は様々な物質の成り立ちを量子力学に基づき定量性を持って理解するために、非常に有効な手段です。
本チュートリアルコースでは密度汎関数理論、数値局在基底法、擬ポテンシャル法に基づく汎用第一原理電子状態計算ソフトウエア
OpenMXの使用方法を解説し、また簡単な実習計算を体験して頂きます。

またすでに計算を始めておられる方で、疑問に感じる点をお持ちの方のご質問も歓迎致します。
チュートリアル後に、皆様の御研究に活用して頂けますと幸いに存じます。

OpenMXの概要•

密度汎関数理論に基づく第一原理電子状態計算ソフトウエア•

ノルム保存型擬ポテンシャルを使用•

最適化数値局在基底により、Kohn-Sham軌道を展開

基本機能

拡張機能

•

LDA, GGA, LDA+UによるSCF計算•

オーダーN, 低次スケーリング法•

全エネルギー、力の計算•

ノンコリニア磁性•

バンド分散、状態密度•

スピン軌道相互作用•

構造最適化(BFGS,RF,EF法等)•

非平衡グリーン関数法による電気伝導•

有効電荷解析(マリケン,ボロノイ,ESP)•

Berry位相による電気分極•

NEV,NVT分子動力学法•

最局在化Wannier関数•

電荷ドーピング•

有効遮蔽体(ESM)法による電場印加•

フェルミ面解析•

NEB法による反応経路探索•

電子密度等のCubeデータ•

バンドunfolding法•

擬ポテンシャル&基底のデータベース 等•

TH法によるSTM解析 等

OpenMXの歴史
2000年 開発開始
2003年 公開 (GNU-GPL)
2003年 共同研究：
産総研
物材機構
ソウル大
金沢大
北陸先端大
日産自動車
富士通 等
2013年

１7回のリリース、Ver. 3.7に至る

http://www.openmx-square.org

OpenMXの開発者• T. Ozaki (U.Tokyo)• T. V. Truong Duy (U.Tokyo)• H. Kino (NIMS)• C.-C. Lee (JAIST)• J. Yu (SNU)• Y. Okuno (Fuji FILM)• M. J. Han (KAIST)• Yang Xiao (NUAA)• M. Ohfuti (Fujitsu)• Y. Gohda (TIT)• F. Ishii (Kanazawa Univ.)• T. Ohwaki (Nissan)• H. Weng (CAS)• M. Toyoda (Osaka Univ.)• H. Kim (SNU)• P. Pou (UAM)対象としている物質
実験との共同研究によりZrB2上のシリセン構造を初めて実証
A. Fleurence et al., Phys. Rev. Lett. 108, 245501 (2012).

実験との共同研究によりIr酸化物のJeff=1/2Mott状態を初めて実証

これまでの実績

B.J. Kim et al., Phys. Rev. Lett. 101, 076402 (2008).

トポロジカル絶縁体の理論的予測
C.-H. Kim et al., Phys. Rev. Lett. 108, 106401 (2012).
H. Weng et al., Phy. Rev. X 4, 011002 (2014).

リチウムイオン電池の電極界面の第一原理分子動力学計算
T. Ohwaki et al., J. Chem. Phys. 136, 134101 (2012).
T. Ohwaki et al., J. Chem. Phys. 140, 244105 (2014).

磁石材料の磁気異方性の第一原理計算
Z. Torbatian et al., Appl. Phys. Lett. 104, 242403 (2014).
I. Kitagawa et al., Phys. Rev. B 81, 214408 (2010).

表面酸化Si基板上のグラフェンナノリボンの大規模電気伝導計算
H. Jippo et al., Appl. Phys. Express 7, 025101 (2014).
M. Ohfuchi et al., Appl. Phys. Express 4, 095101 (2011).

シリセン、グラフェン
カーボンナノチューブ
遷移金属酸化物
金属酸化物
分子磁性体
希土類磁石材料
金属間化合物
リチウムイオン電池材料、鉄鋼材料
金属(酸化物)アモルファス 等

bcc-鉄中の金属炭化物析出物の界面構造予測
H. Sawada et al., Modelling Simul. Mater. Sci. Eng. 21, 045012 (2013).

アモルファス金属酸化物を特徴付ける普遍的な中距離秩序構造を提案
K. Nishio et al., Phys. Rev. Lett. 340, 155502 (2013).

関連論文
およそ240報

ZrB2上のシリセン構造を同定
日本物理学会2013年 第5号

表紙より

シリコン原子からなる蜂の巣格子構造(シリセン)は電子デバイス応用に有用な電子構造を持つことが理論的に知られており、実験による実現が期待されていた。
JAISTマテリアルサイエンス・高村(由)グルー
プはZrB2上にシリコンのネットワーク構造が形成することを発見。密度汎関数理論による大規模第一原理電子状態計算と実験事実を詳細に比較することで、最終的にZrB2上のシリセン構造の同定に初めて成功。

Si 2pの化学シフト
理論
実験

第一原理計算によって得られたZrB2上のシリセン構造

A. Fleurence et al., Phys. Rev. Lett. 108,245501 (2012).

グラフェンナノリボンデバイスの電気伝導計算
表面酸化Si基板上でのグラフェンナノリボンの電気伝導特性を非平衡グリーン関数による計算。 ゲート電圧による電流のon/off比は基板上で105程度で、実験値と整合。


電荷移動の解析か
らナノリボンから
基板表面に電荷移
動 → ナノリボンにホールドープ ショットキーバリ
アを形成し、チャ
ネル中央の価電子帯のトップがEFにピン止め
 ホール状態は酸素準位と混成し、ギャップ中に僅かに染み出す

 染み出したホール
準位を介してトン
ネル電流の増加

M. Ohfuchi et al., Appl. Phys. Express 4,095101 (2011).

BCC-鉄と金属炭化物の界面構造
鉄鋼材料の硬度・靱性の制御に有効な手段は金属炭化物を析出させることである。第一原理計算とモデルポテンシャル法を用いて析出物の構造予測を行った。

H. Sawada et al., Modelling
Simul. Mater. Sci. Eng. 21,
045012 (2013).

HRTEM image
実験値
2.2nm

析出物:

TiC, VC, NbC

整合析出
母体(Fe)
準整合析出

実習内容 (1)
 ログイン
phiにログインする (後で説明)。
 コンパイル体験
openmx3.7/sourceに移動してmakefileのCC, LIB, FCの設定を確認後、以下を実行。
% make clean% make install テスト計算
openmx/workに移動して、以下のコマンドでメタン分子の一点計算を実行。
% mpirun –np 5 openmx Methane.dat –nt 2 | tee met.std
 runtestの実行
正常にコンパイルが出来ていることを確認するために、以下の計算を実行。
% mpirun –np 8 openmx –runtest –nt 2

実際の計算イメージ
1. VESTA等を使用して構造モデルを作成
2. 入力ファイルをeditorを作成
3. 計算機にPBSシステム等を使用してジョブを投げる
4. XCrySDen等を使用して結果を可視化
Editorを用いた入力ファイル作成

XCrySDenを用いて差電子密度を解析

コンパイルの実行
OpenMXの実行ファイルを生成するためにopenmx3.7/sourceに移動してmakefileのLIB, CC, FCを計算機環境に応じて書き直す必要がある。適切な設定が出来たならばopenmx3.7/source中で$ make clean$ make installを実行する。コンパイルを完了するために7分程度を要する。
(今回の実習ではLIB, CC, FCは適切に設定されているため、変更の必要はありません)
正常にコンパイルが終了した場合には、openmx3.7/work中に実行ファイル openmx が生成される。
LIB, CC, FCの設定方法に関してはDuy氏による講義を参照のこと。
またマニュアルの9pも参照のこと。

テスト計算の実行
openmx/workに移動して、以下のコマンドでメタン分子の一点計算を実行。
$ mpirun –np 5 openmx Methane.dat –nt 2 | tee met.std-npによってMPIプロセス数を、-ntによってOpenMPスレッド数を指定する。
5秒程度で計算が終了し、11個のファイルと 一つのディレクトリが生成される。
マニュアルの15pを参照のこと。

Methane.datの中身(1) Keywordの後に値を与える。
(2) Keywordの順番は任意で良い。
(3) コメント文は #を文頭に入れる。

いくつかのkeywords:
システム名：この名前で様々なファイルを保存系に含まれる原子数
原子種の定義原子の座標

マニュアルの25pを参照のこと。

met.outの中身
SCFの履歴全エネルギー

KS固有値
マリケン電荷マニュアルの15pを参照のこと。

出力ファイル

詳細な説明がマニュアル42ページに記載されていますので、ご参照下さい。

可視化に関して
 ・*.tden.cube, *.sden.cube, *.dden.cube等のcubeファイルは様々な ソフトウエアで可視化できる。
XCrySDen
VESTA

等が使いやすい。

 また*.mdファイルはxyz形式で保存されており、これもまた分子動力学シミュレーションの可視化ソフトで構造最適化過程や分子動力学を可視化できる。
phiにはXCrSDenがinstallされており、$ xcrysden
として使用できる。
 その他、DOSやバンド図等は標準ではgnuplotで可視化する 。

runtestの実行
コンパイルが終了した後に、正常にコンパイルが出来ていることを確認するために、以下の計算を実行する。
$ mpirun –np 8 openmx –runtest –nt 2
計算が終了するとruntest.resultが生成する。

予め計算してある全エネルギーと力と比較し、diffを表示する。
マニュアルの22pを参照のこと。

OpenMXにおける局在基底法の理論と実装• 密度汎関数理論• KS方程式の数学的構造• LCPAO法の特徴• 全エネルギーの表式• 各項の積分方法• 力の計算• FFTグリッド• プリミティブ基底• 変分最適化基底• 最適化基底のデータベース• Δ因子による精度の評価

第一原理計算手法の分類
波動関数理論

計算の複雑さe.g., 配置間相互作用法 (CI)O(eN)
密度汎関数理論特徴
高精度
基底・励起状態
高コスト

O(N3)

中精度
基底状態
低コスト

O(N3～)

高精度
基底状態
高コスト
高並列性

O(N3～)

中精度
基底・励起状態
高コスト

量子モンテカルロ法

多体グリーン関数法

密度汎関数理論
基底状態のエネルギーは密度の汎関数で厳密に記述できる。 (Hohenberg and Kohn, 1964)

基底状態の多体問題は有効ポテンシャルに対する一体問題に帰着される。(Kohn-Sham, 1965)
W.Kohn (1923-)

一般化勾配近似(GGA)の成功事例と失敗事例
成功事例:
1. 精度:
凝集エネルギー:
結合長:
Bulk modulus:
エネルギーバリア:

平均誤差

0.3 eV
1%程度の過大評価
5%程度の過小評価
30％程度の過小評価

2. 高精度な水素結合
3. 磁気秩序の記述 (e.g., bcc Fe)

失敗事例:
1. バンドギャップ:
2. vdW相互作用:
3. 強相関系:

40%程度の過小評価
No binding in many cases
軌道分極が記述できない(d-,f-states)

KS方程式の数学的構造
3次元の連立非線形偏微分方程式を自己無撞着に解く必要がある。
OpenMX: LCPAO法OpenMX: 平面波-FFT法

入力電子密度 = 出力電子密度 → 自己無撞着条件

OpenMXの計算手法に関する論文等
2000 block BOP法の開発
PRB 59, 16061, PRB 61, 7972
2001 ab initio recursion 法の開発 PRB 64 195126
2001 recursion法に基づくO(N)逆行列計算法の開発 PRB 64, 195110
2003 変分基底法の開発 PRB 67, 155108, PRB 69, 195113, JCP 121, 10879
2005 LCPAO法の高精度実装法の開発 PRB 72, 045121
2005 O(N) LDA+U法の開発 PRB 73, 045110
2006 O(N) Krylov 部分空間法の開発 PRB 74, 245101
2007 高速・高精度Green関数積分法の開発 PRB 75, 035123
2008 スピン方位制約条件付きノンコリニア密度汎関数法の開発 web notes
2009 高速二電子積分法の開発 JCP 130, 124114
2009 高速球ベッセル変換法の開発 CPC 181, 277
2010 高精度非平衡グリーン関数法の開発 PRB 81, 035116
2011 高並列化O(N2~)法の開発
PRB 82, 075131
2011 O(N)厳密交換汎関数法の開発 PRA 83, 032515
2012 O(N)法+有効媒質法の開発 JCP 136, 134101
2014 O(N)法の超並列化手法の開発 CPC 185, 777 (2014)2014 3D-FFTの超並列化手法の開発 CPC 185, 153 (2014)2014 大規模NBO解析手法の開発 JCP 140, 244105 (2014)

LCPAO法(Linear-Combination of Pseudo Atomic Orbital Method)
一電子KS軌道を擬原子軌道の線形結合で表現する方法

特徴:
・原子様基底を用いているため、計算結果の物理的・化学的解釈が容易。
・基底関数系は厳密には完全系ではないため、完全な収束解を得ることは難しい。
しかし物理的起源に基づく基底のため、注意深く構成された場合には比較的少ない基底関数で実質的に収束解を得ることが可能。
・基底が局在しているため、硬い擬ポテンシャルの使用が可。高精度な擬ポテ
ンシャルを作成することが可能。
・基底の局在性のため、ハミルトニアン行列の計算量と保存メモリ量がO(N)となる。
・様々なO(N)法のアイデアと整合性が良い。

実装: 全エネルギー (1)全エネルギーは５つの項の和で与えられる。数値的に高精度な実装を行うためにはクーロンエネルギーの再定式化が重要である。

運動エネルギー
外場とのクーロン相互作用エネルギー

Hatreeエネルギー

交換相関エネルギー
核間クーロンエネルギー
TO and H.Kino, PRB 72, 045121 (2005)

実装: 全エネルギー (2)
三つのクーロンエネルギーの和を再配分することで、厳密に三つの新しい項を導入できる。
s
中性原子エネルギー
短距離相互作用でかつ2中心積分に分解できる。

差電子ハートリーエネルギー
長距離相互作用であるが、小さな寄与である。
遮蔽核間反発エネルギー
短距離相互作用でかつ２中心相互作用。
差電子密度

中性原子ポテンシャル

実装: 全エネルギー (3)したがって、全エネルギーは次式で与えられる。

それぞれの項の計算は計算精度と効率を考慮し、適切な数値積分法が適用される。

}

フーリエ変換を用いて、q空間で積分}

等間隔グリッド(FFTグリッド)を用いて実空間で積分

細グリッドを用いて実空間で積分

2中心積分の計算方法
基底関数のフーリエ変換

例えば、重なり積分は次式の様に計算される。

角度方向は解
析的に積分できるため、動径方向のみ数値積分すれば良い。

力の計算

Easy calc.
See the left

数値グリッドを用いているにも関わらず、任意のグリッド幅において、解析的に力の計算が実行できる。

FFTグリッドに対するカットオフエネルギー
二つのエネルギー項 Eδee + Exc はFFTグリッドを用いて実空間上で計算される。グリッドの細かさは平面波のカットオフエネ
ルギーで与えられる。

カットオフエネルギーは次式でグリッド細かさと関連付けられる。

プリミティブ擬原子基底関数1.

閉じ込めポテンシャルの下で原子のKS方程式
を自己無撞着に解く。

s-orbital of oxygen

2. 次にこの閉じ込めポテンシャル下でノルム保存擬ポテンシャルを作成。
3. 各角運動量成分lに対して基底状態と励起状態をstep2で得られた擬ポテン
シャルに対して数値的に求める。

多くの系において、精度と計算効率は二つのパラメーターで制御される。

カットオフ半径基底関数の数PRB 67, 155108 (2003)
PRB 69, 195113 (2004)

基底関数に対する収束性
基底のカットオフ半径と基底を増やすとエネルギーと構造パラメーターは系統的に収束している。数学的には厳密ではないが、二つのパラメーターは変分パラメーターと見なし得る。

分子固体

二原子分子の基底状態の計算
擬原子基底による二原子分子の基底状態計算結果

局所密度近似(LDA)における全ての“成功”と“失敗”がDNP程度の基底サイズで再現できる。DNP=二つの価電子基底+一つの分極関数基底数値擬原子基底関数の変分最適化

LCPAO法の計算精度・計算速度は基本的に基底のカットオフ半径
基底の数で制御される。しかし基底関数の形自体を最適化出来るならば、計算量を増やすことなしに、計算精度を向上させることが出来る。

Ozaki, PRB 67, 155108 (2003)

基底関数の変分最適化
1電子KS波動関数

縮約軌道

縮約係数aを固定し、cに関して変分をとれば

a: 縮約係数→
変数cがaに対する従属変数であると見なし、さらに変数cの観点でKS方程式
が自己無撞着に解けていると仮定すると、次式が得られる。

Ozaki, PRB 67, 155108 (2003)

基底関数の最適化手順
1. 典型的な化学的環境を選択

2. 動径波動関数を変分最適化

3. 部分空間内で回転し冗長成分を取り除いた後に統合

××

Δ-factor
FLAPW+LO法であるWien2kの体積-エネルギー曲線の収束解を参照として、そこからのエネルギーのずれΔによって各コードの数値
精度を見積もる。(V0から±6%の範囲のエネルギー計算からV0,
B0, B1を求め、Birch-Murnaghan状態方程式を利用し、Δ-factor
を計算)

Lejaeghere et al., Critical Reviews in Solid State and Materials Sciences 39, 1-24 (2014).

Δ-factorによる各コードの比較

http://molmod.ugent.be/deltacodesdft より引用。

OpenMXによる計算の実際:計算上で留意すべき点• カットオフエネルギーの選択• エネルギーカーブの計算• SCF計算の理論と実際• 構造最適化の理論と実際• 基底関数の選択• 表面状態やfloating状態• Overcompleteness• リスタート• バイナリモードでの出力カットオフエネルギーの選択
scf.energycutoff

200# default=150 Ryd

実空間を離散化するFFTグリッドはscf.energycutoffで指定する。FFTグリッド
はPoisson方程式の解法、Eδee, Excの計算に使用する。
多くの系で200 Ryd程度で収束する。
擬ポテンシャルにセミコア状態を含む場合には、大きめのカットオフを取った方が安全。
300-400 Ryd. 程度。

メモリの増加量:O(E3/2)メタン分子の全エネルギーカットオフエネルギーの選択
H2Oの最適化構造
水分子の最適化構造のscf.energycutoffに対する依存性。
180Ryd.でほぼ収束していることが分かる

エネルギーカーブの計算

格子定数を変化させてエネルギー曲線を計算するには以下のkeywordが使用できる。
MD.Type
MD.EvsLC.Step
MD.maxIter

EvsLC
0.4
32

#
# default=0.4%# default=1

scf.energycutoffによってグリッドを決定すると、格子定数の増加に応じて、不連続にグ
リッド数が変化する。このグリッド数の突然の変化はエネル
ギーカーブに不連続なジャンプを生じさせる。これを避けるた
めにはグリッド数を固定して計算を行うべきである。次のkeywordで明示的にグリッド数が指定できる。
scf.Ngrid32 32 32

それぞれの数字はa-, b-, and c-axesに対するグリッド数。scf.energycutoffとscf.Ngridをどちらも指定した場合にはscf.Ngridが優先される。

MD.EvsLC.flag
1 1 0# a, b, c-軸に対して1: 格子定数を変化, 0: 格子定数を変化させない。

SCF計算: 電子密度の混合
KS有効ポテンシャルはρの関数である。
しかし本来、ρはKS方程式から求められる。

単純混合法
入力と出力の電子密度を単純に足し合わせて、次の入力電子密度とする。

大きなバンドギャップのある系、小さい系に有効な方法

OpenMXで利用可能なSCF収束法
Simple mixing (Simple)
関連するkeywords: scf.Init.Mixing.Weight,
scf.Min.Mixing.Weight, scf.Max.Mixing.Weight
Residual minimization method in the direct
inversion iterative subspace (RMM-DIIS)
関連するkeywords: scf.Init.Mixing.Weight,
scf.Min.Mixing.Weight, scf.Max.Mixing.Weight,
scf.Mixing.History, scf.Mixing.StartPulay
Guaranteed reduction Pulay method (GR-Pulay)
関連するkeywords: scf.Init.Mixing.Weight,
scf.Min.Mixing.Weight, scf.Max.Mixing.Weight,
scf.Mixing.History, scf.Mixing.StartPulay
Kerker mixing (Kerker)
関連するkeywords: scf.Init.Mixing.Weight,
scf.Min.Mixing.Weight, scf.Max.Mixing.Weight,
scf.Kerker.factor
RMM-DIIS with Kerker metric (RMM-DIISK)
関連するkeywords: scf.Init.Mixing.Weight,
scf.Min.Mixing.Weight, scf.Max.Mixing.Weight,
scf.Mixing.History, scf.Mixing.StartPulay,
scf.Mixing.EveryPulay, scf.Kerker.factor• 5つのSCF収束法が利用可能。
• 多くの系でRMM-DIISKがベスト。
• RMM-DIIS、RMM-DIISKはどち
らもPulay混合法であるが、前者は密度行列、後者は電子密度のフーリエ変換に対してPulay混合を適用する。
• 非常に収束が難しい系はKerker
がベター。
• plus Uの場合にはRMM-DIISの方がRMM-DIISKより収束性に優れることもしばしばある。

RMM-DIISKによるSCF収束
アイデア
過去の(フーリエ空間での)残差ベクトルの線形和のノルムを最小化する。

KerkerfactorFの最小化条件から次式が得られる。

|q|が小さい成分(長波長成分)に重きを置く。

“最適”な入力電子密度は次式で与えられる。

G.Kresse and J. Furthmeuller,
PRB 54, 11169 (1996).

様々な電荷密度混合法の比較
シアリック酸

Pt13クラスター

Pt63クラスター

いずれの場合もRMM-DIISもしくはRMM-DIISKが高速である。

RMM-DIISKのkeywords
RMM-DIISKは以下のkeywordsで制御される。
scf.Init.Mixing.Weightscf.Min.Mixing.Weightscf.Max.Mixing.Weightscf.Mixing.History
scf.Mixing.StartPulay
scf.Mixing.EveryPulay
scf.Kerker.factor
•
•
•
•

•
•

0.01
0.01
0.2
40
20
1
5.0

scf.Mixing.StartPulayはRMM-DIISKをスタートするSCFステップを指定する。
それ以前のステップはKerker法が適用される。RMM-DIISK法は準ニュートン法であるから、ある程度の収束が達成された後に、乗り換えるべきである。さもないと収束しない場合がある。
RMM-DIISKにおける履歴数はscf.Mixing.Historyで指定する。
scf.Init.Mixing.Weight, scf.Min.Mixing.Weight, scf.Max.Mixing.Weight はKerker法を適用される際の混合比の初期値、最小値、最大値である。収束状況をモニターし、この範囲で自動的に調整される。
Kerker因子はscf.Kerker.factorで与えらえる。指定しない場合には自動的に最適値が推定され、使用される。
RMM-DIISK法に乗り換えた後に、RMM-DIISKによる混合を実施する頻度を
scf.Mixing.EveryPulayで指定する。1の場合には毎回、実施することに対応。

構造最適化
分子や固体の性質を精密に調べるためにはエネルギーの最小値を与える構造(最適構造)を計算する必要がある。

最急降下法(Steepest decent method)
エネルギーの勾配方向に少しづつ移動する。

1

5

3
2

4

6

0
あまり効率的ではない。

ニュートン法による構造最適化
エネルギーをTaylor展開する

座標の微分を計算

これをΔxについて解けば、もしHessian行列Hが求められるならば、この方法は効率的である。
しかし一般の場合に、Hを求めることは困難である。

準ニュートン法による構造最適化
OpenMXにおける構造最適化は準ニュートン法に基づいている。
OpenMX Ver. 3.7においては４つの最適化方法が実装されている。

近似ヘシアン行列Hの計算手法
DIIS

BFGS

RF(rational function)

H=I

BFGS

BFGS+RFEF(eigenvector following)
BFGS plus monitoringof eigenvalues of H

Broyden-Fletcher-Goldfarb-Shanno (BFGS) 法

赤枠で囲まれた内積が正であるならば、Hの正定値性は保たれる。

構造最適化法の比較
分子固体• ベンチマーク計算からはEF法とRF法の性能が良いことが分かる。
• 構造に自由度が大きい分子系は収束しにくい。
このベンチマーク計算の入力ファイルと出力ファイルは以下のディレクトリにある。
openmx3.7/work/geoopt_example".

基底関数の指定基底関数は以下の様に指定する。

• O7.0はO7.0.paoに対応する。
• -s2p2d1は動径関数として、s軌道を二つ、p軌道を2つ、d軌道を1つ選択することを意味する。
• この時、酸素Oに対して2×1+2×3+1×5=13基底が割り当てられることになる。

O7.0.paoやO_PBE13.vpsの保存先は以下のkeywordで指定する。
DATA.PATH../DFT_DATA13 # default=../DFT_DATA13

Default値は../DFT_DATA13である。

基底関数の選択: H2Oの場合
Database(2013)においてH7.0.paoとO7.0.paoをクリックすると以下の情報が得られる。
http://www.jaist.ac.jp/~t-ozaki/vps_pao2013/H/index.html

H7.0.pao3

http://www.jaist.ac.jp/~t-ozaki/vps_pao2013/O/index.html

O7.0.pao

1

6

4

5

2

7

8

擬原子状態で低エネルギーの軌道から選択するとH7.0-s2p1, O7.0-s2p2d1となる。

基底関数の選択: Siの場合 (1)Si7.0.pao
14
8Si7.0.pao中の固有エネルギーの低い方から順次、基底に含める。

25

構造特性に関しては、
Si7.0-s2p2d1で十分
である。

3
7

Si7.0-s3p2d2f1とSi8.0-s3p2d2f1の比較から、カットオフ
半径は7.0(a.u.)で収束していることが分かる。

6

基底関数の選択: Siの場合 (2)バンド構造の比較からもSi7.0-s2p2d1でほぼ収束していることが分かる。

より完全な収束解を得るためにはSi7.0-s3p2d2f1(Si7.0-s3p3d2f1)を用いた方が良いが、大規模系を取り扱う場合や、広範囲な構造探索を行うためにはSi7.0s2p2d1は精度とコストのバランスの取れた基底関数である。

フローティング状態 in 3C-SiC
C6.0-s2p2
Si8.0-s2p2

C6.0-s2p2d1
Si8.0-s2p2d1

C6.0-s3p3d2f1
Si8.0-s3p3d2f1

Wien2k (FLAPW)

 d軌道を含めることが、バンド構造を再現する上で重要である。
 少なくとも5eVまでのバンド構造はs2p2d1基底で再現されている。

基底関数重なり誤差 (BSSE)

ベンチマーク計算から分子系に対するBSSEは～0.5 kcal/mol

仕事関数の計算
fcc Al (111) surface
真空層にゴースト原子を置くことで仕事関数の高精度計算が可能である。
(本計算はDatabase(2006)の基底関数を用いて行われたため、Database(2013)では状況が改善されている可能性がある。)

仕事関数
Exp : 4.24 eV
PW : 4.09 eV

4.08 eV
4.07 eV

3.15 eVY. Morikawa et al.,
PRB 69,041403(R) (2004);
TABLE I.

2

1

By Jippo and Ohfuti (Fujitsu)

1

2

Al(111) fcc 1x1 6L, GGA-PBE, Fix to the equilibrium lattice parameter
Al6.0-s1p2d1, Al_PBE(Database ver. 2006)
E6.0-s1p2d1, E(Database ver. 2006)

基底関数のOvercompletenessfcc-Agの全エネルギー

fcc, hcp, bcc等の密な構造に対して、多数の基底関数を使用した際に、overcompletenessと呼ばれる数値不安定性が発生する場合がある。

Overcompletenessの原因
s軌道の基底を持つ水素分子を考える。
原子間距離をdとしてd→0とする。
dH’の計算の際にv1のために桁落ちが起きる。その誤差が1/√λ1で増幅される。

Overcompletenessへの便宜的な対策
Overcompletenessが発生した際に、中性ポテン
シャルのプロジェクター展開をオフにすることで、この問題を回避できることが多い。
scf.ProExpn.VNA

off# on|off, default = on

KSハミルトニアン自身の精度が関連していると考えられるが、この理由に関しては、完全解明には至っていない。

リスタート• 計算が終了したにも関わらずSCFが収束していない場
合、以下のkeywordを含めることで最終の電子密度分布を用いてSCF計算をリスタートすることが出来る。
(ただしLDA+U法の場合には機能しない)
scf.restarton# on|off,default=off• 構造最適化や分子動力学計算が途中で終了した場合、最終座標や最終速度を含んだ入力ファイル(*.dat#)が利用できる。このリスタート用の入力ファイルを用いて、計算を途中からリスタート出来る。

バイナリモードでの出力• 大規模な計算を実行し、複数のcubeファイルを出力する
場合、ファイルシステムの性能によっては長時間を要す
る場合がある。この場合、以下のキーワードを用いてバ
イナリモードで出力することで出力時間を低減することが可能である。
OutData.bin.flag

on# default=off, on|off• バイナリファイルはbin2txt.cを用いて、テキストファイ
ルに変換できる。
コンパイル方法: gcc bin2txt.c -lm -o bin2txt
使用方法: ./bin2txt *.bin大規模計算に関して
512コアを使用して1000原子系の通常の対角法による(DNP基底を用いた)計算が可能である。runtestL2のベンチマーク計算の結果を以下に示す。
$ mpirun -np 128 openmx -runtestL2 -nt 4

ファイル名の数字は系に含まれる原子数を示す。Elapsed timeは1 MDあたりの計算時間 (秒)である。一日程度で構造最適化が可能。
マニュアル96pも参照のこと。

「京」上でのO(N)法の並列化効率

系: 131072原子からなるダイヤモンド格子

13万1072コアを使用しておよそ68％の並列化効率実習内容 (2)1. 構造最適化
Methane2.datを用いて、構造最適化を行う。
マニュアル63pを参照のこと。

2. 状態密度の計算
Cdia.datを用いて、状態密度を計算する。
マニュアル76pを参照のこと。

全ての入力ファイル*.datはopenmx/work/
以下にあります。

余裕のある方は以下の計算にもチャレンジして下さい。

3. Nudged elastic band (NEB)法による反応障壁の計算
C2H4_NEB.datを用いて、2個のエチレン分子の環状付加反応によるシクロブタンの生成過程における反応障壁を求める。
マニュアル168pを参照のこと。

4. 炭素鎖の電子透過率の計算
Lead-Chain.dat, NEGF-Chain.datを用いて、電子透過率を計算する。マニュアル129pを参照のこと。

構造最適化 (Methane2.datを用いて)メタン分子の炭素原子のx 座
標を右のように0.3 Å に変更
し、これを初期構造として構
造最適化する。
ジョブの投入方法$ mpirun –np 5 openmx Methane2.dat –nt 2

この計算は最急降下法 (SD)を用いた場
合で、右の履歴で構造が収束する。
XCrySDenを用いるとmet2.mdを可視化できる。その際にはxyzファイルとして読み込むこと。
マニュアル63pも参照のこと。
構造最適化ステップ

状態密度の計算 (Cdia.datを用いて)
Kohn-Sham 固有値の状態密度（DOS）は次の２つの段階を経由して計算される。
(1) SCF計算
ダイヤモンドのSCF計算を実行する。
ジョブの投入方法$ mpirun –np 2 openmx Cdia.dat –nt 2(2) DOS計算(2.a) openmx/source中でDosMainをコンパイル$ make DosMain(2.b) workに移動してDosMainを使ってDOS計算$ ./DosMain cdia.Dos.Val cdia.Dos.Vec
以下、マニュアル76pを参照のこと。