From 49d2be4982ef72e5ec2b801c182ac6d885b1775d Mon Sep 17 00:00:00 2001 From: Daniel Pozsar Date: Fri, 8 Nov 2024 11:34:19 +0100 Subject: [PATCH] restructuring code --- src/grogu_magn/__init__.py | 13 ++ src/grogu_magn/core.py | 263 ++++++++++++++++++++++++++++++++++++ src/grogu_magn/grogu.py | 19 +++ src/grogu_magn/io.py | 219 ++++++++++++++++++++++++++++++ src/grogu_magn/magnetism.py | 88 ++++++++++++ src/grogu_magn/utils.py | 232 +------------------------------ 6 files changed, 607 insertions(+), 227 deletions(-) create mode 100644 src/grogu_magn/core.py create mode 100644 src/grogu_magn/io.py create mode 100644 src/grogu_magn/magnetism.py diff --git a/src/grogu_magn/__init__.py b/src/grogu_magn/__init__.py index 93d0ca4..8548e29 100644 --- a/src/grogu_magn/__init__.py +++ b/src/grogu_magn/__init__.py @@ -17,3 +17,16 @@ # 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. + +from grogu_magn.core import * +from grogu_magn.io import * +from grogu_magn.magnetism import * +from grogu_magn.utils import * + +try: + from tqdm import tqdm + + tqdm_imported = True +except: + print("Please install tqdm for nice progress bar.") + tqdm_imported = False diff --git a/src/grogu_magn/core.py b/src/grogu_magn/core.py new file mode 100644 index 0000000..c294e06 --- /dev/null +++ b/src/grogu_magn/core.py @@ -0,0 +1,263 @@ +import numpy as np +from numpy.linalg import inv + +from grogu_magn.utils import blow_up_orbindx, commutator, parse_magnetic_entity + + +def parallel_Gk(HK, SK, eran, eset): + return inv(SK * eran.reshape(eset, 1, 1) - HK) + + +def sequential_GK(HK, SK, eran, eset): + Gk = np.zeros(shape=(eset, HK.shape[0], HK.shape[1]), dtype="complex128") + for j in range(eset): + Gk[j] = inv(SK * eran[j] - HK) + + return Gk + + +def calc_Vu(H, Tu): + """_summary_ + + Args: + H (_type_): _description_ + Tu (_type_): _description_ + + Returns: + _type_: _description_ + """ + + Vu1 = 1j / 2 * commutator(H, Tu) # equation 100 + Vu2 = 1 / 8 * commutator(commutator(Tu, H), Tu) # equation 100 + + return Vu1, Vu2 + + +def remove_clutter_for_save(pairs, magnetic_entities): + """_summary_ + + Args: + pairs (_type_): _description_ + magnetic_entities (_type_): _description_ + + Returns: + _type_: _description_ + """ + # remove clutter from magnetic entities and pair information + for pair in pairs: + del pair["Gij"] + del pair["Gij_tmp"] + del pair["Gji"] + del pair["Gji_tmp"] + for mag_ent in magnetic_entities: + del mag_ent["Gii"] + del mag_ent["Gii_tmp"] + del mag_ent["Vu1"] + del mag_ent["Vu2"] + + return pairs, magnetic_entities + + +def build_hh_ss(dh): + """_summary_ + + Args: + dh (_type_): _description_ + + Returns: + _type_: _description_ + """ + NO = dh.no # shorthand for number of orbitals in the unit cell + + # preprocessing Hamiltonian and overlap matrix elements + h11 = dh.tocsr(dh.M11r) + h11 += dh.tocsr(dh.M11i) * 1.0j + h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") + + h22 = dh.tocsr(dh.M22r) + h22 += dh.tocsr(dh.M22i) * 1.0j + h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") + + h12 = dh.tocsr(dh.M12r) + h12 += dh.tocsr(dh.M12i) * 1.0j + h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") + + h21 = dh.tocsr(dh.M21r) + h21 += dh.tocsr(dh.M21i) * 1.0j + h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype("complex128") + + sov = ( + dh.tocsr(dh.S_idx) + .toarray() + .reshape(NO, dh.n_s, NO) + .transpose(0, 2, 1) + .astype("complex128") + ) + + # Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation + U = np.vstack( + [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])] + ) + # This is the permutation that transforms ud1ud2 to u12d12 + # That is this transforms FROM SPIN BOX to ORBITAL BOX => U + # the inverse transformation is U.T u12d12 to ud1ud2 + # That is FROM ORBITAL BOX to SPIN BOX => U.T + + # From now on everything is in SPIN BOX!! + hh, ss = np.array( + [ + U.T + @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) + @ U + for i in range(dh.lattice.nsc.prod()) + ] + ), np.array( + [ + U.T + @ np.block( + [[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]] + ) + @ U + for i in range(dh.lattice.nsc.prod()) + ] + ) + + return hh, ss, NO + + +def setup_pairs_and_magnetic_entities( + magnetic_entities, pairs, dh, simulation_parameters +): + """_summary_ + + Args: + magnetic_entities (_type_): _description_ + pairs (_type_): _description_ + dh (_type_): _description_ + simulation_parameters (_type_): _description_ + + Returns: + _type_: _description_ + """ + # for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions + for mag_ent in magnetic_entities: + parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes + mag_ent["orbital_indeces"] = parsed + mag_ent["spin_box_indeces"] = blow_up_orbindx( + parsed + ) # calculate spin box indexes + # if orbital is not set use all + if "l" not in mag_ent.keys(): + mag_ent["l"] = "all" + if isinstance(mag_ent["atom"], int): + mag_ent["tags"] = [ + f"[{mag_ent['atom']}]{dh.atoms[mag_ent['atom']].tag}({mag_ent['l']})" + ] + mag_ent["xyz"] = [dh.xyz[mag_ent["atom"]]] + if isinstance(mag_ent["atom"], list): + mag_ent["tags"] = [] + mag_ent["xyz"] = [] + # iterate over atoms + for atom_idx in mag_ent["atom"]: + mag_ent["tags"].append( + f"[{atom_idx}]{dh.atoms[atom_idx].tag}({mag_ent['l']})" + ) + mag_ent["xyz"].append(dh.xyz[atom_idx]) + + # calculate size for Greens function generation + spin_box_shape = len(mag_ent["spin_box_indeces"]) + + mag_ent["energies"] = ( + [] + ) # we will store the second order energy derivations here + + # These will be the perturbed potentials from eq. 100 + mag_ent["Vu1"] = [] # so they are independent in memory + mag_ent["Vu2"] = [] + + mag_ent["Gii"] = [] # Greens function + mag_ent["Gii_tmp"] = [] # Greens function for parallelization + for i in simulation_parameters["ref_xcf_orientations"]: + # Rotations for every quantization axis + mag_ent["Vu1"].append([]) + mag_ent["Vu2"].append([]) + # Greens functions for every quantization axis + mag_ent["Gii"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape, spin_box_shape), + dtype="complex128", + ) + ) + mag_ent["Gii_tmp"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape, spin_box_shape), + dtype="complex128", + ) + ) + + # for every site we have to store 2x3 Greens function (and the associated _tmp-s) + # in the 3 reference directions, because G_ij and G_ji are both needed + for pair in pairs: + # calculate distance + xyz_ai = magnetic_entities[pair["ai"]]["xyz"] + xyz_aj = magnetic_entities[pair["aj"]]["xyz"] + xyz_aj = xyz_aj + pair["Ruc"] @ simulation_parameters["cell"] + pair["dist"] = np.linalg.norm(xyz_ai - xyz_aj) + + # calculate size for Greens function generation + spin_box_shape_i = len(magnetic_entities[pair["ai"]]["spin_box_indeces"]) + spin_box_shape_j = len(magnetic_entities[pair["aj"]]["spin_box_indeces"]) + pair["tags"] = [] + for mag_ent in [magnetic_entities[pair["ai"]], magnetic_entities[pair["aj"]]]: + tag = "" + # get atoms of magnetic entity + atoms_idx = mag_ent["atom"] + orbitals = mag_ent["l"] + + # if magnetic entity contains one atoms + if isinstance(atoms_idx, int): + tag += f"[{atoms_idx}]{dh.atoms[atoms_idx].tag}({orbitals})" + + # if magnetic entity contains more than one atoms + if isinstance(atoms_idx, list): + # iterate over atoms + atom_group = "{" + for atom_idx in atoms_idx: + atom_group += f"[{atom_idx}]{dh.atoms[atom_idx].tag}({orbitals})--" + # end {} of the atoms in the magnetic entity + tag += atom_group[:-2] + "}" + pair["tags"].append(tag) + pair["energies"] = [] # we will store the second order energy derivations here + + pair["Gij"] = [] # Greens function + pair["Gji"] = [] + pair["Gij_tmp"] = [] # Greens function for parallelization + pair["Gji_tmp"] = [] + for i in simulation_parameters["ref_xcf_orientations"]: + # Greens functions for every quantization axis + pair["Gij"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape_i, spin_box_shape_j), + dtype="complex128", + ) + ) + pair["Gij_tmp"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape_i, spin_box_shape_j), + dtype="complex128", + ) + ) + pair["Gji"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape_j, spin_box_shape_i), + dtype="complex128", + ) + ) + pair["Gji_tmp"].append( + np.zeros( + (simulation_parameters["eset"], spin_box_shape_j, spin_box_shape_i), + dtype="complex128", + ) + ) + + return pairs, magnetic_entities diff --git a/src/grogu_magn/grogu.py b/src/grogu_magn/grogu.py index e69de29..93d0ca4 100644 --- a/src/grogu_magn/grogu.py +++ b/src/grogu_magn/grogu.py @@ -0,0 +1,19 @@ +# 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. diff --git a/src/grogu_magn/io.py b/src/grogu_magn/io.py new file mode 100644 index 0000000..f9a5b47 --- /dev/null +++ b/src/grogu_magn/io.py @@ -0,0 +1,219 @@ +from argparse import ArgumentParser +from pickle import dump, load + +import numpy as np + +default_args = dict( + infile=None, + outfile=None, + scf_xcf_orientation=np.array([0, 0, 1]), + ref_xcf_orientations=[ + dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]), + dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]), + dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]), + ], + kset=2, + kdirs="xyz", + ebot=None, + eset=42, + esetp=1000, +) + +# parser = ArgumentParser() +# parser.add_argument('--input' , dest = 'infile' , default=None , help = 'Input file name') +# parser.add_argument('--output' , dest = 'outfile', default=None , help = 'Output file name') + +# parser.add_argument('--kset' , dest = 'kset' , default = 2 , type=int , help = 'k-space resolution of Jij calculation') +# parser.add_argument('--kdirs' , dest = 'kdirs' , default = 'xyz' , help = 'Definition of k-space dimensionality') +# parser.add_argument('--ebot' , dest = 'ebot' , default = None , type=float, help = 'Bottom energy of the contour') +# parser.add_argument('--eset' , dest = 'eset' , default = 42 , type=int , help = 'Number of energy points on the contour') +# parser.add_argument('--eset-p' , dest = 'esetp' , default = 1000 , type=int , help = 'Parameter tuning the distribution on the contour') +# cmd_line_args = parser.parse_args() + + +def save_pickle(outfile, data): + """_summary_ + + Args: + outfile (_type_): _description_ + data (_type_): _description_ + """ + # save dictionary + with open(outfile, "wb") as output_file: + dump(data, output_file) + + +def load_pickle(infile, data): + """_summary_ + + Args: + infile (_type_): _description_ + data (_type_): _description_ + + Returns: + _type_: _description_ + """ + with open(infile, "wb") as input_file: + data = load(data, input_file) + + return data + + +def print_parameters(simulation_parameters): + """_summary_ + + Args: + simulation_parameters (_type_): _description_ + """ + print( + "================================================================================================================================================================" + ) + print("Input file: ") + print(simulation_parameters["infile"]) + print("Output file: ") + print(simulation_parameters["outfile"]) + print( + "Number of nodes in the parallel cluster: ", + simulation_parameters["parallel_size"], + ) + print( + "================================================================================================================================================================" + ) + print("Cell [Ang]: ") + print(simulation_parameters["cell"]) + print( + "================================================================================================================================================================" + ) + print("DFT axis: ") + print(simulation_parameters["scf_xcf_orientation"]) + print("Quantization axis and perpendicular rotation directions:") + for ref in simulation_parameters["ref_xcf_orientations"]: + print(ref["o"], " --» ", ref["vw"]) + print( + "================================================================================================================================================================" + ) + print("Parameters for the contour integral:") + print("Number of k points: ", simulation_parameters["kset"]) + print("k point directions: ", simulation_parameters["kdirs"]) + print("Ebot: ", simulation_parameters["ebot"]) + print("Eset: ", simulation_parameters["eset"]) + print("Esetp: ", simulation_parameters["esetp"]) + print( + "================================================================================================================================================================" + ) + + +def print_atoms_and_pairs(magnetic_entities, pairs): + """_summary_ + + Args: + magnetic_entities (_type_): _description_ + pairs (_type_): _description_ + """ + print("Atomic information: ") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print( + "[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz" + ) + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + # iterate over magnetic entities + for mag_ent in magnetic_entities: + # iterate over atoms + for tag, xyz in zip(mag_ent["tags"], mag_ent["xyz"]): + # coordinates and tag + print(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}") + print("") + + print( + "================================================================================================================================================================" + ) + print("Anisotropy [meV]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print("Magnetic entity x [Ang] y [Ang] z [Ang]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + # iterate over magnetic entities + for mag_ent in magnetic_entities: + # iterate over atoms + for tag, xyz in zip(mag_ent["tags"], mag_ent["xyz"]): + # coordinates and tag + print(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}") + print("Consistency check: ", mag_ent["K_consistency"]) + print("Anisotropy diag: ", mag_ent["K"]) + print("") + + print( + "================================================================================================================================================================" + ) + print("Exchange [meV]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print("Magnetic entity1 Magnetic entity2 [i j k] d [Ang]") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + # iterate over pairs + for pair in pairs: + # print pair parameters + print( + f"{pair['tags'][0]} {pair['tags'][1]} {pair['Ruc']} d [Ang] {pair['dist']}" + ) + # print magnetic parameters + print("Isotropic: ", pair["J_iso"]) + print("DMI: ", pair["D"]) + print("Symmetric-anisotropy: ", pair["J_S"]) + print("J: ", pair["J"].flatten()) + print("Energies for debugging: ") + print(np.array(pair["energies"])) + print( + "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)" + ) + o1, o2, o3 = pair["energies"] + print(np.array([o2[-1], o3[0], o1[0]])) + print("Test J_xx = E(y,z) = E(z,y)") + print(o2[-1], o3[-1]) + print("") + + print( + "================================================================================================================================================================" + ) + + +def print_runtime_information(times): + """_summary_ + + Args: + times (_type_): _description_ + """ + print("Runtime information: ") + print(f"Total runtime: {times['end_time'] - times['start_time']} s") + print( + "----------------------------------------------------------------------------------------------------------------------------------------------------------------" + ) + print(f"Initial setup: {times['setup_time'] - times['start_time']} s") + print( + f"Hamiltonian conversion and XC field extraction: {times['H_and_XCF_time'] - times['setup_time']:.3f} s" + ) + print( + f"Pair and site datastructure creatrions: {times['site_and_pair_dictionaries_time'] - times['H_and_XCF_time']:.3f} s" + ) + print( + f"k set cration and distribution: {times['k_set_time'] - times['site_and_pair_dictionaries_time']:.3f} s" + ) + print( + f"Rotating XC potential: {times['reference_rotations_time'] - times['k_set_time']:.3f} s" + ) + print( + f"Greens function inversion: {times['green_function_inversion_time'] - times['reference_rotations_time']:.3f} s" + ) + print( + f"Calculate energies and magnetic components: {times['end_time'] - times['green_function_inversion_time']:.3f} s" + ) diff --git a/src/grogu_magn/magnetism.py b/src/grogu_magn/magnetism.py new file mode 100644 index 0000000..38a3e14 --- /dev/null +++ b/src/grogu_magn/magnetism.py @@ -0,0 +1,88 @@ +import numpy as np + + +def calculate_anisotropy_tensor(mag_ent): + """_summary_ + + Args: + mag_ent (_type_): _description_ + + Returns: + _type_: _description_ + """ + + energies = mag_ent["energies"] + Kxx = energies[1, 1] - energies[1, 0] + Kyy = energies[0, 1] - energies[0, 0] + Kzz = 0 + + calculated_diff = Kyy - Kxx + expected_diff = energies[2, 0] - energies[2, 1] + consistency_check = abs(calculated_diff - expected_diff) + + return Kxx, Kyy, Kzz, consistency_check + + +def calculate_exchange_tensor(pair): + """_summary_ + + Args: + pair (_type_): _description_ + + Returns: + _type_: _description_ + """ + energies = pair["energies"] + # Initialize output arrays + J = np.zeros((3, 3)) + D = np.zeros(3) + + # J matrix calculations + # J(1,1) = mean([DEij(2,2,2), DEij(2,2,3)]) + J[0, 0] = np.mean([energies[1, 3], energies[2, 3]]) + + # J(1,2) = -mean([DEij(1,2,3), DEij(2,1,3)]) + J[0, 1] = -np.mean([energies[2, 1], energies[2, 2]]) + J[1, 0] = J[0, 1] + + # J(1,3) = -mean([DEij(1,2,2), DEij(2,1,2)]) + J[0, 2] = -np.mean([energies[1, 1], energies[1, 2]]) + J[2, 0] = J[0, 2] + + # J(2,2) = mean([DEij(2,2,1), DEij(1,1,3)]) + J[1, 1] = np.mean([energies[0, 3], energies[2, 0]]) + + # J(2,3) = -mean([DEij(1,2,1), DEij(2,1,1)]) + J[1, 2] = -np.mean([energies[0, 1], energies[0, 2]]) + J[2, 1] = J[1, 2] + + # J(3,3) = mean([DEij(1,1,1), DEij(1,1,2)]) + J[2, 2] = np.mean([energies[0, 0], energies[1, 0]]) + + # D vector calculations + # D(1) = mean([DEij(1,2,1), -DEij(2,1,1)]) + D[0] = np.mean([energies[0, 1], -energies[0, 2]]) + + # D(2) = mean([DEij(2,1,2), -DEij(1,2,2)]) + D[1] = np.mean([energies[1, 2], -energies[1, 1]]) + + # D(3) = mean([DEij(1,2,3), -DEij(2,1,3)]) + D[2] = np.mean([energies[2, 1], -energies[2, 2]]) + + J_iso = np.trace(J) / 3 + J_S = (J - J_iso * np.eye(3)).flatten() + + return J_iso, J_S, D, J + + +def int_de_ke(traced, we): + """_summary_ + + Args: + traced (_type_): _description_ + we (_type_): _description_ + + Returns: + _type_: _description_ + """ + return np.trapz(-1 / np.pi * np.imag(traced * we)) diff --git a/src/grogu_magn/utils.py b/src/grogu_magn/utils.py index c719ef3..dd2b13e 100644 --- a/src/grogu_magn/utils.py +++ b/src/grogu_magn/utils.py @@ -18,11 +18,9 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE # SOFTWARE. -from itertools import permutations, product -from pprint import pprint - import numpy as np from scipy.special import roots_legendre +from sisl.io.siesta import eigSileSiesta # Pauli matrices tau_x = np.array([[0, 1], [1, 0]]) @@ -306,235 +304,15 @@ def blow_up_orbindx(orb_indices): return np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten() -def calculate_anisotropy_tensor(mag_ent): - """_summary_ - - Args: - mag_ent (_type_): _description_ - - Returns: - _type_: _description_ - """ - - energies = mag_ent["energies"] - Kxx = energies[1, 1] - energies[1, 0] - Kyy = energies[0, 1] - energies[0, 0] - Kzz = 0 - - calculated_diff = Kyy - Kxx - expected_diff = energies[2, 0] - energies[2, 1] - consistency_check = abs(calculated_diff - expected_diff) - - return Kxx, Kyy, Kzz, consistency_check - - -def calculate_exchange_tensor(pair): +def read_siesta_emin(eigfile): """_summary_ Args: - pair (_type_): _description_ + eigfile (_type_): _description_ Returns: _type_: _description_ """ - energies = pair["energies"] - # Initialize output arrays - J = np.zeros((3, 3)) - D = np.zeros(3) - - # J matrix calculations - # J(1,1) = mean([DEij(2,2,2), DEij(2,2,3)]) - J[0, 0] = np.mean([energies[1, 3], energies[2, 3]]) - - # J(1,2) = -mean([DEij(1,2,3), DEij(2,1,3)]) - J[0, 1] = -np.mean([energies[2, 1], energies[2, 2]]) - J[1, 0] = J[0, 1] - - # J(1,3) = -mean([DEij(1,2,2), DEij(2,1,2)]) - J[0, 2] = -np.mean([energies[1, 1], energies[1, 2]]) - J[2, 0] = J[0, 2] - - # J(2,2) = mean([DEij(2,2,1), DEij(1,1,3)]) - J[1, 1] = np.mean([energies[0, 3], energies[2, 0]]) - - # J(2,3) = -mean([DEij(1,2,1), DEij(2,1,1)]) - J[1, 2] = -np.mean([energies[0, 1], energies[0, 2]]) - J[2, 1] = J[1, 2] - - # J(3,3) = mean([DEij(1,1,1), DEij(1,1,2)]) - J[2, 2] = np.mean([energies[0, 0], energies[1, 0]]) - - # D vector calculations - # D(1) = mean([DEij(1,2,1), -DEij(2,1,1)]) - D[0] = np.mean([energies[0, 1], -energies[0, 2]]) - - # D(2) = mean([DEij(2,1,2), -DEij(1,2,2)]) - D[1] = np.mean([energies[1, 2], -energies[1, 1]]) - - # D(3) = mean([DEij(1,2,3), -DEij(2,1,3)]) - D[2] = np.mean([energies[2, 1], -energies[2, 2]]) - - J_iso = np.trace(J) / 3 - J_S = (J - J_iso * np.eye(3)).flatten() - - return J_iso, J_S, D, J - - -def print_parameters(simulation_parameters): - """_summary_ - - Args: - simulation_parameters (_type_): _description_ - """ - print( - "================================================================================================================================================================" - ) - print("Input file: ") - print(simulation_parameters["path"]) - print("Output file: ") - print(simulation_parameters["outpath"]) - print( - "Number of nodes in the parallel cluster: ", - simulation_parameters["parallel_size"], - ) - print( - "================================================================================================================================================================" - ) - print("Cell [Ang]: ") - print(simulation_parameters["cell"]) - print( - "================================================================================================================================================================" - ) - print("DFT axis: ") - print(simulation_parameters["scf_xcf_orientation"]) - print("Quantization axis and perpendicular rotation directions:") - for ref in simulation_parameters["ref_xcf_orientations"]: - print(ref["o"], " --» ", ref["vw"]) - print( - "================================================================================================================================================================" - ) - print("Parameters for the contour integral:") - print("Number of k points: ", simulation_parameters["kset"]) - print("k point directions: ", simulation_parameters["kdirs"]) - print("Ebot: ", simulation_parameters["ebot"]) - print("Eset: ", simulation_parameters["eset"]) - print("Esetp: ", simulation_parameters["esetp"]) - print( - "================================================================================================================================================================" - ) - - -def print_atoms_and_pairs(magnetic_entities, pairs): - """_summary_ - - Args: - magnetic_entities (_type_): _description_ - pairs (_type_): _description_ - """ - print("Atomic information: ") - print( - "----------------------------------------------------------------------------------------------------------------------------------------------------------------" - ) - print( - "[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz" - ) - print( - "----------------------------------------------------------------------------------------------------------------------------------------------------------------" - ) - # iterate over magnetic entities - for mag_ent in magnetic_entities: - # iterate over atoms - for tag, xyz in zip(mag_ent["tags"], mag_ent["xyz"]): - # coordinates and tag - print(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}") - print("") - - print( - "================================================================================================================================================================" - ) - print("Anisotropy [meV]") - print( - "----------------------------------------------------------------------------------------------------------------------------------------------------------------" - ) - print("Magnetic entity x [Ang] y [Ang] z [Ang]") - print( - "----------------------------------------------------------------------------------------------------------------------------------------------------------------" - ) - # iterate over magnetic entities - for mag_ent in magnetic_entities: - # iterate over atoms - for tag, xyz in zip(mag_ent["tags"], mag_ent["xyz"]): - # coordinates and tag - print(f"{tag} {xyz[0]} {xyz[1]} {xyz[2]}") - print("Consistency check: ", mag_ent["K_consistency"]) - print("Anisotropy diag: ", mag_ent["K"]) - print("") - - print( - "================================================================================================================================================================" - ) - print("Exchange [meV]") - print( - "----------------------------------------------------------------------------------------------------------------------------------------------------------------" - ) - print("Magnetic entity1 Magnetic entity2 [i j k] d [Ang]") - print( - "----------------------------------------------------------------------------------------------------------------------------------------------------------------" - ) - # iterate over pairs - for pair in pairs: - # print pair parameters - print( - f"{pair['tags'][0]} {pair['tags'][1]} {pair['Ruc']} d [Ang] {pair['dist']}" - ) - # print magnetic parameters - print("Isotropic: ", pair["J_iso"]) - print("DMI: ", pair["D"]) - print("Symmetric-anisotropy: ", pair["J_S"]) - print("J: ", pair["J"].flatten()) - print("Energies for debugging: ") - pprint(np.array(pair["energies"])) - print( - "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)" - ) - o1, o2, o3 = pair["energies"] - pprint(np.array([o2[-1], o3[0], o1[0]])) - print("Test J_xx = E(y,z) = E(z,y)") - print(o2[-1], o3[-1]) - print("") - - print( - "================================================================================================================================================================" - ) + eigs = eigSileSiesta(eigfile).read_data() - -def print_runtime_information(times): - """_summary_ - - Args: - times (_type_): _description_ - """ - print("Runtime information: ") - print(f"Total runtime: {times['end_time'] - times['start_time']} s") - print( - "----------------------------------------------------------------------------------------------------------------------------------------------------------------" - ) - print(f"Initial setup: {times['setup_time'] - times['start_time']} s") - print( - f"Hamiltonian conversion and XC field extraction: {times['H_and_XCF_time'] - times['setup_time']:.3f} s" - ) - print( - f"Pair and site datastructure creatrions: {times['site_and_pair_dictionaries_time'] - times['H_and_XCF_time']:.3f} s" - ) - print( - f"k set cration and distribution: {times['k_set_time'] - times['site_and_pair_dictionaries_time']:.3f} s" - ) - print( - f"Rotating XC potential: {times['reference_rotations_time'] - times['k_set_time']:.3f} s" - ) - print( - f"Greens function inversion: {times['green_function_inversion_time'] - times['reference_rotations_time']:.3f} s" - ) - print( - f"Calculate energies and magnetic components: {times['end_time'] - times['green_function_inversion_time']:.3f} s" - ) + return eigs.min()