🐸

石の表面のペプチドをGROMACSで

2025/05/09に公開

あらすじ

前回までの記事で、INTERFACE-MDという~~おもちゃ~~ 無機結晶用力場を見つけたので、せっかくなのでこれを使って何かチュートリアルにないことをしたいと思います。

ということで、pyrophyllite という石の結晶力場がありましたので、これをGROMACSで利用、ついでに小さいペプチドで有名な chignolin を乗せてみようと思います。

環境

Ubuntu 24.04.2 LTS
GROMACS 2025.1 GPU版

準備

今回はAmberToolsを使いません。
代わりに、INTERFACE-MDでも紹介されているGROMACSに変換するツールを使います。 $HOGE を適当なディレクトリとします。

cd $HOGE
git clone https://github.com/kolmank/interfaceff2gro
cd interfaceff2gro/msi2lmp/src
make
cd $HOGE

またこれは MDAnalysis というライブラリを使います。今回ｳﾁは pip を使いましたが、人によっては conda とか mamba を使った方がいいかもしれません。

pip install --upgrade MDAnalysis

前回の記事でも使ったパラメータを準備します。

unzip interface_ff_1_5.zip

これで、 HOGE/INTERFACE_FF_1_5 ができます。

トポロジー作成

作業ディレクトリを作成、以降このディレクトリ中にファイルを準備、各種コマンドを実行していきます。

mkdir pyrophyllite_chignolin
cd pyrophyllite_chignolin

石スラブ

ファイル変換。interfaceff2gro を使って INTERFACE-MD で提供されているモデルをGROMACS用に変換します。

$HOGE/interfaceff2gro/interfaceff2gro.py $HOGE/INTERFACE_FF_1_5/MODEL_DATABASE/CLAY_MINERALS/pyrophyllite15_single_layer.car
cp ../INTERFACE_FF_1_5/MODEL_DATABASE/CLAY_MINERALS/pyrophyllite15_single_layer.* .

あと system.top が作業ディレクトリにできていると思いますが使いません。

editconf を使って周期境界条件を調整します。

gmx editconf -f pyrophyllite15_single_layer.gro -o slab.gro -box 2.58000 2.68980 3.6

slab.gro はこんな感じ。

ペプチド

この論文

のチュートリアル例２がChignolinのフォールディングエネルギー曲面の計算です。途中までこの例に従います。

wget https://files.rcsb.org/download/1UAO.pdb
head -n 281 1UAO.pdb > 1UAO_model1.pdb

1UAO_model1.pdb はこんな感じ。

力場の種類を選択します。INTERFACE-MD に従いCHARMM-27、TIP4Pとします。

gmx pdb2gmx -f 1UAO_model1.pdb -ignh -o 1uao.gro

端末

Select the Force Field:

From '/usr/local/gromacs/share/gromacs/top':
 1: AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24, 1999-2012, 2003)
 2: AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
 3: AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29, 461-469, 1996)
 4: AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21, 1049-1074, 2000)
 5: AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65, 712-725, 2006)
 6: AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al., Proteins 78, 1950-58, 2010)
 7: AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
 8: CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
 9: GROMOS96 43a1 force field
10: GROMOS96 43a2 force field (improved alkane dihedrals)
11: GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
12: GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
13: GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
14: GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856, DOI: 10.1007/s00249-011-0700-9)
15: OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
8 ← 8を入力してEnter

Using the Charmm27 force field in directory charmm27.ff
Opening force field file /usr/local/gromacs/share/gromacs/top/charmm27.ff/watermodels.dat

Select the Water Model:
 1: TIP3P   TIP 3-point, recommended
 2: TIP4P   TIP 4-point
 3: TIPS3P  CHARMM TIP 3-point with LJ on H's
 4: TIP5P   TIP 5-point (see https://gitlab.com/gromacs/gromacs/-/issues/1348 for issues)
 5: SPC     simple point charge
 6: SPC/E   extended simple point charge
 7: None
2 ← 2を入力してEnter

続いて、石スラブにchignolinを追加します。 slab.gro は周期境界長さを定義しているので、こちらを基とします。

gmx insert-molecules -f slab.gro -ci 1uao.gro -nmol 1 -o slab_chignolin.gro

これでこんな感じの系になります。

topol.top の最初に石スラブのitpファイルをincludeするように追記するのと、最後の方にある molecules ブロックに石スラブを追加します。

topol.top 21行目～

#include "charmm27.ff/forcefield.itp"
#include "pyrophyllite15_single_layer.itp"

[ moleculetype ]
; Name            nrexcl
Protein_chain_A     3

topol.top 最後の方

[ molecules ]
; Compound        #mols
SURF 1
Protein_chain_A     1

水で満たす

水溶液にします。

gmx solvate -cp slab_chignolin.gro -cs tip4p.gro -p topol.top -o solv.gro

solv.gro はこんな感じになります。

拘束条件

石スラブの拘束条件を追加できるようにします。今思うとイオンを加えた後でもよかった。

gmx make_ndx -f solv.gro #  q <Enter> として、そのまま index.ndx を作成
gmx genrestr -f solv.gro -o posre_SURF.itp

端末

Select group to position restrain
Group     0 (         System) has  2934 elements
Group     1 (          Other) has   600 elements
Group     2 (           SURF) has   600 elements
Group     3 (        Protein) has   138 elements
Group     4 (      Protein-H) has    77 elements
Group     5 (        C-alpha) has    10 elements
Group     6 (       Backbone) has    30 elements
Group     7 (      MainChain) has    39 elements
Group     8 (   MainChain+Cb) has    46 elements
Group     9 (    MainChain+H) has    50 elements
Group    10 (      SideChain) has    88 elements
Group    11 (    SideChain-H) has    38 elements
Group    12 (    Prot-Masses) has   138 elements
Group    13 (    non-Protein) has  2796 elements
Group    14 (          Water) has  2196 elements
Group    15 (            SOL) has  2196 elements
Group    16 (      non-Water) has   738 elements
Select a group: 2
Selected 2: 'SURF'

pyrophyllite15_single_layer.itp を編集します。位置拘束条件を有効にするフラグ変数名がchignolinのものとかぶるので、こちらを POSRES_SURF と変更します。

pyrophyllite15_single_layer.itp 15056行付近

; Include Position restraint file
#ifdef POSRES_SURF
#include "posre_SURF.itp"
#endif

イオン追加で中性化

Chignolin はいくつかカルボン酸を持っていて、このままだと系が電荷をもっています。この状況は特に周期境界条件を考えると物理的にもおかしいことになりますので、カウンターイオンを加えて中性化します。

コンパイル済みトポロジーファイルが欲しいので、罪のないインプットファイル ions.mdp を準備します。

ions.mdp

define      = -DPOSRES_SURF

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
cutoff-scheme   = Verlet    ; Buffered neighbor searching
coulombtype     = cutoff    ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

coulombtype がEwald和を使うようなやつだと系が中性でないと文句を言うので cutoff にします。

gmx grompp -f ions.mdp -c solv.gro -p topol.top  -r solv.gro -o ions.tpr 
gmx genion -s ions.tpr -o solv_ions.gro -p topol.top -pname NA -nname CL -neutral

genionの出力

Will try to add 2 NA ions and 0 CL ions.
Select a continuous group of solvent molecules
Group     0 (         System) has  2934 elements
Group     1 (          Other) has   600 elements
Group     2 (           SURF) has   600 elements
Group     3 (        Protein) has   138 elements
Group     4 (      Protein-H) has    77 elements
Group     5 (        C-alpha) has    10 elements
Group     6 (       Backbone) has    30 elements
Group     7 (      MainChain) has    39 elements
Group     8 (   MainChain+Cb) has    46 elements
Group     9 (    MainChain+H) has    50 elements
Group    10 (      SideChain) has    88 elements
Group    11 (    SideChain-H) has    38 elements
Group    12 (    Prot-Masses) has   138 elements
Group    13 (    non-Protein) has  2796 elements
Group    14 (          Water) has  2196 elements
Group    15 (            SOL) has  2196 elements
Group    16 (      non-Water) has   738 elements
Select a group: 15 
Selected 15: 'SOL'

solv_inos.gro はこんな感じに出来上がりました。

周期境界条件のイメージセルも描いてみるとこんな感じになります。

これで準備完了です。

MD実行

エネルギー最小化

em.mdp

define      = -DPOSRES_SURF

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator  = steep         ; Algorithm (steep = steepest descent minimization)
emtol       = 1000.0        ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep      = 0.01          ; Minimization step size
nsteps      = 50000         ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
cutoff-scheme   = Verlet    ; Buffered neighbor searching
coulombtype     = PME       ; Treatment of long range electrostatic interactions
rcoulomb        = 1.0       ; Short-range electrostatic cut-off
rvdw            = 1.0       ; Short-range Van der Waals cut-off
pbc             = xyz       ; Periodic Boundary Conditions in all 3 dimensions

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

緩和された構造 em.gro が得られます。

見ただけではちょっとよくわかりませんね。

スラブに付着

人工力でchignolinを石スラブに押し付けてみます。物理的におかしい構造になる可能性を否定しきれませんが、まあこのあと平衡化MDをしっかり流せば大丈夫でしょう。知らんけど。

pull.mdp

define    = -DPOSRES_SURF

; RUN CONTROL PARAMETERS
integrator   = md
dt           = 0.002
nsteps       = 100000
;
nstxout      = 500
nstlog       = 500
nstenergy    = 500
;
;nstlist      = 1
pbc          = xyz
rlist        = 1.0
;
coulombtype  = PME
rcoulomb     = 1.0
;
vdwtype      = cut-off
rvdw         = 0.8
dispcorr     = enerpres
;
tcoupl       = v-rescale
tc_grps      = System
tau_t        = 0.1
ref_t        = 300.0 
;
pcoupl       = no
;
gen_vel      = yes
gen_temp    = 300
gen_seed    = -1
;
constraints  = h-bonds
constraint_algorithm = LINCS

pull         = yes
pull-ngroups = 2
pull-ncoords = 1

pull-group1-name = SURF
pull-group1-pbcatom = -1
pull-group2-name = Protein

pull-coord1-type = umbrella
pull-coord1-geometry = distance
pull-coord1-groups = 1 2
pull-coord1-dim = N N Y
pull-coord1-k = 500

このファイルをgromppでコンパイル、 mdrun で実行します。 pull-group1-name で指定するのがGROMACSが自動で定義される原子グループなので、別に index.ndx を渡さなくてもよいようです。

gmx grompp -f pull.mdp -c em.gro -p topol.top  -r solv_ions.gro -o pull.tpr
gmx mdrun -deffnm pull

軌跡をVMDなどで見やすいよう、chignolinを中心に持ってきつつ(-center オプション、このあと対話処理で中心に持ってくるグループを選択できる)、周期境界条件を超えた分子を戻してやります(-pbc molオプション)。 xtc ファイルなのは、軌跡を動画で見るだけなら十分な精度でファイルが非常に軽くなるので xtc ファイルにしています。

gmx trjconv -f pull.trr -s pull.tpr -o pull.xtc -pbc mol -center

ここで変換中の出力に

 There were 126 inconsistent shifts. Check your topology

という警告？がいっぱい出るかと思いますが、これは石スラブが周期境界条件を超えた結合をもつためのようです。

参考： https://gromacs.bioexcel.eu/t/inconsistent-shifts-during-trjconv/3238

動画はこんな感じ。

思ったより強く押し付けているようです。

平衡化

NPT平衡化を行います。以下のファイルを準備します。

npt.mdp

define    = -DPOSRES_SURF

; RUN CONTROL PARAMETERS
integrator   = md
dt           = 0.002
nsteps       = 100000
;
nstxout      = 500
nstlog       = 500
nstenergy    = 500
;
;nstlist      = 1
pbc          = xyz
rlist        = 1.0
;
coulombtype  = PME
rcoulomb     = 1.0
;
vdwtype      = cut-off
rvdw         = 0.8
dispcorr     = enerpres
;
tcoupl       = v-rescale
tc_grps      = System
tau_t        = 0.1
ref_t        = 300.0 
;
pcoupl       = c-rescale
pcoupltype   = semiisotropic
tau_p        = 5.0
compressibility = 0 4.5e-5
ref_p        = 1.0 1.0
refcoord_scaling = all
;
gen_vel      = yes
gen_temp    = 300
gen_seed    = -1
;
constraints  = h-bonds
constraint_algorithm = LINCS

ポイント：

石スラブは動かないように（ define = -DPOSRES_SURF ）
圧力の変化はｚ軸のみの変動で行います。 (pcoupltype = semiisotropic と compressibility = 0 4.5e-5)

これをコンパイル、実行します。

gmx grompp -f npt.mdp -p topol.top -c pull.gro -r solv_ions.gro -o npt.tpr
gmx mdrun -deffnm npt

これもまたＶＭＤなどで見やすいように調整します。

gmx trjconv -f npt.trr -s npt.tpr -pbc mol -center -o npt.xtc

これをVMDで上手く描画してやるとこうなります。

逆側から映していてわかりにくくて申し訳ないのですが、両端のチロシン・トリプトファンの芳香環をずっとスラブにつけてもぞもぞ動いている感じです。

スラブからはみ出しているように見えますが、周期境界条件を考えるとこうなりますので大丈夫です。

サンプリング

これは演習問題とする。

計算速度

ちなみに全部で2928座標（2381原子＋547仮想サイト）で、GTX 1050 Ti、 AMD Ryzen 5 3600 でNPTを実行したときのパフォーマンスはこんな感じでした。

                 (ns/day)    (hour/ns)
Performance:      344.662        0.070

その他

より大きい系をシミュレーションするには、この結晶構造をコピーして並べる必要がありますが、残念ながら周期境界条件を超えて連続した結合を考慮して並べる方法が分かりませんでした。なので今回はここまで。

結晶上で付着した状態ではフォールディングの自由エネルギー曲面が水溶液中とどう違うかとか面白いんじゃないかなと思います。

追記

mdpファイルで設定できるオプションに

periodic-molecules = yes

があったのですが、設定をしてもぱっと見違いは見当たりませんでした。

あらすじ

環境

準備

トポロジー作成

石スラブ

ペプチド

水で満たす

拘束条件

イオン追加で中性化

MD実行

エネルギー最小化

スラブに付着

平衡化

サンプリング

計算速度

その他

追記

Discussion