almost working

class-solution
Daniel Pozsar 3 months ago
parent b11ab970d6
commit 7631cd38b6

@ -2,12 +2,12 @@
# TODO # TODO
- Definition of magnetic entities: [x] Definition of magnetic entities:
* Through simple sequence o forbitals in the unit cell * Through simple sequence o forbitals in the unit cell
* Through atom specification * Through atom specification
* Through atom and orbital specification * Through atom and orbital specification
- Separation of TR and TRB components of the Hamiltonian, Identification of the exchange field. [x] Separation of TR and TRB components of the Hamiltonian, Identification of the exchange field.
- Definition of commutator expressions, old projection matrix elements [x] Definition of commutator expressions, old projection matrix elements
- Efficient calculation of Green's functions [x] Efficient calculation of Green's functions
- Calculation of energy and momentum resolved derivatives [] Calculation of energy and momentum resolved derivatives
- Parallel BZ and serial energy integral [] Parallel BZ and serial energy integral

@ -23,7 +23,28 @@ from itertools import permutations, product
import numpy as np import numpy as np
from scipy.special import roots_legendre from scipy.special import roots_legendre
# 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]])
# define some useful functions # define some useful functions
def hsk(H, ss, sc_off, k=(0, 0, 0)):
"""
One way to speed up Hk and Sk generation
"""
k = np.asarray(k, np.float64) # this two conversion lines
k.shape = (-1,) # are from the sisl source
# this generates the list of phases
phases = np.exp(-1j * 2 * np.pi * k @ sc_off.T)
HK = np.einsum("abc,a->bc", H, phases)
SK = np.einsum("abc,a->bc", ss, phases)
return HK, SK
def make_contour(emin=-20, emax=0.0, enum=42, p=150): def make_contour(emin=-20, emax=0.0, enum=42, p=150):
@ -80,14 +101,7 @@ def make_kset(dirs="xyz", NUMK=20):
return kset return kset
# Pauli matrices def commutator(a, b):
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]])
def comm(a, b):
"Shorthand for commutator" "Shorthand for commutator"
return a @ b - b @ a return a @ b - b @ a
@ -104,7 +118,7 @@ def tau_u(u):
def crossM(u): def crossM(u):
""" """
Definition for the cross-product matrix. Definition for the cross-product matrix.
Acting as a crossproduct with vector u. Acting as a cross product with vector u.
""" """
return np.array([[0, -u[2], u[1]], [u[2], 0, -u[0]], [-u[1], u[0], 0]]) return np.array([[0, -u[2], u[1]], [u[2], 0, -u[0]], [-u[1], u[0], 0]])
@ -152,6 +166,7 @@ def spin_tracer(M):
M12 = M[0::2, 1::2] M12 = M[0::2, 1::2]
M21 = M[1::2, 0::2] M21 = M[1::2, 0::2]
M22 = M[1::2, 1::2] M22 = M[1::2, 1::2]
M_o = dict() M_o = dict()
M_o["x"] = M12 + M21 M_o["x"] = M12 + M21
M_o["y"] = 1j * (M12 - M21) M_o["y"] = 1j * (M12 - M21)
@ -197,3 +212,15 @@ def blow_up_orbindx(orb_indices):
Function to blow up orbital indeces to make SPIN BOX indices. Function to blow up orbital indeces to make SPIN BOX indices.
""" """
return np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten() return np.array([[2 * o, 2 * o + 1] for o in orb_indices]).flatten()
def calculate_exchange_tensor(pair):
o1, o2, o3 = pair["energies"] # o1=x, o2=y, o3=z
# 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])]),
J_ii = np.array([o2[-1], o3[0], o1[0]]) # xx, yy, zz
J_S = -0.5 * np.array([o3[1] + o3[2], o2[1] + o2[1], o1[1] + o1[2]]) # yz, zx, xy
D = 0.5 * np.array([o1[1] - o1[2], o2[2] - o2[1], o3[1] - o3[2]]) # x, y, z
return J_ii.sum() / 3, D, np.concatenate([J_ii[:2] - J_ii.sum()/3, J_S]).flatten()

@ -2,7 +2,7 @@
"cells": [ "cells": [
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 10, "execution_count": 14,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
@ -11,26 +11,12 @@
"os.environ[\"OPENBLAS_NUM_THREADS\"] = \"4\" # export OPENBLAS_NUM_THREADS=4 \n", "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"4\" # export OPENBLAS_NUM_THREADS=4 \n",
"os.environ[\"MKL_NUM_THREADS\"] = \"4\" # export MKL_NUM_THREADS=6\n", "os.environ[\"MKL_NUM_THREADS\"] = \"4\" # export MKL_NUM_THREADS=6\n",
"os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"4\" # export VECLIB_MAXIMUM_THREADS=4\n", "os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"4\" # export VECLIB_MAXIMUM_THREADS=4\n",
"os.environ[\"NUMEXPR_NUM_THREADS\"] = \"4\" # export NUMEXPR_NUM_THREADS=6" "os.environ[\"NUMEXPR_NUM_THREADS\"] = \"4\""
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.linalg import eigh\n",
"#import matplotlib.pyplot as plt\n",
"import sisl\n",
"import warnings\n",
"from grogu.useful import *"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 12, "execution_count": 15,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [
{ {
@ -39,21 +25,35 @@
"'0.14.3'" "'0.14.3'"
] ]
}, },
"execution_count": 12, "execution_count": 15,
"metadata": {}, "metadata": {},
"output_type": "execute_result" "output_type": "execute_result"
} }
], ],
"source": [ "source": [
"np.set_printoptions(linewidth=2*75)\n", "import numpy as np\n",
"sisl.__version__" "import sisl\n",
"from grogu.useful import *\n",
"from mpi4py import MPI\n",
"from numpy.linalg import inv\n",
"import warnings\n",
"\n",
"sisl.__version__\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 13, "execution_count": 16,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of nodes in the parallel cluster: 1\n"
]
}
],
"source": [ "source": [
"# this cell mimicks an input file\n", "# this cell mimicks an input file\n",
"fdf = sisl.get_sile('/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf')\n", "fdf = sisl.get_sile('/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf')\n",
@ -62,7 +62,7 @@
"# list of reference directions for around which we calculate the derivatives\n", "# list of reference directions for around which we calculate the derivatives\n",
"# o is the quantization axis, v and w are two axes perpendicular to it\n", "# o is the quantization axis, v and w are two axes perpendicular to it\n",
"# at this moment the user has to supply o,v,w on the input. \n", "# at this moment the user has to supply o,v,w on the input. \n",
"# we can have some default sfor this\n", "# we can have some default for this\n",
"ref_xcf_orientations=[dict(o=np.array([1,0,0]),vw=[np.array([0,1,0]),np.array([0,0,1])]),\n", "ref_xcf_orientations=[dict(o=np.array([1,0,0]),vw=[np.array([0,1,0]),np.array([0,0,1])]),\n",
" dict(o=np.array([0,1,0]),vw=[np.array([1,0,0]),np.array([0,0,1])]),\n", " dict(o=np.array([0,1,0]),vw=[np.array([1,0,0]),np.array([0,0,1])]),\n",
" dict(o=np.array([0,0,1]),vw=[np.array([1,0,0]),np.array([0,1,0])]),]\n", " dict(o=np.array([0,0,1]),vw=[np.array([1,0,0]),np.array([0,1,0])]),]\n",
@ -80,60 +80,51 @@
" dict(ai=0,aj=2,Ruc=np.array([-1,0,0])),\n", " dict(ai=0,aj=2,Ruc=np.array([-1,0,0])),\n",
" dict(ai=1,aj=2,Ruc=np.array([-1,0,0]))]\n", " dict(ai=1,aj=2,Ruc=np.array([-1,0,0]))]\n",
"\n", "\n",
"# Bz sampling\n", "# Brilloun zone sampling and Green function contour integral\n",
"nk=np.array([10,10,0])" "kset = 20\n",
"kdirs = \"xy\"\n",
"ebot = -40\n",
"eset = 50\n",
"esetp = 10000\n",
"\n",
"\n",
"# MPI parameters\n",
"comm = MPI.COMM_WORLD\n",
"size = comm.Get_size()\n",
"rank = comm.Get_rank()\n",
"root_node = 0\n",
"if rank == root_node:\n",
" print('Number of nodes in the parallel cluster: ',size)\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 14, "execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# digestion of the input\n", "# digestion of the input\n",
"# read in geometry and hamiltonian\n", "# read in hamiltonian\n",
"geo = fdf.read_geometry()\n",
"dh = fdf.read_hamiltonian()\n", "dh = fdf.read_hamiltonian()\n",
"\n", "\n",
"# unit cell index\n",
"uc_in_sc_idx=dh.lattice.sc_index([0,0,0])\n", "uc_in_sc_idx=dh.lattice.sc_index([0,0,0])\n",
"\n", "\n",
"# get indecies for orbitals of the magnetic entities\n", "\n",
"for i,d in enumerate(magnetic_entities):\n",
" parsed = parse_magnetic_entity(dh,**d)\n",
" magnetic_entities[i]['orbital_indeces'] = parsed\n",
" magnetic_entities[i]['spin_box_indeces'] = blow_up_orbindx(parsed)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"-5.82448514"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# WE WILL NOT NEED THIS!!\n", "# WE WILL NOT NEED THIS!!\n",
"eigfile=sisl.io.siesta.eigSileSiesta('/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.EIG')\n", "eigfile=sisl.io.siesta.eigSileSiesta('/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.EIG')\n",
"EF=eigfile.read_fermi_level()\n", "EF=eigfile.read_fermi_level()\n",
"EF" "\n",
] "\n",
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"NO=dh.no # shorthand for number of orbitals in the unit cell\n", "NO=dh.no # shorthand for number of orbitals in the unit cell\n",
"\n",
"# preprocessing Hamiltonian and overlap matrix elements\n", "# preprocessing Hamiltonian and overlap matrix elements\n",
"h11 = dh.tocsr(dh.M11r)\n", "h11 = dh.tocsr(dh.M11r)\n",
"h11 += dh.tocsr(dh.M11i)*1.0j\n", "h11 += dh.tocsr(dh.M11i)*1.0j\n",
@ -153,20 +144,15 @@
"\n", "\n",
"sov = dh.tocsr(dh.S_idx).toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n", "sov = dh.tocsr(dh.S_idx).toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n",
"\n", "\n",
"# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n",
"\n", "\n",
"# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n",
"U=np.vstack([np.kron(np.eye(NO,dtype=int),[1,0]),np.kron(np.eye(NO,dtype=int),[0,1])])\n", "U=np.vstack([np.kron(np.eye(NO,dtype=int),[1,0]),np.kron(np.eye(NO,dtype=int),[0,1])])\n",
"# This is the permutation that transforms ud1ud2 to u12d12\n", "# This is the permutation that transforms ud1ud2 to u12d12\n",
"# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n", "# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n",
"\n",
"# the inverse transformation is U.T u12d12 to ud1ud2\n", "# the inverse transformation is U.T u12d12 to ud1ud2\n",
"# That is FROM ORBITAL BOX to SPIN BOX => U.T\n", "# That is FROM ORBITAL BOX to SPIN BOX => U.T\n",
"\n", "\n",
"# this is a test\n", "# From now on everything is in SPIN BOX!!\n",
"# u12d12=np.array([np.kron([1,1],np.arange(NO)),np.kron([0,1],np.ones(NO))],dtype=int).T\n",
"# ud1ud2=np.array([np.kron(np.arange(NO),[1,1]),np.kron(np.ones(NO),[0,1])],dtype=int).T\n",
"# np.allclose(u12d12,U@ud1ud2)\n",
"\n",
"hh,ss = \\\n", "hh,ss = \\\n",
"np.array([U.T@np.block([[h11[:,:,i],h12[:,:,i]],\n", "np.array([U.T@np.block([[h11[:,:,i],h12[:,:,i]],\n",
" [h21[:,:,i],h22[:,:,i]]])@U for i in range(dh.lattice.nsc.prod())]), \\\n", " [h21[:,:,i],h22[:,:,i]]])@U for i in range(dh.lattice.nsc.prod())]), \\\n",
@ -189,516 +175,341 @@
"hTRB = (hh-hTR)/2\n", "hTRB = (hh-hTR)/2\n",
"\n", "\n",
"# extracting the exchange field\n", "# extracting the exchange field\n",
"traced=[spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())]\n", "traced=[spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77\n",
"\n",
"XCF = np.array([np.array([f['x'] for f in traced]),\n", "XCF = np.array([np.array([f['x'] for f in traced]),\n",
" np.array([f['y'] for f in traced]),\n", " np.array([f['y'] for f in traced]),\n",
" np.array([f['z'] for f in traced])])\n", " np.array([f['z'] for f in traced])]) # equation 77\n",
"\n",
"# Check if exchange field has scalar part\n", "# Check if exchange field has scalar part\n",
"max_xcfs=abs(np.array(np.array([f['c'] for f in traced]))).max()\n", "max_xcfs=abs(np.array(np.array([f['c'] for f in traced]))).max()\n",
"if max_xcfs > 1e-12:\n", "if max_xcfs > 1e-12:\n",
" warnings.warn(f\"Exchange field has non negligeble scalr part. Largest value is {max_xcfs}\")" " warnings.warn(f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\")\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 17, "execution_count": 18,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"# collecting information for all proposed rotations, this might generate a lot of unnceessary data\n", "# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n",
"# information being collected: \n", "for i, mag_ent in enumerate(magnetic_entities):\n",
"# - matrices rotating the direction of the exchange field to the proposed reference direction\n", " parsed = parse_magnetic_entity(dh,**mag_ent) # parse orbital indexes\n",
"# - rotated xc field \n", " magnetic_entities[i]['orbital_indeces'] = parsed\n",
"# - single commutators of the unit cell localized exchange field for the two variations perpendicular to the reference direction\n", " magnetic_entities[i]['spin_box_indeces'] = blow_up_orbindx(parsed) # calculate spin box indexes\n",
"# - double commutators of the unit cell localized exchange field for the two variations perpendicular to the reference direction\n", " spin_box_shape = len(mag_ent[\"spin_box_indeces\"]) # calculate size for Greens function generation\n",
"\n",
"\n",
"# Transformation matrices for rotating the exchange field\n",
"for i,orient in enumerate(ref_xcf_orientations):\n",
"\n", "\n",
" # obtain rotated exchange field \n", " mag_ent['energies'] = [] # we will store the second order energy derivations here\n",
" R=RotMa2b(scf_xcf_orientation,orient['o'])\n",
" \n", " \n",
" rot_XCF = np.einsum('ij,jklm->iklm',R, XCF)\n", " mag_ent['Gii'] = [] # Greens function\n",
" rot_H_XCF = sum([np.kron(rot_XCF[i],tau) for i,tau in enumerate([tau_x,tau_y,tau_z])])\n", " mag_ent['Gii_tmp'] = [] # Greens function for parallelization\n",
" rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx]\n", " mag_ent[\"Vu1\"] = [ list([]) for _ in range(len(ref_xcf_orientations))] # These will be the perturbed potentials from eq. 100\n",
" \n", " mag_ent[\"Vu2\"] = [ list([]) for _ in range(len(ref_xcf_orientations))]\n",
" # obtain total Hamiltonian with the rotated exchange field\n", " for i in ref_xcf_orientations:\n",
" rot_H = hTRS+rot_H_XCF\n", " mag_ent['Gii'].append(np.zeros((eset, spin_box_shape, spin_box_shape),dtype='complex128')) # Greens functions for every quantization axis\n",
" mag_ent['Gii_tmp'].append(np.zeros((eset, spin_box_shape, spin_box_shape),dtype='complex128'))\n",
"\n",
"# 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\n",
"for pair in pairs:\n",
" spin_box_shape_i, spin_box_shape_j = len(magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]), len(magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]) # calculate size for Greens function generation\n",
" \n", " \n",
" # Update magnetic entities with commutator data for each \n", " pair['energies'] = [] # we will store the second order energy derivations here\n",
" for mag_ent in magnetic_entities:\n",
" mag_ent['Commutator_Data'] = []\n",
"\n", "\n",
" for u in orient['vw']:\n", " pair['Gij'] = [] # Greens function\n",
" \n", " pair['Gji'] = []\n",
" Tu = np.kron(np.eye(NO,dtype=int),tau_u(u))\n", " pair['Gij_tmp'] = [] # Greens function for parallelization\n",
" pair['Gji_tmp'] = []\n",
"\n", "\n",
" Vu1 = 1j/2*comm(rot_H_XCF_uc,Tu)\n", " pair[\"Vij\"] = [list([]) for _ in range(len(ref_xcf_orientations))] # These will be the perturbed potentials from eq. 100\n",
" Vu2 = 1/8*comm(comm(Tu,rot_H_XCF_uc),Tu)\n", " pair[\"Vji\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n",
"\n", "\n",
" for mag_ent in magnetic_entities:\n", " for i in ref_xcf_orientations: \n",
" idx = mag_ent['spin_box_indeces']\n", " pair['Gij'].append(np.zeros((eset,spin_box_shape_i,spin_box_shape_j),dtype='complex128'))\n",
" mag_ent['Commutator_Data'].append(dict(Vu1=Vu1[idx,:][:,idx],\n", " pair['Gij_tmp'].append(np.zeros((eset,spin_box_shape_i,spin_box_shape_j),dtype='complex128')) # Greens functions for every quantization axis\n",
" Vu2=Vu2[idx,:][:,idx]))\n", " pair['Gji'].append(np.zeros((eset,spin_box_shape_j,spin_box_shape_i),dtype='complex128'))\n",
"\n", " pair['Gji_tmp'].append(np.zeros((eset,spin_box_shape_j,spin_box_shape_i),dtype='complex128'))\n"
"\n",
"# TODO STORE ABOVE INFO!!!!!"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 19,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"k = np.array([0.11441,0.234432,0.0])\n", "kset = make_kset(dirs=kdirs,NUMK=kset) # generate k space sampling\n",
"\n" "wkset = np.ones(len(kset)) / len(kset) # generate weights for k points\n",
"kpcs = np.array_split(kset,size) # split the k points based on MPI size\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": 19, "execution_count": 20,
"metadata": {}, "metadata": {},
"outputs": [ "outputs": [],
{
"ename": "NameError",
"evalue": "name 'magn_ent_idices' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[19], line 20\u001b[0m\n\u001b[1;32m 17\u001b[0m GK\u001b[38;5;241m=\u001b[39mvecs\u001b[38;5;129m@np\u001b[39m\u001b[38;5;241m.\u001b[39mdiag(\u001b[38;5;241m1\u001b[39m\u001b[38;5;241m/\u001b[39m(ze\u001b[38;5;241m-\u001b[39mvals))\u001b[38;5;129m@vecs\u001b[39m\u001b[38;5;241m.\u001b[39mconj()\u001b[38;5;241m.\u001b[39mT\n\u001b[1;32m 19\u001b[0m \u001b[38;5;66;03m#phase=np.exp(1j*np.dot(np.dot(k,dh.rcell),offset))\u001b[39;00m\n\u001b[0;32m---> 20\u001b[0m Gii\u001b[38;5;241m=\u001b[39mGK[\u001b[43mmagn_ent_idices\u001b[49m[ai],:][:,magn_ent_idices[ai]]\n\u001b[1;32m 21\u001b[0m Gjj\u001b[38;5;241m=\u001b[39mGK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]\n\u001b[1;32m 22\u001b[0m Gij\u001b[38;5;241m=\u001b[39mGK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]\u001b[38;5;241m*\u001b[39mnp\u001b[38;5;241m.\u001b[39mexp(\u001b[38;5;241m-\u001b[39mk\u001b[38;5;129m@offset\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m*\u001b[39mnp\u001b[38;5;241m.\u001b[39mpi)\n",
"\u001b[0;31mNameError\u001b[0m: name 'magn_ent_idices' is not defined"
]
}
],
"source": [ "source": [
"#This goes inside the for loop for k \n", "# collecting information for all proposed rotations, this might generate a lot of unnecessary data\n",
"k = np.array([0,0.1,0.0])\n", "# information being collected: \n",
"# - matrices rotating the direction of the exchange field to the proposed reference direction\n",
"# - rotated xc field \n",
"# - single commutators of the unit cell localized exchange field for the two variations perpendicular to the reference direction\n",
"# - double commutators of the unit cell localized exchange field for the two variations perpendicular to the reference direction\n",
"\n", "\n",
"# this will contain all the data needed to calculate the energy variations upon rotation\n",
"hamiltonians = []\n",
"\n", "\n",
"phases = np.exp(-2j*np.pi*k@dh.lattice.sc_off.T)\n", "# iterate over the reference directions (quantization axes)\n",
"for i,orient in enumerate(ref_xcf_orientations):\n",
"\n", "\n",
"np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell),dh.lattice.sc_off.T))\n", " # obtain rotated exchange field \n",
"ss.shape,rot_H_XCF.shape,hTRS.shape,phases.shape\n", " R=RotMa2b(scf_xcf_orientation,orient['o'])\n",
" rot_XCF = np.einsum('ij,jklm->iklm',R, XCF)\n",
" rot_H_XCF = sum([np.kron(rot_XCF[i],tau) for i,tau in enumerate([tau_x,tau_y,tau_z])])\n",
" rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx]\n",
" \n",
" # obtain total Hamiltonian with the rotated exchange field\n",
" rot_H = hTRS+rot_H_XCF # equation 76\n",
"\n", "\n",
"HK=np.einsum('abc,a->bc',hTRS+rot_H_XCF,phases)\n", " hamiltonians.append(dict(orient=orient['o'], H=rot_H, rotations=[])) # store orientation and rotated Hamiltonian\n",
"SK=np.einsum('abc,a->bc',ss,phases)\n", " \n",
" for u in orient['vw']: # these are the infinitezimal rotations (for now) perpendicular to the quantization axis\n",
" Tu = np.kron(np.eye(NO,dtype=int),tau_u(u)) # section 2.H\n",
"\n", "\n",
"vals,vecs=eigh(HK,SK)\n", " Vu1 = 1j/2*commutator(rot_H_XCF_uc,Tu) # equation 100\n",
" Vu2 = 1/8*commutator(commutator(Tu,rot_H_XCF_uc),Tu) # equation 100\n",
"\n", "\n",
"#This goes inside the for loop for ze\n", " for mag_ent in magnetic_entities:\n",
"ze = -.1+0.000j\n", " 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\n",
"GK=vecs@np.diag(1/(ze-vals))@vecs.conj().T\n", " mag_ent[\"Vu2\"][i].append(Vu2[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :])\n",
"\n", "\n",
"#phase=np.exp(1j*np.dot(np.dot(k,dh.rcell),offset))\n", " for pair in pairs:\n",
"Gii=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]\n", " ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"] # get the pair orbital sizes from the magnetic entities\n",
"Gjj=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]\n", " aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n",
"Gij=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]*np.exp(-k@offset*2*np.pi)\n", " pair[\"Vij\"][i].append(Vu1[:, ai][aj, :]) # fill up the perturbed potentials (for now) based on the on-site projections\n",
"Gji=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]*np.exp( k@offset*2*np.pi)" " pair[\"Vji\"][i].append(Vu1[:, aj][ai, :])\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 21,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [
"source": [
"np.allclose(np.linalg.inv(ze*SK-HK),GK), abs(np.linalg.inv(ze*SK-HK)-GK).max()"
]
},
{ {
"cell_type": "code", "name": "stdout",
"execution_count": null, "output_type": "stream",
"metadata": {}, "text": [
"outputs": [], "Number of magnetic entities being calculated: 4\n",
"source": [ "We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.\n",
"def hsk(dh,k=(0,0,0)):\n", "The shape of the Hamiltonian and the Greens function is 84x84.\n"
" '''\n",
" One way to speed up Hk and Sk generation\n",
" '''\n",
" k = np.asarray(k, np.float64) # this two conversion lines\n",
" k.shape = (-1,) # are from the sisl source\n",
"\n",
" # this generates the list of phases\n",
" phases = np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell),\n",
" dh.lattice.sc_off.T))\n",
"\n",
" HK11 = np.einsum('abc,c->ab',dh.h11,phases)\n",
" HK12 = np.einsum('abc,c->ab',dh.h12,phases)\n",
" HK21 = np.einsum('abc,c->ab',dh.h21,phases)\n",
" HK22 = np.einsum('abc,c->ab',dh.h22,phases)\n",
"\n",
" SK = np.einsum('abc,c->ab',dh.sov,phases)\n",
"\n",
" return HK11,HK12,HK21,HK22,SK\n"
] ]
}, }
{ ],
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [ "source": [
"if rank == root_node:\n",
" print('Number of magnetic entities being calculated: ',len(magnetic_entities))\n",
" print('We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.')\n",
" print(f'The shape of the Hamiltonian and the Greens function is {NO}x{NO}.')\n",
"comm.Barrier()\n",
"#----------------------------------------------------------------------\n",
"\n",
"# make energy contour \n",
"# we are working in eV now !\n",
"# and sisil shifts E_F to 0 !\n",
"cont = make_contour(emin=ebot,enum=eset,p=esetp)\n",
"eran = cont.ze\n",
"\n",
"#---------------------------------------------------------------------- \n",
"# sampling the integrand on the contour and the BZ\n",
"for k in kpcs[rank]:\n",
" wk = wkset[rank] # weight of k point in BZ integral\n",
" for i, hamiltonian_orientation in enumerate(hamiltonians): # iterate over reference directions\n",
" # calculate Greens function\n",
" H = hamiltonian_orientation[\"H\"]\n",
" HK,SK = hsk(H, ss, dh.sc_off, k)\n",
" Gk = inv(SK * eran.reshape(eset, 1, 1) - HK)\n",
" \n", " \n",
" # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
" for mag_ent in magnetic_entities:\n",
" mag_ent[\"Gii_tmp\"][i] += Gk[:,mag_ent[\"spin_box_indeces\"]][...,mag_ent[\"spin_box_indeces\"]] * wk\n",
"\n", "\n",
"rot_xcf_mats=[]\n", " for pair in pairs:\n",
"for i,rot in enumerate(rots):\n", " # add phase shift based on the cell difference\n",
" R=RotM(*rot) # rotation matrix to go from scf to reference direction\n", " phase=np.exp(1j * 2 * np.pi * k @ pair[\"Ruc\"].T)\n",
" R@scf_xcf_orientation # exchange field orientation \n", " \n",
"\n", " # get the pair orbital sizes from the magnetic entities\n",
"\n", " ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]\n",
"rot_xcf_mats = [RotM(*rot) for rot in rots] \n", " aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n",
"# directions to ratate around \n",
"dirs_xcf = [R@scf_xcf_orientation for R in rot_xcf_mats]\n",
"dirs_deriv=[]\n",
"for o in dirs_xcf:\n",
" x,y,z = np.eye(3)\n",
" R = RotMa2b(z,o)\n",
" dirs_deriv.append([R@x,R@y])\n",
"\n",
"\n",
"# geting the exchange field part of the reference Hamiltonians for the rotated exchange field\n",
"dh.hxc_rots=[]\n",
"for mtx in rot_xcf_mats:\n",
" rotated_xcf = np.einsum('ij,jklm->iklm',mtx, XCF)\n",
" dh.hxc_rots.append(sum([np.kron(rotated_xcf[i],tau) \n",
" for i,tau in enumerate([tau_x,tau_y,tau_z]) ]))\n",
"\n",
"\n", "\n",
" # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
" pair['Gij_tmp'][i] += Gk[:,ai][..., aj] * phase * wk\n",
" pair['Gji_tmp'][i] += Gk[:,aj][..., ai] * phase * wk\n",
"\n", "\n",
"cummutators=[]\n", "# summ reduce partial results of mpi nodes\n",
"for i in range(len(hamiltonians)):\n",
" for mag_ent in magnetic_entities:\n",
" comm.Reduce(mag_ent[\"Gii_tmp\"][i], mag_ent[\"Gii\"][i], root=root_node)\n",
"\n", "\n",
"\n" " for pair in pairs:\n",
" comm.Reduce(pair[\"Gij_tmp\"][i], pair[\"Gij\"][i], root=root_node)\n",
" comm.Reduce(pair[\"Gji_tmp\"][i], pair[\"Gji\"][i], root=root_node)\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 22,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"i=0\n", "# split magnetic entities and pairs over MPI nodes\n",
"mag_ent_parallel_set = np.array_split(magnetic_entities,size)\n",
"pairs_parallel_set = np.array_split(pairs,size)\n",
"\n", "\n",
"dh.hxc_rots[i][uc_in_sc_idx]" "# iterate over the magnetic entities\n",
] "for tracker, mag_ent in enumerate(mag_ent_parallel_set[rank]):\n",
}, " # iterate over the quantization axes\n",
{ " for i, Gii in enumerate(mag_ent[\"Gii\"]):\n",
"cell_type": "code", " storage = []\n",
"execution_count": null, " # iterate over the first and second order local perturbations\n",
"metadata": {}, " for Vu1, Vu2 in zip(mag_ent[\"Vu1\"][i], mag_ent[\"Vu2\"][i]):\n",
"outputs": [], " # The Szunyogh-Lichtenstein formula\n",
"source": [] " traced = np.trace((Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii),axis1=1,axis2=2)\n",
}, " # evaluation of the contour integral\n",
{ " storage.append(np.trapz(-1/np.pi * np.imag(traced * cont.we)))\n",
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dh.hTRS.shape"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"np.kron(np.random.rand(3,3),tau_x)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dh.hh"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Sanity check for the exchange field orientations\n",
"xcf_test=np.array([[abs(xcfield[i].c).max(),\n",
" abs(xcfield[i].x).max(),\n",
" abs(xcfield[i].y).max(),\n",
" abs(xcfield[i].z).max()] for i in range(dh.lattice.nsc.prod())])\n",
"# Check how the exchange field components behave\n",
"plt.plot(np.linalg.norm(dh.lattice.sc_off,axis=1),xcf_test[:,0],'o',label='scalar')\n",
"plt.plot(np.linalg.norm(dh.lattice.sc_off,axis=1),xcf_test[:,1],'o',label='x')\n",
"plt.plot(np.linalg.norm(dh.lattice.sc_off,axis=1),xcf_test[:,2],'o',label='y')\n",
"plt.plot(np.linalg.norm(dh.lattice.sc_off,axis=1),xcf_test[:,3],'o',label='z')\n",
"plt.yscale('log')\n",
"plt.ylim(1e-22,1e1)\n",
"plt.legend()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"i=41\n",
"DAT=np.vstack([np.hstack([dh.h11[:,:,i]+EF*dh.sov[:,:,i],dh.h12[:,:,i]]),\n",
" np.hstack([dh.h21[:,:,i], dh.h22[:,:,i]]+EF*dh.sov[:,:,i])])\n",
"dh.lattice.sc_off[i]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"for i,v in enumerate(dh.lattice.sc_off):\n",
" if np.allclose(v,np.array([1,0,0])):\n",
" print(i)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.pcolor((U.T@DAT@U)[:6,:6].imag)\n",
"plt.colorbar()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dh.h11 = dh.tocsr(dh.M11r)\n",
"dh.h11 += dh.tocsr(dh.M11i)*1.0j\n",
"dh.h11 = dh.h11.toarray().astype('complex128')\n",
"dh.h22 = dh.tocsr(dh.M22r)\n",
"dh.h22 += dh.tocsr(dh.M22i)*1.0j\n",
"dh.h22 = dh.h22.toarray().astype('complex128')\n",
"\n",
"dh.h12 = dh.tocsr(dh.M12r)\n",
"dh.h12 += dh.tocsr(dh.M12i)*1.0j\n",
"dh.h12 = dh.h12.toarray().astype('complex128')\n",
"dh.h21 = dh.tocsr(dh.M21r)\n",
"dh.h21 += dh.tocsr(dh.M21i)*1.0j\n",
"dh.h21 = dh.h21.toarray().reshape(dh.no,dh.n_s,dh.no).transpose(0,2,1).astype('complex128')\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dh.tocsr(dh.S_idx"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def hsk(dh,k=(0,0,0)):\n",
" '''\n",
" One way to speed up Hk and Sk generation\n",
" '''\n",
" k = np.asarray(k, np.float64) # this two conversion lines\n",
" k.shape = (-1,) # are from the sisl source\n",
"\n", "\n",
" # this generates the list of phases\n", " # fill up the magnetic entities dictionary with the energies\n",
" phases = np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell),\n", " idx = np.array([len(mag_ent_parallel_set[i]) for i in range(rank)]).sum() + tracker\n",
" dh.lattice.sc_off.T))\n", " magnetic_entities[int(idx)][\"energies\"].append(storage)\n",
"\n", "\n",
" HK11 = np.einsum('abc,c->ab',dh.h11,phases)\n", "# iterate over the pairs\n",
" HK12 = np.einsum('abc,c->ab',dh.h12,phases)\n", "for tracker, pair in enumerate(pairs_parallel_set[rank]):\n",
" HK21 = np.einsum('abc,c->ab',dh.h21,phases)\n", " # iterate over the quantization axes\n",
" HK22 = np.einsum('abc,c->ab',dh.h22,phases)\n", " for i, (Gij, Gji) in enumerate(zip(pair[\"Gij\"], pair[\"Gji\"])):\n",
" site_i = magnetic_entities[pair[\"ai\"]]\n",
" site_j = magnetic_entities[pair[\"aj\"]]\n",
"\n", "\n",
" SK = np.einsum('abc,c->ab',dh.sov,phases)\n", " storage = []\n",
" # iterate over the first order local perturbations in all possible orientations for the two sites\n",
" for Vui in site_i[\"Vu1\"][i]:\n",
" for Vuj in site_j[\"Vu1\"][i]:\n",
" # The Szunyogh-Lichtenstein formula\n",
" traced = np.trace((Vui @ Gij @ Vuj @ Gji),axis1=1,axis2=2)\n",
" # evaluation of the contour integral\n",
" storage.append(np.trapz(-1/np.pi * np.imag(traced * cont.we)))\n",
"\n", "\n",
" return HK11,HK12,HK21,HK22,SK\n" " # fill up the pairs dictionary with the energies\n",
" idx = np.array([len(pairs_parallel_set[i]) for i in range(rank)]).sum() + tracker\n",
" pairs[int(idx)][\"energies\"].append(storage)\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 23,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [],
"source": [ "source": [
"myHk=(lambda x:np.vstack([np.hstack([x[0],x[1]]),\n", "def calculate_exchange_tensor(pair):\n",
" np.hstack([x[2],x[3]])]))(hsk(dh))\n", " Eo, Ev, Ew = pair[\"energies\"]\n",
"mySk=(lambda x:np.vstack([np.hstack([x[-1],0*x[-1]]),\n", " J = np.array([Ew[-1], Ev[-1], Eo[0]]) # xx, yy, zz\n",
" np.hstack([0*x[-1],x[-1]])]))(hsk(dh))\n", " JS = -0.5 * np.array([Eo[1] + Eo[2], Ev[1] + Ev[2], Ew[1] + Ew[2]]) # yz, zx, xy\n",
"\n", " D = 0.5 * np.array([Eo[1] - Eo[2], Ev[2] - Ev[1], Ew[1] - Ew[2]]) # x, y, z\n",
"myHkshift=myHk+EF*mySk" " return J.sum()/3 * 1000"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 24,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [
"source": [
"H0=np.vstack( [np.hstack([dh.h11[:,:,0],dh.h12[:,:,0]]),\n",
" np.hstack([dh.h21[:,:,0],dh.h22[:,:,0]])])\n",
"\n",
"H0shifted=np.vstack( [np.hstack([dh.h11[:,:,0]+EF*dh.sov[:,:,0],dh.h12[:,:,0]]),\n",
" np.hstack([dh.h21[:,:,0], dh.h22[:,:,0]+EF*dh.sov[:,:,0]])])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"(U.T@H0shifted@U)[:3,:3]"
]
},
{ {
"cell_type": "code", "name": "stdout",
"execution_count": null, "output_type": "stream",
"metadata": {}, "text": [
"outputs": [], "-60.72359053548613\n",
"source": [ "-60.531975234450115\n",
"plt.pcolor((U.T@H0shifted@U)[:5,:5])\n", "-60.524676226428824\n",
"plt.colorbar()" "-6.55042989834691\n",
"-6.047933492978864\n"
] ]
}, }
{ ],
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [ "source": [
"np.linalg.eigvalsh(myHk+EF*mySk)" "print(calculate_exchange_tensor(pairs[0])) # isotropic should be -82 meV\n",
"print(calculate_exchange_tensor(pairs[1])) # these should all be around -41.9 in the isotropic part\n",
"print(calculate_exchange_tensor(pairs[2])) # these should all be around -41.9 in the isotropic part\n",
"print(calculate_exchange_tensor(pairs[3])) # these should all be around -41.9 in the isotropic part\n",
"print(calculate_exchange_tensor(pairs[4])) # these should all be around -41.9 in the isotropic part\n"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 25,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [
"source": [
"dh.Hk().toarray()[:3,:3]"
]
},
{ {
"cell_type": "code", "data": {
"execution_count": null, "text/plain": [
"metadata": {}, "-6.043716409664797"
"outputs": [],
"source": [
"plt.matshow(np.abs(dh.Hk().toarray()-myHk))\n",
"plt.colorbar()"
] ]
}, },
{ "execution_count": 25,
"cell_type": "code",
"execution_count": null,
"metadata": {}, "metadata": {},
"outputs": [], "output_type": "execute_result"
}
],
"source": [ "source": [
"plt.plot(np.linalg.eigvalsh(dh.Hk().toarray())-np.linalg.eigvalsh(myHk))" "-61.33097171216109\n",
"-60.52198328932686\n",
"-60.51657719027764\n",
"-6.545208546361317\n",
"-6.043716409664797"
] ]
}, },
{ {
"cell_type": "code", "cell_type": "code",
"execution_count": null, "execution_count": 26,
"metadata": {}, "metadata": {},
"outputs": [], "outputs": [
"source": []
},
{ {
"cell_type": "code", "ename": "SyntaxError",
"execution_count": null, "evalue": "invalid syntax (1876172784.py, line 5)",
"metadata": {}, "output_type": "error",
"outputs": [], "traceback": [
"source": [ "\u001b[0;36m Cell \u001b[0;32mIn[26], line 5\u001b[0;36m\u001b[0m\n\u001b[0;31m ========================================\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
"magnetic_entities = [parse_magnetic_entity(dh,atom=[0,1],l=2),parse_magnetic_entity(dh,atom=2,l=2)]"
] ]
}, }
{ ],
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [ "source": [
"# TODO\n", "# symmetrizing Hamiltonian and overlap matrix to make them hermitian \n",
"\n", "# Check if exchange field has scalar part\n",
"- Definition of magnetic entities:\n", "# parallel over integrals\n",
" * Through simple sequence o forbitals in the unit cell\n", "\n",
" * Through atom specification\n", "========================================\n",
" * Through atom and orbital specification\n", " \n",
"- Separation of TR and TRB components of the Hamiltonian, Identification of the exchange field. \n", "Atom Angstrom\n",
"- Definition of commutator expressions, old projection matrix elements\n", "# Label, x y z Sx Sy Sz #Q Lx Ly Lz Jx Jy Jz\n",
"- Efficient calculation of Green's functions\n", "--------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"- Calculation of energy and momentum resolved derivatives\n", "Te1 1.8955 1.0943 13.1698 -0.0000 0.0000 -0.1543 # 5.9345 -0.0000 0.0000 -0.0537 -0.0000 0.0000 -0.2080 \n",
"- Parallel BZ and serial energy integral" "Te2 1.8955 1.0943 7.4002 0.0000 -0.0000 -0.1543 # 5.9345 0.0000 -0.0000 -0.0537 0.0000 -0.0000 -0.2080 \n",
"Ge3 -0.0000 2.1887 10.2850 0.0000 0.0000 -0.1605 # 3.1927 -0.0000 0.0000 0.0012 0.0000 0.0000 -0.1593 \n",
"Fe4 -0.0000 0.0000 11.6576 0.0001 -0.0001 2.0466 # 8.3044 0.0000 -0.0000 0.1606 0.0001 -0.0001 2.2072 \n",
"Fe5 -0.0000 0.0000 8.9124 -0.0001 0.0001 2.0466 # 8.3044 -0.0000 0.0000 0.1606 -0.0001 0.0001 2.2072 \n",
"Fe6 1.8955 1.0944 10.2850 0.0000 0.0000 1.5824 # 8.3296 -0.0000 -0.0000 0.0520 -0.0000 0.0000 1.6344 \n",
"==================================================================================================================================\n",
" \n",
"Exchange meV\n",
"--------------------------------------------------------------------------------\n",
"# at1 at2 i j k # d (Ang)\n",
"--------------------------------------------------------------------------------\n",
"Fe4 Fe5 0 0 0 # 2.7452\n",
"Isotropic -82.0854\n",
"DMI 0.12557 -0.00082199 6.9668e-08\n",
"Symmetric-anisotropy -0.60237 -0.83842 -0.00032278 -1.2166e-05 -3.3923e-05\n",
"--------------------------------------------------------------------------------\n",
"Fe4 Fe6 0 0 0 # 2.5835\n",
"Isotropic -41.9627\n",
"DMI 1.1205 -1.9532 0.0018386\n",
"Symmetric-anisotropy 0.26007 -0.00013243 0.12977 -0.069979 -0.042066\n",
"--------------------------------------------------------------------------------\n"
] ]
},
{
"cell_type": "markdown",
"metadata": {},
"source": []
} }
], ],
"metadata": { "metadata": {

@ -0,0 +1,451 @@
import os
from sys import stdout
from tqdm import tqdm
from timeit import default_timer as timer
os.environ["OMP_NUM_THREADS"] = "1" # export OMP_NUM_THREADS=4
os.environ["OPENBLAS_NUM_THREADS"] = "1" # export OPENBLAS_NUM_THREADS=4
os.environ["MKL_NUM_THREADS"] = "1" # export MKL_NUM_THREADS=6
os.environ["VECLIB_MAXIMUM_THREADS"] = "1" # export VECLIB_MAXIMUM_THREADS=4
os.environ["NUMEXPR_NUM_THREADS"] = "1" # export NUMEXPR_NUM_THREADS=6
import numpy as np
import sisl
from grogu.useful import *
from mpi4py import MPI
from numpy.linalg import inv
import warnings
start_time = timer()
# this cell mimicks an input file
fdf = sisl.get_sile(
"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/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 = [
dict(atom=3, l=2),
dict(atom=4, l=2),
dict(atom=5, l=2),
# dict(atom=[3, 4]),
]
# pair information
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])),
]
# 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(path="Not yet specified.",
scf_xcf_orientation=scf_xcf_orientation,
ref_xcf_orientations=ref_xcf_orientations,
kset=kset,
kdirs=kdirs,
ebot=ebot,
eset=eset,
esetp=esetp,
parallel_size=size)
# digestion of the input
# read in hamiltonian
dh = fdf.read_hamiltonian()
try:
simulation_parameters["geom"] = fdf.read_geometry()
except:
print("Error reading geometry.")
# unit cell index
uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])
setup_time = timer()
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())
]
)
# 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"] for f in traced]),
np.array([f["y"] for f in traced]),
np.array([f["z"] for f in traced]),
]
) # 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:
warnings.warn(
f"Exchange field has non negligible scalar part. Largest value is {max_xcfs}"
)
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):
parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes
magnetic_entities[i]["orbital_indeces"] = parsed
magnetic_entities[i]["spin_box_indeces"] = blow_up_orbindx(
parsed
) # calculate spin box indexes
spin_box_shape = len(
mag_ent["spin_box_indeces"]
) # calculate size for Greens function generation
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
mag_ent["Vu1"] = [
list([]) for _ in range(len(ref_xcf_orientations))
] # These will be the perturbed potentials from eq. 100
mag_ent["Vu2"] = [list([]) for _ in range(len(ref_xcf_orientations))]
for i in ref_xcf_orientations:
mag_ent["Gii"].append(
np.zeros((eset, spin_box_shape, spin_box_shape), dtype="complex128")
) # Greens functions for every quantization axis
mag_ent["Gii_tmp"].append(
np.zeros((eset, spin_box_shape, spin_box_shape), dtype="complex128")
)
# for every site we have to store 2x3 Greens function (and the associated _tmp-s)
# in the 3 reference directions, because G_ij and G_ji are both needed
for pair in pairs:
spin_box_shape_i, spin_box_shape_j = len(
magnetic_entities[pair["ai"]]["spin_box_indeces"]
), len(
magnetic_entities[pair["aj"]]["spin_box_indeces"]
) # calculate size for Greens function generation
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"] = []
pair["Vij"] = [
list([]) for _ in range(len(ref_xcf_orientations))
] # These will be the perturbed potentials from eq. 100
pair["Vji"] = [list([]) for _ in range(len(ref_xcf_orientations))]
for i in ref_xcf_orientations:
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")
)
pair["Gji_tmp"].append(
np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype="complex128")
)
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)
k_set_time = timer()
# this will contain all the data needed to calculate the energy variations upon rotation
hamiltonians = []
# iterate over the reference directions (quantization axes)
for i, orient in enumerate(ref_xcf_orientations):
# obtain rotated exchange field
R = RotMa2b(scf_xcf_orientation, orient["o"])
rot_XCF = np.einsum("ij,jklm->iklm", R, XCF)
rot_H_XCF = sum(
[np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])]
)
rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx]
# obtain total Hamiltonian with the rotated exchange field
rot_H = hTRS + rot_H_XCF # equation 76
hamiltonians.append(
dict(orient=orient["o"], H=rot_H, rotations=[])
) # store orientation and rotated Hamiltonian
for u in orient[
"vw"
]: # these are the infinitezimal rotations (for now) perpendicular to the quantization axis
Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H
Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100
Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100
for mag_ent in magnetic_entities:
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"], :]
)
for pair in pairs:
ai = magnetic_entities[pair["ai"]][
"spin_box_indeces"
] # get the pair orbital sizes from the magnetic entities
aj = magnetic_entities[pair["aj"]]["spin_box_indeces"]
pair["Vij"][i].append(
Vu1[:, ai][aj, :]
) # fill up the perturbed potentials (for now) based on the on-site projections
pair["Vji"][i].append(Vu1[:, aj][ai, :])
reference_rotations_time = timer()
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]:
wk = wkset[rank] # weight of k point in BZ integral
for i, hamiltonian_orientation in enumerate(
hamiltonians
): # iterate over reference directions
# calculate Greens function
H = hamiltonian_orientation["H"]
HK, SK = hsk(H, ss, dh.sc_off, k)
Gk = inv(SK * eran.reshape(eset, 1, 1) - HK)
# store the Greens function slice of the magnetic entities (for now) based on the on-site projections
for mag_ent in magnetic_entities:
mag_ent["Gii_tmp"][i] += (
Gk[:, mag_ent["spin_box_indeces"]][..., mag_ent["spin_box_indeces"]]
* wk
)
for pair in pairs:
# add phase shift based on the cell difference
phase = np.exp(1j * 2 * np.pi * k @ pair["Ruc"].T)
# get the pair orbital sizes from the magnetic entities
ai = magnetic_entities[pair["ai"]]["spin_box_indeces"]
aj = magnetic_entities[pair["aj"]]["spin_box_indeces"]
# store the Greens function slice of the magnetic entities (for now) based on the on-site projections
pair["Gij_tmp"][i] += Gk[:, ai][..., aj] * phase * wk
pair["Gji_tmp"][i] += Gk[:, aj][..., ai] * phase * wk
# summ reduce partial results of mpi nodes
for i in range(len(hamiltonians)):
for mag_ent in magnetic_entities:
comm.Reduce(mag_ent["Gii_tmp"][i], mag_ent["Gii"][i], root=root_node)
for pair in pairs:
comm.Reduce(pair["Gij_tmp"][i], pair["Gij"][i], root=root_node)
comm.Reduce(pair["Gji_tmp"][i], pair["Gji"][i], root=root_node)
green_function_inversion_time = timer()
if rank == root_node:
# iterate over the magnetic entities
for tracker, mag_ent in enumerate(magnetic_entities):
# iterate over the quantization axes
for i, Gii in enumerate(mag_ent["Gii"]):
storage = []
# iterate over the first and second order local perturbations
for Vu1, Vu2 in zip(mag_ent["Vu1"][i], mag_ent["Vu2"][i]):
# The Szunyogh-Lichtenstein formula
traced = np.trace((Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2)
# evaluation of the contour integral
storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))
# fill up the magnetic entities dictionary with the energies
mag_ent["energies"].append(storage)
# iterate over the pairs
for tracker, pair in enumerate(pairs):
# iterate over the quantization axes
for i, (Gij, Gji) in enumerate(zip(pair["Gij"], pair["Gji"])):
site_i = magnetic_entities[pair["ai"]]
site_j = magnetic_entities[pair["aj"]]
storage = []
# iterate over the first order local perturbations in all possible orientations for the two sites
for Vui in site_i["Vu1"][i]:
for Vuj in site_j["Vu1"][i]:
# The Szunyogh-Lichtenstein formula
traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2)
# evaluation of the contour integral
storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))
# fill up the pairs dictionary with the energies
pairs[tracker]["energies"].append(storage)
end_time = timer()
print("############################### GROGU OUTPUT ###################################")
print("================================================================================")
print("Input file: ")
print(simulation_parameters["path"])
print("Number of nodes in the parallel cluster: ", simulation_parameters["parallel_size"])
print("================================================================================")
try:
print("Cell [Ang]: ")
print(simulation_parameters["geom"].cell)
except:
print("Geometry could not be read.")
print("================================================================================")
print("DFT axis: ")
print(simulation_parameters["scf_xcf_orientation"])
print("Quantization axis and perpendicular rotation directions:")
for ref in ref_xcf_orientations:
print(ref["o"], " --» ", ref["vw"])
print("================================================================================")
print("number of k points: ", simulation_parameters["kset"])
print("k point directions: ", simulation_parameters["kdirs"])
print("================================================================================")
print("Parameters for the contour integral:")
print("Ebot: ", simulation_parameters["ebot"])
print("Eset: ", simulation_parameters["eset"])
print("Esetp: ", simulation_parameters["esetp"])
print("================================================================================")
print("Atomic informations: ")
print("")
print("")
print("Not yet specified.")
print("")
print("")
print("================================================================================")
print("Exchange [meV]")
print("--------------------------------------------------------------------------------")
print("Atom1 Atom2 [i j k] d [Ang]")
print("--------------------------------------------------------------------------------")
for pair in pairs:
J_iso, J_S, D = calculate_exchange_tensor(pair)
J_iso = J_iso * sisl.unit_convert("eV", "meV")
J_S = J_S * sisl.unit_convert("eV", "meV")
D = D * sisl.unit_convert("eV", "meV")
print(pair["ai"], pair["aj"], pair["Ruc"], "distance")
print("Isotropic: ", J_iso)
print("DMI: ", D)
print("Symmetric-anisotropy: ", J_S)
print("")
print("================================================================================")
print("Runtime information: ")
print("Total runtime: ", end_time - start_time)
print("--------------------------------------------------------------------------------")
print("Initial setup: ", setup_time - start_time)
print(f"Hamiltonian conversion and XC field extraction: {H_and_XCF_time - setup_time:.3f} s")
print(f"Pair and site datastructure creatrions: {site_and_pair_dictionaries_time - H_and_XCF_time:.3f} s")
print(f"k set cration and distribution: {k_set_time - site_and_pair_dictionaries_time:.3f} s")
print(f"Rotating XC potential: {reference_rotations_time - k_set_time:.3f} s")
print(f"Greens function inversion: {green_function_inversion_time - reference_rotations_time:.3f} s")
print(f"Calculate energies and magnetic components: {end_time - green_function_inversion_time:.3f} s")
Loading…
Cancel
Save