DCDFTBMDDivide-and-Conquer Density Functional Tight-Binding Molecular Dynamics

ユーザーマニュアルバージョン 2.0
2019年12月

Nakai Group
Department of Chemistry and Biochemistry,
School of Advanced Science and Engineering, Waseda University

概要
DCDFTBMD は、分割統治型密度汎関数強束縛 (DC-DFTB) 法に基づく量子分子動
力学 (MD) シミュレーションを実行するためのプログラムである。DC-DFTB 計
算は線形に近いスケーリングを達成し、大規模分子の取り扱いに適している。また、MPIとOpenMPを用いたハイブリッド並列実装により、高速な量子MDシミュレーションが可能である。さらに、本プログラムは、プロパティ計算、構造最
適化計算、および振動数計算の機能を持つ。

DCDFTBMDウェブサイト
開発背景、ダウンロード、マニュアル、および論文一覧を含むDCDFTBMDに関する一般的な情報は、以下のウェブサイトをご参照ください。
http://www.chem.waseda.ac.jp/dcdftbmd/

更新履歴
バージョン2.0での機能強化


メタダイナミクス計算

Quick-minあるいはFIRE法による構造最適化



LSQ法による初期Mulliken電荷予測

数値微分によるHessian計算の中断と再開

数値微分によるエネルギー勾配計算



拡張パラメータファイルとf電子の取り扱い

SCC計算におけるDFTB+U補正



SCC計算におけるonsite補正D3H5法のF原子への拡張

以下の機能は、MPI/OpenMP版のみ利用可能

MPI通信を用いたBroyden法

Orbitally resolved SCC計算の高速化

Orbitally resolved SCC計算におけるspline補間法

バージョン1.0からのデフォルト設定の変更


温度一定MD計算における温度の許容振れ幅
変更前: 100 K → 変更後: 10000 K

DFTB電荷非依存項の微分計算における数値微分刻み幅
変更前: 固定値 (0.01 atomic unit) → 変更後: オプション化

D3のHessian計算における数値微分刻み幅
変更前: 固定値 (0.0001 atomic unit) → 変更後: オプション化

Damped γ-function計算におけるH原子識別のための参照Hubbard値
変更前: 固定値 (0.4195 atomic unit) → 変更後: オプション化

以下の変更は、serial, OpenMP版のみ適用


中間データに対するディスクI/Oの使用の無効化

バージョン1.0からの不具合修正


DFTB3計算における原子エネルギー分割

開殻系における全原子中性とみなす初期電荷推測

以下の修正は、MPI/OpenMP版のみ適用


Orbitally resolved SCC計算におけるHubbard derivativeパラメータ値の確認

目次
1 DCDFTBMDの入力

1

1.1 Keywordセクション2

1.1.1 概要2

1.1.2 SCC

3

基本設定

3

電荷密度混合3

初期電荷推測4

固有値問題ソルバー

4

3次補正

5

修正γ関数

5

開殻系とスピン分極補正

6

外場効果

6

DFTB+U補正

7

Onsite補正

7

電子プロパティ評価

7

MPI/OpenMP版のみの追加オプション

8

1.1.3 NCC

10固有値問題ソルバー

10

開殻系

10

電子プロパティ評価

10

1.1.4 DC

11

基本設定

11

詳細オプション

11

MPI/OpenMP版のみの追加オプション

12

1.1.5 DISP/DISPERSION

14

基本設定

14

Grimme’s DFT-D2の詳細オプション

15i

Slater-Kirkwoodの詳細オプション

15

Lennard-Jonesの詳細オプション

15

Grimme’s DFT-D3の詳細オプション

15Řezáč’s D3H4の詳細オプション

16

DFT-ulgの詳細オプション

17

dDMCの詳細オプション

17Řezáč’s D3H5の詳細オプション

17

Kubillus’s D3Xの詳細オプション

18

振動数計算に対する詳細オプション

19

MPI/OpenMP版のみの追加オプション

19

1.1.6 PBC

21

基本設定

21

Ewald和の詳細オプション

21

多重極展開の詳細オプション

22

MPI/OpenMP版のみの追加オプション

22

1.1.7 OPT/OPTIMIZE

23

基本設定

23

制約付き最適化24

格子ベクトル最適化

24

SCC収束の高速化

25

BFGSの詳細オプション

25

Steepest descentの詳細オプション

25

Conjugate gradientの詳細オプション

26

Quick-minの詳細オプション

26

FIREの詳細オプション

26

1.1.8 FREQ/FREQUENCY

27

基本設定

27

詳細オプション

28

MPI/OpenMP版のみの追加オプション

29

1.1.9 MD

30
ii

基本設定

30

エネルギー一定アンサンブル30

温度一定アンサンブル31

温度制御の詳細オプション32

圧力一定アンサンブル32

圧力制御の詳細オプション

33

乱数生成

33

ソフトポテンシャル

34

RATTLE法を用いた拘束動力学

36

SCC収束の高速化

36

アニーリング37

メタダイナミクス

37

メタダイナミクスの詳細オプション

38

1.1.10 MISC

40

追加解析とプロパティ評価

40

グローバルオプション

41

Serial, OpenMP版のみの追加オプション

42

1.2 Titleセクション43

1.3 Parameterセクション44

1.3.1 概要

44

1.3.2 定数パラメータの追加指定

45

1.3.3 定数パラメータ追加指定ファイルの形式

47

スピン定数ファイル

47

U–Jパラメータファイル

47

Onsiteパラメータファイル

48

Slater-Kirkwood分散力補正用パラメータファイル

48

Lennard-Jones分散力補正用パラメータファイル

49

1.4 Geometryセクション50
iii

1.4.1 概要

50

1.4.2 追加情報の指定

52

1.4.3 DC-DFTBにおける部分系の手動指定

53

1.5 補助インプットファイル

54

1.5.1 バイナリ形式ファイル

54

1.5.2 ASCII形式ファイル

56

SCC計算における初期電荷の読み込み

56

振動数計算におけるHessian行列の読み込み

58

MD計算における初期速度の読み込み

58

RATTLE法における結合距離拘束の手動指定

59

アニーリングにおける温度スケジュールの設定

60

メタダイナミクスにおけるCVの設定

61

メタダイナミクスにおけるバイアスポテンシャルの読み込み
メタダイナミクスにおけるポテンシャル障壁の設定

74
75

2 DCDFTBMDの実行

76

3 DCDFTBMDの出力

77

3.1 標準アウトプットファイル

77

3.2 詳細データファイル

80

3.3 その他の出力

85

3.3.1 DFTB計算における出力85

3.3.2 振動数計算における出力

86

3.3.3 MD計算における出力87

3.3.4 メタダイナミクスにおける出力

89

参考文献

91
iv

目次

1 DCDFTBMDの入力
DCDFTBMDの実行には、インプットファイルとDFTBのパラメータ (Slater-Kosterファイル) が必要である。インプットファイルは、 “dftb.inp” をデフォルトのファイル名とし、主として以下の4つのセクションから構成される。
Keywordセクション: 計算条件の指定 (1.1節)



Titleセクション: 計算についてのコメントを記述 (1.2節)

Parameterセクション: DFTBパラメータの指定 (1.3節)



Geometryセクション: 分子構造の指定 (1.4節)

各セクションの終わりには空行が挿入される (本マニュアルでは、記号で示す) 。以下は、水分子に対するSCC-DFTBエネルギー一点計算のインプットファ
イルの例である。
Keywordセクション

SCC=TRUE˽DC=FALSE

Titleセクション

TITLE
Parameterセクション2
O˽2
oo.spl˽oh.spl
H˽1
ho.spl˽hh.spl


Geometryセクション3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
H0˽−1.01910555˽−0.08185946˽−0.00000000
H0˽−0.08185946˽−1.01910555˽−0.00000000


DC計算で系を手動分割する場合は、Geometryセクションの後に追加の情報を入力する (1.4.3節を参照) 。空白領域 (本マニュアルでは、記号˽で示す) は、半角スペースで入力する必要がある (タブは不可) 。

1目次

1.1 Keywordセクション
1.1.1 概要
SCC・DCキーワードの有効化・

SCC=TRUE˽DC=FALSE

無効化

この例では、1行に2つのキーワードが指定される。以下のように、2行に分けた指定も可能である。
SCC=TRUESCCキーワードの有効化

DC=FALSE

DCキーワードの無効化

現在利用可能なキーワードは以降の節で説明される。いずれのキーワードも未指定の場合、SCCとDCのキーワードが有効化される。その他のキーワードは以下の形式で有効化できる。
keyword=TRUE

この指定方法では、波かっこ内に記述されたデフォルト値が各オプションに用いられる。いくつかのキーワードは、計算条件を指定するためのオプションをもつ。丸かっこ内の斜体での記述は、オプションのデータ型を表す。オプションは、
以下の形式で指定できる。
keyword=(option1=value1˽option2=value2˽…)

option1・option2を1行で指定

keyword=(option1=value1)

option1・option2を2行に分けて

keyword=(option2=value2˽…)

指定

オプションが赤字もしくは青字で記述されている場合、そのオプションはserial,
OpenMP版もしくはMPI/OpenMP版でのみ機能する。

2

目次

1.1.2 SCC
SCCキーワードは、SCC計算 (Mulliken電荷を自己無撞着に決定) するか否かを指定する。既に記述した通り、デフォルトはSCC=TRUEである。すなわち、 SCCDFTB (DFTB2) [1] 計算が実行される。

基本設定

MAXITER (integer) {200}
SCCサイクルの最大数



ENERGYCONV/ECONV (real) {1.0e−09}
エネルギー変化の収束条件 (Hartree)

DENSITYCONV/DCONV (real) {1.0e−06}
密度変化の収束条件 (atomic unit)

電荷密度混合

MIXER (integer) {1}
電荷密度混合法
解析的Hessian計算におけるcoupled perturbed DFTB計算ではBroyden法
を使用







1: Modified Broyden [2]2: Simple mixing
3: Modified Anderson [3]4: DIIS [4]

ALMIX (real) {2.0e−01}
電荷密度混合割合0より大きく1以下の値を指定



BROYITER/MAXBROYDEN (integer) {70}
Broyden法による電荷密度混合の最大数

NGEN (integer) {4}Anderson, DIIS法を開始するサイクル数
この値に到達するまではsimple mixingを使用
3

目次

初期電荷推測READCHARGE (method) {NONE}
電荷密度ファイルの読み込み
 NONE: 読み込みなし (価電子数から電荷密度を決定)
 BINARY: 計算時に生成される ”chrgfile” ファイルの使用
ファイルが見つからない場合は価電子数から電荷密度を決定
 ASCII: “charge.dat” ファイルによる手動指定 (形式は1.5.2節を参照)
ファイルが見つからない場合は価電子数から電荷密度を決定



ZEROCHARGE (logical) {FALSE}
SCC計算の初期電荷を常に中性原子の値に指定
ZEROCHARGE=FALSEの場合、前ステップで収束した電荷を初期電荷に使用固有値問題ソルバー
SOLVER (integer) {1}
Hamiltonian行列に対する固有値と固有ベクトルの計算 (対角化) 法 1: DSYGVD
 2: DSYGVPSEUDODIAG/PDIAG (logical) {FALSE}
擬対角化 [5]SCC収束後に通常の対角化を一度実行

PDIAGTHRESH (real) {1.0e−02}
擬対角化を開始する密度変化の閾値PDIAGKAPPA (real) {4.0e−02}
擬対角化におけるパラメータ

4

目次

3次補正


THIRDDIAG (logical) {FALSE}SCC Hamiltonianに対する3次補正 (対角成分のみ) [6]
Hubbard derivativeをParameterセクションに指定 (1.3.2節を参照)



THIRDFULL (logical) {FALSE}SCC Hamiltonianに対する3次補正 [7]Hubbard derivativeをParameterセクションに指定 (1.3.2節を参照)
DFTB3モデル [7] では、DAMPXH=TRUEを追加指定

修正γ関数


DAMPXH (logical) {FALSE}
X–Hペアに対するSCC相互作用の短距離でのダンピング [6–8]

DAMPXHZETA (real) {4.0e+00}
DAMPXH=TRUEでのダンピング関数を制御するパラメータ

DAMPXHREFU (real) {4.195e−01}
DAMPXH=TRUEにおけるH原子識別のための参照Hubbard値



UDERIVGAUSS (logical) {FALSE}THIRDDIAG=TRUEにおけるガウス関数によるHubbard derivativeへの電荷依存項の導入 [9]

UDERIVGAUSSD0 (real) {−9.0e−02}
UDERIVGAUSS=TRUEにおけるガウス関数中のパラメータ (高さ)

UDERIVGAUSSG0 (real) {1.61e+01}
UDERIVGAUSS=TRUEにおけるガウス関数中のパラメータ (幅)UDERIVGAUSSQ0 (real) {7.5e−01}
UDERIVGAUSS=TRUEにおけるガウス関数中のパラメータ (中心位置)

5

目次



GAMMASPLINE (logical) {FALSE}
3次spline補間を用いたSCC相互作用の短距離項の計算 (詳細は参考文献[10] を参照)
MPI/OpenMP版では、RUNTYPE=SHELL0との組み合わせ不可GAMMASPLINEGRID (real) {4.0e−02}
GAMMASPLINE=TRUEにおけるグリッド点の間隔 (atomic unit)

開殻系とスピン分極補正


SPIN (logical) {FALSE}
スピン分極補正の適用 [11, 12]
スピン定数を含むファイルをParameterセクションに指定 (1.3.2, 1.3.3節
を参照)

SPINRELAX (logical) {FALSE}
スピン分極における全スピンの最適化



NELECTHRESH (real) {1.0e−06}
SPINRELAX=TRUEにおける全電子数保存の閾値

UNPAIREDELEC (real) {0.0e+00}
不対電子数 (正の実数で指定可能)
正の値の場合はGeometryセクションで指定されたスピン多重度を棄却(1.4.1節を参照)



RESTRICT (logical) {FALSE}
制限開殻アプローチによる二重項の計算
RESTRICT=FALSEは非制限開殻アプローチによる計算

外場効果
POINTCHARGE (logical) {FALSE}
点電荷の導入点電荷の座標と値は原子座標の後に “PC” から始まる行で指定 (1.4.1節
を参照)

6

目次



NUMPC (integer) {0}点電荷の導入数

EFIELD (logical) {FALSE}
外部静電場の導入
クラスターモデルのみ計算可能
電場ベクトルは原子座標の後に “EF” から始まる行で指定 (1.4.1節を参照)DFTB+U補正

PLUSU (logical) {FALSE}DFTB+U補正の適用 [13]U–Jパラメータを含むファイルをParameterセクションに指定 (1.3.2,
1.3.3節を参照)



PLUSUTYPE (integer) {1}
DFTB+U補正の方法
 1: Fully localized limit (FLL) 2: Pseudo self-interaction correction (pSIC)

Onsite補正

ONSITE (logical) {FALSE}
Onsite補正の適用 [14]Onsiteパラメータを含むファイルをParameterセクションに指定 (1.3.2,
1.3.3節を参照)



STOREBLOCKPOP (logical) {TRUE}
PLUSU=TRUE, ONSITE=TRUEで収束した電荷密度をメモリに格納

電子プロパティ評価


MAYER (logical) {FALSE}
Mayer結合次数の解析 [15]

7

目次



CM3 (logical) {FALSE}
CM3電荷の計算 [16, 17]



ORBITALPOP/ORBPOP (logical) {FALSE}
Mulliken密度解析の結果を各原子minimal basis毎に出力



POLAR (logical) {FALSE}
外部静電場を用いた双極子モーメントと静的分極率計算
クラスターモデルのみ計算可能



POLEFSTR (real) {2.0e−03}
POLAR=TRUEにおける外部静電場の強さ

MPI/OpenMP版のみの追加オプション


RUNTYPE (method) {ATOM}DFTB計算の解き方のベース
Serial, OpenMP版はRUNTYPE=ATOMに相当
 SHELL0: 原子価軌道
THIRDFULL=TRUE, THIRDDIAG=TRUEの指定不可
 SHELL: 原子価軌道
THIRDDIAG=TRUEの指定不可
 ATOM: 原子
ORSCC=TRUEの指定不可

ORSCC (logical) {FALSE}
異なる軌道角運動量に対して異なるHubbardパラメータを使用 [12, 18]
RUNTYPE=SHELLのみ指定可

BROYPARALLEL (logical) {FALSE}
Broyden法に対するMPI通信の使用



MIXERSAVEMEMORY (logical) {FALSE}RUNTYPE=ATOMかつ閉殻系の計算において電荷密度混合計算時のメ
モリ使用量を削減

8

目次



GSSAVEMEMORY (logical) {FALSE}GAMMASPLINE=TRUEかつ周期系の計算においてSCC相互作用の短距
離項計算時のメモリ使用量を削減



GSDOMAIN (logical) {FALSE}GAMMASPLINE=TRUEにおける連結リスト [19, 20] を用いたSCC相互作用の短距離項の計算



GSDOMAINBLOCKXYZ/GSDOMAINBLOCK (integer) {5}
GSDOMAIN=TRUEで系をセル分割する際のx, y, z方向のブロック数



GSDOMAINBLOCKX (integer) {5}
GSDOMAIN=TRUEで系をセル分割する際のx方向のブロック数



GSDOMAINBLOCKY (integer) {5}
GSDOMAIN=TRUEで系をセル分割する際のy方向のブロック数



GSDOMAINBLOCKZ (integer) {5}
GSDOMAIN=TRUEで系をセル分割する際のz方向のブロック数



GSDOMAINNEIGHBOR (integer) {1}
GSDOMAIN=TRUEでリンクするセルの隣接参照数9

目次

1.1.3 NCC
NCCキーワードは、non-SCC (NCC-DFTBあるいはDFTB1) [21, 22] 計算するか否かを指定する。NCC=TRUEはSCC=FALSEによる指定も可能である。

固有値問題ソルバー
SOLVER (integer) {1}
Hamiltonian行列に対する固有値と固有ベクトルの計算 (対角化) 法 1: DSYGVD
 2: DSYGV

開殻系
RESTRICT (logical) {FALSE}
制限開殻アプローチによる二重項の計算
RESTRICT=FALSEは非制限開殻アプローチによる計算

電子プロパティ評価


PRINTMULLIKEN (logical) {FALSE}
Mulliken電荷密度解析の実行



MAYER (logical) {FALSE}
Mayer結合次数の解析 [15]PRINTMULLIKEN=TRUEの指定が必要

CM3 (logical) {FALSE}
CM3電荷の計算 [16, 17]
PRINTMULLIKEN=TRUEの指定が必要
DFTB2で用いられるパラメータを計算に使用



ORBITALPOP/ORBPOP (logical) {FALSE}
Mulliken密度解析の結果を各原子minimal basis毎に出力
PRINTMULLIKEN=TRUEの指定が必要

10

目次

1.1.4 DC
DCキーワードは、分割統治電子状態計算 [23, 24] を適用するか否かを指定する。
SCCキーワードと同様、このキーワードは未指定の場合有効化される。DC-DFTB
計算を無効化したい場合、DC=FALSEと指定する必要がある。

基本設定


BETA (real) {8.0e+02}
Fermi関数における温度の逆数のパラメータ (atomic unit)

BUFRAD/BUFFERRADIUS (real) {5.0e+00}
球状バッファ領域の半径 (Å)



SUBTYPE (method) {AUTO}
部分系の構築法 ATOM: 1原子 SPECIFY: 手動指定 (1.4.3節を参照)
 SYSTEM: 全原子
 AUTO: 格子状空間に分割し自動決定
デフォルトは一辺5.0 Åの立方体領域
 AUTOXH: SUBTYPE=AUTOで自動決定する部分系においてH原子を含む共有結合を切断しないように調整
 SEMIAUTO: 一部の部分系を手動指定し、残りを自動決定 (1.4.3節を参照)
 SEMIAUTOXH: SUBTYPE=SEMIAUTOで自動決定する部分系においてH原子を含む共有結合を切断しないように調整

詳細オプション


DELTARXYZ/DELTAR (real) {5.0e+00}
SUBTYPE=AUTOで分割する際のグリッド一辺の長さ (Å)



DELTARX (real) {5.0e+00}
SUBTYPE=AUTOで分割する際のグリッドのx成分の長さ (Å)

11

目次



DELTARY (real) {5.0e+00}
SUBTYPE=AUTOで分割する際のグリッドのy成分の長さ (Å)



DELTARZ (real) {5.0e+00}
SUBTYPE=AUTOで分割する際のグリッドのz成分の長さ (Å)



SUBXHSCALE (real) {1.1e+00}
SUBTYPE=AUTOXHとSUBTYPE=SEMIAUTOXHでH原子を含む共有結合を判定する際の共有結合半径 [25] の和のスケーリングファクターTRANSCOORD (logical) {TRUE}SUBTYPE=AUTO,SUBTYPE=SEMIAUTO,SUBTYPE=AUTOXH,SUBTYPE=SEMIAUTOXHで系を自動分割する際の座標系
 TRUE: 標準配向に変換
 FALSE: 入力構造をそのまま使用
NOFERMI (logical) {FALSE}
共通のFermi準位計算を省略

MPI/OpenMP版のみの追加オプション

FERMITYPE (method) {PARALLEL}
Fermi準位計算の方法 (詳細は参考文献 [26] を参照)
 PARALLEL: CRI (compute-reduce-interpolate) アルゴリズム
 SERIAL: GSS (gather-sort-solve) アルゴリズム



DOMAIN (logical) {FALSE}
連結リスト [19, 20] を用いた部分系の局在化領域への割り当て



DOMAINBLOCKXYZ/DOMAINBLOCK (integer) {5}
DOMAIN=TRUEで系をセル分割する際のx, y, z方向のブロック数

DOMAINBLOCKX (integer) {5}
DOMAIN=TRUEで系をセル分割する際のx方向のブロック数

DOMAINBLOCKY (integer) {5}
DOMAIN=TRUEで系をセル分割する際のy方向のブロック数

12

目次



DOMAINBLOCKZ (integer) {5}
DOMAIN=TRUEで系をセル分割する際のz方向のブロック数

DOMAININDEXTYPE (method) {ATOM}DOMAIN=TRUEでリンクする粒子の定義
 SUBSYSTEM: 部分系の全構成原子が同一領域に存在する時にインデックスを付与 ATOM: 全原子についてインデックスを付与



DOMAINNEIGHBOR (integer) {1}
DOMAIN=TRUEでリンクするセルの隣接参照数



SINGLEPREC (logical) {FALSE}
局在化領域決定における単精度浮動小数点の使用



XHDOMAIN (logical) {FALSE}SUBTYPE=AUTOXHでH原子を含む共有結合の判定に連結リスト [19,20] を使用XHDOMAINBLOCKXYZ/XHDOMAINBLOCK (integer) {5}
XHDOMAIN=TRUEで系をセル分割する際のx, y, z方向のブロック数

XHDOMAINBLOCKX (integer) {5}
XHDOMAIN=TRUEで系をセル分割する際のx方向のブロック数

XHDOMAINBLOCKY (integer) {5}
XHDOMAIN=TRUEで系をセル分割する際のy方向のブロック数

XHDOMAINBLOCKZ (integer) {5}
XHDOMAIN=TRUEで系をセル分割する際のz方向のブロック数

XHDOMAINNEIGHBOR (integer) {1}
XHDOMAIN=TRUEでリンクするセルの隣接参照数

13

目次

1.1.5 DISP/DISPERSIONDISP/DISPERSIONキーワードは、 非共有結合性相互作用 (分散力・水素結合・ハロゲン結合) を経験的に補正するか否かを指定する。

基本設定


DISPTYPE (integer) {1}非共有結合性相互作用補正法
 1: Grimme’s DFT-D2 [27]
 2: Slater-Kirkwood [28]計算に必要なパラメータを含むファイルをParameterセクションに
指定 (1.3.2, 1.3.3節を参照)
 3: Lennard-Jones [29]デフォルトのLennard-JonesパラメータはUFF力場 [30] の値
READUFF=TRUEによる手動指定も可能
 4: Grimme’s DFT-D3 (zero-damping) [31]
 5: Grimme’s DFT-D3(BJ) (Becke and Johnson damping) [31, 32]
 6: Řezáč’s D3H4 [33]水 素 結 合 を 経 験 的 に 補 正 す る た め 、 SCC キ ー ワ ー ド に お い て
DAMPXH=FALSEの指定を想定 (例えば、参考文献 [34] を参照)
 7: DFT-ulg [35]
 8: dDMC [36]
 9: Řezáč’s D3H5 [37]H5補正はもとのSCC相互作用項を変更するため、SCCキーワードに
おいてDAMPXH=FALSEの指定を想定
NCC=TRUE (あるいはSCC=FALSE) との組み合わせ不可
 10: Kubillus’s D3X [38]

14

目次

Grimme’s DFT-D2 (DISPTYPE=1)の詳細オプション
S6 (real) {1.0e+00}
Grimme’s DFT-D2におけるスケーリングファクター

Slater-Kirkwood (DISPTYPE=2)の詳細オプション


UPDATESKPARAMS (logical) {FALSE}
Slater-Kirkwood分散力パラメータをDFTB計算毎に更新
UPDATESKPARAMS=FALSEの場合、初期構造から得られる値を常に使用

Lennard-Jones (DISPTYPE=3)の詳細オプション
READUFF (logical) {FALSE}
Lennard-Jonesパラメータを手動指定 (1.3.2, 1.3.3節を参照)

Grimme’s DFT-D3 (DISPTYPE=4, 5)の詳細オプション


ZEROS8 (real) {6.73e−01}
Grimme’s DFT-D3におけるr−8項のスケーリングファクターZERORS6 (real) {1.235e+00}
Grimme’s DFT-D3における減衰関数中のパラメータ

BJS8 (real) {5.883e−01}
Grimme’s DFT-D3(BJ)におけるr−8項のスケーリングファクター



BJA1 (real) {5.719e−01}
Grimme’s DFT-D3(BJ)におけるパラメータ

BJA2 (real) {3.6017e+00}
Grimme’s DFT-D3(BJ)におけるパラメータ

THREEBODY (logical) {FALSE}
Grimme’s DFT-D3とDFT-D3(BJ)における3体補正項の取り込み

15

目次

Řezáč’s D3H4 (DISPTYPE=6)の詳細オプション


D3H4S6 (real) {1.0e+00}
Řezáč’s D3H4におけるD3分散力に対するパラメータ

D3H4RS6 (real) {1.215e+00}
Řezáč’s D3H4におけるD3分散力に対するパラメータ

D3H4ALPHA (real) {3.0e+01}Řezáč’s D3H4におけるD3分散力に対するパラメータ

D3H4SHH (real) {3.0e−01}Řezáč’s D3H4におけるH–H反発に対するパラメータ

D3H4EHH (real) {1.431e+01}Řezáč’s D3H4におけるH–H反発に対するパラメータ

D3H4R0HH (real) {2.35e+00}
Řezáč’s D3H4におけるH–H反発に対するパラメータ

D3H4CNN (real) {2.01e+00}
Řezáč’s D3H4におけるH4補正に対するパラメータ

D3H4CNO (real) {8.0e−01}Řezáč’s D3H4におけるH4補正に対するパラメータ

D3H4CON (real) {2.58e+00}
Řezáč’s D3H4におけるH4補正に対するパラメータ

D3H4COO (real) {1.11e+00}
Řezáč’s D3H4におけるH4補正に対するパラメータ

D3H4CWAT (real) {1.32e+00}
Řezáč’s D3H4におけるH4補正に対するパラメータ

D3H4CSNHR3 (real) {2.33e+00}
Řezáč’s D3H4におけるH4補正に対するパラメータ

D3H4CSCOO (real) {1.22e+00}
Řezáč’s D3H4におけるH4補正に対するパラメータ16

目次

DFT-ulg (DISPTYPE=7)の詳細オプション
ULGS (real) {7.012e−01}
DFT-ulgにおけるパラメータ
デフォルトはPBE汎関数に対する値

ULGB (real) {6.966e−01}
DFT-ulgにおけるパラメータ
デフォルトはPBE汎関数に対する値

dDMC (DISPTYPE=8)の詳細オプション


DDMCA (real) {1.018e+00}
dDMCにおけるパラメータDDMCB0 (real) {1.857e+00}
dDMCにおけるパラメータ



DDMCS (real) {4.6e+01}
dDMCにおけるパラメータŘezáč’s D3H5 (DISPTYPE=9)の詳細オプション


D3H5S6 (real) {1.0e+00}
Řezáč’s D3H5におけるD3分散力に対するパラメータ

D3H5S8 (real) {4.9e−01}Řezáč’s D3H5におけるD3分散力に対するパラメータ

D3H5RS6 (real) {1.25e+00}
Řezáč’s D3H5におけるD3分散力に対するパラメータ

D3H5ALPHA (real) {2.961e+01}Řezáč’s D3H5におけるD3分散力に対するパラメータ

D3H5SHH (real) {3.0e−01}Řezáč’s D3H5におけるH–H反発に対するパラメータ

17

目次



D3H5EHH (real) {1.431e+01}Řezáč’s D3H5におけるH–H反発に対するパラメータ

D3H5R0HH (real) {2.35e+01}Řezáč’s D3H5におけるH–H反発に対するパラメータ

D3H5SR (real) {7.14e−01}Řezáč’s D3H5におけるH5補正に対するパラメータ

D3H5SW (real) {2.5e−01}Řezáč’s D3H5におけるH5補正に対するパラメータ

D3H5KNH (real) {1.8e−01}Řezáč’s D3H5におけるH5補正に対するパラメータ

D3H5KOH (real) {6.0e−02}Řezáč’s D3H5におけるH5補正に対するパラメータ

D3H5KFH (real) {6.0e−02}Řezáč’s D3H5におけるH5補正に対するパラメータ

D3H5KSH (real) {2.1e−01}Řezáč’s D3H5におけるH5補正に対するパラメータ

Kubillus’s D3X (DISPTYPE=10)の詳細オプション


D3XS8 (real) {3.209e+00}
Kubillus’s D3XにおけるD3分散力に対するパラメータ

D3XA1 (real) {7.46e−01}
Kubillus’s D3XにおけるD3分散力に対するパラメータ

D3XA2 (real) {4.191e+00}
Kubillus’s D3XにおけるD3分散力に対するパラメータ

D3XC1 (real) {7.761e+00}
Kubillus’s D3XにおけるX補正に対するパラメータ

D3XC2 (real) {5.0e−02}
Kubillus’s D3XにおけるX補正に対するパラメータ
18

目次



D3XC3 (real) {4.518e+00}
Kubillus’s D3XにおけるX補正に対するパラメータ



D3XDNCL (real) {1.526e+00}
Kubillus’s D3XにおけるX補正に対するパラメータ



D3XDNBR (real) {1.349e+00}
Kubillus’s D3XにおけるX補正に対するパラメータ



D3XDNI (real) {1.521e+00}
Kubillus’s D3XにおけるX補正に対するパラメータ

D3XDOCL (real) {1.237e+00}
Kubillus’s D3XにおけるX補正に対するパラメータ

D3XDOBR (real) {1.099e+00}
Kubillus’s D3XにおけるX補正に対するパラメータ

D3XDOI (real) {1.313e+00}
Kubillus’s D3XにおけるX補正に対するパラメータ

振動数計算に対する詳細オプション


D3STEPSIZE (real) {1.0e−04}
解 析 的 Hessian 計 算(FREQ/FREQUENCY キ ー ワ ー ド に お け る

FREQTYPE=1) でGrimme’s DFT-D3 (DISPTYPE=4–6, 9, 10) の解析的1次微分を用いたHessian計算における数値微分刻み幅 (atomic unit)
MPI/OpenMP版のみの追加オプション


H4IDOMAIN (logical) {FALSE}Řezáč’s D3H4 (DISPTYPE=6) のH4補正計算でH原子の索引に連結リス
ト [19, 20] を使用



H4IDOMAINBLOCKXYZ/H4IDOMAINBLOCK (integer) {5}
H4IDOMAIN=TRUEで系をセル分割する際のx, y, z方向のブロック数



H4IDOMAINBLOCKX (integer) {5}
H4IDOMAIN=TRUEで系をセル分割する際のx方向のブロック数
19

目次



H4IDOMAINBLOCKY (integer) {5}
H4IDOMAIN=TRUEで系をセル分割する際のy方向のブロック数



H4IDOMAINBLOCKZ (integer) {5}
H4IDOMAIN=TRUEで系をセル分割する際のz方向のブロック数



H4IDOMAINNEIGHBOR (integer) {1}
H4IDOMAIN=TRUEでリンクするセルの隣接参照数



H4VDOMAIN (logical) {FALSE}Řezáč’s D3H4 (DISPTYPE=6) のH4補正計算で原子価の計算に連結リス
ト [19, 20] を使用



H4VDOMAINBLOCKXYZ/H4VDOMAINBLOCK (integer) {5}H4VDOMAIN=TRUEで系をセル分割する際のx, y, z方向のブロック数



H4VDOMAINBLOCKX (integer) {5}H4VDOMAIN=TRUEで系をセル分割する際のx方向のブロック数



H4VDOMAINBLOCKY (integer) {5}H4VDOMAIN=TRUEで系をセル分割する際のy方向のブロック数



H4VDOMAINBLOCKZ (integer) {5}H4VDOMAIN=TRUEで系をセル分割する際のz方向のブロック数



H4VDOMAINNEIGHBOR (integer) {1}
H4VDOMAIN=TRUEでリンクするセルの隣接参照数

20

目次

1.1.6 PBC
PBCキーワードは、周期境界条件に関する計算条件を指定する。Geometryセクシ
ョンで3つの並進ベクトルが検出された場合 (1.4.1節を参照) 、系は周期系とみなされ、PBC=TRUEの明記がなくても本キーワードは有効化される。

基本設定
STRESS (logical) {FALSE}
応力テンソルと格子ベクトルの微分計算



COULOMBTYPE (integer) {1}
Coulomb相互作用の計算法
 1: Ewald和 2: 多重極展開 [39] (詳細は参考文献 [10] を参照)
以下の制限が存在 応力テンソル計算は不可 (STRESS=FALSE)
 点電荷ありの計算は不可(SCC キ ー ワ ー ド に お い て

POINTCHARGE=FALSE) 全エネルギーの原子分割解析は不可 (MISCキーワードにおいて
PRINTATOME=FALSE) MPI/OpenMP版で原子価軌道ベースの計算は不可 (SCCキーワ
ードおいてRUNTYPE=ATOM) 立方体セルのみ動作確認
Ewald和 (COULOMBTYPE=1)の詳細オプション


EWALDALPHA (real) {0.0e+00}
Ewald和における実空間と逆格子空間の計算割合0より大きく1より小さい場合、指定した値を (一点計算でのみ) 使用
それ以外の場合、自動決定される値を使用 (デフォルトかつ推奨)

MAXGAMMA (integer) {7}
Ewald和におけるユニットセルの繰り返し回数

21

目次

多重極展開 (COULOMBTYPE=2)の詳細オプション
FMMP (integer) {8}
多重極展開における展開次数

FMMN (integer) {3}
多重極展開におけるグループ化されるセル数

FMMK (integer) {5}
多重極展開における階層数

MPI/OpenMP版のみの追加オプション


EWALDSAVEMEMORY (logical) {FALSE}
Ewald和逆格子空間計算時のメモリ使用量を削減
計算コストは増加



EWALDMEMORYGB (real) {8.0e+00}
Ewald和逆格子空間計算時に使用する一時配列の最大メモリ使用量 (GB)指定した値以上のメモリ使用量が見込まれる場合は、
EWALDSAVEMEMORY=FALSEの指定によらずメモリ使用量を削減する方法で計算を実行

FMMSAVEMEMORY (logical) {FALSE}
多重極展開体球調和関数計算時のメモリ使用量を削減
計算コストは増加

22

目次

1.1.7 OPT/OPTIMIZE
OPT/OPTIMIZEキーワードは、構造最適化計算を実行するか否かを指定する。

基本設定
MAXITER (integer) {50}
構造最適化サイクルの最大数



GRADCONV/GCONV (real) {1.0e−04}
エネルギー勾配変化の収束条件 (atomic unit)

OPTTYPE (integer) {1}
構造最適化の方法








1: Broyden-Fletcher-Goldfarb-Shanno (BFGS)
2: Steepest descent
3: Conjugate gradient
4: Quick-min [40]5: Fast inertial relaxation engine (FIRE) [41]MAXMOVE (real) {3.0e−01}
OPTTYPE=2–5における原子の最大移動サイズ (Å)

PRINTCOORD (logical) {TRUE}
構造最適化各ステップの原子座標を標準アウトプットファイルに出力
PRINTCOORD=FALSEの場合、最終ステップの原子座標のみ標準アウ
トプットファイルに出力

PRINTGRAD (logical) {TRUE}
構造最適化各ステップのエネルギー勾配を標準アウトプットファイルに出力
PRINTGRAD=FALSEの場合、各ステップの力の最大成分と根二乗平均
値のみ標準アウトプットファイルに出力

23

目次

制約付き最適化CONSTRAINT (logical) {FALSE}
x, y, z成分に対する制約つき構造最適化 (指定方法は1.4.2節を参照)

格子ベクトル最適化


LATTICEOPT (logical) {FALSE}
格子ベクトルの最適化
格子ベクトルの最適化 (CONSTRAINT=TRUE) および点電荷の使用(SCCキーワードにおけるPOINTCHARGE=TRUE in SCC keyword) との組み合わせ不可OPTTYPE=1 (BFGS) の場合、conjugate gradientのルーチンを使用

EXTPRESSURE (real) {0.0e+00}
格子ベクトル最適化における外圧 (atomic unit)
値が0以外の場合、系のエンタルピー (E + pV) を追加出力

ISOSCALING (logical) {FALSE}
等方的なセルサイズ変化による格子ベクトル最適化
FIXANGLE=TRUEと同時指定不可



FIXANGLE (logical) {FALSE}
セルの角度を固定した状態での格子ベクトル最適化
ISOSCALING=TRUEと同時指定不可



FIXANGLE1 (logical) {FALSE}FIXANGLE=TRUEにおいて1つ目の格子ベクトルを固定



FIXANGLE2 (logical) {FALSE}FIXANGLE=TRUEにおいて2つ目の格子ベクトルを固定



FIXANGLE3 (logical) {FALSE}FIXANGLE=TRUEにおいて3つ目の格子ベクトルを固定

24

目次



LATTICEONLY (logical) {FALSE}
原子座標は固定し格子ベクトルのみ最適化
全ての原子座標後に記号*を指定することでも指定可能 (1.4.2節を参照)

SCC収束の高速化


LSMULLIKEN (logical) {FALSE}
原子座標の最小二乗フィッティングによるMulliken電荷初期値の予測[42, 43]
DFTB+U補正 (SCCキーワードにおけるPLUSU=TRUE) およびonsite補
正 (SCCキーワードにおけるONSITE=TRUE) を含む場合は非対応



LSGEN (integer) {5}
LSMULLIKEN=TRUEにおけるデータ数の最大値



LSFITGEOM (logical) {FALSE}LSMULLIKEN=TRUEで過去のステップの原子座標の予測するステップの原子座標への重ね合わせ
クラスターモデルのみ使用可能
本オプションの有効化は不要

BFGS (OPTTYPE=1)の詳細オプションBFGSSOLVER (integer) {1}
BFGSのHessian行列に対する固有値と固有ベクトルの計算 (対角化) 法 1: DSYEVD
 2: DSYEVSteepest descent (OPTTYPE=2)の詳細オプション
STEPSIZE/STEPSTDC (real) {1.0e+02}
Steepest descentにおける変位を制御するパラメータ (atomic unit)

25

目次

Conjugate gradient (OPTTYPE=3)の詳細オプション


STEPCGRATOM (real) {1.0e−01}
Conjugate gradient における 原子 位置 の変 位を制 御す るパ ラメ ータ(atomic unit)



STEPCGRLATTICE (real) {1.0e−01}
Conjugate gradientにおける格子ベクトルの変位を制御するパラメータ(atomic unit)

Quick-min (OPTTYPE=4)の詳細オプション


DELTATQMIN/DELTASTEPQMIN (real) {1.0e−15}
Quick-minにおける時間刻み (s)FIRE (OPTTYPE=5)の詳細オプション


DELTATFIRE/DELTASTEPFIRE (real) {1.0e−15}
FIREにおける時間刻み (s)FIREMAXDELTAT/FIREMAXDELTASTEP (real) {1.0e−14}
FIREにおける最大時間刻み (s)FIRENMIN (integer) {5}
FIREにおけるパラメータ

FIREFINC (real) {1.1e+00}
FIREにおけるパラメータ

FIREFDEC (real) {5.0e−01}
FIREにおけるパラメータ



FIREALPHASTART (real) {1.0e−01}
FIREにおけるパラメータ

FIREFALPHA (real) {9.9e−01}
FIREにおけるパラメータ

26

目次

1.1.8 FREQ/FREQUENCY
FREQ/FREQUENCYキーワードは、調和振動解析を実行するか否かを指定する。
構造最適化計算 (OPT=TRUEあるいはOPTIMIZE=TRUE) と同時に有効化された場合、構造最適化最終構造に対する振動数計算がなされる。この計算は、クラスターモデルのみ実行可能である。

基本設定
FREQTYPE (integer) {1}
Hessian行列の構築法
 1: 解析的2次微分
以下の場合は使用不可
 開殻系の計算
 点 電 荷 あ り の 計 算(SCC キ ー ワ ー ド に お い て

POINTCHARGE=TRUE) 3次spline補間を用いたSCC相互作用の短距離項の計算 (SCCキーワードにおいてGAMMASPLINE=TRUE) DFTB+U 補 正 を 含 む 計 算(SCC キ ー ワ ー ド に お い て

PLUSU=TRUE)
 Onsite 補 正 を 含 む 計 算(SCC キ ー ワ ー ド に お い て

ONSITE=TRUE)
 MPI/OpenMP版で原子価軌道ベースの計算 (SCCキーワードに
おいてRUNTYPE=SHELL, SHELL0とORSCC=TRUE) 2: 1次微分の数値微分


STEPSIZE (real) {1.0e−04}
FREQTYPE=2における数値微分刻み幅 (atomic unit)

PRINTHESS (logical) {FALSE}
Hessian行列の下三角部分を標準アウトプットファイルに出力
リスタートファイル “restart_hess” の出力

PRECFORMAT (logical) {FALSE}
通常出力に加えて高精度フォーマットでも結果を出力
27

目次



READHESS (method) {NONE}
Hessian行列ファイルの読み込み
 NONE: 読み込みなし (解析的あるいは数値的に計算)
 BINARY: PRINTHESS=TRUEで生成される “restart_hess” ファイルの使用
 ASCII: “hess.dat” ファイルによる手動指定 (形式は1.5.2節を参照)

詳細オプションHESSSOLVER (integer) {1}
Hessian行列に対する固有値と固有ベクトルの計算 (対角化) 法 1: DSYEVD
 2: DSYEV

CPCONV (real) {1.0e−06}
SCCを用いた解析的Hessian計算におけるcoupled perturbed DFTB計算の収束条件 (atomic unit)

EFSTR (real) {2.0e−03}
SCCを用いたIR, Raman強度計算における外部静電場の強さ



PROJECTION (logical) {TRUE}
並進・回転モードを振動解析から分離して計算
PROJECTION=FALSEの場合、Hessian行列をそのまま対角化



READISOTOPE (logical) {FALSE}
同位体質量の指定 (詳細は1.4.2節を参照)



THERMOTEMP (real) {2.9815e+02}
熱化学解析における温度 (K)THERMOPRES (real) {1.0e+00}
熱化学解析における圧力 (atm)

THERMOSCALE (real) {1.0e+00}
熱化学解析における振動数のスケーリングファクター

28

目次



INTERRUPTCP (logical) {FALSE}SCCを用いた解析的Hessian計算におけるcoupled perturbed DFTB計算の中断
部分的に構築されたHessian行列は “hess.dat” ファイルに出力CPSTART (integer) {1}
INTERRUPTCP=TRUEにおける部分的なHessian行列構築を開始する原子インデックス

CPEND (integer) {1}
INTERRUPTCP=TRUEにおける部分的なHessian行列構築を終了する原子インデックス

INTERRUPTSEMINUM (logical) {FALSE}FREQTYPE=2におけるHessian行列構築の中断
部分的に構築されたHessian行列は “hess.dat” ファイルに出力

SEMINUMSTART (integer) {1}
INTERRUPTSEMINUM=TRUEにおける部分的なHessian行列構築を開始する原子インデックス



SEMINUMEND (integer) {1}
INTERRUPTSEMINUM=TRUEにおける部分的なHessian行列構築を終了する原子インデックス

MPI/OpenMP版のみの追加オプション


SAVEMEMORYANA (logical) {FALSE}SCCを用いた解析的Hessian計算時のメモリ使用量を削減
計算コストは増加

29

目次

1.1.9 MD
MDキーワードは、速度Verlet法 [44] に基づく分子動力学 (MD) 計算を実行するか否かを指定する。

基本設定


NSTEP/MAXSTEP (integer) {10000}
MD計算のステップ数

DELTAT/DELTASTEP (real) {1.0e−15}
MD計算における時間刻み (s)PRINT (integer) {1}
MD 計 算 中 の 原 子座標 、 速 度 、 Mulliken 電 荷 の “traject”, “velocity”,“mulliken” ファイルへの出力頻度



CHECKPOINT (integer) {100}
第2リスタートファイル “restart_chk” の出力頻度RESTART (logical) {FALSE}
リスタートファイル “restart” を読み込んでMD計算を再開



READVELOCITY (logical) {FALSE}
初期速度ファイルの読み込み
 TRUE: “veloc.dat” ファイルによる手動指定 (形式は1.5.2節を参照)
 FALSE: 読み込みなし

エネルギー一定アンサンブル

NVE (logical) {TRUE}
NVEアンサンブル計算
NVE=FALSEの場合、NVTアンサンブル計算



INITTEMP (real) {2.9815e+02}
初期速度生成における目標温度 (K)初期速度はMaxwell-Boltzmann分布に従ってランダムに生成

30

目次

温度一定アンサンブル

NVT (logical) {FALSE}
NVTアンサンブル計算
NVT=FALSEの場合、NVEアンサンブル計算



BATHTEMP/TEMPERATURE (real) {2.9815e+02}
温度一定MD計算における熱浴温度 (K)

NVTTYPE (integer) {3}
温度制御法 1: Velocity scaling [45]
 2: Nosé-Hoover chain [46]
 3: Nosé-Hoover chain [46]NVTTYPE=2と比較して、いくつかの熱浴設定を変更可能
 4: Berendsen [47]
 5: Andersen [48]

ERRORTEMP (real) {1.0e+04}
BATHTEMP/TEMPERATUREで指定した熱浴温度からの瞬間温度の許
容振れ幅 (K)NVTTYPE=2–5において振れ幅を超えると速度スケーリングを実行



RMROTATION (logical) {FALSE}
重心の回転運動を除去
クラスターモデルのみ使用可能



RMSTEP (integer) {10000}
NVTTYPE=1–4において重心の並進 (RMROTATION=TRUEの場合、回転) 運動を除去する頻度



CALCPRESSURE (logical) {FALSE}
NVEおよびNVTアンサンブルにおける圧力計算
周期境界条件を課したモデルのみ計算可能

31

目次

温度制御の詳細オプション

SCALESTEP (integer) {20}Velocity scaling (NVTTYPE=1) における速度スケーリングの頻度

NOSECLEAN (logical) {FALSE}
Nosé-Hoover chain (NVTTYPE=2, 3) においてMD計算再開時に熱浴パ
ラメータを初期化



NBATH (integer) {3}Nosé-Hoover chain (NVTTYPE=3) における熱浴の数



COUPLINGSTR (real) {1.5e+03}
Nosé-Hoover chain (NVTTYPE=3) におけるカップリングの強さ (cm−1)



TAUTEMP (real) {1.0e−13}
Berendsen (NVTTYPE=4) における緩和時間 (s)

COLLISIONFREQ (real) {1.0e−02}
Andersen (NVTTYPE=5) における仮想熱浴との衝突頻度 (fs−1)



MASSIVE (logical) {FALSE}
Andersen (NVTTYPE=5) における速度スケーリングの方法
 TRUE: 各原子の速度をある確率で同時にスケーリング
 FALSE: 原子毎に速度をある確率でスケーリング

圧力一定アンサンブル

NPH (logical) {FALSE}
NPHアンサンブル計算
NPH=FALSEの場合、NVEアンサンブル計算



NPT (logical) {FALSE}
NPTアンサンブル計算
NPT=FALSEの場合、NVEアンサンブル計算

32

目次



EXTPRESSURE (real) {1.01325e+05}
圧力一定MD計算における目標圧力 (Pa)

NPHTYPE (integer) {1}
圧力制御法
 1: Berendsen [47]

NPTTYPE (integer) {1}(等方的) NPTアンサンブル計算法
 1: Hoover-Andersen [46]ISOSCALING (logical) {FALSE}
等方的なセルサイズ変化による圧力一定MD計算
デフォルトは異方的なセルサイズ変化

圧力制御の詳細オプションTAUPRESS (real) {2.0e−12}
Berendsen (NPHTYPE=1) における緩和時間 (s)

COUPLINGSTRLAT (real) {1.5e+01}
Hoover-Andersen (NPTTYPE=1) における熱浴のカップリングの強さ(cm−1)

乱数生成
RANDTYPE (integer) {1}
乱数生成方法 1: Fortran組み込み手続き
 2: 参考文献 [49] の乱数生成アルゴリズムに従った手続き



SEEDTYPE (integer) {1}
速度のための初期乱数の指定
 1: 初期構造に依存した乱数
 2: 現在時刻から乱数を作成
 3: RANDOMSEEDオプションにより任意の乱数を作成

33

目次



RANDOMSEED (integer) {0}SEEDTYPE=3における変位のための乱数作成のための番号 (正の整数)

ソフトポテンシャル

SOFT/SOFTPOT (logical) {FALSE}
ソフトポテンシャルの使用



SOFTK (real) {1.5e+00}
ソフトポテンシャルの硬さ (kcal/(mol∙Å2))



SOFTCENTERTYPE (method) {ATOM}
ソフトポテンシャルの中心位置
 ATOM: 手動指定した原子 SPECIFY: 手動指定した座標
 COC: 系の座標中心
 COM: 系の重心

SOFTATOMNO (integer) {1}
SOFTCENTERTYPE=ATOMにおけるソフトポテンシャルの中心となる原子のラベル

SOFTPOINTXYZ/SOFTPOINT (real) {0.0e+00}
SOFTCENTERTYPE=SPECIFYにおけるソフトポテンシャルの中心となるxyz座標 (x = y = z) (Å)



SOFTPOINTX (real) {0.0e+00}
SOFTCENTERTYPE=SPECIFYにおけるソフトポテンシャルの中心となるx座標 (Å)



SOFTPOINTY (real) {0.0e+00}
SOFTCENTERTYPE=SPECIFYにおけるソフトポテンシャルの中心となるy座標 (Å)



SOFTPOINTZ (real) {0.0e+00}
SOFTCENTERTYPE=SPECIFYにおけるソフトポテンシャルの中心となるz座標 (Å)
34

目次



SOFTSHAPETYPE (method) {SPHERE}
ソフトポテンシャルの形状
 SPHERE: 球 SQUARE: 立方体あるいは直方体
 ELLIPSOID: 楕円体



SOFTRANGE (real) {1.0e+01}
SOFTSHAPETYPE=SPHEREにおけるソフトポテンシャルの半径 (Å)



SOFTSIDEXYZ/SOFTSIDE (real) {1.0e+01}
SOFTSHAPETYPE=SQUAREにおけるソフトポテンシャルのグリッド
一辺の長さ (Å)



SOFTSIDEX (real) {1.0e+01}
SOFTSHAPETYPE=SQUAREにおけるソフトポテンシャルのグリッドのx成分の長さ (Å)



SOFTSIDEY (real) {1.0e+01}
SOFTSHAPETYPE=SQUAREにおけるソフトポテンシャルのグリッドのy成分の長さ (Å)



SOFTSIDEZ (real) {1.0e+01}
SOFTSHAPETYPE=SQUAREにおけるソフトポテンシャルのグリッドのz成分の長さ (Å)



SOFTSEMIAXISX (real) {1.0e+01}
SOFTSHAPETYPE=ELLIPSOIDにおけるソフトポテンシャルの径のx成分の長さ (Å)



SOFTSEMIAXISY (real) {1.0e+01}
SOFTSHAPETYPE=ELLIPSOIDにおけるソフトポテンシャルの径のy成分の長さ (Å)



SOFTSEMIAXISZ (real) {1.0e+01}
SOFTSHAPETYPE=ELLIPSOIDにおけるソフトポテンシャルの径のz成分の長さ (Å)
35

目次

RATTLE法を用いた拘束動力学


RATTLE (method) {NONE}RATTLE法 [50] に基づく原子間距離に対する拘束動力学
 NONE: 拘束条件なし
 COVALENTXH: X–Hペアの原子間距離を共有結合半径 [25] の和に拘束
 INITGEOMXH: X–Hペアの原子間距離を初期構造の値に拘束 SPECIFY: “rattle.dat” ファイルによる手動指定 (形式は1.5.2節を参照)

RATTLEXH (real) {1.2e+00}
RATTLE=COVALENTXHとRATTLE=INITGEOMXHにおけるX–Hペアに対する共有結合の判定基準 (Å)



RATTLEITER (integer) {100}RATTLE法における拘束条件を満たすための反復サイクルの最大数



RATTLECONV (real) {1.0e−06}
RATTLE法における拘束条件に対する収束条件 (atomic unit)

SCC収束の高速化


LIMULLIKEN (logical) {FALSE}
Lagrange補間によるMulliken電荷初期値の予測 [43, 51]
DFTB+U補正 (SCCキーワードにおけるPLUSU=TRUE) およびonsite補
正 (SCCキーワードにおけるONSITE=TRUE) を含む場合は非対応



LIDEGREE (integer) {0}
LIMULLIKEN=TRUEにおけるLagrange補間多項式の次数デフォルトは前ステップの値のみから予測

36

目次



LIMINDEG (logical) {FALSE}
LIMULLIKEN=TRUEでLagrange補間多項式の次数を毎ステップ最適化
0以上LIDEGREEオプションの値以下の整数値が探索範囲
Nステップ目における収束電荷と予測電荷の絶対誤差が最小となる次数をN + 1ステップ目に使用



LSMULLIKEN (logical) {FALSE}
原子座標の最小二乗フィッティングによるMulliken電荷初期値の予測[42, 43].
DFTB+U補正 (SCCキーワードにおけるPLUSU=TRUE) およびonsite補
正 (SCCキーワードにおけるONSITE=TRUE) を含む場合は非対応



LSGEN (integer) {5}
LSMULLIKEN=TRUEにおけるデータ数の最大値



LSFITGEOM (logical) {FALSE}LSMULLIKEN=TRUEで過去のステップの原子座標の予測するステップの原子座標への重ね合わせ
クラスターモデルのみ使用可能
本オプションの有効化は不要

アニーリング

ANNEAL (logical) {FALSE}
アニーリング計算の実行
温度のスケジュールは “anneal.dat” ファイルで最大99個まで手動指定(形式は1.5.2節を参照)

メタダイナミクスMETAD/METADYNAMICS (logical) {FALSE}
メタダイナミクス計算 [52] の実行
集団変数 (CV) は “metacv.dat” ファイルで指定 (形式と実装済みCVは
1.5.2節を参照)

37

目次



METAMAXCV (integer) {10}メタダイナミクス計算におけるCVの最大数



METAFREQ (integer) {100}メタダイナミクス計算におけるガウス関数の堆積頻度



METAMAXGAUSS (integer) {10000}
メタダイナミクス計算におけるガウス関数の最大堆積数



METAHEIGHT (real) {4.0e−04}
メタダイナミクス計算におけるガウス関数の高さ (Hartree)

METARESTART (logical) {FALSE}
リスタートファイル “restart_meta” を読み込んでメタダイナミクス計
算を再開



METAREADBIAS (logical) {FALSE}
メタダイナミクス計算におけるガウス関数堆積履歴ファイルの読み込み TRUE: “bias.dat” ファイルによる手動指定 (形式は1.5.2節を参照)
 FALSE: 読み込みなし

メタダイナミクスの詳細オプション
WTMETAD/WTMETADYNAMICS (logical) {FALSE}
Well-temperedメタダイナミクス計算 [53] の実行



WTMETABIASTEMP (real) {1.2e+03}
Well-temperedメタダイナミクス計算におけるバイアス温度 (K)

METAPBCMIC (logical) {TRUE}
周 期 系 の メ タ ダ イ ナ ミ ク ス 計 算 で CV の 原 子 間 距 離 計 算 に 対 し て
minimum image conventionを適用
系の座標中心と重心および一部のCVの計算には不適用 (詳細は1.5.2節
を参照)

38

目次



METAPRINTFES (logical) {TRUE}
メタダイナミクス計算において自由エネルギー面のデータをガウス関数堆積毎に “fes.dat” ファイルに出力
自由エネルギー面のグリッド情報を “metacv.dat” ファイルに指定 (詳細は1.5.2節を参照)
Well-temperedメタダイナミクス計算ではスケーリングファクター [53]適用後の結果を出力



METAWALL (logical) {FALSE}
メタダイナミクス計算においてCVのサンプリング範囲を拘束するポテ
ンシャル障壁の適用
ポテンシャル障壁の定義を “metacv.dat” ファイルに指定 (詳細は1.5.2
節を参照)

39

目次

1.1.10 MISCMISCキーワードは、プログラム実行に関するその他の条件を指定する。

追加解析とプロパティ評価


PRINTFORCE/FORCE (logical) {FALSE}
エネルギーと力の計算
点電荷ありの計算 (SCCキーワードにおいてPOINTCHARGE=TRUE)
では点電荷の力も出力
応力テンソル計算 (PBCキーワードにおいてSTRESS=TRUE) では応力テンソル、格子ベクトルの微分、体積、圧力も出力



PRINTRESOLVEDFORCE (logical) {FALSE}
全力に対する各成分の寄与 (electronic, repulsive, dispersionが含まれる
場合はdispersion項) を出力
応力テンソル計算 (PBCキーワードにおいてSTRESS=TRUE) では全応力テンソルに対する各成分の寄与も出力
全力は標準アウトプットファイル、各成分の寄与は詳細データファイルに出力

PRINTMO (logical) {FALSE}
分子軌道係数を出力
部分系の数が1の場合のみ出力可能



PRINTDENSITY (logical) {FALSE}
密度行列の下三角部分を標準アウトプットファイルに出力

PRINTATOME (logical) {FALSE}
全エネルギーに対する各原子の寄与 (electronic, repulsive, dispersionが含まれる場合はdispersion項) を出力
各成分の寄与は詳細データファイルに出力

40

目次

グローバルオプション


FORCETYPE (integer) {1}エネルギー勾配の計算方法
 1: 解析的1次微分
 2: 全エネルギーの数値微分



FORCESTEPSIZE (real) {1.0e−04}
FORCETYPE=2における数値微分刻み幅 (atomic unit)

HSSTEPSIZE (real) {1.0e−02}
DFTB電荷非依存項の微分計算における数値微分刻み幅 (atomic unit)

WRITEDAT (method) {APPEND}
詳細データファイルの書き込み
 NONE: 書き込みなし
 APPEND: 各DFTBステップの情報を追加書き込み
 REPLACE: 各DFTBステップの情報を新規書き込み



WRITECHRGFILE (logical) {TRUE}
電荷密度ファイルの読み込み(SCC キ ー ワ ー ド に お け る

READCHARGE=BINARY) で用いる “chrgfile” ファイルの書き込み


PRINTLAPTIME (logical) {FALSE}DFTB計算の各手続きとプログラム実行にかかった経過時間の出力

PRINTLOCALREGION (logical) {FALSE}
局在化領域に含まれる部分系、原子、基底関数の数を標準アウトプット
ファイルに出力
出力後はプログラムが直ちに終了



PRINTHS (logical) {FALSE}DFTB電荷非依存な0次Hamiltonianと重なり行列の下三角部分を標準ア
ウトプットファイルに出力
部分系の数が1の場合のみ出力可能
出力後はプログラムが直ちに終了
41目次

Serial, OpenMP版のみの追加オプション


USESCRATCH (logical) {FALSE}
中間データに対するディスクI/Oの使用
USESCRATCH=FALSEの場合、中間データをメモリに格納

42

目次

1.2 Titleセクション
コメント行

TITLE

Titleセクションは、化合物名、分子の対称性、電子状態といった計算に関係する情報を付記するために用いられる。1行あたり128文字まで入力可能である。全
コメント行は標準アウトプットファイルにそのまま出力される (3.1節を参照) 。

43

目次

1.3 Parameterセクション
1.3.1 概要2

原子の種類の数O˽2

原子1の名前と原子価軌道のサイズ

oo.spl˽oh.spl

原子1–原子1と原子1–原子2のパラメータH˽1

原子2の名前と原子価軌道のサイズ

ho.spl˽hh.spl

原子2–原子1と原子2–原子2のパラメータ

Parameterセクションは、入力構造に含まれる原子の種類の数の指定から始まる。
水分子の例では、O原子とH原子の合計で2となる。次の行では、1つ目の原子の種類 (原子1) について原子の種類の記号と原子価軌道のサイズを指定する。後者の指定で用いる整数値1, 2, 3, 4は順にs, p, d, f軌道まで考慮することに対応する。改行後、原子1を含むパラメータファイルのパスをこのセクションで指定する原子の種類の順番を厳守して列挙する。すなわち、原子1–原子1のパラメータ、
原子1–原子2 (2つ目の原子の種類) のパラメータ、原子1–原子3 (3つ目の原子の種類) のパラメータ、 … のように指定する。原子1についての入力後、同様の入力を原子2、原子3、 … について繰り返す。この例では、4つのパラメータファイル (oo.spl, oh.spl, ho.spl, hh.spl) がインプットファイルと同じディレクトリ
に格納されていることを想定している。カレントディレクトリを明記した入力
は以下のようになる。
2O˽2./oo.spl˽./oh.spl
H˽1
./ho.spl˽./hh.spl

f電子を取り扱う場合、関連するパラメータファイルは拡張フォーマットである
ことを確実にする必要がある。

44

目次

1.3.2 定数パラメータの追加指定
いくつかのオプションは、プログラムのソースコードあるいはパラメータファ
イルに記述されていない原子の種類に依存する定数パラメータを必要とする。
定数パラメータは、以下の2つの方法で指定される。


原子価軌道のサイズの後に定数パラメータを直接指定
この方法は、SCC Hamiltonianに対する3次補正 (SCCキーワードにおける THIRDFULL=TRUE あ る い は THIRDDIAG=TRUE) で Hubbard
derivativeパラメータを指定する場合に必要となる。MPI/OpenMP版でorbitally resolved SCC計算 (SCCキーワードにおけるORSCC=TRUE) の場合、s, p, d, f原子価軌道に対する値を別々に指定する。

O˽2˽−0.1575

原子ベースの場合 (通常の指定)Cu˽3˽−0.20˽−0.0575˽−0.0575

原子価軌道ベースの場合 (順番は(f), d, p, s)

定数パラメータを含むファイルを追加指定
この方法は、以下の5つで定数パラメータを指定する場合に必要となる。
 スピン分極を含む場 合の スピン 定数 (SCC キーワードにおけ る
SPIN=TRUE) DFTB+U補正を含む場合のU–Jパラメータ (SCCキーワードにおけるPLUSU=TRUE)
 Onsite補正を含む場合のonsiteパラメータ (SCCキーワードにおけるONSITE=TRUE)
 Slater-Kirkwood分散力補正を含む場合の計算に必要となるパラメー
タ (DISP/DISPERSIONキーワードにおけるDISPTYPE=2)
 Lennard-Jones分散力補正を含む場合のLennard-Jonesパラメータ(DISP/DISPERSION キ ー ワ ー ド に お け る DISPTYPE=3 か つ
READUFF=TRUE)
追加指定するファイルのパスは、原子価軌道のサイズの後 (Hubbard
derivativeパラメータがある場合はその後) に記述する。複数のファイルを指定する場合、分散力補正用パラメータ、スピン定数、U–Jパラメータ、
そしてonsiteパラメータを含むファイルの順に記述する。
45

目次

O原子に対する指定例を以下に示す。
O˽2˽o.wll

スピン定数ファイルのみ指定

O˽2˽o.dispsk

分散力補正用パラメータファイルのみ指定

O˽2˽o.uj

U–Jパラメータファイルのみ指定

O˽2˽o.oc

Onsiteパラメータファイルのみ指定

O˽2˽o.dispsk˽o.wll

スピン定数と分散力補正用パラメータを含むファイルを指定O˽2˽o.dispsk˽o.uj

分散力補正用パラメータとU–Jパラメータを含むファイルを指定

O˽2˽o.dispsk˽o.oc

分散力補正用パラメータとonsiteパラメータを含むファイルを指定

O˽2˽o.wll˽o.uj

スピン定数とU–Jパラメータを含むファイルを指定

O˽2˽o.wll˽o.oc

スピン定数とonsiteパラメータを含むファイルを指定

O˽2˽o.uj˽o.oc

U–Jパラメータとonsiteパラメータを含むファ

イルを指定
O˽2˽o.dispsk˽o.wll˽o.uj

分散力補正用パラメータ、スピン定数、U–Jパ
ラメータを含むファイルを指定

O˽2˽o.dispsk˽o.wll˽o.oc

分散力補正用パラメータ、スピン定数、onsite
パラメータを含むファイルを指定

O˽2˽o.dispsk˽o.ud˽o.oc

分散力補正用パラメータ、 U–J パラメータ、
onsiteパラメータを含むファイルを指定

O˽2˽o.wll˽o.ud˽o.oc

スピン定数、U–Jパラメータ、onsiteパラメータを含むファイルを指定

O˽2˽o.dispsk˽o.wll˽o.ud˽o.oc

分散力補正用パラメータ、スピン定数、U–Jパ
ラメータ、onsiteパラメータを含むファイルを
全て指定

46

目次

1.3.3 定数パラメータ追加指定ファイルの形式
1.3.2節に記述した定数パラメータを指定するためのファイルは、それぞれ異なる形式を持つ。C原子に対する指定例を以下に示す。

スピン定数ファイル−0.031
−0.025
−0.025
−0.023スピン定数ファイルの行数は、原子価軌道のサイズに依存する (s, p, d, f軌道に対して1, 4, 9, 16) 。各行では、原子のスピン定数 (atomic unit) を以下の順番で
指定する。



s: ss
p: ss, sp, ps, pp

d: ss, sp, sd, ps, pp, pd, ds, dp, dd
f: ss, sp, sd, sf, ps, pp, pd, pf, ds, dp, dd, df, fs, fp, fd, ff

参考までに、DFTB+のマニュアル [54] にまとめられたH, C, N, O, S, Fe, Ni原子に対するPBEの値を含む拡張子 “.wll” のファイルがslkoディレクトリに格納されている。
U–Jパラメータファイル
0.178371875˽1
0.178371875˽1

U–Jパラメータファイルの行数は、原子価軌道のサイズに依存する (s, p, d, f軌道に対して1, 2, 3, 4) 。各行は2列からなり、第1列はon-site Coulomb相互作用に関するパラメータ (atomic unit)、第2列は原子価軌道の属するグループのインデッ
クスを指定する。例えば、s軌道とp軌道を異なるグループとして取り扱う場合、第2行第2列は “1” から “2” へ変更する。 ある原子の種類に対してDFTB+U補正を適用しない場合であっても、第1列の値を0としてこのファイルを準備する必
要がある。
47

目次

Onsiteパラメータファイル−0.00000˽0.00000
−0.04973˽0.10512
−0.04973˽0.10512
−0.01203˽0.02643

Onsiteパラメータファイルの行数は、原子価軌道のサイズに依存する (s, p, d, f
軌道に対して1, 4, 9, 16) 。各行は2列からなり、第1列は同じスピンに対するパ
ラメータ (atomic unit)、第2列は異なるスピンに対するパラメータ (atomic unit)
を以下の順番で指定する。





s: ss
p: ss, sp, ps, pp
d: ss, sp, sd, ps, pp, pd, ds, dp, dd
f: ss, sp, sd, sf, ps, pp, pd, pf, ds, dp, dd, df, fs, fp, fd, ff

参考までに、参考文献 [14] あるいはDFTB+のマニュアル [54] にまとめられた
H, C, N, O, S原子に対するPBEの値を含む拡張子 “.oc” のファイルがslkoディレ
クトリに格納されている。
Slater-Kirkwood分散力補正用パラメータファイル−1.382˽1.382˽1.382˽1.064˽3.8˽3.8˽3.8˽3.8˽2.5

このファイルは1行9列からなり、第1–4列は原子分極率 (Å3)、第5–8列はvan der
Waals半径 (Å)、第9列は有効電荷を指定する。原子分極率とvan der Waals半径の値は、配位数0, 1, 2, 3に対して順番に指定される。参考までに、DFTB+のマニュアル [54] にまとめられたH, C, N, O, S原子に対する分散力補正用パラメータを含む拡張子 “.dispsk” のファイルがslkoディレクトリに格納されている。

48

目次

Lennard-Jones分散力補正用パラメータファイル−3.851˽0.105

このファイルは1行2列からなり、第1列は非結合性距離 (Å)、第2列は非結合性エネルギー (kcal/mol) を指定する。参考までに、参考文献 [30] にまとめられたH,
C, N, O, S, Fe, Ni原子に対する分散力補正用パラメータを含む拡張子 “.displj” のファイルがslkoディレクトリに格納されている。

49

目次

1.4 Geometryセクション
1.4.1 概要
原子数、電荷、スピン

3˽0˽1

多重度
O0˽−0.06275391˽−0.06275391˽−0.00000000

原子ラベル1の座標

H0˽−1.01910555˽−0.08185946˽−0.00000000

原子ラベル2の座標H0˽−0.08185946˽−1.01910555˽−0.00000000

原子ラベル3の座標

Geometryセクションは、入力構造の原子数、電荷、スピン多重度という3つの整数値の指定から始まる。改行後、原子の位置を小数表記のデカルト座標 (Å) で指定する (整数表記は不可) 。
入力構造の全ての原子の位置を指定した後、ラベルTVを用いて3つの並進ベクトルを指定すると周期境界条件を課したモデルとなる (TVは大文字) 。以下は、ダイヤモンド結晶構造 (基本単位格子) に対する指定例である。
2˽0˽1
C0˽−0.00000000˽−0.00000000˽−0.00000000
C0˽−0.89169750˽−0.89169750˽−0.89169750
TV˽−0.00000000˽−1.78339500˽−1.78339500
TV˽−1.78339500˽−0.00000000˽−1.78339500
TV˽−1.78339500˽−1.78339500˽−0.00000000

50

目次

外部静電場を導入する場合 (SCCキーワードにおけるEFIELD=TRUE) 、入力構
造の全ての原子の位置を指定した後、ラベルEFを用いて電場ベクトルを指定する (EFは大文字) 。以下は、水分子に対して外部静電場 (−0.001, −0.001, −0.000)を与える例である。
3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
H0˽−1.01910555˽−0.08185946˽−0.00000000
H0˽−0.08185946˽−1.01910555˽−0.00000000
EF˽−0.00100000˽−0.00100000˽−0.00000000

外部静電場を導入する場合 (SCCキーワードにおけるPOINTCHARGE=TRUE) 、
入力構造の全ての原子の位置を指定した後、ラベルPCを用いて点電荷の座標と値を指定する (PCは大文字) 。 以下は、水分子に対して電荷+1の点電荷を位置(5.0, 5.0, 5.0) (Å) に与える例である。
3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
H0˽−1.01910555˽−0.08185946˽−0.00000000
H0˽−0.08185946˽−1.01910555˽−0.00000000
PC˽−5.00000000˽−5.00000000˽−5.00000000˽−1.00000000

51

目次

1.4.2 追加情報の指定
以下のようにシミュレーション中に原子位置に制約をかけたい場合や、調和振動解析で同位体質量を指定したい場合には、原子のデカルト座標の後に列を追加する。


構造最適化あるいはMD計算において原子位置を固定
アスタリスク記号*を指定すると、原子位置は固定される。



制約付き構造最適化(OPT/OPTIMIZE キ ー ワ ー ド に お け る

CONSTRAINT=TRUE)
x, y, z成分に対する制約は、値0.0 (変位不可) あるいは1.0 (変位可) を用いて指定する。例えば、制約ベクトル “0.0˽0.0˽1.0” を課された原子は、z座標のみ値が更新される。制約ベクトル “0.0˽0.0˽0.0” は、アスタリス
ク記号の指定と同義である。


調和振動解析における原子質量の変更 (FREQ/FREQUENCYキーワード
におけるREADISOTOPE=TRUE)
パラメ ータ ファ イル に記述 され た原 子質 量の初 期設 定値 は、 形式“ISO=value” により値valueへ変更される (ISO=とvalueの間の空白領域
はなし) 。質量を変更する原子の位置に制約をかける場合、アスタリス
ク記号あるいは制約ベクトルの後にこの形式を指定する。

以下は、水分子に対して上記の指定方法が混在した例である。
3˽0˽1
O0˽−0.06275391˽−0.06275391˽−0.00000000
H0˽−1.01910555˽−0.08185946˽−0.00000000˽0.0˽1.0˽0.0
H0˽−0.08185946˽−1.01910555˽−0.00000000˽*˽ISO=2.014102

制約付き構造最適化の場合、原子ラベル2のH原子はy座標のみ値が更新される。
また、原子ラベル3のH原子の位置はシミュレーション中に更新されない。
FREQ=(READISOTOPE=TRUE)が追加指定された場合、原子ラベル3のH原子の質量は2.014102に変更される。すなわち、HOD分子の調和振動数と熱力学量が
得られる。
52

目次

1.4.3 DC-DFTBにおける部分系の手動指定
本 節 で は 、 DC キ ー ワ ー ド に お い て SUBTYPE=SPECIFY あ る い は
SUBTYPE=SEMIAUTOを用いて部分系を指定する方法を補足する。Geometryセクションの終わりに挿入された空行の後で部分系のインデックスを列挙する。
以下は、22原子からなる入力構造に対してSUBTYPE=SPECIFYを用いて5つの部分系を指定する例である。
1˽1˽1˽1˽1

原子ラベル1, 2, 3, 4, 5に対する部分系のインデックス2˽2˽2˽2

原子ラベル6, 7, 8, 9に対する部分系のインデックス3˽3˽3˽3

原子ラベル10, 11, 12, 13に対する部分系のインデックス

4˽4˽4˽4

原子ラベル14, 15, 16, 17に対する部分系のインデックス

5˽5˽5˽5˽5

原子ラベル18, 19, 20, 21, 22に対する部分系のインデックス

この例の1行目は、原子ラベル1–5の原子は部分系1に属することを指定する。2,3行目では、原子ラベル6–9, 10−13の原子は部分系2, 3に属することが指定され、同様の指定がさらに続く。指定に用いる行数は部分系の数と同じである必要はない。プログラムは、指定された部分系のインデックスの数が入力構造の原子数に到達するまで新しい行の読み込みを試みる。
SUBTYPE=SEMIAUTOの場合、部分系を自動決定させる原子のインデックスとして “0” を用いる。以下は、22原子からなる入力構造に対する指定例である (上記との比較のため、指定形式を揃えている) 。
1˽1˽1˽1˽1

原子ラベル1, 2, 3, 4, 5に対する部分系のインデックス0˽0˽0˽0

原子ラベル6, 7, 8, 9に対する部分系のインデックス0˽0˽0˽0

原子ラベル10, 11, 12, 13に対する部分系のインデックス0˽0˽0˽0

原子ラベル14, 15, 16, 17に対する部分系のインデックス2˽2˽2˽2˽2

原子ラベル18, 19, 20, 21, 22に対する部分系のインデックス

この例では、原子ラベル1–5, 18–22の原子は部分系1, 2に属する。原子ラベル6–17の原子は、部分系3–Xに自動的に割り当てられる (Xの値は入力構造と自動決定で用いるグリッドの長さに依存する) 。全てのインデックスが0あるいは正の整数の場合、SUBTYPE=AUTOあるいはSUBTYPE=SPECIFYと同義である。
53

目次

1.5 補助インプットファイル
インプットファイル、パラメータファイル、および1.3.3節に示した定数パラメー
タ追加指定ファイルに加えて、いくつかの補助インプットファイルが計算条件に依存して必要となる。各補助インプットファイルは、固定ファイル名でインプットファイルと同じディレクトリに格納して使用する。その方法には、(i) 前回のプログラム実行で出力されたバイナリファイルの使用と (ii) 固定形式のASCII
ファイルの準備という2種類がある。
1.5.1 バイナリ形式ファイル
バイナリ形式ファイルは、計算再開時の情報を提供する目的で用いられる。各フ
ァイルの内容については、3.3節あるいは関連するオプションの説明に譲る。バイナリ形式ファイルは、異なるバージョン間で互換性を持たない ((i) OpenMP(serial) 版とMPI/OpenMP版や (ii) バージョン1.0と2.0の間での流用不可) 。
“chrgfile”
SCC=(READCHARGE=BINARY)が指定された場合に必要
前 回 の プ ロ グ ラ ム 実 行 で SCC=TRUE と MISC キ ー ワ ー ド に お け る
WRITECHRGFILE=TRUEの指定が必要

“restart_hess”
FREQ=(RESTARTHESS=BINARY)が指定された場合に必要
前 回 の プ ロ グ ラ ム 実 行 で FREQ/FREQUENCY キ ー ワ ー ド に お け る
PRINTHESS=TRUEの指定が必要

“restart”
MD=(RESTART=TRUE)が指定された場合に必要“restart” ファイルに加えて “restart_chk” ファイルをリネーム後に使
用可能

54

目次



“restart_meta”
MD=(METARESTART=TRUE)が指定された場合に必要“restart_meta” ファイルに加えて “restart_meta_chk” ファイルをリネーム後に使用可能
前回と今回のプログラム実行でMDキーワードにおけるMETAD=TRUE,METADYNAMICS=TRUE,WTMETAD=TRUE,WTMETADYNAMICS=TRUEの指定が必要

55

あ

る

い

は

目次

1.5.2 ASCII形式ファイル
ASCIIファイル形式は、主に入力構造に依存する詳細情報を提供する目的で用いられる。
SCC計算における初期電荷の読み込み“charge.dat”
SCC=(READCHARGE=ASCII)が指定された場合に必要
前回のプログラム実行でSCC=TRUEとSCC=TRUEの指定が必要
以下は、O原子の原子ラベルが1、H原子の原子ラベルが2, 3の場合の水分子に対する指定例である。
Serial, OpenMP版では、 “charge.dat” ファイルの1行目は “ALLATOM”
か “ATOMSPECIES” のいずれかである。
 ALLATOMでは、2行目以降に正味電荷の値を原子毎に指定する。
ALLATOM−0.6

O原子 (原子ラベル1) のMulliken電子密度: 6.6−0.3

H原子 (原子ラベル2) のMulliken電子密度: 0.7−0.3

H原子 (原子ラベル3) のMulliken電子密度: 0.7

 ATOMSPECIESでは、2行目以降の第1列に原子の種類の記号、第2列に正味電荷の値を指定する。プログラムは、空行を検出するまで新
しい行の読み込みを試みる。
ATOMSPECIES
O˽−0.6O原子 (原子ラベル1) のMulliken電子密度: 6.6

H˽−0.3

H原子 (原子ラベル2, 3) のMulliken電子密度: 0.7

56

目次

MPI/OpenMP版では、 “charge.dat” ファイルの1行目は “ALLSHELL” か“SHELLSPECIES” のいずれかである。
 ALLSHELLでは、2行目以降に正味電荷の値を原子価軌道毎に指定する。
ALLSHELL−0.2

O原子 (原子ラベル1) s軌道のMulliken電子密度: 1.8−0.8O原子 (原子ラベル1) p軌道のMulliken電子密度: 4.8−0.3

H原子 (原子ラベル2) s軌道のMulliken電子密度: 0.7−0.3

H原子 (原子ラベル3) s軌道のMulliken電子密度: 0.7

 SHELLSPECIESでは、2行目以降の第1列に原子の種類の記号、第2列に原子価軌道の記号 (s/p/d/f) 、第3列に正味電荷の値を指定する。プログラムは、空行を検出するまで新しい行の読み込みを試みる。開殻系の場合、α/βスピンを区別する記号 (a/b) を第4列に追加指定する必要がある。
SHELLSPECIESO˽s˽−0.2O原子 (原子ラベル1) s軌道のMulliken電子密度: 1.8O˽p˽−0.8O原子 (原子ラベル1) p軌道のMulliken電子密度: 4.8H˽s˽−0.3

H 原子 ( 原子ラベル 2, 3) s 軌道の Mulliken 電子密度 :
0.7

57

目次

振動数計算におけるHessian行列の読み込み“hess.dat”
FREQ=(RESTARTHESS=ASCII)が指定された場合に必要
前 回 の プ ロ グ ラ ム 実 行 で 解 析 的 あ る い は 数 値 的 Hessian 計 算(FREQ/FREQUENCY キ ー ワ ー ド に お け る FREQTYPE=1 あ る い は
FREQTYPE=2)に お い て INTERRUPTCP=TRUE あ る い は

INTERRUPTSEMINUM=TRUEの指定が必要
各行は3列からなり、第1, 2列はHessian行列要素のインデックス、第3列
はHessian行列要素の値を指定する。非対角要素は、ファイル中で1度だけ指定する。すなわち、プログラムはN原子からなる入力構造に対して3N(3N + 1)/2行のファイル読み込みを想定している。
1˽1˽−0.000171

Hessian行列要素 (H)11: −0.000171

2˽1˽−0.000000

Hessian行列要素 (H)21 と (H)12: 0.000000

3˽1˽−0.000083

Hessian行列要素 (H)31 と (H)13: 0.000083…

MD計算における初期速度の読み込み“veloc.dat”
MD=(READVELOCITY=TRUE)が指定された場合に必要
各行は3列からなり、速度のx, y, z成分 (Å/fs) を原子毎に指定する (原子の種類の記号は不要) 。
−0.00154193˽−0.00142221˽−0.00058536

原子ラベル1の速度成分−0.02998618˽−0.02026292˽−0.00706185

原子ラベル2の速度成分−0.00549584˽−0.00232603˽−0.00223541

原子ラベル3の速度成分

…

58

目次

RATTLE法における結合距離拘束の手動指定“rattle.dat”
MD=(RATTLE=SPECIFY)が指定された場合に必要各行は最低3列からなり、第1, 2列は距離を拘束する原子ラベル、第3列は結合距離の拘束方法を表すアルファベット大文字を指定する。プログ
ラムは、ファイルの最終行まで空行なくデータがあることを想定している。
 C: 結合距離を共有結合半径 [25] の和に拘束 I: 結合距離を初期構造の値に拘束 S: 結合距離を第4列で指定する値 (Å) に拘束
第4列が未指定の場合、共有結合半径の和に拘束
1˽2˽C

原子ラベル 1–2間の結合距離を共有結合半径の和に
拘束2˽3˽I

原子ラベル2–3間の結合距離を初期構造の値に拘束1˽3˽S˽0.95

原子ラベル1–3間の結合距離を0.95 Åに拘束

59

目次

アニーリングにおける温度スケジュールの設定“anneal.dat”
MD=(ANNEAL=TRUE)が指定された場合に必要各行は最低3列からなり、第1列は温度の変化方法を表すアルファベット
大文字、第2列は当該スケジュールを使用するステップ数、第3列は当該
スケジュールの最終ステップにおける目標温度 (K) を指定する。プログ
ラムは、ファイルの最終行まで空行なくデータがあることを想定している。
 C: 温度は一定値を維持 L: 温度は前スケジュールの値から線形に変化
 G: 温度は前スケジュールの値から等比級数的に変化 E: 温度は前スケジュールの値から指数関数的に変化指数関数の形状は第4列に指定する値により変更可能
第4列が未指定の場合、値5.0を使用
 S: 温度は前スケジュールの値からシグモイド関数的に変化
シグモイド関数の形状は第4列に指定する値により変更可能
第4列が未指定の場合、値3.0を使用
C˽1001˽0300.0

最初のステップを300 Kに指定
INITTEMPオプションの値を上書き

L˽1000˽1000.0

1000ステップかけて300 Kから1000 Kへ線形に昇温C˽1000˽1000.0

1000ステップの間1000 Kを維持 (平衡化)E˽1000˽0500.0˽3.0

1000ステップかけて1000 Kから500 Kへ減衰因子3.0を用いて指数関数的に降温C˽1000˽0500.0

1000ステップの間500 Kを維持 (平衡化)S˽1000˽0300.0˽2.0

1000ステップかけて 500 Kから 300 Kへゲイン 2.0を用いてシグモイド関数的に降温

60

目次

メタダイナミクスにおけるCVの設定
“metacv.dat”
MD=(METAD=TRUE),

MD=(METADYNAMICS=TRUE),

MD=(WTMETAD=TRUE), MD=(WTMETADYNAMICS=TRUE)が指定された場合に必要
第1行では、第1列にCVの種類を指定した後にガウス関数の幅と付随するオプションを指定する。自由エネルギー面のデータを計算中に出力する
場合 (MDキーワードにおけるMETAPRINTFES=TRUE) 、グリッド情報
を下限、上限、間隔の順番で3列追加指定する (単位はガウス関数の幅と同じである) 。いくつかのCVでは、原子集団のインデックスを第2, 3行に指定する。複数のCVを指定する場合、2番目以降のCVは1つ前の指定の次行から指定する。
実装済みCVの種類と原子数10の場合の指定例を以下に示す。
 ABSOLUTEPOSITION
原子あるいは原子集団の中心の絶対位置
周期系では、METAPBCMIC=TRUEの指定によらずminimum image
conventionは不適用
ABSOLUTEPOSITION˽0.1˽1˽Z˽LABELS˽COC

原子ラベル5のz座標をCVに指定

5

第1行第2–6列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 原子集団の原子数
 第4列: CVとする座標方向 (X/Y/Z)
 第5列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番
 第6列: 原子集団の中心の指定方法
COC ··· 座標中心
COM ··· 重心
61

目次 BONDDISTANCE
2原子で定義される距離 (原子1–原子2)
原子ラベル1–2の距離をCVに指定

BONDDISTANCE˽0.1˽1˽2

第2–4列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 原子1の原子ラベル
 第4列: 原子2の原子ラベル
 BONDANGLE
3原子で定義される角度 (原子1–原子2–原子3)
原子ラベル1–2–3の角度をCVに指定

BONDANGLE˽0.1˽1˽2˽3

第2–5列は以下の値を指定
 第2列: ガウス関数の幅 (deg) 第3列: 原子1の原子ラベル
 第4列: 原子2の原子ラベル
 第5列: 原子3の原子ラベル
 BONDDIHEDRAL
4原子で定義される二面角 (原子1–原子2–原子3–原子4)BONDDIHEDRAL˽0.1˽1˽2˽3˽4

原子ラベル 1–2–3–4 の二面角を CV に指定

第2–6列は以下の値を指定
 第2列: ガウス関数の幅 (deg) 第3列: 原子1の原子ラベル
 第4列: 原子2の原子ラベル
 第5列: 原子3の原子ラベル
 第6列: 原子4の原子ラベル

62

目次

 BONDDISTANCEDIFFERENCE
2原子で定義される距離の差

BONDDISTANCEDIFFERENCE˽0.1˽1˽2˽3˽4

距離r12 − r34をCVに指定

第2–6列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 原子1の原子ラベル
 第4列: 原子2の原子ラベル
 第5列: 原子3の原子ラベル
 第6列: 原子4の原子ラベル
 BONDDISTANCEADDITION
2原子で定義される距離の和

BONDDISTANCEADDITION˽0.1˽1˽2˽3˽4

第2–6列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 原子1の原子ラベル
 第4列: 原子2の原子ラベル
 第5列: 原子3の原子ラベル
 第6列: 原子4の原子ラベル63

距離r12 + r34をCVに指定

目次

 MINIMUMDISTANCE
以下の式で定義される原子間距離の最小値

MINIMUMDISTANCE˽0.1˽3˽2˽500˽LABELS1˽2˽3

原子集団1は原子ラベル1, 2, 3の3原子から構成

4˽5

原子集団2は原子ラベル4, 5の2原子から構成

第1行第2–6列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 原子集団1の原子数
 第4列: 原子集団2の原子数
 第5列: 正の実数パラメータ
 第6列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番

64

目次

 MEANDISTANCE
以下の式で定義される2つの原子集団間の平均距離

MEANDISTANCE˽0.1˽3˽2˽6˽LABELS1˽2˽3

原子集団1は原子ラベル1, 2, 3の3原子から構成

4˽5

原子集団2は原子ラベル4, 5の2原子から構成

第1行第2–6列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 原子集団1の原子数
 第4列: 原子集団2の原子数
 第5列: 正の整数パラメータ
 第6列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番

65

目次

 ATOMPOINTPLANEDISTANCE
1原子で定義される点 (原子1) と3原子で定義される平面 (原子2,原子3, 原子4) 間の距離
周期系では、METAPBCMIC=TRUEの指定によらずminimum image
conventionは不適用
ATOMPOINTPLANEDISTANCE˽0.1˽1˽2˽3˽4

原子ラベル 1 の点と原子
ラベル2, 3, 4の平面間の距離をCVに指定

第2–6列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 点を指定する原子1の原子ラベル
 第4列: 平面を指定する原子2の原子ラベル
 第5列: 平面を指定する原子3の原子ラベル
 第6列: 平面を指定する原子4の原子ラベル

66

目次

 ATOMPLANEPLANEANGLE
3原子で定義される平面 (原子1, 原子2, 原子3と原子4, 原子5, 原子6) 間の角度
周期系では、METAPBCMIC=TRUEの指定によらずminimum image
conventionは不適用
ATOMPLANEPLANEANGLE˽0.1˽1˽2˽3˽4˽5˽6

原子ラベル1, 2, 3と4, 5, 6の平面間の角度を CV に指定

第2–8列は以下の値を指定
 第2列: ガウス関数の幅 (deg) 第3列: 1番目の平面を指定する原子1の原子ラベル
 第4列: 1番目の平面を指定する原子2の原子ラベル
 第5列: 1番目の平面を指定する原子3の原子ラベル
 第6列: 2番目の平面を指定する原子4の原子ラベル
 第7列: 2番目の平面を指定する原子5の原子ラベル
 第8列: 2番目の平面を指定する原子6の原子ラベル

67

目次

 MIDPOINTDISTANCE
原子集団の中心間距離
MIDPOINTDISTANCE˽0.1˽3˽3˽LABELS˽COC
1˽2˽3

1番目の点は原子ラベル1, 2, 3の3原子を用いて計算

4˽5˽6

2番目の点は原子ラベル4, 5, 6の3原子を用いて計算

第1行第2–6列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 1番目の点の指定に用いる原子数 第4列: 2番目の点の指定に用いる原子数 第5列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番
 第6列: 原子集団の中心の指定方法
COC ··· 座標中心
COM ··· 重心

68

目次

 GYRATIONRADIUS
原子集団の断面回転半径
GYRATIONRADIUS˽0.1˽10˽LABELS˽COC

10原子全てを用いて断面回転半径を計算

第2–5列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 原子集団の原子数
全原子を用いる場合、第2行の指定は不要
 第4列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番
 第5列: 原子集団の中心の指定方法
COC ··· 座標中心、以下の式で定義される非質量荷重断面回転半径を計算 ( は原子座標)

COM ··· 重心、以下の式で定義される質量荷重断面回転半径を計
算 ( は原子質量)

69

目次

 INERTIAMOMENT
原子集団の慣性モーメント
INERTIAMOMENT˽0.1˽10˽LABELS˽COC

10原子全てを用いて慣性モーメントを計算

第2–5列は以下の値を指定
 第2列: ガウス関数の幅
 第3列: 原子集団の原子数
全原子を用いる場合、第2行の指定は不要
 第4列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番
 第5列: 原子集団の中心の指定方法
COC ··· 座標中心、以下の式で定義される非質量荷重慣性モーメ
ントを計算 ( は原子座標)

ガウス関数の幅はÅ2単位
COM ··· 重心、以下の式で定義される質量荷重慣性モーメントを計算 ( は原子質量)

ガウス関数の幅はAMU∙Å2単位

70

目次

 RATIONALCOORDINATIONNUMBER
以下の式で定義される有利関数形の配位数

RATIONALCOORDINATIONNUMBER˽0.1˽3˽2˽6˽12˽2.0˽LABELS˽AVERAGE
1˽2˽3

原子集団1は原子ラベル1, 2, 3の3原子から構成

4˽5

原子集団2は原子ラベル4, 5の2原子から構成

第1行第2–9列は以下の値を指定
 第2列: ガウス関数の幅
配位数は無次元量
 第3列: 原子集団1の原子数
 第4列: 原子集団2の原子数
 第5列: 正の整数パラメータ
 第6列: 正の整数パラメータ
 第7列: 正の実数パラメータ(Å)

 第8列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番
 第9列: 規格化因子
SUM ···AVERAGE ···

71

目次

 FERMICOORDINATIONNUMBER
以下の式で定義されるFermi関数形の配位数

FERMICOORDINATIONNUMBER˽0.1˽3˽2˽8.0˽2.0˽LABELS˽AVERAGE
1˽2˽3

原子集団1は原子ラベル1, 2, 3の3原子から構成

4˽5

原子集団2は原子ラベル4, 5の2原子から構成

第1行第2–8列は以下の値を指定
 第2列: ガウス関数の幅
配位数は無次元量
 第3列: 原子集団1の原子数
 第4列: 原子集団2の原子数
 第5列: 正の実数パラメータ
 第6列: 正の実数パラメータ(Å)

 第7列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番
 第8列: 規格化因子
SUM ···AVERAGE ···

72

目次

 RMSD
原子集団の平均二乗偏差 (RMSD)
RMSD˽0.1˽10˽LABELS˽COC˽FIT˽ref.xyz

ref.xyz に含まれる構造に対するRMSDをCVに指定

第2–7列は以下の値を指定
 第2列: ガウス関数の幅 (Å)
 第3列: 原子集団の原子数
全原子を用いる場合、第2行の指定は不要
 第4列: 原子の指定方法
LABELS ··· 原子ラベル
SPECIES ··· Parameterセクションにおける原子の種類の指定順番
 第5列: 原子集団の中心の指定方法
COC ··· 座標中心、以下の式で定義される非質量荷重RMSDを計
算 ( は原子座標)

COM ··· 重心、以下の式で定義される質量荷重RMSDを計算 ( は原子質量)

 第6列: RMSDの計算方法
AS-IS ··· 原子座標をそのまま使用
FIT ··· 現在のステップの原子座標の参照する原子座標への重ね合わせクラスターモデルのみ使用可能 第7列: 参照構造を含むxyzファイル名

73

目次

メタダイナミクスにおけるバイアスポテンシャルの読み込み“bias.dat”
MD=(METAREADBIAS=TRUE)が指定された場合に必要
前回と今回のプログラム実行でMDキーワードにおけるMETAD=TRUE,METADYNAMICS=TRUE,WTMETAD=TRUE,

あ

る

い

は

WTMETADYNAMICS=TRUEの指定が必要
ファイル形式は、メタダイナミクス計算におけるガウス関数堆積の履歴を含む “biaspot” ファイルと同じである (3.3.4節も参照) 。
以下は、距離をCVとする場合の出力例である。
GAUSSIAN˽BIAS˽POTENTIAL:˽˽˽˽˽˽˽˽˽1***˽AT˽T=˽˽˽˽˽˽˽˽˽10.00˽FSEC,˽THIS˽RUN’S˽STEP˽NO.=˽˽˽˽˽˽20
˽˽˽˽Gaussian˽height˽˽˽˽˽=˽˽˽˽˽˽˽˽0.0040000000˽a.u.
˽˽˽˽Collective˽variable˽=˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽˽1˽˽˽˽˽˽Coordinate˽˽˽˽˽˽˽˽=˽˽˽˽˽˽˽˽1.6146539213˽Angstrom˽˽˽˽˽˽Gaussian˽width˽˽˽˽=˽˽˽˽˽˽˽˽0.1000000000˽Angstrom…

以下の点に注意 ガウス関数のインデックス (第1行) はMETAMAXGAUSSオプショ
ンで指定した値を超えないこと
 シミュレーション時間やステップ数 (第2行) はファイル読み込み時に無視される情報
 ガウス関数の高さ (第3行) は普通のメタダイナミクスでは一定値(METAHEIGHTオプションで指定した値) 、well-temperedメタダイ
ナミクスでは値が徐々に減少
 “metacv.dat” ファイルで複数のCVが指定された場合、CVのインデックスと値、ガウス関数の幅 (第4–6行) は繰り返し出力

74

目次

メタダイナミクスにおけるポテンシャル障壁の設定
MD=(METAWALL=TRUE)が指定された場合に必要
以下の偶関数多項式をCVのサンプリング範囲を拘束するポテンシャル障壁として適用可能

式中のパラメータは、ポテンシャル障壁を課す方向とともに “metacv.dat” ファイルで指定する。以下は、2次元のメタダイナミクス計算において1番目のCVの値が2.0より大きくなる場合に2次のポテンシャル障壁を適用する例である。
<Description of first CV>

詳細はメタダイナミクスにおけるCVの設定を参照<Description of second CV>

詳細はメタダイナミクスにおけるCVの設定を参照



CVとポテンシャル障壁の指定を分離する空行

U˽1˽1.0˽2.0˽2

ポテンシャル障壁に関する記述

最終行第1–5列は以下の値を指定


第1列: ポテンシャル障壁を課す方向
L ··· 下限
U ··· 上限



第2列: CVのインデックス



第3列: 係数



第4列: CVの上限・下限値(kcal/(mol∙Å2))ガウス関数の幅と同じ単位

第5列: 正の偶数パラメータ

複数のポテンシャル障壁を課す場合、2番目以降のポテンシャル障壁に関する記
述は1つ前の指定の次行から指定する。

75

目次

2 DCDFTBMDの実行
インプットファイル “dftb.inp” を作成したら、パラメータファイルや補助イン
プットファイルといった計算の実行に必要なファイルを適切なディレクトリに格納する。例えば、1節に示した水分子に対するSCC-DFTBエネルギー一点計算では、4つのパラメータファイル (oo.spl, oh.spl, ho.spl, hh.spl) をインプットフ
ァイルと同じディレクトリに準備する。実行バイナリファイルは引数なしで使用する。
Serial版<installation path>/dftb_serial.00.x

OpenMP版<installation path>/dftb_openmp.00.x

MPI/OpenMP版

mpiexec˽–np˽(numberofprocesses)˽<installation

path>/dftb_mpiomp.00.x

インプットファイルの拡張子部分以外のファイル名は、第1引数を用いて任意に指定できる。以下は、 “file.inp” ファイルを用いる例である。
Serial版<installation path>/dftb_serial.00.x file

OpenMP版<installation path>/dftb_openmp.00.x file

MPI/OpenMP版

mpiexec˽–np˽(numberofprocesses)˽<installation

path>/dftb_mpiomp.00.x file

OpenMP並列計算では、使用するスレッド数を計算の実行前に設定する。
bash

export˽OMP_NUM_THREADS=(number of threads)

csh, tcsh

setenv˽OMP_NUM_THREADS˽(number of threads)

大規模計算が止まる場合、スタックサイズの増加により解決する場合がある。
bash

ulimit˽–s˽unlimited
export˽OMP_STACKSIZE=(stack size value)

csh, tcsh

unlimit
setenv˽OMP_STACKSIZE˽(stack size value)

Serial, OpenMP版における基底関数の数の上限は16384である。

76

目次

3 DCDFTBMDの出力計算が正常に終了した場合、標準アウトプットファイルと詳細データファイル
が通常出力される。入力に “dftb.inp” ファイルを用いた場合、標準アウトプッ
ト フ ァ イルと詳細 デ ータファイルのファ イル名はそれぞれ “dftb.out” と“dftb.dat” となる。拡張子部分以外のファイル名を変更した場合 (例: 2節の“file.inp” ) 、標準アウトプットファイルと詳細データファイルのファイル名も適切に変更される (例: “file.out” と “file.dat”) 。以下では、1節に示した水分子に対するSCC-DFTBエネルギー一点計算のserial, OpenMP版の標準アウトプット
ファイルと詳細データファイルを概説する。部分系の情報やMulliken電荷の出力形式が多少異なるものの、MPI/OpenMP版の対応するファイルはserial, OpenMP版の出力に類似する。なお、改行 () と空白領域 (˽) の記号は今後省略される。

3.1 標準アウトプットファイル
標準アウトプットファイルは、プログラムのヘッダーから始まる。バージョン番
号は公開中の実行バイナリファイルに依存する。プログラム実行開始時刻の出
力後、Titleセクションに記述されたコメント行が表示される (1.2節を参照) 。
*********************************************************************************************************************
DDDDDD

CCCCCCDDDDDD

FFFFFFFF TTTTTTTTT BBBBBBBDDDDDDDD

CCCCCCCC DDDDDDDD FFFFFFFF TTTTTTTTT BBB

MMMM DDDDDBBB MMMMMM DDDDDDD

DDD

DDD CCCDDD

DDD FFTTTBBBBBB MM M M MM DDD

DDDDDD CCC

DDD

DDD

FFFFFFTTTBBBBBBBB MMDDDDDD CCC

DDD

DDD

FFFFFFTTTBBBMM DDDDDD

BBB MM

MM DDDDDD

BB MM

MM DDDDDDD

DDDDDDDD

CCCCCCCC DDDDDDDD

FF

TTT

BBB

DDDDDD

CCCCCCFFTTTBBBBBBBBDDDDDDMM

M

DDD

MM DDDDD

VERSION 2.0 (Dec 2019)*********************************************************************************************************************
Execution of DCDFTBMD begun %a %b %d %H:%M:%S %Y-------TITLE--------77目次

次に、パラメータファイルのパスと計算条件の要約が出力される。
2 kinds of elements will be used
Pair of element 1 and element 1 : {path}/oo.spl
Pair of element 1 and element 2 : {path}/oh.spl
Pair of element 2 and element 1 : {path}/ho.spl
Pair of element 2 and element 2 : {path}/hh.spl
SCC= True------------------------Max cycle= 200

Energy convergence= 1.0E-09

Density convergence= 1.0E-06

Charge mixing method= Broyden

Broyden iteration= 70

Mixing parameter= 2.0E-01

Eigensolver= DSYGVD…
DC= False------------------------DISP= False------------------------PBC= False------------------------OPT= False------------------------FREQ= False------------------------MD= False------------------------MISC= False-------------------------

78

目次

次に、系に依存する量が出力される。
Total number of basis set shells = 4
Number of basis functions=6

Number of electrons=8

Charge of system=0

Spin multiplicity=1

Number of occupied orbitals=4

Total number of atoms=3

原子数、電荷、スピン多重度はGeometryセクションで指定した値に対応する(1.4.1節を参照) 。原子価軌道数、基底関数の数、電子数は水分子に対して以下のように計算される。


原子価軌道数: O(s) + O(p) + H(s) + H(s) = 4



基底関数の数: O(s) + O(px) + O(py) + O(pz) + H(s) + H(s) = 6



電子数: O(s2p4) + H(s1) + H(s1) = 8

これらの量は使用するパラメータファイルに依存するため、minimal basisと価電子を考慮する化学的直感から異なった値となる可能性がある。
最後に、SCC-DFTBエネルギー一点計算の結果が出力される。
*** Start single point energy calculation ************************************** Start SCC-DFTB calculation ***********************************
Iter. Total energy

Energy diff.

Density diff.

Fermi level------------------------------------------------------------------------------------------------------1-0.0726464067-0.0726464067-2.0000000000-0.0025000000

2-4.0712418245-4.1438882313-0.0200187809-0.0025000000

3-4.0741988428-0.0029570183-0.0590590622-0.3564583333

4-4.0779379555-0.0037391127-0.0004103232-0.3525000000

5-4.0779381355-0.0000001800-0.0000024570-0.3525000000

6-4.0779381355-0.0000000000-0.0000000001-0.3562500000

Final SCC-DFTB Energy = -4.0779381355 Eh after 6 iterations
Execution of DCDFTBMD terminated normally %a %b %d %H:%M:%S %Y

SCC収束後の全エネルギーは−4.0779381355 Hartreeである。
79

目次

3.2 詳細データファイル
詳細データファイルは、入力構造に関連するプロパティ情報を含む。始めに、構築された部分系と局在化領域が出力される。この計算はDC=FALSEであるため、3原子全てが部分系1に属する。
SUBSYSTEM : 1------------------------1

2

3

LOCALIZATION REGION : 1----------------------------------------1

MPI/OpenMP版では、MISCキーワードにおいてPRINTLOCALREGION=TRUEの場合を除き簡潔な表示となる。
SUBSYSTEM--------------------111

次に、分子軌道のエネルギー (Hartree) と電子占有数が第1, 2列に出力される。
EIGENVALUES [Eh] AND OCCUPANCIES------------------------------------------------------------0.8489865909

2.0000000000-0.4143390336

2.0000000000-0.3137537414

2.0000000000-0.2591755518

2.0000000000-0.3992760730

0.0000000000-0.5583952222

0.0000000000

80

目次

MISCキーワードにおいてPRINTMO=TRUEの場合、分子軌道係数が追加出力される。占有軌道と仮想軌道はそれぞれOCC.とVIRT.で区別される。
EIGENVALUES [Eh], OCCUPANCIES, AND MO COEFFICIENTS-----------------------------------------------------------------------------------------1

OCC.
1

O

-0.8489865909

S

0 -0.8522740053

P-10.0017546133

P00.0000000000

P

10.0017546133

2

H

S

0 -0.1530727458

3

H

S

0 -0.1530727458

2

OCC.
1

O

-0.4143390336

S

0 -0.0000000000

P-10.4626963766

P

00.0000000000

P

1 -0.4626963766

2

H

S

0 -0.3908041098

3

H

S

0

2.0000000000

2.0000000000

0.3908041098…

次に、原子毎のMulliken電荷密度と正味電荷の値が第2, 3列に出力される。
MULLIKEN POPULATIONS AND NET ATOMIC CHARGES---------------------------------------------------------------------------------1

6.5926170634-0.5926170634

20.7036914683-0.2963085317

30.7036914683-0.2963085317

81

目次

SCCあるいはNCCキーワードにおいてORBITALPOP=TRUEの場合、原子軌道毎のMulliken電荷密度が追加出力される。
MULLIKEN POPULATIONS AND NET ATOMIC CHARGES---------------------------------------------------------------------------------1

2

O

6.5926170634

S

0 1.7342141199

P-1 1.4292014718

P

0 2.0000000000

P

1 1.4292014718

H

0.7036914683
S

3-0.5926170634-0.2963085317

0 0.7036914683

H

0.7036914683
S-0.2963085317

0 0.7036914683

MPI/OpenMP版では、原子価軌道毎のMulliken電荷密度が出力される。
MULLIKEN POPULATIONS--------------------------------------1.73421412

4.85840294

0.70369147

0.70369147

次に、Mulliken電荷を用いて計算された双極子モーメントのx, y, z成分とノルムがatomic unitとDebye単位で表示される。周期境界条件を課したモデルではこの情報は出力されない。
DIPOLE MOMENT [a.u.]----------------------------------xy0.4545263289

0.4545263289

NORM :

z0.0000000000

0.6427972988

DIPOLE MOMENT [Debye]--------------------------------------xyz

1.1552906835

1.1552906835

0.0000000000

NORM :

1.6338277531

82

目次

SCCあるいはNCCキーワードにおいてMAYER=TRUEの場合、原子価とMayer結合次数が行列の対角と非対角要素として表示される。
ATOMIC VALENCIES AND MAYER BOND ORDERS-------------------------------------------------------------------------1

2

3

1

1.796085

0.898042

0.898042

2

0.898042

0.912201

0.014159

3

0.898042

0.014159

0.912201

SCCあるいはNCCキーワードにおいてCM3=TRUEの場合、CM3電荷とそれを用いて計算された双極子モーメントが追加出力される。
CM3 NET ATOMIC CHARGES-------------------------------------------1-0.5833492665

2-0.2916746332

3-0.2916746332

CM3 DIPOLE MOMENT [a.u.]-----------------------------------------xyz

0.4474181000

0.4474181000

0.0000000000

NORM :

0.6327447451

CM3 DIPOLE MOMENT [Debye]---------------------------------------------xy1.1372233679
NORM :

1.1372233679

z0.0000000000

1.6082767104

83

目次

最 後 に 、 全 エ ネ ル ギ ー に 対 す る 各 成 分 の 寄 与 (Hartree) が 表 示 さ れ る 。
MPI/OpenMP版では、MISCキーワードにおけるPRINTATOME=TRUEの指定が対応する情報を取得するために必要である。
ENERGY CONTRIBUTION [Eh]--------------------------------------------ENERGY H0= -4.1689461692

ENERGY SCC= -0.0183616273

TOTAL ELECTRONIC ENERGY = -4.1505845420
REPULSIVE ENERGY= -0.0726464065

TOTAL ENERGY= -4.0779381355

MISC キ ー ワ ー ド に お い て PRINTATOME=TRUE の 場 合 、 全 エ ネ ル ギ ー 、electronicエネルギー、repulsiveエネルギーに対する各原子の寄与 (Hartree) が追加出力される。
ATOM RESOLVED ENERGY CONTRIBUTION [Eh]-----------------------------------------------------------------------TOTALELECTRONICREPULSIVE

1-3.3763986686-3.4127218718-0.0363232032

2-0.3507697335-0.3689313351-0.0181616016

3-0.3507697335-0.3689313351-0.0181616016

84

目次

3.3 その他の出力
その他の出力は、計算条件に依存して追加生成される。
3.3.1 DFTB計算における出力“chrgfile”
SCC=TRUEかつMISC=(WRITECHRGFILE=TRUE)が指定された場合に生成
Mulliken密度解析の結果をバイナリ形式で格納
SCCキーワードにおいてREADCHARGE=BINARYの場合、初期電荷設定に使用 (1.5.1節を参照)

以下のファイルは、serial, OpenMP版のみ生成
“zacaz”, “zbcbz”, “zzszz”
MISC=(USESCRATCH=TRUE)が指定された場合に生成局在化領域の軌道係数 (zacaz, zbcbzはα, βスピンに対応) と重なり行列(zzszz) のデータをバイナリ形式で格納
一時ファイルであり、解析不要

85

目次

3.3.2 振動数計算における出力“restart_hess”
FREQ=(PRINTHESS=TRUE)が指定された場合に生成
Hessian行列要素の計算結果をバイナリ形式で格納
FREQ/FREQUENCYキーワードにおいてRESTARTHESS=BINARYの場合、Hessian行列設定に使用 (1.5.1節を参照)

“hess.dat”
解析的あるいは数値的Hessian計算 (FREQ/FREQUENCYキーワードにおける

FREQTYPE=1

あ

る

い

FREQ=(INTERRUPTCP=TRUE)

はあFREQTYPE=2)

で

る

は

い

FREQ=(INTERRUPTSEMINUM=TRUE)が指定された場合に生成
Hessian行列要素のインデックスと値を出力 (形式は1.5.2節を参照)
FREQ/FREQUENCYキーワードにおいてRESTARTHESS=ASCIIの場合、
Hessian行列設定に使用 (1.5.2節を参照)
“dftb_vib.out”
FREQ=TRUEが指定された場合に生成
調和振動解析の結果をGaussian出力形式で格納

86

目次

3.3.3 MD計算における出力“restart”
MD=TRUEが指定された場合に生成MD計算の直近のステップに関する情報をバイナリ形式で格納
MDキーワードにおいてRESTART=TRUEの場合、計算の再開に使用(1.5.1節を参照)

“restart_chk”
MD=TRUEが指定された場合に生成MD計算の直近のチェックポイントステップに関する情報をバイナリ形
式で格納
MDキーワードにおいてRESTART=TRUEの場合、“restart” にファイル名をリネーム後計算の再開に使用可能 (1.5.1節を参照)

“traject”
MD=TRUEが指定された場合に生成MD計算における原子座標をステップ毎にxyzファイル形式で出力

“velocity”
MD=TRUEが指定された場合に生成MD計算における速度 (Å/fs) をステップ毎にxyzファイル形式で出力

87

目次



“mulliken”
MD=TRUEが指定された場合に生成MD計算におけるMulliken正味電荷をステップ毎に以下の形式で出力
 Serial, OpenMP版
 第1行: 原子数 (Natom)
 第2行: コメント (シミュレーション時間とステップ数)
 第3–(Natom + 2)行: 原子ラベル (第1列) 、原子の種類の記号 (第2列) 、原子のMulliken正味電荷 (第3列) MPI/OpenMP版
 第1行: 原子価軌道数 (Nshell, 第1列) とNatom (第2列) 第2行: コメント (シミュレーション時間とステップ数)
 第3–(Nshell + 2)行: 原子ラベル (第1列) 、原子の種類の記号 (第2列) 、原子価軌道の記号 (第3列) 、原子価軌道のMulliken正味電荷 (第4列)88目次

3.3.4 メタダイナミクスにおける出力
MD計算における出力 (3.3.3節を参照) からの追加出力を列挙する。
“restart_meta”
MD=(METAD=TRUE),

MD=(METADYNAMICS=TRUE),
あMD=(WTMETAD=TRUE),

る

い

は

MD=(WTMETADYNAMICS=TRUE)が指定された場合に生成
メタダイナミクス計算におけるガウス関数堆積の履歴をバイナリ形式
で格納
MDキーワードにおいてMETARESTART=TRUEの場合、計算の再開に使用 (1.5.1節を参照)
“restart_meta_chk”
MD=(METAD=TRUE),

MD=(METADYNAMICS=TRUE),
あMD=(WTMETAD=TRUE),

る

い

は

MD=(WTMETADYNAMICS=TRUE)が指定された場合に生成
メタダイナミクス計算の直近のチェックポイントステップでのガウス関数堆積の履歴をバイナリ形式で格納
MDキーワードにおいてMETARESTART=TRUEの場合、 “restart_meta”にファイル名をリネーム後計算の再開に使用可能 (1.5.1節を参照)
“biaspot”
MD=(METAD=TRUE),

MD=(METADYNAMICS=TRUE),
あMD=(WTMETAD=TRUE),

る

い

は

MD=(WTMETADYNAMICS=TRUE)が指定された場合に生成
メタダイナミクス計算におけるガウス関数堆積の履歴を出力 (形式は1.5.2節を参照)
MDキーワードにおいてMETAREADBIAS=TRUEの場合、 “bias.dat” にファイル名をリネーム後計算の再開に使用可能 (1.5.2節を参照)

89

目次



“fes.dat”
MD=(METAPRINTFES=TRUE)が指定された場合に生成
自由エネルギー面のデータを以下の形式で格納
 第1–(CV数)列: CVのグリッド点 (単位はガウス関数の幅と同じ)
 最終列: 自由エネルギー (Hartree)

90

目次

参考文献
[1] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, T. Frauenheim, S. Suhai, G.
Seifert, Phys. Rev. B 58, 7260 (1998).
[2] D. D. Johnson, Phys. Rev. B 38, 12807 (1988).
[3] V. Eyert, J. Comput. Phys. 124, 271 (1996).
[4] P. Pulay, J. Comput. Chem. 3, 556 (1982).
[5] J. J. P. Stewart, P. Császár, P. Pulay, J. Comput. Chem. 3, 227 (1982).
[6] Y. Yang, H. Yu, D. York, Q. Cui, M. Elstner, J. Phys. Chem. A 111, 10861 (2007).
[7] M. Gaus, Q. Cui, M. Elstner, J. Chem. Theory Comput. 7, 931 (2011).
[8] M. Elstner, J. Phys. Chem. A 111, 5614 (2007).
[9] Y. Yang, H. Yu, D. York, M. Elstner, Q. Cui, J. Chem. Theory Comput. 4, 2067 (2008).
[10] Y. Nishimura, H. Nakai, J. Comput. Chem. 39, 105 (2018).
[11] C. Köhler, G. Seifert, U. Gerstmann, M. Elstner, H. Overhof, T. Frauenheim, Phys. Chem.
Chem. Phys. 3, 5109 (2001).
[12] C. Köhler, G. Seifert, T. Frauenheim, Chem. Phys. 309, 23 (2005).
[13] B. Hourahine, S. Sanna, B. Aradi, C. Köhler, T. Niehaus, T. Frauenheim, J. Phys. Chem. A
111, 5671 (2007).
[14] A. Domínguez, T. A. Niehaus, T. Frauenheim, J. Phys. Chem. A 119, 3535 (2015).
[15] I. Mayer, Chem. Phys. Lett. 97, 270 (1983).
[16] J. A. Kalinowski, B. Lesyng, J. D. Thompson, C. J. Cramer, D. G. Truhlar, J. Phys. Chem. A
108, 2545 (2004).
[17] S. Kaminski, M. Gaus, M. Elstner, J. Phys. Chem. A 116, 11927 (2012).
[18] M. Gaus, H. Jin, D. Demapan, A. S. Christensen, P. Goyal, M. Elstner, Q. Cui, J. Chem. Theory
Comput. 11, 4205 (2015).
[19] B. Quentrec, C. Brot, J. Comput. Phys. 13, 430 (1973).
[20] W. Mattson, B. M. Rice, Comput. Phys. Commun. 119, 135 (1999).
[21] D. Porezag, T. Frauenheim, T. Köhler, G. Seifert, R. Kaschner, Phys. Rev. B 51, 12947 (1995).
[22] G. Seifert, D. Porezag, T. Frauenheim, Int. J. Quantum Chem. 58, 185 (1996).
[23] W. Yang, T.-S. Lee, J. Chem. Phys. 103, 5674 (1995).
[24] M. Kobayashi, H. Nakai, in Linear-Scaling Techniques in Computational Chemistry andPhysics, R. Zalesny, M. G. Papadopoulos, P. G. Mezey, J. Leszczynski, Eds.; Springer:
Dordrecht, pp. 97–127 (2011).
[25] B. Cordero, V. Gómez, A. E. Platero-Prats, M. Revés, J. Echeverría, E. Cremades, F. Barragán,
S. Alvarez, Dalton Trans. 21, 2832 (2008).

91

目次

[26] H. Nishizawa, Y. Nishimura, M. Kobayashi, S. Irle, H. Nakai, J. Comput. Chem. 37, 1983(2016).
[27] S. Grimme, J. Comput. Chem. 27, 1787 (2006).
[28] M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, E. Kaxiras, J. Chem. Phys. 114, 5149 (2001).
[29] L. Zhechkov, T. Heine, S. Patchkovskii, G. Seifert, H. A. Duarte, J. Chem. Theory Comput. 1,
841 (2005).
[30] A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard III, W. M. Skiff, J. Am. Chem. Soc.
114, 10024 (1992).
[31] S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys. 132, 154104 (2010).
[32] S. Grimme, S. Ehrlich, L. Goerigk, J. Comput. Chem. 32, 1456 (2011).
[33] J. Řezáč, P. Hobza, J. Chem. Theory Comput. 8, 141 (2012).
[34] V. M. Miriyala, J. Řezáč, J. Comput. Chem. 38, 688 (2017).
[35] H. Kim, J.-M. Choi, W. A. Goddard III, J. Phys. Chem. Lett. 3, 360 (2012).
[36] R. Petraglia, S. N. Steinmann, C. Corminboeuf, Int. J. Quantum Chem. 115, 1265 (2015).
[37] J. Řezáč, J. Chem. Theory Comput. 13, 4804 (2017).
[38] M. Kubillus, T. Kubař, M. Gaus, J. Řezáč, M. Elstner, J. Chem. Theory Comput. 11, 332 (2015).
[39] C. G. Lambert, T. A. Darden, J. A. Board Jr., J. Comput. Phys. 126, 274 (1996).
[40] H. Jónsson, G. Mills, K. W. Jacobsen, in Classical and Quantum Dynamics in Condensed
Phase Simulations, B. J. Berne, G. Ciccotti, D. F. Coker, Eds.; World Scientific: Singapore,pp. 385–404 (1998).
[41] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumbsch, Phys. Rev. Lett. 97, 170201(2006).
[42] T. Atsumi, H. Nakai, Chem. Phys. Lett. 490, 102 (2010).
[43] Y. Nishimura, H. Nakai, J. Comput. Chem. 40, 1538 (2019).
[44] W. C. Swope, H. C. Andersen, P. H. Berens, K. R. Wilson, J. Chem. Phys. 76, 637 (1982).
[45] L. V. Woodcock, Chem. Phys. Lett. 10, 257 (1971).
[46] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, M. L. Klein, Mol. Phys. 87, 1117 (1996).
[47] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, J. R. Haak, J. Chem.
Phys. 81, 3684 (1984).
[48] H. C. Andersen, J. Chem. Phys. 72, 2384 (1980).
[49] G. Marsaglia, A. Zaman, W. W. Tsang, Stat. Probab. Lett. 8, 35 (1990).
[50] H. C. Andersen, J. Comput. Phys. 52, 24 (1983).
[51] T. Atsumi, H. Nakai, J. Chem. Phys. 128, 094101 (2008).
[52] A. Laio, M. Parrinello, Proc. Natl. Acad. Sci. USA 99, 12562 (2002).
[53] A. Barducci, G. Bussi, M. Parrinello, Phys. Rev. Lett. 100, 020603 (2008).
[54] http://www.dftbplus.org/fileadmin/DFTBPLUS/public/dftbplus/latest/manual.pdf

92