液体のベンゼンの系を作るときの苦労の記録

これは何

ふとベンゼンの液体の系を作ろうかと思い立ち、いろいろツールを使ってやってみたのですが、落とし穴というか手間取ったところとかあったのでメモしておこうかと思います。
今後何か簡単な分子の液体系を作るときに必要なコマンド・入力ファイルをコピペしてちょっと編集すればできるようにしておこうかと。

あらすじ

• openbabelで生成
• Ambertoolsでパラメータ決め
• Gromacsで系を作り、MD
  • エネルギー最小化
  • deform
  • NPT

分子生成とパラメータ決め

分子構造の生成

ベンゼン程度であればすべてコマンドで完結させてしまいたいところです。ベンゼンといえば c1ccccc1 です。この8文字を入力するだけで3次元の構造ファイルができたらいいなあ、ということで。

echo c1ccccc1 | obabel -ismi - -o pdb -Obenzene.pdb --gen3D

PDBファイルができました。

分子種決め

続いて、パラメータ決めします。

antechamber -fi pdb -i benzene.pdb -fo mol2 -o benzene.mol2 -at gaff2 -c bcc -pf y

これで GAFF2での原子タイプと電荷の情報を持った mol2 ファイルができます。

落とし穴

しかしここで落とし穴ポイントが一つ。原子タイプが正しく設定されていないことがあります。
benzene.mol2 の中身を見てみます。

benzene.mol2 の 8行目～

@<TRIPOS>ATOM
      1 C            1.3830    -0.2220     0.0050 c3         1 UNL      -0.090700
      2 C1           0.5070    -1.3070    -0.0080 c3         1 UNL      -0.090700

炭素原子が c3 、つまり sp3 炭素として認識されています。（ ca だったら問題ないですので、「結合の補正ファイル」の節までお進みください。）
このまま進むと結合情報の設定で失敗しますし、うまくいってもシクロヘキサンみたいに曲がってしまいます。
なおこの落とし穴は経験済みです。

https://x.com/MeguruChem/status/1759510328713560266

対策ですが、とりあえず benzene.pdb ファイル中で CONECT から始まる行を削除したもので antechamber を実行しますと

benzene.mol2 の 8行目～

@<TRIPOS>ATOM
      1 C            1.3830    -0.2220     0.0050 ca         1 UNL      -0.130000
      2 C1           0.5070    -1.3070    -0.0080 ca         1 UNL      -0.130000

と、ちゃんと芳香環炭素原子である ca となっています。これで次のステップに進みます。

ちなみに確認のため別環境でやってみると無編集のPDBファイルでも普通に芳香環炭素原子として認識されていました。 openbabel か ambertools のバージョンの問題・・・？よくわかりませんが、結果ファイルの確認はちゃんとしましょうということで。

結合の補正ファイル

あと、AmberToolsのGaff定義ファイルに若干足りない結合があるようなので、補正ファイルを作っておきます。

parmchk2 -i benzene.mol2 -f mol2 -o benzene.frcmod

パラメータファイルづくり

Amber

まず、 tleap でAmber用の top ファイルと crd ファイルを作ります。

cat << EOS | tleap -f -
source leaprc.gaff2
MOL = loadmol2 benzene.mol2
loadamberparams benzene.frcmod
saveamberparm MOL benzene.prmtop benzene.inpcrd
EOS

いや、 インプットファイルにして渡してもいいんですが、なんとなく端末のヒストリーからたどれるようにしたくて。

Gromacs

そしてできたAmber用パラメータファイルをGromacs用に変換します。

cat << EOS | python
import parmed as pmd
mol = pmd.load_file("benzene.prmtop", "benzene.inpcrd")
mol.save("benzene.top", format="gromacs")
mol.save("benzene.gro")
EOS

いや、pythonスクリプトファイルにしてもいいんですが、（略

溶液モデル準備

箱の用意

周期境界を設定し、その中に多数のベンゼン分子を入れます。

gmx editconf -f benzene.gro -box 10
gmx insert-molecules -f out.gro -ci benzene.gro -nmol 1000

周期境界を 10 nm立方とし、その中にベンゼン分子を1000個追加します。
全部で1001個のベンゼン分子が 10x10x10 nm の箱の中にいる状態です。

先の節で生成したgromacs用トポロジーファイルを編集し、分子数を修正します。

benzene.top の最終セクション

[ molecules ]
; Compound       #mols
UNL                  1001

ほぐす

あとはエネルギー最小化しておきます。構造最適化した分子を離して置いているからそんな影響ないだろと思ったのですが、意外とこのステップがないと安定しません。

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

なお、インプットファイルはこんな感じ：

em.mdp

integrator   = steep
nsteps       = 5000
nstlog       = 1
emstep = 0.0005
dt           = 0.0002
constraints              = h-bonds
constraint-algorithm     = Lincs
continuation             = no
Lincs-iter = 2
Lincs-order = 4
nstlist      = 1
ns_type      = grid
pbc          = xyz
rlist        = 1.0
coulombtype  = PME
rcoulomb     = 1.0
rvdw         = 1.0
cutoff-scheme = Verlet

これでできた em.gro をこれから圧縮していきます。

液体へ圧縮

さて、 10x10x10 nm立方に1001分子というのは結構薄くて、 0.130 g/ml 程度です。一方ベンゼンの液体は 0.876 g/ml あります。この液体とも気体とも言えない系を圧縮して液体にしてやります。

これくらいスペースをとってやらないと分子を追加するときに分子同士が衝突して上手くいかないことが多く、こういう風に広めにとって分子を並べ、圧縮して目的の凝縮系を得るというのはよく見る手順です。

deform

1001分子のベンゼン (分子量78.11) の系の密度を 0.876 g/ml としたいとき、一辺 5.292 nm の立方体にすればよいです。

$n$ を分子数、 $M_W$ を分子量、 $\rho$ が設定したい密度、 $N_\mathrm{A}$ をアボガドロ定数とすると

$$L = \sqrt[3]{\frac{nM_W}{\rho N_\mathrm{A}}}\times 単位変換$$

という式を考えると、

$$\sqrt[3]{\frac {78.11 \times 1001}{6.02214076\times 10^{23} }\times \frac{1}{0.876} } \times 10 \times 1000 \times 1000 = 5.29210641994 \mathrm {nm}$$

となります。三乗根の中身はベンゼン1001分子が与えられた密度でどれだけの体積を占めるかを計算し、末尾の $\times 10 \times 1000 \times 1000$ で cm → mm → μm → nm の変換をしています。

1 ns のあいだに 10 nm から 5.292 nm に縮めるとき、 4.708 nm 縮めますので、 速度は 0.004708 nm/ps です。これを deform に設定します。

compress.mdp

nsteps       = 500000
nstxout      = 2500
nstlog       = 1000
dt           = 0.002
constraint-algorithm   = lincs
constraints = h-bonds
continuation = no
cutoff-scheme = Verlet
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       = no
gen-vel      = yes
gen-temp     = 300

deform  = -4.708e-3 -4.708e-3 -4.708e-3 0 0 0
deform-init-flow = yes

これでMDを実行します。

gmx grompp -f compress.mdp -c em.gro -p benzene.top -o compress.tpr
gmx mdrun -deffnm compress

マシンにもよりますが、10分程度で終わるかなと思います。圧縮の様子を見てみましょう：

上手くできているようです。

常温常圧

300K 1atmのNPTで安定するか試してみます。
圧縮率はこちらを参考に設定しました。
https://faq.an.shimadzu.co.jp/faq/show/2971?site_domain=default

npt.mdp

integrator   = md
nsteps       = 1000000
nstxout      = 5000
dt           = 0.002
constraint-algorithm   = lincs
constraints = h-bonds
continuation = no

cutoff-scheme = Verlet
nstlog       = 1000
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        = 10.0
ref-p        = 1.0
pcoupltype   = isotropic
compressibility = 1.0e-4
gen-vel      = no

これで実行します。

gmx grompp -f npt.mdp -c compress.gro -p benzene.top -o npt.tpr
gmx mdrun -deffnm npt
gmx energy -f npt.edr # Density を選択
xmgrace energy.xvg

よさそうですね。

系の見た目はこんな感じです。分子がぎっしり。

NPTで縮められる？

素直な発想として、常温常圧で液体なら deform しなくても、薄い系から圧力熱浴をつけてNPT MDを実行しても凝縮して液体の密度付近に落ち着くのでは？とお思いかもしれません。しかしやってみると：

むしろ気体側に行ってしまいました。なんか蒸気圧が過飽和でも核がないので液体にならない、みたいな状態ですね。まあMDのタイムスケールだと核が急にできても熱浴でうまく冷やしてやらないと破綻しそうな気がします。

じゃあ100気圧とか高い圧力をかけるとどうかという感じですが、以下の点からお勧めしません：

• 壁の動きが大きく破綻しやすい。
• 破綻を避けるため、圧縮率を小さくするとか、tau-p を長くするとかの恣意的な工夫が必要。
• 圧力をどれくらいにすれば蒸気圧に打ち勝って圧縮できるかは系によるので、試行錯誤のステップが要る。圧力を大きくしすぎるとこれも破綻しやすい。
• 圧縮できたところで ref-p 周りで振動するので、目標とする密度に設定するときに苦労する。

正直申しますと、昔は溶液系を作るのに初期配置を目標の密度にできない場合、高い圧力で無理やり縮める方法を使ってたのですが、最近 deform という手段に気づきまして、それでこのメモ記事が生えました。