Source code for pyscf.gw.gw_exact

#!/usr/bin/env python
# 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.
#
# Authors: Timothy Berkelbach <tim.berkelbach@gmail.com>
#          Qiming Sun <osirpt.sun@gmail.com>
#

'''
Spin-restricted G0W0 approximation with exact frequency integration
'''


from functools import reduce
import numpy
import numpy as np
from scipy.optimize import newton

from pyscf import lib
from pyscf.lib import logger
from pyscf import ao2mo
from pyscf import dft
from pyscf.mp.mp2 import get_nocc, get_nmo, get_frozen_mask, _mo_without_core
from pyscf import __config__

einsum = lib.einsum

BLKMIN = getattr(__config__, 'cc_ccsd_blkmin', 4)
MEMORYMIN = getattr(__config__, 'cc_ccsd_memorymin', 2000)

[docs] def kernel(gw, mo_energy, mo_coeff, td_e, td_xy, eris=None, orbs=None, verbose=logger.NOTE): '''GW-corrected quasiparticle orbital energies Returns: A list : converged, mo_energy, mo_coeff ''' # mf must be DFT; for HF use xc = 'hf' mf = gw._scf assert (isinstance(mf, (dft.rks.RKS , dft.uks.UKS, dft.roks.ROKS , dft.uks.UKS, dft.rks_symm.RKS , dft.uks_symm.UKS, dft.rks_symm.ROKS, dft.uks_symm.UKS))) assert (gw.frozen == 0 or gw.frozen is None) if eris is None: eris = gw.ao2mo(mo_coeff) if orbs is None: orbs = range(gw.nmo) v_mf = mf.get_veff() - mf.get_j() v_mf = reduce(numpy.dot, (mo_coeff.T, v_mf, mo_coeff)) nocc = gw.nocc nmo = gw.nmo nvir = nmo-nocc vk_oo = -np.einsum('piiq->pq', eris.oooo) vk_ov = -np.einsum('iqpi->pq', eris.ovoo) vk_vv = -np.einsum('ipqi->pq', eris.ovvo).conj() vk = np.block([[vk_oo, vk_ov],[vk_ov.T, vk_vv]]) nexc = len(td_e) # factor of 2 for normalization, see tdscf/rhf.py td_xy = 2*np.asarray(td_xy) # (nexc, 2, nocc, nvir) td_z = np.sum(td_xy, axis=1).reshape(nexc,nocc,nvir) tdm_oo = einsum('via,iapq->vpq', td_z, eris.ovoo) tdm_ov = einsum('via,iapq->vpq', td_z, eris.ovov) tdm_vv = einsum('via,iapq->vpq', td_z, eris.ovvv) tdm = [] for oo,ov,vv in zip(tdm_oo,tdm_ov,tdm_vv): tdm.append(np.block([[oo, ov],[ov.T, vv]])) tdm = np.asarray(tdm) conv = True mo_energy = np.zeros_like(gw._scf.mo_energy) for p in orbs: tdm_p = tdm[:,:,p] if gw.linearized: ep = gw._scf.mo_energy[p] sigma = get_sigma_element(gw, ep, tdm_p, tdm_p, td_e).real dsigma_dw = get_sigma_deriv_element(gw, ep, tdm_p, tdm_p, td_e).real zn = 1.0/(1-dsigma_dw) mo_energy[p] = ep + zn*(sigma.real + vk[p,p] - v_mf[p,p]) else: def quasiparticle(omega): sigma = get_sigma_element(gw, omega, tdm_p, tdm_p, td_e) return omega - gw._scf.mo_energy[p] - (sigma.real + vk[p,p] - v_mf[p,p]) try: mo_energy[p] = newton(quasiparticle, gw._scf.mo_energy[p], tol=1e-6, maxiter=100) except RuntimeError: conv = False mo_energy[p] = gw._scf.mo_energy[p] logger.warn(gw, 'Root finding for GW eigenvalue %s did not converge. ' 'Setting it equal to the reference MO energy.'%(p)) mo_coeff = gw._scf.mo_coeff if gw.verbose >= logger.DEBUG: numpy.set_printoptions(threshold=nmo) logger.debug(gw, ' GW mo_energy =\n%s', mo_energy) numpy.set_printoptions(threshold=1000) return conv, mo_energy, mo_coeff
[docs] def get_sigma_element(gw, omega, tdm_p, tdm_q, td_e, eta=None, vir_sgn=1): if eta is None: eta = gw.eta nocc = gw.nocc evi = lib.direct_sum('v-i->vi', td_e, gw._scf.mo_energy[:nocc]) eva = lib.direct_sum('v+a->va', td_e, gw._scf.mo_energy[nocc:]) sigma = np.sum( tdm_p[:,:nocc]*tdm_q[:,:nocc]/(omega + evi - 1j*eta) ) sigma += np.sum( tdm_p[:,nocc:]*tdm_q[:,nocc:]/(omega - eva + vir_sgn*1j*eta) ) return sigma
[docs] def get_sigma_deriv_element(gw, omega, tdm_p, tdm_q, td_e, eta=None, vir_sgn=1): if eta is None: eta = gw.eta nocc = gw.nocc evi = lib.direct_sum('v-i->vi', td_e, gw._scf.mo_energy[:nocc]) eva = lib.direct_sum('v+a->va', td_e, gw._scf.mo_energy[nocc:]) dsigma = -np.sum( tdm_p[:,:nocc]*tdm_q[:,:nocc]/(omega + evi - 1j*eta)**2 ) dsigma += -np.sum( tdm_p[:,nocc:]*tdm_q[:,nocc:]/(omega - eva + vir_sgn*1j*eta)**2 ) return dsigma
[docs] def get_g(omega, mo_energy, mo_occ, eta): sgn = mo_occ - 1 return 1.0/(omega - mo_energy + 1j*eta*sgn)
[docs] class GWExact(lib.StreamObject): '''non-relativistic restricted GW Saved results mo_energy : Orbital energies mo_coeff Orbital coefficients ''' eta = getattr(__config__, 'gw_gw_GW_eta', 1e-8) linearized = getattr(__config__, 'gw_gw_GW_linearized', False) _keys = { 'eta', 'linearized', 'mol', 'frozen', 'mo_energy', 'mo_coeff', 'mo_occ', } def __init__(self, mf, frozen=None, tdmf=None): self.mol = mf.mol self._scf = mf self._tdscf = tdmf self.verbose = self.mol.verbose self.stdout = self.mol.stdout self.max_memory = mf.max_memory self.frozen = frozen ################################################## # don't modify the following attributes, they are not input options self._nocc = None self._nmo = None self.mo_energy = None self.mo_coeff = mf.mo_coeff self.mo_occ = mf.mo_occ
[docs] def dump_flags(self, verbose=None): log = logger.new_logger(self, verbose) log.info('') log.info('******** %s ********', self.__class__) log.info('method = %s', self.__class__.__name__) nocc = self.nocc nvir = self.nmo - nocc log.info('GW nocc = %d, nvir = %d', nocc, nvir) if self.frozen is not None: log.info('frozen = %s', self.frozen) logger.info(self, 'use perturbative linearized QP eqn = %s', self.linearized) return self
@property def nocc(self): return self.get_nocc() @nocc.setter def nocc(self, n): self._nocc = n @property def nmo(self): return self.get_nmo() @nmo.setter def nmo(self, n): self._nmo = n get_nocc = get_nocc get_nmo = get_nmo get_frozen_mask = get_frozen_mask
[docs] def get_g0(self, omega, eta=None): if eta is None: eta = self.eta return get_g(omega, self._scf.mo_energy, self.mo_occ, eta)
[docs] def get_g(self, omega, eta=None): if eta is None: eta = self.eta return get_g(omega, self.mo_energy, self.mo_occ, eta)
[docs] def kernel(self, mo_energy=None, mo_coeff=None, td_e=None, td_xy=None, eris=None, orbs=None): if mo_coeff is None: mo_coeff = self._scf.mo_coeff if mo_energy is None: mo_energy = self._scf.mo_energy if self._tdscf is None: from pyscf import tdscf self._tdscf = tdscf.dRPA(self._scf) nocc, nvir = self.nocc, self.nmo-self.nocc self._tdscf.nstates = nocc*nvir self._tdscf.verbose = 0 self._tdscf.kernel() if td_e is None: td_e = self._tdscf.e if td_xy is None: td_xy = self._tdscf.xy cput0 = (logger.process_clock(), logger.perf_counter()) self.dump_flags() self.converged, self.mo_energy, self.mo_coeff = \ kernel(self, mo_energy, mo_coeff, td_e, td_xy, eris=eris, orbs=orbs, verbose=self.verbose) logger.timer(self, 'GW', *cput0) return self.mo_energy
[docs] def reset(self, mol=None): if mol is not None: self.mol = mol self._scf.reset(mol) self._tdscf.reset(mol) return self
[docs] def ao2mo(self, mo_coeff=None): nmo = self.nmo nao = self.mo_coeff.shape[0] nmo_pair = nmo * (nmo+1) // 2 nao_pair = nao * (nao+1) // 2 mem_incore = (max(nao_pair**2, nmo**4) + nmo_pair**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.mol.incore_anyway): return _make_eris_incore(self, mo_coeff) elif getattr(self._scf, 'with_df', None): logger.warn(self, 'GW Exact detected DF being used in the HF object. ' 'MO integrals are computed based on the DF 3-index tensors.\n' 'Developer TODO: Write dfgw.GWExact for the ' 'DF-GW calculations') raise NotImplementedError #return _make_df_eris_outcore(self, mo_coeff) else: return _make_eris_outcore(self, mo_coeff)
class _ChemistsERIs: '''(pq|rs) Identical to rccsd _ChemistsERIs except no vvvv.''' def __init__(self, mol=None): self.mol = mol self.mo_coeff = None self.nocc = None self.fock = None self.oooo = None self.ovoo = None self.oovv = None self.ovvo = None self.ovov = None self.ovvv = None def _common_init_(self, mycc, mo_coeff=None): if mo_coeff is None: mo_coeff = mycc.mo_coeff self.mo_coeff = mo_coeff = _mo_without_core(mycc, mo_coeff) # Note: Recomputed fock matrix since SCF may not be fully converged. dm = mycc._scf.make_rdm1(mycc.mo_coeff, mycc.mo_occ) fockao = mycc._scf.get_hcore() + mycc._scf.get_veff(mycc.mol, dm) self.fock = reduce(numpy.dot, (mo_coeff.conj().T, fockao, mo_coeff)) self.nocc = mycc.nocc self.mol = mycc.mol mo_e = self.fock.diagonal() try: gap = abs(mo_e[:self.nocc,None] - mo_e[None,self.nocc:]).min() if gap < 1e-5: logger.warn(mycc, 'HOMO-LUMO gap %s too small for GW', gap) except ValueError: # gap.size == 0 pass return self def _make_eris_incore(mycc, mo_coeff=None, ao2mofn=None): cput0 = (logger.process_clock(), logger.perf_counter()) eris = _ChemistsERIs() eris._common_init_(mycc, mo_coeff) nocc = eris.nocc nmo = eris.fock.shape[0] if callable(ao2mofn): eri1 = ao2mofn(eris.mo_coeff).reshape([nmo]*4) else: eri1 = ao2mo.incore.full(mycc._scf._eri, eris.mo_coeff) eri1 = ao2mo.restore(1, eri1, nmo) eris.oooo = eri1[:nocc,:nocc,:nocc,:nocc].copy() eris.ovoo = eri1[:nocc,nocc:,:nocc,:nocc].copy() eris.ovov = eri1[:nocc,nocc:,:nocc,nocc:].copy() eris.oovv = eri1[:nocc,:nocc,nocc:,nocc:].copy() eris.ovvo = eri1[:nocc,nocc:,nocc:,:nocc].copy() eris.ovvv = eri1[:nocc,nocc:,nocc:,nocc:].copy() logger.timer(mycc, 'GW integral transformation', *cput0) return eris def _make_eris_outcore(mycc, mo_coeff=None): cput0 = (logger.process_clock(), logger.perf_counter()) log = logger.Logger(mycc.stdout, mycc.verbose) eris = _ChemistsERIs() eris._common_init_(mycc, mo_coeff) mol = mycc.mol mo_coeff = eris.mo_coeff nocc = eris.nocc nao, nmo = mo_coeff.shape nvir = nmo - nocc eris.feri1 = lib.H5TmpFile() eris.oooo = eris.feri1.create_dataset('oooo', (nocc,nocc,nocc,nocc), 'f8') eris.ovoo = eris.feri1.create_dataset('ovoo', (nocc,nvir,nocc,nocc), 'f8', chunks=(nocc,1,nocc,nocc)) eris.ovov = eris.feri1.create_dataset('ovov', (nocc,nvir,nocc,nvir), 'f8', chunks=(nocc,1,nocc,nvir)) eris.ovvo = eris.feri1.create_dataset('ovvo', (nocc,nvir,nvir,nocc), 'f8', chunks=(nocc,1,nvir,nocc)) eris.ovvv = eris.feri1.create_dataset('ovvv', (nocc,nvir,nvir,nvir), 'f8') eris.oovv = eris.feri1.create_dataset('oovv', (nocc,nocc,nvir,nvir), 'f8', chunks=(nocc,nocc,1,nvir)) max_memory = max(MEMORYMIN, mycc.max_memory-lib.current_memory()[0]) ftmp = lib.H5TmpFile() ao2mo.full(mol, mo_coeff, ftmp, max_memory=max_memory, verbose=log) eri = ftmp['eri_mo'] nocc_pair = nocc*(nocc+1)//2 tril2sq = lib.square_mat_in_trilu_indices(nmo) oo = eri[:nocc_pair] eris.oooo[:] = ao2mo.restore(1, oo[:,:nocc_pair], nocc) oovv = lib.take_2d(oo, tril2sq[:nocc,:nocc].ravel(), tril2sq[nocc:,nocc:].ravel()) eris.oovv[:] = oovv.reshape(nocc,nocc,nvir,nvir) oo = oovv = None tril2sq = lib.square_mat_in_trilu_indices(nmo) blksize = min(nvir, max(BLKMIN, int(max_memory*1e6/8/nmo**3/2))) for p0, p1 in lib.prange(0, nvir, blksize): q0, q1 = p0+nocc, p1+nocc off0 = q0*(q0+1)//2 off1 = q1*(q1+1)//2 buf = lib.unpack_tril(eri[off0:off1]) tmp = buf[ tril2sq[q0:q1,:nocc] - off0 ] eris.ovoo[:,p0:p1] = tmp[:,:,:nocc,:nocc].transpose(1,0,2,3) eris.ovvo[:,p0:p1] = tmp[:,:,nocc:,:nocc].transpose(1,0,2,3) eris.ovov[:,p0:p1] = tmp[:,:,:nocc,nocc:].transpose(1,0,2,3) eris.ovvv[:,p0:p1] = tmp[:,:,nocc:,nocc:].transpose(1,0,2,3) buf = tmp = None log.timer('GW integral transformation', *cput0) return eris if __name__ == '__main__': from pyscf import gto mol = gto.Mole() mol.verbose = 5 mol.atom = [ [8 , (0. , 0. , 0.)], [1 , (0. , -0.757 , 0.587)], [1 , (0. , 0.757 , 0.587)]] mol.basis = 'cc-pvdz' mol.build() mf = dft.RKS(mol) mf.xc = 'hf' mf.kernel() gw = GWExact(mf) gw.kernel() print(gw.mo_energy) # [-20.10555941 -1.2264133 -0.68160937 -0.53066324 -0.44679866 # 0.17291986 0.24457082 0.74758225 0.80045126 1.11748749 # 1.15083528 1.19081826 1.40406946 1.43593671 1.63324755 # 1.79839248 1.88459324 2.31461977 2.48839545 3.26047431 # 3.32486673 3.49601314 3.77699489 4.14575936] nocc = mol.nelectron//2 gw.linearized = True gw.kernel(orbs=[nocc-1,nocc]) print(gw.mo_energy[nocc-1] - -0.44684106) print(gw.mo_energy[nocc] - 0.17292032)