#!/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>
#
"""
PBC gamma-point spin-restricted G0W0 Greens function
This implementation has N^4 scaling, and is faster than GW-CD (N^4)
and analytic GW (N^6) methods.
GW-AC is recommended for valence states only, and is inaccurate for core states.
Method:
T. Zhu and G.K.-L. Chan, J. Chem. Theory. Comput. 17, 727-741 (2021)
Compute self-energy on imaginary frequency with density fitting,
then analytically continued to real frequency
Other useful references:
J. Chem. Theory Comput. 12, 3623-3635 (2016)
New J. Phys. 14, 053020 (2012)
"""
from functools import reduce
import numpy as np
from pyscf.ao2mo._ao2mo import nr_e2
from pyscf.lib import current_memory, temporary_env
from pyscf.pbc import df, scf
from pyscf.pbc.df.fft_ao2mo import _format_kpts
from pyscf.pbc.df.df_ao2mo import warn_pbc2d_eri
from fcdmft.gw.mol.gw_ac import GWAC as GWAC_mol
[docs]
class GWAC(GWAC_mol):
def __init__(self, mf, frozen=None, auxbasis=None):
if abs(mf.kpt).max() > 1e-9:
raise NotImplementedError
warn_pbc2d_eri(mf)
GWAC_mol.__init__(self, mf, frozen=frozen, auxbasis=auxbasis)
return
[docs]
def initialize_df(self, auxbasis=None):
"""Initialize density fitting.
Parameters
----------
auxbasis : str, optional
name of auxiliary basis set, by default None
"""
if getattr(self._scf, 'with_df', None):
self.with_df = self._scf.with_df
else:
self.with_df = df.DF(self._scf.mol)
if auxbasis is not None:
self.with_df.auxbasis = auxbasis
else:
try:
self.with_df.auxbasis = df.make_auxbasis(self._scf.mol, mp2fit=True)
except RuntimeError:
self.with_df.auxbasis = df.make_auxbasis(self._scf.mol, mp2fit=False)
self._keys.update(['with_df'])
return
[docs]
def ao2mo(self, mo_coeff=None):
"""Transform density-fitting integral from AO to MO.
Parameters
----------
mo_coeff : double 2d array, optional
coefficient from AO to MO, by default None
Returns
-------
Lpq : double 3d array
three-center density-fitting matrix in MO
"""
if mo_coeff is None:
mo_coeff = self.mo_coeff
nmo = mo_coeff.shape[1]
nao = self.mo_coeff.shape[0]
naux = self.with_df.get_naoaux()
kpts = self._scf.with_df.kpts
max_memory = max(2000, self._scf.max_memory - current_memory()[0] - nao**2 * naux * 8 / 1e6)
mo = np.asarray(mo_coeff, order='F')
ijslice = (0, nmo, 0, nmo)
kptijkl = _format_kpts(kpts)
eri_3d = []
for LpqR, _, _ in self._scf.with_df.sr_loop(kptijkl[:2], max_memory=0.3 * max_memory, compact=False):
Lpq = None
Lpq = nr_e2(LpqR.reshape(-1, nao, nao), mo, ijslice, aosym='s1', mosym='s1', out=Lpq)
eri_3d.append(Lpq)
eri_3d = np.vstack(eri_3d).reshape(-1, nmo, nmo)
return eri_3d
[docs]
def loop_ao2mo(self, mo_coeff=None, ijslice=None):
"""Transform density-fitting integral from AO to MO by block.
Parameters
----------
mo_coeff : double 2d array, optional
coefficient from AO to MO, by default None
ijslice : tuple, optional
tuples for (1st idx start, 1st idx end, 2nd idx start, 2nd idx end), by default None
Returns
-------
eri_3d : double 3d array
three-center density-fitting matrix in MO in a block
"""
if mo_coeff is None:
mo_coeff = self.mo_coeff
nmo = mo_coeff.shape[1]
nao = self.mo_coeff.shape[0]
naux = self.with_df.get_naoaux()
kpts = self._scf.with_df.kpts
max_memory = max(2000, self._scf.max_memory - current_memory()[0] - nao**2 * naux * 8 / 1e6)
mo = np.asarray(mo_coeff, order='F')
if ijslice is None:
ijslice = (0, nmo, 0, nmo)
nislice = ijslice[1] - ijslice[0]
njslice = ijslice[3] - ijslice[2]
kptijkl = _format_kpts(kpts)
eri_3d = []
for LpqR, _, _ in self._scf.with_df.sr_loop(kptijkl[:2], max_memory=0.2 * max_memory, compact=False):
Lpq = None
Lpq = nr_e2(LpqR.reshape(-1, nao, nao), mo, ijslice, aosym='s1', mosym='s1', out=Lpq)
eri_3d.append(Lpq)
eri_3d = np.vstack(eri_3d).reshape(-1, nislice, njslice)
return eri_3d
[docs]
def get_sigma_exchange(self, mo_coeff):
"""Get exchange self-energy (EXX).
Parameters
----------
mo_coeff : double 2d array
orbital coefficient
Returns
-------
vk : double 2d array
exchange self-energy
"""
dm = self._scf.make_rdm1()
if isinstance(self._scf.with_df, df.GDF):
rhf = scf.RHF(self.mol).density_fit()
elif isinstance(self._scf.with_df, df.RSDF):
rhf = scf.RHF(self.mol).rs_density_fit()
if hasattr(self._scf, 'sigma'):
rhf = scf.addons.smearing_(rhf, sigma=self._scf.sigma, method=self._scf.smearing_method)
rhf.exxdiv = None
rhf.with_df = self.with_df
vk = rhf.get_veff(self.mol, dm) - rhf.get_j(self.mol, dm)
vk = reduce(np.matmul, (mo_coeff.T, vk, mo_coeff))
return vk
if __name__ == '__main__':
import os
from pyscf.pbc import gto, scf, df, dft, tools
from pyscf.pbc.lib import chkfile
ucell = gto.Cell()
ucell.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 = 64000,
verbose = 5,
pseudo = 'gth-pade',
basis='gth-szv',
precision=1e-12)
kmesh = [3, 1, 1]
cell = tools.super_cell(ucell, kmesh)
cell.verbose = 5
gdf = df.RSDF(cell)
gdf_fname = 'diamond_sc_gdf_ints.h5'
gdf._cderi_to_save = gdf_fname
if not os.path.isfile(gdf_fname):
gdf.build()
chkfname = 'diamond_sc.chk'
if os.path.isfile(chkfname):
kmf = dft.RKS(cell).rs_density_fit()
kmf.xc = 'pbe'
kmf.exxdiv = None
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
kmf.conv_tol = 1e-12
data = chkfile.load(chkfname, 'scf')
kmf.__dict__.update(data)
else:
kmf = dft.RKS(cell).rs_density_fit()
kmf.xc = 'pbe'
kmf.exxdiv = None
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
kmf.conv_tol = 1e-12
kmf.chkfile = chkfname
kmf.kernel()
gw = GWAC(kmf)
omega = np.linspace(0.3, 1.3, 101)
gw.kernel()
assert abs(gw.mo_energy[5] - 0.52637379) < 1e-4
assert abs(gw.mo_energy[10] - 0.62044176) < 1e-4
assert abs(gw.mo_energy[12] - 0.96572544) < 1e-4
assert abs(gw.mo_energy[15] - 1.0751724) < 1e-4