Quickstart: pdb2reaction all --scan-lists

pdb2reaction all --scan-lists (-s) builds a reaction path from one structure. It drives the distances you choose to target values under harmonic restraints (a scan) and uses the scan endpoints to search the minimum energy path (MEP). With --tsopt, it continues to transition-state (TS) optimization and an intrinsic reaction coordinate (IRC) calculation. You define the reaction coordinate yourself, so the rest of the model changes less and the intended TS is easier to reach. The starting structure can be a reactant, an intermediate, or a product.

The commands below use 1.R.pdb from the bundled examples/; run them in that directory.

What it is for

  • No product structure: drive the forming and breaking bonds to make the product (P) from the reactant (R).

  • Two-step reactions: drive one step after another, such as a methyl transfer followed by a proton transfer.

  • MEP and TS from the scan: go on from the scanned path to the MEP and the TS in the same run.

Choosing a scan command

Goal

Command

Restrained structures and a scan trajectory only

pdb2reaction scan

Go on from the scan to the MEP, and optionally to TS and IRC

pdb2reaction all -s ...

A 2D or 3D energy map over two or three coordinates

pdb2reaction scan2d / scan3d

Examples

1. Check the input first

--dry-run checks the input, the model extraction, the charge and spin parity (electron count vs. multiplicity), and the mapping of the scan atoms to the model, and then stops without any calculation.

pdb2reaction all -i 1.R.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' \
  -s '[(4360, 4419, 1.60)]' --dry-run

The check has passed when the console shows [all] --dry-run parity check OK: ... and [all] Planned stages: extract -> scan -> path_opt., and ends with [Dry run] --dry-run completed. Input command is valid. When the check fails, the run stops with an error such as --dry-run parity check failed; look up the message in Troubleshooting.

2. One stage

Drive the distance between SAM CS1 (atom 4360) and GPP C7 (atom 4419; C6 in IUPAC numbering) to 1.60 Å. Give the atoms by number or by name; the two commands are the same scan.

# Atom numbers (1-based by default)
pdb2reaction all -i 1.R.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' \
 -s '[(4360, 4419, 1.60)]' -o ./result_scan

# Atom names
pdb2reaction all -i 1.R.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' \
 -s '[("SAM,320,CS1", "GPP,321,C7", 1.60)]' -o ./result_scan

The run succeeded when the ====== Pipeline summary ====== block near the end of the console shows Scientific status: success; summary.json holds the same value in scientific_status.

3. Two stages

Each bracketed list after -s is one stage, called a literal (a list written out as text), and the stages run in order:

# Stage 1: drive the methyl-transfer distance to 1.60 Å
# Stage 2: then drive the proton transfer to 0.90 Å
pdb2reaction all -i 1.R.pdb -c 'SAM,GPP,MG' -l 'SAM:1,GPP:-3' -s \
  '[("SAM,320,CS1","GPP,321,C7",1.60)]' \
  '[("GPP,321,H11","GLU,186,OE2",0.90)]' \
  -o ./result_scan

This short example drives only the bonds that form. In your own reaction, put the bonds that break and the H atoms that move in the same stage too; see Decide how to split the reaction.

Writing --scan-lists

  • Tuples: each literal is a list of tuples: distance (atom1, atom2, target_Å), angle (atom1, atom2, atom3, target_deg), or dihedral (atom1, atom2, atom3, atom4, target_deg). Distances are in ångströms, and angles and dihedrals in degrees.

  • Atoms: give 1-based integer atom numbers in the order of the full input structure (--scan-zero-based for 0-based), or atom names in double quotes. With -c, both are mapped to the extracted model.

  • Stages: the tuples in one literal move together, and each stage starts from the relaxed result of the previous one. Give -s once and list the literals after it; all stops with an error when -s is repeated.

  • Quoting: wrap each literal in single quotes so that the shell leaves the parentheses and spaces alone.

  • Bundled PDBs: they have an empty chain column, so use the name or number forms such as "SAM,320,CS1", not the forms with a chain ID.

For the full syntax of atom names and quoting, see Scan-list spec.

Output files

The commands above write:

result_scan/
├── summary.log
├── summary.json                 # Results, with scientific_status
├── mep_trj.pdb                  # MEP over all segments
├── energy_diagram_MEP.png       # MEP energy profile
└── _work/                       # Intermediate files, including the HEI (TS candidate); kept after the run
    ├── scan/
    │   ├── preopt/              # Optimized starting structure
    │   ├── stage_01/            # Scan stage 1
    │   │   ├── result.{xyz,pdb} # Restrained endpoint (optimized without restraints only with --scan-endopt)
    │   │   ├── scan_trj.xyz     # Scan trajectory
    │   │   └── scan.pdb
    │   ├── stage_02/            # Scan stage 2 (two-stage run)
    │   └── result.json          # Scan results of each stage
    └── path_opt/                # MEP search (path_search/ with --refine-path, the recursive MEP search)
        └── hei_seg_01.{xyz,pdb} # Highest-energy image of segment 1

These commands stop after the MEP search and do not create segments/. Add --tsopt for the R/TS/P structures and IRC of each reactive segment in segments/seg_NN/, and --thermo for freq/.

Checking the result

  1. Completion: scientific_status is success when every requested stage converged; otherwise it is partial or failed, with the reasons in scientific_status_reasons. With --tsopt, two checks are left for you: that the imaginary mode moves the bonds that form or break, and that the endpoints are the intended R and P.

  2. Scan: open _work/scan/stage_01/scan_trj.xyz in a viewer and check that the distances change as intended. At the end of each stage the console prints [stage 1] Covalent-bond changes (start vs final): Yes or No, and _work/scan/result.json records it in stages[].bond_changes.

  3. MEP: open mep_trj.pdb and _work/path_opt/hei_seg_01.pdb, the highest-energy image (HEI) and TS candidate, and check that energy_diagram_MEP.png shows a clear barrier.

  4. TS (with --tsopt): a successful TS optimization gives one imaginary mode along the reaction coordinate. The console then prints [tsopt] Converged (n_imag=1)., and summary.json records the count in post_segments[].tsopt.n_imaginary_modes. Open segments/seg_01/ts/vib/imag_*_trj.xyz in a viewer and check that the mode moves the bonds that form or break.

  5. Endpoints (with --tsopt): open segments/seg_01/irc/finished_irc_trj.xyz and the optimized endpoints segments/seg_01/reactant.pdb and product.pdb, and check that they are the intended R and P. Even if the IRC does not converge, the result is usable when the endpoint optimizations reach the intended R and P.

Notes

  • Input: PDB/mmCIF, XYZ, or GJF. -c (PDB/mmCIF only) cuts out a cluster model from the full system; without it the structure is used as is.

  • Defaults of all and scan: the two commands share the scan engine, but the option names and the defaults for the optimizations before and after the scan differ:

    Command

    Step / restraint

    Relaxation limit

    Optimization before / after the scan

    pdb2reaction all

    --scan-max-step-size 0.20 Å; --scan-restraint-k 300 eV/Ų

    --scan-relax-max-cycles 100000

    before: on (--preopt/--no-preopt); after: off (--scan-endopt/--no-scan-endopt)

    pdb2reaction scan

    --max-step-size 0.20 Å; --restraint-k 300 eV/Ų

    --relax-max-cycles 100000

    before: off (--preopt/--no-preopt); after: off (--endopt/--no-endopt)

    The restraint strength can also be set with YAML bias.k.

  • Bond changes and --refine-path: every scan endpoint goes on to the MEP search, whether or not its stage shows a bond change. --refine-path alone decides whether the MEP is refined recursively (path_search/).

  • Scan endpoints: a finished scan gives restrained structures. They are not minima or transition states until an optimization without restraints, or a TS optimization and IRC, confirms them.

  • scan on its own: the standalone scan also accepts a YAML or JSON spec file and scan ranges; all -s takes inline target tuples only.

  • scan --dry-run: it checks the input, the charge and spin parity, and the parse of --scan-lists, but not the extraction or atom mapping of all. Use all --dry-run for those.

Next steps