grogu/test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Mac:54919] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Mac.501/jf.0/778436608/sm_segment.Mac.501.2e660000.0 could be created.\n"
     ]
    }
   ],
   "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 src.grogu_magn.useful import *\n",
    "from mpi4py import MPI\n",
    "from numpy.linalg import inv\n",
    "import warnings\n",
    "\n",
    "# runtime information\n",
    "times = dict()\n",
    "times[\"start_time\"] = timer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "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(\"./lat3_791/Fe3GeTe2.fdf\")\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",
    "\n",
    "# human readable definition of magnetic entities\n",
    "magnetic_entities = [\n",
    "    dict(atom=3, l=2),\n",
    "    dict(atom=4, l=2),\n",
    "    dict(atom=5, l=2),\n",
    "    dict(\n",
    "        atom=[3, 4],\n",
    "    ),\n",
    "]\n",
    "\n",
    "# pair information\n",
    "# these should all be around -41.9 in the isotropic part\n",
    "# isotropic should be -82 meV\n",
    "pairs = [\n",
    "    dict(ai=0, aj=3, Ruc=np.array([0, 0, 0])),\n",
    "    dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),\n",
    "    dict(ai=1, aj=0, 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, 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 = -30\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",
    "\n",
    "simulation_parameters = dict(\n",
    "    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",
    "\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",
    "times[\"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",
    "\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",
    "# 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\"] 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",
    "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",
    "times[\"H_and_XCF_time\"] = timer()"
   ]
  },
  {
   "cell_type": "code",
   "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",
    "    # calculate spin box indexes\n",
    "    magnetic_entities[i][\"spin_box_indeces\"] = blow_up_orbindx(parsed)\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",
    "    mag_ent[\"Gii\"] = []  # Greens function\n",
    "    mag_ent[\"Gii_tmp\"] = []  # Greens function for parallelization\n",
    "    # These will be the perturbed potentials from eq. 100\n",
    "    mag_ent[\"Vu1\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n",
    "    mag_ent[\"Vu2\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n",
    "    for i in ref_xcf_orientations:\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 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",
    "\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",
    "times[\"site_and_pair_dictionaries_time\"] = timer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "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[root_node] = tqdm(kpcs[root_node], desc=\"k loop\", file=stdout)\n",
    "\n",
    "times[\"k_set_time\"] = timer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this will contain the three hamiltonians in the reference directions 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",
    "    # 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(\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 = (\n",
    "        hTRS + rot_H_XCF\n",
    "    )  # equation 76 #######################################################################################\n",
    "\n",
    "    hamiltonians.append(\n",
    "        dict(orient=orient[\"o\"], H=rot_H)\n",
    "    )  # store orientation and rotated Hamiltonian\n",
    "\n",
    "    # these are the infinitezimal 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",
    "            # fill up the perturbed potentials (for now) based on the on-site projections\n",
    "            mag_ent[\"Vu1\"][i].append(\n",
    "                Vu1[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n",
    "            )\n",
    "            mag_ent[\"Vu2\"][i].append(\n",
    "                Vu2[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n",
    "            )\n",
    "\n",
    "times[\"reference_rotations_time\"] = timer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "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",
      "k loop: 100%|██████████| 400/400 [19:35<00:00,  2.94s/it]   \n"
     ]
    }
   ],
   "source": [
    "if rank == root_node:\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",
    "# 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",
    "    # 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",
    "        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] += (\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 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",
    "times[\"green_function_inversion_time\"] = timer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "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",
      "[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",
      "[3]Fe(all)        -7.339158738013707e-06        4.149278510690423e-06        11.657585837928032\n",
      "[4]Fe(all)        -7.326987662162937e-06        4.158274523275774e-06        8.912422537596708\n",
      "\n",
      "================================================================================================================================================================\n",
      "Exchange [meV]\n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "Magnetic entity1        Magnetic entity2       [i  j  k]       d [Ang]\n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "[3]Fe(2)        {[3]Fe(all)--[4]Fe(all)}             [0 0 0]    d [Ang] Not yet.\n",
      "Isotropic:  -7404.692892728154\n",
      "DMI:  [-9.31336133e-01 -7.86063119e-04  1.83575526e-06]\n",
      "Symmetric-anisotropy:  [ 2.03042884e+00  1.04458862e+00 -1.40962466e-03 -6.58426045e-04\n",
      " -6.97969537e-02]\n",
      "\n",
      "[3]Fe(2)        [4]Fe(2)             [0 0 0]    d [Ang] Not yet.\n",
      "Isotropic:  -60.89330922309355\n",
      "DMI:  [-9.32966923e-01 -8.92579299e-04 -2.04258659e-06]\n",
      "Symmetric-anisotropy:  [-5.95741551e+00  7.27737654e+00  6.90431275e-04 -8.11057566e-04\n",
      " -5.49031203e-06]\n",
      "\n",
      "[4]Fe(2)        [3]Fe(2)             [0 0 0]    d [Ang] Not yet.\n",
      "Isotropic:  -60.893309223093574\n",
      "DMI:  [9.32966923e-01 8.92579299e-04 2.04258658e-06]\n",
      "Symmetric-anisotropy:  [-5.95741551e+00  7.27737654e+00  6.90431275e-04  9.74101032e-04\n",
      " -5.49031203e-06]\n",
      "\n",
      "[3]Fe(2)        [5]Fe(2)             [0 0 0]    d [Ang] Not yet.\n",
      "Isotropic:  -60.55651225519789\n",
      "DMI:  [ 3.78506176e+00 -6.13838308e+00  3.59037036e-03]\n",
      "Symmetric-anisotropy:  [-0.34565796  0.65919757  0.07106945 -6.23149871 -0.0424978 ]\n",
      "\n",
      "[4]Fe(2)        [5]Fe(2)             [0 0 0]    d [Ang] Not yet.\n",
      "Isotropic:  -60.54974989595536\n",
      "DMI:  [-3.79945963e+00  6.15244494e+00  3.58990840e-03]\n",
      "Symmetric-anisotropy:  [-0.33638052  0.65239116  0.07106826  6.24185638  0.03636701]\n",
      "\n",
      "[3]Fe(2)        [5]Fe(2)             [-1  0  0]    d [Ang] Not yet.\n",
      "Isotropic:  -6.6253295899433216\n",
      "DMI:  [5.95251705 7.64859703 6.50501652]\n",
      "Symmetric-anisotropy:  [-0.65822877  0.72396528 -0.031302    7.69961304  0.03239586]\n",
      "\n",
      "[4]Fe(2)        [5]Fe(2)             [-1  0  0]    d [Ang] Not yet.\n",
      "Isotropic:  -6.123086646494101\n",
      "DMI:  [6.19414647 4.23019689 6.50504332]\n",
      "Symmetric-anisotropy:  [ 0.32693117  0.22187887 -0.03129943  4.24610256 -0.09833472]\n",
      "\n",
      "================================================================================================================================================================\n",
      "Runtime information: \n",
      "Total runtime: 118.107833833 s\n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "Initial setup: 0.12060429199999945 s\n",
      "Hamiltonian conversion and XC field extraction: 0.558 s\n",
      "Pair and site datastructure creatrions: 0.011 s\n",
      "k set cration and distribution: 0.018 s\n",
      "Rotating XC potential: 0.214 s\n",
      "Greens function inversion: 117.061 s\n",
      "Calculate energies and magnetic components: 0.127 s\n"
     ]
    }
   ],
   "source": [
    "if rank == root_node:\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",
    "\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",
    "\n",
    "    times[\"end_time\"] = timer()\n",
    "    print_output(simulation_parameters, magnetic_entities, pairs, dh, times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (3105939143.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[9], 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": [
    "import sisl.viz\n",
    "\n",
    "dh.geometry.tile(2, 1).plot(axes=\"xy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "xyz[-3:]: red, green, blue\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABNoAAAHACAYAAAB0/gUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCsElEQVR4nO3df5hWdZ0//ucMyAwWM2rK8GMnf5b4E0iTMFu1Rln1YmU/5aLtAktqW9l+VLYfUiqWJqlpbEayaUbq9vG35ie9IGUjM1lNlNZMTRSElEH9qjOICjJzf//g47QToAOcmXtmeDyu61x5n/u87/v1PkznNfdzzn1ORalUKgUAAAAA2CqV5S4AAAAAAHoDQRsAAAAAFEDQBgAAAAAFELQBAAAAQAEEbQAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEAB+pa7gO6otbU1zz//fAYMGJCKiopylwPQ45VKpaxatSpDhgxJZaW/8egzAMXSZzak1wAUq6O9RtC2Ec8//3zq6+vLXQZAr7N8+fL81V/9VbnLKDt9BqBz6DN/ptcAdI536zWCto0YMGBAkvU7r6ampszVAPR8zc3Nqa+vbzu+buv0GYBi6TMb0msAitXRXiNo24i3T62uqanRlAAK5Ksr6+kzAJ1Dn/kzvQagc7xbr3EBAwAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAAChA33IX0KusWJHMnp388Y/JgAHJCSckhx2WVFSUuzIAeoM330xuvDG59971veXww5NPfSqpri53ZQAA0C2tW5fccUcyd+76//7wh5N/+If1sU1nELQV5fLLkzPPTEqlpLLyz+uOOCK5/faktrac1QHQ0/32t8lxxyUvvpj0/X/t+6qrki99KbnrruRDHypvfQAA0M0sXpyMGZM888yff4X+8Y+TL385ueWW5Oiji39PXx0twq23Jv/7fyctLUlr6/qIdN269c/9+tfJ+PHlrQ+Anm3FiuSoo5KXX17/+H/2mZdeShoakhdeKF99AADQzbz+evLxjyfPPrv+8du/QpdK658bOzb5wx+Kf19B29YqlZJvfGPTXw9taVl/fuKiRV1aFgC9yKxZyWuvre8pf6mlJWlqWn92GwAAkCS54YZk+fKN/wrd2rp++bd/K/59BW1b67nnkv/+7/WB26b07bv+66MAsCVuumnjvyG8rbV1/TYAAECS9V8+rHyH1GvduvWXPy6aoG1rvf76u29TUZG88Ubn1wJA77R6dTHbAADANmL16vV/j34nb75Z/PsK2rZWfX3ynve88zZvvZXsv3/X1ANA7zNy5J+v3roxffuu3wYAAEiSDB/+zr9CV1Z2TlQjaNta/fsnJ5+c9Omz8ecrKpIddkhOOKFLywKgF/nCF/5884ONWbcu+fznu64eAADo5v75n9/5V+jW1uSLXyz+fQVtRfjmN5NhwzYM2/r0Wb/8x38k1dXlqQ2Anu+oo5LTTlv/3//z5jtvX3RiypTkiCO6vCwAAOiuhg1LLr54/X//z2u1VVSsX/7X/0r+8R+Lf19BWxFqa5Pf/CY566xkp53Wr6usXH+v2PvvT449trz1AdCzVVQkl1+e/PjHyb77/nn9/vsn116bfOc75asNAAC6qS9/OfnZz5KPfOTP63bfff3dRm+8cdNfTtwaFaXSO90uc9vU3Nyc2traNDU1paamZvMGt7YmTU3J9tsnVVWdUyBAD7NVx9VeaKv3R3Pz+vBtwIDiiwPogfSZDdknAO2tXr3+q6Q1Ne2/JNJRHT2uvsNl4dgilZXJjjuWuwoAejMfmAAAYLO8230si+KrowAAAABQAEEbAAAAABRA0AYAAAAABRC0AQAAAEABBG0AAAAAUABBGwAAAAAUQNAGAAAAAAUQtAEAAABAAQRtAAAAAFCAsgZt9957b8aOHZshQ4akoqIit99+e7vnb7311hx99NF53/vel4qKiixatOhdX3P27NmpqKhot1RXV3fOBADo1vQZAACgK5U1aFu9enWGDx+emTNnbvL5ww47LBdddNFmvW5NTU1WrFjRtjz77LNFlAtAD6PPAAAAXalvOd/8mGOOyTHHHLPJ5ydMmJAkWbp06Wa9bkVFRQYNGrQ1pQHQC+gzAABAV+qV12h77bXXsuuuu6a+vj7HH398HnvssXfcfs2aNWlubm63AMCm6DMAAMDG9Lqgbe+9987VV1+dn/3sZ7nuuuvS2tqaQw89NH/60582OWb69Ompra1tW+rr67uwYgB6En0GAADYlF4XtI0ePToTJ07MiBEjcvjhh+fWW2/NLrvskn//93/f5JipU6emqampbVm+fHkXVgxAT6LPAAAAm1LWa7R1he222y4jR47M4sWLN7lNVVVVqqqqurAqAHoLfQYAAHhbrzuj7S+1tLTk0UcfzeDBg8tdCgC9kD4DAAC8raxntL322mvtzgBYsmRJFi1alJ122invf//78/LLL2fZsmV5/vnnkyRPPvlkkmTQoEFtd3ubOHFihg4dmunTpydJvvnNb+YjH/lI9tprr7z66qu55JJL8uyzz+aUU07p4tkBUG76DAAA0JXKGrQ99NBDOfLII9seT5kyJUkyadKkzJ49O3fccUcmT57c9vyJJ56YJJk2bVrOO++8JMmyZctSWfnnE/NeeeWVnHrqqWlsbMyOO+6Ygw46KPfff3/23XffLpgRAN2JPgMAAHSlilKpVCp3Ed1Nc3Nzamtr09TUlJqamnKXA9DjOa62Z38AFKu7H1fvvffeXHLJJVm4cGFWrFiR2267LePGjXvHMfPnz8+UKVPy2GOPpb6+PmeffXb+6Z/+qcPv2d33CUBP09Hjaq+/RhsAAEA5rV69OsOHD8/MmTM7tP2SJUty3HHH5cgjj8yiRYtyxhln5JRTTsncuXM7uVIAtlavv+soAABAOR1zzDE55phjOrz9rFmzsvvuu+fSSy9Nkuyzzz6577778t3vfjdjxozprDIBKIAz2gAAALqRBQsWpKGhod26MWPGZMGCBZscs2bNmjQ3N7dbAOh6gjYAAIBupLGxMXV1de3W1dXVpbm5OW+88cZGx0yfPj21tbVtS319fVeUCsBfELQBAAD0cFOnTk1TU1Pbsnz58nKXBLBNco02AACAbmTQoEFZuXJlu3UrV65MTU1N+vfvv9ExVVVVqaqq6oryAHgHzmgDAADoRkaPHp158+a1W3f33Xdn9OjRZaoIgI4StAEAAHSi1157LYsWLcqiRYuSJEuWLMmiRYuybNmyJOu/9jlx4sS27T/3uc/lmWeeyVe+8pU88cQT+cEPfpAbb7wxZ555ZjnKB2AzCNoAAAA60UMPPZSRI0dm5MiRSZIpU6Zk5MiROffcc5MkK1asaAvdkmT33XfPnXfembvvvjvDhw/PpZdemquuuipjxowpS/0AdJxrtAEAAHSiI444IqVSaZPPz549e6NjHnnkkU6sCoDO4Iw2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAAgjaAAAAAKAAgjYAAAAAKICgDQAAAAAKIGgDAAAAgAII2gAAAACgAII2AAAAACiAoA0AAAAACiBoAwAAAIACCNoAAAAAoACCNgAAAAAogKANAAAAAAogaAMAAACAApQ1aLv33nszduzYDBkyJBUVFbn99tvbPX/rrbfm6KOPzvve975UVFRk0aJFHXrdm266KcOGDUt1dXUOOOCA3HXXXcUXD0C3p88AAABdqaxB2+rVqzN
      "text/plain": [
       "<Figure size 1500x500 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "coords = dh.xyz[-3:]\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.xlabel(\"x\")\n",
    "plt.ylabel(\"z\")\n",
    "plt.subplot(132)\n",
    "plt.scatter(coords[:, 1], coords[:, 2], color=[\"r\", \"g\", \"b\"])\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.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "print(\"xyz[-3:]: red, green, blue\")"
   ]
  },
  {
   "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
}