grogu/test.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.14.3'"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import os\n",
    "from tqdm import tqdm\n",
    "from sys import getsizeof\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",
    "import pickle\n",
    "from numpy.linalg import inv\n",
    "import warnings\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": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================================================================================================================================================================\n",
      "Input file: \n",
      "Not yet specified.\n",
      "Output file: \n",
      "test.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:  1000\n",
      "k point directions:  x\n",
      "Ebot:  -13\n",
      "Eset:  500\n",
      "Esetp:  10000\n",
      "================================================================================================================================================================\n",
      "Setup done. Elapsed time: 5447.39860525 s\n",
      "================================================================================================================================================================\n"
     ]
    }
   ],
   "source": [
    "# this cell mimicks an input file\n",
    "fdf = sisl.get_sile(\n",
    "    \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\"\n",
    ")  # ./Jij_for_Marci_6p45ang/CrBr.fdf\n",
    "\n",
    "outfile = \"test\"\n",
    "if not outfile.endswith(\".pickle\"):\n",
    "    outfile += \".pickle\"\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",
    "\n",
    "# human readable definition of magnetic entities ./lat3_791/Fe3GeTe2.fdf\n",
    "magnetic_entities = [\n",
    "    dict(atom=3, l=2),\n",
    "    dict(atom=4, l=2),\n",
    "    dict(atom=5, l=2),\n",
    "]\n",
    "# pair information ./lat3_791/Fe3GeTe2.fdf\n",
    "pairs = [\n",
    "    # isotropic should be -82 meV\n",
    "    dict(ai=0, aj=1, 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([0, 0, 0])),\n",
    "]\n",
    "\n",
    "\"\"\"\n",
    "# human readable definition of magnetic entities ./Jij_for_Marci_6p45ang/CrBr.fdf\n",
    "magnetic_entities = [\n",
    "    dict(atom=0, l=2),\n",
    "    dict(atom=1, l=2),\n",
    "    dict(atom=2, l=2),\n",
    "]\n",
    "# pair information ./Jij_for_Marci_6p45ang/CrBr.fdf\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=1, Ruc=np.array([1, 0, 0])),\n",
    "    dict(ai=0, aj=2, Ruc=np.array([1, 0, 0])),\n",
    "    dict(ai=0, aj=1, Ruc=np.array([-1, 0, 0])),\n",
    "    dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])),\n",
    "    dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])),\n",
    "    dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])),\n",
    "    dict(ai=0, aj=1, Ruc=np.array([0, 1, 0])),\n",
    "    dict(ai=0, aj=2, Ruc=np.array([0, 1, 0])),\n",
    "]\n",
    "\"\"\"\n",
    "\n",
    "# Brilloun zone sampling and Green function contour integral\n",
    "kset = 1000\n",
    "kdirs = \"x\"\n",
    "ebot = -13\n",
    "eset = 500\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",
    "\n",
    "simulation_parameters = dict(\n",
    "    path=\"Not yet specified.\",\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 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": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-12.806739\n",
      "-0.01254111\n",
      "xyz[-3:]: red, green, blue\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/danielpozsar/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.9/site-packages/matplotlib/cbook.py:1762: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return math.isfinite(val)\n",
      "/Users/danielpozsar/Documents/oktatás/elte/phd/grogu_project/.venv/lib/python3.9/site-packages/matplotlib/cbook.py:1398: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return np.asarray(x, float)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABLwAAAGsCAYAAADXMb4GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3r0lEQVR4nO3deXxU9b3/8ffsM9lJgCwQdhAXXFFErNWKoLVWr3axtS1Vq+0VW5V7W6W32p9Wpdpby9Vaqf1Ze71XbeuvVm1rrUjdRUAQFUH2HZIAIZlkktnP74+ZTGayDEmY5Ewmr+fjwSMz55w5+c7XEU7e+Xw/x2IYhiEAAAAAAAAgR1jNHgAAAAAAAACQSQReAAAAAAAAyCkEXgAAAAAAAMgpBF4AAAAAAADIKQReAAAAAAAAyCkEXgAAAAAAAMgpBF4AAAAAAADIKXazB9BRNBrVvn37VFhYKIvFYvZwAADAIGEYhpqamlRVVSWrld/pZSOu8wAAQF/05Tov6wKvffv2qbq62uxhAACAQWr37t0aPXq02cNAF7jOAwAAR6M313lZF3gVFhZKir2JoqKijJ8/FArp5Zdf1pw5c+RwODJ+/sGO+UmP+UmP+UmP+UmP+UmP+UkvFArpueee07e+9a3EtQSyT39f5wEAgNzk9XpVXV3dq+u8rAu82srbi4qK+i3wysvLU1FRET8wdIH5SY/5SY/5SY/5SY/5SY/5Sa9tfiSxVC6L9fd1HgAAyG29uc6jwQUAAAAAAAByCoEXAAAAAAAAcgqBFwAAAAAAAHIKgRcAAAAAAAByCoEXAAAAAAAAcgqBFwAAAAAAAHIKgRcAAAAAAAByCoEXAAAAAAAAcgqBFwAAAAAAAHIKgRcAAAAAAAByCoEXAAAAAAAAcgqBFwAAAAAAAHIKgRcAAAAAAAByit3sAQAAgNxnGIaWbzskp82q6eNKzR4OctSWumbVNPo1pjRPY8ryzB4OAAAwERVeAACgX0Sjhmq9fv3fN7fp/J+/rq/+ZoX+8+WNZg8LOex372zX1x5boT+t2WP2UAAAgMmo8AIAABnhC4T10roaPf/BPm2qadLB5oDCUSOxP99p08QRBQpHorLb+J0bAAAA+g+BFwAAOGoPLdusX722Va2hSKd900YV6ytnjNHnT65SgYtLDwAAAPQ/rjoBAMBRe/TNbWoNRTR+eL7+5ZRROmfKCJUXuVSW75LTTjUXAAAABhaBFwAAOGr+eGXXk9+aoaoSj8mjAQAAwFDXq1+5RiIR3X777Ro/frw8Ho8mTpyon/zkJzKM9v4chmHojjvuUGVlpTwej2bPnq3NmzdnfOAAACA7RKOGQpHYtYCLai4AAABkgV5dld5333165JFH9Mtf/lIbNmzQfffdp/vvv18PPfRQ4pj7779fDz74oJYsWaIVK1YoPz9fc+fOld/vz/jgAQCA+YKRaOKxy2EzcSQAAABATK+WNL7zzju69NJLdfHFF0uSxo0bp6efflorV66UFKvuWrx4sX70ox/p0ksvlSQ98cQTKi8v13PPPacrr7yy0zkDgYACgUDiudfrlSSFQiGFQqG+vas02s7ZH+fOBcxPesxPesxPesxPesxPetk8P82t7WOyGhGFQkaao/tHNs4LAAAAzNOrwOuss87So48+qk2bNmnKlCn64IMP9NZbb+mBBx6QJG3fvl01NTWaPXt24jXFxcWaMWOGli9f3mXgtWjRIt15552dtr/88svKy8vr7fvpsaVLl/bbuXMB85Me85Me85Me85Me85NeNs5PY1CS7LLI0MsvvSSLxewRAQAAYKjrVeB12223yev1aurUqbLZbIpEIrrnnnt01VVXSZJqamokSeXl5SmvKy8vT+zraOHChVqwYEHiudfrVXV1tebMmaOioqJevZmeCIVCWrp0qS644AI5HI6Mn3+wY37SY37SY37SY37SY37Sy+b52XO4VVr9plwOmy6+eK4pYwiFQnr++edN+d4AAADIPr0KvP74xz/qySef1FNPPaXjjz9ea9eu1c0336yqqirNmzevTwNwuVxyuVydtjscjn69oO/v8w92zE96zE96zE96zE96zE962Tg/EcVaE7jstqwbGwAAAIamXgVe3//+93XbbbclliZOmzZNO3fu1KJFizRv3jxVVFRIkmpra1VZWZl4XW1trU4++eTMjRoAAGSNQDgiiTs0AgAAIHv06sq0paVFVmvqS2w2m6LR2N2Zxo8fr4qKCi1btiyx3+v1asWKFZo5c2YGhgsAALJNIBy7DnA5CLwAAACQHXpV4XXJJZfonnvu0ZgxY3T88cfr/fff1wMPPKBrrrlGkmSxWHTzzTfr7rvv1uTJkzV+/Hjdfvvtqqqq0mWXXdYf4wcAACYLhGKBl9NG4AUAAIDs0KvA66GHHtLtt9+uG264QXV1daqqqtK3v/1t3XHHHYljfvCDH8jn8+n6669XQ0ODzj77bL300ktyu90ZHzwAADBfMBKv8LLbTB4JAAAAENOrX8UWFhZq8eLF2rlzp1pbW7V161bdfffdcjqdiWMsFovuuusu1dTUyO/365VXXtGUKVMyPnAAAJAdAqF4Dy+WNA5qkUhEt99+u8aPHy+Px6OJEyfqJz/5iQzDSBxjGIbuuOMOVVZWyuPxaPbs2dq8ebOJowYAAOgaV6YAAOCoJHp40bR+ULvvvvv0yCOP6Je//KU2bNig++67T/fff78eeuihxDH333+/HnzwQS1ZskQrVqxQfn6+5s6dK7/fb+LIAQAAOuvVkkYAAICO2gIvJ0saB7V33nlHl156qS6++GJJ0rhx4/T0009r5cqVkmLVXYsXL9aPfvQjXXrppZKkJ554QuXl5XruuecSd/FOFggEFAgEEs+9Xu8AvBMAAAAqvAAAwFEKUuGVE8466ywtW7ZMmzZtkiR98MEHeuutt3TRRRdJkrZv366amhrNnj078Zri4mLNmDFDy5cv7/KcixYtUnFxceJPdXV1/78RAAAAUeEFAACOUiAc7+FF4DWo3XbbbfJ6vZo6dapsNpsikYjuueceXXXVVZKkmpoaSVJ5eXnK68rLyxP7Olq4cKEWLFiQeO71egm9AADAgCDwAgAAR6W9hxdLGgezP/7xj3ryySf11FNP6fjjj9fatWt18803q6qqSvPmzevTOV0ul1wuV4ZHCgAAcGQEXgAA4KgEQvHAi7s0Dmrf//73ddtttyV6cU2bNk07d+7UokWLNG/ePFVUVEiSamtrVVlZmXhdbW2tTj75ZDOGDAAA0C2uTAEAwFFpW9LotHFZMZi1tLTIak39b2iz2RSNxgLN8ePHq6KiQsuWLUvs93q9WrFihWbOnDmgYwUAADgSKrwAAMBRSTStp8JrULvkkkt0zz33aMyYMTr++OP1/vvv64EHHtA111wjSbJYLLr55pt19913a/LkyRo/frxuv/12VVVV6bLLLjN38AAAAB0QeAEAgKNCD6/c8NBDD+n222/XDTfcoLq6OlVVVenb3/627rjjjsQxP/jBD+Tz+XT99deroaFBZ599tl566SW53W4TRw4AANAZgRcAADgq3KUxNxQWFmrx4sVavHhxt8dYLBbddddduuuuuwZuYAAAAH3AlSkAADgq7RVeXFYAAAAgO3BlCgAAjkqQwAsAAABZhitTAABwVOjhBQAAgGxD4AUAAI5KoocXd2kEAABAluDKFAAAHJVAiCWNAAAAyC5cmQIAgKMSjMQCLyeBFwAAALIEV6YAAOCotFd40cMLAAAA2YHACwAAHJVEDy8qvAAAAJAluDIFAABHhbs0AgAAINsQeAEAgKPSFnjRwwsAAADZgitTAABwVIJh7tIIAACA7MKVKQAAOCqJHl4OLisAAACQHbgyBQAAfRaJGgpFDEn08AIAAED2IPACAAB91racUWJJIwAAALIHV6YAAKDPkgMvmtYDAAAgW3BlCgAA+qytf5fVItmtFpNHAwAAAMQQeAEAgD4LJO7QaJPFQuAFAACA7EDgBQAA+ow7NAIAACAbcXUKAAD6zB+KVXg5bVxSAAAAIHtwdQoAAPosGIkvaaTCCwAAAFmEq1M
      "text/plain": [
       "<Figure size 1500x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "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",
    "plt.figure(figsize=(15, 5))\n",
    "plt.subplot(121)\n",
    "plt.plot(np.sort(dh.eig()))\n",
    "plt.grid()\n",
    "plt.subplot(122)\n",
    "DOS = sisl.physics.electron.DOS(np.linspace(-15, 85, 50), dh.eig())\n",
    "plt.plot(DOS, np.linspace(-15, 85, 50))\n",
    "print(np.real(dh.eig()).min())\n",
    "print(np.imag(dh.eig()).min())\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": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Hamiltonian and exchange field rotated. Elapsed time: 5448.302713666 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",
    "# 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",
    "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": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Site and pair dictionaries created. Elapsed time: 5448.330635 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 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",
    "    # 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",
    "    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",
    "    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": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "k loop:   0%|          | 0/1000 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k set created. Elapsed time: 5448.339522166 s\n",
      "================================================================================================================================================================\n"
     ]
    }
   ],
   "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\")\n",
    "\n",
    "if rank == root_node:\n",
    "    times[\"k_set_time\"] = timer()\n",
    "    print(f\"k set created. Elapsed time: {times['k_set_time']} s\")\n",
    "    print(\n",
    "        \"================================================================================================================================================================\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rotations done perpendicular to quantization axis. Elapsed time: 5448.55170175 s\n",
      "================================================================================================================================================================\n"
     ]
    }
   ],
   "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 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",
    "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": 48,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting matrix inversions\n",
      "Total number of k points: 1000\n",
      "Number of energy samples per k point: 500\n",
      "Total number of directions: 3\n",
      "Total number of matrix inversions: 1500000\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: 750.0 KB\n",
      "Expected memory usage on root node: 732.421875 MB\n",
      "================================================================================================================================================================\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "k loop: 100%|██████████| 1000/1000 [39:28<00:00,  2.37s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculated Greens functions. Elapsed time: 7816.996438041 s\n",
      "================================================================================================================================================================\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\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 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",
    "        # 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",
    "        # 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",
    "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": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Magnetic entities integrated.\n",
      "Pairs integrated.\n",
      "Magnetic parameters calculated.\n",
      "##################################################################### GROGU OUTPUT #############################################################################\n",
      "================================================================================================================================================================\n",
      "Input file: \n",
      "Not yet specified.\n",
      "Output file: \n",
      "test.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:  1000\n",
      "k point directions:  x\n",
      "Ebot:  -13\n",
      "Eset:  500\n",
      "Esetp:  10000\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] Not yet.\n",
      "Isotropic:  -29.367092812942317\n",
      "DMI:  [-8.89792781e-01  4.16580540e+00 -3.59948517e-06]\n",
      "Symmetric-anisotropy:  [-4.56670448  1.17667365  0.33235705  4.15578143 -0.03089926]\n",
      "Energies for debugging: \n",
      "array([[-0.02819042, -0.00085889,  0.00092069, -0.02597706],\n",
      "       [-0.03690207, -0.00415578,  0.00417583, -0.0339338 ],\n",
      "       [-0.03160427, -0.00033236, -0.00033235, -0.03122101]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.0339338 , -0.03160427, -0.02819042])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.03393379729617979 -0.03122101431812066\n",
      "\n",
      "[4]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
      "Isotropic:  -74.98944973858333\n",
      "DMI:  [-6.82212886e+00  1.12891691e+01  5.77778279e-04]\n",
      "Symmetric-anisotropy:  [-1.03411109  0.63716653 -0.16031601 10.67952816 -0.32649278]\n",
      "Energies for debugging: \n",
      "array([[-0.07435228, -0.00649564,  0.00714862, -0.07459251],\n",
      "       [-0.07593353, -0.01067953,  0.01189881, -0.07602356],\n",
      "       [-0.0795701 ,  0.00016089,  0.00015974, -0.07975458]])\n",
      "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
      "array([-0.07602356, -0.0795701 , -0.07435228])\n",
      "Test J_xx = E(y,z) = E(z,y)\n",
      "-0.07602356083253432 -0.07975458355080561\n",
      "\n",
      "================================================================================================================================================================\n",
      "Runtime information: \n",
      "Total runtime: 2369.9501148750005 s\n",
      "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
      "Initial setup: 0.17246508400057792 s\n",
      "Hamiltonian conversion and XC field extraction: 0.904 s\n",
      "Pair and site datastructure creatrions: 0.028 s\n",
      "k set cration and distribution: 0.009 s\n",
      "Rotating XC potential: 0.212 s\n",
      "Greens function inversion: 2368.445 s\n",
      "Calculate energies and magnetic components: 0.180 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",
    "    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",
    "\n",
    "    print(\"Pairs integrated.\")\n",
    "\n",
    "    # calculate magnetic parameters\n",
    "    for pair in pairs:\n",
    "        J_iso, J_S, D = 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",
    "\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Magnetic entities integrated.\n",
    "Pairs integrated.\n",
    "Magnetic parameters calculated.\n",
    "##################################################################### GROGU OUTPUT #############################################################################\n",
    "================================================================================================================================================================\n",
    "Input file: \n",
    "Not yet specified.\n",
    "Output file: \n",
    "test.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:  1000\n",
    "k point directions:  x\n",
    "Ebot:  -13\n",
    "Eset:  500\n",
    "Esetp:  10000\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] Not yet.\n",
    "Isotropic:  -29.367092812942317\n",
    "DMI:  [-8.89792781e-01  4.16580540e+00 -3.59948517e-06]\n",
    "Symmetric-anisotropy:  [-4.56670448  1.17667365  0.33235705  4.15578143 -0.03089926]\n",
    "Energies for debugging: \n",
    "array([[-0.02819042, -0.00085889,  0.00092069, -0.02597706],\n",
    "       [-0.03690207, -0.00415578,  0.00417583, -0.0339338 ],\n",
    "       [-0.03160427, -0.00033236, -0.00033235, -0.03122101]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.0339338 , -0.03160427, -0.02819042])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.03393379729617979 -0.03122101431812066\n",
    "\n",
    "[4]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -74.98944973858333\n",
    "DMI:  [-6.82212886e+00  1.12891691e+01  5.77778279e-04]\n",
    "Symmetric-anisotropy:  [-1.03411109  0.63716653 -0.16031601 10.67952816 -0.32649278]\n",
    "Energies for debugging: \n",
    "array([[-0.07435228, -0.00649564,  0.00714862, -0.07459251],\n",
    "       [-0.07593353, -0.01067953,  0.01189881, -0.07602356],\n",
    "       [-0.0795701 ,  0.00016089,  0.00015974, -0.07975458]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.07602356, -0.0795701 , -0.07435228])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.07602356083253432 -0.07975458355080561\n",
    "\n",
    "================================================================================================================================================================\n",
    "Runtime information: \n",
    "Total runtime: 2369.9501148750005 s\n",
    "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
    "Initial setup: 0.17246508400057792 s\n",
    "Hamiltonian conversion and XC field extraction: 0.904 s\n",
    "Pair and site datastructure creatrions: 0.028 s\n",
    "k set cration and distribution: 0.009 s\n",
    "Rotating XC potential: 0.212 s\n",
    "Greens function inversion: 2368.445 s\n",
    "Calculate energies and magnetic components: 0.180 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (1614865498.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[50], line 1\u001b[0;36m\u001b[0m\n\u001b[0;31m    Magnetic entities integrated.\u001b[0m\n\u001b[0m             ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "Magnetic entities integrated.\n",
    "Pairs integrated.\n",
    "Magnetic parameters calculated.\n",
    "##################################################################### GROGU OUTPUT #############################################################################\n",
    "================================================================================================================================================================\n",
    "Input file: \n",
    "Not yet specified.\n",
    "Output file: \n",
    "test.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:  1000\n",
    "k point directions:  x\n",
    "Ebot:  -13\n",
    "Eset:  100\n",
    "Esetp:  10000\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] Not yet.\n",
    "Isotropic:  -29.36445059067917\n",
    "DMI:  [-8.91784089e-01  4.16579134e+00 -3.59901207e-06]\n",
    "Symmetric-anisotropy:  [-4.5685362   1.18301202  0.3323615   4.15576769 -0.02915828]\n",
    "Energies for debugging: \n",
    "array([[-0.02818144, -0.00086263,  0.00092094, -0.02597893],\n",
    "       [-0.03690122, -0.00415577,  0.00417581, -0.03393299],\n",
    "       [-0.03160383, -0.00033237, -0.00033236, -0.03122058]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.03393299, -0.03160383, -0.02818144])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.03393298678958339 -0.031220578858898954\n",
    "\n",
    "[4]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -74.99083333178646\n",
    "DMI:  [-6.82222254e+00  1.12891474e+01  5.95795483e-04]\n",
    "Symmetric-anisotropy:  [-1.03214008  0.64046775 -0.16031723 10.67951504 -0.32709547]\n",
    "Energies for debugging: \n",
    "array([[-0.07435037, -0.00649513,  0.00714932, -0.07459916],\n",
    "       [-0.07593293, -0.01067952,  0.01189878, -0.07602297],\n",
    "       [-0.0795695 ,  0.00016091,  0.00015972, -0.07975399]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.07602297, -0.0795695 , -0.07435037])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.07602297341278706 -0.07975398680085131\n",
    "\n",
    "================================================================================================================================================================\n",
    "Runtime information: \n",
    "Total runtime: 501.2763705420002 s\n",
    "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
    "Initial setup: 0.12531970899999578 s\n",
    "Hamiltonian conversion and XC field extraction: 0.888 s\n",
    "Pair and site datastructure creatrions: 0.062 s\n",
    "k set cration and distribution: 0.024 s\n",
    "Rotating XC potential: 0.245 s\n",
    "Greens function inversion: 499.875 s\n",
    "Calculate energies and magnetic components: 0.056 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (1097899960.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[40], line 1\u001b[0;36m\u001b[0m\n\u001b[0;31m    Magnetic entities integrated.\u001b[0m\n\u001b[0m             ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "Magnetic entities integrated.\n",
    "Pairs integrated.\n",
    "Magnetic parameters calculated.\n",
    "##################################################################### GROGU OUTPUT #############################################################################\n",
    "================================================================================================================================================================\n",
    "Input file: \n",
    "Not yet specified.\n",
    "Output file: \n",
    "test.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:  1000\n",
    "k point directions:  x\n",
    "Ebot:  -13\n",
    "Eset:  50\n",
    "Esetp:  10000\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] Not yet.\n",
    "Isotropic:  -29.359042196157713\n",
    "DMI:  [-8.93813379e-01  4.16574774e+00 -3.59758916e-06]\n",
    "Symmetric-anisotropy:  [-4.57147523  1.18770348  0.33237646  4.15572505 -0.02850636]\n",
    "Energies for debugging: \n",
    "array([[-0.02817134, -0.00086531,  0.00092232, -0.02597527],\n",
    "       [-0.03689861, -0.00415573,  0.00417577, -0.03393052],\n",
    "       [-0.03160252, -0.00033238, -0.00033237, -0.03121927]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.03393052, -0.03160252, -0.02817134])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.03393051742760642 -0.031219274782673133\n",
    "\n",
    "[4]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -74.98947011735306\n",
    "DMI:  [-6.82126017e+00  1.12890803e+01  6.51055432e-04]\n",
    "Symmetric-anisotropy:  [-1.03170277  0.64211178 -0.16032097 10.67947464 -0.32714492]\n",
    "Energies for debugging: \n",
    "array([[-0.07434736, -0.00649412,  0.00714841, -0.07459988],\n",
    "       [-0.07593108, -0.01067947,  0.01189869, -0.07602117],\n",
    "       [-0.07956766,  0.00016097,  0.00015967, -0.07975216]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.07602117, -0.07956766, -0.07434736])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.07602117288317486 -0.07975215648287909\n",
    "\n",
    "[3]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -75.04936117965987\n",
    "DMI:  [ 6.85701206e+00 -1.13051160e+01  6.39732773e-04]\n",
    "Symmetric-anisotropy:  [ -1.05318057   0.68103129  -0.16032    -10.70364791   0.31624841]\n",
    "Energies for debugging: \n",
    "array([[-0.07436833,  0.00654076, -0.00717326, -0.07467721],\n",
    "       [-0.07601601,  0.01070365, -0.01190658, -0.07610254],\n",
    "       [-0.07956796,  0.00016096,  0.00015968, -0.07975245]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.07610254, -0.07956796, -0.07436833])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.07610254174486064 -0.07975245001820208\n",
    "\n",
    "================================================================================================================================================================\n",
    "Runtime information: \n",
    "Total runtime: 261.890385875 s\n",
    "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
    "Initial setup: 0.1216103749998183 s\n",
    "Hamiltonian conversion and XC field extraction: 1.045 s\n",
    "Pair and site datastructure creatrions: 0.011 s\n",
    "k set cration and distribution: 0.006 s\n",
    "Rotating XC potential: 0.224 s\n",
    "Greens function inversion: 260.450 s\n",
    "Calculate energies and magnetic components: 0.033 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (1463926646.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[30], line 1\u001b[0;36m\u001b[0m\n\u001b[0;31m    Magnetic entities integrated.\u001b[0m\n\u001b[0m             ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "Magnetic entities integrated.\n",
    "Pairs integrated.\n",
    "Magnetic parameters calculated.\n",
    "##################################################################### GROGU OUTPUT #############################################################################\n",
    "================================================================================================================================================================\n",
    "Input file: \n",
    "Not yet specified.\n",
    "Output file: \n",
    "test.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:  3000\n",
    "k point directions:  x\n",
    "Ebot:  -13\n",
    "Eset:  50\n",
    "Esetp:  10000\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] Not yet.\n",
    "Isotropic:  -29.437239263932703\n",
    "DMI:  [-8.91096571e-01  4.16386016e+00 -3.68076797e-06]\n",
    "Symmetric-anisotropy:  [-4.68204753  1.22980695  0.33079224  4.15200152 -0.03180686]\n",
    "Energies for debugging: \n",
    "array([[-0.02820743, -0.00085929,  0.0009229 , -0.025985  ],\n",
    "       [-0.03707978, -0.004152  ,  0.00417572, -0.03411929],\n",
    "       [-0.03167345, -0.0003308 , -0.00033079, -0.03129215]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.03411929, -0.03167345, -0.02820743])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.03411928679618277 -0.03129215156448653\n",
    "\n",
    "[4]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -74.98086703681366\n",
    "DMI:  [-6.81158085e+00  1.13265508e+01  5.68487147e-04]\n",
    "Symmetric-anisotropy:  [-1.03582622  0.63764125 -0.15959679 10.72110977 -0.32326625]\n",
    "Energies for debugging: \n",
    "array([[-0.07434323, -0.00648831,  0.00713485, -0.07458268],\n",
    "       [-0.07592777, -0.01072111,  0.01193199, -0.07601669],\n",
    "       [-0.07959289,  0.00016017,  0.00015903, -0.07977653]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.07601669, -0.07959289, -0.07434323])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.07601669325607426 -0.07977653333847946\n",
    "\n",
    "================================================================================================================================================================\n",
    "Runtime information: \n",
    "Total runtime: 808.34989875 s\n",
    "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
    "Initial setup: 0.1594020829999998 s\n",
    "Hamiltonian conversion and XC field extraction: 0.995 s\n",
    "Pair and site datastructure creatrions: 0.021 s\n",
    "k set cration and distribution: 0.021 s\n",
    "Rotating XC potential: 0.232 s\n",
    "Greens function inversion: 806.884 s\n",
    "Calculate energies and magnetic components: 0.037 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Magnetic entities integrated.\n",
    "Pairs integrated.\n",
    "Magnetic parameters calculated.\n",
    "##################################################################### GROGU OUTPUT #############################################################################\n",
    "================================================================================================================================================================\n",
    "Input file: \n",
    "Not yet specified.\n",
    "Output file: \n",
    "test.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:  2000\n",
    "k point directions:  x\n",
    "Ebot:  -13\n",
    "Eset:  50\n",
    "Esetp:  10000\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] Not yet.\n",
    "Isotropic:  -29.313466612791018\n",
    "DMI:  [-8.85106686e-01  4.16843252e+00 -3.50602459e-06]\n",
    "Symmetric-anisotropy:  [-4.7530289   1.27630805  0.34469396  4.15760367 -0.03005581]\n",
    "Energies for debugging: \n",
    "array([[-0.02803716, -0.00085505,  0.00091516, -0.02583675],\n",
    "       [-0.03702114, -0.0041576 ,  0.00417926, -0.0340665 ],\n",
    "       [-0.03178903, -0.0003447 , -0.00034469, -0.03139132]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.0340665 , -0.03178903, -0.02803716])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.03406649551707385 -0.031391324924774755\n",
    "\n",
    "[4]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -74.97640334609967\n",
    "DMI:  [-6.84069197e+00  1.12809486e+01  8.75098125e-04]\n",
    "Symmetric-anisotropy:  [-1.03102826  0.63307456 -0.16123618 10.67285449 -0.32447457]\n",
    "Energies for debugging: \n",
    "array([[-0.07434333, -0.00651622,  0.00716517, -0.07457845],\n",
    "       [-0.07591065, -0.01067285,  0.01188904, -0.07600743],\n",
    "       [-0.07959609,  0.00016211,  0.00016036, -0.07978159]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.07600743, -0.07959609, -0.07434333])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.07600743160328334 -0.07978158717380361\n",
    "\n",
    "================================================================================================================================================================\n",
    "Runtime information: \n",
    "Total runtime: 534.5927785000001 s\n",
    "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
    "Initial setup: 0.14542879100008577 s\n",
    "Hamiltonian conversion and XC field extraction: 0.903 s\n",
    "Pair and site datastructure creatrions: 0.016 s\n",
    "k set cration and distribution: 0.014 s\n",
    "Rotating XC potential: 0.242 s\n",
    "Greens function inversion: 533.223 s\n",
    "Calculate energies and magnetic components: 0.050 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (577185453.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[10], line 1\u001b[0;36m\u001b[0m\n\u001b[0;31m    Magnetic entities integrated.\u001b[0m\n\u001b[0m             ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "Magnetic entities integrated.\n",
    "Pairs integrated.\n",
    "Magnetic parameters calculated.\n",
    "##################################################################### GROGU OUTPUT #############################################################################\n",
    "================================================================================================================================================================\n",
    "Input file: \n",
    "Not yet specified.\n",
    "Output file: \n",
    "test.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:  1000\n",
    "k point directions:  x\n",
    "Ebot:  -13\n",
    "Eset:  50\n",
    "Esetp:  10000\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] Not yet.\n",
    "Isotropic:  -29.359042196157713\n",
    "DMI:  [-8.93813379e-01  4.16574774e+00 -3.59758916e-06]\n",
    "Symmetric-anisotropy:  [-4.57147523  1.18770348  0.33237646  4.15572505 -0.02850636]\n",
    "Energies for debugging: \n",
    "array([[-0.02817134, -0.00086531,  0.00092232, -0.02597527],\n",
    "       [-0.03689861, -0.00415573,  0.00417577, -0.03393052],\n",
    "       [-0.03160252, -0.00033238, -0.00033237, -0.03121927]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.03393052, -0.03160252, -0.02817134])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.03393051742760642 -0.031219274782673133\n",
    "\n",
    "[4]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -74.98947011735306\n",
    "DMI:  [-6.82126017e+00  1.12890803e+01  6.51055432e-04]\n",
    "Symmetric-anisotropy:  [-1.03170277  0.64211178 -0.16032097 10.67947464 -0.32714492]\n",
    "Energies for debugging: \n",
    "array([[-0.07434736, -0.00649412,  0.00714841, -0.07459988],\n",
    "       [-0.07593108, -0.01067947,  0.01189869, -0.07602117],\n",
    "       [-0.07956766,  0.00016097,  0.00015967, -0.07975216]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.07602117, -0.07956766, -0.07434736])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.07602117288317486 -0.07975215648287909\n",
    "\n",
    "================================================================================================================================================================\n",
    "Runtime information: \n",
    "Total runtime: 272.760337666 s\n",
    "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
    "Initial setup: 0.13726933300000033 s\n",
    "Hamiltonian conversion and XC field extraction: 0.970 s\n",
    "Pair and site datastructure creatrions: 0.013 s\n",
    "k set cration and distribution: 0.017 s\n",
    "Rotating XC potential: 0.217 s\n",
    "Greens function inversion: 271.373 s\n",
    "Calculate energies and magnetic components: 0.034 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (2508083104.py, line 2)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  Cell \u001b[0;32mIn[10], line 2\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": [
    "##################################################################### GROGU OUTPUT #############################################################################\n",
    "================================================================================================================================================================\n",
    "Input file: \n",
    "Not yet specified.\n",
    "Output file: \n",
    "test.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:  20\n",
    "k point directions:  xy\n",
    "Ebot:  -15\n",
    "Eset:  100\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] Not yet.\n",
    "Isotropic:  -63.51146633376187\n",
    "DMI:  [-9.32966002e-01 -8.89800564e-04 -2.02862680e-06]\n",
    "Symmetric-anisotropy:  [-3.33721683e+00  1.29672703e+00  6.03687008e-04 -8.09695186e-04\n",
    " -6.00702131e-06]\n",
    "Energies for debugging: \n",
    "array([[-6.22147393e-02, -9.32959995e-04,  9.32972009e-04,\n",
    "        -6.14709765e-02],\n",
    "       [-6.75755955e-02,  8.09695186e-07, -9.69905943e-07,\n",
    "        -6.68486832e-02],\n",
    "       [-5.36624493e-02, -6.05715634e-07, -6.01658381e-07,\n",
    "        -5.36632256e-02]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.06684868, -0.05366245, -0.06221474])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.06684868316242108 -0.053663225627735525\n",
    "\n",
    "[4]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -60.97776480541008\n",
    "DMI:  [-3.79946502e+00  6.15248803e+00  3.52644581e-03]\n",
    "Symmetric-anisotropy:  [0.09164177 0.11139034 0.07107136 6.24190815 0.03636875]\n",
    "Energies for debugging: \n",
    "array([[-6.08663745e-02, -3.83583377e-03,  3.76309627e-03,\n",
    "        -6.11807969e-02],\n",
    "       [-6.06347990e-02, -6.24190815e-03,  6.06306790e-03,\n",
    "        -6.08861230e-02],\n",
    "       [-5.98977554e-02, -6.75449092e-05, -7.45978008e-05,\n",
    "        -5.98154343e-02]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.06088612, -0.05989776, -0.06086637])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.060886123032493376 -0.05981543426889457\n",
    "\n",
    "[3]Fe(2)        [5]Fe(2)        [0 0 0]        d [Ang] Not yet.\n",
    "Isotropic:  -60.9832595805832\n",
    "DMI:  [ 3.78505688e+00 -6.13845903e+00  3.52688440e-03]\n",
    "Symmetric-anisotropy:  [ 0.08120798  0.11261323  0.07107254 -6.23157554 -0.04249828]\n",
    "Energies for debugging: \n",
    "array([[-6.08706463e-02,  3.82755516e-03, -3.74255860e-03,\n",
    "        -6.11770808e-02],\n",
    "       [-6.06458648e-02,  6.23157554e-03, -6.04534251e-03,\n",
    "        -6.09020516e-02],\n",
    "       [-5.98977110e-02, -6.75456589e-05, -7.45994277e-05,\n",
    "        -5.98153841e-02]])\n",
    "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
    "array([-0.06090205, -0.05989771, -0.06087065])\n",
    "Test J_xx = E(y,z) = E(z,y)\n",
    "-0.060902051599344476 -0.05981538413498572\n",
    "\n",
    "================================================================================================================================================================\n",
    "Runtime information: \n",
    "Total runtime: 203.083672208 s\n",
    "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
    "Initial setup: 0.1483790410000001 s\n",
    "Hamiltonian conversion and XC field extraction: 0.861 s\n",
    "Pair and site datastructure creatrions: 0.017 s\n",
    "k set cration and distribution: 0.022 s\n",
    "Rotating XC potential: 0.247 s\n",
    "Greens function inversion: 201.715 s\n",
    "Calculate energies and magnetic components: 0.073 s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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"
   ]
  }
 ],
 "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
}