lamm-mit/qmmm-adaptive

QM/MM Adaptive Dynamics Simulations

Overview

QM/MM (Quantum Mechanics/Molecular Mechanics) hybrid simulations enable accurate computational investigation of enzyme mechanisms, chemical reactions in biological environments, and reaction pathways. This skill combines quantum mechanical accuracy for reactive regions with classical force fields for the surrounding protein and solvent, enabling study of phenomena inaccessible to pure classical or pure quantum approaches.

Adaptive sampling techniques (metadynamics, umbrella sampling, accelerated MD) overcome energy barriers and efficiently explore reaction coordinates, making QM/MM practical for computational drug discovery and enzyme engineering.

Core Capabilities

1. QM/MM Setup and Equilibration

Hybrid System Preparation:

Build QM/MM systems with defined quantum and classical regions:

from qmmm import QMMMSystem

# Define QM region (reaction center)
qm_atoms = ["C1", "C2", "N3", "O4"]  # Reactive atoms
qm_method = "MOPAC"  # or ORCA, TeraChem for higher accuracy

# Define MM region (protein + solvent)
mm_forcefield = "AMBER14"

# Create hybrid system
system = QMMMSystem(
    structure="enzyme_complex.pdb",
    qm_atoms=qm_atoms,
    qm_method=qm_method,
    mm_forcefield=mm_forcefield,
    qm_charge=0,
    qm_multiplicity=1,
    buffer_distance=15  # Å from QM center
)

# Equilibrate in stages
system.equilibrate_qm_region(steps=5000)     # Fix MM atoms
system.equilibrate_hybrid(steps=10000)       # Allow coupling
system.equilibrate_full(steps=50000)         # Full system

Key Parameters:

• QM region size (3-30 atoms typically)
• Link atom treatment (capping vs. effective)
• Electrostatic embedding (mechanical vs. electronic)
• Boundary handling

2. Metadynamics - Explore Reaction Coordinate

Free Energy Landscape Calculation:

Map reaction pathways and identify transition states:

# Define collective variable (reaction coordinate)
cv = ReactionCoordinate(
    atoms=["C1", "C2"],
    type="distance"  # or angle, dihedral, custom
)

# Run metadynamics
metad = Metadynamics(
    system=qmmm_system,
    cv=cv,
    sigma=0.1,           # Width of Gaussian hills
    height=5.0,          # Energy height (kcal/mol)
    stride=100,          # Add hill every N steps
    temperature=300,     # K
    timestep=1.0         # fs
)

trajectory = metad.run(steps=100000)  # 100 ps with ~1000 hills

# Get free energy profile
free_energy = trajectory.get_free_energy()
barriers = free_energy.analyze_barriers()

Outputs:

• Free energy profile along reaction coordinate
• Transition state geometry and energy
• Multiple transition pathways (if present)
• Diffusion coefficient along coordinate

3. Umbrella Sampling - Precise Barrier Calculation

Constrained MD for PMF Calculation:

Systematically sample reaction coordinate windows:

# Define windows along reaction coordinate
windows = [(1.5 + i*0.1) for i in range(20)]  # 1.5 Å to 3.5 Å

umbrella = UmbrellaSampling(
    system=qmmm_system,
    cv=cv,
    windows=windows,
    force_constant=100.0,  # kcal/mol/Ų
    timestep=1.0,
    temperature=300
)

# Run each window (can be parallelized)
umbrella.run_windows(
    equilibration=5000,   # 5 ps equilibration per window
    production=20000,     # 20 ps production per window
    save_frequency=100,   # Save every 0.1 ps
    n_parallel=4          # Run 4 windows in parallel
)

# Analyze using WHAM
free_energy = umbrella.analyze_wham()

Advantages:

• Higher resolution along coordinate
• Precise barrier heights (±0.5 kcal/mol)
• Error estimation via bootstrap
• Umbrella integration for rate constants

4. Accelerated MD (aMD) - Escape Barriers

Enhanced Sampling via Potential Energy Boost:

Accelerate sampling by reducing energy barriers:

amd = AcceleratedMD(
    system=qmmm_system,
    boost_type="dual",        # Boost both dihedral + total energy
    e_threshold=10.0,         # Boost below E+10 kcal/mol
    alpha_dihedral=5.0,       # Dihedral boost strength
    alpha_total=10.0,         # Total energy boost strength
    timestep=2.0,
    temperature=300
)

trajectory = amd.run(steps=200000)  # 400 ps with acceleration

# Reweight trajectory to remove boost
free_energy = trajectory.reweight_by_boost_energy()

Applications:

• Escape local energy minima
• Sample multiple conformational states
• Identify transient intermediates
• Transition pathway discovery

5. Transition State Finding

Locate and Characterize TS Structures:

# From metadynamics trajectory
ts_finder = TransitionStateFinder(
    metad_trajectory=metad_traj,
    free_energy_surface=fe_profile
)

# Find TS structure (highest energy on MEP)
ts_structure = ts_finder.find_ts(
    refine=True,
    opt_method="L-BFGS"
)

# Characterize TS
ts_analysis = ts_structure.analyze()
print(f"TS Energy: {ts_analysis['energy']:.2f} kcal/mol")
print(f"Barrier height: {ts_analysis['barrier_forward']:.2f} kcal/mol")
print(f"Barrier height (reverse): {ts_analysis['barrier_reverse']:.2f} kcal/mol")
print(f"TS geometry verified: {ts_analysis['verified']}")

6. Enzyme Mechanism Studies

Multi-Step Reaction Mechanisms:

Study sequential bond breaking/forming:

# Example: Serine protease reaction mechanism
# Step 1: Nucleophilic attack (Ser195 on substrate)
# Step 2: Tetrahedral intermediate formation
# Step 3: Acyl-enzyme formation
# Step 4: Water activation and acyl hydrolysis

steps = [
    {
        "name": "Nucleophilic attack",
        "cv": Distance(["Ser195-OG", "C1-substrate"]),
        "expected_barrier": 12.0  # kcal/mol
    },
    {
        "name": "Proton transfer",
        "cv": Distance(["His57-NE", "Proton"]),
        "expected_barrier": 8.0
    },
    {
        "name": "Acylation",
        "cv": Distance(["Ser195-C", "Peptide-N"]),
        "expected_barrier": 5.0
    },
    {
        "name": "Deacylation",
        "cv": Distance(["Water-O", "Acyl-C"]),
        "expected_barrier": 15.0
    }
]

# Run each step with QM/MM
for step in steps:
    result = qmmm.study_step(step)
    print(f"{step['name']}: {result['barrier']:.1f} kcal/mol")

QM Methods (QM Region)

Fast (for MD)

• MOPAC PM6-D3H4X - Semi-empirical, very fast
• DFTB+ - Density functional tight binding
• GFN2-xTB - Extended tight binding

Balanced (for TS finding)

• ORCA B3LYP/6-31G(d) - DFT, good accuracy
• TeraChem B3LYP - GPU-accelerated DFT
• MOPAC PM7 - Improved semi-empirical

High Accuracy (for benchmarks)

• ORCA DLPNO-CCSD(T) - Wave function theory
• MOLPRO - High-level quantum chemistry

MM Methods (MM Region)

• AMBER14 - Standard biomolecular FF
• CHARMM36 - High accuracy for proteins
• OPLS-AA - Good for organics

Sampling Techniques

| Technique | Speed | Accuracy | Use Case | |-----------|-------|----------|----------| | Metadynamics | Fast | Good | Explore landscape, find TS guess | | Umbrella Sampling | Medium | Excellent | Precise barriers, PMF convergence | | aMD | Very Fast | Medium | Escape local minima, sample states | | Targeted MD | Fast | Good | Pull along specific pathway | | Steered MD | Very Fast | Medium | Non-equilibrium pulling experiments |

Use Cases

Enzyme Mechanism Elucidation:

• Identify transition states
• Calculate barrier heights
• Propose mechanistic models
• Explain selectivity and reactivity

Drug Metabolism Prediction:

• Predict metabolic pathways
• Identify reactive metabolites
• Calculate oxidation barriers
• Assess bioactivation risk

Protein Engineering:

• Design improved catalysts
• Predict activity mutations
• Optimize substrate selectivity
• Lower activation barriers

Chemical Reactivity in Proteins:

• Metal coordination mechanisms (metalloproteases, P450 oxidations)
• Radical chemistry (peroxidases, LOX)
• Photochemistry (retinal isomerization, photolyase)
• Hydrolysis mechanisms (proteases, esterases)

Integration with Other Skills

Input:

• Protein structures from pdb skill
• Ligand/substrate SMILES from pubchem
• Enzyme sequences from uniprot
• Literature mechanistic insights from pubmed

Output:

• TS structures for further analysis by ase
• Barriers for publication/comparison
• Mechanistic models for hypothesis testing
• Reactive intermediates for openmm MD validation

Performance Characteristics

Computational Cost:

• QM/MM equilibration: ~1-4 hours (500 atoms, 100 ps)
• Metadynamics (100 ps): ~8-24 hours
• Umbrella sampling (20 windows × 20 ps): ~40-120 hours
• aMD (400 ps): ~4-12 hours

Hardware:

• CPU: Excellent for MOPAC/DFTB based QM/MM
• GPU: Required for TeraChem or faster convergence
• Parallelization: Windows/replicas parallelize perfectly

Limitations

• QM region size limited (~30 atoms for DFT)
• Boundary effects at QM/MM interface
• Charge transfer between QM/MM regions approximated
• Polarization effects sometimes underestimated
• Requires quality starting structure

Example Workflow

# 1. Prepare QM/MM system
python qmmm_setup.py --pdb enzyme_complex.pdb \
    --qm-atoms "catalytic_residues.txt" \
    --qm-method MOPAC --force-field AMBER14

# 2. Equilibrate
python qmmm_equilibrate.py --system prepared_qmmm.inp \
    --temperature 300 --steps 50000

# 3. Run metadynamics to explore
python qmmm_metad.py --system equilibrated.inp \
    --cv "distance C1 C2" --duration 100 \
    --output metad_landscape.txt

# 4. Find transition state
python qmmm_ts_finder.py --metad-traj metadynamics.dcd \
    --refine true --output ts_structure.pdb

# 5. Run umbrella sampling for precise barrier
python qmmm_umbrella.py --system equilibrated.inp \
    --cv "distance C1 C2" --windows 20 \
    --window-width 0.1 --output pmf.txt

# 6. Analyze mechanism
python qmmm_mechanism.py --ts-structure ts_structure.pdb \
    --umbrella-pmf pmf.txt --report mechanism.md

References

• QM/MM review: Senn & Thiel (2007) Angew. Chem.
• Metadynamics: Laio & Parrinello (2002) PNAS
• Umbrella sampling: Torrie & Valleau (1977) J. Comp. Phys
• aMD: Hamelberg et al. (2004) J. Chem. Phys
• ORCA: Neese (2012) Wiley Interdiscip. Rev. Comput. Mol. Sci.
• TeraChem: Ufimtsev & Martinez (2008) J. Chem. Phys