- Generated by
1.16.1
|
votca 2026-dev
|
Provides a means for comparing floating point numbers. More...
Namespaces | |
| namespace | csg |
| namespace | tools |
| namespace | xtp |
| Charge transport classes. | |
Classes | |
| struct | Log |
Typedefs | |
| using | Index = Eigen::Index |
Provides a means for comparing floating point numbers.
ewdbgpol-equivalent calculator, built entirely on the new Ewald machinery (EwaldRegistry, EwaldRealSpaceSum, EwaldReciprocalSpaceSum, EwaldShapeCorrection, EwaldPeriodicDipoleOperator) developed alongside this class, rather than the legacy xtp/ewald/ code – registered under a different Identify() ("ewaldbackground") so the two can coexist and be compared directly.
Calculates three electron repulsion integrals for GW and DFT.
Small wrapper for a segment id and the corresponding QMState or filename.
base class to derive regions from
defines a qm region and runs dft and gwbse calculations
Container to define fragments of QMmolecules, containing atomindices, no pointers to atoms, it also handles the parsing of strings etc.. Values should have own destructor.
Takes a list of atoms, and the corresponding density matrix and puts out a table of wavefunction partial charges.
defines a polar region and of interacting electrostatic and induction segments
Takes a list of atoms, and the corresponding density and overlap matrices and puts out a table of partial charges.
Class to set up the topology, e.g division of molecules into different regions for a specific job.
Small container for convergence checking, stores current and last value and gives the diff, or just the current value in the first iteration.
Takes a list of atoms, and creates CHELPG grid.
Parameters for Extended Hueckel calculation.
Iterative solvers for the periodic induced-dipole equation A*mu = F_perm, with A = EwaldPeriodicDipoleOperator.
Shape/surface term of a 3D Ewald sum: the uniform depolarizing field arising from the boundary condition applied at the edge of an infinite periodic sum, expressed via the total system dipole moment.
Central, index-based store of PolarSegment objects, keyed by segment id and charge state.
Reciprocal-space term of a 3D Ewald sum: the field at a target position from the Gaussian-screened periodic image of every registered charge and (static or induced) dipole.
Real-space term of a 3D Ewald sum: the total erfc(alpha*r)-screened field at a target position from every registered segment, summed over periodic images, converged by adaptive radial shells.
Real-space, erfc(alpha*r)-screened multipole interaction tensor for use in a 3D Ewald summation.
Matrix-free periodic induced-dipole self-consistent-field operator, for use with Eigen::ConjugateGradient.
Takes a list of atoms, and the corresponding density matrix and puts out a table of partial charges.
Screening of a GW-BSE calculation by a classical polarizable environment, folded into the screened interaction rather than iterated against it.
Small container for the individual energy terms in a polar region.
Writes an orbital file to a .cube file.
Small container to keep occupation of BSE states for each atom.
Basic Container for QMAtoms,PolarSites and Atoms.
Parser to read strings containing indexes in the format "1 2:7 8" and returning a sorted expanded std::vector<Index> with only unique entries It can also do the opposite.
A graph visitor class for creating graph visitor objects.
A graph visitor determines the graph topology.
This file is a compilation of graph related algorithms.
A breadth first (DF) graph visitor.
A breadth first (BF) graph visitor.
Implements relative method - do not use for comparing with zero use this most of the time, tolerance needs to be meaningful in your context
Function taken from https://stackoverflow.com/questions/17333/what-is-the-most-effective-way-for-float-and-double-comparison user ShitalShal
This graph visitor will explore the vertices closest to the starting node first and proceed outwards.
These algorithms require the interplay of the graph and graph visitor classes and thus cannot be made methods of either. In most cases a graph visitor is specified explores the graph to determine some metric
This visitor will calculate the distance of each node from the starting node it is built on top of the graph breadth first visitor. As the visitor moves through the graph it adds a 'Dist' attribute to each graph node with an integer value corresponding to how far it is removed from the starting node.
E.g.
0 - 1 - 2 - 3
If vertex 1 is the starting vertex than the graph node associated with vertex 1 will have a distance of 0. Vertices 0 and 2 a distance of 1 and vertex 3 a distnace of 2.
This class serves as the base class for creating graph visitor objects it contains the framework by which all derived graph visitors should follow. The function of the graph visitor is to explore a graph, this can be done in different ways such as a breadth first or depth first algorithm or others.
The general flow of the class is as follows:
Follows Li, D'Avino, Duchemin, Beljonne and Blase, Phys. Rev. B (2018) (arXiv:1801.01755), with a Thole induced-dipole environment in place of their charge-response model. The bare Coulomb interaction of the QM subsystem is replaced by
v -> v + v_reac, v_reac = v_12 chi^(2) v_21,
the static response of the environment to a change in the QM charge density, felt back by the QM subsystem. For inducible point dipoles, mu = A^-1 E with A = alpha^-1 + T_Thole – the same operator PolarRegion solves – so
v_reac = -F A^-1 F^T,
with F mapping a charge density to the field it produces at the sites. It is negative semidefinite: it screens.
This class provides the pieces, each separately testable:
AuxFieldAtPoints F, per auxiliary function ReactionFieldKernel B = -F A^-1 F^T for the explicit Thole region ShellKernel the cheap tail beyond it (legacy radial_dielectric) SymmetrizedReactionField R = T^T B T, in the metric of the stored three-centre integrals
B lives between auxiliary functions taken as charge densities. In the RI representation used by TCMatrix_gwbse, M = (mn|Q) T with T = TCMatrix_gwbse::InvSqrt(), a density rho_mn has fit coefficients c = V^-1 (Q|mn), and
(mn|v_reac|kl) = c_mn^T B c_kl = M_mn (T^T B T) M_kl^T,
using V^-1 = T T^T. That is what R = T^T B T is for: it acts on M exactly as the bare interaction does, which is the identity in that metric.
This mirrors DipoleDipoleInteraction's own existing interface (see dipoledipoleinteraction.h) exactly – same required Eigen::EigenBase typedefs/constants, same InnerIterator/operator()/multiply() shape – so it plugs into Eigen::ConjugateGradient the same way PolarRegion::CalcInducedDipolesViaPCG already uses DipoleDipoleInteraction, with the periodic real+reciprocal field (EwaldRealSpaceSum + EwaldReciprocalSpaceSum) standing in for DipoleDipoleInteraction's own aperiodic eeInteractor:: FillTholeInteraction wherever the "interaction between two sites" building block is needed.
Segments in ids may have any number of sites each (this was not always true of this class – an earlier version required exactly one site per segment; generalizing it was deliberately sequenced after that simpler version was already validated, so that any new bug introduced by the generalization itself would be easy to isolate from bugs in the underlying physics, which by this point had already been separately confirmed). The vector layout any v/x this operator is applied to must follow is: for each id in ids (in order), 3 entries per site of that segment (xyz), in the order EwaldRegistry::Get(id, Neutral)'s own site iteration gives them – i.e. a variable-width, offset-table ("CSR-style") layout, not a fixed 3-per-segment one.
multiply(v) computes A*v, where A is the induced-dipole self-consistent operator restricted to ids: A*v = D*v + (periodic field that induced dipoles v produce at every ids site, from every OTHER ids site's periodic images, at every translation including zero – see below for why zero-translation same-segment pairs are no longer excluded here)
The intramolecular term above was NOT present in an earlier version of this class, based on a mistaken belief that DipoleDipoleInteraction excludes intra-segment coupling the same way EwaldRealSpaceSum does for the periodic sum (see that class's own documentation) – checking DipoleDipoleInteraction::multiply() directly shows this is wrong: its inner loop couples every pair of sites in its own flat sites_ list via FillTholeInteraction, with no segment-awareness at all beyond operator()'s own same-site (not same-segment) diagonal check. Both legacy PolarBackground (see its own explicitly-labeled "Real-space, intramolecular contribution" loop, inside its own induction-iteration section specifically – NOT its permanent-field section, which has no intramolecular loop at all) and modern PolarRegion/ DipoleDipoleInteraction therefore include intramolecular induced- induced coupling as standard practice; omitting it here was an actual correctness gap in this class, not a deliberate scoping choice, caught only once a real full-polarizability comparison against legacy surfaced a large, systematic, atom-specific discrepancy that a static- field-only (Tier 1) comparison could never have exposed, since this term contributes nothing until induction is actually active.
This class's own damping choice for that term went through TWO mistaken revisions before landing on Thole-damped, worth recording in full since both were confident and both were wrong in different directions:
A separate, related check was also worth doing directly rather than by analogy: is legacy's OWN Thole damping model here (its L3()/L5(), called via UpdateAllBls from FU12_ERFC_At_By) the same exponential (Thole/van Duijnen-Swart) form ComputeThole already uses, or a different (e.g. linear) one? An earlier investigation concluded "linear" – but that conclusion was traced from XInteractor, which PolarBackground constructs and uses ONLY inside its own explicitly- labeled "CUTOFF TREATMENT" branches (coulombmethod.method="cutoff", not the "ewald" default this class is actually compared against). Tracing EwdInteractor's own L3()/L5() directly instead – confirmed as the class PolarBackground actually uses by default, via its own _ewdactor member and the "ewald" default in ewdbgpol.xml's own schema – shows L3()=1-exp(-ta1*tu3), L5()=1-(1+ta1*tu3)*exp(-ta1*tu3): the SAME exponential form ComputeThole already implements, with ta1 traced back to the same polarmethod.aDamp/thole_a parameter. So this class's own Thole model was never the mismatch; AddIntraSegmentCoupling calling ComputeThole now (with the real thole_a – see this class's own constructor documentation) is a genuine fix, not a model-form change.
A FOURTH, related mismatch was found in EwaldReciprocalSpaceSum: an earlier version of that class excluded the target's own segment from its own reciprocal-space sum entirely, modeled on EwaldRealSpaceSum's own real-space exclusion. Tracing legacy's own reciprocal-space code (PolarBackground::KThread::SP_SFactorCalc/FP_KFieldCalc for the permanent case, SU_SFactorCalc/FU_KFieldCalc for the induced case) shows it never excludes anything – every registered site, including the target's own segment, contributes to the one global structure factor applied back to every site. This matters specifically because of the Ewald-split identity erfc(a*r)/r + [reciprocal-space contribution] = 1/r EXACTLY: if EwaldRealSpaceSum's own real-space sum legitimately excludes a same-segment, zero-translation pair (which it still correctly does, unchanged), the reciprocal-space contribution to that SAME pair must still be included for the two Ewald halves to combine into the correct total wherever they're meant to – e.g. here, where the erfc-screened intramolecular term above is only the real- space HALF of the intended full-strength (Thole-damped) intramolecular interaction; the other half arrives automatically once EwaldReciprocalSpaceSum stops excluding anything at all (see that class's own documentation for the fuller account). This is why the periodic-field term's own documentation above no longer says "excluding that site's own segment's zero-translation copy" the way an earlier revision of this comment did – only EwaldRealSpaceSum's own real-space exclusion remains; EwaldReciprocalSpaceSum's own contribution is unconditional now.
The periodic-field term is computed by temporarily writing v into registry_'s induced-dipole slots and evaluating EwaldRealSpaceSum:: AddFieldAt / EwaldReciprocalSpaceSum::AddFieldAtMany at every ids site; confirmed symmetric for this combination when restricted to ids-only coupling (see test_ewaldperiodicdipolesymmetry.cc – written against the earlier, segment-excluding version of EwaldReciprocalSpaceSum; worth re-confirming this still holds now that exclusion is gone entirely, though the underlying argument – a symmetric physical kernel stays symmetric whether or not a term is added to both sides of a pair uniformly – does not depend on that exclusion having existed), a prerequisite for ConjugateGradient. The intramolecular term is added/transposed into both sites' own result the same block/ block.transpose() pairing DipoleDipoleInteraction::multiply() itself uses, which keeps the combined operator symmetric overall. This intramolecular contribution is linear in v and vanishes at v=0, so it does not affect baseline_ (see below) at all.
The sign/scale convention relative to DipoleDipoleInteraction's own established one has since been checked against a physically unambiguous case – a single polarizable site in one fixed external charge's field must develop an induced dipole aligned with that field – and matches (see test_ewaldperiodicdipoleoperator.cc); this does not by itself prove correctness for larger, genuinely self-consistent (multi-site, multi-iteration) cases, which is worth validating separately before trusting this for anything real – and is worth re-validating again specifically for the genuinely-multi-site-per-segment case this generalization introduces, since neither existing test exercises that case at all. The intramolecular-coupling addition above is itself NOT YET covered by either existing test (both use single-site segments where the term is identically zero) – validating it directly, e.g. against a small system where PolarRegion's own DipoleDipoleInteraction can be run for comparison, is a genuine open task, not something this comment should be read as already having confirmed.
A genuine bug was caught and fixed here during that same validation: AddFieldAt/AddFieldAtMany sum the field from every OTHER registered segment, not just other members of ids – so a naive multiply(v) would silently include a constant contribution from any registered segment outside ids (e.g. a fixed external background), independent of v. That breaks the linearity ConjugateGradient requires (multiply(0) must be 0). The fix: baseline_ (computed once at construction, as the same raw computation evaluated at v=0) captures exactly that external leak, and multiply(v) subtracts it, leaving only the genuine ids-internal coupling. Segments outside ids still correctly influence the physics overall – through the permanent field baked into b before this operator is ever invoked – just not redundantly inside A itself.
IMPORTANT, and easy to get backwards: the right-hand side b to pass to ConjugateGradient::solveWithGuess is b = +V_permanent, not b = -V, even though PolarRegion::CalcInducedDipolesViaPCG itself uses b = -V for the (legacy-derived) DipoleDipoleInteraction operator. The two are not interchangeable conventions to copy blindly: eeInteractor's own V() stores the negative of the physical field (confirmed via its own energy expression, e = q*phi + mu.V, which only matches the standard U = q*phi - mu.E if V = -E), whereas EwaldRealSpaceInteractor/ EwaldRealSpaceSum/EwaldReciprocalSpaceSum's own V() stores the genuine, un-negated physical field (confirmed via test_ewaldrealspaceinteractor .cc's own explicit CalcStaticEnergy formula, q*phi - mu.E). This sign mismatch was caught the hard way – as a factor-of-exactly-(-1) discrepancy in test_ewaldperiodicdipoleoperator.cc, after the double-counting bug above was already fixed and could no longer explain it – rather than reasoned out from first principles, which is exactly why it is worth stating this explicitly here rather than trusting a future caller to rediscover it.
operator()(i,j) is used by Eigen's preconditioner setup (via InnerIterator) to extract the literal diagonal only; DiagonalPrecon- ditioner never reads an off-diagonal value even though it is visited during setup. DipoleDipoleInteraction's own operator()(i,j) computes the real interaction tensor for every cross-site pair regardless (its own comment notes this is "not a fast method"), which is tolerable there because that tensor is O(1) per pair (aperiodic). Here, the periodic per-pair tensor would cost a full real+reciprocal evaluation, so cross-site entries return 0.0 unconditionally rather than computing (and discarding) the real value – same-site entries still return the genuine getPInv() block, exactly as DipoleDipoleInteraction does, since that is both cheap and the only part of operator()'s output the preconditioner actually keeps. "Same-site" here means the literal same PolarSite, not merely the same segment – two different sites of the same multi-site segment are cross-site entries here too (0.0), even though multiply() now does couple them via the intramolecular term described above: that coupling is never read by DiagonalPreconditioner either (same reasoning as the periodic term), so operator()'s own 0.0 for it is equally safe to leave uncomputed, not an oversight relative to multiply()'s own, more complete behavior.
getPInv() alone was NOT actually the true diagonal for a real interval spanning this class's own reciprocal-space-no-exclusion change and the self_field_matrix_ fix below: with the former in place but not yet the latter, a site's own true diagonal was getPInv() + (that site's own spurious reciprocal-space self-field contribution, see self_field_matrix_'s own documentation) – a real, if modest-looking (~4% relative to getPInv() for the specific system this was first noticed on), mismatch between what operator() reported and what multiply() actually used, silently degrading DiagonalPreconditioner's own effectiveness (observed as PCG iteration counts roughly quadrupling for one real system, alongside the separate, larger reciprocal-space-related slowdown from that same change). Adding self_field_matrix_ back into multiply()'s own diagonal (see below) cancels that same self-field term there too, restoring getPInv() as the genuine true diagonal and operator()'s own reported value as correct again – not a coincidence: both fixes address the same underlying leak, just in two different places it needed correcting.
Registered segments not included in ids (the set this operator solves for) are treated as a fixed, non-reacting external background: their own sites still contribute to the periodic field seen by the segments in ids (via the same intermolecular-only exclusion already implemented in EwaldRealSpaceSum/EwaldReciprocalSpaceSum), but their own induced dipoles are never themselves solved for or perturbed by multiply(). This is what lets this operator solve for the induced dipoles of, e.g., a growing/adaptive subset of sites without needing to re-solve the whole periodic system's own dipoles at once.
registry, real_sum, and recip_sum are all stored by reference and must outlive this object; real_sum and recip_sum must have been constructed against the same registry passed here.
Units: bohr, Hartree atomic units throughout.
This mediates pairwise interactions between charge (rank 0) and dipole (rank 1) sites, screened for use as the real-space term of an Ewald split. Quadrupoles (rank 2) are out of scope for this first version; see the project notes for why.
This is an independent implementation, not a port of the legacy xtp/ewald/ code: it operates directly on PolarSite/StaticSite (bohr, atomic units, and their actual stored multipole ordering) rather than the legacy APolarSite representation (nm, different ordering). The B-function recursion below was cross-checked against the legacy EwdInteractor's own UpdateAllBls for correctness, but the surrounding structure, naming, and the field/energy expressions built on top of it were derived fresh.
Damping conventions:
Units and sign convention: bohr, Hartree atomic units throughout, matching PolarSite/StaticSite. Fields point from source to target and are added into the target's V() / V_noE() accumulators, matching eeInteractor's own convention.
This is a fresh implementation of the algorithm identified in the legacy xtp/ewald/ code's RThread::FP_FieldCalc (real-space contribution grows shell-by-shell until the per-shell contribution becomes negligible), not a port of it. In particular the convergence criterion is re-derived: the legacy code's own "test dipole" / RMS-count construction was never fully understood even by the people who wrote it, so this class instead tracks the plain magnitude of each shell's own field contribution and stops once that drops below an absolute tolerance (with a minimum radius enforced throughout, so an accidentally-small early shell can't trigger early stopping). This is a considered simplification, not a proven-equivalent reformulation – see the project notes on validating its own convergence behaviour (does the total field actually stabilize as the tolerance tightens?) before trusting it on a production system.
Scope: this class computes the intermolecular real-space field only. Interactions between two sites in the same registered segment (whether at zero translation, i.e. genuinely intramolecular, or at a nonzero translation, i.e. that segment's own periodic image) are excluded entirely and left to whatever already handles intramolecular coupling – this mirrors the existing split in eeInteractor between ApplyStaticField/ApplyInducedField (intermolecular) and Cholesky_IntraSegment/FillTholeInteraction (intramolecular), rather than inventing a new boundary.
Performance: AddFieldAt caches the real, geometry-determined neighbor set (which (source segment, periodic-translation) pairs actually contributed to a given target's own converged sum) per (target site, source charge state) pair, keyed by that site's own address together with source_state. The first call for a given (target, state) pair runs the full shell-by-shell convergence search (scanning every registered segment each shell); every subsequent call for that SAME pair reuses the cached (source, translation) list directly, skipping the search entirely. This mirrors legacy PolarBackground's own RThread:: PolarNbs() mechanism (built once, at its own first SOR iteration, then reused for every later one) rather than being a novel optimization – a direct comparison against a legacy log for the same real system found this class's own repeated full-system rescan, once per PCG iteration, to be the dominant cause of a real ~8x wall-clock slowdown relative to legacy's own SOR solve at a comparable iteration count (with thread count and iteration count both directly ruled out as explanations first). This caching is safe specifically because site POSITIONS are fixed for the lifetime of a single PCG solve (only induced dipole VALUES change between iterations, which this cache never stores) – it would NOT be safe to reuse across a call sequence where target or source positions themselves change.
Units: bohr, Hartree atomic units throughout, matching PolarSite and EwaldRealSpaceInteractor.
This is derived from scratch, not ported from the legacy xtp/ewald/ code. Ewald's trick surrounds each point multipole with a canceling Gaussian charge distribution; that periodic Gaussian density's own potential satisfies Poisson's equation in Fourier space (phi_k = 4*pi*rho_k/k^2), which gives, for a charge+dipole source set:
S(k) = sum_j [ q_j - i*(k.mu_j) ] * exp(-i*k.r_j) (structure factor) phi(r) = (4*pi/V) * sum_{k!=0} (1/k^2) * exp(-k^2/4a^2) * S(k) * exp(i*k.r) E(r) = (4*pi/V) * sum_{k!=0} (k/k^2) * exp(-k^2/4a^2) * Im[ S(k) * exp(i*k.r) ]
summed over the full set of k!=0 (not a k/-k half-space with an extra factor of 2) – pairing k and -k terms is exactly what makes this sum real by construction, and is worth stating explicitly since getting it wrong (e.g. by introducing a spurious factor of 2, or its inverse) would not be caught by compilation, only by a numerical check against a known reference. See the project notes on validating this class for exactly that reason.
k-vector truncation is a simple |k| < k_max spherical cutoff, not the legacy code's per-vector "graded" adaptive truncation. That grading is self-contained (depends only on structure-factor magnitude, not on any real-space convergence state – confirmed by tracing the legacy code directly), so it is addable later as a genuinely separable enhancement to the truncation strategy, on top of this simple-cutoff implementation, without needing to revisit the underlying k-vector definitions or field formula above.
Scope: unlike EwaldRealSpaceSum, this includes EVERY registered site – no exclusion of the target's own segment at all, even for the intramolecular case. This was not always true of this class: an earlier version excluded the target's own segment entirely, modeled on EwaldRealSpaceSum's own real-space exclusion, based on a mistaken assumption that the two sums should treat exclusion the same way. Tracing legacy PolarBackground's own reciprocal-space code directly (PolarBackground::KThread::SP_SFactorCalc/FP_KFieldCalc for the permanent case, SU_SFactorCalc/FU_KFieldCalc for the induced case) shows it never excludes anything from its own global structure factor – every registered site contributes, and the resulting field is applied back to every site with no segment-awareness at all. This matters specifically because of how the Ewald split works: erfc(a*r)/r
Performance note: S(k) is a sum over every registered site, and is identical for every target sharing the same source_state – no per-target correction of any kind is needed now (there is nothing target-specific left to correct for), so AddFieldAtMany's own benefit over repeated AddFieldAt calls is now simply computing S(k) once per k-vector and reusing it across every target in the batch, rather than once per (k-vector, target) pair. This is an O(N_k * N) batch, rather than the O(N_targets * N_k * N) that repeated single-target AddFieldAt calls would cost – the difference that matters once this is driven from an iterative solve calling it once per site per iteration. AddFieldAt (single target) is kept for convenience/testing and simply delegates to AddFieldAtMany with a one-element batch; it does not itself get the batching benefit, since there's nothing to share it with.
Units: bohr, Hartree atomic units throughout.
This is the modern replacement for the legacy Ewald code's PolarTop / PolarSeg / APolarSite machinery. Where the legacy design shared raw PolarSeg* pointers across multiple containers (BGN/MGN/FGN/FGC, QM0/MM1/MM2) and reassigned them at runtime with manually-maintained ownership flags, EwaldRegistry instead holds every PolarSegment by value, in one place, and everything else (background machinery, regions, ...) refers to a segment by its stable (id, charge state) key rather than by pointer. "Moving" a segment between roles becomes a change in which key is queried, never a change in who owns the underlying data.
A segment that only ever needs one representation (the common case: the neutral background) is registered once, under EwaldChargeState::Neutral. A segment that additionally needs an excited-state representation (the foreground molecule under study) is registered a second time, under EwaldChargeState::Electron or EwaldChargeState::Hole, alongside its existing neutral entry. Both entries remain independently addressable; neither is preferred or "active" by default. Callers explicitly say which state they want when they query.
EwaldRegistry deliberately does not know how to build a PolarSegment from a Topology; that is a separate concern (see SegmentMapper). This class is only ever a store.
Derived from scratch (De Leeuw, Perram & Smith 1980's standard result for the order-of-summation-dependent "surface term" of a conditionally convergent dipole lattice sum), not ported from the legacy code, though the resulting depolarization factors were cross-checked against it and do match: 4*pi/(3*V) for an isotropic (cube/sphere) boundary, or 4*pi/V acting on the z-component only for a slab (2D-periodic, z-normal) boundary – both under the vacuum-boundary convention (surrounding dielectric constant = 1; a general surrounding-medium factor is not implemented here since the legacy code's own use of it was always at this vacuum default).
M = sum_j (q_j * r_j + mu_j) (total system dipole moment) E = -(4*pi / (3*V)) * M (cube/sphere) E = -(4*pi / V) * M_z * z_hat (slab)
Unlike EwaldRealSpaceSum/EwaldReciprocalSpaceSum, this is not an intermolecular-only sum: M includes every registered segment, including the target's own. The depolarizing field is a mean-field property of the whole sample's own surface polarization acting back on itself (the same way a uniformly polarized dielectric sphere experiences its own depolarizing field), not a pairwise molecule-molecule interaction, so there is no analogous self-interaction to exclude.
Units: bohr, Hartree atomic units throughout.
Extracted from EwaldBackground (where they began life as file-local helpers in a calculator header) for two reasons: the MM/MM and QM/MM calculators need the same solvers without including a calculator, and a private calculator header cannot be unit-tested – the test target does not put xtp/src/libxtp on its include path, so the solvers went untested while every other Ewald layer had its own test.
Two solvers are provided, deliberately:
SolveWithJOR – Jacobi over-relaxation, reproducing legacy PolarBackground's own SOR scheme (same criterion, same omega) so the two codes can be compared iteration by iteration. SolveWithIndefinitenessCheck – preconditioned conjugate gradient, far fewer iterations on the same system, with a per-iteration curvature check (see below).
Both converge to the same fixed point; that equivalence is asserted directly in test_ewaldsolvers.cc rather than assumed.
How energies are evaluated depends critically on the id of the region. a) Lower ids means being evaluated later. So expensive regions should have low ids, as they can then already incorporate partially converged results from higher id regions b) The energy of a region, includes all interactions with regions of higher ids. i.e. E(region0)=E0+E01+E02, whereas E(region1)=E1+E12 and not E10 The reason is that DFT codes only return E0+ all interaction energies and not the individual terms, this has to be changed once we want to have multiple QM regions.
PATHWAYS TO A NEW THREADED CALCULATOR ... 1 Define 'JobContainer' (needs to define iterator), 'pJob' ( = *iterator) ... 2 Derive new calculator as ': public ParallelXJobCalc<JobContainer,pJob>' ... 3 Specialize XJOBS_FROM_TABLE< JobContainer, pJob> in xjob.cc ... 4 Register new calculator (see end of parallelxjobcalc.cc) REQUIRED MEMBERS FOR pJob pJob::JobResult (struct)
Converges the self-consistent, fully polarized (neutral-state) background of the whole periodic system, and writes it to an HDF5 checkpoint file via EwaldRegistry::WriteToCpt.
User-facing length-related options (coulombmethod.alpha, coulombmethod.k_max, realspace.r_min) are specified in nm/nm^-1, not this class's own internal bohr/bohr^-1 units – converted once in ParseOptions, immediately on read. This is deliberate: everything downstream of ParseOptions (member variables, logging, the actual EwaldRealSpaceSum/EwaldReciprocalSpaceSum construction) stays in the same bohr-based units PolarSite/EwaldRegistry already use natively, so only the options-parsing boundary itself needs to know about the user-facing unit choice.
Segments may have any number of sites each – this now matches EwaldPeriodicDipoleOperator's own generalized (offset-table) vector layout directly, built the same way here for the same reason: an earlier version of both this class and the operator required exactly one site per segment, and that limitation was lifted from the operator first, validated there in isolation (see test_ewaldperiodicdipoleoperatormultisite.cc), before being carried through here – so that a bug in this generalization, if there were one, would be easy to isolate from the operator's own already-tested multi-site behavior.
The outer "loop until 2nd-order fields converged" structure the legacy PolarBackground::Polarize() itself uses (see the project notes on that class) is not reproduced here at all: EwaldPeriodicDipoleOperator's own ConjugateGradient solve already is that convergence loop, in one call, rather than a separate SOR iteration this class would need to drive itself.
For an earlier history see ctp repo commit 77795ea591b29e664153f9404c8655ba28dc14e9
For earlier commit history see ctp commit 77795ea591b29e664153f9404c8655ba28dc14e9
| using votca::Index = Eigen::Index |
1.16.1