QM/MM with the Built-in OpenMM/Amber Driver

For QM/MM beyond the simple water-only case, TeraChem links OpenMM and uses an Amber-formatted topology + coordinate set to define the MM region. This is the most flexible TeraChem-driven QM/MM mode: any system that Amber's tleap can parameterize is fair game.

Activating the driver

The combination of keywords that activates the OpenMM/Amber driver is:

prmtop      waters.prmtop
coordinates waters.rst7
qmindices   qmregion.dat

prmtop is a standard Amber parameter file (typically generated by tleap), and coordinates is the matching Amber inpcrd / rst7 coordinate file. Together they fully describe the MM region's connectivity, charges, and geometry. The qmindices file then carves out which atoms TeraChem treats quantum-mechanically.

Optionally:

velocities  waters.rst7

points TeraChem at a file containing per-atom velocities (usually the same rst7 as the coordinates, when Amber has written one). Without this keyword TeraChem generates random velocities at tinit. Velocities are only used by run md; other run modes ignore them. A specific random seed for velocity generation can be requested with the seed keyword.

The `qmindices` file

The qmindices file lists one QM atom per line. The first column is the 0-based index of the atom in the prmtop file, and an optional second column overrides the element label that TeraChem would otherwise infer:

2water.qmregion

0
1
2

Here atoms 0, 1, and 2 of the prmtop (the first water molecule) form the QM region. Custom atom-type labels are useful when you want to use a non-standard basis set for a specific atom:

custom_types.qmregion

0 Oa
1 Ha
2 Ha
3 O
4 H
5 H

In this example two waters are in the QM region. The first water is given custom labels Oa/Ha so a custom basis set can be assigned to those atoms; the second water uses the standard element labels.

Index base inconsistency with atomic constraints

QM atom indices in the qmindices file start at 0. But if you later add atomic constraints (e.g. for the optimizer), those use 1-based atom numbering. This off-by-one mismatch is a well-known foot-gun.

Restrictions on the prmtop

The OpenMM/Amber driver imposes two constraints on the parameter file:

No periodic boundary conditions. TeraChem does not currently support PBC in the prmtop. To run a "bulk-like" simulation, build a non-periodic spherical droplet and apply TeraChem's spherical BCs to confine it.
Element types must be recognizable to TeraChem, OR the per-atom type must be overridden in the second column of the qmindices file.

Optimization and NEB

For optimizations on QM/MM systems set up this way, Cartesian coordinates are recommended unless a substantial fraction of the system is frozen:

run             minimize
min_coordinates cartesian

Internal coordinates (DLC/HDLC/TRIC) often struggle with large solvent shells.

For nudged-elastic-band runs, the coordinates file is the first endpoint (neb1.rst7 by convention) and a second file neb2.rst7 in the same directory is the other endpoint (alternatively, nebcrd <file> overrides the second-endpoint name).

Examples

Single-point gradient on a QM-water/MM-water pair

waters_openmm_gradient.in

prmtop          2water.prmtop
qmindices       2water.qmregion
coordinates     2water.rst7

basis           3-21g
method          hf
charge          0
spinmult        1
threall         1e-14
convthre        1e-7

run             gradient
end

Short MD trajectory using Amber-supplied initial velocities

waters_openmm_md.in

prmtop          2water.prmtop
qmindices       2water.qmregion
coordinates     2water.rst7
velocities      2water.rst7

basis           3-21g
method          hf
charge          0
spinmult        1

run             md
nstep           10
timestep        1.0
seed            100

end

CASSCF QM/MM on propene

CAS-type methods work seamlessly with the OpenMM driver:

propene_qmmm_casscf.in

prmtop                 propene.parm7
coordinates            propene.rst7
qmindices              propene.qmregion

basis                  6-31g
method                 HF
run                    gradient
charge                 0
spinmult               1

casscf                 yes
closed                 7
active                 2
cassinglets            2
castarget              1
castargetmult          1
alphacas               yes
alpha                  0.64

end

NEB between two Amber endpoints

propene_neb.in

prmtop          propene.parm7
qmindices       propene.qmregion
coordinates     neb1.rst             # other endpoint is neb2.rst in the same directory

basis           6-31g
method          b3lyp
precision       mixed
charge          0
spinmult        1
threall         1e-14
convthre        1e-10

run             ts
ts_method       neb_free
min_method      lbfgs
min_image       8
min_tolerance   1e-4
min_tolerance_e 1e-4
nstep           400

end

Summary of keywords

Keyword	Type	Required?	Description
`prmtop`	filename	yes (with Amber coords)	Amber parameter file (parm7/prmtop) describing the full QM+MM system
`coordinates`	filename	yes	Amber-format `inpcrd` or `rst7` coordinate file
`qmindices`	filename	yes (with prmtop)	File listing 0-based indices of atoms in the QM region, one per line. Optional second column overrides element label.
`velocities`	filename	no	File containing initial atomic velocities (only used by `run md`)
`seed`	int	no	Random seed for initial-velocity generation
`nebcrd`	filename	no	Path to the second NEB endpoint (defaults to `neb2.rst7`)