From 65702e972506dc26a76ef159cff865cef9d161d8 Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Thu, 31 Oct 2024 09:14:12 +0900 Subject: [PATCH] add documentation in source codes --- src/hwave/dos.py | 79 +++++- src/hwave/qlms.py | 10 + src/hwave/qlmsio/read_input.py | 60 +++++ src/hwave/qlmsio/read_input_k.py | 38 +++ src/hwave/qlmsio/wan90.py | 52 +++- src/hwave/solver/base.py | 53 ++++ src/hwave/solver/perf.py | 38 ++- src/hwave/solver/rpa.py | 379 +++++++++++++++++++++++--- src/hwave/solver/uhfk.py | 414 +++++++++++++++++++++++++---- src/hwave/solver/uhfr.py | 439 ++++++++++++++++++++++++++++--- 10 files changed, 1447 insertions(+), 115 deletions(-) diff --git a/src/hwave/dos.py b/src/hwave/dos.py index 2f8854b..7b62ae5 100644 --- a/src/hwave/dos.py +++ b/src/hwave/dos.py @@ -9,12 +9,37 @@ class DoS: + """ + A class to represent the Density of States (DoS). + + Attributes + ---------- + dos : np.ndarray + The density of states array. + ene : np.ndarray + The energy levels array. + ene_num : int + The number of energy levels. + norb : int + The number of orbitals. + """ + dos: np.ndarray ene: np.ndarray ene_num: int norb: int def __init__(self, ene: np.ndarray, dos: np.ndarray): + """ + Initialize the DoS object. + + Parameters + ---------- + ene : np.ndarray + The energy levels array. + dos : np.ndarray + The density of states array. + """ assert ene.shape[0] == dos.shape[1] self.ene = ene self.dos = dos @@ -22,6 +47,16 @@ def __init__(self, ene: np.ndarray, dos: np.ndarray): self.norb = dos.shape[0] def plot(self, filename: str = "", verbose: bool = False): + """ + Plot the density of states. + + Parameters + ---------- + filename : str, optional + The filename to save the plot (default is ""). + verbose : bool, optional + If True, print additional information (default is False). + """ try: import matplotlib.pyplot as plt except ImportError: @@ -45,6 +80,16 @@ def plot(self, filename: str = "", verbose: bool = False): plt.close() def write_dos(self, output: str, verbose: bool = False): + """ + Write the density of states to a file. + + Parameters + ---------- + output : str + The output filename. + verbose : bool, optional + If True, print additional information (default is False). + """ if verbose: print("Writing DOS to file: ", output) total_dos = np.sum(self.dos, axis=0) @@ -60,21 +105,51 @@ def write_dos(self, output: str, verbose: bool = False): fw.write("{:15.8f} ".format(self.dos[j, i])) fw.write("\n") - def __read_geom(file_name="./dir-model/zvo_geom.dat"): + """ + Read the geometry from a file. + + Parameters + ---------- + file_name : str, optional + The filename to read the geometry from (default is "./dir-model/zvo_geom.dat"). + + Returns + ------- + np.ndarray + A 3x3 array representing the geometry. + """ with open(file_name, "r") as fr: uvec = np.zeros((3, 3)) for i, line in enumerate(itertools.islice(fr, 3)): # take first 3 lines uvec[i, :] = np.array(line.split()) return uvec - def calc_dos( input_dict: dict, ene_window: list | None = None, ene_num: int = 101, verbose: bool = False, ) -> DoS: + """ + Calculate the density of states (DoS). + + Parameters + ---------- + input_dict : dict + Dictionary containing input parameters and file paths. + ene_window : list, optional + List containing the energy window [ene_low, ene_high]. If None, defaults to [ene_min - 0.2, ene_max + 0.2]. + ene_num : int, optional + Number of energy points (default is 101). + verbose : bool, optional + If True, print additional information (default is False). + + Returns + ------- + DoS + An instance of the DoS class containing the calculated density of states. + """ try: import libtetrabz except ImportError: diff --git a/src/hwave/qlms.py b/src/hwave/qlms.py index c7b5068..447d0eb 100644 --- a/src/hwave/qlms.py +++ b/src/hwave/qlms.py @@ -16,6 +16,16 @@ from requests.structures import CaseInsensitiveDict def run(*, input_dict: Optional[dict] = None, input_file: Optional[str] = None): + """ + Run the main process with the given input dictionary or input file. + + Parameters: + input_dict (Optional[dict]): A dictionary containing input parameters. + input_file (Optional[str]): A path to a TOML file containing input parameters. + + Raises: + RuntimeError: If neither input_dict nor input_file is provided, or if both are provided. + """ if input_dict is None: if input_file is None: raise RuntimeError("Neither input_dict nor input_file are passed") diff --git a/src/hwave/qlmsio/read_input.py b/src/hwave/qlmsio/read_input.py index cbf8d30..513d4bf 100644 --- a/src/hwave/qlmsio/read_input.py +++ b/src/hwave/qlmsio/read_input.py @@ -7,8 +7,20 @@ class QLMSInput(): + """ + Class to handle QLMS input files and parameters. + """ + valid_namelist = ["trans", "coulombinter", "coulombintra", "pairhop", "hund", "exchange", "ising", "pairlift", "interall", "initial", "onebodyg"] + def __init__(self, file_name_list, solver_type="UHFr"): + """ + Initialize the QLMSInput class. + + Parameters: + file_name_list (dict): Dictionary of file names. + solver_type (str): Type of solver to use. Default is "UHFr". + """ self.file_names = file_name_list self.ham_param = CaseInsensitiveDict() self.ham_param["Transfer"] = self._read_ham("trans", value_type="complex") @@ -28,6 +40,15 @@ def __init__(self, file_name_list, solver_type="UHFr"): self.green["OneBodyG"] = self._read_green("onebodyg") def get_param(self, key): + """ + Get parameters based on the key. + + Parameters: + key (str): Key to identify the parameter type. + + Returns: + dict or None: Returns the corresponding parameter dictionary or None if the key is invalid. + """ if key == "mod" or key == "parameter": #return self.mod_param return None @@ -41,6 +62,16 @@ def get_param(self, key): return None def _read_para(self, file_key, start_line=5): + """ + Read parameters from a file. + + Parameters: + file_key (str): Key to identify the file. + start_line (int): Line number to start reading from. Default is 5. + + Returns: + CaseInsensitiveDict: Dictionary of parameters read from the file. + """ file_name = self.file_names[file_key] value = CaseInsensitiveDict() with open(file_name, "r") as f: @@ -54,6 +85,16 @@ def _read_para(self, file_key, start_line=5): return value def _read_ham(self, file_key, value_type="real"): + """ + Read Hamiltonian parameters from a file. + + Parameters: + file_key (str): Key to identify the file. + value_type (str): Type of values to read ("real" or "complex"). Default is "real". + + Returns: + dict or None: Dictionary of Hamiltonian parameters or None if the file key is not found. + """ if file_key in self.file_names: file_name = self.file_names[file_key] return self._load_text(file_name, value_type) @@ -61,6 +102,15 @@ def _read_ham(self, file_key, value_type="real"): return None def _read_green(self, file_key): + """ + Read Green's function data from a file. + + Parameters: + file_key (str): Key to identify the file. + + Returns: + numpy.ndarray or None: Array of Green's function data or None if the file key is not found. + """ if file_key in self.file_names: file_name = self.file_names[file_key] data = np.loadtxt(file_name, skiprows = 5) @@ -69,6 +119,16 @@ def _read_green(self, file_key): return data def _load_text(self, file_name, value_type): + """ + Load text data from a file. + + Parameters: + file_name (str): Name of the file to read. + value_type (str): Type of values to read ("real" or "complex"). + + Returns: + dict: Dictionary of data read from the file. + """ if value_type == "real": value_width = 1 _make_value = lambda v: float(v[0]) diff --git a/src/hwave/qlmsio/read_input_k.py b/src/hwave/qlmsio/read_input_k.py index fd5ae51..4c7e4f3 100644 --- a/src/hwave/qlmsio/read_input_k.py +++ b/src/hwave/qlmsio/read_input_k.py @@ -10,9 +10,19 @@ class QLMSkInput(): + """ + Class to handle QLMS input files and parameters. + """ valid_namelist = [s.lower() for s in ["path_to_input", "Geometry", "Transfer", "CoulombIntra", "CoulombInter", "Hund", "Ising", "PairLift", "Exchange", "PairHop", "Extern"]] def __init__(self, info_inputfile, solver_type="UHFk"): + """ + Initialize the QLMSkInput object. + + Parameters: + info_inputfile (dict): Dictionary containing input file information. + solver_type (str): Type of solver to use. Default is "UHFk". + """ logger.debug(">>> QLMSkInput init") # [file.input] @@ -91,6 +101,16 @@ def __init__(self, info_inputfile, solver_type="UHFk"): self.green["onebodyg_uhf"] = self._read_green(file_name) def _read_data(self, file_name, value_type="real"): + """ + Read data from a file. + + Parameters: + file_name (str): Name of the file to read. + value_type (str): Type of values in the file ("real" or "complex"). Default is "real". + + Returns: + dict: Dictionary containing the data read from the file. + """ info = {} try: data = np.loadtxt(file_name, skiprows = 5) @@ -110,6 +130,15 @@ def _read_data(self, file_name, value_type="real"): return info def _read_green(self, file_name): + """ + Read green function data from a file. + + Parameters: + file_name (str): Name of the file to read. + + Returns: + numpy.ndarray: Array containing the green function data. + """ try: _data = np.loadtxt(file_name, dtype=np.int32, skiprows = 5) except FileNotFoundError: @@ -118,6 +147,15 @@ def _read_green(self, file_name): return _data def get_param(self, key): + """ + Get parameters based on the provided key. + + Parameters: + key (str): Key to specify which parameters to return ("mod", "ham", "output", etc.). + + Returns: + dict or None: Dictionary containing the requested parameters or None if the key is invalid. + """ if key == "mod" or key == "parameter": return None elif key == "ham" or key == "hamiltonian": diff --git a/src/hwave/qlmsio/wan90.py b/src/hwave/qlmsio/wan90.py index 183e8e0..6fc0511 100644 --- a/src/hwave/qlmsio/wan90.py +++ b/src/hwave/qlmsio/wan90.py @@ -7,6 +7,16 @@ logger = logging.getLogger("qlms").getChild("wan90") def read_geom(name_in): + """ + Reads the geometry from a file. + + Parameters: + name_in (str): The name of the input file. + + Returns: + dict: A dictionary containing the number of orbitals (`norb`), + the real space vectors (`rvec`), and the center positions (`center`). + """ try: with open(name_in, 'r') as f: # skip header @@ -30,7 +40,18 @@ def read_geom(name_in): data = {"norb": norb, "rvec": rvec, "center": center} return data + def read_geometry(name_in): + """ + Reads the geometry information from a file. + + Parameters: + name_in (str): The name of the input file. + + Returns: + dict: A dictionary containing unit vectors, degree, cell vectors, + site to vector mapping, and number of orbitals. + """ info_geometry = {} unit_vec_line_end = 3 cell_vec_line_start = 4 @@ -61,7 +82,15 @@ def read_geometry(name_in): info_geometry["n_orb"] = norb return info_geometry + def write_geom(name_out, info_geometry): + """ + Writes the geometry information to a file. + + Parameters: + name_out (str): The name of the output file. + info_geometry (dict): A dictionary containing the geometry information. + """ with open(name_out, "w") as fw: #Unit_cell for vec in info_geometry["unit_vec"]: @@ -74,7 +103,17 @@ def write_geom(name_out, info_geometry): pos = idx*1.0/(n_orb+1) fw.write("{} {} {}\n".format(pos, pos, pos)) + def read_w90(name_in): + """ + Reads the Wannier90 data from a file. + + Parameters: + name_in (str): The name of the input file. + + Returns: + dict: A dictionary containing the interaction data. + """ try: with open(name_in, 'r') as f: # skip header @@ -103,7 +142,17 @@ def read_w90(name_in): return data + def write_w90(name_in, info_int, info_geometry, interact_shape): + """ + Writes the Wannier90 data to a file. + + Parameters: + name_in (str): The name of the output file. + info_int (dict): A dictionary containing the interaction information. + info_geometry (dict): A dictionary containing the geometry information. + interact_shape (tuple): A tuple containing the shape of the interaction. + """ with open(name_in, "w") as fw: fw.write("# wannier90 format\n") norb = info_int["n_orb"] @@ -122,10 +171,11 @@ def write_w90(name_in, info_int, info_geometry, interact_shape): int_value = value[1] fw.write("{} {} {} {} {} {} {}\n".format(site_org, site[0], site[1], site[2]+1, site[3]+1, int_value.real, int_value.imag)) + if __name__ == "__main__": path_to_sample = "../../sample/dir-model" import os data_geom = read_geom(os.path.join(path_to_sample, "zvo_geom.dat")) data_hr = read_w90(os.path.join(path_to_sample,"zvo_hr.dat")) print(data_geom) - print(data_hr) + print(data_hr) \ No newline at end of file diff --git a/src/hwave/solver/base.py b/src/hwave/solver/base.py index 71780eb..6200d55 100644 --- a/src/hwave/solver/base.py +++ b/src/hwave/solver/base.py @@ -6,7 +6,18 @@ logger = logging.getLogger("qlms").getChild("solver") class solver_base(): + """ + Base class for the solver. + """ def __init__(self, param_ham, info_log, info_mode, param_mod=None): + """ + Initialize the solver base. + + :param param_ham: Hamiltonian parameters. + :param info_log: Information log. + :param info_mode: Information mode. + :param param_mod: Optional parameter modifications. + """ logger = logging.getLogger("qlms").getChild(self.name) self.param_mod = param_mod @@ -118,6 +129,12 @@ def __init__(self, param_ham, info_log, info_mode, param_mod=None): # pprint(self.param_mod) def _check_info_mode(self, info_mode): + """ + Check the information mode for validity. + + :param info_mode: Information mode to check. + :return: Exit code (0 if valid, non-zero if invalid). + """ logger = logging.getLogger("qlms").getChild(self.name) fix_list = { @@ -136,6 +153,12 @@ def _check_info_mode(self, info_mode): return exit_code def _check_param_mod(self, param_mod): + """ + Check the parameter modifications for validity. + + :param param_mod: Parameter modifications to check. + :return: Error code (0 if valid, non-zero if invalid). + """ logger = logging.getLogger("qlms").getChild(self.name) error_code = 0 @@ -149,6 +172,13 @@ def _check_param_mod(self, param_mod): return error_code def _check_param_range(self, param_mod, range_list): + """ + Check if the parameters are within the specified ranges. + + :param param_mod: Parameter modifications to check. + :param range_list: List of valid ranges for parameters. + :return: Error code (0 if valid, non-zero if invalid). + """ logger = logging.getLogger("qlms").getChild(self.name) error_code = 0 @@ -166,6 +196,13 @@ def _check_param_range(self, param_mod, range_list): return error_code def _round_to_int(self, val, mode): + """ + Round a value to an integer based on the specified mode. + + :param val: Value to round. + :param mode: Rounding mode (e.g., "round", "round-up"). + :return: Rounded integer value. + """ import math mode = mode.lower() # case-insensitive if mode == "as-is": @@ -197,10 +234,26 @@ def _round_to_int(self, val, mode): return ret def solve(self, path_to_output): + """ + Solve the problem and output the results. + + :param path_to_output: Path to the output file. + """ pass def get_results(self): + """ + Get the results of the solver. + + :return: Tuple containing physics and Green's function results. + """ return (self.physics, self.Green) def save_results(self, info_outputfile, green_info): + """ + Save the results to the specified output file. + + :param info_outputfile: Information about the output file. + :param green_info: Information about the Green's function. + """ pass diff --git a/src/hwave/solver/perf.py b/src/hwave/solver/perf.py index b7c9677..8331ea9 100644 --- a/src/hwave/solver/perf.py +++ b/src/hwave/solver/perf.py @@ -2,10 +2,19 @@ import time class PerfDB: + """ + A class to store and manage performance data for function calls. + """ def __init__(self): + """ + Initializes the PerfDB instance with empty dictionaries for count and value. + """ self._db_count = {} self._db_value = {} def __del__(self): + """ + Prints the performance statistics when the PerfDB instance is deleted. + """ if len(self._db_count) == 0: return print("--------------------------------------------------------------------------------") @@ -21,14 +30,40 @@ def __del__(self): )) print("--------------------------------------------------------------------------------") def put(self, name, value): + """ + Adds a performance entry to the database. + + Args: + name (str): The name of the function. + value (float): The elapsed time for the function call. + """ self._db_count[name] = self._db_count.get(name, 0) + 1 self._db_value[name] = self._db_value.get(name, 0) + value _perf_db_data = PerfDB() def do_profile(func): + """ + A decorator that profiles the execution time of a function and stores it in PerfDB. + + Args: + func (function): The function to be profiled. + + Returns: + function: The wrapped function with profiling. + """ @wraps(func) def wrapper(*args, **kwargs): + """ + Wrapper function that profiles the execution time of the decorated function. + + Args: + *args: Variable length argument list for the decorated function. + **kwargs: Arbitrary keyword arguments for the decorated function. + + Returns: + The result of the decorated function. + """ # start time _start = time.perf_counter() @@ -41,8 +76,7 @@ def wrapper(*args, **kwargs): # elapsed time _elapsed = _end - _start - #print("{}: elapsed time {:.3f} ms".format(func.__name__, _elapsed * 1000)) + # Store the elapsed time in PerfDB _perf_db_data.put(func.__module__ + '.' + func.__name__, _elapsed) return result - return wrapper diff --git a/src/hwave/solver/rpa.py b/src/hwave/solver/rpa.py index 5657cec..c6aff92 100644 --- a/src/hwave/solver/rpa.py +++ b/src/hwave/solver/rpa.py @@ -51,15 +51,30 @@ class Lattice: Note: Bx, By, Bz must be chosen so that Lx, Ly, Lz are divisable by Bx, By, Bz, respectively. otherwise, initialization fails. - """ + def __init__(self, params): + """ + Initialize the Lattice object. + + Args: + params (dict): Dictionary containing lattice parameters. + """ logger.debug(">>> Lattice.__init__") self._init_lattice(params) self._show_params() def _init_lattice(self, params): + """ + Initialize lattice parameters from the given dictionary. + + Args: + params (dict): Dictionary containing lattice parameters. + + Raises: + SystemExit: If required parameters are missing or invalid. + """ logger.debug(">>> Lattice._init_lattice") if "CellShape" not in params: @@ -118,6 +133,9 @@ def _init_lattice(self, params): self.nvol = nvol def _show_params(self): + """ + Log the lattice parameters. + """ logger.info("Lattice parameters:") logger.info(" CellShape = {}".format(self.cellshape)) logger.info(" cell volume = {}".format(self.cellvol)) @@ -133,6 +151,14 @@ class Interaction: Construct Hamiltonian from input """ def __init__(self, lattice, param_ham, info_mode): + """ + Initialize the Interaction object. + + Args: + lattice (Lattice): The lattice object. + param_ham (dict): Dictionary containing Hamiltonian parameters. + info_mode (dict): Dictionary containing mode information. + """ logger.debug(">>> Interaction.__init__") self.lattice = lattice @@ -143,7 +169,6 @@ def __init__(self, lattice, param_ham, info_mode): self._has_interaction_exchange = False self._has_interaction_pairhop = False - # mode options self.enable_spin_orbital = info_mode.get("enable_spin_orbital", False) @@ -159,18 +184,45 @@ def __init__(self, lattice, param_ham, info_mode): pass def has_interaction(self): + """ + Check if there is any interaction. + + Returns: + bool: True if there is any interaction, False otherwise. + """ return self._has_interaction def has_interaction_coulomb(self): + """ + Check if there is Coulomb interaction. + + Returns: + bool: True if there is Coulomb interaction, False otherwise. + """ return self._has_interaction_coulomb def has_interaction_exchange(self): + """ + Check if there is exchange interaction. + + Returns: + bool: True if there is exchange interaction, False otherwise. + """ return self._has_interaction_exchange def has_interaction_pairhop(self): + """ + Check if there is pair hopping interaction. + + Returns: + bool: True if there is pair hopping interaction, False otherwise. + """ return self._has_interaction_pairhop def _init_interaction(self): + """ + Initialize interaction parameters. + """ logger.debug(">>> Interaction._init_interaction") # reinterpret interaction coefficient on sublattice @@ -194,6 +246,15 @@ def _init_interaction(self): pass def _reshape_geometry(self, geom): + """ + Reshape the geometry parameters for the sublattice. + + Args: + geom (dict): Dictionary containing geometry parameters. + + Returns: + dict: Reshaped geometry parameters. + """ logger.debug(">>> Interaction._reshape_geometry") Bx,By,Bz = self.lattice.subshape @@ -219,6 +280,16 @@ def _reshape_geometry(self, geom): return geom_new def _reshape_interaction(self, ham, enable_spin_orbital): + """ + Reshape the interaction parameters for the sublattice. + + Args: + ham (dict): Dictionary containing Hamiltonian parameters. + enable_spin_orbital (bool): Flag to enable spin-orbital interaction. + + Returns: + dict: Reshaped Hamiltonian parameters. + """ logger.debug(">>> Interaction._reshape_interaction") Bx,By,Bz = self.lattice.subshape @@ -248,7 +319,7 @@ def _round(x, n): for bz,by,bx in itertools.product(range(Bz),range(By),range(Bx)): - # original cell index of both endes + # original cell index of both ends # x0 -> x1=x0+r x0,y0,z0 = bx, by, bz x1,y1,z1 = x0 + rx, y0 + ry, z0 + rz @@ -273,6 +344,13 @@ def _round(x, n): return ham_new def _export_interaction(self, type, file_name): + """ + Export the interaction parameters to a file. + + Args: + type (str): Type of interaction. + file_name (str): Name of the file to save the interaction parameters. + """ logger.debug(">>> Interaction._export_interaction") intr = self.param_ham[type] @@ -307,6 +385,16 @@ def _export_interaction(self, type, file_name): )) def _make_ham_trans(self): + """ + Create the Hamiltonian transfer terms and perform Fourier transform. + + This method initializes the Hamiltonian transfer terms from the parameters, + reshapes them, and performs a Fourier transform to obtain the Hamiltonian + in momentum space. + + Raises: + Warning: If 'Transfer' key is not found in the Hamiltonian parameters. + """ logger.debug(">>> Interaction._make_ham_trans") nx,ny,nz = self.lattice.shape @@ -323,7 +411,7 @@ def _make_ham_trans(self): self.ham_extern_q = None return - if self.enable_spin_orbital == True: + if self.enable_spin_orbital: # assume orbital index includes spin index tab_r = np.zeros((nx,ny,nz,nd,nd), dtype=np.complex128) @@ -383,9 +471,20 @@ def _make_ham_trans(self): else: self.ham_extern_r = None self.ham_extern_q = None + def _make_ham_inter(self): + """ + Create the interaction Hamiltonian terms and perform Fourier transform. + This method initializes the interaction Hamiltonian terms from the parameters, + reshapes them, and performs a Fourier transform to obtain the Hamiltonian + in momentum space. - def _make_ham_inter(self): + The interaction Hamiltonian includes Coulomb, Hund, Ising, PairLift, Exchange, + and PairHop interactions. + + Raises: + Warning: If interaction types are not found in the Hamiltonian parameters. + """ logger.debug(">>> Interaction._make_ham_inter") nx,ny,nz = self.lattice.shape @@ -395,7 +494,7 @@ def _make_ham_inter(self): nd = norb * ns # Interaction Hamiltonian W[r,b,bp,a,ap] - # H = W(r)^{\beta\beta^\prime\alpha\alpah^\prime} + # H = W(r)^{\beta\beta^\prime\alpha\alpha^\prime} # * c_{i\alpha}^\dagger c_{i\alpha^\prime} c_{j\beta^\prime}^\dagger c_{j\beta} ham_r = np.zeros((nx,ny,nz,*(ns,norb)*4), dtype=np.complex128) @@ -407,12 +506,17 @@ def _make_ham_inter(self): 'Ising': { (0,0,0,0): 1, (1,1,1,1): 1, (0,0,1,1): -1, (1,1,0,0): -1 }, 'PairLift': { (0,1,0,1): 1, (1,0,1,0): 1 }, 'Exchange': { (0,1,1,0): -1, (1,0,0,1): -1 }, - #-- 'PairHop': { (0,0,1,1): 1, (1,1,0,0): 1 }, } # coulomb-type interactions def _append_inter(type): + """ + Append Coulomb-type interactions to the Hamiltonian. + + Args: + type (str): Type of Coulomb interaction. + """ logger.debug("_append_inter {}".format(type)) spins = spin_table[type] for (irvec,orbvec), v in self.param_ham[type].items(): @@ -429,11 +533,17 @@ def _append_inter(type): # c_{i\alpha\down}^\dagger c_{j\alpha^\prime\down} # + (up <-> down) def _append_pairhop(type): + """ + Append PairHop-type interactions to the Hamiltonian. + + Args: + type (str): Type of PairHop interaction. + """ spins = spin_table[type] for (irvec,orbvec), v in self.param_ham[type].items(): # take account of same-site interaction only - if (irvec == (0,0,0)): - a, b = orbvec + if irvec == (0, 0, 0): + a, b = orbvec for spinvec, w in spins.items(): s1,s2,s3,s4 = spinvec # beta beta' alpha alpha' @@ -496,6 +606,14 @@ class RPA: """ @do_profile def __init__(self, param_ham, info_log, info_mode): + """ + Initialize the RPA object. + + Args: + param_ham (dict): Dictionary containing Hamiltonian parameters. + info_log (dict): Dictionary containing log information. + info_mode (dict): Dictionary containing mode information. + """ logger.debug(">>> RPA.__init__") self.param_ham = param_ham @@ -513,7 +631,14 @@ def __init__(self, param_ham, info_log, info_mode): pass def _set_scheme(self, info_mode): - # handle calc_scheme: must be called after setting up interactions + """ + Set the calculation scheme based on the provided mode information. + + This method handles the `calc_scheme` parameter, which must be called after setting up interactions. + + Args: + info_mode (dict): Dictionary containing mode information. + """ self.calc_scheme = info_mode.get("calc_scheme", "auto") @@ -542,6 +667,13 @@ def _set_scheme(self, info_mode): self.enable_reduced = self.calc_scheme.lower() in ["reduced", "squashed"] def _init_param(self): + """ + Initialize parameters for the RPA object. + + This method sets up various parameters such as temperature, energy cutoff, + number of Matsubara frequencies, number of orbitals, spin degrees of freedom, + and coefficients for tail and external fields. + """ logger.debug(">>> RPA._init_param") self.T = self.param_mod.get("T", 0.0) @@ -618,6 +750,16 @@ def _init_param(self): pass def _round_to_int(self, val, mode): + """ + Round a value to an integer based on the specified mode. + + Args: + val (float): The value to be rounded. + mode (str): The rounding mode. Can be one of "as-is", "round", "round-up", or "round-down". + + Returns: + int: The rounded integer value. + """ import math mode = mode.lower() # case-insensitive if mode == "as-is": @@ -649,6 +791,9 @@ def _round_to_int(self, val, mode): return ret def _show_params(self): + """ + Log the RPA parameters. + """ logger.debug(">>> RPA._show_params") logger.info("RPA parameters:") logger.info(" norbit = {}".format(self.norb)) @@ -678,6 +823,13 @@ def _show_params(self): @do_profile def solve(self, green_info, path_to_output): + """ + Solve the RPA calculations. + + Args: + green_info (dict): Dictionary containing Green's function information. + path_to_output (str): Path to the output directory. + """ logger.info("Start RPA calculations") beta = 1.0/self.T @@ -861,6 +1013,13 @@ def solve(self, green_info, path_to_output): @do_profile def save_results(self, info_outputfile, green_info): + """ + Save the RPA results to the specified output file. + + Args: + info_outputfile (dict): Dictionary containing output file information. + green_info (dict): Dictionary containing Green's function information. + """ logger.info("Save RPA results") path_to_output = info_outputfile["path_to_output"] @@ -892,7 +1051,15 @@ def save_results(self, info_outputfile, green_info): pass def _init_wavevec(self): - # wave vectors on sublatticed geometry + """ + Initialize wave vectors on sublatticed geometry. + + This method calculates the reciprocal lattice vectors and the wave number table + for the sublatticed geometry. + + The wave vectors are stored in `self.kvec` and the wave number table is stored + in `self.wavenum_table` and `self.wave_table`. + """ def _klist(n): return np.roll( (np.arange(n)-(n//2)), -(n//2) ) @@ -901,7 +1068,8 @@ def _klist(n): omg = np.dot(rvec[0], np.cross(rvec[1], rvec[2])) kvec = np.array([ np.cross(rvec[(i+1)%3], rvec[(i+2)%3])/omg * 2*np.pi/self.lattice.shape[i] - for i in range(3) ]) + for i in range(3) + ]) self.kvec = kvec # store reciprocal lattice vectors @@ -917,15 +1085,23 @@ def _klist(n): v = kvec[0] * kx + kvec[1] * ky + kvec[2] * kz wtable[ix,iy,iz] = v self.wave_table = wtable.reshape(nvol,3) - def _find_index_range(self, freq_range): - # decode matsubara frequency index list - # 1. single index - # 2. min, max, step - # 3. keyword - # note index n in [0 .. Nmat-1] corresponds to - # w_n = (2*n-Nmat) * pi / beta + """ + Decode Matsubara frequency index list. + + This method decodes the Matsubara frequency index list based on the provided + frequency range. The frequency range can be specified as a single index, a + range with min and max, a range with min, max, and step, or a keyword. + + Args: + freq_range (int, list, str): The frequency range specification. + Returns: + list: The list of frequency indices. + + Raises: + ValueError: If the frequency range specification is invalid. + """ nmat = self.nmat if type(freq_range) == int: @@ -958,8 +1134,22 @@ def _find_index_range(self, freq_range): return freq_index + @do_profile @do_profile def read_init(self, info_inputfile): + """ + Read initial configurations for the RPA calculations. + + This method reads the initial configurations for the RPA calculations from + the specified input files. The configurations include chi0q, trans_mod, and + green_init. + + Args: + info_inputfile (dict): Dictionary containing input file information. + + Returns: + dict: Dictionary containing the read configurations. + """ logger.info("RPA read initial configs") info = {} @@ -980,6 +1170,21 @@ def read_init(self, info_inputfile): return info def _read_chi0q(self, file_name): + """ + Read chi0q data from a file. + + This method reads the chi0q data from the specified file and checks its size + and shape based on the calculation scheme. + + Args: + file_name (str): Name of the file to read the chi0q data from. + + Returns: + numpy.ndarray: The read chi0q data. + + Raises: + SystemExit: If reading the file fails or the data shape is unexpected. + """ logger.debug(">>> RPA._read_chi0q") try: @@ -1060,6 +1265,15 @@ def _read_chi0q(self, file_name): return chi0q def _read_trans_mod(self, file_name): + """ + Read the transformation module from a file and perform FFT. + + Args: + file_name (str): Name of the file to read the transformation module from. + + Returns: + numpy.ndarray: The transformed module in k-space. + """ logger.debug(">>> RPA._read_trans_mod") try: @@ -1084,6 +1298,15 @@ def _read_trans_mod(self, file_name): return tab_k def _read_green(self, file_name): + """ + Read the Green's function from a file. + + Args: + file_name (str): Name of the file to read the Green's function from. + + Returns: + numpy.ndarray: The Green's function reshaped to the appropriate dimensions. + """ logger.debug(">>> RPA._read_green") try: logger.info("read green from {}".format(file_name)) @@ -1097,10 +1320,17 @@ def _read_green(self, file_name): nvol = self.lattice.nvol nd = self.nd return green.reshape(nvol,nd,nd) - def _reshape_green(self, green_): - # convert green function into sublattice + def _reshape_green(self, green_): + """ + Convert the Green's function into sublattice format. + Args: + green_ (numpy.ndarray): The original Green's function. + + Returns: + numpy.ndarray: The Green's function reshaped for the sublattice. + """ Lx,Ly,Lz = self.lattice.cellshape Lvol = Lx * Ly * Lz Bx,By,Bz = self.lattice.subshape @@ -1153,6 +1383,15 @@ def _unpack_index(x, n): return green_sub def _calc_trans_mod(self, g0): + """ + Calculate the transformation module. + + Args: + g0 (numpy.ndarray): The Green's function. + + Returns: + numpy.ndarray: The transformed Hamiltonian in k-space. + """ logger.debug(">>> RPA._calc_trans_mod") nx,ny,nz = self.lattice.shape @@ -1180,6 +1419,19 @@ def _calc_trans_mod(self, g0): @do_profile def _calc_epsilon_k(self, green_info): + """ + Calculate the epsilon\_k values for the RPA calculations. + + This method calculates the epsilon\_k values based on the provided Green's function information. + It determines the transfer term H0(k) and performs necessary transformations based on the spin-orbital + interactions and external fields. + + Args: + green_info (dict): Dictionary containing Green's function information. + + Returns: + None + """ logger.debug(">>> RPA._calc_epsilon_k") nvol = self.lattice.nvol @@ -1262,6 +1514,20 @@ def _calc_epsilon_k(self, green_info): @do_profile def _find_mu(self, Ncond, T): + """ + Find the chemical potential (mu) for the given number of electrons and temperature. + + This method calculates the chemical potential (mu) such that the number of electrons + matches the specified value (Ncond) at the given temperature (T). It uses the eigenvalues + and eigenvectors of the Hamiltonian to perform the calculation. + + Args: + Ncond (float): The number of electrons. + T (float): The temperature. + + Returns: + tuple: A tuple containing the distribution (dist) and the chemical potential (mu). + """ logger.debug(">>> RPA._find_mu") from scipy import optimize @@ -1278,6 +1544,17 @@ def _find_mu(self, Ncond, T): occupied_number = Ncond def _fermi(t, mu, ev): + """ + Calculate the Fermi-Dirac distribution. + + Args: + t (float): Temperature. + mu (float): Chemical potential. + ev (numpy.ndarray): Eigenvalues. + + Returns: + numpy.ndarray: Fermi-Dirac distribution values. + """ w_ = (ev - mu) / t mask_ = w_ < ene_cutoff w1_ = np.where( mask_, w_, 0.0 ) @@ -1286,11 +1563,20 @@ def _fermi(t, mu, ev): return v_ def _calc_delta_n(mu): + """ + Calculate the difference between the number of electrons and the target number. + + Args: + mu (float): Chemical potential. + + Returns: + float: Difference between the calculated and target number of electrons. + """ ff = _fermi(T, mu, w) nn = np.einsum('gkal,gkl,gkal->', np.conjugate(v), ff, v) return nn.real - occupied_number - # find mu s.t. (mu) = N0 + # Find mu such that (mu) = N0 is_converged = False if (_calc_delta_n(ev[0]) * _calc_delta_n(ev[-1])) < 0.0: logger.debug("RPA._find_mu: try bisection") @@ -1313,10 +1599,14 @@ def _calc_delta_n(mu): @do_profile def _calc_green(self, beta, mu): """ - ew: eigenvalues ew[g,k,i] i-th eigenvalue of wavenumber k in block g - ev: eigenvectors ev[g,k,a,i] i-th eigenvector corresponding to ew[g,k,i] - beta: inverse temperature - mu: chemical potential + Calculate the Green's function for the RPA calculations. + + Args: + beta (float): Inverse temperature. + mu (float): Chemical potential. + + Returns: + tuple: A tuple containing the Green's function and the tail correction. """ logger.debug(">>> RPA._calc_green") @@ -1347,15 +1637,19 @@ def _calc_green(self, beta, mu): green_tail = np.einsum('gkaj,gkbj,gl->glkab', ev, np.conj(ev), np.tile(np.ones(nmat), (nblock,1))) * aa * 0.5 * beta return green, green_tail - @do_profile def _calc_chi0q(self, green_kw, green0_tail, beta): - # green function - # green_kw[g,l,k,a,b]: G_ab(k,ie) in block g - # l : matsubara freq i\epsilon_l = i\pi(2*l+1-N)/\beta for l=0..N-1 - # k : wave number kx,ky,kz - # a,b : orbital and spin index a=(s,alpha) b=(t,beta) + """ + Calculate the chi0q values for the RPA calculations. + + Args: + green_kw (numpy.ndarray): Green's function in Matsubara frequency space. + green0_tail (numpy.ndarray): Tail correction for the Green's function. + beta (float): Inverse temperature. + Returns: + numpy.ndarray: The chi0q values in Matsubara frequency space. + """ logger.debug(">>> RPA._calc_chi0q") nx,ny,nz = self.lattice.shape @@ -1409,9 +1703,18 @@ def _calc_chi0q(self, green_kw, green0_tail, beta): axis=1).reshape(nblock,nmat,nvol,*nd_shape) * (-1.0/beta) return chi0_qw - @do_profile def _solve_rpa(self, chi0q, ham): + """ + Solve the RPA equations. + + Args: + chi0q (numpy.ndarray): The chi0q values. + ham (numpy.ndarray): The Hamiltonian in k-space. + + Returns: + numpy.ndarray: The solution to the RPA equations. + """ logger.debug(">>> RPA._solve_rpa") nvol = self.lattice.nvol @@ -1435,6 +1738,16 @@ def _solve_rpa(self, chi0q, ham): def run(*, input_dict: Optional[dict] = None, input_file: Optional[str] = None): + """ + Run the RPA calculations based on the provided input dictionary and file. + + Args: + input_dict (Optional[dict]): Dictionary containing input parameters. + input_file (Optional[str]): Path to the input file (not used in the current implementation). + + Raises: + RuntimeError: If input_dict is not provided. + """ logger = logging.getLogger("run") if input_dict is None: diff --git a/src/hwave/solver/uhfk.py b/src/hwave/solver/uhfk.py index ffbbdc7..fc30929 100644 --- a/src/hwave/solver/uhfk.py +++ b/src/hwave/solver/uhfk.py @@ -10,6 +10,15 @@ class UHFk(solver_base): @do_profile def __init__(self, param_ham, info_log, info_mode, param_mod=None): + """ + Initialize the UHFk solver. + + Parameters: + param_ham (dict): Hamiltonian parameters. + info_log (str): Log information. + info_mode (dict): Mode information. + param_mod (dict, optional): Additional parameters. Defaults to None. + """ self.name = "uhfk" super().__init__(param_ham, info_log, info_mode, param_mod) @@ -47,11 +56,20 @@ def __init__(self, param_ham, info_log, info_mode, param_mod=None): self.ham = np.zeros((nvol,nd,nd), dtype=np.complex128) def _init_mode(self, param): + """ + Initialize mode parameters. + + Parameters: + param (dict): Mode parameters. + """ self.iflag_fock = param.get('flag_fock', True) self.enable_spin_orbital = param.get('enable_spin_orbital', False) @do_profile def _init_param(self): + """ + Initialize and check parameters. + """ # check and store parameters # - cell shape @@ -99,6 +117,9 @@ def _init_param(self): @do_profile def _init_lattice(self): + """ + Initialize lattice parameters. + """ Lx,Ly,Lz = self.param_mod.get("CellShape") self.cellshape = (Lx,Ly,Lz) self.cellvol = Lx * Ly * Lz @@ -128,6 +149,9 @@ def _init_lattice(self): @do_profile def _init_orbit(self): + """ + Initialize orbital parameters. + """ norb = self.param_ham["Geometry"]["norb"] ns = 2 # spin dof @@ -140,6 +164,9 @@ def _init_orbit(self): @do_profile def _show_param(self): + """ + Log the parameters. + """ logger.info("Show parameters") logger.info(" Enable Fock = {}".format(self.iflag_fock)) logger.info(" Cell Shape = {}".format(self.cellshape)) @@ -169,6 +196,15 @@ def _show_param(self): @do_profile def _reshape_geometry(self, geom): + """ + Reshape the geometry for the supercell. + + Parameters: + geom (dict): Geometry information. + + Returns: + dict: Reshaped geometry. + """ Bx,By,Bz = self.subshape bvol = Bx * By * Bz @@ -193,6 +229,16 @@ def _reshape_geometry(self, geom): @do_profile def _reshape_interaction(self, ham, enable_spin_orbital): + """ + Reshape the interaction Hamiltonian for the supercell. + + Parameters: + ham (dict): Interaction Hamiltonian. + enable_spin_orbital (bool): Flag to enable spin-orbital coupling. + + Returns: + dict: Reshaped interaction Hamiltonian. + """ Bx,By,Bz = self.subshape nx,ny,nz = self.shape @@ -243,11 +289,18 @@ def _round(x, n): ham_new[(ir, ov)] = v return ham_new - + @do_profile @do_profile def _reshape_green(self, green): - # convert green function into sublattice + """ + Convert the Green's function into sublattice. + + Parameters: + green (numpy.ndarray): Original Green's function. + Returns: + numpy.ndarray: Reshaped Green's function for the sublattice. + """ Lx,Ly,Lz = self.cellshape Lvol = Lx * Ly * Lz Bx,By,Bz = self.subshape @@ -312,7 +365,15 @@ def _unpack_site(idx, n): @do_profile def _deflate_green(self, green_sub): - # convert green function back to original lattice + """ + Convert the Green's function back to the original lattice. + + Parameters: + green_sub (numpy.ndarray): Green's function for the sublattice. + + Returns: + numpy.ndarray: Green's function for the original lattice. + """ Lx,Ly,Lz = self.cellshape Lvol = Lx * Ly * Lz Bx,By,Bz = self.subshape @@ -370,6 +431,10 @@ def _unpack_site(idx, n): @do_profile def _init_interaction(self): + """ + Initialize interaction coefficients on the sublattice. + If the system has a sublattice, the interaction coefficients are reinterpreted. + """ # reinterpret interaction coefficient on sublattice if self.has_sublattice: # backup @@ -389,8 +454,10 @@ def _init_interaction(self): else: tbl = self._reshape_interaction(self.param_ham[type], False) self.param_ham[type] = tbl - def _init_wavevec(self): + """ + Initialize wave vectors on the sublatticed geometry. + """ # wave vectors on sublatticed geometry def _klist(n): return np.roll( (np.arange(n)-(n//2)), -(n//2) ) @@ -416,9 +483,11 @@ def _klist(n): v = kvec[0] * kx + kvec[1] * ky + kvec[2] * kz wtable[ix,iy,iz] = v self.wave_table = wtable.reshape(nvol,3) - @do_profile def _dump_param_ham(self): + """ + Dump the parameters of the Hamiltonian. + """ if logger.getEffectiveLevel() < logging.INFO: return @@ -433,14 +502,18 @@ def _dump_param_ham(self): for (irvec,orbvec), v in self.param_ham[type].items(): print("\t",irvec,orbvec," = ",v) - @do_profile def _check_interaction(self): + """ + Check the validity of the interaction parameters. + """ self._check_cellsize() self._check_orbital_index() self._check_hermite() - def _check_cellsize(self): + """ + Check if the interaction range exceeds the cell shape. + """ err = 0 for k in self.param_ham.keys(): if k in ["Geometry", "Initial"]: @@ -478,8 +551,10 @@ def _check_cellsize(self): else: logger.error(msg) exit(1) - def _check_orbital_index(self): + """ + Check if the orbital indices are valid. + """ norb = self.param_ham["Geometry"]["norb"] err = 0 for k in self.param_ham.keys(): @@ -524,6 +599,16 @@ def _check_orbital_index(self): exit(1) def _check_hermite(self): + """ + Check the Hermiticity of the transfer parameters. + + This function verifies that the transfer parameters in the Hamiltonian + are Hermitian. If the parameters are not Hermitian, it logs an error + message and exits the program if strict Hermiticity is enforced. + + Attributes: + max_print (int): Maximum number of errors to print. + """ max_print = 32 if "Transfer" in self.param_ham: @@ -557,6 +642,16 @@ def _check_hermite(self): @do_profile def solve(self, green_info, path_to_output): + """ + Solve the UHFk equations. + + This function performs the main loop of the UHFk calculations, iterating + until convergence is achieved or the maximum number of iterations is reached. + + Parameters: + green_info (dict): Information about the Green's function. + path_to_output (str): Path to the output directory. + """ print_level = self.info_log["print_level"] print_check = self.info_log.get("print_check", None) @@ -612,6 +707,19 @@ def solve(self, green_info, path_to_output): @do_profile def _initial_green(self, green_info): + """ + Initialize the Green's function. + + This function initializes the Green's function based on the provided + information. It can load the Green's function from a file, initialize it + in coordinate space, or generate it randomly. + + Parameters: + green_info (dict): Information about the Green's function. + + Returns: + numpy.ndarray: The initialized Green's function. + """ logger.debug(">>> _initial_green") _data = None @@ -632,6 +740,12 @@ def _initial_green(self, green_info): return _data def _initial_green_random_reshape(self): + """ + Initialize the Green's function with random numbers and reshape it for the sublattice. + + Returns: + numpy.ndarray: The initialized and reshaped Green's function. + """ lx,ly,lz = self.cellshape lvol = self.cellvol norb = self.norb_orig if self.has_sublattice else self.norb @@ -639,10 +753,8 @@ def _initial_green_random_reshape(self): nd = ns * norb np.random.seed(self.param_mod["RndSeed"]) - #rand = np.random.rand(nvol * nd * nd).reshape(nvol,ns,norb,ns,norb) - #green = 0.01 * (rand - 0.5) - # G[(s,i,a),(t,j,b)] -> G[ij,s,a,t,b] + # Generate random Green's function and reshape x1 = np.random.rand(lvol * nd * lvol * nd).reshape(ns,lvol,norb,ns,lvol,norb) x2 = np.transpose(x1,(1,4,0,2,3,5)) x3 = np.average(x2,axis=0).reshape(lvol,ns,norb,ns,norb) @@ -655,16 +767,20 @@ def _initial_green_random_reshape(self): return green def _initial_green_random(self): + """ + Initialize the Green's function with random numbers. + + Returns: + numpy.ndarray: The initialized Green's function. + """ nvol = self.nvol norb = self.norb ns = self.ns nd = ns * norb np.random.seed(self.param_mod["RndSeed"]) - #rand = np.random.rand(nvol * nd * nd).reshape(nvol,ns,norb,ns,norb) - #green = 0.01 * (rand - 0.5) - # G[(s,i,a),(t,j,b)] -> G[ij,s,a,t,b] + # Generate random Green's function and reshape x1 = np.random.rand(nvol * nd * nvol * nd).reshape(ns,nvol,norb,ns,nvol,norb) x2 = np.transpose(x1,(1,4,0,2,3,5)) x3 = np.average(x2,axis=0).reshape(nvol,ns,norb,ns,norb) @@ -674,6 +790,16 @@ def _initial_green_random(self): return green def _initial_green_uhf(self, info, geom): + """ + Initialize the Green's function from UHF data. + + Parameters: + info (dict): UHF Green's function data. + geom (dict): Geometry information. + + Returns: + numpy.ndarray: The initialized Green's function. + """ lx,ly,lz = self.cellshape lvol = self.cellvol norb = self.norb @@ -701,6 +827,15 @@ def _pack_site(ix,iy,iz): @do_profile def _make_ham_trans(self): + """ + Construct the Hamiltonian in the momentum space. + + This function constructs the Hamiltonian in the momentum space by performing + a Fourier transform on the real-space Hamiltonian parameters. + + Returns: + None + """ logger.debug(">>> _make_ham_trans") nx,ny,nz = self.shape @@ -708,28 +843,28 @@ def _make_ham_trans(self): norb = self.norb nd = self.nd - if self.enable_spin_orbital == True: + if self.enable_spin_orbital: tab_r = np.zeros((nx,ny,nz,nd,nd), dtype=np.complex128) for (irvec,orbvec), v in self.param_ham["Transfer"].items(): tab_r[(*irvec, *orbvec)] = v - # fourier transform + # Fourier transform tab_k = np.fft.ifftn(tab_r, axes=(0,1,2), norm='forward') ham = tab_k else: - # data structure of T_ab(r): T(rx,ry,rz,a,b) + # Data structure of T_ab(r): T(rx, ry, rz, a, b) tab_r = np.zeros((nx,ny,nz,norb,norb), dtype=np.complex128) for (irvec,orbvec), v in self.param_ham["Transfer"].items(): if orbvec[0] < norb and orbvec[1] < norb: tab_r[(*irvec, *orbvec)] = v else: - pass # skip spin dependence + pass # Skip spin dependence - # fourier transform + # Fourier transform tab_k = np.fft.ifftn(tab_r, axes=(0,1,2), norm='forward') # 2x2 unit matrix to introduce spin index @@ -740,11 +875,18 @@ def _make_ham_trans(self): 'kab, st -> ksatb', tab_k.reshape(nvol,norb,norb), spin ).reshape(nvol,nd,nd) - # store + # Store self.ham_trans = ham @do_profile def _make_ham_inter(self): + """ + Construct the interaction Hamiltonian. + + This function constructs the interaction Hamiltonian by reshaping the + interaction parameters for the sublattice. It also initializes the + interaction and spin tables. + """ logger.debug(">>> _make_ham_inter") nx,ny,nz = self.shape @@ -985,6 +1127,13 @@ def _make_ham_inter(self): @do_profile def _make_ham(self): + """ + Construct the Hamiltonian matrix. + + This function constructs the Hamiltonian matrix by combining the + transfer and interaction Hamiltonians. It uses the Green's function + to update the Hamiltonian matrix elements. + """ logger.debug(">>> _make_ham") nx,ny,nz = self.shape @@ -1066,6 +1215,12 @@ def _make_ham(self): @do_profile def _diag(self): + """ + Diagonalize the Hamiltonian matrix. + + This function reshapes the Hamiltonian matrix and performs + eigenvalue decomposition to obtain the eigenvalues and eigenvectors. + """ logger.debug(">>> _diag") nx,ny,nz = self.shape @@ -1091,6 +1246,11 @@ def _diag(self): @do_profile def _green(self): + """ + Update the Green's function. + + This function updates the Green's function based on the current Hamiltonian. + """ logger.debug(">>> _green") nx,ny,nz = self.shape @@ -1132,13 +1292,44 @@ def _green(self): self.Green = gab_r.reshape(nvol,ns,norb,ns,norb) def _find_dist(self, w, v, ncond): + """ + Determine the distribution of eigenvalues. + + Parameters: + w (numpy.ndarray): Eigenvalues. + v (numpy.ndarray): Eigenvectors. + ncond (int): Number of occupied states. + + Returns: + tuple: Distribution of eigenvalues and chemical potential. + """ if self.T == 0: return self._find_dist_zero_t(w, v, ncond) else: return self._find_dist_nonzero_t(w, v, ncond) def _find_dist_zero_t(self, w, v, ncond): + """ + Determine the distribution of eigenvalues at zero temperature. + + Parameters: + w (numpy.ndarray): Eigenvalues. + v (numpy.ndarray): Eigenvectors. + ncond (int): Number of occupied states. + + Returns: + tuple: Distribution of eigenvalues and chemical potential. + """ def _ksq_table(width): + """ + Generate a table of squared wave vectors. + + Parameters: + width (int): Width of the table. + + Returns: + numpy.ndarray: Table of squared wave vectors. + """ nx,ny,nz = self.shape nvol = self.nvol @@ -1164,12 +1355,34 @@ def _ksq_table(width): return dist.reshape(w.shape), 0.0 def _find_dist_nonzero_t(self, w, v, ncond): + """ + Determine the distribution of eigenvalues at non-zero temperature. + + Parameters: + w (numpy.ndarray): Eigenvalues. + v (numpy.ndarray): Eigenvectors. + ncond (int): Number of occupied states. + + Returns: + tuple: Distribution of eigenvalues and chemical potential. + """ from scipy import optimize ev = np.sort(w.flatten()) occupied_number = ncond def _fermi(t, mu, ev): + """ + Fermi-Dirac distribution function. + + Parameters: + t (float): Temperature. + mu (float): Chemical potential. + ev (numpy.ndarray): Eigenvalues. + + Returns: + numpy.ndarray: Fermi-Dirac distribution values. + """ w_ = (ev - mu) / t mask_ = w_ < self.ene_cutoff w1_ = np.where( mask_, w_, 0.0 ) @@ -1178,11 +1391,20 @@ def _fermi(t, mu, ev): return v_ def _calc_delta_n(mu): + """ + Calculate the difference in the number of occupied states. + + Parameters: + mu (float): Chemical potential. + + Returns: + float: Difference in the number of occupied states. + """ ff = _fermi(self.T, mu, w) nn = np.einsum('kal,kl,kal->', np.conjugate(v), ff, v) return nn.real - occupied_number - # find mu s.t. (mu) = N0 + # Find mu such that (mu) = N0 is_converged = False if (_calc_delta_n(ev[0]) * _calc_delta_n(ev[-1])) < 0.0: logger.debug("+++ find mu: try bisection") @@ -1198,9 +1420,15 @@ def _calc_delta_n(mu): dist = _fermi(self.T, mu, w) return dist, mu - @do_profile def _calc_phys(self): + """ + Calculate physical quantities. + + This function calculates the expectation values of the number of + occupied states (NCond) and the spin projection (Sz). It also + calculates the residue and updates the Green's function. + """ logger.debug(">>> _calc_phys") nx,ny,nz = self.shape @@ -1209,39 +1437,43 @@ def _calc_phys(self): norb = self.norb ns = self.ns - # expectation value of Ncond + # Expectation value of NCond gab_r = self.Green.reshape(nvol,nd,nd) - n = np.sum(np.diagonal(gab_r[0])) * nvol self.physics["NCond"] = n.real - logger.debug("ncond = {}".format(n)) - # expectation value of Sz + # Expectation value of Sz gab_r = self.Green - sigma_z = np.array(np.array([[1,0],[0,-1]])) + sigma_z = np.array([[1, 0], [0, -1]]) sz = np.sum(np.diagonal(np.einsum('satb,st->ab', gab_r[0], sigma_z))) * nvol self.physics["Sz"] = 0.5 * sz.real - logger.debug("sz = {}".format(sz)) - # residue + # Residue rest = np.linalg.norm(self.Green - self.Green_prev) self.physics["Rest"] = rest / np.size(self.Green) * 2 - logger.debug("rest = {}".format(rest)) - # update + # Update mix = self.param_mod["Mix"] - g_new = (1.0 - mix) * self.Green_prev + mix * self.Green - g_new[np.where(abs(g_new) < self.threshold)] = 0.0 - self.Green = g_new @do_profile def _calc_energy(self): + """ + Calculate the total energy of the system. + + This function calculates the total energy of the system by summing the band energy + and interaction energies. It handles both zero and non-zero temperature cases. + + The results are stored in the `self.physics["Ene"]` dictionary. + + Returns: + None + """ logger.debug(">>> _calc_energy") nx,ny,nz = self.shape @@ -1278,12 +1510,22 @@ def _calc_energy(self): mu = mus[k] def _fermi(t, mu, ev): + """ + Fermi-Dirac distribution function. + + Args: + t (float): Temperature. + mu (float): Chemical potential. + ev (numpy.ndarray): Eigenvalues. + + Returns: + numpy.ndarray: Fermi-Dirac distribution values. + """ w = (ev - mu) / t mask_ = w < self.ene_cutoff w1 = np.where( mask_, w, 0.0 ) v1 = 1.0 / (1.0 + np.exp(w1)) v = np.where( mask_, v1, 0.0 ) - #v = np.where( w > self.ene_cutoff, 0.0, 1.0 / (1.0 + np.exp(w)) ) return v wt = -(w - mu) / T @@ -1329,7 +1571,6 @@ def _fermi(t, mu, ev): energy_total += energy[type].real logger.debug("energy: {} = {}".format(type, energy[type])) else: - # logger.info("energy: {} skip".format(type)) pass energy["Total"] = energy_total @@ -1337,10 +1578,31 @@ def _fermi(t, mu, ev): @do_profile def get_results(self): + """ + Get the results of the calculations. + + This function returns the physical quantities and the Green's function. + + Returns: + tuple: A tuple containing the physics dictionary and the Green's function. + """ return (self.physics, self.Green) @do_profile def save_results(self, info_outputfile, green_info): + """ + Save the results to files. + + This function saves the calculated energies, eigenvalues, eigenvectors, and Green's function + to the specified output files. + + Args: + info_outputfile (dict): Information about the output files. + green_info (dict): Information about the Green's function. + + Returns: + None + """ path_to_output = info_outputfile["path_to_output"] if "energy" in info_outputfile.keys(): @@ -1364,7 +1626,6 @@ def save_results(self, info_outputfile, green_info): logger.info("save_results: save energy in file {}".format(file_name)) if "eigen" in info_outputfile.keys(): - # eigenvalue[spin_block,k,eigen_index] -> [k,spin_block*eigen_index] eg = self._green_list["eigenvalue"] egs = eg.shape @@ -1374,11 +1635,6 @@ def save_results(self, info_outputfile, green_info): evs = ev.shape evv = np.einsum('skab,st->ksatb', ev, np.eye(evs[0])).reshape(evs[1],evs[0]*evs[2],evs[0]*evs[3]) - # # wavevec[k,eigen_index] = \vec(k) - # wv = self.wave_table - # wvs = wv.shape - # wvv = np.transpose(np.broadcast_to(wv, ((egs[0]*egs[2]),wvs[0],wvs[1])), (1,0,2)) - file_name = os.path.join(path_to_output, info_outputfile["eigen"]) np.savez(file_name, eigenvalue = egg, @@ -1413,14 +1669,31 @@ def save_results(self, info_outputfile, green_info): @do_profile def _read_green(self, file_name): + """ + Read the Green's function from a file. + + Parameters: + file_name (str): The name of the file to read the Green's function from. + + Returns: + numpy.ndarray or None: The Green's function data if the file is found, otherwise None. + """ try: v = np.load(file_name) except FileNotFoundError: logger.info("_read_green: file {} not found".format(file_name)) return None return self._read_green_from_data(v) - def _read_green_from_data(self, ginfo): + """ + Read the Green's function from the provided data. + + Parameters: + ginfo (numpy.lib.npyio.NpzFile): The data containing the Green's function. + + Returns: + numpy.ndarray or None: The Green's function data if found, otherwise None. + """ if self.has_sublattice: if "green_sublattice" in ginfo.files: logger.debug("_read_green: read green_sublattice") @@ -1439,9 +1712,14 @@ def _read_green_from_data(self, ginfo): logger.debug("_read_green: no data found") data = None return data - @do_profile def _save_green(self, file_name): + """ + Save the Green's function to a file. + + Parameters: + file_name (str): The name of the file to save the Green's function to. + """ if self.has_sublattice: green_orig = self._deflate_green(self.Green) np.savez(file_name, green = green_orig, green_sublattice = self.Green) @@ -1451,6 +1729,19 @@ def _save_green(self, file_name): @do_profile def _save_greenone(self, file_name, green_info): + """ + Save the one-body Green's function to a file. + + This function saves the one-body Green's function to a specified file. It requires + the `onebodyg_uhf` and `geometry_uhf` keys in the `green_info` dictionary. + + Parameters: + file_name (str): The name of the file to save the Green's function to. + green_info (dict): Information about the Green's function, including `onebodyg_uhf` and `geometry_uhf`. + + Returns: + None + """ if not "onebodyg_uhf" in green_info or not "geometry_uhf" in green_info: logger.error("_save_greenone: onebodyg_uhf and geometry_uhf are required") return None @@ -1492,8 +1783,20 @@ def _pack_site(ix,iy,iz): )) logger.info("save_results: save greenone to file {}".format(file_name)) - def _save_trans_mod(self, file_name): + """ + Save the transformed Hamiltonian to a file. + + This function saves the transformed Hamiltonian to a specified file. It performs + a Fourier transform on the Hamiltonian matrix and saves the result. If the system + has a sublattice, it also saves the deflated Green's function. + + Parameters: + file_name (str): The name of the file to save the transformed Hamiltonian to. + + Returns: + None + """ nx,ny,nz = self.shape nvol = self.nvol nd = self.nd @@ -1519,6 +1822,12 @@ def _save_trans_mod(self, file_name): logger.info("save_results: save trans_mod to file {}".format(file_name)) def _export_geometry(self, file_name): + """ + Export the geometry information to a file. + + Parameters: + file_name (str): The name of the file to save the geometry information to. + """ geom = self.param_ham["Geometry"] with open(file_name, "w") as fw: # write unit vector @@ -1533,6 +1842,13 @@ def _export_geometry(self, file_name): fw.write("{:.8f} {:.8f} {:.8f}\n".format(center[k][0], center[k][1], center[k][2])) def _export_interaction(self, type, file_name): + """ + Export the interaction parameters to a file. + + Parameters: + type (str): The type of interaction to export. + file_name (str): The name of the file to save the interaction parameters to. + """ intr = self.param_ham[type] min_r = [0,0,0] @@ -1563,8 +1879,14 @@ def _export_interaction(self, type, file_name): fw.write("{:3} {:3} {:3} {:3} {:3} {:.12f} {:.12f}\n".format( *irvec, *orbvec, v.real, v.imag )) - def _export_hamiltonian(self, path_to_output, prefix): + """ + Export the Hamiltonian parameters to files. + + Parameters: + path_to_output (str): The directory path to save the output files. + prefix (str): The prefix for the output file names. + """ for type in self.param_ham.keys(): if type in ["Geometry"]: file_name = os.path.join(path_to_output, prefix + type + ".dat") diff --git a/src/hwave/solver/uhfr.py b/src/hwave/solver/uhfr.py index a2db826..bdf0f52 100644 --- a/src/hwave/solver/uhfr.py +++ b/src/hwave/solver/uhfr.py @@ -9,18 +9,40 @@ logger = logging.getLogger("qlms").getChild("uhfr") class Interact_UHFr_base(): + """ + Base class for interaction Hamiltonians in the UHFr solver. + """ + def __init__(self, ham_info, Nsize): + """ + Initialize the interaction Hamiltonian. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + """ self.Nsize = Nsize self.Ham_tmp = np.zeros(tuple([(2 * self.Nsize) for i in range(4)]), dtype=complex) self.Ham_trans_tmp = np.zeros(tuple([(2 * self.Nsize) for i in range(2)]), dtype=complex) self.param_ham = self._transform_interall(ham_info) self._check_range() - #Change interaction to interall type def _transform_interall(self, ham_info): + """ + Transform interaction to interall type. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ return ham_info def _calc_hartree(self): + """ + Calculate the Hartree term of the Hamiltonian. + """ site = np.zeros(4, dtype=np.int32) for site_info, value in self.param_ham.items(): for i in range(4): @@ -33,7 +55,10 @@ def _calc_hartree(self): pass def _calc_fock(self): - site = np.zeros(4,dtype=np.int32) + """ + Calculate the Fock term of the Hamiltonian. + """ + site = np.zeros(4, dtype=np.int32) for site_info, value in self.param_ham.items(): for i in range(4): site[i] = site_info[2 * i] + site_info[2 * i + 1] * self.Nsize @@ -43,13 +68,24 @@ def _calc_fock(self): pass def get_ham(self, type): + """ + Get the Hamiltonian. + + Args: + type (str): Type of Hamiltonian to calculate ("hartree" or "hartreefock"). + + Returns: + tuple: Hamiltonian and transformed Hamiltonian. + """ self._calc_hartree() if type == "hartreefock": self._calc_fock() return self.Ham_tmp, self.Ham_trans_tmp - # check input def _check_range(self): + """ + Check the range of the Hamiltonian parameters. + """ err = 0 for site_info, value in self.param_ham.items(): for i in range(4): @@ -62,6 +98,13 @@ def _check_range(self): exit(1) def _check_hermite(self, strict_hermite, tolerance): + """ + Check the Hermiticity of the Hamiltonian. + + Args: + strict_hermite (bool): Whether to enforce strict Hermiticity. + tolerance (float): Tolerance for Hermiticity check. + """ err = 0 for site_info, value in self.param_ham.items(): list = tuple([site_info[i] for i in [6,7,4,5,2,3,0,1]]) @@ -78,13 +121,39 @@ def _check_hermite(self, strict_hermite, tolerance): logger.warn(msg) def check_hermite(self, strict_hermite=False, tolerance=1.0e-8): - self._check_hermite(strict_hermite, tolerance) + """ + Public method to check the Hermiticity of the Hamiltonian. + Args: + strict_hermite (bool): Whether to enforce strict Hermiticity. + tolerance (float): Tolerance for Hermiticity check. + """ + self._check_hermite(strict_hermite, tolerance) class CoulombIntra_UHFr(Interact_UHFr_base): + """ + Class for Coulomb Intra interaction in the UHFr solver. + """ def __init__(self, ham_info, Nsize): + """ + Initialize the CoulombIntra_UHFr class. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + """ self.__name__ = "CoulombIntra" super().__init__(ham_info, Nsize) + def _transform_interall(self, ham_info): + """ + Transform interaction to interall type for Coulomb Intra. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ param_tmp = {} for site_info, value in ham_info.items(): sinfo = tuple([site_info[0], 0, site_info[0], 0, site_info[0], 1, site_info[0], 1]) @@ -92,10 +161,30 @@ def _transform_interall(self, ham_info): return param_tmp class CoulombInter_UHFr(Interact_UHFr_base): + """ + Class for Coulomb Inter interaction in the UHFr solver. + """ def __init__(self, ham_info, Nsize): + """ + Initialize the CoulombInter_UHFr class. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + """ self.__name__ = "CoulombInter" super().__init__(ham_info, Nsize) + def _transform_interall(self, ham_info): + """ + Transform interaction to interall type for Coulomb Inter. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ param_tmp = {} for site_info, value in ham_info.items(): for spin_i, spin_j in itertools.product([0,1], repeat=2): @@ -104,10 +193,30 @@ def _transform_interall(self, ham_info): return param_tmp class Hund_UHFr(Interact_UHFr_base): + """ + Class for Hund interaction in the UHFr solver. + """ def __init__(self, ham_info, Nsize): + """ + Initialize the Hund_UHFr class. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + """ self.__name__ = "Hund" super().__init__(ham_info, Nsize) + def _transform_interall(self, ham_info): + """ + Transform interaction to interall type for Hund. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ param_tmp = {} for site_info, value in ham_info.items(): for spin_i in range(2): @@ -116,10 +225,30 @@ def _transform_interall(self, ham_info): return param_tmp class PairHop_UHFr(Interact_UHFr_base): + """ + Class for Pair Hop interaction in the UHFr solver. + """ def __init__(self, ham_info, Nsize): + """ + Initialize the PairHop_UHFr class. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + """ self.__name__ = "PairHop" super().__init__(ham_info, Nsize) + def _transform_interall(self, ham_info): + """ + Transform interaction to interall type for Pair Hop. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ param_tmp = {} for site_info, value in ham_info.items(): sinfo = tuple([site_info[0], 0, site_info[1], 0, site_info[0], 1, site_info[1], 1]) @@ -129,10 +258,30 @@ def _transform_interall(self, ham_info): return param_tmp class Exchange_UHFr(Interact_UHFr_base): + """ + Class for Exchange interaction in the UHFr solver. + """ def __init__(self, ham_info, Nsize): + """ + Initialize the Exchange_UHFr class. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + """ self.__name__ = "Exchange" super().__init__(ham_info, Nsize) + def _transform_interall(self, ham_info): + """ + Transform interaction to interall type for Exchange. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ param_tmp = {} for site_info, value in ham_info.items(): sinfo = tuple([site_info[0], 0, site_info[1], 0, site_info[1], 1, site_info[0], 1]) @@ -140,12 +289,31 @@ def _transform_interall(self, ham_info): sinfo = tuple([site_info[0], 1, site_info[1], 1, site_info[1], 0, site_info[0], 0]) param_tmp[sinfo] = value return param_tmp - class Ising_UHFr(Interact_UHFr_base): + """ + Class for Ising interaction in the UHFr solver. + """ def __init__(self, ham_info, Nsize): + """ + Initialize the Ising_UHFr class. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + """ self.__name__ = "Ising" super().__init__(ham_info, Nsize) + def _transform_interall(self, ham_info): + """ + Transform interaction to interall type for Ising. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ param_tmp = {} for site_info, value in ham_info.items(): for spin_i, spin_j in itertools.product([0,1], repeat=2): @@ -154,10 +322,30 @@ def _transform_interall(self, ham_info): return param_tmp class PairLift_UHFr(Interact_UHFr_base): + """ + Class for Pair Lift interaction in the UHFr solver. + """ def __init__(self, ham_info, Nsize): + """ + Initialize the PairLift_UHFr class. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + """ self.__name__ = "PairLift" super().__init__(ham_info, Nsize) + def _transform_interall(self, ham_info): + """ + Transform interaction to interall type for Pair Lift. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ param_tmp = {} for site_info, value in ham_info.items(): sinfo = tuple([site_info[0], 0, site_info[0], 1, site_info[1], 0, site_info[1], 1]) @@ -167,21 +355,61 @@ def _transform_interall(self, ham_info): return param_tmp class InterAll_UHFr(Interact_UHFr_base): + """ + Class for InterAll interaction in the UHFr solver. + """ def __init__(self, ham_info, Nsize, strict_hermite, tolerance): + """ + Initialize the InterAll_UHFr class. + + Args: + ham_info (dict): Hamiltonian information. + Nsize (int): Size of the system. + strict_hermite (bool): Whether to enforce strict Hermiticity. + tolerance (float): Tolerance for Hermiticity check. + """ self.__name__ = "InterAll" super().__init__(ham_info, Nsize) self._check_hermite(strict_hermite, tolerance) + def _transform_interall(self, ham_info): + """ + Transform interaction to interall type for InterAll. + + Args: + ham_info (dict): Hamiltonian information. + + Returns: + dict: Transformed Hamiltonian information. + """ return ham_info class Term_base: + """ + Base class for terms in the UHFr solver. + """ + def __init__(self, term_info, Nsize, coeff=1.0): + """ + Initialize the Term_base class. + + Args: + term_info (dict): Term information. + Nsize (int): Size of the system. + coeff (float): Coefficient for the term. + """ self.Nsize = Nsize self.term_info = term_info self.coeff = coeff self._check_range() def get_data(self): + """ + Get the term data. + + Returns: + np.ndarray: Term data as a complex numpy array. + """ data = np.zeros(tuple([(2 * self.Nsize) for i in range(2)]), dtype=complex) for site_info, value in self.term_info.items(): # set value @@ -191,6 +419,9 @@ def get_data(self): return data def _check_range(self): + """ + Check the range of the term parameters. + """ err = 0 for site_info, value in self.term_info.items(): for i in range(2): @@ -203,6 +434,13 @@ def _check_range(self): exit(1) def _check_hermite(self, strict_hermite, tolerance): + """ + Check the Hermiticity of the term. + + Args: + strict_hermite (bool): Whether to enforce strict Hermiticity. + tolerance (float): Tolerance for Hermiticity check. + """ err = 0 for site_info, value in self.term_info.items(): list = tuple([site_info[i] for i in [2,3,0,1]]) @@ -219,24 +457,62 @@ def _check_hermite(self, strict_hermite, tolerance): logger.warn(msg) def check_hermite(self, strict_hermite=False, tolerance=1.0e-8): - self._check_hermite(strict_hermite, tolerance) + """ + Public method to check the Hermiticity of the term. + Args: + strict_hermite (bool): Whether to enforce strict Hermiticity. + tolerance (float): Tolerance for Hermiticity check. + """ + self._check_hermite(strict_hermite, tolerance) class Transfer_UHFr(Term_base): + """ + Class for Transfer term in the UHFr solver. + """ def __init__(self, term_info, Nsize): + """ + Initialize the Transfer_UHFr class. + + Args: + term_info (dict): Term information. + Nsize (int): Size of the system. + """ self.__name__ = "Transfer" super().__init__(term_info, Nsize, coeff=-1.0) class Green_UHFr(Term_base): + """ + Class for Green term in the UHFr solver. + """ def __init__(self, term_info, Nsize): + """ + Initialize the Green_UHFr class. + + Args: + term_info (dict): Term information. + Nsize (int): Size of the system. + """ self.__name__ = "Green" super().__init__(term_info, Nsize, coeff=1.0) - from .base import solver_base class UHFr(solver_base): + """ + Class for the Unrestricted Hartree-Fock-Roothaan (UHFr) solver. + """ + @do_profile def __init__(self, param_ham, info_log, info_mode, param_mod=None): + """ + Initialize the UHFr solver. + + Args: + param_ham (dict): Hamiltonian parameters. + info_log (dict): Logging information. + info_mode (dict): Mode information. + param_mod (dict, optional): Additional parameters for the solver. + """ self.name = "uhfr" super().__init__(param_ham, info_log, info_mode, param_mod) self.physics = {"Ene": 0, "NCond": 0, "Sz": 0, "Rest": 1.0} @@ -275,6 +551,13 @@ def __init__(self, param_ham, info_log, info_mode, param_mod=None): @do_profile def solve(self, green_info, path_to_output): + """ + Solve the UHFr problem. + + Args: + green_info (dict): Information about the Green's function. + path_to_output (str): Path to the output directory. + """ print_level = self.info_log["print_level"] print_step = self.info_log["print_step"] print_check = self.info_log.get("print_check", None) @@ -327,6 +610,15 @@ def solve(self, green_info, path_to_output): @do_profile def _initial_G(self, green_info): + """ + Initialize the Green's function. + + Args: + green_info (dict): Information about the Green's function. + + Returns: + np.ndarray: Initialized Green's function. + """ _green_list = self.green_list green = np.zeros((2 * self.Nsize, 2 * self.Nsize), dtype=complex) if green_info["Initial"] is not None: @@ -364,6 +656,9 @@ def _initial_G(self, green_info): @do_profile def _makeham_const(self): + """ + Create the constant part of the Hamiltonian. + """ self.Ham_trans = np.zeros((2 * self.Nsize, 2 * self.Nsize), dtype=complex) # Transfer integrals trans = Transfer_UHFr(self.param_ham["Transfer"], self.Nsize) @@ -372,6 +667,9 @@ def _makeham_const(self): @do_profile def _makeham(self): + """ + Create the Hamiltonian matrix. + """ import time self.Ham = np.zeros((2 * self.Nsize, 2 * self.Nsize), dtype=complex) self.Ham = self.Ham_trans.copy() @@ -381,6 +679,9 @@ def _makeham(self): @do_profile def _makeham_mat(self): + """ + Create the Hamiltonian matrix with interaction terms. + """ # TODO Add Hund, Exchange, Ising, PairHop, and PairLift self.Ham_local = np.zeros(tuple([(2 * self.Nsize) for i in range(4)]), dtype=complex) if self.iflag_fock is True: @@ -417,9 +718,19 @@ def _makeham_mat(self): self.Ham_trans += Ham_trans_tmp self.Ham_local = self.Ham_local.reshape((2 * self.Nsize) ** 2, (2 * self.Nsize) ** 2) - @do_profile def _diag(self): + """ + Diagonalize the Hamiltonian matrix and store the eigenvalues and eigenvectors. + + This method iterates over the blocks of the Green's function list, constructs + the Hamiltonian matrix for each block, and then diagonalizes it using numpy's + `linalg.eigh` function. The eigenvalues and eigenvectors are stored in the + `green_list` dictionary. + + Attributes: + _green_list (dict): Dictionary containing the Green's function blocks. + """ _green_list = self.green_list for k, block_g_info in _green_list.items(): g_label = block_g_info["label"] @@ -437,6 +748,16 @@ def _diag(self): @do_profile def _fermi(self, mu, eigenvalue): + """ + Calculate the Fermi-Dirac distribution for a given chemical potential and eigenvalues. + + Args: + mu (float): Chemical potential. + eigenvalue (np.ndarray): Array of eigenvalues. + + Returns: + np.ndarray: Array of Fermi-Dirac distribution values. + """ fermi = np.zeros(eigenvalue.shape) for idx, value in enumerate(eigenvalue): if (value - mu) / self.T > self.ene_cutoff: @@ -444,9 +765,19 @@ def _fermi(self, mu, eigenvalue): else: fermi[idx] = 1.0 / (np.exp((value - mu) / self.T) + 1.0) return fermi - @do_profile def _green(self): + """ + Update the Green's function. + + This method updates the Green's function based on the current temperature. + If the temperature is zero, it uses the eigenvalues and eigenvectors to + construct the Green's function. For finite temperatures, it calculates the + chemical potential and uses the Fermi-Dirac distribution. + + Attributes: + _green_list (dict): Dictionary containing the Green's function blocks. + """ _green_list = self.green_list # R_SLT = U^{*} in _green # L_SLT = U^T in _green @@ -472,9 +803,18 @@ def _green(self): from scipy import optimize def _calc_delta_n(mu): + """ + Calculate the difference in electron number for a given chemical potential. + + Args: + mu (float): Chemical potential. + + Returns: + float: Difference in electron number. + """ n_eigen = np.einsum("ij, ij -> j", np.conjugate(eigenvec), eigenvec).real fermi = self._fermi(mu, eigenvalue) - return np.dot(n_eigen, fermi)-occupied_number + return np.dot(n_eigen, fermi) - occupied_number self.Green = np.zeros((2 * self.Nsize, 2 * self.Nsize), dtype=complex) for k, block_g_info in _green_list.items(): @@ -483,7 +823,7 @@ def _calc_delta_n(mu): eigenvec = self.green_list[k]["eigenvector"] occupied_number = block_g_info["occupied"] - #determine mu + # determine mu is_converged = False if (_calc_delta_n(eigenvalue[0]) * _calc_delta_n(eigenvalue[-1])) < 0.0: @@ -502,9 +842,17 @@ def _calc_delta_n(mu): for idx1, org_site1 in enumerate(g_label): for idx2, org_site2 in enumerate(g_label): self.Green[org_site1][org_site2] += tmp_green[idx1][idx2] - @do_profile def _calc_energy(self): + """ + Calculate the energy of the system. + + This method calculates the band energy and interaction energy of the system. + The results are stored in the `self.physics["Ene"]` dictionary. + + Attributes: + _green_list (dict): Dictionary containing the Green's function blocks. + """ _green_list = self.green_list Ene = {} Ene["band"] = 0 @@ -514,7 +862,6 @@ def _calc_energy(self): occupied_number = block_g_info["occupied"] Ene["band"] += np.sum(eigenvalue[:occupied_number]) else: - for k, block_g_info in _green_list.items(): eigenvalue = self.green_list[k]["eigenvalue"] eigenvec = self.green_list[k]["eigenvector"] @@ -528,11 +875,11 @@ def _calc_energy(self): else: ln_Ene[idx] = -(value - mu) / self.T tmp_n = np.einsum("ij, j, ij -> i", np.conjugate(eigenvec), fermi, eigenvec) - Ene["band"] += mu*np.sum(tmp_n) - self.T * np.sum(ln_Ene) + Ene["band"] += mu * np.sum(tmp_n) - self.T * np.sum(ln_Ene) Ene["InterAll"] = 0 green_local = self.Green.reshape((2 * self.Nsize) ** 2) - Ene["InterAll"] -= np.dot(green_local.T, np.dot(self.Ham_local, green_local))/2.0 + Ene["InterAll"] -= np.dot(green_local.T, np.dot(self.Ham_local, green_local)) / 2.0 ene = 0 for value in Ene.values(): @@ -541,8 +888,16 @@ def _calc_energy(self): self.physics["Ene"] = Ene logger.debug(Ene) + @do_profile def _calc_phys(self): + """ + Calculate physical properties of the system. + + This method calculates the number of conduction electrons (NCond), + the spin polarization (Sz), and the residual error (Rest). + The results are stored in the `self.physics` dictionary. + """ n = 0 for site in range(2 * self.Nsize): n += self.Green[site][site] @@ -557,21 +912,36 @@ def _calc_phys(self): rest = 0.0 mix = self.param_mod["mix"] for site1, site2 in itertools.product(range(2 * self.Nsize), range(2 * self.Nsize)): - rest += abs(self.Green[site1][site2] - self.Green_old[site1][site2])**2 + rest += abs(self.Green[site1][site2] - self.Green_old[site1][site2]) ** 2 self.Green[site1][site2] = self.Green_old[site1][site2] * (1.0 - mix) + mix * self.Green[site1][site2] self.physics["Rest"] = np.sqrt(rest) / (2.0 * self.Nsize * self.Nsize) self.Green[np.where(abs(self.Green) < self.threshold)] = 0 + @do_profile def get_results(self): - return (self.physics, self.Green) + """ + Get the results of the UHFr calculation. + + Returns: + tuple: A tuple containing the physics dictionary and the Green's function. + """ + return self.physics, self.Green + @do_profile def save_results(self, info_outputfile, green_info): + """ + Save the results of the UHFr calculation to files. + + Args: + info_outputfile (dict): Information about the output files. + green_info (dict): Information about the Green's function. + """ path_to_output = info_outputfile["path_to_output"] if "energy" in info_outputfile.keys(): - output_str = "Energy_total = {}\n".format(self.physics["Ene"]["Total"].real) + output_str = "Energy_total = {}\n".format(self.physics["Ene"]["Total"].real) output_str += "Energy_band = {}\n".format(self.physics["Ene"]["band"].real) output_str += "Energy_interall = {}\n".format(self.physics["Ene"]["InterAll"].real) output_str += "NCond = {}\n".format(self.physics["NCond"]) @@ -583,7 +953,7 @@ def save_results(self, info_outputfile, green_info): for key, _green_list in self.green_list.items(): eigenvalue = _green_list["eigenvalue"] eigenvector = _green_list["eigenvector"] - np.savez(os.path.join(path_to_output, key+"_"+info_outputfile["eigen"]), eigenvalue = eigenvalue, eigenvector=eigenvector) + np.savez(os.path.join(path_to_output, key + "_" + info_outputfile["eigen"]), eigenvalue=eigenvalue, eigenvector=eigenvector) if "green" in info_outputfile.keys(): if "OneBodyG" in green_info and green_info["OneBodyG"] is not None: @@ -591,25 +961,25 @@ def save_results(self, info_outputfile, green_info): _green_info = np.array(_green_info, dtype=np.int32) output_str = "" for info in _green_info: - nsite1 = info[0] + self.Nsize*info[1] - nsite2 = info[2] + self.Nsize*info[3] + nsite1 = info[0] + self.Nsize * info[1] + nsite2 = info[2] + self.Nsize * info[3] output_str += "{} {} {} {} {} {}\n".format(info[0], info[1], info[2], info[3], self.Green[nsite1][nsite2].real, self.Green[nsite1][nsite2].imag) with open(os.path.join(path_to_output, info_outputfile["green"]), "w") as fw: fw.write(output_str) else: logger.error("OneBodyG is required to output green function.") - + if "initial" in info_outputfile.keys(): output_str = "===============================\n" - green_nonzero = self.Green[np.where(abs(self.Green)>0)] + green_nonzero = self.Green[np.where(abs(self.Green) > 0)] ncisajs = len(green_nonzero) output_str += "NCisAjs {}\n".format(ncisajs) output_str += "===============================\n" output_str += "===============================\n" output_str += "===============================\n" for nsite1, nsite2 in itertools.product(range(2 * self.Nsize), range(2 * self.Nsize)): - if abs(self.Green[nsite1][nsite2])>0: - output_str += "{} {} {} {} {} {}\n".format(nsite1%self.Nsize, nsite1//self.Nsize, nsite2%self.Nsize, nsite2//self.Nsize, + if abs(self.Green[nsite1][nsite2]) > 0: + output_str += "{} {} {} {} {} {}\n".format(nsite1 % self.Nsize, nsite1 // self.Nsize, nsite2 % self.Nsize, nsite2 // self.Nsize, self.Green[nsite1][nsite2].real, self.Green[nsite1][nsite2].imag) with open(os.path.join(path_to_output, info_outputfile["initial"]), "w") as fw: fw.write(output_str) @@ -630,12 +1000,12 @@ def save_results(self, info_outputfile, green_info): output_str = "" for i, j in itertools.product(range(Ns), range(Ns)): fij = 0.0 - for n in range(Ne//2): - fij += ev[i][n*2] * ev[j][n*2+1] - ev[i][n*2+1] * ev[j][n*2] + for n in range(Ne // 2): + fij += ev[i][n * 2] * ev[j][n * 2 + 1] - ev[i][n * 2 + 1] * ev[j][n * 2] output_str += " {:3} {:3} {:.12f} {:.12f}\n".format( i, j, np.real(fij), np.imag(fij)) - with open(os.path.join(path_to_output, key+"_"+info_outputfile["fij"]), "w") as fw: + with open(os.path.join(path_to_output, key + "_" + info_outputfile["fij"]), "w") as fw: fw.write(output_str) else: TwoSz = self.param_mod["2Sz"] @@ -650,16 +1020,23 @@ def save_results(self, info_outputfile, green_info): output_str = "" for i, j in itertools.product(range(Ns), range(Ns)): fij = 0.0 - for n in range(Ne//2): + for n in range(Ne // 2): fij += ev_up[i][n] * ev_dn[j][n] output_str += " {:3} {:3} {:.12f} {:.12f}\n".format( i, j, np.real(fij), np.imag(fij)) - with open(os.path.join(path_to_output, "Sz0_"+info_outputfile["fij"]), "w") as fw: + with open(os.path.join(path_to_output, "Sz0_" + info_outputfile["fij"]), "w") as fw: fw.write(output_str) else: logger.warning("NOT IMPLEMENTED: Sz even and Sz != 0: this case will be implemented in near future") + @do_profile def get_Ham(self): - return self.Ham + """ + Get the Hamiltonian matrix. + + Returns: + np.ndarray: The Hamiltonian matrix. + """ + return self.Ham \ No newline at end of file