Skip to content

SOC-NAMD-QMMM

SOC-NAMD-QMMM is nonadiabatic molecular dynamics that combines three capabilities in one propagation loop:

  • NAMD — Tully fewest-switches surface hopping (FSSH) on MRSF-TDDFT states;
  • SOC — spin-orbit-coupled intersystem crossing (ISC), so hops occur on the spin-adiabatic manifold of mixed singlet/triplet states;
  • QM/MM — the MRSF-TDDFT chromophore is embedded in a classical (OpenMM) MM environment through the ESPF operator.

It is dispatched by runtype=namd with [md] soc=true and [input] qmmm_flag=true, and is configured through the [md] and [qmmm] sections. This page introduces the method and gives a complete, runnable input deck.

Development preview

This workflow documents the implementation branch in OpenQP PR #205. It is not part of OpenQP 1.2.0; use that source branch or a later release that includes runtype=namd.

Overview and theory

FSSH nonadiabatic dynamics

Classical nuclei propagate on one active potential energy surface with velocity-Verlet while the electronic amplitudes evolve; stochastic hops between surfaces reproduce population transfer. OpenQP implements FSSH generalized to an arbitrary number of states on MRSF-TDDFT surfaces, with energy-based decoherence (EDC), finite-difference time-derivative couplings (TDC), and trivial-crossing (diabatic) following. See Tully FSSH.

SOC-NAMD (intersystem crossing)

With [md] soc=true, the spin-orbit MRSF-TDDFT driver (soc_mrsf) builds and diagonalizes the spin-orbit Hamiltonian

H = diag(E_MCH) + H_SOC

whose eigenvectors define the spin-adiabatic states. Here the MCH (molecular Coulomb Hamiltonian) basis is the set of spin-pure MRSF singlet and triplet states, and H_SOC are the spin-orbit matrix elements that couple them. For [tdhf] nstate singlets and the same number of triplets, the manifold has ns + 3*nt SOC states (each triplet contributes three Ms sublevels).

Surface hopping is carried out on this spin-mixed manifold, so a hop between a predominantly singlet and a predominantly triplet spin-adiabatic state is an intersystem-crossing event. Continuity of the eigenvectors from step to step is maintained by U-phase tracking (the phase/ordering of the eigenvector matrix U), and the electronic amplitudes are propagated with local diabatization.

Active-surface force. soc_basis selects the force basis. The default soc_basis=adiabatic propagates on spin-adiabatic SOC eigenstates and uses the weighted-MCH diagonal gradient: each contributing spin-pure (MCH) component carries its own gradient, weighted by its population in the active spin-adiabatic state, with components below grad_wthr dropped so the force stays continuous through strong spin mixing. The soc_du_dt_corr and soc_tdc_grad_corr flags are optional approximate corrections for this adiabatic force path.

For production trajectories, soc_basis=mch propagates in the spin-pure MCH basis and uses exact active-root MCH gradients; with QM/MM this selects the NAMD_SOC_MCH_QMMM driver. At an ISC hop, velocities are rescaled to conserve the total energy (including the ESPF embedding energy change).

ESPF QM/MM embedding

The MRSF-TDDFT QM region is embedded in the OpenMM MM environment through the electrostatic potential-fitted (ESPF) operator: the MM point charges enter the QM core Hamiltonian (polarizing the QM density), and the QM density reacts on the MM atoms through ESPF-fitted atomic charges. All electronic quantities — the reference SCF, the singlet and triplet MRSF states, and the SOC matrix — inherit the same embedding. The default full-ESPF electrostatics (embedding=electrostatic) give a finite-difference-exact analytic gradient and energy-conserving dynamics. Rigid-water SHAKE/RATTLE constraints are applied to the MM region inside the velocity-Verlet loop, so a normal (~0.5 fs) timestep can be used; QM atoms are never constrained. A periodic water box uses particle-mesh Ewald (cutoff=PME) with ESPF-PME electrostatics. See the [qmmm] page and References.

Scope and limitations

SOC-NAMD-QMMM builds the QM molecule from [qmmm] qm_atoms only, so whole-molecule QM regions are supported. Covalent QM/MM boundaries (hydrogen link atoms) in nonadiabatic dynamics are not yet available — single-point QM/MM and ground-state QM/MM MD do handle covalent boundaries (see Link atoms). Use a solvated chromophore in a periodic (PME) water box, with the whole chromophore in the QM region.

How the driver is selected

runfunc.compute_namd picks the surface-hopping class from qmmm_flag, soc, and soc_basis:

[input] qmmm_flag [md] soc [md] soc_basis Class
false false ignored NAMD (gas-phase FSSH)
false true adiabatic NAMD_SOC (gas-phase spin-adiabatic SOC-NAMD)
false true mch NAMD_SOC_MCH (gas-phase MCH-basis SOC-NAMD)
true false ignored NAMD_QMMM (FSSH + ESPF QM/MM)
true true adiabatic NAMD_SOC_QMMM (spin-adiabatic SOC-NAMD + ESPF QM/MM)
true true mch NAMD_SOC_MCH_QMMM (MCH-basis SOC-NAMD + ESPF QM/MM)

Example deck: SOC-NAMD-QMMM in a periodic water box

A complete deck for a chromophore solvated in a periodic TIP3P water box. The whole chromophore is the QM region (qm_atoms); the water is MM.

[input]
runtype    = namd
qmmm_flag  = true
method     = tdhf
functional = bhhlyp
basis      = 6-31g*

[scf]
type         = rohf
multiplicity = 3

[tdhf]
type   = mrsf
nstate = 3

[md]
soc        = true
soc_basis  = mch
active     = 1
init_state = S1
nstep      = 200
dt         = 0.5
thrshe     = 0.1
init_temp  = 300.0
grad_wthr  = 0.001

[qmmm]
pdb_file         = chromophore_water.pdb
forcefield_files = amber14-all.xml,amber14/tip3p.xml
qm_atoms         = 0-14
cutoff           = PME
embedding        = electrostatic
rigidwater       = True

The same job in the compact Python API — job.qmmm(...) enables ESPF QM/MM and job.workflow.namd(...) selects the surface-hopping run (see Run from Python):

from oqp.openqp import OpenQP

job = OpenQP("chromophore_socnamd_qmmm", silent=1)

# QM geometry + atom selection from the PDB (whole-molecule QM region)
job.molecule("chromophore_water.pdb 0-14", basis="6-31g*")
job.theory.mrsf(functional="bhhlyp", nstate=3)   # ROHF triplet reference + MRSF

# ESPF QM/MM embedding in a periodic TIP3P water box
job.qmmm(
    pdb_file="chromophore_water.pdb",
    forcefield=["amber14-all.xml", "amber14/tip3p.xml"],
    qm_atoms="0-14",
    cutoff="PME",
    embedding="electrostatic",
    rigidwater=True,
)

# SOC-NAMD on the spin-adiabatic manifold, exact-gradient (MCH) force basis
job.workflow.namd(
    soc=True,
    soc_basis="mch",
    init_state="S1",
    nstep=200,
    dt=0.5,
    thrshe=0.1,
    init_temp=300.0,
)

mol = job.run()

Notes on the deck:

  • Reference / states. SOC-NAMD requires an MRSF-TDDFT setup: a high-spin ROHF reference ([scf] type=rohf multiplicity=3) and [tdhf] type=mrsf. With [tdhf] nstate=3 the spin-adiabatic manifold has ns + 3*nt = 3 + 9 = 12 states, so [md] active may range 1..12.
  • Initial surface. init_state=S1 starts on the state of dominant S1 character and overrides active.
  • Force basis. soc_basis=mch uses exact active-root MCH gradients and selects NAMD_SOC_MCH_QMMM. Use soc_basis=adiabatic to test the spin-adiabatic weighted-gradient path and its optional correction flags.
  • Gap gate. thrshe=0.1 is the recommended SOC-NAMD value; the large default would allow spurious S0 hops at the Franck-Condon geometry.
  • QM region. qm_atoms must be a whole molecule (see Scope and limitations).
  • Periodicity. cutoff=PME selects the periodic ESPF-PME branch for a solvated box; use NoCutoff for an isolated cluster.

Outputs and energy conservation

The run writes a trajectory log with, per nuclear step, the time, the active state, the MCH/spin-adiabatic state energies, the total energy E_tot (potential + kinetic, including the ESPF embedding energy), and hopping events.

Under NVE dynamics with the default full-ESPF electrostatics, E_tot should be flat (no systematic drift) — this is the main check that the trajectory is physically meaningful. To verify, plot E_tot versus time from the trajectory log and confirm it fluctuates around a constant with no trend.

If a SOC force path produces a slow drift, econs=true rescales velocities each step to conserve E_tot as a temporary stabilizer. Prefer soc_basis=mch for production SOC-NAMD-QMMM trajectories because it uses exact active-root MCH gradients; see the force-basis note above.

References