You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1229 lines
58 KiB
1229 lines
58 KiB
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"[Mac:75276] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Mac.501/jf.0/1847132160/sm_segment.Mac.501.6e190000.0 could be created.\n"
|
|
]
|
|
},
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"'0.14.3'"
|
|
]
|
|
},
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"from tqdm import tqdm\n",
|
|
"from sys import getsizeof\n",
|
|
"from timeit import default_timer as timer\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"import sisl\n",
|
|
"import sisl.viz\n",
|
|
"from src.grogu_magn.utils import *\n",
|
|
"from mpi4py import MPI\n",
|
|
"import pickle\n",
|
|
"from numpy.linalg import inv\n",
|
|
"import warnings\n",
|
|
"\n",
|
|
"\"\"\" \n",
|
|
"# Some input parsing\n",
|
|
"parser = argparse.ArgumentParser()\n",
|
|
"parser.add_argument('--kset' , dest = 'kset' , default = 2 , type=int , help = 'k-space resolution of Jij calculation')\n",
|
|
"parser.add_argument('--kdirs' , dest = 'kdirs' , default = 'xyz' , help = 'Definition of k-space dimensionality')\n",
|
|
"parser.add_argument('--eset' , dest = 'eset' , default = 42 , type=int , help = 'Number of energy points on the contour')\n",
|
|
"parser.add_argument('--eset-p' , dest = 'esetp' , default = 10 , type=int , help = 'Parameter tuning the distribution on the contour')\n",
|
|
"parser.add_argument('--input' , dest = 'infile' , required = True , help = 'Input file name')\n",
|
|
"parser.add_argument('--output' , dest = 'outfile', required = True , help = 'Output file name')\n",
|
|
"parser.add_argument('--Ebot' , dest = 'Ebot' , default = -20.0 , type=float, help = 'Bottom energy of the contour')\n",
|
|
"parser.add_argument('--npairs' , dest = 'npairs' , default = 1 , type=int , help = 'Number of unitcell pairs in each direction for Jij calculation')\n",
|
|
"parser.add_argument('--adirs' , dest = 'adirs' , default = False , help = 'Definition of pair directions')\n",
|
|
"parser.add_argument('--use-tqdm', dest = 'usetqdm', default = 'not' , help = 'Use tqdm for progressbars or not')\n",
|
|
"parser.add_argument('--cutoff' , dest = 'cutoff' , default = 100.0 , type=float, help = 'Real space cutoff for pair generation in Angs')\n",
|
|
"parser.add_argument('--pairfile', dest = 'pairfile', default = False , help = 'File to read pair information')\n",
|
|
"args = parser.parse_args()\n",
|
|
"\"\"\"\n",
|
|
"# runtime information\n",
|
|
"times = dict()\n",
|
|
"times[\"start_time\"] = timer()\n",
|
|
"########################\n",
|
|
"# it works if data is in downloads folder\n",
|
|
"########################\n",
|
|
"sisl.__version__"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"file = open(\"Fe3GeTe2_benchmark_on_15k_300eset.pickle\", \"rb\")\n",
|
|
"object_file = pickle.load(file)\n",
|
|
"object_file"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"dat = sisl.io.siesta.eigSileSiesta(\n",
|
|
" \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.EIG\"\n",
|
|
")\n",
|
|
"siesta_eigs = dat.read_data()\n",
|
|
"siesta_eigs.min()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Ef = dat.read_fermi_level()\n",
|
|
"Ef"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"================================================================================================================================================================\n",
|
|
"Input file: \n",
|
|
"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\n",
|
|
"Output file: \n",
|
|
"./Fe3GeTe2_notebook.pickle\n",
|
|
"Number of nodes in the parallel cluster: 1\n",
|
|
"================================================================================================================================================================\n",
|
|
"Cell [Ang]: \n",
|
|
"[[ 3.79100000e+00 0.00000000e+00 0.00000000e+00]\n",
|
|
" [-1.89550000e+00 3.28310231e+00 0.00000000e+00]\n",
|
|
" [ 1.25954923e-15 2.18160327e-15 2.05700000e+01]]\n",
|
|
"================================================================================================================================================================\n",
|
|
"DFT axis: \n",
|
|
"[0 0 1]\n",
|
|
"Quantization axis and perpendicular rotation directions:\n",
|
|
"[1 0 0] --» [array([0, 1, 0]), array([0, 0, 1])]\n",
|
|
"[0 1 0] --» [array([1, 0, 0]), array([0, 0, 1])]\n",
|
|
"[0 0 1] --» [array([1, 0, 0]), array([0, 1, 0])]\n",
|
|
"================================================================================================================================================================\n",
|
|
"Parameters for the contour integral:\n",
|
|
"Number of k points: 10\n",
|
|
"k point directions: xy\n",
|
|
"Ebot: -13\n",
|
|
"Eset: 300\n",
|
|
"Esetp: 1000\n",
|
|
"================================================================================================================================================================\n",
|
|
"Setup done. Elapsed time: 4.909605666 s\n",
|
|
"================================================================================================================================================================\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"################################################################################\n",
|
|
"#################################### INPUT #####################################\n",
|
|
"################################################################################\n",
|
|
"path = (\n",
|
|
" \"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\"\n",
|
|
")\n",
|
|
"outfile = \"./Fe3GeTe2_notebook\"\n",
|
|
"\n",
|
|
"# this information needs to be given at the input!!\n",
|
|
"scf_xcf_orientation = np.array([0, 0, 1]) # z\n",
|
|
"# list of reference directions for around which we calculate the derivatives\n",
|
|
"# o is the quantization axis, v and w are two axes perpendicular to it\n",
|
|
"# at this moment the user has to supply o,v,w on the input.\n",
|
|
"# we can have some default for this\n",
|
|
"ref_xcf_orientations = [\n",
|
|
" dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),\n",
|
|
" dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),\n",
|
|
" dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),\n",
|
|
"]\n",
|
|
"magnetic_entities = [\n",
|
|
" dict(atom=3, l=2),\n",
|
|
" dict(atom=4, l=2),\n",
|
|
" dict(atom=5, l=2),\n",
|
|
"]\n",
|
|
"pairs = [\n",
|
|
" dict(ai=0, aj=1, Ruc=np.array([0, 0, 0])),\n",
|
|
" dict(ai=0, aj=2, Ruc=np.array([0, 0, 0])),\n",
|
|
" dict(ai=1, aj=2, Ruc=np.array([0, 0, 0])),\n",
|
|
" dict(ai=0, aj=2, Ruc=np.array([-1, -1, 0])),\n",
|
|
" dict(ai=1, aj=2, Ruc=np.array([-1, -1, 0])),\n",
|
|
" dict(ai=0, aj=2, Ruc=np.array([-1, 0, 0])),\n",
|
|
" dict(ai=1, aj=2, Ruc=np.array([-1, 0, 0])),\n",
|
|
" dict(ai=1, aj=2, Ruc=np.array([-2, 0, 0])),\n",
|
|
" dict(ai=1, aj=2, Ruc=np.array([-3, 0, 0])),\n",
|
|
"]\n",
|
|
"# Brilloun zone sampling and Green function contour integral\n",
|
|
"kset = 10\n",
|
|
"kdirs = \"xy\"\n",
|
|
"ebot = -13\n",
|
|
"eset = 300\n",
|
|
"esetp = 1000\n",
|
|
"################################################################################\n",
|
|
"#################################### INPUT #####################################\n",
|
|
"################################################################################\n",
|
|
"\n",
|
|
"# MPI parameters\n",
|
|
"comm = MPI.COMM_WORLD\n",
|
|
"size = comm.Get_size()\n",
|
|
"rank = comm.Get_rank()\n",
|
|
"root_node = 0\n",
|
|
"\n",
|
|
"# rename outfile\n",
|
|
"if not outfile.endswith(\".pickle\"):\n",
|
|
" outfile += \".pickle\"\n",
|
|
"\n",
|
|
"simulation_parameters = dict(\n",
|
|
" path=path,\n",
|
|
" outpath=outfile,\n",
|
|
" scf_xcf_orientation=scf_xcf_orientation,\n",
|
|
" ref_xcf_orientations=ref_xcf_orientations,\n",
|
|
" kset=kset,\n",
|
|
" kdirs=kdirs,\n",
|
|
" ebot=ebot,\n",
|
|
" eset=eset,\n",
|
|
" esetp=esetp,\n",
|
|
" parallel_size=size,\n",
|
|
")\n",
|
|
"\n",
|
|
"# digestion of the input\n",
|
|
"# read sile\n",
|
|
"fdf = sisl.get_sile(path)\n",
|
|
"# read in hamiltonian\n",
|
|
"dh = fdf.read_hamiltonian()\n",
|
|
"simulation_parameters[\"cell\"] = fdf.read_geometry().cell\n",
|
|
"\n",
|
|
"# unit cell index\n",
|
|
"uc_in_sc_idx = dh.lattice.sc_index([0, 0, 0])\n",
|
|
"\n",
|
|
"if rank == root_node:\n",
|
|
" print_parameters(simulation_parameters)\n",
|
|
" times[\"setup_time\"] = timer()\n",
|
|
" print(f\"Setup done. Elapsed time: {times['setup_time']} s\")\n",
|
|
" print(\n",
|
|
" \"================================================================================================================================================================\"\n",
|
|
" )"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n",
|
|
"plt.figure(figsize=(15, 5))\n",
|
|
"plt.subplot(121)\n",
|
|
"plt.plot(np.sort(dh.eigh()), marker=\"o\", linestyle=\" \", label=\"sisl\")\n",
|
|
"plt.plot(np.sort(siesta_eigs[0, 0]), marker=\"o\", linestyle=\" \", label=\"siesta\")\n",
|
|
"plt.ylim(None, 10)\n",
|
|
"plt.xlim(None, 75)\n",
|
|
"plt.legend()\n",
|
|
"plt.grid()\n",
|
|
"plt.subplot(122)\n",
|
|
"DOS = sisl.physics.electron.DOS(np.linspace(-15, 85, 1000), dh.eigh())\n",
|
|
"plt.plot(DOS, np.linspace(-15, 85, 1000))\n",
|
|
"DOS = sisl.physics.electron.DOS(np.linspace(-15, 85, 1000), siesta_eigs[0, 0])\n",
|
|
"plt.plot(DOS, np.linspace(-15, 85, 1000))\n",
|
|
"plt.ylim(None, 10)\n",
|
|
"\n",
|
|
"coords = dh.xyz[-3:]\n",
|
|
"\n",
|
|
"shift = np.array([-1, 0, 0]) @ simulation_parameters[\"cell\"]\n",
|
|
"\n",
|
|
"\n",
|
|
"plt.figure(figsize=(15, 5))\n",
|
|
"plt.subplot(131)\n",
|
|
"plt.scatter(coords[:, 0], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n",
|
|
"plt.scatter(\n",
|
|
" (coords + shift)[:, 0], (coords + shift)[:, 2], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n",
|
|
")\n",
|
|
"plt.xlabel(\"x\")\n",
|
|
"plt.ylabel(\"z\")\n",
|
|
"plt.subplot(132)\n",
|
|
"plt.scatter(coords[:, 1], coords[:, 2], color=[\"r\", \"g\", \"b\"])\n",
|
|
"plt.scatter(\n",
|
|
" (coords + shift)[:, 1], (coords + shift)[:, 2], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n",
|
|
")\n",
|
|
"plt.xlabel(\"y\")\n",
|
|
"plt.ylabel(\"z\")\n",
|
|
"plt.subplot(133)\n",
|
|
"plt.scatter(coords[:, 0], coords[:, 1], color=[\"r\", \"g\", \"b\"])\n",
|
|
"plt.scatter(\n",
|
|
" (coords + shift)[:, 0], (coords + shift)[:, 1], color=[\"r\", \"g\", \"b\"], marker=\"x\"\n",
|
|
")\n",
|
|
"plt.xlabel(\"x\")\n",
|
|
"plt.ylabel(\"y\")\n",
|
|
"print(\"xyz[-3:]: red, green, blue\")\n",
|
|
"\n",
|
|
"print(np.linalg.norm(coords[0] - coords[1]))\n",
|
|
"print(np.linalg.norm(coords[0] - coords[2]))\n",
|
|
"print(np.linalg.norm(coords[2] - coords[1]))\n",
|
|
"print(np.linalg.norm(coords[0] - (coords + shift)[2]))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Hamiltonian and exchange field rotated. Elapsed time: 5.473272791 s\n",
|
|
"================================================================================================================================================================\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"NO = dh.no # shorthand for number of orbitals in the unit cell\n",
|
|
"\n",
|
|
"# preprocessing Hamiltonian and overlap matrix elements\n",
|
|
"h11 = dh.tocsr(dh.M11r)\n",
|
|
"h11 += dh.tocsr(dh.M11i) * 1.0j\n",
|
|
"h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
|
|
"\n",
|
|
"h22 = dh.tocsr(dh.M22r)\n",
|
|
"h22 += dh.tocsr(dh.M22i) * 1.0j\n",
|
|
"h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
|
|
"\n",
|
|
"h12 = dh.tocsr(dh.M12r)\n",
|
|
"h12 += dh.tocsr(dh.M12i) * 1.0j\n",
|
|
"h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
|
|
"\n",
|
|
"h21 = dh.tocsr(dh.M21r)\n",
|
|
"h21 += dh.tocsr(dh.M21i) * 1.0j\n",
|
|
"h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
|
|
"\n",
|
|
"sov = (\n",
|
|
" dh.tocsr(dh.S_idx)\n",
|
|
" .toarray()\n",
|
|
" .reshape(NO, dh.n_s, NO)\n",
|
|
" .transpose(0, 2, 1)\n",
|
|
" .astype(\"complex128\")\n",
|
|
")\n",
|
|
"\n",
|
|
"\n",
|
|
"# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n",
|
|
"U = np.vstack(\n",
|
|
" [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])]\n",
|
|
")\n",
|
|
"# This is the permutation that transforms ud1ud2 to u12d12\n",
|
|
"# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n",
|
|
"# the inverse transformation is U.T u12d12 to ud1ud2\n",
|
|
"# That is FROM ORBITAL BOX to SPIN BOX => U.T\n",
|
|
"\n",
|
|
"# From now on everything is in SPIN BOX!!\n",
|
|
"hh, ss = np.array(\n",
|
|
" [\n",
|
|
" U.T @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) @ U\n",
|
|
" for i in range(dh.lattice.nsc.prod())\n",
|
|
" ]\n",
|
|
"), np.array(\n",
|
|
" [\n",
|
|
" U.T\n",
|
|
" @ np.block([[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]])\n",
|
|
" @ U\n",
|
|
" for i in range(dh.lattice.nsc.prod())\n",
|
|
" ]\n",
|
|
")\n",
|
|
"\n",
|
|
"\n",
|
|
"# symmetrizing Hamiltonian and overlap matrix to make them hermitian\n",
|
|
"for i in range(dh.lattice.sc_off.shape[0]):\n",
|
|
" j = dh.lattice.sc_index(-dh.lattice.sc_off[i])\n",
|
|
" h1, h1d = hh[i], hh[j]\n",
|
|
" hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2\n",
|
|
" s1, s1d = ss[i], ss[j]\n",
|
|
" ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2\n",
|
|
"\n",
|
|
"\n",
|
|
"###################################################################################\n",
|
|
"# either this is shit\n",
|
|
"\n",
|
|
"# identifying TRS and TRB parts of the Hamiltonian\n",
|
|
"TAUY = np.kron(np.eye(NO), tau_y)\n",
|
|
"hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())])\n",
|
|
"hTRS = (hh + hTR) / 2\n",
|
|
"hTRB = (hh - hTR) / 2\n",
|
|
"\n",
|
|
"# extracting the exchange field\n",
|
|
"traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77\n",
|
|
"XCF = np.array(\n",
|
|
" [\n",
|
|
" np.array([f[\"x\"] / 2 for f in traced]),\n",
|
|
" np.array([f[\"y\"] / 2 for f in traced]),\n",
|
|
" np.array([f[\"z\"] / 2 for f in traced]),\n",
|
|
" ]\n",
|
|
") # equation 77\n",
|
|
"\n",
|
|
"###################################################################################\n",
|
|
"\n",
|
|
"\n",
|
|
"# Check if exchange field has scalar part\n",
|
|
"max_xcfs = abs(np.array(np.array([f[\"c\"] / 2 for f in traced]))).max()\n",
|
|
"if max_xcfs > 1e-12:\n",
|
|
" warnings.warn(\n",
|
|
" f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\"\n",
|
|
" )\n",
|
|
"\n",
|
|
"if rank == root_node:\n",
|
|
" times[\"H_and_XCF_time\"] = timer()\n",
|
|
" print(\n",
|
|
" f\"Hamiltonian and exchange field rotated. Elapsed time: {times['H_and_XCF_time']} s\"\n",
|
|
" )\n",
|
|
" print(\n",
|
|
" \"================================================================================================================================================================\"\n",
|
|
" )"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"R [[1. 0. 0.]\n",
|
|
" [0. 1. 0.]\n",
|
|
" [0. 0. 1.]]\n",
|
|
"Z rot doesnt change anything in XCF: True\n",
|
|
"MAX(myhk - sislhk) 2.974862582050264e-06 this is zero as should be\n",
|
|
"MAX(myhk_rot - sislhk) 2.974862582050264e-06 rotation fucks it up\n",
|
|
"MAX(rot_H - sislhk) 8.881784197001252e-16 rotation fucks it up\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"myhk, sk = hsk(hh, ss, dh.sc_off, np.array([0.5, 0, 0]))\n",
|
|
"sislhk = dh.Hk(np.array([0.5, 0, 0])).toarray()\n",
|
|
"\n",
|
|
"R = RotMa2b(scf_xcf_orientation, ref_xcf_orientations[2][\"o\"])\n",
|
|
"print(\"R\", R)\n",
|
|
"rot_XCF = np.einsum(\"ij,jklm->iklm\", R, XCF)\n",
|
|
"print(\"Z rot doesnt change anything in XCF: \", np.allclose(XCF, rot_XCF))\n",
|
|
"\n",
|
|
"###################################################################################\n",
|
|
"# either this is shit\n",
|
|
"rot_H_XCF = sum(\n",
|
|
" [np.kron(rot_XCF[i], tau) for i, tau in enumerate([tau_x, tau_y, tau_z])]\n",
|
|
")\n",
|
|
"rot_H = hTRS + rot_H_XCF # equation 76\n",
|
|
"###################################################################################\n",
|
|
"\n",
|
|
"\n",
|
|
"myhk_rot, sk = hsk(rot_H, ss, dh.sc_off, np.array([0.5, 0, 0]))\n",
|
|
"\n",
|
|
"print(\"MAX(myhk - sislhk)\", np.max(abs(myhk - sislhk)), \"this is zero as should be\")\n",
|
|
"print(\"MAX(myhk_rot - sislhk)\", np.max(abs(myhk_rot - sislhk)), \"rotation fucks it up\")\n",
|
|
"print(\"MAX(rot_H - sislhk)\", np.max(abs(rot_H - hh)), \"rotation fucks it up\")"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Site and pair dictionaries created. Elapsed time: 5.680114458 s\n",
|
|
"================================================================================================================================================================\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n",
|
|
"for mag_ent in magnetic_entities:\n",
|
|
" parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes\n",
|
|
" mag_ent[\"orbital_indeces\"] = parsed\n",
|
|
" mag_ent[\"spin_box_indeces\"] = blow_up_orbindx(parsed) # calculate spin box indexes\n",
|
|
" # if orbital is not set use all\n",
|
|
" if \"l\" not in mag_ent.keys():\n",
|
|
" mag_ent[\"l\"] = \"all\"\n",
|
|
" if isinstance(mag_ent[\"atom\"], int):\n",
|
|
" mag_ent[\"tags\"] = [\n",
|
|
" f\"[{mag_ent['atom']}]{dh.atoms[mag_ent['atom']].tag}({mag_ent['l']})\"\n",
|
|
" ]\n",
|
|
" mag_ent[\"xyz\"] = [dh.xyz[mag_ent[\"atom\"]]]\n",
|
|
" if isinstance(mag_ent[\"atom\"], list):\n",
|
|
" mag_ent[\"tags\"] = []\n",
|
|
" mag_ent[\"xyz\"] = []\n",
|
|
" # iterate over atoms\n",
|
|
" for atom_idx in mag_ent[\"atom\"]:\n",
|
|
" mag_ent[\"tags\"].append(\n",
|
|
" f\"[{atom_idx}]{dh.atoms[atom_idx].tag}({mag_ent['l']})\"\n",
|
|
" )\n",
|
|
" mag_ent[\"xyz\"].append(dh.xyz[atom_idx])\n",
|
|
"\n",
|
|
" # calculate size for Greens function generation\n",
|
|
" spin_box_shape = len(mag_ent[\"spin_box_indeces\"])\n",
|
|
"\n",
|
|
" mag_ent[\"energies\"] = [] # we will store the second order energy derivations here\n",
|
|
"\n",
|
|
" # These will be the perturbed potentials from eq. 100\n",
|
|
" mag_ent[\"Vu1\"] = [] # so they are independent in memory\n",
|
|
" mag_ent[\"Vu2\"] = []\n",
|
|
"\n",
|
|
" mag_ent[\"Gii\"] = [] # Greens function\n",
|
|
" mag_ent[\"Gii_tmp\"] = [] # Greens function for parallelization\n",
|
|
" for i in ref_xcf_orientations:\n",
|
|
" # Rotations for every quantization axis\n",
|
|
" mag_ent[\"Vu1\"].append([])\n",
|
|
" mag_ent[\"Vu2\"].append([])\n",
|
|
" # Greens functions for every quantization axis\n",
|
|
" mag_ent[\"Gii\"].append(\n",
|
|
" np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n",
|
|
" )\n",
|
|
" mag_ent[\"Gii_tmp\"].append(\n",
|
|
" np.zeros((eset, spin_box_shape, spin_box_shape), dtype=\"complex128\")\n",
|
|
" )\n",
|
|
"\n",
|
|
"# for every site we have to store 2x3 Greens function (and the associated _tmp-s)\n",
|
|
"# in the 3 reference directions, because G_ij and G_ji are both needed\n",
|
|
"for pair in pairs:\n",
|
|
" # calculate distance\n",
|
|
" xyz_ai = magnetic_entities[pair[\"ai\"]][\"xyz\"]\n",
|
|
" xyz_aj = magnetic_entities[pair[\"aj\"]][\"xyz\"]\n",
|
|
" xyz_aj = xyz_aj + pair[\"Ruc\"] @ simulation_parameters[\"cell\"]\n",
|
|
" pair[\"dist\"] = np.linalg.norm(xyz_ai - xyz_aj)\n",
|
|
"\n",
|
|
" # calculate size for Greens function generation\n",
|
|
" spin_box_shape_i = len(magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"])\n",
|
|
" spin_box_shape_j = len(magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"])\n",
|
|
" pair[\"tags\"] = []\n",
|
|
" for mag_ent in [magnetic_entities[pair[\"ai\"]], magnetic_entities[pair[\"aj\"]]]:\n",
|
|
" tag = \"\"\n",
|
|
" # get atoms of magnetic entity\n",
|
|
" atoms_idx = mag_ent[\"atom\"]\n",
|
|
" orbitals = mag_ent[\"l\"]\n",
|
|
"\n",
|
|
" # if magnetic entity contains one atoms\n",
|
|
" if isinstance(atoms_idx, int):\n",
|
|
" tag += f\"[{atoms_idx}]{dh.atoms[atoms_idx].tag}({orbitals})\"\n",
|
|
"\n",
|
|
" # if magnetic entity contains more than one atoms\n",
|
|
" if isinstance(atoms_idx, list):\n",
|
|
" # iterate over atoms\n",
|
|
" atom_group = \"{\"\n",
|
|
" for atom_idx in atoms_idx:\n",
|
|
" atom_group += f\"[{atom_idx}]{dh.atoms[atom_idx].tag}({orbitals})--\"\n",
|
|
" # end {} of the atoms in the magnetic entity\n",
|
|
" tag += atom_group[:-2] + \"}\"\n",
|
|
" pair[\"tags\"].append(tag)\n",
|
|
" pair[\"energies\"] = [] # we will store the second order energy derivations here\n",
|
|
"\n",
|
|
" pair[\"Gij\"] = [] # Greens function\n",
|
|
" pair[\"Gji\"] = []\n",
|
|
" pair[\"Gij_tmp\"] = [] # Greens function for parallelization\n",
|
|
" pair[\"Gji_tmp\"] = []\n",
|
|
" for i in ref_xcf_orientations:\n",
|
|
" # Greens functions for every quantization axis\n",
|
|
" pair[\"Gij\"].append(\n",
|
|
" np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n",
|
|
" )\n",
|
|
" pair[\"Gij_tmp\"].append(\n",
|
|
" np.zeros((eset, spin_box_shape_i, spin_box_shape_j), dtype=\"complex128\")\n",
|
|
" )\n",
|
|
" pair[\"Gji\"].append(\n",
|
|
" np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n",
|
|
" )\n",
|
|
" pair[\"Gji_tmp\"].append(\n",
|
|
" np.zeros((eset, spin_box_shape_j, spin_box_shape_i), dtype=\"complex128\")\n",
|
|
" )\n",
|
|
"\n",
|
|
"if rank == root_node:\n",
|
|
" times[\"site_and_pair_dictionaries_time\"] = timer()\n",
|
|
" print(\n",
|
|
" f\"Site and pair dictionaries created. Elapsed time: {times['site_and_pair_dictionaries_time']} s\"\n",
|
|
" )\n",
|
|
" print(\n",
|
|
" \"================================================================================================================================================================\"\n",
|
|
" )"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"k loop: 0%| | 0/100 [00:00<?, ?it/s]"
|
|
]
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"k set created. Elapsed time: 5.716477833 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": 7,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Rotations done perpendicular to quantization axis. Elapsed time: 6.074103166 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 = hTRS + rot_H_XCF # equation 76\n",
|
|
"\n",
|
|
" hamiltonians.append(\n",
|
|
" dict(\n",
|
|
" orient=orient[\"o\"],\n",
|
|
" H=rot_H,\n",
|
|
" GS=np.zeros((eset, rot_H.shape[1], rot_H.shape[2]), dtype=\"complex128\"),\n",
|
|
" GS_tmp=np.zeros((eset, rot_H.shape[1], rot_H.shape[2]), dtype=\"complex128\"),\n",
|
|
" )\n",
|
|
" ) # store orientation and rotated Hamiltonian\n",
|
|
"\n",
|
|
" # these are the rotations (for now) perpendicular to the quantization axis\n",
|
|
" for u in orient[\"vw\"]:\n",
|
|
" Tu = np.kron(np.eye(NO, dtype=int), tau_u(u)) # section 2.H\n",
|
|
"\n",
|
|
" Vu1 = 1j / 2 * commutator(rot_H_XCF_uc, Tu) # equation 100\n",
|
|
" Vu2 = 1 / 8 * commutator(commutator(Tu, rot_H_XCF_uc), Tu) # equation 100\n",
|
|
"\n",
|
|
" for mag_ent in magnetic_entities:\n",
|
|
" idx = mag_ent[\"spin_box_indeces\"]\n",
|
|
" # fill up the perturbed potentials (for now) based on the on-site projections\n",
|
|
" mag_ent[\"Vu1\"][i].append(Vu1[:, idx][idx, :])\n",
|
|
" mag_ent[\"Vu2\"][i].append(Vu2[:, idx][idx, :])\n",
|
|
"\n",
|
|
"if rank == root_node:\n",
|
|
" times[\"reference_rotations_time\"] = timer()\n",
|
|
" print(\n",
|
|
" f\"Rotations done perpendicular to quantization axis. Elapsed time: {times['reference_rotations_time']} s\"\n",
|
|
" )\n",
|
|
" print(\n",
|
|
" \"================================================================================================================================================================\"\n",
|
|
" )"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Starting matrix inversions\n",
|
|
"Total number of k points: 100\n",
|
|
"Number of energy samples per k point: 300\n",
|
|
"Total number of directions: 3\n",
|
|
"Total number of matrix inversions: 90000\n",
|
|
"The shape of the Hamiltonian and the Greens function is 84x84=7056\n",
|
|
"Memory taken by a single Hamiltonian is: 0.015625 KB\n",
|
|
"Expected memory usage per matrix inversion: 0.5 KB\n",
|
|
"Expected memory usage per k point for parallel inversion: 450.0 KB\n",
|
|
"Expected memory usage on root node: 43.9453125 MB\n",
|
|
"================================================================================================================================================================\n"
|
|
]
|
|
},
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"k loop: 100%|██████████| 100/100 [03:17<00:00, 1.97s/it]"
|
|
]
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Calculated Greens functions. Elapsed time: 203.218047166 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 sisl shifts E_F to 0 !\n",
|
|
"cont = make_contour(emin=ebot, enum=eset, p=esetp)\n",
|
|
"eran = cont.ze\n",
|
|
"\n",
|
|
"# ----------------------------------------------------------------------\n",
|
|
"# sampling the integrand on the contour and the BZ\n",
|
|
"for k in kpcs[rank]:\n",
|
|
" wk = wkset[rank] # weight of k point in BZ integral\n",
|
|
" # iterate over reference directions\n",
|
|
" for i, hamiltonian_orientation in enumerate(hamiltonians):\n",
|
|
" # calculate Greens function\n",
|
|
" H = hamiltonian_orientation[\"H\"]\n",
|
|
" HK, SK = hsk(H, ss, dh.sc_off, k)\n",
|
|
"\n",
|
|
" # solve Greens function sequentially for the energies, because of memory bound\n",
|
|
" Gk = np.zeros(shape=(eset, HK.shape[0], HK.shape[1]), dtype=\"complex128\")\n",
|
|
" for j in range(eset):\n",
|
|
" Gk[j] = inv(SK * eran[j] - HK)\n",
|
|
"\n",
|
|
" # saving this for total charge\n",
|
|
" hamiltonian_orientation[\"GS_tmp\"] += Gk @ SK * wk\n",
|
|
"\n",
|
|
" # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
|
|
" for mag_ent in magnetic_entities:\n",
|
|
" mag_ent[\"Gii_tmp\"][i] += (\n",
|
|
" Gk[:, mag_ent[\"spin_box_indeces\"], :][:, :, mag_ent[\"spin_box_indeces\"]]\n",
|
|
" * wk\n",
|
|
" )\n",
|
|
"\n",
|
|
" for pair in pairs:\n",
|
|
" # add phase shift based on the cell difference\n",
|
|
" phase = np.exp(1j * 2 * np.pi * k @ pair[\"Ruc\"].T)\n",
|
|
"\n",
|
|
" # get the pair orbital sizes from the magnetic entities\n",
|
|
" ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]\n",
|
|
" aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n",
|
|
"\n",
|
|
" # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
|
|
" pair[\"Gij_tmp\"][i] += Gk[:, ai][..., aj] * phase * wk\n",
|
|
" pair[\"Gji_tmp\"][i] += Gk[:, aj][..., ai] / phase * wk\n",
|
|
"\n",
|
|
"# summ reduce partial results of mpi nodes\n",
|
|
"for i in range(len(hamiltonians)):\n",
|
|
" # for total charge\n",
|
|
" comm.Reduce(hamiltonians[i][\"GS_tmp\"], hamiltonians[i][\"GS\"], root=root_node)\n",
|
|
"\n",
|
|
" for mag_ent in magnetic_entities:\n",
|
|
" comm.Reduce(mag_ent[\"Gii_tmp\"][i], mag_ent[\"Gii\"][i], root=root_node)\n",
|
|
"\n",
|
|
" for pair in pairs:\n",
|
|
" comm.Reduce(pair[\"Gij_tmp\"][i], pair[\"Gij\"][i], root=root_node)\n",
|
|
" comm.Reduce(pair[\"Gji_tmp\"][i], pair[\"Gji\"][i], root=root_node)\n",
|
|
"\n",
|
|
"if rank == root_node:\n",
|
|
" times[\"green_function_inversion_time\"] = timer()\n",
|
|
" print(\n",
|
|
" f\"Calculated Greens functions. Elapsed time: {times['green_function_inversion_time']} s\"\n",
|
|
" )\n",
|
|
" print(\n",
|
|
" \"================================================================================================================================================================\"\n",
|
|
" )"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Total charge: 39.857798089095745\n",
|
|
"Total charge: 39.85779799083869\n",
|
|
"Total charge: 39.95778482001014\n",
|
|
"Magnetic entities integrated.\n",
|
|
"Pairs integrated.\n",
|
|
"Magnetic parameters calculated.\n",
|
|
"##################################################################### GROGU OUTPUT #############################################################################\n",
|
|
"================================================================================================================================================================\n",
|
|
"Input file: \n",
|
|
"/Users/danielpozsar/Downloads/nojij/Fe3GeTe2/monolayer/soc/lat3_791/Fe3GeTe2.fdf\n",
|
|
"Output file: \n",
|
|
"./Fe3GeTe2_notebook.pickle\n",
|
|
"Number of nodes in the parallel cluster: 1\n",
|
|
"================================================================================================================================================================\n",
|
|
"Cell [Ang]: \n",
|
|
"[[ 3.79100000e+00 0.00000000e+00 0.00000000e+00]\n",
|
|
" [-1.89550000e+00 3.28310231e+00 0.00000000e+00]\n",
|
|
" [ 1.25954923e-15 2.18160327e-15 2.05700000e+01]]\n",
|
|
"================================================================================================================================================================\n",
|
|
"DFT axis: \n",
|
|
"[0 0 1]\n",
|
|
"Quantization axis and perpendicular rotation directions:\n",
|
|
"[1 0 0] --» [array([0, 1, 0]), array([0, 0, 1])]\n",
|
|
"[0 1 0] --» [array([1, 0, 0]), array([0, 0, 1])]\n",
|
|
"[0 0 1] --» [array([1, 0, 0]), array([0, 1, 0])]\n",
|
|
"================================================================================================================================================================\n",
|
|
"Parameters for the contour integral:\n",
|
|
"Number of k points: 10\n",
|
|
"k point directions: xy\n",
|
|
"Ebot: -13\n",
|
|
"Eset: 300\n",
|
|
"Esetp: 1000\n",
|
|
"================================================================================================================================================================\n",
|
|
"Atomic information: \n",
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
|
|
"[atom index]Element(orbitals) x [Ang] y [Ang] z [Ang] Sx Sy Sz Q Lx Ly Lz Jx Jy Jz\n",
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
|
|
"[3]Fe(2) -7.339158738013707e-06 4.149278510690423e-06 11.657585837928032\n",
|
|
"\n",
|
|
"[4]Fe(2) -7.326987662162937e-06 4.158274523275774e-06 8.912422537596708\n",
|
|
"\n",
|
|
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
|
|
"\n",
|
|
"================================================================================================================================================================\n",
|
|
"Anisotropy [meV]\n",
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
|
|
"Magnetic entity x [Ang] y [Ang] z [Ang]\n",
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
|
|
"[3]Fe(2) -7.339158738013707e-06 4.149278510690423e-06 11.657585837928032\n",
|
|
"Consistency check: 6.818420305920725e-05\n",
|
|
"Anisotropy diag: [-0.07482326 -0.00552873 0. ]\n",
|
|
"\n",
|
|
"[4]Fe(2) -7.326987662162937e-06 4.158274523275774e-06 8.912422537596708\n",
|
|
"Consistency check: 7.686025895869975e-05\n",
|
|
"Anisotropy diag: [ 0.0004411 -0.07750451 0. ]\n",
|
|
"\n",
|
|
"[5]Fe(2) 1.8954667088117545 1.0943913231921656 10.285002698393109\n",
|
|
"Consistency check: 5.268649241663148e-05\n",
|
|
"Anisotropy diag: [ 0.01120506 -0.0414785 0. ]\n",
|
|
"\n",
|
|
"================================================================================================================================================================\n",
|
|
"Exchange [meV]\n",
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
|
|
"Magnetic entity1 Magnetic entity2 [i j k] d [Ang]\n",
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
|
|
"[3]Fe(2) [4]Fe(2) [0 0 0] d [Ang] 2.745163300331324\n",
|
|
"Isotropic: -84.46399731420344\n",
|
|
"DMI: [ 1.88620081e-01 -7.64018484e-04 -4.37922496e-07]\n",
|
|
"Symmetric-anisotropy: [ 2.38296852e+00 1.42209663e-04 -6.68866203e-05 1.42209663e-04\n",
|
|
" 2.40029619e+00 -8.14152579e-06 -6.68866203e-05 -8.14152579e-06\n",
|
|
" -4.78326471e+00]\n",
|
|
"J: [-8.20810288e+01 1.42209663e-04 -6.68866203e-05 1.42209663e-04\n",
|
|
" -8.20637011e+01 -8.14152579e-06 -6.68866203e-05 -8.14152579e-06\n",
|
|
" -8.92472620e+01]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-8.92399658e-02, 1.88628223e-04, -1.88611940e-04,\n",
|
|
" -8.91429459e-02],\n",
|
|
" [-8.92545582e-02, 8.30905104e-07, -6.97131863e-07,\n",
|
|
" -8.91774451e-02],\n",
|
|
" [-7.49844563e-02, -1.42647586e-07, -1.41771741e-07,\n",
|
|
" -7.49846125e-02]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.08917745, -0.07498446, -0.08923997])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.08917744514035186 -0.07498461245605947\n",
|
|
"\n",
|
|
"[3]Fe(2) [5]Fe(2) [0 0 0] d [Ang] 2.5835033632437767\n",
|
|
"Isotropic: -41.629965239513076\n",
|
|
"DMI: [ 1.16104440e+00 -2.02989607e+00 1.42007213e-04]\n",
|
|
"Symmetric-anisotropy: [-0.3583508 0.2300684 -0.0464857 0.2300684 -0.60812243 -0.03183621\n",
|
|
" -0.0464857 -0.03183621 0.96647323]\n",
|
|
"J: [-4.19883160e+01 2.30068404e-01 -4.64857039e-02 2.30068404e-01\n",
|
|
" -4.22380877e+01 -3.18362116e-02 -4.64857039e-02 -3.18362116e-02\n",
|
|
" -4.06634920e+01]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-0.04075812, 0.00119288, -0.00112921, -0.04080062],\n",
|
|
" [-0.04056886, 0.00207638, -0.00198341, -0.04056711],\n",
|
|
" [-0.04367555, -0.00022993, -0.00023021, -0.04340952]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.04056711, -0.04367555, -0.04075812])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.04056710918296237 -0.04340952290351278\n",
|
|
"\n",
|
|
"[4]Fe(2) [5]Fe(2) [0 0 0] d [Ang] 2.583501767937866\n",
|
|
"Isotropic: -41.63818539049319\n",
|
|
"DMI: [-1.15569804e+00 2.01131866e+00 1.42813172e-04]\n",
|
|
"Symmetric-anisotropy: [-0.35474671 0.23007014 0.03254717 0.23007014 -0.61181222 0.01758968\n",
|
|
" 0.03254717 0.01758968 0.96655893]\n",
|
|
"J: [-4.19929321e+01 2.30070141e-01 3.25471741e-02 2.30070141e-01\n",
|
|
" -4.22499976e+01 1.75896770e-02 3.25471741e-02 1.75896770e-02\n",
|
|
" -4.06716265e+01]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-0.04076979, -0.00117329, 0.00113811, -0.04082375],\n",
|
|
" [-0.04057347, -0.00204387, 0.00197877, -0.04057565],\n",
|
|
" [-0.04367625, -0.00022993, -0.00023021, -0.04341022]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.04057565, -0.04367625, -0.04076979])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.04057564848089523 -0.043410215712362775\n",
|
|
"\n",
|
|
"[3]Fe(2) [5]Fe(2) [-1 -1 0] d [Ang] 2.5834973202859075\n",
|
|
"Isotropic: -41.66159765171782\n",
|
|
"DMI: [-2.33447929e+00 5.88920633e-05 -4.25773045e-05]\n",
|
|
"Symmetric-anisotropy: [-7.33142752e-01 -1.52756996e-04 1.46307682e-04 -1.52756996e-04\n",
|
|
" -2.16761810e-01 4.07813099e-02 1.46307682e-04 4.07813099e-02\n",
|
|
" 9.49904562e-01]\n",
|
|
"J: [-4.23947404e+01 -1.52756996e-04 1.46307682e-04 -1.52756996e-04\n",
|
|
" -4.18783595e+01 4.07813099e-02 1.46307682e-04 4.07813099e-02\n",
|
|
" -4.07116931e+01]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-4.05161871e-02, -2.37526060e-03, 2.29369798e-03,\n",
|
|
" -4.04786366e-02],\n",
|
|
" [-4.09071991e-02, -2.05199745e-07, -8.74156186e-08,\n",
|
|
" -4.09798629e-02],\n",
|
|
" [-4.32780824e-02, 1.10179692e-07, 1.95334301e-07,\n",
|
|
" -4.38096179e-02]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.04097986, -0.04327808, -0.04051619])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.040979862867913205 -0.043809617939982115\n",
|
|
"\n",
|
|
"[4]Fe(2) [5]Fe(2) [-1 -1 0] d [Ang] 2.583495745338251\n",
|
|
"Isotropic: -41.651354364546386\n",
|
|
"DMI: [ 2.33450755e+00 -3.01207900e-04 -4.35448493e-05]\n",
|
|
"Symmetric-anisotropy: [-7.24739149e-01 -1.52908277e-04 1.43933638e-04 -1.52908277e-04\n",
|
|
" -2.27695138e-01 -4.07805826e-02 1.43933638e-04 -4.07805826e-02\n",
|
|
" 9.52434287e-01]\n",
|
|
"J: [-4.23760935e+01 -1.52908277e-04 1.43933638e-04 -1.52908277e-04\n",
|
|
" -4.18790495e+01 -4.07805826e-02 1.43933638e-04 -4.07805826e-02\n",
|
|
" -4.06989201e+01]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-4.05168794e-02, 2.37528813e-03, -2.29372696e-03,\n",
|
|
" -4.04793245e-02],\n",
|
|
" [-4.08809608e-02, 1.57274262e-07, -4.45141537e-07,\n",
|
|
" -4.09418727e-02],\n",
|
|
" [-4.32787745e-02, 1.09363428e-07, 1.96453126e-07,\n",
|
|
" -4.38103143e-02]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.04094187, -0.04327877, -0.04051688])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.04094187270410039 -0.04381031432397496\n",
|
|
"\n",
|
|
"[3]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] 2.583541444641373\n",
|
|
"Isotropic: -41.622961610451405\n",
|
|
"DMI: [ 1.15553716e+00 2.03004764e+00 -1.00662008e-04]\n",
|
|
"Symmetric-anisotropy: [-0.35259845 -0.22988479 0.04629996 -0.22988479 -0.61302608 -0.01757494\n",
|
|
" 0.04629996 -0.01757494 0.96562453]\n",
|
|
"J: [-4.19755601e+01 -2.29884795e-01 4.62999578e-02 -2.29884795e-01\n",
|
|
" -4.22359877e+01 -1.75749425e-02 4.62999578e-02 -1.75749425e-02\n",
|
|
" -4.06573371e+01]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-0.04075738, 0.00117311, -0.00113796, -0.04081079],\n",
|
|
" [-0.04055729, -0.00207635, 0.00198375, -0.0405554 ],\n",
|
|
" [-0.04366118, 0.00022978, 0.00022999, -0.04339572]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.0405554 , -0.04366118, -0.04075738])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.040555403759627504 -0.043395716365454126\n",
|
|
"\n",
|
|
"[4]Fe(2) [5]Fe(2) [-1 0 0] d [Ang] 2.5835398672184064\n",
|
|
"Isotropic: -41.62008807379083\n",
|
|
"DMI: [-1.16093102e+00 -2.01118580e+00 -9.97040031e-05]\n",
|
|
"Symmetric-anisotropy: [-0.36016123 -0.22988672 -0.03265224 -0.22988672 -0.60539869 0.03182552\n",
|
|
" -0.03265224 0.03182552 0.96555993]\n",
|
|
"J: [-4.19802493e+01 -2.29886721e-01 -3.26522393e-02 -2.29886721e-01\n",
|
|
" -4.22254868e+01 3.18255164e-02 -3.26522393e-02 3.18255164e-02\n",
|
|
" -4.06545281e+01]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-0.04074714, -0.00119276, 0.00112911, -0.04078908],\n",
|
|
" [-0.04056191, 0.00204384, -0.00197853, -0.04056407],\n",
|
|
" [-0.04366189, 0.00022979, 0.00022999, -0.04339642]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.04056407, -0.04366189, -0.04074714])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.04056407415889504 -0.043396424451279274\n",
|
|
"\n",
|
|
"[4]Fe(2) [5]Fe(2) [-2 0 0] d [Ang] 5.951322298958084\n",
|
|
"Isotropic: -1.7294071924601657\n",
|
|
"DMI: [-0.08510436 0.24159561 1.44684598]\n",
|
|
"Symmetric-anisotropy: [ 0.09622056 0.08175665 0.02144594 0.08175665 -0.06049024 -0.0414623\n",
|
|
" 0.02144594 -0.0414623 -0.03573033]\n",
|
|
"J: [-1.63318663 0.08175665 0.02144594 0.08175665 -1.78989743 -0.0414623\n",
|
|
" 0.02144594 -0.0414623 -1.76513752]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-1.72606525e-03, -4.36420623e-05, 1.26566667e-04,\n",
|
|
" -1.82716801e-03],\n",
|
|
" [-1.80420979e-03, -2.63041556e-04, 2.20149672e-04,\n",
|
|
" -1.73934458e-03],\n",
|
|
" [-1.75262685e-03, 1.36508933e-03, -1.52860262e-03,\n",
|
|
" -1.52702868e-03]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.00173934, -0.00175263, -0.00172607])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.0017393445796559598 -0.0015270286786416\n",
|
|
"\n",
|
|
"[4]Fe(2) [5]Fe(2) [-3 0 0] d [Ang] 9.638732176310562\n",
|
|
"Isotropic: -0.4535195863892767\n",
|
|
"DMI: [ 0.04836847 0.05429183 -0.15906095]\n",
|
|
"Symmetric-anisotropy: [ 0.14297234 -0.0154665 0.00345288 -0.0154665 0.19102043 0.01270769\n",
|
|
" 0.00345288 0.01270769 -0.33399277]\n",
|
|
"J: [-0.31054725 -0.0154665 0.00345288 -0.0154665 -0.26249915 0.01270769\n",
|
|
" 0.00345288 0.01270769 -0.78751236]\n",
|
|
"Energies for debugging: \n",
|
|
"array([[-7.80081681e-04, 3.56607876e-05, -6.10761579e-05,\n",
|
|
" -7.42232776e-04],\n",
|
|
" [-7.94943031e-04, -5.77447108e-05, 5.08389441e-05,\n",
|
|
" -8.34582436e-04],\n",
|
|
" [ 2.17234470e-04, -1.43594444e-04, 1.74527450e-04,\n",
|
|
" 2.13487935e-04]])\n",
|
|
"J_ii for debugging: (check if this is the same as in calculate_exchange_tensor)\n",
|
|
"array([-0.00083458, 0.00021723, -0.00078008])\n",
|
|
"Test J_xx = E(y,z) = E(z,y)\n",
|
|
"-0.0008345824358667467 0.0002134879350515751\n",
|
|
"\n",
|
|
"================================================================================================================================================================\n",
|
|
"Runtime information: \n",
|
|
"Total runtime: 198.828617583 s\n",
|
|
"----------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
|
|
"Initial setup: 0.12195141599999992 s\n",
|
|
"Hamiltonian conversion and XC field extraction: 0.564 s\n",
|
|
"Pair and site datastructure creatrions: 0.207 s\n",
|
|
"k set cration and distribution: 0.036 s\n",
|
|
"Rotating XC potential: 0.358 s\n",
|
|
"Greens function inversion: 197.144 s\n",
|
|
"Calculate energies and magnetic components: 0.398 s\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"if rank == root_node:\n",
|
|
" # Calculate total charge\n",
|
|
" for hamiltonian in hamiltonians:\n",
|
|
" GS = hamiltonian[\"GS\"]\n",
|
|
" traced = np.trace((GS), axis1=1, axis2=2)\n",
|
|
" integral = np.trapz(-1 / np.pi * np.imag(traced * cont.we))\n",
|
|
" print(\"Total charge: \", integral)\n",
|
|
"\n",
|
|
" # iterate over the magnetic entities\n",
|
|
" for tracker, mag_ent in enumerate(magnetic_entities):\n",
|
|
" # iterate over the quantization axes\n",
|
|
" for i, Gii in enumerate(mag_ent[\"Gii\"]):\n",
|
|
" storage = []\n",
|
|
" # iterate over the first and second order local perturbations\n",
|
|
" for Vu1, Vu2 in zip(mag_ent[\"Vu1\"][i], mag_ent[\"Vu2\"][i]):\n",
|
|
" # The Szunyogh-Lichtenstein formula\n",
|
|
" traced = np.trace((Vu2 @ Gii + 0.5 * Gii @ Vu1 @ Gii), axis1=1, axis2=2)\n",
|
|
" # evaluation of the contour integral\n",
|
|
" storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n",
|
|
"\n",
|
|
" # fill up the magnetic entities dictionary with the energies\n",
|
|
" magnetic_entities[tracker][\"energies\"].append(storage)\n",
|
|
" # convert to np array\n",
|
|
" magnetic_entities[tracker][\"energies\"] = np.array(\n",
|
|
" magnetic_entities[tracker][\"energies\"]\n",
|
|
" )\n",
|
|
" print(\"Magnetic entities integrated.\")\n",
|
|
"\n",
|
|
" # iterate over the pairs\n",
|
|
" for tracker, pair in enumerate(pairs):\n",
|
|
" # iterate over the quantization axes\n",
|
|
" for i, (Gij, Gji) in enumerate(zip(pair[\"Gij\"], pair[\"Gji\"])):\n",
|
|
" site_i = magnetic_entities[pair[\"ai\"]]\n",
|
|
" site_j = magnetic_entities[pair[\"aj\"]]\n",
|
|
"\n",
|
|
" storage = []\n",
|
|
" # iterate over the first order local perturbations in all possible orientations for the two sites\n",
|
|
" for Vui in site_i[\"Vu1\"][i]:\n",
|
|
" for Vuj in site_j[\"Vu1\"][i]:\n",
|
|
" # The Szunyogh-Lichtenstein formula\n",
|
|
" traced = np.trace((Vui @ Gij @ Vuj @ Gji), axis1=1, axis2=2)\n",
|
|
" # evaluation of the contour integral\n",
|
|
" storage.append(np.trapz(-1 / np.pi * np.imag(traced * cont.we)))\n",
|
|
" # fill up the pairs dictionary with the energies\n",
|
|
" pairs[tracker][\"energies\"].append(storage)\n",
|
|
" # convert to np array\n",
|
|
" pairs[tracker][\"energies\"] = np.array(pairs[tracker][\"energies\"])\n",
|
|
"\n",
|
|
" print(\"Pairs integrated.\")\n",
|
|
"\n",
|
|
" # calculate magnetic parameters\n",
|
|
" for mag_ent in magnetic_entities:\n",
|
|
" Kxx, Kyy, Kzz, consistency = calculate_anisotropy_tensor(mag_ent)\n",
|
|
" mag_ent[\"K\"] = np.array([Kxx, Kyy, Kzz]) * sisl.unit_convert(\"eV\", \"meV\")\n",
|
|
" mag_ent[\"K_consistency\"] = consistency\n",
|
|
"\n",
|
|
" for pair in pairs:\n",
|
|
" J_iso, J_S, D, J = calculate_exchange_tensor(pair)\n",
|
|
" pair[\"J_iso\"] = J_iso * sisl.unit_convert(\"eV\", \"meV\")\n",
|
|
" pair[\"J_S\"] = J_S * sisl.unit_convert(\"eV\", \"meV\")\n",
|
|
" pair[\"D\"] = D * sisl.unit_convert(\"eV\", \"meV\")\n",
|
|
" pair[\"J\"] = J * sisl.unit_convert(\"eV\", \"meV\")\n",
|
|
"\n",
|
|
" print(\"Magnetic parameters calculated.\")\n",
|
|
"\n",
|
|
" times[\"end_time\"] = timer()\n",
|
|
" print(\n",
|
|
" \"##################################################################### GROGU OUTPUT #############################################################################\"\n",
|
|
" )\n",
|
|
"\n",
|
|
" print_parameters(simulation_parameters)\n",
|
|
" print_atoms_and_pairs(magnetic_entities, pairs)\n",
|
|
" print_runtime_information(times)\n",
|
|
"\n",
|
|
" # remove clutter from magnetic entities and pair information\n",
|
|
" for pair in pairs:\n",
|
|
" del pair[\"Gij\"]\n",
|
|
" del pair[\"Gij_tmp\"]\n",
|
|
" del pair[\"Gji\"]\n",
|
|
" del pair[\"Gji_tmp\"]\n",
|
|
" for mag_ent in magnetic_entities:\n",
|
|
" del mag_ent[\"Gii\"]\n",
|
|
" del mag_ent[\"Gii_tmp\"]\n",
|
|
" del mag_ent[\"Vu1\"]\n",
|
|
" del mag_ent[\"Vu2\"]\n",
|
|
" # create output dictionary with all the relevant data\n",
|
|
" results = dict(\n",
|
|
" parameters=simulation_parameters,\n",
|
|
" magnetic_entities=magnetic_entities,\n",
|
|
" pairs=pairs,\n",
|
|
" runtime=times,\n",
|
|
" )\n",
|
|
" # save dictionary\n",
|
|
" with open(outfile, \"wb\") as output_file:\n",
|
|
" pickle.dump(results, output_file)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"ename": "SyntaxError",
|
|
"evalue": "invalid syntax (3105939143.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 ========================================\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"========================================\n",
|
|
" \n",
|
|
"Atom Angstrom\n",
|
|
"# Label, x y z Sx Sy Sz #Q Lx Ly Lz Jx Jy Jz\n",
|
|
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
|
|
"Te1 1.8955 1.0943 13.1698 -0.0000 0.0000 -0.1543 # 5.9345 -0.0000 0.0000 -0.0537 -0.0000 0.0000 -0.2080 \n",
|
|
"Te2 1.8955 1.0943 7.4002 0.0000 -0.0000 -0.1543 # 5.9345 0.0000 -0.0000 -0.0537 0.0000 -0.0000 -0.2080 \n",
|
|
"Ge3 -0.0000 2.1887 10.2850 0.0000 0.0000 -0.1605 # 3.1927 -0.0000 0.0000 0.0012 0.0000 0.0000 -0.1593 \n",
|
|
"Fe4 -0.0000 0.0000 11.6576 0.0001 -0.0001 2.0466 # 8.3044 0.0000 -0.0000 0.1606 0.0001 -0.0001 2.2072 \n",
|
|
"Fe5 -0.0000 0.0000 8.9124 -0.0001 0.0001 2.0466 # 8.3044 -0.0000 0.0000 0.1606 -0.0001 0.0001 2.2072 \n",
|
|
"Fe6 1.8955 1.0944 10.2850 0.0000 0.0000 1.5824 # 8.3296 -0.0000 -0.0000 0.0520 -0.0000 0.0000 1.6344 \n",
|
|
"==================================================================================================================================\n",
|
|
" \n",
|
|
"Exchange meV\n",
|
|
"--------------------------------------------------------------------------------\n",
|
|
"# at1 at2 i j k # d (Ang)\n",
|
|
"--------------------------------------------------------------------------------\n",
|
|
"Fe4 Fe5 0 0 0 # 2.7452\n",
|
|
"Isotropic -82.0854\n",
|
|
"DMI 0.12557 -0.00082199 6.9668e-08\n",
|
|
"Symmetric-anisotropy -0.60237 -0.83842 -0.00032278 -1.2166e-05 -3.3923e-05\n",
|
|
"--------------------------------------------------------------------------------\n",
|
|
"Fe4 Fe6 0 0 0 # 2.5835\n",
|
|
"Isotropic -41.9627\n",
|
|
"DMI 1.1205 -1.9532 0.0018386\n",
|
|
"Symmetric-anisotropy 0.26007 -0.00013243 0.12977 -0.069979 -0.042066\n",
|
|
"--------------------------------------------------------------------------------\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": ".venv",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 3
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython3",
|
|
"version": "3.9.6"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|