実空間差分法に基づく第一原理計算プログラム RSPACE の使い方

(電子状態計算編)

平成 29 年 4 月 1 日 初版作成

筑波大学計算科学研究センター
v8e

はじめに

1

本手引書は，筑波大学計算科学研究センター小野研究室で開発した「実空間差分法に基づく第一原理計算プログラム RSPACE」の基本的な使い方をまとめたものである．本書で取り上げている，例題の計算はすべて Intel Fortran Compiler と Sun Grid Engine をインストールした
PC クラスターで行うことを前提に記述してある．1.1 節に挙げる他の計算機で計算を行う場合，
または 1.1 節に載っていない計算機で計算を行う場合は 1.2 節の連絡先まで連絡してください．

1.1

動作確認機種

本プログラムは以下のマシンに関して動作確認を行っている．
・京コンピュータ
・筑波大学計算科学研究センター COMA・東京大学物性研究所 sekirei

1.2

連絡先

本プログラムの不具合，及び本手引書の誤植に関する連絡先〒 305-8577 茨城県つくば市天王台 1-1-1
筑波大学計算科学研究センター量子物性研究部門
小野倫也
e-mail:ono@ccs.tsukuba.ac.jp
TEL 029-853-4273 FAX 029-853-4273

2

本プログラムの概要

本章では，「実空間差分法に基づく第一原理計算プログラム RSPACE」を用いてシミュレー
ションを行うにあたり，パラメータ決定の際に必要となるフローチャートと，本プログラム固有の記述形式についてまとめている．第一原理計算で，一般的に用いられている単語や記述形
式については省略しているため，そのような単語などで分からないものがあれば，他の文献を参照のこと．

2.1

単位系について本プログラムで扱う単位は全て原子単位系 (a.u. または Hartree) である．原子単位系は，以下のように定義されている．
電子の質量
電荷素量
プランク定数真空の誘電率

1

m
e2
h̄
ε0

→
→
→
→

1
1
1
1
4π

したがって次の物理量は，原子単位系では下のようになる．
長さ
質量
時間
エネルギー

2.2

1 (a.u.)1 (a.u.)1 (a.u.)1 (a.u.)→ 5.2918 × 10−11 (m)→ 9.1094 × 10−31 (kg)→ 2.4189 × 10−17 (sec)→27.212(eV)

境界条件について本プログラムは，x, y, z すべて非周期系，x は周期系で y, z は非周期系，x, y は周期系で z は非周期系，x, y, z すべて周期系の 4 種類の境界条件の設定が可能である．計算するモデルを考慮し，適切な境界条件を選ぶこと．なお，アスキー入出力データは，すべての境界条件で互換性があるが，バイナリ入出力データは，x, y, z すべて非周期系のみ波動関数の虚数部が無いため，
他の境界条件との互換性がない．

2.3

フローチャート

本プログラムのフローチャートを，図 1 に示す．図に出てくる単語は，パラメータ設定の際に必要となるので，参考にするとよい．

2.4

擬ポテンシャルについて

本プログラムは，価電子と内核電子の相互作用を，主に Projector augmented-wave 法で取り扱っているが, 一部の元素でノルム保存型擬ポテンシャルを用いる．このノルム保存型擬ポテン
シャルは，物質・材料研究機構 小林一昭氏の作製した擬ポテンシャルデータベース NCPS97 を本プログラム用に入力ファイルを改編したものである．擬ポテンシャルのデータベース NCPS
シリーズについては，http://www.geocities.co.jp/Technopolis/4765/INFO/guidee.html を参照
のこと．

2.5

並列処理の表記について

本プログラムは，大規模計算を高速に行うことを目的に開発されている．そのためには，複数の計算ノードを用いて並列計算を行う．本節では，並列化方法の詳細については説明を省略するが，以後の説明のため，本書におけるの並列化に関する言葉の定義を行う．
並列方式は，以下にあげる 2 つに大別できる．

1. プロセス並列 計算を複数のプロセスに分けて並列化する．各プロセスが使用するメモリ空間に重複は無く，ほとんど独立に計算を行う．プロセス間でデータの相互参照が必要な場合は，データをネットワークケーブル等を経由して他のプロセスに転送する．

2

開始
読み込み
計算条件(read.inp)
原子座標(atom.xyz)
初期電荷分布(rho.*****)
初期波動関数(wfc.*****)

原子核からのポテンシャルを計算

静電電子間ポテンシャルを計算(Poisson方程式をCG法で解く)
N
収束？
Y
交換相関ポテンシャルを計算

Y

N
部分空間で
対角化する？

NN
固有値の並べ替えする？

Yself consistent loop

CG or DIIS法を繰り返す？

Y
固有値の
並べ替え

Rayleigh-Littz法

N
電荷分布
が収束？

電荷分布
を更新

Y
原子に働く力を計算
し原子を移動させる(原子構造最適化)

Y
構造最適化を続ける？
N
計算結果の出力
原子座標(atom.xyz)
全エネルギー etc. (mdresult.dat)
電荷分布(rho.*****)
波動関数(wfc.*****)

終了

図 1: フローチャート

3

structure optimization loop

Kohn-Sham方程式をCG or DIIS法で解く

2. スレッド並列 1 つのプロセスが，複数の演算器 (コア) を用いて計算を行う．メモリ空間は，プロセス単位で確保するため，各スレッドでメモリ空間を共有する．したがって，スレッド間でのデータの相互参照が必要な場合でも，メモリ空間を共有しているためプロセス並列に比べ容易であるが，全計算で使用可能なスレッド数やメモリ量は，ハードウエアに搭載されている演算器数 (コア数) やメモリ量に制限される．
本プログラムは両方の並列処理が行えるようにコーディングされている．ただし，並列処理の方法を明示的に区別するため，実行条件を記述する際の例を以下に記す．
プロセス並列化する．· · · 1. の方式で並列化を行う．
スレッド並列化する．· · · 2. の方式で並列化を行う．
並列化する．· · · 1. と 2. の両方の方式で並列化を行う．
m プロセスで並列化する．· · · 1. の方式で m 個のプロセスで並列化する．
n スレッドで並列化する．· · · 2. の方式で n 個のスレッドで並列化する．
n×m で並列化する．· · ·1 プロセスあたり n 個のスレッドで，m プロセス使用する．合計 n×m個の演算器 (コア) を使用して並列化する．このタイプの並列化をハイブリッド並列化という．
本文中，特に断りのない限りは上の例に従う．

3

ソースプログラムのコンパイル3.1

ソースプログラムのコピー

RSPACEyymmdd.zip を計算機上の適当なフォルダに展開する．ここで yymmdd は，プログラムを配布した年月日である．展開したディレクトリの中に v8combi と pspaw というディレクトリ
が作成される．前者はソースファイルが格納されたディレクトリ，後者は練習用の擬ポテンシャ
ルが格納されたディレクトリである．

3.2

makefile の設定

ソースプログラムをコンパイルするにあたり，v8combi の中の makefile を変更する．展開したフォルダの中に makefile がある．初期値は，Intel Fortran Compiler をインストールした PC
クラスター用に設定しており，makefile の 4 行目がsystem = xeonになっている．Intel Fortran Compiler をインストールした PC クラスターを使うのであれば，
このまま変更する必要はない．1.1 で紹介した他のコンピュータを使うには，表 1 のように設定
する．

3.3

コンパイル

ディレクトリ v8combi の下で,>make1
1“>”は，プロンプトである．入力しないように．

4

表 1: makefile の設定
計算機
京コンピュータ
筑波大学計算科学研究センター COMA
東京大学物性研究所 sekirei

オプション
kei
coma
sekireiと入力する．コンパイル後，kukan8，ele，ldos，orbcharge というファイルが作成されていることを確認する．

4

簡単な分子の計算を行う

まず最初に，本節では，簡単な例題として，1×1 で非周期境界条件下における CO 分子の電子状態計算を行う．

4.1

ディレクトリの作成

まず，各自のホームディレクトリ (ここでは ∼/とする) の下に実行用ディレクトリ ∼/CO を作り，さらに ∼/CO の下にディレクトリ pspaw(∼/CO/pspaw) を作成する．そして, ソースファイル
をコンパイルしたディレクトリ v8combi から ∼/CO に，kukan8，read.inp，atom.xyz, job.sh をコピーする．次に，RSPACEyymmdd.zip を展開したディレクトリの pspaw ディレクトリから，炭素の擬ポテンシャルファイル paw.006 と酸素の擬ポテンシャルファイル paw.008 を，∼/CO/pspawにコピーする．

4.2

パラメータファイル read.inp の作成∼/CO の中の read.inp を，次の例のように書き換える．赤字で書いた部分が，基本計算で変更する可能性のあるパラメータである．基本計算用の各パラメータの意味は，後で説明する．
read.inp ファイルの作成例 (初めの 69 行のみ)
1: ![ndisp]
2: 66
3: ![nprocx],[nprocy],[nprocz],[nprock]
4: 1 1 1 1
5: ![xmax],[ymax],[zmax]
6: 10.0d0 10.0d0 10.0d0
7: ![nxmax],[nymax],[nzmax]
8: 30 30 30
9: ! [npxmax],[npymax],[npzmax]
10: -8 8 26
11: ![nso],[socang(1)],[socang(2)]12: 0 0.0d0 0.0d0
13: ![nsym]
14: 0
15: ![neigmx]
16: 8
17: ![natom]
18: 2
19: ![num atcell]
20: -1
21: ![num ppcell]
22: -1
23: ![num ppcell d]
24: -1

5

25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:

![gmax]20.0d0![nperi],[npopshow]
0 0![numkx],[numky],[numkz],[ksym],[kband]1 1 1 0 0
!k point (write all k points)![X], [Y], [Z], [weight]
0.00d0 0.00d0 0.00d0 1![cexco]
VWN![nspv]1![epsvh],[epssd],[ratio diis],[eps scf]
1.0D-8 1.0D-8 0.3d0 1.0D-8![ncgmin],[ncgmax],[ncgres]1 800 801![nprecon cg],[nprecon diis]
1 0![nsdmax],[nkscg]
4 1![ndiismax],[ncgscf],[nretcg]
4 200 2000![nrrz],[nchange]
1 1![eta],[looplimit],[nbrydn],[etamag(1)],[etamag(2)]0.20D0 200 10 0.5d0 2.0d0![tmstep],[nmd start],[nmd end],[ngdiis],[sconst]
0.0D0 0 1 1 0.0d0![biasx],[biasy],[biasz]
0.0D0 0.0D0 0.0D0![tf],[tfmin],[tfmax]
0.003D0 -10.0D0 30.0D0![nmdxyz]
0![chrgd],[npolcon],[polconocc]
0.0D0 0 2.7d0![endjel],[chrjel]
0.0d0 0.0d0![fcut]1.0D-4![npre],[nevhist]
1 3![northo]1

プログラムのバージョンによっては，サンプルが変更されている可能性もあるので，すべてのパラメータをチェックすること．また，各パラメータの意味を以下に記す．

6

nproc{x,y,z,k} 並列計算制御用．5 章参照．本章の計算では，全て 1 にしておく．
{x,y,z}max
計算領域設定用．-{x,y,z}max〜+{x,y,z}max が計算領域になる．非周期境界を設定する場合には，atom.xyz(4.3 参照) 節
で決めたすべての原子座標が境界面から 9〜15(a.u.) 以上の距離を保つように，計算領域を設定する．
n{x,y,z}max
グリッド数．必ず整数を選ぶこと．2×nxmax が x 方向のグリッ
ド数，xmax/nxmax がグリッド幅 dx になる．グリッド幅の設定については Appendix を参照のこと．
neigmx
計算する波動関数の数．遷移金属を含まない場合は原子数 ×4，含む場合は原子数 ×9 を推奨する．
natom
原子数．
nperi
境界条件の設定用．非周期境界条件:0，x, y 方向非周期，z 方
向周期境界条件:1，x 方向非周期，y, z 方向周期境界条件:2，x, y, z 方向周期境界条件:3．
numk {x,y,z}各方向の k 点数．x, y, z すべての方向に非周期境界条件を課した場合は 1 にしておく．周期境界条件で複数の k 点を計算する場合は，33 行目以降の行に k 点の情報を並べて書く．また，基礎編の計算においては，numkx × numky × numkz の積を numkx に入れ，他を 1 としても良い．
nspv
交換相関項ポテンシャルにスピンの効果を無視:1，考慮:2．
eta
電子状態の自己無撞着計算における電子密度の更新比．
looplimit
電子状態の自己無撞着計算における反復計算の最大回数．
nbrydn電子密度混合時の Broyden 法の履歴数．1 に設定すると単純な電子密度混合になる．
tmstep原子構造最適化時の時間ステップ．
nmd end
原子構造最適化時のステップ数．
ngdiis
原子構造最適化時における Pulay 法の履歴数．1 に設定すると単純に力の方向に沿って原子を動かす．
npre
初期電子状態の読み込み制御．読み込む:0，新規作成:1．

4.3

原子座標ファイルの作成

次に，原子座標などの情報を記述したファイル atom.xyz を作成する．本章の計算では，以下にのように設定する．各数字の意味は，左から原子の x 座標，y 座標，z 座標，原子番号，x 方向の構造最適化のスイッチ，y 方向の構造最適化のスイッチ，z 方向の構造最適化のスイッチ，原子の質量 [原子量 × 陽子の質量 (1838 a.u.)] である．構造最適化のスイッチが 0 の場合，その方向へは原子を動かさない．本章では電子状態のみを計算するので，すべてのスイッチを 0 にしておく．原子が複数個ある場合は，原子の個数分だけ，下に行を連ねていく．したがって原子が N 個あれば，このファイルは N 行必要である．
atom.xyz ファイルの作成例 (初めの 3 行のみ)
1:! [x], [y], [z], [atom number], switch [x], [y], [z], [weight], switches [soc], [pp], [na]2:0.5500000000000000D+00 0.5500000000000000D+00 0.5500000000000000D+00 6 0 0 0 22056.00 00 1 1a7

3:-0.5500000000000000D+00 -0.5500000000000000D+00 -0.5500000000000000D+00 8 0 0 0 29408.00 00 1 2a

CO 分子の平衡原子間距離は 2.13〜2.14(a.u.) であるが，本書では，後に原子構造最適化計算
により平衡原子間距離を求めるため，実験値よりも短めに設定しておく．

4.4

実行スクリプトの作成

プログラム実行用に，実行スクリプトを作成する．ここでは，Sun Grid Engine によるジョブ
管理システムをインストールした PC クラスタでの実行を想定したジョブスクリプトについて紹介する．他の計算機システムで実行する際には，計算機システムに合ったスクリプトファイ
ルを準備すること．
このジョブスクリプトでは，2 行目で使用するコアの総数 (2.5 節の n × m)，3 行目で 1 プロ
セス当たりのスレッド数 (2.5 節の n) を入力する．プロセス数 (2.5 節の m) については，スクリプトの中で自動的に計算する．
job.sh ファイルの作成例
1: #!/bin/csh
2: #$ -cwd -pe mpi 1
3: setenv OMP NUM THREADS 1
4: mv $USER.$PE.$JOB ID.1 hosts
5: set TEST = ""6: foreach aaaa (‘cat hosts‘)
7: @ I++
8: set TEST = ($TEST ‘echo $aaaa‘)
9: end10: rm hosts
11: set i=0
12: while ( $i < $NSLOTS )13: @ i = $i + $OMP NUM THREADS
14: echo $TEST[$i] >> hosts
15: end
16: @ i = $NSLOTS / $OMP NUM THREADS
17:
18: mpirun -np $i -machinefile hosts ./kukan84.5電子状態計算プログラムの実行

プログラムの実行は qsub というコマンドで，ジョブを投入する．
>qsub job.sh
順番待ちさえなければ，電子状態計算プログラムが即時に実行される．電子状態計算プログラ
ムが終了すると，rho.00000，vht.00000，wfc.00000，mdresult.dat が作成される．
>cat mdresult.datと，入力すると計算結果を確認できる．以下に mdresult.dat の内容を示す．
computed by kukan8e 03/29/2017-01========== data from read.inp ==========
ndisp
:
66
nprocx:1
nprocy:1nprocz:1
xmax:
10.0000000000000
ymax:
10.0000000000000
zmax:
10.0000000000000
nxmax
:
30
nymax
:
30
nzmax
:
30
npxmax:-8
npymax:88

npzmax:nso:
socang:nsym:
neigmx:natom:
num_atcell :num_ppcell :num_ppcell_d:
gmaxps:nperi:
npopshow
:
numkx
:
numky
:
numkz
:
ksym
:
kband
:
cexco
:
nspv
:
epsvh
:
epssd
:
ratio_diis :eps_scf:
ncgmin:ncgmax:ncgres:nprecon_cg :nprecon_diis:
nsdmax:nkscg:
ndiismax:
ncgscf:nretcg:nrrz:nchange:eta:
etamag:looplimit:
nbrydn:tmstep:nmd_start:
nmd_end:ngdiis:sconst
:
biasx
:
biasy
:
biasz
:
tf
:
tfmin
:
tfmax
:
nmdxyz:chrgd:
npolcon:polconocc:
endjel:chrjel
:
fcut
:
npre
:
nevhist:northo:eps:
eps_eig_diis:
alambda_diis:
alambda_min :alambda_max :nradmx:nrprjmx:
nprjmx:lsphel:nlmax:
lrhomx:nfiltyp:26
00.000E+00 0.000E+00
0
8
2
2
2
2
20.0000000000000
0
0
1
1
1
0
0
vwn
1
1.000000000000000E-008
1.000000000000000E-008
0.300000000000000
1.000000000000000E-008
1
800
801
1
0
4
1
4
200
2000
1
1
0.200000000000000
0.5000E+00 0.2000E+01200
100.000000000000000E+000
01
10.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
3.000000000000000E-003-10.0000000000000
30.0000000000000
00.000000000000000E+000
02.70000000000000
0.000000000000000E+000
0.000000000000000E+000
1.000000000000000E-004
1
3
1
1.000000000000000E-016
1.000000000000000E-014
0.500000000000000
0.100000000000000
1.00000000000000
1502
6
18
8
4
25
1

0.000E+00

9psctoff:
1.05000000000000
nqmx
:
400
psftrad:10.0000000000000
psctrat:1.00000000000000
psext:
2.00000000000000
filpp: 2.000000000000000E-002
rctpcc:2.00000000000000
nmesh
:
2
npmesh:2veta: 0.200000000000000
new_pwx:0
new_pwy:0
new_pwz:0
new_rsx:0
new_rsy:0
new_rsz:0
nint1dmax
:
1000
ndec
:
12
nf
:
4
nfdg
:
4
nfh
:
4
nrasy
:
0
nrasz
:
0
nrepx
:
1
zs_pre: 0.500000000000000
pol_pre:1.50000000000000==========================================
Filtering of pseudopotentials is implemented by Fermi distribution.
Filtering of pseudopotentials is implemented by Fermi distribution.
npxmax, npymax, and npzmax are updated below.
========== data computed by listvecdim ==========
npxmax:8
npymax:8
npzmax:8====================================================
cell size [lx,ly,lz (a.u.)]
0.200000000D+02
0.200000000D+02
0.200000000D+02
coarse grid spacing [hx,hy,hz (a.u.)]
0.333333333D+00
0.333333333D+00
0.333333333D+00
dense grid spacing [hx,hy,hz (a.u.)]
0.166666667D+00
0.166666667D+00
0.166666667D+00
cutoff energy of wave functions [gmax_coarse,gmax_dense (Ry)]
0.888264396D+02
0.355305758D+03
number of electron
chargespin
0.100000000D+02
0.000000000D+00
0.100000000D+01
electric field [ex,ey,ez (a.u./a.u.)]
0.000000000D+00
0.000000000D+00
0.000000000D+00
time step (a.u.) (fs)0.000000000D+00
0.000000000D+00
thermo (k)kt (a.u.)
0.947268709D+03
0.300000000D-02
cufoff of force acting on atoms (a.u.)
0.100000000D-03
charge mixing ratio of scf
0.200000000D+00
nrasy,nrasz
00
Wavefunctions & electron density are generated automatically.
total charge
10.0000000000000(electrons)
fermi level -0.156250000000000(hartree)atomic charge(electrons)
1
3.96711149641430
2
6.03288852488762
k,spin,
1
1
1
1
1
1
1
1
1
1
1
1
1
1

band ,
1
2
3
4
5
6
7

eigen value,-0.1146306028E+01-0.5206152117E+00-0.4881698896E+00-0.4881698895E+00-0.3228998209E+00-0.3496017757E-01-0.3496017755E-01

residual norm ,0.8007059172E-05
0.5328887496E-05
0.1228950837E-06
0.2098904900E-06
0.4285043306E-05
0.2643743139E-04
0.2643696682E-04

10occupation
0.1000000000E+01
0.1000000000E+01
0.1000000000E+01
0.1000000000E+01
0.1000000000E+01
0.2763756664E-17
0.2763756649E-17

11
80.1501474936E-01
0.6926060496E-02
0.1610230826E-24
Scf= 55dp=
0.8743034034D-08
one electron energy-5.93232167908051(hartree)exchange correlation energy -5.93120017284993(hartree)hartreeenergy-29.8579787407655(hartree)ewald
energy12.5967331459555(hartree)field correction
0.000000000000000E+000 (hartree)
B_con correction
0.000000000000000E+000 (hartree)
helmholtz free energy-29.1247674467404(hartree)Totalenergy-29.1247674467404(hartree)atomic force
1
0.2581769902112212D+00
0.2581769902114270D+00
0.2581769900183384D+00
2 -0.2581769902112212D+00 -0.2581769902114270D+00 -0.2581769900183384D+00
mdstep=1
atom coordinate (input)
0.5500000000000000D+00
0.5500000000000000D+00
0.5500000000000000D+00
6-0.5500000000000000D+00 -0.5500000000000000D+00 -0.5500000000000000D+00
8atom coordinate (output)
0.5500000000000000D+00
0.5500000000000000D+00
0.5500000000000000D+00
6-0.5500000000000000D+00 -0.5500000000000000D+00 -0.5500000000000000D+00
8
real time450.494811928016(sec)com. time 0.126425336726243(sec)com. vol. 0.000000000000000E+000 (byte)

0
0

0
0

0
0

22056.00
29408.00

00
00

1
1

1a
2a

0
0

0
0

0
0

22056.00
29408.00

00
00

1
1

1a
2a

mdresult.dat の total energy と atom force が，上の内容と 10 桁程度一致していれば，計算は成功である．コンパイラによっても計算順序が異なるため，10 桁以降の部分については，上の結果と異なることがある．

5

計算を並列化する

2.5 節で述べたように，並列化にはプロセス並列化とスレッド並列化がある．本章では，両方の並列化方法について紹介する．なお，本章の例題の mdresult.dat は，実行時間を除いて 4 章と同じである．ただし，並列化による演算順序の変更により，演算結果が 4 章の結果とわずかに異なることがある．

5.1

4×1 でスレッド並列化する

まず初めに，設定が簡単なスレッド並列化を行う．read.inp, atom.xyz ファイルは 4 章と同じであり，実行スクリプト job.sh のみ異なる．下記に 4 章からの変更部分を青字で示す．
job.sh ファイルの作成例
1: #!/bin/csh
2: #$ -cwd -pe mpi 4
3: setenv OMP NUM THREADS 4
4: mv $USER.$PE.$JOB ID.1 hosts
5: set TEST = ""6: foreach aaaa (‘cat hosts‘)
7: @ I++
8: set TEST = ($TEST ‘echo $aaaa‘)
9: end10: rm hosts
11: set i=0
12: while ( $i < $NSLOTS )13: @ i = $i + $OMP NUM THREADS
14: echo $TEST[$i] >> hosts
15: end
16: @ i = $NSLOTS / $OMP NUM THREADS
17:
18: mpirun -np $i -machinefile hosts ./kukan8

本書で使用を想定している PC クラスタでは，ジョブの投入方法も同じである．ただし，並列計算は多くの計算資源を使用するため，計算機システムによっては，並列計算のジョブ投入方法が11

異なる場合がある．使用する計算機システムの説明書を参照のこと．計算結果は，mdresult.dat
で確認できる．出力される mdresult.dat の total energy と atom force が，4 章のものと 10 桁程
度一致していればよい．また，mdredult.dat の下から 3 行目の real time が，4 章の結果と比べて減少していることを確認する．

5.2

1×8 でプロセス並列化する

8 プロセスで並列化する方法を説明する．atom.xyz ファイルは 4 章から変更はない．まず，read.inp の並列化制御行を下記のように書き換える．ここでは，計算領域を x, y, z 方向に 2 分割し，各サブ領域に 1 プロセスを割り当てる．このとき，各方向のグリッド数が，サブ領域数の倍数になっていなければならない (2×nxmax が nprocx で割り切れねばならない．) ことに注意する．下記に，4 章からの変更点を青字で示す．
read.inp ファイルの作成例 (初めの 69 行のみ)
1: ![ndisp]
2: 66
3: ![nprocx],[nprocy],[nprocz],[nprock]
4: 2 2 2 1
5: ![xmax],[ymax],[zmax]
6: 10.0d0 10.0d0 10.0d0
7: ![nxmax],[nymax],[nzmax]
8: 30 30 30
9: ! [npxmax],[npymax],[npzmax]
10: -8 8 26
11: ![nso],[socang(1)],[socang(2)]12: 0 0.0d0 0.0d0
13: ![nsym]
14: 0
15: ![neigmx]
16: 8
17: ![natom]
18: 2
19: ![num atcell]
20: -1
21: ![num ppcell]
22: -1
23: ![num ppcell d]
24: -1
25: ![gmax]
26: 20.0d0
27: ![nperi],[npopshow]
28: 0 0
29: ![numkx],[numky],[numkz],[ksym],[kband]
30: 1 1 1 0 0
31: !k point (write all k points)
32: ![X], [Y], [Z], [weight]
33: 0.00d0 0.00d0 0.00d0 1
34: ![cexco]
35: VWN
36: ![nspv]
37: 1
38: ![epsvh],[epssd],[ratio diis],[eps scf]
39: 1.0D-8 1.0D-8 0.3d0 1.0D-8
40: ![ncgmin],[ncgmax],[ncgres]41: 1 800 801
42: ![nprecon cg],[nprecon diis]
43: 1 044: ![nsdmax],[nkscg]
45: 4 1
46: ![ndiismax],[ncgscf],[nretcg]
47: 4 200 2000
48: ![nrrz],[nchange]
49: 1 1
50: ![eta],[looplimit],[nbrydn],[etamag(1)],[etamag(2)]51: 0.20D0 200 10 0.5d0 2.0d0
52: ![tmstep],[nmd start],[nmd end],[ngdiis],[sconst]
53: 0.0D0 0 1 1 0.0d0
54: ![biasx],[biasy],[biasz]
55: 0.0D0 0.0D0 0.0D0
56: ![tf],[tfmin],[tfmax]
57: 0.003D0 -10.0D0 30.0D0
58: ![nmdxyz]
59: 0
60: ![chrgd],[npolcon],[polconocc]

12

61:
62:
63:
64:
65:
66:
67:
68:
69:

0.0D0 0 2.7d0![endjel],[chrjel]
0.0d0 0.0d0![fcut]1.0D-4![npre],[nevhist]
1 3![northo]1

また，下記に実行用のスクリプトファイル job.sh の 4 章からの変更部分を青字で示す．
job.sh ファイルの作成例
1: #!/bin/csh
2: #$ -cwd -pe mpi 8
3: setenv OMP NUM THREADS 1
4: mv $USER.$PE.$JOB ID.1 hosts
5: set TEST = ""6: foreach aaaa (‘cat hosts‘)
7: @ I++
8: set TEST = ($TEST ‘echo $aaaa‘)
9: end10: rm hosts
11: set i=0
12: while ( $i < $NSLOTS )13: @ i = $i + $OMP NUM THREADS
14: echo $TEST[$i] >> hosts
15: end
16: @ i = $NSLOTS / $OMP NUM THREADS
17:
18: mpirun -np $i -machinefile hosts ./kukan8

ジョブの投入方法，ならびに計算結果の確認方法は，5.1 節と同じである．

5.3

4×8 でハイブリッド並列化する

4×8 で並列化する方法を説明する．read.inp, atom.xyz ファイルは 5.1 節と同じであり，実行スクリプト job.sh のみ異なる．下記に 4 章からの変更部分を青字で示す．
job.sh ファイルの作成例
1: #!/bin/csh
2: #$ -cwd -pe mpi 32
3: setenv OMP NUM THREADS 4
4: mv $USER.$PE.$JOB ID.1 hosts
5: set TEST = ""6: foreach aaaa (‘cat hosts‘)
7: @ I++
8: set TEST = ($TEST ‘echo $aaaa‘)
9: end10: rm hosts
11: set i=0
12: while ( $i < $NSLOTS )13: @ i = $i + $OMP NUM THREADS
14: echo $TEST[$i] >> hosts
15: end
16: @ i = $NSLOTS / $OMP NUM THREADS
17:
18: mpirun -np $i -machinefile hosts ./kukan8

ジョブの投入方法，ならびに計算結果の確認方法は，5.1 節と同じである．

6

原子構造の最適化

6.1

パラメータファイルの作成

原子構造最適化を行うには，read.inp と atom.xyz を書き換える．以下に，read.inp の例を示
す．5.2 節からの変更点を青字にしている．
read.inp ファイルの作成例 (初めの 69 行のみ)

13

1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:

![ndisp]
66![nprocx],[nprocy],[nprocz],[nprock]2 2 2 1![xmax],[ymax],[zmax]
10.0d0 10.0d0 10.0d0![nxmax],[nymax],[nzmax]
30 30 30
! [npxmax],[npymax],[npzmax]-8 8 26![nso],[socang(1)],[socang(2)]0 0.0d0 0.0d0![nsym]0![neigmx]
8![natom]2![num atcell]-1![num ppcell]-1![num ppcell d]-1![gmax]20.0d0![nperi],[npopshow]
0 0![numkx],[numky],[numkz],[ksym],[kband]1 1 1 0 0
!k point (write all k points)![X], [Y], [Z], [weight]
0.00d0 0.00d0 0.00d0 1![cexco]
VWN![nspv]1![epsvh],[epssd],[ratio diis],[eps scf]
1.0D-8 1.0D-8 0.3d0 1.0D-8![ncgmin],[ncgmax],[ncgres]1 800 801![nprecon cg],[nprecon diis]
1 0![nsdmax],[nkscg]
4 1![ndiismax],[ncgscf],[nretcg]
4 200 2000![nrrz],[nchange]
1 1![eta],[looplimit],[nbrydn],[etamag(1)],[etamag(2)]0.20D0 200 10 0.5d0 2.0d0![tmstep],[nmd start],[nmd end],[ngdiis],[sconst]
80.0D0 0 20 5 0.0d0![biasx],[biasy],[biasz]
0.0D0 0.0D0 0.0D0![tf],[tfmin],[tfmax]
0.003D0 -10.0D0 30.0D0![nmdxyz]
0![chrgd],[npolcon],[polconocc]
0.0D0 0 2.7d0![endjel],[chrjel]
0.0d0 0.0d0![fcut]1.0D-4![npre],[nevhist]
1 3![northo]1

次に atom.xyz の例を示す．ここでは，原子構造を最適化するので，各方向の構造最適化のス
イッチを 1 にする．4 章からの変更点を青字にしている．
atom.xyz の作成例 (初めの 3 行のみ)
1:! [x], [y], [z], [atom number], switch [x], [y], [z], [weight], switches [soc], [pp], [na]2: 0.5500000000000000D+00 0.5500000000000000D+00 0.5500000000000000D+00 6 1 1 1 22056.00 00 1 1a
3:-0.5500000000000000D+00 -0.5500000000000000D+00 -0.5500000000000000D+00 8 1 1 1 29408.00 00 1 2a

14

6.2

プログラムの実行

パラメータファイルの変更が終了したら，プログラムを実行する．プログラムの実行は，5.3
節で用いた実行スクリプト job.sh を用いて，ジョブを投入する．
プログラムが終了すると，mdresult.dat に各原子構造最適化ステップごとの電子状態の情報が出力される．mdresult.dat の最後の部分 mdstep=20 のところで，total energy が約-29.170556(a.u.)，
atom coordinate (output) から計算した原子間距離 2 つの原子の距離が 2.14(a.u.) 程度になっていれば成功である．

7

スピン分極した系の電子状態を計算する

ほとんどの原子やイオン，O2 分子，S2 分子などはスピン分極しており，このようなものを扱
うにはスピンを考慮した計算を行う必要がある．ここでは，O2 分子の電子状態計算を例に説明
する．

7.1

ディレクトリの準備

4 章と同様にディレクトリの作成とファイルのコピーを行う．擬ポテンシャルを格納するディレクトリ pspaw には，酸素用の擬ポテンシャル paw.008 をコピーする．

7.2

パラメータファイルの作成

O2 分子を計算するには read.inp と atom.xyz を書き換える．以下に，read.inp の例を示す．5.2
章からの変更点を青字にしている．
read.inp ファイルの作成例 (初めの 69 行のみ)
1: ![ndisp]
2: 66
3: ![nprocx],[nprocy],[nprocz],[nprock]
4: 2 2 2 1
5: ![xmax],[ymax],[zmax]
6: 10.0d0 10.0d0 10.0d0
7: ![nxmax],[nymax],[nzmax]
8: 30 30 30
9: ! [npxmax],[npymax],[npzmax]
10: -8 8 26
11: ![nso],[socang(1)],[socang(2)]12: 0 0.0d0 0.0d0
13: ![nsym]
14: 0
15: ![neigmx]
16: 8
17: ![natom]
18: 2
19: ![num atcell]
20: -1
21: ![num ppcell]
22: -1
23: ![num ppcell d]
24: -1
25: ![gmax]
26: 20.0d0
27: ![nperi],[npopshow]
28: 0 0
29: ![numkx],[numky],[numkz],[ksym],[kband]
30: 1 1 1 0 0
31: !k point (write all k points)
32: ![X], [Y], [Z], [weight]

15

33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:

0.00d0 0.00d0 0.00d0 1![cexco]
VWN![nspv]2![epsvh],[epssd],[ratio diis],[eps scf]
1.0D-8 1.0D-8 0.3d0 1.0D-8![ncgmin],[ncgmax],[ncgres]1 800 801![nprecon cg],[nprecon diis]
1 0![nsdmax],[nkscg]
4 1![ndiismax],[ncgscf],[nretcg]
4 200 2000![nrrz],[nchange]
1 1![eta],[looplimit],[nbrydn],[etamag(1)],[etamag(2)]0.20D0 200 10 0.5d0 2.0d0![tmstep],[nmd start],[nmd end],[ngdiis],[sconst]
0.0D0 0 1 1 0.0d0![biasx],[biasy],[biasz]
0.0D0 0.0D0 0.0D0![tf],[tfmin],[tfmax]
0.003D0 -10.0D0 30.0D0![nmdxyz]
0![chrgd],[npolcon],[polconocc]
0.0D0 0 2.7d0![endjel],[chrjel]
0.0d0 0.0d0![fcut]1.0D-4![npre],[nevhist]
1 3![northo]1

次に atom.xyz の例を示す．
atom.xyz の作成例 (初めの 3 行のみ)
1:! [x], [y], [z], [atom number], switch [x], [y], [z], [weight], switches [soc], [pp], [na]2: 0.1100000000000000D+01 0.0000000000000000D+00 0.0000000000000000D+00 8 1 1 1 29408.00 00 1 1a
3:-0.1100000000000000D+01 0.0000000000000000D+00 0.0000000000000000D+00 8 1 1 1 29408.00 00 1 2a

7.3

プログラムの実行

パラメータファイルの変更が終了したら，プログラムを実行する．プログラムの実行は，5.3
節で用いた実行スクリプト job.sh を用いて，ジョブを投入する．

7.4

結果の確認

計算が終了すると mdresult.dat ファイルが出力される．以下に，mdresult.dat の出力例を示す．
computed by kukan8e 03/29/2017-01========== data from read.inp ==========
ndisp
:
66
nprocx:2
nprocy:2nprocz:2
xmax:
10.0000000000000
ymax:
10.0000000000000
zmax:
10.0000000000000
nxmax
:
30
nymax
:
30
nzmax
:
30
npxmax:-8
npymax:8
npzmax
:
26
nso
:
0

16

socang:
nsym:
neigmx:natom:
num_atcell :num_ppcell :num_ppcell_d:
gmaxps:nperi:
npopshow
:
numkx
:
numky
:
numkz
:
ksym
:
kband
:
cexco
:
nspv
:
epsvh
:
epssd
:
ratio_diis :eps_scf:
ncgmin:ncgmax:ncgres:nprecon_cg :nprecon_diis:
nsdmax:nkscg:
ndiismax:
ncgscf:nretcg:nrrz:nchange:eta:
etamag:looplimit:
nbrydn:tmstep:nmd_start:
nmd_end:ngdiis:sconst
:
biasx
:
biasy
:
biasz
:
tf
:
tfmin
:
tfmax
:
nmdxyz:chrgd:
npolcon:polconocc:
endjel:chrjel
:
fcut
:
npre
:
nevhist:northo:eps:
eps_eig_diis:
alambda_diis:
alambda_min :alambda_max :nradmx:nrprjmx:
nprjmx:lsphel:nlmax:
lrhomx:nfiltyp:
psctoff:
nqmx:

0.000E+00 0.000E+00
0
8
2
2
2
2
20.0000000000000
0
0
1
1
1
0
0
vwn
2
1.000000000000000E-008
1.000000000000000E-008
0.300000000000000
1.000000000000000E-008
1
800
801
1
0
4
1
4
200
2000
1
1
0.200000000000000
0.5000E+00 0.2000E+01200
100.000000000000000E+000
01
10.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
3.000000000000000E-003-10.0000000000000
30.0000000000000
00.000000000000000E+000
02.70000000000000
0.000000000000000E+000
0.000000000000000E+000
1.000000000000000E-004
1
3
1
1.000000000000000E-016
1.000000000000000E-014
0.500000000000000
0.100000000000000
1.00000000000000
1502
6
18
8
4
25
1
1.05000000000000
400

0.000E+00

17psftrad:10.0000000000000
psctrat:1.00000000000000
psext:
2.00000000000000
filpp: 2.000000000000000E-002
rctpcc:2.00000000000000
nmesh
:
2
npmesh:2veta: 0.200000000000000
new_pwx:0
new_pwy:0
new_pwz:0
new_rsx:0
new_rsy:0
new_rsz:0
nint1dmax
:
1000
ndec
:
12
nf
:
4
nfdg
:
4
nfh
:
4
nrasy
:
0
nrasz
:
0
nrepx
:
1
zs_pre: 0.500000000000000
pol_pre:1.50000000000000==========================================
Filtering of pseudopotentials is implemented by Fermi distribution.
npxmax, npymax, and npzmax are updated below.
========== data computed by listvecdim ==========
npxmax:8
npymax:8
npzmax:8====================================================
cell size [lx,ly,lz (a.u.)]
0.200000000D+02
0.200000000D+02
0.200000000D+02
coarse grid spacing [hx,hy,hz (a.u.)]
0.333333333D+00
0.333333333D+00
0.333333333D+00
dense grid spacing [hx,hy,hz (a.u.)]
0.166666667D+00
0.166666667D+00
0.166666667D+00
cutoff energy of wave functions [gmax_coarse,gmax_dense (Ry)]
0.888264396D+02
0.355305758D+03
number of electron
chargespin
0.120000000D+02
0.000000000D+00
0.200000000D+01
electric field [ex,ey,ez (a.u./a.u.)]
0.000000000D+00
0.000000000D+00
0.000000000D+00
time step (a.u.) (fs)0.000000000D+00
0.000000000D+00
thermo (k)kt (a.u.)
0.947268709D+03
0.300000000D-02
cufoff of force acting on atoms (a.u.)
0.100000000D-03
charge mixing ratio of scf
0.200000000D+00
mag. mix. rat., weight of mag. in F_broyd
0.500000000D+00
0.200000000D+01
nrasy,nrasz
00
Wavefunctions & electron density are generated automatically.
total charge
12.0000000000000(electrons)
fermi level -0.200971118756570(hartree)spin polarization -1.99998396516831(\mu_B)
atomic charge(electrons)
1
6.00000133848120
25.99999819350704
k,spin,
1
1
1
1
1
1
1
1
1
1
1
1
1
1

band ,
1
2
3
4
5
6
7

eigen value,-0.1208809610E+01-0.6768339336E+00-0.4707114513E+00-0.4570551337E+00-0.4570551337E+00-0.1636900129E+00-0.1636900128E+00

residual norm ,0.1983694557E-06
0.7265701623E-07
0.8016383386E-07
0.7569689472E-06
0.7644773082E-06
0.8111198396E-05
0.8111483469E-05

18

occupation
0.1000000000E+01
0.1000000000E+01
0.1000000000E+01
0.1000000000E+01
0.1000000000E+01
0.4008717049E-05
0.4008717023E-05

11
80.4081082372E-01
0.2043701728E-01
0.9965035554E-35
12
1-0.1254553948E+01
0.1165788014E-06
0.1000000000E+01
12
2-0.7390928147E+00
0.8173880198E-07
0.1000000000E+01
12
3-0.5193707662E+00
0.2067187376E-05
0.1000000000E+01
12
4-0.5193707662E+00
0.2070250955E-05
0.1000000000E+01
12
5-0.5132248426E+00
0.6074156678E-07
0.1000000000E+01
12
6-0.2382522248E+00
0.3582800592E-05
0.9999959913E+00
12
7-0.2382522247E+00
0.3585528225E-05
0.9999959913E+00
12
80.3827315151E-01
0.1039612093E-01
0.2321905455E-34
Scf= 30dp=
0.9568422665D-08
one electron energy-7.29258225155358(hartree)exchange correlation energy -6.53042442690354(hartree)hartreeenergy-43.0468963687653(hartree)ewald
energy16.3636363636364(hartree)field correction
0.000000000000000E+000 (hartree)
B_con correction
0.000000000000000E+000 (hartree)
helmholtz free energy-40.5062673294884(hartree)Totalenergy-40.5062666835861(hartree)atomic force
1
0.7867502558242023D-01
0.9811362833289650D-09
0.4125588759507082D-12
2 -0.7867502558242023D-01 -0.9811362833289650D-09 -0.4125588759507082D-12
mdstep=1
atom coordinate (input)
0.1100000000000000D+01
0.0000000000000000D+00
0.0000000000000000D+00
8-0.1100000000000000D+01
0.0000000000000000D+00
0.0000000000000000D+00
8atom coordinate (output)
0.1100000000000000D+01
0.0000000000000000D+00
0.0000000000000000D+00
8-0.1100000000000000D+01
0.0000000000000000D+00
0.0000000000000000D+00
8
real time
115.700261442922(sec)com. time
24.1898171515204(sec)com. vol.
13922409984.0000(byte)

0
0

0
0

0
0

29408.00
29408.00

00
00

1
1

1a
2a

0
0

0
0

0
0

29408.00
29408.00

00
00

1
1

1a
2a

mdresult.dat の total energy と atom force が，上の内容と 10 桁程度一致していれば，計算は成功である．
引き続き，原子構造最適化を行うには，同じく 6 章の要領で行う．

8

結晶や表面の電子状態計算

本書のこれまでの計算は，孤立したクラスターの電子状態を求める計算であった．本章では周期構造をもった結晶や表面の計算方法を説明する．

8.1

スーパーセル

周期性を持つ結晶や表面を計算する場合には，その 1 周期分を 1 つのユニットとしたスーパー
セルを用いる．例えば，図 2 のように規則的に原子の並んだ，一層の薄膜を考えよう．図において，白抜きの丸は原子である．薄膜は図中央の正方形を最小単位として周期構造をなしている．このような場合，すべての原子を計算しなくても，正方形内だけを計算すればよい．ここで注意しなくてはならないことは，計算対象となる原子の数は 5 個ではなく 2 個である．なぜなら，原子 A を中心に考えると B，C，D は原子 A のコピーである．このことは，スーパーセ
ルを 1 周期分上下左右に平行移動させれば，原子 A が B，C，D に重なることを考えれば簡単に理解できる．このように，スーパーセルを用いると計算しなくてはならない原子の数が，見かけより少ない場合がよくある．モデルを組むときに，スーパーセルを 1 周期分平行移動させると原子が重なってしまうことがないように注意する必要がある．

19

B
C

E

D
A

図 2: スーパーセル

8.2

リチウム結晶の電子状態計算を行う

結晶は，x，y ，z のすべての方向に周期性を持っている．したがって，read.inp の境界条件を決定するパラメータ nperi は 3 である．
4 章と同様にディレクトリの作成とファイルのコピーを行う．擬ポテンシャルを格納するディレクトリ pspaw には，リチウム用の擬ポテンシャル paw.003 をコピーする．
まず，read.inp を以下のように書き換える．リチウムの格子定数は 6.64(a.u.) なので，計算領
域 2×{x,y,z}max が 6.64(a.u.) になるように {x,y,z}max を 3.32 とする．以下に，read.inp の例を示す．4 章からの変更点を青字にしている．
read.inp ファイルの作成例 (初めの 69 行のみ)
1: ![ndisp]
2: 66
3: ![nprocx],[nprocy],[nprocz],[nprock]
4: 1 1 1 1
5: ![xmax],[ymax],[zmax]
6: 3.32d0 3.32d0 3.32d0
7: ![nxmax],[nymax],[nzmax]
8: 10 10 10
9: ! [npxmax],[npymax],[npzmax]
10: -8 8 26
11: ![nso],[socang(1)],[socang(2)]12: 0 0.0d0 0.0d0
13: ![nsym]
14: 0
15: ![neigmx]
16: 8
17: ![natom]
18: 2
19: ![num atcell]
20: -1
21: ![num ppcell]
22: -1
23: ![num ppcell d]
24: -1

20

25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:

![gmax]20.0d0![nperi],[npopshow]
3 0![numkx],[numky],[numkz],[ksym],[kband]1 1 1 0 0
!k point (write all k points)![X], [Y], [Z], [weight]
0.00d0 0.00d0 0.00d0 1![cexco]
VWN![nspv]1![epsvh],[epssd],[ratio diis],[eps scf]
1.0D-8 1.0D-8 0.3d0 1.0D-8![ncgmin],[ncgmax],[ncgres]1 800 801![nprecon cg],[nprecon diis]
1 0![nsdmax],[nkscg]
4 1![ndiismax],[ncgscf],[nretcg]
4 200 2000![nrrz],[nchange]
1 1![eta],[looplimit],[nbrydn],[etamag(1)],[etamag(2)]0.20D0 200 10 0.5d0 2.0d0![tmstep],[nmd start],[nmd end],[ngdiis],[sconst]
0.0D0 0 1 1 0.0d0![biasx],[biasy],[biasz]
0.0D0 0.0D0 0.0D0![tf],[tfmin],[tfmax]
0.003D0 -10.0D0 30.0D0![nmdxyz]
0![chrgd],[npolcon],[polconocc]
0.0D0 0 2.7d0![endjel],[chrjel]
0.0d0 0.0d0![fcut]1.0D-4![npre],[nevhist]
1 3![northo]
1

次に，atom.xyz を以下のように書き換える．8.1 節で述べたように，リチウムは体心立方晶なので最小単位セルに必要な原子数は 2 個であり，座標ファイルは次の例のようになる．
atom.xyz の作成例 (初めの 3 行のみ)
1:! [x], [y], [z], [atom number], switch [x], [y], [z], [weight], switches [soc], [pp], [na]2:0.3320000000000000D+01 0.3320000000000000D+01 0.3320000000000000D+01 3 0 0 0 12866.00 00 1 1a
3:0.0000000000000000D+00 0.0000000000000000D+00 0.0000000000000000D+00 3 0 0 0 12866.00 00 1 2a

ここで，1 行目の原子の座標の正負は任意だが，プログラムの中で自動的に上のように変更されるため，mdresult.dat などの出力ファイルには全て変更された後の値が出力される 2 ．5.1 章で作成した実行スクリプト job.sh を使い，プログラムを実行する．以下に，電子状態計算の実行結果 mdresult.dat を示す．total energy が次の内容と 10 桁程度一致していれば，計算は成功である．
mdresult.dat の出力例
computed by kukan8e 03/29/2017-01========== data from read.inp ==========
ndisp
:
66
nprocx:1
nprocy:1nprocz:1
xmax:
3.32000000000000
ymax:
3.32000000000000
zmax:
3.32000000000000
nxmax
:
10
nymax
:
10
nzmax
:
10
npxmax:-8

周期境界条件の場合は，プログラム中で-xmax+dx/2 ≤ x ≤ +xmax+dx/2 の範囲内に原子が納まるように変更
される．
2

21

npymax:
npzmax:nso:
socang:nsym:
neigmx:natom:
num_atcell :num_ppcell :num_ppcell_d:
gmaxps:nperi:
npopshow
:
numkx
:
numky
:
numkz
:
ksym
:
kband
:
cexco
:
nspv
:
epsvh
:
epssd
:
ratio_diis :eps_scf:
ncgmin:ncgmax:ncgres:nprecon_cg :nprecon_diis:
nsdmax:nkscg:
ndiismax:
ncgscf:nretcg:nrrz:nchange:eta:
etamag:looplimit:
nbrydn:tmstep:nmd_start:
nmd_end:ngdiis:sconst
:
biasx
:
biasy
:
biasz
:
tf
:
tfmin
:
tfmax
:
nmdxyz:chrgd:
npolcon:polconocc:
endjel:chrjel
:
fcut
:
npre
:
nevhist:northo:eps:
eps_eig_diis:
alambda_diis:
alambda_min :alambda_max :nradmx:nrprjmx:
nprjmx:lsphel:nlmax:
lrhomx:826
00.000E+00 0.000E+00
0
8
2
2
2
2
20.0000000000000
3
0
1
1
1
0
0
vwn
1
1.000000000000000E-008
1.000000000000000E-008
0.300000000000000
1.000000000000000E-008
1
800
801
1
0
4
1
4
200
2000
1
1
0.200000000000000
0.5000E+00 0.2000E+01200
100.000000000000000E+000
01
10.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
3.000000000000000E-003-10.0000000000000
30.0000000000000
00.000000000000000E+000
02.70000000000000
0.000000000000000E+000
0.000000000000000E+000
1.000000000000000E-004
1
3
1
1.000000000000000E-016
1.000000000000000E-014
0.500000000000000
0.100000000000000
1.00000000000000
1502
6
18
8
4
25

0.000E+00

22nfiltyp:1
psctoff:
1.05000000000000
nqmx
:
400
psftrad:10.0000000000000
psctrat:1.00000000000000
psext:
2.00000000000000
filpp: 2.000000000000000E-002
rctpcc:2.00000000000000
nmesh
:
2
npmesh:2veta: 0.200000000000000
new_pwx:
11
new_pwy:11
new_pwz:
11
new_rsx:
11
new_rsy:11
new_rsz:
11
nint1dmax
:
1000
ndec
:
12
nf
:
4
nfdg
:
4
nfh
:
4
nrasy
:
0
nrasz
:
0
nrepx
:
1
zs_pre: 0.500000000000000
pol_pre:1.50000000000000==========================================
Filtering of pseudopotentials is implemented by Fermi distribution.
npxmax, npymax, and npzmax are updated below.
========== data computed by listvecdim ==========
npxmax:13
npymax:
13
npzmax:13====================================================
cell size [lx,ly,lz (a.u.)]
0.664000000D+01
0.664000000D+01
0.664000000D+01
coarse grid spacing [hx,hy,hz (a.u.)]
0.332000000D+00
0.332000000D+00
0.332000000D+00
dense grid spacing [hx,hy,hz (a.u.)]
0.166000000D+00
0.166000000D+00
0.166000000D+00
cutoff energy of wave functions [gmax_coarse,gmax_dense (Ry)]0.895413376D+02
0.358165351D+03
sample k point [skx,sky,skz (2pi/l)] weight
0.000000000D+00
0.000000000D+00
0.000000000D+00
1number of electron
chargespin
0.200000000D+01
0.000000000D+00
0.100000000D+01
electric field [ex,ey,ez (a.u./a.u.)]
0.000000000D+00
0.000000000D+00
0.000000000D+00
time step (a.u.) (fs)0.000000000D+00
0.000000000D+00
thermo (k)kt (a.u.)
0.947268709D+03
0.300000000D-02
cufoff of force acting on atoms (a.u.)
0.100000000D-03
charge mixing ratio of scf
0.200000000D+00
nrasy,nrasz
00
Wavefunctions & electron density are generated automatically.
total charge
2.00000000000000(electrons)
fermi level 0.000000000000000E+000 (hartree)k,spin,
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

band ,
1
2
3
4
5
6
7
8

eigen value,-0.1167088446E+00
0.2003272717E+00
0.2003272717E+00
0.2003272717E+00
0.4749743482E+00
0.5117631824E+00
0.5117631824E+00
0.5117631824E+00

residual norm ,0.1730416405E-06
0.9706159597E-07
0.1144501253E-06
0.2579767073E-06
0.6641748123E-05
0.7903519841E-07
0.1145198747E-06
0.1304499639E-06

23

occupation
0.1000000000E+01
0.9992107771E-29
0.9992107757E-29
0.9992107757E-29
0.0000000000E+00
0.0000000000E+00
0.0000000000E+00
0.0000000000E+00

Scf= 21dp=
0.9339998342D-08
one electron energy-0.233417689152553(hartree)exchange correlation energy -2.94987053657499(hartree)hartreeenergy-2.729244276725520E-003 (hartree)ewald
energy-0.548077326733234(hartree)field correction
0.000000000000000E+000 (hartree)
B_con correction
0.000000000000000E+000 (hartree)energy offset of ps. pot. -8.985242223158289E-004 (hartree)energy jellium/ions
0.000000000000000E+000 (hartree)
helmholtz free energy-3.73499332095982(hartree)Totalenergy-3.73499332095982(hartree)atomic force
10.1799008914857649D-17 -0.6247067097839938D-17 -0.2711183739934840D-16
2 -0.1799008914857646D-17
0.6247067097839936D-17
0.2711183739934840D-16
mdstep=1
atom coordinate (input)
0.3320000000000000D+01
0.3320000000000000D+01
0.3320000000000000D+01
30.0000000000000000D+00
0.0000000000000000D+00
0.0000000000000000D+00
3
atom coordinate (output)
0.3320000000000000D+01
0.3320000000000000D+01
0.3320000000000000D+01
30.0000000000000000D+00
0.0000000000000000D+00
0.0000000000000000D+00
3
real time
9.41127430897905(sec)com. time 0.790023893845500(sec)com. vol. 0.000000000000000E+000 (byte)

8.3

0
0

0
0

0
0

12866.00
12866.00

00
00

1
1

1a
2a

0
0

0
0

0
0

12866.00
12866.00

00
00

1
1

1a
2a

リチウム表面の電子状態計算を行う

次に，表面の電子状態計算を行う．第一原理電子状態計算手法で，表面のように面垂直方向の一方は真空，反対方向には結晶が無限に続く物体の計算方法は，現在のところ実用化されていない．したがって，表面の計算を行うには，原子層を何層か積んだ薄膜を表面に見立てて計算している．ここでは，リチウム (001)4 原子層薄膜の電子状態計算を行う．薄膜は，x，y 方向に周期性を持っている．したがって，read.inp の境界条件を決定するパラメータ nperi は 2 である．ディレクトリの作成とファイルのコピーは，8.2 節と全く同じである．作業後，8.2 節と同じく，read.inp と atom.xyz を以下のように編集する．read.inp は，8.2 節からの変更点を青字で示す．リチウムの格子定数は 6.64(a.u.) なので，計算領域 2×{x,y,z }max が 6.64(a.u.) になるように {x,y,z }max を 3.32 とする．また，周期性のない z 方向は，境界面を原子から離しておかなければならないので，長めに設定しておく．また，4 原子層なので，最小単位セルに必
要な原子数は 4 個である．
read.inp ファイルの作成例 (初めの 69 行のみ)
1: ![ndisp]
2: 66
3: ![nprocx],[nprocy],[nprocz],[nprock]
4: 1 1 1 1
5: ![xmax],[ymax],[zmax]
6: 3.32d0 3.32d0 13.28d0
7: ![nxmax],[nymax],[nzmax]
8: 10 10 40
9: ! [npxmax],[npymax],[npzmax]
10: -8 8 26
11: ![nso],[socang(1)],[socang(2)]12: 0 0.0d0 0.0d0
13: ![nsym]
14: 0
15: ![neigmx]
16: 16
17: ![natom]
18: 4
19: ![num atcell]
20: -1
21: ![num ppcell]
22: -1
23: ![num ppcell d]
24: -1
25: ![gmax]
26: 20.0d0

24

27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:

![nperi],[npopshow]
2 0![numkx],[numky],[numkz],[ksym],[kband]1 1 1 0 0
!k point (write all k points)![X], [Y], [Z], [weight]
0.00d0 0.00d0 0.00d0 1![cexco]
VWN![nspv]1![epsvh],[epssd],[ratio diis],[eps scf]
1.0D-8 1.0D-8 0.3d0 1.0D-8![ncgmin],[ncgmax],[ncgres]1 800 801![nprecon cg],[nprecon diis]
1 0![nsdmax],[nkscg]
4 1![ndiismax],[ncgscf],[nretcg]
4 200 2000![nrrz],[nchange]
1 1![eta],[looplimit],[nbrydn],[etamag(1)],[etamag(2)]0.20D0 200 10 0.5d0 2.0d0![tmstep],[nmd start],[nmd end],[ngdiis],[sconst]
0.0D0 0 1 1 0.0d0![biasx],[biasy],[biasz]
0.0D0 0.0D0 0.0D0![tf],[tfmin],[tfmax]
0.003D0 -10.0D0 30.0D0![nmdxyz]
0![chrgd],[npolcon],[polconocc]
0.0D0 0 2.7d0![endjel],[chrjel]
0.0d0 0.0d0![fcut]1.0D-4![npre],[nevhist]
1 3![northo]1

atom.xyz の作成例 (初めの 5 行のみ)
1:! [x], [y], [z], [atom number], switch [x], [y], [z], [weight], switches [soc], [pp], [na]2: 0.3320000000000000D+01 0.3320000000000000D+01 -0.4980000000000000D+01 3 0 0 0 12866.00 00 1 1a
3: 0.0000000000000000D+00 0.0000000000000000D+00 -0.1660000000000000D+01 3 0 0 0 12866.00 00 1 2a
4: 0.3320000000000000D+01 0.3320000000000000D+01 0.1660000000000000D+01 3 0 0 0 12866.00 00 1 3a
5: 0.0000000000000000D+00 0.0000000000000000D+00 0.4980000000000000D+01 3 0 0 0 12866.00 00 1 4a5.1 章で作成した実行スクリプト job.sh を使い，プログラムを実行する．以下に，電子状態計
算の実行結果 mdresult.dat を示す．total energy が，次の内容と 10 桁程度一致していれば，計算は成功である．
mdresult.dat の出力例
computed by kukan8e 03/29/2017-01========== data from read.inp ==========
ndisp
:
66
nprocx:1
nprocy:1nprocz:1
xmax:
3.32000000000000
ymax:
3.32000000000000
zmax:
13.2800000000000
nxmax
:
10
nymax
:
10
nzmax
:
40
npxmax:-8
npymax:8
npzmax:26
nso:
0
socang:0.000E+00 0.000E+00 0.000E+00
nsym:
0
neigmx
:
16
natom
:
4
num_atcell :425num_ppcell :num_ppcell_d:
gmaxps:nperi:
npopshow
:
numkx
:
numky
:
numkz
:
ksym
:
kband
:
cexco
:
nspv
:
epsvh
:
epssd
:
ratio_diis :eps_scf:
ncgmin:ncgmax:ncgres:nprecon_cg :nprecon_diis:
nsdmax:nkscg:
ndiismax:
ncgscf:nretcg:nrrz:nchange:eta:
etamag:looplimit:
nbrydn:tmstep:nmd_start:
nmd_end:ngdiis:sconst
:
biasx
:
biasy
:
biasz
:
tf
:
tfmin
:
tfmax
:
nmdxyz:chrgd:
npolcon:polconocc:
endjel:chrjel
:
fcut
:
npre
:
nevhist:northo:eps:
eps_eig_diis:
alambda_diis:
alambda_min :alambda_max :nradmx:nrprjmx:
nprjmx:lsphel:nlmax:
lrhomx:nfiltyp:
psctoff:
nqmx:
psftrad:psctrat
:
psext
:
filpp
:
rctpcc:4
420.0000000000000
2
0
1
1
1
0
0
vwn
1
1.000000000000000E-008
1.000000000000000E-008
0.300000000000000
1.000000000000000E-008
1
800
801
1
0
4
1
4
200
2000
1
1
0.200000000000000
0.5000E+00 0.2000E+01200
100.000000000000000E+000
01
10.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
3.000000000000000E-003-10.0000000000000
30.0000000000000
00.000000000000000E+000
02.70000000000000
0.000000000000000E+000
0.000000000000000E+000
1.000000000000000E-004
1
3
1
1.000000000000000E-016
1.000000000000000E-014
0.500000000000000
0.100000000000000
1.00000000000000
1502
6
18
8
4
25
1
1.05000000000000
400
10.0000000000000
1.00000000000000
2.00000000000000
2.000000000000000E-002
2.00000000000000

26

nmesh:2npmesh:2veta: 0.200000000000000
new_pwx:
11
new_pwy:11
new_pwz:0
new_rsx:
11
new_rsy:11
new_rsz:0
nint1dmax
:
1000
ndec
:
12
nf
:
4
nfdg
:
4
nfh
:
4
nrasy
:
0
nrasz
:
0
nrepx
:
1
zs_pre: 0.500000000000000
pol_pre:1.50000000000000==========================================
Filtering of pseudopotentials is implemented by Fermi distribution.
npxmax, npymax, and npzmax are updated below.
========== data computed by listvecdim ==========
npxmax:13
npymax:
13
npzmax:13====================================================
cell size [lx,ly,lz (a.u.)]
0.664000000D+01
0.664000000D+01
0.265600000D+02
coarse grid spacing [hx,hy,hz (a.u.)]
0.332000000D+00
0.332000000D+00
0.332000000D+00
dense grid spacing [hx,hy,hz (a.u.)]
0.166000000D+00
0.166000000D+00
0.166000000D+00
cutoff energy of wave functions [gmax_coarse,gmax_dense (Ry)]0.895413376D+02
0.358165351D+03
sample k point [skx,sky,skz (2pi/l)] weight
0.000000000D+00
0.000000000D+00
0.000000000D+00
1number of electron
chargespin
0.400000000D+01
0.000000000D+00
0.100000000D+01
electric field [ex,ey,ez (a.u./a.u.)]
0.000000000D+00
0.000000000D+00
0.000000000D+00
time step (a.u.) (fs)0.000000000D+00
0.000000000D+00
thermo (k)kt (a.u.)
0.947268709D+03
0.300000000D-02
cufoff of force acting on atoms (a.u.)
0.100000000D-03
charge mixing ratio of scf
0.200000000D+00
nrasy,nrasz
00
Wavefunctions & electron density are generated automatically.
total charge
4.00000000000000(electrons)
fermi level -0.150688094499856(hartree)atomic charge(electrons)
1 0.966936137479371
2
1.03306424821140
3
1.03306323851557
4 0.966936375746787
k,spin,
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

band ,
1
2
3
4
5
6
7
8
9
10
11

eigen value,-0.2176215299E+00-0.1755650933E+00-0.1258110980E+00-0.7437630773E-01-0.1130978723E-01
0.8103562302E-01
0.9697083340E-01
0.9697083340E-01
0.1215480039E+00
0.1475229732E+00
0.1475229733E+00

residual norm ,0.1319010351E-06
0.3407204571E-06
0.1192467683E-05
0.5275046797E-05
0.2239496850E-04
0.3499588687E-03
0.3846284026E-06
0.4321257975E-06
0.6505736683E-03
0.4424579051E-05
0.4431719605E-05

27

occupation
0.9999999998E+00
0.9997496332E+00
0.2503670419E-03
0.8968862531E-11
0.6651558559E-20
0.2848105360E-33
0.1405070360E-35
0.1405070360E-35
0.3888554921E-39
0.6753459162E-43
0.6753459131E-4311
120.2221597812E+00
0.3024374259E-04
0.0000000000E+00
1
1
13
0.2221597813E+00
0.3029742479E-04
0.0000000000E+00
1
1
14
0.2396237192E+00
0.6079475767E-02
0.0000000000E+00
1
1
15
0.3130861322E+00
0.1915124964E-02
0.0000000000E+00
1
1
16
0.3130862711E+00
0.1913407736E-02
0.0000000000E+00
Scf= 35dp=
0.7924186543D-08
one electron energy-0.786348332780760(hartree)exchange correlation energy -5.90841330064037(hartree)hartreeenergy5.33757666721861(hartree)ewald
energy-5.96650333694422(hartree)field correction
0.000000000000000E+000 (hartree)
B_con correction
0.000000000000000E+000 (hartree)
helmholtz free energy-7.32371622143901(hartree)Totalenergy-7.32368830314674(hartree)atomic force
10.8454094657409363D-11
0.8454009593672319D-11 -0.3023884784036666D-02
2 -0.2125844702972279D-11 -0.2125766453801245D-11
0.1471990085794482D-02
30.1697195248769126D-11
0.1697289820595744D-11 -0.1471870953382727D-02
4 -0.8025445203206208D-11 -0.8025532960466819D-11
0.3023765651624911D-02
mdstep=1
atom coordinate (input)
0.3320000000000000D+01
0.3320000000000000D+01 -0.4980000000000000D+01
30.0000000000000000D+00
0.0000000000000000D+00 -0.1660000000000000D+01
3
0.3320000000000000D+01
0.3320000000000000D+01
0.1660000000000000D+01
30.0000000000000000D+00
0.0000000000000000D+00
0.4980000000000000D+01
3
atom coordinate (output)
0.3320000000000000D+01
0.3320000000000000D+01 -0.4980000000000000D+01
30.0000000000000000D+00
0.0000000000000000D+00 -0.1660000000000000D+01
3
0.3320000000000000D+01
0.3320000000000000D+01
0.1660000000000000D+01
30.0000000000000000D+00
0.0000000000000000D+00
0.4980000000000000D+01
3
real time
73.7137713429984(sec)com. time
13.5961111705983(sec)com. vol. 0.000000000000000E+000 (byte)

8.4

0
0
0
0

0
0
0
0

0
0
0
0

12866.00
12866.00
12866.00
12866.00

00
00
00
00

1
1
1
1

1a
2a
3a
4a

0
0
0
0

0
0
0
0

0
0
0
0

12866.00
12866.00
12866.00
12866.00

00
00
00
00

1
1
1
1

1a
2a
3a
4a

サンプリング k メッシュを細かくした計算を行う

結晶の中には，あらゆる波数ベクトル (周期) を持った電子が存在している．すなわち，波動関数は格子と完全に同じ周期のもの (波数ベクトルの成分が整数) は，ブリルアン領域内の波動関数のひとつにすぎない．8.2 章や 8.3 章で行った計算は，数ある波動関数の中から波数ベクトルの成分が整数のものだけを扱った．しかし，厳密な議論では，ブリルアン領域内のすべての波動関数を扱わなければならない．k 点を増やすということは，ブリルアン領域をメッシュで分割し，すべての波動関数を考慮に入れようというものであり，メッシュが細かくなれば細かくなるほど精度が増す．しかし，分割数が多くなればなるほど計算に要する時間も増すので，計算するモデルの大きさに等も考慮して，k 点の数を決めなければならない．
ここでは，リチウム結晶の電子状態を，第一ブリルアン領域を 4 × 4 × 4 に分割した 64k 点で計算する．それぞれの k 点のブリルアン領域での座標は，cx = L2πx ，cy = L2πy ，cz = L2πz とすると

28

( 1
)
1
1
(− 21 cx , − 21 cy , − 21 cz ) ,(− 21 cx , − 4 c1y , −)2 cz ,− 2 cz ,(− 21 cx , 0,)
− 2 cx , 14 cy , − 12 cz ,
)
( 1
1
1
(− 21 cx , − 21 cy , − 41 cz ) ,(− 21 cx , − 4 c1y , −)4 cz ,− 4 cz ,)(− 21 cx , 0,
− 2 cx , 14 cy , − 14 cz ,( 1)
1(− 21 cx , − 21 cy , 0) ,(− 21 cx , − 4 )cy , 0 ,0 ,)(− 21 cx , 0,
− 2 cx , 14 cy , 0 ,
( 1
)
1
1
(− 21 cx , − 21 cy , 41 cz ) ,(− 21 cx , − 41cy ,)4 cz ,
c ,
4 z
)
(− 21 cx , 0,1
1
− 2 cx , 4 cy , 4 cz ,

( 1
)
1
1
(− 41 cx , − 12 cy , − 12 cz ) ,(− 41 cx , − 4 c1y , −)2 cz ,− 2 cz ,(− 14 cx , 0,)
− 4 cx , 14 cy , − 12 cz ,
( 1
)
1
1
(− 41 cx , − 12 cy , − 14 cz ) ,(− 41 cx , − 4 c1y , −)4 cz ,− 4 cz ,(− 14 cx , 0,)
− 4 cx , 14 cy , − 14 cz ,( 1)
1(− 41 cx , − 12 cy , 0) ,(− 41 cx , − 4 )cy , 0 ,0 ,(− 14 cx , 0,)
− 4 cx , 14 cy , 0 ,
( 1
)
1
1
(− 14 cx , − 12 cy , 41 cz ) ,(− 14 cx , − 41cy ,)4 cz ,c ,4 z(− 41 cx , 0,)
1
1
− 4 cx , 4 cy , 4 cz ,()
1
1(0, − 21 cy , − 12 cz ) ,(0, − 4 c1y , −)2 cz ,− 2 cz ,(0, 0,)
0, 14 cy , − 12 cz ,()
1
1(0, − 21 cy , − 14 cz ) ,(0, − 4 c1y , −)4 cz ,− 4 cz ,(0, 0,)
0, 14 cy , − 14 cz ,()
1(0, − 21 cy , 0) ,0, − 4 cy , 0 ,(0, 0, 0) ,( 1)
(0, 4 c1y , 0 1, )(0, − 21 cy , 41 cz ) ,(0, − 41cy ,)4 cz ,c ,4 z
(0, 0,)
1
1
0, 4 cy , 4 cz ,

(1
)
1
1
( 21 cx , − 21 cy , − 21 cz )( 21 cx , − 4 c1y , −)2 cz− 2 cz( 21 cx , 0,)cx , 14 cy , − 12 cz
2
(1
)
1
1
( 21 cx , − 21 cy , − 41 cz )( 21 cx , − 4 c1y , −)4 cz− 4 cz( 21 cx , 0,)cx , 14 cy , − 14 cz
2(1)
1( 21 cx , − 21 cy , 0)
( 12 cx , − 4 )cy , 00( 21 cx , 0,)cx , 14 cy , 0
2
(1
)
1
1
( 21 cx , − 21 cy , 41 cz )( 21 cx , − 41cy ,)4 cz
c
4 z( 21 cx , 0,)11
c , c , c
2 x 4 y 4 zになる．これらすべての点を計算するわけだが，ここで (kx , ky , kz ) と (−kx , −ky , −kz ) の波動
関数は複素共役の関係にあることを用いれば，計算するのは (kx , ky , kz ) だけでよい．よって実際に計算に必要な点は，上の例の場合で 64 点から 36 点に減る．
ブリルアン領域内の座標のうち，係数 cx , cy , cz を除いたものを，read.inp に記述していく．波動関数の複素共役の関係を利用して計算を省いた点は，weight を 2 にしておく．
以下に変更した read.inp を示す．8.2 節からの変更点を，青字で示す．
read.inp ファイルの作成例 (初めの 69 行のみ)
1: ![ndisp]
2: 66
3: ![nprocx],[nprocy],[nprocz],[nprock]
4: 1 1 1 1
5: ![xmax],[ymax],[zmax]
6: 3.32d0 3.32d0 3.32d0
7: ![nxmax],[nymax],[nzmax]
8: 10 10 10
9: ! [npxmax],[npymax],[npzmax]
10: -8 8 26
11: ![nso],[socang(1)],[socang(2)]12: 0 0.0d0 0.0d0
13: ![nsym]
14: 0
15: ![neigmx]
16: 8
17: ![natom]
18: 2
19: ![num atcell]
20: -1
21: ![num ppcell]
22: -1
23: ![num ppcell d]
24: -1
25: ![gmax]
26: 20.0d0
27: ![nperi],[npopshow]
28: 3 0
29: ![numkx],[numky],[numkz],[ksym],[kband]
30: 36 1 1 0 0
31: !k point (write all k points)
32: ![X], [Y], [Z], [weight]
33: -0.50D0 -0.50D0 -0.50D0 1
34: -0.25D0 -0.50D0 -0.50D0 2
35:
0.00D0 -0.50D0 -0.50D0 1
36: -0.50D0 -0.25D0 -0.50D0 2
37: -0.25D0 -0.25D0 -0.50D0 2
38:
0.00D0 -0.25D0 -0.50D0 2
39:
0.25D0 -0.25D0 -0.50D0 2
40: -0.50D0 0.00D0 -0.50D0 1
41: -0.25D0 0.00D0 -0.50D0 2
42:
0.00D0 0.00D0 -0.50D0 1
43: -0.50D0 -0.50D0 -0.25D0 2
44: -0.25D0 -0.50D0 -0.25D0 2
45:
0.00D0 -0.50D0 -0.25D0 2
46:
0.25D0 -0.50D0 -0.25D0 2
47: -0.50D0 -0.25D0 -0.25D0 2

29

48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:

-0.25D0 -0.25D0 -0.25D0 2
0.00D0 -0.25D0 -0.25D0 2
0.25D0 -0.25D0 -0.25D0 2-0.50D0 0.00D0 -0.25D0 2-0.25D0 0.00D0 -0.25D0 2
0.00D0 0.00D0 -0.25D0 2
0.25D0 0.00D0 -0.25D0 2-0.50D0 0.25D0 -0.25D0 2-0.25D0 0.25D0 -0.25D0 2
0.00D0 0.25D0 -0.25D0 2
0.25D0 0.25D0 -0.25D0 2-0.50D0 -0.50D0 0.00D0 1-0.25D0 -0.50D0 0.00D0 2
0.00D0 -0.50D0 0.00D0 1-0.50D0 -0.25D0 0.00D0 2-0.25D0 -0.25D0 0.00D0 2
0.00D0 -0.25D0 0.00D0 2
0.25D0 -0.25D0 0.00D0 2-0.50D0 0.00D0 0.00D0 1-0.25D0 0.00D0 0.00D0 2
0.00D0 0.00D0 0.00D0 1![cexco]
VWN![nspv]1![epsvh],[epssd],[ratio diis],[eps scf]
1.0D-8 1.0D-8 0.3d0 1.0D-8![ncgmin],[ncgmax],[ncgres]1 800 801![nprecon cg],[nprecon diis]
1 0![nsdmax],[nkscg]
4 1![ndiismax],[ncgscf],[nretcg]
4 200 2000![nrrz],[nchange]
1 1![eta],[looplimit],[nbrydn],[etamag(1)],[etamag(2)]0.20D0 200 10 0.5d0 2.0d0![tmstep],[nmd start],[nmd end],[ngdiis],[sconst]
0.0D0 0 1 1 0.0d0![biasx],[biasy],[biasz]
0.0D0 0.0D0 0.0D0![tf],[tfmin],[tfmax]
0.003D0 -10.0D0 30.0D0![nmdxyz]
0![chrgd],[npolcon],[polconocc]
0.0D0 0 2.7d0![endjel],[chrjel]
0.0d0 0.0d0![fcut]1.0D-4![npre],[nevhist]
1 3![northo]1

原子座標ファイル atom.xyz は 8.2 章と全く同じでよい．実行スクリプトは，8.2 節と同じである．以下に，電子状態計算結果 mdresult.dat の出力例を記す．total energy が，次の内容と 10桁程度一致していれば，計算は成功である．
mdresult.dat の出力例
computed by kukan8e 03/29/2017-01========== data from read.inp ==========
ndisp
:
66
nprocx:1
nprocy:1nprocz:1
xmax:
3.32000000000000
ymax:
3.32000000000000
zmax:
3.32000000000000
nxmax
:
10
nymax
:
10
nzmax
:
10
npxmax:-8
npymax:8
npzmax
:
26
nso
:
0

30

socang:
nsym:
neigmx:natom:
num_atcell :num_ppcell :num_ppcell_d:
gmaxps:nperi:
npopshow
:
numkx
:
numky
:
numkz
:
ksym
:
kband
:
cexco
:
nspv
:
epsvh
:
epssd
:
ratio_diis :eps_scf:
ncgmin:ncgmax:ncgres:nprecon_cg :nprecon_diis:
nsdmax:nkscg:
ndiismax:
ncgscf:nretcg:nrrz:nchange:eta:
etamag:looplimit:
nbrydn:tmstep:nmd_start:
nmd_end:ngdiis:sconst
:
biasx
:
biasy
:
biasz
:
tf
:
tfmin
:
tfmax
:
nmdxyz:chrgd:
npolcon:polconocc:
endjel:chrjel
:
fcut
:
npre
:
nevhist:northo:eps:
eps_eig_diis:
alambda_diis:
alambda_min :alambda_max :nradmx:nrprjmx:
nprjmx:lsphel:nlmax:
lrhomx:nfiltyp:
psctoff:
nqmx:

0.000E+00 0.000E+00
0
8
2
2
2
2
20.0000000000000
3
0
36
1
1
0
0
vwn
1
1.000000000000000E-008
1.000000000000000E-008
0.300000000000000
1.000000000000000E-008
1
800
801
1
0
4
1
4
200
2000
1
1
0.200000000000000
0.5000E+00 0.2000E+01200
100.000000000000000E+000
01
10.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
0.000000000000000E+000
3.000000000000000E-003-10.0000000000000
30.0000000000000
00.000000000000000E+000
02.70000000000000
0.000000000000000E+000
0.000000000000000E+000
1.000000000000000E-004
1
3
1
1.000000000000000E-016
1.000000000000000E-014
0.500000000000000
0.100000000000000
1.00000000000000
1502
6
18
8
4
25
1
1.05000000000000
400

0.000E+00

31psftrad:10.0000000000000
psctrat:1.00000000000000
psext:
2.00000000000000
filpp: 2.000000000000000E-002
rctpcc:2.00000000000000
nmesh
:
2
npmesh:2veta: 0.200000000000000
new_pwx:
11
new_pwy:11
new_pwz:
11
new_rsx:
11
new_rsy:11
new_rsz:
11
nint1dmax
:
1000
ndec
:
12
nf
:
4
nfdg
:
4
nfh
:
4
nrasy
:
0
nrasz
:
0
nrepx
:
1
zs_pre: 0.500000000000000
pol_pre:1.50000000000000==========================================
Filtering of pseudopotentials is implemented by Fermi distribution.
npxmax, npymax, and npzmax are updated below.
========== data computed by listvecdim ==========
npxmax:13
npymax:
13
npzmax:13====================================================
cell size [lx,ly,lz (a.u.)]
0.664000000D+01
0.664000000D+01
0.664000000D+01
coarse grid spacing [hx,hy,hz (a.u.)]
0.332000000D+00
0.332000000D+00
0.332000000D+00
dense grid spacing [hx,hy,hz (a.u.)]
0.166000000D+00
0.166000000D+00
0.166000000D+00
cutoff energy of wave functions [gmax_coarse,gmax_dense (Ry)]0.895413376D+02
0.358165351D+03
sample k point [skx,sky,skz (2pi/l)] weight-0.500000000D+00 -0.500000000D+00 -0.500000000D+00
1-0.250000000D+00 -0.500000000D+00 -0.500000000D+00
2
0.000000000D+00 -0.500000000D+00 -0.500000000D+00
1-0.500000000D+00 -0.250000000D+00 -0.500000000D+00
2-0.250000000D+00 -0.250000000D+00 -0.500000000D+00
2
0.000000000D+00 -0.250000000D+00 -0.500000000D+00
2
0.250000000D+00 -0.250000000D+00 -0.500000000D+00
2-0.500000000D+00
0.000000000D+00 -0.500000000D+00
1-0.250000000D+00
0.000000000D+00 -0.500000000D+00
2
0.000000000D+00
0.000000000D+00 -0.500000000D+00
1-0.500000000D+00 -0.500000000D+00 -0.250000000D+00
2-0.250000000D+00 -0.500000000D+00 -0.250000000D+00
2
0.000000000D+00 -0.500000000D+00 -0.250000000D+00
2
0.250000000D+00 -0.500000000D+00 -0.250000000D+00
2-0.500000000D+00 -0.250000000D+00 -0.250000000D+00
2-0.250000000D+00 -0.250000000D+00 -0.250000000D+00
2
0.000000000D+00 -0.250000000D+00 -0.250000000D+00
2
0.250000000D+00 -0.250000000D+00 -0.250000000D+00
2-0.500000000D+00
0.000000000D+00 -0.250000000D+00
2-0.250000000D+00
0.000000000D+00 -0.250000000D+00
2
0.000000000D+00
0.000000000D+00 -0.250000000D+00
2
0.250000000D+00
0.000000000D+00 -0.250000000D+00
2-0.500000000D+00
0.250000000D+00 -0.250000000D+00
2-0.250000000D+00
0.250000000D+00 -0.250000000D+00
2
0.000000000D+00
0.250000000D+00 -0.250000000D+00
2
0.250000000D+00
0.250000000D+00 -0.250000000D+00
2-0.500000000D+00 -0.500000000D+00
0.000000000D+00
1-0.250000000D+00 -0.500000000D+00
0.000000000D+00
2
0.000000000D+00 -0.500000000D+00
0.000000000D+00
1-0.500000000D+00 -0.250000000D+00
0.000000000D+00
2-0.250000000D+00 -0.250000000D+00
0.000000000D+00
232

0.000000000D+00 -0.250000000D+00
0.000000000D+00
2
0.250000000D+00 -0.250000000D+00
0.000000000D+00
2-0.500000000D+00
0.000000000D+00
0.000000000D+00
1-0.250000000D+00
0.000000000D+00
0.000000000D+00
2
0.000000000D+00
0.000000000D+00
0.000000000D+00
1number of electron
chargespin
0.200000000D+01
0.000000000D+00
0.100000000D+01
electric field [ex,ey,ez (a.u./a.u.)]
0.000000000D+00
0.000000000D+00
0.000000000D+00
time step (a.u.) (fs)0.000000000D+00
0.000000000D+00
thermo (k)kt (a.u.)
0.947268709D+03
0.300000000D-02
cufoff of force acting on atoms (a.u.)
0.100000000D-03
charge mixing ratio of scf
0.200000000D+00
nrasy,nrasz
00
Wavefunctions & electron density are generated automatically.
total charge
2.00000000000000(electrons)
fermi level 4.990448640675993E-003 (hartree)k,spin,
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

band ,
1
2
3
4
5
6
7
8

eigen value,
0.1276914927E+00
0.1276914948E+00
0.1276914965E+00
0.1276914990E+00
0.1276915002E+00
0.1276915020E+00
0.2960142531E+00
0.2960142645E+00

residual norm ,0.2427870669E-06
0.6582638555E-07
0.5501522944E-07
0.1562261997E-06
0.5639430819E-07
0.4607476703E-07
0.4601822030E-06
0.5889270274E-06

occupation
0.1726650930E-17
0.1726649728E-17
0.1726648729E-17
0.1726647329E-17
0.1726646589E-17
0.1726645584E-17
0.7412918183E-42
0.7412890130E-42····································
36
1
1-0.1171182197E+00
0.1348294474E-06
0.1000000000E+0136
12
0.1998661808E+00
0.1491183997E-06
0.6149620142E-28
36
1
30.1998661837E+00
0.4256397938E-07
0.6149614221E-28
36
140.1998661867E+00
0.1011528362E-06
0.6149608010E-2836150.4733945145E+00
0.4014575081E-03
0.0000000000E+00
36
160.5106322628E+00
0.9608676918E-07
0.0000000000E+00
36
170.5106322637E+00
0.9334456635E-07
0.0000000000E+00
36
180.5106322655E+00
0.7253040098E-07
0.0000000000E+00
Scf= 26dp=
0.5979090750D-08
one electron energy-7.767003258683454E-002 (hartree)
exchange correlation energy -2.95120089971380(hartree)hartreeenergy-1.176486457492671E-003 (hartree)ewald
energy-0.548077326733235(hartree)field correction
0.000000000000000E+000 (hartree)
B_con correction
0.000000000000000E+000 (hartree)energy offset of ps. pot. -3.460192461495514E-004 (hartree)energy jellium/ions
0.000000000000000E+000 (hartree)
helmholtz free energy-3.57990357878086(hartree)Totalenergy-3.57847076473751(hartree)atomic force
1 -0.1117944411968533D-16 -0.1230101590793667D-16
0.3809163105032855D-17
2
0.1117944411968532D-16
0.1230101590793667D-16 -0.3809163105032954D-17
mdstep=1
atom coordinate (input)
0.3320000000000000D+01
0.3320000000000000D+01
0.3320000000000000D+01
30.0000000000000000D+00
0.0000000000000000D+00
0.0000000000000000D+00
3
atom coordinate (output)
0.3320000000000000D+01
0.3320000000000000D+01
0.3320000000000000D+01
30.0000000000000000D+00
0.0000000000000000D+00
0.0000000000000000D+00
3
real time
151.834109208983(sec)com. time
26.8924054023810(sec)com. vol. 0.000000000000000E+000 (byte)

33

0
0

0
0

0
0

12866.00
12866.00

00
00

1
1

1a
2a

0
0

0
0

0
0

12866.00
12866.00

00
00

1
1

1a
2a

A

グリッド幅の決定方法について

グリッド幅 dx は，セルに含まれる原子の中で，最も小さなグリッド幅を要求されるものよりもさらに小さな値を用いる．各原子に必要なグリッド幅は以下のとおりである．

DX(a.u.)
0.27
0.30
0.33
0.35
0.50
0.60

元素
Cu
Fe,W
C,N,F,Cl
O,S
Li,Na
H,Al,Si

なお，上の表に記載されていない元素に関しては連絡先まで確認のこと．

34