@ -42,16 +42,16 @@ from mpi4py import MPI
try :
try :
from tqdm import tqdm
from tqdm import tqdm
tqdm_imported = True
tqdm_imported : bool = True
except :
except :
print ( " Please install tqdm for nice progress bar. " )
print ( " Please install tqdm for nice progress bar. " )
tqdm_imported = False
tqdm_imported : bool = False
from grogupy import *
from grogupy import *
def parse_command_line ( ) :
def parse_command_line ( ) - > dict :
""" This function can read input from the command line. """
""" This function can read input from the command line. """
parser = ArgumentParser ( )
parser = ArgumentParser ( )
@ -59,12 +59,19 @@ def parse_command_line():
parser . add_argument (
parser . add_argument (
" -i " ,
" -i " ,
" --input " ,
" --input " ,
dest = " infile" ,
dest = " grogupy_ infile" ,
default = None ,
default = None ,
type = str ,
type = str ,
help = " Input file name " ,
help = " Input file name for the Grogupy fdf file " ,
required = True ,
required = True ,
)
)
parser . add_argument (
" --siesta-input " ,
dest = " infile " ,
default = None ,
type = str ,
help = " Input file name for the Siesta fdf file " ,
)
parser . add_argument (
parser . add_argument (
" -o " ,
" -o " ,
" --output " ,
" --output " ,
@ -129,16 +136,16 @@ def parse_command_line():
return cmd_line_args
return cmd_line_args
def main ( simulation_parameters , magnetic_entities , pairs ):
def main ( simulation_parameters : dict , magnetic_entities : list , pairs : list ) - > None :
# runtime information
# runtime information
times = dict ( )
times : dict = dict ( )
times [ " start_time " ] = timer ( )
times [ " start_time " ] = timer ( )
# MPI parameters
# MPI parameters
comm = MPI . COMM_WORLD
comm = MPI . COMM_WORLD
size = comm . Get_size ( )
size : int = comm . Get_size ( )
rank = comm . Get_rank ( )
rank : int = comm . Get_rank ( )
root_node = 0
root_node : int = 0
if rank == root_node :
if rank == root_node :
# include parallel size in simulation parameters
# include parallel size in simulation parameters
@ -172,8 +179,8 @@ def main(simulation_parameters, magnetic_entities, pairs):
# add third direction for off-diagonal anisotropy
# add third direction for off-diagonal anisotropy
for orientation in simulation_parameters [ " ref_xcf_orientations " ] :
for orientation in simulation_parameters [ " ref_xcf_orientations " ] :
o1 = orientation [ " vw " ] [ 0 ]
o1 : np . array = orientation [ " vw " ] [ 0 ]
o2 = orientation [ " vw " ] [ 1 ]
o2 : np . array = orientation [ " vw " ] [ 1 ]
orientation [ " vw " ] . append ( ( o1 + o2 ) / np . sqrt ( 2 ) )
orientation [ " vw " ] . append ( ( o1 + o2 ) / np . sqrt ( 2 ) )
# read sile
# read sile
@ -186,7 +193,7 @@ def main(simulation_parameters, magnetic_entities, pairs):
simulation_parameters [ " cell " ] = fdf . read_geometry ( ) . cell
simulation_parameters [ " cell " ] = fdf . read_geometry ( ) . cell
# unit cell index
# unit cell index
uc_in_sc_idx = dh . lattice . sc_index ( [ 0 , 0 , 0 ] )
uc_in_sc_idx : int = dh . lattice . sc_index ( [ 0 , 0 , 0 ] )
if rank == root_node :
if rank == root_node :
print ( " \n \n \n \n \n " )
print ( " \n \n \n \n \n " )
@ -205,9 +212,9 @@ def main(simulation_parameters, magnetic_entities, pairs):
" ================================================================================================================================================================ "
" ================================================================================================================================================================ "
)
)
NO = dh . no # shorthand for number of orbitals in the unit cell
NO : int = dh . no # shorthand for number of orbitals in the unit cell
# reformat Ham ltonian and Overlap matrix for manipulations
# reformat Ham i ltonian and Overlap matrix for manipulations
hh , ss = build_hh_ss ( dh )
hh , ss = build_hh_ss ( dh )
# symmetrizing Hamiltonian and Overlap matrix to make them hermitian
# symmetrizing Hamiltonian and Overlap matrix to make them hermitian
@ -219,14 +226,18 @@ def main(simulation_parameters, magnetic_entities, pairs):
ss [ i ] , ss [ j ] = ( s1 + s1d . T . conj ( ) ) / 2 , ( s1d + s1 . T . conj ( ) ) / 2
ss [ i ] , ss [ j ] = ( s1 + s1d . T . conj ( ) ) / 2 , ( s1d + s1 . T . conj ( ) ) / 2
# identifying TRS and TRB parts of the Hamiltonian
# identifying TRS and TRB parts of the Hamiltonian
TAUY = np . kron ( np . eye ( NO ) , TAU_Y )
TAUY : np . array = np . kron ( np . eye ( NO ) , TAU_Y )
hTR = np . array ( [ TAUY @ hh [ i ] . conj ( ) @ TAUY for i in range ( dh . lattice . nsc . prod ( ) ) ] )
hTR : np . array = np . array (
hTRS = ( hh + hTR ) / 2
[ TAUY @ hh [ i ] . conj ( ) @ TAUY for i in range ( dh . lattice . nsc . prod ( ) ) ]
hTRB = ( hh - hTR ) / 2
)
hTRS : np . array = ( hh + hTR ) / 2
hTRB : np . array = ( hh - hTR ) / 2
# extracting the exchange field
# extracting the exchange field
traced = [ spin_tracer ( hTRB [ i ] ) for i in range ( dh . lattice . nsc . prod ( ) ) ] # equation 77
traced : np . array = [
XCF = np . array (
spin_tracer ( hTRB [ i ] ) for i in range ( dh . lattice . nsc . prod ( ) )
] # equation 77
XCF : np . array = np . array (
[
[
np . array ( [ f [ " x " ] / 2 for f in traced ] ) ,
np . array ( [ f [ " x " ] / 2 for f in traced ] ) ,
np . array ( [ f [ " y " ] / 2 for f in traced ] ) ,
np . array ( [ f [ " y " ] / 2 for f in traced ] ) ,
@ -235,7 +246,7 @@ def main(simulation_parameters, magnetic_entities, pairs):
)
)
# 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 : float = 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 (
f " Exchange field has non negligible scalar part. Largest value is { max_xcfs } "
f " Exchange field has non negligible scalar part. Largest value is { max_xcfs } "
@ -270,10 +281,10 @@ def main(simulation_parameters, magnetic_entities, pairs):
)
)
# generate weights for k points
# generate weights for k points
wkset = np . ones ( len ( kset ) ) / len ( kset )
wkset : np . array = np . ones ( len ( kset ) ) / len ( kset )
# split the k points based on MPI size
# split the k points based on MPI size
kpcs = np . array_split ( kset , size )
kpcs : np . array = np . array_split ( kset , size )
# use progress bar if available
# use progress bar if available
if rank == root_node and tqdm_imported :
if rank == root_node and tqdm_imported :
@ -289,20 +300,20 @@ def main(simulation_parameters, magnetic_entities, pairs):
# this will contain the three Hamiltonian in the
# this will contain the three Hamiltonian in the
# reference directions needed to calculate the energy
# reference directions needed to calculate the energy
# variations upon rotation
# variations upon rotation
hamiltonians = [ ]
hamiltonians : list = [ ]
# iterate over the reference directions (quantization axes)
# iterate over the reference directions (quantization axes)
for i , orient in enumerate ( simulation_parameters [ " ref_xcf_orientations " ] ) :
for i , orient in enumerate ( simulation_parameters [ " ref_xcf_orientations " ] ) :
# obtain rotated exchange field and Hamiltonian
# obtain rotated exchange field and Hamiltonian
R = RotMa2b ( simulation_parameters [ " scf_xcf_orientation " ] , orient [ " o " ] )
R : np . array = RotMa2b ( simulation_parameters [ " scf_xcf_orientation " ] , orient [ " o " ] )
rot_XCF = np . einsum ( " ij,jklm->iklm " , R , XCF )
rot_XCF : np . array = np . einsum ( " ij,jklm->iklm " , R , XCF )
rot_H_XCF = sum (
rot_H_XCF : np . array = sum (
[ np . kron ( rot_XCF [ i ] , tau ) for i , tau in enumerate ( [ TAU_X , TAU_Y , TAU_Z ] ) ]
[ 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 ]
rot_H_XCF_uc : np . array = rot_H_XCF [ uc_in_sc_idx ]
# obtain total Hamiltonian with the rotated exchange field
# obtain total Hamiltonian with the rotated exchange field
rot_H = hTRS + rot_H_XCF # equation 76
rot_H : np . array = hTRS + rot_H_XCF # equation 76
# store the relevant information of the Hamiltonian
# store the relevant information of the Hamiltonian
hamiltonians . append ( dict ( orient = orient [ " o " ] , H = rot_H ) )
hamiltonians . append ( dict ( orient = orient [ " o " ] , H = rot_H ) )
@ -310,11 +321,11 @@ def main(simulation_parameters, magnetic_entities, pairs):
# these are the rotations (for now) perpendicular to the quantization axis
# these are the rotations (for now) perpendicular to the quantization axis
for u in orient [ " vw " ] :
for u in orient [ " vw " ] :
# section 2.H
# section 2.H
Tu = np . kron ( np . eye ( NO , dtype = int ) , tau_u ( u ) )
Tu : np . array = np . kron ( np . eye ( NO , dtype = int ) , tau_u ( u ) )
Vu1 , Vu2 = calc_Vu ( rot_H_XCF_uc , Tu )
Vu1 , Vu2 = calc_Vu ( rot_H_XCF_uc , Tu )
for mag_ent in magnetic_entities :
for mag_ent in magnetic_entities :
idx = mag_ent [ " spin_box_indices " ]
idx : int = mag_ent [ " spin_box_indices " ]
# fill up the perturbed potentials (for now) based on the on-site projections
# 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 [ " Vu1 " ] [ i ] . append ( onsite_projection ( Vu1 , idx , idx ) )
mag_ent [ " Vu2 " ] [ i ] . append ( onsite_projection ( Vu2 , idx , idx ) )
mag_ent [ " Vu2 " ] [ i ] . append ( onsite_projection ( Vu2 , idx , idx ) )
@ -347,7 +358,7 @@ def main(simulation_parameters, magnetic_entities, pairs):
)
)
# https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv
# https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv
# memory is O(64 n**2) for complex matrices
# memory is O(64 n**2) for complex matrices
memory_size = getsizeof ( hamiltonians [ 0 ] [ " H " ] . base ) / 1024
memory_size : int = getsizeof ( hamiltonians [ 0 ] [ " H " ] . base ) / 1024
print (
print (
f " Memory taken by a single Hamiltonian is: { getsizeof ( hamiltonians [ 0 ] [ ' H ' ] . base ) / 1024 } KB "
f " Memory taken by a single Hamiltonian is: { getsizeof ( hamiltonians [ 0 ] [ ' H ' ] . base ) / 1024 } KB "
)
)
@ -376,12 +387,12 @@ def main(simulation_parameters, magnetic_entities, pairs):
# sampling the integrand on the contour and the BZ
# sampling the integrand on the contour and the BZ
for k in kpcs [ rank ] :
for k in kpcs [ rank ] :
# weight of k point in BZ integral
# weight of k point in BZ integral
wk = wkset [ rank ]
wk : float = wkset [ rank ]
# iterate over reference directions
# iterate over reference directions
for i , hamiltonian_orientation in enumerate ( hamiltonians ) :
for i , hamiltonian_orientation in enumerate ( hamiltonians ) :
# calculate Hamiltonian and Overlap matrix in a given k point
# calculate Hamiltonian and Overlap matrix in a given k point
H = hamiltonian_orientation [ " H " ]
H : np . array = hamiltonian_orientation [ " H " ]
HK , SK = hsk ( H , ss , dh . sc_off , k )
HK , SK = hsk ( H , ss , dh . sc_off , k )
if simulation_parameters [ " parallel_solver_for_Gk " ] :
if simulation_parameters [ " parallel_solver_for_Gk " ] :
@ -392,22 +403,22 @@ def main(simulation_parameters, magnetic_entities, pairs):
# store the Greens function slice of the magnetic entities
# store the Greens function slice of the magnetic entities
for mag_ent in magnetic_entities :
for mag_ent in magnetic_entities :
idx = mag_ent [ " spin_box_indices " ]
idx : int = mag_ent [ " spin_box_indices " ]
mag_ent [ " Gii_tmp " ] [ i ] + = onsite_projection ( Gk , idx , idx ) * wk
mag_ent [ " Gii_tmp " ] [ i ] + = onsite_projection ( Gk , idx , idx ) * wk
for pair in pairs :
for pair in pairs :
# add phase shift based on the cell difference
# add phase shift based on the cell difference
phase = np . exp ( 1 j * 2 * np . pi * k @ pair [ " Ruc " ] . T )
phase : np . array = np . exp ( 1 j * 2 * np . pi * k @ pair [ " Ruc " ] . T )
# get the pair orbital sizes from the magnetic entities
# get the pair orbital sizes from the magnetic entities
ai = magnetic_entities [ pair [ " ai " ] ] [ " spin_box_indices " ]
ai : int = magnetic_entities [ pair [ " ai " ] ] [ " spin_box_indices " ]
aj = magnetic_entities [ pair [ " aj " ] ] [ " spin_box_indices " ]
aj : int = magnetic_entities [ pair [ " aj " ] ] [ " spin_box_indices " ]
# store the Greens function slice of the magnetic entities
# store the Greens function slice of the magnetic entities
pair [ " Gij_tmp " ] [ i ] + = onsite_projection ( Gk , ai , aj ) * phase * wk
pair [ " Gij_tmp " ] [ i ] + = onsite_projection ( Gk , ai , aj ) * phase * wk
pair [ " Gji_tmp " ] [ i ] + = onsite_projection ( Gk , aj , ai ) / phase * wk
pair [ " Gji_tmp " ] [ i ] + = onsite_projection ( Gk , aj , ai ) / phase * wk
# sum m reduce partial results of mpi nodes
# sum reduce partial results of mpi nodes
for i in range ( len ( hamiltonians ) ) :
for i in range ( len ( hamiltonians ) ) :
for mag_ent in magnetic_entities :
for mag_ent in magnetic_entities :
comm . Reduce ( mag_ent [ " Gii_tmp " ] [ i ] , mag_ent [ " Gii " ] [ i ] , root = root_node )
comm . Reduce ( mag_ent [ " Gii_tmp " ] [ i ] , mag_ent [ " Gii " ] [ i ] , root = root_node )
@ -430,11 +441,11 @@ def main(simulation_parameters, magnetic_entities, pairs):
for tracker , mag_ent in enumerate ( magnetic_entities ) :
for tracker , mag_ent in enumerate ( magnetic_entities ) :
# iterate over the quantization axes
# iterate over the quantization axes
for i , Gii in enumerate ( mag_ent [ " Gii " ] ) :
for i , Gii in enumerate ( mag_ent [ " Gii " ] ) :
storage = [ ]
storage : list = [ ]
# iterate over the first and second order local perturbations
# iterate over the first and second order local perturbations
for Vu1 , Vu2 in zip ( mag_ent [ " Vu1 " ] [ i ] , mag_ent [ " Vu2 " ] [ i ] ) :
for Vu1 , Vu2 in zip ( mag_ent [ " Vu1 " ] [ i ] , mag_ent [ " Vu2 " ] [ i ] ) :
# The Szunyogh-Lichtenstein formula
# The Szunyogh-Lichtenstein formula
traced = np . trace (
traced : np . array = np . trace (
( Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii ) , axis1 = 1 , axis2 = 2
( Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii ) , axis1 = 1 , axis2 = 2
) # this is the on site projection
) # this is the on site projection
# evaluation of the contour integral
# evaluation of the contour integral
@ -450,16 +461,16 @@ def main(simulation_parameters, magnetic_entities, pairs):
for tracker , pair in enumerate ( pairs ) :
for tracker , pair in enumerate ( pairs ) :
# iterate over the quantization axes
# iterate over the quantization axes
for i , ( Gij , Gji ) in enumerate ( zip ( pair [ " Gij " ] , pair [ " Gji " ] ) ) :
for i , ( Gij , Gji ) in enumerate ( zip ( pair [ " Gij " ] , pair [ " Gji " ] ) ) :
site_i = magnetic_entities [ pair [ " ai " ] ]
site_i : int = magnetic_entities [ pair [ " ai " ] ]
site_j = magnetic_entities [ pair [ " aj " ] ]
site_j : int = magnetic_entities [ pair [ " aj " ] ]
storage = [ ]
storage : list = [ ]
# iterate over the first order local perturbations in all possible orientations for the two sites
# iterate over the first order local perturbations in all possible orientations for the two sites
# actually all possible orientations without the orientation for the off-diagonal anisotropy
# actually all possible orientations without the orientation for the off-diagonal anisotropy
# that is why we only take the first two of each Vu1
# that is why we only take the first two of each Vu1
for Vui in site_i [ " Vu1 " ] [ i ] [ : 2 ] :
for Vui in site_i [ " Vu1 " ] [ i ] [ : 2 ] :
for Vuj in site_j [ " Vu1 " ] [ i ] [ : 2 ] :
for Vuj in site_j [ " Vu1 " ] [ i ] [ : 2 ] :
# The Szunyogh-Lichtenstein formula
# The Szunyogh-Lichtenstein formula
traced = np . trace (
traced : np . array = np . trace (
( Vui @ Gij @ Vuj @ Gji ) , axis1 = 1 , axis2 = 2
( Vui @ Gij @ Vuj @ Gji ) , axis1 = 1 , axis2 = 2
) # this is the on site projection
) # this is the on site projection
# evaluation of the contour integral
# evaluation of the contour integral
@ -499,7 +510,7 @@ def main(simulation_parameters, magnetic_entities, pairs):
# remove unwanted stuff before saving
# remove unwanted stuff before saving
pairs , magnetic_entities = remove_clutter_for_save ( pairs , magnetic_entities )
pairs , magnetic_entities = remove_clutter_for_save ( pairs , magnetic_entities )
# create output dictionary with all the relevant data
# create output dictionary with all the relevant data
results = dict (
results : dict = dict (
parameters = simulation_parameters ,
parameters = simulation_parameters ,
magnetic_entities = magnetic_entities ,
magnetic_entities = magnetic_entities ,
pairs = pairs ,
pairs = pairs ,
@ -514,28 +525,25 @@ def main(simulation_parameters, magnetic_entities, pairs):
if __name__ == " __main__ " :
if __name__ == " __main__ " :
# loading parameters
# loading parameters
# it is not clear how to give grogu.fdf path...
command_line_arguments = parse_command_line ( )
# command_line_arguments = parse_command_line()
fdf_arguments , magnetic_entities , pairs = read_grogupy_fdf (
# fdf_arguments, magnetic_entities, pairs = read_fdf(command_line_arguments["infile"])
command_line_arguments [ " grogupy_infile " ]
)
# right now we do not use any of these input, but it shows
# the order of priority when processing arguments
# combine the inputs to a single dictionary
default_arguments = False
simulation_parameters : Final [ dict ] = process_input_args (
fdf_arguments = False
DEFAULT_ARGUMENTS , fdf_arguments , command_line_arguments , ACCEPTED_INPUTS
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 #
# This is the input file for now #
####################################################################################################
####################################################################################################
magnetic_entities = [
magnetic_entities : list = [
dict ( atom = 3 , l = 2 ) ,
dict ( atom = 3 , l = 2 ) ,
dict ( atom = 4 , l = 2 ) ,
dict ( atom = 4 , l = 2 ) ,
dict ( atom = 5 , l = 2 ) ,
dict ( atom = 5 , l = 2 ) ,
]
]
pairs = [
pairs : list = [
dict ( ai = 0 , aj = 1 , Ruc = np . array ( [ 0 , 0 , 0 ] ) ) ,
dict ( ai = 0 , aj = 1 , Ruc = np . array ( [ 0 , 0 , 0 ] ) ) ,
dict ( ai = 0 , aj = 2 , 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 = 1 , aj = 2 , Ruc = np . array ( [ 0 , 0 , 0 ] ) ) ,
@ -546,7 +554,7 @@ if __name__ == "__main__":
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 ] ) ) ,
]
]
simulation_parameters = dict ( )
simulation_parameters : dict = dict ( )
simulation_parameters [ " infile " ] = (
simulation_parameters [ " infile " ] = (
" /Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf "
" /Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf "
)
)