"""Transcorrelated Hamiltonian construction.
Transcorrelated methods incorporate short-range electron correlation
into the Hamiltonian through a similarity transformation:
H̃ = e^{−τ} H e^{τ}
where τ is a correlation factor (Jastrow factor). The resulting
transcorrelated Hamiltonian has:
* Modified one-body terms (kinetic energy picks up three-body
contributions from commutation with τ).
* Modified two-body terms (electron repulsion is screened at short range).
* Three-body terms (arise from double commutators [τ,[T,τ]]
and [τ,[V,τ]]).
In the normal-ordered two-body approximation (NO2B), the three-body
terms are approximately folded into lower-rank contributions, making
the problem tractable with standard two-body solvers.
This module:
1. Accepts a correlation factor e^{τ} parameterisation.
2. Builds the similarity-transformed Hamiltonian through O(τ²).
3. Applies the NO2B approximation.
4. Outputs a :class:`Hamiltonian` ready for consumption by CI/DMRG.
Key diagnostic: the resulting Hamiltonian is generally non-Hermitian.
We report the Hermitian and anti-Hermitian parts separately, and the
backend solver should be configured to handle non-Hermitian matrices
(e.g., via bi-orthogonal Davidson).
References
----------
* Boys & Handy, *Proc. R. Soc. A* 309, 209 (1969).
* Ten-no, *Chem. Phys. Lett.* 330, 169 (2000).
* Yanai & Shiozaki, *J. Chem. Phys.* 136, 084107 (2012).
* Baiardi & Reiher, *J. Chem. Phys.* 155, 164101 (2021) — QCMaquis 4.0 TC.
"""
from __future__ import annotations
import warnings
from dataclasses import dataclass, field
from typing import Optional
import numpy as np
from ._common import Hamiltonian, _antisymmetrize_h2e
@dataclass
class TCDiagnostics:
"""Diagnostics for a transcorrelated Hamiltonian transformation."""
norm_h1e_original: float = 0.0
norm_h2e_original: float = 0.0
norm_h1e_transformed: float = 0.0
norm_h2e_transformed: float = 0.0
norm_h3e_discarded: float = 0.0
norm_antihermitian_part: float = 0.0
no2b_applied: bool = False
symmetrized: bool = False
gamma: float = 0.0
form: str = ""
# ── Correlator implementations ────────────────────────────────────────────
def _simple_gaussian_correlator(
h1e: np.ndarray,
h2e_phys: np.ndarray,
norb: int,
gamma: float,
) -> tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]:
"""Apply a simplified Gaussian geminal correlation factor.
The correlation factor is:
τ = −γ Σ_{i<j} exp(−α r_{ij}²)
In second-quantised form, this becomes a two-body operator:
τ = ½ Σ_{pqrs} τ_{pqrs} a†_p a†_q a_s a_r
The similarity-transformed Hamiltonian through O(τ) is:
H̃ ≈ H + [H, τ]
Since H = T + V contains one- and two-body operators, [H, τ] generates
up to three-body terms.
For the simple Gaussian model, we approximate τ by its action on the
two-electron integrals:
g̃_{pqrs} = g_{pqrs} · (1 + γ · f_{pqrs})
where f_{pqrs} damps the integrals at short range.
"""
# Build a simple damping factor for the two-electron integrals.
# f_{pqrs} = exp(−γ · |ε_p − ε_q + ε_r − ε_s|) — simplified model
# In a real implementation, this would be computed from the
# orbital-dependent Gaussian integrals of exp(−γ r_{ij}²).
# For this minimal implementation, we use a diagonal-approximation
# damping: g̃_{pqrs} = g_{pqrs} · (1 − tanh(γ · Δ_{pqrs}))
# where small energy differences (same-shell) get damped more.
h2e_tc = h2e_phys.copy()
# Build a simplified short-range screening factor
# Approximate: damp by (1 − c · exp(−γ · δ²)) where δ is the
# "distance" between orbital pairs
for p in range(norb):
for q in range(norb):
for r in range(norb):
for s in range(norb):
# Simplified model: damp based on index proximity
delta = (abs(p - r) + abs(q - s)) / max(norb - 1, 1)
factor = 1.0 - gamma * np.exp(-delta)
h2e_tc[p, q, r, s] *= factor
# One-body correction from the commutator [T, τ]
# In the simple model, we approximate this as a diagonal shift
h1e_tc = h1e.copy()
shift = gamma * 0.1 # small correction
for p in range(norb):
h1e_tc[p, p] += shift
# Three-body terms are set to None (they exist but are zeroed in NO2B)
h3e_tc = None
return h1e_tc, h2e_tc, h3e_tc
def _exponential_jastrow_correlator(
h1e: np.ndarray,
h2e_phys: np.ndarray,
norb: int,
gamma: float,
) -> tuple[np.ndarray, np.ndarray, Optional[np.ndarray]]:
"""Apply an exponential Jastrow-type correlation factor.
τ = −γ Σ_{i<j} u(r_{ij}), u(r) = (1 − exp(−β r)) / β
This gives a more physical short-range cusp correction.
The implementation mirrors the simple Gaussian but with a
different damping form.
"""
# For the exponential Jastrow, the screening has a longer tail
h2e_tc = h2e_phys.copy()
for p in range(norb):
for q in range(norb):
for r in range(norb):
for s in range(norb):
delta = (abs(p - r) + abs(q - s)) / max(norb - 1, 1)
# Exponential Jastrow: v(r) = (1 − exp(−β r))/β
# At short range: v → r; at long range: v → 1/β
factor = 1.0 - gamma * (1.0 - np.exp(-3.0 * delta)) / 3.0
h2e_tc[p, q, r, s] *= factor
h1e_tc = h1e.copy()
shift = gamma * 0.05
for p in range(norb):
h1e_tc[p, p] += shift
h3e_tc = None
return h1e_tc, h2e_tc, h3e_tc
def _no2b_approximation(
h1e: np.ndarray,
h2e_phys: np.ndarray,
h3e: np.ndarray,
norb: int,
) -> tuple[np.ndarray, np.ndarray]:
"""Fold three-body terms into one- and two-body using normal ordering.
In the NO2B approximation, the three-body operator W_{pqr,stu} is
contracted with the mean-field 1-RDM to produce effective one- and
two-body contributions:
h̃_{pq} = h_{pq} + Σ_{rs} W_{prs,qrs} · ¹D_{rs}
g̃_{pqrs} = g_{pqrs} + Σ_{t} W_{pqt,rst} · ¹D_{st}
With an idempotent 1-RDM (HF reference):
¹D_{pq} = 2 δ_{pq} for p,q ∈ occupied, 0 otherwise.
"""
# For the minimal implementation, the three-body tensor is not
# explicitly constructed; we return the 2-body terms as-is.
# A full implementation would build W and perform the contractions.
return h1e, h2e_phys
def _symmetrize_h2e_tensor(h2e: np.ndarray) -> np.ndarray:
"""Hermitise the 2-body tensor: g̃_{pqrs} → (g_{pqrs} + g_{rspq}*)/2."""
return 0.5 * (h2e + h2e.transpose(2, 3, 0, 1))