🧪

Gromacs with CP2K でQM/MM計算チュートリアル

2025/06/24に公開

はじめに

前回の記事でCP2KをリンクしたGROMACSのビルドができました。
そのテストのため、こちらのチュートリアルをやってみます。

インプットファイルや構造ファイルだけでなく、解説資料も用意されているチュートリアルです。

これのトップディレクトリにパワポによる解説もありますが、より詳細なコマンドもあるHandbook を参考にしてみます。

このチュートリアルには3つの例：

小分子のみQM計算で動かす例
スチルベンの2面角ポテンシャルを求める例
GFPタンパク質のスペクトル

があります。このうち、1. と 3. を今回やってみようと思います。

準備

まずこのチュートリアルリポジトリをクローンします。

git clone https://github.com/bioexcel/gromacs-2022-cp2k-tutorial.git

これで必要なPDBファイルやインプット用の mdp ファイルが全て揃います。

単分子

簡単な分子全体をQM領域とするパターンです。

系は全部でこの12原子だけです。

クローンしたリポジトリ中の nma ディレクトリに入ります。

cd /path/to/gromacs-2022-cp2k-tutorial/nma

まずはパラメータを作ります。

gmx_d pdb2gmx -f nma.pdb  # 1, 1 を選択

分子に設定する力場の種類と水の種類を聞かれるので、両方で１を入力します。
できた力場で、まずはエネルギー最小化します。

gmx_d grompp -f em.mdp -p topol.top -c conf.gro -o nma-em.tpr
gmx_d mdrun -s nma-em.tpr

続いてそこからNVTシミュレーションをやります。

gmx_d grompp -f md-nvt.mdp -p topol.top -c conf.gro -o nma-nvt.tpr
gmx_d mdrun -deffnm nma-nvt

MD自体はこれだけです。まあ応用するときはこのNVTインプットを参考にいろいろやっていく感じかなと。
とりあえずどんな感じの結果になったかを見ておきます。温度がどういう遷移をしたかを出力してみます。

echo 7 | gmx energy -f nma-nvt.edr
xmgrace energy.xvg

ハンドブックの結果と見た感じ異なりますが、300K 付近で揺れているというのはまあ再現できているかなと思います。

VMDでトラジェクトリを見るとこんな感じです。

緑色蛍光タンパク質の色素

緑色蛍光タンパク質（4eul）について、その中の色素分子をQM領域としてMDとスペクトル計算を行う例です。

なんとなく見返して気付いたのですが、この節ではGROMACSのコマンドが gmx だったり gmx_d だったりしますね。まあ editconf みたいな入出力がフォーマットファイルだったりするようなのはバイナリが単精度か倍精度かなんて関係ないと思いますし、 energy も読み込んだ edr ファイルが単精度か倍精度かで処理を変えて対応しているようです。一方 mdrun はさすがに CP2Kとリンクした倍精度版でないといけないでしょう。
ということで、どっちでもよさそうなやつは gmx で、倍精度版でやっておいたよさそうなやつは gmx_d としました。

Tutorial-Handbookからの変更点

QM領域を CRO 残基全体としました。チュートリアル通りに指定すると、 grompp でQM領域の電荷が -0.23 なんとかになっててインプットファイルで指定された電荷－1 と齟齬があるというWarningが出るためです。チュートリアルの指定原子は発光するチロシンの六員環と変形してできた5員環、またその橋掛けのメチレンと6員環についた酸素原子だけというほぼ最小構成でした。チュートリアルの原子指定が電荷が等しくならないのは意図があってかどうかわかりませんが、この記事ではGromacsのコマンドの整合性を優先しました。
スペクトル計算時の NVT 接続間隔を 10ステップごとから毎ステップにしました。エネルギーが階段状になるのがなんとなく嫌だったので。
TDDFTのnstate数を 5 から 20 に増やしました。なんか一度5のままやってみたのですが、励起エネルギーの最大エネルギーが3 eV 程度で、この色素の特徴的な波長には達してないというか、サンプルにあるような短波長側のサブピークが見えませんでした。サンプルからQM領域の原子を増やしたせいでしょうか？
スペクトル情報収集時のファイル名、awkコマンドがHandbookそのままだと上手く動かないと思われます。たぶんCP2Kのバージョンアップかなんかで、ファイル名とか出力フォーマットが変わったんじゃないかという気がします。

系の準備

まずパラメータ決めをします。

gmx pdb2gmx -f 4eul.pdb

これで分子の力場には 1 AMBER03 を、水の力場には 1 TIP3P を選択します。

続いて溶媒の水を追加します。

gmx solvate -cp conf.gro -o conf.gro -p topol.top -shell 10

で、次は中性化のためイオンを追加します。ただこれには tpr ファイルが必要なので、いったん適当な mdp インプットファイルでコンパイルします。一番エラーが出なさそうなエネルギー最小化のインプットを使います。

gmx grompp -f em.mdp -p topol.top -c conf.gro -o egfp-genions.tpr -maxwarn 1
gmx genion -s egfp-genions.tpr -p topol.top -o conf.gro -neutral

2番目のコマンドで、どの分子をイオンに置き換えるか聞かれますので、 13: SOL を選びます。

ここまで来たらいったん古典力場で最小化します。

gmx grompp -f em.mdp -p topol.top -c conf.gro -o egfp-em.tpr
gmx mdrun -deffnm egfp-em

これで古典力場によりエネルギー最小化された構造が egfp-em.gro ファイルに出力されます。

あと、QM原子団を設定した index.ndx を作っておきます。

gmx make_ndx -f conf.gro

make_ndx では次のように入力し、 CRO 残基を QMatoms と名前を付けます。 17 という数字はたぶん変わらないと思いますが、念のためグループ一覧を表示して今回作った r_CRO グループが該当する番号が何かを確認しておきましょう。

r CRO
name 17 QMatoms
q

これで指定された残基はこのような感じです。

NVTによる平衡化

まず古典力場により平衡化します。

gmx grompp -f md-mm-nvt.mdp -p topol.top -c egfp-em.gro -t egfp-em.trr -o egfp-mm-nvt.tpr
gmx mdrun -deffnm egfp-mm-nvt

なお、この古典力場MDをさぼって最小化構造から直接QMMMのNVT計算をするとQMMMの系が崩壊していました。

続いて、いよいよQMMM計算によるNVT MD計算です。といっても非常に重たいので、平衡化の手順と考えるとやらないよりマシ、レベルの時間長さです。

gmx_d grompp -f md-qmmm-nvt.mdp -p topol.top -c egfp-mm-nvt.gro -n index.ndx -o egfp-qmmm-nvt.tpr
gmx_d mdrun -deffnm egfp-qmmm-nvt -ntmpi 1

なんかここからはMPIスレッド数を１としまして、並列化をOpenMPで行うようにしました。これで2時間くらい待ちました。

このMDでの温度変化の様子です。

gmx energy -f egfp-qmmm-nvt.edr

このコマンドでどれを選ぶか聞かれますので、 14 Temperature を選びます。これで energy.xvg ファイルができます。これはgraceというソフトウェアでプロットできます。

xmgrace energy.xvg

NVT結合間隔が10ステップごとである様子が良く分かります。

スペクトル

つづいて、GROMACSのQM計算の応用方法の一つ、テンプレートを渡して自分なりの特殊な計算をさせるバージョンです。

先ほどのNVT計算中に、 egfp-qmmm-nvt_cp2k.inp というファイルができています。GROMACSの内部処理としてはこのファイルをCP2Kに渡して実行、結果をパースしているようです。
このファイルを基にテンプレートファイルを作ります。

cp egfp-qmmm-nvt_cp2k.inp egfp-qmmm-spec.inp
vi egfp-qmmm-spec.inp

$END DFT と &QMMM の間に、以下のように追記します。

  &END DFT

  &PROPERTIES
    &TDDFPT
       NSTATES     20
       MAX_ITER    10
       CONVERGENCE [eV] 1.0e-3
    &END TDDFPT
  &END PROPERTIES

  &QMMM

これは「Tutorial-Handbookからの変更点」でも書いたとおり、チュートリアルのものから変えています。 NSTATES はメインのピークの高周波側にあるサブピークがあるあたりの 4 eV がカバーできる程度の数としました。ほんとはもっと増やして 5 eV までカバーしたかったのですが、Gaussian のTDDFTにおけるState数よりも強く計算時間に効いてくるので、最低限の状態数としました。

さらにインプットファイルの修正を行います。まず md-qmmm-spec.mdp を修正します。

nsttcouple               = 1

これで準備ができました。

実行

コンパイルと実行です。

gmx_d grompp -f md-qmmm-spec.mdp -qmi egfp-qmmm-spec.inp -p topol.top -c egfp-qmmm-nvt.gro -n index.ndx -o egfp-qmmm-spec.tpr

テンプレートファイルの指定は -qmi egfp-qmmm-spec.inp で行っています。つまり mdp ファイル中ではないんですね。

ちなみに core-i7 14700F と RTX 5060 Ti でのパフォーマンスはこんな感じ：

               Core t (s)   Wall t (s)        (%)
       Time:  1416366.355    50584.513     2800.0
                         14h03:04
                 (ns/day)    (hour/ns)
Performance:        0.000   139121.323

NVTと同じ100ステップですが、今度は半日以上もかかっています。マシンパワーによりますが、ゲーミングPCとかなら寝る前に実行し、次の日の仕事が終わってから見ると良いと思います。

あんまり時間がかかるようなら、ある程度のステップ数進んでいたら中断してもよいと思います。次の結果を処理する計算は別にサンプル数がそろっていなくてもよさそうですし。

スペクトルの計算

スペクトル計算のTDDFPTはGROMACSのあずかり知らぬ処理です。なので、MDの構造の各ステップにおけるスペクトルを統計処理するところは自分でやる必要があります。

とりあえず、 QM計算方法に INPUT を指定していると、CP2Kの出力をすべてファイルに書き出しているようです。今回は egfp-qmmm-spec_cp2k.out にありますので、ここからスペクトルに相当するところ、エネルギー（eV）とその強度を抽出します。

grep " TDDFPT|" egfp-qmmm-spec_cp2k.out | awk '{if(NF == 7){print $3 " " $7 }}' > excitations

そして、これからスペクトルを描写するプログラムをチュートリアル作者が書いてくれています。

./conv.py excitations 0.1 2 5

これで、 spec.xvg というファイルができますので、 xmgrace spec.xvg としてみます。

横軸は eV です。縦軸はたぶん任意単位の強度かなと思います。
チュートリアルの結果にくらべ 4.2 eV くらいに余計なピークができてる気がする。たぶんチュートリアルにはない原子があるからかな？と思います。

そのほか

実はCUDA版をインストールしてまして、特に使わんだろうと思ってたんですが、

 SUBROUTINE                       CALLS  ASD         SELF TIME        TOTAL TIME
                                MAXIMUM       AVERAGE  MAXIMUM  AVERAGE  MAXIMUM
 pw_gpu_r3dc1d_3d                111965 15.0  3275.69  3275.69  3275.69  3275.69

という項目があったので、無駄ではなかったようです。・・・平面波？

あとスペクトル計算のところ、普通のQMMM MD計算して細かくトラジェクトリを出しておいて、そのトラジェクトリからTDDFT計算を実行するCP2K用インプットファイルを大量に作りクラスタで並行計算、統計処理した方が正直効率がいい気がする。なんなら、古典MDの軌跡からTDDFT計算していって統計処理しても問題ないことも多い気がする。