@ -2,35 +2,29 @@
"cells": [
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"import os\n",
"os.environ[\"OMP_NUM_THREADS\"] = \"4\" # export OMP_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[\"VECLIB_MAXIMUM_THREADS\"] = \"4\" # export VECLIB_MAXIMUM_THREADS=4\n",
"os.environ[\"NUMEXPR_NUM_THREADS\"] = \"4\""
]
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.14.3'"
"name": "stderr",
"output_type": "stream",
"text": [
"[Daniels-MacBook-Air.local:09850] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Daniels-MacBook-Air.501/jf.0/320733184/sm_segment.Daniels-MacBook-Air.501.131e0000.0 could be created.\n"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import os\n",
"from sys import stdout\n",
"from tqdm import tqdm\n",
"from timeit import default_timer as timer\n",
"\n",
"os.environ[\"OMP_NUM_THREADS\"] = \"1\" # export OMP_NUM_THREADS=4\n",
"os.environ[\"OPENBLAS_NUM_THREADS\"] = \"1\" # export OPENBLAS_NUM_THREADS=4\n",
"os.environ[\"MKL_NUM_THREADS\"] = \"1\" # export MKL_NUM_THREADS=6\n",
"os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"1\" # export VECLIB_MAXIMUM_THREADS=4\n",
"os.environ[\"NUMEXPR_NUM_THREADS\"] = \"1\" # export NUMEXPR_NUM_THREADS=6\n",
"\n",
"import numpy as np\n",
"import sisl\n",
"from grogu.useful import *\n",
@ -38,12 +32,12 @@
"from numpy.linalg import inv\n",
"import warnings\n",
"\n",
"sisl.__version__ \n"
"start_time = timer() \n"
]
},
{
"cell_type": "code",
"execution_count": 16 ,
"execution_count": 2 ,
"metadata": {},
"outputs": [
{
@ -56,34 +50,42 @@
],
"source": [
"# 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(\n",
" \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\"\n",
")\n",
"# this information needs to be given at the input!!\n",
"scf_xcf_orientation = np.array([0, 0, 1]) # z\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",
"# at this moment the user has to supply o,v,w on the input.\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 = [\n",
" 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",
" dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),\n",
"]\n",
"\n",
"# human readable definition of magnetic entities\n",
"magnetic_entities=[dict(atom=3 ,l=2),\n",
"magnetic_entities = [\n",
" dict(atom=3, l=2),\n",
" dict(atom=4, l=2),\n",
" dict(atom=5, l=2),\n",
" dict(atom=[3,4],)]\n",
"# dict(atom=[3, 4]),\n",
"]\n",
"\n",
"# pair information\n",
"pairs=[dict(ai=0,aj=1,Ruc=np.array([0,0,0])), # isotropic should be -82 meV\n",
"pairs = [\n",
" dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV\n",
" dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])), # these should all be around -41.9 in the isotropic part\n",
" dict(ai=1,aj=2,Ruc=np.array([0,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([0, 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",
"]\n",
"\n",
"# Brilloun zone sampling and Green function contour integral\n",
"kset = 20\n",
"kdirs = \"xy\"\n",
"ebot = -4 0\n",
"ebot = -3 0\n",
"eset = 50\n",
"esetp = 10000\n",
"\n",
@ -94,70 +96,89 @@
"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": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
" print(\"Number of nodes in the parallel cluster: \", size)\n",
"\n",
"simulation_parameters = dict(path=\"Not yet specified.\",\n",
" scf_xcf_orientation=scf_xcf_orientation, \n",
" ref_xcf_orientations=ref_xcf_orientations,\n",
" kset=kset,\n",
" kdirs=kdirs, \n",
" ebot=ebot,\n",
" eset=eset, \n",
" esetp=esetp,\n",
" parallel_size=size)\n",
"\n",
"# digestion of the input\n",
"# read in hamiltonian\n",
"dh = fdf.read_hamiltonian()\n",
"try:\n",
" simulation_parameters[\"geom\"] = fdf.read_geometry()\n",
"except:\n",
" print(\"Error reading geometry.\")\n",
"\n",
"# unit cell index\n",
"uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])\n",
"\n",
"\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",
"\n",
"\n",
"setup_time = timer()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"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",
"h11 = h11.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128' )\n",
"h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\" )\n",
"\n",
"h22 = dh.tocsr(dh.M22r)\n",
"h22 += dh.tocsr(dh.M22i) * 1.0j\n",
"h22 = h22.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128' )\n",
"h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\" )\n",
"\n",
"h12 = dh.tocsr(dh.M12r)\n",
"h12 += dh.tocsr(dh.M12i) * 1.0j\n",
"h12 = h12.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128' )\n",
"h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\" )\n",
"\n",
"h21 = dh.tocsr(dh.M21r)\n",
"h21 += dh.tocsr(dh.M21i) * 1.0j\n",
"h21 = h21.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128' )\n",
"h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\" )\n",
"\n",
"sov = dh.tocsr(dh.S_idx).toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n",
"sov = (\n",
" dh.tocsr(dh.S_idx)\n",
" .toarray()\n",
" .reshape(NO, dh.n_s, NO)\n",
" .transpose(0, 2, 1)\n",
" .astype(\"complex128\")\n",
")\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(\n",
" [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])]\n",
")\n",
"# This is the permutation that transforms ud1ud2 to u12d12\n",
"# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n",
"# the inverse transformation is U.T u12d12 to ud1ud2\n",
"# That is FROM ORBITAL BOX to SPIN BOX => U.T\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",
"np.array([U.T@np.block([[sov[:,:,i] ,sov[:,:,i]*0],\n",
" [sov[:,:,i]*0,sov[:,:,i] ]])@U for i in range(dh.lattice.nsc.prod())])\n",
"hh, ss = np.array(\n",
" [\n",
" U.T @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) @ U\n",
" for i in range(dh.lattice.nsc.prod())\n",
" ]\n",
"), np.array(\n",
" [\n",
" U.T\n",
" @ np.block([[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]])\n",
" @ U\n",
" for i in range(dh.lattice.nsc.prod())\n",
" ]\n",
")\n",
"\n",
"\n",
"# symmetrizing Hamiltonian and overlap matrix to make them hermitian\n",
@ -176,138 +197,195 @@
"\n",
"# extracting the exchange field\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])]) # equation 77\n",
"XCF = np.array(\n",
" [\n",
" 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",
" ]\n",
") # 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",
"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 negligible scalar part. Largest value is {max_xcfs}\")\n"
" warnings.warn(\n",
" f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\"\n",
" )\n",
"\n",
"H_and_XCF_time = timer()"
]
},
{
"cell_type": "code",
"execution_count": 18 ,
"execution_count": 4 ,
"metadata": {},
"outputs": [],
"source": [
"# 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",
" mag_ent['energies'] = [] # we will store the second order energy derivations here\n",
" \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",
" magnetic_entities[i][\"orbital_indeces\"] = parsed\n",
" magnetic_entities[i][\"spin_box_indeces\"] = blow_up_orbindx(\n",
" parsed\n",
" ) # calculate spin box indexes\n",
" spin_box_shape = len(\n",
" mag_ent[\"spin_box_indeces\"]\n",
" ) # calculate size for Greens function generation\n",
"\n",
" mag_ent[\"energies\"] = [] # we will store the second order energy derivations here\n",
"\n",
" mag_ent[\"Gii\"] = [] # Greens function\n",
" mag_ent[\"Gii_tmp\"] = [] # Greens function for parallelization\n",
" mag_ent[\"Vu1\"] = [\n",
" list([]) for _ in range(len(ref_xcf_orientations))\n",
" ] # 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",
" mag_ent[\"Gii\"].append(\n",
" np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n",
" ) # Greens functions for every quantization axis\n",
" mag_ent[\"Gii_tmp\"].append(\n",
" np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n",
" )\n",
"\n",
"# for every site we have to store 2x3 Greens function (and the associated _tmp-s)\n",
"# 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",
" pair['Gij'] = [] # Greens function\n",
" pair['Gji'] = []\n",
" pair['Gij_tmp'] = [] # Greens function for parallelization\n",
" pair['Gji_tmp'] = []\n",
"\n",
" pair[\"Vij\"] = [list([]) for _ in range(len(ref_xcf_orientations))] # These will be the perturbed potentials from eq. 100\n",
" spin_box_shape_i, spin_box_shape_j = len(\n",
" magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]\n",
" ), len(\n",
" magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n",
" ) # calculate size for Greens function generation\n",
"\n",
" pair[\"energies\"] = [] # we will store the second order energy derivations here\n",
"\n",
" pair[\"Gij\"] = [] # Greens function\n",
" pair[\"Gji\"] = []\n",
" pair[\"Gij_tmp\"] = [] # Greens function for parallelization\n",
" pair[\"Gji_tmp\"] = []\n",
"\n",
" pair[\"Vij\"] = [\n",
" list([]) for _ in range(len(ref_xcf_orientations))\n",
" ] # These will be the perturbed potentials from eq. 100\n",
" pair[\"Vji\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\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"
" pair[\"Gij\"].append(\n",
" np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n",
" )\n",
" pair[\"Gij_tmp\"].append(\n",
" np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n",
" ) # Greens functions for every quantization axis\n",
" pair[\"Gji\"].append(\n",
" np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n",
" )\n",
" pair[\"Gji_tmp\"].append(\n",
" np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n",
" )\n",
"\n",
"site_and_pair_dictionaries_time = timer()\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"execution_count": 5 ,
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"k loop: 0%| | 0/400 [00:00<?, ?it/s]"
]
}
],
"source": [
"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"
"kpcs = np.array_split(kset, size) # split the k points based on MPI size\n",
"kpcs[root_node] = tqdm(kpcs[root_node], desc='k loop', file=stdout)\n",
"\n",
"k_set_time = timer()"
]
},
{
"cell_type": "code",
"execution_count": 20 ,
"execution_count": 6 ,
"metadata": {},
"outputs": [],
"source": [
"# 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",
"# this will contain all the data needed to calculate the energy variations upon rotation\n",
"hamiltonians = []\n",
"\n",
"# iterate over the reference directions (quantization axes)\n",
"for i, orient in enumerate(ref_xcf_orientations):\n",
"\n",
" # 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",
" R = RotMa2b(scf_xcf_orientation, orient[\"o\"])\n",
" rot_XCF = np.einsum(\"ij,jklm->iklm\", R, XCF)\n",
" rot_H_XCF = sum(\n",
" [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])]\n",
" )\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",
" hamiltonians.append(dict(orient=orient['o'], H=rot_H, rotations=[])) # store orientation and rotated Hamiltonian\n",
" hamiltonians.append(\n",
" dict(orient=orient[\"o\"], H=rot_H, rotations=[])\n",
" ) # 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",
" for u in orient[\n",
" \"vw\"\n",
" ]: # 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",
" 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",
" mag_ent[\"Vu1\"][i].append(\n",
" Vu1[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n",
" ) # fill up the perturbed potentials (for now) based on the on-site projections\n",
" mag_ent[\"Vu2\"][i].append(\n",
" 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",
" ai = magnetic_entities[pair[\"ai\"]][\n",
" \"spin_box_indeces\"\n",
" ] # 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"
" pair[\"Vij\"][i].append(\n",
" Vu1[:, ai][aj, :]\n",
" ) # fill up the perturbed potentials (for now) based on the on-site projections\n",
" pair[\"Vji\"][i].append(Vu1[:, aj][ai, :])\n",
"\n",
"reference_rotations_time = timer()"
]
},
{
"cell_type": "code",
"execution_count": 21 ,
"execution_count": 7 ,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of magnetic entities being calculated: 4 \n",
"Number of magnetic entities being calculated: 3 \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"
"The shape of the Hamiltonian and the Greens function is 84x84.\n",
"k loop: 100%|██████████| 400/400 [03:05<00:00, 2.15it/s]\n"
]
}
],
"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",
" print(\"Number of magnetic entities being calculated: \", len(magnetic_entities))\n",
" print(\n",
" \"We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.\"\n",
" )\n",
" print(f\"The shape of the Hamiltonian and the Greens function is {NO}x{NO}.\")\n",
"comm.Barrier()\n",
"# ----------------------------------------------------------------------\n",
"\n",
@ -321,7 +399,9 @@
"# 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",
" for i, hamiltonian_orientation in enumerate(\n",
" hamiltonians\n",
" ): # iterate over reference directions\n",
" # calculate Greens function\n",
" H = hamiltonian_orientation[\"H\"]\n",
" HK, SK = hsk(H, ss, dh.sc_off, k)\n",
@ -329,7 +409,10 @@
"\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",
" mag_ent[\"Gii_tmp\"][i] += (\n",
" Gk[:, mag_ent[\"spin_box_indeces\"]][..., mag_ent[\"spin_box_indeces\"]]\n",
" * wk\n",
" )\n",
"\n",
" for pair in pairs:\n",
" # add phase shift based on the cell difference\n",
@ -340,8 +423,8 @@
" 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",
" 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",
@ -350,21 +433,86 @@
"\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"
" comm.Reduce(pair[\"Gji_tmp\"][i], pair[\"Gji\"][i], root=root_node)\n",
"\n",
"green_function_inversion_time = timer()\n"
]
},
{
"cell_type": "code",
"execution_count": 22 ,
"execution_count": 8 ,
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"############################### GROGU OUTPUT ###################################\n",
"================================================================================\n",
"Input file: \n",
"Not yet specified.\n",
"Number of nodes in the parallel cluster: 1\n",
"================================================================================\n",
"Cell [Ang]: \n",
"[[ 3.79100000e+00 0.00000000e+00 0.00000000e+00]\n",
" [-1.89550000e+00 3.28310231e+00 0.00000000e+00]\n",
" [ 1.25954923e-15 2.18160327e-15 2.05700000e+01]]\n",
"================================================================================\n",
"DFT axis: \n",
"[0 0 1]\n",
"Quantization axis and perpendicular rotation directions:\n",
"[1 0 0] --» [array([0, 1, 0]), array([0, 0, 1])]\n",
"[0 1 0] --» [array([1, 0, 0]), array([0, 0, 1])]\n",
"[0 0 1] --» [array([1, 0, 0]), array([0, 1, 0])]\n",
"================================================================================\n",
"number of k points: 20\n",
"k point directions: xy\n",
"================================================================================\n",
"Parameters for the contour integral:\n",
"Ebot: -30\n",
"Eset: 50\n",
"Esetp: 10000\n",
"================================================================================\n",
"Atomic informations: \n",
"\n",
"\n",
"Not yet specified.\n",
"\n",
"\n",
"================================================================================\n",
"Exchange [meV]\n",
"--------------------------------------------------------------------------------\n",
"Atom1 Atom2 [i j k] d [Ang]\n",
"--------------------------------------------------------------------------------\n",
"0 1 [0 0 0] distance\n",
"Isotropic: -60.893309223088636\n",
"DMI: [-5.95741551e+00 7.27737654e+00 6.90431275e-04 -8.11057566e-04\n",
" -5.49031203e-06]\n",
"Symmetric-anisotropy: [-9.32966923e-01 -8.92579299e-04 -2.04258660e-06]\n",
"\n",
"0 2 [0 0 0] distance\n",
"Isotropic: -60.556512255197724\n",
"DMI: [-0.34565796 0.65919757 0.07106945 -6.23149871 -0.0424978 ]\n",
"Symmetric-anisotropy: [ 3.78506176e+00 -6.13838308e+00 3.59037036e-03]\n",
"\n",
"================================================================================\n",
"Runtime information: \n",
"Total runtime: 187.525066541\n",
"--------------------------------------------------------------------------------\n",
"Initial setup: 0.275851458\n",
"Hamiltonian conversion and XC field extraction: 1.232 s\n",
"Pair and site datastructure creatrions: 0.025 s\n",
"k set cration and distribution: 0.029 s\n",
"Rotating XC potential: 0.366 s\n",
"Greens function inversion: 185.537 s\n",
"Calculate energies and magnetic components: 0.060 s\n"
]
}
],
"source": [
"# 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",
"if rank == root_node:\n",
" # iterate over the magnetic entities\n",
"for tracker, mag_ent in enumerate(mag_ent_parallel_set[rank]):\n",
" for tracker, mag_ent in enumerate(magnetic_entities):\n",
" # iterate over the quantization axes\n",
" for i, Gii in enumerate(mag_ent[\"Gii\"]):\n",
" storage = []\n",
@ -376,11 +524,10 @@
" storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n",
"\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",
" mag_ent[\"energies\"].append(storage)\n",
"\n",
" # iterate over the pairs\n",
"for tracker, pair in enumerate(pairs_parallel_set[rank] ):\n",
" for tracker, pair in enumerate(pairs):\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",
@ -396,38 +543,96 @@
" storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n",
"\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",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"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"
" pairs[tracker][\"energies\"].append(storage)\n",
"\n",
" end_time = timer()\n",
"\n",
" print(\"############################### GROGU OUTPUT ###################################\")\n",
" print(\"================================================================================\")\n",
" print(\"Input file: \")\n",
" print(simulation_parameters[\"path\"])\n",
" print(\"Number of nodes in the parallel cluster: \", simulation_parameters[\"parallel_size\"])\n",
" print(\"================================================================================\")\n",
" try:\n",
" print(\"Cell [Ang]: \")\n",
" print(simulation_parameters[\"geom\"].cell)\n",
" except:\n",
" print(\"Geometry could not be read.\")\n",
" print(\"================================================================================\")\n",
" print(\"DFT axis: \")\n",
" print(simulation_parameters[\"scf_xcf_orientation\"])\n",
" print(\"Quantization axis and perpendicular rotation directions:\")\n",
" for ref in ref_xcf_orientations:\n",
" print(ref[\"o\"], \" --» \", ref[\"vw\"])\n",
" print(\"================================================================================\")\n",
" print(\"number of k points: \", simulation_parameters[\"kset\"])\n",
" print(\"k point directions: \", simulation_parameters[\"kdirs\"])\n",
" print(\"================================================================================\")\n",
" print(\"Parameters for the contour integral:\")\n",
" print(\"Ebot: \", simulation_parameters[\"ebot\"])\n",
" print(\"Eset: \", simulation_parameters[\"eset\"])\n",
" print(\"Esetp: \", simulation_parameters[\"esetp\"])\n",
" print(\"================================================================================\")\n",
" print(\"Atomic informations: \")\n",
" print(\"\")\n",
" print(\"\")\n",
" print(\"Not yet specified.\")\n",
" print(\"\")\n",
" print(\"\")\n",
" print(\"================================================================================\")\n",
" print(\"Exchange [meV]\")\n",
" print(\"--------------------------------------------------------------------------------\")\n",
" print(\"Atom1 Atom2 [i j k] d [Ang]\")\n",
" print(\"--------------------------------------------------------------------------------\")\n",
" for pair in pairs:\n",
" J_iso, J_S, D = calculate_exchange_tensor(pair)\n",
" J_iso = J_iso * sisl.unit_convert(\"eV\", \"meV\")\n",
" J_S = J_S * sisl.unit_convert(\"eV\", \"meV\")\n",
" D = D * sisl.unit_convert(\"eV\", \"meV\")\n",
" \n",
" print(pair[\"ai\"], pair[\"aj\"], pair[\"Ruc\"], \"distance\")\n",
" print(\"Isotropic: \", J_iso)\n",
" print(\"DMI: \", D)\n",
" print(\"Symmetric-anisotropy: \", J_S)\n",
" print(\"\")\n",
" \n",
" print(\"================================================================================\")\n",
" print(\"Runtime information: \")\n",
" print(\"Total runtime: \", end_time - start_time)\n",
" print(\"--------------------------------------------------------------------------------\")\n",
" print(\"Initial setup: \", setup_time - start_time)\n",
" print(f\"Hamiltonian conversion and XC field extraction: {H_and_XCF_time - setup_time:.3f} s\")\n",
" print(f\"Pair and site datastructure creatrions: {site_and_pair_dictionaries_time - H_and_XCF_time:.3f} s\")\n",
" print(f\"k set cration and distribution: {k_set_time - site_and_pair_dictionaries_time:.3f} s\")\n",
" print(f\"Rotating XC potential: {reference_rotations_time - k_set_time:.3f} s\")\n",
" print(f\"Greens function inversion: {green_function_inversion_time - reference_rotations_time:.3f} s\")\n",
" print(f\"Calculate energies and magnetic components: {end_time - green_function_inversion_time:.3f} s\")\n"
]
},
{
"cell_type": "code",
"execution_count": 24,
"execution_count": 9 ,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-60.72359053548613\n",
"-60.531975234450115\n",
"-60.524676226428824\n",
"-6.55042989834691\n",
"-6.047933492978864\n"
"(-0.06089330922308864, array([-9.32966923e-04, -8.92579299e-07, -2.04258660e-09]), array([-5.95741551e-03, 7.27737654e-03, 6.90431275e-07, -8.11057566e-07,\n",
" -5.49031203e-09]))\n",
"(-0.060556512255197724, array([ 3.78506176e-03, -6.13838308e-03, 3.59037036e-06]), array([-3.45657964e-04, 6.59197575e-04, 7.10694456e-05, -6.23149871e-03,\n",
" -4.24978017e-05]))\n"
]
},
{
"ename": "IndexError",
"evalue": "list index out of range",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[9], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(pairs[\u001b[38;5;241m0\u001b[39m])) \u001b[38;5;66;03m# isotropic should be -82 meV\u001b[39;00m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(pairs[\u001b[38;5;241m1\u001b[39m])) \u001b[38;5;66;03m# these should all be around -41.9 in the isotropic part\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(\u001b[43mpairs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[43m]\u001b[49m)) \u001b[38;5;66;03m# these should all be around -41.9 in the isotropic part\u001b[39;00m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(pairs[\u001b[38;5;241m3\u001b[39m])) \u001b[38;5;66;03m# these should all be around -41.9 in the isotropic part\u001b[39;00m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(pairs[\u001b[38;5;241m4\u001b[39m])) \u001b[38;5;66;03m# these should all be around -41.9 in the isotropic part\u001b[39;00m\n",
"\u001b[0;31mIndexError\u001b[0m: list index out of range"
]
}
],
@ -456,6 +661,8 @@
}
],
"source": [
"These are reasonably converged:\n",
"\n",
"-61.33097171216109\n",
"-60.52198328932686\n",
"-60.51657719027764\n",
