[qmmm]

The [qmmm] section configures hybrid quantum-mechanics/molecular-mechanics (QM/MM) calculations. A quantum region described by an OpenQP method is embedded in a classical force-field environment handled by OpenMM. The established .inp surface includes single-point QM/MM energy and ground-state QM/MM molecular dynamics, while NAMD uses the separate embedded surface-hopping driver.

Ground-state md(...) qmmm(...) always selects the OpenMM QMMM_MD engine, including nonperiodic cutoff=NoCutoff calculations. Periodic methods (PME, Ewald, and CutoffPeriodic) require that OpenMM path, and any future spherical-boundary QM/MM implementation must use it as well. Only an all-QM gas-phase md(...) calculation without qmmm(...) uses the native OpenQP velocity-Verlet driver. In QM/MM NAMD, OpenMM supplies the environmental force model and constraints, while OpenQP continues to perform electronic propagation, surface-hop decisions, and hop momentum adjustment. See Dynamics Backend Selection.

In .oqp, qmmm(...) is accepted with energy, opt, ground-state md, and embedded namd; writing it also enables qmmm_flag. A QM/MM opt(...) minimizes the QM/MM energy with respect to the QM atoms and, optionally, the MM residues within opt(qmmm_radius=...) Å of the QM region (see examples/QMMM/ala-dipeptide_RHF-QMMM-OPT-linkatom.oqp). Other drivers, such as a standalone grad, are rejected with qmmm(...).

dft/pbe0/def2-svp qmmm()
geom="ala.pdb 9 10 17 18 19"

Available in OpenQP 1.3.0

The QM/MM implementation and its covalent-boundary frontier_scheme controls are included in OpenQP 1.3.0.

Background

The QM subsystem is polarized by the MM electrostatic potential through the electrostatic potential-fitted (ESPF) operator: the MM point charges enter the QM core Hamiltonian, and the reaction of the QM density on the MM atoms is represented by ESPF-fitted atomic charges. This gives an analytic, energy-conserving QM/MM gradient. See References for the ESPF operator and its periodic (particle-mesh Ewald) extension.

Two ways of defining the QM region are supported, matching the driver paths:

• Single-point QM/MM energy reads the geometry and QM selection from [input] system = file.pdb <indices> (see [input] system). Dangling covalent bonds crossing the QM/MM boundary are capped automatically (see Link atoms).
• Ground-state QM/MM molecular dynamics and NAMD-QMMM read the PDB, force field, and QM selection from the [qmmm] keys pdb_file, forcefield_files, and qm_atoms below. NAMD additionally uses runtype=namd and the [md] section.

Give the QM region once. qm_atoms (0-based) and the indices after the PDB name in [input] system (1-based) describe the same atoms, and OpenQP derives whichever one is missing. When both are given they must select the same atoms; otherwise the input is rejected. qm_atoms="0-2" corresponds to system = file.pdb 1 2 3.

Minimal QM/MM Example

Single-point QM/MM energy in .oqp:

dft/bhhlyp/6-31g* qmmm()
geom="ala.pdb 9 10 17 18 19"

Python:

from oqp.openqp import OpenQP

job = OpenQP("ala_qmmm", silent=1)
job.molecule("ala.pdb 9 10 17 18 19", basis="6-31g*")
job.theory("hf", functional="bhhlyp")
job.qmmm(embedding="electrostatic")
mol = job.run()

Legacy .inp (QM selection inline in [input] system):

[input]
qmmm_flag  = true
runtype    = energy
method     = hf
functional = bhhlyp
basis      = 6-31g*
system     = ala.pdb 9 10 17 18 19

[scf]
type = rhf

Embedded nonadiabatic molecular dynamics in .oqp:

mrsf(nstate=5)/bhhlyp/6-31g scf(conv=1e-8) tdhf(conv=1e-8)
namd(S1,scheme=TDC_NAC) md(dt=0.5,nstep=400,velocity="water_box.vel")
qmmm(forcefield_files="amber14-all.xml,amber14/tip3p.xml",qm_atoms="0-2",cutoff=PME)
geom="water_box.pdb"

For TDC_NAC and NAC, the current QM/MM implementation evaluates the analytic vector on the QM centres while holding the instantaneous embedding contribution fixed. Embedding-operator and MM-coordinate derivative components are omitted, and only QM velocities are rescaled. The QM region must contain no link atoms. TDC_NAC still propagates the electronic amplitudes with the phase-tracked overlap TDC, so the environment's change along the actual trajectory remains present in electronic propagation; the QM-region analytic vector is requested only for a selected hop direction.

The legacy .inp spelling is:

[input]
qmmm_flag  = true
runtype    = namd
method     = tdhf
functional = bhhlyp
basis      = 6-31g

[scf]
type = rohf
multiplicity = 3
conv = 1.0e-8

[tdhf]
type   = mrsf
nstate = 5
conv   = 1.0e-8

[md]
active     = 1
tdc        = npi
rescale    = hop_analytic_nac
thrshe     = 0.367493
frustrated = reflect

[qmmm]
pdb_file         = water_box.pdb
forcefield_files = amber14-all.xml,amber14/tip3p.xml
qm_atoms         = 0-2
cutoff           = PME
embedding        = electrostatic

Python API

In the compact OpenQP Python API, job.qmmm(...) enables QM/MM: it sets [input] qmmm_flag=true and the [qmmm] section in one call. forcefield is an alias for forcefield_files; a list is joined into the comma-separated string OpenQP expects, and qm_atoms accepts a string ("0-2") or a list of indices. Any other [qmmm] keyword can be passed through as a keyword argument.

from oqp.openqp import OpenQP

# Single-point QM/MM energy (QM selection inline in job.molecule)
job = OpenQP("ala_qmmm", silent=1)
job.molecule("ala.pdb 9 10 17 18 19", basis="6-31g*", charge=0)
job.theory("hf", functional="bhhlyp")
job.qmmm(embedding="electrostatic")
mol = job.run()

# QM/MM molecular dynamics: PDB, force field, and QM atoms in job.qmmm(...)
job = OpenQP("water_box_qmmm", silent=1)
job.molecule("water_box.pdb", basis="6-31g")
job.theory.mrsf(functional="bhhlyp", nstate=5)
job.qmmm(
    pdb_file="water_box.pdb",
    forcefield=["amber14-all.xml", "amber14/tip3p.xml"],
    qm_atoms="0-2",
    cutoff="PME",
    embedding="electrostatic",
)
job.workflow.namd(
    scheme="custom", tdc="npi", rescale="isotropic",
    thrshe=0.367493, frustrated="reflect",
)  # add soc=True for SOC
job.workflow.md(nstep=200, dt=0.5)
mol = job.run()

See the SOC-NAMD-QMMM workflow and Run from Python for the full nonadiabatic QM/MM setup.

Keywords

pdb_file

Field	Value
Type	string (path)
Default	(empty)
Used by	QM/MM molecular dynamics and SOC-NAMD-QMMM

Path to the PDB file that defines the full QM+MM system (coordinates and topology) for ground-state runtype=md, nonadiabatic runtype=namd and QM/MM geometry optimisation (runtype=optimize, see the optimize page). Single-point QM/MM energy instead takes its PDB path and QM indices together from [input] system = file.pdb <indices>.

forcefield_files

Field	Value
Type	string (comma- or space-separated list)
Default	(empty)
Used by	QM/MM molecular dynamics and SOC-NAMD-QMMM

OpenMM force-field XML files applied to the MM region, e.g. amber14-all.xml,amber14/tip3p.xml for a protein/water system or tip3p.xml for a pure water box. Multiple files are combined by OpenMM in order.

forcefield

Field	Value
Type	string list
Default	amber14-all.xml,amber14/tip3p.xml
Used by	ground-state QM/MM setup

Default force-field list used when a driver builds the MM system without an explicit forcefield_files value. New QM/MM-MD decks should set forcefield_files explicitly.

qm_atoms

Field	Value
Type	string (index list)
Default	(empty)
Used by	QM/MM molecular dynamics and SOC-NAMD-QMMM

Zero-based indices of the atoms placed in the QM region, as individual indices and/or ranges, e.g. 0 1 2 or 0-2 or 0-8 12 15. Either qm_atoms or the 1-based selection after the PDB path in [input] system is sufficient; the other is derived. When both are given they must select the same atoms. The selection is normalised to ascending (topology) order by every driver, so the order you write does not matter. Whole-molecule QM selections (e.g. a solute in a solvent box) are the common case. A selection that cuts a covalent bond is capped with a hydrogen link atom and the MM frontier charge is treated per frontier_scheme; in NAMD this is limited to schemes without analytic NAC (see Link atoms).

active_atoms

Field	Value
Type	string (index list, ranges, or name: groups)
Default	(empty)
Used by	QM/MM geometry optimisation, QM/MM MD and NAMD-QMMM

The atoms the calculation is allowed to move. The spelling follows ORCA's %qmmm ActiveAtoms {0:5 16 21:30} end, so a selection can be copied from either side without editing: indices are zero-based (as qm_atoms is) and a range may be written first:last (ORCA) or first-last. Braces, commas, semicolons and spaces all separate groups, so {0:5 16}, 0-5,16 and 0:5 16 select the same atoms.

A group written name: selects by PDB atom name, matched case-insensitively over the whole system — name:P,OP1,OP2,O5' is a nucleic-acid backbone in one line. The QM region is always active and does not need to be listed.

Leaving this empty keeps each driver's own default: an optimisation moves the QM region alone, dynamics propagates every atom. See Active atoms for what "held" means physically.

frozen_atoms

Field	Value
Type	string (index list, ranges, or name: groups)
Default	(empty)
Used by	QM/MM geometry optimisation, QM/MM MD and NAMD-QMMM

Atoms taken back out of the active set, in the same syntax as active_atoms. This is the one-line way to hold a protein or nucleic-acid backbone that active_radius would otherwise set free:

qmmm(active_radius=4.0,frozen_atoms="name:P,OP1,OP2,O5',C5',C4',C3',O3'")

Used on its own it means "everything else moves", which is usually what a dynamics deck wants. It may not name a QM atom — the QM region is what the calculation moves — and that is reported as an error rather than silently ignored.

active_radius

Field	Value
Type	float (Å)
Default	0.0
Used by	QM/MM geometry optimisation, QM/MM MD and NAMD-QMMM

Adds every MM residue with an atom within this distance of a QM atom to the active set — ORCA's ActiveCore_Extension. 0 adds nothing. Whole residues are selected, so waters and side chains stay intact; in a periodic box the distance is the minimum-image distance. The residue that carries the QM atoms is the exception: a covalent cut runs through it, so its MM atoms are selected one by one.

active_from_pdb

Field	Value
Type	boolean
Default	False
Used by	QM/MM geometry optimisation, QM/MM MD and NAMD-QMMM

Reads the active set from the B-factor column of pdb_file — 1 marks an active atom, 0 an atom to hold — which is ORCA's Use_Active_InfoFromPDB. This keeps the selection with the structure, so a prepared system carries its own active region and no index list has to be maintained alongside it. A file whose B-factor column is zero everywhere is reported as an error rather than run as an empty selection.

A runnable deck is examples/QMMM/H2CO-water_BHHLYP-MRSF-NAMD-QMMM-active.inp.

cutoff

Field	Value
Type	string
Default	NoCutoff
Values	NoCutoff, PME, Ewald, CutoffNonPeriodic, CutoffPeriodic
Used by	QM/MM molecular dynamics and SOC-NAMD-QMMM

OpenMM nonbonded method for the MM region. NoCutoff (and CutoffNonPeriodic) treat the system as an isolated cluster with a direct Coulomb sum, even when the PDB carries a CRYST1 record. PME (particle-mesh Ewald), Ewald and CutoffPeriodic select a periodic (orthorhombic) box: the MM-MM part is evaluated by OpenMM with that method, and the QM/MM electrostatics of the electrostatic embedding are evaluated by Ewald summation — the MM potential and its gradient at every QM centre, the forces on the MM atoms from the ESPF charges and their images, and the interaction of the QM charges with their own periodic images, which is made self-consistent with the ESPF charges outside the SCF — in NAMD with the relaxed charges of the propagated state — (the run stops with an error if that loop does not converge). The periodic branch is available to runtype = md, to the QM/MM runtype = optimize, and to same-spin FSSH (runtype = namd, [md] soc = false); SOC-NAMD requires NoCutoff. See ewald_tol, lj_switch, h_lj and mm_charge_width for the related controls. These controls belong to the drivers built on the ESPF QM/MM force (runtype = md, namd and the QM/MM optimize); a single-point QM/MM energy uses the non-periodic point-charge path and the input checker rejects them there rather than ignoring them.

embedding

Field	Value
Type	string
Default	electrostatic
Values	electrostatic, mechanical
Used by	QM/MM electrostatic coupling

Selects how the MM environment couples to the QM subsystem.

Value	Meaning
electrostatic	Full ESPF electrostatic embedding. The MM charges polarize the QM density through the ESPF operator, and the QM density reacts on the MM atoms via ESPF-fitted charges. This is the production value and gives the analytic, energy-conserving QM/MM gradient.
mechanical	No electrostatic coupling into the QM Hamiltonian: the QM subsystem runs a gas-phase SCF (zero embedding field), and OpenMM evaluates the QM–MM electrostatics with the fixed force-field charges of the QM atoms (intra-QM non-bonded pairs excluded). The MM energy is then differentiated at the charges it was built with, so the force is the exact gradient of the energy and the scheme conserves energy in dynamics; it simply omits the polarisation of the QM density by the environment.

Legacy spellings such as espf and split appear in older decks; new inputs should use electrostatic.

frontier_scheme

Field	Value
Type	string
Default	none
Values	none, rcd, rc, z1
Used by	ESPF electrostatics at a covalent QM/MM boundary

When the QM/MM partition cuts a covalent bond, the MM host atom (M1, the MM end of the severed bond) sits ~1.5 Å from the QM density. frontier_scheme selects how that frontier charge is treated in the ESPF embedding. It is a no-op for whole-molecule QM regions (no cut bond). Covalent QM/MM boundaries are handled by the single-point, optimization, ground-state QM/MM MD, and NAMD paths; NAMD across a covalent boundary requires a scheme without analytic NAC.

Value	Meaning
none	Full-field embedding: the QM density sees the complete MM charge set. This is the default and the validated ESPF baseline. ESPF couples the MM potential to the QM atomic-charge operators (h += Σ_A φ_A Q̂_A), not to the raw density via 1/\|r−R_M\| integrals, which structurally suppresses the electron spill-out that motivates frontier redistribution in density-based embedding — so the ESPF method uses the full MM charges even at a covalent boundary.
rcd	Redistributed charge and dipole: delete M1's charge and place virtual point charges at the M1–M2 bond midpoints (2·q₁/N) plus −q₁/N on each MM neighbour M2, conserving the total charge and the dipole about M1. The virtual charges sit at bond midpoints (linear in the real atom positions), so the analytic gradient stays exact.
rc	Redistributed charge: midpoint charges only (conserves the total charge, not the dipole).
z1	Charge deletion: remove M1's charge (conserves neither; provided for comparison).

rcd/rc/z1 are optional refinements, not the ESPF default. See References for the ESPF charge-operator formulation and the redistribution schemes.

rigidwater

Field	Value
Type	boolean
Default	False
Used by	MM constraint setup

Apply rigid-water (SHAKE/RATTLE) constraints to MM water molecules. QM atoms are never constrained. Rigid water removes the stiff O-H stretch from the MM region and allows a normal MD timestep (~0.5-1 fs). The nonadiabatic velocity-Verlet loop always constrains MM rigid water; QM atoms move under the QM forces.

ewald_tol

Field	Value
Type	float
Default	unset (OpenMM default, 5e-4)
Used by	periodic QM/MM MD, NAMD-QMMM and QM/MM optimisation (cutoff = PME / Ewald)

Error tolerance of the OpenMM PME/Ewald evaluation of the MM-MM electrostatics. OpenMM's default gives a force error of about 0.5 kJ/mol/nm, which is the floor of any finite-difference force check or NVE energy-conservation test on a periodic box; use 1e-6 for such validation work. Ignored for non-periodic cutoffs. Must be a positive number.

lj_switch

Field	Value
Type	boolean
Default	False
Used by	QM/MM MD, NAMD-QMMM and QM/MM optimisation with a finite cutoff

Switch the MM Lennard-Jones interactions smoothly to zero over the last 15 % of the nonbonded cutoff (OpenMM setUseSwitchingFunction). A plain truncation makes the MM energy discontinuous when pairs cross the cutoff, which shows up as a drift of the total energy in NVE dynamics (about +10 kJ/mol/ps for a solvated dipeptide with the OpenMM defaults). No effect under NoCutoff.

h_lj

Field	Value
Type	boolean
Default	False
Used by	QM/MM MD, NAMD-QMMM and QM/MM optimisation (electrostatic embedding)

Give Lennard-Jones parameters to MM hydrogens that have none in the force field (TIP3P water hydrogens), using the CHARMM TIP3P values (R_min/2 = 0.2245 Å, ε = 0.046 kcal/mol). Without them nothing but the QM density keeps a water hydrogen from approaching a QM oxygen too closely, since a point charge has no Pauli wall. QM atoms are never modified. On its own this is a weak barrier; see mm_charge_width.

mm_charge_width

Field	Value
Type	float (Å)
Default	unset (point charges)
Used by	QM/MM MD, NAMD-QMMM and QM/MM optimisation (electrostatic embedding)

Represent the MM charges as Gaussians of this width in the QM/MM electrostatics (full-ESPF embedding = electrostatic only; other embeddings reject it). The QM-MM pair potential 1/r becomes erf(μr)/r with μ = 1/(√2 w), consistently in the energy, the QM and MM forces, the direct sum and the real-space part of the Ewald sum; the MM-MM interactions and the QM-image term are untouched. This is the erf damping of the reference Tinker ESPF implementation (its ERFMU keyword); 0.7 Å corresponds to its default μ = 1/Å. It removes the 1/r singularity a bare point charge presents to the QM density, which otherwise lets a solvent hydrogen collapse onto a QM carbonyl oxygen through a polarisation runaway of the ESPF charges, after which the MRSF Z-vector equations diverge. Must be a positive number; 0 or unset means point charges. In a periodic box the damping correction is summed over the minimum image inside the real-space cutoff (half the shortest box edge) only, so the width must be short-ranged on that scale: w ≤ L_min / (9 √2) (1.26 Å for a 16 Å box); larger values are rejected with a message giving the bound.

nonbondedmethod

Field	Value
Type	string
Default	NoCutoff
Used by	ground-state QM/MM setup

OpenMM nonbonded method for the ground-state QM/MM path. The nonadiabatic and newer MD paths use cutoff instead.

constraints

Field	Value
Type	string
Default	None
Used by	ground-state QM/MM MD

OpenMM bond/angle constraint setting for the MM region in the ground-state MD path (for example, constrain X-H bonds).

temperature

Field	Value
Type	float (K)
Default	300.0
Used by	ground-state QM/MM MD

Legacy target/initial temperature for ground-state QM/MM MD. New concise input uses md(temperature=...) for both ground-state and nonadiabatic dynamics.

timestep

Field	Value
Type	number (fs)
Default	1
Used by	ground-state QM/MM MD

Legacy QM/MM MD timestep. New concise input uses md(dt=...) for all dynamics.

nsteps

Field	Value
Type	integer
Default	1
Used by	legacy static QM/MM bookkeeping

Retained for compatibility with the older static-driver configuration. The current ground-state OpenMM MD engine also reads n_steps. New concise input uses md(nstep=...) for all dynamics; qmmm(nsteps=N) remains a compatibility alias.

istate

Field	Value
Type	integer
Default	0
Used by	ground-state QM/MM

Obsolete numeric selector from the disconnected legacy libopenmm path. It is retained in the sectioned schema for compatibility but is reserved in .oqp; do not write qmmm(istate=...). Ground-state canonical md uses S0, and state-aware drivers own their physical state labels.

Legacy Ground-state OpenMM MD Keywords

These [qmmm] propagation keys remain available to old sectioned inputs. New concise inputs put propagation in md(...) and keep qmmm(...) limited to the QM/MM structure, force field, embedding, boundary, and periodic settings. The ground-state QM/MM command must run without MPI.

n_steps

Field	Value
Type	integer
Default	1000
Used by	number of ground-state QM/MM-MD integration steps

This is the preferred MD spelling. Concise .oqp accepts nsteps as an alias for this key, while a traditional sectioned .inp keeps nsteps as separate legacy bookkeeping with its own default of one.

ensemble

Field	Value
Type	string
Default	nve
Supported values	nve, nvt

This is a legacy [qmmm] propagation key. New .oqp and Python inputs use md(ensemble=nve|nvt) instead. The current unified md(...) surface supports NVE and Langevin NVT. NPT is a defined ensemble and the schema retains pressure and barostat_interval, but the current QM/MM runtime rejects it: the OpenMM barostat trial would otherwise use a linearized external force instead of a newly evaluated QM/MM energy for the proposed box.

friction

Field	Value
Type	float, ps^-1
Default	1.0
Used by	Langevin integration for NVT/NPT

pressure

Field	Value
Type	float, bar
Default	1.0
Used by	NPT barostat target

barostat_interval

Field	Value
Type	integer steps
Default	25
Used by	NPT barostat update interval

NPT requires a periodic nonbonded method rather than NoCutoff.

trajectory_format

Field	Value
Type	string
Default	pdb
Values	pdb, dcd

trajectory_file

Field	Value
Type	path
Schema default	empty
Used by	coordinate trajectory output

When omitted from the active MD configuration, the driver uses qmmm_trajectory.<trajectory_format>.

log_file

Field	Value
Type	path
Schema default	empty
Used by	OpenMM state-data output

The active MD driver uses qmmm_trajectory.dat when no path is supplied.

report_interval

Field	Value
Type	integer steps
Default	1
Used by	trajectory, log, and energy reporting frequency

energy_file

Field	Value
Type	path
Schema default	empty
Used by	saved QM, MM, kinetic, and total-energy arrays

The active MD driver uses total_energy.npz when no path is supplied.

qm_atoms_xyz

Field	Value
Type	path
Default	empty
Used by	optional replacement coordinates for the QM atoms

qm_list

Field	Value
Type	integer list
Default	empty
Used by	mapping QM atoms to rows of qm_atoms_xyz

When qm_atoms_xyz is present and qm_list is omitted, rows are used in QM atom order. If supplied, qm_list must have the same length as the QM atom selection and contain valid XYZ indices.

Active atoms (what moves)

A solvated or biomolecular QM/MM system is mostly environment, and propagating all of it is what makes a long trajectory or a large optimisation expensive. The [qmmm] keys active_atoms, frozen_atoms, active_radius and active_from_pdb select which atoms actually move. They follow ORCA's %qmmm ActiveAtoms naming and syntax, and they are read by the QM/MM optimiser, ground-state QM/MM MD and the NAMD-QMMM drivers alike.

The active set is assembled as

QM region  +  active_radius shell  +  active_atoms  −  frozen_atoms

and the QM region is always in it.

Held is not absent

Freezing an atom changes what moves, never the physics. A held atom keeps its MM charge, its contribution to the ESPF embedding field that polarizes the QM density, and the force it exerts on the atoms that do move. It is simply not integrated — so no work is done on it, and energy conservation still applies to the atoms that move. What changes is the cost and the reported number of degrees of freedom.

Each driver holds an atom in the way that is natural for it:

Path	How a held atom is held
runtype=optimize	its coordinates never enter the optimisation vector
runtype=namd (FSSH, SOC, SOC-MCH)	its rows are skipped in the velocity-Verlet update, so it keeps zero velocity; the Maxwell draw, the centre-of-mass correction and the degree-of-freedom count all run over the atoms that move
runtype=md	OpenMM gives the particle zero mass, which is how OpenMM holds an atom in place

Rigid-water constraints are resolved against the selection, because SHAKE/RATTLE cannot satisfy a constraint between a moving atom and a fixed one — it would drag the fixed atom. A constrained pair that straddles the boundary is made whole: held if either atom was frozen, moving otherwise. Constraints inside the held set are then dropped, since those atoms never move.

Defaults preserve every existing deck: with none of these keys set, an optimisation moves the QM region alone and dynamics propagates every atom. A NAMD restart is bound to the selection so a restart cannot silently change what moves, but only once a selection is actually requested — checkpoints written before these keys existed keep validating.

The optimiser also accepts the [optimize] spellings qmmm_radius — the released name of the movable shell — together with qmmm_active and qmmm_freeze, as aliases of active_radius, active_atoms and frozen_atoms; [qmmm] wins if both are given, and only the [qmmm] spelling reaches the dynamics drivers.

Covalent QM/MM boundaries

Most QM/MM calculations put a whole molecule (or a set of whole molecules) in the QM region — a chromophore in a solvent box, a substrate in a binding pocket. When the QM/MM partition instead cuts through a covalent bond (e.g. carving a side chain out of a protein backbone), two problems appear at the frontier, and OpenQP handles each automatically.

1. The dangling bond → hydrogen link atom

Cutting a bond leaves the QM frontier atom with an unsatisfied valence. OpenQP caps it with a hydrogen link atom on the broken QM–MM bond at the scaled (IMOMM) position R_L = R_QM + g·(R_MM − R_QM), with g = (r_H + r_QM)/(r_QM + r_MM) from covalent radii. Because the cap position is a fixed linear function of its two real host atoms, the force computed on the link atom is redistributed onto them by the chain rule — no extra degrees of freedom, and the analytic gradient stays exact. Details in Link atoms.

2. The frontier MM charge → ESPF is robust; frontier_scheme refines it

The MM host atom (M1, the MM end of the severed bond) sits only ~1.5 Å from the QM density. In a density-based electrostatic embedding, embedding M1's full point charge that close over-polarizes the QM density (electron spill-out), which is why such methods delete or shift the frontier charge.

OpenQP's embedding is ESPF, which is structurally more robust here: the MM potential couples to the QM atomic-charge operators, h += Σ_A φ_A Q̂_A (Huix-Rotllant & Ferré; see Background), not to the raw density via 1/|r−R_M| integrals. Because the coupling runs through bounded atomic-charge operators evaluated at the QM nuclei, the spill-out that plagues density-based embedding is suppressed — so the full MM charges can be used "with no scaling" even at a covalent boundary. That is the default, frontier_scheme = none.

For a still-more-faithful boundary field, the optional redistributions delete the raw M1 charge and spread it over virtual point charges at the M1–M2 bond midpoints (M2 = the MM neighbours of M1):

• rcd — 2·q₁/N at each of the N midpoints plus −q₁/N on each M2, conserving both the total charge and the dipole about M1.
• rc — q₁/N at each midpoint (total charge only).
• z1 — delete M1 (conserves neither; for comparison).

Each virtual charge sits at a bond midpoint — a linear function of the real M1/M2 positions — so the electrostatic force on it redistributes onto its hosts by the chain rule, exactly like the link-atom cap, keeping the analytic gradient exact.

Gradient accuracy

The full analytic QM/MM gradient (link-atom projection + ESPF density response + frontier redistribution) matches a finite-difference of the self-consistent QM/MM energy to ~1e-3 au — the ESPF grid-derivative floor, a property of the embedding independent of the boundary treatment. The redistribution does not degrade this; at the deleted host M1 the redistributed schemes are in fact more accurate than full-field, because they remove the raw close-range point charge.

Which path supports it

Covalent-boundary QM/MM is available in the single-point, optimization, ground-state QM/MM MD (runtype=md), and NAMD (runtype=namd) paths. NAMD across a covalent boundary is limited to schemes that do not use analytic NAC (BaeckAn, Overlap, or custom with tdc=npi and isotropic rescaling); TDC_NAC and NAC are rejected when the QM region has link atoms, because the link-centre NAC direction is not yet projected onto its QM and MM hosts, and a default rescale=auto falls back to isotropic rescaling. A NAMD deck is examples/QMMM/ala-dipeptide_BHHLYP-MRSF-NAMD-QMMM-linkatom.oqp. A runnable ground-state deck is examples/QMMM/ala-dipeptide_BHHLYP-QMMM-MD-RCD.inp (alanine dipeptide, QM = the C-terminal amide cutting the ALA C–CA bond, frontier_scheme=rcd); see QM/MM examples.

Link atoms

When the QM region defined through [input] system = file.pdb <indices> cuts a covalent bond, OpenQP caps each dangling bond automatically with a hydrogen link atom — there is no keyword to enable it; the boundary bonds are detected from the PDB topology. The link hydrogen is placed along the broken QM-MM bond at a scaled distance set by the covalent-radius factor

g = (r_H + r_QM) / (r_QM + r_MM)

where r_H, r_QM, and r_MM are the covalent radii of hydrogen and of the QM and MM boundary atoms. Only hydrogen capping is currently supported. The link-atom energy gradient is redistributed onto its two real host atoms by the chain rule of this scaled position, so no extra degrees of freedom are added.

Across a covalent boundary the MM frontier-host charge is treated per frontier_scheme (default none = full-field ESPF), and the ESPF grid switching width is selected automatically — see ESPF grid switching below.

Covalent boundaries in nonadiabatic dynamics

Automatic link-atom capping also applies to QM/MM NAMD, but only with schemes that do not use analytic NAC. TDC_NAC and NAC require a whole-molecule QM region.

The ESPF gradient accounts for the fitting grid moving with the atoms (the translational counter-terms of every grid point are charged to its parent atom), so the force is the exact derivative of the embedded energy and QM/MM molecular dynamics conserves the total energy in vacuum and in a periodic box, with or without a link atom.

ESPF grid switching

The ESPF embedding gradient carries a small finite-difference residual (the grid-derivative floor, ~1e-3 a.u.). At a covalent boundary that residual does not spread over the QM region: it concentrates as an over-shoot on the MM host atom, and it is governed by the smooth-switching width ESPF_SWSCALE.

The shipped whole-molecule value, 1.8, over-smooths at a covalent boundary, and the tighter 1.5 is about 20% worse without one — neither is a safe global default. OpenQP therefore selects it per system, from whether the QM/MM partition actually cut a bond:

QM/MM boundary	ESPF_SWSCALE
covalent (a bond is cut; link atoms present)	1.5, selected automatically
whole-molecule (no bond cut)	1.8

The choice is made by the QM/MM driver when it builds the partition, because that is where it is known whether a bond was cut: the native ESPF gradient is handed a grid and a density and cannot tell. A covalent-boundary run reports the selection:

[QM/MM] 1 covalent boundary bond(s); frontier-charge embedding = ...
[QM/MM] covalent boundary detected; ESPF_SWSCALE=1.5 (whole-molecule default
        is 1.8). Set ESPF_SWSCALE to override.

Setting ESPF_SWSCALE yourself always wins and the automatic choice stands down entirely, so a parameter sweep or a reproduction is never silently overridden:

export ESPF_SWSCALE=1.8      # force the whole-molecule width at a boundary

ESPF environment variables

These are runtime environment variables, not [qmmm] keywords, read by the native ESPF module:

Variable	Default	Meaning
ESPF_SMOOTH	on	enable smooth grid switching at all
ESPF_SWSCALE	selected (1.5 / 1.8)	switching width scale, see above
ESPF_SWDELTA	0.7	switching offset
ESPF_WSCALE	1.0	grid weight scale; keep at 1.0
ESPF_WDERIV	on	include the grid weight-derivative gradient term

Applies to every electrostatically embedded QM/MM path

Grid switching is a property of baseline ESPF, so it affects every path that runs embedding = electrostatic, including SOC-NAMD-QMMM, and it is independent of frontier_scheme. Under embedding = mechanical there is no electrostatic coupling into the QM Hamiltonian and therefore no ESPF grid at all, so none of these variables — including the automatic ESPF_SWSCALE selection — has any effect.

Notes

• Set [input] qmmm_flag=true to activate any QM/MM path; without it the [qmmm] section is ignored.
• For nonadiabatic QM/MM dynamics, combine this section with the [md] section and see the SOC-NAMD-QMMM workflow.
• Use cutoff = PME (or Ewald) with a solvated periodic water box for production QM/MM-MD; NoCutoff is for isolated clusters. For energy conservation in solution set ewald_tol = 1e-6 and lj_switch = true; for excited-state (MRSF) dynamics in water add mm_charge_width = 0.7 (and optionally h_lj = true).