diff --git a/ebcc/cc/base.py b/ebcc/cc/base.py index f01488b7..70c8d49e 100644 --- a/ebcc/cc/base.py +++ b/ebcc/cc/base.py @@ -30,7 +30,7 @@ from ebcc.ext.eccc import BaseExternalCorrection from ebcc.ext.tcc import BaseTailor from ebcc.ham.base import BaseElectronBoson, BaseFock - from ebcc.opt.base import BaseBruecknerEBCC + from ebcc.opt.base import BaseBruecknerEBCC, BaseOptimisedEBCC from ebcc.util import Namespace T = floating @@ -81,6 +81,7 @@ class BaseEBCC(ABC): CDERIs: type[BaseERIs] ElectronBoson: type[BaseElectronBoson] Brueckner: type[BaseBruecknerEBCC] + Optimised: type[BaseOptimisedEBCC] ExternalCorrection: type[BaseExternalCorrection] Tailor: type[BaseTailor] @@ -294,11 +295,11 @@ def kernel(self, eris: Optional[ERIsInputType] = None) -> float: converged = converged_e and converged_t if converged: self.log.debug("") - self.log.output(f"{ANSI.g}Converged.{ANSI.R}") + self.log.output(f"{ANSI.g}Converged{ANSI.R}.") break else: self.log.debug("") - self.log.warning(f"{ANSI.r}Failed to converge.{ANSI.R}") + self.log.warning(f"{ANSI.r}Failed to converge{ANSI.R}.") # Include perturbative correction if required: if self.ansatz.has_perturbative_correction: @@ -385,11 +386,11 @@ def solve_lambda( # Check for convergence: if converged: self.log.debug("") - self.log.output(f"{ANSI.g}Converged.{ANSI.R}") + self.log.output(f"{ANSI.g}Converged{ANSI.R}.") break else: self.log.debug("") - self.log.warning(f"{ANSI.r}Failed to converge.{ANSI.R}") + self.log.warning(f"{ANSI.r}Failed to converge{ANSI.R}.") self.log.debug("") self.log.debug("Time elapsed: %s", timer.format_time(timer())) @@ -451,6 +452,21 @@ def brueckner(self, *args: Any, **kwargs: Any) -> float: bcc = self.Brueckner(self, *args, **kwargs) return bcc.kernel() + def optimised(self, *args: Any, **kwargs: Any) -> float: + """Run an optimised coupled cluster calculation. + + The coupled cluster object will be updated in-place. + + Args: + *args: Arguments to pass to the optimised object. + **kwargs: Keyword arguments to pass to the optimised object. + + Returns: + Correlation energy. + """ + opt = self.Optimised(self, *args, **kwargs) + return opt.kernel() + def external_correction(self, *args: Any, **kwargs: Any) -> float: """Run an externally corrected coupled cluster calculation. @@ -625,7 +641,7 @@ def energy( eris=eris, amplitudes=amplitudes, ) - res: float = ensure_scalar(func(**kwargs)).real + res: float = np.real(ensure_scalar(func(**kwargs))) return astype(res, float) def energy_perturbative( @@ -650,7 +666,7 @@ def energy_perturbative( amplitudes=amplitudes, lambdas=lambdas, ) - res: float = ensure_scalar(func(**kwargs)).real + res: float = np.real(ensure_scalar(func(**kwargs))) return res @abstractmethod diff --git a/ebcc/cc/gebcc.py b/ebcc/cc/gebcc.py index 42f4ac6d..10771db8 100644 --- a/ebcc/cc/gebcc.py +++ b/ebcc/cc/gebcc.py @@ -18,7 +18,7 @@ from ebcc.ham.eris import GERIs from ebcc.ham.fock import GFock from ebcc.ham.space import Space -from ebcc.opt.gbrueckner import BruecknerGEBCC +from ebcc.opt.gopt import BruecknerGEBCC, OptimisedGEBCC if TYPE_CHECKING: from typing import Any, Optional, TypeAlias, Union @@ -47,6 +47,7 @@ class GEBCC(BaseEBCC): Fock = GFock ElectronBoson = GElectronBoson Brueckner = BruecknerGEBCC + Optimised = OptimisedGEBCC ExternalCorrection = ExternalCorrectionGEBCC Tailor = TailorGEBCC diff --git a/ebcc/cc/rebcc.py b/ebcc/cc/rebcc.py index 4563368e..9e841093 100644 --- a/ebcc/cc/rebcc.py +++ b/ebcc/cc/rebcc.py @@ -16,7 +16,7 @@ from ebcc.ham.eris import RERIs from ebcc.ham.fock import RFock from ebcc.ham.space import Space -from ebcc.opt.rbrueckner import BruecknerREBCC +from ebcc.opt.ropt import BruecknerREBCC, OptimisedREBCC if TYPE_CHECKING: from typing import Any, Optional, TypeAlias, Union @@ -43,6 +43,7 @@ class REBCC(BaseEBCC): CDERIs = RCDERIs ElectronBoson = RElectronBoson Brueckner = BruecknerREBCC + Optimised = OptimisedREBCC ExternalCorrection = ExternalCorrectionREBCC Tailor = TailorREBCC diff --git a/ebcc/cc/uebcc.py b/ebcc/cc/uebcc.py index 99383c0f..5e138c10 100644 --- a/ebcc/cc/uebcc.py +++ b/ebcc/cc/uebcc.py @@ -16,7 +16,7 @@ from ebcc.ham.eris import UERIs from ebcc.ham.fock import UFock from ebcc.ham.space import Space -from ebcc.opt.ubrueckner import BruecknerUEBCC +from ebcc.opt.uopt import BruecknerUEBCC, OptimisedUEBCC if TYPE_CHECKING: from typing import Any, Optional, TypeAlias, Union @@ -45,6 +45,7 @@ class UEBCC(BaseEBCC): CDERIs = UCDERIs ElectronBoson = UElectronBoson Brueckner = BruecknerUEBCC + Optimised = OptimisedUEBCC ExternalCorrection = ExternalCorrectionUEBCC Tailor = TailorUEBCC diff --git a/ebcc/eom/base.py b/ebcc/eom/base.py index d2c4c0ef..1a739200 100644 --- a/ebcc/eom/base.py +++ b/ebcc/eom/base.py @@ -312,11 +312,11 @@ def davidson( # Check for convergence: if all(converged): self.log.debug("") - self.log.output(f"{ANSI.g}Converged.{ANSI.R}") + self.log.output(f"{ANSI.g}Converged{ANSI.R}.") else: self.log.debug("") self.log.warning( - f"{ANSI.r}Failed to converge {sum(not c for c in converged)} roots.{ANSI.R}" + f"{ANSI.r}Failed to converge {sum(not c for c in converged)} roots{ANSI.R}." ) # Update attributes: diff --git a/ebcc/ham/space.py b/ebcc/ham/space.py index 8360c6d3..11dcdd5b 100644 --- a/ebcc/ham/space.py +++ b/ebcc/ham/space.py @@ -168,6 +168,17 @@ def vmask(self, char: str) -> NDArray[B]: """ return self.mask(char)[self.virtual] + def xmask(self, char: str) -> NDArray[B]: + """Like `mask`, but returns only a mask into only the correlated sector. + + Args: + char: Character to convert. + + Returns: + Mask of the space. + """ + return self.mask(char)[self.correlated] + def oslice(self, char: str) -> _slice: """Like `slice`, but returns only a slice into only the occupied sector. @@ -194,6 +205,19 @@ def vslice(self, char: str) -> _slice: nocc = self.nocc return slice(max(s.start, nocc) - nocc, s.stop - nocc) + def xslice(self, char: str) -> _slice: + """Like `slice`, but returns only a slice into only the correlated sector. + + Args: + char: Character to convert. + + Returns: + Slice of the space. + """ + s = self.slice(char) + nfocc = self.nfocc + return slice(s.start - nfocc, s.stop - nfocc) + # Full space: @property diff --git a/ebcc/opt/__init__.py b/ebcc/opt/__init__.py index bfd9a84d..3c5703d5 100644 --- a/ebcc/opt/__init__.py +++ b/ebcc/opt/__init__.py @@ -1,5 +1,5 @@ """Orbital-optimised coupled cluster approaches.""" -from ebcc.opt.gbrueckner import BruecknerGEBCC -from ebcc.opt.rbrueckner import BruecknerREBCC -from ebcc.opt.ubrueckner import BruecknerUEBCC +from ebcc.opt.gopt import BruecknerGEBCC +from ebcc.opt.ropt import BruecknerREBCC +from ebcc.opt.uopt import BruecknerUEBCC diff --git a/ebcc/opt/base.py b/ebcc/opt/base.py index 1d699b78..3c425e85 100644 --- a/ebcc/opt/base.py +++ b/ebcc/opt/base.py @@ -6,7 +6,7 @@ from dataclasses import dataclass from typing import TYPE_CHECKING -import numpy +import numpy # FIXME for pyscfad from pyscf import lib from ebcc import numpy as np @@ -22,18 +22,19 @@ from numpy import floating from numpy.typing import NDArray - from ebcc.cc.base import BaseEBCC, SpinArrayType + from ebcc.cc.base import BaseEBCC, ERIsInputType, SpinArrayType from ebcc.core.damping import BaseDamping from ebcc.util import Namespace T = floating # FIXME Custom versions of PySCF functions +# FIXME Care with precision for PySCF @dataclass class BaseOptions(_BaseOptions): - """Options for Brueckner-orbital calculations. + """Options for Brueckner-orbital and orbital-optimised calculations. Args: e_tol: Threshold for converged in the correlation energy. @@ -45,15 +46,15 @@ class BaseOptions(_BaseOptions): """ e_tol: float = 1e-8 - t_tol: float = 1e-8 + t_tol: float = 1e-7 max_iter: int = 20 diis_space: int = 9 diis_min_space: int = 1 damping: float = 0.0 -class BaseBruecknerEBCC(ABC): - """Base class for Brueckner-orbital coupled cluster.""" +class _BaseOptimisedEBCC(ABC): + """Base class for orbital-optimised coupled cluster approaches.""" # Types Options: type[BaseOptions] = BaseOptions @@ -68,7 +69,7 @@ def __init__( options: Optional[BaseOptions] = None, **kwargs: Any, ) -> None: - r"""Initialise the Brueckner EBCC object. + r"""Initialise the orbital-optimised EBCC object. Args: cc: Parent `EBCC` object. @@ -105,18 +106,130 @@ def __init__( cc.log.info(f" > damping: {ANSI.y}{self.options.damping}{ANSI.R}") cc.log.debug("") + @property + @abstractmethod + def name(self) -> str: + """Get the name of the method.""" + pass + @property def spin_type(self) -> str: """Get the spin type.""" return self.cc.spin_type + @abstractmethod + def kernel(self) -> float: + """Run the orbital-optimised coupled cluster calculation. + + Returns: + Correlation energy. + """ + pass + + @abstractmethod + def transform_amplitudes( + self, + u: SpinArrayType, + amplitudes: Optional[Namespace[SpinArrayType]] = None, + is_lambda: bool = False, + ) -> Namespace[SpinArrayType]: + """Transform the amplitudes into the orbital-optimised basis. + + Args: + u: Rotation matrix. + amplitudes: Cluster amplitudes. + is_lambda: Whether the amplitudes are Lambda amplitudes. + + Returns: + Transformed cluster amplitudes. + """ + pass + + @abstractmethod + def mo_to_correlated(self, mo_coeff: Any) -> Any: + """Transform the MO coefficients into the correlated basis. + + Args: + mo_coeff: MO coefficients. + + Returns: + Correlated slice of MO coefficients. + """ + pass + + @abstractmethod + def mo_update_correlated(self, mo_coeff: Any, mo_coeff_corr: Any) -> Any: + """Update the correlated slice of a set of MO coefficients. + + Args: + mo_coeff: MO coefficients. + mo_coeff_corr: Correlated slice of MO coefficients. + + Returns: + Updated MO coefficients. + """ + pass + + @abstractmethod + def update_coefficients( + self, u_tot: SpinArrayType, mo_coeff_new: Any, mo_coeff_ref: Any + ) -> Any: + """Update the MO coefficients. + + Args: + u_tot: Total rotation matrix. + mo_coeff_new: New MO coefficients. + mo_coeff_ref: Reference MO coefficients. + + Returns: + Updated MO coefficients. + """ + pass + + +class BaseBruecknerEBCC(_BaseOptimisedEBCC): + """Base class for Brueckner-orbital coupled cluster.""" + @property def name(self) -> str: """Get the name of the method.""" return f"{self.spin_type}B{self.cc.ansatz.name}" + @abstractmethod + def get_t1_norm(self, amplitudes: Optional[Namespace[SpinArrayType]] = None) -> T: + """Get the norm of the T1 amplitude. + + Args: + amplitudes: Cluster amplitudes. + + Returns: + Norm of the T1 amplitude. + """ + pass + + @abstractmethod + def get_rotation_matrix( + self, + u_tot: Optional[SpinArrayType] = None, + damping: Optional[BaseDamping] = None, + t1: Optional[SpinArrayType] = None, + ) -> tuple[SpinArrayType, SpinArrayType]: + """Update the rotation matrix. + + Also returns the total rotation matrix. + + Args: + u_tot: Total rotation matrix. + damping: Damping object. + t1: T1 amplitudes. + + Returns: + Rotation matrix and total rotation matrix. + """ + pass + def kernel(self) -> float: - """Run the Bruckner-orbital coupled cluster calculation. + """Run the Brueckner orbital coupled cluster calculation. Returns: Correlation energy. @@ -124,9 +237,10 @@ def kernel(self) -> float: timer = util.Timer() # Make sure the initial CC calculation is converged: + eris = self.cc.get_eris() if not self.cc.converged: with lib.temporary_env(self.cc, log=NullLogger()): - self.cc.kernel() + self.cc.kernel(eris=eris) # Set up DIIS: damping = self.Damping(options=self.options) @@ -137,7 +251,7 @@ def kernel(self) -> float: mo_coeff_ref = self.mo_to_correlated(mo_coeff_ref) u_tot = None - self.cc.log.output("Solving for Brueckner orbitals.") + self.cc.log.output(f"Solving for {ANSI.m}Brueckner orbitals{ANSI.R}.") self.cc.log.debug("") self.log.info( f"{ANSI.B}{'Iter':>4s} {'Energy (corr.)':>16s} {'Energy (tot.)':>18s} " @@ -154,7 +268,7 @@ def kernel(self) -> float: converged = False for niter in range(1, self.options.max_iter + 1): # Update rotation matrix: - u, u_tot = self.get_rotation_matrix(u_tot=u_tot, damping=damping) + u, u_tot = self.get_rotation_matrix(u_tot=u_tot, damping=damping, t1=self.cc.t1) # Update MO coefficients: mo_coeff_new = self.update_coefficients(u_tot, mo_coeff_new, mo_coeff_ref) @@ -179,7 +293,8 @@ def kernel(self) -> float: options=self.cc.options, ) self.cc.amplitudes = amplitudes - self.cc.kernel() + eris = self.cc.get_eris() + self.cc.kernel(eris=eris) de = abs(e_prev - self.cc.e_tot) dt = self.get_t1_norm() @@ -202,11 +317,11 @@ def kernel(self) -> float: converged = converged_e and converged_t if converged: self.log.debug("") - self.log.output(f"{ANSI.g}Converged.{ANSI.R}") + self.log.output(f"{ANSI.g}Converged{ANSI.R}.") break else: self.log.debug("") - self.log.warning(f"{ANSI.r}Failed to converge.{ANSI.R}") + self.log.warning(f"{ANSI.r}Failed to converge{ANSI.R}.") self.cc.log.debug("") self.cc.log.output("E(corr) = %.10f", self.cc.e_corr) @@ -217,91 +332,189 @@ def kernel(self) -> float: return self.cc.e_corr - @abstractmethod - def get_rotation_matrix( - self, - u_tot: Optional[SpinArrayType] = None, - damping: Optional[BaseDamping] = None, - t1: Optional[SpinArrayType] = None, - ) -> tuple[SpinArrayType, SpinArrayType]: - """Update the rotation matrix. - - Also returns the total rotation matrix. - Args: - u_tot: Total rotation matrix. - damping: Damping object. - t1: T1 amplitude. +class BaseOptimisedEBCC(_BaseOptimisedEBCC): + """Base class for optimised coupled cluster.""" - Returns: - Rotation matrix and total rotation matrix. - """ - pass + @property + def name(self) -> str: + """Get the name of the method.""" + return f"{self.spin_type}O{self.cc.ansatz.name}" @abstractmethod - def transform_amplitudes( - self, u: SpinArrayType, amplitudes: Optional[Namespace[SpinArrayType]] = None - ) -> Namespace[SpinArrayType]: - """Transform the amplitudes into the Brueckner orbital basis. + def get_grad_norm(self, u: SpinArrayType) -> T: + """Get the norm of the gradient. Args: u: Rotation matrix. - amplitudes: Cluster amplitudes. Returns: - Transformed cluster amplitudes. + Norm of the gradient. """ pass @abstractmethod - def get_t1_norm(self, amplitudes: Optional[Namespace[SpinArrayType]] = None) -> T: - """Get the norm of the T1 amplitude. + def energy( + self, + eris: Optional[ERIsInputType] = None, + rdm1: Optional[SpinArrayType] = None, + rdm2: Optional[SpinArrayType] = None, + ) -> float: + """Calculate the energy. Args: - amplitudes: Cluster amplitudes. + eris: Electron repulsion integrals. + rdm1: One-particle reduced density matrix. + rdm2: Two-particle reduced density matrix. Returns: - Norm of the T1 amplitude. + Total energy. """ pass @abstractmethod - def mo_to_correlated(self, mo_coeff: Any) -> Any: - """Transform the MO coefficients into the correlated basis. + def get_rotation_matrix( + self, + u_tot: Optional[SpinArrayType] = None, + damping: Optional[BaseDamping] = None, + eris: Optional[ERIsInputType] = None, + rdm1: Optional[SpinArrayType] = None, + rdm2: Optional[SpinArrayType] = None, + ) -> tuple[SpinArrayType, SpinArrayType]: + """Update the rotation matrix. + + Also returns the total rotation matrix. Args: - mo_coeff: MO coefficients. + u_tot: Total rotation matrix. + damping: Damping object. + eris: Electronic repulsion integrals. + rdm1: One-particle reduced density matrix. + rdm2: Two-particle reduced density matrix. Returns: - Correlated slice of MO coefficients. + Rotation matrix and total rotation matrix. """ pass - @abstractmethod - def mo_update_correlated(self, mo_coeff: Any, mo_coeff_corr: Any) -> Any: - """Update the correlated slice of a set of MO coefficients. - - Args: - mo_coeff: MO coefficients. - mo_coeff_corr: Correlated slice of MO coefficients. + def kernel(self) -> float: + """Run the optimised coupled cluster calculation. Returns: - Updated MO coefficients. + Correlation energy. """ - pass + timer = util.Timer() - @abstractmethod - def update_coefficients( - self, u_tot: SpinArrayType, mo_coeff_new: Any, mo_coeff_ref: Any - ) -> Any: - """Update the MO coefficients. + # Make sure the initial CC calculation is converged: + eris = self.cc.get_eris() + with lib.temporary_env(self.cc, log=NullLogger()): + if not self.cc.converged: + self.cc.kernel(eris=eris) + if not self.cc.converged_lambda: + self.cc.solve_lambda(eris=eris) - Args: - u_tot: Total rotation matrix. - mo_coeff_new: New MO coefficients. - mo_coeff_ref: Reference MO coefficients. + # Set up DIIS: + damping = self.Damping(options=self.options) - Returns: - Updated MO coefficients. - """ - pass + # Initialise coefficients: + mo_coeff_new: NDArray[T] = np.copy(np.asarray(self.cc.mo_coeff, dtype=types[float])) + mo_coeff_ref: NDArray[T] = np.copy(np.asarray(self.cc.mo_coeff, dtype=types[float])) + mo_coeff_ref = self.mo_to_correlated(mo_coeff_ref) + u_tot = None + + # Initialise density matrices: + rdm1 = self.cc.make_rdm1_f(eris=eris) + rdm2 = self.cc.make_rdm2_f(eris=eris) + assert np.allclose( + self.cc.e_tot, self.energy(rdm1=rdm1, rdm2=rdm2, eris=eris) + ) # FIXME remove + + self.cc.log.output(f"Solving for {ANSI.m}optimised orbitals{ANSI.R}.") + self.cc.log.debug("") + self.log.info( + f"{ANSI.B}{'Iter':>4s} {'Energy (corr.)':>16s} {'Energy (tot.)':>18s} " + f"{'Δ(Energy)':>13s} {'|Gradient|':>13s}{ANSI.R}" + ) + self.log.info( + "%4d %16.10f %18.10f", + 0, + self.cc.e_corr, + self.cc.e_tot, + ) + + converged = False + for niter in range(1, self.options.max_iter + 1): + # Update rotation matrix: + u, u_tot = self.get_rotation_matrix( + u_tot=u_tot, damping=damping, eris=eris, rdm1=rdm1, rdm2=rdm2 + ) + + # Update MO coefficients: + mo_coeff_new = self.update_coefficients(u_tot, mo_coeff_new, mo_coeff_ref) + + # Transform mean-field and amplitudes: + self.mf.mo_coeff = numpy.asarray(mo_coeff_new) + self.mf.e_tot = self.mf.energy_tot() + amplitudes = self.transform_amplitudes(u) + lambas = self.transform_amplitudes(u, is_lambda=True) + + # Update CC calculation: + e_prev = self.cc.e_tot + with lib.temporary_env(self.cc, log=NullLogger()): + self.cc.__class__.__init__( + self.cc, + self.mf, + log=self.cc.log, + ansatz=self.cc.ansatz, + space=self.cc.space, + omega=self.cc.omega, + g=self.cc.bare_g, + G=self.cc.bare_G, + options=self.cc.options, + ) + self.cc.amplitudes = amplitudes + self.cc.lambdas = lambas + eris = self.cc.get_eris() + + # Update the density matrices: + rdm1 = self.cc.make_rdm1_f(eris=eris) + rdm2 = self.cc.make_rdm2_f(eris=eris) + + # Update the energy: + e_tot = self.energy(rdm1=rdm1, rdm2=rdm2, eris=eris) + self.cc.e_corr = e_tot - self.cc.mf.e_tot + + # Log the iteration: + de = abs(e_prev - self.cc.e_tot) + grad = self.get_grad_norm(u) + converged_e = bool(de < self.options.e_tol) + converged_grad = bool(grad < self.options.t_tol) + self.log.info( + f"%4s %16.10f %18.10f" + f" {[ANSI.r, ANSI.g][int(converged_e)]}%13.3e{ANSI.R}" + f" {[ANSI.r, ANSI.g][int(converged_grad)]}%13.3e{ANSI.R}", + niter, + self.cc.e_corr, + self.cc.e_tot, + de, + grad, + ) + + # Check for convergence: + converged = converged_e and converged_grad + if converged: + self.log.debug("") + self.log.output(f"{ANSI.g}Converged{ANSI.R}.") + break + else: + self.log.debug("") + self.log.warning(f"{ANSI.r}Failed to converge{ANSI.R}.") + + self.cc.log.debug("") + self.cc.log.output("E(corr) = %.10f", self.cc.e_corr) + self.cc.log.output("E(tot) = %.10f", self.cc.e_tot) + self.cc.log.debug("") + self.cc.log.debug("Time elapsed: %s", timer.format_time(timer())) + self.cc.log.debug("") + + return self.cc.e_corr diff --git a/ebcc/opt/gbrueckner.py b/ebcc/opt/gbrueckner.py deleted file mode 100644 index 18c5d66d..00000000 --- a/ebcc/opt/gbrueckner.py +++ /dev/null @@ -1,174 +0,0 @@ -"""Generalised Brueckner-orbital coupled cluster.""" - -from __future__ import annotations - -from typing import TYPE_CHECKING - -import scipy.linalg - -from ebcc import numpy as np -from ebcc import util -from ebcc.backend import _put -from ebcc.core.precision import types -from ebcc.opt.base import BaseBruecknerEBCC - -if TYPE_CHECKING: - from typing import Optional, Union - - from numpy import floating - from numpy.typing import NDArray - - from ebcc.cc.gebcc import GEBCC, SpinArrayType - from ebcc.core.damping import BaseDamping - from ebcc.util import Namespace - - T = floating - - -class BruecknerGEBCC(BaseBruecknerEBCC): - """Generalised Brueckner-orbital coupled cluster.""" - - # Attributes - cc: GEBCC - - def get_rotation_matrix( - self, - u_tot: Optional[SpinArrayType] = None, - damping: Optional[BaseDamping] = None, - t1: Optional[SpinArrayType] = None, - ) -> tuple[SpinArrayType, SpinArrayType]: - """Update the rotation matrix. - - Also returns the total rotation matrix. - - Args: - u_tot: Total rotation matrix. - damping: Damping object. - t1: T1 amplitude. - - Returns: - Rotation matrix and total rotation matrix. - """ - if t1 is None: - t1 = self.cc.t1 - if u_tot is None: - u_tot = np.eye(self.cc.space.ncorr, dtype=types[float]) - - zocc = np.zeros((self.cc.space.ncocc, self.cc.space.ncocc)) - zvir = np.zeros((self.cc.space.ncvir, self.cc.space.ncvir)) - t1_block: NDArray[T] = np.block([[zocc, -t1], [np.transpose(t1), zvir]]) - - u = scipy.linalg.expm(t1_block) - - u_tot = u_tot @ u - if np.linalg.det(u_tot) < 0: - u_tot = _put(u_tot, np.ix_(np.arange(u_tot.shape[0]), np.array([0])), -u_tot[:, 0]) - - a: NDArray[T] = np.asarray(np.real(scipy.linalg.logm(u_tot)), dtype=types[float]) - if damping is not None: - a = damping(a, error=t1) - - u_tot = scipy.linalg.expm(a) - - return u, u_tot - - def transform_amplitudes( - self, u: SpinArrayType, amplitudes: Optional[Namespace[SpinArrayType]] = None - ) -> Namespace[SpinArrayType]: - """Transform the amplitudes into the Brueckner orbital basis. - - Args: - u: Rotation matrix. - amplitudes: Cluster amplitudes. - - Returns: - Transformed cluster amplitudes. - """ - if not amplitudes: - amplitudes = self.cc.amplitudes - - nocc = self.cc.space.ncocc - ci = u[:nocc, :nocc] - ca = u[nocc:, nocc:] - - # Transform T amplitudes: - for name, key, n in self.cc.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type): - args: list[Union[tuple[int, ...], NDArray[T]]] = [ - self.cc.amplitudes[name], - tuple(range(n * 2)), - ] - for i in range(n): - args += [ci, (i, i + n * 2)] - for i in range(n): - args += [ca, (i + n, i + n * 3)] - args += [tuple(range(n * 2, n * 4))] - self.cc.amplitudes[name] = util.einsum(*args) - - # Transform S amplitudes: - for name, key, n in self.cc.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type): - raise util.ModelNotImplemented # TODO - - # Transform U amplitudes: - for name, key, nf, nb in self.cc.ansatz.coupling_cluster_ranks(spin_type=self.spin_type): - raise util.ModelNotImplemented # TODO - - return self.cc.amplitudes - - def get_t1_norm(self, amplitudes: Optional[Namespace[SpinArrayType]] = None) -> T: - """Get the norm of the T1 amplitude. - - Args: - amplitudes: Cluster amplitudes. - - Returns: - Norm of the T1 amplitude. - """ - if not amplitudes: - amplitudes = self.cc.amplitudes - weight: T = np.linalg.norm(amplitudes["t1"]) - return weight - - def mo_to_correlated(self, mo_coeff: NDArray[T]) -> NDArray[T]: - """Transform the MO coefficients into the correlated basis. - - Args: - mo_coeff: MO coefficients. - - Returns: - Correlated slice of MO coefficients. - """ - return mo_coeff[:, self.cc.space.correlated] - - def mo_update_correlated(self, mo_coeff: NDArray[T], mo_coeff_corr: NDArray[T]) -> NDArray[T]: - """Update the correlated slice of a set of MO coefficients. - - Args: - mo_coeff: MO coefficients. - mo_coeff_corr: Correlated slice of MO coefficients. - - Returns: - Updated MO coefficients. - """ - mo_coeff = _put( - mo_coeff, - np.ix_(np.arange(mo_coeff.shape[0]), self.cc.space.correlated), # type: ignore - mo_coeff_corr, - ) - return mo_coeff - - def update_coefficients( - self, u_tot: SpinArrayType, mo_coeff: NDArray[T], mo_coeff_ref: NDArray[T] - ) -> NDArray[T]: - """Update the MO coefficients. - - Args: - u_tot: Total rotation matrix. - mo_coeff: New MO coefficients. - mo_coeff_ref: Reference MO coefficients. - - Returns: - Updated MO coefficients. - """ - mo_coeff_new_corr = util.einsum("pi,ij->pj", mo_coeff_ref, u_tot) - mo_coeff_new = self.mo_update_correlated(mo_coeff, mo_coeff_new_corr) - return mo_coeff_new diff --git a/ebcc/opt/gopt.py b/ebcc/opt/gopt.py new file mode 100644 index 00000000..24ece314 --- /dev/null +++ b/ebcc/opt/gopt.py @@ -0,0 +1,305 @@ +"""Generalised orbital-optimised coupled cluster.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import scipy.linalg + +from ebcc import numpy as np +from ebcc import util +from ebcc.backend import _inflate, _put, ensure_scalar +from ebcc.core.precision import astype, types +from ebcc.opt.base import BaseBruecknerEBCC, BaseOptimisedEBCC, _BaseOptimisedEBCC + +if TYPE_CHECKING: + from typing import Optional, Union + + from numpy import floating + from numpy.typing import NDArray + + from ebcc.cc.gebcc import GEBCC, ERIsInputType, SpinArrayType + from ebcc.core.damping import BaseDamping + from ebcc.util import Namespace + + T = floating + + +class _OptimisedGEBCC(_BaseOptimisedEBCC): + """Generalised orbital-optimised coupled cluster.""" + + # Attributes + cc: GEBCC + + def transform_amplitudes( + self, + u: SpinArrayType, + amplitudes: Optional[Namespace[SpinArrayType]] = None, + is_lambda: bool = False, + ) -> Namespace[SpinArrayType]: + """Transform the amplitudes into the orbital-optimised basis. + + Args: + u: Rotation matrix. + amplitudes: Cluster amplitudes. + is_lambda: Whether the amplitudes are Lambda amplitudes. + + Returns: + Transformed cluster amplitudes. + """ + if not amplitudes and not is_lambda: + amplitudes = self.cc.amplitudes + elif not amplitudes: + amplitudes = self.cc.lambdas + amplitudes = amplitudes.copy() + + nocc = self.cc.space.ncocc + ci = u[:nocc, :nocc] + ca = u[nocc:, nocc:] + if is_lambda: + ci, ca = ca, ci + + # Transform T amplitudes: + for name, key, n in self.cc.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type): + if is_lambda: + name = name.replace("t", "l") + args: list[Union[tuple[int, ...], NDArray[T]]] = [ + amplitudes[name], + tuple(range(n * 2)), + ] + for i in range(n): + args += [ci, (i, i + n * 2)] + for i in range(n): + args += [ca, (i + n, i + n * 3)] + args += [tuple(range(n * 2, n * 4))] + amplitudes[name] = util.einsum(*args) + + # Transform S amplitudes: + for name, key, n in self.cc.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type): + raise util.ModelNotImplemented # TODO + + # Transform U amplitudes: + for name, key, nf, nb in self.cc.ansatz.coupling_cluster_ranks(spin_type=self.spin_type): + raise util.ModelNotImplemented # TODO + + return amplitudes + + def mo_to_correlated(self, mo_coeff: NDArray[T]) -> NDArray[T]: + """Transform the MO coefficients into the correlated basis. + + Args: + mo_coeff: MO coefficients. + + Returns: + Correlated slice of MO coefficients. + """ + return mo_coeff[:, self.cc.space.correlated] + + def mo_update_correlated(self, mo_coeff: NDArray[T], mo_coeff_corr: NDArray[T]) -> NDArray[T]: + """Update the correlated slice of a set of MO coefficients. + + Args: + mo_coeff: MO coefficients. + mo_coeff_corr: Correlated slice of MO coefficients. + + Returns: + Updated MO coefficients. + """ + mo_coeff = _put( + mo_coeff, + np.ix_(np.arange(mo_coeff.shape[0]), self.cc.space.correlated), # type: ignore + mo_coeff_corr, + ) + return mo_coeff + + def update_coefficients( + self, u_tot: SpinArrayType, mo_coeff: NDArray[T], mo_coeff_ref: NDArray[T] + ) -> NDArray[T]: + """Update the MO coefficients. + + Args: + u_tot: Total rotation matrix. + mo_coeff: New MO coefficients. + mo_coeff_ref: Reference MO coefficients. + + Returns: + Updated MO coefficients. + """ + mo_coeff_new_corr = util.einsum("pi,ij->pj", mo_coeff_ref, u_tot) + mo_coeff_new = self.mo_update_correlated(mo_coeff, mo_coeff_new_corr) + return mo_coeff_new + + +class BruecknerGEBCC(BaseBruecknerEBCC, _OptimisedGEBCC): + """Generalised Brueckner-orbital coupled cluster.""" + + def get_t1_norm(self, amplitudes: Optional[Namespace[SpinArrayType]] = None) -> T: + """Get the norm of the T1 amplitude. + + Args: + amplitudes: Cluster amplitudes. + + Returns: + Norm of the T1 amplitude. + """ + if not amplitudes: + amplitudes = self.cc.amplitudes + weight: T = np.linalg.norm(amplitudes["t1"]) + return weight + + def get_rotation_matrix( + self, + u_tot: Optional[SpinArrayType] = None, + damping: Optional[BaseDamping] = None, + t1: Optional[SpinArrayType] = None, + ) -> tuple[SpinArrayType, SpinArrayType]: + """Update the rotation matrix. + + Also returns the total rotation matrix. + + Args: + u_tot: Total rotation matrix. + damping: Damping object. + t1: T1 amplitudes. + + Returns: + Rotation matrix and total rotation matrix. + """ + if u_tot is None: + u_tot = np.eye(self.cc.space.ncorr, dtype=types[float]) + if t1 is None: + t1 = self.cc.t1 + + # Get the T1 blocks + zocc = np.zeros((self.cc.space.ncocc, self.cc.space.ncocc)) + zvir = np.zeros((self.cc.space.ncvir, self.cc.space.ncvir)) + t1_block: NDArray[T] = np.block([[zocc, -t1], [np.transpose(t1), zvir]]) + + # Get the orbital gradient + u = scipy.linalg.expm(t1_block) + + u_tot = u_tot @ u + if np.linalg.det(u_tot) < 0: + u_tot = _put(u_tot, np.ix_(np.arange(u_tot.shape[0]), np.array([0])), -u_tot[:, 0]) + + a: NDArray[T] = np.asarray(np.real(scipy.linalg.logm(u_tot)), dtype=types[float]) + if damping is not None: + a = damping(a, error=t1) + + u_tot = scipy.linalg.expm(a) + + return u, u_tot + + +class OptimisedGEBCC(BaseOptimisedEBCC, _OptimisedGEBCC): + """Generalised optimised coupled cluster.""" + + def get_grad_norm(self, u: SpinArrayType) -> T: + """Get the norm of the gradient. + + Args: + u: Rotation matrix. + + Returns: + Norm of the gradient. + """ + grad: T = np.linalg.norm(scipy.linalg.logm(u)) + return grad + + def get_rotation_matrix( + self, + u_tot: Optional[SpinArrayType] = None, + damping: Optional[BaseDamping] = None, + eris: Optional[ERIsInputType] = None, + rdm1: Optional[SpinArrayType] = None, + rdm2: Optional[SpinArrayType] = None, + ) -> tuple[SpinArrayType, SpinArrayType]: + """Update the rotation matrix. + + Also returns the total rotation matrix. + + Args: + u_tot: Total rotation matrix. + damping: Damping object. + eris: Electronic repulsion integrals. + rdm1: One-particle reduced density matrix. + rdm2: Two-particle reduced density matrix. + + Returns: + Rotation matrix and total rotation matrix. + """ + if u_tot is None: + u_tot = np.eye(self.cc.space.ncorr, dtype=types[float]) + if rdm1 is None: + rdm1 = self.cc.make_rdm1_f(eris=eris) + if rdm2 is None: + rdm2 = self.cc.make_rdm2_f(eris=eris) + eris = self.cc.get_eris(eris=eris) + + # Get one-particle contribution to generalised Fock matrix + mo_coeff = self.cc.mo_coeff[:, self.cc.space.correlated] + h1e_ao: NDArray[T] = np.asarray(self.cc.mf.get_hcore(), dtype=types[float]) + h1e = util.einsum("pq,pi,qj->ij", h1e_ao, mo_coeff, mo_coeff) + f = np.dot(h1e, rdm1) * 0.5 + + # Get two-particle contribution to generalised Fock matrix + h2e = eris.xxxx + f += util.einsum("prst,stqr->pq", h2e, rdm2) * 0.5 + + # Get the orbital gradient + o, v = self.cc.space.xslice("o"), self.cc.space.xslice("v") + x_vo = (f - np.transpose(f))[v, o] / self.cc.energy_sum("vo") + x = _inflate( + (self.cc.space.ncorr, self.cc.space.ncorr), + np.ix_(self.cc.space.xmask("v"), self.cc.space.xmask("o")), # type: ignore + x_vo, + ) + u = scipy.linalg.expm(x - np.transpose(x)) + + u_tot = u_tot @ u + if np.linalg.det(u_tot) < 0: + u_tot = _put(u_tot, np.ix_(np.arange(u_tot.shape[0]), np.array([0])), -u_tot[:, 0]) + + a: NDArray[T] = np.asarray(np.real(scipy.linalg.logm(u_tot)), dtype=types[float]) + if damping is not None: + a = damping(a) + + u_tot = scipy.linalg.expm(a) + + return u, u_tot + + def energy( + self, + eris: Optional[ERIsInputType] = None, + rdm1: Optional[SpinArrayType] = None, + rdm2: Optional[SpinArrayType] = None, + ) -> float: + """Calculate the energy. + + Args: + eris: Electron repulsion integrals. + rdm1: One-particle reduced density matrix. + rdm2: Two-particle reduced density matrix. + + Returns: + Total energy. + """ + if rdm1 is None: + rdm1 = self.cc.make_rdm1_f(eris=eris) + if rdm2 is None: + rdm2 = self.cc.make_rdm2_f(eris=eris) + eris = self.cc.get_eris(eris=eris) + + # Get the Hamiltonian + h1e = util.einsum( + "pq,pi,qj->ij", self.cc.mf.get_hcore(), self.cc.mo_coeff, self.cc.mo_coeff + ) + h2e = eris.xxxx + + # Get the energy + e_tot = util.einsum("pq,qp->", h1e, rdm1) + e_tot += util.einsum("pqrs,pqrs->", h2e, rdm2) * 0.5 + e_tot += self.cc.mf.energy_nuc() + + res: float = np.real(ensure_scalar(e_tot)) + return astype(res, float) diff --git a/ebcc/opt/rbrueckner.py b/ebcc/opt/rbrueckner.py deleted file mode 100644 index dc16fb34..00000000 --- a/ebcc/opt/rbrueckner.py +++ /dev/null @@ -1,174 +0,0 @@ -"""Restricted Brueckner-orbital coupled cluster.""" - -from __future__ import annotations - -from typing import TYPE_CHECKING - -import scipy.linalg - -from ebcc import numpy as np -from ebcc import util -from ebcc.backend import _put -from ebcc.core.precision import types -from ebcc.opt.base import BaseBruecknerEBCC - -if TYPE_CHECKING: - from typing import Optional, Union - - from numpy import floating - from numpy.typing import NDArray - - from ebcc.cc.rebcc import REBCC, SpinArrayType - from ebcc.core.damping import BaseDamping - from ebcc.util import Namespace - - T = floating - - -class BruecknerREBCC(BaseBruecknerEBCC): - """Restricted Brueckner-orbital coupled cluster.""" - - # Attributes - cc: REBCC - - def get_rotation_matrix( - self, - u_tot: Optional[SpinArrayType] = None, - damping: Optional[BaseDamping] = None, - t1: Optional[SpinArrayType] = None, - ) -> tuple[SpinArrayType, SpinArrayType]: - """Update the rotation matrix. - - Also returns the total rotation matrix. - - Args: - u_tot: Total rotation matrix. - damping: Damping object. - t1: T1 amplitude. - - Returns: - Rotation matrix and total rotation matrix. - """ - if t1 is None: - t1 = self.cc.t1 - if u_tot is None: - u_tot = np.eye(self.cc.space.ncorr, dtype=types[float]) - - zocc = np.zeros((self.cc.space.ncocc, self.cc.space.ncocc)) - zvir = np.zeros((self.cc.space.ncvir, self.cc.space.ncvir)) - t1_block: NDArray[T] = np.block([[zocc, -t1], [np.transpose(t1), zvir]]) - - u = scipy.linalg.expm(t1_block) - - u_tot = u_tot @ u - if np.linalg.det(u_tot) < 0: - u_tot = _put(u_tot, np.ix_(np.arange(u_tot.shape[0]), np.array([0])), -u_tot[:, 0]) - - a: NDArray[T] = np.asarray(np.real(scipy.linalg.logm(u_tot)), dtype=types[float]) - if damping is not None: - a = damping(a, error=t1) - - u_tot = scipy.linalg.expm(a) - - return u, u_tot - - def transform_amplitudes( - self, u: SpinArrayType, amplitudes: Optional[Namespace[SpinArrayType]] = None - ) -> Namespace[SpinArrayType]: - """Transform the amplitudes into the Brueckner orbital basis. - - Args: - u: Rotation matrix. - amplitudes: Cluster amplitudes. - - Returns: - Transformed cluster amplitudes. - """ - if not amplitudes: - amplitudes = self.cc.amplitudes - - nocc = self.cc.space.ncocc - ci = u[:nocc, :nocc] - ca = u[nocc:, nocc:] - - # Transform T amplitudes: - for name, key, n in self.cc.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type): - args: list[Union[SpinArrayType, tuple[int, ...]]] = [ - self.cc.amplitudes[name], - tuple(range(n * 2)), - ] - for i in range(n): - args += [ci, (i, i + n * 2)] - for i in range(n): - args += [ca, (i + n, i + n * 3)] - args += [tuple(range(n * 2, n * 4))] - self.cc.amplitudes[name] = util.einsum(*args) - - # Transform S amplitudes: - for name, key, n in self.cc.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type): - raise util.ModelNotImplemented # TODO - - # Transform U amplitudes: - for name, key, nf, nb in self.cc.ansatz.coupling_cluster_ranks(spin_type=self.spin_type): - raise util.ModelNotImplemented # TODO - - return self.cc.amplitudes - - def get_t1_norm(self, amplitudes: Optional[Namespace[SpinArrayType]] = None) -> T: - """Get the norm of the T1 amplitude. - - Args: - amplitudes: Cluster amplitudes. - - Returns: - Norm of the T1 amplitude. - """ - if not amplitudes: - amplitudes = self.cc.amplitudes - weight: T = np.linalg.norm(amplitudes["t1"]) - return weight - - def mo_to_correlated(self, mo_coeff: NDArray[T]) -> NDArray[T]: - """Transform the MO coefficients into the correlated basis. - - Args: - mo_coeff: MO coefficients. - - Returns: - Correlated slice of MO coefficients. - """ - return mo_coeff[:, self.cc.space.correlated] - - def mo_update_correlated(self, mo_coeff: NDArray[T], mo_coeff_corr: NDArray[T]) -> NDArray[T]: - """Update the correlated slice of a set of MO coefficients. - - Args: - mo_coeff: MO coefficients. - mo_coeff_corr: Correlated slice of MO coefficients. - - Returns: - Updated MO coefficients. - """ - mo_coeff = _put( - mo_coeff, - np.ix_(np.arange(mo_coeff.shape[0]), self.cc.space.correlated), # type: ignore - mo_coeff_corr, - ) - return mo_coeff - - def update_coefficients( - self, u_tot: SpinArrayType, mo_coeff: NDArray[T], mo_coeff_ref: NDArray[T] - ) -> NDArray[T]: - """Update the MO coefficients. - - Args: - u_tot: Total rotation matrix. - mo_coeff: New MO coefficients. - mo_coeff_ref: Reference MO coefficients. - - Returns: - Updated MO coefficients. - """ - mo_coeff_new_corr = util.einsum("pi,ij->pj", mo_coeff_ref, u_tot) - mo_coeff_new = self.mo_update_correlated(mo_coeff, mo_coeff_new_corr) - return mo_coeff_new diff --git a/ebcc/opt/ropt.py b/ebcc/opt/ropt.py new file mode 100644 index 00000000..255713db --- /dev/null +++ b/ebcc/opt/ropt.py @@ -0,0 +1,305 @@ +"""Restricted orbital-optimised coupled cluster.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import scipy.linalg + +from ebcc import numpy as np +from ebcc import util +from ebcc.backend import _inflate, _put, ensure_scalar +from ebcc.core.precision import astype, types +from ebcc.opt.base import BaseBruecknerEBCC, BaseOptimisedEBCC, _BaseOptimisedEBCC + +if TYPE_CHECKING: + from typing import Optional, Union + + from numpy import floating + from numpy.typing import NDArray + + from ebcc.cc.rebcc import REBCC, ERIsInputType, SpinArrayType + from ebcc.core.damping import BaseDamping + from ebcc.util import Namespace + + T = floating + + +class _OptimisedREBCC(_BaseOptimisedEBCC): + """Restricted orbital-optimised coupled cluster.""" + + # Attributes + cc: REBCC + + def transform_amplitudes( + self, + u: SpinArrayType, + amplitudes: Optional[Namespace[SpinArrayType]] = None, + is_lambda: bool = False, + ) -> Namespace[SpinArrayType]: + """Transform the amplitudes into the orbital-optimised basis. + + Args: + u: Rotation matrix. + amplitudes: Cluster amplitudes. + is_lambda: Whether the amplitudes are Lambda amplitudes. + + Returns: + Transformed cluster amplitudes. + """ + if not amplitudes and not is_lambda: + amplitudes = self.cc.amplitudes + elif not amplitudes: + amplitudes = self.cc.lambdas + amplitudes = amplitudes.copy() + + nocc = self.cc.space.ncocc + ci = u[:nocc, :nocc] + ca = u[nocc:, nocc:] + if is_lambda: + ci, ca = ca, ci + + # Transform T amplitudes: + for name, key, n in self.cc.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type): + if is_lambda: + name = name.replace("t", "l") + args: list[Union[SpinArrayType, tuple[int, ...]]] = [ + amplitudes[name], + tuple(range(n * 2)), + ] + for i in range(n): + args += [ci, (i, i + n * 2)] + for i in range(n): + args += [ca, (i + n, i + n * 3)] + args += [tuple(range(n * 2, n * 4))] + amplitudes[name] = util.einsum(*args) + + # Transform S amplitudes: + for name, key, n in self.cc.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type): + raise util.ModelNotImplemented # TODO + + # Transform U amplitudes: + for name, key, nf, nb in self.cc.ansatz.coupling_cluster_ranks(spin_type=self.spin_type): + raise util.ModelNotImplemented # TODO + + return amplitudes + + def mo_to_correlated(self, mo_coeff: NDArray[T]) -> NDArray[T]: + """Transform the MO coefficients into the correlated basis. + + Args: + mo_coeff: MO coefficients. + + Returns: + Correlated slice of MO coefficients. + """ + return mo_coeff[:, self.cc.space.correlated] + + def mo_update_correlated(self, mo_coeff: NDArray[T], mo_coeff_corr: NDArray[T]) -> NDArray[T]: + """Update the correlated slice of a set of MO coefficients. + + Args: + mo_coeff: MO coefficients. + mo_coeff_corr: Correlated slice of MO coefficients. + + Returns: + Updated MO coefficients. + """ + mo_coeff = _put( + mo_coeff, + np.ix_(np.arange(mo_coeff.shape[0]), self.cc.space.correlated), # type: ignore + mo_coeff_corr, + ) + return mo_coeff + + def update_coefficients( + self, u_tot: SpinArrayType, mo_coeff: NDArray[T], mo_coeff_ref: NDArray[T] + ) -> NDArray[T]: + """Update the MO coefficients. + + Args: + u_tot: Total rotation matrix. + mo_coeff: New MO coefficients. + mo_coeff_ref: Reference MO coefficients. + + Returns: + Updated MO coefficients. + """ + mo_coeff_new_corr = util.einsum("pi,ij->pj", mo_coeff_ref, u_tot) + mo_coeff_new = self.mo_update_correlated(mo_coeff, mo_coeff_new_corr) + return mo_coeff_new + + +class BruecknerREBCC(BaseBruecknerEBCC, _OptimisedREBCC): + """Restricted Brueckner-orbital coupled cluster.""" + + def get_t1_norm(self, amplitudes: Optional[Namespace[SpinArrayType]] = None) -> T: + """Get the norm of the T1 amplitude. + + Args: + amplitudes: Cluster amplitudes. + + Returns: + Norm of the T1 amplitude. + """ + if not amplitudes: + amplitudes = self.cc.amplitudes + weight: T = np.linalg.norm(amplitudes["t1"]) + return weight + + def get_rotation_matrix( + self, + u_tot: Optional[SpinArrayType] = None, + damping: Optional[BaseDamping] = None, + t1: Optional[SpinArrayType] = None, + ) -> tuple[SpinArrayType, SpinArrayType]: + """Update the rotation matrix. + + Also returns the total rotation matrix. + + Args: + u_tot: Total rotation matrix. + damping: Damping object. + t1: T1 amplitudes. + + Returns: + Rotation matrix and total rotation matrix. + """ + if u_tot is None: + u_tot = np.eye(self.cc.space.ncorr, dtype=types[float]) + if t1 is None: + t1 = self.cc.t1 + + # Get the T1 blocks + zocc = np.zeros((self.cc.space.ncocc, self.cc.space.ncocc)) + zvir = np.zeros((self.cc.space.ncvir, self.cc.space.ncvir)) + t1_block: NDArray[T] = np.block([[zocc, -t1], [np.transpose(t1), zvir]]) + + # Get the orbital gradient + u = scipy.linalg.expm(t1_block) + + u_tot = u_tot @ u + if np.linalg.det(u_tot) < 0: + u_tot = _put(u_tot, np.ix_(np.arange(u_tot.shape[0]), np.array([0])), -u_tot[:, 0]) + + a: NDArray[T] = np.asarray(np.real(scipy.linalg.logm(u_tot)), dtype=types[float]) + if damping is not None: + a = damping(a, error=t1) + + u_tot = scipy.linalg.expm(a) + + return u, u_tot + + +class OptimisedREBCC(BaseOptimisedEBCC, _OptimisedREBCC): + """Restricted optimised coupled cluster.""" + + def get_grad_norm(self, u: SpinArrayType) -> T: + """Get the norm of the gradient. + + Args: + u: Rotation matrix. + + Returns: + Norm of the gradient. + """ + grad: T = np.linalg.norm(scipy.linalg.logm(u)) + return grad + + def get_rotation_matrix( + self, + u_tot: Optional[SpinArrayType] = None, + damping: Optional[BaseDamping] = None, + eris: Optional[ERIsInputType] = None, + rdm1: Optional[SpinArrayType] = None, + rdm2: Optional[SpinArrayType] = None, + ) -> tuple[SpinArrayType, SpinArrayType]: + """Update the rotation matrix. + + Also returns the total rotation matrix. + + Args: + u_tot: Total rotation matrix. + damping: Damping object. + eris: Electronic repulsion integrals. + rdm1: One-particle reduced density matrix. + rdm2: Two-particle reduced density matrix. + + Returns: + Rotation matrix and total rotation matrix. + """ + if u_tot is None: + u_tot = np.eye(self.cc.space.ncorr, dtype=types[float]) + if rdm1 is None: + rdm1 = self.cc.make_rdm1_f(eris=eris) + if rdm2 is None: + rdm2 = self.cc.make_rdm2_f(eris=eris) + eris = self.cc.get_eris(eris=eris) + + # Get one-particle contribution to generalised Fock matrix + mo_coeff = self.cc.mo_coeff[:, self.cc.space.correlated] + h1e_ao: NDArray[T] = np.asarray(self.cc.mf.get_hcore(), dtype=types[float]) + h1e = util.einsum("pq,pi,qj->ij", h1e_ao, mo_coeff, mo_coeff) + f = np.dot(h1e, rdm1) * 0.5 + + # Get two-particle contribution to generalised Fock matrix + h2e = eris.xxxx + f += util.einsum("prst,stqr->pq", h2e, rdm2) * 0.5 + + # Get the orbital gradient + o, v = self.cc.space.xslice("o"), self.cc.space.xslice("v") + x_vo = 0.5 * (f - np.transpose(f))[v, o] / self.cc.energy_sum("vo") + x = _inflate( + (self.cc.space.ncorr, self.cc.space.ncorr), + np.ix_(self.cc.space.xmask("v"), self.cc.space.xmask("o")), # type: ignore + x_vo, + ) + u = scipy.linalg.expm(x - np.transpose(x)) + + u_tot = u_tot @ u + if np.linalg.det(u_tot) < 0: + u_tot = _put(u_tot, np.ix_(np.arange(u_tot.shape[0]), np.array([0])), -u_tot[:, 0]) + + a: NDArray[T] = np.asarray(np.real(scipy.linalg.logm(u_tot)), dtype=types[float]) + if damping is not None: + a = damping(a) # FIXME can x be the error? + + u_tot = scipy.linalg.expm(a) + + return u, u_tot + + def energy( + self, + eris: Optional[ERIsInputType] = None, + rdm1: Optional[SpinArrayType] = None, + rdm2: Optional[SpinArrayType] = None, + ) -> float: + """Calculate the energy. + + Args: + eris: Electron repulsion integrals. + rdm1: One-particle reduced density matrix. + rdm2: Two-particle reduced density matrix. + + Returns: + Total energy. + """ + if rdm1 is None: + rdm1 = self.cc.make_rdm1_f(eris=eris) + if rdm2 is None: + rdm2 = self.cc.make_rdm2_f(eris=eris) + eris = self.cc.get_eris(eris=eris) + + # Get the Hamiltonian + h1e = util.einsum( + "pq,pi,qj->ij", self.cc.mf.get_hcore(), self.cc.mo_coeff, self.cc.mo_coeff + ) + h2e = eris.xxxx + + # Get the energy + e_tot = util.einsum("pq,qp->", h1e, rdm1) + e_tot += util.einsum("pqrs,pqrs->", h2e, rdm2) * 0.5 + e_tot += self.cc.mf.energy_nuc() + + res: float = np.real(ensure_scalar(e_tot)) + return astype(res, float) diff --git a/ebcc/opt/ubrueckner.py b/ebcc/opt/ubrueckner.py deleted file mode 100644 index b4db762d..00000000 --- a/ebcc/opt/ubrueckner.py +++ /dev/null @@ -1,216 +0,0 @@ -"""Unrestricted Brueckner-orbital coupled cluster.""" - -from __future__ import annotations - -from typing import TYPE_CHECKING - -import scipy.linalg - -from ebcc import numpy as np -from ebcc import util -from ebcc.backend import _put -from ebcc.core.precision import types -from ebcc.opt.base import BaseBruecknerEBCC - -if TYPE_CHECKING: - from typing import Optional - - from numpy import floating - from numpy.typing import NDArray - - from ebcc.cc.uebcc import UEBCC, SpinArrayType - from ebcc.core.damping import BaseDamping - from ebcc.util import Namespace - - T = floating - - -class BruecknerUEBCC(BaseBruecknerEBCC): - """Unrestricted Brueckner-orbital coupled cluster.""" - - # Attributes - cc: UEBCC - - def get_rotation_matrix( - self, - u_tot: Optional[SpinArrayType] = None, - damping: Optional[BaseDamping] = None, - t1: Optional[SpinArrayType] = None, - ) -> tuple[SpinArrayType, SpinArrayType]: - """Update the rotation matrix. - - Also returns the total rotation matrix. - - Args: - u_tot: Total rotation matrix. - damping: Damping object. - t1: T1 amplitude. - - Returns: - Rotation matrix and total rotation matrix. - """ - if t1 is None: - t1 = self.cc.t1 - if u_tot is None: - u_tot = util.Namespace( - aa=np.eye(self.cc.space[0].ncorr, dtype=types[float]), - bb=np.eye(self.cc.space[1].ncorr, dtype=types[float]), - ) - - t1_block: Namespace[NDArray[T]] = util.Namespace() - zocc = np.zeros((self.cc.space[0].ncocc, self.cc.space[0].ncocc)) - zvir = np.zeros((self.cc.space[0].ncvir, self.cc.space[0].ncvir)) - t1_block.aa = np.block([[zocc, -t1.aa], [np.transpose(t1.aa), zvir]]) - zocc = np.zeros((self.cc.space[1].ncocc, self.cc.space[1].ncocc)) - zvir = np.zeros((self.cc.space[1].ncvir, self.cc.space[1].ncvir)) - t1_block.bb = np.block([[zocc, -t1.bb], [np.transpose(t1.bb), zvir]]) - - u = util.Namespace(aa=scipy.linalg.expm(t1_block.aa), bb=scipy.linalg.expm(t1_block.bb)) - - u_tot.aa = u_tot.aa @ u.aa - u_tot.bb = u_tot.bb @ u.bb - if np.linalg.det(u_tot.aa) < 0: - u_tot.aa = _put( - u_tot.aa, np.ix_(np.arange(u_tot.aa.shape[0]), np.array([0])), -u_tot.aa[:, 0] - ) - if np.linalg.det(u_tot.bb) < 0: - u_tot.bb = _put( - u_tot.bb, np.ix_(np.arange(u_tot.aa.shape[0]), np.array([0])), -u_tot.bb[:, 0] - ) - - a = np.concatenate( - [np.ravel(scipy.linalg.logm(u_tot.aa)), np.ravel(scipy.linalg.logm(u_tot.bb))], axis=0 - ) - a: NDArray[T] = np.asarray(np.real(a), dtype=types[float]) - if damping is not None: - xerr = np.concatenate([t1.aa.ravel(), t1.bb.ravel()]) - a = damping(a, error=xerr) - - u_tot.aa = scipy.linalg.expm(np.reshape(a[: u_tot.aa.size], u_tot.aa.shape)) - u_tot.bb = scipy.linalg.expm(np.reshape(a[u_tot.aa.size :], u_tot.bb.shape)) - - return u, u_tot - - def transform_amplitudes( - self, u: SpinArrayType, amplitudes: Optional[Namespace[SpinArrayType]] = None - ) -> Namespace[SpinArrayType]: - """Transform the amplitudes into the Brueckner orbital basis. - - Args: - u: Rotation matrix. - amplitudes: Cluster amplitudes. - - Returns: - Transformed cluster amplitudes. - """ - if not amplitudes: - amplitudes = self.cc.amplitudes - - nocc = (self.cc.space[0].ncocc, self.cc.space[1].ncocc) - ci = {"a": u.aa[: nocc[0], : nocc[0]], "b": u.bb[: nocc[1], : nocc[1]]} - ca = {"a": u.aa[nocc[0] :, nocc[0] :], "b": u.bb[nocc[1] :, nocc[1] :]} - - # Transform T amplitudes: - for name, key, n in self.cc.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type): - for comb in util.generate_spin_combinations(n, unique=True): - args = [getattr(self.cc.amplitudes[name], comb), tuple(range(n * 2))] - for i in range(n): - args += [ci[comb[i]], (i, i + n * 2)] - for i in range(n): - args += [ca[comb[i + n]], (i + n, i + n * 3)] - args += [tuple(range(n * 2, n * 4))] - setattr(self.cc.amplitudes[name], comb, util.einsum(*args)) - - # Transform S amplitudes: - for name, key, n in self.cc.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type): - raise util.ModelNotImplemented # TODO - - # Transform U amplitudes: - for name, key, nf, nb in self.cc.ansatz.coupling_cluster_ranks(spin_type=self.spin_type): - raise util.ModelNotImplemented # TODO - - return self.cc.amplitudes - - def get_t1_norm(self, amplitudes: Optional[Namespace[SpinArrayType]] = None) -> T: - """Get the norm of the T1 amplitude. - - Args: - amplitudes: Cluster amplitudes. - - Returns: - Norm of the T1 amplitude. - """ - if not amplitudes: - amplitudes = self.cc.amplitudes - weight_a = np.linalg.norm(amplitudes["t1"].aa) - weight_b = np.linalg.norm(amplitudes["t1"].bb) - weight: T = (weight_a**2 + weight_b**2) ** 0.5 - return weight - - def mo_to_correlated( - self, mo_coeff: tuple[NDArray[T], NDArray[T]] - ) -> tuple[NDArray[T], NDArray[T]]: - """Transform the MO coefficients into the correlated basis. - - Args: - mo_coeff: MO coefficients. - - Returns: - Correlated slice of MO coefficients. - """ - return ( - mo_coeff[0][:, self.cc.space[0].correlated], - mo_coeff[1][:, self.cc.space[1].correlated], - ) - - def mo_update_correlated( - self, - mo_coeff: tuple[NDArray[T], NDArray[T]], - mo_coeff_corr: tuple[NDArray[T], NDArray[T]], - ) -> tuple[NDArray[T], NDArray[T]]: - """Update the correlated slice of a set of MO coefficients. - - Args: - mo_coeff: MO coefficients. - mo_coeff_corr: Correlated slice of MO coefficients. - - Returns: - Updated MO coefficients. - """ - space = self.cc.space - mo_coeff = ( - _put( - mo_coeff[0], - np.ix_(np.arange(mo_coeff[0].shape[0]), space[0].correlated), # type: ignore - mo_coeff_corr[0], - ), - _put( - mo_coeff[1], - np.ix_(np.arange(mo_coeff[1].shape[0]), space[1].correlated), # type: ignore - mo_coeff_corr[1], - ), - ) - return mo_coeff - - def update_coefficients( - self, - u_tot: SpinArrayType, - mo_coeff: tuple[NDArray[T], NDArray[T]], - mo_coeff_ref: tuple[NDArray[T], NDArray[T]], - ) -> tuple[NDArray[T], NDArray[T]]: - """Update the MO coefficients. - - Args: - u_tot: Total rotation matrix. - mo_coeff: New MO coefficients. - mo_coeff_ref: Reference MO coefficients. - - Returns: - Updated MO coefficients. - """ - mo_coeff_new_corr = ( - util.einsum("pi,ij->pj", mo_coeff_ref[0], u_tot.aa), - util.einsum("pi,ij->pj", mo_coeff_ref[1], u_tot.bb), - ) - mo_coeff_new = self.mo_update_correlated(mo_coeff, mo_coeff_new_corr) - return mo_coeff_new diff --git a/ebcc/opt/uopt.py b/ebcc/opt/uopt.py new file mode 100644 index 00000000..7b02c59f --- /dev/null +++ b/ebcc/opt/uopt.py @@ -0,0 +1,405 @@ +"""Unrestricted orbital-optimised coupled cluster.""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import scipy.linalg + +from ebcc import numpy as np +from ebcc import util +from ebcc.backend import _inflate, _put, ensure_scalar +from ebcc.core.precision import astype, types +from ebcc.opt.base import BaseBruecknerEBCC, BaseOptimisedEBCC, _BaseOptimisedEBCC + +if TYPE_CHECKING: + from typing import Optional + + from numpy import floating + from numpy.typing import NDArray + + from ebcc.cc.uebcc import UEBCC, ERIsInputType, SpinArrayType + from ebcc.core.damping import BaseDamping + from ebcc.util import Namespace + + T = floating + + +class _OptimisedUEBCC(_BaseOptimisedEBCC): + """Unrestricted orbital-optimised coupled cluster.""" + + # Attributes + cc: UEBCC + + def transform_amplitudes( + self, + u: SpinArrayType, + amplitudes: Optional[Namespace[SpinArrayType]] = None, + is_lambda: bool = False, + ) -> Namespace[SpinArrayType]: + """Transform the amplitudes into the orbital-optimised basis. + + Args: + u: Rotation matrix. + amplitudes: Cluster amplitudes. + is_lambda: Whether the amplitudes are Lambda amplitudes. + + Returns: + Transformed cluster amplitudes. + """ + if not amplitudes and not is_lambda: + amplitudes = self.cc.amplitudes + elif not amplitudes: + amplitudes = self.cc.lambdas + amplitudes = amplitudes.copy() + + nocc = (self.cc.space[0].ncocc, self.cc.space[1].ncocc) + ci = {"a": u.aa[: nocc[0], : nocc[0]], "b": u.bb[: nocc[1], : nocc[1]]} + ca = {"a": u.aa[nocc[0] :, nocc[0] :], "b": u.bb[nocc[1] :, nocc[1] :]} + if is_lambda: + ci, ca = ca, ci + + # Transform T amplitudes: + for name, key, n in self.cc.ansatz.fermionic_cluster_ranks(spin_type=self.spin_type): + if is_lambda: + name = name.replace("t", "l") + for comb in util.generate_spin_combinations(n, unique=True): + args = [getattr(amplitudes[name], comb), tuple(range(n * 2))] + for i in range(n): + args += [ci[comb[i]], (i, i + n * 2)] + for i in range(n): + args += [ca[comb[i + n]], (i + n, i + n * 3)] + args += [tuple(range(n * 2, n * 4))] + setattr(amplitudes[name], comb, util.einsum(*args)) + + # Transform S amplitudes: + for name, key, n in self.cc.ansatz.bosonic_cluster_ranks(spin_type=self.spin_type): + raise util.ModelNotImplemented # TODO + + # Transform U amplitudes: + for name, key, nf, nb in self.cc.ansatz.coupling_cluster_ranks(spin_type=self.spin_type): + raise util.ModelNotImplemented # TODO + + return amplitudes + + def mo_to_correlated( + self, mo_coeff: tuple[NDArray[T], NDArray[T]] + ) -> tuple[NDArray[T], NDArray[T]]: + """Transform the MO coefficients into the correlated basis. + + Args: + mo_coeff: MO coefficients. + + Returns: + Correlated slice of MO coefficients. + """ + return ( + mo_coeff[0][:, self.cc.space[0].correlated], + mo_coeff[1][:, self.cc.space[1].correlated], + ) + + def mo_update_correlated( + self, + mo_coeff: tuple[NDArray[T], NDArray[T]], + mo_coeff_corr: tuple[NDArray[T], NDArray[T]], + ) -> tuple[NDArray[T], NDArray[T]]: + """Update the correlated slice of a set of MO coefficients. + + Args: + mo_coeff: MO coefficients. + mo_coeff_corr: Correlated slice of MO coefficients. + + Returns: + Updated MO coefficients. + """ + space = self.cc.space + mo_coeff = ( + _put( + mo_coeff[0], + np.ix_(np.arange(mo_coeff[0].shape[0]), space[0].correlated), # type: ignore + mo_coeff_corr[0], + ), + _put( + mo_coeff[1], + np.ix_(np.arange(mo_coeff[1].shape[0]), space[1].correlated), # type: ignore + mo_coeff_corr[1], + ), + ) + return mo_coeff + + def update_coefficients( + self, + u_tot: SpinArrayType, + mo_coeff: tuple[NDArray[T], NDArray[T]], + mo_coeff_ref: tuple[NDArray[T], NDArray[T]], + ) -> tuple[NDArray[T], NDArray[T]]: + """Update the MO coefficients. + + Args: + u_tot: Total rotation matrix. + mo_coeff: New MO coefficients. + mo_coeff_ref: Reference MO coefficients. + + Returns: + Updated MO coefficients. + """ + mo_coeff_new_corr = ( + util.einsum("pi,ij->pj", mo_coeff_ref[0], u_tot.aa), + util.einsum("pi,ij->pj", mo_coeff_ref[1], u_tot.bb), + ) + mo_coeff_new = self.mo_update_correlated(mo_coeff, mo_coeff_new_corr) + return mo_coeff_new + + +class BruecknerUEBCC(BaseBruecknerEBCC, _OptimisedUEBCC): + """Unrestricted Brueckner-orbital coupled cluster.""" + + def get_t1_norm(self, amplitudes: Optional[Namespace[SpinArrayType]] = None) -> T: + """Get the norm of the T1 amplitude. + + Args: + amplitudes: Cluster amplitudes. + + Returns: + Norm of the T1 amplitude. + """ + if not amplitudes: + amplitudes = self.cc.amplitudes + weight_a = np.linalg.norm(amplitudes["t1"].aa) + weight_b = np.linalg.norm(amplitudes["t1"].bb) + weight: T = (weight_a**2 + weight_b**2) ** 0.5 + return weight + + def get_rotation_matrix( + self, + u_tot: Optional[SpinArrayType] = None, + damping: Optional[BaseDamping] = None, + t1: Optional[SpinArrayType] = None, + ) -> tuple[SpinArrayType, SpinArrayType]: + """Update the rotation matrix. + + Also returns the total rotation matrix. + + Args: + u_tot: Total rotation matrix. + damping: Damping object. + t1: T1 amplitudes. + + Returns: + Rotation matrix and total rotation matrix. + """ + if u_tot is None: + u_tot = util.Namespace( + aa=np.eye(self.cc.space[0].ncorr, dtype=types[float]), + bb=np.eye(self.cc.space[1].ncorr, dtype=types[float]), + ) + if t1 is None: + t1 = self.cc.t1 + + # Get the T1 blocks + t1_block: Namespace[NDArray[T]] = util.Namespace() + zocc = np.zeros((self.cc.space[0].ncocc, self.cc.space[0].ncocc)) + zvir = np.zeros((self.cc.space[0].ncvir, self.cc.space[0].ncvir)) + t1_block.aa = np.block([[zocc, -t1.aa], [np.transpose(t1.aa), zvir]]) + zocc = np.zeros((self.cc.space[1].ncocc, self.cc.space[1].ncocc)) + zvir = np.zeros((self.cc.space[1].ncvir, self.cc.space[1].ncvir)) + t1_block.bb = np.block([[zocc, -t1.bb], [np.transpose(t1.bb), zvir]]) + + # Get the orbital gradient + u = util.Namespace(aa=scipy.linalg.expm(t1_block.aa), bb=scipy.linalg.expm(t1_block.bb)) + + u_tot.aa = u_tot.aa @ u.aa + u_tot.bb = u_tot.bb @ u.bb + if np.linalg.det(u_tot.aa) < 0: + u_tot.aa = _put( + u_tot.aa, np.ix_(np.arange(u_tot.aa.shape[0]), np.array([0])), -u_tot.aa[:, 0] + ) + if np.linalg.det(u_tot.bb) < 0: + u_tot.bb = _put( + u_tot.bb, np.ix_(np.arange(u_tot.aa.shape[0]), np.array([0])), -u_tot.bb[:, 0] + ) + + a = np.concatenate( + [np.ravel(scipy.linalg.logm(u_tot.aa)), np.ravel(scipy.linalg.logm(u_tot.bb))], axis=0 + ) + a: NDArray[T] = np.asarray(np.real(a), dtype=types[float]) + if damping is not None: + xerr = np.concatenate([t1.aa.ravel(), t1.bb.ravel()]) + a = damping(a, error=xerr) + + u_tot.aa = scipy.linalg.expm(np.reshape(a[: u_tot.aa.size], u_tot.aa.shape)) + u_tot.bb = scipy.linalg.expm(np.reshape(a[u_tot.aa.size :], u_tot.bb.shape)) + + return u, u_tot + + +class OptimisedUEBCC(BaseOptimisedEBCC, _OptimisedUEBCC): + """Unrestricted optimised coupled cluster.""" + + def get_grad_norm(self, u: SpinArrayType) -> T: + """Get the norm of the gradient. + + Args: + u: Rotation matrix. + + Returns: + Norm of the gradient. + """ + grad_a = np.linalg.norm(scipy.linalg.logm(u.aa)) + grad_b = np.linalg.norm(scipy.linalg.logm(u.bb)) + grad: T = (grad_a**2 + grad_b**2) ** 0.5 + return grad + + def get_rotation_matrix( + self, + u_tot: Optional[SpinArrayType] = None, + damping: Optional[BaseDamping] = None, + eris: Optional[ERIsInputType] = None, + rdm1: Optional[SpinArrayType] = None, + rdm2: Optional[SpinArrayType] = None, + ) -> tuple[SpinArrayType, SpinArrayType]: + """Update the rotation matrix. + + Also returns the total rotation matrix. + + Args: + u_tot: Total rotation matrix. + damping: Damping object. + eris: Electronic repulsion integrals. + rdm1: One-particle reduced density matrix. + rdm2: Two-particle reduced density matrix. + + Returns: + Rotation matrix and total rotation matrix. + """ + if u_tot is None: + u_tot = util.Namespace( + aa=np.eye(self.cc.space[0].ncorr, dtype=types[float]), + bb=np.eye(self.cc.space[1].ncorr, dtype=types[float]), + ) + if rdm1 is None: + rdm1 = self.cc.make_rdm1_f(eris=eris) + if rdm2 is None: + rdm2 = self.cc.make_rdm2_f(eris=eris) + eris = self.cc.get_eris(eris=eris) + + # Get one-particle contribution to generalised Fock matrix + mo_coeff = util.Namespace( + a=self.cc.mo_coeff[0][:, self.cc.space[0].correlated], + b=self.cc.mo_coeff[1][:, self.cc.space[1].correlated], + ) + h1e_ao: NDArray[T] = np.asarray(self.cc.mf.get_hcore(), dtype=types[float]) + h1e = util.Namespace( + aa=util.einsum("pq,pi,qj->ij", h1e_ao, mo_coeff.a, mo_coeff.a), + bb=util.einsum("pq,pi,qj->ij", h1e_ao, mo_coeff.b, mo_coeff.b), + ) + f = util.Namespace( + aa=np.dot(h1e.aa, rdm1.aa), + bb=np.dot(h1e.bb, rdm1.bb), + ) + + # Get two-particle contribution to generalised Fock matrix + h2e = util.Namespace( + aaaa=eris.aaaa.xxxx, + aabb=eris.aabb.xxxx, + bbbb=eris.bbbb.xxxx, + ) + f.aa += util.einsum("prst,stqr->pq", h2e.aaaa, rdm2.aaaa) + f.aa += util.einsum("prst,qrst->pq", h2e.aabb, rdm2.aabb) + f.bb += util.einsum("prst,stqr->pq", h2e.bbbb, rdm2.bbbb) + f.bb += util.einsum("stpr,stqr->pq", h2e.aabb, rdm2.aabb) + + # Get the orbital gradient + oa, va = self.cc.space[0].xslice("o"), self.cc.space[0].xslice("v") + ob, vb = self.cc.space[1].xslice("o"), self.cc.space[1].xslice("v") + x_vo = util.Namespace( + aa=(f.aa - np.transpose(f.aa))[va, oa] / self.cc.energy_sum("vo", "aa"), + bb=(f.bb - np.transpose(f.bb))[vb, ob] / self.cc.energy_sum("vo", "bb"), + ) + x = util.Namespace( + aa=_inflate( + (self.cc.space[0].ncorr, self.cc.space[0].ncorr), + np.ix_(self.cc.space[0].xmask("v"), self.cc.space[0].xmask("o")), # type: ignore + x_vo.aa, + ), + bb=_inflate( + (self.cc.space[1].ncorr, self.cc.space[1].ncorr), + np.ix_(self.cc.space[1].xmask("v"), self.cc.space[1].xmask("o")), # type: ignore + x_vo.bb, + ), + ) + u = util.Namespace( + aa=scipy.linalg.expm(x.aa - np.transpose(x.aa)), + bb=scipy.linalg.expm(x.bb - np.transpose(x.bb)), + ) + + u_tot.aa = u_tot.aa @ u.aa + u_tot.bb = u_tot.bb @ u.bb + if np.linalg.det(u_tot.aa) < 0: + u_tot.aa = _put( + u_tot.aa, np.ix_(np.arange(u_tot.aa.shape[0]), np.array([0])), -u_tot.aa[:, 0] + ) + if np.linalg.det(u_tot.bb) < 0: + u_tot.bb = _put( + u_tot.bb, np.ix_(np.arange(u_tot.aa.shape[0]), np.array([0])), -u_tot.bb[:, 0] + ) + + a = np.concatenate( + [np.ravel(scipy.linalg.logm(u_tot.aa)), np.ravel(scipy.linalg.logm(u_tot.bb))], axis=0 + ) + a: NDArray[T] = np.asarray(np.real(a), dtype=types[float]) + if damping is not None: + a = damping(a) + + u_tot.aa = scipy.linalg.expm(np.reshape(a[: u_tot.aa.size], u_tot.aa.shape)) + u_tot.bb = scipy.linalg.expm(np.reshape(a[u_tot.aa.size :], u_tot.bb.shape)) + + return u, u_tot + + def energy( + self, + eris: Optional[ERIsInputType] = None, + rdm1: Optional[SpinArrayType] = None, + rdm2: Optional[SpinArrayType] = None, + ) -> float: + """Calculate the energy. + + Args: + eris: Electron repulsion integrals. + rdm1: One-particle reduced density matrix. + rdm2: Two-particle reduced density matrix. + + Returns: + Total energy. + """ + if rdm1 is None: + rdm1 = self.cc.make_rdm1_f(eris=eris) + if rdm2 is None: + rdm2 = self.cc.make_rdm2_f(eris=eris) + eris = self.cc.get_eris(eris=eris) + + # Get the Hamiltonian + h1e = util.Namespace( + aa=util.einsum( + "pq,pi,qj->ij", self.cc.mf.get_hcore(), self.cc.mo_coeff[0], self.cc.mo_coeff[0] + ), + bb=util.einsum( + "pq,pi,qj->ij", self.cc.mf.get_hcore(), self.cc.mo_coeff[1], self.cc.mo_coeff[1] + ), + ) + h2e = util.Namespace( + aaaa=eris.aaaa.xxxx, + aabb=eris.aabb.xxxx, + bbbb=eris.bbbb.xxxx, + ) + + # Get the energy + e_tot = util.einsum("pq,qp->", h1e.aa, rdm1.aa) + e_tot += util.einsum("pq,qp->", h1e.bb, rdm1.bb) + e_tot += util.einsum("pqrs,pqrs->", h2e.aaaa, rdm2.aaaa) * 0.5 + e_tot += util.einsum("pqrs,pqrs->", h2e.aabb, rdm2.aabb) + e_tot += util.einsum("pqrs,pqrs->", h2e.bbbb, rdm2.bbbb) * 0.5 + e_tot += self.cc.mf.energy_nuc() + + res: float = np.real(ensure_scalar(e_tot)) + return astype(res, float) diff --git a/examples/17-occsd.py b/examples/17-occsd.py new file mode 100644 index 00000000..9c7fce90 --- /dev/null +++ b/examples/17-occsd.py @@ -0,0 +1,25 @@ +""" +Example of an orbital-optimised calculation using a CCSD reference. +""" + +import numpy as np +from pyscf import gto, scf + +from ebcc import EBCC + +# Define the molecule using PySCF +mol = gto.Mole() +mol.atom = "H 0 0 0; F 0 0 1.1" +mol.basis = "cc-pvdz" +mol.build() + +# Run a RHF calculation using PySCF +mf = scf.RHF(mol) +mf.kernel() + +# Run a CCSD calculation +ccsd = EBCC(mf, ansatz="CCSD") +ccsd.kernel() + +# Run a orbital-optimised calculation using the CCSD reference +ccsd.optimised(e_tol=1e-6, t_tol=1e-5) diff --git a/tutorials/01-Introduction.ipynb b/tutorials/01-Introduction.ipynb index 2321ab96..5842a902 100644 --- a/tutorials/01-Introduction.ipynb +++ b/tutorials/01-Introduction.ipynb @@ -45,8 +45,8 @@ "text": [ "System: uname_result(system='Linux', node='ollie-desktop', release='6.8.0-49-generic', version='#49~22.04.1-Ubuntu SMP PREEMPT_DYNAMIC Wed Nov 6 17:42:15 UTC 2', machine='x86_64') Threads 1\n", "Python 3.10.12 (main, Nov 6 2024, 20:22:13) [GCC 11.4.0]\n", - "numpy 1.26.1 scipy 1.14.0 h5py 3.10.0\n", - "Date: Sun Nov 24 12:52:20 2024\n", + "numpy 2.2.0 scipy 1.14.0 h5py 3.12.1\n", + "Date: Wed Dec 11 09:25:27 2024\n", "PySCF version 2.6.2\n", "PySCF path /home/ollie/.local/lib/python3.10/site-packages/pyscf\n", "\n", @@ -68,7 +68,7 @@ "number of NR cGTOs = 32\n", "basis = aug-cc-pvdz\n", "ecp = {}\n", - "CPU time: 0.48\n", + "CPU time: 0.50\n", "\n", "\n", "******** ********\n", @@ -85,8 +85,8 @@ "SCF max_cycles = 50\n", "direct_scf = True\n", "direct_scf_tol = 1e-13\n", - "chkfile to save SCF result = /tmp/tmpj0td81lu\n", - "max_memory 4000 MB (current use 125 MB)\n", + "chkfile to save SCF result = /tmp/tmpgm7qpzzt\n", + "max_memory 4000 MB (current use 129 MB)\n", "Set gradient conv threshold to 3.16228e-05\n", "Initial guess from minao.\n", "init E= -99.6650091730287\n", @@ -114,7 +114,7 @@ { "data": { "text/plain": [ - "-100.00114790766195" + "np.float64(-100.00114790766195)" ] }, "execution_count": 1, @@ -186,15 +186,15 @@ " / _ \\| '_ \\ / __| / __|\n", " | __/| |_) || (__ | (__\n", " \\___||_.__/ \\___| \\___|\n", - " \u001b[1m1.5.0\u001b[m\u001b[m\n", + " \u001b[1m1.6.0\u001b[m\u001b[m\n", "numpy:\n", - " > Version: 1.26.1\n", + " > Version: 2.2.0\n", " > Git hash: N/A\n", "pyscf:\n", " > Version: 2.6.2\n", " > Git hash: N/A\n", "ebcc:\n", - " > Version: 1.5.0\n", + " > Version: 1.6.0\n", " > Git hash: N/A\n", "OMP_NUM_THREADS = 1\n", "\n", @@ -232,12 +232,12 @@ " 12 -0.2374225560 -100.2385704637 \u001b[32m 2.400e-09\u001b[m \u001b[31m 2.683e-08\u001b[m\n", " 13 -0.2374225572 -100.2385704648 \u001b[32m 1.178e-09\u001b[m \u001b[32m 9.110e-09\u001b[m\n", "\n", - "\u001b[32mConverged.\u001b[m\n", + "\u001b[32mConverged\u001b[m.\n", "\n", "E(corr) = -0.2374225572\n", "E(tot) = -100.2385704648\n", "\n", - "Time elapsed: 222 ms\n", + "Time elapsed: 197 ms\n", "\n", "Correlation energy: -0.2374225572\n", "Total energy: -100.2385704648\n" @@ -289,9 +289,9 @@ " 12 \u001b[31m 1.510e-08\u001b[m\n", " 13 \u001b[32m 2.888e-09\u001b[m\n", "\n", - "\u001b[32mConverged.\u001b[m\n", + "\u001b[32mConverged\u001b[m.\n", "\n", - "Time elapsed: 282 ms\n", + "Time elapsed: 259 ms\n", "\n", "\n", "Error in approximate 1RDM: 0.00842\n" @@ -338,18 +338,18 @@ "\n", "Solving for IP excitations using the Davidson solver.\n", "\n", - "\u001b[32mConverged.\u001b[m\n", + "\u001b[32mConverged\u001b[m.\n", "\n", "\u001b[1mRoot Energy Weight Conv.\u001b[m\n", - " 0 0.5573504810 0.95419 \u001b[32m True\u001b[m\n", - " 1 0.5573553537 0.95419 \u001b[32m True\u001b[m\n", - " 2 0.6780342604 0.96114 \u001b[32m True\u001b[m\n", - " 3 1.3574272581 0.00053347 \u001b[32m True\u001b[m\n", - " 4 1.3695781840 0.6841 \u001b[32m True\u001b[m\n", + " 0 0.5573552294 0.95419 \u001b[32m True\u001b[m\n", + " 1 0.5573553699 0.95419 \u001b[32m True\u001b[m\n", + " 2 0.6780342414 0.96114 \u001b[32m True\u001b[m\n", + " 3 1.3574267676 0.00053338 \u001b[32m True\u001b[m\n", + " 4 1.3574296809 0.00053366 \u001b[32m True\u001b[m\n", "\n", - "Time elapsed: 185 ms\n", + "Time elapsed: 537 ms\n", "\n", - "First ionisation potential: -0.55735048\n", + "First ionisation potential: -0.55735523\n", "Singles weight: 0.9542\n" ] } @@ -367,7 +367,7 @@ "id": "7dc63f4e-d40c-470e-8a25-e9523d01629c", "metadata": {}, "source": [ - "Brueckner orbital self-consistency is also generally available to methods. The `brueckner` method automatically calls the `kernel` itself this time, updating in-place the coefficients and amplitudes. For finer control the solvers can be used from `ebcc.opt`." + "Brueckner orbital self-consistency is also generally available to methods, as long as they have a T1 amplitude. The `brueckner` method automatically calls the `kernel` itself this time, updating in-place the coefficients and amplitudes. For finer control the solvers can be used from `ebcc.opt`." ] }, { @@ -386,31 +386,30 @@ "\n", "\u001b[1mOptions\u001b[m:\n", " > e_tol: \u001b[33m1e-08\u001b[m\n", - " > t_tol: \u001b[33m1e-08\u001b[m\n", + " > t_tol: \u001b[33m1e-07\u001b[m\n", " > max_iter: \u001b[33m20\u001b[m\n", " > diis_space: \u001b[33m9\u001b[m\n", " > diis_min_space: \u001b[33m1\u001b[m\n", " > damping: \u001b[33m0.0\u001b[m\n", "\n", - "Solving for Brueckner orbitals.\n", + "Solving for \u001b[35mBrueckner orbitals\u001b[m.\n", "\n", "\u001b[1mIter Energy (corr.) Energy (tot.) Conv. Δ(Energy) |T1|\u001b[m\n", " 0 -0.2374225572 -100.2385704648 \u001b[32m True\u001b[m\n", " 1 -0.2405424771 -100.2377076113 \u001b[31m True\u001b[m \u001b[31m 3.120e-03\u001b[m \u001b[31m 7.248e-03\u001b[m\n", " 2 -0.2396041405 -100.2378515978 \u001b[31m True\u001b[m \u001b[31m 1.440e-04\u001b[m \u001b[31m 1.930e-03\u001b[m\n", " 3 -0.2397546740 -100.2378203372 \u001b[31m True\u001b[m \u001b[31m 3.126e-05\u001b[m \u001b[31m 1.479e-05\u001b[m\n", - " 4 -0.2397542965 -100.2378205140 \u001b[31m True\u001b[m \u001b[31m 1.767e-07\u001b[m \u001b[31m 4.722e-07\u001b[m\n", - " 5 -0.2397542967 -100.2378205119 \u001b[31m True\u001b[m \u001b[32m 2.028e-09\u001b[m \u001b[31m 2.912e-08\u001b[m\n", - " 6 -0.2397542975 -100.2378205119 \u001b[31m True\u001b[m \u001b[32m 4.846e-12\u001b[m \u001b[32m 6.940e-09\u001b[m\n", + " 4 -0.2397542965 -100.2378205140 \u001b[31m True\u001b[m \u001b[31m 1.767e-07\u001b[m \u001b[31m 4.723e-07\u001b[m\n", + " 5 -0.2397542967 -100.2378205119 \u001b[31m True\u001b[m \u001b[32m 2.029e-09\u001b[m \u001b[32m 2.912e-08\u001b[m\n", "\n", - "\u001b[32mConverged.\u001b[m\n", + "\u001b[32mConverged\u001b[m.\n", "\n", - "E(corr) = -0.2397542975\n", + "E(corr) = -0.2397542967\n", "E(tot) = -100.2378205119\n", "\n", - "Time elapsed: 829 ms\n", + "Time elapsed: 660 ms\n", "\n", - "T1 norm: 6.94e-09\n" + "T1 norm: 2.91e-08\n" ] } ], @@ -420,6 +419,102 @@ "print(f\"T1 norm: {np.linalg.norm(ccsd.t1):.3g}\")" ] }, + { + "cell_type": "markdown", + "id": "a5d50be0-9910-4be0-b5a0-f08e25c02d0d", + "metadata": {}, + "source": [ + "Orbital-optimised variants can also be used for any ansatz that has fermionic reduced density matrices implemented, and behaves similarly to the Brueckner orbital self-consistency." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "da607b73-c658-4444-aa4c-f0bef1c0b49b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Solving for excitation amplitudes.\n", + "\n", + "\u001b[1mIter Energy (corr.) Energy (tot.) Δ(Energy) Δ(Ampl.)\u001b[m\n", + " 0 -0.2426301777 -100.2406963929\n", + " 1 -0.2324846749 -100.2305508901 \u001b[31m 1.015e-02\u001b[m \u001b[31m 2.026e-02\u001b[m\n", + " 2 -0.2399011124 -100.2379673276 \u001b[31m 7.416e-03\u001b[m \u001b[31m 5.390e-03\u001b[m\n", + " 3 -0.2390696180 -100.2371358332 \u001b[31m 8.315e-04\u001b[m \u001b[31m 2.505e-03\u001b[m\n", + " 4 -0.2397551480 -100.2378213633 \u001b[31m 6.855e-04\u001b[m \u001b[31m 1.406e-03\u001b[m\n", + " 5 -0.2397699213 -100.2378361365 \u001b[31m 1.477e-05\u001b[m \u001b[31m 1.397e-04\u001b[m\n", + " 6 -0.2397557461 -100.2378219613 \u001b[31m 1.418e-05\u001b[m \u001b[31m 6.782e-05\u001b[m\n", + " 7 -0.2397530628 -100.2378192781 \u001b[31m 2.683e-06\u001b[m \u001b[31m 2.423e-05\u001b[m\n", + " 8 -0.2397542098 -100.2378204251 \u001b[31m 1.147e-06\u001b[m \u001b[31m 8.766e-06\u001b[m\n", + " 9 -0.2397544044 -100.2378206197 \u001b[31m 1.946e-07\u001b[m \u001b[31m 1.565e-06\u001b[m\n", + " 10 -0.2397542810 -100.2378204962 \u001b[31m 1.235e-07\u001b[m \u001b[31m 3.180e-07\u001b[m\n", + " 11 -0.2397542895 -100.2378205047 \u001b[32m 8.490e-09\u001b[m \u001b[31m 7.512e-08\u001b[m\n", + " 12 -0.2397542964 -100.2378205117 \u001b[32m 6.975e-09\u001b[m \u001b[31m 3.000e-08\u001b[m\n", + " 13 -0.2397542958 -100.2378205110 \u001b[32m 6.721e-10\u001b[m \u001b[32m 6.867e-09\u001b[m\n", + "\n", + "\u001b[32mConverged\u001b[m.\n", + "\n", + "E(corr) = -0.2397542958\n", + "E(tot) = -100.2378205110\n", + "\n", + "Time elapsed: 185 ms\n", + "\n", + "\n", + "\u001b[1m\u001b[4mROCCSD\u001b[m\n", + "\u001b[1m******\u001b[m\n", + "\n", + "\u001b[1mOptions\u001b[m:\n", + " > e_tol: \u001b[33m1e-08\u001b[m\n", + " > t_tol: \u001b[33m1e-07\u001b[m\n", + " > max_iter: \u001b[33m20\u001b[m\n", + " > diis_space: \u001b[33m9\u001b[m\n", + " > diis_min_space: \u001b[33m1\u001b[m\n", + " > damping: \u001b[33m0.0\u001b[m\n", + "\n", + "Solving for \u001b[35moptimised orbitals\u001b[m.\n", + "\n", + "\u001b[1mIter Energy (corr.) Energy (tot.) Δ(Energy) |Gradient|\u001b[m\n", + " 0 -0.2397542958 -100.2378205110\n", + " 1 -0.2394091346 -100.2378565469 \u001b[31m 3.452e-04\u001b[m \u001b[31m 3.625e-03\u001b[m\n", + " 2 -0.2392961209 -100.2378623248 \u001b[31m 5.778e-06\u001b[m \u001b[31m 1.280e-03\u001b[m\n", + " 3 -0.2391843024 -100.2378654312 \u001b[31m 3.106e-06\u001b[m \u001b[31m 7.171e-04\u001b[m\n", + " 4 -0.2391603116 -100.2378656344 \u001b[31m 2.032e-07\u001b[m \u001b[31m 1.792e-04\u001b[m\n", + " 5 -0.2391636350 -100.2378656519 \u001b[31m 1.743e-08\u001b[m \u001b[31m 4.991e-05\u001b[m\n", + " 6 -0.2391635241 -100.2378656520 \u001b[32m 1.256e-10\u001b[m \u001b[31m 3.744e-06\u001b[m\n", + " 7 -0.2391635481 -100.2378656520 \u001b[32m 2.757e-12\u001b[m \u001b[31m 4.906e-07\u001b[m\n", + " 8 -0.2391635251 -100.2378656520 \u001b[32m 3.411e-13\u001b[m \u001b[31m 1.193e-07\u001b[m\n", + " 9 -0.2391635354 -100.2378656520 \u001b[32m 1.990e-13\u001b[m \u001b[32m 6.743e-08\u001b[m\n", + "\n", + "\u001b[32mConverged\u001b[m.\n", + "\n", + "E(corr) = -0.2391635354\n", + "E(tot) = -100.2378656520\n", + "\n", + "Time elapsed: 881 ms\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "np.float64(-0.23916353543478408)" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ccsd.amplitudes = ccsd.init_amps()\n", + "ccsd.kernel() # restore\n", + "\n", + "ccsd.optimised()" + ] + }, { "cell_type": "markdown", "id": "4d476a97-196c-44bf-910f-ecbd5e3aca44", @@ -430,7 +525,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 8, "id": "dd86aaae-82e4-47e0-8da9-2a6ec2349e8b", "metadata": {}, "outputs": [ @@ -439,11 +534,11 @@ "output_type": "stream", "text": [ "RHF:\n", - " RCC2, RCC3, RCCD, RCCSD, RCCSD(T), RCCSDT, RCCSD_SD_1_1, RCCSD_SD_1_2, RCCSD_S_1_1, RCCSDt, RCCSDt, RDCD, RDCSD, RDFCC2, RDFCCD, RDFCCSD, RDFDCD, RDFDCSD, RDFQCISD, RMP2, RMP3, RQCISD\n", + " RCC2, RCC3, RCCD, RCCSD, RCCSD(T), RCCSDT, RCCSD_SD_1_1, RCCSD_SD_1_2, RCCSD_S_1_1, RCCSDt, RCCSDt, RDCD, RDCSD, RDFCC2, RDFCCD, RDFCCSD, RDFDCD, RDFDCSD, RDFQCISD, RMP2, RMP3, RQCISD, RecCC\n", "UHF:\n", - " UCC2, UCC3, UCCD, UCCSD, UCCSD(T), UCCSDT, UCCSD_SD_1_1, UCCSD_SD_1_2, UCCSD_S_1_1, UCCSDt, UCCSDt, UDCD, UDCSD, UDFCC2, UDFCCD, UDFCCSD, UDFDCD, UDFDCSD, UDFQCISD, UMP2, UMP3, UQCISD\n", + " UCC2, UCC3, UCCD, UCCSD, UCCSD(T), UCCSDT, UCCSD_SD_1_1, UCCSD_SD_1_2, UCCSD_S_1_1, UCCSDt, UCCSDt, UDCD, UDCSD, UDFCC2, UDFCCD, UDFCCSD, UDFDCD, UDFDCSD, UDFQCISD, UMP2, UMP3, UQCISD, UecCC\n", "GHF:\n", - " GCC2, GCC3, GCCD, GCCSD, GCCSD(T), GCCSDT, GCCSDTQ, GCCSD_SD_1_1, GCCSD_SD_1_2, GCCSD_S_1_1, GCCSDt, GCCSDt, GMP2, GMP3, GQCISD\n" + " GCC2, GCC3, GCCD, GCCSD, GCCSD(T), GCCSDT, GCCSDTQ, GCCSD_SD_1_1, GCCSD_SD_1_2, GCCSD_S_1_1, GCCSDt, GCCSDt, GMP2, GMP3, GQCISD, GecCC\n" ] } ], @@ -463,7 +558,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 9, "id": "92d01dd2-19c8-41f3-910b-5f2e40d6d9b4", "metadata": {}, "outputs": [