{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "[Daniels-MacBook-Air.local:67834] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Daniels-MacBook-Air.501/jf.0/2346057728/sm_segment.Daniels-MacBook-Air.501.8bd60000.0 could be created.\n" ] }, { "data": { "text/plain": [ "'0.14.3'" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from tqdm import tqdm\n", "from sys import getsizeof\n", "from timeit import default_timer as timer\n", "\n", "import numpy as np\n", "import sisl\n", "import sisl.viz\n", "from src.grogu_magn.useful import *\n", "from mpi4py import MPI\n", "import pickle\n", "from numpy.linalg import inv\n", "import warnings\n", "\n", "\"\"\" \n", "# Some input parsing\n", "parser = argparse.ArgumentParser()\n", "parser.add_argument('--kset' , dest = 'kset' , default = 2 , type=int , help = 'k-space resolution of Jij calculation')\n", "parser.add_argument('--kdirs' , dest = 'kdirs' , default = 'xyz' , help = 'Definition of k-space dimensionality')\n", "parser.add_argument('--eset' , dest = 'eset' , default = 42 , type=int , help = 'Number of energy points on the contour')\n", "parser.add_argument('--eset-p' , dest = 'esetp' , default = 10 , type=int , help = 'Parameter tuning the distribution on the contour')\n", "parser.add_argument('--input' , dest = 'infile' , required = True , help = 'Input file name')\n", "parser.add_argument('--output' , dest = 'outfile', required = True , help = 'Output file name')\n", "parser.add_argument('--Ebot' , dest = 'Ebot' , default = -20.0 , type=float, help = 'Bottom energy of the contour')\n", "parser.add_argument('--npairs' , dest = 'npairs' , default = 1 , type=int , help = 'Number of unitcell pairs in each direction for Jij calculation')\n", "parser.add_argument('--adirs' , dest = 'adirs' , default = False , help = 'Definition of pair directions')\n", "parser.add_argument('--use-tqdm', dest = 'usetqdm', default = 'not' , help = 'Use tqdm for progressbars or not')\n", "parser.add_argument('--cutoff' , dest = 'cutoff' , default = 100.0 , type=float, help = 'Real space cutoff for pair generation in Angs')\n", "parser.add_argument('--pairfile', dest = 'pairfile', default = False , help = 'File to read pair information')\n", "args = parser.parse_args()\n", "\"\"\"\n", "# runtime information\n", "times = dict()\n", "times[\"start_time\"] = timer()\n", "########################\n", "# it works if data is in downloads folder\n", "########################\n", "sisl.__version__" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "dat = sisl.io.siesta.eigSileSiesta(\n", " \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.EIG\"\n", ")\n", "siesta_eigs = dat.read_data()\n", "siesta_eigs.min()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ef = dat.read_fermi_level()\n", "Ef" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "================================================================================================================================================================\n", "Input file: \n", "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\n", "Output file: \n", "./Fe3GeTe2_notebook.pickle\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", "Parameters for the contour integral:\n", "Number of k points: 15\n", "k point directions: xy\n", "Ebot: -13\n", "Eset: 300\n", "Esetp: 1000\n", "================================================================================================================================================================\n", "Setup done. Elapsed time: 10.633095125 s\n", "================================================================================================================================================================\n" ] } ], "source": [ "################################################################################\n", "#################################### INPUT #####################################\n", "################################################################################\n", "path = (\n", " \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\"\n", ")\n", "outfile = \"./Fe3GeTe2_notebook\"\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 = [\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", "]\n", "magnetic_entities = [\n", " dict(atom=3, l=2),\n", " dict(atom=4, l=2),\n", " dict(atom=5, l=2),\n", "]\n", "pairs = [\n", " dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),\n", " dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])),\n", " dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])),\n", " dict(ai=0, aj=2, Ruc=np.array([-1, -1, 0])),\n", " dict(ai=1, aj=2, Ruc=np.array([-1, -1, 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([-2, 0, 0])),\n", " dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])),\n", "]\n", "# Brilloun zone sampling and Green function contour integral\n", "kset = 15\n", "kdirs = \"xy\"\n", "ebot = -13\n", "eset = 300\n", "esetp = 1000\n", "################################################################################\n", "#################################### INPUT #####################################\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", "\n", "# rename outfile\n", "if not outfile.endswith(\".pickle\"):\n", " outfile += \".pickle\"\n", "\n", "simulation_parameters = dict(\n", " path=path,\n", " outpath=outfile,\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", "\n", "# digestion of the input\n", "# read sile\n", "fdf = sisl.get_sile(path)\n", "# read in hamiltonian\n", "dh = fdf.read_hamiltonian()\n", "simulation_parameters[\"cell\"] = fdf.read_geometry().cell\n", "\n", "# unit cell index\n", "uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])\n", "\n", "if rank == root_node:\n", " print_parameters(simulation_parameters)\n", " times[\"setup_time\"] = timer()\n", " print(f\"Setup done. Elapsed time: {times['setup_time']} s\")\n", " print(\n", " \"================================================================================================================================================================\"\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "import matplotlib.pyplot as plt\n", "\n", "plt.figure(figsize=(15, 5))\n", "plt.subplot(121)\n", "plt.plot(np.sort(dh.eigh()), marker=\"o\", linestyle=\" \", label=\"sisl\")\n", "plt.plot(np.sort(siesta_eigs[0, 0]), marker=\"o\", linestyle=\" \", label=\"siesta\")\n", "plt.ylim(None, 10)\n", "plt.xlim(None, 75)\n", "plt.legend()\n", "plt.grid()\n", "plt.subplot(122)\n", "DOS = sisl.physics.electron.DOS(np.linspace(-15, 85, 1000), dh.eigh())\n", "plt.plot(DOS, np.linspace(-15, 85, 1000))\n", "DOS = sisl.physics.electron.DOS(np.linspace(-15, 85, 1000), siesta_eigs[0, 0])\n", "plt.plot(DOS, np.linspace(-15, 85, 1000))\n", "plt.ylim(None, 10)\n", "\n", "coords = dh.xyz[-3:]\n", "\n", "shift = np.array([-1, 0, 0]) @ simulation_parameters[\"cell\"]\n", "\n", "\n", "plt.figure(figsize=(15, 5))\n", "plt.subplot(131)\n", "plt.scatter(coords[:, 0], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n", "plt.scatter(\n", " (coords + shift)[:, 0], (coords + shift)[:, 2], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n", ")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"z\")\n", "plt.subplot(132)\n", "plt.scatter(coords[:, 1], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n", "plt.scatter(\n", " (coords + shift)[:, 1], (coords + shift)[:, 2], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n", ")\n", "plt.xlabel(\"y\")\n", "plt.ylabel(\"z\")\n", "plt.subplot(133)\n", "plt.scatter(coords[:, 0], coords[:, 1], color=[\"r\", \"g\", \"b\"])\n", "plt.scatter(\n", " (coords + shift)[:, 0], (coords + shift)[:, 1], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n", ")\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "print(\"xyz[-3:]: red, green, blue\")\n", "\n", "print(np.linalg.norm(coords[0] - coords[1]))\n", "print(np.linalg.norm(coords[0] - coords[2]))\n", "print(np.linalg.norm(coords[2] - coords[1]))\n", "print(np.linalg.norm(coords[0] - (coords + shift)[2]))" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Hamiltonian and exchange field rotated. Elapsed time: 11.290487791 s\n", "================================================================================================================================================================\n" ] } ], "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", "\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", "\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", "\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", "\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(\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 = 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", "for i in range(dh.lattice.sc_off.shape[0]):\n", " j = dh.lattice.sc_index(-dh.lattice.sc_off[i])\n", " h1, h1d = hh[i], hh[j]\n", " hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2\n", " s1, s1d = ss[i], ss[j]\n", " ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2\n", "\n", "\n", "###################################################################################\n", "# either this is shit\n", "\n", "# identifying TRS and TRB parts of the Hamiltonian\n", "TAUY = np.kron(np.eye(NO), tau_y)\n", "hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())])\n", "hTRS = (hh + hTR) / 2\n", "hTRB = (hh - hTR) / 2\n", "\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(\n", " [\n", " np.array([f[\"x\"] / 2 for f in traced]),\n", " np.array([f[\"y\"] / 2 for f in traced]),\n", " np.array([f[\"z\"] / 2 for f in traced]),\n", " ]\n", ") # equation 77\n", "\n", "###################################################################################\n", "\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(\n", " f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\"\n", " )\n", "\n", "if rank == root_node:\n", " times[\"H_and_XCF_time\"] = timer()\n", " print(\n", " f\"Hamiltonian and exchange field rotated. Elapsed time: {times['H_and_XCF_time']} s\"\n", " )\n", " print(\n", " \"================================================================================================================================================================\"\n", " )" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "R [[1. 0. 0.]\n", " [0. 1. 0.]\n", " [0. 0. 1.]]\n", "Z rot doesnt change anything in XCF: True\n", "MAX(myhk - sislhk) 2.974862582050264e-06 this is zero as should be\n", "MAX(myhk_rot - sislhk) 2.974862582050264e-06 rotation fucks it up\n", "MAX(rot_H - sislhk) 8.881784197001252e-16 rotation fucks it up\n" ] } ], "source": [ "myhk, sk = hsk(hh, ss, dh.sc_off, np.array([0.5, 0, 0]))\n", "sislhk = dh.Hk(np.array([0.5, 0, 0])).toarray()\n", "\n", "R = RotMa2b(scf_xcf_orientation, ref_xcf_orientations[2][\"o\"])\n", "print(\"R\", R)\n", "rot_XCF = np.einsum(\"ij,jklm->iklm\", R, XCF)\n", "print(\"Z rot doesnt change anything in XCF: \", np.allclose(XCF, rot_XCF))\n", "\n", "###################################################################################\n", "# either this is shit\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 = hTRS + rot_H_XCF # equation 76\n", "###################################################################################\n", "\n", "\n", "myhk_rot, sk = hsk(rot_H, ss, dh.sc_off, np.array([0.5, 0, 0]))\n", "\n", "print(\"MAX(myhk - sislhk)\", np.max(abs(myhk - sislhk)), \"this is zero as should be\")\n", "print(\"MAX(myhk_rot - sislhk)\", np.max(abs(myhk_rot - sislhk)), \"rotation fucks it up\")\n", "print(\"MAX(rot_H - sislhk)\", np.max(abs(rot_H - hh)), \"rotation fucks it up\")" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Site and pair dictionaries created. Elapsed time: 11.50304025 s\n", "================================================================================================================================================================\n" ] } ], "source": [ "# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n", "for mag_ent in magnetic_entities:\n", " parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes\n", " mag_ent[\"orbital_indeces\"] = parsed\n", " mag_ent[\"spin_box_indeces\"] = blow_up_orbindx(parsed) # calculate spin box indexes\n", " # if orbital is not set use all\n", " if \"l\" not in mag_ent.keys():\n", " mag_ent[\"l\"] = \"all\"\n", " if isinstance(mag_ent[\"atom\"], int):\n", " mag_ent[\"tags\"] = [\n", " f\"[{mag_ent['atom']}]{dh.atoms[mag_ent['atom']].tag}({mag_ent['l']})\"\n", " ]\n", " mag_ent[\"xyz\"] = [dh.xyz[mag_ent[\"atom\"]]]\n", " if isinstance(mag_ent[\"atom\"], list):\n", " mag_ent[\"tags\"] = []\n", " mag_ent[\"xyz\"] = []\n", " # iterate over atoms\n", " for atom_idx in mag_ent[\"atom\"]:\n", " mag_ent[\"tags\"].append(\n", " f\"[{atom_idx}]{dh.atoms[atom_idx].tag}({mag_ent['l']})\"\n", " )\n", " mag_ent[\"xyz\"].append(dh.xyz[atom_idx])\n", "\n", " # calculate size for Greens function generation\n", " spin_box_shape = len(mag_ent[\"spin_box_indeces\"])\n", "\n", " mag_ent[\"energies\"] = [] # we will store the second order energy derivations here\n", "\n", " # These will be the perturbed potentials from eq. 100\n", " mag_ent[\"Vu1\"] = [] # so they are independent in memory\n", " mag_ent[\"Vu2\"] = []\n", "\n", " mag_ent[\"Gii\"] = [] # Greens function\n", " mag_ent[\"Gii_tmp\"] = [] # Greens function for parallelization\n", " for i in ref_xcf_orientations:\n", " # Rotations for every quantization axis\n", " mag_ent[\"Vu1\"].append([])\n", " mag_ent[\"Vu2\"].append([])\n", " # Greens functions for every quantization axis\n", " mag_ent[\"Gii\"].append(\n", " np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n", " )\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", " # calculate distance\n", " xyz_ai = magnetic_entities[pair[\"ai\"]][\"xyz\"]\n", " xyz_aj = magnetic_entities[pair[\"aj\"]][\"xyz\"]\n", " xyz_aj = xyz_aj + pair[\"Ruc\"] @ simulation_parameters[\"cell\"]\n", " pair[\"dist\"] = np.linalg.norm(xyz_ai - xyz_aj)\n", "\n", " # calculate size for Greens function generation\n", " spin_box_shape_i = len(magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"])\n", " spin_box_shape_j = len(magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"])\n", " pair[\"tags\"] = []\n", " for mag_ent in [magnetic_entities[pair[\"ai\"]], magnetic_entities[pair[\"aj\"]]]:\n", " tag = \"\"\n", " # get atoms of magnetic entity\n", " atoms_idx = mag_ent[\"atom\"]\n", " orbitals = mag_ent[\"l\"]\n", "\n", " # if magnetic entity contains one atoms\n", " if isinstance(atoms_idx, int):\n", " tag += f\"[{atoms_idx}]{dh.atoms[atoms_idx].tag}({orbitals})\"\n", "\n", " # if magnetic entity contains more than one atoms\n", " if isinstance(atoms_idx, list):\n", " # iterate over atoms\n", " atom_group = \"{\"\n", " for atom_idx in atoms_idx:\n", " atom_group += f\"[{atom_idx}]{dh.atoms[atom_idx].tag}({orbitals})--\"\n", " # end {} of the atoms in the magnetic entity\n", " tag += atom_group[:-2] + \"}\"\n", " pair[\"tags\"].append(tag)\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", " for i in ref_xcf_orientations:\n", " # Greens functions for every quantization axis\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", " )\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", "if rank == root_node:\n", " times[\"site_and_pair_dictionaries_time\"] = timer()\n", " print(\n", " f\"Site and pair dictionaries created. Elapsed time: {times['site_and_pair_dictionaries_time']} s\"\n", " )\n", " print(\n", " \"================================================================================================================================================================\"\n", " )" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "k loop: 0%| | 0/225 [00:00iklm\", 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(\n", " dict(\n", " orient=orient[\"o\"],\n", " H=rot_H,\n", " GS=np.zeros((eset, rot_H.shape[1], rot_H.shape[2]), dtype=\"complex128\"),\n", " GS_tmp=np.zeros((eset, rot_H.shape[1], rot_H.shape[2]), dtype=\"complex128\"),\n", " )\n", " ) # store orientation and rotated Hamiltonian\n", "\n", " # these are the rotations (for now) perpendicular to the quantization axis\n", " for u in orient[\"vw\"]:\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", " idx = mag_ent[\"spin_box_indeces\"]\n", " # fill up the perturbed potentials (for now) based on the on-site projections\n", " mag_ent[\"Vu1\"][i].append(Vu1[:, idx][idx, :])\n", " mag_ent[\"Vu2\"][i].append(Vu2[:, idx][idx, :])\n", "\n", "if rank == root_node:\n", " times[\"reference_rotations_time\"] = timer()\n", " print(\n", " f\"Rotations done perpendicular to quantization axis. Elapsed time: {times['reference_rotations_time']} s\"\n", " )\n", " print(\n", " \"================================================================================================================================================================\"\n", " )" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Starting matrix inversions\n", "Total number of k points: 225\n", "Number of energy samples per k point: 300\n", "Total number of directions: 3\n", "Total number of matrix inversions: 202500\n", "The shape of the Hamiltonian and the Greens function is 84x84=7056\n", "Memory taken by a single Hamiltonian is: 0.015625 KB\n", "Expected memory usage per matrix inversion: 0.5 KB\n", "Expected memory usage per k point for parallel inversion: 450.0 KB\n", "Expected memory usage on root node: 98.876953125 MB\n", "================================================================================================================================================================\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "k loop: 100%|██████████| 225/225 [09:42<00:00, 2.59s/it]\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Calculated Greens functions. Elapsed time: 593.774521541 s\n", "================================================================================================================================================================\n" ] } ], "source": [ "if rank == root_node:\n", " print(\"Starting matrix inversions\")\n", " print(f\"Total number of k points: {kset.shape[0]}\")\n", " print(f\"Number of energy samples per k point: {eset}\")\n", " print(f\"Total number of directions: {len(hamiltonians)}\")\n", " print(\n", " f\"Total number of matrix inversions: {kset.shape[0] * len(hamiltonians) * eset}\"\n", " )\n", " print(f\"The shape of the Hamiltonian and the Greens function is {NO}x{NO}={NO*NO}\")\n", " # https://stackoverflow.com/questions/70746660/how-to-predict-memory-requirement-for-np-linalg-inv\n", " # memory is O(64 n**2) for complex matrices\n", " memory_size = getsizeof(hamiltonians[0][\"H\"].base) / 1024\n", " print(\n", " f\"Memory taken by a single Hamiltonian is: {getsizeof(hamiltonians[0]['H'].base) / 1024} KB\"\n", " )\n", " print(f\"Expected memory usage per matrix inversion: {memory_size * 32} KB\")\n", " print(\n", " f\"Expected memory usage per k point for parallel inversion: {memory_size * len(hamiltonians) * eset * 32} KB\"\n", " )\n", " print(\n", " f\"Expected memory usage on root node: {len(np.array_split(kset, size)[0]) * memory_size * len(hamiltonians) * eset * 32 / 1024} MB\"\n", " )\n", " print(\n", " \"================================================================================================================================================================\"\n", " )\n", "\n", "comm.Barrier()\n", "# ----------------------------------------------------------------------\n", "\n", "# make energy contour\n", "# we are working in eV now !\n", "# and sisl 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", " # iterate over reference directions\n", " for i, hamiltonian_orientation in enumerate(hamiltonians):\n", " # calculate Greens function\n", " H = hamiltonian_orientation[\"H\"]\n", " HK, SK = hsk(H, ss, dh.sc_off, k)\n", "\n", " # solve Greens function sequentially for the energies, because of memory bound\n", " Gk = np.zeros(shape=(eset, HK.shape[0], HK.shape[1]), dtype=\"complex128\")\n", " for j in range(eset):\n", " Gk[j] = inv(SK * eran[j] - HK)\n", "\n", " # saving this for total charge\n", " hamiltonian_orientation[\"GS_tmp\"] += Gk @ SK * wk\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] += (\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", " 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 total charge\n", " comm.Reduce(hamiltonians[i][\"GS_tmp\"], hamiltonians[i][\"GS\"], root=root_node)\n", "\n", " for mag_ent in magnetic_entities:\n", " comm.Reduce(mag_ent[\"Gii_tmp\"][i], mag_ent[\"Gii\"][i], root=root_node)\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", "\n", "if rank == root_node:\n", " times[\"green_function_inversion_time\"] = timer()\n", " print(\n", " f\"Calculated Greens functions. Elapsed time: {times['green_function_inversion_time']} s\"\n", " )\n", " print(\n", " \"================================================================================================================================================================\"\n", " )" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Total charge: 39.907801636668175\n", "Total charge: 39.90780156034552\n", "Total charge: 39.992244200462636\n", "Magnetic entities integrated.\n", "Pairs integrated.\n", "Magnetic parameters calculated.\n", "##################################################################### GROGU OUTPUT #############################################################################\n", "================================================================================================================================================================\n", "Input file: \n", "/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\n", "Output file: \n", "./Fe3GeTe2_notebook.pickle\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", "Parameters for the contour integral:\n", "Number of k points: 15\n", "k point directions: xy\n", "Ebot: -13\n", "Eset: 300\n", "Esetp: 1000\n", "================================================================================================================================================================\n", "Atomic information: \n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz\n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "[3]Fe(2) -7.339158738013707e-06 4.149278510690423e-06 11.657585837928032\n", "\n", "[4]Fe(2) -7.326987662162937e-06 4.158274523275774e-06 8.912422537596708\n", "\n", "[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n", "\n", "================================================================================================================================================================\n", "Exchange [meV]\n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] 2.745163300331324\n", "Isotropic: -87.36625034925267\n", "DMI: [ 3.21770364e-02 -1.00594555e-03 -3.51968374e-07]\n", "Symmetric-anisotropy: [ 2.12625625e+00 7.29980046e-05 -4.82064787e-05 7.29980046e-05\n", " 2.14479369e+00 -9.50239206e-06 -4.82064787e-05 -9.50239206e-06\n", " -4.27104994e+00]\n", "J: [-8.52399941e+01 7.29980046e-05 -4.82064787e-05 7.29980046e-05\n", " -8.52214567e+01 -9.50239206e-06 -4.82064787e-05 -9.50239206e-06\n", " -9.16373003e+01]\n", "Energies for debugging: \n", "array([[-9.16002498e-02, 3.21865388e-05, -3.21675340e-05,\n", " -9.17228673e-02],\n", " [-9.16743508e-02, 1.05415203e-06, -9.57739074e-07,\n", " -9.17598493e-02],\n", " [-7.87200461e-02, -7.33499730e-08, -7.26460363e-08,\n", " -7.87201389e-02]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.09175985, -0.07872005, -0.09160025])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.09175984926880505 -0.07872013892300818\n", "\n", "[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] 2.5835033632437767\n", "Isotropic: -41.52730925965426\n", "DMI: [ 1.17785417e+00 -2.06225136e+00 -4.69189233e-05]\n", "Symmetric-anisotropy: [ 0.02746399 0.15418338 -0.07648137 0.15418338 -0.1972997 -0.04062499\n", " -0.07648137 -0.04062499 0.16983571]\n", "J: [-4.14998453e+01 1.54183382e-01 -7.64813730e-02 1.54183382e-01\n", " -4.17246090e+01 -4.06249865e-02 -7.64813730e-02 -4.06249865e-02\n", " -4.13574736e+01]\n", "Energies for debugging: \n", "array([[-0.04145412, 0.00121848, -0.00113723, -0.04156529],\n", " [-0.04126082, 0.00213873, -0.00198577, -0.04129418],\n", " [-0.04188393, -0.00015423, -0.00015414, -0.04170551]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.04129418, -0.04188393, -0.04145412])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.04129418024412361 -0.041705510285739864\n", "\n", "[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] 2.583501767937866\n", "Isotropic: -41.532465661632486\n", "DMI: [-1.16839004e+00 2.04641313e+00 -4.67852015e-05]\n", "Symmetric-anisotropy: [ 0.02908454 0.15418348 0.06718218 0.15418348 -0.19939635 0.03271419\n", " 0.06718218 0.03271419 0.17031181]\n", "J: [-4.15033811e+01 1.54183477e-01 6.71821756e-02 1.54183477e-01\n", " -4.17318620e+01 3.27141946e-02 6.71821756e-02 3.27141946e-02\n", " -4.13621538e+01]\n", "Energies for debugging: \n", "array([[-0.0414585 , -0.0012011 , 0.00113568, -0.04157913],\n", " [-0.0412658 , -0.0021136 , 0.00197923, -0.04130058],\n", " [-0.0418846 , -0.00015423, -0.00015414, -0.04170618]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.04130058, -0.0418846 , -0.0414585 ])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.04130058163838069 -0.04170618060677052\n", "\n", "[3]Fe(2) [5]Fe(2) [-1 -1 0] d [Ang] 2.5834973202859075\n", "Isotropic: -41.516776604415945\n", "DMI: [-2.41730078e+00 8.52711009e-05 3.85154429e-05]\n", "Symmetric-anisotropy: [-3.24863466e-01 -1.57805287e-04 1.53264979e-04 -1.57805287e-04\n", " 1.46547000e-01 7.51504636e-02 1.53264979e-04 7.51504636e-02\n", " 1.78316466e-01]\n", "J: [-4.18416401e+01 -1.57805287e-04 1.53264979e-04 -1.57805287e-04\n", " -4.13702296e+01 7.51504636e-02 1.53264979e-04 7.51504636e-02\n", " -4.13384601e+01]\n", "Energies for debugging: \n", "array([[-4.11335642e-02, -2.49245124e-03, 2.34215032e-03,\n", " -4.11229167e-02],\n", " [-4.15433560e-02, -2.38536080e-07, -6.79938781e-08,\n", " -4.17094547e-02],\n", " [-4.16175425e-02, 1.96320730e-07, 1.19289844e-07,\n", " -4.19738255e-02]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.04170945, -0.04161754, -0.04113356])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.04170945467884016 -0.041973825462496055\n", "\n", "[4]Fe(2) [5]Fe(2) [-1 -1 0] d [Ang] 2.583495745338251\n", "Isotropic: -41.514139744541524\n", "DMI: [ 2.41735105e+00 -1.91201970e-04 3.71831534e-05]\n", "Symmetric-anisotropy: [-3.21261199e-01 -1.57863970e-04 1.40654575e-04 -1.57863970e-04\n", " 1.43245335e-01 -7.51485430e-02 1.40654575e-04 -7.51485430e-02\n", " 1.78015864e-01]\n", "J: [-4.18354009e+01 -1.57863970e-04 1.40654575e-04 -1.57863970e-04\n", " -4.13708944e+01 -7.51485430e-02 1.40654575e-04 -7.51485430e-02\n", " -4.13361239e+01]\n", "Energies for debugging: \n", "array([[-4.11342275e-02, 2.49249959e-03, -2.34220250e-03,\n", " -4.11235752e-02],\n", " [-4.15380203e-02, 5.05473958e-08, -3.31856545e-07,\n", " -4.16963053e-02],\n", " [-4.16182136e-02, 1.95047123e-07, 1.20680816e-07,\n", " -4.19744966e-02]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.04169631, -0.04161821, -0.04113423])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.04169630526113004 -0.04197449662603875\n", "\n", "[3]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] 2.583541444641373\n", "Isotropic: -41.518970343111015\n", "DMI: [1.16817703e+00 2.06220317e+00 7.21860779e-06]\n", "Symmetric-anisotropy: [ 0.03010538 -0.15396904 0.07631411 -0.15396904 -0.20084252 -0.03269774\n", " 0.07631411 -0.03269774 0.17073714]\n", "J: [-4.14888650e+01 -1.53969044e-01 7.63141143e-02 -1.53969044e-01\n", " -4.17198129e+01 -3.26977448e-02 7.63141143e-02 -3.26977448e-02\n", " -4.13482332e+01]\n", "Energies for debugging: \n", "array([[-0.04144688, 0.00120087, -0.00113548, -0.04156696],\n", " [-0.04124959, -0.00213852, 0.00198589, -0.04128287],\n", " [-0.04187266, 0.00015398, 0.00015396, -0.04169486]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.04128287, -0.04187266, -0.04144688])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.04128287386394945 -0.041694856070965985\n", "\n", "[4]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] 2.5835398672184064\n", "Isotropic: -41.51853887711384\n", "DMI: [-1.17770380e+00 -2.04619919e+00 9.44971395e-06]\n", "Symmetric-anisotropy: [ 0.02607311 -0.15396929 -0.0672924 -0.15396929 -0.19537854 0.0406109\n", " -0.0672924 0.0406109 0.16930543]\n", "J: [-4.14924658e+01 -1.53969287e-01 -6.72923993e-02 -1.53969287e-01\n", " -4.17139174e+01 4.06108984e-02 -6.72923993e-02 4.06108984e-02\n", " -4.13492334e+01]\n", "Energies for debugging: \n", "array([[-0.04144387, -0.00121831, 0.00113709, -0.04155449],\n", " [-0.04125459, 0.00211349, -0.00197891, -0.04128939],\n", " [-0.04187334, 0.00015398, 0.00015396, -0.04169554]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.04128939, -0.04187334, -0.04144387])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.041289392429914126 -0.04169553910557238\n", "\n", "[4]Fe(2) [5]Fe(2) [-2 0 0] d [Ang] 5.951322298958084\n", "Isotropic: -1.7091250393227355\n", "DMI: [ 0.03576957 0.26364426 -0.18258214]\n", "Symmetric-anisotropy: [ 0.06980213 0.03737695 0.02629686 0.03737695 -0.14019341 -0.03727655\n", " 0.02629686 -0.03727655 0.07039128]\n", "J: [-1.63932291 0.03737695 0.02629686 0.03737695 -1.84931845 -0.03727655\n", " 0.02629686 -0.03727655 -1.63873376]\n", "Energies for debugging: \n", "array([[-1.60572733e-03, 7.30461189e-05, 1.50698510e-06,\n", " -1.77828935e-03],\n", " [-1.67174019e-03, -2.89941114e-04, 2.37347403e-04,\n", " -1.62742800e-03],\n", " [-1.92034755e-03, -2.19959089e-04, 1.45205198e-04,\n", " -1.65121781e-03]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.00162743, -0.00192035, -0.00160573])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.0016274280031076719 -0.0016512178077259976\n", "\n", "[4]Fe(2) [5]Fe(2) [-3 0 0] d [Ang] 9.638732176310562\n", "Isotropic: -0.09189981541370944\n", "DMI: [ 0.00859892 0.0105977 -0.36714056]\n", "Symmetric-anisotropy: [-0.05448044 -0.03796364 -0.01728604 -0.03796364 0.04913009 -0.00640036\n", " -0.01728604 -0.00640036 0.00535035]\n", "J: [-0.14638025 -0.03796364 -0.01728604 -0.03796364 -0.04276972 -0.00640036\n", " -0.01728604 -0.00640036 -0.08654947]\n", "Energies for debugging: \n", "array([[-4.10410032e-05, 1.49992875e-05, -2.19856035e-06,\n", " -1.34816950e-05],\n", " [-1.32057931e-04, 6.68833577e-06, 2.78837344e-05,\n", " -1.85874459e-04],\n", " [-7.20577538e-05, -3.29176918e-04, 4.05104201e-04,\n", " -1.06886051e-04]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-1.85874459e-04, -7.20577538e-05, -4.10410032e-05])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.00018587445917353333 -0.00010688605065147818\n", "\n", "================================================================================================================================================================\n", "Runtime information: \n", "Total runtime: 583.827944292 s\n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "Initial setup: 0.10415358400000052 s\n", "Hamiltonian conversion and XC field extraction: 0.657 s\n", "Pair and site datastructure creatrions: 0.213 s\n", "k set cration and distribution: 0.023 s\n", "Rotating XC potential: 0.264 s\n", "Greens function inversion: 581.984 s\n", "Calculate energies and magnetic components: 0.582 s\n" ] } ], "source": [ "if rank == root_node:\n", " # Calculate total charge\n", " for hamiltonian in hamiltonians:\n", " GS = hamiltonian[\"GS\"]\n", " traced = np.trace((GS), axis1=1, axis2=2)\n", " integral = np.trapz(-1 / np.pi * np.imag(traced * cont.we))\n", " print(\"Total charge: \", integral)\n", "\n", " # iterate over the magnetic entities\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", " # 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", " # fill up the magnetic entities dictionary with the energies\n", " magnetic_entities[tracker][\"energies\"].append(storage)\n", " # convert to np array\n", " magnetic_entities[tracker][\"energies\"] = np.array(\n", " magnetic_entities[tracker][\"energies\"]\n", " )\n", " print(\"Magnetic entities integrated.\")\n", "\n", " # iterate over the pairs\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", " site_j = magnetic_entities[pair[\"aj\"]]\n", "\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", " # fill up the pairs dictionary with the energies\n", " pairs[tracker][\"energies\"].append(storage)\n", " # convert to np array\n", " pairs[tracker][\"energies\"] = np.array(pairs[tracker][\"energies\"])\n", "\n", " print(\"Pairs integrated.\")\n", "\n", " # calculate magnetic parameters\n", " for pair in pairs:\n", " J_iso, J_S, D, J = calculate_exchange_tensor(pair)\n", " pair[\"J_iso\"] = J_iso * sisl.unit_convert(\"eV\", \"meV\")\n", " pair[\"J_S\"] = J_S * sisl.unit_convert(\"eV\", \"meV\")\n", " pair[\"D\"] = D * sisl.unit_convert(\"eV\", \"meV\")\n", " pair[\"J\"] = J * sisl.unit_convert(\"eV\", \"meV\")\n", "\n", " print(\"Magnetic parameters calculated.\")\n", "\n", " times[\"end_time\"] = timer()\n", " print(\n", " \"##################################################################### GROGU OUTPUT #############################################################################\"\n", " )\n", "\n", " print_parameters(simulation_parameters)\n", " print_atoms_and_pairs(magnetic_entities, pairs)\n", " print_runtime_information(times)\n", "\n", " # remove clutter from magnetic entities and pair information\n", " for pair in pairs:\n", " del pair[\"Gij\"]\n", " del pair[\"Gij_tmp\"]\n", " del pair[\"Gji\"]\n", " del pair[\"Gji_tmp\"]\n", " for mag_ent in magnetic_entities:\n", " del mag_ent[\"Gii\"]\n", " del mag_ent[\"Gii_tmp\"]\n", " del mag_ent[\"Vu1\"]\n", " del mag_ent[\"Vu2\"]\n", " # create output dictionary with all the relevant data\n", " results = dict(\n", " parameters=simulation_parameters,\n", " magnetic_entities=magnetic_entities,\n", " pairs=pairs,\n", " runtime=times,\n", " )\n", " # save dictionary\n", " with open(outfile, \"wb\") as output_file:\n", " pickle.dump(results, output_file)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (3105939143.py, line 1)", "output_type": "error", "traceback": [ "\u001b[0;36m Cell \u001b[0;32mIn[13], line 1\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": [ "========================================\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": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": ".venv", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.6" } }, "nbformat": 4, "nbformat_minor": 2 }