Constrained QM/MM with FlexiBLE

QM/MM is normally applied to systems where the QM and MM regions hold chemically distinct species — e.g. a QM solute in an MM solvent, or a QM enzymatic active site embedded in an MM protein. A harder problem arises when QM and MM contain chemically identical particles, such as a QM solute plus its first solvation shell surrounded by an MM solvent bath. In ordinary AIMD nothing prevents the molecules from diffusing across the QM/MM boundary, so after a few thousand steps the partition is destroyed.

TeraChem ships an implementation of the Flexible Boundary Layer using Exchange (FlexiBLE) method to keep the partition intact during dynamics.¹ Unlike adaptive QM/MM schemes that let the count of QM particles fluctuate while fixing the volume, FlexiBLE fixes the count of QM particles and lets the volume of the QM region fluctuate. FlexiBLE rigorously preserves canonical ensemble averages and reproduces dynamical properties well on sub-diffusional timescales for particles in the inner QM region.

FlexiBLE is layered on the OpenMM/Amber driver — you need a prmtop + rst7 + qmindices setup as for any Amber-driven TeraChem QM/MM job — and is enabled by setting qmmmbp flexible.

Bias potential

FlexiBLE adds a classical bias to the system Hamiltonian:

\[V^{\mathrm{bias}}_{ijk\cdots} = -k_B T \log f_{ijk\cdots},\]

where the bias function \(f_{ijk\cdots}\) is normalized over all permutations of identical particles between QM and MM:

\[f_{ijk\cdots} = \frac{h_{ijk\cdots}}{\sum_L \hat{P}_L(h_{ijk\cdots})},\]

with \(\hat{P}_L\) permuting the QM/MM identity assignments. The unnormalized \(h_{ijk\cdots}\) is a product of pair functions \(g_{ik}\) between each QM and MM particle:

\[g_{ik} = \begin{cases} 1, & x_i < x_k \\ \exp\!\left(-\dfrac{\alpha^3 (x_i - x_k)^3}{1 + \alpha(x_i - x_k)}\right), & x_i \geq x_k, \end{cases}\]

where \(x_i\) and \(x_k\) are radial distances from the system origin and \(\alpha\) controls how strongly the bias penalizes QM/MM mixing. Larger \(\alpha\) gives sharper separation. The functional form ensures that the bias penalizes any configuration in which a QM particle is farther from the origin than an MM particle.

The sum over permutations would naively involve \((N_{\mathrm{QM}} + N_{\mathrm{MM}})! / (N_{\mathrm{QM}}! N_{\mathrm{MM}}!)\) terms, but with a tree algorithm and a screening threshold, only a handful of permutations contribute and the cost stays manageable.

Configuring FlexiBLE

The required keyword to turn FlexiBLE on is:

qmmmbp flexible

The behavior is then controlled by a handful of flexible_* keywords:

Keyword	Purpose
`flexible_alpha`	Exponent \(\alpha\) in Å\(^{-1}\). Default and recommended value: 15. Larger values improve QM/MM separation but require a smaller timestep.
`flexible_thresh`	Screening threshold for discarding negligible permutations. Default 0.1. Smaller values include more permutations at increased cost.
`flexible_scale`	Per-iteration multiplicative tightening of the threshold. Default 0.5.
`flexible_maxit`	Maximum number of tree-sweep iterations. The job aborts if convergence is not reached. Default 10.
`flexible_target`	Target-atom assignment (per molecule group; use `-1` for solvent-only).
`flexible_type`	Boundary geometry: `spherical` (default) or `capsule`
`flexible_usecom`	`yes` (default) places the origin at the total system center of mass; `no` uses an explicit origin set by `flexible_center`
`flexible_center`	Origin coordinates (Å), used when `flexible_usecom no`
`molecule_info`	File describing how atoms are grouped into molecules
`capsule_length`	Length of the capsule (Å), used when `flexible_type capsule`
`capsule_vector`	Direction of the capsule's long axis (only direction matters; magnitude is renormalized)

When flexible_usecom yes (the default), TeraChem also propagates gradient contributions from the center of mass so that total classical energy and momentum are rigorously conserved. When combining FlexiBLE with a spherical droplet boundary condition (mdbc spherical), set the same origin convention for both — use flexible_usecom no together with mdbc_no_centering yes so the FlexiBLE bias and the droplet potential share an origin.

Boundary geometries

Spherical (default)

In the default spherical mode, QM particles are biased to stay inside a sphere around the origin and MM particles to stay outside. Critically, the radius is not fixed — the number of QM particles is what is conserved, while the sphere's effective volume fluctuates with the dynamics.

flexible_sphere.in

prmtop          spcfwat.prmtop
coordinates     geomvels.rst7
velocities      geomvels.rst7
qmindices       md_qmindex.txt
molecule_info   moleculeinfo.txt

basis           mini_[scaled]
method          rhf
charge          0
spinmult        1

run             md
nstep           40
timestep        0.5
thermostat      nhc
t0              300

mdbc                spherical
mdbc_no_centering   yes
md_r1               25
md_k1               10

qmmmbp              flexible
flexible_type       spherical
flexible_alpha      15
flexible_thresh     0.1
flexible_scale      0.5
flexible_maxit      10
flexible_target     0
flexible_usecom     no
flexible_center     0.0 0.0 0.0

precision           double
threall             1.0e-16
convthre            1.0e-8
end

Capsule

For QM regions that are extended in one direction (e.g. a conjugated oligomer or peptide chain), a sphere wastes solvent on the short axes. FlexiBLE supports a capsule geometry — a cylinder capped with hemispheres at each end — whose long axis is set by capsule_vector and whose length is set by capsule_length:

FlexiBLE QM/MM capsule boundary geometry

flexible_capsule.in

prmtop          ch3oh_water.prmtop
coordinates     init.rst7
velocities      init.rst7
qmindices       qmindices.txt
molecule_info   moleculeinfo.txt

basis           sto-3g
method          rhf
charge          0
spinmult        1

run             md
nstep           30

mdbc                spherical
md_r1               25
md_k1               10

qmmmbp              flexible
flexible_type       capsule
flexible_usecom     no
flexible_center     0.0 0.0 0.0
capsule_vector      0.0 1.0 0.0
capsule_length      4.0

# Two QM groups, separate parameters per group:
flexible_alpha      15 15
flexible_thresh     0.1 0.1
flexible_scale      0.5 0.5
flexible_maxit      10 10
flexible_target     -1 0

precision           double
threall             1.0e-16
convthre            1.0e-8
end

When the capsule center is fixed at the origin (flexible_usecom no plus flexible_center 0 0 0), it is the user's responsibility to ensure that the initial coordinates place the long axis of the QM region along the direction set by capsule_vector. FlexiBLE will not rotate the system to fit the capsule.

Solvent restrictions

The current FlexiBLE implementation supports homogeneous solvents whose molecules have a single heavy atom — water, ammonia, methane, and the like. Mixed solvents and electrolytes are not yet supported. In addition to solvent molecules, the bias acts on the solute by treating it as identical to the solvent, which can introduce artifacts if the solute reaches the QM/MM boundary. The recommended workaround is to fix one atom of the solute at the system origin and use a solvent shell large enough to contain at least one complete coordination sphere.

If the initial coordinates have too much QM/MM mixing, FlexiBLE will fail with an error message rather than producing biased dynamics — adjust the qmindices selection or pre-equilibrate the system to clean up the initial partition.

J. Chem. Phys. 155, 224112 (2021). ↩