QM/MM in AMBER with QUICK¶

This tutorial describes how to run minimization, equilibration, and production with AMBER files generated by CHARMM-GUI, including QM/MM simulations with QUICK and semiempirical methods.

Assumed file names:

• step3_input.parm7
• step3_input.rst7
• step4.0_minimization.mdin
• step4.1_equilibration.mdin
• step5_production.mdin

If your files use different names, substitute them in the commands below.

1. Minimization¶

Run minimization with sander on CPU. If dihe.restraint exists, update the restraint force placeholder first:

grep -q "FC" dihe.restraint && sed -e "s/FC/1.0/g" dihe.restraint > step4.0_minimization.rest

Run minimization:

sander -O -i step4.0_minimization.mdin -p step3_input.parm7 -c step3_input.rst7 -o step4.0_minimization.mdout -r step4.0_minimization.rst7 -inf step4.0_minimization.mdinfo -ref step3_input.rst7

For GPU execution, use pmemd.cuda with the same arguments when the calculation type supports it:

pmemd.cuda -O ...

2. Equilibration¶

If dihe.restraint exists:

sed -e "s/FC/1.0/g" dihe.restraint > step4.1_equilibration.rest

Run equilibration:

sander -O -i step4.1_equilibration.mdin -p step3_input.parm7 -c step4.0_minimization.rst7 -o step4.1_equilibration.mdout -r step4.1_equilibration.rst7 -inf step4.1_equilibration.mdinfo -ref step3_input.rst7 -x step4.1_equilibration.nc

GPU form:

pmemd.cuda -O ...

3. Production¶

sander -O -i step5_production.mdin -p step3_input.parm7 -c step4.1_equilibration.rst7 -o step5.mdout -r step5.rst7 -inf step5.mdinfo -x step5.nc

For MPI-based execution:

mpirun -np 2 sander.MPI -O -i step5_production.mdin -p step3_input.parm7 -c step4.1_equilibration.rst7 -o step5.mdout -r step5.rst7 -inf step5.mdinfo -x step5.nc

Semiempirical QM/MM with AM1 is not accelerated through pmemd.cuda. To run QM/MM where the quantum region is evaluated on GPU, use QUICK through AmberTools.

4. QM/MM Dynamics with QUICK¶

SANDER can use different QUICK installations for QM/MM:

• serial
• MPI parallel
• CUDA accelerated
• HIP accelerated
• CUDA + MPI
• HIP + MPI

QUICK can be used through two main interfaces:

• FBI, the file-based interface
• API, the in-memory application programming interface

5. File-Based Interface¶

With the file-based interface, sander can call QUICK executables such as:

• quick
• quick.MPI
• quick.cuda
• quick.hip
• quick.cuda.MPI
• quick.hip.MPI

Before running QUICK simulations, load the AMBER environment:

source $AMBERHOME/amber.sh

or, for C shell:

source $AMBERHOME/amber.csh

FBI Example with Mechanical Embedding¶

The example below places the first two residues in the QM region and runs B3LYP/def2-SVP with quick.cuda.MPI on two GPUs:

&cntrl
  ...
  ifqnt = 1,
/
&qmmm
  qmmask = ':1-2',
  qm_theory = 'extern',
  qmmm_int = 5,
/
&quick
  method = 'B3LYP',
  basis = 'def2-svp',
  executable = 'quick.cuda.MPI',
  do_parallel = 'mpirun -np 2',
/

Main flags:

• ifqnt = 1: enables QM/MM
• qmmask: selects atoms in the QM region
• qm_theory = 'extern': uses an external program through FBI
• qmmm_int = 5: mechanical embedding
• executable: selects the QUICK executable
• do_parallel: command prefix for MPI QUICK executables

In many FBI workflows, use serial sander rather than sander.MPI, because system calls inside MPI programs can be limited. This means the MM region may run on one CPU core even when the QM calculation uses multiple GPUs.

6. QUICK API¶

The API communicates directly in memory and is usually much faster than FBI.

Example for electrostatic embedding:

&cntrl
  ...
  ifqnt = 1,
/
&qmmm
  qmmask = ':1-2',
  qm_theory = 'quick',
  qmmm_int = 1,
  qm_ewald = 0,
/
&quick
  method = 'B3LYP',
  basis = 'def2-svp',
/

Important points:

• qm_theory = 'quick' activates the QUICK API.
• qmmm_int = 1 uses electrostatic embedding.
• The same input can be used with sander, sander.MPI, sander.quick.cuda, or sander.quick.cuda.MPI.
• QUICK output is written to quick.out unless changed.

7. &quick Namelist Summary¶

Variable	Type	Interface	Description
method	String	API/FBI	Method: HF, B3LYP, wB97X-D, PBE, etc.; default is BLYP.
basis	String	API/FBI	Basis set: 6-31G*, def2-SVP, cc-pVDZ, etc.; default is 6-31G.
executable	String	FBI only	QUICK executable, such as quick.cuda.MPI or quick.hip.
do_parallel	String	FBI only	Prefix such as mpirun -np 4 or srun --gpus=4.
scf_cyc	Integer	API/FBI	Maximum SCF cycles; default is 200.
reuse_dmx	Integer	API only	Reuse previous density matrix; 1 is default and accelerates MD.
denserms	Float	API only	Density-matrix convergence criterion.
intcutoff	Float	API only	Integral cutoff.
xccutoff	Float	API only	Exchange-correlation pruning threshold.
basiscutoff	Float	API only	Cutoff for insignificant basis functions.
gradcutoff	Float	API only	Gradient convergence criterion.
export	String	API only	Export orbitals, for example molden.
keywords	String	API only	Full QUICK keyword line, for example B3LYP BASIS=cc-pVDZ CHARGE=0 MULT=1 GRADIENT EXTCHARGES.
outfprefix	String	API only	Prefix for QUICK output; default creates quick.out.
debug	Integer	API/FBI	Debug level.
use_template	Integer	FBI only	Use a QUICK template file.

Practical rule:

• For maximum speed and electrostatic embedding, use the API with qm_theory='quick' and sander.quick.cuda or sander.quick.cuda.MPI.
• Use FBI only when a specific QUICK option is not exposed through the API.

8. FBI vs API¶

Feature	FBI, qm_theory='extern'	API, qm_theory='quick'
Communication	Through files	Directly in memory
Speed in MD	Usually much slower	Fastest option
MM part	Usually serial	Can use MPI/OpenMP depending on executable
Multi-GPU QM	Supported with quick.cuda.MPI	Supported with sander.quick.cuda.MPI
SANDER executable	Usually serial sander	sander.quick.cuda or sander.quick.cuda.MPI
Electrostatic embedding	Supported	Supported and recommended
Mechanical embedding	Supported	Supported
QUICK output	Many separate files	Usually one quick.out
Typical use	Rare compatibility cases	Default choice for performance

9. Example QUICK Commands¶

Single GPU:

sander.quick.cuda -O -i md.in -o md.out -p topo.prmtop -c coord.rst

Multi-GPU with parallel MM:

mpirun -np 8 sander.quick.cuda.MPI -O -i md.in -o md.out -p topo.prmtop -c coord.rst

Example .mdin fragment:

&qmmm
  qm_theory = 'quick',
  qmmask = ':LIG | :120,130,145',
/
&quick
  method = 'PM6',
  basis = '6-31G*',
/

10. Study-Specific QUICK Production Flags¶

Example production input with explicit iqmatoms:

&qmmm
  iqmatoms=593, 594, 595, 596, 597, 598, 3722, 3723, 3724, 3725, 3726, 3727, 3728, 3729, 3730, 3731, 45573, 45574, 45575, 45624, 45625, 45626, 56283, 56284, 56285, 58281, 58282, 58283, 58680, 58681, 58682, 58608, 58609, 58610, 58863, 58864, 58865, 9066, 9067, 9068, 9069, 9070, 9071, 9072, 9073, 9074, 9075, 9076, 9077, 9078, 9079, 9080, 9081, 9082, 9083, 9084, 9085, 9086, 9087, 9088,
  qmcharge=-1,
  qm_theory='PM3',
  qmcut=12.0,
  qmshake=0,
  adjust_q=1
/
&quick
  method = 'B3LYP',
  basis = 'def2-svp',
/

The production simulation was run with:

mpirun -np 2 sander.quick.cuda.MPI -O -i step5_production.mdin -p step3_input.parm7 -c step4.1_equilibration.rst7 -o step5.mdout -r step5.rst7 -inf step5.mdinfo -x step5.nc

In this case, -np was limited to 2 because using more CPU ranks caused immediate GPU memory overflow on two RTX 4090 GPUs.

11. Odd Electron Count Error¶

If MPI_ABORT occurs and quick.out contains:

QMMM: System specified with odd number of electrons (107)
QMMM: but odd spin (1). You most likely have the charge of
QMMM: QM region (qmcharge) set incorrectly.

The selected QM region has an odd number of electrons. Review qmcharge and spin/multiplicity settings so they are chemically consistent with the selected QM region.

Example input after adding explicit charge and spin settings:

&qmmm
  iqmatoms=3729, 45573, 45574, 45575, 3728, 3730, 3731, 598, 596, 597, 45624, 45625, 45626, 56283, 56284, 56285, 58281, 58282, 58283, 58590, 58591, 58592, 9067, 9068, 9066, 9089, 9106, 9109, 9110, 9107, 9108,
  qmcharge=-1,
  spin=1,
  qm_theory='quick',
  qmcut=12.0,
  qmshake=0,
  adjust_q=1
/
&quick
  method = 'B3LYP',
  basis = 'def2-svp',
/

12. Semiempirical AM1 Fallback¶

In the reference workflow, the B3LYP/def2-SVP QUICK run was estimated at about 6,123.8 hours. Because of this cost, the QM/MM MD was run with the semiempirical method qm_theory='AM1'.

Production flags:

&qmmm
  iqmatoms=3729, 45573, 45574, 45575, 3728, 3730, 3731, 598, 596, 597, 45624, 45625, 45626, 56283, 56284, 56285, 58281, 58282, 58283, 58590, 58591, 58592, 9067, 9068, 9066, 9089, 9106, 9109, 9110, 9107, 9108,
  qmcharge=-1,
  spin=1,
  qm_theory='AM1',
  qmcut=12.0,
  qmshake=1,
  adjust_q=1,
  qm_ewald=1, qm_pme=1,
/

Run on CPU with 28 cores:

mpirun -np 28 sander.MPI -O -i step5_production.mdin -p step3_input.parm7 -c step4.1_equilibration.rst7 -o step5.mdout -r step5.rst7 -inf step5.mdinfo -x step5.nc

13. Error: QM Region + Cutoff Larger than Box¶

During QM/MM with AMBER, the following error can appear:

****************************************************
ERROR: QM region + cutoff larger than box dimension:
QM-MM Cutoff = 35.0000
 Coord   Lower     Upper    Size    Radius of largest sphere inside unit cell
   X   -40.743    39.788    80.531    46.761
   Y   -40.855    52.676    93.531    46.761
   Z   -42.214    43.047    85.261    46.761
****************************************************
SANDER BOMB in subroutine QM_CHECK_PERIODIC<qm_mm.f>
QM region + cutoff larger than box
cannot continue, need larger box.

AMBER checks whether the QM region fits inside the simulation box using:

• radius of the QM region
• QM/MM cutoff, qmcut
• smallest box dimension

The required condition is:

QM_radius + qmcut < half_box_dimension

where:

half_box_dimension = smallest_box_dimension / 2

Example:

• qmcut = 35 angstrom
• QM_radius = 8.147 angstrom
• smallest box dimension = 80.531 angstrom

Then:

QM_radius + qmcut = 43.147 angstrom
half_box_dimension = 80.531 / 2 = 40.265 angstrom

Since 43.147 > 40.265, the QM region plus cutoff does not fit in the box.

Reducing qmcut to 30 gives:

QM_radius + qmcut = 38.147 angstrom
half_box_dimension = 40.265 angstrom

Now 38.147 < 40.265, so the simulation can proceed.

14. Check the QM Radius Automatically¶

Save the following script as calc_qm_radius.py:

import sys
import math

fn = sys.argv[1]
qmcut = float(sys.argv[2])
min_box_dimension = float(sys.argv[3])

iqm_list = [
    3729, 45573, 45574, 45575, 3728, 3730, 3731, 598, 596, 597,
    45624, 45625, 45626, 56283, 56284, 56285, 58281, 58282, 58283,
    58590, 58591, 58592, 9067, 9068, 9066, 9089, 9106, 9109,
    9110, 9107, 9108,
]

atoms = []
with open(fn) as f:
    for line in f:
        if line.startswith("ATOM") or line.startswith("HETATM"):
            serial = int(line[6:11].strip())
            x = float(line[30:38].strip())
            y = float(line[38:46].strip())
            z = float(line[46:54].strip())
            resseq = line[22:26].strip()
            atoms.append((serial, int(resseq) if resseq.isdigit() else None, x, y, z))

coords = [(x, y, z) for (serial, resseq, x, y, z) in atoms if serial in iqm_list]
if not coords:
    coords = [(x, y, z) for (serial, resseq, x, y, z) in atoms if resseq in iqm_list]

if not coords:
    print("No IQM atom found. Check whether iqm_list contains atom serials or residue IDs.")
    print("Total atoms in PDB:", len(atoms))
    sys.exit(1)

cx = sum(x for x, y, z in coords) / len(coords)
cy = sum(y for x, y, z in coords) / len(coords)
cz = sum(z for x, y, z in coords) / len(coords)

dists = [
    math.sqrt((x - cx) ** 2 + (y - cy) ** 2 + (z - cz) ** 2)
    for x, y, z in coords
]

qm_radius = max(dists)
required_half_box = qm_radius + qmcut
required_box_size = 2 * required_half_box
half_box_dimension = min_box_dimension / 2

print("QM atom count:", len(coords))
print("QM centroid: {:.3f} {:.3f} {:.3f}".format(cx, cy, cz))
print("QM radius: {:.3f} angstrom".format(qm_radius))
print("Required half box: {:.3f} angstrom".format(required_half_box))
print("Required box size: {:.3f} angstrom".format(required_box_size))

if required_half_box > half_box_dimension:
    print(
        "WARNING: QM region + qmcut ({:.2f} angstrom) > half box ({:.2f} angstrom)".format(
            required_half_box, half_box_dimension
        )
    )
    print("Solution: increase the box size or reduce qmcut.")
else:
    print(
        "OK: QM region + qmcut = {:.2f} angstrom fits in a box of {:.2f} angstrom".format(
            required_half_box, min_box_dimension
        )
    )

Run:

python calc_qm_radius.py step3_input.pdb 35 80.531
python calc_qm_radius.py step3_input.pdb 30 80.531

Use the output to choose a qmcut that fits inside the simulation box.