vibeqc.run_pbc_bipole_rhf

vibeqc.run_pbc_bipole_rhf(system, basis, kmesh, options=None, *, linear_dep_threshold=1e-07, canonical_orth_normalize_diag_first=True, level_shift_schedule=None, use_mom=False, use_oda=False, oda_trust_lambda_max=1.0, use_incremental_fock=False, use_ewald_j_split=None, ewald_omega=None, ewald_precision=1e-08, v_ne_grid_options=None, use_multipole_diag=False, use_multipole_far_field=None, multipole_l_max=2, use_exchange_ewald_split=None, exchange_exxdiv='ewald', use_fock_symmetry=None, use_fock_symmetry_reduce=False, progress=None, verbose=None, initial_density=None, dft_plus_u=None)[source]

Multi-k closed-shell RHF via the CRYSTAL-gauge BIPOLE scaffold.

dft_plus_u: optional list of HubbardSite. When set, the Dudarev rotationally-invariant per-spin V_U is added to every per-k Fock matrix using the same per-spin Bloch-summed convention as run_pbc_bipole_uhf() (closed-shell: P_σ = P_total / 2, E_U_total = 2 × E_σ). The +U energy lands on result.e_dft_plus_u.

Algorithm (real-space two-electron / bielectronic build):

  1. Real-space one-electron integrals S(g), T(g), V_ne(g) at opts.lattice_opts.cutoff_bohr. For 3D systems V_ne uses the same Ewald α as E_nn.

  2. Bloch-sum to S(k), Hcore(k) per k-point; canonical-orth X(k).

  3. Initial guess via opts.initial_guess (default SAD).

  4. SCF iter: a. Build F^{2e}(g). With use_ewald_j_split=True this is

    J_SR(g;ω) + J_LR(g;ω) + V_bg·S(g) - ½K(g) where the exchange convention depends on use_exchange_ewald_split (below). With the flag off, use the legacy direct-only build_fock_2e_real_space scaffold.

    1. Bloch-sum F^{2e}(g) → F(k); add Hcore(k).

    2. Energy: E_elec = Σ_g tr[D(g)Hcore(g)] + ½Σ_g tr[D(g)F²e(g)] in real-space block form (real-space lattice-sum convention).

    3. Optional DIIS extrapolation of F(k) via [F,DS] errors.

    4. Optional LEVSHIFT shift on F(k).

    5. Diagonalise F(k) → C(k), ε(k).

    6. Optional MOM reorder of occupied subspace.

    7. Rebuild D_real via real_space_density_from_kpoints.

    8. Optional ODA mixing on density.

  5. E_total = E_elec + E_nuc.

use_ewald_j_split defaults to None. In that mode the driver automatically uses the CRYSTAL-gauge Ewald-J split for 3D systems and keeps the old direct-only path for dim < 3 diagnostic runs. Pass False explicitly only when you want the legacy direct-only F²e scaffold for debugging. (Passing True on a dim < 3 system raises — the Ewald split needs a 3D reciprocal lattice.)

use_exchange_ewald_split (2026-06-10 energy-assembly redesign; multi-k q≠0 channels 2026-06-11; multi-k default flip 2026-06-13) defaults to None = auto: ON for any 3D run under the Ewald J split (Γ AND multi-k), OFF otherwise. When ON, the exchange uses the Ewald split convention (module docstring of vibeqc.bipole_fock_ewald):

K(k) = K_SR(erfc ω, direct) + K_LR(erf ω, reciprocal, q+G≠0)
       + (ξ_M − π/(V_sc·ω²))·S(k)·D(k)·S(k)

with ξ_M the probe-charge Ewald (Madelung) constant of the BvK supercell (= the unit cell at Γ; V_sc = n_k·V) when exchange_exxdiv='ewald' (the default; PySCF-equivalent) or 0 when 'none'. At multi-k the LR term couples every k-point pair through the momentum-transfer channels q = k − k′ (see vibeqc.bipole_fock_ewald.compute_K_long_range_at_k()). In this mode the SCF density is the full Bloch fold (the Γ-locality projection P(g≠0)=0 is not applied), and the EXT EL-SPHEROPOLE term is omitted from the total — at the corrected gauge it is a double-count (MgO Γ fixed-density audit, 2026-06-10: the reassembled total matches PySCF GDF RHF to truncation with no spheropole term). When OFF (explicit False), the legacy convention is kept: full-Coulomb direct-space K, Γ-locality projection at n_k = 1, spheropole term added — known to mis-state absolute energies on tight ionic cells (kept only for the legacy-gauge analytic gradient + parity diagnostics). The corrected multi-k gauge needs a Monkhorst-Pack BlochKMesh carrying its mesh metadata; under the auto default an ad-hoc k-list (band path / explicit list) at multi-k falls back to the legacy gauge with a log note, and an explicit True with such a mesh raises.

For dim < 3 the whole one- and two-electron Coulomb gauge falls back to DIRECT_TRUNCATED (no Ewald, no reciprocal sum), and the EXT EL-SPHEROPOLE correction — a 3D-Ewald reciprocal-space (K=0 limit) term — is identically zero, so it is omitted and e_ext_el_spheropole is None. The resulting energy is the direct-truncated value: vacuum-independent and equal to the molecular RHF energy in the isolated-cell limit (see tests/test_pbc_bipole_dim_lt3.py).

For 3D systems the default V_ne implementation is analytic: erfc-screened nuclear attraction from libint plus a reciprocal-space AO-pair Fourier-transform sum. Passing v_ne_grid_options opts into the older grid-quadrature long-range V_ne path for diagnostics.

Parameters:

  • system (PeriodicSystem)

  • basis (BasisSet)

  • kmesh (BlochKMesh)

  • linear_dep_threshold (float)

  • canonical_orth_normalize_diag_first (bool)

  • level_shift_schedule (Optional['LevelShiftSchedule'])

  • use_mom (bool)

  • use_oda (bool)

  • oda_trust_lambda_max (float)

  • use_incremental_fock (bool)

  • use_ewald_j_split (Optional[bool])

  • ewald_omega (Optional[float])

  • ewald_precision (float)

  • v_ne_grid_options (Optional[GridOptions])

  • use_multipole_diag (bool)

  • use_multipole_far_field (Optional[bool])

  • multipole_l_max (int)

  • use_exchange_ewald_split (Optional[bool])

  • exchange_exxdiv (str)

  • use_fock_symmetry (Optional[bool])

  • use_fock_symmetry_reduce (bool)

  • progress (Union[bool, ProgressLogger, None])

  • verbose (Optional[int])

  • initial_density (Optional[Sequence[np.ndarray]])

  • dft_plus_u (Optional[List['HubbardSite']])

Return type:

PBCBipoleRHFResult