クラスターモデルの組み方

pdb2reaction は、基質の周りの残基をタンパク質から切り出し、切った結合を水素（キャップ水素）で埋めたクラスターモデルで計算します。

早見表

目的

すること

節

デフォルトのモデルを組む

-c に基質・補因子・金属を並べる

デフォルトのモデルを組む

重い計算の前に原子数と電荷を見る

extract を単独で実行するか、all --dry-run

モデルを確かめる

計算を軽くする

--exclude-backbone、--no-include-h2o、-r 0 --selected-resn、PyMOL で手で削る

モデルを削る

反応に関わる残基が入っていない、または切れ端だけ入っている

-r を広げる、-c に足す

モデルを広げる

自分で組んだモデルを使う

-c を省き、境界の原子を --freeze-atoms で固定する

自分で組んだモデルを使う

境界を手で決める・確かめる

チェックリストに沿って確かめる

手でモデルを組む・確かめるとき

切り出さずに酵素全体で経路を見る

path-opt・path-search（MEP の探索まで）

切り出さずに酵素全体で計算する

原子を固定する・距離を拘束する

--freeze-links（デフォルトで有効）、--freeze-atoms、--distance-restraint

原子の固定と距離の拘束

デフォルトのモデルを組む

-c に基質を書き、補因子と金属も同じリストに入れてください。リストの各残基から、それぞれ距離での探索が始まります。

pdb2reaction extract -i 1.R.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' -o model.pdb --out-json

これは examples/ の同梱例です。モデルは model.pdb に書き出され、同じディレクトリの result.json に、キャップ水素を含まない原子数（n_atoms_extracted）、キャップ水素の数（n_link_hydrogens）、電荷（total_charge）が入ります。モデルの原子数は n_atoms_extracted + n_link_hydrogens です。

デフォルトのモデルは次のように組まれます。

  • 半径：-c の原子から -r（デフォルト 2.6 Å）以内に原子が 1 つでもある残基を入れます。その原子が主鎖の原子なら、両隣の残基も入ります。水も入ります（--include-h2o）。

  • 切り方：つながったアミノ酸の並びは主鎖を残し、両端を CA で切ります。前後の残基がモデルに入らなかった残基は CB で切り、側鎖だけを残します。-c に入れたアミノ酸は、ペプチド結合でつながった隣の残基が半径に入るので、全部の原子を残します。

  • キャップ水素：CA か CB で切った結合ごとに、その炭素から 1.09 Å の位置に水素を置きます（残基 LKH の原子 HL）。

  • 複数の構造：反応物と生成物を 1 回で渡すと、構造ごとに選ばれた残基の和集合をすべての構造に使うので、どのモデルも同じ境界になります。入力では、同じ原子が同じ順に並んでいる必要があります。同梱例では、1.R.pdb だけなら 632 原子・電荷 +1、1.R.pdb と 3.P.pdb を一緒に渡すと 669 原子・電荷 0 になります。3.P.pdb では Cys167 の主鎖の水素が半径に入り、Cys167・Val166・Asp168（−1）がその実行のすべてのモデルに加わるためです。

モデルを確かめる

重い計算の前に、原子数と電荷を確かめてください。上のように extract を単独で実行するか、all に --dry-run を付けます。

pdb2reaction all -i 1.R.pdb 3.P.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' --dry-run

-c があると、--dry-run は一時ディレクトリでモデルを切り出し、電荷と、電子数が偶数か奇数かを確かめます。計算は行わず、出力も作りません。端末では次の行を見てください。

  • 原子数：Atoms after truncation: N（入力が複数なら [extract:multi] Atoms after truncation (model k): N）は、キャップ水素を足す前の原子数です。これにキャップ水素の数を足します。キャップ水素の数は、入力が 1 つなら [extract] Link-H to add: M、入力が複数なら [extract:multi] link-H targets common across models: M に出ます。上の 1.R.pdb と 3.P.pdb なら N + M = 669 です。

  • 電荷：Total active site model charge の行です。--dry-run では [all] --dry-run extract: model total_charge = … も出ます。

  • 警告：キャップ水素の付かない所で結合を切ったとき（all はここで止まります）と、アミノ酸らしい残基が 非標準の名前 のときに、extract は警告を出します。境界・キャップ・電荷を確かめてください。

-c を付けた all を実際に実行すると、出力ディレクトリの _work/models/model_<入力の名前>.pdb にモデルが残ります。ビューアで開き、反応に関わる残基が入っているかを確かめてください。

モデルを削る

DFT で最適化できるのは、多くても 300 原子くらいまでです。上の -c と -l で 1.R.pdb だけを切り出したときの、オプションごとの原子数と電荷を下の表に示します。原子数はキャップ水素を含みます。

モデルの作り方

オプション

原子数

電荷

デフォルトの切り出し

—

632

+1

主鎖の原子を除く

--exclude-backbone

397

+1

さらに水も除く

--exclude-backbone --no-include-h2o

367

+1

残基を自分で選ぶ

-r 0 --selected-resn '44,63,186'

129

−1

手で編集する

切り出した PDB を PyMOL で開き、要らない部分を消して保存する

—

—

  • --exclude-backbone：アミノ酸から主鎖の原子（N・CA・C・O とその水素）を除き、側鎖を CB で切ります。プロリンは環を残し、ペプチド結合でつながった -c のアミノ酸どうしの間では主鎖を残します。主鎖の原子しか中心の近くにない残基は選ばれません。

  • --no-include-h2o：水を除きます。

  • -r 0 --selected-resn：距離で残基を足さず、-c の残基と、指定した残基だけでモデルを作ります。同梱例の 44・63・186 番は、SAM のメチル炭素（CS1）にいちばん近い 3 残基です。自分の系では、反応に関わる残基を選んでください。

  • 手で削る：PyMOL などのビューアで要らない残基を消して保存し、自分で組んだモデルを使うの手順で使います。

all でも同じオプションを使えるので、最初から小さなモデルで MLIP の探索ができます。例は MLIP の TS を DFT で確かめる の例 1 にあります。残基を除くと電荷が変わるので、毎回ログを確かめてください。

モデルを広げる

反応に関わる残基がモデルに入っていないとき、または切れ端だけが入っているときは、モデルを広げてください。

  • -r を広げる（デフォルト 2.6 Å）：たとえば 3.5 Å や 4.0 Å にします。

  • --radius-het2het（デフォルト 0 で無効）：中心の側も相手の側も C・H 以外の原子どうしだけで測る、2 つ目の距離のしきい値を加えます。全体の半径を広げずに、近くの N・O の相手を拾えます。

  • 残基を -c に足す：残基を丸ごと残したいときに使います。-c のアミノ酸からも距離での探索が始まり、-r が 0 より大きければペプチド結合でつながった隣の残基が加わるので、--exclude-backbone なしでは主鎖を含む全部の原子が残ります。--selected-resn で足した残基は中心にならないので、モデルに隣の残基が無ければ側鎖だけになります。

  • 相手の鎖の残基が半径の外にあるときは、-c 'A:SAM,A:GPP,A:MG,B:MET:38' のように chain を付けて -c に足してください。

モデルを広げると計算は重くなり、精度が大きさとともに上がるとは限りません。モデルを広げても障壁が変わらないことを、系ごとに確かめてください。タンパク質の全体を環境に入れたいときは、ML/MM のツールキット mlmm-toolkit を使ってください。

自分で組んだモデルを使う

-c を省くと、all は切り出しを行わず、入力をそのまま使います。電荷は、全体を -q で渡すか、PDB/mmCIF の入力なら残基名ごとに -l で渡してください。

pdb2reaction all -i cluster_R.pdb cluster_P.pdb -l 'LIG:+1' --tsopt --thermo

切り出した PDB を削ったときは、残したキャップ水素の親原子はそのまま固定されます。消した残基のキャップ水素も消してください。親原子の無いキャップ水素があると isolated LKH/HL で止まります。新しく切った所には、下のチェックリストのとおりにキャップ水素を付けます。extract を使わずに組んだモデルには、自分で付けない限り LKH/HL のキャップ水素が無いので、--freeze-links は何も固定しません。境界の原子は --freeze-atoms で固定してください。

手でモデルを組む・確かめるとき

境界を自分で決めるときと、extract が組んだモデルを確かめるときは、次の点を見てください。

  • 切る結合は、距離だけでなく化学で選びます。

  • 主鎖の断片を残すときは、両端を CA でそろえ、キャップ水素で原子価を満たします。

  • 側鎖・リガンド・補因子の境界は、なるべく極性の無い C–C 単結合で切ります（ふつうは CA–CB）。ペプチドの C–N、極性の C–O・C–N、芳香環や共役の結合、S–S 結合、金属への配位結合は、半径の外に出ても切りません。相手をモデルに入れるか、境界を動かします。

  • キャップ水素は、切った結合の残る側の原子に 1 つずつ、元の結合の向きに沿って 1.09 Å の位置へ置きます。名前を残基 LKH（chain L）の原子 HL にすると、その親原子が自動で固定されます。

  • 境界の原子が意図した原子価を持ち、切った結合ごとにキャップ水素が 1 つずつあることを確かめます。

  • 反応物・中間体・生成物では、同じ原子を同じ順に並べます。選び方とキャップは 1 回決めてすべての状態に当て、状態ごとに別々に切り出しません。

  • 切った後は電荷と多重度を数え直し、最小エネルギー経路（MEP）や Hessian の計算の前に境界を目で見ます。

  • 共有結合した補因子、修飾残基、金属の部位、上の決まりに合わない境界は、自動の切り出しの後に手で直します。

切り出さずに酵素全体で計算する

クラスターモデルを切り出さずに、酵素全体の構造で反応経路を見ることもできます。このときは path-opt・path-search の GSM・DMF による最小エネルギー経路（MEP）の探索までで、虚振動 1 つの遷移状態（TS）は求められません。

原子の固定と距離の拘束

固定は、原子の力を 0 にしてその場から動かさないことです。拘束は、2 原子の距離を調和ばねで目標の値へ引くことです。切り出したモデルでは、キャップ水素の親原子を固定して、構造最適化・経路探索・固有反応座標（IRC）の計算の間に境界が崩れないようにします。

原子を固定する 3 つの方法

  • --freeze-links/--no-freeze-links（デフォルトで有効）：extract が付けたキャップ水素の親原子を固定します。PDB と mmCIF の入力で効きます。XYZ と GJF の入力には LKH の残基が無いので、同じ原子の PDB を --ref-pdb で渡すと効きます。座標は -i のものを使います。sp はこの原子を自動では固定しません。

  • --freeze-atoms 'i,j,k'：1 始まりの原子のインデックスで、どの形式の入力にも使えます。-c を付けた all では元の入力の番号で書き、切り出したモデルの番号に自動で写されます。ほかのコマンドでは、渡した構造の番号です。

  • YAML の geom.freeze_atoms（--config で渡す）：リストが長いときや、ほかの設定と一緒に残したいときに使います。

pdb2reaction extract -i complex.pdb -c 'A:GPP:301,A:SAM:302' -l 'GPP:-3,SAM:1' -o model.pdb
pdb2reaction opt -i model.pdb -l 'GPP:-3,SAM:1' -m 1   # キャップ水素の親原子が固定される
pdb2reaction tsopt -i ts_candidate.xyz -q 0 -m 1 --freeze-atoms '12,15,28,29,42'

geom:
  freeze_atoms: [12, 15, 28, 29, 42]   # 1 始まり

実行では 3 つの指定の和集合を固定します。3 つとも all・opt・tsopt・freq・irc・path-opt・path-search・scan・scan2d・scan3d で使えます。sp で使えるのは --freeze-atoms と YAML だけです。

固定の効果

  • opt・tsopt・scan・freq・irc と GSM（Growing String Method）の経路探索では、固定した原子の力を 0 にします。

  • 固定した原子を Hessian から外すので、freq は動ける原子で部分 Hessian 振動解析（PHVA）を行います。取り除く剛体運動は 固定境界での剛体モード にあります。

  • --mep-mode dmf の DMF（Direct Max Flux）の経路探索は、力を 0 にする代わりに調和拘束（k = 300 eV/Ų）で固定した原子を留めるので、座標が少し動くことがあります（path-opt の注意点）。

距離を拘束する

--distance-restraint は 2 原子のあいだに調和の拘束を足します。(i, j, 目標の距離)（Å）を渡すか、(i, j) で最初の距離を保ちます。

pdb2reaction opt -i input.pdb -q 0 -m 1 \
    --distance-restraint '[(1,5,2.0)]' --restraint-k 20.0 --out-dir ./result_opt_rest

  • 原子のインデックスは 1 始まりで、--zero-based を付けると 0 始まりになります。--restraint-k は力の定数です（デフォルト 300 eV/Ų）。

  • scan・scan2d・scan3d では、--restraint-k が各段のばねの強さです（距離は eV/Ų、角度は eV/rad²）。YAML の bias.k より CLI の値が優先されます。各段はその段の組だけを拘束します。all では --scan-restraint-k で同じ値を渡します。スキャンリストの書き方は スキャンリスト仕様 にあります。

  • opt の result.json（--out-json）とスキャンの結果のエネルギーは、拘束のエネルギーを除いた値です。

使用上の注意点

  • 300 原子は DFT の経験的な目安で、コードの上限ではありません。GPU のメモリは原子数だけでは決まらないので、同じ設定の短い試行で GPU メモリ を測ってください。

  • chain の欄：同梱の PDB は chain の欄が空です。chain が空の PDB では、残基を名前か番号で指定してください。chain のある PDB では、-c 'A:TYR:44' のように chain:残基名:番号 で書きます。

  • キャップ水素は CA と CB だけ：extract がキャップ水素を付けるのは、アミノ酸を CA か CB で切った所だけです。ほかの切断にはキャップ水素が付きません。

  • どの入力も同じ原子：原子の数や並びが違う入力は [multi] Atom count mismatch か [multi] Atom order mismatch で止まります。

  • --no-freeze-links は、境界をわざと緩める診断の実行のときだけ使ってください。本番の計算では --freeze-links を有効のままにします。

関連ドキュメント