ASEでMOPACを動かす

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

通常のやり方
1. エラーの内容
2. 回避策（失敗）
無理やり動かす
計算実行
まとめ

通常のやり方

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

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

エラーの内容

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」を探しますが、それがないのでエラーになります。

追記 2024.6.20
ASEのアップデートによりエラーは解消しました。

回避策（失敗）

それなら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を無理やり動かしましたが、ASEのバージョンアップによってこんな面倒なことをやらなくてもおそらく動くようになっていると思います。

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