- 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 on EwaldRegistry, EwaldRealSpaceSum, EwaldReciprocalSpaceSum, EwaldShapeCorrection and EwaldPeriodicDipoleOperator instead of the legacy xtp/ewald/ code. Registered as "ewaldbackground" so that both can coexist.
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 iterative solvers.
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.
Follows the Eigen::EigenBase interface of DipoleDipoleInteraction (dipoledipoleinteraction.h), with the periodic real- plus reciprocal-space field (EwaldRealSpaceSum, EwaldReciprocalSpaceSum, EwaldShapeCorrection) in place of the aperiodic eeInteractor::FillTholeInteraction.
Vector layout: for each id in ids (in order), 3 entries (xyz) per site of that segment, in the site order of EwaldRegistry::Get(id, Neutral). Segments may have any number of sites, so the layout is variable-width (see offsets_).
The induced dipoles solve mu = P*(F_perm + C*mu), i.e. A*mu = (P^-1 - C)*mu = F_perm, where C*mu is the field produced by the induced dipoles of the ids sites: the periodic real-space sum (excluding zero-translation same-segment pairs), the reciprocal-space sum (no exclusions), the shape/surface term, and the intramolecular Thole-damped coupling between distinct sites of the same segment (AddIntraSegmentCoupling). Like legacy PolarBackground and DipoleDipoleInteraction, intramolecular induced-induced coupling is included; it is Thole damped with the same exponential (van Duijnen-Swart) form as legacy EwdInteractor::L3()/L5() (FU12_ERFC_At_By damps intra- and intermolecular pairs alike).
Because the reciprocal-space sum excludes nothing, it also contains each site's spurious self-field -M*v_i, which is cancelled by self_field_matrix_ (as in legacy FU12_ERF_At_By with identical sites). With that correction getPInv() is the exact diagonal block of A.
Sign convention: b = +F_perm, not b = -V as in PolarRegion::CalcInducedDipolesViaPCG. eeInteractor::V() stores the negative of the field (e = q*phi + mu.V), whereas the Ewald classes' V() stores the physical field E (energy q*phi - mu.E).
The field sums include every registered segment, not only ids, so RawMultiply(v) contains a v-independent contribution from segments outside ids. multiply(v) = RawMultiply(v) - baseline_, with baseline_ = RawMultiply(0), which restores linearity (multiply(0) = 0). Segments outside ids thus act as a fixed, non-responding background; their influence enters only through the permanent field in b.
The operator is symmetric when restricted to ids (see test_ewaldperiodicdipolesymmetry.cc); the intramolecular term is added as a block/block.transpose() pair to preserve this. The sign/scale convention is tested on a single polarizable site in a point-charge field (test_ewaldperiodicdipoleoperator.cc), multi-site segments in test_ewaldperiodicdipoleoperatormultisite.cc.
operator()(i,j) is only used by Eigen's preconditioner setup, which keeps the diagonal. Same-site entries return the getPInv() block; all other entries (including distinct sites of one segment) return 0.0, since a periodic pair tensor would require a full Ewald evaluation.
registry, real_sum, recip_sum and shape are stored by reference and must outlive this object; the sums must be built on the same registry.
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.
Follows the shell-by-shell scheme of legacy RThread::FP_FieldCalc, but with a simpler convergence criterion: the sum stops once the magnitude of a shell's field contribution drops below an absolute tolerance, provided a minimum radius r_min has been reached (legacy uses a "test dipole" RMS construction instead).
Scope: the target's own segment at zero translation (the intramolecular pairs, including the r = 0 self-pair) is excluded; its other periodic images are included. Intramolecular coupling is handled elsewhere, as in the split of eeInteractor between ApplyStaticField/ApplyInducedField and Cholesky_IntraSegment/FillTholeInteraction.
Performance: the (source segment, translation) pairs that contribute to a target's converged sum are cached per (target site address, source charge state), so only the first call runs the shell search, as with legacy RThread::PolarNbs(). Without the cache the repeated search made each solver iteration ~8x slower than legacy. The cache assumes that site positions do not change during its lifetime; only induced-dipole values may 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 a factor of 2); pairing k and -k makes the sum real. A wrong factor of 2 would only show up in a numerical check against a reference.
k-vectors are truncated with a spherical cutoff |k| < k_max, not the graded per-vector truncation of the legacy code. That grading depends only on the structure-factor magnitude and could be added later without changing the k-vectors or the field formula.
Scope: unlike EwaldRealSpaceSum, every registered site contributes; the target's own segment is not excluded, as in legacy PolarBackground (KThread::SP_SFactorCalc/FP_KFieldCalc, SU_SFactorCalc/FU_KFieldCalc). Since erfc(a*r)/r plus the reciprocal-space term is exactly 1/r, a pair excluded from the real-space sum (same segment, zero translation) must still be included here for the two halves to combine correctly; see EwaldPeriodicDipoleOperator::AddIntraSegmentCoupling.
Performance: S(k) is the same for all targets with the same source_state, so AddFieldAtMany computes it once per k-vector, O(N_k * N) per batch instead of O(N_targets * N_k * N) for repeated AddFieldAt calls. AddFieldAt delegates to AddFieldAtMany with a one-element batch.
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.
Kept in a public header (rather than in the EwaldBackground calculator) so that both the MM/MM and QM/MM calculators can use them and so that they can be unit-tested.
SolveWithJOR – Jacobi over-relaxation, reproducing legacy PolarBackground's SOR scheme (same criterion, same omega) for iteration-by-iteration comparison. SolveWithIndefinitenessCheck – preconditioned conjugate gradient with a per-iteration curvature check.
Both converge to the same fixed point (tested in test_ewaldsolvers.cc).
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.
Length-related options (coulombmethod.alpha, coulombmethod.k_max, realspace.r_min) are given in nm or nm^-1 and converted to bohr or bohr^-1 in ParseOptions; everything downstream uses atomic units.
Segments may have any number of sites; b and x use the variable-width layout of EwaldPeriodicDipoleOperator.
The outer "loop until 2nd-order fields converged" of legacy PolarBackground::Polarize() is not reproduced: the iterative solve (PCG or JOR) is that loop.
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