Direct Reaction Field (DRF)

Standard QM/MM electrostatic embedding lets a QM region polarize in the static field of MM point charges, but the MM atoms themselves remain unresponsive: they cannot polarize in return when the QM electron density changes (e.g. on photoexcitation). The direct reaction field method of Thole and van Duijnen¹ repairs this by giving each MM atom an isotropic atomic polarizability \(\alpha_i\), then adding the induced-dipole response directly to the QM Hamiltonian — both to the one-electron core operator and to the two-electron Coulomb/exchange operators. The result is a single QM/MM Hamiltonian whose ground- and excited-state eigenstates are all orthogonal and whose state crossings are well-defined, without the root-flipping and excitation-energy errors that plague linear-response and state-specific polarizable continuum schemes.

TeraChem's implementation is the integral-exact DRF (also called DRF-2e-Pol or IEDRF) developed by Liu, Humeniuk, Glover² and extended to analytical gradients, TDDFT/CIS, and conical-intersection optimization by Humeniuk and Glover³. All polarization integrals are evaluated exactly through a dedicated GPU-accelerated kernel library — earlier DRF implementations expanded these integrals in Taylor series around atomic centers, which broke smoothness and precluded analytic gradients.

DRF can be combined with HF, DFT, TDDFT/RPA, CIS, CASCI, and CASSCF, and supports analytic gradients and non-adiabatic couplings for all of them.

Build requirement

DRF is only present in TeraChem builds compiled with the DRF2EPOL preprocessor flag. Without it, setting direct_reaction_field yes aborts with "TeraChem has not been compiled with DRF2EPOL". The official binary releases enable it by default.

OpenMM/Amber driver only

DRF reads atomic polarizabilities from an Amber prmtop file and only works with the OpenMM/Amber QM/MM driver (see OpenMM / Amber Driver). The internal-water and Amber File drivers are not supported. Setting direct_reaction_field yes without prmtop/qmindices aborts.

Theory

Polarization energy and the DRF Hamiltonian

The total charge distribution is split into electronic \(\hat\rho(\mathbf{r}) = \sum_a(-1)\delta(\mathbf{r}-\hat{\mathbf{r}}_a)\), QM nuclear, and classical MM-point-charge contributions. Each polarizable MM atom \(i\) at position \(\mathbf{R}_i\) carries an isotropic atomic polarizability \(\alpha_i\) and develops an induced dipole

\[ \mathbf{p}_i \;=\; \alpha_i\,\mathbf{E}(\mathbf{R}_i), \]

driven by the total electric field at \(\mathbf{R}_i\). That field has contributions from the free QM charges (electrons \(\hat{f}^{(e)}_i\) and nuclei + MM point charges \(f^{(n)}_i\)) and from every other induced dipole, which couples the dipoles through the dipole-field tensor \(\mathbf{T}^{(ij)}\):

\[ \mathbf{E}(\mathbf{R}_i) \;=\; f_i[\rho] \;-\; \sum_{j\neq i}\mathbf{T}^{(ij)}\mathbf{p}_j. \]

Collecting the fields into a \(3N_{\text{pol}}\)-component supervector \(\mathbf{f}\), the effective polarizability tensor of the whole MM region becomes

\[ \mathbf{A} \;=\; \big(\boldsymbol{\alpha}^{-1} + \mathbf{T}\big)^{-1}, \]

and the classical-form polarization energy \(U_{\text{pol}}^{\text{classical}} = -\tfrac{1}{2}\mathbf{f}^\top\mathbf{A}\mathbf{f}\) becomes, upon promoting \(\mathbf{f}\) to its quantum-mechanical operator counterpart \(\hat{\mathbf{f}} = \hat{\mathbf{f}}^{(e)} + \mathbf{f}^{(n)}\), the polarization Hamiltonian:

\[ \hat{H}_{\text{pol}} \;=\; -\tfrac{1}{2}\hat{\mathbf{f}}^{(e)\,\top}\mathbf{A}\,\hat{\mathbf{f}}^{(e)} \;-\; \mathbf{f}^{(n)\,\top}\mathbf{A}\,\hat{\mathbf{f}}^{(e)} \;-\; \tfrac{1}{2}\mathbf{f}^{(n)\,\top}\mathbf{A}\,\mathbf{f}^{(n)}. \]

Because the field \(\hat{\mathbf{f}}^{(e)}\) is linear in the electronic positions, \(\hat H_{\text{pol}}\) decomposes cleanly into a constant ("zero-electron") term, a one-electron term, and a two-electron term:

\[ \hat H_{\text{pol}} \;=\; \tfrac{1}{2}\sum_{a\neq b}\hat h^{(2)}(a,b) \;+\; \sum_a \hat h^{(1)}(a) \;+\; h^{(0)}, \]

with

\[ \begin{aligned} \hat h^{(2)}(a,b) &= -\sum_{i,j} \hat{f}^{(e)\,\top}_{ia}\,\mathbf{A}_{ij}\,\hat{f}^{(e)}_{jb}, \\[2pt] \hat h^{(1)}(a) &= \sum_{i,j,n} Q_n\,f^{(n)\,\top}_{in}\,\mathbf{A}_{ij}\,\hat{f}^{(e)}_{ja} \;-\; \tfrac{1}{2}\sum_{i,j} \hat{f}^{(e)\,\top}_{ia}\,\mathbf{A}_{ij}\,f^{(e)}_j, \\[2pt] h^{(0)} &= -\tfrac{1}{2}\sum_{m,n}\sum_{i,j} Q_m Q_n\,f^{(n)}_{im}\,\mathbf{A}_{ij}\,f^{(n)}_{jn}. \end{aligned} \]

Each piece adds to the corresponding piece of the bare QM Hamiltonian: \(h^{(0)}\) to nuclear repulsion, \(\hat h^{(1)}\) to the core Hamiltonian, and \(\hat h^{(2)}\) to the two-electron Coulomb/exchange operators (modulated by same_site_exact — see below). The "2e-Pol" in the name emphasizes the two-electron contribution that distinguishes IEDRF from mean-field flavors of polarizable embedding.

Damped polarization operator

The bare polarization integral \(\langle\mu|(\mathbf{r}-\mathbf{R}_i)/|\mathbf{r}-\mathbf{R}_i|^3|\nu\rangle\) diverges at \(r=R_i\) because the MM atoms enter as point-induced dipoles, when physically they should have a finite charge distribution. TeraChem follows the Thole short-range damping prescription, replacing the field by

\[ \mathbf{F}^{(e)}_{\mu\nu}(\mathbf{R}_i) \;=\; \Big\langle\mu\Big|\frac{\mathbf{r}-\mathbf{R}_i}{|\mathbf{r}-\mathbf{R}_i|^3}\, C\big(|\mathbf{r}-\mathbf{R}_i|\big)\Big|\nu\Big\rangle, \qquad C(r) \;=\; \big[1 - e^{-\alpha r^2}\big]^q. \]

The exponent \(\alpha\) is set by cutoff_exponent (units of bohr\(^{-2}\); default 4.0) and the power \(q\) is set by cutoff_power (default 2). \(\alpha\) must be positive; \(q\) must be \(\geq 2\). The same damping is applied to the fields from nuclei and MM point charges so that the short-range contributions cancel for a neutral molecule.

Same-site polarization integrals

When the same polarizable atom \(i\) appears on both sides of the two-electron polarization integral \(\langle\mu\nu||\hat h^{(2)}|\lambda\sigma\rangle\), the term reduces to a sum of one-electron "core-polarization" integrals \(\sum_{\alpha\beta}A_{i,\alpha\beta}I^{\alpha\beta}_{\mu\nu}(\mathbf{R}_i)\). TeraChem supports two ways of handling this same-site block:

same_site_exact yes (default): evaluate \(I^{\alpha\beta}_{\mu\nu}\) exactly.
same_site_exact no: insert a resolution-of-the-identity approximation in the AO basis. This roughly halves the integral count and is appropriate for large MM regions, but is less accurate.

Electron spill-out and MM-ECP

A QM electron sees only an attractive \({-1}/r\) tail from MM point charges and atomic-dipole induced fields, so a small fraction of the wave function leaks ("spills out") onto polarizable cation centers — an unphysical artifact that DRF aggravates relative to plain electrostatic embedding. The standard cure is to place a repulsive effective core potential (ECP) on each MM atom, so that the QM density sees a Pauli-like wall at short range. Enable this with the mm_ecp keyword, pointing it at a basis-set file whose ECP parameters cover the MM elements:

mm_ecp    mm_ecp_marefat2020

The bundled mm_ecp_marefat2020 parameters cover H through Ar and are taken from Marefat Khah et al., J. Chem. Theory Comput. 16, 1373 (2020). Only the ECP block of the basis file is read; the orbital basis on MM atoms is ignored. By default ECPs are removed from MM atoms that are bonded to the QM region (the "link atom" side); set ecps_on_mm_link_atoms yes to keep them.

Recommended pairing

direct_reaction_field should almost always be used with mm_ecp for production work. Without ECPs the QM density can polarize unphysically onto highly polarizable MM atoms, especially with diffuse basis sets.

Worked examples

RHF gradient in a polarizable solvent shell

Single-point energy + gradient on a QM water surrounded by five MM waters, each polarizable, at HF/aug-cc-pVDZ:

DRF_rhf_gradient.in

prmtop                   solvated.prmtop
qmindices                solvated.qmregion
coordinates              solvated.rst7

basis                    aug-cc-pvdz
method                   hf
charge                   0
spinmult                 1

direct_reaction_field    yes
same_site_exact          yes
cutoff_exponent          4.0

run                      gradient
precision                double
threall                  1.0e-15

CASSCF excited-state gradient with DRF

S\(_1\) analytic gradient with CASSCF(2,2) and DRF embedding:

DRF_casscf_S1_gradient.in

prmtop                   solvated.prmtop
qmindices                solvated.qmregion
coordinates              solvated.rst7

basis                    sto-3g
method                   hf
charge                   0
spinmult                 1

casscf                   yes
closed                   4
active                   2
cassinglets              2
castarget                1            # S1

cphfiter                 1000
cphftol                  1.0e-14

direct_reaction_field    yes
same_site_exact          yes
cutoff_exponent          4.0

run                      gradient
precision                double
threall                  1.0e-15
convthre                 1.0e-12

Saving induced dipoles

DRF_with_dipoles.in

prmtop                   LiF_C60.prmtop
qmindices                LiF_C60.qmregion
coordinates              LiF_C60.rst7

basis                    sto-3g
method                   hf
charge                   0
spinmult                 1

direct_reaction_field    yes
same_site_exact          yes
cutoff_exponent          4.0
save_induced_dipoles     yes

run                      energy

This writes one file per electronic state — scr/induced_dipoles_I.dat, where I = 0 is the ground state — with the per-site induced dipole in Debye. Useful for visualizing the polarization response or diagnosing unphysical spill-out.

Non-adiabatic coupling with CASCI

DRF couplings between two CASCI states, used in conical-intersection optimization:

DRF_nac.in

prmtop                   solvated.prmtop
qmindices                solvated.qmregion
coordinates              solvated.rst7

basis                    sto-3g
method                   hf
casci                    yes
closed                   4
active                   2
cassinglets              2

nacstate1                0
nacstate2                1
tvarnac                  yes

direct_reaction_field    yes
same_site_exact          yes
cutoff_exponent          4.0

cphfiter                 1000
cphftol                  1.0e-14
run                      coupling

Many more worked inputs (RPA/TDDFT, CASCI, link atoms, diffuse basis sets, diamond-with-link-atoms DFT) are bundled under tests/tests/DRF2ePol/ in the TeraChem source.

Keyword reference

Keyword	Type	Default	Description
`direct_reaction_field`	bool	`no`	Master switch. Enables polarizable embedding with the IEDRF Hamiltonian. Requires the OpenMM driver (`prmtop` + `qmindices`) and a `DRF2EPOL` build.
`same_site_exact`	bool	`yes`	Evaluate same-site two-electron polarization integrals exactly. Set to `no` for resolution-of-identity (faster, less accurate).
`cutoff_exponent`	float	`4.0`	Damping exponent \(\alpha\) (bohr\(^{-2}\), must be > 0) in \(C(r) = [1-\exp(-\alpha r^2)]^q\).
`cutoff_power`	int	`2`	Damping power \(q\) (must be \(\geq 2\)). On the GPU only the compiled-in values of \(q\) are accepted; CPU evaluation supports arbitrary \(q\).
`cpp_integrals`	string	`gpu`	Where to evaluate polarization integrals: `gpu` (default, fast) or `cpu` (slower, supports arbitrary `cutoff_power`).
`save_induced_dipoles`	bool	`no`	Write the induced dipole on each polarizable MM site (Debye) for every electronic state to `scr/induced_dipoles_I.dat`.
`mm_ecp`	string	not set	Filename of the basis file whose ECP block is read and placed on the MM atoms. The bundled set `mm_ecp_marefat2020` covers H–Ar. Strongly recommended for DRF.
`ecps_on_mm_link_atoms`	bool	`no`	Keep MM-ECPs on MM atoms that are bonded to the QM region. Off by default because the ECPs would raise the orbital energy on the bonded QM atom and the H capping atom.

B. T. Thole and P. T. van Duijnen, Theor. Chim. Acta 55, 307 (1980). doi:10.1007/BF00549429. ↩
X. Liu, A. Humeniuk, W. J. Glover, J. Chem. Theory Comput. 18, 6826 (2022). doi:10.1021/acs.jctc.2c00662. Introduces the integral-exact DRF reformulation (QM/MM-IEDRF). ↩
A. Humeniuk and W. J. Glover, J. Chem. Theory Comput. 20, 2111 (2024). doi:10.1021/acs.jctc.3c01018. Multistate IEDRF with analytic gradients and conical-intersection optimization; this is the implementation in TeraChem and the source of the equations summarized above (Sec. 2.1–2.2). ↩