QM/MM

Run a quantum region inside a classical environment: a ligand in a protein pocket, a reacting site in an enzyme, a solute in explicit solvent. The QM atoms get a real wavefunction; everything else becomes point charges that polarize it.

What you get: QM region selection (by index, by radius, or by whole residue), link atoms across covalent cuts with four boundary-charge schemes, Gaussian-smeared charges, Thole polarizable embedding, an AMBER-form MM crate, analytic QM and MM forces, a full-structure gradient across the cut, geometry optimization.

What is missing: no AMBER prmtop reader (go through OpenMM), and no periodic boundary conditions.

Before you start: relax the geometry with xtb

DFT on an unrelaxed structure wastes hours and often fails.

Measured on a 27-atom carbocation built by hand:

hand-built    C7-H24 = 0.902 Å      SCF: Stalled@37, never converged
xtb GFN2      C7-H24 = 1.084 Å      SCF: converged
                                    545 kcal/mol lower

An SCF has no reason to converge on a structure with a 0.9 Å C–H bond. xtb fixes it in seconds; DFT spends hours failing on it.

xtb input.xyz --opt --gfn 2 --chrg 1        # seconds
ferric single_point.toml                     # then the expensive step

Use it for the initial geometry of anything you did not get from a crystal structure or a previous optimization: hand-built ligands, docked poses, edited residues, anything from a SMILES string.

Build note: xtb built with -O3 miscompiles gradients on this platform. Build with -Doptimization=2.

The smallest working example

import ferric

# The full structure: every atom, QM and MM alike.
symbols = ["O", "H", "H",   "O", "H", "H"]     # two waters
coords  = [[0.0, 0.0, 0.0], [0.76, 0.59, 0.0], [-0.76, 0.59, 0.0],
           [0.0, 0.0, 3.0], [0.76, 0.59, 3.0], [-0.76, 0.59, 3.0]]
charges = [-0.834, 0.417, 0.417,  -0.834, 0.417, 0.417]   # MM partial charges (e)

# QM = the first water; the second becomes point charges.
sys = ferric.QmmmSystem(symbols, coords, charges, qm_indices=[0, 1, 2])

res = ferric.run_qmmm(sys, "cc-pvdz", method="rhf")
print(res.energy)

charges are ignored for atoms that end up in the QM region, so you can pass a full force-field charge set and let the selection decide.

Selecting the QM region

Three ways:

# 1. Explicit indices — full control.
ferric.QmmmSystem(symbols, coords, charges, qm_indices=[0, 1, 2])

# 2. Everything within a radius of some seed atoms.
ferric.QmmmSystem(symbols, coords, charges,
                  qm_seeds=[0], qm_radius_angstrom=5.0)

# 3. Whole residues — a residue joins if any of its atoms is in range.
#    Prefer this for proteins: a radius cut alone will slice through a
#    residue and leave you with a chemically meaningless fragment.
ferric.QmmmSystem(symbols, coords, charges,
                  qm_seeds=[0], qm_radius_angstrom=5.0,
                  residue_ids=residue_ids)

residue_ids together with qm_indices is an error: an explicit index list and a residue expansion are contradictory instructions.

Cutting a covalent bond

If the QM/MM boundary crosses a bond, you need a link atom and a scheme for what to do with the host's charge:

sys = (ferric.QmmmSystem(symbols, coords, charges, qm_indices=qm)
       .with_link_atoms(bonds)                        # bonds crossing the cut
       .with_boundary_charges(bonds, "rcd"))          # keep | delete-host | rc | rcd

The parser is strict: an unknown string is an error.

Energies, forces, and optimization

res = ferric.run_qmmm(sys, "cc-pvdz", method="uhf")
# KS-DFT: method="rks" or "uks" and xc, e.g. method="rks", xc="PBE".
# xc without rks/uks (or rks/uks without xc) raises ValueError.
res.energy
res.qm_gradient()      # dE/dR on the QM atoms
res.mm_forces()        # forces on the MM sites
res.full_gradient()    # across the cut, link rows folded back onto the frontier

opt = ferric.run_optimize_qmmm(sys, "cc-pvdz", method="rhf", move_mm="none")

move_mm chooses which MM atoms relax: "none" (the default, where only QM atoms move), ("within", radius_angstrom), ("residues", [ids]) (the system must have been built with residue_ids=), or "all". Any value other than "none" requires mm_topology=: without a force field, nothing holds the moving MM atoms' geometry together.

Polarizable embedding

Thole polarizable sites, if fixed charges are not enough:

sys = ferric.QmmmSystem(symbols, coords, charges,
                        qm_indices=qm,
                        polarizabilities_angstrom3=alphas)
res = ferric.run_qmmm(sys, "cc-pvdz", method="rhf", thole_a=2.1304)
res.e_pol                  # the polarization energy
res.induced_dipoles()

What is validated, and how far

Checked against pyscf.qmmm.mm_charge: the energy shift agrees to <1e-8 Ha, MM forces to 3e-10. An empty or all-zero MM region is bit-identical to a gas phase calculation.

The QM–MM Lennard-Jones term across a covalent cut is checked against OpenMM's own exclusion and 1-4 handling (ethanol cut at C–C, realistic LJ): energy and gradient agree to ~1e-18 Ha. With no bond path between the regions it is bit-identical to the plain every-pair sum.

Not validated: periodic boundary conditions (absent).

From the CLI

A [qmmm] section runs QM/MM from a TOML file, with no Python needed. The QM region becomes the molecule that is solved, and the MM region becomes the external potential it is solved in. A complete, runnable input is examples/water-qmmm.toml (repository, not the wheel):

[molecule]
xyz = "testdata/molecules/water.xyz"   # required by the parser, not read with [qmmm]
charge = 0                             # apply to the QM region
multiplicity = 1

[basis]
name = "sto-3g"

[method]
kind = "rhf"
task = "energy"

[qmmm]
pqr = "testdata/molecules/water_na.pqr"
qm_indices = [0, 1, 2]            # zero-based; or: qm_seeds = [0], qm_radius_angstrom = 1.5
# link_bonds = [[0, 3]]           # required when the cut crosses a covalent bond
# boundary_scheme = "delete-host" # default; also "keep", "rc", "rcd"

Measured on that file (water QM, one Na⁺ 4 Å away as MM, STO-3G): vacuum −74.9629466809, embedded −74.9653197421 Ha (−1.489 kcal/mol). That matches ferric.run_rhf(point_charges=...) to all ten printed digits.

The geometry and the MM charges both come from the PQR. An xyz has no partial charges, and an MM region without charges contributes nothing. [molecule] xyz must still be present, because the parser requires it, but it is not read when [qmmm] is present. That matters when you compute a vacuum reference: delete [qmmm] and ferric falls back to the xyz. If the xyz holds a different geometry, you are comparing two different molecules. In the example above, the xyz is an optimized water that gives −74.9631468000, which is 0.13 kcal/mol of error in the difference. Compare at one geometry.

boundary_scheme defaults to "delete-host", not "keep". Keeping the host charge across a covalent cut puts a bare point charge inside the link atom's bond length (measured 0.443 Å on an ethane C–C cut). A geometry optimization in that field diverges instead of failing with an error. "keep" is still available and is the right choice when the cut isn't covalent. An unknown scheme name is an error.

[qmmm] and [external_potential] together are refused. The MM region is an external potential, so combining them would count a contribution twice.

Scope of the CLI section: it reads only a PQR, because it needs charges and geometry together, and does electrostatic embedding only. Polarizable sites, smeared charges, the MM force field and MM relaxation are Python-only.

Known limits

• No AMBER prmtop reader. Go through OpenMM, or build the arrays.
• No periodic boundary conditions.
• The pocket field is fixed unless you ask otherwise. move_mm="none" is the default in run_optimize_qmmm; the MM sites do not relax with the QM region until you widen it.
• QM/MM dispersion is the force field's Lennard-Jones only. D3/D4/XDM/VV10 are QM-atom-pairwise and never see the MM atoms. With mm_topology= set, the MM force field supplies Lennard-Jones between every QM real atom and every MM atom (link atoms and boundary charges carry none). Without it there is no QM–MM van der Waals term at all.
• QM–MM Lennard-Jones follows the force field's exclusion rules across the cut. Pairs 1-2 or 1-3 apart through the topology's bond list (the cut bond included) get no LJ, and 1-4 pairs are scaled by the topology's scale_lj_14 (AMBER default 0.5). The MM–MM terms use the same bond graph, so the cut bond has to be in the topology's bonds. QM–MM electrostatics do not use these exclusions or scale_coul_14. The embedding puts every MM charge in the QM Hamiltonian, and the boundary scheme (keep/delete-host/rc/rcd) handles the charges next to the cut.
• The MM crate is AMBER-form (harmonic bonds and angles, periodic torsions, Lennard-Jones, Coulomb), validated against OpenMM.