Source code for pyscf.cc.momgfccsd

#!/usr/bin/env python
# Copyright 2014-2021 The PySCF Developers. All Rights Reserved.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Author: Oliver Backhouse <olbackhouse@gmail.com>
#

"""
GF-CCSD solver via moment constraints.

See reference: Backhouse, Booth, arXiv:2206.13198 (2022).
"""

from collections import defaultdict

import numpy as np
import scipy.linalg

from pyscf import lib, cc, ao2mo
from pyscf.lib import logger
from pyscf.agf2 import mpi_helper


[docs] def kernel( gfccsd, hole_moments=None, part_moments=None, t1=None, t2=None, l1=None, l2=None, eris=None, imds=None, verbose=None, ): if gfccsd.verbose >= logger.WARN: gfccsd.check_sanity() gfccsd.dump_flags() log = logger.new_logger(gfccsd, verbose) if (l1 is None and gfccsd._cc.l1 is None) or (l2 is None and gfccsd._cc.l2 is None): raise ValueError( "Lambda amplitudes must be set for %s. This " "can be done by calling solve_lambda on the " "CC method or by setting l1, l2 attributes. " % gfccsd.__class__.__name__ ) if (hole_moments is None or part_moments is None) and imds is None: ip = hole_moments is None ea = part_moments is None imds = gfccsd.make_imds(eris=eris, ip=ip, ea=ea) if hole_moments is None: log.info("Building hole moments:") hole_moments = gfccsd.build_hole_moments(t1=t1, t2=t2, l1=l1, l2=l2, imds=imds) else: log.info("Hole moments passed by argument.") if part_moments is None: log.info("Building particle moments:") part_moments = gfccsd.build_part_moments(t1=t1, t2=t2, l1=l1, l2=l2, imds=imds) else: log.info("Particle moments passed by argument.") if gfccsd.hermi_moments: hole_moments = 0.5 * (hole_moments + hole_moments.swapaxes(1, 2).conj()) part_moments = 0.5 * (part_moments + part_moments.swapaxes(1, 2).conj()) if gfccsd.hermi_solver: solver = block_lanczos_symm eig = eigh_block_tridiagonal else: solver = block_lanczos_nosymm eig = eig_block_tridiagonal log.info("Solving for the hole moments.") blocks = solver(gfccsd, hole_moments) orth = mat_sqrt(hole_moments[0], hermi=gfccsd.hermi_solver) eh, vh = eig(gfccsd, *blocks, orth=orth) log.info("Solving for the particle moment.") blocks = solver(gfccsd, part_moments) orth = mat_sqrt(part_moments[0], hermi=gfccsd.hermi_solver) ep, vp = eig(gfccsd, *blocks, orth=orth) # Check the moments if gfccsd.niter[0] is not None: for n in range(2*gfccsd.niter[0]+2): a = lib.einsum("xk,yk,k->xy", vh[0], vh[1].conj(), eh**n) a /= np.max(np.abs(a)) b = hole_moments[n] / np.max(np.abs(hole_moments[n])) err = np.max(np.abs(a - b)) (logger.debug1 if err < 1e-8 else logger.warn)( gfccsd, "Error in hole moment %d: %10.6g", n, err) if gfccsd.niter[0] is not None: for n in range(2*gfccsd.niter[1]+2): a = lib.einsum("xk,yk,k->xy", vp[0], vp[1].conj(), ep**n) a /= np.max(np.abs(a)) b = part_moments[n] / np.max(np.abs(part_moments[n])) err = np.max(np.abs(a - b)) (logger.debug1 if err < 1e-8 else logger.warn)( gfccsd, "Error in particle moment %d: %10.6g", n, err) mask = np.argsort(eh.real) eh, vh = eh[mask], (vh[0][:, mask], vh[1][:, mask]) mask = np.argsort(ep.real) ep, vp = ep[mask], (vp[0][:, mask], vp[1][:, mask]) return eh, vh, ep, vp
[docs] def mat_sqrt(m, hermi=False): """Return the square root of a matrix. """ if hermi: w, v = np.linalg.eigh(m) mask = w >= 0 w, v = w[mask], v[:, mask] out = np.dot(v * w[None]**0.5, v.T.conj()) else: w, v = np.linalg.eig(m) out = np.dot(v * w[None]**(0.5+0j), np.linalg.inv(v)) return out
[docs] def mat_isqrt(m, tol=1e-16, hermi=False): """Return the inverse square root of a matrix. """ if hermi: w, v = np.linalg.eigh(m) mask = w > tol w, v = w[mask], v[:, mask] out = np.dot(v * w[None]**-0.5, v.T.conj()) else: w, v = np.linalg.eig(m) mask = np.abs(w) >= tol vinv = np.linalg.inv(v)[mask] w, v = w[mask], v[:, mask] out = np.dot(v * w[None]**(-0.5+0j), vinv) return out
[docs] def build_block_tridiagonal(a, b, c=None): """Construct a block tridiagonal matrix from a list of on-diagonal and off-diagonal blocks. """ z = np.zeros_like(a[0], dtype=a[0].dtype) if c is None: c = [x.T.conj() for x in b] h = np.block([[ a[i] if i == j else b[j] if j == i-1 else c[i] if i == j-1 else z for j in range(len(a))] for i in range(len(a))] ) return h
[docs] def eig_block_tridiagonal(gfccsd, a, b, c, orth=None): """Diagonalise a non-Hermitian block-tridiagonal Hamiltonian and transform its eigenvectors appropriately. """ h_tri = build_block_tridiagonal(a, b, c) e, u = np.linalg.eig(h_tri) if orth is not None: vl = np.dot(orth, u[:gfccsd.nmo]) vr = np.dot(np.linalg.inv(u)[:, :gfccsd.nmo], orth).T.conj() else: vl = u[:gfccsd.nmo] vr = np.linalg.inv(u)[:, :gfccsd.nmo].T.conj() return e, (vl, vr)
[docs] def eigh_block_tridiagonal(gfccsd, a, b, orth=None): """Diagonalise a Hermitian block-tridiagonal Hamiltonian and transform its eigenvectors appropriately. """ h_tri = build_block_tridiagonal(a, b) e, u = np.linalg.eigh(h_tri) if orth is not None: v = np.dot(orth, u[:gfccsd.nmo]) else: v = u[:gfccsd.nmo] return e, (v, v)
def _matrix_info(x, hermi=False): norm = np.abs(np.einsum("pq,qp->", x, x)) eigvals = np.linalg.eigvals(x) mineig = np.min(np.abs(eigvals)) maxeig = np.max(np.abs(eigvals)) return norm, mineig, maxeig
[docs] def block_lanczos_symm(gfccsd, moments, verbose=None): """Hermitian block Lanczos solver, returns a set of poles that best reproduce the inputted moments. Args: gfccsd : MomGFCCSD GF-CCSD object moments : ndarray (2*niter+2, n, n) Array of moments with which the resulting poles should be consistent with. Kwargs: verbose : int Level of verbosity. Returns: a : ndarray (niter+1, n, n) On-diagonal blocks of the block tridiagonal Hamiltonian. b : ndarray (niter, n, n) Off-diagonal blocks of the block tridiagonal Hamiltonian. """ log = logger.new_logger(gfccsd, verbose) log.debug1("block_lanczos_symm: %d moments", len(moments)) nmo = gfccsd.nmo niter = (len(moments) - 2) // 2 dtype = np.complex128 a = np.zeros((niter+1, nmo, nmo), dtype=dtype) b = np.zeros((niter, nmo, nmo), dtype=dtype) t = np.zeros((len(moments), nmo, nmo), dtype=dtype) v = defaultdict(lambda: np.zeros((nmo, nmo), dtype=dtype)) v[0, 0] = np.eye(nmo).astype(dtype) orth = mat_isqrt(moments[0], hermi=True) for i in range(len(moments)): t[i] = np.linalg.multi_dot((orth, moments[i], orth)) a[0] = t[1] log.debug1("Raw moments:") log.debug1(" %4s %12s %12s %12s", "N", "norm", "min(|eig|)", "max(|eig|)") for i in range(len(moments)): log.debug1(" %4d %12.6g %12.6g %12.6g", i, *_matrix_info(moments[i], hermi=True)) log.debug1("Orthogonalised moments:") log.debug1(" %4s %12s %12s %12s", "N", "norm", "min(|eig|)", "max(|eig|)") for i in range(len(moments)): log.debug1(" %4d %12.6g %12.6g %12.6g", i, *_matrix_info(t[i], hermi=True)) for i in range(niter): log.info("Iteration %d", i) b2 = np.zeros((nmo, nmo), dtype=dtype) for j in range(i+2): for l in range(i+1): b2 += np.linalg.multi_dot((v[i, l].T.conj(), t[j+l+1], v[i, j-1])) b2 -= np.dot(a[i], a[i]) if i: b2 -= np.dot(b[i-1], b[i-1]) b[i] = mat_sqrt(b2, hermi=True) binv = mat_isqrt(b2, hermi=True) for j in range(i+2): r = ( + v[i, j-1] - np.dot(v[i, j], a[i]) - np.dot(v[i-1, j], b[i-1]) ) v[i+1, j] = np.dot(r, binv) for j in range(i+2): for l in range(i+2): a[i+1] += np.linalg.multi_dot((v[i+1, l].T.conj(), t[j+l+1], v[i+1, j])) log.debug1(" %4s %12s %12s %12s", "mat", "norm", "min(|eig|)", "max(|eig|)") log.debug1(" %4s %12.6g %12.6g %12.6g", "B^2", *_matrix_info(b2)) log.debug1(" %4s %12.6g %12.6g %12.6g", "B", *_matrix_info(b[i])) log.debug1(" %4s %12.6g %12.6g %12.6g", "B^-1", *_matrix_info(binv)) log.debug1(" %4s %12.6g %12.6g %12.6g", "A", *_matrix_info(a[i+1])) biorth_error = 0.0 for j in range(i+2): x = np.zeros_like(v[0, 0]) for k in range(i+2): for l in range(i+2): x += np.linalg.multi_dot((v[i+1, l].T.conj(), t[k+l], v[j, k])) biorth_error = max(biorth_error, np.max(np.abs(x - np.eye(nmo)*((i+1)==j)))) log.info(" Error in biorthogonality: %12.6g", biorth_error) return a, b
[docs] def block_lanczos_nosymm(gfccsd, moments, verbose=None): """Non-Hermitian block Lanczos solver, returns a set of poles that best reproduce the inputted moments. Args: gfccsd : MomGFCCSD GF-CCSD object moments : ndarray (2*niter+2, n, n) Array of moments with which the resulting poles should be consistent with. Kwargs: verbose : int Level of verbosity. Returns: a : ndarray (niter+1, n, n) On-diagonal blocks of the block tridiagonal Hamiltonian. b : ndarray (niter, n, n) Upper off-diagonal blocks of the block tridiagonal Hamiltonian. c : ndarray (niter, n, n) Lower off-diagonal blocks of the block tridiagonal Hamiltonian. """ log = logger.new_logger(gfccsd, verbose) log.debug1("block_lanczos_nosymm: %d moments", len(moments)) nmo = gfccsd.nmo niter = (len(moments) - 2) // 2 dtype = np.complex128 a = np.zeros((niter+1, nmo, nmo), dtype=dtype) b = np.zeros((niter, nmo, nmo), dtype=dtype) c = np.zeros((niter, nmo, nmo), dtype=dtype) t = np.zeros((len(moments), nmo, nmo), dtype=dtype) v = defaultdict(lambda: np.zeros((nmo, nmo), dtype=dtype)) w = defaultdict(lambda: np.zeros((nmo, nmo), dtype=dtype)) v[0, 0] = np.eye(nmo).astype(dtype) w[0, 0] = np.eye(nmo).astype(dtype) orth = mat_isqrt(moments[0]) for i in range(len(moments)): t[i] = np.linalg.multi_dot((orth, moments[i], orth)) a[0] = t[1] log.debug1("Raw moments:") log.debug1(" %4s %12s %12s %12s", "N", "norm", "min(|eig|)", "max(|eig|)") for i in range(len(moments)): log.debug1(" %4d %12.6g %12.6g %12.6g", i, *_matrix_info(moments[i])) log.debug1("Orthogonalised moments:") log.debug1(" %4s %12s %12s %12s", "N", "norm", "min(|eig|)", "max(|eig|)") for i in range(len(moments)): log.debug1(" %4d %12.6g %12.6g %12.6g", i, *_matrix_info(t[i])) for i in range(niter): log.info("Iteration %d", i) b2 = np.zeros((nmo, nmo), dtype=dtype) c2 = np.zeros((nmo, nmo), dtype=dtype) for j in range(i+2): for l in range(i+1): b2 += np.linalg.multi_dot((w[i, l], t[j+l+1], v[i, j-1])) c2 += np.linalg.multi_dot((w[i, j-1], t[j+l+1], v[i, l])) b2 -= np.dot(a[i], a[i]) c2 -= np.dot(a[i], a[i]) if i: b2 -= np.dot(c[i-1], c[i-1]) c2 -= np.dot(b[i-1], b[i-1]) b[i] = mat_sqrt(b2) c[i] = mat_sqrt(c2) binv = mat_isqrt(b2) cinv = mat_isqrt(c2) for j in range(i+2): r = ( + v[i, j-1] - np.dot(v[i, j], a[i]) - np.dot(v[i-1, j], b[i-1]) ) v[i+1, j] = np.dot(r, cinv) s = ( + w[i, j-1] - np.dot(a[i], w[i, j]) - np.dot(c[i-1], w[i-1, j]) ) w[i+1, j] = np.dot(binv, s) for j in range(i+2): for l in range(i+2): a[i+1] += np.linalg.multi_dot((w[i+1, l], t[j+l+1], v[i+1, j])) log.debug1(" %4s %12s %12s %12s", "mat", "norm", "min(|eig|)", "max(|eig|)") log.debug1(" %4s %12.6g %12.6g %12.6g", "B^2", *_matrix_info(b2)) log.debug1(" %4s %12.6g %12.6g %12.6g", "B", *_matrix_info(b[i])) log.debug1(" %4s %12.6g %12.6g %12.6g", "B^-1", *_matrix_info(binv)) log.debug1(" %4s %12.6g %12.6g %12.6g", "C^2", *_matrix_info(c2)) log.debug1(" %4s %12.6g %12.6g %12.6g", "C", *_matrix_info(c[i])) log.debug1(" %4s %12.6g %12.6g %12.6g", "C^-1", *_matrix_info(cinv)) log.debug1(" %4s %12.6g %12.6g %12.6g", "A", *_matrix_info(a[i+1])) biorth_error = 0.0 for j in range(i+2): x = np.zeros_like(v[0, 0]) y = np.zeros_like(v[0, 0]) for k in range(i+2): for l in range(i+2): x += np.linalg.multi_dot((w[i+1, l], t[k+l], v[j, k])) y += np.linalg.multi_dot((w[j, l], t[k+l], v[i+1, k])) biorth_error = max(biorth_error, np.max(np.abs(x - np.eye(nmo)*((i+1)==j)))) biorth_error = max(biorth_error, np.max(np.abs(y - np.eye(nmo)*((i+1)==j)))) log.info(" Error in biorthogonality: %12.6g", biorth_error) return a, b, c
def _kd(n, i): v = np.zeros((n,)) v[i] = 1.0 return v
[docs] def contract_ket_hole(gfccsd, eom, t1, t2, v, orb): r"""Contract a vector with \bar{a}^\dagger_p |\Psi>. """ nocc, nvir = t1.shape if orb < nocc: return v[orb] else: b1 = t1[:, orb-nocc] b2 = t2[:, :, orb-nocc] b = eom.amplitudes_to_vector(b1, b2) return np.dot(v, b)
[docs] def build_ket_hole(gfccsd, eom, t1, t2, orb): r"""Build \bar{a}^\dagger_p |\Psi>. """ nocc, nvir = t1.shape if orb < nocc: b1 = np.eye(nocc)[orb] b2 = np.zeros((nocc, nocc, nvir)) else: b1 = t1[:, orb-nocc] b2 = t2[:, :, orb-nocc] return eom.amplitudes_to_vector(b1, b2)
[docs] def build_bra_hole(gfccsd, eom, t1, t2, l1, l2, orb): """Get the first- and second-order contributions to the left-hand transformed vector for a given orbital for the hole part of the Green's function. """ nocc, nvir = t1.shape if orb < nocc: e1 = _kd(nocc, orb) e1 -= lib.einsum("ie,e->i", l1, t1[orb]) tmp = t2[orb] * 2.0 tmp -= t2[orb].swapaxes(1, 2) e1 -= lib.einsum("imef,mef->i", l2, tmp) tmp = -lib.einsum("ijea,e->ija", l2, t1[orb]) e2 = 2.0 * tmp e2 -= tmp.swapaxes(0, 1) tmp = lib.einsum("ja,i->ija", l1, _kd(nocc, orb)) e2 += tmp * 2.0 e2 -= tmp.swapaxes(0, 1) else: e1 = l1[:, orb-nocc].copy() e2 = l2[:, :, orb-nocc] * 2.0 e2 -= l2[:, :, :, orb-nocc] return eom.amplitudes_to_vector(e1, e2)
[docs] def contract_ket_part(gfccsd, eom, t1, t2, v, orb): r"""Contract a vector with \bar{a}_p |\Psi>. """ nocc, nvir = t1.shape if orb < nocc: b1 = t1[orb] b2 = t2[orb] b = eom.amplitudes_to_vector(b1, b2) return np.dot(v, b) else: return -v[orb-nocc]
[docs] def build_ket_part(gfccsd, eom, t1, t2, orb): r"""Build \bar{a}_p |\Psi>. """ nocc, nvir = t1.shape if orb < nocc: b1 = t1[orb] b2 = t2[orb] else: b1 = -np.eye(nvir)[orb-nocc] b2 = np.zeros((nocc, nvir, nvir)) return eom.amplitudes_to_vector(b1, b2)
[docs] def build_bra_part(gfccsd, eom, t1, t2, l1, l2, orb): """Get the first- and second-order contributions to the left-hand transformed vector for a given orbital for the particle part of the Green's function. """ nocc, nvir = t1.shape if orb < nocc: e1 = -l1[orb] e2 = -l2[orb] * 2.0 e2 += l2[:, orb] else: e1 = _kd(nvir, orb-nocc) e1 -= lib.einsum("mb,m->b", l1, t1[:, orb-nocc]) tmp = t2[:, :, :, orb-nocc] * 2.0 tmp -= t2[:, :, orb-nocc] e1 -= lib.einsum("kmeb,kme->b", l2, tmp) tmp = -lib.einsum("ikba,k->iab", l2, t1[:, orb-nocc]) e2 = tmp * 2.0 e2 -= tmp.swapaxes(1, 2) tmp = lib.einsum("ib,a->iab", l1, _kd(nvir, orb-nocc)) e2 += tmp * 2.0 e2 -= tmp.swapaxes(1, 2) return eom.amplitudes_to_vector(e1, e2)
[docs] class MomGFCCSD(lib.StreamObject): """Green's function coupled cluster singles and doubles using the moment-resolved solver. Attributes: verbose : int Print level. Default value equals to :class:`Mole.verbose`. niter : tuple of (int, int) Number of block Lanczos iterations for occupied and virtual sectors. If either are `None` then said sector will not be computed. weight_tol : float Threshold for weight in the physical space to consider a pole an ionisation or removal event. Default value is 1e-1. hermi_moments : bool Whether to Hermitise the moments, default value is False. hermi_solver : obol Whether to use the real-valued, symmetric block Lanczos solver, default value is False. Results: eh : ndarray Energies of the compressed hole Green's function vh : tuple of ndarray Left- and right-hand transition amplitudes of the compressed hole Green's function ep : ndarray Energies of the compressed particle Green's function vp : tuple of ndarray Left- and right-hand transition amplitudes of the compressed particle Green's function """ _keys = { 'verbose', 'stdout', 'niter', 'weight_tol', 'hermi_moments', 'hermi_solver', 'eh', 'ep', 'vh', 'vp', 'chkfile', } def __init__(self, mycc, niter=(2, 2)): self._cc = mycc self.verbose = mycc.verbose self.stdout = mycc.stdout if isinstance(mycc, cc.uccsd.UCCSD): raise NotImplementedError("MomGFCCSD for unrestricted CCSD") if isinstance(niter, int): self.niter = (niter, niter) else: self.niter = niter self.weight_tol = 1e-1 self.hermi_moments = False self.hermi_solver = False self.eh = None self.ep = None self.vh = None self.vp = None self._t1 = None self._t2 = None self._l1 = None self._l2 = None self.chkfile = self._cc.chkfile
[docs] def dump_flags(self, verbose=None): log = logger.new_logger(self, verbose) log.info("") log.info("******** %s ********", self.__class__) log.info("niter = %s", self.niter) log.info("nmo = %s", self.nmo) log.info("nocc = %s", self.nocc) log.info("weight_tol = %s", self.weight_tol) log.info("hermi_moments = %s", self.hermi_moments) log.info("hermi_solver = %s", self.hermi_solver) log.info("chkfile = %s", self.chkfile)
def _finalize(self): self.ipgfccsd() self.eagfccsd() return self
[docs] def reset(self, mol=None): self._cc.reset(mol) return self
@property def eomip_method(self): return self._cc.eomip_method() @property def eomea_method(self): return self._cc.eomea_method() build_bra_hole = build_bra_hole build_bra_part = build_bra_part contract_ket_hole = contract_ket_hole contract_ket_part = contract_ket_part
[docs] def make_imds(self, eris=None, ip=True, ea=True): """Build EOM intermediates. """ imds = cc.eom_rccsd._IMDS(self._cc, eris=eris) if ip: imds.make_ip() if ea: imds.make_ea() return imds
[docs] def build_hole_moments(self, t1=None, t2=None, l1=None, l2=None, imds=None, niter=None): """Build moments of the hole (IP-EOM-CCSD) Green's function. """ if t1 is None: t1 = self._cc.t1 if t2 is None: t2 = self._cc.t2 if l1 is None: l1 = self._cc.l1 if l2 is None: l2 = self._cc.l2 if niter is None: niter = self.niter[0] nmom = 2 * niter + 2 moments = np.zeros((nmom, self.nmo, self.nmo)) cput0 = (logger.process_clock(), logger.perf_counter()) eom = self.eomip_method() if imds is None: imds = self.make_imds(ea=False) diag = eom.get_diag(imds) for p in mpi_helper.nrange(self.nmo): ket = self.build_bra_hole(eom, t1, t2, l1, l2, p) for n in range(nmom): for q in range(self.nmo): moments[n, q, p] += self.contract_ket_hole(eom, t1, t2, ket, q) if (n+1) != nmom: ket = -eom.l_matvec(ket, imds, diag) mpi_helper.barrier() moments = mpi_helper.allreduce(moments) logger.timer(self, "IP-EOM-CCSD moments", *cput0) return moments
[docs] def build_part_moments(self, t1=None, t2=None, l1=None, l2=None, imds=None, niter=None): """Build moments of the particle (EA-EOM-CCSD) Green's function. """ if t1 is None: t1 = self._cc.t1 if t2 is None: t2 = self._cc.t2 if l1 is None: l1 = self._cc.l1 if l2 is None: l2 = self._cc.l2 if niter is None: niter = self.niter[1] nmom = 2 * niter + 2 moments = np.zeros((nmom, self.nmo, self.nmo)) cput0 = (logger.process_clock(), logger.perf_counter()) eom = self.eomea_method() if imds is None: imds = self.make_imds(ip=False) diag = eom.get_diag(imds) for p in mpi_helper.nrange(self.nmo): ket = self.build_bra_part(eom, t1, t2, l1, l2, p) for n in range(nmom): for q in range(self.nmo): moments[n, q, p] -= self.contract_ket_part(eom, t1, t2, ket, q) if (n+1) != nmom: ket = eom.l_matvec(ket, imds, diag) mpi_helper.barrier() moments = mpi_helper.allreduce(moments) logger.timer(self, "EA-EOM-CCSD moments", *cput0) return moments
[docs] def make_rdm1(self, ao_repr=False, eris=None, imds=None): """Build the first-order reduced density matrix at the CCSD level using the zeroth-order moment of the hole part of the CCSD Green's function. """ if imds is None: imds = self.make_imds(eris=eris, ea=False) dm1 = self.build_hole_moments(imds=imds, niter=0)[0] dm1 = dm1 + dm1.T.conj() if ao_repr: mo = self._cc.mo_coeff dm1 = np.linalg.multi_dot((mo, dm1, mo.T.conj())) return dm1
[docs] def kernel(self, **kwargs): eh, vh, ep, vp = kernel(self, **kwargs) self.eh = eh self.vh = vh self.ep = ep self.vp = vp self._finalize() return eh, vh, ep, vp
[docs] def dump_chk(self, chkfile=None, key="gfccsd"): if chkfile is None: chkfile = self.chkfile lib.chkfile.dump(chkfile, key+"/eh", self.eh) lib.chkfile.dump(chkfile, key+"/vh_left", self.vh[0]) lib.chkfile.dump(chkfile, key+"/vh_right", self.vh[1]) lib.chkfile.dump(chkfile, key+"/ep", self.ep) lib.chkfile.dump(chkfile, key+"/vp_left", self.vp[0]) lib.chkfile.dump(chkfile, key+"/vp_right", self.vp[1]) lib.chkfile.dump(chkfile, key+"/niter", np.array(self.niter)) return self
[docs] def update_from_chk_(self, chkfile=None, key="gfccsd"): if chkfile is None: chkfile = self.chkfile self.eh = lib.chkfile.load(chkfile, key+"/eh") self.vh = ( lib.chkfile.load(chkfile, key+"/vh_left"), lib.chkfile.load(chkfile, key+"/vh_right"), ) self.ep = lib.chkfile.load(chkfile, key+"/ep") self.vp = ( lib.chkfile.load(chkfile, key+"/vp_left"), lib.chkfile.load(chkfile, key+"/vp_right"), ) self.niter = tuple(lib.chkfile.load(chkfile, key+"/niter"))
update = update_from_chk = update_from_chk_
[docs] def ipgfccsd(self, nroots=5): """Print and return ionisation potentials. """ eh, (vh, uh) = self.eh, self.vh mask = np.abs(np.sum(vh * uh.conj(), axis=0)) > self.weight_tol mask = np.arange(mask.size)[mask][::-1] e_ip = -eh[mask] v_ip, u_ip = vh[:, mask], uh[:, mask] nroots = min(nroots, len(e_ip)) logger.note(self, " %s %s %16s %10s", "", "", "Energy", "Weight") for n in range(nroots): qpwt = np.abs(np.sum(v_ip[:, n] * u_ip[:, n].conj())).real warn = "" if np.abs(e_ip[n].imag) > 1e-8: warn += "(Warning: imag part: %.6g)" % e_ip[n].imag logger.note(self, " %2s %2d %16.10g %10.6g %s" % ("IP", n, e_ip[n].real, qpwt, warn)) if nroots == 1: return e_ip[0].real, v_ip[:, 0], u_ip[:, 0] else: return e_ip.real, v_ip, u_ip
[docs] def eagfccsd(self, nroots=5): """Print and return electron affinities. """ ep, (vp, up) = self.ep, self.vp mask = np.abs(np.sum(vp * up.conj(), axis=0)) > self.weight_tol e_ea = ep[mask] v_ea, u_ea = vp[:, mask], up[:, mask] nroots = min(nroots, len(e_ea)) logger.note(self, " %s %s %16s %10s", "", "", "Energy", "Weight") for n in range(nroots): qpwt = np.abs(np.sum(v_ea[:, n] * u_ea[:, n].conj())).real warn = "" if np.abs(e_ea[n].imag) > 1e-8: warn += "(Warning: imag part: %.6g)" % e_ea[n].imag logger.note(self, " %2s %2d %16.10g %10.6g %s" % ("EA", n, e_ea[n].real, qpwt, warn)) if nroots == 1: return e_ea[0].real, v_ea[:, 0], u_ea[:, 0] else: return e_ea.real, v_ea, u_ea
@property def nmo(self): return self._cc.nmo @property def nocc(self): return self._cc.nocc
if __name__ == "__main__": from pyscf import gto, scf mol = gto.M( #atom="O 0 0 0; O 0 0 1", atom="N 0 0 0; N 0 0 1", basis="cc-pvdz", verbose=0, ) mf = scf.RHF(mol) mf = mf.run() ccsd = cc.CCSD(mf) ccsd = ccsd.run() ccsd.solve_lambda() niter = 5 gfcc = MomGFCCSD(ccsd, (niter, niter)) gfcc.kernel() ip1, vip1 = ccsd.ipccsd(nroots=8) ip2, vip2, uip2 = gfcc.ipgfccsd(nroots=8) ea1, vea1 = ccsd.eaccsd(nroots=8) ea2, vea2, uea2 = gfcc.eagfccsd(nroots=8) print(" %12s %12s %12s" % ("EOM", "GF", "Error")) print("IP1 %12.8f %12.8f %12.8f" % (ip1[0],ip2[0],np.abs(ip1[0]-ip2[0]))) print("IP2 %12.8f %12.8f %12.8f" % (ip1[1],ip2[1],np.abs(ip1[1]-ip2[1]))) print("EA1 %12.8f %12.8f %12.8f" % (ea1[0],ea2[0],np.abs(ea1[0]-ea2[0]))) print("EA2 %12.8f %12.8f %12.8f" % (ea1[1],ea2[1],np.abs(ea1[1]-ea2[1])))