{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "from sys import stdout\n", "from tqdm import tqdm\n", "from timeit import default_timer as timer\n", "\n", "os.environ[\"OMP_NUM_THREADS\"] = \"1\" # export OMP_NUM_THREADS=4\n", "os.environ[\"OPENBLAS_NUM_THREADS\"] = \"1\" # export OPENBLAS_NUM_THREADS=4\n", "os.environ[\"MKL_NUM_THREADS\"] = \"1\" # export MKL_NUM_THREADS=6\n", "os.environ[\"VECLIB_MAXIMUM_THREADS\"] = \"1\" # export VECLIB_MAXIMUM_THREADS=4\n", "os.environ[\"NUMEXPR_NUM_THREADS\"] = \"1\" # export NUMEXPR_NUM_THREADS=6\n", "\n", "import numpy as np\n", "import sisl\n", "from src.grogu_magn.useful import *\n", "from mpi4py import MPI\n", "from numpy.linalg import inv\n", "import warnings\n", "\n", "# runtime information\n", "times = dict()\n", "times[\"start_time\"] = timer()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this cell mimicks an input file\n", "fdf = sisl.get_sile(\"./lat3_791/Fe3GeTe2.fdf\") # ./Jij_for_Marci_6p45ang/CrBr.fdf\n", "# this information needs to be given at the input!!\n", "scf_xcf_orientation = np.array([0, 0, 1]) # z\n", "# list of reference directions for around which we calculate the derivatives\n", "# o is the quantization axis, v and w are two axes perpendicular to it\n", "# at this moment the user has to supply o,v,w on the input.\n", "# we can have some default for this\n", "ref_xcf_orientations = [\n", " dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),\n", " dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),\n", " dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),\n", "]\n", "\n", "\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", " # dict(atom=[3, 4],),\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", " # these should all be around -41.9 in the isotropic part\n", " # dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])),\n", " # dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])),\n", " # dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])),\n", " # dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])),\n", "]\n", "\n", "\"\"\"\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 = 20\n", "kdirs = \"xy\"\n", "ebot = -30\n", "eset = 50\n", "esetp = 10000\n", "\n", "\n", "# MPI parameters\n", "comm = MPI.COMM_WORLD\n", "size = comm.Get_size()\n", "rank = comm.Get_rank()\n", "root_node = 0\n", "if rank == root_node:\n", " print(\"Number of nodes in the parallel cluster: \", size)\n", "\n", "simulation_parameters = dict(\n", " path=\"Not yet specified.\",\n", " scf_xcf_orientation=scf_xcf_orientation,\n", " ref_xcf_orientations=ref_xcf_orientations,\n", " kset=kset,\n", " kdirs=kdirs,\n", " ebot=ebot,\n", " eset=eset,\n", " esetp=esetp,\n", " parallel_size=size,\n", ")\n", "\n", "# digestion of the input\n", "# read in hamiltonian\n", "dh = fdf.read_hamiltonian()\n", "try:\n", " simulation_parameters[\"geom\"] = fdf.read_geometry()\n", "except:\n", " print(\"Error reading geometry.\")\n", "\n", "# unit cell index\n", "uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])\n", "\n", "times[\"setup_time\"] = timer()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "NO = dh.no # shorthand for number of orbitals in the unit cell\n", "\n", "# preprocessing Hamiltonian and overlap matrix elements\n", "h11 = dh.tocsr(dh.M11r)\n", "h11 += dh.tocsr(dh.M11i) * 1.0j\n", "h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n", "\n", "h22 = dh.tocsr(dh.M22r)\n", "h22 += dh.tocsr(dh.M22i) * 1.0j\n", "h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n", "\n", "h12 = dh.tocsr(dh.M12r)\n", "h12 += dh.tocsr(dh.M12i) * 1.0j\n", "h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n", "\n", "h21 = dh.tocsr(dh.M21r)\n", "h21 += dh.tocsr(dh.M21i) * 1.0j\n", "h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n", "\n", "sov = (\n", " dh.tocsr(dh.S_idx)\n", " .toarray()\n", " .reshape(NO, dh.n_s, NO)\n", " .transpose(0, 2, 1)\n", " .astype(\"complex128\")\n", ")\n", "\n", "\n", "# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n", "U = np.vstack(\n", " [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])]\n", ")\n", "# This is the permutation that transforms ud1ud2 to u12d12\n", "# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n", "# the inverse transformation is U.T u12d12 to ud1ud2\n", "# That is FROM ORBITAL BOX to SPIN BOX => U.T\n", "\n", "# From now on everything is in SPIN BOX!!\n", "hh, ss = np.array(\n", " [\n", " U.T @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) @ U\n", " for i in range(dh.lattice.nsc.prod())\n", " ]\n", "), np.array(\n", " [\n", " U.T\n", " @ np.block([[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]])\n", " @ U\n", " for i in range(dh.lattice.nsc.prod())\n", " ]\n", ")\n", "\n", "\n", "# symmetrizing Hamiltonian and overlap matrix to make them hermitian\n", "for i in range(dh.lattice.sc_off.shape[0]):\n", " j = dh.lattice.sc_index(-dh.lattice.sc_off[i])\n", " h1, h1d = hh[i], hh[j]\n", " hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2\n", " s1, s1d = ss[i], ss[j]\n", " ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2\n", "\n", "# identifying TRS and TRB parts of the Hamiltonian\n", "TAUY = np.kron(np.eye(NO), tau_y)\n", "hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())])\n", "hTRS = (hh + hTR) / 2\n", "hTRB = (hh - hTR) / 2\n", "\n", "# extracting the exchange field\n", "traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77\n", "XCF = np.array(\n", " [\n", " np.array([f[\"x\"] for f in traced]),\n", " np.array([f[\"y\"] for f in traced]),\n", " np.array([f[\"z\"] for f in traced]),\n", " ]\n", ") # equation 77\n", "\n", "# Check if exchange field has scalar part\n", "max_xcfs = abs(np.array(np.array([f[\"c\"] for f in traced]))).max()\n", "if max_xcfs > 1e-12:\n", " warnings.warn(\n", " f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\"\n", " )\n", "\n", "times[\"H_and_XCF_time\"] = timer()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n", "for i, mag_ent in enumerate(magnetic_entities):\n", " parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes\n", " magnetic_entities[i][\"orbital_indeces\"] = parsed\n", " # calculate spin box indexes\n", " magnetic_entities[i][\"spin_box_indeces\"] = blow_up_orbindx(parsed)\n", " # calculate size for Greens function generation\n", " spin_box_shape = len(mag_ent[\"spin_box_indeces\"])\n", "\n", " mag_ent[\"energies\"] = [] # we will store the second order energy derivations here\n", "\n", " mag_ent[\"Gii\"] = [] # Greens function\n", " mag_ent[\"Gii_tmp\"] = [] # Greens function for parallelization\n", " # These will be the perturbed potentials from eq. 100\n", " mag_ent[\"Vu1\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n", " mag_ent[\"Vu2\"] = [list([]) for _ in range(len(ref_xcf_orientations))]\n", " for i in ref_xcf_orientations:\n", " # Greens functions for every quantization axis\n", " mag_ent[\"Gii\"].append(\n", " np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n", " )\n", " mag_ent[\"Gii_tmp\"].append(\n", " np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n", " )\n", "\n", "# for every site we have to store 2x3 Greens function (and the associated _tmp-s)\n", "# in the 3 reference directions, because G_ij and G_ji are both needed\n", "for pair in pairs:\n", " # calculate size for Greens function generation\n", " spin_box_shape_i = len(magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"])\n", " spin_box_shape_j = len(magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"])\n", "\n", " pair[\"energies\"] = [] # we will store the second order energy derivations here\n", "\n", " pair[\"Gij\"] = [] # Greens function\n", " pair[\"Gji\"] = []\n", " pair[\"Gij_tmp\"] = [] # Greens function for parallelization\n", " pair[\"Gji_tmp\"] = []\n", " for i in ref_xcf_orientations:\n", " # Greens functions for every quantization axis\n", " pair[\"Gij\"].append(\n", " np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n", " )\n", " pair[\"Gij_tmp\"].append(\n", " np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n", " )\n", " pair[\"Gji\"].append(\n", " np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n", " )\n", " pair[\"Gji_tmp\"].append(\n", " np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n", " )\n", "\n", "times[\"site_and_pair_dictionaries_time\"] = timer()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling\n", "wkset = np.ones(len(kset)) / len(kset) # generate weights for k points\n", "kpcs = np.array_split(kset, size) # split the k points based on MPI size\n", "kpcs[root_node] = tqdm(kpcs[root_node], desc=\"k loop\", file=stdout)\n", "\n", "times[\"k_set_time\"] = timer()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# this will contain the three hamiltonians in the reference directions needed to calculate the energy variations upon rotation\n", "hamiltonians = []\n", "\n", "# iterate over the reference directions (quantization axes)\n", "for i, orient in enumerate(ref_xcf_orientations):\n", " # obtain rotated exchange field\n", " R = RotMa2b(scf_xcf_orientation, orient[\"o\"])\n", " rot_XCF = np.einsum(\"ij,jklm->iklm\", R, XCF)\n", " rot_H_XCF = sum(\n", " [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])]\n", " )\n", " rot_H_XCF_uc = rot_H_XCF[uc_in_sc_idx]\n", "\n", " # obtain total Hamiltonian with the rotated exchange field\n", " rot_H = (\n", " hTRS + rot_H_XCF\n", " ) # equation 76 #######################################################################################\n", "\n", " hamiltonians.append(\n", " dict(orient=orient[\"o\"], H=rot_H)\n", " ) # store orientation and rotated Hamiltonian\n", "\n", " # these are the infinitezimal rotations (for now) perpendicular to the quantization axis\n", " for u in orient[\"vw\"]:\n", " Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H\n", "\n", " Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100\n", " Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100\n", "\n", " for mag_ent in magnetic_entities:\n", " # fill up the perturbed potentials (for now) based on the on-site projections\n", " mag_ent[\"Vu1\"][i].append(\n", " Vu1[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n", " )\n", " mag_ent[\"Vu2\"][i].append(\n", " Vu2[:, mag_ent[\"spin_box_indeces\"]][mag_ent[\"spin_box_indeces\"], :]\n", " )\n", "\n", "times[\"reference_rotations_time\"] = timer()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if rank == root_node:\n", " print(\"Number of magnetic entities being calculated: \", len(magnetic_entities))\n", " print(\n", " \"We have to calculate the Greens function for three reference direction and we are going to calculate 15 energy integrals per site.\"\n", " )\n", " print(f\"The shape of the Hamiltonian and the Greens function is {NO}x{NO}.\")\n", "comm.Barrier()\n", "# ----------------------------------------------------------------------\n", "\n", "# make energy contour\n", "# we are working in eV now !\n", "# and sisil shifts E_F to 0 !\n", "cont = make_contour(emin=ebot, enum=eset, p=esetp)\n", "eran = cont.ze\n", "\n", "# ----------------------------------------------------------------------\n", "# sampling the integrand on the contour and the BZ\n", "for k in kpcs[rank]:\n", " wk = wkset[rank] # weight of k point in BZ integral\n", " # iterate over reference directions\n", " for i, hamiltonian_orientation in enumerate(hamiltonians):\n", " # calculate Greens function\n", " H = hamiltonian_orientation[\"H\"]\n", " HK, SK = hsk(H, ss, dh.sc_off, k)\n", " Gk = inv(SK * eran.reshape(eset, 1, 1) - HK)\n", "\n", " # 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", "times[\"green_function_inversion_time\"] = timer()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if rank == root_node:\n", " # iterate over the magnetic entities\n", " for tracker, mag_ent in enumerate(magnetic_entities):\n", " # iterate over the quantization axes\n", " for i, Gii in enumerate(mag_ent[\"Gii\"]):\n", " storage = []\n", " # iterate over the first and second order local perturbations\n", " for Vu1, Vu2 in zip(mag_ent[\"Vu1\"][i], mag_ent[\"Vu2\"][i]):\n", " # The Szunyogh-Lichtenstein formula\n", " traced = np.trace((Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2)\n", " # evaluation of the contour integral\n", " storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n", "\n", " # fill up the magnetic entities dictionary with the energies\n", " magnetic_entities[tracker][\"energies\"].append(storage)\n", "\n", " # iterate over the pairs\n", " for tracker, pair in enumerate(pairs):\n", " # iterate over the quantization axes\n", " for i, (Gij, Gji) in enumerate(zip(pair[\"Gij\"], pair[\"Gji\"])):\n", " site_i = magnetic_entities[pair[\"ai\"]]\n", " site_j = magnetic_entities[pair[\"aj\"]]\n", "\n", " storage = []\n", " # iterate over the first order local perturbations in all possible orientations for the two sites\n", " for Vui in site_i[\"Vu1\"][i]:\n", " for Vuj in site_j[\"Vu1\"][i]:\n", " # The Szunyogh-Lichtenstein formula\n", " traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2)\n", " # evaluation of the contour integral\n", " storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n", " # fill up the pairs dictionary with the energies\n", " pairs[tracker][\"energies\"].append(storage)\n", "\n", " times[\"end_time\"] = timer()\n", " print_output(simulation_parameters, magnetic_entities, pairs, dh, times)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "##################################################################### GROGU OUTPUT #############################################################################\n", "================================================================================================================================================================\n", "Input file: \n", "Not yet specified.\n", "Number of nodes in the parallel cluster: 1\n", "================================================================================================================================================================\n", "Cell [Ang]: \n", "[[ 3.79100000e+00 0.00000000e+00 0.00000000e+00]\n", " [-1.89550000e+00 3.28310231e+00 0.00000000e+00]\n", " [ 1.25954923e-15 2.18160327e-15 2.05700000e+01]]\n", "================================================================================================================================================================\n", "DFT axis: \n", "[0 0 1]\n", "Quantization axis and perpendicular rotation directions:\n", "[1 0 0] --» [array([0, 1, 0]), array([0, 0, 1])]\n", "[0 1 0] --» [array([1, 0, 0]), array([0, 0, 1])]\n", "[0 0 1] --» [array([1, 0, 0]), array([0, 1, 0])]\n", "================================================================================================================================================================\n", "number of k points: 20\n", "k point directions: xy\n", "================================================================================================================================================================\n", "Parameters for the contour integral:\n", "Ebot: -30\n", "Eset: 50\n", "Esetp: 10000\n", "================================================================================================================================================================\n", "Atomic informations: \n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz\n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "[3]Fe(2) -7.339158738013707e-06 4.149278510690423e-06 11.657585837928032\n", "\n", "[4]Fe(2) -7.326987662162937e-06 4.158274523275774e-06 8.912422537596708\n", "\n", "[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n", "\n", "[3]Fe(all) -7.339158738013707e-06 4.149278510690423e-06 11.657585837928032\n", "[4]Fe(all) -7.326987662162937e-06 4.158274523275774e-06 8.912422537596708\n", "\n", "================================================================================================================================================================\n", "Exchange [meV]\n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n", "Isotropic: -63.511236909297565\n", "DMI: [-9.32966923e-01 -8.92579299e-04 -2.04258659e-06]\n", "Symmetric-anisotropy: [-3.33948782e+00 1.29796666e+00 6.90431275e-04 -8.11057566e-04\n", " -5.49031203e-06]\n", "Energies for debugging: \n", "array([[-6.22132703e-02, -9.32961433e-04, 9.32972414e-04,\n", " -6.14697157e-02],\n", " [-6.75757803e-02, 8.11057566e-07, -9.74101032e-07,\n", " -6.68507247e-02],\n", " [-5.36159327e-02, -6.92473862e-07, -6.88388688e-07,\n", " -5.36168088e-02]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.06685072, -0.05361593, -0.06221327])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.0668507247298092 -0.05361680878933305\n", "####################################################################################################################\n", "# This is from a different run\n", "[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] Not yet.\n", "Isotropic: -63.511236909297565\n", "DMI: [-9.32966923e-01 -8.92579299e-04 -2.04258659e-06]\n", "Symmetric-anisotropy: [-3.33948782e+00 1.29796666e+00 6.90431275e-04 -8.11057566e-04\n", " -5.49031203e-06]\n", "Energies for debugging: \n", "array([[-6.22132703e-02, -9.32961433e-04, 9.32972414e-04,\n", " -6.14697157e-02],\n", " [-6.75757803e-02, 8.11057566e-07, -9.74101032e-07,\n", " -6.68507247e-02],\n", " [-5.36159327e-02, -6.92473862e-07, -6.88388688e-07,\n", " -5.36168088e-02]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.06685072, -0.05361593, -0.06221327])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.0668507247298092 -0.05361680878933305\n", "####################################################################################################################\n", "[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n", "Isotropic: -60.982930301608\n", "DMI: [ 3.78506176e+00 -6.13838308e+00 3.59037036e-03]\n", "Symmetric-anisotropy: [ 0.08076008 0.11287844 0.07106945 -6.23149871 -0.0424978 ]\n", "Energies for debugging: \n", "array([[-6.08700519e-02, 3.82755956e-03, -3.74256396e-03,\n", " -6.11765688e-02],\n", " [-6.06451611e-02, 6.23149871e-03, -6.04526746e-03,\n", " -6.09021702e-02],\n", " [-5.98973147e-02, -6.74790753e-05, -7.46598160e-05,\n", " -5.98149915e-02]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.06090217, -0.05989731, -0.06087005])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.06090217021953539 -0.059814991464384966\n", "####################################################################################################################\n", "# This is from a different run\n", "[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n", "Isotropic: -60.982930301608\n", "DMI: [ 3.78506176e+00 -6.13838308e+00 3.59037036e-03]\n", "Symmetric-anisotropy: [ 0.08076008 0.11287844 0.07106945 -6.23149871 -0.0424978 ]\n", "Energies for debugging: \n", "array([[-6.08700519e-02, 3.82755956e-03, -3.74256396e-03,\n", " -6.11765688e-02],\n", " [-6.06451611e-02, 6.23149871e-03, -6.04526746e-03,\n", " -6.09021702e-02],\n", " [-5.98973147e-02, -6.74790753e-05, -7.46598160e-05,\n", " -5.98149915e-02]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.06090217, -0.05989731, -0.06087005])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.06090217021953539 -0.059814991464384966\n", "####################################################################################################################\n", "[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] Not yet.\n", "Isotropic: -60.97738429102214\n", "DMI: [-3.79945963e+00 6.15244494e+00 3.58990840e-03]\n", "Symmetric-anisotropy: [0.09125387 0.11162376 0.07106826 6.24185638 0.03636701]\n", "Energies for debugging: \n", "array([[-6.08657605e-02, -3.83582664e-03, 3.76309262e-03,\n", " -6.11802619e-02],\n", " [-6.06341379e-02, -6.24185638e-03, 6.06303351e-03,\n", " -6.08861304e-02],\n", " [-5.98973587e-02, -6.74783547e-05, -7.46581715e-05,\n", " -5.98150413e-02]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.06088613, -0.05989736, -0.06086576])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.06088613041916854 -0.05981504131171201\n", "\n", "[3]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] Not yet.\n", "Isotropic: -6.828322238590878\n", "DMI: [5.95251705 7.64859703 6.50501652]\n", "Symmetric-anisotropy: [-0.45523613 0.13725614 -0.031302 7.69961304 0.03239586]\n", "Energies for debugging: \n", "array([[-0.00669107, 0.00592012, -0.00598491, -0.00651034],\n", " [-0.00748301, -0.00769961, 0.00759758, -0.00728356],\n", " [-0.00590136, 0.00653632, -0.00647371, -0.00567313]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.00728356, -0.00590136, -0.00669107])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.007283558363709989 -0.005673134248977463\n", "\n", "[4]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] Not yet.\n", "Isotropic: -6.332622110203909\n", "DMI: [6.19414647 4.23019689 6.50504332]\n", "Symmetric-anisotropy: [ 0.53646663 -0.33927457 -0.03129943 4.24610256 -0.09833472]\n", "Energies for debugging: \n", "array([[-0.0066719 , 0.00629248, -0.00609581, -0.00652981],\n", " [-0.00599877, -0.0042461 , 0.00421429, -0.00579616],\n", " [-0.00590121, 0.00653634, -0.00647374, -0.00567299]])\n", "J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n", "array([-0.00579616, -0.00590121, -0.0066719 ])\n", "Test J_xx = E(y,z) = E(z,y)\n", "-0.005796155475608501 -0.00567298926736884\n", "\n", "================================================================================================================================================================\n", "Runtime information: \n", "Total runtime: 110.461674709 s\n", "----------------------------------------------------------------------------------------------------------------------------------------------------------------\n", "Initial setup: 0.15987975 s\n", "Hamiltonian conversion and XC field extraction: 0.601 s\n", "Pair and site datastructure creatrions: 0.010 s\n", "k set cration and distribution: 0.016 s\n", "Rotating XC potential: 0.289 s\n", "Greens function inversion: 109.307 s\n", "Calculate energies and magnetic components: 0.078 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" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sisl.viz\n", "\n", "dh.geometry.tile(2, 1).plot(axes=\"xy\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "\n", "coords = dh.xyz[-3:]\n", "\n", "\n", "plt.figure(figsize=(15, 5))\n", "plt.subplot(131)\n", "plt.scatter(coords[:, 0], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"z\")\n", "plt.subplot(132)\n", "plt.scatter(coords[:, 1], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n", "plt.xlabel(\"y\")\n", "plt.ylabel(\"z\")\n", "plt.subplot(133)\n", "plt.scatter(coords[:, 0], coords[:, 1], color=[\"r\", \"g\", \"b\"])\n", "plt.xlabel(\"x\")\n", "plt.ylabel(\"y\")\n", "print(\"xyz[-3:]: red, green, blue\")" ] } ], "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 }