#!/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: Tianyu Zhu <zhutianyu1991@gmail.com>
#
"""
Spin-restricted Quasiparticle Self-consistent MP2
"""
import time, h5py
from functools import reduce
import numpy, scipy
import numpy as np
from pyscf import lib
from pyscf.lib import logger
from pyscf import ao2mo
from pyscf import df, dft, scf
from pyscf.mp.mp2 import MP2, get_nocc, get_nmo, get_frozen_mask, _mo_energy_without_core, _mo_without_core, _mem_usage
from pyscf import __config__
from fcdmft.gw.mol.gw_ac import get_g0
einsum = lib.einsum
[docs]
def kernel(mp2, mo_energy, mo_coeff):
"""QSMP2-corrected quasiparticle orbital energies
Parameters
----------
mp2 : MP2 object
mo_energy : double array
molecular orbital energies
mo_coeff : double ndarray
molecular orbital coefficients
Returns
-------
conv : bool
whether the quasiparticle equation is converged
mo_energy : double array
MP2-corrected quasiparticle orbital energies
mo_coeff : double ndarray
MP2-corrected quasiparticle orbital coefficients
"""
mf = mp2._scf
assert mp2.frozen is None
nmo = mp2.nmo
nocc = mp2.nocc
nvir = nmo - nocc
if mp2.ncore is None:
ncore = 0
else:
ncore = mp2.ncore
eta = mp2.eta
mp2.mo_energy = mp2._scf.mo_energy
# get hcore and ovlp
hcore = mf.get_hcore()
ovlp = mf.get_ovlp()
mp2_diis = lib.diis.DIIS(mp2, mp2.diis_file)
mp2_diis.space = mp2.diis_space
diis_start_cycle = mp2.diis_start_cycle
qsmp2_conv = False
cycle = 0
conv_pts = [mo_energy[nocc - 1], mo_energy[nocc]]
gf = get_g0(conv_pts, mo_energy, eta)
# cutoff for high virtuals
if mp2.ecut is None:
ncut = nmo
else:
ncut = len(mo_energy[mo_energy < mo_energy[nocc] + mp2.ecut])
while not qsmp2_conv and cycle < max(1, mp2.max_cycle):
# density matrix from updated QSMP2 density
dm_iter = 2.0 * np.dot(mo_coeff[:, :nocc], mo_coeff[:, :nocc].T)
# update eris
eris = mp2.ao2mo(mo_coeff)
# update veff
veff = mf.get_veff(mp2.mol, dm=dm_iter)
# compute QSMP2 self-energy
eia = mo_energy[:nocc, None] - mo_energy[None, nocc:]
sigma = np.zeros((nmo, nmo), dtype=np.complex128)
for p in range(ncore, ncut):
for q in range(ncore, ncut):
if p == q:
sigma[p, q] = get_sigma_mp2(mp2, eris, p, q, mo_energy[p] + 1j * eta, eia)
elif q > p:
sigma[p, q] += 0.5 * get_sigma_mp2(mp2, eris, p, q, mo_energy[p] + 1j * eta, eia)
sigma[p, q] += 0.5 * get_sigma_mp2(mp2, eris, p, q, mo_energy[q] + 1j * eta, eia)
sigma[q, p] = sigma[p, q]
# obtain vsig in AO basis
CS = np.dot(mo_coeff.T, ovlp)
vsig = reduce(np.dot, (CS.T, sigma.real, CS))
# complete Hamiltonian through DIIS
ham = veff + vsig
ham = mp2.run_diis(ham, cycle, mp2_diis)
ham = hcore + ham
# diagonalize
mo_energy, mo_coeff = scipy.linalg.eigh(ham, ovlp)
mp2.mo_energy = mo_energy
mp2.mo_coeff = mo_coeff
gf_old = gf.copy()
gf = get_g0(conv_pts, mo_energy, eta)
norm_dgf = np.linalg.norm(gf - gf_old) / nmo**2 / 2.0
if norm_dgf < mp2.conv_tol:
qsmp2_conv = True
logger.info(mp2, 'QSMP2 cycle= %d |dgf|= %4.3g', cycle + 1, norm_dgf)
logger.info(
mp2, ' HOMO = %.2f eV, LUMO = %.2f eV', mo_energy[nocc - 1] * 27.211386, mo_energy[nocc] * 27.211386
)
with np.printoptions(threshold=nmo):
logger.debug(mp2, ' QSMP2 mo_energy =\n%s', mo_energy)
cycle += 1
return qsmp2_conv, mo_energy, mo_coeff
[docs]
def get_sigma_mp2(mp2, eris, p, q, freq, eia):
"""Compute MP2 self-energy in MO basis on real axis
Parameters
----------
mp2 : QSMP2 object
eris : double ndarray
ERI in MO basis, nmo**4
p : int
row index of the self-energy matrix
q : int
column index of the self-energy matrix
freq : complex
frequency point
eia : double ndarray
occupied and virtual orbital energy difference
Returns
-------
sigma : double
self-energy matrix element
"""
mo_energy = mp2.mo_energy
iapb = eris.ovnv[:, :, p, :]
iaqb = eris.ovnv[:, :, q, :]
iajp = eris.ovon[:, :, :, p]
iajq = eris.ovon[:, :, :, q]
nocc = mp2.nocc
ewa = freq - mo_energy[nocc:]
eiw = mo_energy[:nocc] - freq
e_iab = lib.direct_sum('ia+b->iab', eia, ewa)
e_iaj = -lib.direct_sum('ia+j->iaj', eia, eiw)
t_iab = iapb / e_iab
t_iaj = iajp / e_iaj
sigma = 2 * einsum('iab,iab', t_iab, iaqb.conj())
sigma -= einsum('iab,iba', t_iab, iaqb.conj())
sigma += 2 * einsum('iaj,iaj', t_iaj, iajq.conj())
sigma -= einsum('iaj,jai', t_iaj, iajq.conj())
return sigma
# TO BE MURGED WITH PySCF _ChemistsERIs
class _ChemistsERIs:
def __init__(self, mol=None):
self.mol = mol
self.mo_coeff = None
self.nocc = None
self.fock = None
self.e_hf = None
self.orbspin = None
self.ovov = None
self.ovon = None
self.ovnv = None
def _common_init_(self, mp, mo_coeff=None):
if mo_coeff is None:
mo_coeff = mp.mo_coeff
if mo_coeff is None:
raise RuntimeError(
'mo_coeff, mo_energy are not initialized.\n' 'You may need to call mf.kernel() to generate them.'
)
self.mo_coeff = _mo_without_core(mp, mo_coeff)
self.mol = mp.mol
if mo_coeff is mp._scf.mo_coeff and mp._scf.converged:
# The canonical MP2 from a converged SCF result. Rebuilding fock
# and e_hf can be skipped
self.mo_energy = _mo_energy_without_core(mp, mp._scf.mo_energy)
self.fock = numpy.diag(self.mo_energy)
self.e_hf = mp._scf.e_tot
else:
dm = mp._scf.make_rdm1(mo_coeff, mp.mo_occ)
vhf = mp._scf.get_veff(mp.mol, dm)
fockao = mp._scf.get_fock(vhf=vhf, dm=dm)
self.fock = self.mo_coeff.conj().T.dot(fockao).dot(self.mo_coeff)
self.e_hf = mp._scf.energy_tot(dm=dm, vhf=vhf)
self.mo_energy = self.fock.diagonal().real
return self
def _make_eris(mp, mo_coeff=None, ao2mofn=None, verbose=None):
log = logger.new_logger(mp, verbose)
time0 = (time.process_time(), time.perf_counter())
eris = _ChemistsERIs()
eris._common_init_(mp, mo_coeff)
mo_coeff = eris.mo_coeff
nocc = mp.nocc
nmo = mp.nmo
nvir = nmo - nocc
mem_incore, mem_outcore, mem_basic = _mem_usage(nocc, nvir)
mem_now = lib.current_memory()[0]
max_memory = max(0, mp.max_memory - mem_now)
if max_memory < mem_basic:
log.warn(
'Not enough memory for integral transformation. ' 'Available mem %s MB, required mem %s MB',
max_memory,
mem_basic,
)
co = numpy.asarray(mo_coeff[:, :nocc], order='F')
cv = numpy.asarray(mo_coeff[:, nocc:], order='F')
if mp.mol.incore_anyway or (mp._scf._eri is not None and mem_incore < max_memory):
log.debug('transform (ia|jb) incore')
if callable(ao2mofn):
eris.ovon = ao2mofn((co, cv, co, mo_coeff))
eris.ovnv = ao2mofn((co, cv, mo_coeff, cv))
else:
eris.ovon = ao2mo.general(mp._scf._eri, (co, cv, co, mo_coeff)).reshape(nocc, nvir, nocc, nmo)
eris.ovnv = ao2mo.general(mp._scf._eri, (co, cv, mo_coeff, cv)).reshape(nocc, nvir, nmo, nvir)
elif getattr(mp._scf, 'with_df', None):
# To handle the PBC or custom 2-electron with 3-index tensor.
# Call dfmp2.MP2 for efficient DF-MP2 implementation.
log.warn(
'DF-HF is found. (ia|jb) is computed based on the DF '
'3-tensor integrals.\n'
'You can switch to dfmp2.MP2 for better performance'
)
log.debug('transform (ia|jb) with_df')
eris.ovon = mp._scf.with_df.ao2mo((co, cv, co, mo_coeff)).reshape(nocc, nvir, nocc, nmo)
eris.ovnv = mp._scf.with_df.ao2mo((co, cv, mo_coeff, cv)).reshape(nocc, nvir, nmo, nvir)
else:
raise NotImplementedError
log.debug('transform (ia|jb) outcore')
eris.feri = lib.H5TmpFile()
# ao2mo.outcore.general(mp.mol, (co,cv,co,cv), eris.feri,
# max_memory=max_memory, verbose=log)
# eris.ovov = eris.feri['eri_mo']
eris.ovov = _ao2mo_ovov(mp, co, cv, eris.feri, max(2000, max_memory), log)
time1 = log.timer('Integral transformation', *time0)
return eris
[docs]
class QSMP2(MP2):
# DIIS parameters for QSMP2
diis_space = 10
diis_start_cycle = 1
diis_file = None
def __init__(self, mf, frozen=None, ncore=None, ecut=None):
MP2.__init__(self, mf, frozen=frozen)
self.ncore = ncore
self.ecut = ecut
self.max_cycle = 20
self.conv_tol = 1e-5
self.eta = 5e-3
self.mo_energy = None
self.sigma = None
[docs]
def dump_flags(self):
log = logger.Logger(self.stdout, 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('QSMP2 nocc = %d, nvir = %d', nocc, nvir)
log.info('max_memory %d MB (current use %d MB)', self.max_memory, lib.current_memory()[0])
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 kernel(self, mo_energy=None, mo_coeff=None):
"""
Args:
mo_energy : 1D array (nmo), mean-field mo energy
mo_coeff : 2D array (nmo, nmo), mean-field mo coefficient
Returns:
self.mo_energy : QSMP2 quasiparticle energy
self.mo_coeff : QSMP2 MO coefficients (new orbitals)
"""
if mo_coeff is None:
mo_coeff = _mo_without_core(self, self._scf.mo_coeff)
if mo_energy is None:
mo_energy = _mo_energy_without_core(self, self._scf.mo_energy)
cput0 = (time.process_time(), time.perf_counter())
self.dump_flags()
self.converged, self.mo_energy, self.mo_coeff = kernel(self, mo_energy, mo_coeff)
logger.timer(self, 'QSMP2', *cput0)
return self.mo_energy, self.mo_coeff
[docs]
def ao2mo(self, mo_coeff=None):
return _make_eris(self, mo_coeff, verbose=self.verbose)
[docs]
def run_diis(self, ham, istep, adiis):
if adiis and istep >= self.diis_start_cycle:
ham = adiis.update(ham)
logger.debug1(self, 'DIIS for step %d', istep)
return ham
if __name__ == '__main__':
from pyscf import gto, dft, scf
mol = gto.Mole()
mol.verbose = 5
mol.atom = [[8, (0.0, 0.0, 0.0)], [1, (0.0, -0.7571, 0.5861)], [1, (0.0, 0.7571, 0.5861)]]
mol.basis = 'def2-svp'
mol.build()
mf = scf.RHF(mol)
mf.kernel()
nocc = mol.nelectron // 2
mp2 = QSMP2(mf, ncore=0, ecut=2.0)
mp2.kernel()