diff --git a/README.md b/README.md index 186d5bc..a36b258 100644 --- a/README.md +++ b/README.md @@ -3,21 +3,14 @@ More on the theoretical background can be seen on [arXiv](https://arxiv.org/abs/ # TODO ## Testing -- Calculate total charge for all Green's function --> - Run tests on different magnetic materials and compare it to Grogu Matlab --> ran on Jij_for_Marci_6p45ang, but I could not compare data -- Change order in: - * magnetic entities --> - * pairs --> -- Compair single pair run to multiple pair run --> it is the same -- Check behavior by distance to test symmetry and Green function phase --> it is symmetric and decays by distance in an AFM way -- Check behavior by orbital indexing --> there is no unwanted mixing and otherwise consistent with symmetry and distance ## Developing -- Add documentation to `useful` functions +- Add documentation to `useful` functions and more... - Check the symmetrization of the Hamiltonian and overlap matrix to make them hermitian - Check if exchange field has scalar part -- Output dictionary and rewrite output printing with the information - Add more tests!! +- io stuff # Building wheel See detailed documentation on [PYPI](https://packaging.python.org/en/latest/tutorials/packaging-projects/). diff --git a/src/grogu_magn/grogu.py b/src/grogu_magn/grogu.py index 93d0ca4..b01b08c 100644 --- a/src/grogu_magn/grogu.py +++ b/src/grogu_magn/grogu.py @@ -17,3 +17,480 @@ # LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. + +# Copyright (c) [2024] [Daniel Pozsar] +# +# Permission is hereby granted, free of charge, to any person obtaining a copy +# of this software and associated documentation files (the "Software"), to deal +# in the Software without restriction, including without limitation the rights +# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +# copies of the Software, and to permit persons to whom the Software is +# furnished to do so, subject to the following conditions: +# +# The above copyright notice and this permission notice shall be included in all +# copies or substantial portions of the Software. +# +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +# SOFTWARE. + + +import os + +os.environ["OMP_NUM_THREADS"] = "1" # export OMP_NUM_THREADS=1 +os.environ["OPENBLAS_NUM_THREADS"] = "1" # export OPENBLAS_NUM_THREADS=1 +os.environ["MKL_NUM_THREADS"] = "1" # export MKL_NUM_THREADS=1 +os.environ["VECLIB_MAXIMUM_THREADS"] = "1" # export VECLIB_MAXIMUM_THREADS=1 +os.environ["NUMEXPR_NUM_THREADS"] = "1" # export NUMEXPR_NUM_THREADS=1 + + +from timeit import default_timer as timer + +# runtime information +times = dict() +times["start_time"] = timer() + + +import warnings +from sys import getsizeof + +import sisl +from mpi4py import MPI + +from src.grogu_magn import * + +# input output stuff +###################################################################### +###################################################################### +###################################################################### + +infile = ( + "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf" +) +outfile = "./Fe3GeTe2_notebook" + +magnetic_entities = [ + dict(atom=3, l=2), + dict(atom=4, l=2), + dict(atom=5, l=2), +] +pairs = [ + dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])), + dict(ai=0, aj=2, Ruc=np.array([-1, -1, 0])), + dict(ai=1, aj=2, Ruc=np.array([-1, -1, 0])), + dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-2, 0, 0])), + dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])), +] +simulation_parameters = default_args +simulation_parameters["infile"] = infile +simulation_parameters["outfile"] = outfile +simulation_parameters["kset"] = 20 +simulation_parameters["kdirs"] = "xy" +simulation_parameters["eset"] = 600 +simulation_parameters["esetp"] = 10000 + +fdf = sisl.io.fdfSileSiesta("input.fdf") +fdf.get("XCF_Rotation") + +###################################################################### +###################################################################### +###################################################################### + + +# MPI parameters +comm = MPI.COMM_WORLD +size = comm.Get_size() +rank = comm.Get_rank() +root_node = 0 + +# include parallel size in simulation parameters +simulation_parameters["parallel_size"] = size + +# check versions for debugging +if rank == root_node: + try: + print("sisl version: ", sisl.__version__) + except: + print("sisl version unknown.") + try: + print("numpy version: ", np.__version__) + except: + print("numpy version unknown.") + +# rename outfile +if not simulation_parameters["outfile"].endswith(".pickle"): + simulation_parameters["outfile"] += ".pickle" + +# if ebot is not given put it 0.1 eV under the smallest energy +if simulation_parameters["ebot"] is None: + try: + eigfile = simulation_parameters["infile"][:-3] + "EIG" + simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 0.1 + simulation_parameters["automatic_ebot"] = True + except: + print("Could not determine ebot.") + print("Parameter was not given and .EIG file was not found.") + +# read sile +fdf = sisl.get_sile(simulation_parameters["infile"]) + +# read in hamiltonian +dh = fdf.read_hamiltonian() + +# read unit cell vectors +simulation_parameters["cell"] = fdf.read_geometry().cell + +# unit cell index +uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0]) + +if rank == root_node: + print("\n\n\n\n\n") + print( + "#################################################################### JOB INFORMATION ###########################################################################" + ) + print_job_description(simulation_parameters) + print( + "################################################################################################################################################################" + ) + print("\n\n\n\n\n") + + times["setup_time"] = timer() + print(f"Setup done. Elapsed time: {times['setup_time']} s") + print( + "================================================================================================================================================================" + ) + + +# reformat Hamltonian and Overlap matrix for manipulations +hh, ss, NO = build_hh_ss(dh) + +# symmetrizing Hamiltonian and Overlap matrix to make them hermitian +for i in range(dh.lattice.sc_off.shape[0]): + j = dh.lattice.sc_index(-dh.lattice.sc_off[i]) + h1, h1d = hh[i], hh[j] + hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2 + s1, s1d = ss[i], ss[j] + ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2 + +# identifying TRS and TRB parts of the Hamiltonian +TAUY = np.kron(np.eye(NO), tau_y) +hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())]) +hTRS = (hh + hTR) / 2 +hTRB = (hh - hTR) / 2 + +# extracting the exchange field +traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77 +XCF = np.array( + [ + np.array([f["x"] / 2 for f in traced]), + np.array([f["y"] / 2 for f in traced]), + np.array([f["z"] / 2 for f in traced]), + ] +) + +# check if exchange field has scalar part +max_xcfs = abs(np.array(np.array([f["c"] / 2 for f in traced]))).max() +if max_xcfs > 1e-12: + warnings.warn( + f"Exchange field has non negligible scalar part. Largest value is {max_xcfs}" + ) + +if rank == root_node: + times["H_and_XCF_time"] = timer() + print( + f"Hamiltonian and exchange field rotated. Elapsed time: {times['H_and_XCF_time']} s" + ) + print( + "================================================================================================================================================================" + ) + +# initialize pairs and magnetic entities based on input information +pairs, magnetic_entities = setup_pairs_and_magnetic_entities( + magnetic_entities, pairs, dh, simulation_parameters +) + +if rank == root_node: + times["site_and_pair_dictionaries_time"] = timer() + print( + f"Site and pair dictionaries created. Elapsed time: {times['site_and_pair_dictionaries_time']} s" + ) + print( + "================================================================================================================================================================" + ) + +# generate k space sampling +kset = make_kset( + dirs=simulation_parameters["kdirs"], NUMK=simulation_parameters["kset"] +) + +# generate weights for k points +wkset = np.ones(len(kset)) / len(kset) + +# split the k points based on MPI size +kpcs = np.array_split(kset, size) + +# use progress bar if available +if rank == root_node and tqdm_imported: + kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop") + +if rank == root_node: + times["k_set_time"] = timer() + print(f"k set created. Elapsed time: {times['k_set_time']} s") + print( + "================================================================================================================================================================" + ) + +# this will contain the three Hamiltonian in the +# reference directions needed to calculate the energy +# variations upon rotation +hamiltonians = [] + +# iterate over the reference directions (quantization axes) +for i, orient in enumerate(simulation_parameters["ref_xcf_orientations"]): + # obtain rotated exchange field and Hamiltonian + R = RotMa2b(simulation_parameters["scf_xcf_orientation"], orient["o"]) + rot_XCF = np.einsum("ij,jklm->iklm", R, XCF) + rot_H_XCF = sum( + [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])] + ) + rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx] + + # obtain total Hamiltonian with the rotated exchange field + rot_H = hTRS + rot_H_XCF # equation 76 + + # store the relevant information of the Hamiltonian + hamiltonians.append(dict(orient=orient["o"], H=rot_H)) + + if simulation_parameters["calculate_charge"]: + hamiltonians[-1]["GS"] = np.zeros( + (simulation_parameters["eset"], rot_H.shape[1], rot_H.shape[2]), + dtype="complex128", + ) + hamiltonians[-1]["GS_tmp"] = np.zeros( + (simulation_parameters["eset"], rot_H.shape[1], rot_H.shape[2]), + dtype="complex128", + ) + + # these are the rotations (for now) perpendicular to the quantization axis + for u in orient["vw"]: + # section 2.H + Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) + Vu1, Vu2 = calc_Vu(rot_H_XCF_uc, Tu) + + for mag_ent in magnetic_entities: + idx = mag_ent["spin_box_indices"] + # fill up the perturbed potentials (for now) based on the on-site projections + mag_ent["Vu1"][i].append(onsite_projection(Vu1, idx, idx)) + mag_ent["Vu2"][i].append(onsite_projection(Vu2, idx, idx)) + +if rank == root_node: + times["reference_rotations_time"] = timer() + print( + f"Rotations done perpendicular to quantization axis. Elapsed time: {times['reference_rotations_time']} s" + ) + print( + "================================================================================================================================================================" + ) + +# provide helpful information to estimate the runtime and memory +# requirements of the Greens function calculations +if rank == root_node: + print("Starting matrix inversions.") + if simulation_parameters["padawan_mode"]: + print("Padawan mode: ") + print(f"Total number of k points: {kset.shape[0]}") + print(f"Number of energy samples per k point: {simulation_parameters['eset']}") + print(f"Total number of directions: {len(hamiltonians)}") + print( + f"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * simulation_parameters['eset']}" + ) + print( + f"The shape of the Hamiltonian and the Greens function is {NO}x{NO}={NO*NO}" + ) + # https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv + # memory is O(64 n**2) for complex matrices + memory_size = getsizeof(hamiltonians[0]["H"].base) / 1024 + print( + f"Memory taken by a single Hamiltonian is: {getsizeof(hamiltonians[0]['H'].base) / 1024} KB" + ) + print(f"Expected memory usage per matrix inversion: {memory_size * 32} KB") + print( + f"Expected memory usage per k point for parallel inversion: {memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32} KB" + ) + print( + f"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * simulation_parameters['eset'] * 32 / 1024} MB" + ) + print( + "================================================================================================================================================================" + ) + +# MPI barrier +comm.Barrier() + +# make energy contour +cont = make_contour( + emin=simulation_parameters["ebot"], + enum=simulation_parameters["eset"], + p=simulation_parameters["esetp"], +) +eran = cont.ze + +# sampling the integrand on the contour and the BZ +for k in kpcs[rank]: + # weight of k point in BZ integral + wk = wkset[rank] + + # iterate over reference directions + for i, hamiltonian_orientation in enumerate(hamiltonians): + # calculate Hamiltonian and Overlap matrix in a given k point + H = hamiltonian_orientation["H"] + HK, SK = hsk(H, ss, dh.sc_off, k) + + if simulation_parameters["parallel_solver_for_Gk"]: + Gk = parallel_Gk(HK, SK, eran, simulation_parameters["eset"]) + else: + # solve Greens function sequentially for the energies, because of memory bound + Gk = sequential_GK(HK, SK, eran, simulation_parameters["eset"]) + + # saving this for total charge + if simulation_parameters["calculate_charge"]: + hamiltonian_orientation["GS_tmp"] += Gk @ SK * wk + + # store the Greens function slice of the magnetic entities + for mag_ent in magnetic_entities: + idx = mag_ent["spin_box_indices"] + mag_ent["Gii_tmp"][i] += onsite_projection(Gk, idx, idx) * wk + + for pair in pairs: + # add phase shift based on the cell difference + phase = np.exp(1j * 2 * np.pi * k @ pair["Ruc"].T) + + # get the pair orbital sizes from the magnetic entities + ai = magnetic_entities[pair["ai"]]["spin_box_indices"] + aj = magnetic_entities[pair["aj"]]["spin_box_indices"] + + # store the Greens function slice of the magnetic entities + pair["Gij_tmp"][i] += onsite_projection(Gk, ai, aj) * phase * wk + pair["Gji_tmp"][i] += onsite_projection(Gk, aj, ai) / phase * wk + +# summ reduce partial results of mpi nodes +for i in range(len(hamiltonians)): + # for total charge calculation + if simulation_parameters["calculate_charge"]: + comm.Reduce(hamiltonians[i]["GS_tmp"], hamiltonians[i]["GS"], root=root_node) + + for mag_ent in magnetic_entities: + comm.Reduce(mag_ent["Gii_tmp"][i], mag_ent["Gii"][i], root=root_node) + + for pair in pairs: + comm.Reduce(pair["Gij_tmp"][i], pair["Gij"][i], root=root_node) + comm.Reduce(pair["Gji_tmp"][i], pair["Gji"][i], root=root_node) + +if rank == root_node: + times["green_function_inversion_time"] = timer() + print( + f"Calculated Greens functions. Elapsed time: {times['green_function_inversion_time']} s" + ) + print( + "================================================================================================================================================================" + ) + + +if rank == root_node: + # Calculate total charge + if simulation_parameters["calculate_charge"]: + for hamiltonian in hamiltonians: + GS = hamiltonian["GS"] + traced = np.trace((GS), axis1=1, axis2=2) + simulation_parameters["charges"].append(int_de_ke(traced, cont.we)) + + # iterate over the magnetic entities + for tracker, mag_ent in enumerate(magnetic_entities): + # iterate over the quantization axes + for i, Gii in enumerate(mag_ent["Gii"]): + storage = [] + # iterate over the first and second order local perturbations + for Vu1, Vu2 in zip(mag_ent["Vu1"][i], mag_ent["Vu2"][i]): + # The Szunyogh-Lichtenstein formula + traced = np.trace( + (Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2 + ) # this is the on site projection + # evaluation of the contour integral + storage.append(int_de_ke(traced, cont.we)) + # fill up the magnetic entities dictionary with the energies + magnetic_entities[tracker]["energies"].append(storage) + # convert to np array + magnetic_entities[tracker]["energies"] = np.array( + magnetic_entities[tracker]["energies"] + ) + + # iterate over the pairs + for tracker, pair in enumerate(pairs): + # iterate over the quantization axes + for i, (Gij, Gji) in enumerate(zip(pair["Gij"], pair["Gji"])): + site_i = magnetic_entities[pair["ai"]] + site_j = magnetic_entities[pair["aj"]] + storage = [] + # iterate over the first order local perturbations in all possible orientations for the two sites + for Vui in site_i["Vu1"][i]: + for Vuj in site_j["Vu1"][i]: + # The Szunyogh-Lichtenstein formula + traced = np.trace( + (Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2 + ) # this is the on site projection + # evaluation of the contour integral + storage.append(int_de_ke(traced, cont.we)) + # fill up the pairs dictionary with the energies + pairs[tracker]["energies"].append(storage) + # convert to np array + pairs[tracker]["energies"] = np.array(pairs[tracker]["energies"]) + + # calculate magnetic parameters + for mag_ent in magnetic_entities: + Kxx, Kyy, Kzz, consistency = calculate_anisotropy_tensor(mag_ent) + mag_ent["K"] = np.array([Kxx, Kyy, Kzz]) * sisl.unit_convert("eV", "meV") + mag_ent["K_consistency"] = consistency + + for pair in pairs: + J_iso, J_S, D, J = calculate_exchange_tensor(pair) + pair["J_iso"] = J_iso * sisl.unit_convert("eV", "meV") + pair["J_S"] = J_S * sisl.unit_convert("eV", "meV") + pair["D"] = D * sisl.unit_convert("eV", "meV") + pair["J"] = J * sisl.unit_convert("eV", "meV") + + times["end_time"] = timer() + print("\n\n\n\n\n") + print( + "##################################################################### GROGU OUTPUT #############################################################################" + ) + + print_parameters(simulation_parameters) + print_atoms_and_pairs(magnetic_entities, pairs) + print( + "################################################################################################################################################################" + ) + print_runtime_information(times) + print("") + + # remove unwanted stuff before saving + pairs, magnetic_entities = remove_clutter_for_save(pairs, magnetic_entities) + # create output dictionary with all the relevant data + results = dict( + parameters=simulation_parameters, + magnetic_entities=magnetic_entities, + pairs=pairs, + runtime=times, + ) + + # save results + save_pickle(simulation_parameters["outfile"], results) + + print("\n\n\n\n\n")