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 two public method keywords. method=correction refines parameters for an existing mol2 topology. method=abinitio reads a PDB target and automatically routes to NCAA or MetalAA according to the target residue kind.

Parmfit calls AmberTools utilities such as parmchk2, antechamber, prepgen, resp, espgen, and tleap to handle parameter assignment and template generation.

Rebuild topology with parmed

Parmfit reassigns atom types (MAPLE Z0/Z1, M1/Y1, etc.) on deployed residues. After type reassignment, tleap may not correctly populate the %FLAG BONDS_INC_HYDROGEN topology section for H-bond constraints. Use parmed to rebuild the topology:

parmed <base>.prmtop
outparm <base>_parmed.prmtop
quit

parmed automatically determines hydrogen-mass constraint groups from atomic masses in the topology.

Routes

Route	Public Input	Target	Main Products
Correction	`#parmfit(method=correction,mol2=...)`	Ligand, cofactor, or small molecule with existing mol2 topology	Refined Amber mol2/frcmod and GROMACS top/gro
NCAA	`#parmfit(method=abinitio,target=...)`	Protein residue target in a PDB file	Residue prepin/frcmod, processed PDB, tleap validation files
MetalAA	`#parmfit(method=abinitio,target=...)`	Ion target in a PDB metal site	Metal-site mol2, metal frcmod, tleap deployment files

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.

Parameter	Default	Description
`qm_engine`	`g16`	External QM engine: `g16`, `g09`, or `orca`
`qm_nproc`	`8`	Number of processors for QM jobs
`qm_mem`	`24`	Memory in GB for QM jobs; adjust this according to the available memory and the size of the selected QM reference model.
`opt_max_iter`	`256`	Maximum LBFGS iterations for MLIP geometry optimization. In Correction, use `max_iter`
`opt_max_step`	`0.2`	Maximum LBFGS step size in Å for MLIP geometry optimization. In Correction, use `max_step`

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 uses the accepted parameters from the previous round to refresh the MM reference profiles, then re-runs Stage 1 and Stage 2. A cycle is accepted only if total loss improves; otherwise the workflow rolls back and stops. Setting torsion_refine_rounds=0 runs Stage 1 only.

TorsionFit applies to Correction and NCAA routes. Center-bond selection is route-specific: Correction auto-selects rotatable single bonds from the mol2 topology; NCAA restricts fitting to sidechain torsions and disables torsion ensemble terms internally.

Parameter	Default	Description
`torsionfit`	`true`	Enable proper torsion fitting after the bond/angle stage
`torsion_bonds`	Auto	Explicit center bond atom pairs (e.g. `6-8 12-14`); otherwise auto-selected from mol2 topology
`torsion_steps`	`36`	Number of scan increments per torsion
`torsion_refine_rounds`	`3`	Maximum Stage 1→Stage 2 fast-MM cycles
`torsion_refine_max_iter`	`256`	Maximum iterations for the Stage 2 global optimizer
`torsion_refine_tol`	`1.0e-6`	Convergence tolerance for Stage 2
`torsion_ensemble`	`false`	Add rigid rotor ensemble frames as extra Stage 2 constraints in Correction; NCAA disables this internally
`torsion_ensemble_ratio`	`0.3`	Target ensemble frame ratio relative to center bonds × steps in Correction
`torsion_ensemble_weight`	`0.5`	Weight of ensemble loss relative to scan loss in Correction Stage 2
`backend`	`cgbs`	TorsionFit scan optimizer backend
`constraint_mode`	`projected`	Scan constraint mode: `projected` or `fixinternals`

QM Reference (iqm)

iqm controls whether external QM jobs provide reference-quality geometry and Hessian 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.

iqm does not affect RESP

Abinitio routes always run Gaussian ESP jobs for RESP charge fitting regardless of the iqm setting. iqm controls only the reference geometry, Hessian, and torsion-scan data, not the RESP charge pipeline.

Route	`iqm=false`	`iqm=true`
Correction	MLIP opt + MLIP Hessian (+ MLIP TorsionFit if enabled)	QM opt_frequency → QM geometry + QM Hessian (+ QM-referenced TorsionFit via `qm_mode`). If `bonded=none` but `torsionfit=true`, only a QM geometry optimization is run—no Hessian is computed.
NCAA	MLIP capped-model opt + MLIP Hessian	QM opt_frequency on the capped model → QM geometry + QM Hessian, but only when `bonded≠none` or `torsionfit=true`. If both are disabled, `iqm=true` has no effect.
MetalAA	MLIP large-model opt + MLIP Hessian	QM opt_frequency after MLIP opt → QM geometry + QM Hessian for bonded fitting. No TorsionFit, no `qm_mode`.

Levels of Theory

Parameter	Correction	NCAA	MetalAA
`chg_level`	—	RESP ESP (always active)	RESP ESP (always active)
`opt_level`	QM opt/opt_freq/constrained_opt when `iqm=true`	QM opt/opt_freq when `iqm=true`	QM opt_freq when `iqm=true`
`sp_level`	QM single-point via `qm_mode` when `iqm=true`; defaults to `opt_level`	QM single-point via `qm_mode` when `iqm=true`; defaults to `opt_level`	Unused (no TorsionFit)

TorsionFit qm_mode

qm_mode controls how QM reference data is used inside TorsionFit (Correction and NCAA). It has no effect on MetalAA.

Mode	Description
1	MLIP drives the torsion scan, then QM single-point energy and forces are calculated at each frame
2	MLIP pre-optimization, then QM constrained torsion optimization at each frame
3	QM engine evaluates energy and forces at each scan step; the built-in optimizer of MAPLE drives the scan

qm_compare (Correction only)

When iqm=true, qm_compare=true, and refinement is enabled, Correction runs an additional MLIP-only comparison branch alongside the QM path. This produces _mlip-suffixed output files that let users compare MLIP-driven parameters against QM-referenced ones without re-running the workflow.

Input Patterns

Correction

Correction uses the inline coordinate block as the structure to optimize and fit. The mol2 file supplies the starting atom types, charges, and parameterization context.

#model=uma(size=uma-s-1p1,task=omol)
#parmfit(method=correction,mol2=toluene_ff.mol2)
#device=gpu0

C   0.000000  1.520000  0.000000
...

Abinitio

Abinitio routes require exactly one PDB block. The target selects a residue in that PDB; MAPLE classifies it as a protein residue or ion and dispatches to NCAA or MetalAA.

#model=uma(size=uma-s-1p1)
#parmfit(method=abinitio,target=HIF89,rn=HF1,cmo=0 1)
#device=gpu0

PDB MaDA-1.pdb

External Parmfit Config

For longer 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 # or ; comments is accepted. The config file accepts parmfit parameters only; do not place global MAPLE keys in it.

#parmfit(method=abinitio,input=NCAA.parmfit)

PDB MaDA-1.pdb

Note

An external config must not set pdb=, input=, or frcmod=. Provide PDB structures through the PDB <path> block of the input file; correction frcmod files are generated automatically.

Output Layout

Parmfit writes a main .out file and a <base>_work/ directory. The excerpts below illustrate key sections of the .out file for each route.

Correction

The setup header confirms the input mol2, topology size, bonded refinement method, and TorsionFit configuration:

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

The mSeminario bond/angle changes table lists each modified BOND and ANGLE term with its atom-pair type code and old → new force constants:

------------------------------------------------------------------------------------------------------------
                                                   BONDS
------------------------------------------------------------------------------------------------------------
(1,2)            c3-ca             k=243.910000  r=1.515000           |  k=259.042910  r=1.501978
(1,8)            c3-hc             k=345.250000  r=1.096000           |  k=357.895134  r=1.088538
(2,3)            ca-ca             k=354.250000  r=1.399000           |  k=386.871870  r=1.388923

------------------------------------------------------------------------------------------------------------
                                                   ANGLES
------------------------------------------------------------------------------------------------------------
(2,1,8)          ca-c3-hc                k=43.610000  theta=110.6300  |  k=53.192800  theta=110.8773
(2,1,9)          ca-c3-hc                k=43.610000  theta=110.6300  |  k=55.257488  theta=110.9240
(1,2,3)          c3-ca-ca                k=60.740000  theta=120.8300  |  k=97.095122  theta=120.8023

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=0.000000  n=1  phase=0.0000   |  k=0.326863  n=1  phase=0.0094
(3,2,1,8)         ca-ca-c3-hc       term=2 k=NA  n=NA  phase=NA            |  k=0.339285  n=2  phase=-179.9939
(3,2,1,8)         ca-ca-c3-hc       term=3 k=NA  n=NA  phase=NA            |  k=0.292462  n=3  phase=107.5106
(3,2,1,8)         ca-ca-c3-hc       term=4 k=NA  n=NA  phase=NA            |  k=0.223706  n=4  phase=-0.1078

The torsion energy trace gives per-center-bond angle-vs-energy profiles, ending with MAE/RMSE. The columns orig_ref (original GAFF2) and stage2_ref (final fitted) should be compared against MLIP_ref to confirm fitting improvement:

center bond (1, 2):
  angle_deg   MLIP_ref   orig_ref   stage0_ref   stage1_ref   stage2_ref
     0.0035   0.014655   0.082197    0.021962    0.019633    0.019663
    10.0035   0.015672   0.066117    0.012318    0.011285    0.010432
     ...

The final PARMFIT CORRECTION RESULT block lists all output files, stage timing, and a summary of how many bonds, angles, and dihedrals were changed.

NCAA

The abinitio header identifies the target residue, chirality, and charge, and reports progress through each stage:

Target: _89:HIF (protein)
  [NCAA] MLIP model preparation + reference optimization ...
  NCAA target selector: HIF89
  residue name: HF1
  chirality: L
  charge/mult: 0 1
  protein model: ff14SB
  bonded refinement: mseminario
  QM ESP method: HF/6-31G(d)
  [NCAA] multiconformer RESP ...
  [NCAA] AmberTools template build ...
  [NCAA] MLIP Hessian + mSeminario ...
  [NCAA] TorsionFit skipped ...
  [NCAA] writing refined templates + tleap input ...
  [NCAA] tleap validation ...

The PARMFIT ABINITIO RESULT block confirms the output files, stage timing, tleap status, and the exact tleap command to run:

Status: completed
Route: NCAA
Target: _89:HIF
Residue name: HF1
Chirality: L
Representative conformer: ref
RESP conformers: alpha, beta

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

Stage timing:
  multiconformer RESP              9m 21s   slowest
  model preparation + reference    1m 29s
  Hessian + mSeminario             1m 12s
  tleap validation                     3s

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

Tleap status:
  Errors: 0
  Warnings: 3
  Notes: 0

MetalAA

The abinitio header prints the metal target, site parameters, large-model charge, and per-stage logs:

Target: A462:FE (ion)
  cluster_cutoff: 3.00 A
  donor_cutoff: 2.70 A
  metal site charge/mult: -1 6
  metal oxidation: 3
  large model charge/mult: -2 6
  chgmod: 1
  [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] running tleap validation ...

The renamed atom types table shows how MAPLE reassigns types locally: the metal receives M1, direct donors receive Y1/Y2/..., 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    Y1    -0.307250
  A461:HEM     NA    nd    Y2    -0.065018
  A461:HEM     NB    nc    Y3    -0.127778
  A461:HEM     ND    nd    Y4    -0.092271
  A461:HEM     NC    nc    Y5    -0.056100
  A462:FE      FE    FE    M1     0.293950

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

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

Stage timing:
  large RESP/Gaussian ESP       30m 35s   slowest
  Hessian + Seminario               55s
  large optimization                25s
  tleap validation                   4s

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:
  Environment contains potentially charged standard residues
  outside the site core: _3:LYS, _4:GLU, ...