et209d
/
grogu
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
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.
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '21ab0dca84'
${ noResults }
grogu/test.ipynb

706 lines
32 KiB
Raw Blame History

{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"[Daniels-Air:17787] shmem: mmap: an error occurred while determining whether or not /var/folders/yh/dx7xl94n3g52ts3td8qcxjcc0000gn/T//ompi.Daniels-Air.501/jf.0/2288320512/sm_segment.Daniels-Air.501.88650000.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",
"from mpi4py import MPI\n",
"from numpy.linalg import inv\n",
"import warnings\n",
"\n",
"start_time = timer()\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of nodes in the parallel cluster: 1\n"
]
}
],
"source": [
"# this cell mimicks an input file\n",
"fdf = sisl.get_sile(\n",
" \"./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",
"# list of reference directions for around which we calculate the derivatives\n",
"# o is the quantization axis, v and w are two axes perpendicular to it\n",
"# at this moment the user has to supply o,v,w on the input.\n",
"# we can have some default for this\n",
"ref_xcf_orientations = [\n",
" dict(o=np.array([1, 0, 0]), vw=[np.array([0, 1, 0]), np.array([0, 0, 1])]),\n",
" dict(o=np.array([0, 1, 0]), vw=[np.array([1, 0, 0]), np.array([0, 0, 1])]),\n",
" dict(o=np.array([0, 0, 1]), vw=[np.array([1, 0, 0]), np.array([0, 1, 0])]),\n",
"]\n",
"\n",
"# human readable definition of magnetic entities\n",
"magnetic_entities = [\n",
" dict(atom=3, l=2),\n",
" dict(atom=4, l=2),\n",
" dict(atom=5, l=2),\n",
"# dict(atom=[3, 4]),\n",
"]\n",
"\n",
"# pair information\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 = -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(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": 3,
"metadata": {},
"outputs": [],
"source": [
"NO = dh.no # shorthand for number of orbitals in the unit cell\n",
"\n",
"# preprocessing Hamiltonian and overlap matrix elements\n",
"h11 = dh.tocsr(dh.M11r)\n",
"h11 += dh.tocsr(dh.M11i) * 1.0j\n",
"h11 = h11.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
"\n",
"h22 = dh.tocsr(dh.M22r)\n",
"h22 += dh.tocsr(dh.M22i) * 1.0j\n",
"h22 = h22.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
"\n",
"h12 = dh.tocsr(dh.M12r)\n",
"h12 += dh.tocsr(dh.M12i) * 1.0j\n",
"h12 = h12.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
"\n",
"h21 = dh.tocsr(dh.M21r)\n",
"h21 += dh.tocsr(dh.M21i) * 1.0j\n",
"h21 = h21.toarray().reshape(NO, dh.n_s, NO).transpose(0, 2, 1).astype(\"complex128\")\n",
"\n",
"sov = (\n",
" dh.tocsr(dh.S_idx)\n",
" .toarray()\n",
" .reshape(NO, dh.n_s, NO)\n",
" .transpose(0, 2, 1)\n",
" .astype(\"complex128\")\n",
")\n",
"\n",
"\n",
"# Reorganization of Hamiltonian and overlap matrix elements to SPIN BOX representation\n",
"U = np.vstack(\n",
" [np.kron(np.eye(NO, dtype=int), [1, 0]), np.kron(np.eye(NO, dtype=int), [0, 1])]\n",
")\n",
"# This is the permutation that transforms ud1ud2 to u12d12\n",
"# That is this transforms FROM SPIN BOX to ORBITAL BOX => U\n",
"# the inverse transformation is U.T u12d12 to ud1ud2\n",
"# That is FROM ORBITAL BOX to SPIN BOX => U.T\n",
"\n",
"# From now on everything is in SPIN BOX!!\n",
"hh, ss = np.array(\n",
" [\n",
" U.T @ np.block([[h11[:, :, i], h12[:, :, i]], [h21[:, :, i], h22[:, :, i]]]) @ U\n",
" for i in range(dh.lattice.nsc.prod())\n",
" ]\n",
"), np.array(\n",
" [\n",
" U.T\n",
" @ np.block([[sov[:, :, i], sov[:, :, i] * 0], [sov[:, :, i] * 0, sov[:, :, i]]])\n",
" @ U\n",
" for i in range(dh.lattice.nsc.prod())\n",
" ]\n",
")\n",
"\n",
"\n",
"# symmetrizing Hamiltonian and overlap matrix to make them hermitian\n",
"for i in range(dh.lattice.sc_off.shape[0]):\n",
" j = dh.lattice.sc_index(-dh.lattice.sc_off[i])\n",
" h1, h1d = hh[i], hh[j]\n",
" hh[i], hh[j] = (h1 + h1d.T.conj()) / 2, (h1d + h1.T.conj()) / 2\n",
" s1, s1d = ss[i], ss[j]\n",
" ss[i], ss[j] = (s1 + s1d.T.conj()) / 2, (s1d + s1.T.conj()) / 2\n",
"\n",
"# identifying TRS and TRB parts of the Hamiltonian\n",
"TAUY = np.kron(np.eye(NO), tau_y)\n",
"hTR = np.array([TAUY @ hh[i].conj() @ TAUY for i in range(dh.lattice.nsc.prod())])\n",
"hTRS = (hh + hTR) / 2\n",
"hTRB = (hh - hTR) / 2\n",
"\n",
"# extracting the exchange field\n",
"traced = [spin_tracer(hTRB[i]) for i in range(dh.lattice.nsc.prod())] # equation 77\n",
"XCF = np.array(\n",
" [\n",
" np.array([f[\"x\"] for f in traced]),\n",
" np.array([f[\"y\"] for f in traced]),\n",
" np.array([f[\"z\"] for f in traced]),\n",
" ]\n",
") # equation 77\n",
"\n",
"# Check if exchange field has scalar part\n",
"max_xcfs = abs(np.array(np.array([f[\"c\"] for f in traced]))).max()\n",
"if max_xcfs > 1e-12:\n",
" warnings.warn(\n",
" f\"Exchange field has non negligible scalar part. Largest value is {max_xcfs}\"\n",
" )\n",
"\n",
"H_and_XCF_time = timer()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# for every site we have to store 3 Greens function (and the associated _tmp-s) in the 3 reference directions\n",
"for i, mag_ent in enumerate(magnetic_entities):\n",
" parsed = parse_magnetic_entity(dh, **mag_ent) # parse orbital indexes\n",
" magnetic_entities[i][\"orbital_indeces\"] = parsed\n",
" 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(\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(\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(\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": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"k loop: 0%| | 0/400 [00:00<?, ?it/s]"
]
}
],
"source": [
"kset = make_kset(dirs=kdirs, NUMK=kset) # generate k space sampling\n",
"wkset = np.ones(len(kset)) / len(kset) # generate weights for k points\n",
"kpcs = np.array_split(kset, size) # split the k points based on MPI size\n",
"kpcs[root_node] = tqdm(kpcs[root_node], desc='k loop', file=stdout)\n",
"\n",
"k_set_time = timer()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# this will contain all the data 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(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",
"\n",
" for mag_ent in magnetic_entities:\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\"]][\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(\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": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"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",
"k loop: 27%|██▋ | 107/400 [00:29<01:20, 3.62it/s]\n"
]
},
{
"ename": "KeyboardInterrupt",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[7], line 26\u001b[0m\n\u001b[1;32m 24\u001b[0m H \u001b[38;5;241m=\u001b[39m hamiltonian_orientation[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mH\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n\u001b[1;32m 25\u001b[0m HK, SK \u001b[38;5;241m=\u001b[39m hsk(H, ss, dh\u001b[38;5;241m.\u001b[39msc_off, k)\n\u001b[0;32m---> 26\u001b[0m Gk \u001b[38;5;241m=\u001b[39m \u001b[43minv\u001b[49m\u001b[43m(\u001b[49m\u001b[43mSK\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43meran\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mreshape\u001b[49m\u001b[43m(\u001b[49m\u001b[43meset\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mHK\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 28\u001b[0m \u001b[38;5;66;03m# store the Greens function slice of the magnetic entities (for now) based on the on-site projections\u001b[39;00m\n\u001b[1;32m 29\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m mag_ent \u001b[38;5;129;01min\u001b[39;00m magnetic_entities:\n",
"File \u001b[0;32m<__array_function__ internals>:200\u001b[0m, in \u001b[0;36minv\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
"File \u001b[0;32m~/Downloads/grogu/.venv/lib/python3.9/site-packages/numpy/linalg/linalg.py:538\u001b[0m, in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m 536\u001b[0m signature \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mD->D\u001b[39m\u001b[38;5;124m'\u001b[39m \u001b[38;5;28;01mif\u001b[39;00m isComplexType(t) \u001b[38;5;28;01melse\u001b[39;00m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124md->d\u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m 537\u001b[0m extobj \u001b[38;5;241m=\u001b[39m get_linalg_error_extobj(_raise_linalgerror_singular)\n\u001b[0;32m--> 538\u001b[0m ainv \u001b[38;5;241m=\u001b[39m \u001b[43m_umath_linalg\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43minv\u001b[49m\u001b[43m(\u001b[49m\u001b[43ma\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msignature\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43msignature\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mextobj\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mextobj\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 539\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m wrap(ainv\u001b[38;5;241m.\u001b[39mastype(result_t, copy\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m))\n",
"\u001b[0;31mKeyboardInterrupt\u001b[0m: "
]
}
],
"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",
" 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",
" Gk = inv(SK * eran.reshape(eset, 1, 1) - HK)\n",
"\n",
" # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
" for mag_ent in magnetic_entities:\n",
" mag_ent[\"Gii_tmp\"][i] += (\n",
" Gk[:, mag_ent[\"spin_box_indeces\"]][..., mag_ent[\"spin_box_indeces\"]]\n",
" * wk\n",
" )\n",
"\n",
" for pair in pairs:\n",
" # add phase shift based on the cell difference\n",
" phase = np.exp(1j * 2 * np.pi * k @ pair[\"Ruc\"].T)\n",
"\n",
" # get the pair orbital sizes from the magnetic entities\n",
" ai = magnetic_entities[pair[\"ai\"]][\"spin_box_indeces\"]\n",
" aj = magnetic_entities[pair[\"aj\"]][\"spin_box_indeces\"]\n",
"\n",
" # store the Greens function slice of the magnetic entities (for now) based on the on-site projections\n",
" pair[\"Gij_tmp\"][i] += Gk[:, ai][..., aj] * phase * wk\n",
" pair[\"Gji_tmp\"][i] += Gk[:, aj][..., ai] * phase * wk\n",
"\n",
"# summ reduce partial results of mpi nodes\n",
"for i in range(len(hamiltonians)):\n",
" for mag_ent in magnetic_entities:\n",
" comm.Reduce(mag_ent[\"Gii_tmp\"][i], mag_ent[\"Gii\"][i], root=root_node)\n",
"\n",
" for pair in pairs:\n",
" comm.Reduce(pair[\"Gij_tmp\"][i], pair[\"Gij\"][i], root=root_node)\n",
" comm.Reduce(pair[\"Gji_tmp\"][i], pair[\"Gji\"][i], root=root_node)\n",
"\n",
"green_function_inversion_time = timer()\n"
]
},
{
"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",
" 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(print_atomic_indices(pair, magnetic_entities, dh))\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": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"import sisl.viz\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"dh.geometry.plot(axes=\"xy\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"coords = dh.xyz[-3:]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plt.figure(figsize=(15,5))\n",
"plt.subplot(131)\n",
"plt.scatter(coords[:,0], coords[:,2], color=[\"r\", \"g\", \"b\"])\n",
"plt.xlabel(\"x\")\n",
"plt.ylabel(\"z\")\n",
"plt.subplot(132)\n",
"plt.scatter(coords[:,1], coords[:,2], color=[\"r\", \"g\", \"b\"])\n",
"plt.xlabel(\"y\")\n",
"plt.ylabel(\"z\")\n",
"plt.subplot(133)\n",
"plt.scatter(coords[:,0], coords[:,1], color=[\"r\", \"g\", \"b\"])\n",
"plt.xlabel(\"x\")\n",
"plt.ylabel(\"y\")\n",
"print(\"xyz[-3:]: red, green, blue\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"print(calculate_exchange_tensor(pairs[0])[0]) # isotropic should be -82 meV\n",
"print(calculate_exchange_tensor(pairs[1])[0]) # these should all be around -41.9 in the isotropic part\n",
"#print(calculate_exchange_tensor(pairs[2])) # these should all be around -41.9 in the isotropic part\n",
"#print(calculate_exchange_tensor(pairs[3])) # these should all be around -41.9 in the isotropic part\n",
"#print(calculate_exchange_tensor(pairs[4])) # these should all be around -41.9 in the isotropic part\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"These are reasonably converged:\n",
"\n",
"-61.33097171216109\n",
"-60.52198328932686\n",
"-60.51657719027764\n",
"-6.545208546361317\n",
"-6.043716409664797"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# symmetrizing Hamiltonian and overlap matrix to make them hermitian \n",
"# Check if exchange field has scalar part\n",
"# parallel over integrals\n",
"\n",
"========================================\n",
" \n",
"Atom Angstrom\n",
"# Label, x y z Sx Sy Sz #Q Lx Ly Lz Jx Jy Jz\n",
"--------------------------------------------------------------------------------------------------------------------------------------------------------------------------\n",
"Te1 1.8955 1.0943 13.1698 -0.0000 0.0000 -0.1543 # 5.9345 -0.0000 0.0000 -0.0537 -0.0000 0.0000 -0.2080 \n",
"Te2 1.8955 1.0943 7.4002 0.0000 -0.0000 -0.1543 # 5.9345 0.0000 -0.0000 -0.0537 0.0000 -0.0000 -0.2080 \n",
"Ge3 -0.0000 2.1887 10.2850 0.0000 0.0000 -0.1605 # 3.1927 -0.0000 0.0000 0.0012 0.0000 0.0000 -0.1593 \n",
"Fe4 -0.0000 0.0000 11.6576 0.0001 -0.0001 2.0466 # 8.3044 0.0000 -0.0000 0.1606 0.0001 -0.0001 2.2072 \n",
"Fe5 -0.0000 0.0000 8.9124 -0.0001 0.0001 2.0466 # 8.3044 -0.0000 0.0000 0.1606 -0.0001 0.0001 2.2072 \n",
"Fe6 1.8955 1.0944 10.2850 0.0000 0.0000 1.5824 # 8.3296 -0.0000 -0.0000 0.0520 -0.0000 0.0000 1.6344 \n",
"==================================================================================================================================\n",
" \n",
"Exchange meV\n",
"--------------------------------------------------------------------------------\n",
"# at1 at2 i j k # d (Ang)\n",
"--------------------------------------------------------------------------------\n",
"Fe4 Fe5 0 0 0 # 2.7452\n",
"Isotropic -82.0854\n",
"DMI 0.12557 -0.00082199 6.9668e-08\n",
"Symmetric-anisotropy -0.60237 -0.83842 -0.00032278 -1.2166e-05 -3.3923e-05\n",
"--------------------------------------------------------------------------------\n",
"Fe4 Fe6 0 0 0 # 2.5835\n",
"Isotropic -41.9627\n",
"DMI 1.1205 -1.9532 0.0018386\n",
"Symmetric-anisotropy 0.26007 -0.00013243 0.12977 -0.069979 -0.042066\n",
"--------------------------------------------------------------------------------\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": ".venv",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Reference in new issue View Git Blame Copy Permalink