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.

Under active development

SOC-NAMD and SOC-NAMD-QMMM are included in OpenQP 1.3.0 but remain experimental. A runnable input does not imply that the present force and coupling approximations are ready for production scientific conclusions. Use this workflow for method development and numerical verification, and report the selected basis and the limitations below with every result.

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 current development and numerical verification, soc_basis=mch is the preferred path because it 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). This preference does not change the experimental status of SOC-NAMD.

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 (particle-mesh Ewald, cutoff=PME, ESPF-PME electrostatics) is available to the same-spin FSSH driver; the SOC-NAMD drivers run isolated clusters (cutoff=NoCutoff). See the [qmmm] page and References.

Scope and limitations

The validated SOC-NAMD-QMMM setting is a whole-molecule QM region: a solvated chromophore with the whole chromophore in the QM region. Covalent QM/MM boundaries (hydrogen link atoms) are also accepted, in same-spin and SOC-NAMD, with schemes that do not use analytic NAC; see Link atoms.

In the periodic branch (cutoff=PME, same-spin FSSH only) the QM-image term is iterated to self-consistency with the relaxed ESPF charges of the propagated state, so the force integrated is the derivative of the reported energy of that state; each iteration is logged as QM-image field, active-state iteration, and after a surface hop the field is re-iterated for the new state before its force is integrated. The relaxed charges inherit the Z-vector residual, so periodic MRSF dynamics should set [tdhf] zvconv = 1e-8 (the input checker warns otherwise): the image loop then converges in two gradient evaluations per step instead of four or five. The spin-adiabatic SOC-NAMD state is a mixture of MCH states without per-iteration relaxed charges, so SOC-NAMD rejects periodic boxes (use cutoff=NoCutoff, the default). The tight-binding (method=dftb/xtb) NAMD path supports non-periodic clusters only.

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 water cluster

A complete deck for a chromophore solvated in a TIP3P water cluster (SOC-NAMD is non-periodic; the same deck with soc=false may use cutoff=PME). The whole chromophore is the QM region (qm_atoms); the water is MM.

.oqp:

mrsf(nstate=3)/bhhlyp/6-31g*
namd(S1,scheme=custom,tdc=npi,rescale=isotropic,thrshe=0.367493,frustrated=reflect,
     soc=true,soc_basis=mch)
md(nstep=200,dt=0.5)
qmmm(forcefield_files="amber14-all.xml,amber14/tip3p.xml",qm_atoms="0-14",rigidwater=true)
geom="chromophore_water.pdb"

The PDB path appears only once. OpenQP derives qmmm.pdb_file from the PDB prefix of geom. Give the QM region once: either the 0-based qm_atoms above or the 1-based selector after the PDB name (geom="chromophore_water.pdb 1-15" selects the same atoms). The Python API likewise derives both qmmm.pdb_file and qmmm.qm_atoms from job.molecule("file.pdb <selector>") when they are omitted from job.qmmm(...).

Python — 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 1-15", basis="6-31g*")
job.theory.mrsf(functional="bhhlyp", nstate=3)   # ROHF triplet reference + MRSF

# ESPF QM/MM embedding in a TIP3P water cluster
job.qmmm(
    forcefield=["amber14-all.xml", "amber14/tip3p.xml"],
    embedding="electrostatic",
    rigidwater=True,
)

# SOC-NAMD on the spin-adiabatic manifold, exact-gradient (MCH) force basis
job.workflow.namd(
    scheme="custom",
    tdc="npi",
    rescale="isotropic",
    thrshe=0.367493,
    frustrated="reflect",
    soc=True,
    soc_basis="mch",
    init_state="S1",
)
job.workflow.md(nstep=200, dt=0.5, temperature=300.0)

mol = job.run()

Legacy .inp:

[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.367493
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
embedding        = electrostatic
rigidwater       = True

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 cutoff. thrshe=0.367493, approximately 10 eV, is the default for both same-spin and SOC dynamics and blocks large-gap transitions outside the intended local crossing region.
  • QM region. qm_atoms is normally a whole molecule; a covalent cut is capped with hydrogen link atoms and needs a scheme without analytic NAC (see Scope and limitations).
  • Periodicity. SOC-NAMD runs isolated clusters (cutoff=NoCutoff). cutoff=PME selects the periodic ESPF-PME branch for a solvated box and is available to the same-spin FSSH driver (soc=false) only.

Trajectory records, restart, and energy conservation

SOC-NAMD writes the same restartable record family as same-spin NAMD:

  • <project>.namd.trj is an appendable packed binary trajectory intended for memory-mapped analysis. Besides geometry, velocity, energy, population, hop, and full-precision random fields, SOC records contain the complex spin-adiabatic overlap and anti-Hermitian TDC as real/imaginary components and identify the active adiabatic or mch representation.
  • <project>.namd.restart.npz is an atomic numerical checkpoint. It preserves the previous SOC eigensystem and singlet/triplet response vectors required for the next overlap and rejects incompatible shapes, dtypes, or non-finite data.
  • <project>.namd.restart.oqp is a directly runnable continuation input. It freezes a date-derived random seed and reconnects the checkpoint and packed trajectory without treating nstep as an additional step count.

The main log prints a concise energy-conservation table. The default nve_policy=warn applies to SOC as well as same-spin NVE trajectories; use nve_policy=error only after tolerances have been calibrated. Under NVE dynamics with the default full-ESPF electrostatics, E_tot should show no systematic drift. Read the packed record without loading it fully:

from oqp.library.namd import read_namd_trajectory

metadata, trj = read_namd_trajectory("chromophore_socnamd_qmmm.namd.trj")
time = trj["time_fs"]
energy = trj["e_tot_hartree"]
soc_tdc = trj["overlap_tdc_au"] + 1j * trj["overlap_tdc_imag_au"]

The independent TD-Baeck–An NACME magnitude diagnostic is defined only for same-spin state pairs. Its global default is therefore contextually disabled for SOC; the SOC record itself retains the complete complex overlap/TDC.

If a SOC force path produces a slow drift, econs=true rescales velocities each step to conserve E_tot as a temporary stabilizer. For development and numerical verification, prefer soc_basis=mch because it uses exact active-root MCH gradients; see the force-basis note above. Do not interpret this as production-readiness of the complete SOC-NAMD method.

References