"""DF-CCSD and CCSD(T) correlation-energy driver for closed-shell RHF references.
Usage (standalone)::
from vibeqc import Molecule, BasisSet, run_rhf, RHFOptions
from vibeqc import run_ccsd, CCSDOptions
mol = Molecule(...)
basis = BasisSet(mol, "cc-pVDZ")
hf = run_rhf(mol, basis, RHFOptions())
opts = CCSDOptions()
opts.aux_basis = "cc-pvdz-ri"
opts.compute_triples = True
result = run_ccsd(mol, basis, hf, opts)
print(f"CCSD(T) = {result.e_total:.12f} Ha")
print(f"CCSD corr = {result.e_ccsd_correlation:.12f} Ha")
print(f"(T) corr = {result.e_t:.12f} Ha")
Or through ``run_job`` with ``method="ccsd(t)"``::
from vibeqc import run_job
run_job(mol, basis="cc-pVDZ", method="ccsd(t)")
"""
from __future__ import annotations
from typing import Optional
from ._vibeqc_core import CCSDOptions as _CCSDOptions
from ._vibeqc_core import CCSDResult as _CCSDResult
from ._vibeqc_core import run_ccsd as _run_ccsd
from .density_fitting import default_aux_basis_for
_OPTION_ATTRIBUTES = (
"density_fit",
"aux_basis",
"max_iter",
"conv_tol_energy",
"conv_tol_residual",
"diis_subspace_size",
"n_frozen_core",
"compute_triples",
"triples_memory_mode",
)
_TRIPLES_NONE = frozenset({"none", "off", "false", "0", "no", "ccsd"})
_TRIPLES_STANDARD_T = frozenset(
{
"(t)",
"t",
"true",
"1",
"yes",
"standard",
"standard(t)",
"raghavachari",
"raghavachari(t)",
"ccsd(t)",
}
)
_TRIPLES_UNIMPLEMENTED = {
"a-ccsd(t)": "A-CCSD(T)",
"accsd(t)": "A-CCSD(T)",
"asymmetric": "A-CCSD(T)",
"ccsd[t]": "CCSD[T]",
"[t]": "CCSD[T]",
"bracket": "CCSD[T]",
"ccsd+t(ccsd)": "CCSD+T(CCSD)",
"+t(ccsd)": "CCSD+T(CCSD)",
}
def _normalise_triples_selector(selector: object) -> str:
if isinstance(selector, bool):
return "(t)" if selector else "none"
if not isinstance(selector, str):
raise ValueError(
"triples= expects 'none', '(t)', or one of the named roadmap "
f"triples variants; got {selector!r}."
)
key = selector.strip().lower().replace(" ", "").replace("_", "-")
if key in _TRIPLES_NONE:
return "none"
if key in _TRIPLES_STANDARD_T:
return "(t)"
if key in _TRIPLES_UNIMPLEMENTED:
variant = _TRIPLES_UNIMPLEMENTED[key]
raise NotImplementedError(
f"triples={selector!r} selects {variant}, which is a v0.14 "
"roadmap variant but is not implemented yet. Currently supported "
"selectors are triples='none' and triples='(t)'."
)
raise ValueError(
"triples= expects 'none', '(t)', or one of the named roadmap triples "
f"variants; got {selector!r}."
)
def resolve_triples_selector(selector: object) -> bool:
"""Return whether a CCSD triples selector requests standard ``(T)``.
The current numerical kernel supports plain CCSD and the standard
Raghavachari perturbative triples correction. Named related
variants (``A-CCSD(T)``, ``CCSD[T]``, ``CCSD+T(CCSD)``) are parsed
so callers get a clear :class:`NotImplementedError` instead of a
silent fallback.
"""
return _normalise_triples_selector(selector) == "(t)"
[docs]
class CCSDOptions(_CCSDOptions):
"""Python-friendly CCSD option struct with aux-basis autodetection."""
[docs]
def __init__(
self,
*,
density_fit: bool = True,
aux_basis: Optional[str] = None,
max_iter: int = 100,
conv_tol_energy: float = 1e-8,
conv_tol_residual: float = 1e-7,
diis_subspace_size: int = 6,
n_frozen_core: int = 0,
compute_triples: bool = True,
triples: Optional[object] = None,
triples_memory_mode: str = "fast",
):
super().__init__()
if triples is not None:
compute_triples = resolve_triples_selector(triples)
self.density_fit = density_fit
if aux_basis is not None:
self.aux_basis = aux_basis
self.max_iter = max_iter
self.conv_tol_energy = conv_tol_energy
self.conv_tol_residual = conv_tol_residual
self.diis_subspace_size = diis_subspace_size
self.n_frozen_core = n_frozen_core
self.compute_triples = compute_triples
self.triples_memory_mode = triples_memory_mode
@property
def triples(self) -> str:
"""Human-readable triples selector: ``"none"`` or ``"(t)"``."""
return "(t)" if bool(self.compute_triples) else "none"
@triples.setter
def triples(self, selector: object) -> None:
self.compute_triples = resolve_triples_selector(selector)
[docs]
def resolve_aux_basis(self, orbital_basis_name: str) -> None:
"""If ``aux_basis`` is empty, auto-detect from the orbital basis name."""
if not self.aux_basis:
self.aux_basis = default_aux_basis_for(orbital_basis_name, kind="ri")
def _copy_ccsd_options(options: _CCSDOptions) -> CCSDOptions:
"""Copy native/Python CCSD options into the Python wrapper type."""
if not isinstance(options, _CCSDOptions):
raise TypeError(
"CCSD options must be a vibeqc.CCSDOptions or native "
f"_vibeqc_core.CCSDOptions instance, got {type(options).__name__}."
)
py_opts = CCSDOptions()
for attr in _OPTION_ATTRIBUTES:
setattr(py_opts, attr, getattr(options, attr))
return py_opts
[docs]
def run_ccsd(molecule, basis, rhf_result, options=None):
"""Run closed-shell DF-CCSD (and optionally CCSD(T)).
Parameters
----------
molecule : Molecule
basis : BasisSet
rhf_result : RHFResult
Converged RHF reference.
options : CCSDOptions, optional
CCSD options. If not provided, defaults are used with
``density_fit=True`` and ``compute_triples=True``.
Returns
-------
CCSDResult
"""
if options is None:
options = CCSDOptions()
elif isinstance(options, _CCSDOptions) and not isinstance(options, CCSDOptions):
# Upcast bare C++ options to Python wrapper for aux resolution.
options = _copy_ccsd_options(options)
if not getattr(rhf_result, "converged", False):
raise RuntimeError("run_ccsd: RHF reference is not converged")
if getattr(molecule, "multiplicity", 1) != 1:
raise ValueError(
"run_ccsd: CCSD requires a closed-shell (singlet) reference; "
f"got multiplicity={getattr(molecule, 'multiplicity', 1)}"
)
if options.density_fit:
options.resolve_aux_basis(basis.name)
return _run_ccsd(molecule, basis, rhf_result, options)
[docs]
def chemical_core_orbital_count(molecule) -> int:
"""Number of chemical-core spatial orbitals for frozen-core CCSD.
Standard noble-gas-core convention per atom: H, He contribute 0;
Li..Ne 1 (1s); Na..Ar 5 (+2s2p); K..Kr 9 (+3s3p); Rb..Xe 18
(+3d4s4p); heavier 27. Matches the usual frozen-core counts of
molecular correlation codes for main-group chemistry.
"""
n = 0
for atom in molecule.atoms:
z = atom.Z
if z <= 2:
continue
if z <= 10:
n += 1
elif z <= 18:
n += 5
elif z <= 36:
n += 9
elif z <= 54:
n += 18
else:
n += 27
return n