et209d
/
grogu
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '11f68b9a43'
${ noResults }
grogu/grgoupy_comondor.py

1851 lines
63 KiB
Raw Normal View History Unescape

valami elromlott menet közben!!!
2 months ago
# 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)