From aa3701c8b1a2f5e523ec92a3bcaf72d25d774654 Mon Sep 17 00:00:00 2001 From: k-yoshimi Date: Mon, 11 Nov 2024 10:49:06 +0900 Subject: [PATCH 1/8] Update docstring format to numpydoc format in rpa.py --- src/hwave/solver/rpa.py | 180 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 172 insertions(+), 8 deletions(-) diff --git a/src/hwave/solver/rpa.py b/src/hwave/solver/rpa.py index 5657cec..e7028be 100644 --- a/src/hwave/solver/rpa.py +++ b/src/hwave/solver/rpa.py @@ -678,6 +678,28 @@ def _show_params(self): @do_profile def solve(self, green_info, path_to_output): + """Solve the RPA equation to calculate susceptibility. + + This is the main method that performs RPA calculations. It either calculates + or loads chi0q, transforms interaction Hamiltonians based on spin state, + and solves the RPA equation. + + Parameters + ---------- + green_info : dict + Dictionary containing Green's function information and calculation parameters. + May include pre-calculated chi0q. + path_to_output : str + Path to directory where output files will be saved. + + Notes + ----- + The calculation flow is: + 1. Calculate/load chi0q + 2. Transform interactions based on spin state (spin-free/spinful/spin-diag) + 3. Transform tensors based on calculation scheme (reduced/squashed/general) + 4. Solve RPA equation to get chiq + """ logger.info("Start RPA calculations") beta = 1.0/self.T @@ -861,6 +883,23 @@ def solve(self, green_info, path_to_output): @do_profile def save_results(self, info_outputfile, green_info): + """Save calculation results to files. + + Parameters + ---------- + info_outputfile : dict + Dictionary containing output file configuration. + green_info : dict + Dictionary containing calculation information and results. + + Notes + ----- + Saves the following information: + - Calculated correlation functions (chi0q/chiq) + - Matsubara frequency indices + - Wave vector unit vectors + - Wave vector indices + """ logger.info("Save RPA results") path_to_output = info_outputfile["path_to_output"] @@ -980,6 +1019,31 @@ def read_init(self, info_inputfile): return info def _read_chi0q(self, file_name): + """Read chi0q from file and validate its shape. + + Parameters + ---------- + file_name : str + Path to the file containing chi0q data. + + Returns + ------- + ndarray + The loaded chi0q array. + + Raises + ------ + AssertionError + If the loaded data doesn't match expected dimensions or format. + + Notes + ----- + Performs checks for: + - Lattice volume consistency + - Number of orbitals consistency + - Spin freedom consistency + - Calculation scheme compatibility + """ logger.debug(">>> RPA._read_chi0q") try: @@ -1060,6 +1124,28 @@ def _read_chi0q(self, file_name): return chi0q def _read_trans_mod(self, file_name): + """Read transfer integral modifications from file. + + Parameters + ---------- + file_name : str + Path to the file containing transfer integral modifications. + + Returns + ------- + ndarray + The transfer integrals in k-space. + + Raises + ------ + RuntimeError + If file reading fails. + + Notes + ----- + - Converts layout if sublattice exists + - Performs Fourier transform to k-space + """ logger.debug(">>> RPA._read_trans_mod") try: @@ -1350,12 +1436,35 @@ def _calc_green(self, beta, mu): @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 bare susceptibility chi0q. + + Parameters + ---------- + green_kw : ndarray + Green's function in k-space and Matsubara frequency. + Shape: [g,l,k,a,b] where: + - g: block index + - l: Matsubara frequency index + - k: wave number index + - a,b: orbital and spin indices + green0_tail : ndarray + High-frequency tail correction for Green's function. + beta : float + Inverse temperature. + + Returns + ------- + ndarray + The calculated chi0q. + + Notes + ----- + Calculation steps: + 1. Fourier transform from Matsubara freq to imaginary time + 2. Transform from k-space to real space + 3. Calculate chi0 in real space and imaginary time + 4. Transform back to k-space and Matsubara frequency + """ logger.debug(">>> RPA._calc_chi0q") nx,ny,nz = self.lattice.shape @@ -1365,7 +1474,7 @@ def _calc_chi0q(self, green_kw, green0_tail, beta): assert nvol == self.lattice.nvol assert nmat == self.nmat - # Fourier transform from matsubara freq to imaginary time + # Fourier transform from Matsubara freq to imaginary time omg = np.exp(-1j * np.pi * (1.0/nmat - 1.0) * np.arange(nmat)) green_kt = np.einsum('gtv,t->gtv', @@ -1376,7 +1485,7 @@ def _calc_chi0q(self, green_kw, green0_tail, beta): # Fourier transform from wave number space to coordinate space green_rt = FFT.ifftn(green_kt.reshape(nblock,nmat,nx,ny,nz,nd*nd), axes=(2,3,4)) - # calculate chi0(r,t)[a,a',b,b'] = G(r,t)[a,b] * G(-r,-t)[b',a'] + # calculate chi0 in real space and imaginary time green_rev = np.flip(np.roll(green_rt, -1, axis=(1,2,3,4)), axis=(1,2,3,4)).reshape(nblock,nmat,nvol,nd,nd) sgn = np.full(nmat, -1) @@ -1412,6 +1521,26 @@ def _calc_chi0q(self, green_kw, green0_tail, beta): @do_profile def _solve_rpa(self, chi0q, ham): + """Solve the RPA equation. + + Parameters + ---------- + chi0q : ndarray + Bare susceptibility. + ham : ndarray + Interaction Hamiltonian. + + Returns + ------- + ndarray + The RPA susceptibility chiq. + + Notes + ----- + Solves the equation: + chiq = [1 + chi0q * W]^(-1) * chi0q + where W is the interaction vertex. + """ logger.debug(">>> RPA._solve_rpa") nvol = self.lattice.nvol @@ -1435,6 +1564,41 @@ def _solve_rpa(self, chi0q, ham): def run(*, input_dict: Optional[dict] = None, input_file: Optional[str] = None): + """Main entry point for running RPA calculations. + + Parameters + ---------- + input_dict : dict, optional + Dictionary containing input parameters. Must be provided if input_file is None. + input_file : str, optional + Path to input file. Not used if input_dict is provided. + + Raises + ------ + RuntimeError + If input_dict is None. + + Notes + ----- + The input dictionary should contain the following sections: + - log: Logging configuration + - print_level: Verbosity level (default: 1) + - print_step: Print frequency (default: 1) + - mode: Calculation mode configuration + - mode: Must be "RPA" for RPA calculations + - file: File paths configuration + - input: Input file paths + - path_to_input: Input directory (default: "") + - output: Output file paths + - path_to_output: Output directory (default: "output") + + The function performs the following steps: + 1. Initializes logging and file paths + 2. Reads interaction definitions + 3. Creates RPA solver instance + 4. Performs RPA calculations + 5. Saves results + """ logger = logging.getLogger("run") if input_dict is None: From 85ae6f1da2ecb964f0c3ee5a6fd8701201e2f94b Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Mon, 11 Nov 2024 12:57:45 +0900 Subject: [PATCH 2/8] Update docstrings to numpydoc format in dos.py, qlms.py, solver/base.py, and solver/perf.py --- src/hwave/dos.py | 90 ++++++++++++++++++++++++++-- src/hwave/qlms.py | 65 ++++++++++++++------ src/hwave/solver/base.py | 124 +++++++++++++++++++++++++++++++++++---- src/hwave/solver/perf.py | 47 +++++++++++++++ 4 files changed, 294 insertions(+), 32 deletions(-) diff --git a/src/hwave/dos.py b/src/hwave/dos.py index db7bd61..0c2eefe 100644 --- a/src/hwave/dos.py +++ b/src/hwave/dos.py @@ -1,3 +1,10 @@ +"""Density of states (DoS) calculation and plotting utilities. + +This module provides functionality for calculating, plotting and analyzing the +density of states from Hartree-Fock calculations. + +""" + from __future__ import annotations import itertools @@ -9,10 +16,26 @@ class DoS: - dos: np.ndarray - ene: np.ndarray - ene_num: int - norb: int + """Class for storing and manipulating density of states data. + + Parameters + ---------- + ene : np.ndarray + Energy grid points + dos : np.ndarray + Density of states values for each orbital at each energy point + + Attributes + ---------- + dos : np.ndarray + Density of states array with shape (norb, nene) + ene : np.ndarray + Energy grid points array with shape (nene,) + ene_num : int + Number of energy points + norb : int + Number of orbitals + """ def __init__(self, ene: np.ndarray, dos: np.ndarray): assert ene.shape[0] == dos.shape[1] @@ -22,6 +45,17 @@ 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. + + Creates a plot showing the total DoS and orbital-resolved DoS. + + Parameters + ---------- + filename : str, optional + If provided, save plot to this file + verbose : bool, optional + If True, print additional output + """ try: import matplotlib.pyplot as plt except ImportError: @@ -45,6 +79,15 @@ def plot(self, filename: str = "", verbose: bool = False): plt.close() def write_dos(self, output: str, verbose: bool = False): + """Write density of states data to file. + + Parameters + ---------- + output : str + Output filename + verbose : bool, optional + If True, print additional output + """ if verbose: print("Writing DOS to file: ", output) total_dos = np.sum(self.dos, axis=0) @@ -62,6 +105,18 @@ def write_dos(self, output: str, verbose: bool = False): def __read_geom(file_name="./dir-model/zvo_geom.dat"): + """Read geometry data from file. + + Parameters + ---------- + file_name : str, optional + Path to geometry file + + Returns + ------- + np.ndarray + Unit cell vectors array with shape (3,3) + """ 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 @@ -75,6 +130,29 @@ def calc_dos( ene_num: int = 101, verbose: bool = False, ) -> DoS: + """Calculate density of states. + + Parameters + ---------- + input_dict : dict + Input parameters dictionary + ene_window : list, optional + Energy window [min, max] for DoS calculation + ene_num : int, optional + Number of energy points + verbose : bool, optional + If True, print additional output + + Returns + ------- + DoS + Calculated density of states object + + Raises + ------ + ImportError + If required libtetrabz package is not installed + """ try: import libtetrabz except ImportError: @@ -130,6 +208,10 @@ def calc_dos( def main(): + """Command-line interface for DoS calculation. + + Parses command line arguments and runs DoS calculation. + """ import tomli import argparse diff --git a/src/hwave/qlms.py b/src/hwave/qlms.py index c7b5068..a463183 100644 --- a/src/hwave/qlms.py +++ b/src/hwave/qlms.py @@ -16,53 +16,77 @@ from requests.structures import CaseInsensitiveDict def run(*, input_dict: Optional[dict] = None, input_file: Optional[str] = None): + """Run Hartree-Fock calculation with given input parameters. + + Parameters + ---------- + input_dict : dict, optional + Dictionary containing input parameters. Cannot be used with input_file. + input_file : str, optional + Path to TOML input file. Cannot be used with input_dict. + + Raises + ------ + RuntimeError + If neither or both input_dict and input_file are provided + """ + # Validate input arguments if input_dict is None: if input_file is None: raise RuntimeError("Neither input_dict nor input_file are passed") + # Load parameters from TOML file with open(input_file, "rb") as f: input_dict = tomli.load(f) else: if input_file is not None: raise RuntimeError("Both input_dict and input_file are passed") - # Initialize information about log + # Initialize logging configuration info_log = input_dict.get("log", {}) - info_log["print_level"] = info_log.get("print_level", 1) - info_log["print_step"] = info_log.get("print_step", 1) + info_log["print_level"] = info_log.get("print_level", 1) # Default print level + info_log["print_step"] = info_log.get("print_step", 1) # Default print step - # Initialize information about mode + # Get calculation mode and file paths info_mode = input_dict.get("mode", {}) info_file = input_dict.get("file", {"input": {}, "output": {}}) - # Initialize information about input files + + # Setup input file paths info_inputfile = info_file.get("input", {}) info_inputfile["path_to_input"] = info_inputfile.get("path_to_input", "") - # Initialize information about output files + # Setup output directory info_outputfile = info_file.get("output", {}) info_outputfile["path_to_output"] = info_outputfile.get("path_to_output", "output") path_to_output = info_outputfile["path_to_output"] os.makedirs(path_to_output, exist_ok=True) + # Configure logging logger = logging.getLogger("qlms") fmt = "%(asctime)s %(levelname)s %(name)s: %(message)s" - # logging.basicConfig(level=logging.DEBUG, format=fmt) logging.basicConfig(level=logging.INFO, format=fmt) + # Validate calculation mode if "mode" not in info_mode: logger.error("mode is not defined in [mode].") exit(1) mode = info_mode["mode"] + + # Initialize solver based on calculation mode if mode == "UHFr": + # Real-space unrestricted Hartree-Fock logger.info("Read def files") file_list = CaseInsensitiveDict() - #interaction files + + # Process interaction file paths for key, file_name in info_inputfile["interaction"].items(): - if key.lower() in ["trans", "coulombinter", "coulombintra", "pairhop", "hund", "exchange", "ising", "pairlift", "interall"]: + if key.lower() in ["trans", "coulombinter", "coulombintra", "pairhop", + "hund", "exchange", "ising", "pairlift", "interall"]: file_list[key] = os.path.join(info_inputfile["path_to_input"], file_name) else: logging.error("Keyword {} is incorrect.".format(key)) exit(1) - #initial and green + + # Process initial state and Green's function files for key, file_name in info_inputfile.items(): if key.lower() == "initial": file_list[key] = os.path.join(info_inputfile["path_to_input"], file_name) @@ -70,34 +94,30 @@ def run(*, input_dict: Optional[dict] = None, input_file: Optional[str] = None): file_list[key] = os.path.join(info_inputfile["path_to_input"], file_name) read_io = qlmsio.read_input.QLMSInput(file_list) + # Read Hamiltonian and Green's function parameters logger.info("Get Hamiltonian information") ham_info = read_io.get_param("ham") logger.info("Get Green function information") green_info = read_io.get_param("green") - # solver = sol_uhf.UHF(ham_info, info_log, info_mode, mod_param_info) solver = sol_uhfr.UHFr(ham_info, info_log, info_mode) elif mode == "UHFk": + # k-space unrestricted Hartree-Fock logger.info("Read definitions from files") read_io = qlmsio.read_input_k.QLMSkInput(info_inputfile) - # logger.info("Get parameter information") - # mod_param_info = info_mode #read_io.get_param("mod") - # pprint.pprint(mod_param_info, width = 1) - logger.info("Get Hamiltonian information") ham_info = read_io.get_param("ham") logger.info("Get Green function information") green_info = read_io.get_param("green") - # pprint.pprint(info_mode, width=1) - solver = sol_uhfk.UHFk(ham_info, info_log, info_mode) elif mode == "RPA": + # Random Phase Approximation logger.info("RPA mode") logger.info("Read interaction definitions from files") @@ -107,12 +127,13 @@ def run(*, input_dict: Optional[dict] = None, input_file: Optional[str] = None): solver = sol_rpa.RPA(ham_info, info_log, info_mode) green_info = read_io.get_param("green") - green_info.update( solver.read_init(info_inputfile) ) + green_info.update(solver.read_init(info_inputfile)) else: logger.warning("mode is incorrect: mode={}.".format(mode)) exit(0) + # Execute calculation logger.info("Start UHF calculation") solver.solve(green_info, path_to_output) logger.info("Save calculation results.") @@ -121,20 +142,28 @@ def run(*, input_dict: Optional[dict] = None, input_file: Optional[str] = None): def main(): + """Command-line interface entry point. + + Parses command line arguments and runs the calculation. + """ import argparse + # Setup argument parser parser = argparse.ArgumentParser(prog='hwave') parser.add_argument('input_toml', nargs='?', default=None, help='input parameter file') parser.add_argument('--version', action='store_true', help='show version') args = parser.parse_args() + # Handle version request if args.version: print('hwave', hwave.__version__) sys.exit(0) + # Validate input file if args.input_toml is None: parser.print_help() sys.exit(1) + # Run calculation run(input_file = args.input_toml) diff --git a/src/hwave/solver/base.py b/src/hwave/solver/base.py index 71780eb..1568e1d 100644 --- a/src/hwave/solver/base.py +++ b/src/hwave/solver/base.py @@ -1,11 +1,43 @@ +"""Base solver class for Hartree-Fock calculations. + +This module provides the base solver class that implements common functionality +for Hartree-Fock calculations, including parameter handling and validation. + +""" import sys from requests.structures import CaseInsensitiveDict -# from pprint import pprint - import logging logger = logging.getLogger("qlms").getChild("solver") + class solver_base(): + """Base solver class for Hartree-Fock calculations. + + Parameters + ---------- + param_ham : dict + Hamiltonian parameters + info_log : dict + Logging configuration + info_mode : dict + Calculation mode parameters + param_mod : dict, optional + Model parameters to override defaults + + Attributes + ---------- + param_mod : CaseInsensitiveDict + Model parameters + param_ham : dict + Hamiltonian parameters + info_log : dict + Logging configuration + threshold : float + Cutoff threshold for Green's function elements + relax_checks : bool + Whether to relax parameter validation checks + """ + def __init__(self, param_ham, info_log, info_mode, param_mod=None): logger = logging.getLogger("qlms").getChild(self.name) @@ -32,12 +64,8 @@ def __init__(self, param_ham, info_log, info_mode, param_mod=None): range_list = { "T": [ 0, None ], - # "2Sz": [ -param_mod["Nsite"], param_mod["Nsite"] ], - # "Nsite": [ 1, None ], - # "Ncond": [ 1, None ], "IterationMax": [ 0, None ], "Mix": [ 0.0, 1.0 ], - # "print_step": [ 1, None ], "EPS": [ 0, None ], } @@ -114,15 +142,23 @@ def __init__(self, param_ham, info_log, info_mode, param_mod=None): # canonicalize self.param_mod["EPS"] = pow(10, -self.param_mod["EPS"]) - # debug - # pprint(self.param_mod) - def _check_info_mode(self, info_mode): + """Check validity of info_mode parameters. + + Parameters + ---------- + info_mode : dict + Mode parameters to validate + + Returns + ------- + int + Number of validation errors found + """ logger = logging.getLogger("qlms").getChild(self.name) fix_list = { "mode": ["UHFr", "UHFk"], - # "flag_fock": [True, False] } exit_code = 0 @@ -136,6 +172,18 @@ def _check_info_mode(self, info_mode): return exit_code def _check_param_mod(self, param_mod): + """Check validity of model parameters. + + Parameters + ---------- + param_mod : dict + Model parameters to validate + + Returns + ------- + int + Number of validation errors found + """ logger = logging.getLogger("qlms").getChild(self.name) error_code = 0 @@ -149,6 +197,20 @@ def _check_param_mod(self, param_mod): return error_code def _check_param_range(self, param_mod, range_list): + """Check if parameters are within valid ranges. + + Parameters + ---------- + param_mod : dict + Model parameters to validate + range_list : dict + Valid ranges for parameters + + Returns + ------- + int + Number of validation errors found + """ logger = logging.getLogger("qlms").getChild(self.name) error_code = 0 @@ -166,6 +228,25 @@ def _check_param_range(self, param_mod, range_list): return error_code def _round_to_int(self, val, mode): + """Round a value to integer according to specified mode. + + Parameters + ---------- + val : float + Value to round + mode : str + Rounding mode to use + + Returns + ------- + int + Rounded integer value + + Raises + ------ + SystemExit + If rounding fails or mode is invalid + """ import math mode = mode.lower() # case-insensitive if mode == "as-is": @@ -197,10 +278,33 @@ def _round_to_int(self, val, mode): return ret def solve(self, path_to_output): + """Solve the Hartree-Fock equations. + + Parameters + ---------- + path_to_output : str + Path to output file + """ pass def get_results(self): + """Get calculation results. + + Returns + ------- + tuple + (physics, Green's function) results + """ return (self.physics, self.Green) def save_results(self, info_outputfile, green_info): + """Save calculation results. + + Parameters + ---------- + info_outputfile : dict + Output file configuration + green_info : dict + Green's function information + """ pass diff --git a/src/hwave/solver/perf.py b/src/hwave/solver/perf.py index b7c9677..5ee0bbc 100644 --- a/src/hwave/solver/perf.py +++ b/src/hwave/solver/perf.py @@ -1,11 +1,34 @@ +"""Performance profiling utilities. + +This module provides utilities for profiling function execution times and +collecting performance statistics. +""" + from functools import wraps import time + class PerfDB: + """Database for storing performance profiling data. + + Stores execution counts and total elapsed time for profiled functions. + Prints summary statistics on deletion. + + Attributes + ---------- + _db_count : dict + Number of calls per function + _db_value : dict + Total elapsed time per function + """ + def __init__(self): + """Initialize empty performance database.""" self._db_count = {} self._db_value = {} + def __del__(self): + """Print summary statistics when object is deleted.""" if len(self._db_count) == 0: return print("--------------------------------------------------------------------------------") @@ -20,13 +43,37 @@ def __del__(self): self._db_count[item] )) print("--------------------------------------------------------------------------------") + def put(self, name, value): + """Add a timing measurement. + + Parameters + ---------- + name : str + Function name + value : float + Elapsed time in seconds + """ 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): + """Decorator for profiling function execution time. + + Parameters + ---------- + func : callable + Function to profile + + Returns + ------- + callable + Wrapped function that measures and records execution time + """ @wraps(func) def wrapper(*args, **kwargs): # start time From 8fb65a2c5038173ff58c0c5cd9f6def9ba380bd8 Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Mon, 11 Nov 2024 13:35:50 +0900 Subject: [PATCH 3/8] Update docstrings to numpydoc format in uhfk.py --- src/hwave/solver/uhfk.py | 193 +++++++++++++++++++++++++++++++++++---- 1 file changed, 177 insertions(+), 16 deletions(-) diff --git a/src/hwave/solver/uhfk.py b/src/hwave/solver/uhfk.py index ffbbdc7..df7f3e8 100644 --- a/src/hwave/solver/uhfk.py +++ b/src/hwave/solver/uhfk.py @@ -10,32 +10,48 @@ class UHFk(solver_base): @do_profile def __init__(self, param_ham, info_log, info_mode, param_mod=None): + """Initialize the UHF solver. + + Parameters + ---------- + param_ham : dict + Hamiltonian parameters + info_log : dict + Logging configuration + info_mode : dict + Solver mode parameters + param_mod : dict, optional + Additional model parameters + """ self.name = "uhfk" super().__init__(param_ham, info_log, info_mode, param_mod) - self._init_mode(info_mode) - self._init_param() - self._init_lattice() - self._init_orbit() + # Initialize solver components + self._init_mode(info_mode) # Set solver modes like Fock term + self._init_param() # Initialize parameters + self._init_lattice() # Setup lattice geometry + self._init_orbit() # Setup orbital structure # self._dump_param_ham() - self._check_interaction() - self._init_interaction() + self._check_interaction() # Validate interactions + self._init_interaction() # Setup interaction terms # self._dump_param_ham() - self._init_wavevec() + self._init_wavevec() # Setup k-vectors - self._show_param() + self._show_param() # Display parameters - # local data + # Initialize physical quantities self.physics = { "Ene": { "Total": 0.0, "Band": 0.0 }, - "NCond": 0.0, - "Sz": 0.0, - "Rest": 1.0 + "NCond": 0.0, # Total particle number + "Sz": 0.0, # Total spin + "Rest": 1.0 # Convergence measure } - nvol = self.nvol - nd = self.nd + # Setup dimensions + nvol = self.nvol # Number of k-points + nd = self.nd # Total number of basis states + # Storage for eigenvalues/vectors self._green_list = { # "eigenvalue": np.zeros((nvol,nd), dtype=np.complex128), # "eigenvector": np.zeros((nvol,nd,nd), dtype=np.complex128) @@ -43,15 +59,34 @@ def __init__(self, param_ham, info_log, info_mode, param_mod=None): #self.Green = np.zeros((nvol,ns,norb,ns,norb), dtype=np.complex128), - # work area + # Initialize Hamiltonian matrix self.ham = np.zeros((nvol,nd,nd), dtype=np.complex128) def _init_mode(self, param): + """Initialize solver modes. + + Parameters + ---------- + param : dict + Mode parameters + """ + # Enable/disable Fock term self.iflag_fock = param.get('flag_fock', True) + + # Enable/disable spin-orbital coupling self.enable_spin_orbital = param.get('enable_spin_orbital', False) @do_profile def _init_param(self): + """Initialize solver parameters. + + Sets up key parameters including: + - Cell shape + - Temperature + - Number of electrons/filling + - Sz constraints + - Numerical precision parameters + """ # check and store parameters # - cell shape @@ -99,6 +134,15 @@ def _init_param(self): @do_profile def _init_lattice(self): + """Initialize lattice geometry. + + Sets up: + - Cell shape and volume + - Sublattice structure + - Consistency checks for sublattice + + Parameters are read from self.param_mod. + """ Lx,Ly,Lz = self.param_mod.get("CellShape") self.cellshape = (Lx,Ly,Lz) self.cellvol = Lx * Ly * Lz @@ -128,6 +172,14 @@ def _init_lattice(self): @do_profile def _init_orbit(self): + """Initialize orbital structure. + + Sets up: + - Number of orbitals + - Spin degrees of freedom + - Total basis dimension + Takes into account supercell structure if present. + """ norb = self.param_ham["Geometry"]["norb"] ns = 2 # spin dof @@ -140,6 +192,15 @@ def _init_orbit(self): @do_profile def _show_param(self): + """Display solver parameters. + + Logs key parameters including: + - Enabled features (Fock term, spin-orbital coupling) + - Lattice dimensions + - Orbital/spin structure + - Physical parameters (filling, temperature) + - Numerical parameters + """ logger.info("Show parameters") logger.info(" Enable Fock = {}".format(self.iflag_fock)) logger.info(" Cell Shape = {}".format(self.cellshape)) @@ -169,6 +230,18 @@ def _show_param(self): @do_profile def _reshape_geometry(self, geom): + """Reshape geometry for sublattice structure. + + Parameters + ---------- + geom : dict + Original geometry definition + + Returns + ------- + dict + Reshaped geometry for sublattice + """ Bx,By,Bz = self.subshape bvol = Bx * By * Bz @@ -193,6 +266,20 @@ def _reshape_geometry(self, geom): @do_profile def _reshape_interaction(self, ham, enable_spin_orbital): + """Reshape interaction terms for sublattice. + + Parameters + ---------- + ham : dict + Original interaction terms + enable_spin_orbital : bool + Whether spin-orbital coupling is enabled + + Returns + ------- + dict + Reshaped interaction terms + """ Bx,By,Bz = self.subshape nx,ny,nz = self.shape @@ -246,6 +333,18 @@ def _round(x, n): @do_profile def _reshape_green(self, green): + """Convert Green function to sublattice basis. + + Parameters + ---------- + green : ndarray + Green function in original basis + + Returns + ------- + ndarray + Green function in sublattice basis + """ # convert green function into sublattice Lx,Ly,Lz = self.cellshape @@ -312,6 +411,18 @@ def _unpack_site(idx, n): @do_profile def _deflate_green(self, green_sub): + """Convert Green function back to original basis. + + Parameters + ---------- + green_sub : ndarray + Green function in sublattice basis + + Returns + ------- + ndarray + Green function in original basis + """ # convert green function back to original lattice Lx,Ly,Lz = self.cellshape Lvol = Lx * Ly * Lz @@ -370,6 +481,11 @@ def _unpack_site(idx, n): @do_profile def _init_interaction(self): + """Initialize interaction terms. + + Processes interaction terms for sublattice structure if needed. + Handles different types of interactions (Coulomb, Hund, etc). + """ # reinterpret interaction coefficient on sublattice if self.has_sublattice: # backup @@ -391,6 +507,13 @@ def _init_interaction(self): self.param_ham[type] = tbl def _init_wavevec(self): + """Initialize wave vectors. + + Sets up: + - k-vectors for sublatticed geometry + - Reciprocal lattice vectors + - Wave number tables + """ # wave vectors on sublatticed geometry def _klist(n): return np.roll( (np.arange(n)-(n//2)), -(n//2) ) @@ -436,6 +559,13 @@ def _dump_param_ham(self): @do_profile def _check_interaction(self): + """Validate interaction terms. + + Performs checks on: + - Cell size compatibility + - Orbital index validity + - Hermiticity of terms + """ self._check_cellsize() self._check_orbital_index() self._check_hermite() @@ -557,6 +687,21 @@ def _check_hermite(self): @do_profile def solve(self, green_info, path_to_output): + """Main solver routine. + + Performs UHF iteration until convergence: + 1. Constructs Hamiltonian + 2. Diagonalizes + 3. Updates Green function + 4. Calculates physical quantities + + Parameters + ---------- + green_info : dict + Initial Green function information + path_to_output : str + Path for output files + """ print_level = self.info_log["print_level"] print_check = self.info_log.get("print_check", None) @@ -1283,7 +1428,7 @@ def _fermi(t, mu, ev): 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)) ) + #v = np.where( w > self.ene_cutoff, 0.0, 1.0 / (1.0 + np.exp(w)) ) ) return v wt = -(w - mu) / T @@ -1341,6 +1486,22 @@ def get_results(self): @do_profile def save_results(self, info_outputfile, green_info): + """Save results to files. + + Parameters + ---------- + info_outputfile : dict + Output file configuration + green_info : dict + Green function information + + Saves: + - Energy and physical quantities + - Eigenvalues/vectors + - Green functions + - One-body quantities + - RPA data if requested + """ path_to_output = info_outputfile["path_to_output"] if "energy" in info_outputfile.keys(): From 1c3512293bab221835c196efa2f8e69e644b6f41 Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Mon, 11 Nov 2024 14:23:21 +0900 Subject: [PATCH 4/8] Update docstrings to numpydoc format in uhfr.py --- src/hwave/solver/uhfr.py | 559 ++++++++++++++++++++++++++++++++++----- 1 file changed, 496 insertions(+), 63 deletions(-) diff --git a/src/hwave/solver/uhfr.py b/src/hwave/solver/uhfr.py index a2db826..03cf0f3 100644 --- a/src/hwave/solver/uhfr.py +++ b/src/hwave/solver/uhfr.py @@ -1,3 +1,10 @@ +"""Unrestricted Hartree-Fock solver with real-space representation. + +This module implements the unrestricted Hartree-Fock (UHF) method for solving +many-body quantum systems in real space. It supports both zero and finite +temperature calculations. +""" + import itertools import logging import numpy as np @@ -8,7 +15,28 @@ logger = logging.getLogger("qlms").getChild("uhfr") -class Interact_UHFr_base(): +class Interact_UHFr_base: + """Base class for interaction terms in UHF calculations. + + Parameters + ---------- + ham_info : dict + Dictionary containing interaction parameters + Nsize : int + System size (number of sites) + + Attributes + ---------- + Nsize : int + System size + Ham_tmp : ndarray + Temporary array for Hamiltonian construction + Ham_trans_tmp : ndarray + Temporary array for transformed Hamiltonian + param_ham : dict + Transformed interaction parameters + """ + def __init__(self, ham_info, Nsize): self.Nsize = Nsize self.Ham_tmp = np.zeros(tuple([(2 * self.Nsize) for i in range(4)]), dtype=complex) @@ -16,40 +44,23 @@ def __init__(self, ham_info, Nsize): self.param_ham = self._transform_interall(ham_info) self._check_range() - #Change interaction to interall type def _transform_interall(self, ham_info): + """Transform interaction parameters to internal format. + + Parameters + ---------- + ham_info : dict + Input interaction parameters + + Returns + ------- + dict + Transformed parameters + """ return ham_info - def _calc_hartree(self): - 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 - # Diagonal Fock term - self.Ham_tmp[site[0]][site[1]][site[2]][site[3]] += value - self.Ham_tmp[site[2]][site[3]][site[0]][site[1]] += value - if site[1] == site[2]: - self.Ham_trans_tmp[site[1]][site[2]] += value - pass - - def _calc_fock(self): - 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 - # OffDiagonal Fock term - self.Ham_tmp[site[0]][site[3]][site[2]][site[1]] -= value - self.Ham_tmp[site[2]][site[1]][site[0]][site[3]] -= value - pass - - def get_ham(self, type): - self._calc_hartree() - if type == "hartreefock": - self._calc_fock() - return self.Ham_tmp, self.Ham_trans_tmp - - # check input def _check_range(self): + """Check that site indices are within valid range.""" err = 0 for site_info, value in self.param_ham.items(): for i in range(4): @@ -62,6 +73,15 @@ def _check_range(self): exit(1) def _check_hermite(self, strict_hermite, tolerance): + """Check Hermiticity of interaction parameters. + + Parameters + ---------- + strict_hermite : bool + If True, exit on Hermiticity violation + 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,9 +98,27 @@ def _check_hermite(self, strict_hermite, tolerance): logger.warn(msg) def check_hermite(self, strict_hermite=False, tolerance=1.0e-8): + """Public interface for Hermiticity check. + + Parameters + ---------- + strict_hermite : bool, optional + If True, exit on Hermiticity violation + tolerance : float, optional + Tolerance for Hermiticity check + """ self._check_hermite(strict_hermite, tolerance) class CoulombIntra_UHFr(Interact_UHFr_base): + """On-site Coulomb interaction term. + + Parameters + ---------- + ham_info : dict + Interaction parameters + Nsize : int + System size + """ def __init__(self, ham_info, Nsize): self.__name__ = "CoulombIntra" super().__init__(ham_info, Nsize) @@ -92,6 +130,15 @@ def _transform_interall(self, ham_info): return param_tmp class CoulombInter_UHFr(Interact_UHFr_base): + """Inter-site Coulomb interaction term. + + Parameters + ---------- + ham_info : dict + Interaction parameters + Nsize : int + System size + """ def __init__(self, ham_info, Nsize): self.__name__ = "CoulombInter" super().__init__(ham_info, Nsize) @@ -104,6 +151,15 @@ def _transform_interall(self, ham_info): return param_tmp class Hund_UHFr(Interact_UHFr_base): + """Hund's coupling term. + + Parameters + ---------- + ham_info : dict + Interaction parameters + Nsize : int + System size + """ def __init__(self, ham_info, Nsize): self.__name__ = "Hund" super().__init__(ham_info, Nsize) @@ -116,6 +172,15 @@ def _transform_interall(self, ham_info): return param_tmp class PairHop_UHFr(Interact_UHFr_base): + """Pair hopping term. + + Parameters + ---------- + ham_info : dict + Interaction parameters + Nsize : int + System size + """ def __init__(self, ham_info, Nsize): self.__name__ = "PairHop" super().__init__(ham_info, Nsize) @@ -129,6 +194,15 @@ def _transform_interall(self, ham_info): return param_tmp class Exchange_UHFr(Interact_UHFr_base): + """Exchange interaction term. + + Parameters + ---------- + ham_info : dict + Interaction parameters + Nsize : int + System size + """ def __init__(self, ham_info, Nsize): self.__name__ = "Exchange" super().__init__(ham_info, Nsize) @@ -142,6 +216,15 @@ def _transform_interall(self, ham_info): return param_tmp class Ising_UHFr(Interact_UHFr_base): + """Ising interaction term. + + Parameters + ---------- + ham_info : dict + Interaction parameters + Nsize : int + System size + """ def __init__(self, ham_info, Nsize): self.__name__ = "Ising" super().__init__(ham_info, Nsize) @@ -154,6 +237,15 @@ def _transform_interall(self, ham_info): return param_tmp class PairLift_UHFr(Interact_UHFr_base): + """Pair lifting term. + + Parameters + ---------- + ham_info : dict + Interaction parameters + Nsize : int + System size + """ def __init__(self, ham_info, Nsize): self.__name__ = "PairLift" super().__init__(ham_info, Nsize) @@ -167,6 +259,19 @@ def _transform_interall(self, ham_info): return param_tmp class InterAll_UHFr(Interact_UHFr_base): + """General interaction term. + + Parameters + ---------- + ham_info : dict + Interaction parameters + Nsize : int + System size + strict_hermite : bool + If True, exit on Hermiticity violation + tolerance : float + Tolerance for Hermiticity check + """ def __init__(self, ham_info, Nsize, strict_hermite, tolerance): self.__name__ = "InterAll" super().__init__(ham_info, Nsize) @@ -175,6 +280,26 @@ def _transform_interall(self, ham_info): return ham_info class Term_base: + """Base class for single-particle terms. + + Parameters + ---------- + term_info : dict + Term parameters + Nsize : int + System size + coeff : float, optional + Overall coefficient + + Attributes + ---------- + Nsize : int + System size + term_info : dict + Term parameters + coeff : float + Overall coefficient + """ def __init__(self, term_info, Nsize, coeff=1.0): self.Nsize = Nsize self.term_info = term_info @@ -182,6 +307,13 @@ def __init__(self, term_info, Nsize, coeff=1.0): self._check_range() def get_data(self): + """Get matrix representation of term. + + Returns + ------- + ndarray + Matrix representation + """ 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 +323,7 @@ def get_data(self): return data def _check_range(self): + """Check that site indices are within valid range.""" err = 0 for site_info, value in self.term_info.items(): for i in range(2): @@ -203,6 +336,15 @@ def _check_range(self): exit(1) def _check_hermite(self, strict_hermite, tolerance): + """Check Hermiticity of term parameters. + + Parameters + ---------- + strict_hermite : bool + If True, exit on Hermiticity violation + 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,14 +361,41 @@ def _check_hermite(self, strict_hermite, tolerance): logger.warn(msg) def check_hermite(self, strict_hermite=False, tolerance=1.0e-8): + """Public interface for Hermiticity check. + + Parameters + ---------- + strict_hermite : bool, optional + If True, exit on Hermiticity violation + tolerance : float, optional + Tolerance for Hermiticity check + """ self._check_hermite(strict_hermite, tolerance) class Transfer_UHFr(Term_base): + """Hopping term. + + Parameters + ---------- + term_info : dict + Hopping parameters + Nsize : int + System size + """ def __init__(self, term_info, Nsize): self.__name__ = "Transfer" super().__init__(term_info, Nsize, coeff=-1.0) class Green_UHFr(Term_base): + """Green's function term. + + Parameters + ---------- + term_info : dict + Green's function parameters + Nsize : int + System size + """ def __init__(self, term_info, Nsize): self.__name__ = "Green" super().__init__(term_info, Nsize, coeff=1.0) @@ -235,6 +404,40 @@ def __init__(self, term_info, Nsize): from .base import solver_base class UHFr(solver_base): + """Unrestricted Hartree-Fock solver. + + Parameters + ---------- + param_ham : dict + Hamiltonian parameters + info_log : dict + Logging parameters + info_mode : dict + Mode parameters + param_mod : dict, optional + Model parameters + + Attributes + ---------- + name : str + Solver name + physics : dict + Physical quantities + iflag_fock : bool + Include Fock terms if True + ene_cutoff : float + Energy cutoff for finite temperature + T : float + Temperature + strict_hermite : bool + Strict Hermiticity check + hermite_tolerance : float + Tolerance for Hermiticity + Nsize : int + System size + Ncond : int + Number of particles + """ @do_profile def __init__(self, param_ham, info_log, info_mode, param_mod=None): self.name = "uhfr" @@ -275,6 +478,15 @@ def __init__(self, param_ham, info_log, info_mode, param_mod=None): @do_profile def solve(self, green_info, path_to_output): + """Solve the UHF equations. + + Parameters + ---------- + green_info : dict + Green's function parameters + path_to_output : str + Output directory path + """ 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 +539,18 @@ def solve(self, green_info, path_to_output): @do_profile def _initial_G(self, green_info): + """Initialize Green's function. + + Parameters + ---------- + green_info : dict + Green's function parameters + + Returns + ------- + ndarray + Initial 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 +588,7 @@ def _initial_G(self, green_info): @do_profile def _makeham_const(self): + """Initialize constant part of 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 +597,7 @@ def _makeham_const(self): @do_profile def _makeham(self): + """Construct full Hamiltonian.""" import time self.Ham = np.zeros((2 * self.Nsize, 2 * self.Nsize), dtype=complex) self.Ham = self.Ham_trans.copy() @@ -381,6 +607,7 @@ def _makeham(self): @do_profile def _makeham_mat(self): + """Construct interaction matrix.""" # 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: @@ -420,6 +647,7 @@ def _makeham_mat(self): @do_profile def _diag(self): + """Diagonalize Hamiltonian.""" _green_list = self.green_list for k, block_g_info in _green_list.items(): g_label = block_g_info["label"] @@ -437,6 +665,20 @@ def _diag(self): @do_profile def _fermi(self, mu, eigenvalue): + """Calculate Fermi-Dirac distribution. + + Parameters + ---------- + mu : float + Chemical potential + eigenvalue : ndarray + Eigenvalues + + Returns + ------- + ndarray + Fermi-Dirac distribution + """ fermi = np.zeros(eigenvalue.shape) for idx, value in enumerate(eigenvalue): if (value - mu) / self.T > self.ene_cutoff: @@ -447,88 +689,231 @@ def _fermi(self, mu, eigenvalue): @do_profile def _green(self): + """Calculate Green's function. + + Parameters + ---------- + None + + Returns + ------- + None + Updates self.Green and self.Green_old attributes + + Notes + ----- + For T=0: + Constructs Green's function from occupied eigenvectors using: + G = U^* U^T where U contains occupied eigenvectors + + For T>0: + 1. Finds chemical potential mu by solving particle number equation + 2. Calculates Fermi-Dirac occupations + 3. Constructs Green's function using eigenvectors and occupations + """ _green_list = self.green_list # R_SLT = U^{*} in _green # L_SLT = U^T in _green R_SLT = np.zeros((2 * self.Nsize, 2 * self.Nsize), dtype=complex) L_SLT = np.zeros((2 * self.Nsize, 2 * self.Nsize), dtype=complex) - # Copy Green_old + # Store previous Green's function for convergence check self.Green_old = self.Green.copy() - if self.T == 0: + + if self.T == 0: # Zero temperature case for k, block_g_info in _green_list.items(): - g_label = block_g_info["label"] - # occupied: return the occupied number - occupied_number = block_g_info["occupied"] - eigenvec = self.green_list[k]["eigenvector"] - eigen_start = self.green_list[k]["eigen_start"] + g_label = block_g_info["label"] # List of orbital indices for this block + occupied_number = block_g_info["occupied"] # Number of occupied states + eigenvec = self.green_list[k]["eigenvector"] # Eigenvectors for this block + eigen_start = self.green_list[k]["eigen_start"] # Starting index for eigenvalues + + # Construct Green's function from occupied eigenvectors + # G = U^* U^T where U contains occupied eigenvectors for eigen_i in range(occupied_number): - evec_i = eigenvec[:, eigen_i] + evec_i = eigenvec[:, eigen_i] # Get i-th eigenvector for idx, org_site in enumerate(g_label): + # Build U^* and U^T matrices R_SLT[org_site][eigen_start + eigen_i] = np.conjugate(evec_i[idx]) L_SLT[eigen_start + eigen_i][org_site] = evec_i[idx] + + # Multiply matrices to get Green's function RMat = np.dot(R_SLT, L_SLT) self.Green = RMat.copy() - else: # for finite temperatures + + else: # Finite temperature case from scipy import optimize + # Function to find chemical potential by solving particle number equation def _calc_delta_n(mu): + # Calculate occupation numbers for each eigenstate n_eigen = np.einsum("ij, ij -> j", np.conjugate(eigenvec), eigenvec).real + # Get Fermi-Dirac distribution fermi = self._fermi(mu, eigenvalue) + # Return difference from target particle number return np.dot(n_eigen, fermi)-occupied_number + # Initialize Green's function self.Green = np.zeros((2 * self.Nsize, 2 * self.Nsize), dtype=complex) + for k, block_g_info in _green_list.items(): - g_label = block_g_info["label"] - eigenvalue = self.green_list[k]["eigenvalue"] - eigenvec = self.green_list[k]["eigenvector"] - occupied_number = block_g_info["occupied"] + g_label = block_g_info["label"] # Orbital indices + eigenvalue = self.green_list[k]["eigenvalue"] # Eigenvalues + eigenvec = self.green_list[k]["eigenvector"] # Eigenvectors + occupied_number = block_g_info["occupied"] # Target particle number - #determine mu + # Find chemical potential using bisection method first is_converged = False - if (_calc_delta_n(eigenvalue[0]) * _calc_delta_n(eigenvalue[-1])) < 0.0: mu, r = optimize.bisect(_calc_delta_n, eigenvalue[0], eigenvalue[-1], full_output=True, disp=False) is_converged = r.converged + + # If bisection fails, try Newton's method if not is_converged: mu, r = optimize.newton(_calc_delta_n, eigenvalue[0], full_output=True) is_converged = r.converged + + # Exit if chemical potential cannot be found if not is_converged: logger.error("find mu: not converged. abort") exit(1) + # Store chemical potential self.green_list[k]["mu"] = mu + + # Calculate Fermi-Dirac occupations fermi = self._fermi(mu, eigenvalue) + + # Construct Green's function using eigenvectors and occupations + # G = U^* f U where f is diagonal matrix of Fermi-Dirac occupations tmp_green = np.einsum("ij, j, kj -> ik", np.conjugate(eigenvec), fermi, eigenvec) + + # Store block of Green's function 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 total energy and its components. + + Parameters + ---------- + None + + Returns + ------- + None + Updates self.physics["Ene"] with dictionary containing: + - "band": Band energy + - "InterAll": Interaction energy + - "Total": Total energy + """ _green_list = self.green_list Ene = {} Ene["band"] = 0 + # Zero temperature case - sum up energies of occupied states if self.T == 0: - for k, block_g_info in _green_list.items(): - eigenvalue = self.green_list[k]["eigenvalue"] - occupied_number = block_g_info["occupied"] - Ene["band"] += np.sum(eigenvalue[:occupied_number]) + Ene["band"] = _calc_zero_temp_energy(_green_list) else: - - for k, block_g_info in _green_list.items(): - eigenvalue = self.green_list[k]["eigenvalue"] - eigenvec = self.green_list[k]["eigenvector"] - mu = self.green_list[k]["mu"] + Ene["band"] = _calc_finite_temp_energy(_green_list) + + def _calc_zero_temp_energy(self, green_list): + """Calculate band energy for zero temperature case. + + Parameters + ---------- + green_list : dict + Dictionary containing Green's function blocks + + Returns + ------- + float + Band energy at zero temperature + + Notes + ----- + At T=0, band energy is simply sum of occupied eigenvalues. + Loops through each block in Green's function and sums eigenvalues + up to the occupation number for that block. + """ + energy = 0 + # Loop through each block (spin up/down or sz-free) + for k, block_g_info in green_list.items(): + eigenvalue = self.green_list[k]["eigenvalue"] # Get eigenvalues for this block + occupied_number = block_g_info["occupied"] # Number of occupied states + # Sum up eigenvalues of occupied states only + energy += np.sum(eigenvalue[:occupied_number]) + return energy + + def _calc_finite_temp_energy(self, green_list): + """Calculate band energy for finite temperature case. + + Parameters + ---------- + green_list : dict + Dictionary containing Green's function blocks + + Returns + ------- + float + Band energy at finite temperature + + Notes + ----- + At finite T, band energy includes: + 1. Chemical potential term: mu * n where n is particle number + 2. Entropy term: -T * sum(ln(1 + exp(-(e-mu)/T))) + Uses Fermi-Dirac distribution and logarithmic terms. + """ + energy = 0 + # Loop through each block (spin up/down or sz-free) + for k, block_g_info in green_list.items(): + eigenvalue = self.green_list[k]["eigenvalue"] # Get eigenvalues + eigenvec = self.green_list[k]["eigenvector"] # Get eigenvectors + mu = self.green_list[k]["mu"] # Chemical potential + + # Calculate Fermi-Dirac occupations fermi = self._fermi(mu, eigenvalue) - - ln_Ene = np.zeros(eigenvalue.shape) - for idx, value in enumerate(eigenvalue): - if -(value - mu) / self.T < self.ene_cutoff: - ln_Ene[idx] = np.log1p(np.exp(-(value - mu) / self.T)) - else: - ln_Ene[idx] = -(value - mu) / self.T + # Calculate logarithmic terms for entropy + ln_Ene = _calc_log_terms(eigenvalue, mu) + # Calculate particle number using eigenvectors and occupations 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) + + # Add mu*N term and entropy term + energy += mu*np.sum(tmp_n) - self.T * np.sum(ln_Ene) + return energy + + def _calc_log_terms(self, eigenvalue, mu): + """Calculate logarithmic terms in grand potential. + + Parameters + ---------- + eigenvalue : ndarray + Array of eigenvalues + mu : float + Chemical potential + + Returns + ------- + ndarray + Array of logarithmic terms + + Notes + ----- + Calculates ln(1 + exp(-(e-mu)/T)) for each eigenvalue. + Uses cutoff to avoid numerical overflow: + - If -(e-mu)/T < cutoff: use log1p for numerical stability + - If -(e-mu)/T >= cutoff: approximate as -(e-mu)/T + """ + ln_Ene = np.zeros(eigenvalue.shape) + for idx, value in enumerate(eigenvalue): + # Check if exponential will overflow + if -(value - mu) / self.T < self.ene_cutoff: + # Use log1p for numerical stability + ln_Ene[idx] = np.log1p(np.exp(-(value - mu) / self.T)) + else: + # For large negative arguments, approximate as -(e-mu)/T + ln_Ene[idx] = -(value - mu) / self.T + return ln_Ene Ene["InterAll"] = 0 green_local = self.Green.reshape((2 * self.Nsize) ** 2) @@ -543,6 +928,24 @@ def _calc_energy(self): @do_profile def _calc_phys(self): + """Calculate physical observables. + + Parameters + ---------- + None + + Returns + ------- + None + Updates self.physics with: + - "NCond": Total particle number + - "Sz": Total spin z component + - "Rest": Convergence measure + + Notes + ----- + Also performs mixing of old and new Green's functions for convergence + """ n = 0 for site in range(2 * self.Nsize): n += self.Green[site][site] @@ -564,10 +967,33 @@ def _calc_phys(self): @do_profile def get_results(self): + """Get calculation results. + + Returns + ------- + tuple + (physics, Green) where: + - physics: Dictionary of physical observables + - Green: Green's function matrix + """ return (self.physics, self.Green) @do_profile def save_results(self, info_outputfile, green_info): + """Save calculation results to files. + + Parameters + ---------- + info_outputfile : dict + Dictionary specifying output files and paths + green_info : dict + Dictionary containing Green's function information + + Returns + ------- + None + Writes results to files specified in info_outputfile + """ path_to_output = info_outputfile["path_to_output"] if "energy" in info_outputfile.keys(): @@ -662,4 +1088,11 @@ def save_results(self, info_outputfile, green_info): @do_profile def get_Ham(self): + """Get Hamiltonian matrix. + + Returns + ------- + ndarray + Hamiltonian matrix + """ return self.Ham From faf5296fbc8fee0e6b36a934ca5749357f95ef34 Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Mon, 11 Nov 2024 14:28:33 +0900 Subject: [PATCH 5/8] fix --- src/hwave/solver/uhfr.py | 71 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/src/hwave/solver/uhfr.py b/src/hwave/solver/uhfr.py index 03cf0f3..000d0e8 100644 --- a/src/hwave/solver/uhfr.py +++ b/src/hwave/solver/uhfr.py @@ -59,6 +59,77 @@ def _transform_interall(self, ham_info): """ return ham_info + def _calc_hartree(self): + """Calculate Hartree terms in the Hamiltonian. + + Calculates the diagonal Hartree terms in the Hamiltonian by iterating through + interaction parameters and adding contributions to Ham_tmp and Ham_trans_tmp. + + Parameters + ---------- + None + + Returns + ------- + None + Updates self.Ham_tmp and self.Ham_trans_tmp arrays + """ + 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 + # Diagonal Fock term + self.Ham_tmp[site[0]][site[1]][site[2]][site[3]] += value + self.Ham_tmp[site[2]][site[3]][site[0]][site[1]] += value + if site[1] == site[2]: + self.Ham_trans_tmp[site[1]][site[2]] += value + pass + + def _calc_fock(self): + """Calculate Fock exchange terms in the Hamiltonian. + + Calculates the off-diagonal Fock exchange terms in the Hamiltonian by iterating + through interaction parameters and adding contributions to Ham_tmp. + + Parameters + ---------- + None + + Returns + ------- + None + Updates self.Ham_tmp array + """ + 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 + # OffDiagonal Fock term + self.Ham_tmp[site[0]][site[3]][site[2]][site[1]] -= value + self.Ham_tmp[site[2]][site[1]][site[0]][site[3]] -= value + pass + + def get_ham(self, type): + """Get the Hamiltonian matrices. + + Calculates and returns the Hartree and Fock terms of the Hamiltonian. + + Parameters + ---------- + type : str + Type of calculation - either "hartree" or "hartreefock" + + Returns + ------- + tuple of ndarray + (Ham_tmp, Ham_trans_tmp) containing the Hamiltonian matrices + """ + self._calc_hartree() + if type == "hartreefock": + self._calc_fock() + return self.Ham_tmp, self.Ham_trans_tmp + + def _check_range(self): """Check that site indices are within valid range.""" err = 0 From a8f64349dad70cc5497bfce46e243c662b84b32c Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Mon, 11 Nov 2024 14:33:33 +0900 Subject: [PATCH 6/8] change order _calc_**_energy functions --- src/hwave/solver/uhfr.py | 84 +++++++++++++++++++++------------------- 1 file changed, 44 insertions(+), 40 deletions(-) diff --git a/src/hwave/solver/uhfr.py b/src/hwave/solver/uhfr.py index 000d0e8..6c65032 100644 --- a/src/hwave/solver/uhfr.py +++ b/src/hwave/solver/uhfr.py @@ -881,12 +881,6 @@ def _calc_energy(self): _green_list = self.green_list Ene = {} Ene["band"] = 0 - # Zero temperature case - sum up energies of occupied states - if self.T == 0: - Ene["band"] = _calc_zero_temp_energy(_green_list) - else: - Ene["band"] = _calc_finite_temp_energy(_green_list) - def _calc_zero_temp_energy(self, green_list): """Calculate band energy for zero temperature case. @@ -915,6 +909,39 @@ def _calc_zero_temp_energy(self, green_list): energy += np.sum(eigenvalue[:occupied_number]) return energy + def _calc_log_terms(self, eigenvalue, mu): + """Calculate logarithmic terms in grand potential. + + Parameters + ---------- + eigenvalue : ndarray + Array of eigenvalues + mu : float + Chemical potential + + Returns + ------- + ndarray + Array of logarithmic terms + + Notes + ----- + Calculates ln(1 + exp(-(e-mu)/T)) for each eigenvalue. + Uses cutoff to avoid numerical overflow: + - If -(e-mu)/T < cutoff: use log1p for numerical stability + - If -(e-mu)/T >= cutoff: approximate as -(e-mu)/T + """ + ln_Ene = np.zeros(eigenvalue.shape) + for idx, value in enumerate(eigenvalue): + # Check if exponential will overflow + if -(value - mu) / self.T < self.ene_cutoff: + # Use log1p for numerical stability + ln_Ene[idx] = np.log1p(np.exp(-(value - mu) / self.T)) + else: + # For large negative arguments, approximate as -(e-mu)/T + ln_Ene[idx] = -(value - mu) / self.T + return ln_Ene + def _calc_finite_temp_energy(self, green_list): """Calculate band energy for finite temperature case. @@ -942,6 +969,7 @@ def _calc_finite_temp_energy(self, green_list): eigenvec = self.green_list[k]["eigenvector"] # Get eigenvectors mu = self.green_list[k]["mu"] # Chemical potential + # Calculate Fermi-Dirac occupations fermi = self._fermi(mu, eigenvalue) # Calculate logarithmic terms for entropy @@ -951,40 +979,16 @@ def _calc_finite_temp_energy(self, green_list): # Add mu*N term and entropy term energy += mu*np.sum(tmp_n) - self.T * np.sum(ln_Ene) - return energy - - def _calc_log_terms(self, eigenvalue, mu): - """Calculate logarithmic terms in grand potential. - - Parameters - ---------- - eigenvalue : ndarray - Array of eigenvalues - mu : float - Chemical potential - - Returns - ------- - ndarray - Array of logarithmic terms - - Notes - ----- - Calculates ln(1 + exp(-(e-mu)/T)) for each eigenvalue. - Uses cutoff to avoid numerical overflow: - - If -(e-mu)/T < cutoff: use log1p for numerical stability - - If -(e-mu)/T >= cutoff: approximate as -(e-mu)/T - """ - ln_Ene = np.zeros(eigenvalue.shape) - for idx, value in enumerate(eigenvalue): - # Check if exponential will overflow - if -(value - mu) / self.T < self.ene_cutoff: - # Use log1p for numerical stability - ln_Ene[idx] = np.log1p(np.exp(-(value - mu) / self.T)) - else: - # For large negative arguments, approximate as -(e-mu)/T - ln_Ene[idx] = -(value - mu) / self.T - return ln_Ene + return energy + + + # Zero temperature case - sum up energies of occupied states + if self.T == 0: + Ene["band"] = _calc_zero_temp_energy(_green_list) + else: + Ene["band"] = _calc_finite_temp_energy(_green_list) + + Ene["InterAll"] = 0 green_local = self.Green.reshape((2 * self.Nsize) ** 2) From bab66fbead6e9d1899ce82a6ffa99520667d439b Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Mon, 11 Nov 2024 14:41:16 +0900 Subject: [PATCH 7/8] bugfix --- src/hwave/solver/uhfr.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/hwave/solver/uhfr.py b/src/hwave/solver/uhfr.py index 6c65032..9985567 100644 --- a/src/hwave/solver/uhfr.py +++ b/src/hwave/solver/uhfr.py @@ -881,7 +881,8 @@ def _calc_energy(self): _green_list = self.green_list Ene = {} Ene["band"] = 0 - def _calc_zero_temp_energy(self, green_list): + + def _calc_zero_temp_energy(green_list): """Calculate band energy for zero temperature case. Parameters @@ -909,7 +910,7 @@ def _calc_zero_temp_energy(self, green_list): energy += np.sum(eigenvalue[:occupied_number]) return energy - def _calc_log_terms(self, eigenvalue, mu): + def _calc_log_terms(eigenvalue, mu): """Calculate logarithmic terms in grand potential. Parameters @@ -942,7 +943,7 @@ def _calc_log_terms(self, eigenvalue, mu): ln_Ene[idx] = -(value - mu) / self.T return ln_Ene - def _calc_finite_temp_energy(self, green_list): + def _calc_finite_temp_energy(green_list): """Calculate band energy for finite temperature case. Parameters From 628b772064a27874d085de249e2669470ec64bcd Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Mon, 11 Nov 2024 14:55:18 +0900 Subject: [PATCH 8/8] Update docstrings to numpydoc format in read_input.py, read_input_k.py, and wan90.py --- src/hwave/qlmsio/read_input_k.py | 103 ++++++++++++++++++++ src/hwave/qlmsio/wan90.py | 158 +++++++++++++++++++++++++++---- 2 files changed, 245 insertions(+), 16 deletions(-) diff --git a/src/hwave/qlmsio/read_input_k.py b/src/hwave/qlmsio/read_input_k.py index fd5ae51..a3cf15f 100644 --- a/src/hwave/qlmsio/read_input_k.py +++ b/src/hwave/qlmsio/read_input_k.py @@ -1,3 +1,15 @@ +"""Input file reader for k-space calculations. + +This module provides functionality to read and parse input files for k-space +calculations in QLMS. It handles geometry, transfer integrals, and various +interaction terms. + +Classes +------- +QLMSkInput + Reads and parses input files for k-space calculations. +""" + import logging import numpy as np import sys @@ -10,6 +22,46 @@ class QLMSkInput(): + """Input file reader for k-space calculations. + + Parameters + ---------- + info_inputfile : dict + Dictionary containing input file information with keys: + - path_to_input : str + Path to input file directory + - interaction : dict + Dictionary specifying interaction files + - initial : str, optional + Filename for initial Green's function in k-space (NPZ format) + - initial_uhf : str, optional + Filename for initial Green's function in real space + - geometry_uhf : str, optional + Filename for geometry in real space + - onebodyg_uhf : str, optional + Filename for one-body Green's function indices + solver_type : str, optional + Type of solver to use (default: "UHFk") + + Attributes + ---------- + valid_namelist : list + List of valid input file section names + ham_param : CaseInsensitiveDict + Dictionary storing Hamiltonian parameters + green : CaseInsensitiveDict + Dictionary storing Green's function data + + Methods + ------- + get_param(key) + Get parameters by key + _read_data(file_name, value_type) + Read data from file into dictionary + _read_green(file_name) + Read Green's function indices from file + """ + 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"): @@ -91,6 +143,25 @@ 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 file into dictionary. + + Parameters + ---------- + file_name : str + Name of file to read + value_type : str, optional + Type of values to read ("real" or "complex") + + Returns + ------- + dict + Dictionary containing data read from file + + Raises + ------ + FileNotFoundError + If input file is not found + """ info = {} try: data = np.loadtxt(file_name, skiprows = 5) @@ -110,6 +181,23 @@ def _read_data(self, file_name, value_type="real"): return info def _read_green(self, file_name): + """Read Green's function indices from file. + + Parameters + ---------- + file_name : str + Name of file to read + + Returns + ------- + ndarray + Array of Green's function indices + + Raises + ------ + FileNotFoundError + If input file is not found + """ try: _data = np.loadtxt(file_name, dtype=np.int32, skiprows = 5) except FileNotFoundError: @@ -118,6 +206,21 @@ def _read_green(self, file_name): return _data def get_param(self, key): + """Get parameters by key. + + Parameters + ---------- + key : str + Key to retrieve parameters for: + - "mod"/"parameter": Returns None + - "ham"/"hamiltonian": Returns Hamiltonian parameters + - "output"/"green": Returns Green's function data + + Returns + ------- + CaseInsensitiveDict or None + Requested parameters or None if 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..5a13889 100644 --- a/src/hwave/qlmsio/wan90.py +++ b/src/hwave/qlmsio/wan90.py @@ -1,3 +1,22 @@ +"""Functions for reading/writing Wannier90 format files. + +This module provides functions to read and write geometry and Hamiltonian data +in Wannier90 format. + +Functions +--------- +read_geom(name_in) + Read geometry data from file. +read_geometry(name_in) + Read extended geometry data from file. +write_geom(name_out, info_geometry) + Write geometry data to file. +read_w90(name_in) + Read Wannier90 Hamiltonian data from file. +write_w90(name_in, info_int, info_geometry, interact_shape) + Write Wannier90 Hamiltonian data to file. +""" + from __future__ import print_function import itertools @@ -7,6 +26,29 @@ logger = logging.getLogger("qlms").getChild("wan90") def read_geom(name_in): + """Read geometry data from file. + + Parameters + ---------- + name_in : str + Input filename + + Returns + ------- + dict + Dictionary containing: + - norb : int + Number of orbitals + - rvec : ndarray + Real space lattice vectors (3x3) + - center : ndarray + Orbital center positions (norb x 3) + + Raises + ------ + SystemExit + If file not found + """ try: with open(name_in, 'r') as f: # skip header @@ -31,6 +73,33 @@ def read_geom(name_in): return data def read_geometry(name_in): + """Read extended geometry data from file. + + Parameters + ---------- + name_in : str + Input filename + + Returns + ------- + dict + Dictionary containing: + - unit_vec : ndarray + Unit cell vectors + - degree : ndarray + Degrees of freedom + - cell_vec : ndarray + Cell vectors + - site2vec : dict + Mapping of sites to vectors + - n_orb : int + Number of orbitals + + Raises + ------ + SystemExit + If file not found + """ info_geometry = {} unit_vec_line_end = 3 cell_vec_line_start = 4 @@ -62,6 +131,15 @@ def read_geometry(name_in): return info_geometry def write_geom(name_out, info_geometry): + """Write geometry data to file. + + Parameters + ---------- + name_out : str + Output filename + info_geometry : dict + Geometry information dictionary + """ with open(name_out, "w") as fw: #Unit_cell for vec in info_geometry["unit_vec"]: @@ -75,6 +153,27 @@ def write_geom(name_out, info_geometry): fw.write("{} {} {}\n".format(pos, pos, pos)) def read_w90(name_in): + """Read Wannier90 Hamiltonian data from file. + + Parameters + ---------- + name_in : str + Input filename + + Returns + ------- + dict + Dictionary mapping (irvec, orbvec) tuples to complex values where: + - irvec : tuple + Real space lattice vector indices + - orbvec : tuple + Orbital indices + + Raises + ------ + SystemExit + If file not found + """ try: with open(name_in, 'r') as f: # skip header @@ -104,23 +203,50 @@ def read_w90(name_in): return data def write_w90(name_in, info_int, info_geometry, interact_shape): + """Write Wannier90 Hamiltonian data to file. + + Parameters + ---------- + name_in : str + Output filename + info_int : dict + Interaction information dictionary + info_geometry : dict + Geometry information dictionary + interact_shape : tuple + Shape of interaction array + """ + # Build output content as a list of strings + output = [] + + # Header information + norb = info_int["n_orb"] + Nlattice = np.prod(interact_shape) + n_total = (Nlattice * norb) ** 2 + + output.append("# wannier90 format") + output.append(str(norb)) + output.append(str(n_total)) + + # Write list of ones, 15 per line + ones = np.ones(n_total, dtype=int) + for i in range(0, n_total, 15): + output.append(" ".join(map(str, ones[i:i+15]))) + + # Interaction data + fmt = "{} {} {} {} {} {} {}" + for idx, (site_org, int_value) in info_int["interaction"].items(): + site = info_geometry["site2vec"][idx] + output.append(fmt.format( + site_org, + site[0], site[1], + site[2]+1, site[3]+1, + int_value.real, int_value.imag + )) + + # Write all output at once with open(name_in, "w") as fw: - fw.write("# wannier90 format\n") - norb = info_int["n_orb"] - fw.write("{}\n".format(norb)) - Nlattice = interact_shape[0]*interact_shape[1]*interact_shape[2] - fw.write("{}\n".format((Nlattice*norb)**2)) - for idx in range((Nlattice*norb)**2): - fw.write("1 ") - if (idx+1)%15 == 0: - fw.write("\n") - if (idx+1)%15 != 0: - fw.write("\n") - for idx, value in info_int["interaction"].items(): - site = info_geometry["site2vec"][idx] - site_org = value[0] - 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)) + fw.write("\n".join(output) + "\n") if __name__ == "__main__": path_to_sample = "../../sample/dir-model"