import argparse
import inspect
import logging
import os
import h5py
import numpy as np
import scipy.linalg as LA
import pyscf.lib.chkfile as chk
from io import StringIO
from numba import jit
from pyscf import gto as mgto
from pyscf import scf as mscf
from pyscf import dft as mdft
from pyscf import df as mdf
from pyscf.pbc import tools, gto, df, scf, dft
from pyscf import __version__ as pyscf_version
from pyscf.pbc.lib import kpts as libkpts
from . import integral_utils as int_utils
from . import kpt_utils
from . import ortho_utils
from .symmetry_utils import get_representation, get_spinor_representation
from ..version import __version__
def _init_guess_from_chk(mf, cell, chkfile):
"""Call mf.init_guess_by_chkfile compatibly across PySCF versions.
PySCF < ~2.1 signature: init_guess_by_chkfile(cell, chkfile_name, ...)
PySCF >= ~2.1 signature: init_guess_by_chkfile(chk=None, ...) (uses self.cell)
GHF signature: init_guess_by_chkfile(chkfile, project)
"""
first_param = next(iter(inspect.signature(mf.init_guess_by_chkfile).parameters))
if first_param in ("chk", "chkfile"):
return mf.init_guess_by_chkfile(chkfile)
else:
return mf.init_guess_by_chkfile(cell, chkfile)
[docs]
def print_high_symmetry_points(args):
"""For given simulation parameters, generate and print list of lattice special points
Parameters
----------
args : map
simulation parameters
"""
import ase.spacegroup
lattice_vectors, symbols, positions = extract_ase_data(args.a, args.atom)
cc = ase.spacegroup.crystal(symbols, positions, cellpar=ase.geometry.cell_to_cellpar(lattice_vectors))
lat = cc.cell.get_bravais_lattice()
special_points = lat.get_special_points()
print("List of special points: {}".format(special_points))
[docs]
def check_high_symmetry_path(args):
"""Check that selected high-symmetry path is correct for the chosen simulation parameters
Parameters
----------
args : map
simulation parameters
"""
if args.high_symmetry_path is None:
return
import ase.spacegroup
lattice_vectors, symbols, positions = extract_ase_data(args.a, args.atom)
logging.info(f"parse: {lattice_vectors} {symbols} {positions}")
cc = ase.spacegroup.crystal(symbols, positions, cellpar=ase.geometry.cell_to_cellpar(lattice_vectors))
lat = cc.cell.get_bravais_lattice()
special_points = lat.get_special_points()
path = args.high_symmetry_path
for sp in special_points.keys():
path = path.replace(sp, "")
path = path.replace(",", "")
if path != "":
raise RuntimeError(("Chosen high symmetry path {} has invalid special points {}. Valid "
"special points are {} ").format(args.high_symmetry_path, path, special_points.keys()))
[docs]
def high_symmetry_path(cell, args):
"""Compute high-symmetry k-path
Parameters
----------
cell : pyscf.pbc.Cell
unit-cell object
args : map
simulation parameters
Returns
-------
tuple
k-points on the chosen high-symmetry path; corresponding
non-interacting Hamiltonian and overlap matrix, and linear k-point
axis and labels (used in band structure plots)
"""
if args.high_symmetry_path is None:
return [None, None, None]
import ase
lattice_vectors, _, _ = extract_ase_data(args.a, args.atom)
kpath = ase.dft.kpoints.bandpath(
args.high_symmetry_path, lattice_vectors,
npoints=args.high_symmetry_path_points
)
kmesh = cell.get_abs_kpts(kpath.kpts)
if args.x2c == 2:
new_mf = args.mean_field(cell, kmesh).density_fit().x2c1e()
elif args.x2c == 1:
new_mf = args.mean_field(cell, kmesh).density_fit().sfx2c1e()
else:
new_mf = args.mean_field(cell, kmesh).density_fit()
H0_hs = new_mf.get_hcore()
Sk_hs = new_mf.get_ovlp()
# get symmetry labels
# lin_kpt_axis = a tuple of (
# linear axis for plotting kpoints,
# special point location along the linear axis,
# symmetry point labels
# )
lin_kpt_axis = kpath.get_linear_kpoint_axis()
return [kmesh, H0_hs, Sk_hs, lin_kpt_axis]
[docs]
def fold_back_to_1stBZ(kpts):
"""Map each k-point from a given list of scaled k-points into the first Brillouin zone
Parameters
----------
kpts : numpy.ndarray
list of k-points
Returns
-------
numpy.ndarray
k-points folded into 1st Brillouin-zone
"""
nkpts = len(kpts)
for i, ik in enumerate(kpts):
kpts[i] = np.array([wrap_1stBZ(kk) for kk in ik])
return kpts
[docs]
def inversion_sym(kmesh_scaled):
"""For a given list of the scaled k-points in the full Brillouin zone select k-points that are
equivalent by the time-reversal symmetry and return them and their corresponding weight and index
Parameters
----------
kpts : numpy.ndarray
list of scaled k-points
Returns
-------
numpy.ndarray
indices of irreducible k-points in the input k-list
numpy.ndarray
inverse index, associating each k-point with its unique equivalent
numpy.ndarray
weights for each irreducible k-point (degeneracy)
numpy.ndarray
truth table for whether a k-point has an irreducible time-reversal equivalent
"""
ind = np.arange(np.shape(kmesh_scaled)[0])
weight = np.zeros(np.shape(kmesh_scaled)[0])
for i, ki in enumerate(kmesh_scaled):
ki = [wrap_1stBZ(l) for l in ki]
kmesh_scaled[i] = ki
# Time-reversal symmetry
Inv = (-1) * np.identity(3)
for i, ki in enumerate(kmesh_scaled):
ki = np.dot(Inv,ki)
ki = [wrap_1stBZ(l) for l in ki]
for l, kl in enumerate(kmesh_scaled[:i]):
if np.allclose(ki,kl):
ind[i] = l
break
uniq = np.unique(ind, return_counts=True)
for i, k in enumerate(uniq[0]):
weight[k] = uniq[1][i]
ir_list = uniq[0]
# Mark down time-reversal-reduced k-points
conj_list = np.zeros(len(kmesh_scaled))
for i, k in enumerate(ind):
if i != k:
conj_list[i] = 1
return ir_list, ind, weight, conj_list
[docs]
def wrap_k(k):
while k < 0 :
k = 1 + k
while (k - 9.9999999999e-1) > 0.0 :
k = k - 1
return k
[docs]
def parse_basis(basis_list):
"""Parse information about chosen basis sets
Parameters
----------
basis_list : str
basis-set information
Returns
-------
list
basis-set information usable in initialization of pyscf Cell object
"""
logging.debug(f"{basis_list}, {len(basis_list) % 2}")
if len(basis_list) % 2 == 0:
b = {}
for atom_i in range(0, len(basis_list), 2):
bas_i = basis_list[atom_i + 1]
if os.path.exists(bas_i) :
with open(bas_i) as bfile:
bas = mgto.parse(bfile.read())
# if basis specified as a standard basis
else:
bas = bas_i
b[basis_list[atom_i]] = bas
return b
else:
return basis_list[0]
[docs]
def parse_geometry(g):
"""Parse geometry of the system
"""
res = ""
if os.path.exists(g) :
with open(g) as gf:
res = gf.read()
else:
res = g
return res
[docs]
def save_data(args, mycell, mf, kmesh, ind, weight, num_ik, ir_list, conj_list, Nk, nk, NQ, F, S, T, hf_dm, madelung, Zs, last_ao):
'''
Save data in Green/WeakCoupling format into a hdf5 file
'''
kptij_idx, kij_conj, kij_trans, kpair_irre_list, num_kpair_stored, kptis, kptjs = int_utils.integrals_grid(mycell, kmesh)
logging.info(f"number of reduced k-pairs: {num_kpair_stored}")
inp_data = h5py.File(args.output_path, "w")
inp_data["symmetry/k/mesh"] = kmesh
inp_data["symmetry/k/mesh_scaled"] = mycell.get_scaled_kpts(kmesh)
inp_data["symmetry/k/bz2ibz"] = ind
inp_data["symmetry/k/weight_ibz"] = weight
inp_data["symmetry/k/ink"] = num_ik
inp_data["symmetry/k/nk"] = nk
inp_data["symmetry/k/ibz2bz"] = ir_list
inp_data["symmetry/k/tr_conj"] = conj_list
# k-point pairs for integrals
inp_data["symmetry/pairs/conj_pairs_list"] = kij_conj
inp_data["symmetry/pairs/trans_pairs_list"] = kij_trans
inp_data["symmetry/pairs/kpair_irre_list"] = kpair_irre_list
inp_data["symmetry/pairs/kpair_idx"] = kptij_idx
inp_data["symmetry/pairs/num_kpair_stored"] = num_kpair_stored
inp_data["HF/Nk"] = Nk
inp_data["HF/nk"] = nk
inp_data["HF/Energy"] = mf.e_tot
inp_data["HF/Energy_nuc"] = mf.energy_nuc()
inp_data["HF/Fock-k"] = F.view(np.float64).reshape(F.shape[0], F.shape[1], F.shape[2], F.shape[3], 2)
inp_data["HF/Fock-k"].attrs["__complex__"] = np.int8(1)
inp_data["HF/S-k"] = S.view(np.float64).reshape(S.shape[0], S.shape[1], S.shape[2], S.shape[3], 2)
inp_data["HF/S-k"].attrs["__complex__"] = np.int8(1)
inp_data["HF/H-k"] = T.view(np.float64).reshape(T.shape[0], T.shape[1], T.shape[2], T.shape[3], 2)
inp_data["HF/H-k"].attrs["__complex__"] = np.int8(1)
inp_data["HF/madelung"] = madelung
inp_data["HF/mo_energy"] = mf.mo_energy
inp_data["HF/mo_coeff"] = mf.mo_coeff
inp_data["mulliken/Zs"] = Zs
inp_data["mulliken/last_ao"] = last_ao
inp_data["params/nao"] = mycell.nao_nr()
inp_data["params/nso"] = S.shape[2]
inp_data["params/ns"] = S.shape[0]
inp_data["params/nel_cell"] = mycell.nelectron
inp_data["params/nk"] = kmesh.shape[0]
nk_arr = np.atleast_1d(np.array(args.nk, dtype=int))
inp_data["symmetry/k/nk_list"] = np.array([nk_arr[0]]*3, dtype=int) if nk_arr.size == 1 else nk_arr
inp_data["params/NQ"] = NQ
inp_data.attrs["__green_version__"] = __version__
inp_data.close()
chk.save(args.output_path, "Cell", mycell.dumps())
inp_data = h5py.File(os.path.join(os.path.dirname(args.output_path),"dm.h5"), "w")
inp_data["HF/dm-k"] = hf_dm.view(np.float64).reshape(hf_dm.shape[0], hf_dm.shape[1], hf_dm.shape[2], hf_dm.shape[3], 2)
inp_data["HF/dm-k"].attrs["__complex__"] = np.int8(1)
inp_data["dm_gamma"] = hf_dm[:, 0, :, :]
inp_data.close()
[docs]
def orthogonalize(mydf, orth, X_k, X_inv_k, F, T, hf_dm, S, mf=None):
'''
Transform Fock-matrix, non-interacting Hamiltonian, density matrix and overlap matrix into an orthogonal basis.
``orth`` selects the per-k transformation:
- "none" : identity (AO basis preserved)
- "lowdin" : symmetric S^{-1/2} orthogonalization
- "mo" : canonical molecular orbitals from ``mf.mo_coeff``
- "natural" : natural orbitals from the mean-field density matrix
``mf`` is required for "mo"; "natural" uses ``hf_dm``. Per-k basis
construction is delegated to ``ortho_utils.{lowdin,mo,natural}_per_k``.
'''
if orth == "mo" and mf is None:
raise ValueError("orthogonalize: mf is required for orth='mo'.")
ns = hf_dm.shape[0]
if orth == "mo" and ns == 2:
# No single C diagonalizes both spin Fock blocks; we fall back to the
# spin-averaged Fock. The resulting MOs are not the eigenstates of
# either F_alpha or F_beta individually.
logging.warning(
"orthogonalize: orth='mo' with ns=2 (UHF/UKS); using "
"spin-averaged MOs (eigenstates of 0.5*(F_alpha+F_beta)), "
"not the canonical alpha/beta MOs."
)
maxdiff = -1
old_shape = [-1, -1]
for ik, k in enumerate(mydf.kpts):
if orth == "none":
X_inv_k.append(np.eye(F.shape[2], dtype=np.complex128))
X_k.append(np.eye(F.shape[2], dtype=np.complex128))
continue
Sk = S[0, ik]
if orth == "lowdin":
x, x_pinv = ortho_utils.lowdin_per_k(Sk)
elif orth == "symmetric_lowdin":
x, x_pinv = ortho_utils.symmetric_lowdin_per_k(Sk)
elif orth == "mo":
# For ns == 2, no single C diagonalizes both F_alpha and F_beta;
# diagonalize the spin-averaged Fock against S to obtain a
# spin-symmetric MO basis. For ns == 1 use mf.mo_coeff directly.
if ns == 2:
F_bar = 0.5 * (F[0, ik] + F[1, ik])
_, C_k = LA.eigh(F_bar, Sk)
else:
C_k = mf.mo_coeff[ik]
x, x_pinv = ortho_utils.mo_per_k(Sk, C_k)
elif orth == "natural":
dmk = 0.5 * (hf_dm[0, ik] + hf_dm[1, ik]) if ns == 2 else hf_dm[0, ik]
x, x_pinv = ortho_utils.natural_per_k(Sk, dmk)
else:
raise ValueError(f"orthogonalize: unknown orth '{orth}'.")
n_ortho, n_nonortho = x.shape
if old_shape[0] >= 0 and n_ortho != old_shape[0] and n_nonortho != old_shape[1]:
raise RuntimeError("Error!!! Different k-point have different number of orthogonal basis.")
old_shape[0] = n_ortho
old_shape[1] = n_nonortho
X_inv_k.append(x_pinv.copy())
X_k.append(x.copy())
diff = np.eye(n_nonortho) - np.dot(x, x_pinv)
diff_max = np.max(np.abs(diff))
maxdiff = max(maxdiff, diff_max)
logging.info(f"max diff from identity {maxdiff}")
if orth == "none":
X_inv_k = np.asarray(X_inv_k).reshape(F.shape[1:])
X_k = np.asarray(X_k).reshape(F.shape[1:])
return X_k, X_inv_k, S, F, T, hf_dm
X_inv_k = np.asarray(X_inv_k).reshape(F.shape[1:])
X_k = np.asarray(X_k).reshape(F.shape[1:])
F = transform(F, X_k, X_inv_k)
T = transform(T, X_k, X_inv_k)
hf_dm = transform(hf_dm, X_inv_k, X_k)
S = np.array([np.eye(F.shape[-1], dtype=np.complex128)] * F.shape[1])
S = np.array([S] * ns)
return X_k, X_inv_k, S, F, T, hf_dm
[docs]
def add_common_params(parser):
'''
Define common command line arguments for Green python module
'''
parser.add_argument("--atom", type=parse_geometry, help="poistions of atoms", required=True)
parser.add_argument("--Nk", type=int, default=0, help="number of plane-waves in each direction for integral evaluation")
parser.add_argument("--basis", type=str, nargs="*", help="basis sets definition. First specify atom then basis for this atom", required=True)
parser.add_argument("--auxbasis", type=str, nargs="*", default=[None], help="auxiliary basis")
parser.add_argument("--ecp", type=str, nargs="*", default=[None], help="effective core potentials")
parser.add_argument("--xc", type=str, nargs="*", default=[None], help="XC functional")
parser.add_argument("--dm0", type=str, default=None, help="initial guess for density matrix")
parser.add_argument("--df_int", type=int, default=1, help="prepare density fitting integrals or not")
parser.add_argument("--int_path", type=str, default="df_int", help="path to store ewald corrected integrals")
parser.add_argument("--hf_int_path", type=str, default="df_hf_int", help="path to store hf integrals")
parser.add_argument("--output_path", type=str, default="input.h5", help="output file with initial data")
parser.add_argument(
"--orth", type=str, default="none",
choices=["none", "lowdin", "symmetric_lowdin", "mo", "natural", "0", "1"],
help=(
"Orbital basis for stored quantities: "
"'none' = keep AO basis (legacy '0'); "
"'lowdin' = canonical Löwdin V·Lambda^{-1/2} (legacy '1'); "
"'symmetric_lowdin' = Hermitian Löwdin S^{-1/2}; "
"'mo' = canonical MOs from mean-field; "
"'natural' = natural orbitals from mean-field density matrix."
),
)
parser.add_argument("--beta", type=float, default=None, help="Emperical parameter for even-Gaussian auxiliary basis")
parser.add_argument("--active_space", type=int, nargs='+', default=None, help="active space orbitals")
parser.add_argument("--spin", type=int, default=0, help="Local spin")
parser.add_argument("--restricted", type=lambda x: (str(x).lower() in ['true','1', 'yes']), default='false', help="Spin restricted calculations.")
parser.add_argument("--memory", type=int, default=700, help="Memory bound for integral chunk in MB")
parser.add_argument("--grid_only", type=lambda x: (str(x).lower() in ['true','1', 'yes']), default='false', help="Only recompute k-grid points")
parser.add_argument("--diffuse_cutoff", type=float, default=0.0, help="Remove the diffused Gaussians whose exponents are less than the cutoff")
parser.add_argument("--damping", type=float, default=0.0, help="Damping factor for mean-field iterations")
parser.add_argument("--max_iter", type=int, default=100, help="Maximum number of iterations in the SCF loop")
parser.add_argument("--keep_cderi", type=lambda x: (str(x).lower() in ['true','1', 'yes']), default='false', help="Keep generated cderi files.")
parser.add_argument("--job", choices=["init", "sym_path", "ewald_corr"], default="init", nargs="+")
parser.add_argument(
"--x2c", type=int, default=0, choices=[0, 1, 2],
help="enable X2C calculations (0: non-rel., 1: sfX2C1e, 2: X2C1e)"
)
advanced = parser.add_argument_group(
"Advanced options",
"Low-level knobs intended for expert users. Default values are appropriate for most calculations."
)
advanced.add_argument(
"--use_j2c_eig_decomposition",
type=lambda x: (str(x).lower() in ['true', '1', 'yes']),
default=False,
help="Use eigenvalue decomposition for j2c factors during DF build. Set false to force Cholesky-based path."
)
[docs]
def add_pbc_params(parser):
'''
Define PBC-specific command line arguments for Green python module
'''
parser.add_argument("--a", type=parse_geometry, help="lattice geometry", required=True)
parser.add_argument("--nk", type=int, nargs='+', help="number of k-points in each direction. Provide 1 value for symmetric mesh or 3 values for anisotropic mesh.", required=True)
parser.add_argument("--pseudo", type=str, nargs="*", default=[None], help="pseudopotential")
parser.add_argument("--shift", type=float, nargs=3, default=[0.0, 0.0, 0.0], help="mesh shift")
parser.add_argument("--center", type=float, nargs=3, default=[0.0, 0.0, 0.0], help="mesh center")
parser.add_argument("--space_symm", type=lambda x: (str(x).lower() in ['true','1', 'yes']), default='true', help="Use space group symmetry")
parser.add_argument("--tr_symm", type=lambda x: (str(x).lower() in ['true','1', 'yes']), default='true', help="Use time-reversal symmetry")
parser.add_argument("--print_high_symmetry_points", default=False, action='store_true', help="Print available high symmetry points for current system and exit.")
parser.add_argument("--high_symmetry_path", type=str, default=None, help="High symmetry path")
parser.add_argument("--high_symmetry_path_points", type=int, default=0, help="Number of points for high symmetry path")
parser.add_argument("--finite_size_kind", choices=["ewald", "gf2", "gw", "gw_s", "coarse_grained"], default="ewald", nargs="+",
help="Two body finite-size correction. Be default computes the second set of integrals that include simple ewald correction.")
_ORTH_ALIASES = {"0": "none", "1": "lowdin"}
[docs]
def init_mol_params(params=None):
'''
Initialize argparse.ArgumentParser for Green/WeakCoupling python module and return a prased parameters map with parameters specific for molecular calculations
'''
parser = argparse.ArgumentParser(description="Green/WeakCoupling initialization script")
add_common_params(parser)
args = parser.parse_args(args=params)
args.orth = _ORTH_ALIASES.get(args.orth, args.orth)
args.basis = parse_basis(args.basis)
args.auxbasis = parse_basis(args.auxbasis)
args.ecp = parse_basis(args.ecp)
args.xc = parse_basis(args.xc)
if args.x2c == 2 and args.restricted:
raise RuntimeError("X2C calculation can not be spin restricted")
if args.xc is not None:
if args.x2c == 2:
args.mean_field = mdft.GKS
else:
args.mean_field = mdft.RKS if args.restricted else mdft.UKS
else:
if args.x2c == 2:
args.mean_field = mscf.GHF
else:
args.mean_field = mscf.RHF if args.restricted else mscf.UHF
args.ns = 1 if args.restricted or args.x2c == 2 else 2
# parameters needed to create empty grid
args.a = [[1,0,0],[0,1,0],[0,0,1]]
args.nk = [1, 1, 1]
args.shift = [0.,0.,0.]
args.center = [0.,0.,0.]
return args
[docs]
def init_pbc_params(params=None):
'''
Initialize argparse.ArgumentParser for Green/WeakCoupling python module and return a prased parameters map with parameters specific for periodic calculations
'''
parser = argparse.ArgumentParser(description="Green/WeakCoupling initialization script")
add_common_params(parser)
add_pbc_params(parser)
args = parser.parse_args(args=params)
args.orth = _ORTH_ALIASES.get(args.orth, args.orth)
if len(args.nk) == 1:
args.nk = [args.nk[0], args.nk[0], args.nk[0]]
elif len(args.nk) != 3:
raise ValueError("--nk must be given 1 or 3 integers, got {}".format(len(args.nk)))
args.basis = parse_basis(args.basis)
args.auxbasis = parse_basis(args.auxbasis)
args.ecp = parse_basis(args.ecp)
args.pseudo = parse_basis(args.pseudo)
args.xc = parse_basis(args.xc)
if args.x2c == 2 and args.restricted:
raise RuntimeError("X2C calculation can not be spin restricted")
if args.xc is not None:
if args.x2c == 2:
args.mean_field = dft.KGKS
else:
args.mean_field = dft.KRKS if args.restricted else dft.KUKS
else:
if args.x2c == 2:
args.mean_field = scf.KGHF
else:
args.mean_field = scf.KRHF if args.restricted else scf.KUHF
args.ns = 1 if args.restricted or args.x2c == 2 else 2
return args
[docs]
def mol_cell(args):
'''
Initialize PySCF unit cell object for a given system
'''
c = mgto.M(
atom = args.atom,
unit = 'A',
basis = args.basis,
ecp = args.ecp,
verbose = 7,
spin = args.spin
)
return c
[docs]
def pbc_cell(args):
'''
Initialize PySCF unit cell object for a given system
'''
spg_symm = args.space_symm
c = gto.M(
a = args.a,
atom = args.atom,
unit = 'A',
basis = args.basis,
ecp = args.ecp,
pseudo = args.pseudo,
verbose = 7,
spin = args.spin,
space_group_symmetry = spg_symm,
# exp_to_discard = args.diffuse_cutoff
)
_a = c.lattice_vectors()
c.exp_to_discard = args.diffuse_cutoff
c.build()
if np.linalg.det(_a) < 0:
raise "Lattice are not in right-handed coordinate system. Please correct your lattice vectors"
return c
[docs]
def wrap_1stBZ(k):
'''
wrap scaled k into [-0.5,0.5) range
:param k: value of k-point at some dimension
'''
while k < -0.5 :
k = k + 1
while (k - 4.9999999999e-1) > 0.0 :
k = k - 1
return k
[docs]
def init_k_mesh(args, mycell):
"""init k-points mesh for GDF
Parameters
----------
args : map
simulation parameters
mycell : pyscf.pbc.Cell
unit cell for simulation
Returns
-------
numpy.ndarray
k-mesh for the Brillouin Zone
numpy.ndarray
k-mesh with unique k-points, forming the irreducible k-mesh
numpy.ndarray
indices of irreducible k-points in the input k-list
numpy.ndarray
truth table for whether a k-point has an irreducible time-reversal equivalent
numpy.ndarray
weights for each irreducible k-point (degeneracy)
numpy.ndarray
inverse index, associating each k-point with its unique equivalent
int
number of irreducible k-points
"""
if args.center is None:
args.center = [0,0,0]
if args.shift is None:
args.shift = [0,0,0]
kstruct = mycell.make_kpts(args.nk, scaled_center=args.center,
space_group_symmetry=args.space_symm, time_reversal_symmetry=args.tr_symm)
if not (args.space_symm or args.tr_symm):
kstruct = libkpts.make_kpts(mycell, kstruct, space_group_symmetry=False, time_reversal_symmetry=False)
kmesh = kstruct.kpts
for i, kk in enumerate(kmesh):
ki = kmesh[i]
ki = mycell.get_scaled_kpts(ki) + args.shift
ki = [wrap_k(l) for l in ki]
kmesh[i] = mycell.get_abs_kpts(ki)
for i, ki in enumerate(kmesh):
ki = mycell.get_scaled_kpts(ki)
ki = [wrap_k(l) for l in ki]
ki = mycell.get_abs_kpts(ki)
kmesh[i] = ki
logging.debug(kmesh)
logging.info(mycell.get_scaled_kpts(kmesh))
logging.info("Compute irreducible k-points")
k_ibz = kstruct.kpts_ibz
nk = kstruct.nkpts
num_ik = kstruct.nkpts_ibz
ir_list = kstruct.ibz2bz
conj_list = kstruct.time_reversal_symm_bz
ind = kstruct.bz2ibz # gives index of ir_list to which each k-point belongs
# but we want ind to be the index in the main list of k-points
ind = ir_list[ind]
# PySCF stores weight in fractions of 1/nk
weight = ind * 0
weight_ir_list = kstruct.weights_ibz
for i, irr_i in enumerate(ir_list):
weight[irr_i] = weight_ir_list[i] * nk
# wrap IBZ k-points into 1st BZ
for i, ki in enumerate(k_ibz):
ki = mycell.get_scaled_kpts(ki)
ki = [wrap_1stBZ(l) for l in ki]
k_ibz[i] = ki
return kmesh, k_ibz, ir_list, conj_list, weight, ind, num_ik, kstruct
[docs]
def init_q_mesh(args, mycell, k_mesh, save_data=True):
"""Initialize q-mesh for GDF
Parameters
----------
mycell : pyscf.pbc.Cell
unit cell for simulation
k_mesh : numpy.ndarray
k-mesh for the Brillouin Zone
Returns
-------
pyscf.pbc.lib.kpts.KPoints
q-mesh struct for the Brillouin Zone
"""
# NOTE: we use the getattr function because "space_symm" and "tr_symm" CLI arguments are only available for periodic systems
tr_symm = bool(getattr(args, "tr_symm", True))
space_symm = bool(getattr(args, "space_symm", True))
# Unlike k-mesh, the presence or absence of relativity doesn't concern q_mesh structure
qstruct = kpt_utils.build_q_struct(mycell, k_mesh, space_symm=space_symm, tr_symm=tr_symm)
# Obtain all info to save
nq = qstruct.nkpts
q_mesh = qstruct.kpts
num_iq = qstruct.nkpts_ibz
ir_list = qstruct.ibz2bz
conj_list = qstruct.time_reversal_symm_bz
ind = qstruct.bz2ibz # gives index of ir_list to which each k-point belongs
# but we want ind to be the index in the main list of k-points
ind = ir_list[ind]
# PySCF stores weight in fractions of 1/nq
weight = ind * 0
weight_ir_list = qstruct.weights_ibz
for i, irr_i in enumerate(ir_list):
weight[irr_i] = weight_ir_list[i] * nq
# Save q-grid info to output file
if save_data:
inp_data = h5py.File(args.output_path, "a")
grid = inp_data["symmetry"]
if "q" not in grid:
grid.create_group("q")
qgrid = grid["q"]
def _write(path, value):
if path in qgrid:
qgrid[path][...] = value
else:
qgrid[path] = value
_write("mesh", q_mesh)
_write("mesh_scaled", mycell.get_scaled_kpts(q_mesh))
_write("bz2ibz", ind)
_write("weight_ibz", weight)
_write("inq", num_iq)
_write("nq", nq)
_write("ibz2bz", ir_list)
_write("tr_conj", conj_list)
inp_data.close()
return qstruct
[docs]
def read_dm(dm0, dm_file):
'''
Read density matrix from smaller kmesh
'''
nao = dm0.shape[-1]
nkpts = dm0.shape[1]
dm = np.zeros((2,nao,nao),dtype=np.complex128)
f = h5py.File(dm_file, 'r')
dm[:,:,:] = f['/dm_gamma'][:]
f.close()
dm_kpts = np.repeat(dm[:,None, :, :], nkpts, axis=1)
return dm_kpts
[docs]
def solve_mean_field(args, mydf, mycell):
'''
Obtain pySCF mean-field solution for a given parameters, unit-cell object and density-fitting object
'''
logging.info("Solve Mean-field")
# prepare and solve DFT
if args.x2c == 0:
mf = args.mean_field(mycell, mydf.kpts).density_fit()
elif args.x2c == 1:
mf = args.mean_field(mycell, mydf.kpts).density_fit().sfx2c1e()
elif args.x2c == 2:
mf = args.mean_field(mycell, mydf.kpts).density_fit().x2c1e()
if args.xc is not None:
mf.xc = args.xc
#mf.max_memory = 10000
mydf._cderi = "cderi.h5"
mf.kpts = mydf.kpts
mf.with_df = mydf
mf.diis_space = 16
mf.damp = args.damping
mf.max_cycle = args.max_iter
mf.chkfile = 'tmp.chk'
if os.path.exists("tmp.chk"):
init_dm = _init_guess_from_chk(mf, mycell, mf.chkfile)
mf.kernel(init_dm)
elif args.dm0 is not None:
init_dm = mf.get_init_guess()
init_dm = read_dm(init_dm, args.dm0)
mf.kernel(init_dm)
else:
mf.kernel()
if args.x2c < 2:
mf.analyze()
return mf
[docs]
def solve_mol_mean_field(args, mydf, mycell):
'''
Obtain pySCF mean-field solution for a given parameters, unit-cell object and density-fitting object
'''
logging.info("Solve LDA")
# prepare and solve DFT
if args.x2c == 0:
mf = args.mean_field(mycell).density_fit()
elif args.x2c == 1:
mf = args.mean_field(mycell).density_fit().sfx2c1e()
elif args.x2c == 2:
mf = args.mean_field(mycell).x2c1e()
if args.xc is not None:
mf.xc = args.xc
# mf.max_memory = 10000
# mydf._cderi = "cderi.h5"
if args.x2c == 2:
tmp_mf = mscf.RHF(mycell).density_fit() if args.restricted else mscf.UHF(mycell).density_fit()
tmp_mf.with_df._cderi_to_save = "cderi_mol.h5"
tmp_mf.with_df.build()
tmp_mf = None
else:
mf.with_df._cderi_to_save = "cderi_mol.h5"
mf.with_df.build()
mf.diis_space = 16
mf.damp = args.damping
mf.max_cycle = args.max_iter
mf.chkfile = 'tmp.chk'
if os.path.exists("tmp.chk"):
init_dm = _init_guess_from_chk(mf, mycell, mf.chkfile)
mf.kernel(init_dm)
elif args.dm0 is not None:
init_dm = mf.get_init_guess()
init_dm = read_dm(init_dm, args.dm0)
mf.kernel(init_dm)
else:
mf.kernel()
if args.x2c < 2:
mf.analyze()
return mf
[docs]
def store_kstruct_ops_info(args, mycell, kmesh, kstruct, X_k=None, X_inv_k=None):
"""Store symmetry operation information for k-points into a hdf5 file in Green'WeakCoupling format
Parameters
----------
args : map
simulation parameters
mycell : pyscf.pbc.Cell
unit cell for simulation
kmesh : numpy.ndarray
k-mesh for the Brillouin Zone
kstruct : pyscf.pbc.symm.KPointsSymmetry
k-point symmetry structure
X_k : numpy.ndarray, optional
Orthogonalization matrices in the full BZ, shape ``(nk, nao, nao)``
(or spin-orbital analog). Required together with ``X_inv_k`` when
``args.orth != "none"`` to rotate AO-space symmetry operators into the
orthogonalized basis.
X_inv_k : numpy.ndarray, optional
Inverse orthogonalization matrices on irreducible k-points, indexed by
``kstruct.bz2ibz`` representatives. Used as
``X_k[k] @ U_ao[k] @ X_inv_k[k_ir]``.
Notes
-----
The function writes/updates the following datasets under
``/symmetry/k`` in ``args.output_path``:
- ``n_stars``: number of k-point stars.
- ``stars/<i>``: indices of BZ points in star ``i``.
- ``k_sym_transform_ao``: one AO-space symmetry transform per BZ k-point,
mapping each point to its representative irreducible k-point.
For ``args.x2c < 2``, ``k_sym_transform_ao`` is built from
``get_representation``. For ``args.x2c == 2``, the full double-group spinor
representation :math:`D^{1/2}(R^{-1}) \\otimes U_\\text{orbital}(R^{-1})` is stored
via :func:`get_spinor_representation`.
Returns
-------
None
Data is written directly to the HDF5 file.
"""
# extract symmetry operation information from kstruct
nk = kmesh.shape[0]
inp_data = h5py.File(args.output_path, "a")
if "symmetry" in inp_data:
symm_grp = inp_data["symmetry"]
else:
symm_grp = inp_data.create_group("symmetry")
if "k" in symm_grp:
symm_grp = symm_grp["k"]
else:
symm_grp.create_group("k")
symm_grp = symm_grp["k"]
stars_ops = kstruct.stars_ops_bz
stars = kstruct.stars
n_stars = len(stars)
# store number of stars for the k-mesh / k-struct
if "n_stars" in symm_grp:
symm_grp["n_stars"][...] = n_stars
else:
symm_grp["n_stars"] = n_stars
# store stars themselves
if "stars" in symm_grp:
for i in range(n_stars):
symm_grp["stars/{}" .format(i)][...] = stars[i]
else:
star_grp = symm_grp.create_group("stars")
for i in range(n_stars):
star_grp["{}" .format(i)] = stars[i]
# construct symmetry operators in AO basis
# NOTE: only one operator per k-point is stored, the one that connects it to the irreducible k-point
# IMPORTANT: In periodic systems, overlap matrices S_k in the Bloch AO basis include lattice phase factors.
# Therefore, S_k matrices at equivalent k-points differ due to these phases. However, the generalized
# eigenproblem (H, S) IS invariant under the symmetry transformation, with eigenvalues matching to
# machine precision. This validates that the stored rotation matrices are correct for physical transformations.
nao = mycell.nao_nr()
if args.x2c < 2:
kspace_orep = np.zeros((nk, nao, nao), dtype=np.complex128)
# Non-relativistic calculations, where nso = nao (i.e. AO space representation is correct)
for ik in range(nk):
iop = stars_ops[ik]
mat_ao = get_representation(ik, iop, mycell, kstruct)
kspace_orep[ik] = mat_ao
else:
# Relativistic (X2C1e): full double-group spinor representation.
# U_spinor(R) = D^{1/2}(R) ⊗ U_orbital(R); SU(2) lifted from kstruct.ops directly.
# For TR k-points the combined operator (U_spinor·Θ)* is stored so that
# the reconstruction X(k) = (Uk @ X_ir @ Uk†)* gives U·Θ·X_ir*·Θ†·U†.
nso = nao * 2
kspace_orep = np.zeros((nk, nso, nso), dtype=np.complex128)
theta = np.kron(np.array([[0, 1], [-1, 0]], dtype=np.complex128), np.eye(nao))
tr_conj_bz = kstruct.time_reversal_symm_bz
for ik in range(nk):
iop = stars_ops[ik]
u_spinor = get_spinor_representation(ik, iop, mycell, kstruct)
kspace_orep[ik] = (u_spinor @ theta).conj() if tr_conj_bz[ik] else u_spinor
kspace_orep = kspace_orep.astype(np.complex128)
# If quantities are saved in an orthogonalized basis, rotate symmetry operators
# to the same basis so U(k) reconstructs H/F/G consistently.
if args.orth != "none":
if (X_k is None or X_inv_k is None):
raise ValueError(
"Cannot transform symmetry operators to orthogonal basis: "
"missing transformation matrices X_k and/or X_inv_k. "
f"(--orth={args.orth} requires --grid_only=false to compute mean-field quantities)"
)
# get mapping from full BZ idx to idx (still in full BZ) of the corresponding irreducible point
bz2ibz = kstruct.ibz2bz[kstruct.bz2ibz]
kspace_orep_orth = np.zeros_like(kspace_orep)
for ik in range(nk):
ik_ir = bz2ibz[ik]
kspace_orep_orth[ik] = X_k[ik] @ kspace_orep[ik] @ X_inv_k[ik_ir]
kspace_orep = kspace_orep_orth
if "k_sym_transform_ao" in symm_grp:
symm_grp["k_sym_transform_ao"][...] = kspace_orep # .view(np.float64).reshape(kspace_orep.shape + (2,))
# symm_grp["k_sym_transform_ao"].attrs["__complex__"] = np.int8(1)
else:
symm_grp["k_sym_transform_ao"] = kspace_orep # .view(np.float64).reshape(kspace_orep.shape + (2,))
# symm_grp["k_sym_transform_ao"].attrs["__complex__"] = np.int8(1)
inp_data.close()
[docs]
def store_auxcell_kstruct_ops_info(args, auxcell, kmesh):
"""Store symmetry operation information for k-points into hdf5 file in Green'WeakCoupling format
for auxcell only case
Parameters
----------
args : map
simulation parameters
auxcell : pyscf.pbc.gto.Cell
auxiliary unit cell for density-fitting
kmesh : numpy.ndarray
k-mesh for the Brillouin Zone
aux_kstruct : pyscf.pbc.symm.KPointsSymmetry
k-point symmetry structure for aux-basis
"""
# generate periodic cell for auxbasis
auxcell.build()
qstruct = init_q_mesh(args, auxcell, kmesh)
irre_q_inds = qstruct.ibz2bz
stars_ops = qstruct.stars_ops_bz
nk = qstruct.nkpts
nao = auxcell.nao_nr()
# star data for qstruct
stars_ops = qstruct.stars_ops_bz
stars = qstruct.stars
n_stars = len(stars)
# read j2c and compute j2c_sqrt and j2c_sqrt_inv for each k-point using lower Cholesky
# decomposition to match the convention used by PySCF when building j3c integrals.
# PySCF computes B = L^{-1} @ eri3c (lower Cholesky, j2c = LL†), so P0_tilde lives
# in the L-basis. Obar = L_bz^{-1} @ mat_ao @ L_irre correctly maps P0_tilde between
# k-points. Using upper Cholesky gives L^T instead of L, producing the wrong Obar.
import scipy.linalg as LA
j2c_data = h5py.File('cderi.h5', 'r')
j2c_decomp = j2c_data['j2c'].attrs['j2c_decomposition']
first_j2c_key = irre_q_inds[0]
nq = j2c_data[f'j2c/{first_j2c_key}'].shape[0]
assert nq == nao, "number of AOs in auxcell and j2c data do not match"
# We will compute j2c_sqrt for all irreducible k-points once and store
# For j2c_sqrt_inv, we will compute it on the fly as we construct kspace_orep_p0
j2c_sqrt_irre = []
for irre_q in irre_q_inds:
j2c_i = j2c_data[f'j2c/{irre_q}'][...]
j2c_i_dagger = j2c_i.conj().T
assert np.allclose(j2c_i, j2c_i_dagger, atol=1e-10, rtol=0), "j2c metric is not Hermitian"
# make it explicitly hermitian
j2c_i = (j2c_i + j2c_i_dagger) / 2
if j2c_decomp == 'cholesky':
j2c_sqrt_i, _ = int_utils.cholesky_decomposed_metric(j2c_i, auxcell, inv=False)
elif j2c_decomp == 'eigenvalue':
j2c_sqrt_i, _ = int_utils.eigenvalue_decomposed_metric(j2c_i, auxcell, inv=False)
else:
raise ValueError("Unsupported j2c decomposition method: {}".format(j2c_decomp))
j2c_sqrt_irre.append(j2c_sqrt_i)
# compute representation in the AO basis for each k-point and each symmetry operation
# NOTE: only one operator per k-point is stored, the one that connects it to the irreducible k-point
kspace_orep_j2c = np.zeros((nk, nao, nao), dtype=np.complex128)
kspace_orep_p0 = np.zeros((nk, nao, nao), dtype=np.complex128)
for ik in range(nk):
# indexing
iop = stars_ops[ik]
irre_q = qstruct.bz2ibz[ik] # index of irreducible k-point in the ibz list
irre_q_bz = qstruct.ibz2bz[irre_q] # index of irreducible k-point in the full bz list
# Short-circuit when ik is its own IBZ representative: the q->q
# transformation must be identity. Computing it via L_bz^{-1} @ mat_ao @ L_irre
# can drift from identity (e.g., when stars_ops[ik] != identity but acts
# trivially on q, or when Cholesky/eigendecomp is recomputed independently
# from the pre-stored sqrt). The downstream GW kernel relies on U=I at IBZ
# reps in eval_p0_bz_from_ibz.
if ik == irre_q_bz:
kspace_orep_p0[ik] = np.eye(nao, dtype=np.complex128)
kspace_orep_j2c[ik] = np.eye(nao, dtype=np.complex128)
continue
# Build transformation operator in the aux-AO basis connecting "ik" with "irre_k"
mat_ao = get_representation(ik, iop, auxcell, qstruct)
# obtain J^{1/2} (q_IBZ) from pre-computed list
j2c_irre_k_sqrt = j2c_sqrt_irre[irre_q]
# compute J^{-1/2} (k_BZ) on the fly from q_irreducible
j2c_irre_i = j2c_data["j2c/{}".format(irre_q_bz)][...]
j2c_i = mat_ao @ j2c_irre_i @ mat_ao.conj().T
if j2c_decomp == 'cholesky':
j2c_ik_sqrt_inv, _ = int_utils.cholesky_decomposed_metric(j2c_i, auxcell, inv=True)
elif j2c_decomp == 'eigenvalue':
j2c_ik_sqrt_inv, _ = int_utils.eigenvalue_decomposed_metric(j2c_i, auxcell, inv=True)
else:
raise ValueError("Unsupported j2c decomposition method: {}".format(j2c_decomp))
# get effective dimensions
ncols = j2c_irre_k_sqrt.shape[1]
nrows = j2c_ik_sqrt_inv.shape[0]
# transform to j2c basis
kspace_orep_p0[ik, :nrows, :ncols] = j2c_ik_sqrt_inv @ mat_ao @ j2c_irre_k_sqrt
kspace_orep_j2c[ik] = mat_ao
# clean up for next iteration
j2c_irre_i = None
j2c_data.close()
# TODO: Integrate the j2c_neg for 2D systems where the metric can be negative for some k-points.
# Save transformed aux kspace_orep_p0 to hdf5 file
inp_data = h5py.File(args.output_path, "a")
if "symmetry" in inp_data:
symm_grp = inp_data["symmetry"]
else:
symm_grp = inp_data.create_group("symmetry")
if "q" in symm_grp:
symm_grp = symm_grp["q"]
else:
symm_grp.create_group("q")
symm_grp = symm_grp["q"]
# Store kspace transformation matrices
kspace_orep_p0 = kspace_orep_p0.astype(np.complex128)
kspace_orep_j2c = kspace_orep_j2c.astype(np.complex128)
if "k_sym_transform_p0" in symm_grp:
symm_grp["k_sym_transform_p0"][...] = kspace_orep_p0 # .view(np.float64).reshape(kspace_orep_p0.shape + (2,))
# symm_grp["k_sym_transform_p0"].attrs["__complex__"] = np.int8(1)
else:
symm_grp["k_sym_transform_p0"] = kspace_orep_p0 # .view(np.float64).reshape(kspace_orep_j2c.shape + (2,))
# symm_grp["k_sym_transform_p0"].attrs["__complex__"] = np.int8(1)
if "k_sym_transform_j2c" in symm_grp:
symm_grp["k_sym_transform_j2c"][...] = kspace_orep_j2c # view(np.float64).reshape(kspace_orep_j2c.shape + (2,))
# symm_grp["k_sym_transform_j2c"].attrs["__complex__"] = np.int8(1)
else:
symm_grp["k_sym_transform_j2c"] = kspace_orep_j2c # .view(np.float64).reshape(kspace_orep_j2c.shape + (2,))
# symm_grp["k_sym_transform_j2c"].attrs["__complex__"] = np.int8(1)
if "n_stars" in symm_grp:
symm_grp["n_stars"][...] = n_stars
else:
symm_grp["n_stars"] = n_stars
# store stars themselves
if "stars" in symm_grp:
for i in range(n_stars):
symm_grp["stars/{}" .format(i)][...] = stars[i]
else:
star_grp = symm_grp.create_group("stars")
for i in range(n_stars):
star_grp["{}" .format(i)] = stars[i]
inp_data.close()
[docs]
def store_mol_symmetry_info(args, mycell, auxcell, kmesh=None):
"""Store trivial symmetry information for molecular calculations.
Molecular cases use a single Gamma-point only, so the k- and q-mesh symmetry
datasets are all one-point identity mappings.
"""
zero_mesh = np.zeros((1, 3), dtype=np.float64) if kmesh is None else np.asarray(kmesh, dtype=np.float64)
point_index = np.array([0], dtype=np.int64)
pair_index = np.array([[0, 0]], dtype=np.int64)
weight = np.array([1.0], dtype=np.float64)
tr_conj = np.array([0], dtype=np.int64)
star0 = np.array([0], dtype=np.int64)
nao = mycell.nao_nr()
nso = 2 * nao if args.x2c == 2 else nao
naux = auxcell.nao_nr()
k_sym_transform_ao = np.eye(nso, dtype=np.complex128).reshape(1, nso, nso)
q_sym_transform_j2c = np.eye(naux, dtype=np.complex128).reshape(1, naux, naux)
q_sym_transform_p0 = np.eye(naux, dtype=np.complex128).reshape(1, naux, naux)
inp_data = h5py.File(args.output_path, "a")
def _write(path, value):
if path in inp_data:
inp_data[path][...] = value
else:
inp_data[path] = value
_write("symmetry/k/mesh", zero_mesh)
_write("symmetry/k/mesh_scaled", zero_mesh)
_write("symmetry/k/bz2ibz", point_index)
_write("symmetry/k/weight_ibz", weight)
_write("symmetry/k/ink", 1)
_write("symmetry/k/nk", 1)
_write("symmetry/k/ibz2bz", point_index)
_write("symmetry/k/tr_conj", tr_conj)
_write("symmetry/k/n_stars", 1)
_write("symmetry/k/k_sym_transform_ao", k_sym_transform_ao)
if "symmetry/k/stars/0" in inp_data:
inp_data["symmetry/k/stars/0"][...] = star0
else:
inp_data.require_group("symmetry/k/stars")
inp_data["symmetry/k/stars/0"] = star0
_write("symmetry/q/mesh", zero_mesh)
_write("symmetry/q/mesh_scaled", zero_mesh)
_write("symmetry/q/bz2ibz", point_index)
_write("symmetry/q/weight_ibz", weight)
_write("symmetry/q/inq", 1)
_write("symmetry/q/nq", 1)
_write("symmetry/q/ibz2bz", point_index)
_write("symmetry/q/tr_conj", tr_conj)
_write("symmetry/q/n_stars", 1)
_write("symmetry/q/k_sym_transform_j2c", q_sym_transform_j2c)
_write("symmetry/q/k_sym_transform_p0", q_sym_transform_p0)
if "symmetry/q/stars/0" in inp_data:
inp_data["symmetry/q/stars/0"][...] = star0
else:
inp_data.require_group("symmetry/q/stars")
inp_data["symmetry/q/stars/0"] = star0
_write("symmetry/pairs/conj_pairs_list", point_index)
_write("symmetry/pairs/trans_pairs_list", point_index)
_write("symmetry/pairs/kpair_irre_list", point_index)
_write("symmetry/pairs/kpair_idx", pair_index)
_write("symmetry/pairs/num_kpair_stored", 1)
inp_data.close()
[docs]
def store_k_grid(args, mycell, kmesh, k_ibz, ir_list, conj_list, weight, ind, num_ik, kstruct=None):
'''
Store reciprocal space information into a hdf5 file in Green'WeakCoupling format
'''
inp_data = h5py.File(args.output_path, "a")
nk = kmesh.shape[0]
kptij_idx, kij_conj, kij_trans, kpair_irre_list, num_kpair_stored, kptis, kptjs = int_utils.integrals_grid(mycell, kmesh)
logging.info(f"number of reduced k-pairs: {num_kpair_stored}")
def _write(path, value):
if path in inp_data:
inp_data[path][...] = value
else:
inp_data[path] = value
# Structured grid layout (preferred).
_write("symmetry/k/mesh", kmesh)
_write("symmetry/k/mesh_scaled", mycell.get_scaled_kpts(kmesh))
_write("symmetry/k/bz2ibz", ind)
_write("symmetry/k/weight_ibz", weight)
_write("symmetry/k/ink", num_ik)
_write("symmetry/k/nk", nk)
_write("symmetry/k/ibz2bz", ir_list)
_write("symmetry/k/tr_conj", conj_list)
# k-point pairs for integrals.
_write("symmetry/pairs/conj_pairs_list", kij_conj)
_write("symmetry/pairs/trans_pairs_list", kij_trans)
_write("symmetry/pairs/kpair_irre_list", kpair_irre_list)
_write("symmetry/pairs/kpair_idx", kptij_idx)
_write("symmetry/pairs/num_kpair_stored", num_kpair_stored)
# Basic params needed by both grid-only and full MF consumers.
_write("params/nk", nk)
nk_arr = np.atleast_1d(np.array(args.nk, dtype=int))
_write("symmetry/k/nk_list", np.array([nk_arr[0]]*3, dtype=int) if nk_arr.size == 1 else nk_arr)
# Store operators for symmetry operations
if kstruct is not None:
store_kstruct_ops_info(args, mycell, kmesh, kstruct)
inp_data.close()
[docs]
def construct_mol_gdf(args, mycell):
'''
Construct Gaussian Density Fitting obejct for a given parameters and unit cell.
We make sure to disable range-separeting implementation
'''
# Use gaussian density fitting to get fitted densities
mydf = df.GDF(mycell, mycell.kpts)
if args.auxbasis is not None:
mydf.auxbasis = args.auxbasis
elif args.beta is not None:
mydf.auxbasis = df.aug_etb(mycell, beta=args.beta)
# Coulomb kernel mesh
if args.Nk > 0:
mydf.mesh = [args.Nk, args.Nk, args.Nk]
return mydf
[docs]
def construct_gdf(args, mycell, kmesh=None):
'''
Construct Gaussian Density Fitting obejct for a given parameters and unit cell.
We make sure to disable range-separeting implementation
'''
# Use gaussian density fitting to get fitted densities
mydf = int_utils.GreenGDF(mycell)
mydf.space_symm = bool(getattr(args, "space_symm", False))
mydf.tr_symm = bool(getattr(args, "tr_symm", False))
mydf.x2c = int(getattr(args, "x2c", 0))
mydf.use_j2c_eig_decomposition = bool(getattr(args, "use_j2c_eig_decomposition", True))
if hasattr(mydf, "_prefer_ccdf"):
mydf._prefer_ccdf = True # Disable RS-GDF switch for new pyscf versions
if args.auxbasis is not None:
mydf.auxbasis = args.auxbasis
elif args.beta is not None:
mydf.auxbasis = df.aug_etb(mycell, beta=args.beta)
# Coulomb kernel mesh
if args.Nk > 0:
mydf.mesh = [args.Nk, args.Nk, args.Nk]
if kmesh is not None:
mydf.kpts = kmesh
return mydf
[docs]
def compute_ewald_correction(args, cell, kmesh, filename, X_k=None):
# Use gaussian density fitting to get fitted densities
mydf = int_utils.GreenGDF(cell)
mydf.space_symm = bool(args.space_symm)
mydf.tr_symm = bool(args.tr_symm)
mydf.x2c = int(args.x2c)
mydf.use_j2c_eig_decomposition = bool(getattr(args, "use_j2c_eig_decomposition", True))
if args.auxbasis is not None:
mydf.auxbasis = args.auxbasis
elif args.beta is not None:
mydf.auxbasis = df.aug_etb(cell, beta=args.beta)
# Coulomb kernel mesh
if args.Nk > 0:
mydf.mesh = [args.Nk, args.Nk, args.Nk]
int_utils.compute_ewald_correction(args, mydf, kmesh, cell.nao_nr(), filename, X_k=X_k)
[docs]
def compute_df_int_dca(args, mycell, kmesh, lattice_kmesh, nao, X_k):
"""Generate density-fitting integrals for correlated methods using q-averaging over the super-lattice points to compensate finite-size error
Parameters
----------
args : map
simulation parameters
mycell : pyscf.pbc.Cell or pyscf.Mol
unit cell object
kmesh : numpy.ndarray
reciprocal space grid
lattice_kmesh : numpy.ndarray
super-lattice k-points
nao : int
number of atomic orbitals in the unit cell
X_k : numpy.ndarray
trasformation matrix for projection onto an orthogonal space
"""
if not bool(args.df_int):
return
mydf = construct_gdf(args, mycell, kmesh)
# Use Ewald for divergence treatment
mydf.exxdiv = 'ewald'
weighted_coulG_old = int_utils.GreenGDF.weighted_coulG
int_utils.GreenGDF.weighted_coulG = int_utils.weighted_coulG_ewald
old_get_coulG = tools.get_coulG
tools.get_coulG = lambda cell, k=np.zeros(3), exx=False, mf=None, mesh=None, Gv=None, wrap_around=True, omega=None, **kwargs: int_utils.get_coarsegrained_coulG(lattice_kmesh, cell, k, exx, mf, mesh, Gv,
wrap_around, omega, **kwargs)
kij_conj, kij_trans, kpair_irre_list, kptij_idx, num_kpair_stored = int_utils.compute_integrals(mycell, mydf, kmesh, nao, X_k, args.int_path, "cderi_ewald_dca.h5", False)
mydf = None
mydf = construct_gdf(args, mycell, kmesh)
int_utils.GreenGDF.weighted_coulG = weighted_coulG_old
int_utils.compute_integrals(mycell, mydf, kmesh, nao, X_k, args.hf_int_path, "cderi_dca.h5", False)
tools.get_coulG = old_get_coulG