et209d
/
grogu
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

WARNING: ran pre commit hooks on these, this may break stuff...

Browse Source
class-solution
Daniel Pozsar 2 months ago
parent 49d2be4982
commit 0db938aa50
30 changed files with 29329 additions and 29929 deletions
Show Stats Download Patch File Download Diff File

BIN
Fe3GeTe2_final_test.pickle
View File

Binary file not shown.

BIN
Fe3GeTe2_notebook.pickle
View File

Binary file not shown.

1
data/lat3_791/Fe.psf
Unescape View File

@ -2829,4 +2829,3 @@
0.159623499658E-47 0.394685762244E-48 0.959152537476E-49 0.229037992853E-49 0.159623499658E-47 0.394685762244E-48 0.959152537476E-49 0.229037992853E-49
0.537289950000E-50 0.123783348782E-50 0.280097688343E-51 0.624177220781E-52 0.537289950000E-50 0.123783348782E-50 0.280097688343E-51 0.624177220781E-52
0.136516118863E-52 0.292975523658E-53 0.616797664740E-54 0.136516118863E-52 0.292975523658E-53 0.616797664740E-54

1
data/lat3_791/Fe3GeTe2.ORB_INDX
Unescape View File

@ -6823,4 +6823,3 @@ spec = Atomic species label
isc = Unit cell indexes to which orbital belongs: isc = Unit cell indexes to which orbital belongs:
center(io) = center(iuo) + sum_(i=1:3) cell_vec(i) * isc(i) center(io) = center(iuo) + sum_(i=1:3) cell_vec(i) * isc(i)
iuo = Equivalent orbital in first unit cell iuo = Equivalent orbital in first unit cell

1
data/lat3_791/Fe3GeTe2.bib
Unescape View File

@ -21,4 +21,3 @@
issue = {086005}, issue = {086005},
doi = {10.1088/0953-8984/24/8/086005}, doi = {10.1088/0953-8984/24/8/086005},
} }

1
data/lat3_791/Fe3GeTe2.times
Unescape View File

@ -72,4 +72,3 @@ Tot.tot: total CPU time in all programs in one node
Nod.avg: average calculation time in one program across nodes Nod.avg: average calculation time in one program across nodes
Nod.max: maximum calculation time in one program across nodes Nod.max: maximum calculation time in one program across nodes
Calculation time: CPU time excluding communications Calculation time: CPU time excluding communications

1
data/lat3_791/TIMES
Unescape View File

@ -72,4 +72,3 @@ Tot.tot: total CPU time in all programs in one node
Nod.avg: average calculation time in one program across nodes Nod.avg: average calculation time in one program across nodes
Nod.max: maximum calculation time in one program across nodes Nod.max: maximum calculation time in one program across nodes
Calculation time: CPU time excluding communications Calculation time: CPU time excluding communications

47
input.fdf
Unescape View File

@ -0,0 +1,47 @@
InputFile /Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf
OutputFile ./Fe3GeTe2_notebook
ScfOrientation [ 0 0 1 ]
%block XCF_Rotation
1 0 0
0 1 0
0 0 1
%endblock XCFRotation
%block XCF_Rotation
0 1 0
1 0 0
0 0 1
%endblock XCFRotation
%block XCF_Rotation
0 0 1
1 0 0
0 1 0
%endblock XCFRotation
%block MagneticEntites # atom index and orbital index
3 2
4 2
5 2
%endblock MagneticEntites
%Pairsblock # MagneticEntites index ai and aj, supercell offset
0 1 0 0 0
0 2 0 0 0
1 2 0 0 0
0 2 -1 -1 0
1 2 -1 -1 0
0 2 -1 0 0
1 2 -1 0 0
1 2 -2 0 0
1 2 -3 0 0
%endPairsblock
INTEGRAL.Kset 3
INTEGRAL.Kdirs xy
INTEGRAL.Ebot -13
INTEGRAL.Eset 300
INTEGRAL.Esetp 1000

799
test.ipynb
View File

File diff suppressed because it is too large Load Diff

533
test.py
Unescape View File

@ -1,61 +1,52 @@
import pickle # Copyright (c) [2024] [Daniel Pozsar]
import warnings #
from sys import getsizeof # 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.
from timeit import default_timer as timer from timeit import default_timer as timer
import numpy as np
import sisl
import sisl.viz
from mpi4py import MPI
from numpy.linalg import inv
from tqdm import tqdm
from grogu_magn.utils import *
"""
# Some input parsing
parser = argparse.ArgumentParser()
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('--eset' , dest = 'eset' , default = 42 , type=int , help = 'Number of energy points on the contour')
parser.add_argument('--eset-p' , dest = 'esetp' , default = 10 , type=int , help = 'Parameter tuning the distribution on the contour')
parser.add_argument('--input' , dest = 'infile' , required = True , help = 'Input file name')
parser.add_argument('--output' , dest = 'outfile', required = True , help = 'Output file name')
parser.add_argument('--Ebot' , dest = 'Ebot' , default = -20.0 , type=float, help = 'Bottom energy of the contour')
parser.add_argument('--npairs' , dest = 'npairs' , default = 1 , type=int , help = 'Number of unitcell pairs in each direction for Jij calculation')
parser.add_argument('--adirs' , dest = 'adirs' , default = False , help = 'Definition of pair directions')
parser.add_argument('--use-tqdm', dest = 'usetqdm', default = 'not' , help = 'Use tqdm for progressbars or not')
parser.add_argument('--cutoff' , dest = 'cutoff' , default = 100.0 , type=float, help = 'Real space cutoff for pair generation in Angs')
parser.add_argument('--pairfile', dest = 'pairfile', default = False , help = 'File to read pair information')
args = parser.parse_args()
"""
# runtime information # runtime information
times = dict() times = dict()
times["start_time"] = timer() times["start_time"] = timer()
################################################################################
#################################### INPUT ##################################### import warnings
################################################################################ from sys import getsizeof
import sisl
from mpi4py import MPI
from src.grogu_magn import *
# input output stuff
######################################################################
######################################################################
######################################################################
path = ( path = (
"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf" "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf"
) )
outfile = "./Fe3GeTe2_benchmark_on_15k_300eset_orb_test3" outfile = "./Fe3GeTe2_notebook"
# this information needs to be given at the input!!
scf_xcf_orientation = np.array([0, 0, 1]) # z
# list of reference directions for around which we calculate the derivatives
# o is the quantization axis, v and w are two axes perpendicular to it
# at this moment the user has to supply o,v,w on the input.
# we can have some default for this
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])]),
]
magnetic_entities = [ magnetic_entities = [
dict(atom=3, l=2), dict(atom=3, l=2),
dict(atom=4, l=2), dict(atom=4, l=2),
dict(atom=5, l=1), dict(atom=5, l=2),
] ]
pairs = [ pairs = [
dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),
@ -68,15 +59,12 @@ pairs = [
dict(ai=1, aj=2, Ruc=np.array([-2, 0, 0])), dict(ai=1, aj=2, Ruc=np.array([-2, 0, 0])),
dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])), dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])),
] ]
# Brilloun zone sampling and Green function contour integral simulation_parameters = default_args
kset = 15
kdirs = "xy" ######################################################################
ebot = -13 ######################################################################
eset = 300 ######################################################################
esetp = 1000
################################################################################
#################################### INPUT #####################################
################################################################################
# MPI parameters # MPI parameters
comm = MPI.COMM_WORLD comm = MPI.COMM_WORLD
@ -84,28 +72,37 @@ size = comm.Get_size()
rank = comm.Get_rank() rank = comm.Get_rank()
root_node = 0 root_node = 0
# check versions for debugging
if rank == root_node:
try:
print(sisl.__version__)
except:
print("sisl version unknown.")
try:
print(np.__version__)
except:
print("numpy version unknown.")
# rename outfile # rename outfile
if not outfile.endswith(".pickle"): if not simulation_parameters["outfile"].endswith(".pickle"):
outfile += ".pickle" simulation_parameters["outfile"] += ".pickle"
simulation_parameters = dict( # if ebot is not given put it 0.1 eV under the smallest energy
path=path, if simulation_parameters["ebot"] is None:
outpath=outfile, try:
scf_xcf_orientation=scf_xcf_orientation, eigfile = simulation_parameters["infile"][:-3] + "EIG"
ref_xcf_orientations=ref_xcf_orientations, simulation_parameters["ebot"] = read_siesta_emin(eigfile) - 0.1
kset=kset, except:
kdirs=kdirs, print("Could not determine ebot.")
ebot=ebot, print("Parameter was not given and .EIG file was not found.")
eset=eset,
esetp=esetp,
parallel_size=size,
)
# digestion of the input
# read sile # read sile
fdf = sisl.get_sile(path) fdf = sisl.get_sile(simulation_parameters["infile"])
# read in hamiltonian # read in hamiltonian
dh = fdf.read_hamiltonian() dh = fdf.read_hamiltonian()
# read unit cell vectors
simulation_parameters["cell"] = fdf.read_geometry().cell simulation_parameters["cell"] = fdf.read_geometry().cell
# unit cell index # unit cell index
@ -119,60 +116,11 @@ if rank == root_node:
"================================================================================================================================================================" "================================================================================================================================================================"
) )
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) # reformat Hamltonian and Overlap matrix for manipulations
h12 += dh.tocsr(dh.M12i) * 1.0j hh, ss, NO = build_hh_ss(dh)
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")
)
# symmetrizing Hamiltonian and Overlap matrix to make them hermitian
# 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())
]
)
# symmetrizing Hamiltonian and overlap matrix to make them hermitian
for i in range(dh.lattice.sc_off.shape[0]): for i in range(dh.lattice.sc_off.shape[0]):
j = dh.lattice.sc_index(-dh.lattice.sc_off[i]) j = dh.lattice.sc_index(-dh.lattice.sc_off[i])
h1, h1d = hh[i], hh[j] h1, h1d = hh[i], hh[j]
@ -194,9 +142,9 @@ XCF = np.array(
np.array([f["y"] / 2 for f in traced]), np.array([f["y"] / 2 for f in traced]),
np.array([f["z"] / 2 for f in traced]), np.array([f["z"] / 2 for f in traced]),
] ]
) # equation 77 )
# Check if exchange field has scalar part # check if exchange field has scalar part
max_xcfs = abs(np.array(np.array([f["c"] / 2 for f in traced]))).max() max_xcfs = abs(np.array(np.array([f["c"] / 2 for f in traced]))).max()
if max_xcfs > 1e-12: if max_xcfs > 1e-12:
warnings.warn( warnings.warn(
@ -212,104 +160,10 @@ if rank == root_node:
"================================================================================================================================================================" "================================================================================================================================================================"
) )
# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions # initialize pairs and magnetic entities based on input information
for mag_ent in magnetic_entities: pairs, magnetic_entities = setup_pairs_and_magnetic_entities(
parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes magnetic_entities, pairs, dh, simulation_parameters
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 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((eset, spin_box_shape, spin_box_shape), dtype="complex128")
)
mag_ent["Gii_tmp"].append(
np.zeros((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 ref_xcf_orientations:
# Greens functions for every quantization axis
pair["Gij"].append(
np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype="complex128")
)
pair["Gij_tmp"].append(
np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype="complex128")
)
pair["Gji"].append(
np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128")
)
pair["Gji_tmp"].append(
np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128")
)
if rank == root_node: if rank == root_node:
times["site_and_pair_dictionaries_time"] = timer() times["site_and_pair_dictionaries_time"] = timer()
@ -320,10 +174,15 @@ if rank == root_node:
"================================================================================================================================================================" "================================================================================================================================================================"
) )
kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling # generate k space sampling
kset = make_kset(
dirs=simulation_parameters["kdirs"], NUMK=simulation_parameters["kset"]
)
wkset = np.ones(len(kset)) / len(kset) # generate weights for k points wkset = np.ones(len(kset)) / len(kset) # generate weights for k points
kpcs = np.array_split(kset, size) # split the k points based on MPI size kpcs = np.array_split(kset, size) # split the k points based on MPI size
kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop")
if tqdm_imported:
kpcs[root_node] = tqdm(kpcs[root_node], desc="k loop")
if rank == root_node: if rank == root_node:
times["k_set_time"] = timer() times["k_set_time"] = timer()
@ -331,225 +190,3 @@ if rank == root_node:
print( print(
"================================================================================================================================================================" "================================================================================================================================================================"
) )
# this will contain the three hamiltonians 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(ref_xcf_orientations):
# obtain rotated exchange field
R = RotMa2b(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 #######################################################################################
hamiltonians.append(
dict(orient=orient["o"], H=rot_H)
) # store orientation and rotated Hamiltonian
# these are the rotations (for now) perpendicular to the quantization axis
for u in orient["vw"]:
Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H
Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100
Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100
for mag_ent in magnetic_entities:
idx = mag_ent["spin_box_indeces"]
# fill up the perturbed potentials (for now) based on the on-site projections
mag_ent["Vu1"][i].append(Vu1[:, idx][idx, :])
mag_ent["Vu2"][i].append(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(
"================================================================================================================================================================"
)
if rank == root_node:
print("Starting matrix inversions")
print(f"Total number of k points: {kset.shape[0]}")
print(f"Number of energy samples per k point: {eset}")
print(f"Total number of directions: {len(hamiltonians)}")
print(
f"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * 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) * eset * 32} KB"
)
print(
f"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * eset * 32 / 1024} MB"
)
print(
"================================================================================================================================================================"
)
comm.Barrier()
# ----------------------------------------------------------------------
# make energy contour
# we are working in eV now !
# and sisl shifts E_F to 0 !
cont = make_contour(emin=ebot, enum=eset, p=esetp)
eran = cont.ze
# ----------------------------------------------------------------------
# sampling the integrand on the contour and the BZ
for k in kpcs[rank]:
wk = wkset[rank] # weight of k point in BZ integral
# iterate over reference directions
for i, hamiltonian_orientation in enumerate(hamiltonians):
# calculate Greens function
H = hamiltonian_orientation["H"]
HK, SK = hsk(H, ss, dh.sc_off, k)
# Gk = inv(SK * eran.reshape(eset, 1, 1) - HK)
# solve Greens function sequentially for the energies, because of memory bound
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)
# store the Greens function slice of the magnetic entities (for now) based on the on-site projections
for mag_ent in magnetic_entities:
mag_ent["Gii_tmp"][i] += (
Gk[:, mag_ent["spin_box_indeces"], :][:, :, mag_ent["spin_box_indeces"]]
* 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_indeces"]
aj = magnetic_entities[pair["aj"]]["spin_box_indeces"]
# store the Greens function slice of the magnetic entities (for now) based on the on-site projections
pair["Gij_tmp"][i] += Gk[:, ai][..., aj] * phase * wk
pair["Gji_tmp"][i] += 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)
# evaluation of the contour integral
storage.append(np.trapz(-1 / np.pi * np.imag(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"]
)
print("Magnetic entities integrated.")
# 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)
# evaluation of the contour integral
storage.append(np.trapz(-1 / np.pi * np.imag(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"])
print("Pairs integrated.")
# 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")
print("Magnetic parameters calculated.")
times["end_time"] = timer()
print(
"##################################################################### GROGU OUTPUT #############################################################################"
)
print_parameters(simulation_parameters)
print_atoms_and_pairs(magnetic_entities, pairs)
print_runtime_information(times)
# 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"]
# create output dictionary with all the relevant data
results = dict(
parameters=simulation_parameters,
magnetic_entities=magnetic_entities,
pairs=pairs,
runtime=times,
)
# save dictionary
with open(outfile, "wb") as output_file:
pickle.dump(results, output_file)

Write Preview
Loading…
Cancel
Save