GrimmeのxTBを利用して、いわゆるScan (reluxed scan)をWindows上のGUIで実行・解析するツールです。1次元、2次元のscanと複数パラメータのconcerted scanができます。
- 2D Scanの鞍点判定をplot時に基準を変えて再実行できるように変更。一部処理を並列化して少しだけ速くしました。
- 2D Scanの鞍点判定をscipyのspline fitを利用する方法に変更。やや処理が重くなり、以前よりまともな判定がでるように。
- リファクタリング + GUIデザインの変更
- cclibを利用した読み込みを追加
- xyzファイルをインプットとして読んだとき、最後の構造を利用するように動作を変更。
- とりあえず完成
- Windows 11
- Python 3.8
- wxpython 4.1.1
- numpy 1.22.2
- scipy 1.8.0
- cclib 1.7.2
- matplotlib 3.5.1
- xtb-6.6.1 for Windows
python -m pip install wxpython cclib matplotlib scipy
で雑に動くと思います。Pythonが古いとエラーがでるかもしれません。
https://github.com/grimme-lab/xtb の右側のreleaseからWindows用のバイナリをダウンロードして適当な場所に置きます。ここでは D:\programs\xtb-6.6.1
においたものとします。
また適当な分子構造ビューワを用意します。
xtbscanフォルダの config.py
をテキストエディタで開き、上のほうにあるいくつかのパラメータを自分の環境に合わせて書き換えます。
VIEWER_PATH = 'D:/programs/jmol/jmol.bat'
XTB_BIN = 'D:/programs/xtb-6.6.1/bin/xtb.exe'
XTB_PARAM_DIR = 'D:/programs/xtb-6.6.1/share/xtb'
また下記をFalseにすると、2D Scanの鞍点の判定にscipyを使わなくなります (以前の適当な方法に戻す)。 その場合scipyはインストールしなくてもよいです。 その下のパラメータは鞍点判定時の厳しさ(小さいほど厳しい)のデフォルト値です。
# Use Scipy to check saddle point in 2d scan
USE_SCIPY = True
CHECK_SADDLE2D_GRAD_TOL = 0.001
xtbscan.pyw
を実行します。
- ファイルの読み込みはファイルをドラッグアンドドロップするか左上の ... ボタンから読み込みます。
- 初期構造ファイルは xyz、Gaussian Job/Log ファイル、もしくはcclib [https://cclib.github.io/] で読み込めるファイルが利用できます。
- 拡張子が gjf/com/log/out のものはGaussianファイルとして、xyz のものは xyzファイルとして処理します。それ以外はclibでそのまま読みこみます。
- 複数構造があるファイルなどの場合はいずれも最後の構造が利用されます。
- xyzファイル以外の場合、対応するxyzファイルが自動で作成されて初期構造ファイルとしてセットされます。
- 複数構造があるxyzファイルの場合、name_last.xyzというファイルが自動で作成されて初期構造ファイルとしてセットされます。
- 上記の場合、同名の古いファイルは上書きされるので注意してください。
- Methodは計算パラメータです。通常はgfn2でよいです。
- chargeは電荷、UHFは不対電子の数で閉殻なら0、ラジカルなら1 etc.です。Gaussianとかの多重度と指定がちょっと違うので注意。
- SolvationはALPBとGBSAがありますが、とりあえずやるならALPBでいいと思います(新しい方)。
- Forceはscanやconstrainの強さです。通常は1.0のままでいいです。完全にパラメータをfixせずに、強いポテンシャルで束縛してscanします。
- 2原子間の距離 (distance)、3原子の角度(angle)、4原子の二面角(dihedral)でscanできます。
- 右に条件を打ち込んでAddボタンを押すと表に追加されます。
- Atomsはスペースかカンマで区切って入力してください(1スタートの番号順です)。
- Start/End/Num Stepは、初期値、最終値、ステップ数です。
- 初期構造ファイルをあらかじめ読み込んでおくと、左上のviewボタンでviewerを開いて作業できるのでよいです。
- 構造を読み込んだ状態でcurrentを押すと現在値がStartに入ります。
- 条件1つだと1次元、2つだと2次元になります。3次元以上はできません。
- concerted mode のチェックを入れるとconcerted scanになります。この場合は複数のパラメータを同期的に変化させて1次元のscanをおこないます。この場合は3つ以上のパラメータの変化もできます。ただし全ての条件でNum Stepを同じ数にする必要があります。
- なおscanを1つも入れないで実行すると単純な構造最適化がおこなわれます。
- Scan以外のパラメータ(距離・角度・二面角)を固定できます。
- 現在値で固定するときは、currentで取得してもいいですが、auto としてもよいです(これはscanのほうではでききません)。
- cpus と memory per cpu を設定して Runを押すと実行されます。
- cpus はパソコンのCPUの物理コア数以下を指定すると良いです (4コア8スレッドなら最大4など)。
- メモリは最高でこの入力 x 利用CPUコア数 になるので、適当に入れてください。
- keep logはxtbの中間ログファイルを残すかどうかです。when failは、エラーがでて失敗したときだけ残します。
- Runを押すと計算がはじまります。実行中にStopを押すとキャンセルされます。途中から続行はできませんので、完全に中止になります。
- 計算が正常に終わるとメッセージボックスが出てきて、結果をロードするか聞かれます(scanでない最適化の場合は聞かれません)
- ロードすると、左下のCSVファイルのところにファイル名が読み込まれます。
- 計算結果はインプットファイル名を
JOB.xyz
とすると、JOB_xtbscan.csv
とJOB_xtbscan.xyz
に出力されます。エネルギーの情報はCSV、構造情報はXYZファイルに書かれています。 - 過去の計算結果も、CSVファイルをロードすれば読み込めます。CSVファイルをドラッグアンドドロップしてもよいです。
- CSVを読み込んだ上で、その下の各ボタンを押すと可視化ができます。
- view (all) は全ての構造のXYZファイルをビューワで開きます。
- show table はCSVの中身をテーブルで表示します。普通のCSVなのでエクセルで読んでも読めます。
- 特定の番号の構造を見るときは、その番号を打ち込んで、viewで見れます。copyを押すと座標情報がクリップボードにコピーされます。
- plot を押すとscanしたパラメータ/エネルギーの図が表示されます。1次元なら普通のグラフ、2次元なら3次元のグラフが表示されます。エネルギーは、最低値を0として相対値(kcal/mol)で表示されます。
- annotationにチェックを入れると、どの点がどの番号かも表示します。上記から番号を指定するときに便利です。
- annotationにチェックを入れていると、「saddleっぽい」点が赤くなります。1次元scanの場合は単に極大値です。2次元の場合、離散値からscipyを利用したspline fittingをおこない、その近似曲面上の一階/二階の微分からsaddle判定をおこないます(USE_SCIPY=Trueのとき)。勾配に関する判定を緩めにしてありますので、間違ったのも拾いますが、そこは人間が判断しましょう。
- 2次元scanの場合、annotationをチェックしてgrad. tol.を入れたうえでplotすると、そのgrad. tol.を基準にして鞍点の判定をやり直します(数秒~十数秒かかります)。赤い点が多すぎるor少なすぎるときは適宜基準を変えてください。そのときに判定された鞍点の番号はログに表示されます。csvに出力されている値は計算時のもので、これは更新しません(default値での判定)。
- 2次元scanの場合、散布図だけでなくsurfaceボタンを押すとPESのような表示もできます。
- なおgrad. tol.の値は物理的に意味のある単位ではなく、適当なスケールのものです。
- XYZファイルをドラッグアンドドロップするとInputファイルとして読み込まれます。結果をみたいときはCSVファイルをドラッグアンドドロップしてください。同名のXYZファイルも同時に読み込まれます。
- 結果のCSVファイルには、各scanパラメータについて、[scan] と [real] の値が書いてあります。[scan] は設定条件からの計算値、[real]は実際の構造から計算した値です。厳密に座標やパラメータを固定せずに強いポテンシャルによる束縛で固定をしているため、多少はずれます。 プロットするときは実際の値[real]が使われています。