Unrestricted electron-boson coupled cluster.
ebcc.cc.uebcc.UEBCC(mf, log=None, ansatz='CCSD', options=None, space=None, omega=None, g=None, G=None, mo_coeff=None, mo_occ=None, fock=None, **kwargs)
Bases: BaseEBCC
Unrestricted electron-boson coupled cluster.
Initialise the EBCC object.
| Parameters: |
|
|---|
Source code in ebcc/cc/base.py
def __init__(
self,
mf: SCF,
log: Optional[Logger] = None,
ansatz: Optional[Union[Ansatz, str]] = "CCSD",
options: Optional[BaseOptions] = None,
space: Optional[SpaceType] = None,
omega: Optional[NDArray[T]] = None,
g: Optional[NDArray[T]] = None,
G: Optional[NDArray[T]] = None,
mo_coeff: Optional[NDArray[T]] = None,
mo_occ: Optional[NDArray[T]] = None,
fock: Optional[BaseFock] = None,
**kwargs: Any,
) -> None:
r"""Initialise the EBCC object.
Args:
mf: PySCF mean-field object.
log: Log to write output to. Default is the global logger, outputting to `stderr`.
ansatz: Overall ansatz.
options: Options for the EBCC calculation.
space: Space containing the frozen, correlated, and active fermionic spaces. Default
assumes all electrons are correlated.
omega: Bosonic frequencies.
g: Electron-boson coupling matrix corresponding to the bosonic annihilation operator
:math:`g_{bpq} p^\dagger q b`. The creation part is assumed to be the fermionic
conjugate transpose to retain Hermiticity in the Hamiltonian.
G: Boson non-conserving term :math:`G_{b} (b^\dagger + b)`.
mo_coeff: Molecular orbital coefficients. Default is the mean-field coefficients.
mo_occ: Molecular orbital occupation numbers. Default is the mean-field occupation.
fock: Fock matrix. Default is the mean-field Fock matrix.
**kwargs: Additional keyword arguments used to update `options`.
"""
# Options:
if options is None:
options = self.Options()
self.options = options
for key, val in kwargs.items():
setattr(self.options, key, val)
# Parameters:
self.log = default_log if log is None else log
self.mf = self._convert_mf(mf)
self._mo_coeff: Optional[NDArray[T]] = (
np.asarray(mo_coeff, dtype=types[float]) if mo_coeff is not None else None
)
self._mo_occ: Optional[NDArray[T]] = (
np.asarray(mo_occ, dtype=types[float]) if mo_occ is not None else None
)
# Ansatz:
if isinstance(ansatz, Ansatz):
self.ansatz = ansatz
elif isinstance(ansatz, str):
self.ansatz = Ansatz.from_string(
ansatz, density_fitting=getattr(self.mf, "with_df", None) is not None
)
else:
raise TypeError("ansatz must be an Ansatz object or a string.")
self._eqns = self.ansatz._get_eqns(self.spin_type)
# Space:
if space is not None:
self.space = space
else:
self.space = self.init_space()
# Boson parameters:
if bool(self.fermion_coupling_rank) != bool(self.boson_coupling_rank):
raise ValueError(
"Fermionic and bosonic coupling ranks must both be zero, or both non-zero."
)
self.omega = np.asarray(omega, dtype=types[float]) if omega is not None else None
self.bare_g = np.asarray(g, dtype=types[float]) if g is not None else None
self.bare_G = np.asarray(G, dtype=types[float]) if G is not None else None
if self.boson_ansatz != "":
self.g = self.get_g()
self.G = self.get_mean_field_G()
else:
assert self.nbos == 0
self.options.shift = False
self.g = None
self.G = None
# Fock matrix:
if fock is None:
self.fock = self.get_fock()
else:
self.fock = fock
# Attributes:
self.e_corr = 0.0
self.amplitudes = util.Namespace()
self.converged = False
self.lambdas = util.Namespace()
self.converged_lambda = False
# Logging:
init_logging(self.log)
self.log.info(f"\n{ANSI.B}{ANSI.U}{self.name}{ANSI.R}")
self.log.debug(f"{ANSI.B}{'*' * len(self.name)}{ANSI.R}")
self.log.debug("")
self.log.info(f"{ANSI.B}Options{ANSI.R}:")
self.log.info(f" > e_tol: {ANSI.y}{self.options.e_tol}{ANSI.R}")
self.log.info(f" > t_tol: {ANSI.y}{self.options.t_tol}{ANSI.R}")
self.log.info(f" > max_iter: {ANSI.y}{self.options.max_iter}{ANSI.R}")
self.log.info(f" > diis_space: {ANSI.y}{self.options.diis_space}{ANSI.R}")
self.log.info(f" > diis_min_space: {ANSI.y}{self.options.diis_min_space}{ANSI.R}")
self.log.info(f" > damping: {ANSI.y}{self.options.damping}{ANSI.R}")
self.log.debug("")
self.log.info(f"{ANSI.B}Ansatz{ANSI.R}: {ANSI.m}{self.ansatz}{ANSI.R}")
self.log.debug("")
self.log.info(f"{ANSI.B}Space{ANSI.R}: {ANSI.m}{self.space}{ANSI.R}")
self.log.debug("")
if self.boson_ansatz != "":
self.log.info(f"{ANSI.B}Bosons{ANSI.R}: {ANSI.m}{self.nbos}{ANSI.R}")
self.log.info(" > Energy shift due to polaritonic basis: %.10f", self.const)
self.log.debug("")
ebcc.cc.uebcc.UEBCC.spin_type: str
property
Get a string representation of the spin type.
ebcc.cc.uebcc.UEBCC.xi: NDArray[T]
property
Get the shift in the bosonic operators to diagonalise the photon Hamiltonian.
| Returns: |
|
|---|
ebcc.cc.uebcc.UEBCC.nmo: int
property
Get the number of molecular orbitals.
| Returns: |
|
|---|
ebcc.cc.uebcc.UEBCC.nocc: tuple[int, int]
property
Get the number of occupied molecular orbitals.
| Returns: |
|
|---|
ebcc.cc.uebcc.UEBCC.nvir: tuple[int, int]
property
Get the number of virtual molecular orbitals.
| Returns: |
|
|---|
ebcc.cc.uebcc.UEBCC.ip_eom(**kwargs)
Get the IP-EOM object.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def ip_eom(self, **kwargs: Any) -> IP_UEOM:
"""Get the IP-EOM object.
Args:
options: Options for the IP-EOM calculation.
**kwargs: Additional keyword arguments.
Returns:
IP-EOM object.
"""
return IP_UEOM(self, **kwargs)
ebcc.cc.uebcc.UEBCC.ea_eom(**kwargs)
Get the EA-EOM object.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def ea_eom(self, **kwargs: Any) -> EA_UEOM:
"""Get the EA-EOM object.
Args:
options: Options for the EA-EOM calculation.
**kwargs: Additional keyword arguments.
Returns:
EA-EOM object.
"""
return EA_UEOM(self, **kwargs)
ebcc.cc.uebcc.UEBCC.ee_eom(**kwargs)
Get the EE-EOM object.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def ee_eom(self, **kwargs: Any) -> EE_UEOM:
"""Get the EE-EOM object.
Args:
options: Options for the EE-EOM calculation.
**kwargs: Additional keyword arguments.
Returns:
EE-EOM object.
"""
return EE_UEOM(self, **kwargs)
ebcc.cc.uebcc.UEBCC.from_rebcc(rcc)
classmethod
Initialise an UEBCC object from an REBCC object.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
@classmethod
def from_rebcc(cls, rcc: REBCC) -> UEBCC:
"""Initialise an `UEBCC` object from an `REBCC` object.
Args:
rcc: Restricted electron-boson coupled cluster object.
Returns:
UEBCC object.
"""
ucc = cls(
rcc.mf,
log=rcc.log,
ansatz=rcc.ansatz,
space=(rcc.space, rcc.space),
omega=rcc.omega,
g=rcc.bare_g,
G=rcc.bare_G,
options=rcc.options,
)
ucc.e_corr = rcc.e_corr
ucc.converged = rcc.converged
ucc.converged_lambda = rcc.converged_lambda
has_amps = bool(rcc.amplitudes)
has_lams = bool(rcc.lambdas)
if has_amps:
amplitudes: Namespace[SpinArrayType] = util.Namespace()
for name, key, n in ucc.ansatz.fermionic_cluster_ranks(spin_type=ucc.spin_type):
amplitudes[name] = util.Namespace()
for comb in util.generate_spin_combinations(n, unique=True):
subscript, _ = util.combine_subscripts(key, comb)
tn = rcc.amplitudes[name]
tn = util.symmetrise(subscript, tn, symmetry="-" * 2 * n)
amplitudes[name][comb] = tn
for name, key, n in ucc.ansatz.bosonic_cluster_ranks(spin_type=ucc.spin_type):
amplitudes[name] = np.copy(rcc.amplitudes[name]) # type: ignore
for name, key, nf, nb in ucc.ansatz.coupling_cluster_ranks(spin_type=ucc.spin_type):
amplitudes[name] = util.Namespace()
for comb in util.generate_spin_combinations(nf, unique=True):
tn = rcc.amplitudes[name]
amplitudes[name][comb] = tn
ucc.amplitudes = amplitudes
if has_lams:
lambdas: Namespace[SpinArrayType] = util.Namespace()
for name, key, n in ucc.ansatz.fermionic_cluster_ranks(spin_type=ucc.spin_type):
lname = name.replace("t", "l")
lambdas[lname] = util.Namespace()
for comb in util.generate_spin_combinations(n, unique=True):
subscript, _ = util.combine_subscripts(key, comb)
tn = rcc.lambdas[lname]
tn = util.symmetrise(subscript, tn, symmetry="-" * 2 * n)
lambdas[lname][comb] = tn
for name, key, n in ucc.ansatz.bosonic_cluster_ranks(spin_type=ucc.spin_type):
lname = "l" + name
lambdas[lname] = np.copy(rcc.lambdas[lname]) # type: ignore
for name, key, nf, nb in ucc.ansatz.coupling_cluster_ranks(spin_type=ucc.spin_type):
lname = "l" + name
lambdas[lname] = util.Namespace()
for comb in util.generate_spin_combinations(nf, unique=True):
tn = rcc.lambdas[lname]
lambdas[lname][comb] = tn
ucc.lambdas = lambdas
return ucc
ebcc.cc.uebcc.UEBCC.init_space()
Initialise the fermionic space.
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def init_space(self) -> SpaceType:
"""Initialise the fermionic space.
Returns:
Fermionic space. All fermionic degrees of freedom are assumed to be correlated.
"""
space = (
Space(
self.mo_occ[0] > 0,
np.zeros(self.mo_occ[0].shape, dtype=np.bool_),
np.zeros(self.mo_occ[0].shape, dtype=np.bool_),
),
Space(
self.mo_occ[1] > 0,
np.zeros(self.mo_occ[1].shape, dtype=np.bool_),
np.zeros(self.mo_occ[1].shape, dtype=np.bool_),
),
)
return space
ebcc.cc.uebcc.UEBCC.init_amps(eris=None)
Initialise the cluster amplitudes.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def init_amps(self, eris: Optional[ERIsInputType] = None) -> Namespace[SpinArrayType]:
"""Initialise the cluster amplitudes.
Args:
eris: Electron repulsion integrals.
Returns:
Initial cluster amplitudes.
"""
eris = self.get_eris(eris)
amplitudes: Namespace[SpinArrayType] = util.Namespace()
# Build T amplitudes
for name, key, n in self.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type):
tn: SpinArrayType = util.Namespace()
for comb in util.generate_spin_combinations(n, unique=True):
if n == 1:
tn[comb] = self.fock[comb][key] / self.energy_sum(key, comb)
elif n == 2:
comb_t = comb[0] + comb[2] + comb[1] + comb[3]
key_t = key[0] + key[2] + key[1] + key[3]
tn[comb] = np.transpose(eris[comb_t][key_t], (0, 2, 1, 3)) / self.energy_sum(
key, comb
)
if comb in ("aaaa", "bbbb"):
# TODO generalise:
tn[comb] = (tn[comb] - np.transpose(tn[comb], (1, 0, 2, 3))) * 0.5
else:
shape = tuple(self.space["ab".index(s)].size(k) for s, k in zip(comb, key))
tn[comb] = np.zeros(shape, dtype=types[float])
amplitudes[name] = tn
# Build S amplitudes:
for name, key, n in self.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type):
if self.omega is None or self.G is None:
raise ValueError("Bosonic parameters not set.")
if n == 1:
amplitudes[name] = -self.G / self.omega # type: ignore
else:
shape = (self.nbos,) * n
amplitudes[name] = np.zeros(shape, dtype=types[float]) # type: ignore
# Build U amplitudes:
for name, key, nf, nb in self.ansatz.coupling_cluster_ranks(spin_type=self.spin_type):
if self.omega is None or self.g is None:
raise ValueError("Bosonic parameters not set.")
if nf != 1:
raise util.ModelNotImplemented
if nb == 1:
tn = util.Namespace(
aa=self.g.aa[key] / self.energy_sum(key, "_aa"),
bb=self.g.bb[key] / self.energy_sum(key, "_aa"),
)
amplitudes[name] = tn
else:
tn = util.Namespace(
aa=np.zeros(
(self.nbos,) * nb + tuple(self.space[0].size(k) for k in key[nb:]),
dtype=types[float],
),
bb=np.zeros(
(self.nbos,) * nb + tuple(self.space[0].size(k) for k in key[nb:]),
dtype=types[float],
),
)
amplitudes[name] = tn
return amplitudes
ebcc.cc.uebcc.UEBCC.init_lams(amplitudes=None)
Initialise the cluster lambda amplitudes.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def init_lams(
self, amplitudes: Optional[Namespace[SpinArrayType]] = None
) -> Namespace[SpinArrayType]:
"""Initialise the cluster lambda amplitudes.
Args:
amplitudes: Cluster amplitudes.
Returns:
Initial cluster lambda amplitudes.
"""
if not amplitudes:
amplitudes = self.amplitudes
lambdas: Namespace[SpinArrayType] = util.Namespace()
# Build L amplitudes:
for name, key, n in self.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type):
lname = name.replace("t", "l")
perm = list(range(n, 2 * n)) + list(range(n))
lambdas[lname] = util.Namespace()
for key in dict(amplitudes[name]).keys():
ln = np.transpose(amplitudes[name][key], perm)
lambdas[lname][key] = ln
# Build LS amplitudes:
for name, key, n in self.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type):
lambdas["l" + name] = amplitudes[name] # type: ignore
# Build LU amplitudes:
for name, key, nf, nb in self.ansatz.coupling_cluster_ranks(spin_type=self.spin_type):
if nf != 1:
raise util.ModelNotImplemented
perm = list(range(nb)) + [nb + 1, nb]
lambdas["l" + name] = util.Namespace()
for key in dict(amplitudes[name]).keys():
ln = np.transpose(amplitudes[name][key], perm)
lambdas["l" + name][key] = ln
return lambdas
ebcc.cc.uebcc.UEBCC.update_amps(eris=None, amplitudes=None)
Update the cluster amplitudes.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def update_amps(
self,
eris: Optional[ERIsInputType] = None,
amplitudes: Optional[Namespace[SpinArrayType]] = None,
) -> Namespace[SpinArrayType]:
"""Update the cluster amplitudes.
Args:
eris: Electron repulsion integrals.
amplitudes: Cluster amplitudes.
Returns:
Updated cluster amplitudes.
"""
amplitudes = self._get_amps(amplitudes=amplitudes)
func, kwargs = self._load_function(
"update_amps",
eris=eris,
amplitudes=amplitudes,
)
res: Namespace[SpinArrayType] = func(**kwargs)
res = util.Namespace(**{key.rstrip("new"): val for key, val in res.items()})
# Divide T amplitudes:
for name, key, n in self.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type):
for comb in util.generate_spin_combinations(n, unique=True):
subscript, _ = util.combine_subscripts(key, comb)
tn = res[name][comb]
tn /= self.energy_sum(key, comb)
tn += amplitudes[name][comb]
tn = util.symmetrise(subscript, tn, symmetry="-" * (2 * n))
res[name][comb] = tn
# Divide S amplitudes:
for name, key, n in self.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type):
res[name] /= self.energy_sum(key, "_" * n) # type: ignore
res[name] += amplitudes[name] # type: ignore
# Divide U amplitudes:
for name, key, nf, nb in self.ansatz.coupling_cluster_ranks(spin_type=self.spin_type):
if nf != 1:
raise util.ModelNotImplemented
tn = res[name].aa
tn /= self.energy_sum(key, "_" * nb + "aa")
tn += amplitudes[name].aa
res[name].aa = tn
tn = res[name].bb
tn /= self.energy_sum(key, "_" * nb + "bb")
tn += amplitudes[name].bb
res[name].bb = tn
return res
ebcc.cc.uebcc.UEBCC.update_lams(eris=None, amplitudes=None, lambdas=None, lambdas_pert=None, perturbative=False)
Update the cluster lambda amplitudes.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def update_lams(
self,
eris: Optional[ERIsInputType] = None,
amplitudes: Optional[Namespace[SpinArrayType]] = None,
lambdas: Optional[Namespace[SpinArrayType]] = None,
lambdas_pert: Optional[Namespace[SpinArrayType]] = None,
perturbative: bool = False,
) -> Namespace[SpinArrayType]:
"""Update the cluster lambda amplitudes.
Args:
eris: Electron repulsion integrals.
amplitudes: Cluster amplitudes.
lambdas: Cluster lambda amplitudes.
lambdas_pert: Perturbative cluster lambda amplitudes.
perturbative: Flag to include perturbative correction.
Returns:
Updated cluster lambda amplitudes.
"""
# TODO active
amplitudes = self._get_amps(amplitudes=amplitudes)
lambdas = self._get_lams(lambdas=lambdas, amplitudes=amplitudes)
func, kwargs = self._load_function(
"update_lams",
eris=eris,
amplitudes=amplitudes,
lambdas=lambdas,
lambdas_pert=lambdas_pert,
)
res: Namespace[SpinArrayType] = func(**kwargs)
res = util.Namespace(**{key.rstrip("new"): val for key, val in res.items()})
# Divide T amplitudes:
for name, key, n in self.ansatz.fermionic_cluster_ranks(
spin_type=self.spin_type, which="l"
):
for comb in util.generate_spin_combinations(n, unique=True):
subscript, _ = util.combine_subscripts(key, comb)
tn = res[name][comb]
tn /= self.energy_sum(key, comb)
tn += lambdas[name][comb]
tn = util.symmetrise(subscript, tn, symmetry="-" * (2 * n))
res[name][comb] = tn
# Divide S amplitudes:
for name, key, n in self.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type, which="l"):
res[name] /= self.energy_sum(key, "_" * n) # type: ignore
res[name] += lambdas[name] # type: ignore
# Divide U amplitudes:
for name, key, nf, nb in self.ansatz.coupling_cluster_ranks(
spin_type=self.spin_type, which="l"
):
if nf != 1:
raise util.ModelNotImplemented
tn = res[name].aa
tn /= self.energy_sum(key, "_" * nb + "aa")
tn += lambdas[name].aa
res[name].aa = tn
tn = res[name].bb
tn /= self.energy_sum(key, "_" * nb + "bb")
tn += lambdas[name].bb
res[name].bb = tn
return res
ebcc.cc.uebcc.UEBCC.make_rdm1_f(eris=None, amplitudes=None, lambdas=None, hermitise=True)
Make the one-particle fermionic reduced density matrix :math:\langle i^+ j \rangle.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def make_rdm1_f(
self,
eris: Optional[ERIsInputType] = None,
amplitudes: Optional[Namespace[SpinArrayType]] = None,
lambdas: Optional[Namespace[SpinArrayType]] = None,
hermitise: bool = True,
) -> SpinArrayType:
r"""Make the one-particle fermionic reduced density matrix :math:`\langle i^+ j \rangle`.
Args:
eris: Electron repulsion integrals.
amplitudes: Cluster amplitudes.
lambdas: Cluster lambda amplitudes.
hermitise: Hermitise the density matrix.
Returns:
One-particle fermion reduced density matrix.
"""
func, kwargs = self._load_function(
"make_rdm1_f",
eris=eris,
amplitudes=amplitudes,
lambdas=lambdas,
)
dm: SpinArrayType = func(**kwargs)
if hermitise:
dm.aa = (dm.aa + np.transpose(dm.aa)) * 0.5
dm.bb = (dm.bb + np.transpose(dm.bb)) * 0.5
return dm
ebcc.cc.uebcc.UEBCC.make_rdm2_f(eris=None, amplitudes=None, lambdas=None, hermitise=True)
Make the two-particle fermionic reduced density matrix :math:\langle i^+j^+lk \rangle.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def make_rdm2_f(
self,
eris: Optional[ERIsInputType] = None,
amplitudes: Optional[Namespace[SpinArrayType]] = None,
lambdas: Optional[Namespace[SpinArrayType]] = None,
hermitise: bool = True,
) -> SpinArrayType:
r"""Make the two-particle fermionic reduced density matrix :math:`\langle i^+j^+lk \rangle`.
Args:
eris: Electron repulsion integrals.
amplitudes: Cluster amplitudes.
lambdas: Cluster lambda amplitudes.
hermitise: Hermitise the density matrix.
Returns:
Two-particle fermion reduced density matrix.
"""
func, kwargs = self._load_function(
"make_rdm2_f",
eris=eris,
amplitudes=amplitudes,
lambdas=lambdas,
)
dm: SpinArrayType = func(**kwargs)
if hermitise:
def transpose1(dm: NDArray[T]) -> NDArray[T]:
dm = (np.transpose(dm, (0, 1, 2, 3)) + np.transpose(dm, (2, 3, 0, 1))) * 0.5
return dm
def transpose2(dm: NDArray[T]) -> NDArray[T]:
dm = (np.transpose(dm, (0, 1, 2, 3)) + np.transpose(dm, (1, 0, 3, 2))) * 0.5
return dm
dm.aaaa = transpose2(transpose1(dm.aaaa))
dm.aabb = transpose2(dm.aabb)
dm.bbbb = transpose2(transpose1(dm.bbbb))
return dm
ebcc.cc.uebcc.UEBCC.make_eb_coup_rdm(eris=None, amplitudes=None, lambdas=None, unshifted=True, hermitise=True)
Make the electron-boson coupling reduced density matrix.
.. math:: \langle b^+ i^+ j \rangle
and
.. math:: \langle b i^+ j \rangle
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def make_eb_coup_rdm(
self,
eris: Optional[ERIsInputType] = None,
amplitudes: Optional[Namespace[SpinArrayType]] = None,
lambdas: Optional[Namespace[SpinArrayType]] = None,
unshifted: bool = True,
hermitise: bool = True,
) -> SpinArrayType:
r"""Make the electron-boson coupling reduced density matrix.
.. math::
\langle b^+ i^+ j \rangle
and
.. math::
\langle b i^+ j \rangle
Args:
eris: Electron repulsion integrals.
amplitudes: Cluster amplitudes.
lambdas: Cluster lambda amplitudes.
unshifted: If `self.options.shift` is `True`, return the unshifted density matrix. Has
no effect if `self.options.shift` is `False`.
hermitise: Hermitise the density matrix.
Returns:
Electron-boson coupling reduced density matrix.
"""
func, kwargs = self._load_function(
"make_eb_coup_rdm",
eris=eris,
amplitudes=amplitudes,
lambdas=lambdas,
)
dm_eb: SpinArrayType = func(**kwargs)
if hermitise:
dm_eb.aa = np.array(
[
(dm_eb.aa[0] + np.transpose(dm_eb.aa[1], (0, 2, 1))) * 0.5,
(dm_eb.aa[1] + np.transpose(dm_eb.aa[0], (0, 2, 1))) * 0.5,
]
)
dm_eb.bb = np.array(
[
(dm_eb.bb[0] + np.transpose(dm_eb.bb[1], (0, 2, 1))) * 0.5,
(dm_eb.bb[1] + np.transpose(dm_eb.bb[0], (0, 2, 1))) * 0.5,
]
)
if unshifted and self.options.shift:
rdm1_f = self.make_rdm1_f(hermitise=hermitise)
shift = util.einsum("x,ij->xij", self.xi, rdm1_f.aa)
dm_eb.aa -= shift[None]
shift = util.einsum("x,ij->xij", self.xi, rdm1_f.bb)
dm_eb.bb -= shift[None]
return dm_eb
ebcc.cc.uebcc.UEBCC.energy_sum(*args, signs_dict=None)
Get a direct sum of energies.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def energy_sum(self, *args: str, signs_dict: Optional[dict[str, str]] = None) -> NDArray[T]:
"""Get a direct sum of energies.
Args:
*args: Energies to sum. Should specify a subscript and spins.
signs_dict: Signs of the energies in the sum. Default sets `("o", "O", "i")` to be
positive, and `("v", "V", "a", "b")` to be negative.
Returns:
Sum of energies.
"""
subscript, spins = args
n = 0
def next_char() -> str:
nonlocal n
if n < 26:
char = chr(ord("a") + n)
else:
char = chr(ord("A") + n)
n += 1
return char
if signs_dict is None:
signs_dict = {}
for k, s in zip("vVaoOib", "---+++-"):
if k not in signs_dict:
signs_dict[k] = s
energies = []
for key, spin in zip(subscript, spins):
factor = 1 if signs_dict[key] == "+" else -1
if key == "b":
assert self.omega is not None
energies.append(self.omega * types[float](factor))
else:
energies.append(np.diag(self.fock[spin + spin][key + key]) * types[float](factor))
subscript = ",".join([next_char() for k in subscript])
energy_sum = util.dirsum(subscript, *energies)
return energy_sum
ebcc.cc.uebcc.UEBCC.amplitudes_to_vector(amplitudes)
Construct a vector containing all of the amplitudes used in the given ansatz.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def amplitudes_to_vector(self, amplitudes: Namespace[SpinArrayType]) -> NDArray[T]:
"""Construct a vector containing all of the amplitudes used in the given ansatz.
Args:
amplitudes: Cluster amplitudes.
Returns:
Cluster amplitudes as a vector.
"""
vectors = []
for name, key, n in self.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type):
for spin in util.generate_spin_combinations(n, unique=True):
tn = amplitudes[name][spin]
subscript, _ = util.combine_subscripts(key, spin)
vectors.append(np.ravel(util.compress_axes(subscript, tn)))
for name, key, n in self.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type):
vectors.append(np.ravel(amplitudes[name])) # type: ignore
for name, key, nf, nb in self.ansatz.coupling_cluster_ranks(spin_type=self.spin_type):
if nf != 1:
raise util.ModelNotImplemented
vectors.append(np.ravel(amplitudes[name].aa))
vectors.append(np.ravel(amplitudes[name].bb))
return np.concatenate(vectors)
ebcc.cc.uebcc.UEBCC.vector_to_amplitudes(vector)
Construct a namespace of amplitudes from a vector.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def vector_to_amplitudes(self, vector: NDArray[T]) -> Namespace[SpinArrayType]:
"""Construct a namespace of amplitudes from a vector.
Args:
vector: Cluster amplitudes as a vector.
Returns:
Cluster amplitudes.
"""
amplitudes: Namespace[SpinArrayType] = util.Namespace()
i0 = 0
sizes: dict[tuple[str, ...], int] = {
(o, s): self.space[i].size(o) for o in "ovOVia" for i, s in enumerate("ab")
}
for name, key, n in self.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type):
amplitudes[name] = util.Namespace()
for spin in util.generate_spin_combinations(n, unique=True):
subscript, csizes = util.combine_subscripts(key, spin, sizes=sizes)
size = util.get_compressed_size(subscript, **csizes)
shape = tuple(self.space["ab".index(s)].size(k) for s, k in zip(spin, key))
tn_tril = vector[i0 : i0 + size]
tn = util.decompress_axes(subscript, tn_tril, shape=shape)
amplitudes[name][spin] = tn
i0 += size
for name, key, n in self.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type):
shape = (self.nbos,) * n
size = self.nbos**n
amplitudes[name] = np.reshape(vector[i0 : i0 + size], shape) # type: ignore
i0 += size
for name, key, nf, nb in self.ansatz.coupling_cluster_ranks(spin_type=self.spin_type):
if nf != 1:
raise util.ModelNotImplemented
amplitudes[name] = util.Namespace()
shape = (self.nbos,) * nb + tuple(self.space[0].size(k) for k in key[nb:])
size = util.prod(shape)
amplitudes[name].aa = np.reshape(vector[i0 : i0 + size], shape)
i0 += size
shape = (self.nbos,) * nb + tuple(self.space[1].size(k) for k in key[nb:])
size = util.prod(shape)
amplitudes[name].bb = np.reshape(vector[i0 : i0 + size], shape)
i0 += size
assert i0 == len(vector)
return amplitudes
ebcc.cc.uebcc.UEBCC.lambdas_to_vector(lambdas)
Construct a vector containing all of the lambda amplitudes used in the given ansatz.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def lambdas_to_vector(self, lambdas: Namespace[SpinArrayType]) -> NDArray[T]:
"""Construct a vector containing all of the lambda amplitudes used in the given ansatz.
Args:
lambdas: Cluster lambda amplitudes.
Returns:
Cluster lambda amplitudes as a vector.
"""
vectors = []
for name, key, n in self.ansatz.fermionic_cluster_ranks(
spin_type=self.spin_type, which="l"
):
for spin in util.generate_spin_combinations(n, unique=True):
tn = lambdas[name][spin]
subscript, _ = util.combine_subscripts(key, spin)
vectors.append(np.ravel(util.compress_axes(subscript, tn)))
for name, key, n in self.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type, which="l"):
vectors.append(np.ravel(lambdas[name])) # type: ignore
for name, key, nf, nb in self.ansatz.coupling_cluster_ranks(
spin_type=self.spin_type, which="l"
):
if nf != 1:
raise util.ModelNotImplemented
vectors.append(np.ravel(lambdas[name].aa))
vectors.append(np.ravel(lambdas[name].bb))
return np.concatenate(vectors)
ebcc.cc.uebcc.UEBCC.vector_to_lambdas(vector)
Construct a namespace of lambda amplitudes from a vector.
| Parameters: |
|
|---|
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def vector_to_lambdas(self, vector: NDArray[T]) -> Namespace[SpinArrayType]:
"""Construct a namespace of lambda amplitudes from a vector.
Args:
vector: Cluster lambda amplitudes as a vector.
Returns:
Cluster lambda amplitudes.
"""
lambdas: Namespace[SpinArrayType] = util.Namespace()
i0 = 0
sizes: dict[tuple[str, ...], int] = {
(o, s): self.space[i].size(o) for o in "ovOVia" for i, s in enumerate("ab")
}
for name, key, n in self.ansatz.fermionic_cluster_ranks(
spin_type=self.spin_type, which="l"
):
lambdas[name] = util.Namespace()
for spin in util.generate_spin_combinations(n, unique=True):
subscript, csizes = util.combine_subscripts(key, spin, sizes=sizes)
size = util.get_compressed_size(subscript, **csizes)
shape = tuple(self.space["ab".index(s)].size(k) for s, k in zip(spin, key))
tn_tril = vector[i0 : i0 + size]
tn = util.decompress_axes(subscript, tn_tril, shape=shape)
lambdas[name][spin] = tn
i0 += size
for name, key, n in self.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type, which="l"):
shape = (self.nbos,) * n
size = self.nbos**n
lambdas[name] = np.reshape(vector[i0 : i0 + size], shape) # type: ignore
i0 += size
for name, key, nf, nb in self.ansatz.coupling_cluster_ranks(
spin_type=self.spin_type, which="l"
):
if nf != 1:
raise util.ModelNotImplemented
lambdas[name] = util.Namespace()
shape = (self.nbos,) * nb + tuple(self.space[0].size(k) for k in key[nb:])
size = util.prod(shape)
lambdas[name].aa = np.reshape(vector[i0 : i0 + size], shape)
i0 += size
shape = (self.nbos,) * nb + tuple(self.space[1].size(k) for k in key[nb:])
size = util.prod(shape)
lambdas[name].bb = np.reshape(vector[i0 : i0 + size], shape)
i0 += size
assert i0 == len(vector)
return lambdas
ebcc.cc.uebcc.UEBCC.get_mean_field_G()
Get the mean-field boson non-conserving term.
| Returns: |
|
|---|
Source code in ebcc/cc/uebcc.py
def get_mean_field_G(self) -> NDArray[T]:
"""Get the mean-field boson non-conserving term.
Returns:
Mean-field boson non-conserving term.
"""
# FIXME should this also sum in frozen orbitals?
assert self.omega is not None
assert self.g is not None
boo: tuple[NDArray[T], NDArray[T]] = (self.g.aa.boo, self.g.bb.boo)
val = util.einsum("Ipp->I", boo[0])
val += util.einsum("Ipp->I", boo[1])
val -= self.xi * self.omega
if self.bare_G is not None:
# Require bare_G to have a spin index for now:
assert self.bare_G.shape == val.shape
val += self.bare_G
return val