実空間差分法に基づく第一原理計算プログラムの使い方(電子状態計算編)

平成 24 年 9 月 5 日
平成 21 年 8 月 8 日
平成 21 年 7 月 26 日
平成 20 年 5 月 30 日
平成 20 年 5 月 28 日

param.f を変更
擬ポテンシャルデータベースに関する記述を追加
v40 − 43c 用に改訂
Li 表面系のコンパイルコマンドに関する記述を変更
初版作成

大阪大学大学院工学研究科精密科学専攻
v40-43c

はじめに

1

本手引書は，大阪大学大学院工学研究科精密科学専攻広瀬研究室で開発した，「実空間差分法に基づく第一原理分子動力学シミュレーションプログラム」の基本的な使い方をまとめたものである．本書で取り上げている，例題の計算はすべて大阪大学大学院工学研究科精密科学専攻
計算物理領域の Pentium4 クラスターで行うことを前提に記述してある．1.1 節に挙げる他の計算機で計算を行う場合，または 1.1 節に載っていない計算機で計算を行う場合は 1.2 節の連絡先まで連絡してください．

1.1

動作確認機種

本プログラムは以下のマシンに関して動作確認を行っている．
・Pentium4 クラスター
・NEC SX7(東北大学情報シナジーセンター)
・HITACHI SR11000(東京大学物性研究所)

1.2

連絡先

本プログラムの不具合，及び本手引書の誤植に関する連絡先〒 565-0871 大阪府吹田市山田丘 2-1
大阪大学大学院工学研究科精密科学専攻
小野倫也
e-mail:ono@prec.eng.osaka-u.ac.jp
TEL 06-6879-7290 FAX 06-6879-7290

2

本プログラムの概要

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

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)

境界条件について本プログラムは，境界条件の設定の違いにより，以下の 4 種類に分けられている．
version40· · ·x, y, z すべて非周期系
version41· · ·x 周期系，y, z は非周期系
version42· · ·x, y 周期系，z は非周期系
version43· · ·x, y, z すべて周期系計算するモデルにより，適切なプログラムを選ぶこと．なお，アスキー入出力データは，すべてのバージョンで互換性があるが，バイナリ入出力データは，version40 のみ波動関数の虚数部分が無いため他のバージョンとの互換性がない．

2.3

フローチャート

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

2.4

擬ポテンシャルについて

本コードの擬ポテンシャルは、物質・材料研究機構 小林一昭氏の作製した擬ポテンシャルデー
タベース NCPS97 を本コード用に入力ファイルを改編したものである。擬ポテンシャルのデータ
ベース NCPS シリーズについては、http://www.geocities.co.jp/Technopolis/4765/INFO/guidee.html
を参照のこと。

2.5

並列処理の表記について

本プログラムは大規模計算を高速に行うことを目的に開発されている．そのためには，複数のプロセッサーを用いて並列処理を行う．本節では，並列化の技法について説明と，本書におけるの並列化に関する言葉の定義を行う．
現在の並列コンピューターは，以下にあげる 2 つに大別できる．

2

コンパイルmake'param.f'は
インプット

コンパイル
make
'param.f'は
インプット

※dec
余分な波動関数を削除

初期電子状態を作成

電子状態計算
構造最適化

結果の解析

pre

kukan

pop

開始

開始

開始
読み込み
計算条件(read.inp)
原子座標(atom.xyz)

初期電子状態
波動関数を計算pre*.f90
計算結果の出力初期電荷分布(ele.0**)
初期波動関数(wfc.0**)

終了

rs.j

読み込み

読み込み

計算条件(read.inp)
原子座標(atom.xyz)
初期電荷分布(ele.0**)
初期波動関数(wfc.0**)

計算条件(read.inp)
原子座標(atom.xyz)
波動関数(wfc.0**)

ポピュレーションを計算
pop3*.f90

原子核からのポテンシャルを計算
pseudo*.f90

計算結果の出力

静電電子間ポテンシャルを計算(Poisson方程式をCG法で解く)
vhcalc*.f90, cg*.f90, kake3d*.f90
収束？
Yatomic, bondポピュレーショ
ン(pop.dat)

N

終了

交換相関ポテンシャルを計算
lda*.f90 or lsda*.f90

NN

RayleighLittz法する？

Y

固有値の
並べ替え

N

電荷分布
が収束？

開始
読み込み

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

波動関数(wfc.0**)

余分な波動関数を削除

構造最適化を続ける？

Y

N
計算結果の出力

計算結果の出力波動関数(wfc.0**)

原子座標(atom.xyz)
全エネルギー etc. (mdresult.dat)
電荷分布(ele.0**)
波動関数(wfc.0**)

終了

N

Y

Rayleigh-Littz法
eig.f or eigch.f※dec

固有値の並べ替えする？

終了

図 1: フローチャート

3

structure optimization loop

Y

CG法を繰り返す？

self consistent loop

Kohn-Sham方程式をC G 法で解く
precon*.f90, kake3dc*.f90, yyinc*.f90

ele
開始
読み込み
電荷分布(ele.0**)

計算結果の出力
アスキーデータの出力(ele.dat)

終了

電荷分布
を更新

1. 共有メモリー型並列コンピューター
1 つのコンピューターに 1 つの大容量メモリーと複数のプロセッサーが載っており，それぞれのプロセッサーはコンピューター内のメモリーを共有している．そのため，あるプロセッサーの演算結果を，同じコンピューター内の他のプロセッサーが瞬時に参照できるというメリットがあるが，1 つのコンピューターに搭載可能なメモリーの容量やプロセッサーの数に限界があり，計算規模が頭打ちしてしまうというデメリットもある．メモリーを共有しているひとつのユニットをノードと呼び，複数のノードを使用すれば，次に述べる 2. のような並列化も可能である．

2. 分散メモリー型並列コンピューター
複数のコンピューター (ノード) がネットワークケーブルでつながれ，クラスターを形成したようなものである．演算はそれぞれのノードが分担して計算を行い，演算結果はネットワーク
ケーブルを経由して他のノードに転送される．初期の並列化は，ほとんどがこのタイプであった．演算量に対しデータ転送が多いプログラムは，通信の部分がボトルネックになってしまいパ
フォーマンスがほとんど向上しないが，データ転送が少ないプログラムであれば，クラスター
内のノードの数に比例してパフォーマンスが向上する．
本プログラムは両方の並列処理が行えるようにコーディングされている．ただし，並列処理の方法により，コンパイルや実行方法に違いがあるため，実行条件を記述する際の例を以下に記す．
複数のプロセッサーで並列化する．· · · 1. の方式で並列化を行う．
複数のノードで並列化する．· · · 2. の方式で並列化を行う1．
並列化する．· · · 1. と 2. の両方の方式で並列化を行う．
m プロセッサーで並列化する．· · · 1. の方式で m 個のプロセッサーで並列化する．
n ノードで並列化する．· · · 2. の方式で n 個のノードで並列化する．
m×n で並列化する．· · ·1 ノードあたり m 個のプロセッサーで，n ノード使用する．合計 m×n
個のプロセッサーを使用して並列化する．本文中，特に断りのない限りは上の例に従う．

簡単な分子の計算を行う

3

本節では，簡単な例題として，すべての方向に非周期境界条件を課したプログラムを使い，1×1 で CO 分子の原子構造最適化計算を行う．

3.1

ディレクトリの作成

まず，各自のホームディレクトリ (ここでは ∼/とする) の下に実行ディレクトリ ∼/CO を作り，さらに ∼/CO の下に source(∼/CO/source) と exe(∼/CO/exe) を作成する．次に ∼/CO/source

1

ほとんどの共有メモリー型コンピュターは，1 ノード内に複数のプロセスを生成して並列処理を行うことができる．この場合，1 ノードしか使用していなくてもこちらの方に含む．

4

にソースファイルを FTP またはファイルコピーによって送り込む．以上の操作ができないとき
は，UNIX の知識のある者に相談する．

3.2

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

コンパイルを行う前に，パラメータファイル param.f を作成する．∼/CO/source の下に param.fのサンプルが作成されているので，サンプルの 1〜24 行目を次の例のように書き換える．プロ
グラムのバージョンによっては，サンプルが変更されている可能性もあるので，全てのパラメー
タをチェックすること．各パラメータの意味を以下に記す．
NPROC(X,Y,Z) · · · 並列計算制御用．4 章参照．本章の計算では，全て 1 にしておく．
(X,Y,Z)MAX· · · -(X,Y,Z)MAX〜+(X,Y,Z)MAX が計算領域になる．非周期境界を設定する場合には，atom.xyz(3.3 参照) で決めた全ての原子座標が境界面から 9〜15(a.u.) 以上の距離を保てるように，計算領域を設定する．
N(X,Y,Z)MAX · · · グリッド数．必ず整数を選ぶこと．2×NXMAX が x 方向のグ
リッド数，XMAX/NXMAX がグリッド幅 DX になる．グリッ
ド幅の設定については Appendix を参照のこと．
NP(X,Y,Z)MAX · · · 擬ポテンシャルの非局所領域に含まれるグリッド数．必ず整数を選ぶこと．(NPXMAX-4)×DX が擬ポテンシャルのカット
オフ半径程度になればよい．擬ポテンシャルのカットオフ半径の設定については Appendix を参照のこと．
NEIGMX· · · 計算する波動関数の数．遷移金属を含まない場合は原子数 ×4，含む場合は原子数 ×9 にしておくこと．
NATOM· · · 原子数．
NUMK
· · · サンプル k ポイントの数．x, y, z すべての方向に非周期境界条件を課した場合は 1 にしておく．
NUMS· · · 交換相関項ポテンシャルにスピンの効果を無視:1，考慮:2．
NMESH· · · Double Grid 法における dense grid の幅 (DX/NMESH)．3 が
望ましい．
NPERI
· · · 境界条件の設定用．非周期境界条件:0，1 方向周期境界条件:1，
2 方向周期境界条件:2，3 方向周期境界条件:3．
NRAS(X,Y,Z)· · · 螺旋境界条件の設定．ここでは説明しないので， 0 にしておく．
NUMK(X,Y,Z) · · · x, y, z 方向のサンプル k ポイントの数．電界印加時の計算にのみ使用するため，本テキストの計算では常に 1 にしておく．
ENDJEL· · · jellium の設定．ここでは説明しないので，0.0D0 にしておく．
CHRJEL· · · jellium の設定．ここでは説明しないので，0.0D0 にしておく．
MACHINE· · · 計算機の種別設定．スカラー:0，ベクトル:1．
param.f ファイルの作成例 (初めの 24 行のみ)
1:!
********** param.f 05/02/2011-01
2:
IMPLICIT REAL*8 (A-H,O-Z)
3:
INCLUDE ’mpif.h’**********

5

4:! number of processors (nodes)
5:
PARAMETER (NPROCX=1,NPROCY=1,NPROCZ=1)
6:! size of supercell
7:
PARAMETER (XMAX=10.0D0,YMAX=10.0D0,ZMAX=10.0D0)
8:! number of grid points
9:
PARAMETER (NXMAX=30,NYMAX=30,NZMAX=30)
10:! number of grid points in core region
11:
PARAMETER (NPXMAX=9,NPYMAX=9,NPZMAX=9)
12:! number of wave functions, number of wave atoms, number grid points for spiral boundary conditions
13:
PARAMETER (NEIGMX=8,NATOM=2,NRASXY=0,NRASXZ=0,NRASYZ=0)
14:! number of sampling k-pints, electron spin
15:
PARAMETER (NUMKX=1,NUMKY=1,NUMKZ=1,NUMK=1,NUMS=1)
16:! number of grid points for double grid technique
17:
PARAMETER (NMESH=3)
18:! switch for periodic boundary conditions
19:
PARAMETER (NPERI=0)
20:! charge of jellium (if you do not use jellium, set CHRJEL to be zero.)
21:
PARAMETER (ENDJEL=0.0D0,CHRJEL=0.0D0)
22:
23:! switch for the type of computer (scalar;0, vector;1)
24:
PARAMETER (MACHINE=0)

コンパイルは，makefile で行う．∼/CO/source の下で>make -f make.dragons lda023と入力する．コンパイルが終了したら ∼/CO/source の下に kukan40，pre，ele，dec，pop ができていることを確認すること．
次に，実行条件入力ファイル read.inp を作成する．FTP などで PC にダウンロードし秀丸等のテキストエディタで変更後，再度 ∼/CO/source に送る．本章の計算では，次の例のように設定
する．4
read.inp の作成例#
CG 法 (P 方程式) の収束判定
CG 法 (K-S 方程式) の収束判定
SCF の収束判定
1.0D-12
1.0D-12
1.0D-8#
CG 法 (P 方程式) 最低反復回数
CG 法 (P 方程式) 最大反復回数
CG 法 (P 方程式) の最出発:N (1 回/N 回)
1
400
500
#
CG 法 (KS 方程式) の最大反復回数 CG 法 (KS 方程式) のスイッチ (0;no 1;yes)
2
1#
SCF における電荷密度の更新比
SCF の最大反復回数
0.10D0200
MD の回数:Ns(start)
MD の回数:Ne(end) (Ne-Ns が MD の回数になる)#
MD の Time step (a.u.)
0.0D0
0
1
#
X 方向の電界 (a.u.)
Y 方向の電界 (a.u.)
Z 方向の電界 (a.u.)
0.0D0
0.0D0
0.0D0
#
k point (個数分縦に並べる)
#
X
Y
Z
ウエイト
0.00D0
0.00D0
0.00D01#
フェルミ温度 (KT) (a.u.)
固有値の最低値 (a.u.)
固有値の最高値 (a.u.)
0.0001D0-1000.0D0
1000.0D0#
Rayleigh-Ritz:N (1 回/N scf 回) 固有値の並べ替え:N (1 回/N scf 回) 動かさない原子の potential の計算を省く (0;no 2;yes)110#
全電子数-全陽子数0.0D0#Force の cutoff
0.5D-22“>”は，プロンプトである．入力しないように．
lda に続く数字はバージョン番号である．version40 の場合は 0，version41 の場合は 1，version42 の場合は 2，
version43 の場合は 3 にする．
4
計算機によっては，日本語で書いた文がエラーの原因になることがある．この場合，コメント文をすべて削除
する．
3

6

3.3

原子座標ファイルの作成

次に原子座標などの情報を記述したファイル atom.xyz を作成する．本章の計算では，以下にのように設定する．各数字の意味は，左から原子の x 座標，y 座標，z 座標，原子番号，x 方向の MD のスイッチ，y 方向の MD のスイッチ，z 方向の MD のスイッチ，原子の質量 [原子量 ×陽子の質量 (1838 a.u.)] である．MD のスイッチが 0 の場合，その方向へは原子を動かさない5．
はじめは電子状態のみを計算するので，全て 0 にしておく．原子が複数個ある場合は，原子の個数分だけ，下に行を連ねていく．したがって原子が N 個あれば，このファイルは N 行必要である．
atom.xyz ファイルの作成例
0.5500000E+00 0.5500000E+00 0.5500000E+00-0.5500000E+00 -0.5500000E+00 -0.5500000E+00

6
8

0
0

0
0

0
0

22056.00
29408.00

CO 分子の平衡原子間距離は 2.1(a.u.) であるが，本章では原子構造最適化計算により電子状
態と平衡原子間距離を求めることを目的としているため，実験値よりも短めに設定しておく．

3.4

実行に必要なファイルのコピー

最後に，実行用ディレクトリ ∼/CO/exe に必要なファイルをコピーする．今，カレントディ
レクトリが ∼/CO/source であるとすると，>cd ..
>cp source/pre exe>cp source/kukan40 exe>cp source/read.inp exe>cp source/atom.xyz exeと入力する．以上で，実行環境は全て整った．

3.5

電子状態計算プログラムの実行

初めて電子状態計算を行うには，あらかじめ入力の波動関数を仮定しておく必要がある6．今，カレントディレクトリが ∼/CO であるとすると，>cd exe>bsub -n 1 -q S -e errfile -o outfile "mpijob mpirun pre"と，入力する．ジョブが実行されると pre により ∼/CO/exe に初期波動関数ファイル ele.000 とwfc.000 が作られる78 ．次に，>bsub -n 1 -q S -e errfile -o outfile "mpijob mpirun kukan40"と入力すると電子状態プログラムが実行される．このプログラムが終了すると，∼/CO/5

現バージョンでは，構造最適化のみ可能である．
計算結果は，仮定した波動関数に依存しない．しかし，本計算は収束計算を行っているため，波動関数の仮定の仕方によっては計算時間が非常に長くなってしまうことがある．ここでは，解に近い波動関数を仮定するようにしている．
7-n, -q, -e, -o の後ろの数字と文字は，それぞれプロセス数，実行キュー名，エラーファイル名，標準出力ファ
イル名である．
8
ジョブの実行確認は “bjobs”で行う．
6

7

exe の下の ele.000 と wfc.000 が更新され，mdresult.dat が作成される．
>cat mdresult.datと，入力すると計算結果を確認できる．以下に mdresult.dat の内容を示す．
mdresult.dat の出力例
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.111111111D+00
0.111111111D+00
0.111111111D+00
nonlocal region [vnlx,vnly,vnlz (a.u.)]
0.600000000D+01
0.600000000D+01
0.600000000D+01
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.315756236D+02
0.100000000D-03
cufoff of the force (a.u.)
0.500000000D-02
mixing ratio of scf
0.100000000D+00
nrasy,nrasz
00
total charge
10.0000000000000
fermi level -0.244140625000000(hartree)atomic charge
13.97034710203453
2
6.02965289796599
1
1
1 eigen value =-0.1138159053D+01
0.1210809104D-09
0.2000000000D+01
1
1
2 eigen value =-0.5232890463D+00
0.9447984641D-11
0.2000000000D+01
1
1
3 eigen value =-0.4875246388D+00
0.7255293815D-11
0.2000000000D+01
1
1
4 eigen value =-0.4875246123D+00
0.2651612204D-10
0.2000000000D+01
1
1
5 eigen value =-0.3252261010D+00
0.1516331236D-09
0.2000000000D+01
1
1
6 eigen value =-0.3528477697D-01
0.1775874064D-09
0.0000000000D+00
1
1
7 eigen value =-0.3528471655D-01
0.8456857513D-08
0.0000000000D+00
1
1
8 eigen value =0.1491283425D-01
0.1573650079D-06
0.0000000000D+00
scf= 107dp=
0.8672886747D-08
one electoron energy-5.92344690344374(hartree)exchange correlation energy
1.57745576380013(hartree)hartreeenergy-29.8137255676793(hartree)ewald
energy12.5967331459555(hartree)field correction
0.000000000000000E+000 (hartree)
helmholtz free energy-21.5629835613674(hartree)totalenergy-21.5629835613674(hartree)mdstep=1
atom coordinate (input)
10.5500000000000000D+00
0.5500000000000000D+00
0.5500000000000000D+00
6
0
2 -0.5500000000000000D+00 -0.5500000000000000D+00 -0.5500000000000000D+00
80atom coordinate (output)
10.5500000000000000D+00
0.5500000000000000D+00
0.5500000000000000D+00
6
0
2 -0.5500000000000000D+00 -0.5500000000000000D+00 -0.5500000000000000D+00
80atomic force
1
0.1995570983038144D+00
0.1995570995934819D+00
0.1995571384730521D+00
2 -0.1995570983038144D+00 -0.1995570995934819D+00 -0.1995571384730521D+00
real time483.160570859909(sec)com. time 3.999900817871094E-002 (sec)com. vol. 0.000000000000000E+000 (byte)======================================================================
This calculation was executed using pseudopotential database NCPS97.
For NCPS database, see K.Kobayashi, Comput. Mater. Sci. 14, 72 (1999).

0
0

0
0

22056.00
29408.00000022056.00
29408.00

mdresult.dat の total energy の部分と atom force の部分が上の内容と一致していれば，計算は成功である.

83.6

原子構造最適化プログラムの実行

電子状態計算の次は，原子構造最適化である．原子構造最適化を行うには atom.xyz，read.inp
を次の例のように書き換える．
read.inp の作成例#
CG 法 (P 方程式) の収束判定
CG 法 (K-S 方程式) の収束判定
SCF の収束判定
1.0D-12
1.0D-12
1.0D-8#
CG 法 (P 方程式) 最低反復回数
CG 法 (P 方程式) 最大反復回数
CG 法 (P 方程式) の最出発:N (1 回/N 回)
1
400
500
#
CG 法 (KS 方程式) の最大反復回数 CG 法 (KS 方程式) のスイッチ (0;no 1;yes)
2
1
SCF の最大反復回数#
SCF における電荷密度の更新比
0.10D0200#
MD の Time step (a.u.)
MD の回数:Ns(start)
MD の回数:Ne(end) (Ne-Ns が MD の回数になる)
80.0D0040
Y 方向の電界 (a.u.)
Z 方向の電界 (a.u.)#
X 方向の電界 (a.u.)
0.0D0
0.0D0
0.0D0
#
k point (個数分縦に並べる)
ウエイト
#
X
Y
Z
0.00D0
0.00D0
0.00D01#
フェルミ温度 (KT) (a.u.)
固有値の最低値 (a.u.)
固有値の最高値 (a.u.)
0.0001D0-1000.0D0
1000.0D0#
Rayleigh-Ritz:N (1 回/N scf 回) 固有値の並べ替え:N (1 回/N scf 回) 動かさない原子の potential の計算を省く (0;no 2;yes)110#
全電子数-全陽子数0.0D0#Force の cutoff
0.5D-2
atom.xyz の作成例
0.6000000E+00 0.6000000E+00 0.6000000E+00-0.6000000E+00 -0.6000000E+00 -0.6000000E+00

6
8

1
1

1
1

1
1

22056.00
29408.00

コピーが終了したら，次に初期波動関数を作成する．しかし，今回は新たに波動関数を仮定しなくても 3.5 節で行った結果の方が解に近いため，前節の結果を利用する．続いて>bsub -n 1 -q S -e errfile -o outfile "mpijob mpirun kukan40"と入力すれば原子構造最適化計算のプログラムが実行される．
プログラムが終了すると，mdresult.dat に各 MD(原子構造最適化) ステップごとの電子状態の情報が出力される．mdresult.dat の最後の部分 MDSTEP=40 のところで，total energy が，
約-21.600(a.u.)，atom coordinate (output) の欄を見て 2 つの原子の距離が約 2.1(a.u.) になっていれば成功である．

4

複数のプロセッサで並列化する

2.5 節で述べたように，並列化にはノード内のみでの並列化と複数のノードを用いた並列化がある．本章では，両方の並列化方法について紹介する．なお，本章の例題の解答は 3 章の解答と同じである．ただし，並列化による演算順序の変更により，演算結果が 3 章の結果と若干異なってしまうこともある．

4.1

2×1 で並列処理する

まず初めに，ひとつのノード内で並列処理を行う．基本的に 3 章とほとんど変わらない．3.2
節の9

>make -f make.dragons lda0の代わりに，>make -f make.dragonp lda0と入力する．これら以外は，3 章と全く同じである．

4.2

1 ノード内で 1×2 で並列処理する

1 ノード内に 2 プロセス発生させて並列処理を行う方法を説明する．実際に 2 ノード以上を使って並列化させる場合は，4.3 節を参照すること．
まず，3 章の 3.1 節までの作業を行った後，カレントディレクトリを ∼/CO/source にする．
param.f をエディタで開き，NPROCX=2 にする．これにより，-XMAX<0 の領域はノード 0 が，
0<-XMAX の領域はノード 1 が受け持つことになる．ここでは行わないが，NPROC(Y,Z) を変
更する場合も同様の分割を行う．トータルのノード数は NPROCX×NPROCY×NPROCZ 個である．並列化を行う場合は，2×N(X,Y,Z)MAX がそれぞれ NPROC(X,Y,Z) で割り切れ，さらに，NPROC(X,Y,Z) は 1 か偶数でなくてはならないという制約がある．
param.f の作成例 (初めの 24 行のみ)
1:!
********** param.f 05/02/2011-01 **********
2:
IMPLICIT REAL*8 (A-H,O-Z)
3:
INCLUDE ’mpif.h’
4:! number of processors (nodes)
5:
PARAMETER (NPROCX=2,NPROCY=1,NPROCZ=1)
6:! size of supercell
7:
PARAMETER (XMAX=10.0D0,YMAX=10.0D0,ZMAX=10.0D0)
8:! number of grid points
9:
PARAMETER (NXMAX=30,NYMAX=30,NZMAX=30)
10:! number of grid points in core region
11:
PARAMETER (NPXMAX=9,NPYMAX=9,NPZMAX=9)
12:! number of wave functions, number of wave atoms, number grid points for spiral boundary conditions
13:
PARAMETER (NEIGMX=8,NATOM=2,NRASXY=0,NRASXZ=0,NRASYZ=0)
14:! number of sampling k-pints, electron spin
15:
PARAMETER (NUMKX=1,NUMKY=1,NUMKZ=1,NUMK=1,NUMS=1)
16:! number of grid points for double grid technique
17:
PARAMETER (NMESH=3)
18:! switch for periodic boundary conditions
19:
PARAMETER (NPERI=0)
20:! charge of jellium (if you do not use jellium, set CHRJEL to be zero.)
21:
PARAMETER (ENDJEL=0.0D0,CHRJEL=0.0D0)
22:
23:! switch for the type of computer (scalar;0, vector;1)
24:
PARAMETER (MACHINE=0)

その他のファイルの内容は 3 章と全く同じでよい．コンパイル，ファイルのコピーは 3 章と同じ作業を行う．∼/CO/exe で>bsub -n 2 -R "span[ptile=2]" -q S -e errfile -o outfile "mpijob mpirun pre"と，入力すると初期波動関数が生成される．ここで，出てくる-n，ptile=の後ろの数字は，それぞれノード数，ノードあたりのプロセッサー数である．したがって，m×n で並列化するときの nがノード数，m がノードあたりのプロセッサー数に相当する9 ．以上の作業が終われば ∼/CO/exeに ele.000，ele.001，wfc.000，wfc.001 の 4 つができているはずである．ここで，ele.00*，wfc.00*の*の数字はノードに割り振った番号である．ele.000 と wfc.000 はプロセッサー 0 が入出力し，ele.001 と wfc.001 はプロセッサー 1 が入出力する．
9

何も指定しない場合は”span[ptile=2]”が仮定される．

10

次に>bsub -n 2 -R "span[ptile=2]" -q S -e errfile -o outfile "mpijob mpirun kukan40"と入力すれば，電子状態計算プログラムが実行される．出力結果は，3 章と同じである．

4.3

2 ノードで 1×2 で並列処理する

今度は，実際に 2 ノード以上を使って並列化させる方法を説明する．それぞれ>bsub -n 2 -R "span[ptile=1]" -q P4 -e errfile -o outfile "mpijob mpirun pre">bsub -n 2 -R "span[ptile=1]" -q P4 -e errfile -o outfile "mpijob mpirun kukan40"と入力すれば，初期波動関数作成，電子状態計算のプログラムが実行される．10

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

5

4 章までの計算は，孤立したクラスターの電子状態を求める計算であった．本章では周期構造をもった結晶や表面の計算方法を説明する．本章の計算は，モデルサイズが小さいため全て
1 プロセッサーで行う．

5.1

スーパーセル

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

5.2

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

結晶は，x，y ，z の全ての方向に周期性を持っている．したがって，結晶の電子状態を計算するときに用いるプログラムは version43 である．ディレクトリの作成は 3 章と全く同じである．
3 章と同じく，まず param.f を作成する．
param.f の作成例 (初めの 24 行のみ)
1:!
********** param.f 05/02/2011-01 **********
2:
IMPLICIT REAL*8 (A-H,O-Z)
3:
INCLUDE ’mpif.h’
4:! number of processors (nodes)
5:
PARAMETER (NPROCX=1,NPROCY=1,NPROCZ=1)
6:! size of supercell
7:
PARAMETER (XMAX=3.32D0,YMAX=3.32D0,ZMAX=3.32D0)
10

2 ノード以上使う場合は，並列計算用のキュー P4 または P8 で計算を行う．

11

*
+
,
)

図 2: スーパーセル

8:! number of grid points
9:
PARAMETER (NXMAX=10,NYMAX=10,NZMAX=10)
10:! number of grid points in core region
11:
PARAMETER (NPXMAX=10,NPYMAX=10,NPZMAX=10)
12:! number of wave functions, number of wave atoms, number grid points for spiral boundary conditions
13:
PARAMETER (NEIGMX=8,NATOM=2,NRASXY=0,NRASXZ=0,NRASYZ=0)
14:! number of sampling k-pints, electron spin
15:
PARAMETER (NUMKX=1,NUMKY=1,NUMKZ=1,NUMK=1,NUMS=1)
16:! number of grid points for double grid technique
17:
PARAMETER (NMESH=3)
18:! switch for periodic boundary conditions
19:
PARAMETER (NPERI=3)
20:! charge of jellium (if you do not use jellium, set CHRJEL to be zero.)
21:
PARAMETER (ENDJEL=0.0D0,CHRJEL=0.0D0)
22:
23:! switch for the type of computer (scalar;0, vector;1)
24:
PARAMETER (MACHINE=0)

リチウムの格子定数は 6.64(a.u.) なので，計算領域 2×(X,Y,Z)MAX が 6.64(a.u.) になるように (X,Y,Z)MAX=3.32 とする．コンパイルは，>make -f make.dragons lda3と入力する．その後，実行用ディレクトリに必要なファイル pre，kukan43 をコピーする．
read.inp の例#
CG 法 (P 方程式) の収束判定
CG 法 (K-S 方程式) の収束判定
SCF の収束判定
1.0D-12
1.0D-12
1.0D-8#
CG 法 (P 方程式) 最低反復回数
CG 法 (P 方程式) 最大反復回数
CG 法 (P 方程式) の最出発:N (1 回/N 回)
1
400
500
#
CG 法 (KS 方程式) の最大反復回数 CG 法 (KS 方程式) のスイッチ (0;no 1;yes)
2
1#
SCF における電荷密度の更新比
SCF の最大反復回数
0.10D0200#
MD の Time step (a.u.)
MD の回数:Ns(start)
MD の回数:Ne(end) (Ne-Ns が MD の回数になる)

12

#
#
#
#
#
#
#

0.0D0
0
1
Y 方向の電界 (a.u.)
Z 方向の電界 (a.u.)
X 方向の電界 (a.u.)
0.0D0
0.0D0
0.0D0
k point (個数分縦に並べる)
X
Y
Z
ウエイト
0.0D0
0.0D0
0.0D0
1
フェルミ温度 (KT) (a.u.)
固有値の最低値 (a.u.)
固有値の最高値 (a.u.)
0.0001D0-1000.0D0
1000.0D0
Rayleigh-Ritz:N (1 回/N scf 回) 固有値の並べ替え:N (1 回/N scf 回) 動かさない原子の potential の計算を省く (0;no 2;yes)110
全電子数-全陽子数0.0D0Force の cutoff
0.5D-2

read.inp を作成したら，次に原子座標ファイル atom.xyz を作成する．5.1 節で述べたように，
リチウムは体心立方晶なので最小単位セルに必要な原子数は 2 個であり，座標ファイルは次の例のようになる．
atom.xyz の作成例
-0.3320000E+01 -0.3320000E+01 -0.3320000E+01
0.0000000E+00 0.0000000E+00 0.0000000E+00

3
3

0
0

0
0

0
0

12866.00
12866.00

ここで，1 行目の原子の座標の正負は任意だが，プログラムの中で自動的に上のように変更されるため，mdresult.dat などの出力ファイルには全て変更された後の値が出力される11 ．ファ
イルのコピー，実行は 3 章と同じであるが， kukan40 が kukan43 に変更されていることに注意する．最後に，電子状態計算の実行結果 mdresult.dat を示す．total energy の部分が次の内容と一致していれば，計算は成功である．
mdresult.dat の出力例
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.110666667D+00
0.110666667D+00
0.110666667D+00
nonlocal region [vnlx,vnly,vnlz (a.u.)]
0.664000000D+01
0.664000000D+01
0.664000000D+01
sample k point [skx,sky,skz (2pi/l)] weight
0.000000000D+00
0.000000000D+00
0.000000000D+00 1
number 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.315756236D+02
0.100000000D-03
cufoff of the force (a.u.)
0.500000000D-02
mixing ratio of scf
0.100000000D+00
nrasy,nrasz
00
total charge
2.00000000000000
fermi level 0.000000000000000E+000 (hartree)
atomic charge
11.00000803902034
2 0.999991960979647
1
1
1 eigen value =-0.1178122856D+00
0.2073869768D-11
1
1
2 eigen value =0.1995609987D+00
0.5303687887D-11
1
1
3 eigen value =0.1995609987D+00
0.9593338376D-11
1
1
4 eigen value =0.1995609987D+00
0.9980901767D-11
1
1
5 eigen value =0.4638158146D+00
0.2989061296D-10
1
1
6 eigen value =0.5167489683D+00
0.2050917370D-10
1
1
7 eigen value =0.5167489683D+00
0.2063266407D-10
1
1
8 eigen value =0.5167489683D+00
0.4127189397D-10
110.2000000000D+01
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00

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

13

scf=

68
dp=
0.9340968942D-08
one electoron energy-0.235624571166907(hartree)exchange correlation energy -5.386069972278303E-003 (hartree)
hartreeenergy-1.621170951704844E-003 (hartree)ewald
energy-0.548077326733235(hartree)field correction
0.000000000000000E+000 (hartree)
helmholtz free energy-0.790709138824125(hartree)totalenergy-0.790709138824125(hartree)mdstep=1
atom coordinate (input)
1 -0.3320000000000000D+01 -0.3320000000000000D+01 -0.3320000000000000D+01
2
0.0000000000000000D+00
0.0000000000000000D+00
0.0000000000000000D+00
atom coordinate (output)
1 -0.3320000000000000D+01 -0.3320000000000000D+01 -0.3320000000000000D+01
2
0.0000000000000000D+00
0.0000000000000000D+00
0.0000000000000000D+00
atomic force
1
0.8770383276775760D-17 -0.1502306554972619D-16 -0.4150725315750125D-16
2 -0.8770383276775760D-17
0.1502306554972619D-16
0.4150725315750125D-16
real time
28.4494550228119(sec)com. time 9.999990463256836E-003 (sec)com. vol. 0.000000000000000E+000 (byte)======================================================================
This calculation was executed using pseudopotential database NCPS97.
For NCPS database, see K.Kobayashi, Comput. Mater. Sci. 14, 72 (1999).

5.3

3
3

0
0

0
0

0
0

12866.00
12866.00

3
3

0
0

0
0

0
0

12866.00
12866.00

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

次に，表面の電子状態計算を行う．第一原理電子状態計算手法で，表面のように面垂直方向の一方は真空，反対方向には結晶が無限に続く物体の計算方法は，現在のところ実用化されていない．したがって，表面の計算を行うには，原子層を何層か積んだ薄膜を表面に見立てて計算している．ここでは，リチウム (001)4 原子層薄膜の電子状態計算を行う．薄膜は，x，y 方向に周期性を持っている．したがって，薄膜の電子状態を計算するときに用いるプログラムは
version42 である．ディレクトリの作成は 3 章と全く同じである．作成後，3 章と同じく，まずparam.f を作成する．
param.f の作成例 (初めの 24 行のみ)
1:!
********** param.f 05/02/2011-01 **********
2:
IMPLICIT REAL*8 (A-H,O-Z)
3:
INCLUDE ’mpif.h’
4:! number of processors (nodes)
5:
PARAMETER (NPROCX=1,NPROCY=1,NPROCZ=1)
6:! size of supercell
7:
PARAMETER (XMAX=3.32D0,YMAX=3.32D0,ZMAX=13.28D0)
8:! number of grid points
9:
PARAMETER (NXMAX=10,NYMAX=10,NZMAX=40)
10:! number of grid points in core region
11:
PARAMETER (NPXMAX=10,NPYMAX=10,NPZMAX=10)
12:! number of wave functions, number of wave atoms, number grid points for spiral boundary conditions
13:
PARAMETER (NEIGMX=16,NATOM=4,NRASXY=0,NRASXZ=0,NRASYZ=0)
14:! number of sampling k-pints, electron spin
15:
PARAMETER (NUMKX=1,NUMKY=1,NUMKZ=1,NUMK=1,NUMS=1)
16:! number of grid points for double grid technique
17:
PARAMETER (NMESH=3)
18:! switch for periodic boundary conditions
19:
PARAMETER (NPERI=2)
20:! charge of jellium (if you do not use jellium, set CHRJEL to be zero.)
21:
PARAMETER (ENDJEL=0.0D0,CHRJEL=0.0D0)
22:
23:! switch for the type of computer (scalar;0, vector;1)
24:
PARAMETER (MACHINE=0)

リチウムの格子定数は 6.64(a.u.) なので，計算領域 2×(X,Y,Z)MAX が 6.64(a.u.) になるように (X,Y,Z)MAX=3.32 とする．また，周期性のない z 方向は，境界面を原子から離しておかな

14ければならないので，長めに設定しておく．read.inp は，次の例のように作成する．
read.inp の例#
CG 法 (P 方程式) の収束判定
CG 法 (K-S 方程式) の収束判定
SCF の収束判定
1.0D-12
1.0D-12
1.0D-8
CG 法 (P 方程式) 最大反復回数
CG 法 (P 方程式) の最出発:N (1 回/N 回)#
CG 法 (P 方程式) 最低反復回数
1
400
500
#
CG 法 (KS 方程式) の最大反復回数 CG 法 (KS 方程式) のスイッチ (0;no 1;yes)
2
1#
SCF における電荷密度の更新比
SCF の最大反復回数
0.05D0200#
MD の Time step (a.u.)
MD の回数:Ns(start)
MD の回数:Ne(end) (Ne-Ns が MD の回数になる)
0.0D0
0
1
Y 方向の電界 (a.u.)
Z 方向の電界 (a.u.)#
X 方向の電界 (a.u.)
0.0D0
0.0D0
0.0D0
#
k point (個数分縦に並べる)
#
X
Y
Z
ウエイト
0.0D0
0.0D0
0.0D0
1
#
フェルミ温度 (KT) (a.u.)
固有値の最低値 (a.u.)
固有値の最高値 (a.u.)
0.0001D0-1000.0D0
1000.0D0#
Rayleigh-Ritz:N (1 回/N scf 回) 固有値の並べ替え:N (1 回/N scf 回) 動かさない原子の potential の計算を省く (0;no 2;yes)110#
全電子数-全陽子数0.0D0#Force の cutoff
0.5D-2

read.inp の作成と下記のコマンド >make -f make.dragons lda2
でコンパイルを行った後，原子座標ファイル atom.xyz を作成する．最小単位セルに必要な原子数は 4 個であり，座標ファイルは次の例のようになる．コンパイル用のコマンドが lda3 ではな
く lda2 となっているところに注意!
atom.xyz の作成例
-0.3320000E+01 -0.3320000E+01 -0.4980000E+01
0.0000000E+00 0.0000000E+00 -0.1660000E+01-0.3320000E+01 -0.3320000E+01 0.1660000E+01
0.0000000E+00 0.0000000E+00 0.4980000E+01

3
3
3
3

0
0
0
0

0
0
0
0

0
0
0
0

12866.00
12866.00
12866.00
12866.00

ファイルのコピー，実行は 3 章と同じであるが，kukan40 が kukan42 に変更されていることに注意する．最後に，電子状態計算の実行結果 mdresult.dat を示す．total energy の部分と atom
force の部分が次の内容と一致していれば，計算は成功である．
mdresult.dat の出力例
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.110666667D+00
0.110666667D+00
0.110666667D+00
nonlocal region [vnlx,vnly,vnlz (a.u.)]
0.664000000D+01
0.664000000D+01
0.664000000D+01
sample k point [skx,sky,skz (2pi/l)] weight
0.000000000D+00
0.000000000D+00
0.000000000D+00 1
number 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.315756236D+02
0.100000000D-03
cufoff of the force (a.u.)
0.500000000D-02
mixing ratio of scf
0.500000000D-01
nrasy,nrasz

15

0
0
total charge
4.00000000000000
fermi level -0.152587890625000(hartree)atomic charge
1 0.967249575076897
2
1.03275143049763
31.03275062842039
4 0.967248366005066
1
1
1 eigen value =-0.2196704018D+00
0.6853601472D-12
0.2000000000D+01
1
1
2 eigen value =-0.1768975089D+00
0.6728713433D-11
0.2000000000D+01
1
1
3 eigen value =-0.1266419235D+00
0.1790864275D-11
0.0000000000D+00
1
1
4 eigen value =-0.7439854043D-01
0.3724121012D-11
0.0000000000D+00
1
1
5 eigen value =-0.1110562288D-01
0.3347725756D-11
0.0000000000D+00
1
1
6 eigen value =0.8068451685D-01
0.1832362815D-09
0.0000000000D+00
1
1
7 eigen value =0.1004363050D+00
0.7532717841D-12
0.0000000000D+00
1
1
8 eigen value =0.1004363050D+00
0.9460792198D-12
0.0000000000D+00
1
1
9 eigen value =0.1219676673D+00
0.2811960464D-10
0.0000000000D+00
1
1 10 eigen value =0.1519722048D+00
0.7256278434D-12
0.0000000000D+00
1
1 11 eigen value =0.1519722048D+00
0.4901169071D-11
0.0000000000D+00
1
1 12 eigen value =0.2269131077D+00
0.5608559071D-12
0.0000000000D+00
1
1 13 eigen value =0.2269131078D+00
0.2327435943D-10
0.0000000000D+00
1
1 14 eigen value =0.2383325693D+00
0.1261513382D-08
0.0000000000D+00
1
1 15 eigen value =0.3192690808D+00
0.1792683991D-10
0.0000000000D+00
1
1 16 eigen value =0.3192690809D+00
0.6374910847D-09
0.0000000000D+00
scf= 142dp=
0.9268376643D-08
one electoron energy-0.793135821317456(hartree)exchange correlation energy -2.020703521583506E-002 (hartree)
hartreeenergy5.33767415324447(hartree)ewald
energy-5.96650333694423(hartree)field correction
0.000000000000000E+000 (hartree)
helmholtz free energy-1.44217204023306(hartree)totalenergy-1.44217204023306(hartree)mdstep=1
atom coordinate (input)
1 -0.3320000000000000D+01 -0.3320000000000000D+01 -0.4980000000000000D+01
3
02
0.0000000000000000D+00
0.0000000000000000D+00 -0.1660000000000000D+01
30
3 -0.3320000000000000D+01 -0.3320000000000000D+01
0.1660000000000000D+01
3
04
0.0000000000000000D+00
0.0000000000000000D+00
0.4980000000000000D+01
30
atom coordinate (output)
1 -0.3320000000000000D+01 -0.3320000000000000D+01 -0.4980000000000000D+01
3
02
0.0000000000000000D+00
0.0000000000000000D+00 -0.1660000000000000D+01
30
3 -0.3320000000000000D+01 -0.3320000000000000D+01
0.1660000000000000D+01
3
04
0.0000000000000000D+00
0.0000000000000000D+00
0.4980000000000000D+01
30atomic force
1 -0.8407493306262627D-16
0.2804451702379245D-16 -0.2554335266887098D-02
2 -0.5458557017826061D-16
0.2985089337772781D-18
0.1253766698119486D-02
3
0.4371815274932037D-17
0.4886340272487307D-16 -0.1253688155347638D-02
40.1342886879659548D-15 -0.7720642868244279D-16
0.2554256724115250D-02
real time
277.704668998718(sec)com. time 3.999137878417969E-002 (sec)com. vol. 0.000000000000000E+000 (byte)======================================================================
This calculation was executed using pseudopotential database NCPS97.
For NCPS database, see K.Kobayashi, Comput. Mater. Sci. 14, 72 (1999).

6

0
0
0
0

0
0
0
0

12866.00
12866.00
12866.00
12866.00

0
0
0
0

0
0
0
0

12866.00
12866.00
12866.00
12866.00

第一原理計算の高精度化

前節までの計算は，最も簡単，言い換えれば粗い近似の第一原理電子状態計算であった．さらに精度を向上させる方法として，スピンの考慮や，周期系モデルにおいてサンプル k ポイン
トの数を増やすといった方法がある．本章では，これらの方法について説明する．

16

6.1

スピンを考慮した電子状態計算を行う

ほとんどの原子やイオン，O2 分子，S2 分子などはスピン分極しており，このようなものを扱
うにはスピンを考慮した計算を行う必要がある．ここでは，O2 分子の電子状態計算を例にして説明する．
param.f，read.inp，atom.xyz を次のように作成する．
param.f の作成例 (初めの 24 行のみ)
1:!
********** param.f 05/02/2011-01 **********
2:
IMPLICIT REAL*8 (A-H,O-Z)
3:
INCLUDE ’mpif.h’
4:! number of processors (nodes)
5:
PARAMETER (NPROCX=1,NPROCY=1,NPROCZ=1)
6:! size of supercell
7:
PARAMETER (XMAX=8.0D0,YMAX=8.0D0,ZMAX=8.0D0)
8:! number of grid points
9:
PARAMETER (NXMAX=24,NYMAX=24,NZMAX=24)
10:! number of grid points in core region
11:
PARAMETER (NPXMAX=9,NPYMAX=9,NPZMAX=9)
12:! number of wave functions, number of wave atoms, number grid points for spiral boundary conditions
13:
PARAMETER (NEIGMX=8,NATOM=2,NRASXY=0,NRASXZ=0,NRASYZ=0)
14:! number of sampling k-pints, electron spin
15:
PARAMETER (NUMKX=1,NUMKY=1,NUMKZ=1,NUMK=1,NUMS=2)
16:! number of grid points for double grid technique
17:
PARAMETER (NMESH=3)
18:! switch for periodic boundary conditions
19:
PARAMETER (NPERI=0)
20:! charge of jellium (if you do not use jellium, set CHRJEL to be zero.)
21:
PARAMETER (ENDJEL=0.0D0,CHRJEL=0.0D0)
22:
23:! switch for the type of computer (scalar;0, vector;1)
24:
PARAMETER (MACHINE=0)
read.inp の作成例#
CG 法 (P 方程式) の収束判定
CG 法 (K-S 方程式) の収束判定
SCF の収束判定
1.0D-12
1.0D-12
1.0D-8#
CG 法 (P 方程式) 最低反復回数
CG 法 (P 方程式) 最大反復回数
CG 法 (P 方程式) の最出発:N (1 回/N 回)
1
400
500
#
CG 法 (KS 方程式) の最大反復回数 CG 法 (KS 方程式) のスイッチ (0;no 1;yes)
2
1#
SCF における電荷密度の更新比
SCF の最大反復回数
0.10D0200#
MD の Time step (a.u.)
MD の回数:Ns(start)
MD の回数:Ne(end) (Ne-Ns が MD の回数になる)
0.0D0
0
1
#
X 方向の電界 (a.u.)
Y 方向の電界 (a.u.)
Z 方向の電界 (a.u.)
0.0D0
0.0D0
0.0D0
#
k point (個数分縦に並べる)
#
X
Y
Z
ウエイト
0.00D0
0.00D0
0.00D01#
フェルミ温度 (KT) (a.u.)
固有値の最低値 (a.u.)
固有値の最高値 (a.u.)
0.0001D0-1000.0D0
1000.0D0#
Rayleigh-Ritz:N (1 回/N scf 回) 固有値の並べ替え:N (1 回/N scf 回) 動かさない原子の potential の計算を省く (0;no 2;yes)110#
全電子数-全陽子数0.0D0#Force の cutoff
0.5D-2
atom.xyz の作成例
0.1100000E+01 0.0000000E+00-0.1100000E+01 0.0000000E+00

0.0000000E+00
0.0000000E+00

8
8

0
0

0
0

0
0

29408.00
29408.00

基本的に，3 章とほとんど変わりはないが，param.f の NUMS が 2 になっていることに注意すること．また，コンパイルは，>make -f make.dragons lsda017と入力する．lda ではなく lsda となっているところに注意! これら以外は，3 章と同じである．
出力結果 mdresult.dat を次に示す．
mdresult.dat の出力例
cell size [lx,ly,lz (a.u.)]
0.160000000D+02
0.160000000D+02
0.160000000D+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.111111111D+00
0.111111111D+00
0.111111111D+00
nonlocal region [vnlx,vnly,vnlz (a.u.)]
0.600000000D+01
0.600000000D+01
0.600000000D+01
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.315756236D+02
0.100000000D-03
cufoff of the force (a.u.)
0.500000000D-02
mixing ratio of scf
0.100000000D+00
nrasy,nrasz
00
total charge
12.0000000000000
fermi level -0.183105468750000(hartree)atomic charge
16.00000124907822
2
5.99999875092187
1
1
1 eigen value =-0.1250449843D+01
0.9853324659D-12
0.1000000000D+01
1
1
2 eigen value =-0.7358900948D+00
0.7557513382D-12
0.1000000000D+01
1
1
3 eigen value =-0.5204972613D+00
0.8099017100D-11
0.1000000000D+01
1
1
4 eigen value =-0.5204972358D+00
0.7842062476D-11
0.1000000000D+01
1
1
5 eigen value =-0.5145541678D+00
0.7139757295D-12
0.1000000000D+01
1
1
6 eigen value =-0.2402934242D+00
0.6102861955D-11
0.1000000000D+01
1
1
7 eigen value =-0.2402933938D+00
0.5418164983D-11
0.1000000000D+01
1
1
8 eigen value =0.6012055962D-01
0.3413178094D-04
0.0000000000D+00
1
2
1 eigen value =-0.1201854523D+01
0.9270068382D-12
0.1000000000D+01
1
2
2 eigen value =-0.6707397602D+00
0.4974411436D-11
0.1000000000D+01
1
2
3 eigen value =-0.4708171441D+00
0.7836487258D-12
0.1000000000D+01
1
2
4 eigen value =-0.4561860939D+00
0.1666281538D-10
0.1000000000D+01
1
2
5 eigen value =-0.4561860124D+00
0.1540088571D-10
0.1000000000D+01
1
2
6 eigen value =-0.1632235824D+00
0.6376329410D-11
0.0000000000D+00
1
2
7 eigen value =-0.1632234930D+00
0.1620357289D-10
0.0000000000D+00
1
2
8 eigen value =0.2828060391D-01
0.1443861832D-03
0.0000000000D+00
scf= 106dp=
0.8882018574D-08
one electoron energy-7.27825895479820(hartree)exchange correlation energy
2.08449247570767(hartree)hartreeenergy-42.8905266088814(hartree)ewald
energy16.3636363636364(hartree)field correction
0.000000000000000E+000 (hartree)
helmholtz free energy-31.7206567243356(hartree)totalenergy-31.7206567243356(hartree)mdstep=1
atom coordinate (input)
1
0.1100000000000000D+01
0.0000000000000000D+00
0.0000000000000000D+00
80
2 -0.1100000000000000D+01
0.0000000000000000D+00
0.0000000000000000D+00
80atom coordinate (output)
1
0.1100000000000000D+01
0.0000000000000000D+00
0.0000000000000000D+00
80
2 -0.1100000000000000D+01
0.0000000000000000D+00
0.0000000000000000D+00
80atomic force
1
0.7898879336295694D-01
0.2620224698811231D-14
0.1005343984723030D-13
2 -0.7898879336295694D-01 -0.2620224698811232D-14 -0.1005343984723030D-13
real time249.685297966003(sec)com. time 2.999567985534668E-002 (sec)com. vol. 0.000000000000000E+000 (byte)======================================================================
This calculation was executed using pseudopotential database NCPS97.
For NCPS database, see K.Kobayashi, Comput. Mater. Sci. 14, 72 (1999).

18

0
0

0
0

29408.00
29408.00000029408.00
29408.00

引き続き，原子構造最適化を行うには，同じく 3 章の要領，read.inp と atom.xyz の書き換えを行った後，kukan40 を実行する．

6.2

サンプル k ポイントを増やす

結晶の中には，あらゆる波数ベクトル (周期) を持った電子が存在していると考えることができる．すなわち，波動関数は格子と同じ周期のもの (波数ベクトルの成分が整数) は，ブリルア
ン領域内の波動関数のひとつにすぎない．5.2 章や 5.3 章で行った計算は，数ある波動関数の中から波数ベクトルの成分が整数のものだけを扱った．しかし，厳密な議論では，ブリルアン領域
内の全ての波動関数を扱わなければならない．サンプル k ポイントを増やすということは，ブ
リルアン領域をメッシュで分割し，全ての波動関数を考慮に入れようというものであり，メッ
シュが細かくなれば細かくなるほど精度が増す．しかし，分割数が多くなればなるほど計算に要する時間も増すので，計算するモデルの大きさに等も考慮して，サンプル k ポイントの数を決めなければならない．ここではリチウム結晶の電子状態を 64k ポイントとって計算する．

read.inp の作成 サンプル k ポイントを増やすには，計算条件設定ファイル read.inp に計算する k ポイントの値を，手で入力する必要がある．本計算では，第一ブリルアン領域を 4 × 4 × 4のメッシュで分割する．それぞれのメッシュ点のブリルアン領域での座標は，cx = L2πx ，cy = L2πy ，cz = L2πz とすると
 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 , 11
− 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,11
− 2 cx , 4 cy , 4 cz ,

 1

1
1
− 41 cx , − 21 cy , − 21 cz  ,− 41 cx , − 4 c1y , −2 cz ,− 2 cz ,− 41 cx , 0,
− 4 cx , 14 cy , − 21 cz ,
 1

1
1
− 41 cx , − 21 cy , − 41 cz  ,− 41 cx , − 4 c1y , −4 cz ,− 4 cz ,− 41 cx , 0,
− 4 cx , 14 cy , − 41 cz , 1
1
− 41 cx , − 21 cy , 0 ,
− 41 cx , − 4 cy , 0 ,0 ,− 41 cx , 0,
− 4 cx , 14 cy , 0 ,
 1

1
1
− 41 cx , − 21 cy , 41 cz  ,− 41 cx , − 41cy ,4 cz ,c ,4 z
− 41 cx , 0,11
− 4 cx , 4 cy , 4 cz ,11
0, − 21 cy , − 21 cz  ,0, − 4 c1y , −2 cz ,− 2 cz ,0, 0,
0, 14 cy , − 21 cz ,11
0, − 21 cy , − 41 cz  ,0, − 4 c1y , −4 cz ,− 4 cz ,0, 0,
0, 14 cy , − 41 cz ,1
0, − 21 cy , 0 ,0, − 4 cy , 0 ,(0, 0, 0) , 1
0, 4 cy , 0 ,11
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,

1
1
c
,
c
x
 21 4 1y , − 2 c1z  21 cx , − 21 cy , − 41 cz  21 cx , − 4 c1y , −4 cz− 4 cz 21 cx , 0,

1
1
c
,
c
x
 21 4 1y , − 4cz
 21 cx , − 21 cy , 0
 21 cx , − 4 cy , 00 21 cx , 0,

1
c
,
c
x
 21 4 1y , 0 1  21 cx , − 21 cy , 41 cz  21 cx , − 41cy ,4 cz
c
4 z 21 cx , 0,1
1
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 に記述していく．波動関数の複素共役の関係を利用して計算を省いた点は，ウエイトのところを 2 にしておく．
####

CG 法 (P 方程式) の収束判定
CG 法 (K-S 方程式) の収束判定
SCF の収束判定
1.0D-12
1.0D-12
1.0D-8
CG 法 (P 方程式) 最低反復回数
CG 法 (P 方程式) 最大反復回数
CG 法 (P 方程式) の最出発:N (1 回/N 回)
1
400
500
CG 法 (KS 方程式) の最大反復回数 CG 法 (KS 方程式) のスイッチ (0;no 1;yes)
2
1
SCF における電荷密度の更新比
SCF の最大反復回数

19

#
#
#
#

#
#
#
#

0.10D0200
MD の回数:Ns(start)
MD の回数:Ne(end) (Ne-Ns が MD の回数になる)
MD の Time step (a.u.)
0.0D0
0
1
X 方向の電界 (a.u.)
Y 方向の電界 (a.u.)
Z 方向の電界 (a.u.)
0.0D0
0.0D0
0.0D0
k point (個数分縦に並べる)
X
Y
Z
ウエイト
-0.50D0-0.50D0-0.50D01-0.25D0-0.50D0-0.50D02
0.00D0-0.50D0-0.50D01-0.50D0-0.25D0-0.50D02-0.25D0-0.25D0-0.50D02
0.00D0-0.25D0-0.50D02
0.25D0-0.25D0-0.50D02-0.50D0
0.00D0-0.50D01-0.25D0
0.00D0-0.50D02
0.00D0
0.00D0-0.50D01-0.50D0-0.50D0-0.25D02-0.25D0-0.50D0-0.25D02
0.00D0-0.50D0-0.25D02
0.25D0-0.50D0-0.25D02-0.50D0-0.25D0-0.25D02-0.25D0-0.25D0-0.25D02
0.00D0-0.25D0-0.25D02
0.25D0-0.25D0-0.25D02-0.50D0
0.00D0-0.25D02-0.25D0
0.00D0-0.25D02
0.00D0
0.00D0-0.25D02
0.25D0
0.00D0-0.25D02-0.50D0
0.25D0-0.25D02-0.25D0
0.25D0-0.25D02
0.00D00.25D0-0.25D02
0.25D00.25D0-0.25D02-0.50D0-0.50D0
0.00D01-0.25D0-0.50D0
0.00D02
0.00D0-0.50D0
0.00D01-0.50D0-0.25D0
0.00D02-0.25D0-0.25D0
0.00D02
0.00D0-0.25D0
0.00D02
0.25D0-0.25D0
0.00D02-0.50D0
0.00D0
0.00D01-0.25D0
0.00D0
0.00D02
0.00D0
0.00D0
0.00D01
フェルミ温度 (KT) (a.u.)
固有値の最低値 (a.u.)
固有値の最高値 (a.u.)
0.0001D0-1000.0D0
1000.0D0
Rayleigh-Ritz:N (1 回/N scf 回) 固有値の並べ替え:N (1 回/N scf 回) 動かさない原子の potential の計算を省く (0;no 2;yes)110
全電子数-全陽子数0.0D0Force の cutoff
0.5D-2

他のファイルの作成 次に，パラメータファイルの作成を行う.5.2 章と異なるところは，9 行目の NUMK が 36 になることだけである. パラメータファイルの作成例を下に記す.
param.f の作成例 (初めの 24 行のみ)
1:!
********** param.f 05/02/2011-01 **********
2:
IMPLICIT REAL*8 (A-H,O-Z)
3:
INCLUDE ’mpif.h’
4:! number of processors (nodes)
5:
PARAMETER (NPROCX=1,NPROCY=1,NPROCZ=1)
6:! size of supercell
7:
PARAMETER (XMAX=3.32D0,YMAX=3.32D0,ZMAX=3.32D0)
8:! number of grid points
9:
PARAMETER (NXMAX=10,NYMAX=10,NZMAX=10)
10:! number of grid points in core region
11:
PARAMETER (NPXMAX=10,NPYMAX=10,NPZMAX=10)
12:! number of wave functions, number of wave atoms, number grid points for spiral boundary conditions
13:
PARAMETER (NEIGMX=8,NATOM=2,NRASXY=0,NRASXZ=0,NRASYZ=0)
14:! number of sampling k-pints, electron spin

20

15:
PARAMETER (NUMKX=1,NUMKY=1,NUMKZ=1,NUMK=36,NUMS=1)
16:! number of grid points for double grid technique
17:
PARAMETER (NMESH=3)
18:! switch for periodic boundary conditions
19:
PARAMETER (NPERI=3)
20:! charge of jellium (if you do not use jellium, set CHRJEL to be zero.)
21:
PARAMETER (ENDJEL=0.0D0,CHRJEL=0.0D0)
22:
23:! switch for the type of computer (scalar;0, vector;1)
24:
PARAMETER (MACHINE=0)

原子座標ファイル atom.xyz は 5.2 章と全く同じでよい．コンパイル，実行も 5.2 章と同じである．以下に電子状態計算結果 mdresult.dat の出力例を記す．
mdresult.dat の出力例
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.110666667D+00
0.110666667D+00
0.110666667D+00
nonlocal region [vnlx,vnly,vnlz (a.u.)]
0.664000000D+01
0.664000000D+01
0.664000000D+01
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 2
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 1
number 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.315756236D+02
0.100000000D-03
cufoff of the force (a.u.)
0.500000000D-02

21

mixing ratio of scf
0.100000000D+00
nrasy,nrasz
00
total charge
1.99999999999993
fermi level 7.149786577779170E-003 (hartree)
atomic charge
1 0.999997434689555
21.00000256531064
1
1
1 eigen value =0.1271181160D+00
1
1
2 eigen value =0.1271183920D+00
1
1
3 eigen value =0.1271198154D+00
1
1
4 eigen value =0.1271200480D+00
1
1
5 eigen value =0.1271217278D+00
1
1
6 eigen value =0.1271219741D+00
1
1
7 eigen value =0.2946099298D+00
1
1
8 eigen value =0.2946100753D+00

0.1083547447D-10
0.1713894363D-10
0.4268511246D-11
0.3571187958D-11
0.1533725279D-10
0.5014411299D-11
0.1855425777D-10
0.6763772503D-10

0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00
0.0000000000D+00····································
36
36
36
36
36
36
36
36
scf=

1
1 eigen value =-0.1175039808D+00
0.1286537827D-10
0.3125000000D-01
1
2 eigen value =0.1996909149D+00
0.1417640513D-10
0.0000000000D+00
1
3 eigen value =0.1996938022D+00
0.8841951083D-11
0.0000000000D+00
1
4 eigen value =0.1996971599D+00
0.7549029704D-11
0.0000000000D+00
1
5 eigen value =0.2923582235D+00
0.6749380587D-02
0.0000000000D+00
1
6 eigen value =0.5164111164D+00
0.9270942044D-10
0.0000000000D+00
1
7 eigen value =0.5164138899D+00
0.4742088553D-10
0.0000000000D+00
1
8 eigen value =0.5164167871D+00
0.6850669262D-11
0.0000000000D+00
68
dp=
0.9186551643D-08
one electoron energy-7.842581501942483E-002 (hartree)
exchange correlation energy -6.650449821272245E-003 (hartree)
hartreeenergy-9.416134072745397E-004 (hartree)ewald
energy-0.548077326733235(hartree)field correction
0.000000000000000E+000 (hartree)
helmholtz free energy-0.634140477571235(hartree)totalenergy-0.634095204981207(hartree)mdstep=1
atom coordinate (input)
1 -0.3320000000000000D+01 -0.3320000000000000D+01 -0.3320000000000000D+01
3
02
0.0000000000000000D+00
0.0000000000000000D+00
0.0000000000000000D+00
30
atom coordinate (output)
1 -0.3320000000000000D+01 -0.3320000000000000D+01 -0.3320000000000000D+01
3
02
0.0000000000000000D+00
0.0000000000000000D+00
0.0000000000000000D+00
30
atomic force
10.2852771071137328D-12 -0.1347960982550435D-12
0.1335356698440086D-12
2 -0.2852771071137320D-12
0.1347960982550435D-12 -0.1335356698440086D-12
real time470.199990034103(sec)com. time 0.109968900680542(sec)com. vol. 0.000000000000000E+000 (byte)======================================================================
This calculation was executed using pseudopotential database NCPS97.
For NCPS database, see K.Kobayashi, Comput. Mater. Sci. 14, 72 (1999).

A

0
0

0
0

12866.00
12866.00

0000

12866.00
12866.00

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

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

12

ソフト型ポテンシャルを用いた場合，C 原子については括弧付きの値が使用可能である．本マニュアルの例題
ではこちらを用いた．

22

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

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

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

B

擬ポテンシャルのカットオフ半径の決定方法

擬ポテンシャルのカットオフ半径の大きさは，セルに含まれる原子の中で，最も大きなカッ
トオフ半径を要求されるものよりもさらに大きな値を用いる．各原子に必要なカットオフ半径の大きさは以下のとおりである．
元素
H
Li
C
N
O
F
Na

カットオフ半径 (a.u.)
0.951
2.800
1.450
1.240
1.519
1.256
4.200

元素
Al
Si
S
Cl
Fe
Cu
W

カットオフ半径 (a.u.)
2.500
2.300
1.750
1.950
2.550
2.600
2.800

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

23