all（一気通貫ワークフロー）

all サブコマンドは、活性部位モデルの抽出と最小エネルギー経路（MEP）の探索を 1 回の実行で行います。指定すれば、各反応段の遷移状態（TS）の最適化と、固有反応座標（IRC）・振動数・DFT の計算まで行います。

--tsopt を付けない場合は、TS 候補までで終わります。TS 候補は各 MEP セグメントで最もエネルギーの高い点（HEI）です。デフォルトの計算バックエンドは、Meta が公開した学習済みの機械学習原子間ポテンシャル（MLIP）である UMA です。

主な用途

与える入力でモードが決まります。

  • R と P から経路とエネルギー図を作る: 反応順に並べた 2 構造以上（反応物、中間体、生成物）を与えると、隣り合う構造の間の MEP を求め、エネルギー図を描く

  • 反応物 1 つから経路を作る: 1 構造と、作る結合・切れる結合を -s で与えると、段階的スキャンで中間体を作り、それらを通る MEP を求める

  • TS 候補 1 つを確かめる（TS-only モード）: 1 構造に --tsopt を付け、-s を付けずに与えると、TS を最適化し、そこから IRC をたどる（n_imag = 1 で、IRC が狙った R と P に着けば TS と確かめられる）

基本的な実行例

例は GPP C6-メチル基転移酵素 BezA（Tsutsumi et al., Angew. Chem. Int. Ed. 2022, 61, e202111217）の系で、スクリプト一式は examples/ にあります。1.R.pdb（反応物）・2.IM.pdb（中間体）・3.P.pdb（生成物）は、水素原子をすべて含む全系の構造です。自分の構造にも水素原子が要ります。例 1〜3 の流れと結果の確かめ方は クイックスタート: all、クイックスタート: --scan-lists、クイックスタート: TS-only モード にあります。

1. MEP を求め、TS 最適化・熱化学・DFT まで計算する

-c で抽出の中心を、-l で非標準残基の電荷を指定します。

pdb2reaction all -i 1.R.pdb 3.P.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' \
    --tsopt --thermo --dft --out-dir ./result_mep

端末に各 TS の [tsopt] Converged (n_imag=1). と、最後の ====== Pipeline summary ====== の下の Scientific status: success が出れば、指定した段はすべて終わっています。result_mep/summary.json にも同じ値が入ります。続けて 実行結果の判定 のとおり端点を確かめてください。最適化した構造は result_mep/segments/seg_NN/ にあります。

2. 反応物から段階的スキャンで経路を作る

ステージ 1 で SAM のメチル炭素（CS1）を GPP の C7 に近づけ（1.60 Å）、ステージ 2 で GPP の H11 を Glu186 の OE2 に移します（0.90 Å）。

pdb2reaction all -i 1.R.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' \
    -s '[("CS1 SAM 320","GPP 321 C7",1.60)]' '[("GPP 321 H11","GLU 186 OE2",0.90)]' \
    --tsopt --thermo --out-dir ./result_scan

1 つのリテラルの中の目標は、同じステージで一緒に動きます。リテラルを並べると順に別のステージとして実行し、各ステージは前のステージの終わりの構造から始まります。各ステージの終わりの構造が MEP 探索の入力になります。-s は 1 回だけ書き、その後にすべてのリテラルを並べてください。反応の分け方は 反応の分け方を決める を参照してください。chain が空の PDB では、原子を残基名・残基番号・原子名の 3 つで順不同に指定します（"CS1 SAM 320"）。chain があるときは A:SAM:320:CS1 と書きます。指定できる形はすべて スキャンリスト仕様 にあります。

3. TS 候補を確かめる（TS-only モード）

入力を 1 つにして --tsopt を付け、-s を付けないと、MEP 探索を省きます。

pdb2reaction all -i TS_candidate.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' \
    --tsopt --thermo --dft

最適化した R・TS・P は result_all/segments/seg_01/ に書き出されます。

4. 途中のセグメントから後処理をやり直す

セグメント N から後処理をやり直すときは、元のコマンドと同じ入力・抽出・経路・計算バックエンドのオプションと同じ --out-dir を指定し、--resume-segment N を足してください。--tsopt-max-cycles などの後処理のオプションは変えられます。

pdb2reaction all -i 1.R.pdb 3.P.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' \
    --tsopt --thermo --dft \
    --resume-segment 1 --out-dir ./result_mep

N より前のセグメントはそのまま残し、セグメント N 以降の後処理と、要約・エネルギー図を書き直します。

処理の仕組みと計算仕様

全系の構造 (PDB / mmCIF / XYZ / GJF)
  ├─ (-c のとき) 活性部位モデルの抽出: extract
  │   └─ 活性部位モデル
  ├─ (1 構造 + -s のとき) 段階的スキャン: scan
  │   └─ 各ステージの終わりの構造 = 中間体
  ├─ MEP 探索: path-opt（デフォルト）または path-search（--refine-path）
  │   └─ mep_trj.xyz と energy_diagram_MEP.png
  └─ (--tsopt のとき) TS 最適化と IRC: tsopt → irc
      ├─ (--thermo のとき) 振動数と熱化学: freq
      └─ (--dft のとき) DFT 一点計算: dft

  1. 入力の準備: altloc（代替配置）を含む PDB では、残基ごとに平均の占有率が最も高いラベルを 1 つ選びます。元素の欄が空の PDB では、元素記号を補います。-c を指定すると、指定した残基のまわりの活性部位モデルを切り出し、切った結合をキャップ水素で埋めます。

  2. 経路の作成: まず入力構造を最適化します（--preopt）。-s を指定すると、段階的スキャンで中間体を作ります。続いて path-opt が、隣り合う構造の間の MEP を GSM（growing string method）か DMF（direct max flux）で求めます。--refine-path では、再帰的な path-search が経路を精密化し、結合が変わる所で反応段に分けます。各反応段の HEI がその段の TS 候補です。

  3. TS の最適化（--tsopt）: 各 HEI をデフォルトでは RS-P-RFO（restricted-step partitioned rational function optimization）で最適化し、最後の Hessian から n_imag を求めます。

  4. IRC の追跡: TS から EulerPC（Euler 予測子–修正子法）で IRC を両方向へたどり、両端を極小まで最適化します。これがそのセグメントの R と P になります。

  5. 熱化学と DFT: --thermo では R・TS・P で freq を実行してギブズエネルギーを求め、--dft では同じ構造で DFT 一点計算を行います。それぞれのエネルギー図も描きます。

all が TS から IRC へ進むのは、TS 最適化が収束し、最後の Hessian を計算でき、n_imag ≥ 1 のときだけです。

TS の結果

all の次の動作

収束、n_imag = 1

IRC を実行し、IRC の両端を最適化します。

収束、n_imag ≥ 2

警告を出し、MEP の方向に最もよく合う虚振動（合うものが無ければ最も低い虚振動）に沿って IRC を実行します。結果は partial です。

収束、n_imag = 0

IRC の前で止まります。

未収束（サイクル上限か --stop-plateau）、--skip-final-freq（最後の Hessian を省く）、Hessian の失敗

IRC の前で止まります。

IRC の前で止まった場合、結果は success になりません。TS のファイルは segments/seg_NN/ts/ に残り、後のセグメントの後処理は行いません。TS 最適化の終わり方の一覧は tsopt の「TS の判定」 にあります。

片方の端点の最適化が収束しない場合、結果は partial になり、segments/seg_NN/endpoint_opt/ を確認用に残します。端点の最適化がエラーで失敗した場合は、エラーを segments/seg_NN/endpoint_opt/failure.json に記録し、そのセグメントは振動数と DFT の段の前で止まります。TS と IRC の構造は残ります。

実行結果の判定

TS 最適化が成功すると、反応モードの虚振動が 1 つ出ます（n_imag = 1）。IRC が収束しなくても、端点の最適化で狙った R と P に着けば、その結果は使えます。

結果は次の 3 か所で確かめます。

  • 端末: 各 TS 最適化は、虚振動 1 つで収束すると [tsopt] Converged (n_imag=1). で終わります。====== Pipeline summary ====== の下に Execution status: と Scientific status: が出ます。結果が success でないときは、RESULT WARNING: の行に理由が出ます。

  • summary.log: ヘッダーに Pipeline mode（MEP、Scan、TS-only）と 2 つのステータスが出ます。[1] は MEP の概要、[2] は各セグメントの MEP 上の障壁 ΔE‡・反応エネルギー ΔE・結合の変化、[3] は各セグメントの後処理で、TS imaginary freq: の下に n_imag が出ます。[4] はエネルギー図の表、[5] は出力のツリーです。

  • summary.json: scientific_status には success・partial・failed のいずれかが、scientific_status_reasons にはその理由が入ります。各 TS の n_imag は post_segments[].tsopt.n_imaginary_modes です。

    • summary.json の障壁: --tsopt のとき、各セグメントの障壁は post_segments[].mlip.barrier_kcal で、最適化した TS と R の MLIP のエネルギー差です。--thermo と --dft のときは、同じ barrier_kcal が gibbs_mlip・dft・gibbs_dft_mlip の下にもあります。segments[].barrier_kcal は TS 最適化の前の MEP 上の障壁で、TS-only モードでは TS − R です。

success は、指定した段がすべて収束したことを示し、--tsopt のときはすべての TS が n_imag = 1 であることも意味します。端点が狙った R と P かは自分で確かめてください。summary.log の [2] の結合の変化と、segments/seg_NN/reactant.*・product.* の構造を、狙った R と P と比べます。n_imag が 1 でないときや、端点が狙いと違うときは TS が取れないとき を参照してください。

主な出力ファイル

all は --out-dir に次のファイルを書き出します。

result_all/
├─ summary.log                  # 結果の要約（テキスト）
├─ summary.json                 # 機械可読な結果（常に出力。all に --out-json はありません）
├─ mep_trj.xyz                  # 全セグメントの MEP 軌跡
├─ mep_trj.pdb                  # 同じ軌跡の PDB
├─ mep_trj.cif                  # 同じ軌跡の mmCIF（mmCIF 入力または大きな PDB 入力のとき）
├─ mep_w_ref.pdb                # MEP を全系の入力に重ねた構造（--write-ref-merge）
├─ energy_diagram_MEP.png       # 全セグメントの MEP のエネルギー
├─ energy_diagram_*_all.png     # 全セグメントの R → TS → P の図（--tsopt、--thermo、--dft）
├─ irc_plot_all.png             # 全セグメントの IRC のエネルギー（--tsopt）
├─ segments/
│  └─ seg_NN/                   # 反応の 1 段: seg_01, seg_02, ...
│     ├─ reactant.*             # 最適化した R・TS・P（入力と同じ形式、--tsopt）
│     ├─ ts.*
│     ├─ product.*
│     ├─ energy_diagram_*.png   # この段の R → TS → P の図
│     ├─ ts/                    # TS 最適化。vib/imag_*_trj.xyz は虚振動のアニメーション
│     ├─ irc/                   # IRC の軌跡と irc_plot.png
│     ├─ endpoint_opt/          # 端点の最適化（--dump のとき、または端点が収束しなかったか失敗したときに残る）
│     ├─ freq/{R,TS,P}/         # 振動数と熱化学（--thermo）
│     └─ dft/{R,TS,P}/          # DFT 一点計算（--dft）
└─ _work/                       # 途中のファイル（TS 候補の HEI を含む）
   ├─ models/                   # 抽出したモデル model_<入力名>.pdb（-c のとき）
   ├─ scan/                     # 段階的スキャン（-s のとき）
   └─ path_opt/                 # MEP 探索と hei_seg_NN.*（--refine-path のときは path_search/）

  • 報告に使う構造: segments/seg_NN/reactant.*・ts.*・product.* を使ってください。seg_NN/ の下の各ディレクトリには、各段の計算のファイルが入っています。

  • TS-only モード: MEP 探索が無いので、MEP のファイルと _work/path_opt/ はありません。R・TS・P は segments/seg_01/ に入ります。

エネルギー図のファイル名は手法を表します。

ファイル名

生成されるとき

内容

energy_diagram_MEP.png

MEP 探索の完了時

全セグメントの MEP のエネルギー

energy_diagram_MLIP.png

--tsopt

R → TS → P、MLIP のエネルギー

energy_diagram_G_MLIP.png

--thermo

R → TS → P、MLIP のギブズエネルギー

energy_diagram_DFT.png

--dft

R → TS → P、MLIP の構造での DFT のエネルギー

energy_diagram_G_DFT_plus_MLIP.png

--dft と --thermo

R → TS → P、DFT のエネルギーに MLIP の熱補正を足した値

energy_diagram_*_all.png

_all の無い図と同じ

全セグメントをまとめた同じ図（出力ディレクトリの直下）

irc_plot.png（seg_NN/irc/ の中）、irc_plot_all.png

--tsopt

1 つのセグメントと全セグメントの IRC のエネルギー

図のエネルギーは、最初の状態（反応物）を基準にした kcal/mol です。

主な CLI オプション

オプション

引数の型

デフォルト

説明

-i, --input

パス（複数可）

（必須）

反応順に並べた 2 構造以上、または -s か --tsopt を付けた 1 構造（.pdb、.cif、.xyz、.gjf）。1 つの -i の後に並べるか、-i を繰り返す

-c, --center

文字列

None

抽出の中心（通常は基質と触媒残基）。残基名（'SAM,GPP'）、残基 ID（'A:123,B:456'）、chain 付きの名前（'A:SAM'、'A:SAM:123'）。省略すると入力全体を使う

-l, --ligand-charge

文字列

None

非標準残基の電荷（例: 'SAM:1,GPP:-3'）、またはその合計（リガンドの総電荷）の数値。PDB/mmCIF 入力のみ

-q, --charge

整数

None

系全体の総電荷。-c のときは抽出したモデルから求める。-c が無いときは、-l、YAML の calc.charge、.gjf 入力のどれも無ければ必須。明示すると求めた値より優先し、警告を出す

-m, --multiplicity

整数

1

スピン多重度（2S+1）

-b, --backend

文字列

uma

計算バックエンド（uma, orb, mace, aimnet2, dft）

-r, --radius

浮動小数点数

2.6

中心原子からの抽出の半径（Å）。0 では -c と --selected-resn の残基だけを残す

--selected-resn

文字列

""

半径による拡張なしで入れる残基（-c と同じ形）

-s, --scan-lists

文字列

None

1 構造の段階的スキャンの目標。1 つのリテラルが 1 ステージ（例: '[("A:SAM:320:CS1","A:GPP:321:C7",1.60)]'。書き方は スキャンリスト仕様）

--tsopt/--no-tsopt

フラグ

False

各セグメントの TS を最適化し、IRC を実行

--thermo/--no-thermo

フラグ

False

R・TS・P の振動数と熱化学（--tsopt が必要）

--dft/--no-dft

フラグ

False

R・TS・P の DFT 一点計算（--tsopt が必要）

--refine-path/--no-refine-path

フラグ

False

隣り合う組ごとの path-opt の代わりに、再帰的な path-search を実行

--mep-mode

gsm / dmf

gsm

MEP の手法: GSM または DMF

--opt-mode

grad / hess

grad

単一構造の最適化とスキャンのオプティマイザ: grad = L-BFGS、hess = RFO

--opt-mode-post

grad / hess

hess（コマンドラインに --opt-mode を書いたときはその値）

TS と IRC 後の端点のオプティマイザ: grad = TS は Dimer・端点は L-BFGS、hess = TS は RS-P-RFO・端点は RFO

--preopt/--no-preopt

フラグ

True

スキャンと MEP 探索の前に入力構造を最適化

--flatten/--no-flatten

フラグ

False

TS 最適化の後に残った余分な虚振動を消す

--stop-plateau/--no-stop-plateau

フラグ

False

収束の前にエネルギーが変わらなくなったら最適化を止める。収束ではなく stalled として報告する

--tsopt-max-cycles

整数

100000

TS 最適化のサイクル上限

--resume-segment

整数

None

--out-dir の MEP を使い、セグメント N から後処理をやり直す（例 4）

--dry-run/--no-dry-run

フラグ

False

計算をせずにオプションを確かめ、計画を表示。-c のときは一時ディレクトリで抽出を行い、電荷を確かめる

-o, --out-dir

パス

./result_all/

出力先ディレクトリ

全オプションは pdb2reaction all --help-advanced か 自動生成のオプションの一覧（英語のみ） を参照してください。

補足: YAML（--config）では、上の表のオプションに無い設定もできます。節とキーの一覧は YAML 設定の一覧 にあります。

使用上の注意点

  • --dft と -b dft: 一緒には使えず、実行の始めにエラーで止まります。-b dft の計算の後に DFT の一点計算を足すときは、別のジョブで pdb2reaction sp -b dft か pdb2reaction dft を実行してください。

  • --dft の費用: 必要なメモリは構造・基底・汎関数・精度・ソフトウェアの構成で変わります。対象の計算ノードで代表的な構造を試し、最大メモリ使用量を見てください。大きなモデルでは、MLIP の計算を先に終え、DFT の一点計算を別のジョブで実行してください。

  • TS-only モードの R と P: IRC のエネルギーの高いほうの端を反応物と呼びます（同じなら左の端）。R・P の名前、ファイル名、障壁、反応エネルギーはこのエネルギーの順に従うもので、化学的に分かった反応の向きではありません。P からの障壁は barrier_kcal − delta_kcal です。summary.json の endpoint_assignment にこの規則が記録され、chemical_direction_known: false となります。

  • TS-only モードの summary.log: [1] は TS と IRC の概要で、[2] は最適化した TS と端点から求めます。

  • 熱化学のファイル: all は熱化学量を thermoanalysis.yaml から読むので、--thermo では --no-dump を指定してもこのファイルを残します。

  • 抽出の半径: -r 0 では半径による拡張を無効にし、-c と --selected-resn で選んだ残基からモデルを組みます。構造上の安全策として、ジスルフィド結合の相手や隣の残基の主鎖が加わることはあります。半径 0 は内部で 0.001 Å として扱います。

  • -c を省いたとき: 抽出を行わず、入力構造の全体を MEP 探索・tsopt・freq・dft に渡します。1 構造のときは、このときも -s か --tsopt が必要です。

  • 入力の形式: -c を使うときは PDB か mmCIF が必要です。-c が無いときは XYZ と GJF も使えます。1 回の実行のすべての構造は、同じ原子を同じ順に持つ必要があります。

  • 電荷と多重度: -c のときの総電荷は抽出したモデルの合計で、アミノ酸・イオン・水は組み込みの値、そのほかの残基は -l の値、-l に無い残基は 0 として数えます。-c が無いときは、入力に -l を当てて求めるか、.gjf のヘッダーから読みます。多重度は -m、無ければ .gjf のヘッダー、それも無ければ 1 です。詳しくは 電荷の指定 を参照してください。

  • 別々に用意した構造: 入力構造を別々に用意すると、反応座標の外の構造の違いも障壁に入ります。障壁を読む前に構造を比べてください。組成が同じ 2 つの機構を比べるときは、両方の経路で共通の原子の集合と順序を使います。

  • --write-ref-merge: 経路を元の全系の入力に重ねた構造を、確認用に書き出します。mep_w_ref* は出力ディレクトリの直下に、hei_w_ref_seg_NN.pdb は _work/path_search/ に入ります。--refine-path、-c、PDB か mmCIF の入力が必要です。

  • --resume-segment: --tsopt・--thermo・--dft のどれかが必要で、--dry-run とは一緒に使えません。保存した入力と MEP がコマンドと合わないときは、エラーで止まります。

変異体と野生型の比較

1 つの経路の中では、すべての構造が同じ原子を同じ順に持ちます。変異体と野生型（WT）では残基が違い、原子数も変わることが多いので、両者の全エネルギーをそのまま差し引くことはできません。代わりに、それぞれの系の中で求めた障壁を比べます。

ΔΔG‡ = (G_TS − G_R)_mutant − (G_TS − G_R)_WT

  • 2 つのモデルで、選ぶ残基の位置と、境界・キャップの決め方をそろえ、狙った変異だけが違うようにします。半径で別々に抽出すると、境界の残基が片方のモデルにだけ入ることがあるので、2 つの選択を比べてください。

  • プロトン化の決め方、電荷の決め方、バックエンドとモデル、精度、拘束、熱化学の条件をそろえます。変異でプロトン化の状態や形式電荷が変わる場合は総電荷も違うので、両方に同じ -q を当てはめないでください。

2 つの実行は、入力と出力先のほかは同じオプションにします。R が化学的な反応物になるよう、それぞれの系の R と P を与えます（Endpoint モード）。G_TS − G_R は post_segments[].gibbs_mlip.barrier_kcal です。

pdb2reaction all -i wt_R.pdb wt_P.pdb -c 'SAM,GPP,MG' -l 'GPP:-3,SAM:1' --tsopt --thermo -o result_wt
pdb2reaction all -i mutant_R.pdb mutant_P.pdb -c 'SAM,GPP,MG' -l 'GPP:-3,SAM:1' --tsopt --thermo -o result_mutant

関連ドキュメント