事例9.1：GROMACSを用いたTIP3P水溶液中の高分子（PEO）MD

事例9.1：GROMACSを用いたTIP3P水溶液中の高分子（PEO）MD

この事例では、fbtk を用いて高分子（PEO）1鎖を大量のTIP3P水分子で溶媒和させたシステムを構築し、最新の OpenFF (Sage 2.2.1) 力場と TIP3P水モデル を適用して GROMACS 形式（.gro, .top）で出力・実行する手順を解説します。

水モデルの取り扱いとシミュレーションの安定性

FBTK + OpenFF + GROMACS の組み合わせは、水モデル（本事例では TIP3P）を極めて合理的かつ堅牢に取り扱います。

1. FBTK による配置: 大量の溶媒分子を高分子の周囲にパッキングする際、最も重要なのは「原子同士の重なり（Overlap）を避ける」ことです。FBTK は UFF ベースの暫定的な形状を用いて、重なりのない高品質な初期配置を高速に作成します。
2. OpenFF による設計図の割り当て: to_openff() を呼び出した際、OpenFF Toolkitは、指定された ["openff-2.2.1.offxml", "tip3p.offxml"] の力場ファイルを統合し、PEO分子にはSage 2.2.1のパラメータを、水分子にはTIP3P水モデルの厳密なパラメータ（O-H距離 0.9572 Å、H-O-H角度 104.52°、電荷など）を設計図として割り当てます。
3. GROMACS による幾何構造の強制と維持: GROMACSの.topファイルには水分子（SOL）に対して[ settles ]ディレクティブが記述されます。これにより、GROMACSでエネルギー最小化 (EM) を実行すると、FBTKが作った暫定構造はTIP3Pの理想形状へと即座に収束します。さらに、MD実行中も水分子はsettlesによって厳密な剛体として固定され、物理的に正しい挙動を示します。PEO分子などの他のH結合に対しては、MD時の constraints = h-bonds 指定により拘束が適用されます。

構築スクリプト (run.py)

import fbtk
from openff.toolkit import ForceField # ForceFieldクラスをインポート

# 1. PEO (20mer) の定義
monomer = fbtk.Molecule.from_smiles("*COC*")
peo = fbtk.Molecule.from_polymer(monomer, degree=20, name="PEO")

# 2. 水分子 (SOL) の定義
water = fbtk.Molecule.from_smiles("O", name="SOL")

# 3. 混合系の構築と緩和
builder = fbtk.Builder(density=1.0)
builder.add_molecule(peo, count=1)
builder.add_molecule(water, count=1000)
system = builder.build()

# 構造緩和（OpenFFへ渡す前の歪み解消）
system.relax()

# 4. OpenFF 形式へ変換し、GROMACS 形式で保存
# OpenFF 2.2.1 と TIP3P 水モデルを両方適用
combined_forcefield = ForceField("openff-2.2.1.offxml", "tip3p.offxml")
interchange = system.to_openff(forcefield=combined_forcefield)
interchange.to_gromacs("peo_water_tip3p") # 出力ファイル名も変更

GROMACS 設定ファイル (.mdp)

計算を再現するために必要な設定ファイルの内容です。

エネルギー最小化 (em.mdp)

integrator  = steep
emtol       = 1000.0
emstep      = 0.01
nsteps      = 5000

cutoff-scheme   = Verlet
nstlist         = 20
vdw-type        = Cut-off
vdw-modifier    = Force-switch
rvdw-switch     = 1.0
rvdw            = 1.2
coulombtype     = PME
rcoulomb        = 1.2

; 最小化中は拘束をかけず、FBTK構造を力場の理想形へ自然に収束させます
constraints     = none
pbc             = xyz

NVT平衡化 MD (md.mdp)

integrator      = md
dt              = 0.002
nsteps          = 10000

nstxout-compressed = 1000
nstlog             = 1000
nstenergy          = 1000

cutoff-scheme   = Verlet
nstlist         = 20
vdw-type        = Cut-off
vdw-modifier    = Force-switch
rvdw-switch     = 1.0
rvdw            = 1.2
coulombtype     = PME
rcoulomb        = 1.2

tcoupl          = v-rescale
tc-grps         = System
tau-t           = 0.1
ref-t           = 300

; Constraints (PEO H-bonds; Water rigidity is handled by 'settles' in .top file)
constraints     = h-bonds
constraint-algorithm = LINCS
pbc             = xyz

実行コマンド

出力された peo_water_tip3p.gro, peo_water_tip3p.top と上記 MDP ファイルを用いて実行します。

# 1. エネルギー最小化
gmx grompp -f em.mdp -c peo_water_tip3p.gro -p peo_water_tip3p.top -o em.tpr
gmx mdrun -v -deffnm em

# 2. MD実行
gmx grompp -f md.mdp -c em.gro -p peo_water_tip3p.top -o md.tpr -maxwarn 1
gmx mdrun -v -deffnm md

# 3. マルチフレーム gro への変換（解析準備）
echo 0 | gmx trjconv -f md.xtc -s md.tpr -o trajectory.gro

解析への応用 (fbtk-analyze)

fbtk-analyze を使うと、インデックスファイルなしで直感的な解析が可能です。

例：PEO の酸素と、水の酸素の間の RDF を計算

fbtk-analyze rdf -q "resname PEO and element O with resname SOL and element O" trajectory.gro -o peo_sol_rdf.dat

実行結果と検証

最新の OpenFF Sage 2.2.1 を用いた GROMACS 計算の結果、以下の安定性が確認されました。

• 収束の速さ: わずか 76 ステップでエネルギー最小化が完了しました。拘束をかけずに EM を行うことで、FBTK のパッキング構造が力場の理想幾何構造へ速やかに Snapping（変形・収束）したことを示しています。
• 熱力学データ: 300K で安定したシミュレーションが継続し、PEO と SOL という残基名が正しく管理されていることが確認されました（Performance: ~78.5 ns/day）。

まとめ

FBTK + OpenFF + GROMACS のツールチェーンにより、最新力場を用いた溶液系の MD 準備から解析までの全工程が大幅に効率化されます。