ASE [https://wiki.fysik.dtu.dk/ase/] を利用してNudged Elastic Band (NEB)法の計算を実行するためのWindows上GUIツールです。 エネルギーと勾配計算の計算ソフトにはXTBとGaussain16が利用できます。
- とりあえず完成（テスト中）。
- Windows 11
- Python 3.8
- ase 3.23.0
- wxpython 4.1.1
- numpy 1.22.2
- matplotlib 3.5.1
- xtb-6.6.1 for Windows
- Gaussian 16
python -m pip install ase wxpython matplotlib
で雑に動くと思います。Pythonが古いとエラーがでるかもしれません。
- https://github.com/grimme-lab/xtb の右側のreleaseからWindows用のバイナリをダウンロードして適当な場所に置きます。ここでは
D:\programs\xtb-6.6.1においたものとします。 - 適当な分子構造ビューワを用意します (Jmolなど)。
- Gaussian for Windowsを利用する場合、インストールします。ここではデフォルトの
C:\G16W以下にインストールしたものとします。
フォルダの config.py をテキストエディタで開き、上のほうにあるいくつかのパラメータを自分の環境に合わせて書き換えます。
以下は利用する外部プログラムのパス等です。G16_SCRATCHはGaussianの一時ファイルの書き込み場所です。
XTB_BIN = 'D:/programs/xtb-6.6.1/bin/xtb.exe'
XTB_PARAM_DIR = 'D:/programs/xtb-6.6.1/share/xtb'
G16_ROOT = 'C:/G16W'
G16_SCRATCH_DIR = 'D:/scratch'
VIEWER_PATH = 'D:/programs/jmol/jmol.bat'
TEXT_EDITOR_PATH = 'C:/Windows/notepad.exe'
下記はASE（の裏で動いているnumpyの裏で動いているBLAS/LAPACK）の使用するCPUやメモリの設定です。
これはXTBやGaussianの設定とは別です。いじらなくても多分大丈夫ですが、万が一メモリ不足で落ちるときは ASE_OMP_STACKSIZE を増やしてください。
ASE_OMP_NUM_THREADS = None # integer or None (auto)
ASE_OMP_STACKSIZE = '256M' # memory per thread (without B)
aseneb_GUI.pyw を実行します。落ちるときは多分外部パッケージが足りていません。
- プロジェクトの作成（一プロジェクト一ディレクトリ推奨）
- 初期構造と最終構造の読み込みと構造最適化・エネルギー計算
- 内挿による初期反応パスの作成
- NEB計算による最適化（メイン）
- 上部のnewボタンからプロジェクトファイル（JSONファイル）を作成してください。このファイルの名前がプロジェクト名、おいた場所がプロジェクトのディレクトリになります。
- 既存のプロジェクトを読み込むときはopenから読み込んでください。
- 画面下部のXTB/G16部分が計算レベルの設定ですので、まずはここを設定してください。タブを選んでいる方が使われます。
- G16の場合、計算手法の設定はテンプレートファイルを利用します。これは通常のGaussianインプットからlink0行（%ではじまる行）を除外し、原子座標部分を @ に置き換えたものを用意します（最後の空行を忘れないように）。外部で作成したものをloadボタンから読み込むか、editボタンでテキストエディタで作成します（雛形が表示されます）。いずれの場合もプロジェクトディレクトリの中に
template.gjfという名前で保存されます。 - G16の右側の部分は通常は使いません（後述）。
- [Input Structures] から読み込みます。ASEで読み込めるファイルなら何でも読めます。複数構造が入っているファイルの場合は、最終構造を利用します。
- ...のボタンもしくはファイル名欄にドラッグ&ドロップでファイル名を読み込んだのち、loadボタンを押すと初期計算が走ります（計算レベルを設定しておく）。XTBは早いですが、G16でDFTならそれなりに掛かります。
- 計算が終わると右側に結果のファイル名（.traj; ASE専用のバイナリファイル）が表示されます。プロジェクトディレクトリには対応するxyzファイルも作成されます。
- viewを押すと構造が確認できます。
- #.imagesに作成する中間構造の数を入れて、Methodを選んで、Runを押します。
- #.imagesは反応経路の長さに合わせて8-20くらいがよいでしょう。Methodはidppが推奨です。
- パスの構造の数は初期・最終構造を含めて#.images+2になります。
- 注意点として、#.imagesを変更した場合、過去のファイル等の可視化もきちんとできなくなるので、ここを変更したときはここから全てやり直しになります。
- Runを押すと計算が走りますが通常は一瞬で終わり、右に結果ファイルが表示されます。
- viewから確認しましょう。
- 構造にはnode番号が振られ、node0 ~ node(#.images+1) の#.images+2個のnodeができます。
- 各種パラメータを設定して、Runを押します。
- Fmax（力の最大値）が小さくなるまで構造最適化が行われます。
- Methodは通常はasenebでよいと思いますが、よくわかりません。
- kは0.1程度が一般的のようです。
- Climbにチェックを入れると、Climbing Image NEB法になります。TS付近の精度が良くなるようです。
- OptimizerはFIREまたは LBFGS-LineSearchがよいらしいです。FIREは堅実に進むが収束が遅めで、Fmaxが小さくなってきた領域での収束が行ききらない場合がり、LBFGS-LineSearchは収束が速くFmaxの小さい領域にも強いが、Fmaxが大きい領域でぶれやすいという問題があります。とりあえずFIREで走らせて、さらに収束させたいときに、そこからLBFGS-LineSearchにする、などがあるようです。compositeを選ぶと、そのような挙動をします。
- Optimizerについては [https://docs.matlantis.com/atomistic-simulation-tutorial/ja/2_3_opt-algorithm.html] などが参考になります。
- Parallelを増やすと、各イメージについて並列でGaussianやXTBが走るようになります（後述）。
- 計算が進むと Calculation Logに経過が表示され、Result Filesに結果のファイルが表示されます。
- 計算を止めたか終わったかしたあと、もう一度Runを押すと、Result Filesの最後のファイルの結果からリスタートし、次のResultファイルが作成されます。やり直したいときはファイルを削除してください。
- Result Filesからファイルを選択し、その右のボタンで結果を確認できます。また対応するxyzファイルがプロジェクトディレクトリの中に作成されています。
- view/plot/infoは、そのファイルの中の最終パスの構造、エネルギー図、数値情報を表示します。
- view(all)/plot(all)/info(all) は最適化中の全てのものを表示します。
- save TSは最終パス中の最高エネルギーの構造をxyzで保存します。
- 主にOpen Shell Singletや3d遷移金属などで、普通にSCF計算をやると駄目なときようです。
- [3] Interpolate までやったあと、NEB計算をする前にG16タブの右側のテキスト欄に波動関数計算に必要な追加キーワードを入れたあとRunを押すと、各構造（node）について、このキーワードを追加した一点計算が走り、そのchkファイルが保存されます。
- G16による計算では、プロジェクトディレクトリ内にg16dataというディレクトリが作成され、その中にログやchkが出力されていきます。このとき既存のchkファイルがあるときは、そのファイルからguess=readで波動関数の初期値が再利用されます。
- そのため上記の方法で各nodeのchkを作ったあと、NEB計算を実行することで波動関数の初期値を変更できます。
- ここの計算でも[4]のParallelの並列化が反映されます。
- NEB計算では構造最適化1回ごとに各イメージのエネルギーと勾配の計算が実行されます。
- Parallelを2以上にすることで、Gaussian等が並列で起動され複数のイメージの計算が並列で実行できます。
- 要求されるCPUコア数は、Parallel * 計算あたりのCPU数（XTB/G16の設定）なので、利用できるCPUコア数を考えて適宜設定してください。
- XTBの場合はXTBの勾配計算よりも並列化のためのオーバヘッドが大きくなるのか、むしろ遅くなるようなのでParallel=1が推奨です。
- Gaussianの場合ある程度分散させたほうが早くなると思います。これはGaussian自体の並列化がサチってくるためです。特にCPUコア数がある場合は、G16のコア数を4や8まで減らしてParallelを増やすべきです。
- ただし手の空いたCPUがでると無駄なので、Parallelはイメージ数の約数にしたほうがよいです。例えば10イメージなら、Paralell=2 or 5とすると、1プロセスが計算するイメージが5つか2つになりますが、3にすると、プロセスAとBは3つのイメージを、プロセスCは4つのイメージを計算することになり、Cが4つ目を計算している間にAとBがが遊んでいることになり無駄です。
- ASENEB以外のMethodでは、初期・最終構造（最初と最後のnode）についてもエネルギーを計算するようです（ただこれは構造が変化しなければ、Gaussianでの計算時は前の結果を再利用するようにしているので計算コストは無視できます。XTBの場合はそもそも計算が速いので無視できます）。
- InitalとFinalのファイルを指定（loadは押さない）した状態で、Run Allを押すと全部の計算がまとめて流れます。
- 画面の状態がおかしくなったら（計算が走ってないのにloadやRunが押せないなど）、再起動するかメニューのEmergency > Resetをしてください。
- 色々デバッグが足りない部分があります。
- プロジェクトファイルは常時上書き保存されます。