ParmFit

Overview

ParmFit is the MAPLE force-field parameterization task. It turns MAPLE calculator results, existing topology files, and optional external QM references into Amber-style parameters for ligand/cofactor correction, noncanonical amino-acid templates, and metal-site parameterization.

ParmFit has three public method keywords. method=refinement refines parameters for an existing mol2 topology. method=ncaa parameterizes noncanonical amino-acid residues and method=metalaa parameterizes a metal site; both read a hydrogenated PDB block.

ParmFit calls AmberTools utilities to handle parameter assignment and template generation; the required tools are listed under Requirements. Before running MD with a ParmFit topology, follow each route’s deployment checks — the tleap validation must report Errors: 0, the fitted charges must match their targets, and MAPLE-typed topologies go through a ParmEd pass before SHAKE.

Routes

The method keyword selects the route; without an explicit method (and without input=), ParmFit defaults to refinement. The input for each route — its header keys and the required structure block (a coordinate block for Refinement, a PDB block for NCAA and MetalAA) — is specified in that route’s manual. The ParmFit Tutorials work through each route on a real system.

Requirements

All routes need a working MAPLE installation and AmberTools on PATH — a run invokes whichever utilities it needs: parmchk2, antechamber, prepgen, espgen, resp, and tleap. Refinement relies on GAFF2 atom typing: supply a pre-generated mol2 (antechamber -at gaff2 ...), or omit mol2= and ParmFit does the typing itself at runtime, building the mol2 from the molecule recorded in the input file. That generated mol2 exists only to carry the GAFF2 types — its charge column holds zero placeholders, and chg_fit therefore defaults to abcg2 (see the Refinement page). In both cases the mol2 atoms must appear in the same order as the coordinate block, hydrogens included — the auto-generated mol2 satisfies this by construction; MAPLE checks only that the atom counts match and then pairs atoms by position. NCAA and MetalAA read exactly one PDB block and need no separate atomic-charge input. Refinement may preserve the charges already stored in the input mol2 when chg_fit=none; the other charge routes derive charges at runtime (see ChargeFit below).

When RESP charges are needed. They always take their ESP from Gaussian: g16 or g09 must be on PATH for Refinement and NCAA runs with chg_fit=resp (the NCAA default) and for every MetalAA run. QM reference jobs (iqm=true) use qm_engine, which defaults to g16; set qm_engine=orca to run them with ORCA instead. A Refinement or NCAA run therefore needs no Gaussian when chg_fit is not resp and the run keeps iqm=false or sets qm_engine=orca; MetalAA always needs Gaussian.

How MAPLE reads PDB structures. The PDB must be explicitly hydrogenated; MAPLE neither adds nor removes hydrogens. While reading the PDB, MAPLE builds a connectivity graph and matches standard protein residues against the active pro_ff Amber templates using element composition, topology, explicit hydrogens, and chirality. Matched residues receive canonical atom names, atom types, partial charges, residue net charges, and Amber connectivity — the atom names are rewritten to the canonical template names before tleap export — so standard input atom names are not required. Cysteine states are resolved in two steps: explicit SG–H connectivity distinguishes protonated CYS from the deprotonated CYM/CYX candidates, and the graph-equivalent CYM/CYX pair is decided by disulfide-bond detection — a cysteine participating in an S–S bond is assigned CYX (net charge 0), otherwise CYM (net charge −1). Alternate conformations resolve by default to the altloc A set, or to the highest-occupancy set when A is absent; keep chooses among them when a listed residue carries several alternate conformations — append the altloc letter (A181B) to build that residue from its altloc-B atoms, while a letterless entry (A181) keeps the default choice. Residues without a complete standard template are reported, not guessed: backbone-only matches (NCAA candidates) and unmatched nonprotein residues, which should be declared via lig_resids/ncaa_resids in MetalAA or parameterized in separate NCAA runs.

Read the tleap validation warnings. Errors: 0 is the gate; the warnings decide whether the topology is safe to use. Most are routine — terminal-name conversions (NTHR -> THR) and an addIons request after neutralization. The tutorials list and explain the warnings of their verified runs.

General Settings

These parameters apply across all ParmFit routes. For iqm behavior, qm_mode, and levels of theory (opt_level / sp_level / chg_level), see the QM Reference section below.

ChargeFit

chg_fit selects the atomic-charge source independently of iqm and qm_compare. Refinement defaults to none; the input mol2 charges are kept as they are. NCAA defaults to resp; none is invalid because the PDB route has no reusable input charges. MetalAA derives site charges through a fixed RESP-only workflow and does not expose chg_fit.

RESP derivation always runs the same two-tool pipeline — Gaussian for the ESP, AmberTools for the RESP stages — whether iqm is enabled or disabled; its Gaussian jobs reuse the configured qm_nproc and qm_mem. The charge-capable MLIP may differ from the main #model. If it matches the main model, ParmFit reuses the active calculator and its configured options. A different charge-capable MLIP is initialized once on the active device with its default options.

TorsionFit

TorsionFit refines proper torsion parameters in two stages. Stage 1 uses a spectral shared-group restrained linear least-squares solver: torsion instances that share the same chemical environment (atom types plus one-hop neighbors) are grouped and fitted together, with a spectral helper providing phase seeds. Stage 2 takes the Stage 1 parameters and performs direct global kPhi/phase optimization against the scan reference profile.

torsion_refine_rounds controls how many Stage 1→Stage 2 fast-MM cycles are attempted. Each cycle refreshes the MM reference profiles — the first from the Stage 0 parameters (the mSeminario-refined starting point), later ones from the previous round’s accepted parameters — and re-runs Stage 1 and Stage 2. A cycle is accepted only if total loss improves; otherwise the workflow rolls back and stops. With torsion_refine_rounds=0, only Stage 1 runs.

TorsionFit applies to Refinement and NCAA routes. The fitted dihedrals differ by route: Refinement auto-selects rotatable single bonds from the mol2 topology; NCAA restricts fitting to sidechain dihedrals and disables torsion ensemble terms internally.

Scan Reuse

Before each scan, TorsionFit looks under <base>_work/torsionfit/ for the scan file that an earlier run wrote for the same target. If that file has the expected number of frames and the same elements, and its scan angles match those of the current starting geometry, TorsionFit reads and reuses it, and the log does not report the reuse. The check ignores the model and the QM level of theory, so delete <base>_work/torsionfit/ before you rerun with either of them changed.

QM Reference (iqm)

iqm controls whether external QM jobs provide reference-quality geometry, Hessian, and torsion data. When iqm=false, all optimization, Hessian evaluation, and torsion reference work use the MAPLE MLIP calculator exclusively. When iqm=true, MAPLE runs QM reference jobs whose scope depends on the route.

Levels of Theory

Level strings are passed to the QM engine verbatim, so write them the way the selected engine expects: B3LYP/def2SVP for Gaussian, B3LYP/def2-SVP for ORCA. A mismatched spelling — for example the hyphenated form with qm_engine=g16 — fails at the first QM job with a route syntax error. RESP chg_level values always follow Gaussian syntax. The built-in defaults of opt_level and sp_level are written B3LYP/def2-SVP, the ORCA spelling, while the default engine is g16. Set opt_level explicitly (for example B3LYP/def2SVP) before enabling iqm with Gaussian.

TorsionFit qm_mode

qm_mode controls how QM reference data is used inside TorsionFit (Refinement and NCAA; mode 1 is the default). It has no effect on MetalAA.

qm_compare (Refinement only)

When iqm=true, qm_compare=true, and refinement is enabled, Refinement runs an additional MLIP-only comparison branch alongside the QM path. The MLIP branch produces _mlip-suffixed output files, so its parameters can be compared against the QM-referenced ones without re-running the workflow. Charge fitting uses the same configured chg_fit method independently on the QM and MLIP branch geometries, so each branch reports its own charged work mol2.

External ParmFit Config

For more complex jobs, use a separate key-value text file and reference it from the input header. Any plain text file with key=value lines and optional # comments is accepted. The config file accepts parmfit parameters only; do not place global MAPLE keys in it. Declare the method inline together with the input= reference, for example #parmfit(method=ncaa,input=...).

#model=uma(size=uma-s-1p1)
#parmfit(method=ncaa,input=NCAA.parmfit)

PDB MaDA.pdb

Output Layout

ParmFit writes a main .out file and a <base>_work/ directory. Next to the .out file, every route also writes a <base>.parmfitrun.in snapshot of the resolved configuration (user entries with defaults filled in); unrecognized keys are reported as warnings in the log. The excerpts below illustrate key sections of the .out file for each route.

Refinement

The setup header echoes the input mol2, bonded refinement method, QM reference state, and TorsionFit configuration:

======================================================================
                       Parmfit Refinement Setup
======================================================================
Input mol2:        toluene_ff.mol2
Initial frcmod:    Toluene_work/Toluene_original.frcmod
Topology atoms:    ...
Topology bonds:    ...
Bonded refinement: mSeminario
QM reference:     disabled
TorsionFit:        enabled
Torsion backend:   cgbs
Torsion bonds:      auto-select non-ring torsion bonds with proper torsions

The mSeminario bond/angle changes table lists each modified BOND and ANGLE term with its type code (atom pairs for BOND, atom triples for ANGLE) and old → new values:

------------------------------------------------------------------------------------------------------------
                                                   BONDS
------------------------------------------------------------------------------------------------------------
(1,2)            c3-ca             k=...  r=...                       |  k=...  r=...
...

------------------------------------------------------------------------------------------------------------
                                                   ANGLES
------------------------------------------------------------------------------------------------------------
(2,1,8)          ca-c3-hc                k=...  theta=...             |  k=...  theta=...
...

The TorsionFit dihedral changes table shows per-torsion kPhi/n/phase changes across multiple term slots; NA entries indicate newly added terms:

------------------------------------------------------------------------------------------------------------
                                                 DIHEDRALS
------------------------------------------------------------------------------------------------------------
(3,2,1,8)         ca-ca-c3-hc       term=1 k=...  n=1  phase=...           |  k=...  n=1  phase=...
(3,2,1,8)         ca-ca-c3-hc       term=2 k=NA   n=NA phase=NA            |  k=...  n=2  phase=...
...

The torsion energy trace gives per-dihedral angle-versus-energy profiles and ends with MAE/RMSE. Its columns are the scan reference (MLIP or QM), the original GAFF2 profile (orig), the mSeminario-refined starting point (stage0), and the refitted profiles after each TorsionFit stage (stage1, stage2). A rejected Stage 2 remains visible as an NA column; the final comparison is labelled stage2_final when Stage 2 is accepted and stage1_final when it is rejected:

torsion bond (1, 2):
  angle_deg   MLIP_ref   orig_ref   stage0_ref   stage1_ref   stage2_ref
       ...          ...        ...          ...          ...           NA
  ...
MLIP_ref vs stage1_final: MAE = ... kcal/mol, RMSE = ... kcal/mol

The final PARMFIT REFINEMENT RESULT block lists the products, charge provenance, stage timing, and parameter-change counts:

Main products:
  Amber mol2:   Toluene_work/Toluene_maple.mol2
  Amber frcmod: Toluene_work/Toluene_maple.frcmod
  ...

Charge fitting:
  method:        resp
  RESP level:    HF/6-31G(d)
  target charge: ...
  actual charge: ...
  work MOL2:     Toluene_work/chargefit/mlip_route_charged.mol2
  final MOL2:    Toluene_work/Toluene_maple.mol2

NCAA

The job header identifies the target residue, chirality, charge method, and progress through each stage:

Target: _89:HIF (protein)
  PDB residue templates: matched=516, backbone_only=1, unmatched=0
    non-standard (matched backbone only): _89:HIF (target)
  PDB target charge: _89:HIF = 0 (from your input)
  [NCAA] MLIP model preparation + reference optimization ...
  NCAA target selector: HIF89
  residue name: HF1
  chirality: L
  net charge / spin multiplicity: 0 1
  protein ff: ff14SB
  bonded refinement: mseminario
  charge fitting: resp
  charge level: HF/6-31G(d)
  [NCAA] atomic charge fitting (resp) ...
  [NCAA] AmberTools template build ...
  [NCAA] MLIP Hessian + mSeminario ...
  [NCAA] TorsionFit skipped ...
  [NCAA] writing refined templates + tleap input ...
  [NCAA] tleap validation ...

The PARMFIT NCAA RESULT block confirms the charge result and lists the output files, tleap status, and the exact tleap command to run:

Status: completed
Route: NCAA
Target: _89:HIF
Residue name: HF1
Chirality: L
Protein ff: ff14SB
Res charge/spin-multi: 0 1
Charge fitting: resp
Charge level: HF/6-31G(d)
Target charge: 0
Actual charge: ...
Charge MOL2: MaDA_work/ncaa/HIF89/MaDA_capped.mol2

Main products:
  refined prepin: MaDA_work/HF1_maple.prepin
  refined frcmod: MaDA_work/HF1_maple.frcmod
  tleap PDB:      MaDA_work/MaDA_ncaa_tleap.pdb
  tleap input:    MaDA_work/MaDA_ncaa_tleap.in

Next step:
  cd MaDA_work
  tleap -s -f MaDA_ncaa_tleap.in |tee MaDA_ncaa_tleap.out

Tleap status:
  Errors: 0
  Warnings: 4
  Notes: 0

MetalAA

The job header prints the metal target, site parameters, ion charge and multiplicity, and per-stage logs:

Target: A462:FE (ion)
  ion_resids: A462
  add_resid: []
  cluster_cutoff: 3.00 A
  donor_cutoff: 2.70 A
  water ff: tip3p
  protein ff: ff14SB
  ion ff: 12_6
  ion_charges: 3
  ion_mults: 6
  ncaa_resids: []
  large model charge/mult: -2 6
  chgmod: 1
  fixchg_resids: []
  QM method: PBE1PBE/def2SV
  [MetalAA] large-model input written.
  [MetalAA] MLIP large-model optimization ...
  [MetalAA] large-model RESP ...
  [MetalAA] site files written.
  [MetalAA] MLIP Hessian + Seminario + final frcmod ...
  [MetalAA] final parameter files written.
  [MetalAA] running tleap validation ...

The renamed atom types table shows how MAPLE reassigns types locally: the metal receives M1, direct donors receive Z1/Z2/..., and other atoms keep their original Amber/GAFF types. Verify that the donor list matches the intended coordination environment:

Renamed atom types:
  residue      atom  old   new   charge
  _400:CYM     SG    SH    Z1    -0.329172
  A461:HEM     NA    nc    Z2    -0.129187
  ...
  A462:FE      FE    FE    M1     0.439642

The PARMFIT METALAA RESULT block lists core residues, output files, tleap status, any environment warnings, and the tleap command:

Status: completed
Route: MetalAA
Target: A462:FE
Protein ff: ff14SB
Ions charge: 3
Large model charge/mult: -2 6
Core residues: _400:CYM, A461:HEM, A462:FE

Main products:
  final frcmod: P450_work/P450_metal.frcmod
  tleap input:  P450_work/P450_metal_tleap.in
  tleap PDB:    P450_work/P450_metal_tleap.pdb
  mol2 files:   P450_work/CM1.mol2, P450_work/HM1.mol2, P450_work/FE1.mol2

Next step:
  cd P450_work
  tleap -s -f P450_metal_tleap.in |tee P450_metal_tleap.out

Tleap status:
  Errors: 0
  Warnings: 8
  Notes: 1

  Warnings:
  Potentially charged standard residues within 8 Å of the metal but outside the large model: _398:ARG

Troubleshooting