Journal of Computer Chemistry, Japan
Online ISSN : 1347-3824
Print ISSN : 1347-1767
ISSN-L : 1347-1767
ソフトウェア・レビュー
DCDFTBMDプログラムの公開
西村 好史吉川 武司中井 浩巳
著者情報
ジャーナル フリー HTML

2018 年 17 巻 5 号 p. A21-A27

詳細
Abstract

The academic-free program DCDFTBMD has been released for performing large-scale quantum chemical calculations based on divide-and-conquer density functional tight-binding (DC-DFTB) method. The DC-DFTB method achieves low computational cost with linear complexity and its massively parallel implementation enables efficient geometry optimization, vibrational frequency analysis, and molecular dynamics simulations. The present article describes DCDFTBMD in terms of download and installation details, structure of input and output, and available features together with associated keywords and options.

1 はじめに

分割統治型密度汎関数強束縛(DC-DFTB)法は,大規模系の量子化学計算を合理的な計算時間で実行可能な手法の一つである.2000年代には,DC法を提案したYangらにより, DC-DFTB法のタンパク質,多糖類,水への適用が研究された [1,2,3].我々は,2011年よりDC-DFTB法に基づく計算を高効率に実行するための理論・プログラム開発を開始し,超並列処理の実装 [4, 5]により数千~数万原子を含む系の全てを量子的に取り扱う分子動力学(MD)計算への展開を可能とした.さらに,開発したDC-DFTB-MD法を凝縮系に適用し,化学反応解析に対する有効性を精査してきた [6,7,8,9,10,11,12,13,14,15].そして,2018年11月よりDCDFTBMD というプログラム名でアカデミックフリーとしての公開を開始した [16].本稿では,DCDFTBMD の公開に関する概要,計算の実行における入出力,およびプログラムに実装されている代表的な機能について報告する.

2 DCDFTBMD の入手と実行準備

アカデミック(官を含む)に所属する利用希望者は,DCDFTBMD を入手するために専用のサイト [16]で利用許諾規約を確認し,登録を行う必要がある.利用許諾規約には,(i)許諾事項,(ii)禁止事項,(iii)無保証,(iv)利用中止,(v)免責事項,(vi)引用,(vii)規約の変更に関する詳細が含まれている.ユーザー登録に必要な情報は,氏名・所属・e-mailである.審査後,ダウンロードサイトの情報が登録したe-mailに送付される.なお,一般の利用希望者は,分子モデリング・可視化ソフトウェアWinmostarと連結したプログラムをX-ability社から入手できる [17].

ダウンロードサイトから入手可能なパッケージは,利用許諾規約,インストーラ,実行バイナリ,実行用シェルスクリプト,日本語・英語マニュアル,サンプル入出力から構成されている.インストーラの起動により,実行環境の構築が以下の手順で行われる.まず,インストールする実行バイナリを1ノード逐次計算に対応するSerial版,ノード内並列計算が可能なOpenMP版,フラットMPIおよびハイブリッド並列計算が可能なMPI+OpenMP版の中から選択する.次に,MPI+OpenMP版を選択した場合,使用するMPIライブラリ(MPICH, OpenMPI, Intel MPI)とそのパスを設定する.最後に実行バイナリのインストールパスを入力する.インストーラの正常終了後,Linuxのsourceコマンドを用いた環境変数の読み込みを行うことで,プログラム実行に必要な準備が完了する.

3 DCDFTBMD の入出力

DCDFTBMD のファイル入出力処理をFigure 1に示す.計算条件や分子構造を指定する入力は,拡張子.inpで与える.ファイル名のデフォルトはdftbであるが,任意の文字列を用いることも可能である.プログラムは,入力とは別にDFTB計算に対する原子の種類のペア毎のパラメータを必要とする.汎用目的・特定目的のパラメータのいくつかは,DFTBのWebサイト [18]からダウンロード可能である.プログラムの実行により,計算条件の要約や分子構造に対するエネルギー値などの情報が拡張子.outを持つ標準出力に出力される.また,DC法における系の分割結果や各原子のMulliken電荷などの情報が拡張子.datを持つ詳細出力に出力される.これらのファイル名には,入力で用いた文字列がそのまま採用される.MD計算の場合,計算を再開するためのファイルや解析のための時系列データが追加で出力される.この他に,計算条件に応じて追加の出力ファイルが適宜作成される.

Figure 1.

 Input/Output structure of DCDFTBMD . The red characters for specifying input and output filenames can be modified arbitrarily.

DCDFTBMD の入力例をFigure 2に示す.入力は,計算条件を指定するkeyword section,計算についてのコメントを記述するtitle section,DFTB計算パラメータを指定するparameter section,分子構造を指定するgeometry sectionの4つに大別され,各sectionの終わりには空行が挿入される.DC計算で系の分割方法を手動指定する場合は,geometry sectionの後に各原子が所属する部分系の番号を入力する.

Figure 2.

 Sample input of DCDFTBMD for DC-DFTB1 single point calculation of carbon nanotube (CNT) system.

Keyword sectionで指定可能なキーワードをTable 1に示す.Entry #1–#5はDC-DFTB計算に関する設定,Entry #6–#8はシミュレーション手法に関する設定,Entry #9は出力制御に関する設定をそれぞれ担当する.このsectionにkeyword=TRUEの形式で入力すると,各キーワードをプログラムのデフォルト設定で使用する.キーワードが未指定の場合,Entry #1と#3のみが有効化され,DFTB2 [19]と呼ばれるレベルのDFTB計算にDC法を適用した一点計算を行う.Figure 2の例では,keyword=FALSEの形式でEntry #1の使用が無効化されている.この場合,Entry #2が有効化され,DFTB1 [20]という最も単純な近似に基づきDC-DFTBエネルギー計算がなされる.各キーワードには計算条件を細かく設定するためのオプションが複数用意されており,キーワードの等号の後に丸括弧を用いて指定する.代表的なオプションのいくつかは4章で述べる.

Table 1.  Available keywords of DCDFTBMD .
Entry Keyword Description
#1 SCC Setup for self-consistent charge (SCC) DFTB calculation
#2 NCC Setup for non-self-consistent charge (NCC) DFTB calculation
#3 DC Setup for DC calculation
#4 DISP Setup for noncovalent interaction correction
#5 PBC Setup for periodic boundary condition
#6 OPT Setup for geometry optimization
#7 FREQ Setup for vibrational frequency analysis
#8 MD Setup for molecular dynamics
#9 MISC Setup for degree of information in output

Title sectionに入力されたコメントは標準出力にそのまま出力される.複数行に渡った入力も可能である.

Parameter sectionでは,まず分子構造に含まれる原子の種類の数を入力する.次に,原子の種類を指定するラベルと原子価軌道のサイズを入力する.原子価軌道のサイズ指定には1, 2, 3という値を用い,順にs, p, d 軌道までを考慮することを意味する.ここで用いるべき値は,次の行で指定するDFTB計算パラメータの開発条件に依存する.DFTB計算パラメータファイル名の入力では,このsectionで指定する原子の種類の順番を厳守する必要がある.また,パラメータファイルを入力とは別のディレクトリに格納して使用する場合,参照パスを適切に指定する必要がある.1つ目の原子の種類について入力が完了したら,同様の入力を原子の種類の数だけ繰り返す.なお,計算条件によっては原子価軌道のサイズの後に追加パラメータを入力することがある.

Geometry sectionは,分子構造に含まれる原子数および電荷とスピン多重度の指定で開始される.分子構造は,Å単位のxyz形式で入力する.構造最適化計算やMD計算で特定の原子の座標を固定する制約を課したい場合や,振動数計算で同位体質量を指定したい場合は,z座標の後に対応する情報を入力する必要がある.

4 DCDFTBMD の機能とオプション

本章では,DCDFTBMD の機能についてTable 1に示した各キーワードのオプションと関連付けて示す.

4.1 DFTB計算手法とプロパティ解析の機能

DC-DFTB計算に関する設定のうち,Entry #1と#2は主にDFTB法の計算レベルやプロパティ解析の有無を制御するものであり,どちらかは必ず有効化する必要がある.3章で述べたように,Figure 2の入力例ではDFTB1という計算レベルが採用される.すなわち,SCC=FALSEあるいはNCC=TRUEを入力で指定すると,以下の式で与えられる全エネルギーを評価する.   

E DFTB 1 = E 0 + E 1 = A > B atom V A B + μ ν AO D μ ν H μ ν 0   (1)

式(1)のVABH0μνはDFTBパラメータを用いて計算される.D μνは密度行列であり,DC法を適用する場合は全系を分割して得られる部分系の密度行列の和で近似される.

SCC=TRUEによりEntry #1が有効化されている時の全エネルギーは,DFTB2レベルのものとなる.   

E DFTB 2 = E DFTB 1 + E 2
  
= E DFTB 1 + 1 2 A B atom γ A B Δ q A Δ q B   (2)

ここで,γΑΒは原子A, BのMulliken電荷ゆらぎΔqA, ΔqB 間の相互作用を原子間距離とDFTBパラメータを用いて記述する.DCDFTBMD では,SCCキーワードにおいてGAMMASPLINEというオプションの有効化により,γΑΒの短距離項の計算コストを削減できる [5].

DFTB2の計算精度を改善するため,電荷ゆらぎによって原子サイズが変化する効果を考慮するDFTB3 [21]が最近提案された.式(3)で表されるDFTB3のエネルギー計算は,THIRDFULLオプションの適用により許可される.   

E DFTB 3 = E DFTB 2 + E 3
  
= E DFTB 2 + 1 3 A B atom Γ A B Δ q A 2 Δ q B   (3)

後方互換性のため,ΓABの非対角成分を0として計算する手法(DFTB3-diag) [22, 23]も実装されている(オプション名THIRDDIAG).

SCC計算において,γΑΒの短距離項に指数関数を掛けて減衰挙動を修正することで水素結合の記述が改善される [22, 23].これは,オプション名DAMPXHの有効化により導入される.文献 [21]にあるDFTB3は,THIRDFULLとDAMPXHのオプションを同時に用いた手法であることに留意する必要がある.

SCC計算では,通常原子間のMulliken電荷ゆらぎの相互作用を評価するが,原子価軌道間の相互作用とすることで電荷ゆらぎの記述の柔軟性が増す.この拡張を実現するORSCCオプションは,DFTB2 [24]とDFTB3 [25]に対して実装されている.この時,式 (2), (3) にあるエネルギー成分E2, E3は以下のようになる.   

E 2 = 1 2 A B atom l A l B γ l A l B Δ q l A Δ q l B   (4)
  
E 3 = 1 3 A B atom l A l B Γ l A l B Δ q l A Δ q l B Δ q A   (5)

式(4), (5)のlAは原子Aの原子価軌道を表す.

開殻系分子のSCC計算において,スピン分極効果を取り込むためにはオプション名SPINが必要となる.スピン分極効果の全エネルギーへの寄与は   

E spin = 1 2 A atom l A l A ' W l A l A ' p l A p l A '   (6)
で与えられる [24,25,26].式(6)の W l A l A ' , p l A はそれぞれDFTBパラメータとMullikenスピン密度である.

SCC計算における外場効果として,点電荷 [27]と静電場 [28, 29]の指定が可能である.前者は,POINTCHARGEオプションの有効化に加えてNUMPCオプションで点電荷の数を入力する必要がある.点電荷の位置座標と電荷の値はgeometry sectionでラベルPCを用いて入力する.後者の計算では,EFIELDオプションを適用する.静電場の各成分はgeometry sectionでラベルEFを用いて入力する.

DCDFTBMD は,デフォルトでは修正Broyden法 [30]を用いてSCC計算の繰り返し手続きの収束性改善を図る.繰り返し計算の最大回数はMAXITERオプションより調節し,収束判定のための全エネルギーとMulliken電荷の閾値は,ECONV,DCONVオプションを通して制御する.

SCC計算において,エネルギーとMulliken電荷は常に算出されるプロパティである.また,孤立系の分子ではMulliken電荷に基づく双極子モーメントが計算される.双極子モーメントの計算精度は,Mayer結合次数計算 [31] (オプション名MAYER)の結果を用いたCM3電荷補正 [32, 33] (オプション名CM3)により改良を試みることができる.さらに,POLARオプションの指定により,有限電場法に基づいて静的分極率を求めることが可能である.

4.2 DC計算手法の機能

DC計算に関する設定は,Entry #3のDCキーワードを介して行う.全系をいくつかの部分系に分割する方法は,オプション名SUBTYPEを用いて以下の項目から選択する: 1原子(ATOM),入力を用いた手動指定(SPECIFY),全原子(SYSTEM),3次元グリッドを用いた自動決定(AUTO, AUTOXH),手動指定と自動決定を組み合わせた方法(SEMIAUTO, SEMIAUTOXH).自動決定法は,各原子の運動および化学反応によって部分系の定義が動的に変換するMD計算で特に重要である.3次元グリッドのサイズは,DELTARオプションを用いて調節できる.AUTOXHでは,水素原子を含む共有結合がグリッド分割によって2つの部分系にまたがらないよう調整する.

DC計算では,部分系の切り出しによる誤差を削減するため,部分系周囲の領域をバッファ領域として考慮する.DCDFTBMD におけるバッファ領域の定義は,部分系に含まれる各原子を中心とするオプション名BUFRADで指定した半径を持つ球状領域に含まれる全ての部分系である.大規模系でのバッファ領域決定を短時間で完了するため,MPI+OpenMP版にはセル分割法 [34]を用いたアルゴリズムが用意されている(オプション名DOMAIN).

4.3 非共有結合相互作用補正に関する機能

DFTB法の全エネルギーの定式化の観点から,弱い相互作用を適切に記述する補正項の導入が大規模系で重要となる.DCDFTBMD には種々の分散力・水素結合・ハロゲン結合補正手法が実装されており,Entry #4のDISPキーワードにあるDISPTYPEオプションから選択する.経験的分散力補正として,Slater-Kirkwood [35], Lennard-Jones [36], DFT-D2 [37], DFT-D3 [38, 39], DFT-ulg [40]の手法を適用可能である.分散力係数をMulliken電荷依存とするdDMC法 [41]の選択も可能である.D3H4 [42], D3H5 [43], D3X [44]と呼ばれる手法は,DFT-D3法に立脚して水素結合やハロゲン結合の補正項を追加する.

4.4 周期境界条件の設定に関する機能

DCDFTBMD は,3次元並進ベクトルをラベルTVとしてgeometry sectionの最終原子の座標の後に入力すると,周期系を取り扱うことができる.周期系のSCC計算では,Mulliken電荷ゆらぎの長距離相互作用を厳密に評価するためにEwald法 [45]を用いる.基本セルの形状が立方体の場合,多重極展開法に基づくアルゴリズム [5]を用いることでこの処理を大幅に高速化できる.アルゴリズムの選択は,Entry #5のPBCキーワード中のCOULOMBTYPEオプションより行う.この他に,応力テンソルの計算を明示的に指定するオプション(オプション名STRESS)などが整備されている.

4.5 シミュレーション手法に関する機能

Entry #6–#8では,それぞれ構造最適化計算,振動数計算,MD計算の条件が指定される.構造最適化計算のアルゴリズムは,OPTキーワードにおいてOPTTYPEオプションよりBFGS法,最急降下法,共役勾配法から選択できる.構造の収束は,エネルギー勾配の絶対値の最大と二乗平均平方根がGCONVオプションで指定した閾値を下回っているかにより判定される.周期系においてLATTICEOPTオプションを用いると,格子ベクトルの最適化が可能となる.

振動数計算では,FREQキーワードのFREQTYPEオプションよりHessian行列の計算を解析的あるいは数値的に行うかを決定する.SCC計算の場合,調和振動数に加えてIR強度とRaman活性が文献 [28, 29]の手続きに従って算出されるため,振動スペクトルを予測できる.熱化学解析はデフォルトで実行され,結果が標準出力に出力される.

MD計算の主な機能は,文献 [6]で既に報告しているため,MDキーワードの主要オプションについてのみ述べる.MD計算のステップ数と時間刻みは,オプション名NSTEPとDELTATより指定する.アンサンブルの選択はNVE, NVT, NPTオプションのいずれかを有効化することで行う.温度制御のアルゴリズムは,NVTTYPEオプションより選択し,熱浴の温度はBATHTEMPオプションで指定する.

5 まとめ

DC-DFTB法に基づく量子化学計算を実行するためのプログラムDCDFTBMD について,プログラムの入手と実行準備の方法,ファイル入出力の概要,およびプログラムに実装されている機能とオプションの一部を報告した.公開済のDCDFTBMD は,シンプルな入力形式の元,大規模系のプロパティ計算,構造最適化計算,振動数計算,およびMD計算に対応している.現在,開発版のプログラムにおいてメタダイナミクスによる自由エネルギー計算 [11]や時間依存DFTB法をDC計算に拡張した励起状態計算 [46]の導入が進んでいる.十分な検証を行った後,これらの機能についても公開版へ組み込むことを予定している.また,GPUアクセラレータに対応したプログラム開発にも着手している.DCDFTBMD のGPU化により,量子MD計算の時空間スケールのさらなる拡張が期待される.

謝辞

本研究は,文部科学省ポスト「京」重点課題5「エネルギーの高効率な創出,変換・貯蔵,利用の新規基盤技術の開発」,早稲田大学理工学術院総合研究所プロジェクト研究「計算化学の社会実装」(研究番号: 17P23),およびJSPS科研費基盤研究(S) JP18H05264の支援を受けて実施された.配布パッケージ作成でご協力いただいた高度情報科学技術研究機構(RIST)の河東田道夫博士に感謝する.

参考文献
 
© 2018 日本コンピュータ化学会
feedback
Top