8. Polarizable Environments: Thole and Charge-Response Models

In a QM/MM job the classical environment of the QM region is divided into regions (staticregion, polarregion, ewaldregion). Static regions carry permanent multipoles only. Polar regions also respond to electric fields, and this chapter describes how. A polar region can contain molecules of two kinds, which can be mixed in one region:

• Thole molecules: inducible point dipoles on the atoms, coupled to each other inside and between molecules with Thole damping [Thole:1981]. This has been XTP’s model so far.

• Charge-response molecules: induced atomic charges from a charge-response matrix \(\Pi\) of the molecule, plus residual induced dipoles, both smeared over Gaussians, with no couplings inside the molecule [Morita:1997] [DAvino:2014].

Both models give the same molecular polarizability tensor when they are parametrized to it. They differ in where on the molecule the response sits, which matters close to the QM region: at van der Waals contacts the QM charge density probes the response atom by atom. In conjugated molecules most of the in-plane polarizability is intramolecular charge flow [Tsiper:2001], which induced charges describe directly and atomic dipoles only effectively.

The same model, with the same smearing, is used in the ground-state QM/MM loop and in the reaction field of environment-screened GW-BSE (Embedded GW-BSE: Screening by a Polarizable Environment).

8.1. Common ground

8.1.1. Sites and moments

Every polar site carries the permanent multipoles of its .mps file (up to quadrupoles) and a polarizability tensor \(\alpha_i\) (the P line, in ų, isotropic or as the six components xx xy xz yy yz zz). The site responds with an induced dipole \(\boldsymbol{\mu}_i\); a charge-response site in addition with an induced charge \(\delta q_i\).

8.1.2. Fields and potentials

During an outer QM/MM iteration every other region adds its field (and, at charge-response sites, its potential) to the sites of a polar region:

• permanent multipoles of static regions and of the other molecules of the same or other polar regions, as point multipoles;

• the induced moments of other polar regions;

• the QM region: the field of its electron density and nuclei;

• the periodic background of an ewaldregion: its field, and on charge-response sites also its potential (see Periodic background).

Two accumulators per site keep apart what is counted in the polar region’s energy (sources in outer regions) and what is counted in the other region’s energy (for example the QM region, whose Hamiltonian contains the induced moments).

8.1.3. Energies

A polar region reports

\[E_\text{static} = E_\text{static-static} + E_\text{static,ext}, \qquad E_\text{polar} = E_\text{indu-indu} + E_\text{indu-stat} + E_\text{internal} + E_\text{polar,ext},\]

with the interactions among induced moments of different sites (\(E_\text{indu-indu}\)), of induced moments with the permanent multipoles of other molecules of the region (\(E_\text{indu-stat}\)), the internal (polarization) energy \(E_\text{internal} = \tfrac12 \sum_i \boldsymbol{\mu}_i^T \alpha_i^{-1} \boldsymbol{\mu}_i\) (plus the charge term below for charge-response molecules), and the interaction with the fields of outer regions. At self-consistency \(E_\text{polar}\) is the polarization energy \(\tfrac12 x^T b\) of the region in the fields it is counted with.

8.2. The Thole model

Each site responds to the total field, external plus that of all other induced dipoles,

\[\boldsymbol{\mu}_i = \alpha_i \Big( \mathbf{E}_i^\text{ext} - \sum_{j \neq i} T_{ij}\, \boldsymbol{\mu}_j \Big) \quad\Leftrightarrow\quad (\alpha^{-1} + T)\, \boldsymbol{\mu} = \mathbf{E}^\text{ext},\]

where the sum runs over all other sites, inside the molecule as well. The dipole tensor is damped at short range,

\[T_{ij} = \frac{\lambda_3}{r^3}\,\mathbb{1} - \frac{3\lambda_5}{r^5}\,\mathbf{r}\mathbf{r}^T, \qquad \lambda_3 = 1 - e^{-a u^3}, \quad \lambda_5 = 1 - (1 + a u^3)\, e^{-a u^3}, \quad a u^3 = \frac{a\, r^3}{\sqrt{\alpha_i \alpha_j}},\]

with \(a\) the option exp_damp (default 0.39) and \(\alpha_i\) the largest eigenvalue of the site’s polarizability tensor. Without the damping, two point dipoles closer than \((4 \alpha_i \alpha_j)^{1/6}\) polarize each other without bound (the polarization catastrophe).

Because the intramolecular couplings are part of the model, the site polarizabilities are not the molecule’s polarizability divided among its atoms: xtp_tools -e molpol adjusts them until the coupled molecule reproduces a target polarizability tensor, usually from a DFT calculation. The target is read from an ORCA log (qmpackage=orca) or computed by xtpdft itself: with dftpackage.xtpdft.polarizability=true the static polarizability is obtained analytically after the SCF (coupled-perturbed Kohn–Sham; closed shells, LDA/GGA/global hybrids) and stored in the orb file, which molpol reads with qmpackage=xtp logfile=system.orb:

xtp_tools -e dftgwbse -c job_name=system tasks=input,dft,parse \
  dftpackage.xtpdft.polarizability=true
xtp_tools -e molpol -c job_name=system qmpackage=xtp logfile=system.orb

Solver. \((\alpha^{-1} + T)\,\boldsymbol{\mu} = \mathbf{E}\) is solved by preconditioned conjugate gradients (tolerance_dipole, max_iter), starting from the solution of the previous outer iteration; a region of one molecule is solved directly.

Coupling with the QM region. The sites see the field of the QM electrons integrated on the DFT grid at the site position and that of the nuclei, both as points. The QM region sees the induced dipoles as point dipoles in its Hamiltonian. In the reaction field of embedded GW-BSE the sites are instead smeared over Gaussians of width site_width \(\cdot\,\alpha^{1/3}\) to keep \(1 + R\) positive definite (Embedded GW-BSE: Screening by a Polarizable Environment).

8.3. The charge-response model

8.3.1. Response of one molecule

A charge-response molecule \(I\) responds to the external potential \(\boldsymbol{\phi}_I\) at its atoms with induced atomic charges and to the field with residual induced dipoles,

\[\delta\mathbf{q}_I = -\Pi_I\, \boldsymbol{\phi}_I, \qquad \boldsymbol{\mu}_i = \alpha_i^\text{res}\, \mathbf{E}_i .\]

\(\Pi_I\) (\(N \times N\), in e²/Hartree) is the charge-response matrix of the molecule type: symmetric, positive semidefinite, and charge conserving (each row sums to zero, so a uniform potential induces nothing and every molecule stays neutral). \(\alpha_i^\text{res}\) carries what induced charges cannot: the out-of-plane response of a planar molecule and a small in-plane residual. In a uniform field \(\mathbf{F}\) the potential is \(\phi_k = -\mathbf{F}\cdot\mathbf{R}_k\), so the molecular polarizability is

\[\alpha_\text{mol} = X^T \Pi\, X + \sum_i \alpha_i^\text{res},\]

with \(X\) the atomic positions. The residual polarizabilities are fitted so that this equals the DFT polarizability of the molecule exactly.

\(\Pi\) is a linear-response property of the whole molecule and already contains its own Coulomb screening. The induced charges and dipoles of one molecule therefore do not interact with each other; they interact only with other molecules and with the QM region.

8.3.2. Smearing

The induced charge and the residual dipole of a charge-response site are Gaussians: the charge density \(\delta q\, g(\mathbf{r})\) with \(g(\mathbf{r}) = (\pi s^2)^{-3/2} e^{-r^2/s^2}\), and the dipole as \(-\boldsymbol{\mu}\cdot\nabla g\). The width \(s\) is the option cr_width (default 0.7 Å). Two Gaussians of widths \(s_i\) and \(s_j\) interact through \(\operatorname{erf}(r/S)/r\) with \(S = \sqrt{s_i^2 + s_j^2}\) and its derivatives,

\[B_0 = \frac{\operatorname{erf}(r/S)}{r}, \quad B_1 = \frac{B_0 - \frac{2}{\sqrt\pi\, S} e^{-r^2/S^2}}{r^2}, \quad B_2 = \frac{3 B_1 - \frac{4}{\sqrt\pi\, S^3} e^{-r^2/S^2}}{r^2},\]

so that, with \(\mathbf{r}\) from the source to the receiving site, a charge \(q\) produces the potential \(q B_0\) and the field \(q B_1 \mathbf{r}\), and a dipole the potential \(B_1\, \boldsymbol{\mu}\cdot\mathbf{r}\) and the field \((B_2\, \mathbf{r}\mathbf{r}^T - B_1 \mathbb{1})\, \boldsymbol{\mu}\). For \(S = 0\) these are the bare Coulomb kernels. Thole sites are points (\(s = 0\)).

• Charge–charge and charge–dipole couplings always use these kernels.

• Dipole–dipole couplings between molecules, when at least one site is a charge-response site, use them as well (cr_dipole_damping gaussian, the default) or the Thole damping above (thole). Between two Thole sites the Thole model is unchanged.

• Permanent multipoles act on charge-response sites as points. The Gaussian correction would be below \(3\cdot 10^{-4}\) relative at the shortest intermolecular distances (1.7 Å for \(s = 0.7\) Å).

The smearing is what makes the induced charges usable at contact distances: a point charge right next to a hydrogen bond over-responds in the same way a point dipole does. The QM/MM loop and the reaction field of embedded GW-BSE use the same width, so the environment responds to the QM density in the same way in both.

8.3.3. Linear system

With the unknowns \(x = (\delta\mathbf{q}, \boldsymbol{\mu})\) (induced charges of all charge-response sites, induced dipoles of all sites), the sources \(c = (-\boldsymbol{\phi}, \mathbf{E})\) and \(M = \operatorname{blockdiag}(\Pi_I, \alpha_i)\), the response of the region is

\[x = M\,(c - T x) \quad\Leftrightarrow\quad (\mathbb{1} + M^{1/2}\, T\, M^{1/2})\, y = M^{1/2} c, \qquad x = M^{1/2} y,\]

where \(T\) contains the couplings above: none inside a charge-response molecule, the Thole couplings inside a Thole molecule, and all couplings between molecules. \(T\) is built so that \((Tx)_q\) is the potential and \((Tx)_\mu\) minus the field of the induced moments. The symmetric operator is positive definite for a stable environment and is solved by conjugate gradients (tolerance_dipole, max_iter). The null space of \(\Pi\) (charge conservation) is handled by \(M^{1/2}\). A non-positive curvature during the iterations stops the job with a message (polarization catastrophe). Regions without charge-response molecules use the Thole solver unchanged.

The internal energy of the induced charges, \(\tfrac12\, \delta\mathbf{q}^T \Pi^{+} \delta\mathbf{q}\) with the pseudo-inverse of \(\Pi\), is added to \(E_\text{internal}\).

8.3.4. Coupling with the QM region

The induced charge and dipole of a charge-response site enter the QM Hamiltonian as the Gaussian charge and dipole above. The matrix elements are three-centre Coulomb integrals of an AO pair with an s- and a p-type primitive Gaussian at the site, normalized to a unit charge and a unit dipole; their interaction with the QM nuclei uses the kernels \(B_0\) and \(B_1\). The potential and field of the QM density at the site, seen by the same Gaussian, are traces of the density matrix with the same integrals, and the nuclei again enter through \(B_0\) and \(B_1\). The coupling is therefore exactly symmetric, without grid errors. Only the xtp DFT package supports this; with other packages a job with charge-response molecules stops at the first outer iteration.

8.3.5. Reaction field of embedded GW-BSE

With charge-response molecules among the explicit polar regions, the reaction-field kernel over the auxiliary basis becomes

\[B = -(M^{1/2} C)^T\, (\mathbb{1} + M^{1/2} T M^{1/2})^{-1}\, (M^{1/2} C), \qquad C = (-\Phi^T, F^T),\]

with \(\Phi\) the potential of each auxiliary function at the charge-response sites (Gaussian-averaged with width cr_width) and \(F\) its field at all sites (Thole sites smeared with site_width \(\cdot\,\alpha^{1/3}\) as before). For Thole molecules alone it reduces to \(-F (\alpha^{-1} + T)^{-1} F^T\). It is formed with a Cholesky factor of the dense operator, as the Thole kernel. Charge-response molecules cannot be part of a dielectric shell region (shell_regions), which responds with \(\alpha/\varepsilon\) per site.

8.4. Parametrization of charge-response molecules

\(\Pi\) is the static linear response of the molecule, obtained once per molecule type from a DFT calculation, like the target polarizability of molpol: natively with xtp_tools -e crfit (xtpdft), or with the python tool xtp_cr_fit (PySCF); see Tools below. The procedure:

1. Perturbations. An external potential localized on atom \(k\) with Hirshfeld weights [Hirshfeld:1977], \(v_k(\mathbf{r}) = w_k(\mathbf{r})\) with \(\sum_k w_k = 1\) (free spherical atoms in the same basis), plus three uniform fields.

2. Response. The first-order density matrices of all \(N + 3\) perturbations from analytic coupled-perturbed Kohn–Sham (crfit, xtp_cr_fit) or central finite differences of SCFs (crfit method=finite_differences); with diffuse functions in the basis (for example def2-SVPD or def2-TZVPD) and the functional of the QM region.

3. Induced charges. Measured with the same Hirshfeld weights: \(\Pi_{jk} = -\langle w_j | \chi | w_k \rangle\), which is symmetric, positive semidefinite and charge conserving by construction. (ESP-fitted charges, as in the original charge-response kernel [Morita:1997], give a non-symmetric, more ill-conditioned matrix with amplified short-wavelength modes; they were less accurate in the tests below.)

4. Residual dipoles. \(\alpha^\text{res} = \alpha_\text{DFT} - X^T \Pi X\) from the field response, distributed over the atoms as one tensor scaled by the trace of the Hirshfeld atomic dipole polarizabilities. \(\alpha^\text{res}\) must be positive definite.

8.4.1. Files

• .crm: the charge-response matrix in the atom order of the molecule type’s .mps file. Lines starting with ! are comments; then the number of atoms \(N\), \(N\) lines element x y z (Å, the geometry of the fit, for information), and \(N\) lines of \(N\) numbers (e²/Hartree). On reading, XTP checks that it is symmetric, charge conserving and positive semidefinite to \(10^{-6}\) relative, and removes the rounding of the file.

• .mps: the permanent multipoles as for any polar molecule, with the residual polarizabilities \(\alpha_i^\text{res}\) on the P lines (six components xx xy xz yy yz zz in ų), in the frame of the template geometry.

• map file: next to <multipoles_n> of the segment,

  <multipoles_n>MP_FILES/water_cr_n.mps</multipoles_n>
  <chargeresponse_n>MP_FILES/water_cr_n.crm</chargeresponse_n>
  

  (and likewise for other states). The matrix is reordered to the order of mpoles in the fragments. Every atom of the molecule must be mapped: a part of \(\Pi\) does not conserve charge. Segments mapped from an explicit .mps file name in the job file get no charge response.

• polar region options: cr_width and cr_dipole_damping (above); they apply to the charge-response molecules of the region.

Residual polarizabilities are in general anisotropic. Use map2md 0: with map2md 1 the sites are placed on the MD atoms without rotating site tensors.

Checkpoints store \(\Pi\), the induced charges and the options with the region; checkpoints without them read as Thole molecules.

8.4.2. Tools

xtp_tools -e crfit is part of XTP and needs nothing else. It reads job_name.xyz and the dftpackage options (xtpdft; closed-shell molecules; LDA, GGA or global hybrid functionals), runs the ground state (converged to at least scf_error, default \(10^{-9}\)) and obtains the first-order density matrices of all perturbations at once by coupled-perturbed Kohn–Sham (method=analytic, the default; see DFTEngine::StaticResponse), solved to a relative residual tolerance (default \(10^{-7}\)). With method=finite_differences it runs instead two SCFs per perturbation at \(\pm\) step (default \(10^{-3}\) Hartree, resp. Hartree/bohr), warm-started from the ground state, and takes \(\Delta P = (P_+ - P_-)/(2s)\) with \(s\) = step. If the geometry has an inversion centre (inversion, default on), only one atom of each pair is perturbed. The free-atom references of the Hirshfeld weights are those of xtpdft’s CDFT (spherically averaged, Hund’s-rule spin, the functional of the calculation). It writes output.crm, output.json (the keys of xtp_cr_fit) and, with mpsfile (a template .mps with the permanent multipoles), output.mps with the residual polarizabilities on the P lines, rotated into the template frame as xtp_cr_mps does. Example:

xtp_tools -e crfit -c job_name=water output=water_cr_n \
  mpsfile=water_n.mps dftpackage.basisset=def2-svpd \
  dftpackage.auxbasisset=aux-def2-svpd \
  dftpackage.functional=XC_HYB_GGA_XC_PBEH

Both methods agree: for pentacene (PBE0/def2-SVPD) \(\Pi\) to \(10^{-6}\); the polarizability differs by \(2 \times 10^{-4}\), the hyperpolarizability contribution of the finite field step (the analytic response agrees with SCFs at step \(10^{-4}\) to \(10^{-7}\)). Compared with xtp_cr_fit at the same functional and basis, the polarizability agrees to a few \(10^{-4}\), and the elements of \(\Pi\) are 3 % (water) to 5–7 % (pentacene) larger; the difference in \(\Pi\) is the choice of free-atom references (xtp_cr_fit uses spherical Hartree–Fock atoms). Pentacene on two cores: 23 min analytic, 44 min with finite differences, 2 h with xtp_cr_fit.

The python tools are installed with XTP when their modules are found at configure time: xtp_cr_fit needs PySCF and numpy, xtp_cr_mps numpy only.

• xtp_cr_fit -x mol.xyz -o mol runs steps 1–4 for one molecule (restricted Kohn–Sham, i.e. closed shell; --charge for closed-shell ions) and writes mol_cr.json (\(\Pi\), \(\alpha_\text{DFT}\), the charge-flow and residual parts, the residual tensors per atom, the settings) and mol.crm. Defaults: --xc pbe0, --basis def2-svpd, density fitting with --auxbasis def2-universal-jkfit; --grid sets the DFT grid and --response-grid (default 1) the coarser grid of the response kernel. An inversion centre of the geometry, if present, is used to halve the number of responses (--no-inversion switches this off). The fit stops with an error if \(\Pi\) is not positive semidefinite or \(\alpha^\text{res}\) not positive definite. Use the same functional as the QM region and a basis with diffuse functions; the geometry should be that of the .mps template (e.g. the molpol input).

• xtp_cr_mps -m template.mps -j mol_cr.json -o mol_cr_n writes the pair of files for the map file: mol_cr_n.mps keeps the permanent multipoles of the template and replaces its P lines by the residual polarizabilities, mol_cr_n.crm holds \(\Pi\). The template must list the same atoms in the same order as the fit; if it is in another frame, the residual tensors are rotated into it (a warning is printed if the geometries differ by more than 0.05 Å rms after alignment).

For water in def2-SVPD either fit takes seconds. The permanent multipoles of other states come from the template .mps files of those states as usual.

8.5. Accuracy and choice of model

The induction energy \(P = \tfrac12 \langle \rho | \chi_\text{env} | \rho \rangle\) of one electron in a frontier orbital of a molecule A next to a molecule B (B the environment), compared with B’s DFT finite-field response (PBE0; charge-response parameters at the same level):

A, B, contact	DFT (eV)	Thole	charge response, \(s = 0.7\) Å
pentacene pairs (bulk, 2.6–2.8 Å)	−0.10 to −0.17	+8 to +9 %	−4 to −10 %
pyridine lone pair + H-bonded water (1.76 Å)	−0.318	+84 %	+5 %
pyridine + 6 waters (coupled cluster of waters)	−0.449	+62 %	+5 %
phenol + 6 waters	−0.178	+8 %	−3 %

For pentacene the Thole model was used with its molpol polarizabilities; scaled to the DFT polarizability it is +15 to +17 %. Isotropic Thole sites cannot follow the anisotropy of an acene (eigenvalues 21/41/68 ų vs. 18/37/91 ų in DFT), while the charge-response model reproduces the DFT tensor by construction. The width matters at hydrogen bonds (point sites: +29 % for the pyridine pair) and hardly at all between π-stacked acenes.

The charge-response model is the better choice for the molecules in direct contact with the QM region when they are conjugated or hydrogen-bonded to it. It costs one DFT linear-response calculation per molecule type; the QM/MM loop itself is about as expensive as with Thole sites. Thole molecules remain appropriate further out and can be mixed into the same region.

8.6. Periodic background

An ewaldregion is a frozen periodic background, converged beforehand by the ewaldbackground calculator with Thole molecules; the molecules of the job’s regions are cut out of it. A polar region in it may hold charge-response molecules: the background acts on them like the static environment of the other regions, with the field on every site and the potential on the charge-response sites, both of its permanent multipoles and its induced dipoles (real-space, reciprocal-space and shape terms, minus the cut-out copies). The background itself stays a Thole model and does not respond to the job.

Two details differ from the coupling between polar regions:

• the field of the background’s induced dipoles is Thole-damped (the Ewald real-space convention), also on charge-response sites with cr_dipole_damping gaussian; and

• the potential is that at the site as a point, not averaged over the Gaussian of width cr_width.

Both matter only for background molecules at contact distance from the polar region; with \(s = 0.7\) Å the Gaussian average of a point source 3 Å away differs from the point value by less than \(10^{-6}\). The potential carries the additive constant of the Ewald sum (the omitted \(k = 0\) term); it cancels, since the induced charges of each molecule sum to zero. The energy terms are those of Thole molecules in the background, with \(\delta q_i \phi_i\) added to the external polar energy.

8.7. Limitations

• Charge-response molecules need the xtp DFT package.

• They cannot be part of a dielectric shell region of embedded GW-BSE.

• An ewaldregion background is always made of Thole molecules (see Periodic background).

• Whole molecules only (see Files), and map2md 0 for anisotropic residual polarizabilities.

• The dense reaction-field kernel of embedded GW-BSE has the dimension \(N_\text{CR} + 3 N_\text{sites}\), slightly more than the \(3 N_\text{sites}\) of the Thole model.