import dataclasses
import numpy as np
import pyscf.cc
from vayesta.core.types import CCSD_WaveFunction
from vayesta.core.util import dot, log_method, log_time, einsum
from vayesta.solver._uccsd_eris import uao2mo
from vayesta.solver.cisd import CISD_Solver
from vayesta.solver.solver import ClusterSolver, UClusterSolver
[docs]class RCCSD_Solver(ClusterSolver):
[docs] @dataclasses.dataclass
class Options(ClusterSolver.Options):
# Convergence
max_cycle: int = 100 # Max number of iterations
conv_tol: float = None # Convergence energy tolerance
conv_tol_normt: float = None # Convergence amplitude tolerance
diis_space: int = None # DIIS space size
diis_start_cycle: int = None # The step to start DIIS, default is 0
iterative_damping: float = None # self-consistent damping parameter
level_shift: float = None # Level shift for virtual orbitals to stabilize ccsd iterations
init_guess: str = "MP2" # ['MP2', 'CISD']
solve_lambda: bool = True # If false, use Lambda=T approximation
n_moments: tuple = None
# Self-consistent mode
sc_mode: int = None
[docs] def kernel(self, t1=None, t2=None, l1=None, l2=None, coupled_fragments=None, t_diagnostic=True):
mf_clus, frozen = self.hamil.to_pyscf_mf(allow_dummy_orbs=True, allow_df=True)
solver_cls = self.get_solver_class(mf_clus)
self.log.debugv("PySCF solver class= %r" % solver_cls)
mycc = solver_cls(mf_clus, frozen=frozen)
if self.opts.max_cycle is not None:
mycc.max_cycle = self.opts.max_cycle
if self.opts.diis_space is not None:
mycc.diis_space = self.opts.diis_space
if self.opts.diis_start_cycle is not None:
mycc.diis_start_cycle = self.opts.diis_start_cycle
if self.opts.iterative_damping is not None:
mycc.iterative_damping = self.opts.iterative_damping
if self.opts.level_shift is not None:
mycc.level_shift = self.opts.level_shift
if self.opts.conv_tol is not None:
mycc.conv_tol = self.opts.conv_tol
if self.opts.conv_tol_normt is not None:
mycc.conv_tol_normt = self.opts.conv_tol_normt
mycc.callback = self.get_callback()
if t1 is None and t2 is None:
t1, t2 = self.generate_init_guess()
self.log.info("Solving CCSD-equations %s initial guess...", "with" if (t2 is not None) else "without")
with log_time(self.log.timing, "Time for T amplitudes: %s"):
mycc.kernel(t1=t1, t2=t2)
self.converged = mycc.converged
if t_diagnostic:
self.t_diagnostic(mycc)
self.print_extra_info(mycc)
if self.opts.solve_lambda:
with log_time(self.log.timing, "Time for lambda amplitudes: %s"):
l1, l2 = mycc.solve_lambda(l1=l1, l2=l2)
self.converged = self.converged and mycc.converged_lambda
else:
self.log.info("Using Lambda=T approximation for Lambda-amplitudes.")
l1, l2 = mycc.t1, mycc.t2
self.wf = CCSD_WaveFunction(self.hamil.mo, mycc.t1, mycc.t2, l1=l1, l2=l2)
# In-cluster Moments
nmom = self.opts.n_moments
if nmom is not None:
try:
from dyson.expressions import CCSD
except ImportError:
self.log.error("Dyson not found - required for moment calculations")
self.log.info("Skipping in-cluster moment calculations")
return
self.log.info("Calculating in-cluster CCSD moments %s" % str(nmom))
with log_time(self.log.timing, "Time for hole moments: %s"):
expr = CCSD["1h"](mf_clus, t1=mycc.t1, t2=mycc.t2, l1=l1, l2=l2)
self.hole_moments = expr.build_gf_moments(nmom[0])
# vecs_bra = expr.build_gf_vectors(nmom[0], left=True)
# amps_bra = [expr.eom.vector_to_amplitudes(amps[n,p], ccm.nmo, ccm.nocc) for p in range(ccm.nmo) for n in range(nmom)]
# vecs_ket = expr.build_gf_vectors(nmom[0], left=False)
# amps_ket = [expr.eom.vector_to_amplitudes(amps[n,p], ccm.nmo, ccm.nocc) for p in range(ccm.nmo) for n in range(nmom)]
# self.ip_moment_amplitudes = (amps_bra, amps_ket)
with log_time(self.log.timing, "Time for hole moments: %s"):
expr = CCSD["1p"](mf_clus, t1=mycc.t1, t2=mycc.t2, l1=l1, l2=l2)
self.particle_moments = expr.build_gf_moments(nmom[1])
# vecs_bra = expr.build_gf_vectors(nmom[0], left=True)
# amps_bra = [expr.eom.vector_to_amplitudes(amps[n,p], ccm.nmo, ccm.nocc) for p in range(ccm.nmo) for n in range(nmom)]
# vecs_ket = expr.build_gf_vectors(nmom[0], left=False)
# amps_ket = [expr.eom.vector_to_amplitudes(amps[n,p], ccm.nmo, ccm.nocc) for p in range(ccm.nmo) for n in range(nmom)]
# self.ea_moment_amplitudes = (amps_bra, amps_ket)
[docs] def get_solver_class(self, mf):
if hasattr(mf, "with_df") and mf.with_df is not None:
return pyscf.cc.dfccsd.RCCSD
return pyscf.cc.ccsd.RCCSD
[docs] def generate_init_guess(self, eris=None):
if self.opts.init_guess in ("default", "MP2"):
# CCSD will build MP2 amplitudes
return None, None
if self.opts.init_guess == "CISD":
cisd = CISD_Solver(self.hamil)
cisd.kernel()
wf = cisd.wf.as_ccsd()
return wf.t1, wf.t2
raise ValueError("init_guess= %r" % self.opts.init_guess)
[docs] @log_method()
def t_diagnostic(self, solver):
self.log.info("T-Diagnostic")
self.log.info("------------")
try:
dg_t1 = solver.get_t1_diagnostic()
dg_d1 = solver.get_d1_diagnostic()
dg_d2 = solver.get_d2_diagnostic()
dg_t1_msg = "good" if dg_t1 <= 0.02 else "inadequate!"
dg_d1_msg = "good" if dg_d1 <= 0.02 else ("fair" if dg_d1 <= 0.05 else "inadequate!")
dg_d2_msg = "good" if dg_d2 <= 0.15 else ("fair" if dg_d2 <= 0.18 else "inadequate!")
fmt = 3 * " %2s= %.3f (%s)"
self.log.info(fmt, "T1", dg_t1, dg_t1_msg, "D1", dg_d1, dg_d1_msg, "D2", dg_d2, dg_d2_msg)
# good: MP2~CCSD~CCSD(T) / fair: use MP2/CCSD with caution
self.log.info(" (T1<0.02: good / D1<0.02: good, D1<0.05: fair / D2<0.15: good, D2<0.18: fair)")
if dg_t1 > 0.02 or dg_d1 > 0.05 or dg_d2 > 0.18:
self.log.info(" at least one diagnostic indicates that CCSD is not adequate!")
except Exception as e:
self.log.error("Exception in T-diagnostic: %s", e)
[docs] def get_callback(self):
return None
def _debug_exact_wf(self, wf):
mo = self.hamil.mo
# Project onto cluster:
ovlp = self.hamil._fragment.base.get_ovlp()
ro = dot(wf.mo.coeff_occ.T, ovlp, mo.coeff_occ)
rv = dot(wf.mo.coeff_vir.T, ovlp, mo.coeff_vir)
t1 = dot(ro.T, wf.t1, rv)
t2 = einsum("Ii,Jj,IJAB,Aa,Bb->ijab", ro, ro, wf.t2, rv, rv)
if wf.l1 is not None:
l1 = dot(ro.T, wf.l1, rv)
l2 = einsum("Ii,Jj,IJAB,Aa,Bb->ijab", ro, ro, wf.l2, rv, rv)
else:
l1 = l2 = None
self.wf = CCSD_WaveFunction(mo, t1, t2, l1=l1, l2=l2)
self.converged = True
def _debug_random_wf(self):
mo = self.hamil.mo
t1 = np.random.rand(mo.nocc, mo.nvir)
l1 = np.random.rand(mo.nocc, mo.nvir)
t2 = np.random.rand(mo.nocc, mo.nocc, mo.nvir, mo.nvir)
l2 = np.random.rand(mo.nocc, mo.nocc, mo.nvir, mo.nvir)
self.wf = CCSD_WaveFunction(mo, t1, t2, l1=l1, l2=l2)
self.converged = True
[docs]class UCCSD_Solver(UClusterSolver, RCCSD_Solver):
[docs] @dataclasses.dataclass
class Options(RCCSD_Solver.Options):
pass
[docs] def get_solver_class(self, mf):
return UCCSD
[docs] def t_diagnostic(self, solver):
"""T diagnostic not implemented for UCCSD in PySCF."""
self.log.info("T diagnostic not implemented for UCCSD in PySCF.")
def _debug_exact_wf(self, wf):
mo = self.hamil.mo
# Project onto cluster:
ovlp = self.hamil._fragment.base.get_ovlp()
roa = dot(wf.mo.coeff_occ[0].T, ovlp, mo.coeff_occ[0])
rob = dot(wf.mo.coeff_occ[1].T, ovlp, mo.coeff_occ[1])
rva = dot(wf.mo.coeff_vir[0].T, ovlp, mo.coeff_vir[0])
rvb = dot(wf.mo.coeff_vir[1].T, ovlp, mo.coeff_vir[1])
t1a = dot(roa.T, wf.t1a, rva)
t1b = dot(rob.T, wf.t1b, rvb)
t2aa = einsum("Ii,Jj,IJAB,Aa,Bb->ijab", roa, roa, wf.t2aa, rva, rva)
t2ab = einsum("Ii,Jj,IJAB,Aa,Bb->ijab", roa, rob, wf.t2ab, rva, rvb)
t2bb = einsum("Ii,Jj,IJAB,Aa,Bb->ijab", rob, rob, wf.t2bb, rvb, rvb)
t1 = (t1a, t1b)
t2 = (t2aa, t2ab, t2bb)
if wf.l1 is not None:
l1a = dot(roa.T, wf.l1a, rva)
l1b = dot(rob.T, wf.l1b, rvb)
l2aa = einsum("Ii,Jj,IJAB,Aa,Bb->ijab", roa, roa, wf.l2aa, rva, rva)
l2ab = einsum("Ii,Jj,IJAB,Aa,Bb->ijab", roa, rob, wf.l2ab, rva, rvb)
l2bb = einsum("Ii,Jj,IJAB,Aa,Bb->ijab", rob, rob, wf.l2bb, rvb, rvb)
l1 = (l1a, l1b)
l2 = (l2aa, l2ab, l2bb)
else:
l1 = l2 = None
self.wf = CCSD_WaveFunction(mo, t1, t2, l1=l1, l2=l2)
self.converged = True
def _debug_random_wf(self):
raise NotImplementedError
# Subclass pyscf UCCSD to enable support of spin-dependent ERIs.
[docs]class UCCSD(pyscf.cc.uccsd.UCCSD):
ao2mo = uao2mo