6. Polarizable Ewald Embedding

A charge carrier in a molecular solid does not sit in vacuum. Its energy is shifted by the electrostatic potential of every other molecule in the sample, and by the polarization those molecules undergo in response to the carrier’s own field. Both contributions are long ranged: the electrostatic one falls off only as \(1/r\) for a charged carrier, and the polarization response is itself driven by that slowly decaying field. A calculation that simply truncates the environment at some radius does not converge to the right answer, it converges to an answer that depends on where the truncation was placed.

The machinery described here computes those shifts for a periodic sample, by combining an Ewald lattice sum over a self-consistently polarized background with an explicit, quantum-mechanical treatment of the molecule carrying the charge. It is used in two ways:

• as a purely classical calculation, in which every molecule including the central one is represented by distributed multipoles, and

• as a three-region QM/MM calculation, in which the central molecule is described by DFT, a shell of neighbours is treated as explicitly polarizable point multipoles, and everything beyond that is the periodic Ewald background.

The two share all of their machinery below the level of the central molecule’s description, which is deliberate: it means the classical calculation is a reference for the QM/MM one, and a disagreement between them is a statement about the description of that molecule rather than about the environment.

6.1. The Ewald decomposition

The quantity to be summed is the electrostatic interaction between a set of multipoles in a periodic cell and all of their periodic images,

(6.1)\[E = \frac{1}{2}\sum_{\mathbf{T}}\ {\sum_{i,j}}' Q_i \, \hat{T}(\mathbf{r}_i - \mathbf{r}_j - \mathbf{T}) \, Q_j ,\]

where \(\mathbf{T}\) runs over lattice translations, \(Q_i\) denotes the multipole moments of site \(i\), \(\hat{T}\) is the interaction tensor, and the prime excludes \(i = j\) in the \(\mathbf{T} = 0\) cell. This sum is only conditionally convergent: its value depends on the order in which the terms are taken, which is another way of saying it depends on the shape of the macroscopic sample and on the boundary condition applied at its surface.

The Ewald construction splits the \(1/r\) kernel with the error function,

(6.2)\[\frac{1}{r} = \underbrace{\frac{\mathrm{erfc}(\alpha r)}{r}}_{\text{short ranged}} + \underbrace{\frac{\mathrm{erf}(\alpha r)}{r}}_{\text{smooth}} ,\]

and evaluates the two pieces in the spaces where each converges quickly. The screened part is summed directly over neighbouring images in real space; the smooth part is summed in reciprocal space, where its Fourier transform decays as \(\exp(-k^2/4\alpha^2)\). The result is assembled from four contributions:

Real space

The erfc-screened interaction, summed over all pairs and lattice translations within a distance cutoff. Implemented by EwaldRealSpaceSum.

Reciprocal space

The erf-screened interaction, as a sum over reciprocal lattice vectors weighted by \(\exp(-k^2/4\alpha^2)/k^2\). Implemented by EwaldReciprocalSpaceSum, which builds a structure factor once per call and replays it against every target. The \(k = 0\) term is omitted; see The k = 0 term and the potential’s gauge.

Shape (surface) term

The uniform depolarizing field that encodes the boundary condition of Eq.6.1, expressed through the total dipole moment \(\mathbf{M} = \sum_j (q_j \mathbf{r}_j + \boldsymbol{\mu}_j)\) of the cell [DeLeeuw:1980]. Under vacuum boundary conditions,

\[\mathbf{E} = -\frac{4\pi}{3V}\mathbf{M} \quad\text{(cube or sphere)}, \qquad \mathbf{E} = -\frac{4\pi}{V}M_z\,\hat{\mathbf{z}} \quad\text{(slab)} .\]

Implemented by EwaldShapeCorrection. This is not a pairwise sum: it is a mean-field property of the whole sample’s surface polarization acting back on itself, so it includes the target’s own segment.

Self-interaction removal

The reciprocal sum, being a sum over the full periodic density, includes each site’s interaction with itself. For a point charge this term vanishes by symmetry. For a static dipole it does not: a dipole’s own erf-screened field at its own position is finite and equals \(\tfrac{4}{3}\alpha^3/\sqrt{\pi}\). Omitting it makes the total depend on \(\alpha\), and this was a real defect in an earlier version of this code.

6.1.1. The \(\alpha\) invariant

The splitting parameter \(\alpha\) in Eq.6.2 is arbitrary. It controls how work is divided between the real- and reciprocal-space sums, and nothing else. The total must therefore be independent of it, and any dependence is a bug rather than a tolerance.

This is the single most useful diagnostic in the whole implementation, because it is sensitive to a missing term rather than to a wrong factor. A term that is absent from one sum but present in the other will leave a residue that grows or shrinks as the split moves, even when every individual number looks plausible. Both defects mentioned above — the dipole self-term, and an omitted erf correction — were found this way, and the current implementation is flat to the tenth digit across a factor of two in \(\alpha\) while the individual channels move by factors of twenty to sixty.

Anyone changing this code should run an \(\alpha\) scan before believing a result. Convergence parameters must be held fixed while \(\alpha\) varies, not scaled with it: a truncation error that tracks \(\alpha\) imitates exactly the failure the test is looking for.

6.2. Segments, fragments and charge states

The classical representation of a molecule is a PolarSegment: a list of PolarSite objects, each carrying a position, a set of permanent multipole moments up to the requested rank, a polarizability tensor, and an induced dipole. Segments are built from the mapping file by the polar mapper, which associates each coarse-grained segment with a .mps file per charge state.

Multipoles are read from .mps files, in the format used by GDMA and related tools. Two conventions in that format are worth stating explicitly because they are easy to get wrong:

• the dipole line is ordered z x y, not x y z;

• multipoles are in \(e\,a_0\) and polarizabilities in \(\mathrm{\AA}^3\), while the code works internally in atomic units throughout (bohr, Hartree).

Each segment is registered for the charge states it is needed in — Neutral, Electron, Hole — and these are held in an EwaldRegistry, keyed by segment id and state. The registry is the single source of truth for the periodic density: it replaces the legacy design in which several containers shared raw pointers to the same segments and their charge state depended on which container had most recently been asked.

A site’s permanent and induced moments are kept separate throughout, and they enter the energy through different channels. This matters for the bookkeeping described in Energy accounting, and for the Thole damping described next, which applies to induced-dipole interactions and not to permanent ones.

6.3. The polarizable background

Before any job is run, the periodic environment must be polarized self-consistently. This is the job of the ewaldbackground calculator, which is run once per frame:

xtp_run -e ewaldbackground -o ewaldbackground.xml -f state.hdf5

Each site \(i\) carries an induced dipole \(\boldsymbol{\mu}_i\) responding to the total field at its position,

(6.3)\[\boldsymbol{\mu}_i = \boldsymbol{\alpha}_i \Big( \mathbf{F}^{\text{perm}}_i + \sum_{j \neq i} \mathbf{T}^{\text{Thole}}_{ij}\,\boldsymbol{\mu}_j \Big) ,\]

which is a linear system in the induced dipoles, solved to self-consistency. The field \(\mathbf{F}^{\text{perm}}\) is the full periodic Ewald field of the permanent multipoles; the coupling \(\mathbf{T}^{\text{Thole}}\) between induced dipoles is damped at short range by the exponential Thole model [Thole:1981], which prevents the polarization catastrophe that an undamped point-dipole model suffers when two polarizable sites approach each other.

The damping factors depend on both sites’ polarizabilities through \(u^3 = a\,r^3 \sqrt{\alpha_i \alpha_j}^{-1}\), so a site with no polarizability is read as complete overlap and damped maximally. That is the correct limit for two real sites and the wrong one for a field evaluation point; see Coupling to the QM region.

Two solvers are available. The default is a preconditioned conjugate gradient; a Jacobi over-relaxation scheme reproducing the legacy SOR iteration is retained for direct comparison. The PCG solver additionally reports a Lanczos estimate of the smallest eigenvalue of the interaction operator, which is a direct algebraic certificate of whether the system is positive definite rather than an inference from residual behaviour.

The result is written to a checkpoint, by default ewaldbackground.hdf5, containing

• the EwaldRegistry — every segment, in every registered charge state, with its converged induced dipoles, and

• an ewald_parameters group holding \(\alpha\), \(k_\text{max}\), \(r_\text{min}\), the field tolerance, the Thole parameter \(a\), the screening factor, the shape, and the simulation box.

The parameters are part of the checkpoint on purpose. See Parameter inheritance.

6.4. The three-region QM/MM setup

A QM/MM job is built from three regions, which must be enumerated in this order:

id	type	role
0	qmregion	The molecule carrying the charge, described by DFT. Holds the segment whose site energy is wanted, in the charge state the job specifies.
1	polarregion	An explicit shell of neighbours, as polarizable point multipoles. These respond to the QM density and to the background, and they polarize each other.
2	ewaldregion	The periodic background, read from the ewaldbackground checkpoint. Frozen: nothing in the job repolarizes it.

An ewaldregion may not be region 0. It owns no segments carved out of the job’s topology — it represents the whole periodic cell — so the recentring that JobTopology applies has no segment of its own to key on.

Because the background is frozen, EwaldRegion is always Converged() and its Reset() is a no-op. Including it can never prevent the inter-region SCF loop from terminating. Its Interactwith* methods return zero, and this is the correct physics rather than a stub: no other region polarizes it. The influence in the other direction is implemented by that region’s own InteractwithEwaldRegion, following the receiver-pull convention used throughout Region::ApplyInfluenceOfOtherRegions.

A two-region qmregion + ewaldregion job is also valid and is useful as a diagnostic, since it removes the classical coupling entirely.

6.4.1. Foreground declaration and suppression

The segments treated explicitly — by the QM region and by the polar region together — are also present in the background, since the background is the whole cell. They must not be counted twice. Removing them is not simply a matter of dropping them from the sums:

• in real space the coincident copy of each explicit segment is suppressed, while its periodic images are kept. Deleting the images too would replace one carved-out cavity with a lattice of vacancies.

• in reciprocal space nothing is held out. A \(k\)-space sum runs over the whole periodic density and cannot have a hole cut in it. Instead, the erf-screened energy of exactly those coincident copies is subtracted afterwards, which removes precisely what the reciprocal sum put back.

The copies to be suppressed are the union over all regions that own segments, and no single region knows that union. JobTopology therefore calls EwaldRegion::RegisterForeground once, after building every region and before evaluating any of them. Without this, a QM/MM job would leave the QM segment’s neutral background copy sitting underneath the QM density — a ghost molecule in every sum, with no symptom to notice it by. The suppressed count is reported per job and checked against the expected total.

The erf correction uses the background’s own multipoles and induced dipoles, not the job’s charge state, because that is what the reciprocal sum actually placed there. Using the job’s state would remove something that was never added.

6.5. Energy accounting

Each pairwise contribution is computed exactly once, and which region books it follows a single rule: the region further inside owns the energy. A site accumulates the field from regions further out in V(), which carries energy, and from regions further in in V_noE(), which does not.

contribution	where it is booked
QM internal	the DFT total energy
QM ↔ polar	the DFT total, via the external multipole matrix (electrons) and ExternalRepulsion (nuclei). PolarRegion::InteractwithQMRegion returns zero by this convention.
QM ↔ background	the DFT total, via the potential on the integration grid (electrons) and a scalar added to \(E_0\) (nuclei)
polar internal (permanent)	E_static_static
polar ↔ background, permanent × permanent and permanent × induced	E_static_ext, returned by EwaldRegion::ApplyFieldTo
polar induced × background	E_polar_ext, as \(\sum \boldsymbol{\mu}\cdot\mathbf{V}\)
polar induction self-energy	E_polar_internal
background internal	not computed — see below

The background’s own internal energy is deliberately absent. It is a constant, independent of the job’s charge state, so it cancels exactly in any charge-state difference. Site energies are therefore correct; an absolute total energy is not, and should not be quoted as one.

One sign convention deserves attention. The Ewald code and the region framework store \(\mathbf{V}\) with opposite signs — the background solver builds its right-hand side as \(+\mathbf{V}\), the polar region as \(-(\mathbf{V} + \mathbf{V}_{\text{noE}})\). ApplyFieldTo negates its own contribution at the boundary so that the polar region’s induction runs in the right direction. Both conventions are internally consistent; only the boundary between them needs care.

6.5.1. Site energies

The quantity of interest is the difference between charge states, referenced to the same molecule in vacuum:

(6.4)\[\Delta E_{\text{env}} = \big[E(h) - E(n)\big]_{\text{embedded}} - \big[E(h) - E(n)\big]_{\text{vacuum}} .\]

Using the same geometries in both brackets makes the intramolecular relaxation cancel exactly, leaving the environment’s response to the charge at its own geometry. In practice this means running the same job file against a qmregion-only setup to obtain the vacuum reference.

6.6. Coupling to the QM region

The background reaches the QM region as a potential, evaluated on the same numerical integration grid the DFT calculation already uses:

\[E_{\text{QM-bg}} = -\int \rho(\mathbf{r})\,\phi(\mathbf{r})\,\mathrm{d}\mathbf{r} + \sum_A Z_A \phi(\mathbf{R}_A) .\]

The first term enters the one-electron Hamiltonian as a matrix built by Ewald_Potential::IntegrateEwald; the second is a scalar added to \(E_0\). The minus sign on the electronic term is the electron charge, and matches what AOMultipole::FillPotential does for the nuclei and for external multipoles.

That sign is worth flagging because getting it wrong is nearly invisible. For a neutral QM region the two terms differ only through the shape of the density, not its total charge, so they very nearly cancel; flipping one turns a near-cancellation into a near-doubling, and the result still looks like a plausible small number. On a neutral methane in a methane background, the difference is \(-13.5\) meV against a correct \(+0.13\) meV.

The grid must be the one DFTEngine uses, since the potential values have to land on exactly the quadrature points that form \(\langle\rho|\phi\rangle\). PrepareEwaldPotentialGrid therefore reads dftpackage.xtpdft.integration_grid, not the grid_for_potential option — the latter governs the opposite direction, the QM density’s influence outward on the classical regions. The coupling between the two grids is a constraint rather than a preference, and it costs nothing: the Ewald contribution is converged at the medium default to below \(10^{-10}\) Ha.

Because the background is frozen, the QM geometry is fixed for the job, and the QM density does not enter \(\phi\), the potential is evaluated once per job and reused across inter-region SCF iterations.

6.6.1. The probe and Thole damping

\(\phi\) at a point is obtained by placing a unit test charge there and calling the same validated energy routines the classical channels use, so that the result is in the same gauge by construction. That probe is a PolarSite, whose constructor unavoidably assigns a polarizability from the element table — which would otherwise feed a meaningless number into the Thole damping of the induced-dipole term.

The field-point path therefore requests an undamped induced-dipole interaction explicitly. This is not a preference: AOMultipole and DFTEngine::ExternalRepulsion already deliver induced dipoles to a QM density undamped, so damping here would give one job two different conventions for the same physical interaction depending on which route the dipole arrived through.

Note that this cannot be expressed by giving the probe no polarizability. The Thole factor would then be formed from \(u^3 = 0\), which the model reads as complete overlap and damps maximally — the opposite of the intended limit.

An analytic alternative to sampling the potential on a grid would hand the background’s moments to the DFT engine as operators, building the AO matrices directly. Most of it was written and is preserved in the git history; it was removed rather than kept, for two reasons.

It was blocked. A rank-1 (induced dipole) source needs operator-centre derivatives that stock libint2 builds do not provide, which is why the real-space part split every background dipole into a pair of point charges instead of using the dipole operator.

It would also not have been faster where it was wanted. Both routes evaluate the same Ewald sum; they differ only in how many places. The grid route does so once per quadrature point, which grows linearly with the size of the QM region; the analytic route once per shell pair, which grows quadratically. So it pays off for small QM regions and loses for large ones — the opposite of the case for adopting it to reach bigger molecules. The cost both routes share is the sum itself, which is addressed by balancing \(\alpha\) (see Choosing \alpha) rather than by changing how the result reaches the Hamiltonian.

6.6.2. The \(k = 0\) term and the potential’s gauge

The reciprocal sum omits the \(k = 0\) term, which corresponds to embedding the cell in a uniform neutralizing background. As a consequence \(\phi\) is defined only up to an additive constant \(\phi_0\), and a region of net charge \(q\) shifts by \(q\,\phi_0\).

For a neutral region this drops out identically: \(\phi_0\) multiplies \(Z_{\text{total}} - N_{\text{electrons}} = 0\) for the QM region, and \(\sum q = 0\) for each neutral classical segment. For a charged region it does not, and the resulting energy carries the usual charged-periodic-cell convention. This is a physical statement about the model, not a numerical artefact, and it is the same convention the classical channel uses — every term in PotentialAt is evaluated through the same energy routines, so the two halves of a QM/MM job share one gauge by construction rather than by coincidence.

6.7. Practical notes

6.7.1. Parameter inheritance

An ewaldregion takes every Ewald parameter from the background checkpoint, and none from the job’s own options. \(\alpha\), \(k_\text{max}\), \(r_\text{min}\), the field tolerance, the Thole parameter, the screening factor, the shape and the box are all read from the ewald_parameters group written by ewaldbackground. The job’s XML cannot override them, and does not offer the option.

This is deliberate. The background’s induced dipoles were converged with a particular splitting and particular cutoffs; the foreground’s interaction with that background is only meaningful if it is evaluated the same way. A job that re-specified \(\alpha\) would be adding a foreground computed under one decomposition to a background computed under another, and the \(\alpha\) invariance of The \alpha invariant would no longer hold — silently, since each half would look internally consistent.

The practical consequence is that changing any Ewald parameter means re-running ewaldbackground. The checkpoint is the unit of configuration, not the job file.

6.7.2. Choosing \(\alpha\)

Since the total is independent of \(\alpha\), the choice is purely one of cost, and ewaldbackground makes it for you. Both cutoffs are fixed multiples of \(\alpha\) – which is what keeps the accuracy \(\alpha\)-independent – so with \(N\) sites in a cell of volume \(V\) the work at one target is

\[N_\text{real} = \frac{N}{V}\,\frac{4\pi}{3} \left(\frac{s_r}{\alpha}\right)^{3}, \qquad N_k = \frac{4\pi}{3}\,\frac{V}{8\pi^{3}}\,(s_k\alpha)^{3} ,\]

falling as \(\alpha^{-3}\) and rising as \(\alpha^{3}\). Minimising the sum gives

\[\alpha_{\text{opt}} = \sqrt{2\pi}\; \frac{N^{1/6}}{V^{1/3}}\;\sqrt{\frac{s_r}{s_k}} .\]

The \(N^{1/6}\) is what makes Ewald \(O(N^{3/2})\): at the optimum the two halves cost the same, and that common value grows only as \(\sqrt{N}\), independent of \(V\).

The two multiples are tied to one accuracy. \(s_r\) is the screening_factor option, leaving a real-space tail of \(\mathrm{erfc}(s_r)\); the reciprocal cutoff is then derived to leave the same tail, \(s_k = 2\sqrt{-\ln \mathrm{erfc}(s_r)}\), which at the default \(s_r = 6\) gives \(k_\text{max} = 12.39\,\alpha\). So screening_factor alone sets how accurate both sums are, and \(\alpha\) follows from it and from the system.

Note

Before this derivation the defaults were \(\alpha = 3/L_\text{min}\) and \(k_\text{max} = 6\alpha\), which together switched the splitting off: in a cubic cell those give \(k_\text{max}/(2\pi/L) = 9/\pi\), a fixed \(\approx 98\) \(k\)-vectors regardless of system size, against a real-space cutoff of \(2L\). Per target on thiophene at experimental density that was 301,691 terms at 1000 segments and 1,508,063 at 5000 – of which 98 were reciprocal in both cases. Results were unaffected, since \(\alpha\) cancels; only the cost was.

Setting alpha or k_max explicitly overrides either half of the derivation. If \(k_\text{max}\) is fixed by hand, the balance changes: reciprocal work stops depending on \(\alpha\) and larger \(\alpha\) becomes monotonically cheaper, until \(k_\text{max}\) is no longer large enough to converge the \(\exp(-k^2/4\alpha^2)\) weight.

6.7.3. Cost

Evaluating the background potential over a DFT integration grid is the dominant cost of a QM/MM job. It scales as

\[N_\text{grid} \times \left( N_\text{real} + N_k \right) ,\]

the number of quadrature points times the work of one Ewald evaluation – the sources inside the real-space cutoff plus the \(k\)-vectors, not times them. The two inner terms are a sum, which is why balancing them against each other (see Choosing \alpha) is what governs the cost, and why leaving either one far larger than the other wastes almost all of the effort on one half of a split whose whole purpose is to avoid that.

All three numbers are logged before the evaluation starts, along with periodic progress so that a long run can be distinguished from a hung one. ewaldbackground additionally prints an Ewald balance line giving \(N_\text{real}\) and \(N_k\) side by side: at a derived \(\alpha\) they come out comparable, and a run where they differ by orders of magnitude is one where \(\alpha\) was set by hand or where the derivation was given a cell it does not suit. It is the first number to look at when a background solve is slower than expected.

The evaluation is parallelized over grid points and happens once per job – the background is frozen, so the potential is laid down on the grid exactly once and reused across every inter-region iteration.

6.7.4. Diagnostics

Two checks are worth running whenever this code is changed.

The \(\alpha\) scan described in The \alpha invariant: run the same job against backgrounds converged at several \(\alpha\), with all other convergence parameters held fixed, and confirm the total does not move while the individual channels do.

The with/without comparison: run a job with and without the ewaldregion, at otherwise identical settings, and difference the QM energies at the first inter-region iteration, before any induced dipoles exist. That difference is the QM–background interaction and nothing else. It is the sharpest instrument available on this code path, and it is what exposed the electron-charge sign described above — which the zeroed-background test, the unit tests, and the \(\alpha\) scan had all missed.