GROMACSで電場を揺らして分子を揺らす

これは何

皆さんはたまによく使う計算ソフトの入力パラメータのドキュメントを眺めたりすると思います。
この間Gromacsのドキュメント見てたら、パルスや電場をかける設定を見かけました。
ちょっと使ってみようと思います。

液体の水の系を用意し、水分子の特定の分子振動だけを揺動してみようと思います。

水にしたのは、振動数が違う分子振動が２つあること、Gromacsに溶媒として準備されているので分子データを準備しなくてもいいことがあり、パルス以外のことで悩むことが少なく簡単にできるかなと思ってのことです。

ファイル準備

トポロジーファイル

力場とかの設定です。 gromacs/share/gromacs/top/amber03.ff/ の中を漁って、TIP3P FLEXIBLEの内容をとってきました。amber03.ff/forcefield.itp と amber03/tip3p.itp を include し、mdpファイル中で define = -DFLEXIBLE とすれば [ system ] 以下だけで終わりそうですが、具体的な top ファイルの構造の簡単なテンプレートが欲しかったので。
[ molecules ] の節は普通は 分子種 分子数 の行が続くのですが、今回はわざと開けています。 gmx solvate で追加された分子数分追加されるので、そこで追記してもらいます。

system.top

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes             1            0.83333333

[ atomtypes ]
; name    at.num    mass    charge ptype  sigma      epsilon
OW           8      16.00    0.0000  A   3.15061e-01  6.36386e-01
HW           1       1.008   0.0000  A   0.00000e+00  0.00000e+00

[ moleculetype ]
; molname       nrexcl
SOL             2

[ atoms ]
; id  at type     res nr  res name  at name  cg nr  charge    mass
  1   OW          1       SOL       OW       1      -0.834    16.00000
  2   HW          1       SOL       HW1      1       0.417     1.00800
  3   HW          1       SOL       HW2      1       0.417     1.00800

[ bonds ]
; i     j       funct   length  force_constant
1       2       1       0.09572 502416.0   0.09572        502416.0
1       3       1       0.09572 502416.0   0.09572        502416.0

[ angles ]
; i     j       k       funct   angle   force_constant
2       1       3       1       104.52  628.02      104.52  628.02

[ system ]
; Name
Generic title in water

[ molecules ]
; Compound       #mols

あとはエネルギー最小化と平衡化用NPTのインプット。OH伸縮を固定してませんので、0.1 fs刻みです。詳細は省きます。

em.mdp

integrator   = steep
nsteps       = 5000
nstlog       = 1
emstep = 0.0005
dt           = 0.0002
constraints              = none
nstlist      = 1
ns_type      = grid
pbc          = xyz
rlist        = 1.0
coulombtype  = PME
rcoulomb     = 1.0
rvdw         = 1.0
cutoff-scheme = Verlet

npt.mdp

integrator   = md
nsteps       = 1000000
nstxout      = 10000
dt           = 0.0001
constraints              = none
cutoff-scheme = Verlet
nstlog       = 10000
nstlist      = 10
ns_type      = grid
pbc          = xyz
rlist        = 1.0
coulombtype  = PME
rcoulomb     = 1.0
rvdw         = 1.0
tcoupl       = V-rescale
tc-grps      = System
tau-t        = 1.0
ref-t        = 300
pcoupl       = c-rescale
tau-p        = 1.0
ref-p        = 1.0
pcoupltype   = isotropic
compressibility = 4.5e-5
gen-vel      = yes
gen-temp     = 300

系の準備

gmx solvate を入力gro無しでやるとboxいっぱいに水を入れてくれます。で、ここでtopファイルを渡していたら最後にSOLの分子数設定をしてくれるので、わざわざ出力見て自分で修正する必要はないです。ただしtopファイルに元からSOLの分子数設定行が書いてあったとしてもそれを上書きせず、追記します。

gmx solvate -box 5 5 5 -p system.top

そのままエネルギー最小化と平衡化のMDを実行します。

gmx grompp -f em.mdp -p system.top -c out.gro -o em.tpr
gmx mdrun -deffnm em
gmx grompp -f npt.mdp -c em.gro -p system.top -o npt.tpr
gmx mdrun -deffnm npt

最後の npt.gro をVMDで見てみます。

何の変哲もないTIP3P FLEXIBLEの水です。これに電場をかけてエネルギーを与えてやります。

実験

電場を揺らす設定のインプットファイル excite.mdp です。20psのNVE、つまり熱浴は付けていません。

excite.mdp

integrator   = md
nsteps       = 200000
nstxout      = 1000
dt           = 0.0001
constraints   = none
cutoff-scheme = Verlet
nstlog       = 1000
nstenergy    = 10
nstlist      = 10
pbc          = xyz
rlist        = 1.0
coulombtype  = PME
rcoulomb     = 1.0
rvdw         = 1.0

tcoupl       = no
pcoupl       = no


electric-field-x = 0.2 390 2 0.5
electric-field-y = 0.2 725 12 0.5

GROMACSの電磁場は、ドキュメントによると

$$E_0\exp\left(-\frac{\left(t-t_0\right)^2}{2\sigma^2}\right)cos\left[\omega\left(t-t_0\right)\right]$$

で表され、 electric-field-x などで与える4つの数字は $E_0\ \omega\ t_0\ \sigma$ に相当します。今回は $E_0$ はMDが破綻しないがエネルギーのグラフで分かりやすいくらい与えられる強さで、 $t_0$ は二つのパルスが十分離れているくらいにして、 $\sigma$ はあまり長すぎず短すぎないようにするくらいは気をつけましたが、だいたいは適当に決めました。そして $\omega$ についてはOH伸縮またはHOH変角の周波数に合わせています。
このインプットでは、 x軸方向に2psの時点でHOH変角に相当する周波数のパルスを、y軸に沿って12psの時点でOH伸縮に相当する周波数のパルスを約 1ps の間適用してます。
でもｳﾁのざっくりとした計算だとOH伸縮は500 ps^{-1} 弱くらいだと思ってたんですが、何かと結合してしまってるんでしょうか？もしくは単位変換ミスった？

これで実行します。core-i7 14700FとRTX5060Tiで1分半くらいでした。

gmx grompp -f excite.mdp -c npt.gro -p system.top -o excite.tpr
gmx mdrun -deffnm excite -field field.xvg

これで実際に水分子が揺らされている様子を確認したいのですが、残念ながら軌跡をそこまで細かくとっていないので、電磁パルスを受けて水分子が激しく動き出す様子は見れません。ストレージをケチったせいです。かわりにエネルギーを細かく出していますので、こちらで確認します。

gmx energy -f excite.edr # 1 2 8

これでまずは結合エネルギー、角度エネルギー、全エネルギーの時系列を得ます。
結合エネルギーはOH伸縮、角度エネルギーはHOH変角に相当します。

これらを並べてプロットしてみます。

このグラフを描くのに実行したgnuplotコマンドです。gnuplotは@と#をコメント文字として設定してやればxvgファイルをプロットできます。

set bmargin 0
set tmargin 0
set lmargin at screen 0.25
set multiplot layout 3,1
set format x ""
p "energy.xvg" u 1:3 w l t "Angle"
set ylabel "Energy/kJ/mol"
p "energy.xvg" u 1:2 w l t "Bond"
unset ylabel
set format x "%g"
set xlabel "time/ps"
set bmargin at screen 0.1
p "energy.xvg" u 1:4 w l t "Total Energy" lw 3

こんな感じで、全エネルギーが2回のパルスそれぞれで増えていますが、結合エネルギーは1回目のパルスでは増えておらず2回目のパルスで増大、一方で変角はその逆である様子がわかります。分子振動の周波数に応じたパルスを当てて励起できています。

よく見ると、変角のパルスで増えたエネルギーが徐々に減っている様子も分かります。

さらに電場の情報がfield.xvgにありますので、これも併せてプロットします。

set datafile commentschars "@#"
p "energy.xvg" u 1:($2*0.0001 - 1) w l t "bond"
rep "energy.xvg" u 1:($3*0.0002 - 1) w l t "angle"
rep "field.xvg" u 1:2 w l t "pulse_x"
rep "field.xvg" u 1:3 w l t "pulse_y"

今回は同じグラフに重ねて出すため、スケールしてシフトしています。

パルス幅とエネルギー上昇区間が似たような感じですね。

余談

これはジャンルとしては非平衡MDです。パルスによりエネルギーを与えるシミュレーションで計算するものとしたら応答関数か熱流かな？あまり詳しくはないですが。

パルスを当てて何がわかるかですが、今回は分子振動の周波数とパルスの周波数が一致もしくは都合がいい関係になると励起する様子を見ていました。ですが分子振動の周波数を知るだけなら普通の平衡MDでなんかの自己相関関数からスペクトルを調べればいいです。それに限らず現象が小さい非平衡でよければ、非平衡の計算をしなくても揺動散逸定理から平衡状態のゆらぎの相関関数で結構いろんなことが説明できます。
ではわざわざ非平衡MDをする動機ですが、思いつくのとしては、例えば変角のエネルギーが徐々に減っていますが、これの時定数だとか、どこの分子振動が代わりに増えるかとか、そういった応答とか散逸過程を具体的に調べられるかなと思います。
また、パルスを当てた直後は系が平衡からずれていて、そのときは平衡状態に比べてパルスを受けたときの応答が変わります。こういう多次元分光のシミュレーションも非平衡状態を与えてやればできるかなと思います。

そういえば、多次元分光のシミュレーションについて、シグナルに熱浴に起因する構造が見えてしまったという話を聞いたことがある気がします。確かにエネルギー散逸を実際に見るようなMDにおいては、熱浴もエネルギーが伝わっている大きな経路の一つなので、さもありなんと思います。今回NVEにしたのもそういう背景があるからです。
でもタンパク質など、熱浴をつけておきたい系もあると思うのですが、そういうときはどうすればいいのでしょうか。。。良く知らないのですが、思いつきとしてはタンパク質には熱浴をつけず、溶媒の水につけるとかでしょうか？あるいは気にせずつけても意外と上手くいくのかもしれません。