diff --git a/README.md b/README.md index 45b4f89..a7063c3 100644 --- a/README.md +++ b/README.md @@ -2,12 +2,12 @@ # TODO -- Definition of magnetic entities: +[x] Definition of magnetic entities: * Through simple sequence o forbitals in the unit cell * Through atom specification * Through atom and orbital specification -- Separation of TR and TRB components of the Hamiltonian, Identification of the exchange field. -- Definition of commutator expressions, old projection matrix elements -- Efficient calculation of Green's functions -- Calculation of energy and momentum resolved derivatives -- Parallel BZ and serial energy integral +[x] Separation of TR and TRB components of the Hamiltonian, Identification of the exchange field. +[x] Definition of commutator expressions, old projection matrix elements +[x] Efficient calculation of Green's functions +[] Calculation of energy and momentum resolved derivatives +[] Parallel BZ and serial energy integral diff --git a/src/grogu/useful.py b/src/grogu/useful.py index df49241..6b02e23 100644 --- a/src/grogu/useful.py +++ b/src/grogu/useful.py @@ -23,7 +23,28 @@ from itertools import permutations, product import numpy as np 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 +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): @@ -80,14 +101,7 @@ def make_kset(dirs="xyz", NUMK=20): return kset -# 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]]) - - -def comm(a, b): +def commutator(a, b): "Shorthand for commutator" return a @ b - b @ a @@ -104,7 +118,7 @@ def tau_u(u): def crossM(u): """ 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]]) @@ -152,6 +166,7 @@ def spin_tracer(M): 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) @@ -197,3 +212,15 @@ def blow_up_orbindx(orb_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() + + +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() diff --git a/test.ipynb b/test.ipynb index c803d0a..70e6293 100644 --- a/test.ipynb +++ b/test.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 10, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ @@ -11,26 +11,12 @@ "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[\"VECLIB_MAXIMUM_THREADS\"] = \"4\" # export VECLIB_MAXIMUM_THREADS=4\n", - "os.environ[\"NUMEXPR_NUM_THREADS\"] = \"4\" # export NUMEXPR_NUM_THREADS=6" - ] - }, - { - "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 *" + "os.environ[\"NUMEXPR_NUM_THREADS\"] = \"4\"" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 15, "metadata": {}, "outputs": [ { @@ -39,21 +25,35 @@ "'0.14.3'" ] }, - "execution_count": 12, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "np.set_printoptions(linewidth=2*75)\n", - "sisl.__version__" + "import numpy as np\n", + "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", - "execution_count": 13, + "execution_count": 16, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of nodes in the parallel cluster: 1\n" + ] + } + ], "source": [ "# this cell mimicks an input file\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", "# 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", - "# 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", " 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", @@ -80,60 +80,51 @@ " dict(ai=0,aj=2,Ruc=np.array([-1,0,0])),\n", " dict(ai=1,aj=2,Ruc=np.array([-1,0,0]))]\n", "\n", - "# Bz sampling\n", - "nk=np.array([10,10,0])" + "# Brilloun zone sampling and Green function contour integral\n", + "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", - "execution_count": 14, + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# digestion of the input\n", - "# read in geometry and hamiltonian\n", - "geo = fdf.read_geometry()\n", + "# read in hamiltonian\n", "dh = fdf.read_hamiltonian()\n", "\n", + "# unit cell index\n", "uc_in_sc_idx=dh.lattice.sc_index([0,0,0])\n", "\n", - "# get indecies for orbitals of the magnetic entities\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": [ + "\n", "# WE WILL NOT NEED THIS!!\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" - ] - }, - { - "cell_type": "code", - "execution_count": 16, - "metadata": {}, - "outputs": [], - "source": [ + "\n", + "\n", "NO=dh.no # shorthand for number of orbitals in the unit cell\n", + "\n", "# preprocessing Hamiltonian and overlap matrix elements\n", "h11 = dh.tocsr(dh.M11r)\n", "h11 += dh.tocsr(dh.M11i)*1.0j\n", @@ -153,20 +144,15 @@ "\n", "sov = dh.tocsr(dh.S_idx).toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n", "\n", - "# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\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", "# This is the permutation that transforms ud1ud2 to u12d12\n", "# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n", - "\n", "# the inverse transformation is U.T u12d12 to ud1ud2\n", "# That is FROM ORBITAL BOX to SPIN BOX => U.T\n", "\n", - "# this is a test\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", + "# From now on everything is in SPIN BOX!!\n", "hh,ss = \\\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", @@ -189,516 +175,341 @@ "hTRB = (hh-hTR)/2\n", "\n", "# extracting the exchange field\n", - "traced=[spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())]\n", - "\n", + "traced=[spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77\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['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", "max_xcfs=abs(np.array(np.array([f['c'] for f in traced]))).max()\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", - "execution_count": 17, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ - "# collecting information for all proposed rotations, this might generate a lot of unnceessary data\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", - "# Transformation matrices for rotating the exchange field\n", - "for i,orient in enumerate(ref_xcf_orientations):\n", - "\n", - " # obtain rotated exchange field \n", - " R=RotMa2b(scf_xcf_orientation,orient['o'])\n", - " \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\n", - " \n", - " # Update magnetic entities with commutator data for each \n", - " for mag_ent in magnetic_entities:\n", - " mag_ent['Commutator_Data'] = []\n", + "# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n", + "for i, mag_ent in enumerate(magnetic_entities):\n", + " parsed = parse_magnetic_entity(dh,**mag_ent) # parse orbital indexes\n", + " magnetic_entities[i]['orbital_indeces'] = parsed\n", + " magnetic_entities[i]['spin_box_indeces'] = blow_up_orbindx(parsed) # calculate spin box indexes\n", + " spin_box_shape = len(mag_ent[\"spin_box_indeces\"]) # calculate size for Greens function generation\n", "\n", - " for u in orient['vw']:\n", + " mag_ent['energies'] = [] # we will store the second order energy derivations here\n", " \n", - " Tu = np.kron(np.eye(NO,dtype=int),tau_u(u))\n", + " mag_ent['Gii'] = [] # Greens function\n", + " mag_ent['Gii_tmp'] = [] # Greens function for parallelization\n", + " mag_ent[\"Vu1\"] = [ list([]) for _ in range(len(ref_xcf_orientations))] # These will be the perturbed potentials from eq. 100\n", + " mag_ent[\"Vu2\"] = [ list([]) for _ in range(len(ref_xcf_orientations))]\n", + " for i in ref_xcf_orientations:\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", + " pair['energies'] = [] # we will store the second order energy derivations here\n", "\n", - " Vu1 = 1j/2*comm(rot_H_XCF_uc,Tu)\n", - " Vu2 = 1/8*comm(comm(Tu,rot_H_XCF_uc),Tu)\n", + " pair['Gij'] = [] # Greens function\n", + " pair['Gji'] = []\n", + " pair['Gij_tmp'] = [] # Greens function for parallelization\n", + " pair['Gji_tmp'] = []\n", "\n", - " for mag_ent in magnetic_entities:\n", - " idx = mag_ent['spin_box_indeces']\n", - " mag_ent['Commutator_Data'].append(dict(Vu1=Vu1[idx,:][:,idx],\n", - " Vu2=Vu2[idx,:][:,idx]))\n", + " pair[\"Vij\"] = [list([]) for _ in range(len(ref_xcf_orientations))] # These will be the perturbed potentials from eq. 100\n", + " pair[\"Vji\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n", "\n", - "\n", - "# TODO STORE ABOVE INFO!!!!!" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": 18, - "metadata": {}, - "outputs": [], - "source": [ - "k = np.array([0.11441,0.234432,0.0])\n", - "\n" + " for i in ref_xcf_orientations: \n", + " pair['Gij'].append(np.zeros((eset,spin_box_shape_i,spin_box_shape_j),dtype='complex128'))\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", + " pair['Gji'].append(np.zeros((eset,spin_box_shape_j,spin_box_shape_i),dtype='complex128'))\n", + " pair['Gji_tmp'].append(np.zeros((eset,spin_box_shape_j,spin_box_shape_i),dtype='complex128'))\n" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, - "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": [ - "#This goes inside the for loop for k \n", - "k = np.array([0,0.1,0.0])\n", - "\n", - "\n", - "phases = np.exp(-2j*np.pi*k@dh.lattice.sc_off.T)\n", - "\n", - "np.exp(-1j * np.dot(np.dot(np.dot(dh.rcell, k), dh.cell),dh.lattice.sc_off.T))\n", - "ss.shape,rot_H_XCF.shape,hTRS.shape,phases.shape\n", - "\n", - "HK=np.einsum('abc,a->bc',hTRS+rot_H_XCF,phases)\n", - "SK=np.einsum('abc,a->bc',ss,phases)\n", - "\n", - "vals,vecs=eigh(HK,SK)\n", - "\n", - "#This goes inside the for loop for ze\n", - "ze = -.1+0.000j\n", - "GK=vecs@np.diag(1/(ze-vals))@vecs.conj().T\n", - "\n", - "#phase=np.exp(1j*np.dot(np.dot(k,dh.rcell),offset))\n", - "Gii=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]\n", - "Gjj=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]\n", - "Gij=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]*np.exp(-k@offset*2*np.pi)\n", - "Gji=GK[magn_ent_idices[ai],:][:,magn_ent_idices[ai]]*np.exp( k@offset*2*np.pi)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, "outputs": [], "source": [ - "np.allclose(np.linalg.inv(ze*SK-HK),GK), abs(np.linalg.inv(ze*SK-HK)-GK).max()" + "kset = make_kset(dirs=kdirs,NUMK=kset) # generate k space sampling\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", - "execution_count": null, + "execution_count": 20, "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", - " # 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", + "# collecting information for all proposed rotations, this might generate a lot of unnecessary data\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", - " 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", + "# this will contain all the data needed to calculate the energy variations upon rotation\n", + "hamiltonians = []\n", "\n", - " SK = np.einsum('abc,c->ab',dh.sov,phases)\n", + "# iterate over the reference directions (quantization axes)\n", + "for i,orient in enumerate(ref_xcf_orientations):\n", "\n", - " return HK11,HK12,HK21,HK22,SK\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ + " # obtain rotated exchange field \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", - "rot_xcf_mats=[]\n", - "for i,rot in enumerate(rots):\n", - " R=RotM(*rot) # rotation matrix to go from scf to reference direction\n", - " R@scf_xcf_orientation # exchange field orientation \n", - "\n", - "\n", - "rot_xcf_mats = [RotM(*rot) for rot in rots] \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", + " hamiltonians.append(dict(orient=orient['o'], H=rot_H, rotations=[])) # store orientation and rotated Hamiltonian\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", + " 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", - "cummutators=[]\n", + " for mag_ent in magnetic_entities:\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", + " mag_ent[\"Vu2\"][i].append(Vu2[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :])\n", "\n", - "\n" + " for pair in pairs:\n", + " ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"] # get the pair orbital sizes from the magnetic entities\n", + " aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n", + " pair[\"Vij\"][i].append(Vu1[:, ai][aj, :]) # fill up the perturbed potentials (for now) based on the on-site projections\n", + " pair[\"Vji\"][i].append(Vu1[:, aj][ai, :])\n" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 21, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Number of magnetic entities being calculated: 4\n", + "We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.\n", + "The shape of the Hamiltonian and the Greens function is 84x84.\n" + ] + } + ], "source": [ - "i=0\n", + "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", + " # 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", + " for pair in pairs:\n", + " # add phase shift based on the cell difference\n", + " phase=np.exp(1j * 2 * np.pi * k @ pair[\"Ruc\"].T)\n", + " \n", + " # get the pair orbital sizes from the magnetic entities\n", + " ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]\n", + " aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\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", + "# 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", - "dh.hxc_rots[i][uc_in_sc_idx]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "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()" + " 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", - "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, + "execution_count": 22, "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", + "# 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", - " # 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", + "# 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", + " storage = []\n", + " # iterate over the first and second order local perturbations\n", + " for Vu1, Vu2 in zip(mag_ent[\"Vu1\"][i], mag_ent[\"Vu2\"][i]):\n", + " # The Szunyogh-Lichtenstein formula\n", + " 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", "\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", + " # fill up the magnetic entities dictionary with the energies\n", + " idx = np.array([len(mag_ent_parallel_set[i]) for i in range(rank)]).sum() + tracker\n", + " magnetic_entities[int(idx)][\"energies\"].append(storage)\n", "\n", - " SK = np.einsum('abc,c->ab',dh.sov,phases)\n", + "# iterate over the pairs\n", + "for tracker, pair in enumerate(pairs_parallel_set[rank]):\n", + " # iterate over the quantization axes\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", - " return HK11,HK12,HK21,HK22,SK\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "myHk=(lambda x:np.vstack([np.hstack([x[0],x[1]]),\n", - " np.hstack([x[2],x[3]])]))(hsk(dh))\n", - "mySk=(lambda x:np.vstack([np.hstack([x[-1],0*x[-1]]),\n", - " np.hstack([0*x[-1],x[-1]])]))(hsk(dh))\n", - "\n", - "myHkshift=myHk+EF*mySk" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "H0=np.vstack( [np.hstack([dh.h11[:,:,0],dh.h12[:,:,0]]),\n", - " np.hstack([dh.h21[:,:,0],dh.h22[:,:,0]])])\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", - "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", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "plt.pcolor((U.T@H0shifted@U)[:5,:5])\n", - "plt.colorbar()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "np.linalg.eigvalsh(myHk+EF*mySk)" + " # 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", - "execution_count": null, + "execution_count": 23, "metadata": {}, "outputs": [], "source": [ - "dh.Hk().toarray()[:3,:3]" + "def calculate_exchange_tensor(pair):\n", + " Eo, Ev, Ew = pair[\"energies\"]\n", + " J = np.array([Ew[-1], Ev[-1], Eo[0]]) # xx, yy, zz\n", + " JS = -0.5 * np.array([Eo[1] + Eo[2], Ev[1] + Ev[2], Ew[1] + Ew[2]]) # yz, zx, xy\n", + " D = 0.5 * np.array([Eo[1] - Eo[2], Ev[2] - Ev[1], Ew[1] - Ew[2]]) # x, y, z\n", + " return J.sum()/3 * 1000" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 24, "metadata": {}, - "outputs": [], - "source": [ - "plt.matshow(np.abs(dh.Hk().toarray()-myHk))\n", - "plt.colorbar()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-60.72359053548613\n", + "-60.531975234450115\n", + "-60.524676226428824\n", + "-6.55042989834691\n", + "-6.047933492978864\n" + ] + } + ], "source": [ - "plt.plot(np.linalg.eigvalsh(dh.Hk().toarray())-np.linalg.eigvalsh(myHk))" + "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", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, + "execution_count": 25, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "-6.043716409664797" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "magnetic_entities = [parse_magnetic_entity(dh,atom=[0,1],l=2),parse_magnetic_entity(dh,atom=2,l=2)]" + "-61.33097171216109\n", + "-60.52198328932686\n", + "-60.51657719027764\n", + "-6.545208546361317\n", + "-6.043716409664797" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", + "execution_count": 26, "metadata": {}, + "outputs": [ + { + "ename": "SyntaxError", + "evalue": "invalid syntax (1876172784.py, line 5)", + "output_type": "error", + "traceback": [ + "\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" + ] + } + ], "source": [ - "# TODO\n", - "\n", - "- Definition of magnetic entities:\n", - " * Through simple sequence o forbitals in the unit cell\n", - " * Through atom specification\n", - " * Through atom and orbital specification\n", - "- Separation of TR and TRB components of the Hamiltonian, Identification of the exchange field. \n", - "- Definition of commutator expressions, old projection matrix elements\n", - "- Efficient calculation of Green's functions\n", - "- Calculation of energy and momentum resolved derivatives\n", - "- Parallel BZ and serial energy integral" + "# symmetrizing Hamiltonian and overlap matrix to make them hermitian \n", + "# Check if exchange field has scalar part\n", + "# parallel over integrals\n", + "\n", + "========================================\n", + " \n", + "Atom Angstrom\n", + "# Label, x y z Sx Sy Sz #Q Lx Ly Lz Jx Jy Jz\n", + "--------------------------------------------------------------------------------------------------------------------------------------------------------------------------\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", + "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": { diff --git a/test.py b/test.py new file mode 100644 index 0000000..29c1380 --- /dev/null +++ b/test.py @@ -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")