# Copyright (c) [2024] [] # # 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 warnings from argparse import ArgumentParser from sys import getsizeof from timeit import default_timer as timer # use numpy number of threads one try: from threadpoolctl import threadpool_info, threadpool_limits user_api = threadpool_info()["user_api"] threadpool_limits(limits=1, user_api=user_api) except: print("Warning: threadpoolctl could not make numpy use single thread!") from pickle import dump, load import numpy as np import sisl from mpi4py import MPI from numpy.linalg import inv from scipy.special import roots_legendre try: from tqdm import tqdm tqdm_imported = True except: print("Please install tqdm for nice progress bar.") tqdm_imported = False # Pauli matrices TAU_X = np.array([[0, 1], [1, 0]]) TAU_Y = np.array([[0, -1j], [1j, 0]]) TAU_Z = np.array([[1, 0], [0, -1]]) TAU_0 = np.array([[1, 0], [0, 1]]) # list of accepted input parameters ACCEPTED_INPUTS = [ "infile", "outfile", "scf_xcf_orientation", "ref_xcf_orientations", "kset", "kdirs", "ebot", "eset", "esetp", "parallel_solver_for_Gk", "padawan_mode", ] DEFAULT_ARGUMENTS = 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, parallel_solver_for_Gk=False, padawan_mode=True, ) def commutator(a, b): """Shorthand for commutator. Commutator of two matrices in the mathematical sense. Args: a : np.array_like The first matrix b : np.array_like The second matrix Returns: np.array_like The commutator of a and b """ return a @ b - b @ a # define some useful functions def hsk(H, ss, sc_off, k=(0, 0, 0)): """Speed up Hk and Sk generation. Calculates the Hamiltonian and the Overlap matrix at a given k point. It is faster that the sisl version. Args: H : np.array_like Hamiltonian in spin box form ss : np.array_like Overlap matrix in spin box form sc_off : list supercell indexes of the Hamiltonian k : tuple, optional The k point where the matrices are set up. Defaults to (0, 0, 0) Returns: np.array_like Hamiltonian at the given k point np.array_like Overlap matrix at the given k point """ # this two conversion lines are from the sisl source k = np.asarray(k, np.float64) k.shape = (-1,) # this generates the list of phases phases = np.exp(-1j * 2 * np.pi * k @ sc_off.T) # phases applied to the hamiltonian HK = np.einsum("abc,a->bc", H, phases) SK = np.einsum("abc,a->bc", ss, phases) return HK, SK def make_kset(dirs="xyz", NUMK=20): """Simple k-grid generator to sample the Brillouin zone. Depending on the value of the dirs argument k sampling in 1,2 or 3 dimensions is generated. If dirs argument does not contain either of x, y or z a kset of a single k-pont at the origin is returned. The total number of k points is the NUMK**(dimensions) Args: dirs : str, optional Directions of the k points in the Brillouin zone. They are the three lattice vectors. Defaults to "xyz" NUMK : int, optional The number of k points in a direction. Defaults to 20 Returns: np.array_like An array of k points that uniformly sample the Brillouin zone in the given directions """ # if there is no xyz in dirs return the Gamma point if not (sum([d in dirs for d in "xyz"])): return np.array([[0, 0, 0]]) kran = len(dirs) * [np.linspace(0, 1, NUMK, endpoint=False)] mg = np.meshgrid(*kran) dirsdict = dict() for d in enumerate(dirs): dirsdict[d[1]] = mg[d[0]].flatten() for d in "xyz": if not (d in dirs): dirsdict[d] = 0 * dirsdict[dirs[0]] kset = np.array([dirsdict[d] for d in "xyz"]).T return kset def make_contour(emin=-20, emax=0.0, enum=42, p=150): """A more sophisticated contour generator. Calculates the parameters for the complex contour integral. It uses the Legendre-Gauss quadrature method. It returns a class that contains the information for the contour integral. Args: emin : int, optional Energy minimum of the contour. Defaults to -20 emax : float, optional Energy maximum of the contour. Defaults to 0.0, so the Fermi level enum : int, optional Number of sample points along the contour. Defaults to 42 p : int, optional Shape parameter that describes the distribution of the sample points. Defaults to 150 Returns: ccont Contains all the information for the contour integral """ x, wl = roots_legendre(enum) R = (emax - emin) / 2 # radius z0 = (emax + emin) / 2 # center point y1 = -np.log(1 + np.pi * p) # lower bound y2 = 0 # upper bound y = (y2 - y1) / 2 * x + (y2 + y1) / 2 phi = (np.exp(-y) - 1) / p # angle parameter ze = z0 + R * np.exp(1j * phi) # complex points for path we = -(y2 - y1) / 2 * np.exp(-y) / p * 1j * (ze - z0) * wl # just an empty container class class ccont: pass cont = ccont() cont.R = R cont.z0 = z0 cont.ze = ze cont.we = we cont.enum = enum return cont def tau_u(u): """Pauli matrix in direction u. Returns the vector u in the basis of the Pauli matrices. Args: u : list or np.array_like The direction Returns: np.array_like Arbitrary direction in the base of the Pauli matrices """ # u is force to be of unit length u = u / np.linalg.norm(u) return u[0] * TAU_X + u[1] * TAU_Y + u[2] * TAU_Z # def crossM(u): """Definition for the cross-product matrix. It acts as a cross product with vector u. Args: u : list or np.array_like The second vector in the cross product Returns: np.array_like The matrix that represents teh cross product with a vector """ return np.array([[0, -u[2], u[1]], [u[2], 0, -u[0]], [-u[1], u[0], 0]]) def RotM(theta, u, eps=1e-10): """Definition of rotation matrix with angle theta around direction u. Args: theta : float The angle of rotation u : np.array_like The rotation axis eps : float, optional Cutoff for small elements in the resulting matrix. Defaults to 1e-10 Returns: np.array_like The rotation matrix """ u = u / np.linalg.norm(u) M = ( np.cos(theta) * np.eye(3) + np.sin(theta) * crossM(u) + (1 - np.cos(theta)) * np.outer(u, u) ) # kill off small numbers M[abs(M) < eps] = 0.0 return M def RotMa2b(a, b, eps=1e-10): """Definition of rotation matrix rotating unit vector a to unit vector b. Function returns array R such that R @ a = b holds. Args: a : np.array_like First vector b : np.array_like Second vector eps : float, optional Cutoff for small elements in the resulting matrix. Defaults to 1e-10 Returns: np.array_like The rotation matrix with the above property """ v = np.cross(a, b) c = a @ b M = np.eye(3) + crossM(v) + crossM(v) @ crossM(v) / (1 + c) # kill off small numbers M[abs(M) < eps] = 0.0 return M def read_siesta_emin(eigfile): """It reads the lowest energy level from the siesta run. It uses the .EIG file from siesta that contains the eigenvalues. Args: eigfile : str The path to the .EIG file Returns: float The energy minimum """ # read the file eigs = eigSileSiesta(eigfile).read_data() return eigs.min() def int_de_ke(traced, we): """It numerically integrates the traced matrix. It is a wrapper from numpy.trapz and it contains the relevant constants to calculate the energy integral from equation 93 or 96. Args: traced : np.array_like The trace of a matrix or a matrix product we : float The weight of a point on the contour Returns: float The energy calculated from the integral formula """ return np.trapz(-1 / np.pi * np.imag(traced * we)) def onsite_projection(matrix, idx1, idx2): """It produces the slices of a matrix for the on site projection. The slicing is along the last two axes as these contains the orbital indexing. Args: matrix : (..., :, :) np.array_like Some matrix idx : np.array_like The indexes of the orbitals Returns: np.array_like Reduced matrix based on the projection """ return matrix[..., idx1, :][..., idx2] def calc_Vu(H, Tu): """Calculates the local perturbation in case of a spin rotation. Args: H : (NO, NO) np.array_like Hamiltonian Tu : (NO, NO) array_like Rotation around u Returns: Vu1 : (NO, NO) np.array_like First order perturbed matrix Vu2 : (NO, NO) np.array_like Second order perturbed matrix """ Vu1 = 1j / 2 * commutator(H, Tu) # equation 100 Vu2 = 1 / 8 * commutator(commutator(Tu, H), Tu) # equation 100 return Vu1, Vu2 def build_hh_ss(dh): """It builds the Hamiltonian and Overlap matrix from the sisl.dh class. It restructures the data in the SPIN BOX representation, where NS is the number of supercells and NO is the number of orbitals. Args: dh : sisl.physics.Hamiltonian Hamiltonian read in by sisl Returns: hh : (NS, NO, NO) np.array_like Hamiltonian in SPIN BOX representation ss : (NS, NO, NO) np.array_like Overlap matrix in SPIN BOX representation """ 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 def setup_pairs_and_magnetic_entities( magnetic_entities, pairs, dh, simulation_parameters ): """It creates the complete structure of the dictionaries and fills some basic data. It creates orbital indexes, spin box indexes, coordinates and tags for magnetic entities. Furthermore it creates the structures for the energies, the perturbed potentials and the Greens function calculation. It dose the same for the pairs. Args: pairs : dict Contains the initial pair information magnetic_entities : dict Contains the initial magnetic entity information dh : sisl.physics.Hamiltonian Hamiltonian read in by sisl simulation_parameters : dict A set of parameters from the simulation Returns: pairs : dict Contains the initial information and the complete structure magnetic_entities : dict Contains the initial information and the complete structure """ # 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_indices"] = parsed mag_ent["spin_box_indices"] = 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" # tag creation for one atom 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"]]] # tag creation for more atoms 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_indices"]) # we will store the second order energy derivations here mag_ent["energies"] = [] # 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 _ 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_indices"]) spin_box_shape_j = len(magnetic_entities[pair["aj"]]["spin_box_indices"]) # tag generation 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 _ 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 def parallel_Gk(HK, SK, eran, eset): """Calculates the Greens function by inversion. It calculates the Greens function on all the energy levels at the same time. Args: HK : (NO, NO), np.array_like Hamiltonian at a given k point SK : (NO, NO), np.array_like Overlap Matrix at a given k point eran : (eset) np.array_like Energy sample along the contour eset : int Number of energy samples along the contour Returns: Gk : (eset, NO, NO), np.array_like Green's function at a given k point """ # Calculates the Greens function on all the energy levels return inv(SK * eran.reshape(eset, 1, 1) - HK) def sequential_GK(HK, SK, eran, eset): """Calculates the Greens function by inversion. It calculates sequentially over the energy levels. Args: HK : (NO, NO), np.array_like Hamiltonian at a given k point SK : (NO, NO), np.array_like Overlap Matrix at a given k point eran : (eset) np.array_like Energy sample along the contour eset : int Number of energy samples along the contour Returns: Gk : (eset, NO, NO), np.array_like Green's function at a given k point """ # creates an empty holder Gk = np.zeros(shape=(eset, HK.shape[0], HK.shape[1]), dtype="complex128") # fills the holder sequentially by the Greens function on a given energy for j in range(eset): Gk[j] = inv(SK * eran[j] - HK) return Gk def remove_clutter_for_save(pairs, magnetic_entities): """Removes unimportant data from the dictionaries. It is used before saving to throw away data that is not needed for post processing. Args: pairs : dict Contains all the pair information magnetic_entities : dict Contains all the magnetic entity information Returns: pairs : dict Contains all the reduced pair information magnetic_entities : dict Contains all the reduced magnetic entity information """ # 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 blow_up_orbindx(orb_indices): """Function to blow up orbital indices to make SPIN BOX indices. Args: orb_indices : np.array_like These are the indices in ORBITAL BOX Returns: orb_indices : np.array_like These are the indices in SPIN BOX """ orb_indices = np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten() return orb_indices def spin_tracer(M): """Spin tracer utility. This takes an operator with the orbital-spin sequence: orbital 1 up, orbital 1 down, orbital 2 up, orbital 2 down, that is in the SPIN-BOX representation, and extracts orbital dependent Pauli traces. Args: M : np.array_like Traceable matrix Returns: dict It contains the traced matrix with "x", "y", "z" and "c" """ M11 = M[0::2, 0::2] M12 = M[0::2, 1::2] M21 = M[1::2, 0::2] M22 = M[1::2, 1::2] M_o = dict() M_o["x"] = M12 + M21 M_o["y"] = 1j * (M12 - M21) M_o["z"] = M11 - M22 M_o["c"] = M11 + M22 return M_o def parse_magnetic_entity(dh, atom=None, l=None, **kwargs): """Function to define orbital indexes of a given magnetic entity. Args: dh : sisl.physics.Hamiltonian Hamiltonian from sisl atom : integer or list of integers, optional Defining atom (or atoms) in the unit cell forming the magnetic entity. Defaults to None l : integer, optional Defining the angular momentum channel. Defaults to None Returns: list The orbital indexes of the given magnetic entity """ # case where we deal with more than one atom defining the magnetic entity if type(atom) == list: dat = [] for a in atom: a_orb_idx = dh.geometry.a2o(a, all=True) if ( type(l) == int ): # if specified we restrict to given l angular momentum channel inside each atom a_orb_idx = a_orb_idx[[o.l == l for o in dh.geometry.atoms[a].orbitals]] dat.append(a_orb_idx) orbital_indeces = np.hstack(dat) # case where we deal with a singel atom magnetic entity elif type(atom) == int: orbital_indeces = dh.geometry.a2o(atom, all=True) if ( type(l) == int ): # if specified we restrict to given l angular momentum channel orbital_indeces = orbital_indeces[ [o.l == l for o in dh.geometry.atoms[atom].orbitals] ] return orbital_indeces # numpy array containing integers labeling orbitals associated to a magnetic entity. def calculate_anisotropy_tensor(mag_ent): """Calculates the renormalized anisotropy tensor from the energies. Args: mag_ent : dict An element from the magnetic entities Returns: K : np.array_like elements of the anisotropy tensor consistency_check : float Absolute value of the difference from the consistency check """ # get the energies energies = mag_ent["energies"] K = np.zeros((3, 3)) # calculate the diagonal tensor elements K[0, 0] = energies[1, 1] - energies[1, 0] K[1, 1] = energies[0, 1] - energies[0, 0] K[2, 2] = 0 # calculate the off-diagonal tensor elements K[0, 1] = (energies[2, 0] + energies[2, 1]) / 2 - energies[2, 2] K[1, 0] = K[0, 1] K[0, 2] = (energies[1, 0] + energies[1, 1]) / 2 - energies[1, 2] K[2, 0] = K[0, 2] K[1, 2] = (energies[0, 0] + energies[0, 1]) / 2 - energies[0, 2] K[2, 1] = K[1, 2] # perform consistency check calculated_diff = K[1, 1] - K[0, 0] expected_diff = energies[2, 0] - energies[2, 1] consistency_check = abs(calculated_diff - expected_diff) return K, consistency_check def calculate_exchange_tensor(pair): """Calculates the exchange tensor from the energies. It produces the isotropic exchange, the relevant elements from the Dzyaloshinskii-Morilla (Dm) tensor, the symmetric-anisotropy and the complete exchange tensor. Args: pair : dict An element from the pairs Returns: J_iso : float Isotropic exchange (Tr[J] / 3) J_S : np.array_like Symmetric-anisotropy (J_S = J - J_iso * I ––> Jxx, Jyy, Jxy, Jxz, Jyz) D : np.array_like DM elements (Dx, Dy, Dz) J : np.array_like Complete exchange tensor flattened (Jxx, Jxy, Jxz, Jyx, Jyy, Jyz, Jzx, Jzy, Jzz) """ # energies from rotations 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 # based on the grogu output pdf # traceless symmetric exchange matrix: # Jxx, Jyy, Jxy, Jxz, Jyz J_S = np.array([J[0, 0] - J_iso, J[1, 1] - J_iso, J[0, 1], J[0, 2], J[1, 2]]) return J_iso, J_S, D, J def read_fdf(path): """It reads the simulation parameters, magnetic entities and pairs from the fdf. Args: path : string The path to the .fdf file Returns: fdf_arguments : dict The read input arguments from the fdf file magnetic_entities : list It contains the dictionaries associated with the magnetic entities pairs : dict It contains the dictionaries associated with the pair information """ # read fdf file fdf = fdfSileSiesta(path) fdf_arguments = dict() InputFile = fdf.get("InputFile") if InputFile is not None: fdf_arguments["infile"] = InputFile OutputFile = fdf.get("OutputFile") if OutputFile is not None: fdf_arguments["outfile"] = OutputFile ScfXcfOrientation = fdf.get("ScfXcfOrientation") if ScfXcfOrientation is not None: fdf_arguments["scf_xcf_orientation"] = np.array(ScfXcfOrientation) XCF_Rotation = fdf.get("XCF_Rotation") if XCF_Rotation is not None: rotations = [] # iterate over rows for rot in XCF_Rotation: # convert row to dictionary dat = np.array(rot.split()[:9], dtype=float) o = dat[:3] vw = dat[3:].reshape(2, 3) rotations.append(dict(o=o, vw=vw)) fdf_arguments["ref_xcf_orientations"] = rotations Kset = fdf.get("INTEGRAL.Kset") if Kset is not None: fdf_arguments["kset"] = Kset Kdirs = fdf.get("INTEGRAL.Kdirs") if Kdirs is not None: fdf_arguments["kdirs"] = Kdirs # This is permitted because it means automatic Ebot definition fdf_arguments["ebot"] = fdf.get("INTEGRAL.Ebot") Eset = fdf.get("INTEGRAL.Eset") if Eset is not None: fdf_arguments["eset"] = Eset Esetp = fdf.get("INTEGRAL.Esetp") if Esetp is not None: fdf_arguments["esetp"] = Esetp ParallelSolver = fdf.get("GREEN.ParallelSolver") if ParallelSolver is not None: fdf_arguments["parallel_solver_for_Gk"] = ParallelSolver PadawanMode = fdf.get("PadawanMode") if PadawanMode is not None: fdf_arguments["padawan_mode"] = PadawanMode Pairs = fdf.get("Pairs") if Pairs is not None: pairs = [] # iterate over rows for fdf_pair in Pairs: # convert data dat = np.array(fdf_pair.split()[:5], dtype=int) # create pair dictionary my_pair = dict(ai=dat[0], aj=dat[1], Ruc=np.array(dat[2:])) pairs.append(my_pair) MagneticEntities = fdf.get("MagneticEntities") if MagneticEntities is not None: magnetic_entities = [] for mag_ent in MagneticEntities: row = mag_ent.split() dat = [] for string in row: if string.find("#") != -1: break dat.append(string) if dat[0] in {"Cluster", "cluster"}: magnetic_entities.append(dict(atom=[int(_) for _ in dat[1:]])) continue elif dat[0] in {"AtomShell", "Atomshell", "atomShell", "atomshell"}: magnetic_entities.append( dict(atom=int(dat[1]), l=[int(_) for _ in dat[2:]]) ) continue elif dat[0] in {"AtomOrbital", "Atomorbital", "tomOrbital", "atomorbital"}: magnetic_entities.append( dict(atom=int(dat[1]), orb=[int(_) for _ in dat[2:]]) ) continue elif dat[0] in {"Orbitals", "orbitals"}: magnetic_entities.append(dict(orb=[int(_) for _ in dat[1:]])) continue else: raise Exception("Unrecognizable magnetic entity in .fdf!") return fdf_arguments, magnetic_entities, pairs def process_input_args( default_arguments, fdf_arguments, command_line_arguments, accepted_inputs=ACCEPTED_INPUTS, ): """It returns the final simulation parameters based on the inputs. The merging is done in the order of priority: 1. command line arguments 2. fdf arguments 3. default arguments Args: default_arguments : dict Default arguments from grogupy fdf_arguments : dict Arguments read from the fdf input file command_line_arguments : dict Arguments from the command line Returns: dict The final simulation parameters """ # iterate over fdf_arguments and update default arguments for key, value in fdf_arguments.values(): if value is not None and key in accepted_inputs: default_arguments[key] = value # iterate over command_line_arguments and update default arguments for key, value in command_line_arguments.values(): if value is not None and key in accepted_inputs: default_arguments[key] = value return default_arguments def save_pickle(outfile, data): """Saves the data in the outfile with pickle. Args: outfile : str Path to outfile data : dict Contains the data """ # save dictionary with open(outfile, "wb") as output_file: dump(data, output_file) def load_pickle(infile): """Loads the data from the infile with pickle. Args: infile : str Path to infile Returns: data : dict A dictionary of data """ # open and read file with open(infile, "wb") as input_file: data = load(data, input_file) return data def print_job_description(simulation_parameters): """It prints the parameters and the description of the job. Args: simulation_parameters : dict It contains the simulations parameters """ 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"], ) if simulation_parameters["parallel_solver_for_Gk"]: print("solver used for Greens function calculation: parallel") else: print("solver used for Greens function calculation: sequential") 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"]) if simulation_parameters["automatic_ebot"]: print( "Ebot: ", simulation_parameters["ebot"], " WARNING: This was automatically determined!", ) else: print("Ebot: ", simulation_parameters["ebot"]) print("Eset: ", simulation_parameters["eset"]) print("Esetp: ", simulation_parameters["esetp"]) print( "================================================================================================================================================================" ) def print_parameters(simulation_parameters): """It prints the simulation parameters for the grogu out. Args: simulation_parameters : dict It contains the simulations parameters """ 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): """It prints the pair and magnetic entity information for the grogu out. Args: magnetic_entities : dict It contains the data on the magnetic entities pairs : dict It contains the data on the pairs """ 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("K: # Kxx, Kxy, Kxz, Kyx, Kyy, Kyz, Kzx, Kzy, Kzz") print(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"], " # Tr[J] / 3") print("") print("DMI: ", pair["D"], " # Dx, Dy, Dz") print("") print( "Symmetric-anisotropy: ", pair["J_S"], " # J_S = J - J_iso * I ––> Jxx, Jyy, Jxy, Jxz, Jyz", ) print("") print("J: # Jxx, Jxy, Jxz, Jyx, Jyy, Jyz, Jzx, Jzy, Jzz") print(pair["J"]) print( "----------------------------------------------------------------------------------------------------------------------------------------------------------------" ) print( "================================================================================================================================================================" ) def print_runtime_information(times): """It prints the runtime information for the grogu out. Args: times : dict It contains the runtime data """ 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" ) def parse_command_line(): """This function can read input from the command line.""" parser = ArgumentParser() parser.add_argument( "-i", "--input", dest="infile", default=None, type=str, help="Input file name", required=True, ) parser.add_argument( "-o", "--output", dest="outfile", default=None, type=str, help="Output file name", ) parser.add_argument( "--kset", dest="kset", default=None, type=int, help="k-space resolution of calculation", ) parser.add_argument( "--kdirs", dest="kdirs", default=None, type=str, 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=None, type=int, help="Number of energy points on the contour", ) parser.add_argument( "--eset-p", dest="esetp", default=None, type=int, help="Parameter tuning the distribution on the contour", ) parser.add_argument( "--parallel-green", dest="parallel_solver_for_Gk", default=None, type=bool, help="Whether to use the parallel or sequential solver for Greens function", ) parser.add_argument( "--padawan-mode", dest="padawan_mode", default=None, type=bool, help="If it is on it turns on extra helpful information for new users", ) cmd_line_args = parser.parse_args() return cmd_line_args def main(simulation_parameters, magnetic_entities, pairs): # runtime information times = dict() times["start_time"] = timer() # MPI parameters comm = MPI.COMM_WORLD size = comm.Get_size() rank = comm.Get_rank() root_node = 0 if rank == root_node: # include parallel size in simulation parameters simulation_parameters["parallel_size"] = size # check versions for debugging 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 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) - 1 simulation_parameters["automatic_ebot"] = True except: print("Could not determine ebot.") print("Parameter was not given and .EIG file was not found.") else: simulation_parameters["automatic_ebot"] = False # add third direction for off-diagonal anisotropy for orientation in simulation_parameters["ref_xcf_orientations"]: o1 = orientation["vw"][0] o2 = orientation["vw"][1] orientation["vw"].append((o1 + o2) / np.sqrt(2)) # 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( "================================================================================================================================================================" ) NO = dh.no # shorthand for number of orbitals in the unit cell # reformat Hamltonian and Overlap matrix for manipulations hh, ss = 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)) # 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"]) # 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 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: # 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: K, consistency = calculate_anisotropy_tensor(mag_ent) mag_ent["K"] = K * 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") if __name__ == "__main__": # loading parameters # it is not clear how to give grogu.fdf path... # command_line_arguments = parse_command_line() # fdf_arguments, magnetic_entities, pairs = read_fdf(command_line_arguments["infile"]) # right now we do not use any of these input, but it shows # the order of priority when processing arguments default_arguments = False fdf_arguments = False command_line_arguments = False # simulation_parameters = process_input_args( # default_arguments, fdf_arguments, command_line_arguments, ACCEPTED_INPUTS # ) #################################################################################################### # This is the input file for now # #################################################################################################### 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 = dict() simulation_parameters["infile"] = "" simulation_parameters["outfile"] = "./" simulation_parameters["scf_xcf_orientation"] = np.array([0, 0, 1]) simulation_parameters["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])]), ] simulation_parameters["kset"] = 0 simulation_parameters["kdirs"] = "" simulation_parameters["ebot"] = None simulation_parameters["eset"] = 0 simulation_parameters["esetp"] = 0 simulation_parameters["parallel_solver_for_Gk"] = False simulation_parameters["padawan_mode"] = True #################################################################################################### # This is the input file for now # #################################################################################################### main(simulation_parameters, magnetic_entities, pairs)