@ -178,9 +178,14 @@ if rank == root_node:
kset = make_kset ( kset = make_kset (
dirs = simulation_parameters [ " kdirs " ] , NUMK = simulation_parameters [ " kset " ]
dirs = simulation_parameters [ " kdirs " ] , NUMK = simulation_parameters [ " kset " ]
) )
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
# 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 tqdm_imported : if tqdm_imported :
kpcs [ root_node ] = tqdm ( kpcs [ root_node ] , desc = " k loop " )
kpcs [ root_node ] = tqdm ( kpcs [ root_node ] , desc = " k loop " )
@ -190,3 +195,237 @@ if rank == root_node:
print (
print (
" ================================================================================================================================================================ "
" ================================================================================================================================================================ "
)
)
# this will contain the three Hamiltonian in the
# reference directions needed to calculate the energy
# variations upon rotation
hamiltonians = [ ]
# iterate over the reference directions (quantization axes)
for i , orient in enumerate ( simulation_parameters [ " ref_xcf_orientations " ] ) :
# obtain rotated exchange field and Hamiltonian
R = RotMa2b ( simulation_parameters [ " scf_xcf_orientation " ] , orient [ " o " ] )
rot_XCF = np . einsum ( " ij,jklm->iklm " , R , XCF )
rot_H_XCF = sum (
[ np . kron ( rot_XCF [ i ] , tau ) for i , tau in enumerate ( [ tau_x , tau_y , tau_z ] ) ]
)
rot_H_XCF_uc = rot_H_XCF [ uc_in_sc_idx ]
# obtain total Hamiltonian with the rotated exchange field
rot_H = hTRS + rot_H_XCF # equation 76
# store the relevant information of the Hamiltonian
hamiltonians . append ( dict ( orient = orient [ " o " ] , H = rot_H ) )
if simulation_parameters [ " calculate_charge " ] :
hamiltonians [ - 1 ] [ " GS " ] = np . zeros (
( simulation_parameters [ " eset " ] , rot_H . shape [ 1 ] , rot_H . shape [ 2 ] ) ,
dtype = " complex128 " ,
)
hamiltonians [ - 1 ] [ " GS_tmp " ] = np . zeros (
( simulation_parameters [ " eset " ] , rot_H . shape [ 1 ] , rot_H . shape [ 2 ] ) ,
dtype = " complex128 " ,
)
# these are the rotations (for now) perpendicular to the quantization axis
for u in orient [ " vw " ] :
# section 2.H
Tu = np . kron ( np . eye ( NO , dtype = int ) , tau_u ( u ) )
Vu1 , Vu2 = calc_Vu ( rot_H_XCF_uc , Tu )
for mag_ent in magnetic_entities :
idx = mag_ent [ " spin_box_indices " ]
# fill up the perturbed potentials (for now) based on the on-site projections
mag_ent [ " Vu1 " ] [ i ] . append ( onsite_projection ( Vu1 , idx ) )
mag_ent [ " Vu2 " ] [ i ] . append ( onsite_projection ( Vu2 , 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. " )
print ( f " Total number of k points: { kset . shape [ 0 ] } " )
print ( f " Number of energy samples per k point: { simulation_parameters [ ' eset ' ] } " )
print ( f " Total number of directions: { len ( hamiltonians ) } " )
print (
f " Total number of matrix inversions: { kset . shape [ 0 ] * len ( hamiltonians ) * simulation_parameters [ ' eset ' ] } "
)
print ( f " The shape of the Hamiltonian and the Greens function is { NO } x { NO } = { NO * NO } " )
# https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv
# memory is O(64 n**2) for complex matrices
memory_size = getsizeof ( hamiltonians [ 0 ] [ " H " ] . base ) / 1024
print (
f " Memory taken by a single Hamiltonian is: { getsizeof ( hamiltonians [ 0 ] [ ' H ' ] . base ) / 1024 } KB "
)
print ( f " Expected memory usage per matrix inversion: { memory_size * 32 } KB " )
print (
f " Expected memory usage per k point for parallel inversion: { memory_size * len ( hamiltonians ) * simulation_parameters [ ' eset ' ] * 32 } KB "
)
print (
f " Expected memory usage on root node: { len ( np . array_split ( kset , size ) [ 0 ] ) * memory_size * len ( hamiltonians ) * simulation_parameters [ ' eset ' ] * 32 / 1024 } MB "
)
print (
" ================================================================================================================================================================ "
)
# MPI barrier
comm . Barrier ( )
# make energy contour
cont = make_contour (
emin = simulation_parameters [ " ebot " ] ,
enum = simulation_parameters [ " eset " ] ,
p = simulation_parameters [ " esetp " ] ,
)
eran = cont . ze
# sampling the integrand on the contour and the BZ
for k in kpcs [ rank ] :
# weight of k point in BZ integral
wk = wkset [ rank ]
# iterate over reference directions
for i , hamiltonian_orientation in enumerate ( hamiltonians ) :
# calculate Hamiltonian and Overlap matrix in a given k point
H = hamiltonian_orientation [ " H " ]
HK , SK = hsk ( H , ss , dh . sc_off , k )
if simulation_parameters [ " parallel_solver_for_Gk " ] :
Gk = parallel_Gk ( HK , SK , eran , simulation_parameters [ " eset " ] )
else :
# solve Greens function sequentially for the energies, because of memory bound
Gk = sequential_GK ( HK , SK , eran , simulation_parameters [ " eset " ] )
# saving this for total charge
if simulation_parameters [ " calculate_charge " ] :
hamiltonian_orientation [ " GS_tmp " ] + = Gk @ SK * wk
# store the Greens function slice of the magnetic entities
for mag_ent in magnetic_entities :
idx = mag_ent [ " spin_box_indices " ]
mag_ent [ " Gii_tmp " ] [ i ] + = onsite_projection ( Gk , idx ) * wk
for pair in pairs :
# add phase shift based on the cell difference
phase = np . exp ( 1 j * 2 * np . pi * k @ pair [ " Ruc " ] . T )
# get the pair orbital sizes from the magnetic entities
ai = magnetic_entities [ pair [ " ai " ] ] [ " spin_box_indices " ]
aj = magnetic_entities [ pair [ " aj " ] ] [ " spin_box_indices " ]
# store the Greens function slice of the magnetic entities
pair [ " Gij_tmp " ] [ i ] + = onsite_projection ( Gk , ai , aj ) * phase * wk
pair [ " Gji_tmp " ] [ i ] + = onsite_projection ( Gk , aj , ai ) / phase * wk
# summ reduce partial results of mpi nodes
for i in range ( len ( hamiltonians ) ) :
# for total charge calculation
if simulation_parameters [ " calculate_charge " ] :
comm . Reduce ( hamiltonians [ i ] [ " GS_tmp " ] , hamiltonians [ i ] [ " GS " ] , root = root_node )
for mag_ent in magnetic_entities :
comm . Reduce ( mag_ent [ " Gii_tmp " ] [ i ] , mag_ent [ " Gii " ] [ i ] , root = root_node )
for pair in pairs :
comm . Reduce ( pair [ " Gij_tmp " ] [ i ] , pair [ " Gij " ] [ i ] , root = root_node )
comm . Reduce ( pair [ " Gji_tmp " ] [ i ] , pair [ " Gji " ] [ i ] , root = root_node )
if rank == root_node :
times [ " green_function_inversion_time " ] = timer ( )
print (
f " Calculated Greens functions. Elapsed time: { times [ ' green_function_inversion_time ' ] } s "
)
print (
" ================================================================================================================================================================ "
)
if rank == root_node :
# Calculate total charge
if simulation_parameters [ " calculate_charge " ] :
for hamiltonian in hamiltonians :
GS = hamiltonian [ " GS " ]
traced = np . trace ( ( GS ) , axis1 = 1 , axis2 = 2 )
simulation_parameters [ " charges " ] . append ( int_de_ke ( traced , cont . we ) )
# iterate over the magnetic entities
for tracker , mag_ent in enumerate ( magnetic_entities ) :
# iterate over the quantization axes
for i , Gii in enumerate ( mag_ent [ " Gii " ] ) :
storage = [ ]
# iterate over the first and second order local perturbations
for Vu1 , Vu2 in zip ( mag_ent [ " Vu1 " ] [ i ] , mag_ent [ " Vu2 " ] [ i ] ) :
# The Szunyogh-Lichtenstein formula
traced = np . trace (
( Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii ) , axis1 = 1 , axis2 = 2
) # this is the on site projection
# evaluation of the contour integral
storage . append ( int_de_ke ( traced , cont . we ) )
# fill up the magnetic entities dictionary with the energies
magnetic_entities [ tracker ] [ " energies " ] . append ( storage )
# convert to np array
magnetic_entities [ tracker ] [ " energies " ] = np . array (
magnetic_entities [ tracker ] [ " energies " ]
)
# iterate over the pairs
for tracker , pair in enumerate ( pairs ) :
# iterate over the quantization axes
for i , ( Gij , Gji ) in enumerate ( zip ( pair [ " Gij " ] , pair [ " Gji " ] ) ) :
site_i = magnetic_entities [ pair [ " ai " ] ]
site_j = magnetic_entities [ pair [ " aj " ] ]
storage = [ ]
# iterate over the first order local perturbations in all possible orientations for the two sites
for Vui in site_i [ " Vu1 " ] [ i ] :
for Vuj in site_j [ " Vu1 " ] [ i ] :
# The Szunyogh-Lichtenstein formula
traced = np . trace (
( Vui @ Gij @ Vuj @ Gji ) , axis1 = 1 , axis2 = 2
) # this is the on site projection
# evaluation of the contour integral
storage . append ( int_de_ke ( traced , cont . we ) )
# fill up the pairs dictionary with the energies
pairs [ tracker ] [ " energies " ] . append ( storage )
# convert to np array
pairs [ tracker ] [ " energies " ] = np . array ( pairs [ tracker ] [ " energies " ] )
# calculate magnetic parameters
for mag_ent in magnetic_entities :
Kxx , Kyy , Kzz , consistency = calculate_anisotropy_tensor ( mag_ent )
mag_ent [ " K " ] = np . array ( [ Kxx , Kyy , Kzz ] ) * sisl . unit_convert ( " eV " , " meV " )
mag_ent [ " K_consistency " ] = consistency
for pair in pairs :
J_iso , J_S , D , J = calculate_exchange_tensor ( pair )
pair [ " J_iso " ] = J_iso * sisl . unit_convert ( " eV " , " meV " )
pair [ " J_S " ] = J_S * sisl . unit_convert ( " eV " , " meV " )
pair [ " D " ] = D * sisl . unit_convert ( " eV " , " meV " )
pair [ " J " ] = J * sisl . unit_convert ( " eV " , " meV " )
times [ " end_time " ] = timer ( )
print (
" ##################################################################### GROGU OUTPUT ############################################################################# "
)
print_parameters ( simulation_parameters )
print_atoms_and_pairs ( magnetic_entities , pairs )
print_runtime_information ( times )
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_pickle ( simulation_parameters [ " outfile " ] , results )
