Ligand Parameterization

Ligand Parameterization

Ligand parameterization is Core 2. The challenge is not the individual steps (Antechamber, parmchk2, RESP fitting are all well-established) but the bookkeeping between them — the same molecule is represented three or four different ways, and any inconsistency in atom identity drops the wrong charges on the wrong atoms. This page walks the pipeline end-to-end and shows where the gates are.

The four representations of one molecule

Representation	Source	What it carries	What it lacks
SDF	Vendor / curator (e.g. PubChem, the medicinal chemist)	bond graph, bond orders, aromaticity, formal charge, ideal geometry, often hydrogens	the bound conformation in your complex
Complex PDB ligand	X-ray / docking	the bound conformation, atom names, chain ID, residue name	bond orders, formal charge (PDB only stores element + coordinates)
Typed mol2	Antechamber output	GAFF2 atom types, force-field-ready geometry	RESP charges (only initial AM1-BCC or Mulliken)
Final mol2	After RESP injection	typed + RESP-charged, ready for tleap	—

The pipeline reconciles these. The non-trivial part is making sure the same physical atom ends up at the same name across all four files.

Pipeline

   SDF ─────┐                                                  
            ▼                                                  
       Antechamber  ──> typed_ligand.mol2     (GAFF2 types, no RESP yet)
            │                                                  
            ▼                                                  
       parmchk2     ──> ligand.frcmod         (bonded + vdW parameters)
            │                                                  
   PDB ─────┤                                                  
            ▼                                                  
       graph_match ──> φ : SDF atom → PDB atom    (identity-safe)
            │                                                  
            ▼                                                  
       coord_embed ──> heavy-atom coords from PDB,                       
                      H coords from SDF/GAFF2 template          
            │                                                  
            ▼                                                  
       PySCF / GPU4PySCF: HF/6-31G* wavefunction               
            │                                                  
            ▼                                                  
       Multiwfn RESP fit ──> q_L (per-atom partial charges)    
            │                                                  
            ▼                                                  
       charge-sum gate: |Σq - Q_formal| ≤ 1e-4                 
            │                                                  
            ▼                                                  
       inject q_L into typed_ligand.mol2 → ligand.mol2         
            │                                                  
            ▼                                                  
       name-based merge μ_L : mol2 ↔ complex-PDB ligand        
            │                                                  
            ▼                                                  
       tleap loadmol2 + frcmod ──> ligand block in topology    

Identity safety: the graph-isomorphism step

This is the single biggest difference from a naive "trust the row order" pipeline.

Strict isomorphism

Build the SDF graph $G_{\mathcal{T}} = (\mathcal{T}, E_{\mathcal{T}})$ from explicit bonds and the PDB graph $G_{\mathcal{C}} = (\mathcal{C}, E_{\mathcal{C}})$ from coordinate-derived bonds (covalent-radius cutoffs). Solve for the edge-preserving map

$$ \phi : \mathcal{T} \to \mathcal{C}, \qquad {(\phi(a), \phi(b)) : (a,b) \in E_{\mathcal{T}}} = E_{\mathcal{C}}, $$

constrained by matching element and graph degree at each vertex. If a unique solution exists, $\phi$ is determined and the pipeline proceeds without ambiguity.

Hypergraph fallback

For highly symmetric ligands (substituted benzenes, symmetric heterocycles), strict isomorphism returns multiple equally valid mappings. CYPForge 1.1.0 then invokes a heavy-atom hypergraph resolver that ranks candidate $\phi$'s by:

1. element match,
2. graph-neighborhood signature (degree + element of neighbors),
3. pairwise heavy-atom distance terms,
4. Kabsch superposition residual against the complex pose.

The resolver returns one of three decisions:

• unique — one candidate clearly wins; admitted.
• equivalent_ok — the remaining ambiguities are indistinguishable under the retained checks (e.g. two equivalent rotations of a phenyl ring), or are explicitly recorded equivalent exchanges; admitted as an audit boundary and recorded as a WARN-level note.
• anything else — rejected; the workflow halts.

equivalent_ok is honest: the implementation cannot prove a strict uniqueness theorem in those cases, so it records the situation rather than pretending. Theorem 11.1 in S3 explicitly notes this is an acceptance boundary, not a proof of uniqueness.

Why this matters

A wrong $\phi$ means the RESP charges (computed in SDF atom order) get injected into the wrong PDB atoms. The charge sum stays correct, but per-atom assignment is scrambled. In a symmetric ligand this can be invisible; in an asymmetric one it produces wrong electrostatics on the substrate's reactive site. Row-index matching is the failure mode CYPForge is built to prevent.

Coordinate embedding

Once $\phi$ is fixed, the typed mol2 is populated:

• Heavy atoms — coordinates copied from the complex PDB via $\phi$. The typed mol2 inherits the bound conformation.
• Hydrogens — rebuilt from the SDF + GAFF2 template geometry. The complex PDB usually doesn't have ligand hydrogens (or has them inconsistently), and rebuilding them from the SDF preserves the chemistry. An optional pre-RESP "H-only" PBE relaxation can be requested (pre-resp-relax pbe-h-only) to clean up hydrogen positions without moving any heavy atom.

After this step the typed mol2 has the same chemistry as the SDF, the same bound conformation as the PDB, and is ready for charge injection. A gate then requires:

$$ \mathrm{RMSD}_{\text{heavy}}(\text{typed\_mol2}, \text{complex-PDB-ligand}) \leq 0.05\ \mathrm{\AA}. $$

If the embedding shifted heavy atoms beyond this, something went wrong — most likely a misidentified $\phi$.

RESP charges

CYPForge uses the standard HF/6-31G* / RESP recipe:

1. Wavefunction — gpu4pyscf (CUDA-accelerated PySCF) computes an HF/6-31G* wavefunction at the embedded geometry. CPU pyscf is supported as a fallback when no GPU is available.
2. ESP grid — written via Multiwfn-compatible format.
3. RESP fit — Multiwfn_noGUI reads the wavefunction, builds the ESP grid, and performs the two-stage RESP fit with the standard hyperbolic restraints.
4. Charge-sum gate — the fitted $q_L$ must satisfy

$$ \left| \sum_{n \in \mathcal{N}_L} q_L(n) - Q_L \right| \leq 10^{-4}, $$

where $Q_L$ is the integer formal charge declared at init time (--formal-charge). A failure here usually means the wavefunction's nuclear charge does not match the declared formal charge, i.e. you declared --formal-charge -1 for a neutral ligand or vice versa.

The charges are then written back into the typed mol2, producing the final ligand.mol2. GAFF2 atom types and bonded parameters are not touched at this stage.

Atom-name merge

When LEaP loads the ligand, it does so by name. The merge map $\mu_L$ identifies the mol2 atom name with the complex-PDB atom name:

$$ \mathcal{N}_{\text{mol2}}^L = \mathcal{N}_{\text{PDB}}^L \implies \mu_L(n) = n. $$

This is unique only if both name sets are identical and duplicate-free. If they differ (which happens when the SDF / Antechamber pipeline invented new H names that the PDB does not have), the workflow records a coordinate/element fallback and emits a WARN; the user is expected to review.

Vina-assisted mode-11

When the ligand pose in your complex PDB came from a docking run rather than an X-ray structure, an additional frame compatibility check is required. Let $H_d$ and $H_c$ be the same-name heme anchor atoms in the docking receptor and the current assembly receptor:

$$ \mathrm{RMSD}(H_d, H_c) = \sqrt{\frac{1}{k} \sum_{i=1}^k |\mathbf{h}_i^d - \mathbf{h}_i^c|^2} \leq 0.25\ \mathrm{\AA}. $$

This gate confirms that the docking pose and the assembly pose are in the same reference frame. It does not validate that the docking pose itself is a chemically reasonable binding mode — that is your responsibility.

What this layer guarantees

If Core 2 reaches PASS:

• the ligand atom identity is either strictly determined (graph isomorphism unique) or admitted under an audit boundary (equivalent_ok recorded),
• per-atom RESP charges are deposited on the correct atoms,
• the charge sum matches the declared formal charge to $10^{-4}$,
• the heavy-atom geometry in the typed mol2 is within $0.05\ \mathrm{\AA}$ of the complex PDB,
• the mol2-to-PDB name merge is independent of file order.

What it does not guarantee:

• that the docking / X-ray pose is a chemically reasonable binding mode,
• that per-atom RESP charges are chemically optimal (only the sum is gated),
• that GAFF2 atom typing is the best choice for your particular chemistry — if GAFF2 fails on your ligand (e.g. organometallic centers), you must supply your own force field.

Where this lives in the code

Step	Source
Antechamber + parmchk2	src/cypforge_core/ligand_gpu4pyscf_esp.py
Graph isomorphism	src/cypforge_core/ligand_mapping_resolver.py
Hypergraph fallback	src/cypforge_core/ligand_heavy_hypergraph_resolver.py
Coordinate embedding	src/cypforge_core/ligand_gpu4pyscf_esp.py
HF/6-31G* + RESP	src/cypforge_core/ligand_gpu4pyscf_esp.py, src/cypforge_core/_qm/qm_restrained_ligand_opt.py
Charge-sum gate	src/cypforge_core/ligand_gpu4pyscf_esp.py
Atom-name merge	src/cypforge_core/ligand_mapping_leapin.py
Vina pose frame check	src/cypforge_core/ligand_pose_frame.py, src/cypforge_core/ligand_vina_pose_selection.py

Next: Gate System and Manifests — how the PASS / WARN / FAIL machinery is implemented.