#!/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>
#
"""
k-point spin-restricted quasiparticle self-consistent MP2
NOTE: this code has bugs in get_sigma_mp2 (off-diag self-energy)
"""
from functools import reduce
import time, h5py
import numpy, scipy
import numpy as np
from scipy.optimize import newton, least_squares
from pyscf import lib
from pyscf.lib import logger
from pyscf.ao2mo import _ao2mo
from pyscf.ao2mo.incore import _conc_mos
from pyscf.pbc import df, dft, scf
from pyscf.pbc.mp.kmp2 import KMP2, get_nocc, get_nmo, get_frozen_mask
from pyscf import __config__
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
nkpts = mp2.nkpts
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
mp2.mo_coeff = mp2._scf.mo_coeff
# cutoff for high virtuals
if mp2.ncut is None:
ncut = nmo
else:
ncut = mp2.ncut + nocc
if ncut >= nmo:
ncut = nmo
# get hcore and ovlp
hcore = np.array(mf.get_hcore())
ovlp = np.array(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
homo = -99.0
lumo = 99.0
for k in range(nkpts):
if homo < mo_energy[k][nocc - 1]:
homo = mo_energy[k][nocc - 1]
if lumo > mo_energy[k][nocc]:
lumo = mo_energy[k][nocc]
conv_pts = [homo, lumo]
gf = get_g0(conv_pts, mo_energy, eta)
qsmp2_conv = False
cycle = 0
while not qsmp2_conv and cycle < max(1, mp2.max_cycle):
# density matrix from updated QSMP2 density
dm_iter = np.zeros((nkpts, nmo, nmo), dtype=np.complex128)
for k in range(nkpts):
dm_iter[k] = 2.0 * np.dot(mo_coeff[k][:, :nocc], mo_coeff[k][:, :nocc].T.conj())
# update veff
veff = np.array(mf.get_veff(dm_iter))
# compute QSMP2 self-energy
sigma = np.zeros((nkpts, nmo, nmo), dtype=np.complex128)
for kp in range(nkpts):
for p in range(ncore, ncut):
for q in range(ncore, ncut):
if p == q:
sigma[kp, p, p] = get_sigma_mp2(mp2, kp, p, p, mo_energy[kp][p] + 1j * eta)
# elif q > p:
# sigma[kp,p,q] = get_sigma_mp2(mp2, kp, p, q, mo_energy[kp][p]+1j*eta,
# eq=mo_energy[kp][q]+1j*eta)
# sigma[kp,q,p] = sigma[kp,p,q]
# obtain vsig in AO basis
vsig = np.zeros((nkpts, nmo, nmo), dtype=np.complex128)
for k in range(nkpts):
CS = np.dot(mo_coeff[k].T.conj(), ovlp[k])
vsig[k] = reduce(np.dot, (CS.T.conj(), sigma[k].real, CS))
# complete Hamiltonian through DIIS
ham = veff + vsig
ham = mp2.run_diis(ham, cycle, mp2_diis)
ham = hcore + ham
# diagonalize
for k in range(nkpts):
mo_energy[k], mo_coeff[k] = scipy.linalg.eigh(ham[k], ovlp[k])
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 / nkpts / 2.0
if norm_dgf < mp2.conv_tol:
qsmp2_conv = True
logger.info(mp2, 'QSMP2 cycle= %d |dgf|= %4.3g', cycle + 1, norm_dgf)
homo = -99.0
lumo = 99.0
for k in range(nkpts):
if homo < mo_energy[k][nocc - 1]:
homo = mo_energy[k][nocc - 1]
if lumo > mo_energy[k][nocc]:
lumo = mo_energy[k][nocc]
logger.info(mp2, ' HOMO = %.2f eV, LUMO = %.2f eV', homo * 27.211386, lumo * 27.211386)
if mp2.verbose >= logger.DEBUG:
numpy.set_printoptions(threshold=nmo)
for k in range(nkpts):
logger.debug(mp2, ' QSMP2 mo_energy @ k%d =\n%s', k, mo_energy[k])
numpy.set_printoptions(threshold=1000)
cycle += 1
return qsmp2_conv, mo_energy, mo_coeff
[docs]
def get_sigma_mp2(mp2, kp, p, q, ep, eq=None):
"""
Compute KMP2 self-energy in MO basis on real axis
"""
mo_energy = mp2.mo_energy
mo_coeff = mp2.mo_coeff
nkpts = mp2.nkpts
nmo = mp2.nmo
nocc = mp2.nocc
nvir = nmo - nocc
kconserv = mp2.khelper.kconserv
fao2mo = mp2._scf.with_df.ao2mo
onvv_ip = np.zeros((nkpts, nocc, 1, nvir, nvir), dtype=mo_coeff[0].dtype)
oovn_ap = np.zeros((nkpts, nocc, nocc, nvir, 1), dtype=mo_coeff[0].dtype)
onvv_iq = np.zeros((nkpts, nocc, 1, nvir, nvir), dtype=mo_coeff[0].dtype)
oovn_aq = np.zeros((nkpts, nocc, nocc, nvir, 1), dtype=mo_coeff[0].dtype)
sigma = 0j
for ki in range(nkpts):
for ka in range(nkpts):
kb = kconserv[ki, ka, kp]
orbo_i = mo_coeff[ki][:, :nocc]
orbn_p = mo_coeff[kp][:, p].reshape(-1, 1)
orbv_a = mo_coeff[ka][:, nocc:]
orbv_b = mo_coeff[kb][:, nocc:]
onvv_ip[ka] = (
fao2mo(
(orbo_i, orbv_a, orbn_p, orbv_b),
(mp2.kpts[ki], mp2.kpts[ka], mp2.kpts[kp], mp2.kpts[kb]),
compact=False,
)
.reshape(nocc, nvir, 1, nvir)
.transpose(0, 2, 1, 3)
)
orbn_q = mo_coeff[kp][:, q].reshape(-1, 1)
onvv_iq[ka] = (
fao2mo(
(orbo_i, orbv_a, orbn_q, orbv_b),
(mp2.kpts[ki], mp2.kpts[ka], mp2.kpts[kp], mp2.kpts[kb]),
compact=False,
)
.reshape(nocc, nvir, 1, nvir)
.transpose(0, 2, 1, 3)
)
for ka in range(nkpts):
kb = kconserv[ki, ka, kp]
eia = mo_energy[ki][:nocc, None] - mo_energy[ka][None, nocc:]
epb = ep - mo_energy[kb][nocc:]
eipab = lib.direct_sum('ia,b->iab', eia, epb)
t2_ipab = onvv_ip[ka][:, 0, :, :] / eipab
wonvv = (2 * onvv_iq[ka][:, 0, :, :] - onvv_iq[kb][:, 0, :, :].transpose(0, 2, 1)).conj()
sigma += einsum('iab,iab', t2_ipab, wonvv) / nkpts
if p != q:
epb = eq - mo_energy[kb][nocc:]
eipab = lib.direct_sum('ia,b->iab', eia, epb)
t2_ipab = onvv_ip[ka][:, 0, :, :] / eipab
sigma += einsum('iab,iab', t2_ipab, wonvv) / nkpts
for ka in range(nkpts):
for ki in range(nkpts):
kj = kconserv[ka, ki, kp]
orbo_i = mo_coeff[ki][:, :nocc]
orbn_p = mo_coeff[kp][:, p].reshape(-1, 1)
orbv_a = mo_coeff[ka][:, nocc:]
orbo_j = mo_coeff[kj][:, :nocc]
oovn_ap[ki] = (
fao2mo(
(orbo_i, orbv_a, orbo_j, orbn_p),
(mp2.kpts[ki], mp2.kpts[ka], mp2.kpts[kj], mp2.kpts[kp]),
compact=False,
)
.reshape(nocc, nvir, nocc, 1)
.transpose(0, 2, 1, 3)
)
orbn_q = mo_coeff[kp][:, q].reshape(-1, 1)
oovn_aq[ki] = (
fao2mo(
(orbo_i, orbv_a, orbo_j, orbn_q),
(mp2.kpts[ki], mp2.kpts[ka], mp2.kpts[kj], mp2.kpts[kp]),
compact=False,
)
.reshape(nocc, nvir, nocc, 1)
.transpose(0, 2, 1, 3)
)
for ki in range(nkpts):
kj = kconserv[ka, ki, kp]
eia = mo_energy[ki][:nocc, None] - mo_energy[ka][None, nocc:]
ejp = mo_energy[kj][:nocc] - ep
eijap = lib.direct_sum('ia,j->ija', eia, ejp)
t2_ijap = oovn_ap[ki][:, :, :, 0] / eijap
woovn = (2 * oovn_aq[ki][:, :, :, 0] - oovn_aq[kj][:, :, :, 0].transpose(1, 0, 2)).conj()
sigma -= einsum('ija,ija', t2_ijap, woovn) / nkpts
if p != q:
ejp = mo_energy[kj][:nocc] - eq
eijap = lib.direct_sum('ia,j->ija', eia, ejp)
t2_ijap = oovn_ap[ki][:, :, :, 0] / eijap
sigma -= einsum('ija,ija', t2_ijap, woovn) / nkpts
if p != q:
sigma = sigma / 2.0
sigma = sigma / nkpts
return sigma
[docs]
def get_g0(omega, mo_energy, eta):
nkpts = len(mo_energy)
nmo = len(mo_energy[0])
nw = len(omega)
gf0 = np.zeros((nkpts, nmo, nmo, nw), dtype=np.complex128)
for iw in range(nw):
for k in range(nkpts):
gf0[k, :, :, iw] = np.diag(1.0 / (omega[iw] + 1j * eta - mo_energy[k]))
return gf0
[docs]
class KQSMP2(KMP2):
# DIIS parameters for QSGW
diis_space = 10
diis_start_cycle = 1
diis_file = None
def __init__(self, mf, ncore=None, ncut=None, frozen=None, mo_coeff=None, mo_occ=None):
KMP2.__init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None)
self.ncore = ncore
self.ncut = ncut
self.max_cycle = 20
self.conv_tol = 1e-5
self.eta = 5e-3
self.mo_energy = None
self.mo_coeff = None
self.sigma = None
@property
def nocc(self):
return self.mol.nelectron // 2
@nocc.setter
def nocc(self, n):
self._nocc = n
@property
def nmo(self):
return len(self._scf.mo_energy[0])
@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 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
nkpts = self.nkpts
log.info('QSMP2 nocc = %d, nvir = %d, nkpts = %d', nocc, nvir, nkpts)
return self
[docs]
def kernel(self, mo_energy=None, mo_coeff=None):
"""
Args:
mo_energy : 2D array (nkpts, nmo), mean-field mo energy
mo_coeff : 3D array (nkpts, nmo, nmo), mean-field mo coefficient
Returns:
self.mo_energy : QSMP2 quasiparticle energy
self.mo_coeff : QSMP2 mo_coeff
"""
if mo_coeff is None:
mo_coeff = np.array(self._scf.mo_coeff)
if mo_energy is None:
mo_energy = np.array(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, 'KQSMP2', *cput0)
return self.mo_energy, self.mo_coeff
[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__':
import os
from pyscf.pbc import gto, scf, mp
from pyscf.pbc.lib import chkfile
cell = gto.Cell()
cell.build(
unit='angstrom',
a="""
0.000000 1.783500 1.783500
1.783500 0.000000 1.783500
1.783500 1.783500 0.000000
""",
atom='C 1.337625 1.337625 1.337625; C 2.229375 2.229375 2.229375',
dimension=3,
max_memory=8000,
verbose=5,
pseudo='gth-pade',
basis='gth-dzv',
precision=1e-10,
)
kpts = cell.make_kpts([3, 1, 1], scaled_center=[0, 0, 0])
gdf = df.GDF(cell, kpts)
gdf_fname = 'gdf_ints_311.h5'
gdf._cderi_to_save = gdf_fname
if not os.path.isfile(gdf_fname):
gdf.build()
chkfname = 'diamond_311.chk'
if os.path.isfile(chkfname):
kmf = scf.KRHF(cell, kpts, exxdiv='ewald')
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
data = chkfile.load(chkfname, 'scf')
kmf.__dict__.update(data)
else:
kmf = scf.KRHF(cell, kpts, exxdiv='ewald')
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
kmf.conv_tol = 1e-12
kmf.chkfile = chkfname
kmf.kernel()
nocc = cell.nelectron // 2
mp2 = KQSMP2(kmf, ncut=4)
mp2.kernel()
print(mp2.mo_energy)