# 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
'''
import numpy as np
import copy
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 mpi_helper, _agf2
from pyscf.agf2 import aux_space as aux
from pyscf.agf2 import chkfile as chkutil
from pyscf.agf2.chempot import binsearch_chempot, minimize_chempot
from pyscf.mp.mp2 import get_frozen_mask as _get_frozen_mask
BLKMIN = getattr(__config__, 'agf2_blkmin', 1)
[docs]
def kernel(agf2, eri=None, gf=None, se=None, verbose=None, dump_chk=True):
log = logger.new_logger(agf2, verbose)
cput1 = cput0 = (logger.process_clock(), logger.perf_counter())
name = agf2.__class__.__name__
if eri is None: eri = agf2.ao2mo()
if gf is None: gf = agf2.gf
if se is None: se = agf2.se
if verbose is None: verbose = agf2.verbose
if gf is None:
gf = agf2.init_gf()
gf_froz = agf2.init_gf(frozen=True)
else:
gf_froz = gf
if se is None:
se = agf2.build_se(eri, gf_froz)
if dump_chk:
agf2.dump_chk(gf=gf, se=se)
if isinstance(agf2.diis, lib.diis.DIIS):
diis = agf2.diis
elif agf2.diis:
diis = lib.diis.DIIS(agf2)
diis.space = agf2.diis_space
diis.min_space = agf2.diis_min_space
else:
diis = None
e_init = agf2.energy_mp2(agf2.mo_energy, se)
log.info('E(init) = %.16g E_corr(init) = %.16g', e_init+eri.e_hf, e_init)
e_1b = eri.e_hf
e_2b = e_init
e_prev = 0.0
se_prev = None
converged = False
for niter in range(1, agf2.max_cycle+1):
if agf2.damping != 0.0:
se_prev = copy.deepcopy(se)
# one-body terms
gf, se, fock_conv = agf2.fock_loop(eri, gf, se)
e_1b = agf2.energy_1body(eri, gf)
# two-body terms
se = agf2.build_se(eri, gf, se_prev=se_prev)
se = agf2.run_diis(se, diis)
e_2b = agf2.energy_2body(gf, se)
if dump_chk:
agf2.dump_chk(gf=gf, se=se)
e_tot = e_1b + e_2b
ip = agf2.get_ip(gf, nroots=1)[0][0]
ea = agf2.get_ea(gf, nroots=1)[0][0]
log.info('cycle = %3d E(%s) = %.15g E_corr(%s) = %.15g dE = %.9g',
niter, name, e_tot, name, e_tot-eri.e_hf, e_tot-e_prev)
log.info('E_1b = %.15g E_2b = %.15g', e_1b, e_2b)
log.info('IP = %.15g EA = %.15g', ip, ea)
cput1 = log.timer('%s iter'%name, *cput1)
if abs(e_tot - e_prev) < agf2.conv_tol:
converged = True
break
e_prev = e_tot
if dump_chk:
agf2.dump_chk(gf=gf, se=se)
log.timer('%s'%name, *cput0)
return converged, e_1b, e_2b, gf, se
[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.
Args:
eri : _ChemistsERIs
Electronic repulsion integrals
gf_occ : GreensFunction
Occupied Green's function
gf_vir : GreensFunction
Virtual 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
Returns:
:class:`SelfEnergy`
'''
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(agf2.stdout, agf2.verbose)
assert type(gf_occ) is aux.GreensFunction
assert type(gf_vir) is aux.GreensFunction
nmo = eri.nmo
tol = agf2.weight_tol
facs = dict(os_factor=os_factor, ss_factor=ss_factor)
ci, ei = gf_occ.coupling, gf_occ.energy
ca, ea = gf_vir.coupling, gf_vir.energy
mem_incore = (gf_occ.nphys*gf_occ.naux**2*gf_vir.naux) * 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, ci, ca))
else:
qeri = _make_qmo_eris_outcore(agf2, eri, (ci, ci, ca))
if isinstance(qeri, np.ndarray):
vv, vev = _agf2.build_mats_ragf2_incore(qeri, ei, ea, **facs)
else:
vv, vev = _agf2.build_mats_ragf2_outcore(qeri, ei, ea, **facs)
e, c = _agf2.cholesky_build(vv, vev)
se = aux.SelfEnergy(e, c, chempot=gf_occ.chempot)
se.remove_uncoupled(tol=tol)
if not (agf2.frozen is None or agf2.frozen == 0):
mask = get_frozen_mask(agf2)
coupling = np.zeros((nmo, se.naux))
coupling[mask] = se.coupling
se = aux.SelfEnergy(se.energy, coupling, chempot=se.chempot)
log.timer('se part', *cput0)
return se
[docs]
def get_jk(agf2, eri, rdm1, with_j=True, with_k=True):
''' Get the J/K matrices.
Args:
eri : ndarray or H5 dataset
Electronic repulsion integrals (NOT as _ChemistsERIs)
rdm1 : 2D array
Reduced density matrix
Kwargs:
with_j : bool
Whether to compute J. Default value is True
with_k : bool
Whether to compute K. Default value is True
Returns:
tuple of ndarrays corresponding to J and K, if either are
not requested then they are set to None.
'''
if isinstance(eri, np.ndarray):
vj, vk = _vhf.incore(eri, rdm1, with_j=with_j, with_k=with_k)
else:
nmo = rdm1.shape[0]
npair = nmo*(nmo+1)//2
vj = vk = None
if with_j:
rdm1_tril = lib.pack_tril(rdm1 + np.tril(rdm1, k=-1))
vj = np.zeros_like(rdm1_tril)
if with_k:
vk = np.zeros_like(rdm1)
blksize = _agf2.get_blksize(agf2.max_memory, (nmo*npair, nmo**3))
blksize = min(1, max(BLKMIN, blksize))
logger.debug1(agf2, 'blksize (ragf2.get_jk) = %d' % blksize)
tril2sq = lib.square_mat_in_trilu_indices(nmo)
for p0, p1 in lib.prange(0, nmo, blksize):
idx = list(np.concatenate(tril2sq[p0:p1]))
eri0 = eri[idx]
# vj built in tril layout with scaled rdm1_tril
if with_j:
vj[idx] = np.dot(eri0, rdm1_tril)
if with_k:
eri0 = lib.unpack_tril(eri0, axis=-1)
eri0 = eri0.reshape(p1-p0, nmo, nmo, nmo)
vk[p0:p1] = lib.einsum('ijkl,jk->il', eri0, rdm1)
if with_j:
vj = lib.unpack_tril(vj)
return vj, vk
[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, vk = agf2.get_jk(eri.eri, rdm1)
fock = eri.h1e + vj - 0.5 * vk
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 : _ChemistERIs
Electronic repulsion integrals
gf : GreensFunction
Auxiliaries of the Green's function
se : SelfEnergy
Auxiliaries of the self-energy
Returns:
:class:`SelfEnergy`, :class:`GreensFunction` and a boolean
indicating whether convergence was successful.
'''
assert type(gf) is aux.GreensFunction
assert type(se) 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
fock = agf2.get_fock(eri, gf)
nelec = eri.nocc * 2
nmo = eri.nmo
naux = se.naux
nqmo = nmo + naux
buf = np.zeros((nqmo, nqmo))
converged = False
opts = dict(tol=agf2.conv_tol_nelec, maxiter=agf2.max_cycle_inner)
rdm1_prev = 0
for niter1 in range(1, agf2.max_cycle_outer+1):
se, opt = minimize_chempot(se, fock, nelec, x0=se.chempot, **opts)
for niter2 in range(1, agf2.max_cycle_inner+1):
w, v = se.eig(fock, chempot=0.0, out=buf)
se.chempot, nerr = binsearch_chempot((w, v), nmo, nelec)
w, v = se.eig(fock, out=buf)
gf = aux.GreensFunction(w, v[:nmo], chempot=se.chempot)
fock = agf2.get_fock(eri, gf)
rdm1 = agf2.make_rdm1(gf)
fock = diis.update(fock, xerr=None)
if niter2 > 1:
derr = np.max(np.absolute(rdm1 - rdm1_prev))
if derr < agf2.conv_tol_rdm1:
break
rdm1_prev = rdm1.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
log.info('fock converged = %s chempot = %.9g dN = %.3g |ddm| = %.3g',
converged, se.chempot, nerr, derr)
log.timer('fock loop', *cput0)
return gf, se, converged
[docs]
def energy_1body(agf2, eri, gf):
''' Calculates the one-body energy according to the RHF form.
Args:
eri : _ChemistsERIs
Electronic repulsion integrals
gf : GreensFunction
Auxiliaries of Green's function
Returns:
One-body energy
'''
assert type(gf) is aux.GreensFunction
rdm1 = agf2.make_rdm1(gf)
fock = agf2.get_fock(eri, gf)
e1b = 0.5 * np.sum(rdm1 * (eri.h1e + fock))
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 : GreensFunction
Auxiliaries of the Green's function
se : SelfEnergy
Auxiliaries of the self-energy
Returns
Two-body energy
'''
assert type(gf) is aux.GreensFunction
assert type(se) is aux.SelfEnergy
gf_occ = gf.get_occupied()
se_vir = se.get_virtual()
e2b = 0.0
for l in mpi_helper.nrange(gf_occ.naux):
vxl = gf_occ.coupling[:,l]
vxk = se_vir.coupling
dlk = gf_occ.energy[l] - se_vir.energy
vv = vxk * vxl[:,None]
e2b += lib.einsum('xk,yk,k->', vv, vv.conj(), 1./dlk)
e2b *= 2
mpi_helper.barrier()
e2b = mpi_helper.allreduce(e2b)
return np.ravel(e2b.real)[0]
[docs]
def energy_mp2(agf2, mo_energy, 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 : GreensFunction
Auxiliaries of the Green's function
se : SelfEnergy
Auxiliaries of the self-energy
Returns
MP2 energy
'''
assert type(se) is aux.SelfEnergy
occ = mo_energy < se.chempot
se_vir = se.get_virtual()
vxk = se_vir.coupling[occ]
dxk = lib.direct_sum('x,k->xk', mo_energy[occ], -se_vir.energy)
emp2 = lib.einsum('xk,xk,xk->', vxk, vxk.conj(), 1./dxk)
return np.ravel(emp2.real)[0]
[docs]
class RAGF2(lib.StreamObject):
''' Restricted 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 : int
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 : SelfEnergy
Auxiliaries of the self-energy
gf : GreensFunction
Auxiliaries of the Green's function
'''
async_io = getattr(__config__, 'agf2_async_io', True)
incore_complete = getattr(__config__, 'agf2_incore_complete', False)
_keys = set((
'async_io', 'mol', 'incore_complete',
'conv_tol', 'conv_tol_rdm1', 'conv_tol_nelec', 'max_cycle',
'max_cycle_outer', 'max_cycle_inner', 'weight_tol', 'fock_diis_space',
'fock_diis_min_space', 'diis', 'diis_space', 'diis_min_space',
'os_factor', 'ss_factor', 'damping',
'mo_energy', 'mo_coeff', 'mo_occ', 'se', 'gf', 'e_1b', 'e_2b', 'e_init',
'frozen', 'converged', 'chkfile',
))
def __init__(self, mf, frozen=None, mo_energy=None, mo_coeff=None, mo_occ=None):
if mo_energy is None: mo_energy = mpi_helper.bcast(mf.mo_energy)
if mo_coeff is None: mo_coeff = mpi_helper.bcast(mf.mo_coeff)
if mo_occ is None: mo_occ = mpi_helper.bcast(mf.mo_occ)
self.mol = mf.mol
self._scf = mf
self.verbose = self.mol.verbose
self.stdout = self.mol.stdout
self.max_memory = mf.max_memory
self.incore_complete = self.incore_complete or self.mol.incore_anyway
self.conv_tol = getattr(__config__, 'agf2_conv_tol', 1e-7)
self.conv_tol_rdm1 = getattr(__config__, 'agf2_conv_tol_rdm1', 1e-8)
self.conv_tol_nelec = getattr(__config__, 'agf2_conv_tol_nelec', 1e-6)
self.max_cycle = getattr(__config__, 'agf2_max_cycle', 50)
self.max_cycle_outer = getattr(__config__, 'agf2_max_cycle_outer', 20)
self.max_cycle_inner = getattr(__config__, 'agf2_max_cycle_inner', 50)
self.weight_tol = getattr(__config__, 'agf2_weight_tol', 1e-11)
self.fock_diis_space = getattr(__config__, 'agf2_diis_space', 6)
self.fock_diis_min_space = getattr(__config__, 'agf2_diis_min_space', 1)
self.diis = getattr(__config__, 'agf2_diis', True)
self.diis_space = getattr(__config__, 'agf2_diis_space', 8)
self.diis_min_space = getattr(__config__, 'agf2_diis_min_space', 1)
self.os_factor = getattr(__config__, 'agf2_os_factor', 1.0)
self.ss_factor = getattr(__config__, 'agf2_ss_factor', 1.0)
self.damping = getattr(__config__, 'agf2_damping', 0.0)
self.mo_energy = mo_energy
self.mo_coeff = mo_coeff
self.mo_occ = mo_occ
self.se = None
self.gf = None
self.e_1b = mf.e_tot
self.e_2b = 0.0
self.e_init = 0.0
self.frozen = frozen
self._nmo = None
self._nocc = None
self.converged = False
self.chkfile = mf.chkfile
energy_1body = energy_1body
energy_2body = energy_2body
fock_loop = fock_loop
build_se_part = build_se_part
get_jk = get_jk
[docs]
def ao2mo(self, mo_coeff=None):
''' Get the electronic repulsion integrals in MO basis.
'''
# happens when e.g. restarting from chkfile
if self._scf._eri is None and self._scf._is_mem_enough():
self._scf._eri = self.mol.intor('int2e', aosym='s8')
mem_incore = ((self.nmo*(self.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):
''' Computes the one-body reduced density matrix in MO basis.
Kwargs:
gf : GreensFunction
Auxiliaries of the Green's function
Returns:
ndarray of density matrix
'''
if gf is None: gf = self.gf
if gf is None: gf = self.init_gf()
return gf.make_rdm1()
[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:
:class:`GreensFunction`, :class:`SelfEnergy`
'''
energy = self.mo_energy
coupling = np.eye(self.nmo)
chempot = binsearch_chempot(np.diag(energy), self.nmo, self.nocc*2)[0]
if frozen:
mask = get_frozen_mask(self)
energy = energy[mask]
coupling = coupling[:,mask]
gf = aux.GreensFunction(energy, coupling, chempot=chempot)
return gf
[docs]
def build_gf(self, eri=None, gf=None, se=None):
''' Builds the auxiliaries of the Green's function by solving
the Dyson equation.
Kwargs:
eri : _ChemistsERIs
Electronic repulsion integrals
gf : GreensFunction
Auxiliaries of the Green's function
se : SelfEnergy
Auxiliaries of the self-energy
Returns:
: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)
fock = self.get_fock(eri, gf)
return se.get_greens_function(fock)
[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 : 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:
: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 = dict(os_factor=os_factor, ss_factor=ss_factor)
gf_occ = gf.get_occupied()
gf_vir = gf.get_virtual()
if gf_occ.naux == 0 or gf_vir.naux == 0:
logger.warn(self, 'Attempting to build a self-energy with '
'no (i,j,a) or (a,b,i) configurations.')
se = aux.SelfEnergy([], [[],]*self.nmo, chempot=gf.chempot)
else:
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 = aux.combine(se_occ, se_vir)
if se_prev is not None and self.damping != 0.0:
se.coupling *= np.sqrt(1.0-self.damping)
se_prev.coupling *= np.sqrt(self.damping)
se = aux.combine(se, se_prev)
se = se.compress(n=(None,0))
return se
[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:
:class:`SelfEnergy`
'''
if diis is None:
return se
se_occ = se.get_occupied()
se_vir = se.get_virtual()
vv_occ = np.dot(se_occ.coupling, se_occ.coupling.T)
vv_vir = np.dot(se_vir.coupling, se_vir.coupling.T)
vev_occ = np.dot(se_occ.coupling * se_occ.energy[None], se_occ.coupling.T)
vev_vir = np.dot(se_vir.coupling * se_vir.energy[None], se_vir.coupling.T)
dat = np.array([vv_occ, vv_vir, vev_occ, vev_vir])
dat = diis.update(dat)
vv_occ, vv_vir, vev_occ, vev_vir = dat
se_occ = aux.SelfEnergy(*_agf2.cholesky_build(vv_occ, vev_occ), chempot=se.chempot)
se_vir = aux.SelfEnergy(*_agf2.cholesky_build(vv_vir, vev_vir), chempot=se.chempot)
se = aux.combine(se_occ, se_vir)
return se
[docs]
def energy_nuc(self):
return self._scf.energy_nuc()
[docs]
def dump_flags(self, verbose=None):
log = logger.new_logger(self, verbose)
log.info('')
log.info('******** %s ********', self.__class__)
log.info('conv_tol = %g', self.conv_tol)
log.info('conv_tol_rdm1 = %g', self.conv_tol_rdm1)
log.info('conv_tol_nelec = %g', self.conv_tol_nelec)
log.info('max_cycle = %g', self.max_cycle)
log.info('max_cycle_outer = %g', self.max_cycle_outer)
log.info('max_cycle_inner = %g', self.max_cycle_inner)
log.info('weight_tol = %g', self.weight_tol)
log.info('diis = %d', self.diis)
log.info('diis_space = %d', self.diis_space)
log.info('diis_min_space = %d', self.diis_min_space)
log.info('fock_diis_space = %d', self.fock_diis_space)
log.info('fock_diis_min_space = %d', self.fock_diis_min_space)
log.info('os_factor = %g', self.os_factor)
log.info('ss_factor = %g', self.ss_factor)
log.info('damping = %g', self.damping)
log.info('nmo = %s', self.nmo)
log.info('nocc = %s', self.nocc)
if self.frozen is not None:
log.info('frozen orbitals = %s', self.frozen)
log.info('max_memory %d MB (current use %d MB)',
self.max_memory, lib.current_memory()[0])
return self
def _finalize(self):
''' Hook for dumping results and clearing up the object.
'''
if self.converged:
logger.info(self, '%s converged', self.__class__.__name__)
else:
logger.note(self, '%s not converged', self.__class__.__name__)
ip = self.get_ip(self.gf, nroots=1)[0][0]
ea = self.get_ea(self.gf, nroots=1)[0][0]
logger.note(self, 'E(%s) = %.16g E_corr = %.16g',
self.__class__.__name__, self.e_tot, self.e_corr)
logger.note(self, 'IP = %.16g EA = %.16g', ip, ea)
logger.note(self, 'Quasiparticle gap = %.16g', ip+ea)
return self
[docs]
def reset(self, mol=None):
if mol is not None:
self.mol = mol
self._scf.reset(mol)
return self
[docs]
def kernel(self, eri=None, gf=None, se=None, dump_chk=True):
if self.verbose >= logger.WARN:
self.check_sanity()
self.dump_flags()
if eri is None: eri = self.ao2mo()
if gf is None: gf = self.gf
if se is None: se = self.se
if gf is None:
gf = self.init_gf()
gf_froz = self.init_gf(frozen=True)
else:
gf_froz = gf
if se is None:
se = self.build_se(eri, gf_froz)
self.converged, self.e_1b, self.e_2b, self.gf, self.se = \
kernel(self, eri=eri, gf=gf, se=se, verbose=self.verbose, dump_chk=dump_chk)
self._finalize()
return self.converged, self.e_1b, self.e_2b, self.gf, self.se
[docs]
def dump_chk(self, chkfile=None, key='agf2', gf=None, se=None,
frozen=None, nmom=None,
mo_energy=None, mo_coeff=None, mo_occ=None):
if chkfile is None:
chkfile = self.chkfile
if not chkfile:
return self
chkutil.dump_agf2(self, chkfile, key,
gf, se, frozen, None,
mo_energy, mo_coeff, mo_occ)
return self
[docs]
def update_from_chk_(self, chkfile=None, key='agf2'):
if chkfile is None:
chkfile = self.chkfile
mol, agf2_dict = chkutil.load_agf2(chkfile, key)
self.__dict__.update(agf2_dict)
return self
update = update_from_chk = update_from_chk_
[docs]
def density_fit(self, auxbasis=None, with_df=None):
from pyscf.agf2 import dfragf2
myagf2 = dfragf2.DFRAGF2(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.get_occupied()
e_ip = list(-gf_occ.energy[-nroots:])[::-1]
v_ip = list(gf_occ.coupling[:,-nroots:].T)[::-1]
return e_ip, v_ip
[docs]
def ipagf2(self, nroots=5):
''' Find the (N-1)-electron charged excitations, corresponding
to the largest :attr:`nroots` poles of the occupied
Green's function.
Kwargs:
nroots : int
Number of roots (poles) requested. Default 1.
Returns:
IP and transition moment (float, 1D array) if :attr:`nroots`
= 1, or array of IPs and moments (1D array, 2D array) if
:attr:`nroots` > 1.
'''
e_ip, v_ip = self.get_ip(self.gf, nroots=nroots)
for n, en, vn in zip(range(nroots), e_ip, v_ip):
qpwt = np.linalg.norm(vn)**2
logger.note(self, 'IP energy level %d E = %.16g QP weight = %0.6g', n, en, qpwt)
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.get_virtual()
e_ea = list(gf_vir.energy[:nroots])
v_ea = list(gf_vir.coupling[:,:nroots].T)
return e_ea, v_ea
[docs]
def eaagf2(self, nroots=5):
''' Find the (N+1)-electron charged excitations, corresponding
to the smallest :attr:`nroots` poles of the virtual
Green's function.
Kwargs:
See ipagf2()
'''
e_ea, v_ea = self.get_ea(self.gf, nroots=nroots)
for n, en, vn in zip(range(nroots), e_ea, v_ea):
qpwt = np.linalg.norm(vn)**2
logger.note(self, 'EA energy level %d E = %.16g QP weight = %0.6g', n, en, qpwt)
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)
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.size
return self._nmo
@nmo.setter
def nmo(self, val):
self._nmo = val
@property
def e_tot(self):
return self.e_1b + self.e_2b
@property
def e_corr(self):
# TODO Should HF energy be recalculated in case DFT orbitals or so were used?
e_hf = mpi_helper.bcast(self._scf.e_tot)
return self.e_tot - e_hf
@property
def qmo_energy(self):
return self.gf.energy
@property
def qmo_coeff(self):
''' Gives the couplings in AO basis '''
return np.dot(self.mo_coeff, self.gf.coupling)
@property
def qmo_occ(self):
coeff = self.gf.get_occupied().coupling
occ = 2.0 * np.linalg.norm(coeff, axis=0) ** 2
vir = np.zeros_like(self.gf.get_virtual().energy)
qmo_occ = np.concatenate([occ, vir])
return qmo_occ
[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.nmo = None
self.nocc = 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()
fock_ao = h1e_ao + agf2._scf.get_veff(agf2.mol, dm)
self.h1e = np.dot(np.dot(mo_coeff.conj().T, h1e_ao), mo_coeff)
self.fock = np.dot(np.dot(mo_coeff.conj().T, fock_ao), mo_coeff)
self.h1e = mpi_helper.bcast(self.h1e)
self.fock = mpi_helper.bcast(self.fock)
self.e_hf = mpi_helper.bcast(agf2._scf.e_tot)
self.nmo = agf2.nmo
self.nocc = agf2.nocc
self.mol = agf2.mol
mo_e = self.fock.diagonal()
gap = abs(mo_e[:self.nocc,None] - mo_e[None,self.nocc:]).min()
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)
eri = ao2mo.incore.full(agf2._scf._eri, eris.mo_coeff, verbose=log)
eri = ao2mo.addons.restore('s4', eri, eris.nmo)
eris.eri = eri
log.timer('MO integral transformation', *cput0)
return eris
def _make_mo_eris_outcore(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)
mol = agf2.mol
mo_coeff = np.asarray(eris.mo_coeff, order='F')
eris.feri = lib.H5TmpFile()
ao2mo.outcore.full(mol, mo_coeff, eris.feri, dataname='mo',
max_memory=agf2.max_memory, verbose=log)
eris.eri = eris.feri['mo']
log.timer('MO integral transformation', *cput0)
return eris
def _make_qmo_eris_incore(agf2, eri, coeffs):
''' Returns ndarray
'''
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(agf2.stdout, agf2.verbose)
cx = np.eye(eri.nmo)
if not (agf2.frozen is None or agf2.frozen == 0):
mask = get_frozen_mask(agf2)
cx = cx[:,mask]
coeffs = (cx,) + coeffs
shape = tuple(x.shape[1] for x in coeffs)
qeri = ao2mo.incore.general(eri.eri, coeffs, compact=False, verbose=log)
qeri = qeri.reshape(shape)
log.timer('QMO integral transformation', *cput0)
return qeri
def _make_qmo_eris_outcore(agf2, eri, coeffs):
''' Returns H5 dataset
'''
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(agf2.stdout, agf2.verbose)
nmo = eri.nmo
ci, cj, ca = coeffs
ni = ci.shape[1]
nj = cj.shape[1]
na = ca.shape[1]
npair = nmo*(nmo+1)//2
mask = get_frozen_mask(agf2)
frozen = np.sum(~mask)
# possible to have incore MO, outcore QMO
if getattr(eri, 'feri', None) is None:
eri.feri = lib.H5TmpFile()
elif 'qmo' in eri.feri:
del eri.feri['qmo']
eri.feri.create_dataset('qmo', (nmo-frozen, ni, nj, na), 'f8')
blksize = _agf2.get_blksize(agf2.max_memory, (nmo*npair, nj*na, npair), (nmo*ni, nj*na))
blksize = min(nmo, max(BLKMIN, blksize))
log.debug1('blksize (ragf2._make_qmo_eris_outcore) = %d', blksize)
tril2sq = lib.square_mat_in_trilu_indices(nmo)
q1 = 0
for p0, p1 in lib.prange(0, nmo, blksize):
if not np.any(mask[p0:p1]):
# block is fully frozen
continue
inds = np.arange(p0, p1)[mask[p0:p1]]
q0, q1 = q1, q1 + len(inds)
idx = list(np.concatenate(tril2sq[inds]))
buf = eri.eri[idx] # (blk, nmo, npair)
buf = buf.reshape((q1-q0)*nmo, -1) # (blk*nmo, npair)
jasym, nja, cja, sja = ao2mo.incore._conc_mos(cj, ca, compact=True)
buf = ao2mo._ao2mo.nr_e2(buf, cja, sja, 's2kl', 's1')
buf = buf.reshape(q1-q0, nmo, nj, na)
buf = lib.einsum('xpja,pi->xija', buf, ci)
eri.feri['qmo'][q0:q1] = np.asarray(buf, order='C')
log.timer('QMO integral transformation', *cput0)
return eri.feri['qmo']
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', verbose=3)
rhf = scf.RHF(mol)
rhf.conv_tol = 1e-11
rhf.run()
ragf2 = RAGF2(rhf, frozen=0)
ragf2.run()
ragf2.ipagf2(nroots=5)
ragf2.eaagf2(nroots=5)
print(mp.MP2(rhf, frozen=ragf2.frozen).run(verbose=0).e_corr)
print(ragf2.e_init)
ragf2 = ragf2.density_fit()
ragf2.run()