量子化学ソフトウェアの導入 first example

密度汎関数法による計算化学ソフトウェアの導入を試みます．
このページで紹介するのは，ドキュメントの最初の例を読んで計算できたソフトウェアです．
計算した際の，入力ファイル，実行コマンド，および出力を記しています．
ソフトウェアの選択基準は，Debian でパッケージ化されており，ウェブサイトに実行例付きチュートリアルが置いてあり，環境設定しなくても計算できる，ということです．

このサイトでは，いくつかの化学ソフトウェアを開発しています．
作成した分子モデルの性質を調べるためには，例えば量子化学ソフトウェアの利用といった，計算科学的手法が必要です．
しかし，私は計算化学は専門外であり，使用頻度も多くありません．そのため，簡便さを重要な選択項目としています．

目次（ページ内リンク）

 候補リスト
 MPQC 2.3.1
 NWChem 7.0.2
 Molcas 22.10
 PSI4 1.3.2

候補リスト

まず，Debian パッケージで分子軌道に関係ありそうなソフトウェアをピックアップし，一覧を作成しました．
この段階では，密度汎関数法に限定していません（後から別目的で使うかもしれないので）．

密度汎関数法ソフトウェアの選択基準

簡単に試せそうな基準として，以下の項目を立ててみました．
Builcule の開発の視点からは，入力ファイルが簡単に作成できることが重要となります．

動かせたソフトウェア

筆者にも動かせたソフトウェアについては，ページ内リンクで情報や実行例を示しています．
これらが利用すべきソフトウェアの候補，というわけです．

 mpqc (2.3.1-22)：Massively Parallel Quantum Chemistry Program
 nwchem (7.0.2-4)：High-performance computational chemistry software (default MPI)
 openmolcas (22.10-1)：Quantum chemistry software package
 psi4 (1:1.3.2+dfsg-5)：Quantum Chemical Program Suite

MPQC 2.3.1

サイト：MPQC

(*) 入力ファイルの形式は 2 種類あると記されています．

以下の実例は，Debian の旧バージョンで作成しました．

ドキュメントの実例：Simple input

入力ファイル

ファイル名は mpqc.in 一択です． 


% B3LYP optimization of water
optimize: yes
method: KS (xc = B3LYP)
basis: 3-21G*
molecule: (angstrom)
    O    0.172   0.000   0.000
    H    0.745   0.000   0.754
    H    0.745   0.000  -0.754

実行と結果

出力ファイル名は -o オプションで指定します．指定しない場合は，結果が標準出力されます．
$ mpqc -o h2o.out
計算終了後，h2o.out，およびバイナリファイルがいくつか生成していました．
h2o.out 以外はバイナリファイルです．

ドキュメントの実例：Object-Oriented Input

入力ファイル


% This input does a Hartree-Fock calculation on water.
molecule<Molecule>: (
  symmetry = C2V
  unit = angstrom
  { atoms geometry } = {
    O     [     0.00000000     0.00000000     0.37000000 ]
    H     [     0.78000000     0.00000000    -0.18000000 ]
    H     [    -0.78000000     0.00000000    -0.18000000 ]
  }
)
basis<GaussianBasisSet>: (
  name = "STO-3G"
  molecule = $:molecule
)
mpqc: (
  mole<CLHF>: (
    molecule = $:molecule
    basis = $:basis
  )
)

実行と結果

上と同じく，入力ファイル名を mpqc.in とします．入力ファイルは -f オプションで指定します．
$ mpqc -f h2o.in -o h2o.out
とすると，h2o.out，およびバイナリファイルがいくつか生成しました．

NWChem 7.0.2

サイト：NWChem

以下で，いくつか実例を試してみます．解説部分は "Getting Started" からの改変した引用です．詳しくはそちらを参照してください．
これらは，Debian の旧バージョンで作成しました．

窒素の構造最適化：Dunning cc-pvdz 基底系を使用した窒素分子の SCF 構造最適化

入力ファイル

入力ファイルを引用します．読みやすくなるかと思い，空行を追加しています．
ファイル名は n2.nw とします．ちなみに，デフォルトのファイル名は "nwchem.nw" だそうです．

入力ファイルはフリーフォーマットで，ディレクティブと呼ばれるコマンドとそのコマンドを適用するデータ（基底系，原子とその座標など）で構成されています．
この例では，TITLE，GEOMETRY，BASIS，TASK の 4 つがディレクティブです．ドキュメントの一部を引用します．

"Getting Started" から Geometry，Basis，Task というページへのリンクが張られており，さらに詳細な情報が得られます． 


 title "Nitrogen cc-pvdz SCF geometry optimization"

 geometry
   n 0 0 0
   n 0 0 1.08
 end

 basis
   n library cc-pvdz
 end

 task scf optimize

実行と結果

結果は標準出力されます．
$ nwchem n2 | tee n2.out
とすれば，結果は標準出力されつつ，file.out というファイルにも出力できます．
他に，バイナリファイルがいろいろ生成しました．

なお，コンパイルのページに並列処理のコマンドを使うことができると記してあります．すなわち，
$ mpirun -np 6 nwchem n2 | tee n2.out
として，プロセス数を指定して実行することも可能です（ここでは 6）．

二重項状態の水：2 次の Møller-Plesset 摂動理論（MP2）を使用した H2O+ の最適化と振動周波数

ドキュメントの抜粋を再編して記します．
入力ファイルには，初期値の設定とそれに続く 3 段階のタスクが記述されています．

  1. 初期値の入力 == ランタイムデータベースの設定
  2. タスク 1：計算コストの低い基底系（STO-3G）を用いて，SCF の幾何学的最適化の予備的な計算を行います
  3. タスク 2：MP2 と分極関数を含む基底系を使用して最適化を終了します
  4. タスク 3：MP2 振動周波数を計算します

まず，初期値の入力部分．
プログラムはこのデータベースから情報を取得するので，異なるブロックや指令の順番を入れ替えても，結果には影響しないそうです．
print level が low に設定されているのは，最初の計算の冗長な出力を避けるためです．
座標，電荷，基底系に続いて，scf ...end ブロックで最初の SCF 計算が設定されています．

次いで，1 番めのタスク．
タスクでは，直前の指示文で設定されたパラメータを用いて，指定された計算をコードに実行させます．
タスクが完了しても，新しい入力ディレクティブで上書きされない限り，データベース内の設定は変更されずに後続のタスクで使用されます．
ここでは，下の内容が記述されています．

  1. キーワード scf と optimize で指定しているように，SCF 計算と構造最適化をおこないます
  2. 次のタスクのために，より適した基底系（6-31G**）を設定し直しています

次いで，2 番めのタスクでは，1 番めのタスクで設定した基底系を使って MP2 構造最適化をおこなっています．
3 番めのタスクではは振動周波数を計算しています．print 文で，SCFとMP2モジュールからの出力を無効化しています． 


start h2o_freq
charge 1
geometry units angstroms
 O       0.0  0.0  0.0
 H       0.0  0.0  1.0
 H       0.0  1.0  0.0
end

basis
 H library sto-3g
 O library sto-3g
end

scf
 uhf; doublet
 print low
end

title "H2O+ : STO-3G UHF geometry optimization"
task scf optimize
basis
 H library 6-31g**
 O library 6-31g**
end

title "H2O+ : 6-31g** UMP2 geometry optimization"
task mp2 optimize
mp2; print none; end
scf; print none; end

title "H2O+ : 6-31g** UMP2 frequencies"
task mp2 freq

実行と結果

入力ファイルの名前を h2o.nw とします．
拡張子を省略せずプロセス数を指定して実行し，標準出力をファイルに保存するなら，
$ mpirun -np 6 nwchem h2o.nw | tee h2o.out
とでもすれば，計算結果が h2o.out というファイルに出力できます．

量子化学ソフトウェアの導入 NWChemでは，NWChem と Gabedit を使って，アンモニアなどの簡単な分子を計算，表示します．

Molcas 22.10

サイト：Molcas / OpenMolcas · GitLab

水の計算

HTML ドキュメントの [3.1.6. Basic Examples] では，まず "3.1.6.1. Simple Calculation on Water" として水の計算例が紹介されています．
ここでは基底状態の計算のみ記していますが，マニュアルの本文では二重項状態と三重項状態の計算もしています．

入力ファイル (1)

まず，水分子のカルテシアン座標（XYZ 座標）を記述したファイルです．ファイル名を "water.xyz" とします． 


  3
Angstrom
  O  0.000000  0.000000  0.000000
  H  0.758602  0.000000  0.504284
  H  0.758602  0.000000  -0.504284

入力ファイル (2)

次は，計算方法を記述したファイルです．ファイル名を "water.input" とします．
カルテシアン座標を記述したファイルを置き換えるだけで，同じ条件で種々の分子の計算ができそうです．
なお，本文同ページの少し後に書いてありますが，coord=water.xyz 行を書き換えれば，一つのファイルにできます． 


&GATEWAY
  coord=water.xyz
  basis=sto-3g
&SEWARD
&SCF

入力ファイルの書式は，プログラムモジュールごとに記していく構造になっています．
& で始まる行が，計算をおこなう「プログラムモジュール」です．それに続く行でモジュールに動作を指示する「キーワード」を設定します．
ここでは，GATEWAY プログラムモジュールに続けて，coord キーワードで上に引用した水の座標 water.xyz を，basis キーワードで基底関数を設定しています．
SEWARD プログラムモジュールで積分を計算し，SCF プログラムモジュールで Hartree-Fock 波動関数の計算を完結させる，とのことです．

実行と結果

"water.xyz" と "water.input" を置いたディレクトリで，
$ pymolcas water.input -f
とすると計算できました．いくつかのファイルが生成しますが，そのうちの 2 ファイルを置きます．

構造最適化

続く "4.6.2 Geometry Optimization" を練習してみます．
新たに EMIL (Extended Molcas Input Language) コマンドという書式が加わります．
EMIL については，[7.3 General input structure. EMIL commands] に記されていますが，モジュールのループ実行，部分実行，変数の変更，ファイル出力などの操作を可能にするコマンドとのことです．

入力ファイル

上記入力ファイル (2) "water.input" を書き換えます．while 文が EMIL コマンドですが，構造が収束するまで繰り返し計算，ということです．
密度汎関数法（DFT/B3LYP）で構造最適化をしています． 


&GATEWAY
  coord=water.xyz
  basis=ANO-S-MB
>>  Do  While
  &SEWARD
  &SCF
    ksdft=b3lyp
  &SLAPAF
>>>  EndDo

実行と結果

[水の計算] で使った "water.xyz" と "water.input" を置いたディレクトリで，
$ pymolcas water.input -f
とすると計算できました．
いくつかのファイルが生成しますが，waterOpt.xyz が構造最適化後のカルテシアン座標のようです．
座標が少し変化しています． 


   3
  -76.354429660502660
 O     0.04257864     0.00000000     0.00000000
 H     0.71602336     0.00000000     0.84129652
 H     0.71602336     0.00000000    -0.84129652

構造最適化については，[5.2 Optimizing geometries] でもう少し詳しく紹介されています．

PSI4 1.3.2

PSI4 は計算負荷の高い部分が C++ で書かれ、C++ クラスのエクスポートやインターフェース部分は Python で書かれているようです．
したがって，PSI4 ｈは，実行ファイルとしても Python モジュールとしても使えるのですが，ここでは実行ファイルとして使います．

サイト：PsiCode
このページにマニュアル PSI4: Open-Source Quantum Chemistry へのリンクが張られています．

cc-pVDZ 基底系を使った，水分子の Hartree–Fock SCF 計算

入力ファイル

下は，"Sample Input Files" の最初のものです．
ファイル名をデフォルトの input.dat として保存します． 


# Any line starting with the # character is a comment line
#! Sample HF/cc-pVDZ H2O computation

memory 600 mb

molecule h2o {
  O
  H 1 0.96
  H 1 0.96 2 104.5
}

set basis cc-pVDZ
energy('scf')

チュートリアルには，
"The program correctly guessed all of these options for us. We can change the default behavior through additional keywords."
と書いてあります． この例から思うに，特に設定しなければ，電荷 0 の三重項状態と解釈されるようです．

実行と結果

$ si4 input.dat output.dat（このようにデフォルトのファイル名を使っている場合は，単に "psi4" だけでも可）
とすれば，結果は output.dat というファイルに出力されます．

電荷，スピン多重度，分子軌道の計算法の設定

水の例の次に，メチレンビラジカルを使い，標記の設定法が記されています．
操作は水の場合と同じなので，入力ファイルのみ引用します．行の途中にコメントを付けてみました． 


#! Sample UHF/6-31G** CH2 computation

molecule ch2 {
  0 3  #電荷とスピン多重度
  C
  H 1 R
  H 1 R 2 A

  R = 1.075
  A = 133.93
}

set basis 6-31G**
set reference uhf  #分子軌道の計算法
energy ('scf')

量子化学ソフトウェアの導入 PSI4では，PSI4 と Gabedit を使って，水などの簡単な分子を計算，表示します