Source code for pyscf.agf2.uagf2

# Copyright 2014-2020 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 J. Backhouse <olbackhouse@gmail.com>
#         George H. Booth <george.booth@kcl.ac.uk>
#

'''
Auxiliary second-order Green's function perturbation theory for
unrestricted references
'''

import time
import numpy as np
from pyscf import lib
from pyscf.lib import logger
from pyscf import __config__
from pyscf import ao2mo
from pyscf.scf import _vhf
from pyscf.agf2 import ragf2, _agf2, mpi_helper
from pyscf.agf2 import aux_space as aux
from pyscf.agf2.chempot import binsearch_chempot, minimize_chempot
from pyscf.mp.ump2 import get_frozen_mask as _get_frozen_mask

BLKMIN = getattr(__config__, 'agf2_blkmin', 1)


[docs] def build_se_part(agf2, eri, gf_occ, gf_vir, os_factor=1.0, ss_factor=1.0): ''' Builds either the auxiliaries of the occupied self-energy, or virtual if :attr:`gf_occ` and :attr:`gf_vir` are swapped, for a single spin. Args: eri : _ChemistsERIs Electronic repulsion integrals gf_occ : tuple of GreensFunction Occupied Green's function for each spin gf_vir : tuple of GreensFunction Virtual Green's function for each spin Kwargs: os_factor : float Opposite-spin factor for spin-component-scaled (SCS) calculations. Default 1.0 ss_factor : float Same-spin factor for spin-component-scaled (SCS) calculations. Default 1.0 Returns: :class:`SelfEnergy` ''' cput0 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(agf2.stdout, agf2.verbose) assert type(gf_occ[0]) is aux.GreensFunction assert type(gf_occ[1]) is aux.GreensFunction assert type(gf_vir[0]) is aux.GreensFunction assert type(gf_vir[1]) is aux.GreensFunction nmo = eri.nmo noa, nob = gf_occ[0].naux, gf_occ[1].naux nva, nvb = gf_vir[0].naux, gf_vir[1].naux tol = agf2.weight_tol facs = {'os_factor': os_factor, 'ss_factor': ss_factor} ci_a, ei_a = gf_occ[0].coupling, gf_occ[0].energy ci_b, ei_b = gf_occ[1].coupling, gf_occ[1].energy ca_a, ea_a = gf_vir[0].coupling, gf_vir[0].energy ca_b, ea_b = gf_vir[1].coupling, gf_vir[1].energy mem_incore = (nmo[0]*noa*(noa*nva+nob*nvb)) * 8/1e6 mem_now = lib.current_memory()[0] if (mem_incore+mem_now < agf2.max_memory) or agf2.incore_complete: qeri = _make_qmo_eris_incore(agf2, eri, (ci_a, ci_a, ca_a), (ci_b, ci_b, ca_b), spin=0) else: qeri = _make_qmo_eris_outcore(agf2, eri, (ci_a, ci_a, ca_a), (ci_b, ci_b, ca_b), spin=0) if isinstance(qeri[0], np.ndarray): vv, vev = _agf2.build_mats_uagf2_incore(qeri, (ei_a, ei_b), (ea_a, ea_b), **facs) else: vv, vev = _agf2.build_mats_uagf2_outcore(qeri, (ei_a, ei_b), (ea_a, ea_b), **facs) e, c = _agf2.cholesky_build(vv, vev) se_a = aux.SelfEnergy(e, c, chempot=gf_occ[0].chempot) se_a.remove_uncoupled(tol=tol) if not (agf2.frozen is None or agf2.frozen == 0): mask = get_frozen_mask(agf2) coupling = np.zeros((nmo[0], se_a.naux)) coupling[mask[0]] = se_a.coupling se_a = aux.SelfEnergy(se_a.energy, coupling, chempot=se_a.chempot) cput0 = log.timer('se part (alpha)', *cput0) mem_incore = (nmo[1]*nob*(nob*nvb+noa*nva)) * 8/1e6 mem_now = lib.current_memory()[0] if (mem_incore+mem_now < agf2.max_memory) or agf2.incore_complete: qeri = _make_qmo_eris_incore(agf2, eri, (ci_a, ci_a, ca_a), (ci_b, ci_b, ca_b), spin=1) else: qeri = _make_qmo_eris_outcore(agf2, eri, (ci_a, ci_a, ca_a), (ci_b, ci_b, ca_b), spin=1) if isinstance(qeri[0], np.ndarray): vv, vev = _agf2.build_mats_uagf2_incore(qeri, (ei_b, ei_a), (ea_b, ea_a), **facs) else: vv, vev = _agf2.build_mats_uagf2_outcore(qeri, (ei_b, ei_a), (ea_b, ea_a), **facs) e, c = _agf2.cholesky_build(vv, vev) se_b = aux.SelfEnergy(e, c, chempot=gf_occ[1].chempot) se_b.remove_uncoupled(tol=tol) if not (agf2.frozen is None or agf2.frozen == 0): mask = get_frozen_mask(agf2) coupling = np.zeros((nmo[1], se_b.naux)) coupling[mask[1]] = se_b.coupling se_b = aux.SelfEnergy(se_b.energy, coupling, chempot=se_b.chempot) cput0 = log.timer('se part (beta)', *cput0) return (se_a, se_b)
[docs] def get_fock(agf2, eri, gf=None, rdm1=None): ''' Computes the physical space Fock matrix in MO basis. If :attr:`rdm1` is not supplied, it is built from :attr:`gf`, which defaults to the mean-field Green's function. Args: eri : _ChemistsERIs Electronic repulsion integrals Kwargs: gf : GreensFunction Auxiliaries of the Green's function rdm1 : 2D array Reduced density matrix Returns: ndarray of physical space Fock matrix ''' if rdm1 is None: rdm1 = agf2.make_rdm1(gf) vj_aa, vk_aa = agf2.get_jk(eri.eri_aa, rdm1=rdm1[0]) vj_bb, vk_bb = agf2.get_jk(eri.eri_bb, rdm1=rdm1[1]) vj_ab = agf2.get_jk(eri.eri_ab, rdm1=rdm1[1], with_k=False)[0] vj_ba = agf2.get_jk(eri.eri_ba, rdm1=rdm1[0], with_k=False)[0] fock_a = eri.h1e[0] + vj_aa + vj_ab - vk_aa fock_b = eri.h1e[1] + vj_bb + vj_ba - vk_bb fock = (fock_a, fock_b) return fock
[docs] def fock_loop(agf2, eri, gf, se): ''' Self-consistent loop for the density matrix via the HF self- consistent field. Args: eri : _ChemistsERIs Electronic repulsion integrals gf : tuple of GreensFunction Auxiliaries of the Green's function for each spin se : tuple of SelfEnergy Auxiliaries of the self-energy for each spin Returns: :class:`SelfEnergy`, :class:`GreensFunction` and a boolean indicating whether convergence was successful. ''' assert type(gf[0]) is aux.GreensFunction assert type(gf[1]) is aux.GreensFunction assert type(se[0]) is aux.SelfEnergy assert type(se[1]) is aux.SelfEnergy cput0 = cput1 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(agf2.stdout, agf2.verbose) diis = lib.diis.DIIS(agf2) diis.space = agf2.fock_diis_space diis.min_space = agf2.fock_diis_min_space focka, fockb = agf2.get_fock(eri, gf) sea, seb = se gfa, gfb = gf nalph, nbeta = agf2.nocc nmoa, nmob = eri.nmo nauxa, nauxb = sea.naux, seb.naux nqmoa, nqmob = nauxa+nmoa, nauxb+nmob bufa, bufb = np.zeros((nqmoa, nqmoa)), np.zeros((nqmob, nqmob)) rdm1a_prev = 0 rdm1b_prev = 0 converged = False opts = {'tol': agf2.conv_tol_nelec, 'maxiter': agf2.max_cycle_inner} for niter1 in range(1, agf2.max_cycle_outer+1): sea, opt = minimize_chempot(sea, focka, nalph, x0=sea.chempot, occupancy=1, **opts) seb, opt = minimize_chempot(seb, fockb, nbeta, x0=seb.chempot, occupancy=1, **opts) for niter2 in range(1, agf2.max_cycle_inner+1): wa, va = sea.eig(focka, chempot=0.0, out=bufa) wb, vb = seb.eig(fockb, chempot=0.0, out=bufb) sea.chempot, nerra = \ binsearch_chempot((wa, va), nmoa, nalph, occupancy=1) seb.chempot, nerrb = \ binsearch_chempot((wb, vb), nmob, nbeta, occupancy=1) nerr = max(nerra, nerrb) wa, va = sea.eig(focka, out=bufa) wb, vb = seb.eig(fockb, out=bufb) gfa = aux.GreensFunction(wa, va[:nmoa], chempot=sea.chempot) gfb = aux.GreensFunction(wb, vb[:nmob], chempot=seb.chempot) gf = (gfa, gfb) focka, fockb = agf2.get_fock(eri, gf) rdm1a, rdm1b = agf2.make_rdm1(gf) focka, fockb = diis.update(np.array((focka, fockb)), xerr=None) if niter2 > 1: derra = np.max(np.absolute(rdm1a - rdm1a_prev)) derrb = np.max(np.absolute(rdm1b - rdm1b_prev)) derr = max(derra, derrb) if derr < agf2.conv_tol_rdm1: break rdm1a_prev = rdm1a.copy() rdm1b_prev = rdm1b.copy() log.debug1('fock loop %d cycles = %d dN = %.3g |ddm| = %.3g', niter1, niter2, nerr, derr) cput1 = log.timer_debug1('fock loop %d'%niter1, *cput1) if derr < agf2.conv_tol_rdm1 and abs(nerr) < agf2.conv_tol_nelec: converged = True break se = (sea, seb) log.info('fock converged = %s' % converged) log.info(' alpha: chempot = %.9g dN = %.3g |ddm| = %.3g', sea.chempot, nerra, derra) log.info(' beta: chempot = %.9g dN = %.3g |ddm| = %.3g', seb.chempot, nerrb, derrb) log.timer('fock loop', *cput0) return gf, se, converged
[docs] def energy_1body(agf2, eri, gf): ''' Calculates the one-body energy according to the UHF form. Args: eri : _ChemistsERIs Electronic repulsion integrals gf : tuple of GreensFunction Auxiliaries of the Green's function for each spin Returns: One-body energy ''' assert type(gf[0]) is aux.GreensFunction assert type(gf[1]) is aux.GreensFunction rdm1 = agf2.make_rdm1(gf) fock = agf2.get_fock(eri, gf) e1b_a = 0.5 * np.sum(rdm1[0] * (eri.h1e[0] + fock[0])) e1b_b = 0.5 * np.sum(rdm1[1] * (eri.h1e[1] + fock[1])) e1b = e1b_a + e1b_b e1b += agf2.energy_nuc() return e1b
[docs] def energy_2body(agf2, gf, se): ''' Calculates the two-body energy using analytically integrated Galitskii-Migdal formula. The formula is symmetric and only one side needs to be calculated. Args: gf : tuple of GreensFunction Auxiliaries of the Green's function for each spin se : tuple of SelfEnergy Auxiliaries of the self-energy for each spin Returns: Two-body energy ''' e2b_a = ragf2.energy_2body(agf2, gf[0], se[0]) e2b_b = ragf2.energy_2body(agf2, gf[1], se[1]) e2b = (e2b_a + e2b_b) * 0.5 return e2b
[docs] def energy_mp2(agf2, gf, se): ''' Calculates the two-body energy using analytically integrated Galitskii-Migdal formula for an MP2 self-energy. Per the definition of one- and two-body partitioning in the Dyson equation, this result is half of :func:`energy_2body`. Args: gf : tuple of GreensFunction Auxiliaries of the Green's function for each spin se : tuple of SelfEnergy Auxiliaries of the self-energy for each spin Returns: MP2 energy ''' emp2_a = ragf2.energy_mp2(agf2, gf[0], se[0]) emp2_b = ragf2.energy_mp2(agf2, gf[1], se[1]) emp2 = (emp2_a + emp2_b) * 0.5 return emp2
[docs] class UAGF2(ragf2.RAGF2): ''' Unrestricted AGF2 with canonical HF reference Attributes: verbose : int Print level. Default value equals to :class:`Mole.verbose` max_memory : float or int Allowed memory in MB. Default value equals to :class:`Mole.max_memory` incore_complete : bool Avoid all I/O. Default is False. conv_tol : float Convergence threshold for AGF2 energy. Default value is 1e-7 conv_tol_rdm1 : float Convergence threshold for first-order reduced density matrix. Default value is 1e-8. conv_tol_nelec : float Convergence threshold for the number of electrons. Default value is 1e-6. max_cycle : int Maximum number of AGF2 iterations. Default value is 50. max_cycle_outer : int Maximum number of outer Fock loop iterations. Default value is 20. max_cycle_inner : int Maximum number of inner Fock loop iterations. Default value is 50. weight_tol : float Threshold in spectral weight of auxiliaries to be considered zero. Default 1e-11. diis : bool or lib.diis.DIIS Whether to use DIIS, can also be a lib.diis.DIIS object. Default value is True. diis_space : int DIIS space size. Default value is 8. diis_min_space : int Minimum space of DIIS. Default value is 1. fock_diis_space : int DIIS space size for Fock loop iterations. Default value is 6. fock_diis_min_space : Minimum space of DIIS. Default value is 1. os_factor : float Opposite-spin factor for spin-component-scaled (SCS) calculations. Default 1.0 ss_factor : float Same-spin factor for spin-component-scaled (SCS) calculations. Default 1.0 damping : float Damping factor for the self-energy. Default value is 0.0 Saved results e_corr : float AGF2 correlation energy e_tot : float Total energy (HF + correlation) e_1b : float One-body part of :attr:`e_tot` e_2b : float Two-body part of :attr:`e_tot` e_init : float Initial correlation energy (truncated MP2) converged : bool Whether convergence was successful se : tuple of SelfEnergy Auxiliaries of the self-energy for each spin gf : tuple of GreensFunction Auxiliaries of the Green's function for each spin ''' energy_1body = energy_1body energy_2body = energy_2body fock_loop = fock_loop build_se_part = build_se_part
[docs] def ao2mo(self, mo_coeff=None): ''' Get the electronic repulsion integrals in MO basis. ''' nmo = max(self.nmo) mem_incore = ((nmo*(nmo+1)//2)**2) * 8/1e6 mem_now = lib.current_memory()[0] if (self._scf._eri is not None and (mem_incore+mem_now < self.max_memory or self.incore_complete)): eri = _make_mo_eris_incore(self, mo_coeff) else: logger.warn(self, 'MO eris are outcore - this may be very ' 'slow for agf2. increasing max_memory or ' 'using density fitting is recommended.') eri = _make_mo_eris_outcore(self, mo_coeff) return eri
[docs] def make_rdm1(self, gf=None): ''' Compute the one-body reduced density matrix in MO basis. Kwargs: gf : tuple of GreensFunction Auxiliaries of the Green's functions for each spin Returns: tuple of ndarray of density matrices ''' if gf is None: gf = self.gf if gf is None: gf = self.init_gf() rdm1_a = gf[0].make_rdm1(occupancy=1) rdm1_b = gf[1].make_rdm1(occupancy=1) return (rdm1_a, rdm1_b)
[docs] def get_fock(self, eri=None, gf=None, rdm1=None): ''' Computes the physical space Fock matrix in MO basis. ''' if eri is None: eri = self.ao2mo() if gf is None: gf = self.gf return get_fock(self, eri, gf=gf, rdm1=rdm1)
[docs] def energy_mp2(self, mo_energy=None, se=None): if mo_energy is None: mo_energy = self.mo_energy if se is None: se = self.build_se(gf=self.gf) self.e_init = energy_mp2(self, mo_energy, se) return self.e_init
[docs] def init_gf(self, frozen=False): ''' Builds the Hartree-Fock Green's function. Returns: tuple of :class:`GreensFunction`, tuple of :class:`SelfEnergy` ''' nmoa, nmob = self.nmo nocca, noccb = self.nocc energy = self.mo_energy coupling = (np.eye(nmoa), np.eye(nmob)) focka = np.diag(energy[0]) fockb = np.diag(energy[1]) cpt_a = binsearch_chempot(focka, nmoa, nocca, occupancy=1)[0] cpt_b = binsearch_chempot(fockb, nmob, noccb, occupancy=1)[1] if frozen: mask = get_frozen_mask(self) energy = (energy[0][mask[0]], energy[1][mask[1]]) coupling = (coupling[0][:,mask[0]], coupling[1][:,mask[1]]) gf_a = aux.GreensFunction(energy[0], coupling[0], chempot=cpt_a) gf_b = aux.GreensFunction(energy[1], coupling[1], chempot=cpt_b) gf = (gf_a, gf_b) return gf
[docs] def build_gf(self, eri=None, gf=None, se=None): ''' Builds the auxiliaries of the Green's functions by solving the Dyson equation for each spin. Kwargs: eri : _ChemistsERIs Electronic repulsion integrals gf : tuple of GreensFunction Auxiliaries of the Green's function for each spin se : tuple of SelfEnergy Auxiliaries of the self-energy for each spin Returns: tuple of :class:`GreensFunction` ''' if eri is None: eri = self.ao2mo() if gf is None: gf = self.gf if gf is None: gf = self.init_gf() if se is None: se = self.build_se(eri, gf) focka, fockb = self.get_fock(eri, gf) gf_a = se[0].get_greens_function(focka) gf_b = se[1].get_greens_function(fockb) return (gf_a, gf_b)
[docs] def build_se(self, eri=None, gf=None, os_factor=None, ss_factor=None, se_prev=None): ''' Builds the auxiliaries of the self-energy. Args: eri : _ChemistsERIs Electronic repulsion integrals gf : tuple of GreensFunction Auxiliaries of the Green's function Kwargs: os_factor : float Opposite-spin factor for spin-component-scaled (SCS) calculations. Default 1.0 ss_factor : float Same-spin factor for spin-component-scaled (SCS) calculations. Default 1.0 se_prev : SelfEnergy Previous self-energy for damping. Default value is None Returns tuple of :class:`SelfEnergy` ''' if eri is None: eri = self.ao2mo() if gf is None: gf = self.gf if gf is None: gf = self.init_gf() if os_factor is None: os_factor = self.os_factor if ss_factor is None: ss_factor = self.ss_factor facs = {'os_factor': os_factor, 'ss_factor': ss_factor} gf_occ = (gf[0].get_occupied(), gf[1].get_occupied()) gf_vir = (gf[0].get_virtual(), gf[1].get_virtual()) se_occ = self.build_se_part(eri, gf_occ, gf_vir, **facs) se_vir = self.build_se_part(eri, gf_vir, gf_occ, **facs) se_a = aux.combine(se_occ[0], se_vir[0]) se_b = aux.combine(se_occ[1], se_vir[1]) if se_prev is not None and self.damping != 0.0: se_a_prev, se_b_prev = se_prev se_a.coupling *= np.sqrt(1.0-self.damping) se_b.coupling *= np.sqrt(1.0-self.damping) se_a_prev.coupling *= np.sqrt(self.damping) se_b_prev.coupling *= np.sqrt(self.damping) se_a = aux.combine(se_a, se_a_prev) se_b = aux.combine(se_b, se_b_prev) se_a = se_a.compress(n=(None,0)) se_b = se_b.compress(n=(None,0)) return (se_a, se_b)
[docs] def run_diis(self, se, diis=None): ''' Runs the direct inversion of the iterative subspace for the self-energy. Args: se : SelfEnergy Auxiliaries of the self-energy diis : lib.diis.DIIS DIIS object Returns: tuple of :class:`SelfEnergy` ''' if diis is None: return se se_occ_a, se_occ_b = (se[0].get_occupied(), se[1].get_occupied()) se_vir_a, se_vir_b = (se[0].get_virtual(), se[1].get_virtual()) vv_occ_a = np.dot(se_occ_a.coupling, se_occ_a.coupling.T) vv_occ_b = np.dot(se_occ_b.coupling, se_occ_b.coupling.T) vv_vir_a = np.dot(se_vir_a.coupling, se_vir_a.coupling.T) vv_vir_b = np.dot(se_vir_b.coupling, se_vir_b.coupling.T) vev_occ_a = np.dot(se_occ_a.coupling * se_occ_a.energy[None], se_occ_a.coupling.T) vev_occ_b = np.dot(se_occ_b.coupling * se_occ_b.energy[None], se_occ_b.coupling.T) vev_vir_a = np.dot(se_vir_a.coupling * se_vir_a.energy[None], se_vir_a.coupling.T) vev_vir_b = np.dot(se_vir_b.coupling * se_vir_b.energy[None], se_vir_b.coupling.T) dat = np.array([vv_occ_a, vv_vir_a, vev_occ_a, vev_vir_a, vv_occ_b, vv_vir_b, vev_occ_b, vev_vir_b]) dat = diis.update(dat) vv_occ_a, vv_vir_a, vev_occ_a, vev_vir_a, \ vv_occ_b, vv_vir_b, vev_occ_b, vev_vir_b = dat se_occ_a = aux.SelfEnergy(*_agf2.cholesky_build(vv_occ_a, vev_occ_a), chempot=se[0].chempot) se_vir_a = aux.SelfEnergy(*_agf2.cholesky_build(vv_vir_a, vev_vir_a), chempot=se[0].chempot) se_occ_b = aux.SelfEnergy(*_agf2.cholesky_build(vv_occ_b, vev_occ_b), chempot=se[1].chempot) se_vir_b = aux.SelfEnergy(*_agf2.cholesky_build(vv_vir_b, vev_vir_b), chempot=se[1].chempot) se = (aux.combine(se_occ_a, se_vir_a), aux.combine(se_occ_b, se_vir_b)) return se
[docs] def density_fit(self, auxbasis=None, with_df=None): from pyscf.agf2 import dfuagf2 myagf2 = dfuagf2.DFUAGF2(self._scf) myagf2.__dict__.update(self.__dict__) if with_df is not None: myagf2.with_df = with_df if auxbasis is not None and myagf2.with_df.auxbasis != auxbasis: myagf2.with_df = myagf2.with_df.copy() myagf2.with_df.auxbasis = auxbasis return myagf2
[docs] def get_ip(self, gf, nroots=5): gf_occ = (gf[0].get_occupied(), gf[1].get_occupied()) spin = np.array([0,]*gf_occ[0].naux + [1,]*gf_occ[1].naux) e_ip = np.concatenate([gf_occ[0].energy, gf_occ[1].energy], axis=0) v_ip = np.concatenate([gf_occ[0].coupling, gf_occ[1].coupling], axis=1) mask = np.argsort(e_ip) spin = list(spin[mask][-nroots:])[::-1] e_ip = list(-e_ip[mask][-nroots:])[::-1] v_ip = list(v_ip[:,mask][:,-nroots:].T)[::-1] return e_ip, v_ip, spin
[docs] def ipagf2(self, nroots=5): e_ip, v_ip, spin = self.get_ip(self.gf, nroots=nroots) for n, en, vn, sn in zip(range(nroots), e_ip, v_ip, spin): qpwt = np.linalg.norm(vn)**2 tag = ['alpha', 'beta'][sn] logger.note(self, 'IP energy level %d E = %.16g QP weight = %0.6g (%s)', n, en, qpwt, tag) if nroots == 1: return e_ip[0], v_ip[0] else: return e_ip, v_ip
[docs] def get_ea(self, gf, nroots=5): gf_vir = (gf[0].get_virtual(), gf[1].get_virtual()) spin = np.array([0,]*gf_vir[0].naux + [1,]*gf_vir[1].naux) e_ea = np.concatenate([gf_vir[0].energy, gf_vir[1].energy], axis=0) v_ea = np.concatenate([gf_vir[0].coupling, gf_vir[1].coupling], axis=1) mask = np.argsort(e_ea) spin = list(spin[mask][:nroots]) e_ea = list(e_ea[mask][:nroots]) v_ea = list(v_ea[:,mask][:,:nroots].T) return e_ea, v_ea, spin
[docs] def eaagf2(self, nroots=5): e_ea, v_ea, spin = self.get_ea(self.gf, nroots=nroots) for n, en, vn, sn in zip(range(nroots), e_ea, v_ea, spin): qpwt = np.linalg.norm(vn)**2 tag = ['alpha', 'beta'][sn] logger.note(self, 'EA energy level %d E = %.16g QP weight = %0.6g (%s)', n, en, qpwt, tag) if nroots == 1: return e_ea[0], v_ea[0] else: return e_ea, v_ea
@property def nocc(self): if self._nocc is None: self._nocc = (np.sum(self.mo_occ[0] > 0), np.sum(self.mo_occ[1] > 0)) return self._nocc @nocc.setter def nocc(self, val): self._nocc = val @property def nmo(self): if self._nmo is None: self._nmo = (self.mo_occ[0].size, self.mo_occ[1].size) return self._nmo @nmo.setter def nmo(self, val): self._nmo = val @property def qmo_energy(self): return (self.gf[0].energy, self.gf[1].energy) @property def qmo_coeff(self): ''' Gives the couplings in AO basis ''' return (np.dot(self.mo_coeff[0], self.gf[0].coupling), np.dot(self.mo_coeff[1], self.gf[1].coupling)) @property def qmo_occ(self): coeff_a = self.gf[0].get_occupied().coupling coeff_b = self.gf[1].get_occupied().coupling occ_a = np.linalg.norm(coeff_a, axis=0) ** 2 occ_b = np.linalg.norm(coeff_b, axis=0) ** 2 vir_a = np.zeros_like(self.gf[0].get_virtual().energy) vir_b = np.zeros_like(self.gf[1].get_virtual().energy) qmo_occ_a = np.concatenate([occ_a, vir_a]) qmo_occ_b = np.concatenate([occ_b, vir_b]) return qmo_occ_a, qmo_occ_b
[docs] def get_frozen_mask(agf2): with lib.temporary_env(agf2, _nocc=None, _nmo=None): return _get_frozen_mask(agf2)
class _ChemistsERIs: ''' (pq|rs) MO integrals stored in s4 symmetry, we only need QMO integrals in low-symmetry tensors and s4 is highest supported by _vhf ''' def __init__(self, mol=None): self.mol = mol self.mo_coeff = None self.nocc = None self.nmo = None self.fock = None self.h1e = None self.eri = None self.e_hf = None def _common_init_(self, agf2, mo_coeff=None): if mo_coeff is None: mo_coeff = agf2.mo_coeff self.mo_coeff = mo_coeff dm = agf2._scf.make_rdm1(agf2.mo_coeff, agf2.mo_occ) h1e_ao = agf2._scf.get_hcore() vhf = agf2._scf.get_veff(agf2.mol, dm) fock_ao = agf2._scf.get_fock(vhf=vhf, dm=dm) self.h1e = (np.dot(np.dot(mo_coeff[0].conj().T, h1e_ao), mo_coeff[0]), np.dot(np.dot(mo_coeff[1].conj().T, h1e_ao), mo_coeff[1])) self.fock = (np.dot(np.dot(mo_coeff[0].conj().T, fock_ao[0]), mo_coeff[0]), np.dot(np.dot(mo_coeff[1].conj().T, fock_ao[1]), mo_coeff[1])) self.h1e = (mpi_helper.bcast(self.h1e[0]), mpi_helper.bcast(self.h1e[1])) self.fock = (mpi_helper.bcast(self.fock[0]), mpi_helper.bcast(self.fock[1])) self.e_hf = mpi_helper.bcast(agf2._scf.e_tot) self.nmo = agf2.nmo nocca, noccb = self.nocc = agf2.nocc self.mol = agf2.mol mo_e = (self.fock[0].diagonal(), self.fock[1].diagonal()) gap_a = abs(mo_e[0][:nocca,None] - mo_e[0][None,nocca:]).min() gap_b = abs(mo_e[1][:noccb,None] - mo_e[1][None,noccb:]).min() gap = min(gap_a, gap_b) if gap < 1e-5: logger.warn(agf2, 'HOMO-LUMO gap %s may be too small for AGF2', gap) return self def _make_mo_eris_incore(agf2, mo_coeff=None): ''' Returns _ChemistsERIs ''' cput0 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(agf2.stdout, agf2.verbose) eris = _ChemistsERIs() eris._common_init_(agf2, mo_coeff) moa, mob = eris.mo_coeff nmoa, nmob = eris.nmo eri_aa = ao2mo.incore.full(agf2._scf._eri, moa, verbose=log) eri_bb = ao2mo.incore.full(agf2._scf._eri, mob, verbose=log) eri_aa = ao2mo.addons.restore('s4', eri_aa, nmoa) eri_bb = ao2mo.addons.restore('s4', eri_bb, nmob) eri_ab = ao2mo.incore.general(agf2._scf._eri, (moa,moa,mob,mob), verbose=log) assert eri_ab.shape == (nmoa*(nmob+1)//2, nmob*(nmob+1)//2) eri_ba = np.transpose(eri_ab) eris.eri_aa = eri_aa eris.eri_ab = eri_ab eris.eri_ba = eri_ba eris.eri_bb = eri_bb eris.eri = ((eri_aa, eri_ab), (eri_ba, eri_bb)) log.timer('MO integral transformation', *cput0) return eris def _make_mo_eris_outcore(agf2, mo_coeff=None): ''' Returns _ChemistsERIs ''' log = logger.Logger(agf2.stdout, agf2.verbose) eris = _ChemistsERIs() eris._common_init_(agf2, mo_coeff) mol = agf2.mol moa = np.asarray(eris.mo_coeff[0], order='F') mob = np.asarray(eris.mo_coeff[1], order='F') nmoa, nmob = eris.nmo eris.feri = lib.H5TmpFile() ao2mo.outcore.full(mol, moa, eris.feri, dataname='mo/aa') ao2mo.outcore.full(mol, mob, eris.feri, dataname='mo/bb') ao2mo.outcore.general(mol, (moa,moa,mob,mob), eris.feri, dataname='mo/ab', verbose=log) ao2mo.outcore.general(mol, (mob,mob,moa,moa), eris.feri, dataname='mo/ba', verbose=log) eris.eri_aa = eris.feri['mo/aa'] eris.eri_ab = eris.feri['mo/ab'] eris.eri_ba = eris.feri['mo/ba'] eris.eri_bb = eris.feri['mo/bb'] eris.eri = ((eris.eri_aa, eris.eri_ab), (eris.eri_ba, eris.eri_bb)) return eris def _make_qmo_eris_incore(agf2, eri, coeffs_a, coeffs_b, spin=None): ''' Returns nested tuple of ndarray spin = None: ((aaaa, aabb), (bbaa, bbbb)) spin = 0: (aaaa, aabb) spin = 1: (bbbb, bbaa) ''' cput0 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(agf2.stdout, agf2.verbose) nmo = eri.nmo nmoa, nmob = nmo cxa = np.eye(nmoa) cxb = np.eye(nmob) if not (agf2.frozen is None or agf2.frozen == 0): mask = get_frozen_mask(agf2) cxa = cxa[:,mask[0]] cxb = cxb[:,mask[1]] # npaira, npairb = nmoa*(nmoa+1)//2, nmob*(nmob+1)//2 cia, cja, caa = coeffs_a cib, cjb, cab = coeffs_b nia, nja, naa = [x.shape[1] for x in coeffs_a] nib, njb, nab = [x.shape[1] for x in coeffs_b] if spin is None or spin == 0: c_aa = (cxa, cia, cja, caa) c_ab = (cxa, cia, cjb, cab) qeri_aa = ao2mo.incore.general(eri.eri_aa, c_aa, compact=False, verbose=log) qeri_ab = ao2mo.incore.general(eri.eri_ab, c_ab, compact=False, verbose=log) qeri_aa = qeri_aa.reshape(cxa.shape[1], nia, nja, naa) qeri_ab = qeri_ab.reshape(cxa.shape[1], nia, njb, nab) if spin is None or spin == 1: c_bb = (cxb, cib, cjb, cab) c_ba = (cxb, cib, cja, caa) qeri_bb = ao2mo.incore.general(eri.eri_bb, c_bb, compact=False, verbose=log) qeri_ba = ao2mo.incore.general(eri.eri_ba, c_ba, compact=False, verbose=log) qeri_bb = qeri_bb.reshape(cxb.shape[1], nib, njb, nab) qeri_ba = qeri_ba.reshape(cxb.shape[1], nib, nja, naa) if spin is None: qeri = ((qeri_aa, qeri_ab), (qeri_ba, qeri_bb)) elif spin == 0: qeri = (qeri_aa, qeri_ab) elif spin == 1: qeri = (qeri_bb, qeri_ba) log.timer('QMO integral transformation', *cput0) return qeri def _make_qmo_eris_outcore(agf2, eri, coeffs_a, coeffs_b, spin=None): ''' Returns nested tuple of H5 dataset spin = None: ((aaaa, aabb), (bbaa, bbbb)) spin = 0: (aaaa, aabb) spin = 1: (bbbb, bbaa) ''' cput0 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(agf2.stdout, agf2.verbose) nmo = eri.nmo nmoa, nmob = nmo mask = get_frozen_mask(agf2) frozena = np.sum(~mask[0]) frozenb = np.sum(~mask[1]) # npaira, npairb = nmoa*(nmoa+1)//2, nmob*(nmob+1)//2 cia, cja, caa = coeffs_a cib, cjb, cab = coeffs_b nia, nja, naa = [x.shape[1] for x in coeffs_a] nib, njb, nab = [x.shape[1] for x in coeffs_b] # possible to have incore MO, outcore QMO if getattr(eri, 'feri', None) is None: eri.feri = lib.H5TmpFile() else: for key in ['aa', 'ab', 'ba', 'bb']: if 'qmo/%s'%key in eri.feri: del eri.feri['qmo/%s'%key] if spin is None or spin == 0: eri.feri.create_dataset('qmo/aa', (nmoa-frozena, nia, nja, naa), 'f8') eri.feri.create_dataset('qmo/ab', (nmoa-frozena, nia, njb, nab), 'f8') blksize = _agf2.get_blksize(agf2.max_memory, (nmoa**3, nmoa*nja*naa), (nmoa*nmob**2, nmoa*njb*nab)) blksize = min(nmoa, max(BLKMIN, blksize)) log.debug1('blksize (uagf2._make_qmo_eris_outcore) = %d', blksize) tril2sq = lib.square_mat_in_trilu_indices(nmoa) q1 = 0 for p0, p1 in lib.prange(0, nmoa, blksize): if not np.any(mask[0][p0:p1]): # block is fully frozen continue inds = np.arange(p0, p1)[mask[0][p0:p1]] q0, q1 = q1, q1 + len(inds) idx = list(np.concatenate(tril2sq[inds])) # aa buf = eri.eri_aa[idx] # (blk, nmoa, npaira) buf = buf.reshape((q1-q0)*nmoa, -1) # (blk*nmoa, npaira) jasym_aa, nja_aa, cja_aa, sja_aa = ao2mo.incore._conc_mos(cja, caa) buf = ao2mo._ao2mo.nr_e2(buf, cja_aa, sja_aa, 's2kl', 's1') buf = buf.reshape(q1-q0, nmoa, nja, naa) buf = lib.einsum('xpja,pi->xija', buf, cia) eri.feri['qmo/aa'][q0:q1] = np.asarray(buf, order='C') # ab buf = eri.eri_ab[idx] # (blk, nmoa, npairb) buf = buf.reshape((q1-q0)*nmob, -1) # (blk*nmoa, npairb) jasym_ab, nja_ab, cja_ab, sja_ab = ao2mo.incore._conc_mos(cjb, cab) buf = ao2mo._ao2mo.nr_e2(buf, cja_ab, sja_ab, 's2kl', 's1') buf = buf.reshape(q1-q0, nmoa, njb, nab) buf = lib.einsum('xpja,pi->xija', buf, cia) eri.feri['qmo/ab'][q0:q1] = np.asarray(buf, order='C') if spin is None or spin == 1: eri.feri.create_dataset('qmo/ba', (nmob-frozenb, nib, nja, naa), 'f8') eri.feri.create_dataset('qmo/bb', (nmob-frozenb, nib, njb, nab), 'f8') max_memory = agf2.max_memory - lib.current_memory()[0] blksize = int((max_memory/8e-6) / max(nmob**3+nmob*njb*nab, nmob*nmoa**2*nja*naa)) blksize = min(nmob, max(BLKMIN, blksize)) log.debug1('blksize (uagf2._make_qmo_eris_outcore) = %d', blksize) tril2sq = lib.square_mat_in_trilu_indices(nmob) q1 = 0 for p0, p1 in lib.prange(0, nmob, blksize): if not np.any(mask[1][p0:p1]): # block is fully frozen continue inds = np.arange(p0, p1)[mask[1][p0:p1]] q0, q1 = q1, q1 + len(inds) idx = list(np.concatenate(tril2sq[inds])) # ba buf = eri.eri_ba[idx] # (blk, nmob, npaira) buf = buf.reshape((q1-q0)*nmob, -1) # (blk*nmob, npaira) jasym_ba, nja_ba, cja_ba, sja_ba = ao2mo.incore._conc_mos(cja, caa) buf = ao2mo._ao2mo.nr_e2(buf, cja_ba, sja_ba, 's2kl', 's1') buf = buf.reshape(q1-q0, nmob, nja, naa) buf = lib.einsum('xpja,pi->xija', buf, cib) eri.feri['qmo/ba'][q0:q1] = np.asarray(buf, order='C') # bb buf = eri.eri_bb[idx] # (blk, nmob, npairb) buf = buf.reshape((q1-q0)*nmob, -1) # (blk*nmob, npairb) jasym_bb, nja_bb, cja_bb, sja_bb = ao2mo.incore._conc_mos(cjb, cab) buf = ao2mo._ao2mo.nr_e2(buf, cja_bb, sja_bb, 's2kl', 's1') buf = buf.reshape(q1-q0, nmob, njb, nab) buf = lib.einsum('xpja,pi->xija', buf, cib) eri.feri['qmo/bb'][q0:q1] = np.asarray(buf, order='C') if spin is None: qeri = ((eri.feri['qmo/aa'], eri.feri['qmo/ab']), (eri.feri['qmo/ba'], eri.feri['qmo/bb'])) elif spin == 0: qeri = (eri.feri['qmo/aa'], eri.feri['qmo/ab']) elif spin == 1: qeri = (eri.feri['qmo/bb'], eri.feri['qmo/ba']) log.timer('QMO integral transformation', *cput0) return qeri if __name__ == '__main__': from pyscf import gto, scf, mp mol = gto.M(atom='O 0 0 0; H 0 0 1; H 0 1 0', basis='cc-pvdz', charge=-1, spin=1, verbose=3) uhf = scf.UHF(mol) uhf.conv_tol = 1e-11 uhf.run() uagf2 = UAGF2(uhf, frozen=0) uagf2.run() uagf2.ipagf2(nroots=5) uagf2.eaagf2(nroots=5) uagf2 = uagf2.density_fit() uagf2.run()