import os
import h5py
import logging
import numpy as np
from pyscf.df import addons
from pyscf.pbc import tools, gto
from . import gdf_s_metric as gdf_S
from . import common_utils as comm
from . import integral_utils as int_utils
from ..pesto import ft
[docs]
class pyscf_init:
'''Initialization class for Green project
Attributes
----------
args : map
simulation parameters
cell : pyscf.pbc.cell
unit cell object
kmesh : numpy.ndarray
Monkhorst-Pack reciprocal space grid
'''
def __init__(self, args):
'''
Initialize PySCF interoperation class
Parameters
----------
args: map
simulation parameters
'''
self.args = args
if self.args.Nk is None:
self.args.Nk = 0
if self.args.spin is None:
self.args.spin = 0
if self.args.damping is None:
self.args.damping = 0
if self.args.max_iter is None:
self.args.max_iter = 100
self.cell = self.cell_object()
[docs]
def compute_df_int(self, nao, X_k):
raise NotImplementedError("Please Implement this method")
[docs]
def mf_object(self, mydf=None):
raise NotImplementedError("Please Implement this method")
[docs]
def df_object(self, mydf=None):
raise NotImplementedError("Please Implement this method")
[docs]
def cell_object(self):
raise NotImplementedError("Please Implement this method")
[docs]
class pyscf_pbc_init (pyscf_init):
'''Initialization class for periodic / solid-state systems for the Green project
'''
def __init__(self, args=None):
super().__init__(comm.init_pbc_params() if args is None else args)
self.kmesh, self.k_ibz, self.ir_list, self.conj_list, self.weight, self.ind, self.num_ik = comm.init_k_mesh(self.args, self.cell)
[docs]
def compute_df_int(self, nao, X_k):
'''
Generate density-fitting integrals for correlated methods
'''
mydf = comm.construct_gdf(self.args, self.cell, self.kmesh)
int_utils.compute_integrals(self.args, self.cell, mydf, self.kmesh, nao, X_k, self.args.hf_int_path, "cderi.h5", True, True)
mydf = None
if 'gf2' in self.args.finite_size_kind or 'gw' in self.args.finite_size_kind or 'gw_s' in self.args.finite_size_kind:
self.compute_twobody_finitesize_correction()
if not self.args.keep_cderi:
os.remove("cderi.h5")
os.system("sync")
return
mydf = comm.construct_gdf(self.args, self.cell, self.kmesh)
# Use Ewald for divergence treatment
mydf.exxdiv = 'ewald'
import importlib.util as iu
new_pyscf = iu.find_spec('pyscf.pbc.df.gdf_builder') is not None
if not new_pyscf :
from pyscf.pbc import df as gdf
weighted_coulG_old = gdf.GDF.weighted_coulG
from pyscf.pbc import df as gdf
import green_igen.df as gggdf
auxcell = gggdf.make_modrho_basis(mydf.cell, mydf.auxbasis,
mydf.exp_to_discard)
kptij_lst = [(ki, ki) for i, ki in enumerate(self.kmesh)]
kptij_lst = np.asarray(kptij_lst)
gdf.GDF.weighted_coulG = int_utils.weighted_coulG_ewald
gggdf._make_j3c(mydf, self.cell, auxcell, kptij_lst, "cderi_ewald.h5")
#weighted_coulG_old = gdf.GDF.weighted_coulG
gdf.GDF.weighted_coulG2 = None
#kij_conj, kij_trans, kpair_irre_list, kptij_idx, num_kpair_stored =
int_utils.compute_integrals(self.args, self.cell, mydf, self.kmesh, nao, X_k, self.args.int_path, "cderi.h5", True, self.args.keep_cderi, cderi_name2="cderi_ewald.h5")
if not new_pyscf :
gdf.GDF.weighted_coulG = weighted_coulG_old
[docs]
def evaluate_high_symmetry_path(self):
if self.args.print_high_symmetry_points:
comm.print_high_symmetry_points(self.args)
return
if self.args.high_symmetry_path is None:
raise RuntimeError("Please specify high-symmetry path")
if self.args.high_symmetry_path is not None:
try:
comm.check_high_symmetry_path(self.args)
except RuntimeError as e:
logging.error("\n\n\n")
logging.error("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
logging.error("!!!!!!!!! Cannot compute high-symmetry path !!!!!!!!!")
logging.error("!! Correct or Disable high-symmetry path evaluation !")
logging.error("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!")
logging.error(e)
exit(-1)
kmesh_hs, Hk_hs, Sk_hs, lin_kpt_axis = comm.high_symmetry_path(
self.cell, self.args
)
xpath, special_points, special_labels = lin_kpt_axis
inp_data = h5py.File(self.args.output_path, "a")
logging.debug(kmesh_hs)
logging.debug(self.cell.get_scaled_kpts(kmesh_hs))
inp_data["high_symm_path/k_mesh"] = self.cell.get_scaled_kpts(kmesh_hs)
inp_data["high_symm_path/r_mesh"] = ft.construct_rmesh(self.args.nk, self.args.nk, self.args.nk)
inp_data["high_symm_path/Hk"] = Hk_hs
inp_data["high_symm_path/Sk"] = Sk_hs
inp_data["high_symm_path/xpath"] = xpath
inp_data["high_symm_path/special_points"] = special_points
inp_data["high_symm_path/special_labels"] = special_labels
[docs]
def compute_twobody_finitesize_correction(self, mydf=None):
if not os.path.exists(self.args.hf_int_path):
os.mkdir(self.args.hf_int_path)
if 'gf2' in self.args.finite_size_kind :
comm.compute_ewald_correction(self.args, self.cell, self.kmesh, self.args.hf_int_path + "/df_ewald.h5")
if 'gw' in self.args.finite_size_kind :
self.evaluate_gw_correction(mydf)
[docs]
def evaluate_gw_correction(self, mydf=None):
if mydf is None:
mydf = comm.construct_gdf(self.args, self.cell, self.kmesh)
mydf.build()
j3c, kptij_lst, j2c_sqrt, uniq_kpts = gdf_S.make_j3c(mydf, self.cell, j2c_sqrt=True, exx=False)
''' Transformation matrix from auxiliary basis to plane-wave '''
AqQ, q_reduced, q_scaled_reduced = gdf_S.transformation_PW_to_auxbasis(mydf, self.cell, j2c_sqrt, uniq_kpts)
q_abs = np.array([np.linalg.norm(qq) for qq in q_reduced])
q_abs = np.array([round(qq, 8) for qq in q_abs])
# Different prefactors for the GW finite-size correction for testing
# In practice, the madelung constant is used, which decays as (1/nk).
X = (6*np.pi**2)/(self.cell.vol*len(self.kmesh))
X = (2.0/np.pi) * np.cbrt(X)
X2 = 2.0 * np.cbrt(1.0/(self.cell.vol*len(self.kmesh)))
f = h5py.File(self.args.hf_int_path + "/AqQ.h5", 'w')
f["AqQ"] = AqQ
f["qs"] = q_reduced
f["qs_scaled"] = q_scaled_reduced
f["q_abs"] = q_abs
f["X"] = X
f["X2"] = X2
f["madelung"] = tools.pbc.madelung(self.cell, self.kmesh)
f.close()
[docs]
def mf_object(self, mydf=None):
return comm.solve_mean_field(self.args, mydf, self.cell)
[docs]
def df_object(self, mydf=None):
return comm.construct_gdf(self.args, self.cell, self.kmesh)
[docs]
def cell_object(self):
return comm.pbc_cell(self.args)
[docs]
class pyscf_mol_init (pyscf_init):
'''Initialization class for molecular systems in the Green project
'''
def __init__(self, args=None):
super().__init__(comm.init_mol_params() if args is None else args)
self.kmesh = np.array([[0.,0.,0.]])
self.k_ibz = np.array([[0.,0.,0.],])
self.ir_list = np.array([0])
self.conj_list= np.array([0])
self.weight= np.array([1.0])
self.ind= np.array([0])
self.num_ik = 1
self.kcell = gto.Cell(verbose=0)
self.kcell.a = [[1,0,0],[0,1,0],[0,0,1]]
self.kcell.atom = self.cell.atom
self.kcell.spin = self.cell.spin
self.kcell.charge = self.cell.charge
self.kcell.unit = 'A'
self.kcell.basis = self.cell.basis
self.kcell.kpts = self.kcell.make_kpts([1, 1, 1])
self.kcell.ecp = self.cell.ecp
self.kcell.build()
[docs]
def compute_df_int(self, nao, X_k):
'''
Generate density-fitting integrals for correlated methods
'''
h_in = h5py.File("cderi_mol.h5", 'r')
h_out = h5py.File("cderi.h5", 'w')
j3c_obj = h_in["/j3c"]
if not isinstance(j3c_obj, h5py.Dataset): # not a dataset
if isinstance(j3c_obj, h5py.Group): # pyscf >= 2.1
h_in.copy(h_in["/j3c"], h_out, "j3c/0")
else:
raise ValueError("Unknown structure of cderi_mol.h5. Perhaps, PySCF upgrade went badly...")
else: # pyscf < 2.1
h_in.copy(h_in["/j3c"], h_out, "j3c/0/0")
kptij = np.zeros((1, 2, 3))
h_out["j3c-kptij"] = kptij
h_in.close()
h_out.close()
mydf = comm.construct_gdf(self.args, self.kcell, self.kmesh)
int_utils.compute_integrals(self.args, self.kcell, mydf, self.kmesh, nao, X_k, "df_hf_int", "cderi.h5", True, self.args.keep_cderi)
mydf = None
[docs]
def df_object(self, mydf=None):
return comm.construct_mol_gdf(self.args, self.kcell)
[docs]
def mf_object(self, mydf=None):
return comm.solve_mol_mean_field(self.args, mydf, self.cell)
[docs]
def cell_object(self):
return comm.mol_cell(self.args)