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

enclosed in main function and correcte d useful import

Browse Source
class-solution
Daniel Pozsar 3 months ago
parent 4b091373d9
commit bd3ec608c4
1 changed files with 359 additions and 344 deletions
Show Stats Download Patch File Download Diff File

313
src/grogu_magn/grogu.py
Unescape View File

@ -7,68 +7,86 @@ import sisl
from mpi4py import MPI
from numpy.linalg import inv
from tqdm import tqdm
from src.grogu_magn.useful import *
def main():
start_time = timer() # runtime information
times = dict()
times["start_time"] = timer()
# this cell mimicks an input file
fdf = sisl.get_sile("./lat3_791/Fe3GeTe2.fdf")
# 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 = [
from useful import *
# runtime information
times = dict()
times["start_time"] = timer()
# this cell mimicks an input file
fdf = sisl.get_sile("./Jij_for_Marci_6p45ang/CrBr.fdf") # ./lat3_791/Fe3GeTe2.fdf
# 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])]),
]
]
# human readable definition of magnetic entities
magnetic_entities = [
"""
# human readable definition of magnetic entities ./lat3_791/Fe3GeTe2.fdf
magnetic_entities = [
dict(atom=3, l=2),
dict(atom=4, l=2),
dict(atom=5, l=2),
dict(
atom=[3, 4],
),
]
# pair information
# these should all be around -41.9 in the isotropic part
# isotropic should be -82 meV
pairs = [
dict(ai=0, aj=3, Ruc=np.array([0, 0, 0])),
]
# pair information ./lat3_791/Fe3GeTe2.fdf
pairs = [
dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV
dict(
ai=0, aj=2, Ruc=np.array([0, 0, 0])
), # these should all be around -41.9 in the isotropic part
dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])),
dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])),
dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])),
] """
# human readable definition of magnetic entities ./Jij_for_Marci_6p45ang/CrBr.fdf
magnetic_entities = [
dict(atom=0, l=2),
dict(atom=1, l=2),
dict(atom=2, l=2),
]
# pair information ./Jij_for_Marci_6p45ang/CrBr.fdf
pairs = [
dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),
dict(ai=1, aj=0, 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=1, Ruc=np.array([1, 0, 0])),
dict(ai=0, aj=2, Ruc=np.array([1, 0, 0])),
dict(ai=0, aj=1, Ruc=np.array([-1, 0, 0])),
dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])),
dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])),
]
# Brilloun zone sampling and Green function contour integral
kset = 20
kdirs = "xy"
ebot = -30
eset = 50
esetp = 10000
# MPI parameters
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
root_node = 0
if rank == root_node:
dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])),
dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])),
dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])),
dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])),
]
# Brilloun zone sampling and Green function contour integral
kset = 20
kdirs = "xy"
ebot = -30
eset = 100
esetp = 10000
# MPI parameters
comm = MPI.COMM_WORLD
size = comm.Get_size()
rank = comm.Get_rank()
root_node = 0
if rank == root_node:
print("Number of nodes in the parallel cluster: ", size)
simulation_parameters = dict(
simulation_parameters = dict(
path="Not yet specified.",
scf_xcf_orientation=scf_xcf_orientation,
ref_xcf_orientations=ref_xcf_orientations,
@ -78,111 +96,109 @@ def main():
eset=eset,
esetp=esetp,
parallel_size=size,
)
)
# digestion of the input
# read in hamiltonian
dh = fdf.read_hamiltonian()
try:
# digestion of the input
# read in hamiltonian
dh = fdf.read_hamiltonian()
try:
simulation_parameters["geom"] = fdf.read_geometry()
except:
except:
print("Error reading geometry.")
# unit cell index
uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])
# unit cell index
uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])
times["setup_time"] = timer()
times["setup_time"] = timer()
NO = dh.no # shorthand for number of orbitals in the unit cell
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")
# 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")
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")
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")
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 = (
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(
# 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
U.T @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) @ U
for i in range(dh.lattice.nsc.prod())
]
), np.array(
), np.array(
[
U.T
@ np.block(
[[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]]
)
@ 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]):
# 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
# 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(
# 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"] for f in traced]),
np.array([f["y"] for f in traced]),
np.array([f["z"] for f in traced]),
]
) # equation 77
) # equation 77
# Check if exchange field has scalar part
max_xcfs = abs(np.array(np.array([f["c"] for f in traced]))).max()
if max_xcfs > 1e-12:
# Check if exchange field has scalar part
max_xcfs = abs(np.array(np.array([f["c"] 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}"
)
times["H_and_XCF_time"] = timer()
times["H_and_XCF_time"] = timer()
# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions
for i, mag_ent in enumerate(magnetic_entities):
# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions
for i, mag_ent in enumerate(magnetic_entities):
parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes
magnetic_entities[i]["orbital_indeces"] = parsed
# calculate spin box indexes
@ -190,9 +206,7 @@ def main():
# 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
mag_ent["energies"] = [] # we will store the second order energy derivations here
mag_ent["Gii"] = [] # Greens function
mag_ent["Gii_tmp"] = [] # Greens function for parallelization
@ -208,9 +222,9 @@ def main():
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:
# 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 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"])
@ -222,12 +236,13 @@ def main():
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")
) # Greens functions for every quantization axis
)
pair["Gji"].append(
np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128")
)
@ -235,20 +250,20 @@ def main():
np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128")
)
times["site_and_pair_dictionaries_time"] = timer()
times["site_and_pair_dictionaries_time"] = timer()
kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling
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[root_node] = tqdm(kpcs[root_node], desc="k loop", file=stdout)
kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling
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[root_node] = tqdm(kpcs[root_node], desc="k loop", file=stdout)
times["k_set_time"] = timer()
times["k_set_time"] = timer()
# this will contain the three hamiltonians in the reference directions needed to calculate the energy variations upon rotation
hamiltonians = []
# 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):
# 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)
@ -274,33 +289,34 @@ def main():
Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100
for mag_ent in magnetic_entities:
# fill up the perturbed potentials (for now) based on the on-site projections
mag_ent["Vu1"][i].append(
Vu1[:, mag_ent["spin_box_indeces"]][mag_ent["spin_box_indeces"], :]
) # fill up the perturbed potentials (for now) based on the on-site projections
)
mag_ent["Vu2"][i].append(
Vu2[:, mag_ent["spin_box_indeces"]][mag_ent["spin_box_indeces"], :]
)
times["reference_rotations_time"] = timer()
times["reference_rotations_time"] = timer()
if rank == root_node:
if rank == root_node:
print("Number of magnetic entities being calculated: ", len(magnetic_entities))
print(
"We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site."
)
print(f"The shape of the Hamiltonian and the Greens function is {NO}x{NO}.")
comm.Barrier()
# ----------------------------------------------------------------------
# make energy contour
# we are working in eV now !
# and sisil 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]:
comm.Barrier()
# ----------------------------------------------------------------------
# make energy contour
# we are working in eV now !
# and sisil 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):
@ -309,6 +325,11 @@ def main():
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] += (
@ -328,8 +349,8 @@ def main():
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)):
# 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)
@ -337,9 +358,9 @@ def main():
comm.Reduce(pair["Gij_tmp"][i], pair["Gij"][i], root=root_node)
comm.Reduce(pair["Gji_tmp"][i], pair["Gji"][i], root=root_node)
times["green_function_inversion_time"] = timer()
times["green_function_inversion_time"] = timer()
if rank == root_node:
if rank == root_node:
# iterate over the magnetic entities
for tracker, mag_ent in enumerate(magnetic_entities):
# iterate over the quantization axes
@ -348,9 +369,7 @@ def main():
# 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
)
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)))
@ -377,7 +396,3 @@ def main():
times["end_time"] = timer()
print_output(simulation_parameters, magnetic_entities, pairs, dh, times)
if __name__ == "__main__":
main()

Write Preview
Loading…
Cancel
Save