# 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 with density fitting
'''
import numpy as np
import ctypes
from pyscf import lib
from pyscf.lib import logger
from pyscf import __config__
from pyscf import ao2mo, df
from pyscf.agf2 import uagf2, dfragf2, mpi_helper, _agf2
from pyscf.agf2 import aux_space as aux
BLKMIN = getattr(__config__, 'agf2_blkmin', 100)
[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[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
nmoa, nmob = eri.nmo
nocca, nvira = gf_occ[0].naux, gf_vir[0].naux
noccb, nvirb = gf_occ[1].naux, gf_vir[1].naux
naux = agf2.with_df.get_naoaux()
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
qeri = _make_qmo_eris_incore(agf2, eri, (ci_a, ci_a, ca_a), (ci_b, ci_b, ca_b))
(qxi_a, qja_a), (qxi_b, qja_b) = qeri
qxi = (qxi_a, qxi_b)
qja = (qja_a, qja_b)
himem_required = naux*(nvira+nmoa) + (nocca*nvira+noccb*nvirb)*(1+2*nmoa) + (2*nmoa**2)
himem_required *= 8e-6
himem_required *= lib.num_threads()
if ((himem_required*1.05 + lib.current_memory()[0]) > agf2.max_memory
and agf2.allow_lowmem_build) or agf2.allow_lowmem_build == 'force':
log.debug('Thread-private memory overhead %.3f exceeds max_memory, using '
'low-memory version.', himem_required)
build_mats_dfuagf2 = _agf2.build_mats_dfuagf2_lowmem
else:
build_mats_dfuagf2 = _agf2.build_mats_dfuagf2_incore
vv, vev = build_mats_dfuagf2(qxi, qja, (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 = uagf2.get_frozen_mask(agf2)
coupling = np.zeros((nmoa, 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)
himem_required = naux*(nvirb+nmob) + (noccb*nvirb+nocca*nvira)*(1+2*nmob) + (2*nmob**2)
himem_required *= 8e-6
himem_required *= lib.num_threads()
if ((himem_required*1.05 + lib.current_memory()[0]) > agf2.max_memory
and agf2.allow_lowmem_build) or agf2.allow_lowmem_build == 'force':
log.debug('Thread-private memory overhead %.3f exceeds max_memory, using '
'low-memory version.', himem_required)
build_mats_dfuagf2 = _agf2.build_mats_dfuagf2_lowmem
else:
build_mats_dfuagf2 = _agf2.build_mats_dfuagf2_incore
rv = np.s_[::-1]
vv, vev = build_mats_dfuagf2(qxi[rv], qja[rv], (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 = uagf2.get_frozen_mask(agf2)
coupling = np.zeros((nmoa, 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]
class DFUAGF2(uagf2.UAGF2):
''' Unrestricted AGF2 with canonical HF reference with density fitting
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.
allow_lowmem_build : bool
Allow the self-energy build to switch to a serially slower
code with lower thread-private memory overhead if needed. One
of True, False or 'force'. Default value is True.
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
'''
_keys = {'_with_df', 'allow_lowmem_build'}
def __init__(self, mf, frozen=None, mo_energy=None, mo_coeff=None, mo_occ=None):
uagf2.UAGF2.__init__(self, mf, frozen=frozen, mo_energy=mo_energy,
mo_coeff=mo_coeff, mo_occ=mo_occ)
if getattr(mf, 'with_df', None) is not None:
self.with_df = mf.with_df
else:
self.with_df = df.DF(mf.mol)
self.with_df.auxbasis = df.make_auxbasis(mf.mol, mp2fit=True)
self.allow_lowmem_build = True
build_se_part = build_se_part
get_jk = dfragf2.get_jk
[docs]
def ao2mo(self, mo_coeff=None):
''' Get the density-fitted electronic repulsion integrals in
MO basis.
'''
eri = _make_mo_eris_incore(self, mo_coeff)
return eri
[docs]
def reset(self, mol=None):
self.with_df.reset(mol)
return uagf2.UAGF2.reset(self, mol)
@property
def with_df(self):
return self._with_df
@with_df.setter
def with_df(self, val):
self._with_df = val
self._with_df.__class__ = dfragf2.DF
class _ChemistsERIs(uagf2._ChemistsERIs):
''' (pq|rs) as (pq|J)(J|rs)
MO tensors are stored in tril from, we only need QMO tensors
in low-symmetry
'''
pass
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)
with_df = agf2.with_df
moa, mob = eris.mo_coeff
nmoa, nmob = moa.shape[1], mob.shape[1]
npaira, npairb = nmoa*(nmoa+1)//2, nmob*(nmob+1)//2
naux = with_df.get_naoaux()
qxy_a = np.zeros((naux, npaira))
qxy_b = np.zeros((naux, npairb))
moa = np.asarray(moa, order='F')
mob = np.asarray(mob, order='F')
sija = (0, nmoa, 0, nmoa)
sijb = (0, nmob, 0, nmob)
sym = {'aosym': 's2', 'mosym': 's2'}
for p0, p1 in with_df.prange():
eri0 = with_df._cderi[p0:p1]
qxy_a[p0:p1] = ao2mo._ao2mo.nr_e2(eri0, moa, sija, out=qxy_a[p0:p1], **sym)
qxy_b[p0:p1] = ao2mo._ao2mo.nr_e2(eri0, mob, sijb, out=qxy_b[p0:p1], **sym)
mpi_helper.barrier()
mpi_helper.allreduce_safe_inplace(qxy_a)
mpi_helper.allreduce_safe_inplace(qxy_b)
eris.eri_a = qxy_a
eris.eri_b = qxy_b
eris.eri_aa = (eris.eri_a, eris.eri_a)
eris.eri_ab = (eris.eri_a, eris.eri_b)
eris.eri_ba = (eris.eri_b, eris.eri_a)
eris.eri_bb = (eris.eri_b, eris.eri_b)
eris.eri = (eris.eri_a, eris.eri_b)
log.timer('MO integral transformation', *cput0)
return eris
def _make_qmo_eris_incore(agf2, eri, coeffs_a, coeffs_b):
''' Returns nested tuple of ndarray
'''
cput0 = (logger.process_clock(), logger.perf_counter())
log = logger.Logger(agf2.stdout, agf2.verbose)
cxa, cxb = np.eye(agf2.nmo[0]), np.eye(agf2.nmo[1])
if not (agf2.frozen is None or agf2.frozen == 0):
mask = uagf2.get_frozen_mask(agf2)
cxa = cxa[:,mask[0]]
cxb = cxb[:,mask[1]]
nmoa, nmob = agf2.nmo
npaira, npairb = nmoa*(nmoa+1)//2, nmob*(nmob+1)//2
with_df = agf2.with_df
naux = with_df.get_naoaux()
cia, cja, caa = coeffs_a
cib, cjb, cab = coeffs_b
xisym_a, nxi_a, cxi_a, sxi_a = ao2mo.incore._conc_mos(cxa, cia, compact=False)
jasym_a, nja_a, cja_a, sja_a = ao2mo.incore._conc_mos(cja, caa, compact=False)
xisym_b, nxi_b, cxi_b, sxi_b = ao2mo.incore._conc_mos(cxb, cib, compact=False)
jasym_b, nja_b, cja_b, sja_b = ao2mo.incore._conc_mos(cjb, cab, compact=False)
sym = {'aosym': 's2', 'mosym': 's1'}
qxi_a = np.zeros((naux, nxi_a))
qxi_b = np.zeros((naux, nxi_b))
qja_a = np.zeros((naux, nja_a))
qja_b = np.zeros((naux, nja_b))
buf = (np.zeros((with_df.blockdim, npaira)), np.zeros((with_df.blockdim, npairb)))
for p0, p1 in mpi_helper.prange(0, naux, with_df.blockdim):
naux0 = p1 - p0
bufa0 = buf[0][:naux0]
bufb0 = buf[1][:naux0]
bufa0[:] = eri.eri[0][p0:p1]
bufb0[:] = eri.eri[1][p0:p1]
qxi_a[p0:p1] = ao2mo._ao2mo.nr_e2(bufa0, cxi_a, sxi_a, out=qxi_a[p0:p1], **sym)
qxi_b[p0:p1] = ao2mo._ao2mo.nr_e2(bufb0, cxi_b, sxi_b, out=qxi_b[p0:p1], **sym)
qja_a[p0:p1] = ao2mo._ao2mo.nr_e2(bufa0, cja_a, sja_a, out=qja_a[p0:p1], **sym)
qja_b[p0:p1] = ao2mo._ao2mo.nr_e2(bufb0, cja_b, sja_b, out=qja_b[p0:p1], **sym)
mpi_helper.barrier()
mpi_helper.allreduce_safe_inplace(qxi_a)
mpi_helper.allreduce_safe_inplace(qxi_b)
mpi_helper.allreduce_safe_inplace(qja_a)
mpi_helper.allreduce_safe_inplace(qja_b)
qxi_a = qxi_a.reshape(naux, -1)
qxi_b = qxi_b.reshape(naux, -1)
qja_a = qja_a.reshape(naux, -1)
qja_b = qja_b.reshape(naux, -1)
log.timer('QMO integral transformation', *cput0)
return ((qxi_a, qja_a), (qxi_b, qja_b))
if __name__ == '__main__':
from pyscf import gto, scf, mp
import pyscf.scf.stability
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).density_fit()
uhf.conv_tol = 1e-11
uhf.run()
for niter in range(1, 11):
stability = scf.stability.uhf_stability(uhf)
if isinstance(stability, tuple):
sint, sext = stability
else:
sint = stability
if np.allclose(sint, uhf.mo_coeff):
break
else:
rdm1 = uhf.make_rdm1(sint, uhf.mo_occ)
uhf.scf(dm0=rdm1)
uagf2 = DFUAGF2(uhf)
uagf2.run()
uagf2.ipagf2(nroots=5)
uagf2.eaagf2(nroots=5)
