ASEでMOPACを動かす

MOPAC

PythonのASE (Atomic Simulation Environment) でMOPACを動かそうとするとエラーが出ますが、ちょっと頑張れば動かすことができます。

スポンサーリンク

通常のやり方

ASE_MOPAC_COMMAND に以下の内容を設定します。

Windowsの場合
"C:\Program Files\MOPAC\bin\mopac.exe" PREFIX.mop > nul 2>&1
Linuxの場合
/path/to/MOPAC.exe PREFIX.mop 2> /dev/null
Windowsの場合

普通はこのように設定すると思いますが、エラーになってしまいます。

エラーの内容

CalculatorにMOPACを指定して実行すると以下のエラーが出ます。

  File "C:\Python311\Lib\site-packages\ase\calculators\calculator.py", line 709, in get_potential_energy
    energy = self.get_property('energy', atoms)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Python311\Lib\site-packages\ase\calculators\calculator.py", line 742, in get_property
    raise PropertyNotImplementedError('{} not present in this '
ase.calculators.calculator.PropertyNotImplementedError: energy not present in this calculation

energyが読めないとのことです。
実は現行のMOPACは出力の中に「TOTAL ENERGY」がありません。
以前のバージョンのMOPACにはそれがあったのですが、現在のMOPACにはなくなっています。
ASEは出力の中から「TOTAL ENERGY」を探しますが、それがないのでエラーになります。

回避策(失敗)

それならDISPキーワードを付ければいいじゃないか!
ということでやってみます。
Calculatorの設定を以下のようにして実行します。

atoms.calc = MOPAC(task="1SCF GRADIENTS DISP")

するとMOPACの入力ファイルは

PM7 1SCF GRADIENTS DISP RELSCF=0.0001 

のような形式で作られます。
計算結果は

          FINAL HEAT OF FORMATION =       -100.09375 KCAL/MOL =    -418.79224 KJ/MOL

          TOTAL ENERGY            =     -20544.23931 KCAL/MOL = ELECTRONIC ENERGY + CORE-CORE REPULSION
          ENERGY OF ATOMS         =      20444.29082 KCAL/MOL
                              SUM =        -99.94849 KCAL/MOL
          DISPERSION ENERGY       =         -0.14526 KCAL/MOL
                              SUM =       -100.09375 KCAL/MOL = FINAL HEAT OF FORMATION



          TOTAL ENERGY            =       -890.88253 EV
          ELECTRONIC ENERGY       =      -2476.20502 EV 

のようになり、TOTAL ENERGYが現れます。
これで行ける!と思ったのですがまたエラーになります。
今度のエラーは

  File "C:\Python311\Lib\site-packages\ase\calculators\mopac.py", line 195, in read_results
    self.final_hof = float(line.split()[5]) * kcal / mol
                     ^^^^^^^^^^^^^^^^^^^^^^
ValueError: could not convert string to float: 'FINAL'

ということで、FINAL HEAT OF FORMATIONが取得できないようです。
MOPACの出力をよく見ると、FINAL HEAT OF FORMATIONは1行目と7行目にもあります。
7行目のほうが邪魔なんだと思います。

スポンサーリンク

無理やり動かす

邪魔なのはここです。

          ENERGY OF ATOMS         =      20444.29082 KCAL/MOL
                              SUM =        -99.94849 KCAL/MOL
          DISPERSION ENERGY       =         -0.14526 KCAL/MOL
                              SUM =       -100.09375 KCAL/MOL = FINAL HEAT OF FORMATION

これを消せばいいので、以下の流れで消します。

ASEは環境変数ASE_MOPAC_COMMANDを実行するので、このコマンドの中に余計な行を消すpythonプログラムを入れます。
入れる方法は何通りかあると思いますが、今回はバッチファイルを書き、mopac実行後にpythonでmopac.outを編集する、という流れで行きます。

バッチファイル(C:\ASE\run_MOPAC_ASE.bat)

@echo off
"C:\Program Files\MOPAC\bin\mopac.exe" PREFIX.mop > nul 2>&1
python C:\ASE\ase_mopac_out.py

Python(C:\ASE\ase_mopac_out.py)

f_in = open("mopac.out", "r")
lines = f_in.readlines()
f_in.close()

f_out = open("mopac.out", "w") 

for line in lines:
    if not "= FINAL HEAT OF FORMATION" in line:
        f_out.write(line)
f_out.close()

環境変数

ASE_MOPAC_COMMAND
C:\ASE\run_MOPAC_ASE.bat
スポンサーリンク

計算実行

酢酸分子を構造最適化してみます。

from ase.build import molecule
from ase.calculators.mopac import MOPAC
from ase.optimize import BFGS

atoms = molecule("CH3COOH")
print(atoms.get_positions())

atoms.calc = MOPAC(task="1SCF GRADIENTS DISP")
optimized_geometry = BFGS(atoms)
optimized_geometry.run()
print(atoms.get_positions())

結果はこうなりました。

[[ 0.        0.15456   0.      ]
 [ 0.166384  1.360084  0.      ]
 [-1.236449 -0.415036  0.      ]
 [-1.867646  0.333582  0.      ]
 [ 1.073776 -0.892748  0.      ]
 [ 2.048189 -0.408135  0.      ]
 [ 0.968661 -1.528353  0.881747]
 [ 0.968661 -1.528353 -0.881747]]
      Step     Time          Energy         fmax
BFGS:    0 00:04:14     -890.882530        1.2871
BFGS:    1 00:04:14     -890.925590        1.2379
BFGS:    2 00:04:14     -890.942710        1.0117
BFGS:    3 00:04:14     -890.956200        0.4009
BFGS:    4 00:04:14     -890.965640        0.2121
BFGS:    5 00:04:14     -890.968700        0.2356
BFGS:    6 00:04:14     -890.970840        0.1548
BFGS:    7 00:04:14     -890.971870        0.0683
BFGS:    8 00:04:14     -890.972330        0.0494
[[ 2.20202201e-02  1.72937150e-01  2.37630992e-05]
 [ 1.62941112e-01  1.36768316e+00  1.24529589e-06]
 [-1.22645303e+00 -3.69728391e-01  3.66907622e-06]
 [-1.94716972e+00  3.09471135e-01 -1.11162937e-06]
 [ 1.06547033e+00 -8.84690234e-01  1.79371699e-05]
 [ 2.07486198e+00 -4.38231106e-01  2.56914400e-05]
 [ 9.84939253e-01 -1.54100530e+00  8.83078508e-01]
 [ 9.84965871e-01 -1.54083540e+00 -8.83149694e-01]]

エラーなく動きました。計算は非常に速いです。
出力からは、BFGSのサイクルが動いて、初期座標と計算後の座標が変わったことは分かると思います。

スポンサーリンク

まとめ

ASEでMOPACを無理やり動かしました。
非常に計算速度が速く、実験値ベースで作られたパラメータに基づいて計算をする半経験MOのほうが便利な場面もありますよね。
この記事がASEでMOPACを動かしたい人の参考になれば幸いです。

コメント

タイトルとURLをコピーしました