Skip to content

SOC-NAMD-QMMM: excited-state dynamics with intersystem crossing in an environment

This tutorial walks you through a spin-orbit-coupled nonadiabatic molecular dynamics simulation with QM/MM embedding — "SOC-NAMD-QMMM" — from the physics behind the acronym to a deck you can run and the numbers it prints.

By the end you will understand:

  • what each ingredient (MRSF-TDDFT, NAMD, SOC, QM/MM) contributes and why you would combine them,
  • how the pieces map onto the [scf], [tdhf], [md], and [qmmm] sections of an OpenQP input,
  • how to run the calculation and read the trajectory it produces, and
  • how to grow the minimal demo into a production run.

The runnable inputs are in inputs/; the reference manual for every keyword is the OpenQP manual, in particular the SOC-NAMD-QMMM workflow page and the [md] / [qmmm] keyword pages.

1. The problem, and why it needs four ingredients

Imagine a molecule that absorbs light, is promoted to an excited singlet state, and then — instead of simply fluorescing back down — crosses into a triplet state and does its chemistry from there. This intersystem crossing (ISC) is at the heart of phosphorescence, photosensitizers, triplet-triplet annihilation, and much of photobiology. Simulating it from first principles, in a realistic environment, needs four things working together. Let's build up the acronym one ingredient at a time.

MRSF-TDDFT — the excited states

You need excited-state energies and, crucially, gradients (forces) to move the nuclei. Ordinary TDDFT struggles near conical intersections and for open-shell/diradical character. MRSF-TDDFT (Mixed-Reference Spin-Flip TDDFT) fixes this: it spin-flips out of a high-spin (triplet) reference and mixes two reference determinants, which removes the spin contamination of plain spin-flip TDDFT and gives balanced singlet and triplet states — exactly what you need when those states are about to cross. In the deck this is [scf] multiplicity=3 type=rohf (the triplet reference) plus [tdhf] type=mrsf.

NAMD — letting the nuclei hop between surfaces

Once a molecule is excited, the Born-Oppenheimer approximation (nuclei on a single fixed electronic surface) breaks down: the wavepacket can switch electronic states as it moves. Nonadiabatic molecular dynamics models this. OpenQP uses Tully's fewest-switches surface hopping (FSSH): the nuclei move classically on one "active" state, an electronic amplitude is propagated alongside, and at each step the trajectory may hop to another state with a probability set by the nonadiabatic coupling. runtype=namd selects this; the [md] section controls the integrator, the active state, and the hopping.

SOC — turning on singlet ↔ triplet transitions

Plain NAMD only hops between states of the same spin (internal conversion, S→S). A singlet cannot become a triplet without the relativistic spin-orbit coupling operator, which mixes spin. Turn it on with [md] soc=True. OpenQP then propagates on a spin-adiabatic ("SHARC-like") manifold built from the singlets and the individual triplet sublevels — for nstate singlets and nt triplets that is ns + 3·nt states (each triplet contributes 3 M_s sublevels). A hop between a singlet block and a triplet block is an intersystem crossing, driven by the computed SOC matrix elements. This is what makes it SOC-NAMD.

QM/MM — putting the molecule in its environment

Photochemistry rarely happens in vacuum. QM/MM treats the interesting molecule quantum-mechanically (QM) and its surroundings — solvent, protein, crystal — with a cheap classical force field (MM). OpenQP couples the two with ESPF electrostatic embedding: the MM point charges polarize the QM density through the electrostatic-potential-fitted (ESPF) operator, and the QM density reacts back on the MM atoms, with an analytic, energy-conserving gradient. qmmm_flag=True plus the [qmmm] section (a PDB, a force field, and which atoms are QM) turn this on; the MM engine is OpenMM.

Put together, SOC-NAMD-QMMM = excited-state surface-hopping dynamics of an MRSF-TDDFT chromophore, with singlet↔triplet intersystem crossing, embedded in an explicit MM environment.

2. The system for this tutorial

We use formaldehyde (H₂CO) in a small cluster of 5 water molecules. It is the smallest system that shows every ingredient:

  • Formaldehyde's low-lying n→π* singlet (S₁) and triplet (T₁) states lie close in energy — a classic ISC playground.
  • 5 waters give a real, polarizing MM environment without being slow.
  • It runs in seconds, so you can iterate.

The QM region is the 4 formaldehyde atoms (qm_atoms=0-3 in the PDB); the 5 waters are MM.

Note. The shipped deck uses [md] nstep=1 — a single nuclear step — so the example finishes instantly and can be used as a smoke test. To see actual dynamics and hops you raise nstep (see §6).

3. Prerequisites

pip install openqp        # the QM engine + CLI
pip install openmm        # the MM backend (required for QM/MM)

Check both import:

python -c "import oqp; import openmm; print('ok')"

If OpenMM is missing, QM/MM decks are reported SKIPPED rather than run.

4. The input, line by line

Open inputs/h2co-water_soc-namd-qmmm.oqp. It is three lines including the geometry. Here is what each one is doing.

The route — the QM electronic model

mrsf(nstate=2)/bhhlyp/6-31g*

mrsf is MRSF-TDDFT, and it carries its own reference: a triplet ROHF determinant. That reference is not the state you care about — it is the mathematical starting point that MRSF spin-flips from to reach the balanced singlet/triplet manifold — and because the model implies it, there is no way to pair MRSF with a reference that cannot support it.

nstate=2 requests two singlet roots (S₀, S₁). With soc=true in the driver the code also forms the triplets and their sublevels automatically, giving the spin-adiabatic manifold described in §1. bhhlyp is a common half-and-half functional for MRSF.

The geometry — QM subsystem and MM environment in one line

geom="formaldehyde_water.pdb 0-3"

The PDB holds all atoms (QM + MM); the trailing 0-3 carves out the QM region by 0-based PDB index. Everything else is MM. One line therefore fixes both the QM geometry and the QM/MM partition — the legacy deck spelled the same thing three times (an inline system block, [qmmm] pdb_file, and [qmmm] qm_atoms), which is three chances for them to disagree.

namd(...) — the dynamics and the hopping

This is the heart of a NAMD run.

namd(nstep=1,dt=0.25,active=5,substep=50,seed=3,trivial=true,soc=true)

It shares the first line with the route in the deck. Only what differs from the defaults is written:

option meaning
nstep=1 number of nuclear steps (raise for a real trajectory)
dt=0.25 nuclear time step (fs)
active=5 initially populated spin-adiabatic state
substep=50 electronic sub-steps per nuclear step
seed=3 fixed RNG seed → reproducible hops
trivial=true detect trivial (weakly-avoided) crossings
soc=true spin-orbit coupling ON → allows ISC

Already the defaults, and therefore not written: initial velocities sampled from a Maxwell distribution at 300 K (init_temp=300, velocity=maxwell), energy-based decoherence (decoherence=edc, Granucci-Persico 2007), and the thrshe=0.1 Hartree gap gate.

Key ideas:

  • active is which state the trajectory starts on. On the spin-adiabatic manifold the states are ordered by energy across singlets and triplet sublevels, so active=5 starts partway up the manifold (a bright singlet, here) rather than the ground state. Because active already names the state completely, the driver takes no separate state label, and the gradient the propagator needs is selected from it each step.
  • substep matters because the electronic amplitude oscillates much faster than the nuclei move; the time-derivative couplings are integrated on this finer grid between nuclear steps.
  • decoherence=edc (the default) corrects a well-known FSSH pathology (over-coherence) so populations relax physically. trivial=true stops the trajectory from missing a hop at a very sharp, weakly-avoided crossing.
  • soc=true is the switch that turns NAMD into SOC-NAMD: without it you get internal conversion only; with it, singlet↔triplet ISC is possible.
  • thrshe is a safety gate: near the Franck-Condon point unrestricted hopping would permit unphysical hops down to S₀; the default 0.1 Hartree is the recommended value.

qmmm(...) — the environment

qmmm(forcefield_files="formaldehyde.xml tip3p.xml")

Supplying qmmm(...) is itself what turns QM/MM on; the PDB named in geom is propagated into the section for you. The MM atoms are parameterized by the forcefield_files. formaldehyde.xml only needs to supply the QM atoms' Lennard-Jones parameters — their electrostatics come from ESPF, not from fixed MM charges. The defaults suit this isolated cluster: cutoff=NoCutoff and full ESPF coupling (embedding=electrostatic). For a solvated periodic box you would add cutoff=PME, which turns on the particle-mesh-Ewald branch of the embedding.

Everything not written — the extended-Huckel starting orbitals, the SCF thresholds, the XC grid — is a default. Exact section calls such as scf(conv=1e-8) remain the escape hatch for every control the format does not name directly.

The QM region here is a whole molecule. Cutting a covalent bond (carving a fragment out of a larger molecule) is a covalent boundary — supported by the single-point and ground-state QM/MM paths, but not by namd(...). See the covalent-boundary docs.

5. The same run, in Python

The compact OpenQP scripting interface builds the identical calculation — job.theory.mrsf(...) sets the reference and states, job.qmmm(...) turns on the embedding, and job.workflow.namd(soc=True, ...) selects runtype=namd and fills the [md] section:

from oqp.openqp import OpenQP

job = OpenQP("h2co-water_soc-namd-qmmm", silent=1)

job.molecule("""
    C   0.000000   0.000000   0.000000
    O   0.000000   0.000000   1.203000
    H   0.000000   0.943000  -0.589000
    H   0.000000  -0.943000  -0.589000
""", charge=0)

job.theory.mrsf(functional="bhhlyp", basis="6-31g*", nstate=2, multiplicity=3)

job.qmmm(
    pdb_file="formaldehyde_water.pdb",
    forcefield=["formaldehyde.xml", "tip3p.xml"],
    qm_atoms="0-3", cutoff="NoCutoff", embedding="electrostatic",
)

job.workflow.namd(
    soc=True, nstep=1, dt=0.25, active=5, substep=50,
    init_temp=300, velocity="maxwell", seed=3,
    decoherence="edc", trivial=True, thrshe=0.1,
)

mol = job.run()

The full script is inputs/h2co-water_soc-namd-qmmm.py.

6. Running it

From the inputs/ folder (so the PDB/force-field files resolve), either style:

cd soc-namd-qmmm/inputs
openqp h2co-water_soc-namd-qmmm.oqp     # input-file style
python h2co-water_soc-namd-qmmm.py      # Python-API style

What OpenQP does at each nuclear step:

  1. Build the MM electrostatic potential at the QM atoms (via OpenMM) and hand it to the ESPF operator.
  2. Run the embedded MRSF-TDDFT calculation → singlet and triplet energies and the SOC matrix elements between them.
  3. Form the spin-adiabatic states and the gradient of the active state.
  4. Integrate the electronic amplitude over substep sub-steps using the time-derivative couplings, and evaluate the fewest-switches hopping probabilities (now including S↔T channels because soc=True).
  5. Possibly hop (with decoherence + trivial-crossing handling), then advance the nuclei by dt under the active-state force (QM atoms by the QM/ESPF force, MM atoms by the MM + coupling force).

The run writes a log (<project>.log) and a trajectory.

7. Reading the output and going further

The .log records, per step: the electronic-state energies, the active state, the hopping probabilities, the SOC couplings, and the kinetic/potential/total energy. A hop shows up as the active-state index changing between steps; an intersystem crossing is a hop whose old and new states belong to different spin blocks.

To turn the smoke test into a real simulation:

Want to… Change
See actual dynamics/hops Raise [md] nstep (e.g. 500–2000)
Start on a different state Change [md] active
Sample many trajectories Vary [md] seed and average
Model a solvated box [qmmm] cutoff=PME with a periodic PDB
More/other excited states Raise [tdhf] nstate
Internal conversion only [md] soc=False

On the SOC gradient. With soc=True the default hops on the spin-adiabatic manifold using a weighted-MCH diagonal gradient. If you need the exact active-root gradient in the molecular-Coulomb-Hamiltonian (MCH) basis, set [md] soc_basis=mch. See the workflow page for the trade-offs.

8. Recap

You combined four ideas into one simulation:

  • MRSF-TDDFT gave balanced singlet/triplet excited states and gradients,
  • NAMD (FSSH) let the nuclei hop between electronic surfaces,
  • SOC opened the singlet↔triplet (intersystem-crossing) channel, and
  • QM/MM (ESPF) embedded the chromophore in an explicit environment.

That is SOC-NAMD-QMMM. From here, the manual gives the full input contract, the periodic (PME) setup, and the compact job.qmmm(...) / job.workflow.namd(...) Python API.

References