#!/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 MP2 quasiparticle energies
"""
from functools import reduce
import time, h5py
import numpy
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, orbs=None, kptlist=None, verbose=logger.NOTE):
"""MP2-corrected quasiparticle orbital energies
Parameters
----------
mp2 : KMP2QP object
mo_energy : double array
molecular orbital energies
mo_coeff : double ndarray
molecular orbital coefficients
orbs : list, optional
a list of orbital indices, by default range(nmo)
kptlist : list, optional
a list of k-point indices, by default range(nkpts)
Returns
-------
conv : bool
whether the quasiparticle equation is converged
mo_energy : double array
MP2-corrected quasiparticle orbital energies
"""
mf = mp2._scf
assert mp2.frozen is None
if orbs is None:
orbs = range(mp2.nmo)
if kptlist is None:
kptlist = range(mp2.nkpts)
nkpts = mp2.nkpts
nmo = mp2.nmo
nocc = mp2.nocc
nvir = nmo - nocc
eta = mp2.eta
mp2.mo_energy = np.array(mf.mo_energy)
mp2.mo_occ = np.array(mf.mo_occ)
mo_occ_1d = np.array(mp2.mo_occ).reshape(-1)
is_metal = False
if np.linalg.norm(np.abs(mo_occ_1d - 1.0) - 1.0) > 1e-5:
is_metal = True
conv = True
mo_energy = np.zeros_like(np.array(mf.mo_energy))
for kp in kptlist:
for p in orbs:
if mp2.linearized:
# linearized QP-MP2
de = 1e-6
ep = mf.mo_energy[kp][p]
if is_metal:
sigma = get_sigma_mp2_metal(mp2, kp, p, ep + 1j * eta).real
dsigma = get_sigma_mp2_metal(mp2, kp, p, ep + de + 1j * eta).real - sigma
else:
sigma = get_sigma_mp2(mp2, kp, p, ep + 1j * eta).real
dsigma = get_sigma_mp2(mp2, kp, p, ep + de + 1j * eta).real - sigma
zn = 1.0 / (1.0 - dsigma / de)
e = ep + zn * sigma
mo_energy[kp, p] = e
else:
# self-consistently solve QP equation
def quasiparticle(omega):
if is_metal:
sigma = get_sigma_mp2_metal(mp2, kp, p, omega + 1j * eta).real
else:
sigma = get_sigma_mp2(mp2, kp, p, omega + 1j * eta).real
return omega - mf.mo_energy[kp][p] - sigma
try:
e = newton(quasiparticle, mf.mo_energy[kp][p], tol=1e-5, maxiter=100)
mo_energy[kp, p] = e
except RuntimeError:
conv = False
if mp2.verbose >= logger.DEBUG:
numpy.set_printoptions(threshold=nmo)
for k in range(nkpts):
logger.debug(mp2, ' MP2 mo_energy @ k%d =\n%s', k, mo_energy[k])
numpy.set_printoptions(threshold=1000)
return conv, mo_energy
[docs]
def get_sigma_mp2(mp2, kp, p, ep):
"""
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)
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)
)
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
# apply regularization
if mp2.regularized:
t2_ipab = t2_ipab * (1 - np.exp(-np.abs(eipab) * mp2.kappa)) ** mp2.norm
wonvv = (2 * onvv_ip[ka][:, 0, :, :] - onvv_ip[kb][:, 0, :, :].transpose(0, 2, 1)).conj()
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)
)
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
# apply regularization
if mp2.regularized:
t2_ijap = t2_ijap * (1 - np.exp(-np.abs(eijap) * mp2.kappa)) ** mp2.norm
woovn = (2 * oovn_ap[ki][:, :, :, 0] - oovn_ap[kj][:, :, :, 0].transpose(1, 0, 2)).conj()
sigma -= einsum('ija,ija', t2_ijap, woovn) / nkpts
sigma = sigma / nkpts
return sigma
[docs]
class KMP2QP(KMP2):
def __init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None):
KMP2.__init__(self, mf, frozen=None, mo_coeff=None, mo_occ=None)
self.eta = 5e-3
self.linearized = False
self.mo_energy = None
self.sigma = None
self.kappa = 1.5
self.regularized = False
self.norm = 1
@property
def nocc(self):
frozen_mask = get_frozen_mask(self)
nkpts = len(self._scf.mo_energy)
nelec = 0.0
for k in range(nkpts):
nelec += np.sum(self._scf.mo_occ[k][frozen_mask[k]])
nelec = int(nelec / nkpts)
return nelec // 2
@nocc.setter
def nocc(self, n):
self._nocc = n
@property
def nmo(self):
frozen_mask = get_frozen_mask(self)
return len(self._scf.mo_energy[0][frozen_mask[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('MP2 nocc = %d, nvir = %d, nkpts = %d', nocc, nvir, nkpts)
if self.frozen is not None:
log.info('frozen orbitals = %s', str(self.frozen))
logger.info(self, 'use perturbative linearized QP eqn = %s', self.linearized)
return self
[docs]
def kernel(self, mo_energy=None, mo_coeff=None, orbs=None, kptlist=None):
"""
Args:
mo_energy : 2D array (nkpts, nmo), mean-field mo energy
mo_coeff : 3D array (nkpts, nmo, nmo), mean-field mo coefficient
orbs: list, orbital indices
kptlist: list, k-point indices
Returns:
self.mo_energy : MP2 quasiparticle energy
"""
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 = kernel(self, mo_energy, mo_coeff, orbs=orbs, kptlist=kptlist)
logger.timer(self, 'KMP2QP', *cput0)
return self.mo_energy
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=4,
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)
kmf.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)
kmf.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 = KMP2QP(kmf)
mp2.linearized = False
mp2.kernel(orbs=range(nocc - 2, nocc + 2))
print(mp2.mo_energy)
assert abs(mp2.mo_energy[0][nocc - 1] - 0.47002092) < 1e-5
assert abs(mp2.mo_energy[0][nocc] - 0.88134227) < 1e-5
# Test on Na (metallic)
cell = gto.Cell()
cell.build(
unit='angstrom',
a="""
-2.11250000000000 2.11250000000000 2.11250000000000
2.11250000000000 -2.11250000000000 2.11250000000000
2.11250000000000 2.11250000000000 -2.11250000000000
""",
atom="""
Na 0.00000 0.00000 0.00000
""",
dimension=3,
max_memory=8000,
verbose=5,
pseudo='gth-pade',
basis='gth-dzvp-molopt-sr',
precision=1e-10,
)
kpts = cell.make_kpts([2, 2, 1], scaled_center=[0, 0, 0])
gdf = df.GDF(cell, kpts)
gdf_fname = 'gdf_ints_221.h5'
gdf._cderi_to_save = gdf_fname
if not os.path.isfile(gdf_fname):
gdf.build()
chkfname = 'na_221.chk'
if os.path.isfile(chkfname):
kmf = scf.KRHF(cell, kpts)
kmf = scf.addons.smearing_(kmf, sigma=5e-3, method='fermi')
kmf.exxdiv = None
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)
kmf = scf.addons.smearing_(kmf, sigma=5e-3, method='fermi')
kmf.exxdiv = None
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
kmf.conv_tol = 1e-12
kmf.chkfile = chkfname
kmf.kernel()
mp2 = KMP2QP(kmf)
mp2.linearized = False
mp2.kernel(kptlist=[0, 1, 2], orbs=range(3, 7))
print(mp2.mo_energy)
assert abs(mp2.mo_energy[0][4] - 0.09788033) < 1e-5
assert abs(mp2.mo_energy[1][4] - 0.22536046) < 1e-5