#!/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>
#
"""
Periodic spin-restricted random phase approximation
(direct RPA/dRPA in chemistry) with N^4 scaling (Gamma only)
Method:
Main routines are based on GW-AC method described in:
T. Zhu and G.K.-L. Chan, J. Chem. Theory. Comput. 17, 727-741 (2021)
X. Ren et al., New J. Phys. 14, 053020 (2012)
"""
import numpy as np
import os
from pyscf.ao2mo import _ao2mo
from pyscf.lib import current_memory, temporary_env
from pyscf.pbc import df, dft, scf
from fcdmft.rpa.mol.rpa import RPA, get_idx_metal, _mo_energy_without_core
[docs]
class RPA(RPA):
[docs]
def get_ehf(self):
"""Get Hartree-Fock energy.
Returns
-------
e_hf : double
Hartree-Fock energy
"""
with temporary_env(self.with_df, verbose=0), temporary_env(self.mol, verbose=0):
dm = self._scf.make_rdm1()
rhf = scf.RHF(self.mol, exxdiv=self._scf.exxdiv)
rhf.verbose = 0
if hasattr(self._scf, 'sigma'):
rhf = scf.addons.smearing_(rhf, sigma=self._scf.sigma, method='fermi')
rhf.with_df = self._scf.with_df
e_hf = rhf.energy_elec(dm=dm)[0]
e_hf += self._scf.energy_nuc()
return e_hf
[docs]
def ao2mo(self, mo_coeff=None, fullmo=False, mo_occ=None):
"""Transform density-fitting integral from AO to MO.
Parameters
----------
mo_coeff : double 1d array, optional
orbital energy, by default None
fullmo : bool, optional
transform AO to MO of all orbitals, by default False
mo_occ : double 1d array, optional
occupation number, by default None
Returns
-------
eri_3d : 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)
is_metal = False
if mo_occ is None:
mo_occ_1d = _mo_energy_without_core(self, self._scf.mo_occ).reshape(-1)
else:
mo_occ_1d = np.array(mo_occ).reshape(-1)
if np.linalg.norm(np.abs(mo_occ_1d - 1.0) - 1.0) > 1e-5:
is_metal = True
mo = np.asarray(mo_coeff, order='F')
if fullmo:
ijslice = (0, nmo, 0, nmo)
elif is_metal:
idx_occ, idx_frac, idx_vir = get_idx_metal(mo_occ_1d)
nocc = len(idx_occ)
nfrac = len(idx_frac)
ijslice = (0, nocc + nfrac, nocc, nmo)
else:
nocc = self.nocc
ijslice = (0, nocc, nocc, nmo)
nislice = ijslice[1] - ijslice[0]
njslice = ijslice[3] - ijslice[2]
print('ijslice', ijslice, flush=True)
from pyscf.pbc.df.fft_ao2mo import _format_kpts
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 = _ao2mo.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 loop_ao2mo(self, mo_coeff=None, fullmo=False, ijslice=None):
if mo_coeff is None:
mo_coeff = self.mo_coeff
nmo = mo_coeff.shape[1]
nao = self.mo_coeff.shape[0]
self.with_df.get_naoaux()
kpts = self._scf.with_df.kpts
is_metal = False
mo_occ_1d = np.array(self._scf.mo_occ).reshape(-1)
if np.linalg.norm(np.abs(mo_occ_1d - 1.0) - 1.0) > 1e-5:
is_metal = True
mo = np.asarray(mo_coeff, order='F')
if ijslice is None:
if fullmo:
ijslice = (0, nmo, 0, nmo)
elif is_metal:
# Assume the ordering of mo_coeff is
# occ, frac, vir
idx_occ, idx_frac, idx_vir = get_idx_metal(mo_occ_1d)
for i in idx_occ:
assert i < idx_frac[0]
for i in idx_frac:
assert i < idx_vir[0]
nocc = len(idx_occ)
nfrac = len(idx_frac)
nvir = len(idx_vir)
ijslice = (0, nocc + nfrac, nocc, nmo)
else:
nocc = self.nocc
ijslice = (0, nocc, nocc, nmo)
nislice = ijslice[1] - ijslice[0]
njslice = ijslice[3] - ijslice[2]
from pyscf.pbc.df.fft_ao2mo import _format_kpts
kptijkl = _format_kpts(kpts)
eri_3d = []
for LpqR, _, sign in self._scf.with_df.sr_loop(
kptijkl[:2], max_memory=0.1 * self._scf.max_memory, compact=False
):
Lpq = None
Lpq = _ao2mo.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
if __name__ == '__main__':
from pyscf.pbc import gto, scf, dft, df, tools
from pyscf.pbc.lib import chkfile
# test on diamond
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=16000,
verbose=5,
pseudo='gth-pbe',
basis='gth-dzv',
precision=1e-12,
)
kmesh = [3, 1, 1]
cell = tools.super_cell(ucell, kmesh)
cell.verbose = 5
gdf = df.RSDF(cell)
gdf_fname = 'gdf_ints.h5'
gdf._cderi_to_save = gdf_fname
if not os.path.isfile(gdf_fname):
gdf.build()
chkfname = 'diamond_hf.chk'
if os.path.isfile(chkfname):
kmf = scf.RHF(cell).density_fit()
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 = scf.RHF(cell).density_fit()
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
kmf.conv_tol = 1e-12
kmf.chkfile = chkfname
kmf.kernel()
rpa = RPA(kmf)
rpa.kernel()
assert abs(rpa.e_corr - -0.5558316165999143) < 1e-6
assert abs(rpa.e_tot - -32.08317615664809) < 1e-6
# Test on Na (metallic)
ucell = gto.Cell()
ucell.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=16000,
verbose=5,
pseudo='gth-pade',
basis='gth-dzvp-molopt-sr',
precision=1e-10,
)
kmesh = [2, 2, 1]
cell = tools.super_cell(ucell, kmesh)
cell.verbose = 5
gdf = df.RSDF(cell)
gdf_fname = 'gdf_ints_na.h5'
gdf._cderi_to_save = gdf_fname
if not os.path.isfile(gdf_fname):
gdf.build()
chkfname = 'na_lda.chk'
if os.path.isfile(chkfname):
kmf = dft.RKS(cell).density_fit()
kmf = scf.addons.smearing_(kmf, sigma=5e-3, method='fermi')
kmf.xc = 'lda'
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).density_fit()
kmf = scf.addons.smearing_(kmf, sigma=5e-3, method='fermi')
kmf.xc = 'lda'
kmf.with_df = gdf
kmf.with_df._cderi = gdf_fname
kmf.conv_tol = 1e-12
kmf.chkfile = chkfname
kmf.kernel()
rpa = RPA(kmf)
rpa.kernel()
assert abs(rpa.e_corr - -0.1612718859071952) < 1e-6 # old value: -0.1612721103400592
assert abs(rpa.e_tot - -190.42773568569532) < 1e-6 # old value: -190.4277327357537
print('passed tests!')