From e742c4fd187391eb5102aecac9a336b532c05016 Mon Sep 17 00:00:00 2001 From: Daniel Pozsar Date: Mon, 28 Oct 2024 20:13:42 +0100 Subject: [PATCH] something is still wrong with the isotropic exchange --- local_operator.ipynb | 228 ++++++++++++++ test.ipynb | 695 ++++++++++++++++++++++++++++--------------- test.py | 2 +- 3 files changed, 680 insertions(+), 245 deletions(-) create mode 100644 local_operator.ipynb diff --git a/local_operator.ipynb b/local_operator.ipynb new file mode 100644 index 0000000..968674a --- /dev/null +++ b/local_operator.ipynb @@ -0,0 +1,228 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/latex": [ + "$\\displaystyle \\operatorname{tr}\\left(\\left(0.5 \\left(AA G_{AA} + 0.5 AB G_{BA} + 0.5 AR G_{RA}\\right) AB + \\left(AA G_{AB} + 0.5 AB G_{BB} + 0.5 AR G_{RB}\\right) BB + 0.5 \\left(AA G_{AR} + 0.5 AB G_{BR} + 0.5 AR G_{RR}\\right) RB\\right) G_{BA} \\right) + \\operatorname{tr}\\left(\\left(0.25 BA G_{AA} AB + 0.5 BA G_{AB} BB + 0.25 BA G_{AR} RB\\right) G_{BB} \\right) + \\operatorname{tr}\\left(\\left(0.25 RA G_{AA} AB + 0.5 RA G_{AB} BB + 0.25 RA G_{AR} RB\\right) G_{BR} \\right) + 0.5 \\operatorname{tr}\\left(\\left(AA G_{AB} + 0.5 AB G_{BB} + 0.5 AR G_{RB}\\right) BA G_{AA} \\right) + 0.5 \\operatorname{tr}\\left(\\left(AA G_{AB} + 0.5 AB G_{BB} + 0.5 AR G_{RB}\\right) BR G_{RA} \\right) + 0.25 \\operatorname{tr}\\left(BA G_{AB} BA G_{AB} \\right) + 0.25 \\operatorname{tr}\\left(BA G_{AB} BR G_{RB} \\right) + 0.25 \\operatorname{tr}\\left(RA G_{AB} BA G_{AR} \\right) + 0.25 \\operatorname{tr}\\left(RA G_{AB} BR G_{RR} \\right)$" + ], + "text/plain": [ + "Trace((0.5*(AA*G_AA + 0.5*AB*G_BA + 0.5*AR*G_RA)*AB + (AA*G_AB + 0.5*AB*G_BB + 0.5*AR*G_RB)*BB + 0.5*(AA*G_AR + 0.5*AB*G_BR + 0.5*AR*G_RR)*RB)*G_BA) + Trace((0.25*BA*G_AA*AB + 0.5*BA*G_AB*BB + 0.25*BA*G_AR*RB)*G_BB) + Trace((0.25*RA*G_AA*AB + 0.5*RA*G_AB*BB + 0.25*RA*G_AR*RB)*G_BR) + 0.5*Trace((AA*G_AB + 0.5*AB*G_BB + 0.5*AR*G_RB)*BA*G_AA) + 0.5*Trace((AA*G_AB + 0.5*AB*G_BB + 0.5*AR*G_RB)*BR*G_RA) + 0.25*Trace(BA*G_AB*BA*G_AB) + 0.25*Trace(BA*G_AB*BR*G_RB) + 0.25*Trace(RA*G_AB*BA*G_AR) + 0.25*Trace(RA*G_AB*BR*G_RR)" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from sympy import MatrixSymbol, BlockMatrix, ZeroMatrix, trace, block_collapse\n", + "\n", + "# Define MatrixSymbols with specific sizes\n", + "# A = 2x2, B = 3x3, R = 4x4\n", + "AA = MatrixSymbol('AA', 2, 2) # 2x2\n", + "AB = MatrixSymbol('AB', 2, 3) # 2x3\n", + "BA = MatrixSymbol('BA', 3, 2) # 3x2\n", + "AR = MatrixSymbol('AR', 2, 4) # 2x4\n", + "RA = MatrixSymbol('RA', 4, 2) # 4x2\n", + "BB = MatrixSymbol('BB', 3, 3) # 3x3\n", + "BR = MatrixSymbol('BR', 3, 4) # 3x4\n", + "RB = MatrixSymbol('RB', 4, 3) # 4x3\n", + "RR = MatrixSymbol('RR', 4, 4) # 4x4\n", + "\n", + "# Create block matrices VA and VB using BlockMatrix\n", + "# We need to use ZeroMatrix with appropriate dimensions for padding\n", + "VA = BlockMatrix([\n", + " [AA, 0.5*AB, 0.5*AR],\n", + " [0.5*BA, ZeroMatrix(3,3), ZeroMatrix(3,4)],\n", + " [0.5*RA, ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", + "])\n", + "\n", + "VB = BlockMatrix([\n", + " [ZeroMatrix(2,2), 0.5*AB, ZeroMatrix(2,4)],\n", + " [0.5*BA, BB, 0.5*BR],\n", + " [ZeroMatrix(4,2), 0.5*RB, ZeroMatrix(4,4)]\n", + "])\n", + "\n", + "# Define G matrix symbols with matching dimensions\n", + "G_AA = MatrixSymbol('G_AA', 2, 2)\n", + "G_AB = MatrixSymbol('G_AB', 2, 3)\n", + "G_BA = MatrixSymbol('G_BA', 3, 2)\n", + "G_AR = MatrixSymbol('G_AR', 2, 4)\n", + "G_RA = MatrixSymbol('G_RA', 4, 2)\n", + "G_BB = MatrixSymbol('G_BB', 3, 3)\n", + "G_BR = MatrixSymbol('G_BR', 3, 4)\n", + "G_RB = MatrixSymbol('G_RB', 4, 3)\n", + "G_RR = MatrixSymbol('G_RR', 4, 4)\n", + "\n", + "# Create G matrix as a BlockMatrix\n", + "G = BlockMatrix([\n", + " [G_AA, G_AB, G_AR],\n", + " [G_BA, G_BB, G_BR],\n", + " [G_RA, G_RB, G_RR]\n", + "])\n", + "\n", + "\n", + "\n", + "# Calculate the trace of VA@G@VB@G\n", + "# First, let's calculate the product step by step\n", + "product = block_collapse(VA @ G @ VB @ G)\n", + "\n", + "# Calculate the trace\n", + "result = trace(product)\n", + "result" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Trace((0.5*(AA*G_AA + 0.5*AB*G_BA)*AB + (AA*G_AB + 0.5*AB*G_BB)*BB)*G_BA) + Trace((0.25*BA*G_AA*AB + 0.5*BA*G_AB*BB)*G_BB) + 0.5*Trace((AA*G_AB + 0.5*AB*G_BB)*BA*G_AA) + 0.25*Trace(BA*G_AB*BA*G_AB)\n" + ] + } + ], + "source": [ + "# Define MatrixSymbols with specific sizes\n", + "AA = MatrixSymbol('AA', 2, 2) # 2x2\n", + "AB = MatrixSymbol('AB', 2, 3) # 2x3\n", + "BA = MatrixSymbol('BA', 3, 2) # 3x2\n", + "BB = MatrixSymbol('BB', 3, 3) # 3x3\n", + "\n", + "# Create block matrices VA and VB using BlockMatrix\n", + "# R-related terms are replaced with zero matrices\n", + "VA = BlockMatrix([\n", + " [AA, 0.5*AB, ZeroMatrix(2,4)],\n", + " [0.5*BA, ZeroMatrix(3,3), ZeroMatrix(3,4)],\n", + " [ZeroMatrix(4,2), ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", + "])\n", + "\n", + "VB = BlockMatrix([\n", + " [ZeroMatrix(2,2), 0.5*AB, ZeroMatrix(2,4)],\n", + " [0.5*BA, BB, ZeroMatrix(3,4)],\n", + " [ZeroMatrix(4,2), ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", + "])\n", + "\n", + "# Define G matrix symbols with matching dimensions (kept original)\n", + "G_AA = MatrixSymbol('G_AA', 2, 2)\n", + "G_AB = MatrixSymbol('G_AB', 2, 3)\n", + "G_BA = MatrixSymbol('G_BA', 3, 2)\n", + "G_AR = MatrixSymbol('G_AR', 2, 4)\n", + "G_RA = MatrixSymbol('G_RA', 4, 2)\n", + "G_BB = MatrixSymbol('G_BB', 3, 3)\n", + "G_BR = MatrixSymbol('G_BR', 3, 4)\n", + "G_RB = MatrixSymbol('G_RB', 4, 3)\n", + "G_RR = MatrixSymbol('G_RR', 4, 4)\n", + "\n", + "# Create G matrix as a BlockMatrix (kept original)\n", + "G = BlockMatrix([\n", + " [G_AA, G_AB, G_AR],\n", + " [G_BA, G_BB, G_BR],\n", + " [G_RA, G_RB, G_RR]\n", + "])\n", + "\n", + "# Calculate the product and trace\n", + "product = block_collapse(VA @ G @ VB @ G)\n", + "result = trace(product)\n", + "print(result)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Trace(AA*G_AB*BB*G_BA)\n" + ] + } + ], + "source": [ + "from sympy import MatrixSymbol, BlockMatrix, ZeroMatrix, trace, block_collapse\n", + "\n", + "# Define MatrixSymbols with specific sizes\n", + "AA = MatrixSymbol('AA', 2, 2) # 2x2\n", + "BB = MatrixSymbol('BB', 3, 3) # 3x3\n", + "\n", + "# Create block matrices VA and VB using BlockMatrix with only AA and BB\n", + "VA = BlockMatrix([\n", + " [AA, ZeroMatrix(2,3), ZeroMatrix(2,4)],\n", + " [ZeroMatrix(3,2), ZeroMatrix(3,3), ZeroMatrix(3,4)],\n", + " [ZeroMatrix(4,2), ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", + "])\n", + "\n", + "VB = BlockMatrix([\n", + " [ZeroMatrix(2,2), ZeroMatrix(2,3), ZeroMatrix(2,4)],\n", + " [ZeroMatrix(3,2), BB, ZeroMatrix(3,4)],\n", + " [ZeroMatrix(4,2), ZeroMatrix(4,3), ZeroMatrix(4,4)]\n", + "])\n", + "\n", + "# Define G matrix symbols with matching dimensions\n", + "G_AA = MatrixSymbol('G_AA', 2, 2)\n", + "G_AB = MatrixSymbol('G_AB', 2, 3)\n", + "G_BA = MatrixSymbol('G_BA', 3, 2)\n", + "G_AR = MatrixSymbol('G_AR', 2, 4)\n", + "G_RA = MatrixSymbol('G_RA', 4, 2)\n", + "G_BB = MatrixSymbol('G_BB', 3, 3)\n", + "G_BR = MatrixSymbol('G_BR', 3, 4)\n", + "G_RB = MatrixSymbol('G_RB', 4, 3)\n", + "G_RR = MatrixSymbol('G_RR', 4, 4)\n", + "\n", + "# Create G matrix as a BlockMatrix\n", + "G = BlockMatrix([\n", + " [G_AA, G_AB, G_AR],\n", + " [G_BA, G_BB, G_BR],\n", + " [G_RA, G_RB, G_RR]\n", + "])\n", + "\n", + "# Calculate the product and simplify\n", + "product = block_collapse(VA @ G @ VB @ G)\n", + "\n", + "# Calculate the trace\n", + "result = trace(product)\n", + "simplified_result = result.simplify()\n", + "\n", + "print(simplified_result)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.6" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/test.ipynb b/test.ipynb index 70e6293..b64d79e 100644 --- a/test.ipynb +++ b/test.ipynb @@ -2,35 +2,29 @@ "cells": [ { "cell_type": "code", - "execution_count": 14, - "metadata": {}, - "outputs": [], - "source": [ - "import os\n", - "os.environ[\"OMP_NUM_THREADS\"] = \"4\" # export OMP_NUM_THREADS=4\n", - "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"4\" # export OPENBLAS_NUM_THREADS=4 \n", - "os.environ[\"MKL_NUM_THREADS\"] = \"4\" # export MKL_NUM_THREADS=6\n", - "os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"4\" # export VECLIB_MAXIMUM_THREADS=4\n", - "os.environ[\"NUMEXPR_NUM_THREADS\"] = \"4\"" - ] - }, - { - "cell_type": "code", - "execution_count": 15, + "execution_count": 1, "metadata": {}, "outputs": [ { - "data": { - "text/plain": [ - "'0.14.3'" - ] - }, - "execution_count": 15, - "metadata": {}, - "output_type": "execute_result" + "name": "stderr", + "output_type": "stream", + "text": [ + "[Daniels-MacBook-Air.local:09850] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Daniels-MacBook-Air.501/jf.0/320733184/sm_segment.Daniels-MacBook-Air.501.131e0000.0 could be created.\n" + ] } ], "source": [ + "import os\n", + "from sys import stdout\n", + "from tqdm import tqdm\n", + "from timeit import default_timer as timer\n", + "\n", + "os.environ[\"OMP_NUM_THREADS\"] = \"1\" # export OMP_NUM_THREADS=4\n", + "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"1\" # export OPENBLAS_NUM_THREADS=4\n", + "os.environ[\"MKL_NUM_THREADS\"] = \"1\" # export MKL_NUM_THREADS=6\n", + "os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"1\" # export VECLIB_MAXIMUM_THREADS=4\n", + "os.environ[\"NUMEXPR_NUM_THREADS\"] = \"1\" # export NUMEXPR_NUM_THREADS=6\n", + "\n", "import numpy as np\n", "import sisl\n", "from grogu.useful import *\n", @@ -38,12 +32,12 @@ "from numpy.linalg import inv\n", "import warnings\n", "\n", - "sisl.__version__\n" + "start_time = timer()\n" ] }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 2, "metadata": {}, "outputs": [ { @@ -56,34 +50,42 @@ ], "source": [ "# this cell mimicks an input file\n", - "fdf = sisl.get_sile('/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf')\n", + "fdf = sisl.get_sile(\n", + " \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\"\n", + ")\n", "# this information needs to be given at the input!!\n", - "scf_xcf_orientation=np.array([0,0,1]) #z\n", + "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", + "# at this moment the user has to supply o,v,w on the input.\n", "# we can have some default for this\n", - "ref_xcf_orientations=[dict(o=np.array([1,0,0]),vw=[np.array([0,1,0]),np.array([0,0,1])]),\n", - " 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", + "ref_xcf_orientations = [\n", + " dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),\n", + " dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),\n", + " dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),\n", + "]\n", "\n", "# human readable definition of magnetic entities\n", - "magnetic_entities=[dict(atom=3 ,l=2),\n", - " dict(atom=4 ,l=2),\n", - " dict(atom=5 ,l=2),\n", - " dict(atom=[3,4],)]\n", + "magnetic_entities = [\n", + " dict(atom=3, l=2),\n", + " dict(atom=4, l=2),\n", + " dict(atom=5, l=2),\n", + "# dict(atom=[3, 4]),\n", + "]\n", "\n", "# pair information\n", - "pairs=[dict(ai=0,aj=1,Ruc=np.array([0,0,0])), # isotropic should be -82 meV\n", - " dict(ai=0,aj=2,Ruc=np.array([0,0,0])), # these should all be around -41.9 in the isotropic part\n", - " dict(ai=1,aj=2,Ruc=np.array([0,0,0])),\n", - " dict(ai=0,aj=2,Ruc=np.array([-1,0,0])),\n", - " dict(ai=1,aj=2,Ruc=np.array([-1,0,0]))]\n", + "pairs = [\n", + " dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])), # isotropic should be -82 meV\n", + " dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])), # these should all be around -41.9 in the isotropic part\n", + "# dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])),\n", + "# dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])),\n", + "# dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])),\n", + "]\n", "\n", "# Brilloun zone sampling and Green function contour integral\n", "kset = 20\n", "kdirs = \"xy\"\n", - "ebot = -40\n", + "ebot = -30\n", "eset = 50\n", "esetp = 10000\n", "\n", @@ -94,254 +96,335 @@ "rank = comm.Get_rank()\n", "root_node = 0\n", "if rank == root_node:\n", - " print('Number of nodes in the parallel cluster: ',size)\n" + " print(\"Number of nodes in the parallel cluster: \", size)\n", + "\n", + "simulation_parameters = dict(path=\"Not yet specified.\",\n", + " scf_xcf_orientation=scf_xcf_orientation, \n", + " ref_xcf_orientations=ref_xcf_orientations,\n", + " kset=kset,\n", + " kdirs=kdirs, \n", + " ebot=ebot,\n", + " eset=eset, \n", + " esetp=esetp,\n", + " parallel_size=size)\n", + "\n", + "# digestion of the input\n", + "# read in hamiltonian\n", + "dh = fdf.read_hamiltonian()\n", + "try:\n", + " simulation_parameters[\"geom\"] = fdf.read_geometry()\n", + "except:\n", + " print(\"Error reading geometry.\")\n", + "\n", + "# unit cell index\n", + "uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])\n", + "\n", + "setup_time = timer()" ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": 17, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ - "# digestion of the input\n", - "# read in hamiltonian\n", - "dh = fdf.read_hamiltonian()\n", - "\n", - "# unit cell index\n", - "uc_in_sc_idx=dh.lattice.sc_index([0,0,0])\n", - "\n", - "\n", - "# WE WILL NOT NEED THIS!!\n", - "eigfile=sisl.io.siesta.eigSileSiesta('/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.EIG')\n", - "EF=eigfile.read_fermi_level()\n", - "\n", - "\n", - "NO=dh.no # shorthand for number of orbitals in the unit cell\n", + "NO = dh.no # shorthand for number of orbitals in the unit cell\n", "\n", "# preprocessing Hamiltonian and overlap matrix elements\n", - "h11 = dh.tocsr(dh.M11r)\n", - "h11 += dh.tocsr(dh.M11i)*1.0j\n", - "h11 = h11.toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\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", + "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", + "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", + "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 = dh.tocsr(dh.S_idx).toarray().reshape(NO,dh.n_s,NO).transpose(0,2,1).astype('complex128')\n", + "sov = (\n", + " dh.tocsr(dh.S_idx)\n", + " .toarray()\n", + " .reshape(NO, dh.n_s, NO)\n", + " .transpose(0, 2, 1)\n", + " .astype(\"complex128\")\n", + ")\n", "\n", "\n", "# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n", - "U=np.vstack([np.kron(np.eye(NO,dtype=int),[1,0]),np.kron(np.eye(NO,dtype=int),[0,1])])\n", + "U = np.vstack(\n", + " [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])]\n", + ")\n", "# This is the permutation that transforms ud1ud2 to u12d12\n", "# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n", "# the inverse transformation is U.T u12d12 to ud1ud2\n", "# That is FROM ORBITAL BOX to SPIN BOX => U.T\n", "\n", "# From now on everything is in SPIN BOX!!\n", - "hh,ss = \\\n", - "np.array([U.T@np.block([[h11[:,:,i],h12[:,:,i]],\n", - " [h21[:,:,i],h22[:,:,i]]])@U for i in range(dh.lattice.nsc.prod())]), \\\n", - "np.array([U.T@np.block([[sov[:,:,i] ,sov[:,:,i]*0],\n", - " [sov[:,:,i]*0,sov[:,:,i] ]])@U for i in range(dh.lattice.nsc.prod())])\n", - "\n", - "\n", - "# symmetrizing Hamiltonian and overlap matrix to make them hermitian \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", + " 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", + "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([np.array([f['x'] for f in traced]),\n", - " np.array([f['y'] for f in traced]),\n", - " np.array([f['z'] for f in traced])]) # equation 77\n", + "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(f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\")\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", + "H_and_XCF_time = timer()" ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n", "for i, mag_ent in enumerate(magnetic_entities):\n", - " parsed = parse_magnetic_entity(dh,**mag_ent) # parse orbital indexes\n", - " magnetic_entities[i]['orbital_indeces'] = parsed\n", - " magnetic_entities[i]['spin_box_indeces'] = blow_up_orbindx(parsed) # calculate spin box indexes\n", - " spin_box_shape = len(mag_ent[\"spin_box_indeces\"]) # calculate size for Greens function generation\n", - "\n", - " mag_ent['energies'] = [] # we will store the second order energy derivations here\n", - " \n", - " mag_ent['Gii'] = [] # Greens function\n", - " mag_ent['Gii_tmp'] = [] # Greens function for parallelization\n", - " mag_ent[\"Vu1\"] = [ list([]) for _ in range(len(ref_xcf_orientations))] # These will be the perturbed potentials from eq. 100\n", - " mag_ent[\"Vu2\"] = [ list([]) for _ in range(len(ref_xcf_orientations))]\n", + " parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes\n", + " magnetic_entities[i][\"orbital_indeces\"] = parsed\n", + " magnetic_entities[i][\"spin_box_indeces\"] = blow_up_orbindx(\n", + " parsed\n", + " ) # calculate spin box indexes\n", + " spin_box_shape = len(\n", + " mag_ent[\"spin_box_indeces\"]\n", + " ) # calculate size for Greens function generation\n", + "\n", + " mag_ent[\"energies\"] = [] # we will store the second order energy derivations here\n", + "\n", + " mag_ent[\"Gii\"] = [] # Greens function\n", + " mag_ent[\"Gii_tmp\"] = [] # Greens function for parallelization\n", + " mag_ent[\"Vu1\"] = [\n", + " list([]) for _ in range(len(ref_xcf_orientations))\n", + " ] # These will be the perturbed potentials from eq. 100\n", + " mag_ent[\"Vu2\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n", " for i in ref_xcf_orientations:\n", - " mag_ent['Gii'].append(np.zeros((eset, spin_box_shape, spin_box_shape),dtype='complex128')) # Greens functions for every quantization axis\n", - " mag_ent['Gii_tmp'].append(np.zeros((eset, spin_box_shape, spin_box_shape),dtype='complex128'))\n", - "\n", - "# for every site we have to store 2x3 Greens function (and the associated _tmp-s) in the 3 reference directions, because G_ij and G_ji are both needed\n", + " mag_ent[\"Gii\"].append(\n", + " np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n", + " ) # Greens functions for every quantization axis\n", + " mag_ent[\"Gii_tmp\"].append(\n", + " np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n", + " )\n", + "\n", + "# for every site we have to store 2x3 Greens function (and the associated _tmp-s)\n", + "# in the 3 reference directions, because G_ij and G_ji are both needed\n", "for pair in pairs:\n", - " spin_box_shape_i, spin_box_shape_j = len(magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]), len(magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]) # calculate size for Greens function generation\n", - " \n", - " pair['energies'] = [] # we will store the second order energy derivations here\n", - "\n", - " pair['Gij'] = [] # Greens function\n", - " pair['Gji'] = []\n", - " pair['Gij_tmp'] = [] # Greens function for parallelization\n", - " pair['Gji_tmp'] = []\n", - "\n", - " pair[\"Vij\"] = [list([]) for _ in range(len(ref_xcf_orientations))] # These will be the perturbed potentials from eq. 100\n", + " spin_box_shape_i, spin_box_shape_j = len(\n", + " magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]\n", + " ), len(\n", + " magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n", + " ) # calculate size for Greens function generation\n", + "\n", + " pair[\"energies\"] = [] # we will store the second order energy derivations here\n", + "\n", + " pair[\"Gij\"] = [] # Greens function\n", + " pair[\"Gji\"] = []\n", + " pair[\"Gij_tmp\"] = [] # Greens function for parallelization\n", + " pair[\"Gji_tmp\"] = []\n", + "\n", + " pair[\"Vij\"] = [\n", + " list([]) for _ in range(len(ref_xcf_orientations))\n", + " ] # These will be the perturbed potentials from eq. 100\n", " pair[\"Vji\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n", "\n", - " for i in ref_xcf_orientations: \n", - " pair['Gij'].append(np.zeros((eset,spin_box_shape_i,spin_box_shape_j),dtype='complex128'))\n", - " pair['Gij_tmp'].append(np.zeros((eset,spin_box_shape_i,spin_box_shape_j),dtype='complex128')) # Greens functions for every quantization axis\n", - " pair['Gji'].append(np.zeros((eset,spin_box_shape_j,spin_box_shape_i),dtype='complex128'))\n", - " pair['Gji_tmp'].append(np.zeros((eset,spin_box_shape_j,spin_box_shape_i),dtype='complex128'))\n" + " for i in ref_xcf_orientations:\n", + " pair[\"Gij\"].append(\n", + " np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n", + " )\n", + " pair[\"Gij_tmp\"].append(\n", + " np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n", + " ) # Greens functions for every quantization axis\n", + " pair[\"Gji\"].append(\n", + " np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n", + " )\n", + " pair[\"Gji_tmp\"].append(\n", + " np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n", + " )\n", + "\n", + "site_and_pair_dictionaries_time = timer()\n" ] }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "k loop: 0%| | 0/400 [00:00iklm',R, XCF)\n", - " rot_H_XCF = sum([np.kron(rot_XCF[i],tau) for i,tau in enumerate([tau_x,tau_y,tau_z])])\n", + "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", + "\n", " # obtain total Hamiltonian with the rotated exchange field\n", - " rot_H = hTRS+rot_H_XCF # equation 76\n", + " rot_H = hTRS + rot_H_XCF # equation 76\n", "\n", - " hamiltonians.append(dict(orient=orient['o'], H=rot_H, rotations=[])) # store orientation and rotated Hamiltonian\n", - " \n", - " for u in orient['vw']: # these are the infinitezimal rotations (for now) perpendicular to the quantization axis\n", - " Tu = np.kron(np.eye(NO,dtype=int),tau_u(u)) # section 2.H\n", + " hamiltonians.append(\n", + " dict(orient=orient[\"o\"], H=rot_H, rotations=[])\n", + " ) # store orientation and rotated Hamiltonian\n", + "\n", + " for u in orient[\n", + " \"vw\"\n", + " ]: # these are the infinitezimal rotations (for now) perpendicular to the quantization axis\n", + " Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H\n", "\n", - " Vu1 = 1j/2*commutator(rot_H_XCF_uc,Tu) # equation 100\n", - " Vu2 = 1/8*commutator(commutator(Tu,rot_H_XCF_uc),Tu) # equation 100\n", + " Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100\n", + " Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100\n", "\n", " for mag_ent in magnetic_entities:\n", - " mag_ent[\"Vu1\"][i].append(Vu1[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]) # fill up the perturbed potentials (for now) based on the on-site projections\n", - " mag_ent[\"Vu2\"][i].append(Vu2[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :])\n", + " mag_ent[\"Vu1\"][i].append(\n", + " Vu1[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n", + " ) # fill up the perturbed potentials (for now) based on the on-site projections\n", + " mag_ent[\"Vu2\"][i].append(\n", + " Vu2[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n", + " )\n", "\n", " for pair in pairs:\n", - " ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"] # get the pair orbital sizes from the magnetic entities\n", + " ai = magnetic_entities[pair[\"ai\"]][\n", + " \"spin_box_indeces\"\n", + " ] # get the pair orbital sizes from the magnetic entities\n", " aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n", - " pair[\"Vij\"][i].append(Vu1[:, ai][aj, :]) # fill up the perturbed potentials (for now) based on the on-site projections\n", - " pair[\"Vji\"][i].append(Vu1[:, aj][ai, :])\n" + " pair[\"Vij\"][i].append(\n", + " Vu1[:, ai][aj, :]\n", + " ) # fill up the perturbed potentials (for now) based on the on-site projections\n", + " pair[\"Vji\"][i].append(Vu1[:, aj][ai, :])\n", + "\n", + "reference_rotations_time = timer()" ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Number of magnetic entities being calculated: 4\n", + "Number of magnetic entities being calculated: 3\n", "We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.\n", - "The shape of the Hamiltonian and the Greens function is 84x84.\n" + "The shape of the Hamiltonian and the Greens function is 84x84.\n", + "k loop: 100%|██████████| 400/400 [03:05<00:00, 2.15it/s]\n" ] } ], "source": [ "if rank == root_node:\n", - " print('Number of magnetic entities being calculated: ',len(magnetic_entities))\n", - " print('We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.')\n", - " print(f'The shape of the Hamiltonian and the Greens function is {NO}x{NO}.')\n", + " print(\"Number of magnetic entities being calculated: \", len(magnetic_entities))\n", + " print(\n", + " \"We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.\"\n", + " )\n", + " print(f\"The shape of the Hamiltonian and the Greens function is {NO}x{NO}.\")\n", "comm.Barrier()\n", - "#----------------------------------------------------------------------\n", + "# ----------------------------------------------------------------------\n", "\n", - "# make energy contour \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", + "cont = make_contour(emin=ebot, enum=eset, p=esetp)\n", "eran = cont.ze\n", "\n", - "#---------------------------------------------------------------------- \n", + "# ----------------------------------------------------------------------\n", "# sampling the integrand on the contour and the BZ\n", "for k in kpcs[rank]:\n", - " wk = wkset[rank] # weight of k point in BZ integral\n", - " for i, hamiltonian_orientation in enumerate(hamiltonians): # iterate over reference directions\n", - " # calculate Greens function\n", + " wk = wkset[rank] # weight of k point in BZ integral\n", + " for i, hamiltonian_orientation in enumerate(\n", + " hamiltonians\n", + " ): # iterate over reference directions\n", + " # calculate Greens function\n", " H = hamiltonian_orientation[\"H\"]\n", - " HK,SK = hsk(H, ss, dh.sc_off, k)\n", + " HK, SK = hsk(H, ss, dh.sc_off, k)\n", " Gk = inv(SK * eran.reshape(eset, 1, 1) - HK)\n", - " \n", + "\n", " # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n", " for mag_ent in magnetic_entities:\n", - " mag_ent[\"Gii_tmp\"][i] += Gk[:,mag_ent[\"spin_box_indeces\"]][...,mag_ent[\"spin_box_indeces\"]] * wk\n", + " 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", + " 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", + " pair[\"Gij_tmp\"][i] += Gk[:, ai][..., aj] * phase * wk\n", + " pair[\"Gji_tmp\"][i] += Gk[:, aj][..., ai] * phase * wk\n", "\n", "# summ reduce partial results of mpi nodes\n", "for i in range(len(hamiltonians)):\n", @@ -350,84 +433,206 @@ "\n", " for pair in pairs:\n", " comm.Reduce(pair[\"Gij_tmp\"][i], pair[\"Gij\"][i], root=root_node)\n", - " comm.Reduce(pair[\"Gji_tmp\"][i], pair[\"Gji\"][i], root=root_node)\n" + " comm.Reduce(pair[\"Gji_tmp\"][i], pair[\"Gji\"][i], root=root_node)\n", + "\n", + "green_function_inversion_time = timer()\n" ] }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 8, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "############################### GROGU OUTPUT ###################################\n", + "================================================================================\n", + "Input file: \n", + "Not yet specified.\n", + "Number of nodes in the parallel cluster: 1\n", + "================================================================================\n", + "Cell [Ang]: \n", + "[[ 3.79100000e+00 0.00000000e+00 0.00000000e+00]\n", + " [-1.89550000e+00 3.28310231e+00 0.00000000e+00]\n", + " [ 1.25954923e-15 2.18160327e-15 2.05700000e+01]]\n", + "================================================================================\n", + "DFT axis: \n", + "[0 0 1]\n", + "Quantization axis and perpendicular rotation directions:\n", + "[1 0 0] --» [array([0, 1, 0]), array([0, 0, 1])]\n", + "[0 1 0] --» [array([1, 0, 0]), array([0, 0, 1])]\n", + "[0 0 1] --» [array([1, 0, 0]), array([0, 1, 0])]\n", + "================================================================================\n", + "number of k points: 20\n", + "k point directions: xy\n", + "================================================================================\n", + "Parameters for the contour integral:\n", + "Ebot: -30\n", + "Eset: 50\n", + "Esetp: 10000\n", + "================================================================================\n", + "Atomic informations: \n", + "\n", + "\n", + "Not yet specified.\n", + "\n", + "\n", + "================================================================================\n", + "Exchange [meV]\n", + "--------------------------------------------------------------------------------\n", + "Atom1 Atom2 [i j k] d [Ang]\n", + "--------------------------------------------------------------------------------\n", + "0 1 [0 0 0] distance\n", + "Isotropic: -60.893309223088636\n", + "DMI: [-5.95741551e+00 7.27737654e+00 6.90431275e-04 -8.11057566e-04\n", + " -5.49031203e-06]\n", + "Symmetric-anisotropy: [-9.32966923e-01 -8.92579299e-04 -2.04258660e-06]\n", + "\n", + "0 2 [0 0 0] distance\n", + "Isotropic: -60.556512255197724\n", + "DMI: [-0.34565796 0.65919757 0.07106945 -6.23149871 -0.0424978 ]\n", + "Symmetric-anisotropy: [ 3.78506176e+00 -6.13838308e+00 3.59037036e-03]\n", + "\n", + "================================================================================\n", + "Runtime information: \n", + "Total runtime: 187.525066541\n", + "--------------------------------------------------------------------------------\n", + "Initial setup: 0.275851458\n", + "Hamiltonian conversion and XC field extraction: 1.232 s\n", + "Pair and site datastructure creatrions: 0.025 s\n", + "k set cration and distribution: 0.029 s\n", + "Rotating XC potential: 0.366 s\n", + "Greens function inversion: 185.537 s\n", + "Calculate energies and magnetic components: 0.060 s\n" + ] + } + ], "source": [ - "# split magnetic entities and pairs over MPI nodes\n", - "mag_ent_parallel_set = np.array_split(magnetic_entities,size)\n", - "pairs_parallel_set = np.array_split(pairs,size)\n", - "\n", - "# iterate over the magnetic entities\n", - "for tracker, mag_ent in enumerate(mag_ent_parallel_set[rank]):\n", - " # iterate over the quantization axes\n", - " for i, Gii in enumerate(mag_ent[\"Gii\"]):\n", - " storage = []\n", - " # iterate over the first and second order local perturbations\n", - " for Vu1, Vu2 in zip(mag_ent[\"Vu1\"][i], mag_ent[\"Vu2\"][i]):\n", - " # The Szunyogh-Lichtenstein formula\n", - " traced = np.trace((Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii),axis1=1,axis2=2)\n", - " # evaluation of the contour integral\n", - " storage.append(np.trapz(-1/np.pi * np.imag(traced * cont.we)))\n", - "\n", - " # fill up the magnetic entities dictionary with the energies\n", - " idx = np.array([len(mag_ent_parallel_set[i]) for i in range(rank)]).sum() + tracker\n", - " magnetic_entities[int(idx)][\"energies\"].append(storage)\n", - "\n", - "# iterate over the pairs\n", - "for tracker, pair in enumerate(pairs_parallel_set[rank]):\n", - " # iterate over the quantization axes\n", - " for i, (Gij, Gji) in enumerate(zip(pair[\"Gij\"], pair[\"Gji\"])):\n", - " site_i = magnetic_entities[pair[\"ai\"]]\n", - " site_j = magnetic_entities[pair[\"aj\"]]\n", - "\n", - " 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", + "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((Vui @ Gij @ Vuj @ Gji),axis1=1,axis2=2)\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 pairs dictionary with the energies\n", - " idx = np.array([len(pairs_parallel_set[i]) for i in range(rank)]).sum() + tracker\n", - " pairs[int(idx)][\"energies\"].append(storage)\n" - ] - }, - { - "cell_type": "code", - "execution_count": 23, - "metadata": {}, - "outputs": [], - "source": [ - "def calculate_exchange_tensor(pair):\n", - " Eo, Ev, Ew = pair[\"energies\"]\n", - " J = np.array([Ew[-1], Ev[-1], Eo[0]]) # xx, yy, zz\n", - " JS = -0.5 * np.array([Eo[1] + Eo[2], Ev[1] + Ev[2], Ew[1] + Ew[2]]) # yz, zx, xy\n", - " D = 0.5 * np.array([Eo[1] - Eo[2], Ev[2] - Ev[1], Ew[1] - Ew[2]]) # x, y, z\n", - " return J.sum()/3 * 1000" + " storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n", + "\n", + " # fill up the magnetic entities dictionary with the energies\n", + " mag_ent[\"energies\"].append(storage)\n", + "\n", + " # iterate over the pairs\n", + " for tracker, pair in enumerate(pairs):\n", + " # iterate over the quantization axes\n", + " for i, (Gij, Gji) in enumerate(zip(pair[\"Gij\"], pair[\"Gji\"])):\n", + " site_i = magnetic_entities[pair[\"ai\"]]\n", + " site_j = magnetic_entities[pair[\"aj\"]]\n", + "\n", + " storage = []\n", + " # iterate over the first order local perturbations in all possible orientations for the two sites\n", + " for Vui in site_i[\"Vu1\"][i]:\n", + " for Vuj in site_j[\"Vu1\"][i]:\n", + " # The Szunyogh-Lichtenstein formula\n", + " traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2)\n", + " # evaluation of the contour integral\n", + " storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n", + "\n", + " # fill up the pairs dictionary with the energies\n", + " pairs[tracker][\"energies\"].append(storage)\n", + "\n", + " end_time = timer()\n", + "\n", + " print(\"############################### GROGU OUTPUT ###################################\")\n", + " print(\"================================================================================\")\n", + " print(\"Input file: \")\n", + " print(simulation_parameters[\"path\"])\n", + " print(\"Number of nodes in the parallel cluster: \", simulation_parameters[\"parallel_size\"])\n", + " print(\"================================================================================\")\n", + " try:\n", + " print(\"Cell [Ang]: \")\n", + " print(simulation_parameters[\"geom\"].cell)\n", + " except:\n", + " print(\"Geometry could not be read.\")\n", + " print(\"================================================================================\")\n", + " print(\"DFT axis: \")\n", + " print(simulation_parameters[\"scf_xcf_orientation\"])\n", + " print(\"Quantization axis and perpendicular rotation directions:\")\n", + " for ref in ref_xcf_orientations:\n", + " print(ref[\"o\"], \" --» \", ref[\"vw\"])\n", + " print(\"================================================================================\")\n", + " print(\"number of k points: \", simulation_parameters[\"kset\"])\n", + " print(\"k point directions: \", simulation_parameters[\"kdirs\"])\n", + " print(\"================================================================================\")\n", + " print(\"Parameters for the contour integral:\")\n", + " print(\"Ebot: \", simulation_parameters[\"ebot\"])\n", + " print(\"Eset: \", simulation_parameters[\"eset\"])\n", + " print(\"Esetp: \", simulation_parameters[\"esetp\"])\n", + " print(\"================================================================================\")\n", + " print(\"Atomic informations: \")\n", + " print(\"\")\n", + " print(\"\")\n", + " print(\"Not yet specified.\")\n", + " print(\"\")\n", + " print(\"\")\n", + " print(\"================================================================================\")\n", + " print(\"Exchange [meV]\")\n", + " print(\"--------------------------------------------------------------------------------\")\n", + " print(\"Atom1 Atom2 [i j k] d [Ang]\")\n", + " print(\"--------------------------------------------------------------------------------\")\n", + " for pair in pairs:\n", + " J_iso, J_S, D = calculate_exchange_tensor(pair)\n", + " J_iso = J_iso * sisl.unit_convert(\"eV\", \"meV\")\n", + " J_S = J_S * sisl.unit_convert(\"eV\", \"meV\")\n", + " D = D * sisl.unit_convert(\"eV\", \"meV\")\n", + " \n", + " print(pair[\"ai\"], pair[\"aj\"], pair[\"Ruc\"], \"distance\")\n", + " print(\"Isotropic: \", J_iso)\n", + " print(\"DMI: \", D)\n", + " print(\"Symmetric-anisotropy: \", J_S)\n", + " print(\"\")\n", + " \n", + " print(\"================================================================================\")\n", + " print(\"Runtime information: \")\n", + " print(\"Total runtime: \", end_time - start_time)\n", + " print(\"--------------------------------------------------------------------------------\")\n", + " print(\"Initial setup: \", setup_time - start_time)\n", + " print(f\"Hamiltonian conversion and XC field extraction: {H_and_XCF_time - setup_time:.3f} s\")\n", + " print(f\"Pair and site datastructure creatrions: {site_and_pair_dictionaries_time - H_and_XCF_time:.3f} s\")\n", + " print(f\"k set cration and distribution: {k_set_time - site_and_pair_dictionaries_time:.3f} s\")\n", + " print(f\"Rotating XC potential: {reference_rotations_time - k_set_time:.3f} s\")\n", + " print(f\"Greens function inversion: {green_function_inversion_time - reference_rotations_time:.3f} s\")\n", + " print(f\"Calculate energies and magnetic components: {end_time - green_function_inversion_time:.3f} s\")\n" ] }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "-60.72359053548613\n", - "-60.531975234450115\n", - "-60.524676226428824\n", - "-6.55042989834691\n", - "-6.047933492978864\n" + "(-0.06089330922308864, array([-9.32966923e-04, -8.92579299e-07, -2.04258660e-09]), array([-5.95741551e-03, 7.27737654e-03, 6.90431275e-07, -8.11057566e-07,\n", + " -5.49031203e-09]))\n", + "(-0.060556512255197724, array([ 3.78506176e-03, -6.13838308e-03, 3.59037036e-06]), array([-3.45657964e-04, 6.59197575e-04, 7.10694456e-05, -6.23149871e-03,\n", + " -4.24978017e-05]))\n" + ] + }, + { + "ename": "IndexError", + "evalue": "list index out of range", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[9], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(pairs[\u001b[38;5;241m0\u001b[39m])) \u001b[38;5;66;03m# isotropic should be -82 meV\u001b[39;00m\n\u001b[1;32m 2\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(pairs[\u001b[38;5;241m1\u001b[39m])) \u001b[38;5;66;03m# these should all be around -41.9 in the isotropic part\u001b[39;00m\n\u001b[0;32m----> 3\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(\u001b[43mpairs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[43m]\u001b[49m)) \u001b[38;5;66;03m# these should all be around -41.9 in the isotropic part\u001b[39;00m\n\u001b[1;32m 4\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(pairs[\u001b[38;5;241m3\u001b[39m])) \u001b[38;5;66;03m# these should all be around -41.9 in the isotropic part\u001b[39;00m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28mprint\u001b[39m(calculate_exchange_tensor(pairs[\u001b[38;5;241m4\u001b[39m])) \u001b[38;5;66;03m# these should all be around -41.9 in the isotropic part\u001b[39;00m\n", + "\u001b[0;31mIndexError\u001b[0m: list index out of range" ] } ], @@ -456,6 +661,8 @@ } ], "source": [ + "These are reasonably converged:\n", + "\n", "-61.33097171216109\n", "-60.52198328932686\n", "-60.51657719027764\n", diff --git a/test.py b/test.py index 29c1380..f020e45 100644 --- a/test.py +++ b/test.py @@ -55,7 +55,7 @@ pairs = [ kset = 20 kdirs = "xy" ebot = -30 -eset = 100 +eset = 50 esetp = 10000