|
votca 2026-dev
|
#include <environmentscreening.h>
Static Public Member Functions | |
| static Eigen::MatrixXd | AuxFieldAtPoints (const AOBasis &auxbasis, const std::vector< Eigen::Vector3d > &points, const std::vector< double > &widths={}) |
| F: the electric field at each point produced by each auxiliary basis function, taken as a charge density. | |
| static std::vector< double > | SiteWidths (const std::vector< PolarSegment > &segments, double site_width) |
| R per site, in bohr: site_width * alpha_iso^(1/3), with alpha_iso = tr(alpha)/3. | |
| static Eigen::MatrixXd | ReactionFieldKernel (const Eigen::MatrixXd &F, const std::vector< PolarSegment > &segments, double exp_damp) |
| B = -F A^-1 F^T for an explicit Thole region. | |
| static Eigen::MatrixXd | ShellKernel (const Eigen::MatrixXd &F, const std::vector< PolarSegment > &segments, double epsilon) |
| The tail beyond the explicit region: legacy's radial_dielectric. | |
| static Eigen::MatrixXd | SymmetrizedReactionField (const Eigen::MatrixXd &B, const Eigen::MatrixXd &T) |
| R = T^T B T, the reaction field in the metric of the stored three-centre integrals, where the bare interaction is the identity. | |
| static Eigen::MatrixXd | DressingMatrix (const Eigen::MatrixXd &R) |
| S = (1 + R)^(1/2), the dressing of the auxiliary index that turns the bare interaction into u = v + v_reac in the RPA and the correlation self-energy (TCMatrix_gwbse::DressAuxIndex). | |
| static Eigen::MatrixXd | Kernel (const AOBasis &auxbasis, const ScreeningEnvironment &env) |
| B for a whole environment: ReactionFieldKernel of the explicit segments plus ShellKernel of the shell, over the auxiliary basis. The two do not couple to each other, as in legacy. | |
| static Eigen::MatrixXd | Metric (const AOBasis &auxbasis) |
| T for an auxiliary basis on its own: the metric TCMatrix_gwbse::Fill folds into the three-centre integrals, built the same way (Pseudo_InvSqrt_GWBSE with the same tolerance), so R from it is the R the GW-BSE run will see. | |
| static ScreeningCheck | Check (const AOBasis &auxbasis, const QMMolecule &atoms, const ScreeningEnvironment &env, const Eigen::MatrixXd &T, Index n_modes=3) |
| Builds R for env and reports whether 1 + R is positive definite, and why not. | |
Definition at line 121 of file environmentscreening.h.
|
static |
F: the electric field at each point produced by each auxiliary basis function, taken as a charge density.
Returns an N_aux x 3 N_points matrix; column 3p+k is field component k at point p. Rows follow the auxiliary basis in exactly the order and normalization TCMatrix_gwbse::Fill uses, because both go through AOBasis::GenerateLibintBasis and getMapToBasisFunctions.
The potential of a single auxiliary function at a point is a libint2 nuclear-attraction integral against Shell::unit(), which needs no derivative support from libint2. The field is its gradient with respect to the point, taken by a four-point central difference with h = 1e-3 bohr. Measured against the closed-form field of a diffuse s-Gaussian (exponent 0.12) at 2 to 15 bohr, that is accurate to about 1e-12 relative; the two-point stencil at the same h is only 1e-8.
The points are assumed to lie outside the auxiliary functions' significant extent, as the method itself assumes – the environment does not overlap the QM subsystem. Inside a function the field is still correct, but the reaction field would no longer mean anything.
Definition at line 61 of file environmentscreening.cc.
|
static |
Builds R for env and reports whether 1 + R is positive definite, and why not.
Everything R depends on – the QM atoms, the auxiliary basis on them, the polar sites and their polarizabilities – is fixed before the QM/MM loop starts, so this is what a job runs up front rather than discovering the problem after the ground state has converged.
The eigenvalues of R are those of B c = lambda V c: for the charge distribution rho = sum_Q c_Q chi_Q, lambda is its reaction energy relative to its own Coulomb energy. For inducible points outside rho that ratio is bounded by the dielectric limit, above -1. What can push it lower is a site inside the tail of rho, where the undamped point-dipole response is no longer bounded – so the report lists, for the n_modes lowest modes, their charge, the QM atoms and auxiliary shells carrying them, and the polar sites carrying their reaction, with the distance of each to the nearest QM atom.
Definition at line 406 of file environmentscreening.cc.
|
static |
S = (1 + R)^(1/2), the dressing of the auxiliary index that turns the bare interaction into u = v + v_reac in the RPA and the correlation self-energy (TCMatrix_gwbse::DressAuxIndex).
Symmetric positive definite. Throws unless 1 + R is.
Definition at line 349 of file environmentscreening.cc.
|
static |
B for a whole environment: ReactionFieldKernel of the explicit segments plus ShellKernel of the shell, over the auxiliary basis. The two do not couple to each other, as in legacy.
Definition at line 323 of file environmentscreening.cc.
|
static |
T for an auxiliary basis on its own: the metric TCMatrix_gwbse::Fill folds into the three-centre integrals, built the same way (Pseudo_InvSqrt_GWBSE with the same tolerance), so R from it is the R the GW-BSE run will see.
Definition at line 363 of file environmentscreening.cc.
|
static |
B = -F A^-1 F^T for an explicit Thole region.
F must be AuxFieldAtPoints at the positions of the sites of segments, in iteration order. A is assembled densely, block for block as DipoleDipoleInteraction::multiply builds it, from the same eeInteractor PolarRegion uses, so the environment responds exactly as the polar region does in the iterative scheme.
Dense on purpose. B needs A^-1 applied to N_aux right-hand sides at once. A matrix-free CG re-evaluates every Thole block on every iteration of every solve; assembling A once and factorizing it costs one pass over the pairs plus O((3N)^3), which is far cheaper for any polar region that fits in memory. Formed as -Y^T Y with Y = L^-1 F^T from the Cholesky factor, so B is symmetric and negative semidefinite by construction rather than up to round-off.
Throws if A is not positive definite – a Thole polarization catastrophe, i.e. sites too close together for their damping – since then there is no stable induced-dipole response to speak of.
Definition at line 222 of file environmentscreening.cc.
|
static |
The tail beyond the explicit region: legacy's radial_dielectric.
Legacy (Ewald3DnD::EvaluateRadialCorrection) took a shell of segments outside the polar cutoff, applied the QM field screened by 1/epsilon, induced directly – no mutual induction – and took the unscreened interaction back with the QM region. That is exactly an effective polarizability alpha_j / epsilon per site with no T coupling, so its susceptibility is block-diagonal and
B_shell = - sum_j F_j (alpha_j / epsilon) F_j^T,
which needs no linear solve. The 1/epsilon stands in for the mutual induction it leaves out. Because the sum runs over real sites, the geometry is whatever the morphology is – a slab simply has no sites in the vacuum – which is why this, rather than an analytic Born term, is the tail treatment.
Definition at line 252 of file environmentscreening.cc.
|
static |
R per site, in bohr: site_width * alpha_iso^(1/3), with alpha_iso = tr(alpha)/3.
Why smear at all: for a charge density rho and inducible points outside it, <rho|v_reac|rho> is bounded by the dielectric limit, above -<rho|v|rho>. A point dipole inside the tail of a diffuse auxiliary function is not bounded, and 1 + R can lose positive definiteness without any unphysical geometry: measured on a production QM/MM job (55 QM atoms, def2-tzvp/aux-def2-tzvp, 3900 Thole sites, closest contact 2.62 A), three carbons (alpha = 18.8 bohr^3) of a neighbouring C60 at 3.05-3.3 A carried 86% of a mode at lambda = -1.056. Smearing each site over a Gaussian of width proportional to alpha^(1/3) – the length scale Thole damping itself uses – removes that: at site_width 0.5 the same job has lambda_min = -0.755, while 1/2 <rho|v_reac|rho> for the HOMO, LUMO and HOMO-LUMO densities moves by 1e-4 relative (0.8: 2e-3, 1.0: 7e-3).
Definition at line 308 of file environmentscreening.cc.
|
static |
R = T^T B T, the reaction field in the metric of the stored three-centre integrals, where the bare interaction is the identity.
Throws unless 1 + R is positive definite. R itself is negative semidefinite, and 1 + R is the effective interaction u = v + v_reac in this metric; an eigenvalue of R at or below -1 would mean the environment screens a charge fluctuation by more than the fluctuation itself, which no stable environment does.
Definition at line 280 of file environmentscreening.cc.