"""Transcorrelated Hamiltonian construction.
Transcorrelated methods incorporate short-range electron correlation
into the Hamiltonian through a similarity transformation:
H̃ = e^{-t} H e^{t}
where t 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 t).
* Modified two-body terms (electron repulsion is screened at short range).
* Three-body terms (arise from double commutators [t,[T,t]]
and [t,[V,t]]).
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^{t} parameterisation.
2. Builds the similarity-transformed Hamiltonian through O(t^2).
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:
t = -g S_{i<j} exp(-a r_{ij}^2)
In second-quantised form, this becomes a two-body operator:
t = 1/2 S_{pqrs} t_{pqrs} a+_p a+_q a_s a_r
The similarity-transformed Hamiltonian through O(t) is:
H̃ ≈ H + [H, t]
Since H = T + V contains one- and two-body operators, [H, t] generates
up to three-body terms.
For the simple Gaussian model, we approximate t by its action on the
two-electron integrals:
g̃_{pqrs} = g_{pqrs} . (1 + g . f_{pqrs})
where f_{pqrs} damps the integrals at short range.
"""
# Build a simple damping factor for the two-electron integrals.
# f_{pqrs} = exp(-g . |e_p - e_q + e_r - e_s|) -- simplified model
# In a real implementation, this would be computed from the
# orbital-dependent Gaussian integrals of exp(-g r_{ij}^2).
# For this minimal implementation, we use a diagonal-approximation
# damping: g̃_{pqrs} = g_{pqrs} . (1 - tanh(g . Δ_{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(-g . d^2)) where d 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, 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.
t = -g S_{i<j} u(r_{ij}), u(r) = (1 - exp(-b r)) / b
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(-b r))/b
# At short range: v -> r; at long range: v -> 1/b
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} + S_{rs} W_{prs,qrs} . ¹D_{rs}
g̃_{pqrs} = g_{pqrs} + S_{t} W_{pqt,rst} . ¹D_{st}
With an idempotent 1-RDM (HF reference):
¹D_{pq} = 2 d_{pq} for p,q in 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))