crystal.crystal
1from __future__ import annotations 2from time import time 3from tree.linalg import MatrixTree 4from tree.network import _Mapper, Tree 5# from crystal.polytope import * 6from polytope import * 7from crystal.types import BoxOfPoints 8from varname.core import varname 9from copy import copy, deepcopy 10import numpy as np 11from torch.linalg import eigh 12from torch import empty, conj, multiply, exp, ones, pi, imag, einsum, zeros, tensor, arange, reshape, sum, complex32, complex64, complex128, moveaxis, dstack, eye, cat, eq, logical_not, flip 13from numpy import conj as scalar_conj 14import scipy.special as sp 15import seaborn as sns 16import matplotlib.pyplot as plt 17import torch 18 19COMPLEX=tensor([1.0j]).dtype 20pauli_z = tensor([[1,0],[0,-1]],dtype=COMPLEX) 21pauli_y = tensor([[0,-1.0j],[1.0j,0]],dtype=COMPLEX) 22pauli_x = tensor([[0,1],[1,0]],dtype=COMPLEX) 23pauli_plus = (pauli_x+1.0j*pauli_y)*(1./2.) 24pauli_minus = (pauli_x-1.0j*pauli_y)*(1./2.) 25eye = eye(2,dtype=COMPLEX) 26pauli_vector = dstack([eye, pauli_x, pauli_y, pauli_z]) 27pauli_vector = moveaxis(pauli_vector,-1,0) 28 29def _set_dict(_dict, structure, indices, alt_indices): 30 try: 31 _dict[structure][0].extend(indices) 32 _dict[structure][1].extend(alt_indices) 33 except: 34 _dict[structure] = [indices, alt_indices] 35 36def _set_structure_dict(_dict, structure_dict, parameter, anomalous): 37 for structure, [indices, alt_indices] in structure_dict.items(): 38 if len(indices)>0: 39 try: 40 _dict[parameter] 41 except: 42 _dict[parameter] = {} 43 try: 44 _dict[parameter][structure] 45 except: 46 _dict[parameter][structure] = {"indices": [], "anomalous": anomalous} 47 try: 48 _dict[parameter][structure]["indices"][0].extend(indices) 49 _dict[parameter][structure]["indices"][1].extend(alt_indices) 50 except: 51 _dict[parameter][structure]["indices"] = [indices, alt_indices] 52 53def split_structure_dict(structure_dict): 54 diagonals, offdiagonals = {}, {} 55 for structure, [indices, alt_indices] in structure_dict.items(): 56 if type(indices)==list: 57 indices = tensor(indices) 58 alt_indices = tensor(alt_indices) 59 _bool = eq(indices, alt_indices) 60 _diagonals = [*indices[_bool].tolist()], [*alt_indices[_bool].tolist()] 61 _not_bool = logical_not(_bool) 62 _offdiagonals = [*indices[_not_bool].tolist()], [*alt_indices[_not_bool].tolist()] 63 diagonals[structure] = _diagonals 64 offdiagonals[structure] = _offdiagonals 65 return diagonals, offdiagonals 66 67class Hamiltonian(): 68 69 @property 70 def dof(self): 71 try: 72 return self.__dict__["_dof"] 73 except: 74 self.__dict__["_dof"] = self.n_leaves 75 return self.__dict__["_dof"] 76 77 def solve(self, overwrite_matrix=False): 78 79 if self.is_anomalous(): 80 self.eigenvalues = empty([2*self.dof]) 81 else: 82 self.eigenvalues = empty([self.dof]) 83 84 if overwrite_matrix: 85 (self.eigenvalues,self.eigenvectors) = eigh(self.matrix, out=(self.eigenvalues, self.matrix)) 86 else: 87 (self.eigenvalues,self.eigenvectors) = eigh(self.matrix) 88 89 # if self.is_anomalous(): 90 # self.eigenvalues[:self.dof] = flip(self.eigenvalues[:self.dof], [0]) 91 # self.eigenvectors[:,:self.dof] = flip(self.eigenvectors[:,:self.dof], [1]) 92 93 return self.eigenvalues, self.eigenvectors 94 95 def _check_eigenvectors(self): 96 if not hasattr(self, "eigenvectors"): 97 raise AttributeError("Self has no attribute eigenvectors!") 98 99class StatisticalMechanics(Hamiltonian): 100 101 def fermi_distribution(self, temperature:float=0): 102 if temperature == 0: 103 if self.anomalous: 104 distribution = ones(2*self.dof, dtype=COMPLEX) 105 else: 106 distribution = ones(self.dof, dtype=COMPLEX) 107 # Depopulate negative energy solutions 108 if self.anomalous: 109 ## Make use of BdG symmetry 110 distribution[..., :self.dof] = 0 111 else: 112 distribution[self.eigenvalues > 0] = 0 113 114 else: 115 distribution = 1/exp((self.eigenvalues/temperature)+1) 116 117 return distribution 118 119 def _get_lattice_resolved_indices(self, site=None): 120 """ 121 Returns indices of specified site as an array. 122 The first self.nd dimensions of the array are the dimensions of the lattice 123 124 Examples 125 -------- 126 127 >>> A = Atom(fractional_coordinates=[0, 0]) 128 >>> B = Atom(fractional_coordinates=[0.5, 0]) 129 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 130 >>> lattice.cut(number_of_cells=3, direction=0) 131 >>> lattice.cut(number_of_cells=3, direction=1) 132 >>> lattice.n_spins = 2 133 134 >>> lattice.__set_indices__() 135 >>> lattice._get_lattice_resolved_indices(site=None) 136 [[[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]], [[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]]] 137 >>> lattice._get_lattice_resolved_indices(site=lattice) 138 [[[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]], [[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]]] 139 >>> lattice._get_lattice_resolved_indices(site=A) 140 [[[[0, 1]], [[4, 5]], [[8, 9]]], [[[12, 13]], [[16, 17]], [[20, 21]]], [[[24, 25]], [[28, 29]], [[32, 33]]]] 141 >>> lattice._get_lattice_resolved_indices(site=lattice[:,:,:,"↑"]) 142 [[[[[0]], [[2]]], [[[4]], [[6]]], [[[8]], [[10]]]], [[[[12]], [[14]]], [[[16]], [[18]]], [[[20]], [[22]]]], [[[[24]], [[26]]], [[[28]], [[30]]], [[[32]], [[34]]]]] 143 >>> lattice._get_lattice_resolved_indices(site=lattice[:,:,A,"↑"]) 144 [[[[[0]]], [[[4]]], [[[8]]]], [[[[12]]], [[[16]]], [[[20]]]], [[[[24]]], [[[28]]], [[[32]]]]] 145 >>> lattice._get_lattice_resolved_indices(site=[A, B]) 146 [[[[0, 1], [4, 5]], [[8, 9], [12, 13]], [[16, 17], [20, 21]]], [[[24, 25], [28, 29]], [[32, 33], [2, 3]], [[6, 7], [10, 11]]], [[[14, 15], [18, 19]], [[22, 23], [26, 27]], [[30, 31], [34, 35]]]] 147 148 """ 149 _slice = [slice(None, None, None) for _ in range(self.nd)] 150 if site==None: 151 site = self 152 _type = type(site) 153 if _type==Lattice: 154 return site[*_slice].indices 155 elif _type==list: 156 indices = [self._get_lattice_resolved_indices(site=site) for site in site] 157 _shape = list(np.shape(indices)) 158 height = site[0].height 159 _shape[-height] = len(site) 160 _shape = _shape[1:] 161 indices = np.reshape(indices, _shape) 162 return indices.tolist() 163 elif _type==_Mapper: 164 return site.indices 165 elif _type==Cell: 166 return site.indices 167 else: 168 lattice_depth = self.height 169 depth = lattice_depth - site.height 170 for _ in range(depth-1): 171 _slice.append(slice(None, None, None)) 172 _slice.append(site) 173 cells = self[*_slice] 174 return cells.indices 175 176 def _get_rolled_indices(self, indices, hop_vector:None): 177 """ 178 indices should be output of self._get_lattice_resolved_indices(...) 179 180 Examples 181 -------- 182 183 >>> A = Atom(fractional_coordinates=[0, 0]) 184 >>> B = Atom(fractional_coordinates=[0.5, 0]) 185 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 186 >>> lattice.cut(number_of_cells=3, direction=0) 187 >>> lattice.cut(number_of_cells=3, direction=1) 188 >>> lattice.n_spins = 2 189 190 >>> lattice.__set_indices__() 191 192 >>> indices = [[[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]], [[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]]] 193 >>> lattice._get_rolled_indices(indices, [1,0]) 194 [[[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]], [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]] 195 196 """ 197 if len(hop_vector)!=self.nd: 198 raise ValueError("Hop must have same length as model dimensions!") 199 for axis in range(self.nd): 200 shift = hop_vector[axis] 201 if self.pbc[axis]: 202 indices = np.roll(indices, shift=-shift, axis=axis) 203 else: 204 indices = np.add(indices, shift) 205 dimension = self.nd[axis] 206 indices[indices >= dimension] = None 207 return indices.tolist() 208 209 def _get_spin_matrix_indices(self, indices, spin_vector, structure=1): 210 """ 211 indices should be output of self._get_lattice_resolved_indices(...) or self._get_rolled_indices(...) 212 213 Examples 214 -------- 215 216 >>> A = Atom(fractional_coordinates=[0, 0]) 217 >>> B = Atom(fractional_coordinates=[0.5, 0]) 218 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 219 >>> lattice.cut(number_of_cells=3, direction=0) 220 >>> lattice.cut(number_of_cells=3, direction=1) 221 >>> lattice.n_spins = 2 222 223 >>> lattice.__set_indices__() 224 225 >>> indices = [[[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]], [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]] 226 >>> lattice._get_spin_matrix_indices(indices, structure=-1) 227 {(-1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]} 228 >>> lattice._get_spin_matrix_indices(indices, [0,1,0,0]) 229 {(1+0j): [[12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 0, 2, 4, 6, 8, 10, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 1, 3, 5, 7, 9, 11], [13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 1, 3, 5, 7, 9, 11, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 0, 2, 4, 6, 8, 10]]} 230 231 """ 232 233 indices = list(flatten(indices)) 234 235 if type(spin_vector)==list: 236 spin_vector = tensor(spin_vector, dtype=COMPLEX) 237 _pauli_vector = einsum('i,ijk->jk', spin_vector, pauli_vector) 238 structure_dict = {} 239 for i in range(2): 240 for j in range(2): 241 _structure = complex(_pauli_vector[i,j] * structure) 242 if abs(_structure)!=0: 243 _set_dict(structure_dict, structure, indices[i::2], indices[j::2]) 244 return structure_dict 245 246 def get_diagonal_indices(self, site, structure, pauli_vector=None): 247 """ 248 Examples 249 -------- 250 251 >>> A = Atom(fractional_coordinates=[0, 0]) 252 >>> B = Atom(fractional_coordinates=[0.5, 0]) 253 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 254 >>> lattice.cut(number_of_cells=3, direction=0) 255 >>> lattice.cut(number_of_cells=3, direction=1) 256 >>> lattice.n_spins = 2 257 258 >>> lattice.__set_indices__() 259 260 >>> lattice.get_diagonal_indices(site=A, structure=-1) 261 {(-1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]} 262 263 >>> lattice.get_diagonal_indices(site=A, structure=-1, pauli_vector=[0,0,0,1]) 264 {(1+0j): [[0, 4, 8, 12, 16, 20, 24, 28, 32, 1, 5, 9, 13, 17, 21, 25, 29, 33], [1, 5, 9, 13, 17, 21, 25, 29, 33, 0, 4, 8, 12, 16, 20, 24, 28, 32]]} 265 """ 266 267 indices = self._get_lattice_resolved_indices(site) 268 269 if type(pauli_vector)!=type(None): 270 structure = self._get_spin_matrix_indices(indices, pauli_vector, structure) 271 else: 272 indices = flatten(list(indices)) 273 structure = {complex(structure): [indices, indices]} 274 275 return structure 276 277 def get_offdiagonal_indices(self, site, neighbour, hop_vector=None, trs=True, structure=1, pauli_vector=None): 278 """ 279 Examples 280 -------- 281 282 >>> lattice.get_offdiagonal_indices(site=A, neighbour=B, hop_vector=[1,0], trs=True, structure=1j) 283 {1j: [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11]], -1j: [[14, 15, 18, 19, 22, 284 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]} 285 >>> lattice.get_offdiagonal_indices(site=A, neighbour=B, hop_vector=[1,0], trs=True, structure=1) 286 {(1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11], [14, 15, 18, 19, 22, 23, 287 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11, 0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 1 288 0, 11]]} 289 290 TODO add parameter to flip spin with Pauli vector 291 292 """ 293 if type(pauli_vector)!=type(None): 294 raise ValueError("Pauli vector not yet implemented!") 295 296 indices = self._get_lattice_resolved_indices(site) 297 alt_indices = self._get_lattice_resolved_indices(neighbour) 298 299 if type(hop_vector)!=type(None): 300 alt_indices = self._get_rolled_indices(alt_indices, hop_vector) 301 302 indices = flatten(list(indices)) 303 alt_indices = flatten(list(alt_indices)) 304 structure_dict = {} 305 if trs: 306 _indices = deepcopy(indices) 307 _alt_indices = deepcopy(alt_indices) 308 309 _set_dict(structure_dict, structure, indices, alt_indices) 310 311 if trs: 312 structure = scalar_conj(structure) 313 _set_dict(structure_dict, structure, _alt_indices, _indices) 314 315 return structure_dict 316 317 def _init_sparse_matrices(self): 318 319 self.__set_indices__() 320 321 self.anomalous = False 322 self.interacting = False 323 324 self.hamiltonian_sparse_matrix = {} 325 self.interaction_sparse_matrix = {} 326 self.hartree_sparse_matrix = {} 327 self.hartree_sparse_matrix = {} 328 self.fock_sparse_matrix = {} 329 self.gorkov_sparse_matrix = {} 330 self._sparse_matrices = [self.hamiltonian_sparse_matrix, self.interaction_sparse_matrix, self.hartree_sparse_matrix, self.fock_sparse_matrix, self.gorkov_sparse_matrix] 331 332 for parameter, list_of_parameters in self._diagonal.items(): 333 # Unpack dictionary 334 interaction = list_of_parameters["interaction"] 335 if interaction: 336 self.interacting = True 337 renormalise = list_of_parameters["renormalise"] 338 for parameters in list_of_parameters["parameters"]: 339 # Unpack dictionary 340 site = parameters["site"] 341 self.__dict__[parameter] = site.__dict__[parameter] 342 anomalous = parameters["anomalous"] 343 if anomalous: 344 self.anomalous = True 345 structure = parameters["structure"] 346 pauli_vector = parameters["pauli_vector"] 347 spin_orbital_matrix = parameters["spin_orbital_matrix"] 348 # Get {structure: [indices, alt indices]} 349 structure_dict = self.get_diagonal_indices(site=site, structure=structure, pauli_vector=pauli_vector) 350 if interaction: 351 _set_structure_dict(self.interaction_sparse_matrix, structure_dict, parameter, anomalous=False) 352 elif renormalise: 353 if anomalous: 354 _set_structure_dict(self.gorkov_sparse_matrix, structure_dict, parameter, anomalous) 355 else: 356 diagonals_dict, offdiagonals_dict = split_structure_dict(structure_dict) 357 _set_structure_dict(self.hartree_sparse_matrix, diagonals_dict, parameter, anomalous) 358 _set_structure_dict(self.fock_sparse_matrix, offdiagonals_dict, parameter, anomalous) 359 else: 360 _set_structure_dict(self.hamiltonian_sparse_matrix, structure_dict, parameter, anomalous) 361 362 for parameter, list_of_parameters in self._offdiagonal.items(): 363 # Unpack dictionary 364 interaction = list_of_parameters["interaction"] 365 if interaction: 366 self.interacting = True 367 renormalise = list_of_parameters["renormalise"] 368 for parameters in list_of_parameters["parameters"]: 369 # Unpack dictionary 370 site = parameters["site"] 371 self.__dict__[parameter] = site.__dict__[parameter] 372 neighbour = parameters["neighbour"] 373 hop_vector = parameters["vector"] 374 anomalous = parameters["anomalous"] 375 if anomalous: 376 self.anomalous = True 377 structure = parameters["structure"] 378 pauli_vector = parameters["pauli_vector"] 379 trs = parameters["trs"] 380 spin_orbital_matrix = parameters["spin_orbital_matrix"] 381 # Get {structure: [indices, alt indices]} 382 structure_dict = self.get_offdiagonal_indices(site=site, neighbour=neighbour, hop_vector=hop_vector, trs=trs, structure=structure, pauli_vector=pauli_vector) 383 if interaction: 384 _set_structure_dict(self.interaction_sparse_matrix, structure_dict, parameter, anomalous=False) 385 elif renormalise: 386 if anomalous: 387 _set_structure_dict(self.gorkov_sparse_matrix, structure_dict, parameter, anomalous) 388 else: 389 _set_structure_dict(self.hartree_sparse_matrix, structure_dict, parameter, anomalous) 390 _set_structure_dict(self.fock_sparse_matrix, structure_dict, parameter, anomalous) 391 raise ValueError("TODO: Split pauli_vector into diagonal and offdiagonal components") 392 else: 393 _set_structure_dict(self.hamiltonian_sparse_matrix, structure_dict, parameter, anomalous) 394 395 for _dict in self._sparse_matrices: 396 for key, items in _dict.items(): 397 value = self.__dict__[key] 398 for structure, map in items.items(): 399 map["update"] = value 400 items[structure]["indices"] = tensor(items[structure]["indices"]) 401 402 def _update_matrix(self): 403 404 if not hasattr(self, "_matrix"): 405 406 self._init_sparse_matrices() 407 408 if self.anomalous: 409 self.__dict__["_matrix"] = zeros([2*self.dof, 2*self.dof], dtype=COMPLEX) 410 else: 411 self.__dict__["_matrix"] = zeros([self.dof, self.dof], dtype=COMPLEX) 412 413 if self.interacting: 414 self.__dict__["interaction"] = zeros([self.dof, self.dof], dtype=COMPLEX) 415 416 for key, items in self.hamiltonian_sparse_matrix.items(): 417 value = self.__dict__[key] 418 for structure, map in items.items(): 419 [indices, alt_indices] = map["indices"] 420 update = map["update"] 421 anomalous = map["anomalous"] 422 if update != 0: 423 value = map["update"] 424 map["update"] = 0 425 _value = value * structure 426 if self.anomalous: 427 # Creates Hamiltonian of the following type: [[H0, Δ], [-H0, Δ*]] 428 if anomalous: 429 self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += _value 430 self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += scalar_conj(_value) 431 else: 432 self.__dict__["_matrix"][[indices, alt_indices]] += _value 433 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= _value 434 else: 435 self.__dict__["_matrix"][[indices, alt_indices]] += _value 436 437 for key, items in self.interaction_sparse_matrix.items(): 438 value = self.__dict__[key] 439 for structure, map in items.items(): 440 [indices, alt_indices] = map["indices"] 441 update = map["update"] 442 anomalous = map["anomalous"] 443 if update != 0: 444 value = map["update"] 445 map["update"] = 0 446 _value = value * structure 447 self.__dict__["interaction"][[indices, alt_indices]] += _value 448 449 for key, items in self.hartree_sparse_matrix.items(): 450 value = self.__dict__[key] 451 for structure, map in items.items(): 452 [indices, alt_indices] = map["indices"] 453 update = map["update"] 454 anomalous = map["anomalous"] 455 if update != 0: 456 value = map["update"] 457 map["update"] = 0 458 _value = value * structure 459 self.__dict__["_matrix"][[indices, alt_indices]] += _value 460 if self.anomalous: 461 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= _value 462 map["values"] = tensor([_value for _ in range(len(indices))]) 463 # elif "values" in map.keys(): 464 # self.__dict__["_matrix"][[indices, alt_indices]] += map["values"] 465 # if self.anomalous: 466 # self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= map["values"] 467 468 for key, items in self.fock_sparse_matrix.items(): 469 value = self.__dict__[key] 470 for structure, map in items.items(): 471 [indices, alt_indices] = map["indices"] 472 update = map["update"] 473 anomalous = map["anomalous"] 474 if update != 0: 475 value = map["update"] 476 map["update"] = 0 477 _value = value * structure 478 self.__dict__["_matrix"][[indices, alt_indices]] -= _value 479 if self.anomalous: 480 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] += _value 481 map["values"] = tensor([_value for _ in range(len(indices))]) 482 # elif "values" in map.keys(): 483 # self.__dict__["_matrix"][[indices, alt_indices]] -= map["values"] 484 # if self.anomalous: 485 # self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] += map["values"] 486 487 for key, items in self.gorkov_sparse_matrix.items(): 488 value = self.__dict__[key] 489 for structure, map in items.items(): 490 [indices, alt_indices] = map["indices"] 491 update = map["update"] 492 anomalous = map["anomalous"] 493 if update != 0: 494 value = map["update"] 495 map["update"] = 0 496 _value = value * structure 497 self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += _value 498 self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += scalar_conj(_value) 499 map["values"] = tensor([_value for _ in range(len(indices))]) 500 # elif "values" in map.keys(): 501 # self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += map["values"] 502 # self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += conj(map["values"]) 503 504 def extract_mean_fields(self): 505 506 # print("matrix:") 507 # print(self._matrix) 508 509 if not hasattr(self, "eigenvectors"): 510 raise Error("Run self.solve() before calling this function!") 511 for items in self.gorkov_sparse_matrix.values(): 512 for map in items.values(): 513 [indices, alt_indices] = map["indices"] 514 # TODO: add Fermi function 515 _alt_indices = alt_indices + self.dof 516 if "values" not in map: 517 map["values"] = 0 518 previous_values = deepcopy(map["values"]) 519 # map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, :self.dof], flip(self.eigenvectors[_alt_indices, self.dof:], [2])) 520 map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, :self.dof], self.eigenvectors[_alt_indices,: self.dof]) 521 values = map["values"] - previous_values 522 # print("##### gorkov #######") 523 # print(map["values"]) 524 # print("#################") 525 self.__dict__["_matrix"][[indices, _alt_indices]] += values 526 self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += conj(values) 527 528 for items in self.fock_sparse_matrix.values(): 529 for map in items.values(): 530 [indices, alt_indices] = map["indices"] 531 if "values" not in map: 532 map["values"] = 0 533 previous_values = deepcopy(map["values"]) 534 # TODO: add Fermi function 535 map["values"] = einsum('i,in,in->i', -self.interaction[indices,alt_indices], self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:])) 536 values = map["values"] - previous_values 537 self.__dict__["_matrix"][[indices, alt_indices]] -= values # Swap sign for anticommutativity!!! 538 if self.anomalous: 539 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] += values 540 541 for items in self.hartree_sparse_matrix.values(): 542 for map in items.values(): 543 [indices, alt_indices] = map["indices"] 544 # TODO: add Fermi function 545 if "values" not in map: 546 map["values"] = 0 547 previous_values = deepcopy(map["values"]) 548 map["values"] = einsum('ij,in,jn->i', self.interaction, self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:])) 549 # print("##### hartree #######") 550 # print(previous_values) 551 # print(map["values"]) 552 # print("#################") 553 # print(self._matrix) 554 values = map["values"] - previous_values 555 self.__dict__["_matrix"][[indices, alt_indices]] += values 556 if self.anomalous: 557 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= values 558 559 # print(self.hartree_sparse_matrix) 560 # print(self.fock_sparse_matrix) 561 # print(self.gorkov_sparse_matrix) 562 print(self._matrix) 563 564 # Normal: 565 # map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:])) 566 567 # def mean_field_step() 568 # for key, items in self.interaction_sparse_matrix.items(): 569 # value = self.__dict__[key] 570 # for structure, map in items.items(): 571 # [indices, alt_indices] = map["indices"] 572 # update = map["update"] 573 # anomalous = map["anomalous"] 574 # if update != 0: 575 # value = map["update"] 576 # map["update"] = 0 577 # _value = value * structure 578 # self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += _value 579 # self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += scalar_conj(_value) 580 # map["values"] = tensor([_value for _ in range(len(indices))]) 581 # else: 582 # self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += map["values"] 583 # self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += conj(map["values"]) 584 585 ######################################################## 586 587 588 589 590 # def _get_indices(self, sites=None, hop=None, spin_flip=False, anomalous=False): 591 # _slice = [slice(None, None, None) for _ in range(self.nd)] 592 # cells = self.cells[*_slice] 593 594 # if type(hop)!=type(None): 595 # if len(hop)!=self.nd: 596 # raise ValueError("Hop must have same length as model dimensions!") 597 # for axis in range(self.nd): 598 # shift = hop[axis] 599 # cells = np.roll(cells, shift=-shift, axis=axis) 600 601 # cells = _Mapper(cells.tolist()) 602 603 # if type(sites)==type(None): 604 # sites = self.cell 605 606 # index = [0 for _ in range(self.nd)] 607 # index.append(0) 608 # test_cell = list(flatten(cells[*index])) 609 # sites_type = type(sites[0]) 610 # for site in sites: 611 # if type(site)!=sites_type: 612 # raise ValueError("All sites in sites must have same type!") 613 614 # def recursive(test_cell): 615 # test_cell = test_cell[0] 616 # structure_type = type(test_cell) 617 # if structure_type!=sites_type: 618 # _slice.append(slice(None, None, None)) 619 # test_cell = recursive(test_cell) 620 # return test_cell 621 622 # test_cell = recursive(test_cell) 623 624 # _slice.append(sites) 625 626 # if self.n_spins == 2: 627 # spin_flip = [0,1] 628 # if spin_flip: 629 # spin_flip = [1,0] 630 631 # _slice.append(spin_flip) 632 633 # cells = cells[*_slice] 634 635 # indices = cells.indices 636 # shape = np.shape(indices)[:-1] 637 # indices = flatten(cells.indices) 638 # if anomalous and self.anomalous: 639 # indices += self.n_leaves 640 641 # return list(indices), shape 642 643 def _slice_eigenvector(self, sites=None, hop=None, spin_flip=False, anomalous=False): 644 indices, shape = self._get_indices(sites=sites, hop=hop, spin_flip=spin_flip, anomalous=anomalous) 645 return self.eigenvectors[indices], shape 646 647 def density_matrix(self, temperature:float=0, sites=None, sites_to=None, hop=None, spin_flip=False, anomalous=False): 648 649 if type(sites)==type(None) and type(sites_to)!=type(None): 650 raise ValueError("If sites = None, then sites_to must also equal None!") 651 elif type(sites)!=type(None) and type(sites_to)==type(None): 652 sites_to = sites 653 elif type(sites)==type(None) and type(sites_to)==type(None): 654 sites = self.cell 655 sites_to = self.cell 656 elif len(sites)!=len(sites_to): 657 raise ValueError("sites and sites_to must have the same length!") 658 Ψ, shape = self._slice_eigenvector(sites=sites) 659 Ψc, shape = self._slice_eigenvector(sites=sites_to, hop=hop, spin_flip=spin_flip, anomalous=anomalous) 660 Ψc = conj(Ψc) 661 dm = multiply(Ψ, Ψc) 662 shape = shape + (Ψ.shape[-1],) 663 dm = reshape(dm, shape) 664 665 # if temperature==0 and not self.anomalous: 666 # return dm 667 # else: 668 return einsum('...n,n->...n', dm, self.fermi_distribution(temperature)) 669 670 def greens_function(self, energy:float or list, resolution:float, temperature:float=0, sites=None, sites_to=None, hop=None, spin_flip=False, anomalous=False, spin_polarised=False): 671 kwargs = {"temperature": temperature, "sites": sites, "sites_to": sites_to, "hop": hop, "spin_flip": spin_flip, "anomalous": anomalous} 672 dm = self.density_matrix(**kwargs) 673 try: 674 energy = energy[:, None] 675 poles = 1/(energy + 1.0j*resolution - self.eigenvalues) 676 ΨΨ = einsum('...n, en -> ...e', dm, poles) 677 except: 678 poles = 1/(energy + 1.0j*resolution - self.eigenvalues) 679 ΨΨ = einsum('...n, n -> ...', dm, poles) 680 if self.n_spins == 2 and not spin_polarised: 681 ΨΨ = sum(ΨΨ, axis=-1) 682 if (ΨΨ.shape)[-1] == 1: 683 ΨΨ = sum(ΨΨ, axis=-1) 684 return -1/pi*imag(ΨΨ) 685 686class BdGPatch(StatisticalMechanics): 687 688 def __setattr__(self, name, value, update = True): 689 """Set attributes down the lineage""" 690 if name=="n_spins": 691 self._set_n_spins(value) 692 self.__dict__[f"{name}"] = value 693 else: 694 try: 695 previous_value = self.__dict__[f"{name}"] 696 for _dict in self._sparse_matrices: 697 items = _dict[name] 698 for structure in items.keys(): 699 items[structure]["update"] = value - previous_value 700 except: 701 pass 702 self.__dict__[f"{name}"] = value 703 [item.__setattr__(name, value, update = False) for item in self] 704 705 # @property 706 # def spins(self): 707 # return _Mapper(self) 708 709 def _set_n_spins(self, n_spins): 710 if n_spins<1 or n_spins>2: 711 raise ValueError("n_spins must be 1 or 2") 712 try: 713 if self.__dict__["_n_spins"]==2 and n_spins==1: 714 self.clear() 715 [leaf.parent.clear() for leaf in self.leaves] 716 except: 717 pass 718 if n_spins==2: 719 [leaf.extend([Spin('↑'), Spin('↓')]) for leaf in self.leaves] 720 721 def add(self, structure=1, pauli_vector:list=None, spin_orbital_matrix=None, anomalous:bool=False, interaction=False, renormalise=False, **kwarg): 722 """Add vertex term as a parameter described by a key-value pair""" 723 724 parameters = {"site": self, "anomalous": anomalous, "structure": structure, "pauli_vector": pauli_vector, "spin_orbital_matrix": spin_orbital_matrix} 725 726 if len(kwarg)==0: 727 raise ValueError("Requires a kwarg!") 728 729 if len(kwarg)>1: 730 raise ValueError("Accepts only one kwarg!") 731 732 key, value = list(kwarg.items())[0] 733 734 if type(value)!=float and type(value)!=int: 735 raise ValueError("Value must be real! (In order for the symmetries to apply correctly, put the imaginary number in structure and the magnitude in the value)") 736 737 if not key in self._diagonal: 738 self._diagonal[key] = {} 739 740 if not "parameters" in self._diagonal[key].keys(): 741 self._diagonal[key]["parameters"] = [] 742 743 self._diagonal[key]["parameters"].append(parameters) 744 self._diagonal[key]["interaction"] = interaction 745 self._diagonal[key]["renormalise"] = renormalise 746 747 self.__setattr__(key, value) 748 749 def hop(self, neighbour: Tree=None, structure=1, pauli_vector:list=None, spin_orbital_matrix=None, vector=None, trs:bool=True, anomalous:bool=False, interaction=False, renormalise=False, **kwarg): 750 """Add directed edge as a neighbour and a parameter (described by a 751 key-value pair) 752 """ 753 parameters = {"site": self, "neighbour": neighbour, "vector": vector, "trs": trs, "anomalous": anomalous, "structure": structure, "pauli_vector": pauli_vector, "spin_orbital_matrix": spin_orbital_matrix} 754 755 if len(kwarg)==0: 756 raise ValueError("Requires a kwarg!") 757 758 if len(kwarg)>1: 759 raise ValueError("Accepts only one kwarg!") 760 761 if type(neighbour)!=type(None): 762 if self.n_leaves!=neighbour.n_leaves: 763 raise ValueError("Self and neighbour have different number of leaves. Hopping is not well-defined!") 764 765 key, value = list(kwarg.items())[0] 766 767 if not key in self._offdiagonal: 768 self._offdiagonal[key] = {} 769 770 if not "parameters" in self._offdiagonal[key].keys(): 771 self._offdiagonal[key]["parameters"] = [] 772 773 self._offdiagonal[key]["parameters"].append(parameters) 774 self._offdiagonal[key]["interaction"] = interaction 775 self._offdiagonal[key]["renormalise"] = renormalise 776 777 self.__setattr__(key, value) 778 779 def is_anomalous(self): 780 if not hasattr(self, "anomalous"): 781 diagonals, offdiagonals = self._diagonal, self._offdiagonal 782 783 self.anomalous = False 784 for items in list(diagonals.values())+list(offdiagonals.values()): 785 for item in items["parameters"]: 786 if item["anomalous"]: 787 self.anomalous = True 788 break 789 790 return self.anomalous 791 792 # def has_interactions(self): 793 # if not hasattr(self, "interacting"): 794 # diagonals, offdiagonals = self._diagonal, self._offdiagonal 795 796 # self.interacting = False 797 # for items in list(diagonals.values())+list(offdiagonals.values()): 798 # if items["interaction"]: 799 # self.interacting = True 800 # break 801 # return self.interacting 802 803 def has_renormalisation_fields(self): 804 if not hasattr(self, "renormalise"): 805 diagonals, offdiagonals = self._diagonal, self._offdiagonal 806 807 self.renormalise = False 808 for items in list(diagonals.values())+list(offdiagonals.values()): 809 if items["renormalise"]: 810 self.renormalise = True 811 break 812 return self.renormalise 813 814 # def _get_spin_matrix_indices(self): 815 # if self.height == 1: 816 # _indices = self.indices 817 # _indices = [[x, y] for y in _indices for x in _indices] 818 # return _indices 819 # elif self.height>1: 820 # _indices = [] 821 # [_indices.extend(site._get_spin_matrix_indices()) for site in self] 822 # print(_indices) 823 # exit() 824 # return _indices 825 826 def read_value(self, list_of_parameters): 827 value = list_of_parameters["value"] 828 try: 829 value -= list_of_parameters["previous_value"] 830 except: 831 pass 832 return value 833 834 def _get_spin_indices(self, pauli_vector): 835 cache = {} 836 if type(pauli_vector)==list: 837 pauli_vector = tensor(pauli_vector, dtype=COMPLEX) 838 pauli_vector = einsum('i,ijk->jk', pauli_vector, pauli_vec) 839 if self.height>0: 840 spin_matrix_indices = self._get_spin_matrix_indices() 841 pauli_vector = pauli_vector.flatten() 842 for index in range(len(spin_matrix_indices)): 843 value = complex(pauli_vector[index%4]) 844 if abs(value)!=0: 845 _indices = spin_matrix_indices[index] 846 try: 847 cache[value][0].append(_indices[0]) 848 cache[value][1].append(_indices[1]) 849 except: 850 cache[value] = [[_indices[0]], [_indices[1]]] 851 return cache 852 853 def anomalous_cache(self, cache, anomalous): 854 if anomalous: 855 _cache = cache.clone().detach() 856 _cache[1,:] += self.n_leaves 857 _alt_cache = cache.clone().detach() 858 _alt_cache[0,:] += self.n_leaves 859 cache = [_cache, _alt_cache, anomalous] 860 else: 861 cache = [cache, cache + self.n_leaves, anomalous] 862 return cache 863 864 #def read_indices(self, list_of_parameters): 865 # try: 866 # return list_of_parameters["cache"] 867 # except: 868 # pass 869 870 # def set_cache(cache, structure, value, anomalous): 871 # try: 872 # cache[structure] = cat([cache[structure], value], axis=-1) 873 # except: 874 # cache[structure] = value 875 876 # if self.is_anomalous(): 877 # cache[structure] = self.anomalous_cache(cache[structure], anomalous) 878 879 # cache = {} 880 881 # for parameters in list_of_parameters["parameters"]: 882 # site = parameters["site"] 883 # structure = parameters["structure"] 884 # pauli_vector = parameters["pauli_vector"] 885 # spin_orbital_matrix = parameters["spin_orbital_matrix"] 886 # anomalous = parameters["anomalous"] 887 # try: 888 # vector = parameters["vector"] 889 # except: 890 # vector = None 891 # try: 892 # neighbour = parameters["neighbour"] 893 # except: 894 # neighbour = None 895 # try: 896 # trs = parameters["trs"] 897 # except: 898 # trs = False 899 900 # ################################################################ 901 # if type(vector)==type(None) and type(neighbour)==type(None): 902 # if type(spin_orbital_matrix)!=type(None): 903 # # indices = site._get_spin_orbital_indices(pauli_vector) 904 # raise ValueError("Not yet implemented!") 905 906 # elif type(pauli_vector)!=type(None): 907 # pauli_vector = tensor(pauli_vector, dtype=COMPLEX) 908 # pauli_vector *= structure 909 # for key, value in site._get_spin_indices(pauli_vector).items(): 910 # structure = key 911 # value = tensor(value) 912 # set_cache(cache, structure, value, anomalous) 913 # else: 914 # value = tensor([site.indices, site.indices]) 915 # set_cache(cache, structure, value, anomalous) 916 # ################################################################ 917 # else: 918 # if neighbour == None: 919 # neighbour = site 920 # vector = parameters["vector"] 921 922 # _slice = [slice(None, None, None) for _ in range(self.nd)] 923 # _neighbour_slice = deepcopy(_slice) 924 925 # if type(site)==Cell: 926 # cells = _Mapper([site]) 927 # neighbour_cells = _Mapper([neighbour]) 928 # elif type(site)!=Lattice: 929 # if site.depth>1: 930 # [_slice.append(slice(None,None,None)) for _ in range(site.height-1)] 931 # [_neighbour_slice.append(slice(None,None,None)) for _ in range(neighbour.height-1)] 932 # _slice.append(site) 933 # _neighbour_slice.append(neighbour) 934 # cells = self[*_slice] 935 # neighbour_cells = self[*_neighbour_slice] 936 # else: 937 # cells = self.cells 938 # neighbour_cells = self.cells 939 940 # if type(vector)!=type(None): 941 # _shape = self.cells.shape 942 # _neighbour_cells = np.empty(_shape, dtype=MatrixTree) 943 # nslice = [slice(None, None, None) for _ in range(self.nd)] 944 # _neighbour_cells[*nslice] = neighbour_cells 945 # for axis in range(self.nd): 946 # shift = vector[axis] 947 # _neighbour_cells = np.roll(_neighbour_cells, shift=shift, axis=axis) 948 # neighbour_cells = _neighbour_cells 949 # # else: 950 # # nslice = [slice(None, None, None) for _ in range(self.nd)] 951 # # cells = _cells[*nslice] 952 # # neighbour_cells = _neighbour_cells[*nslice] 953 954 # if type(pauli_vector)!=type(None): 955 956 # # Neighbour indices: 957 # _indices = [site._get_spin_indices(pauli_vector) for site in flatten(self.cells.tolist())] 958 # _neighbour_indices = [site._get_spin_indices(pauli_vector) for site in flatten(neighbour_cells)] 959 960 # keys = list(_indices[0].keys()) 961 962 # def func(x, y, key): 963 # try: 964 # x[key][0].extend(y[key][0]) 965 # x[key][1].extend(y[key][1]) 966 # except: 967 # x[key] = y[key] 968 969 # __indices = {} 970 # __neighbour_indices = {} 971 # [func(__indices, _index, key) for key in keys for _index in _indices] 972 # [func(__neighbour_indices, _index, key) for key in keys for _index in _neighbour_indices] 973 974 # indices = {} 975 # for key in __indices.keys(): 976 # indices[key] = [__indices[key][0], __neighbour_indices[key][1]] 977 978 # for key, value in indices.items(): 979 # structure = key 980 # value = tensor(value) 981 # set_cache(cache, structure, value, anomalous) 982 983 # else: 984 # try: 985 # indices = flatten(list(_Mapper(cells.tolist()).indices)) 986 # neighbour_indices = flatten(list(_Mapper(neighbour_cells.tolist()).indices)) 987 # except: 988 # indices = flatten(list(cells.indices)) 989 # neighbour_indices = flatten(list(neighbour.indices)) 990 # value = tensor([indices, neighbour_indices]) 991 # set_cache(cache, structure, value, anomalous) 992 # if trs: 993 # value = tensor([neighbour_indices, indices]) 994 # set_cache(cache, scalar_conj(structure), value, anomalous) 995 # # name = neighbour.__name__ 996 # # if type(pauli_vector)==type(None): 997 # # if site.depth == 0: 998 # # neighbour_indices = flatten([cell.indices for cell in cells]) 999 # # elif site.depth == 1: 1000 # # neighbour_indices = flatten([cell[name].indices for cell in cells]) 1001 # # elif site.depth > 1: 1002 # # _slice = [slice(None, None, None) for _ in range(site.depth-1)] 1003 # # _slice.append(name) 1004 # # neighbour_indices = flatten([cell[*_slice].indices for cell in cells]) 1005 1006 # # if type(spin_orbital_matrix)!=type(None): 1007 # # # indices = site._get_spin_orbital_indices(pauli_vector) 1008 # # raise ValueError("Not yet implemented!") 1009 1010 # # elif type(pauli_vector)!=type(None): 1011 # # print(vector) 1012 # # exit() 1013 # # pauli_vector = tensor(pauli_vector, dtype=COMPLEX) 1014 # # pauli_vector *= structure 1015 # # for key, value in site._get_spin_indices(pauli_vector).items(): 1016 # # structure = key 1017 # # value = tensor(value) 1018 # # else: 1019 # # indices = site.indices 1020 # # print(cells) 1021 # # exit() 1022 # # exit() 1023 # # print(neighbour_indices) 1024 # # exit() 1025 # # print(self) 1026 # # print(site) 1027 # # exit() 1028 1029 # list_of_parameters["cache"] = cache 1030 1031 # return list_of_parameters["cache"] 1032 1033 @property 1034 def matrix(self): 1035 return self.get_matrix() 1036 1037 # @property 1038 def get_matrix(self): 1039 """ 1040 Construct the matrix representation of the tree 1041 1042 >>> A = MatrixTree() 1043 >>> B = MatrixTree() 1044 >>> matrix = MatrixTree(A, B) 1045 >>> matrix.add(μ=0.6) 1046 >>> A.hop(B, t=1) 1047 >>> matrix.matrix 1048 array([[0.6, 1. ], 1049 [0. , 0.6]]) 1050 1051 >>> A.t = 3 1052 >>> matrix.matrix 1053 array([[0.6, 3. ], 1054 [0. , 0.6]]) 1055 1056 >>> A.t *= 3 1057 >>> matrix.matrix 1058 array([[0.6, 9. ], 1059 [0. , 0.6]]) 1060 """ 1061 1062 self._update_matrix() 1063 1064 return self._matrix 1065 1066 def self_consistent_step(self): 1067 1068 matrix = self.get_matrix() 1069 1070 print(matrix) 1071 1072 exit() 1073 1074class Spin(BdGPatch, MatrixTree): 1075 """ 1076 A class with a name: ↑ and ↓ spin, belonging 1077 as a child to an Orbital. 1078 1079 Attributes 1080 ---------- 1081 None 1082 1083 Methods 1084 ------- 1085 None 1086 1087 Examples 1088 -------- 1089 >>> spin = Spin('up') 1090 >>> print(spin) 1091 Spin(↑) 1092 1093 >>> spin = Spin('↓') 1094 >>> print(spin) 1095 Spin(↓) 1096 1097 >>> spin = Spin(1) 1098 >>> print(spin) 1099 Spin(↓) 1100 1101 Note 1102 ---- 1103 This class acts as a template for spin indices. 1104 """ 1105 def __init__(self, name='↑'): 1106 """Inherit list class method and add parent attribute to items""" 1107 super().__init__() 1108 if name in ['up', '↑', 0]: 1109 self.__name__ = '↑' 1110 elif name in ['dn', 'down', '↓', 1]: 1111 self.__name__ = '↓' 1112 else: 1113 raise ValueError('Varname must be up or dn!') 1114 self.__dict__["parent"] = None 1115 self.__dict__["_idx"] = None 1116 1117 def __str__(self): 1118 """Returns object class name""" 1119 return f'{self.__class__.__name__}({self.__name__})' 1120 1121 def __repr__(self): 1122 """Returns object class name""" 1123 return f'{self.__class__.__name__}({self.__name__})' 1124 1125class Orbital(BdGPatch, MatrixTree, object): 1126 """ 1127 Orbital representation 1128 1129 Args 1130 ---- 1131 n (int): principal quantum number, determines the energy level and size of the orbital 1132 l (int): azimuthal quantum number, defines the shape of the orbital 1133 m (int): magnetic quantum number, defines the orientation of the orbital 1134 1135 Attributes 1136 ---------- 1137 1138 Methods 1139 ------- 1140 radial_function(self, r) 1141 angular_function(self, θ, φ) 1142 compute_wavefunction_cross_section(self, horizontal_interval, vertical_interval) 1143 compute_probability_density(self) 1144 plot_probability_density_cross_section(self, a0, dark_theme=False, colormap='rocket') 1145 TODO add 3D density plot 1146 1147 Examples 1148 -------- 1149 >>> orbital = Orbital(3,2,0) 1150 >>> print(orbital) 1151 Orbital(3d_{z^2}) 1152 1153 Note 1154 ---- 1155 Integrates over the unspecified axis 1156 1157 Reference 1158 --------- 1159 1. https://ssebastianmag.medium.com/computational-physics-with-python-hydrogen-wavefunctions-electron-density-plots-8fede44b7b12 1160 2. https://medium.com/@bldevries/a-symphony-of-spheres-animating-spherical-harmonics-in-blender-with-python-and-shape-keys-aa67b7ff3d93 1161 1162 """ 1163 def __init__(self, n:int=0, l:int=0, m:int=0): 1164 """Inherit list class method and add parent attribute to items""" 1165 1166 # Quantum numbers validation 1167 if not isinstance(n, int) or n < 1: 1168 raise ValueError('n should be an integer satisfying the condition: n >= 1') 1169 if not isinstance(l, int) or not (0 <= l < n): 1170 raise ValueError('l should be an integer satisfying the condition: 0 <= l < n') 1171 if not isinstance(m, int) or not (-l <= m <= l): 1172 raise ValueError('m should be an integer satisfying the condition: -l <= m <= l') 1173 1174 self.__dict__["quantum_numbers"] = [n,l,m] 1175 self.__dict__["n"] = n 1176 self.__dict__["l"] = l 1177 self.__dict__["m"] = m 1178 1179 super().__init__() 1180 self.__dict__["parent"] = None 1181 self.__dict__["_idx"] = None 1182 1183 def __str__(self): 1184 """Returns object class name""" 1185 return f'{self.__class__.__name__}({self.__name__})' 1186 1187 def __repr__(self): 1188 """Returns object class name""" 1189 return f'{self.__class__.__name__}({self.__name__})' 1190 1191 @property 1192 def __name__(self): 1193 n, l, m = self.n, self.l, self.m 1194 name = n 1195 if l==0: 1196 if m==0: 1197 name = 's' 1198 elif l==1: 1199 if m==0: 1200 name ='p_z' 1201 elif m==1: 1202 name ='p_x' 1203 elif m==-1: 1204 name ='p_y' 1205 elif l==2: 1206 if m==0: 1207 name ='d_{z^2}' 1208 elif m==1: 1209 name ='d_{xz}' 1210 elif m==-1: 1211 name ='d_{yz}' 1212 elif m==2: 1213 name ='d_{xy}' 1214 elif m==-2: 1215 name ='d_{x^2-z^2}' 1216 elif l==3: 1217 if m==0: 1218 name ='f_{z^3}' 1219 elif m==1: 1220 name ='f_{xz^2}' 1221 elif m==-1: 1222 name ='f_{yz^2}' 1223 elif m==2: 1224 name ='f_{xyz}' 1225 elif m==-2: 1226 name ='f_{z(x^2-z^2)}' 1227 elif m==3: 1228 name ='f_{x(x^2-3y^2)}' 1229 elif m==-3: 1230 name ='f_{y(3x^2-y^2)}' 1231 if l<=3: 1232 name=str(n)+name 1233 else: 1234 name = f'({n} {l} {m})' 1235 return name 1236 1237 def radial_function(self, r): 1238 """Compute the normalized radial part of the wavefunction using 1239 Laguerre polynomials and an exponential decay factor. 1240 Args: 1241 r (numpy.ndarray): radial coordinate 1242 Returns: 1243 numpy.ndarray: wavefunction radial component 1244 """ 1245 n = self.n 1246 l = self.l 1247 laguerre = sp.genlaguerre(n-l-1, 2*l+1) 1248 p = 2*r/n 1249 1250 constant_factor = np.sqrt( 1251 ((2/n)**3 * (sp.factorial(n-l-1))) / (2*n*(sp.factorial(n+l))) 1252 ) 1253 return constant_factor * np.exp(-p/2) * (p**l) * laguerre(p) 1254 1255 def compute_radial_density(self, rr): 1256 """ 1257 Compute the radial probability density for the 1258 list of radial points rr. 1259 """ 1260 area = 4*np.pi*np.square(rr) 1261 radius = np.square(self.radial_function(rr)) 1262 self.radial_density = radius * area 1263 return self.radial_density 1264 1265 def plot_radial_density(self, ax, rr): 1266 """ 1267 Compute the radial probability density for the 1268 list of radial points rr. 1269 """ 1270 if not hasattr(self, 'radial_density'): 1271 print('Ψ not computed. Computing...') 1272 self.compute_radial_density(rr) 1273 print('Ψ computed.') 1274 1275 ax.plot(rr, self.radial_density, label=self.__name__) 1276 return ax 1277 1278 def angular_function(self, θ, φ): 1279 """ Compute the normalized angular part of the wavefunction using 1280 Legendre polynomials and a phase-shifting exponential factor. 1281 Args: 1282 θ (numpy.ndarray): polar angle 1283 φ (int): azimuthal angle 1284 Returns: 1285 numpy.ndarray: wavefunction angular component 1286 """ 1287 1288 l = self.l 1289 m = self.m 1290 1291 legendre = sp.lpmv(m, l, np.cos(θ)) 1292 1293 constant_factor = ((-1)**m) * np.sqrt( 1294 ((2*l+1) * sp.factorial(l-np.abs(m))) / 1295 (4*np.pi*sp.factorial(l+np.abs(m))) 1296 ) 1297 return constant_factor * legendre * np.real(np.exp(1.j*m*φ)) 1298 1299 def compute_wavefunction(self, box): 1300 """ Compute the normalized wavefunction as a product 1301 of its radial and angular components. 1302 Args: 1303 interval (int, int, int): (start, end, number of points) 1304 Returns: 1305 numpy.ndarray: wavefunction 1306 """ 1307 1308 if not isinstance(box, BoxOfPoints): 1309 raise ValueError('Arguments should be type BoxOfPoints') 1310 1311 # x-y grid to represent electron spatial distribution 1312 x, y, z = box.meshgrid 1313 1314 # Use epsilon to avoid division by zero during angle calculations 1315 eps = np.finfo(float).eps 1316 1317 # Ψnlm(r,θ,φ) = Rnl(r).Ylm(θ,φ) 1318 r = np.sqrt((x**2 + y**2 + z**2)) 1319 θ = np.arccos(z / (r + eps)) 1320 φ = np.sign(y) * np.arccos(x / np.sqrt(x**2 + y**2 + eps)) 1321 self.Ψ = self.radial_function(r) * self.angular_function(θ,φ) 1322 return self.Ψ 1323 1324 def compute_probability_density(self): 1325 """ Compute the probability density of a given wavefunction. 1326 Args: 1327 psi (numpy.ndarray): wavefunction 1328 Returns: 1329 numpy.ndarray: wavefunction probability density 1330 """ 1331 self.density = np.abs(self.Ψ)**2 1332 return self.density 1333 1334 def plot_probability_density_cross_section(self, box:BoxOfPoints, crop=None, colormap='rocket'): 1335 """ Plot the probability density of the hydrogen 1336 atom's wavefunction for a given quantum state (n,l,m). 1337 Args: 1338 a0 (float): scale factor 1339 dark_theme (bool): If True, uses a dark background for the plot, defaults to False 1340 colormap (str): Seaborn plot colormap, defaults to 'rocket' 1341 """ 1342 1343 _labels = ['x','y','z'] 1344 labels = [] 1345 axes = [] 1346 for label, axis in zip(_labels, box.axes): 1347 if axis.number_of_points!=1: 1348 labels.append(label) 1349 axes.append(axis) 1350 1351 if len(axes)>2: 1352 raise ValueError("ERROR: BoxOfPoints must have one axis with only 1 point") 1353 elif len(axes)<2: 1354 raise ValueError("ERROR: BoxOfPoints has insufficiently many axes") 1355 1356 xlabel = '$' + labels[0] + '/' + axes[0].unit + '$' 1357 ylabel = '$' + labels[1] + '/' + axes[1].unit + '$' 1358 1359 n = self.n 1360 l = self.l 1361 m = self.m 1362 1363 # Colormap validation 1364 try: 1365 sns.color_palette(colormap) 1366 except ValueError: 1367 raise ValueError(f'{colormap} is not a recognized Seaborn colormap.') 1368 1369 # Configure plot aesthetics using matplotlib rcParams settings 1370 plt.rcParams['font.family'] = 'STIXGeneral' 1371 plt.rcParams['mathtext.fontset'] = 'stix' 1372 plt.rcParams['xtick.major.width'] = 4 1373 plt.rcParams['ytick.major.width'] = 4 1374 plt.rcParams['xtick.major.size'] = 15 1375 plt.rcParams['ytick.major.size'] = 15 1376 plt.rcParams['xtick.labelsize'] = 30 1377 plt.rcParams['ytick.labelsize'] = 30 1378 plt.rcParams['axes.linewidth'] = 4 1379 1380 fig, ax = plt.subplots(figsize=(16, 16.5)) 1381 plt.subplots_adjust(top=0.82) 1382 plt.subplots_adjust(right=0.905) 1383 plt.subplots_adjust(left=-0.1) 1384 1385 # Compute and visualize the wavefunction probability density 1386 if not hasattr(self, 'Ψ'): 1387 print('Ψ not computed. Computing...') 1388 self.compute_wavefunction(box) 1389 self.compute_probability_density() 1390 print('Ψ computed.') 1391 1392 prob_density = self.density.squeeze() 1393 1394 [x0,x1,y0,y1]=[axes[0].lower_bound, axes[0].upper_bound,axes[1].lower_bound, axes[1].upper_bound] 1395 1396 if crop!=None: 1397 idx,jdx = axes[0].get_cropped_index(*crop[0]) 1398 prob_density = prob_density[idx:jdx,:] 1399 1400 idy,jdy = axes[1].get_cropped_index(*crop[1]) 1401 prob_density = prob_density[:,idy:jdy] 1402 1403 xx = axes[0].get_cropped_values(*crop[0]) 1404 x0, x1 = np.min(xx), np.max(xx) 1405 yy = axes[1].get_cropped_values(*crop[1]) 1406 y0, y1 = np.min(yy), np.max(yy) 1407 1408 vmax = np.max(prob_density) 1409 # Here we transpose the array to align the calculated z-x plane with Matplotlib's y-x imshow display 1410 im = ax.imshow(prob_density, cmap=sns.color_palette(colormap, as_cmap=True), 1411 extent=[x0,x1,y0,y1], vmin=0, vmax=vmax) 1412 1413 cbar = plt.colorbar(im, fraction=0.046, pad=0.03) 1414 cbar.set_ticks([]) 1415 1416 plt.rcParams['text.color'] = '#000000' 1417 title_color = '#000000' 1418 ax.tick_params(axis='x', colors='#000000') 1419 ax.tick_params(axis='y', colors='#000000') 1420 ax.set_xlabel(xlabel, fontsize=36) 1421 ax.set_ylabel(ylabel, fontsize=36) 1422 1423 ax.set_title('Hydrogen Atom - Orbital Wavefunction', 1424 pad=130, fontsize=44, loc='left', color=title_color) 1425 ax.text(x0, 1.12*y0, ( 1426 r'$|ψ_{n \ell m}(r, θ, φ)|^{2} =' 1427 r' |R_{n\ell}(r) Y_{\ell}^{m}(θ, ψ)|^2$' 1428 ), fontsize=36) 1429 ax.text(0.95*x0, 0.9*y0, r'$({0}, {1}, {2})$'.format(n, l, m), color='#dfdfdf', fontsize=42) 1430 ax.text(1.24*x1, 0.322*y1, 'Probability density', rotation='vertical', fontsize=40) 1431 ax.text(1.08*x1, 1.08*y0, f'{vmax:.3}', fontsize=24) 1432 ax.text(1.12*x1, 1.09*y1, 0, fontsize=24) 1433 ax.text(1.24*x1, 0.82*y1, '−', fontsize=34) 1434 ax.text(1.24*x1, 0.78*y0, '+', fontsize=34, rotation='vertical') 1435 ax.invert_yaxis() 1436 1437 # Save and display the plot 1438 plt.savefig(f'({n},{l},{m})-{labels[0]}{labels[1]}.png', bbox_inches='tight') 1439 1440class Atom(BdGPatch, MatrixTree): 1441 """ 1442 Atomic representation in the real-space basis 1443 1444 Attributes 1445 ---------- 1446 1447 Methods 1448 ------- 1449 1450 Examples 1451 -------- 1452 >>> hydrogen = Atom([0.2,0.4], [Orbital(1,0,0)]) 1453 >>> print(hydrogen) 1454 Atom(hydrogen) 1455 1456 Note 1457 ---- 1458 """ 1459 def __init__(self, fractional_coordinates:list=[0,0,0], orbitals=[]): 1460 """Inherit list class method and add parent attribute to items""" 1461 self.__dict__["fractional_coordinates"] = fractional_coordinates 1462 self._name = varname() 1463 super().__init__(*orbitals) 1464 self.__dict__["parent"] = None 1465 self.__dict__["_idx"] = None 1466 for item in orbitals: 1467 item.__dict__["parent"] = self 1468 1469 @property 1470 def __name__(self): 1471 return self._name 1472 1473 def __str__(self): 1474 """Returns object class name""" 1475 return f'{self.__class__.__name__}({self.__name__})' 1476 1477 def __repr__(self): 1478 """Returns object class name""" 1479 return f'{self.__class__.__name__}({self.__name__})' 1480 1481 @property 1482 def orbitals(self): 1483 return _Mapper(self) 1484 1485class Molecule(BdGPatch, MatrixTree): 1486 """ 1487 Molecule in the real-space basis 1488 1489 Attributes 1490 ---------- 1491 1492 Methods 1493 ------- 1494 1495 Examples 1496 -------- 1497 >>> hydrogen_A = Atom(fractional_coordinates=[0.25, 0.33], orbitals=[Orbital(1, 0, 0)]) 1498 >>> hydrogen_B = Atom(fractional_coordinates=[0.75, 0.33], orbitals=[Orbital(1, 0, 0)]) 1499 >>> oxygen = Atom(fractional_coordinates=[0.5,0.66], orbitals=[Orbital(2, 1, 0), Orbital(2, 1, 1), Orbital(2, 1, -1)]) 1500 >>> water = Molecule(fractional_coordinates=[0.5,0.5], atoms=[hydrogen_A, hydrogen_B, oxygen]) 1501 >>> water.n_spins = 2 1502 >>> water['oxygen',:,:] 1503 [Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑)] 1504 1505 >>> oxygen['2p_z'].n_spins = 1 1506 >>> water['oxygen',:,:] 1507 [Spin(↓), Spin(↑), Spin(↓), Spin(↑)] 1508 1509 >>> water[:,:,:] 1510 [Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑)] 1511 1512 Note 1513 ---- 1514 """ 1515 def __init__(self, fractional_coordinates:list=[0,0,0], atoms=[]): 1516 """Inherit list class method and add parent attribute to items""" 1517 self.fractional_coordinates = fractional_coordinates 1518 self._name = varname() 1519 # name = f'{fractional_coordinates}' 1520 super().__init__(*atoms) 1521 self.__dict__["parent"] = None 1522 self.__dict__["_idx"] = None 1523 for item in atoms: 1524 item.__dict__["parent"] = self 1525 1526 @property 1527 def __name__(self): 1528 return self._name 1529 1530 def __str__(self): 1531 """Returns object class name""" 1532 return f'{self.__class__.__name__}({self.__name__})' 1533 1534 def __repr__(self): 1535 """Returns object class name""" 1536 return f'{self.__class__.__name__}({self.__name__})' 1537 1538 @property 1539 def atoms(self): 1540 return _Mapper(self) 1541 1542# class Cluster(Molecule): 1543# """ 1544# Cluster in the real-space basis 1545 1546# Examples 1547# -------- 1548 1549# Note 1550# ---- 1551# """ 1552# def __init__(self, fractional_coordinates:list=[0,0,0], cell:list=[]): 1553# """Inherit list class method and add parent attribute to items""" 1554# super().__init__(fractional_coordinates, cell) 1555 1556class Cell(BdGPatch, MatrixTree): 1557 """ 1558 Cell in the real-space basis 1559 1560 Examples 1561 -------- 1562 1563 >>> A = Atom(fractional_coordinates=[0, 0]) 1564 >>> B = Atom(fractional_coordinates=[0, 0.5]) 1565 >>> cell = Cell(fractional_coordinates=[0,0], sites=[A, B]) 1566 >>> cell 1567 Cell([0, 0]) 1568 1569 """ 1570 def __init__(self, fractional_coordinates, sites, basis = None): 1571 """Inherit list class method and add parent attribute to items""" 1572 self.fractional_coordinates = list(fractional_coordinates) 1573 if type(basis)!=type(None): 1574 self.__dict__["cartesian_coordinates"] = list(np.matmul(fractional_coordinates, basis)) 1575 super().__init__(*sites) 1576 for site in sites: 1577 site.__dict__["parent"] = self 1578 self.__dict__[site.__name__] = site 1579 1580 @property 1581 def __name__(self): 1582 return str(self.fractional_coordinates) 1583 1584 def __repr__(self): 1585 """Returns object class name""" 1586 return f'{self.__class__.__name__}({self.__name__})' 1587 1588 def __str__(self): 1589 """Returns object class name""" 1590 return f'{self.__class__.__name__}({self.__name__})' 1591 1592class Lattice(BdGPatch, MatrixTree): 1593 """ 1594 Lattice 1595 1596 Examples 1597 -------- 1598 >>> A = Atom(fractional_coordinates=[0, 0]) 1599 >>> B = Atom(fractional_coordinates=[0.5, 0.5]) 1600 >>> my_lattice = Lattice(basis=[[1,0.5],[0,1]], cell=[A, B]) 1601 >>> sublattice_half = my_lattice.cut(number_of_cells=3, direction=0) 1602 >>> sublattice_full = my_lattice.cut(number_of_cells=2, direction=1) 1603 1604 Add lattice parameter 1605 >>> my_lattice.t = 1 1606 >>> my_lattice[0,0].t 1607 1 1608 1609 Add lattice parameter to particular cell 1610 >>> my_lattice[1,0].t = 2 1611 >>> my_lattice[0].t 1612 [1, 1] 1613 >>> my_lattice.t 1614 1 1615 >>> my_lattice[1,0].A 1616 Atom(A) 1617 >>> my_lattice[0].A 1618 [Atom(A), Atom(A)] 1619 >>> my_lattice.A 1620 Atom(A) 1621 1622 # Cutout a sub-region 1623 # >>> square = my_lattice.cutout(boundary=[[-1, -1], [-1, 1], [1, 1], [1, -1]]) 1624 # >>> my_lattice.boundaries["square"] 1625 # {'boundary': [[-1, -1], [-1, 1], [1, 1], [1, -1]], 'cells': [Cell([0, 0]), Cell([0, 1]), Cell([1, 0]), Cell([1, 1])]} 1626 1627 Note 1628 ---- 1629 """ 1630 def __init__(self, basis:list=[[1,0,0],[0,1,0],[0,0,1]], cell=[]): 1631 """Inherit list class method and add parent attribute to items""" 1632 self.__dict__["_name"] = varname() 1633 self.__dict__["basis"] = np.array(basis) 1634 self.__dict__["inv_basis"] = np.linalg.inv(basis) 1635 self.__dict__["cell"] = cell 1636 self.__dict__["nd"] = len(self.basis) 1637 self.__dict__["n_cells"] = [1 for _ in range(self.nd)] 1638 self.__dict__["pbc"] = [True for _ in range(self.nd)] 1639 self.__dict__["boundaries"] = {} 1640 self.__dict__["subregions"] = {} 1641 self.__dict__["n_spins"] = 1 1642 self._initialise_cells() 1643 super().__init__() 1644 1645 def _initialise_cells(self): 1646 """ 1647 Initialises self.cells as an empty array --called by self.__init__() 1648 """ 1649 coord = [0 for _ in range(self.nd)] 1650 sites = deepcopy(self.cell) 1651 cell = Cell(fractional_coordinates = coord, sites = sites, basis = self.basis) 1652 for parent_site in self.cell: 1653 parent_site.clear() 1654 for child_site in cell: 1655 name = parent_site.__name__ 1656 if name == child_site.__name__: 1657 self.__dict__[name] = parent_site 1658 parent_site.extend(child_site) 1659 child_site.__dict__["parent"] = parent_site 1660 array = np.empty(self.n_cells, dtype = Cell) 1661 array[*coord] = cell 1662 self.__dict__["cells"] = array 1663 1664 # def __getattr__(self, name: str): 1665 # try: 1666 # return self.__dict__[name] 1667 # except: 1668 # return _Mapper([getattr(item, name) for item in _Mapper(self.cells.tolist())]) 1669 1670 # def __setattr__(self, name, value, update = True): 1671 # """ 1672 # Set attributes down the lineage, and add attribute to list of 1673 # attributes to be updates 1674 # """ 1675 # if update: 1676 # self.__cache__(name, value) 1677 1678 # self.__dict__[f"{name}"] = value 1679 1680 # [item.__setattr__(name, value, update = False) for item in self] 1681 1682 # def __setattr__(self, name: str, value): 1683 # if name in self.__dict__.keys(): 1684 # self.__dict__[name] = value 1685 # else: 1686 # [setattr(item, name, value) for item in _Mapper(self.cells.tolist())] 1687 1688 @property 1689 def __name__(self): 1690 return self._name 1691 1692 def __str__(self): 1693 """Returns object class name""" 1694 return f'{self.__name__}({[list(x) for x in self.basis]})' 1695 1696 def __repr__(self): 1697 """Returns object class name""" 1698 return f'{self.__class__.__name__}({[list(x) for x in self.basis]})' 1699 1700 def cut(self, number_of_cells: int, direction: int, glue_edges: bool = True): 1701 """ 1702 Returns a deepcopy of self 1703 """ 1704 1705 if not glue_edges: 1706 raise ValueError("TODO: Open boundary conditions not yet implemented!") 1707 1708 _slice = [(slice(0, self.n_cells[i])) for i in range(self.nd)] 1709 self.n_cells[direction] = number_of_cells 1710 cells = np.empty(self.n_cells, dtype=Cell) 1711 cells[*_slice] = self.cells 1712 self.cells = cells 1713 indices = np.moveaxis(np.indices(self.n_cells), 0, -1) 1714 indices = indices[cells == None] 1715 for index in indices: 1716 cell = Cell(fractional_coordinates=index, sites=deepcopy(self.cell), basis=self.basis) 1717 self.cells[*index] = cell 1718 1719 for parent_site in self.cell: 1720 parent_site.clear() 1721 for child_site in cell: 1722 name = parent_site.__name__ 1723 if name == child_site.__name__: 1724 self.__dict__[name] = parent_site 1725 parent_site.extend(child_site) 1726 child_site.__dict__["parent"] = parent_site 1727 cell = deepcopy(self.cell) 1728 cells = [Cell(fractional_coordinates=index, sites=[deepcopy(site) for site in self.cell], basis=self.basis) for index in indices] 1729 indices = np.array(indices).T 1730 self.cells[*indices] = cells 1731 1732 def mapper(parent_site, cell): 1733 parent_site.clear() 1734 for child_site in cell: 1735 name = parent_site.__name__ 1736 if name == child_site.__name__: 1737 self.__dict__[name] = parent_site 1738 parent_site.extend(child_site) 1739 child_site.__dict__["parent"] = parent_site 1740 1741 [[mapper(site, cell) for site in self.cell] for cell in cells] 1742 1743 self.pbc[direction] = glue_edges 1744 1745 self.init_cells() 1746 1747 try: 1748 name = varname() 1749 sublattice = Lattice(basis=copy(self.basis), cell=copy(self.cell)) 1750 sublattice._name = name 1751 sublattice.n_cells = copy(self.n_cells) 1752 sublattice.pbc = copy(self.pbc) 1753 sublattice.cells = deepcopy(self.cells) 1754 sublattice.extend(flatten(sublattice.cells)) 1755 return sublattice 1756 except: 1757 pass 1758 1759 @property 1760 def n_total_cells(self): 1761 return np.prod(self.n_cells) 1762 1763 def cutout(self, boundary, cut_from_infinite_plane = True): 1764 """ 1765 Cutout cells within the boundary. 1766 1767 The boundary is a list of points in the 1768 Cartesian. 1769 """ 1770 try: 1771 name = varname() 1772 except: 1773 name = len(self.boundaries) 1774 1775 if cut_from_infinite_plane: 1776 inv_boundary = np.matmul(boundary, self.inv_basis) 1777 polytope = Polytope(points=inv_boundary) 1778 inside, outside = polytope.cut() 1779 corners = polytope.frac_corners 1780 n_cells = corners[:,1] - corners[:,0] 1781 for d, n in enumerate(n_cells): 1782 self.cut(number_of_cells=int(n+1), direction=d) 1783 1784 polytope = Polytope(points=boundary) 1785 1786 # clip 1787 # ax = polytope.plot() 1788 # cartesian_coordinates = flatten(self.cartesian_coordinates, self.nd-1) 1789 # cartesian_coordinates = clean(cartesian_coordinates) 1790 # cartesian_coordinates = np.transpose(cartesian_coordinates) 1791 # # plt.scatter(*cartesian_coordinates) 1792 1793 # inside = np.matmul(inside, self.basis) 1794 # outside = np.matmul(outside, self.basis) 1795 # ax.scatter(*inside.T) 1796 # ax.scatter(*outside.T, c="r") 1797 # lim=[-5,20] 1798 # ax.set(xlim=lim, ylim=lim) 1799 # plt.show() 1800 # exit() 1801 1802 inside = inside.astype(int) 1803 outside = outside.astype(int) 1804 1805 cells = [cell.cartesian_coordinates for cell in self] 1806 cells = np.array(cells, dtype=float) 1807 1808 cartesian = np.matmul(inside, self.basis) 1809 for cartesian, cell in zip(cartesian, inside): 1810 self.cells[*cell].__dict__["fractional_coordinates"] = cell 1811 self.cells[*cell].__dict__["cartesian_coordinates"] = cartesian 1812 cells = [self.cells[*cell] for cell in inside] 1813 else: 1814 polytope = Polytope(points=boundary, basis=self.basis) 1815 cells = flatten_depth(self.cartesian_coordinates, self.nd-1) 1816 cells = [list(cell) for cell in cells] 1817 cells = np.array(cells, dtype=float) 1818 inside, outside = polytope.cut(cells) 1819 inside = np.matmul(inside, self.inv_basis).astype(int) 1820 outside = np.matmul(outside, self.inv_basis).astype(int) 1821 cells = [self.cells[*cell] for cell in inside] 1822 1823 for cell in outside: 1824 try: 1825 self.cells[*cell] = None 1826 except: 1827 pass 1828 1829 self.init_cells() 1830 1831 self.boundaries[name] = polytope 1832 self.subregions[name] = cells 1833 1834 return self.boundaries[name] 1835 1836 def init_cells(self): 1837 self.clear() 1838 self.extend(flatten(self.cells)) 1839 1840 def __set_indices__(self): 1841 """ 1842 """ 1843 self.initialised = True 1844 def f(l, i): l.__dict__["index"] = i 1845 [f(leaf, i) for i, leaf in enumerate(self.leaves)] 1846 for site in self.cell: 1847 _slice = [slice(None, None, None) for _ in range(self.nd+site.depth)] 1848 _slice.append(site.__name__) 1849 site.clear() 1850 site.extend(flatten(self[*_slice])) 1851 1852 def plt_boundary(self, ax=None): 1853 for boundary in self.boundaries: 1854 ax = boundary.plot(ax) 1855 return ax 1856 1857 def plt_lattice(self, ax=None): 1858 ax.scatter(*self.cartesian().T, c='b') 1859 return ax 1860 1861 def __getitem__(self,keys): 1862 cells = self.cells 1863 _keys = [] 1864 if type(keys)==int: 1865 return _Mapper(self.cells[keys]) 1866 if len(keys)>self.nd: 1867 _keys = keys[self.nd:] 1868 keys = keys[:self.nd] 1869 1870 cells = cells[keys] 1871 if type(cells)==np.ndarray: 1872 cells = _Mapper(cells.tolist()) 1873 1874 _slice = [] 1875 for key in keys: 1876 if key==slice(None,None,None): 1877 _slice.append(slice(None,None,None)) 1878 1879 if len(_keys)>0: 1880 cells = _Mapper(_Mapper(cells)[*_slice,*_keys]) 1881 1882 return cells 1883 1884 def plot(self, ax=None): 1885 if ax==None: 1886 ax = plt.axes() 1887 plt.tight_layout() 1888 cartesian_coordinates = flatten_depth(self.cartesian_coordinates, self.nd-1) 1889 cartesian_coordinates = clean(cartesian_coordinates) 1890 cartesian_coordinates = np.transpose(cartesian_coordinates) 1891 plt.scatter(*cartesian_coordinates) 1892 1893def flatten_depth(items, depth:int): 1894 if depth>0: 1895 return [flatten_depth(item, depth-1) for sublist in items for item in sublist] 1896 elif depth==0: 1897 return items 1898 else: 1899 raise ValueError("Depth must be greater than 0!") 1900 1901def flatten(l, ltypes=(list, tuple, _Mapper)): 1902 if type(l)==np.ndarray: 1903 l=l.tolist() 1904 ltype = type(l) 1905 l = list(l) 1906 i = 0 1907 while i < len(l): 1908 while type(l[i]) in ltypes: 1909 if not l[i]: 1910 l.pop(i) 1911 i -= 1 1912 break 1913 else: 1914 l[i:i + 1] = l[i] 1915 i += 1 1916 return ltype(l) 1917 1918def clean(items): 1919 return [list(item) for item in items if type(item)!=type(None)] 1920 1921def friedel_oscillations(): 1922 torch.set_default_dtype(torch.float32) 1923 1924 A = Atom(fractional_coordinates=[0, 0]) 1925 B = Atom(fractional_coordinates=[0.5, 0.5]) 1926 # A.hop(B, t=1) 1927 # A.hop(B, vector=[1,0], t=1) 1928 # B.hop(B, vector=[1,0], t=1) 1929 # A.hop(A, vector=[0,1], t=1) 1930 # B.hop(B, vector=[0,1], t=1) 1931 # A.add(μ=-3.95) 1932 lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A]) 1933 # lattice = Lattice(basis=[[1, 0],[0, 1]]) 1934 # square = lattice.cutout(boundary=[[-1, -1], [-1, 1], [1, 1], [1, -1]]) 1935 lattice.hop(vector=[1,0], t=1) 1936 lattice.hop(vector=[0,1], t=1) 1937 sublattice_half = lattice.cut(number_of_cells=11, direction=0) 1938 sublattice_half = lattice.cut(number_of_cells=11, direction=1) 1939 # lattice.n_spins = 2 1940 lattice.add(μ=3.95, structure=-1) 1941 # lattice.add(μ=-3.95, structure=-1, pauli_vector=[0,1,1,0]) 1942 # lattice.add(μ=-3.95, structure=1, pauli_vector=[0,0,0,1]) 1943 # lattice[1,0].add(V=-1., structure=-1, pauli_vector=[0,0,0,1]) 1944 lattice[5,5].add(V=1., structure=-1) 1945 # lattice[1,0].hop(lattice[0,0], t=1) 1946 # lattice.hop(vector=[1,0], t=1, pauli_vector=[0,0,0,1]) 1947 # lattice.hop(vector=[1,0], t=1, pauli_vector=[1,0,0,0]) 1948 # exit() 1949 # lattice.add(Δ=5j, anomalous=True, pauli_vector=[0,0,1,0]) 1950 # sublattice_half = my_lattice.cut(number_of_cells=3, direction=0) 1951 # sublattice_full = my_lattice.cut(number_of_cells=3, direction=1) 1952 # print(my_lattice.cells) 1953 # exit() 1954 # square = my_lattice.cutout(boundary=[[-10, -10], [-10, 10], [10, 10], [10, -10]]) 1955 # print(my_lattice[:,:,:,:]) 1956 1957 # print(my_lattice.indices) 1958 # exit() 1959 1960 # matrix = lattice.get_matrix() 1961 1962 1963 # exit() 1964 # matrix = np.real(matrix) 1965 # exit() 1966 # c = plt.imshow(matrix) 1967 # plt.colorbar(c) 1968 # plt.show() 1969 # exit() 1970 1971 w, v = lattice.solve() 1972 1973 # exit() 1974 # print(w.dtype) 1975 # print(v) 1976 # print(w) 1977 # print(np.shape(v)) 1978 # exit() 1979 # dm = my_lattice.density_matrix() 1980 # print(np.shape(dm)) 1981 # fermi = my_lattice.fermi_distribution(temperature=10) 1982 # fermi = my_lattice.thermal_density_matrix(temperature=0) 1983 # print(t1) 1984 # print(t2) 1985 # greens = my_lattice.greens_function(energy=0, resolution=0.1, sites=[A, B], sites_to=[B, A], hop=[1,0], spin_flip=True, anomalous = True) 1986 # greens = my_lattice.greens_function(energy=arange(-4,4,0.1), resolution=0.1, sites=[A, B], sites_to=[B, A], hop=[1,0], spin_flip=True, anomalous = True) 1987 # greens = my_lattice.greens_function(energy=0, resolution=0.1, sites=[A, B], sites_to=[B, A], hop=[1,0], spin_flip=True, anomalous = True, spin_polarised=True) 1988 greens = lattice.greens_function(energy=0, resolution=0.1) 1989 # energy = arange(-4,4,0.1) 1990 # greens = my_lattice.greens_function(energy=energy, resolution=0.1) 1991 # greens = sum(greens, axis=0) 1992 # plt.scatter(energy, greens) 1993 # plt.show() 1994 # exit() 1995 # greens = greens.reshape(11,11) 1996 c = plt.imshow(greens) 1997 plt.colorbar(c) 1998 plt.show() 1999 2000 exit() 2001 2002 # square = my_lattice.cutout(boundary=[[0, 0], [0, 10], [10, 10], [10, 0]]) 2003 2004 def plot(): 2005 ax = plt.axes() 2006 plt.tight_layout() 2007 my_lattice.plot(ax) 2008 # square.plot(ax) 2009 lim=[-1,3] 2010 ax.set(xlim=lim, ylim=lim) 2011 plt.show() 2012 2013def mean_field_theory(): 2014 2015 A = Atom(fractional_coordinates=[0, 0]) 2016 lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A]) 2017 lattice.n_spins = 2 2018 # lattice.hop(vector=[1,0], t=1, structure=-1) 2019 # lattice.hop(vector=[0,1], t=1, structure=-1) 2020 lattice.cut(number_of_cells=1, direction=0) 2021 lattice.cut(number_of_cells=1, direction=1) 2022 # lattice.add(μ=-1.7, structure=-1) 2023 # lattice.add(U=2.8, structure=-1, interaction=True) 2024 # lattice.add(φ=2.40, structure=-1, renormalise=True) 2025 lattice.add(Δ=1.00, anomalous=True, pauli_vector=[1,0,0,0], renormalise=True) 2026 2027 # lattice[5,5].add(V=1., structure=-1) 2028 2029 print(lattice._diagonal) 2030 print(lattice._offdiagonal) 2031 # print(lattice._hartree) 2032 # print(lattice._fock) 2033 print(lattice._gorkov) 2034 # print(lattice._interaction) 2035 exit() 2036 2037 print(lattice.get_matrix()) 2038 exit() 2039 2040 w, v = lattice.solve() 2041 exit() 2042 lattice.extract_fields(w, v) 2043 # lattice.self_consistent_step() 2044 2045 exit() 2046 2047 lattice.add(ρ=3.15, structure=-1, renormalise=True) 2048 lattice.add(φ=2.15, structure=-1, renormalise=True) 2049 lattice.add(Δ=3, anomalous=True, pauli_vector=[0,1j,0,0], renormalise=True) 2050 2051 matrix = lattice.matrix 2052 2053 print(matrix) 2054 2055 2056 exit() 2057 2058 # lattice.add(Δ=5j, anomalous=True, pauli_vector=[0,1,0,0]) 2059 2060 lattice.self_consistent_step() 2061 2062 # lattice.scmft(friction=friction,max_iterations=max_iterations,convergence_factor=convergence_factor) 2063 2064 2065 exit() 2066 2067 w, v = lattice.solve() 2068 2069 print(w) 2070 2071# mean_field_theory() 2072 2073def test_indices(): 2074 2075 A = Atom(fractional_coordinates=[0, 0]) 2076 B = Atom(fractional_coordinates=[0.5, 0]) 2077 lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 2078 lattice.cut(number_of_cells=3, direction=0) 2079 lattice.cut(number_of_cells=3, direction=1) 2080 lattice.n_spins = 2 2081 lattice.hop(t=1, vector=[1,0], structure=-1) 2082 # lattice.hop(t=1, vector=[0,1], structure=-1) 2083 # lattice._initialise_cells 2084 2085 lattice.__set_indices__() 2086 2087 # >>> lattice.get_indices(site=A, anomalous=True, structure=-1, pauli_vector=[0,1,0,0]) 2088 # {(1+0j): [37, 45, 53, 61, 69, 40, 48, 56, 64]} 2089 exit() 2090 2091 print(lattice.get_diagonal_indices(site=A, anomalous=False, structure=-1, pauli_vector=[0,1,0,0])) 2092 exit() 2093 2094 print(lattice.get_offdiagonal_indices(site, neighbour, vector, trs, anomalous, structure, pauli_vector)) 2095 exit() 2096 2097 exit() 2098 2099 2100 # lattice.add(U=2.8, structure=-1, interaction=True) 2101 # lattice.add(φ=2.40, structure=-1, renormalise=True) 2102 # lattice.add(Δ=1.00, anomalous=True, pauli_vector=[1,0,0,0], renormalise=True) 2103 2104 # lattice[5,5].add(V=1., structure=-1) 2105 2106# test_indices() 2107 2108# A = Atom(fractional_coordinates=[0, 0]) 2109# B = Atom(fractional_coordinates=[0.5, 0]) 2110# C = Atom(fractional_coordinates=[0.5, 0]) 2111# A.hop(t=1, neighbour=B, vector=[1,0], structure=-1) 2112# A.hop(t=-1, neighbour=B) 2113# B.hop(t=-1, neighbour=C) 2114# C.hop(t=-1, neighbour=A) 2115# lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B, C]) 2116lattice = Lattice(basis=[[1, 0],[0, 1]]) 2117lattice.cut(number_of_cells=1, direction=0) 2118lattice.cut(number_of_cells=1, direction=1) 2119lattice.n_spins = 1 2120# lattice.add(μ=1, pauli_vector=[0,1,0,0]) 2121# lattice[1,0].add(V=-17.0, structure=-1) 2122# lattice.hop(t=1, vector=[1,0], structure=-1) 2123# lattice.hop(t=1, vector=[0,1], structure=-1) 2124 2125lattice.add(μ=0.1, structure=-1) 2126lattice.add(ρ=0.00, renormalise=True) 2127# lattice.add(Φ=1.00, renormalise=True) 2128lattice.add(Δ=2.28, anomalous=True, renormalise=True) 2129 2130lattice.add(U=2.8, structure=-1, interaction=True) 2131 2132w, v = lattice.solve() 2133lattice.extract_mean_fields() 2134 2135w, v = lattice.solve() 2136lattice.extract_mean_fields() 2137 2138w, v = lattice.solve() 2139lattice.extract_mean_fields() 2140 2141w, v = lattice.solve() 2142lattice.extract_mean_fields() 2143 2144w, v = lattice.solve() 2145lattice.extract_mean_fields() 2146 2147w, v = lattice.solve() 2148lattice.extract_mean_fields() 2149 2150w, v = lattice.solve() 2151lattice.extract_mean_fields() 2152 2153w, v = lattice.solve() 2154lattice.extract_mean_fields() 2155 2156w, v = lattice.solve() 2157lattice.extract_mean_fields() 2158 2159w, v = lattice.solve() 2160lattice.extract_mean_fields() 2161 2162w, v = lattice.solve() 2163lattice.extract_mean_fields() 2164 2165w, v = lattice.solve() 2166lattice.extract_mean_fields() 2167 2168w, v = lattice.solve() 2169lattice.extract_mean_fields() 2170 2171w, v = lattice.solve() 2172lattice.extract_mean_fields() 2173 2174w, v = lattice.solve() 2175lattice.extract_mean_fields() 2176# print(lattice.get_matrix())
COMPLEX =
torch.complex64
pauli_z =
tensor([[ 1.+0.j, 0.+0.j],
[ 0.+0.j, -1.+0.j]])
pauli_y =
tensor([[0.+0.j, -0.-1.j],
[0.+1.j, 0.+0.j]])
pauli_x =
tensor([[0.+0.j, 1.+0.j],
[1.+0.j, 0.+0.j]])
pauli_plus =
tensor([[0.+0.j, 1.+0.j],
[0.+0.j, 0.+0.j]])
pauli_minus =
tensor([[0.+0.j, 0.+0.j],
[1.+0.j, 0.+0.j]])
eye =
tensor([[1.+0.j, 0.+0.j],
[0.+0.j, 1.+0.j]])
pauli_vector =
tensor([[[ 1.+0.j, 0.+0.j],
[ 0.+0.j, 1.+0.j]],
[[ 0.+0.j, 1.+0.j],
[ 1.+0.j, 0.+0.j]],
[[ 0.+0.j, -0.-1.j],
[ 0.+1.j, 0.+0.j]],
[[ 1.+0.j, 0.+0.j],
[ 0.+0.j, -1.+0.j]]])
def
split_structure_dict(structure_dict):
54def split_structure_dict(structure_dict): 55 diagonals, offdiagonals = {}, {} 56 for structure, [indices, alt_indices] in structure_dict.items(): 57 if type(indices)==list: 58 indices = tensor(indices) 59 alt_indices = tensor(alt_indices) 60 _bool = eq(indices, alt_indices) 61 _diagonals = [*indices[_bool].tolist()], [*alt_indices[_bool].tolist()] 62 _not_bool = logical_not(_bool) 63 _offdiagonals = [*indices[_not_bool].tolist()], [*alt_indices[_not_bool].tolist()] 64 diagonals[structure] = _diagonals 65 offdiagonals[structure] = _offdiagonals 66 return diagonals, offdiagonals
class
Hamiltonian:
68class Hamiltonian(): 69 70 @property 71 def dof(self): 72 try: 73 return self.__dict__["_dof"] 74 except: 75 self.__dict__["_dof"] = self.n_leaves 76 return self.__dict__["_dof"] 77 78 def solve(self, overwrite_matrix=False): 79 80 if self.is_anomalous(): 81 self.eigenvalues = empty([2*self.dof]) 82 else: 83 self.eigenvalues = empty([self.dof]) 84 85 if overwrite_matrix: 86 (self.eigenvalues,self.eigenvectors) = eigh(self.matrix, out=(self.eigenvalues, self.matrix)) 87 else: 88 (self.eigenvalues,self.eigenvectors) = eigh(self.matrix) 89 90 # if self.is_anomalous(): 91 # self.eigenvalues[:self.dof] = flip(self.eigenvalues[:self.dof], [0]) 92 # self.eigenvectors[:,:self.dof] = flip(self.eigenvectors[:,:self.dof], [1]) 93 94 return self.eigenvalues, self.eigenvectors 95 96 def _check_eigenvectors(self): 97 if not hasattr(self, "eigenvectors"): 98 raise AttributeError("Self has no attribute eigenvectors!")
def
solve(self, overwrite_matrix=False):
78 def solve(self, overwrite_matrix=False): 79 80 if self.is_anomalous(): 81 self.eigenvalues = empty([2*self.dof]) 82 else: 83 self.eigenvalues = empty([self.dof]) 84 85 if overwrite_matrix: 86 (self.eigenvalues,self.eigenvectors) = eigh(self.matrix, out=(self.eigenvalues, self.matrix)) 87 else: 88 (self.eigenvalues,self.eigenvectors) = eigh(self.matrix) 89 90 # if self.is_anomalous(): 91 # self.eigenvalues[:self.dof] = flip(self.eigenvalues[:self.dof], [0]) 92 # self.eigenvectors[:,:self.dof] = flip(self.eigenvectors[:,:self.dof], [1]) 93 94 return self.eigenvalues, self.eigenvectors
100class StatisticalMechanics(Hamiltonian): 101 102 def fermi_distribution(self, temperature:float=0): 103 if temperature == 0: 104 if self.anomalous: 105 distribution = ones(2*self.dof, dtype=COMPLEX) 106 else: 107 distribution = ones(self.dof, dtype=COMPLEX) 108 # Depopulate negative energy solutions 109 if self.anomalous: 110 ## Make use of BdG symmetry 111 distribution[..., :self.dof] = 0 112 else: 113 distribution[self.eigenvalues > 0] = 0 114 115 else: 116 distribution = 1/exp((self.eigenvalues/temperature)+1) 117 118 return distribution 119 120 def _get_lattice_resolved_indices(self, site=None): 121 """ 122 Returns indices of specified site as an array. 123 The first self.nd dimensions of the array are the dimensions of the lattice 124 125 Examples 126 -------- 127 128 >>> A = Atom(fractional_coordinates=[0, 0]) 129 >>> B = Atom(fractional_coordinates=[0.5, 0]) 130 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 131 >>> lattice.cut(number_of_cells=3, direction=0) 132 >>> lattice.cut(number_of_cells=3, direction=1) 133 >>> lattice.n_spins = 2 134 135 >>> lattice.__set_indices__() 136 >>> lattice._get_lattice_resolved_indices(site=None) 137 [[[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]], [[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]]] 138 >>> lattice._get_lattice_resolved_indices(site=lattice) 139 [[[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]], [[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]]] 140 >>> lattice._get_lattice_resolved_indices(site=A) 141 [[[[0, 1]], [[4, 5]], [[8, 9]]], [[[12, 13]], [[16, 17]], [[20, 21]]], [[[24, 25]], [[28, 29]], [[32, 33]]]] 142 >>> lattice._get_lattice_resolved_indices(site=lattice[:,:,:,"↑"]) 143 [[[[[0]], [[2]]], [[[4]], [[6]]], [[[8]], [[10]]]], [[[[12]], [[14]]], [[[16]], [[18]]], [[[20]], [[22]]]], [[[[24]], [[26]]], [[[28]], [[30]]], [[[32]], [[34]]]]] 144 >>> lattice._get_lattice_resolved_indices(site=lattice[:,:,A,"↑"]) 145 [[[[[0]]], [[[4]]], [[[8]]]], [[[[12]]], [[[16]]], [[[20]]]], [[[[24]]], [[[28]]], [[[32]]]]] 146 >>> lattice._get_lattice_resolved_indices(site=[A, B]) 147 [[[[0, 1], [4, 5]], [[8, 9], [12, 13]], [[16, 17], [20, 21]]], [[[24, 25], [28, 29]], [[32, 33], [2, 3]], [[6, 7], [10, 11]]], [[[14, 15], [18, 19]], [[22, 23], [26, 27]], [[30, 31], [34, 35]]]] 148 149 """ 150 _slice = [slice(None, None, None) for _ in range(self.nd)] 151 if site==None: 152 site = self 153 _type = type(site) 154 if _type==Lattice: 155 return site[*_slice].indices 156 elif _type==list: 157 indices = [self._get_lattice_resolved_indices(site=site) for site in site] 158 _shape = list(np.shape(indices)) 159 height = site[0].height 160 _shape[-height] = len(site) 161 _shape = _shape[1:] 162 indices = np.reshape(indices, _shape) 163 return indices.tolist() 164 elif _type==_Mapper: 165 return site.indices 166 elif _type==Cell: 167 return site.indices 168 else: 169 lattice_depth = self.height 170 depth = lattice_depth - site.height 171 for _ in range(depth-1): 172 _slice.append(slice(None, None, None)) 173 _slice.append(site) 174 cells = self[*_slice] 175 return cells.indices 176 177 def _get_rolled_indices(self, indices, hop_vector:None): 178 """ 179 indices should be output of self._get_lattice_resolved_indices(...) 180 181 Examples 182 -------- 183 184 >>> A = Atom(fractional_coordinates=[0, 0]) 185 >>> B = Atom(fractional_coordinates=[0.5, 0]) 186 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 187 >>> lattice.cut(number_of_cells=3, direction=0) 188 >>> lattice.cut(number_of_cells=3, direction=1) 189 >>> lattice.n_spins = 2 190 191 >>> lattice.__set_indices__() 192 193 >>> indices = [[[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]], [[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]]] 194 >>> lattice._get_rolled_indices(indices, [1,0]) 195 [[[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]], [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]] 196 197 """ 198 if len(hop_vector)!=self.nd: 199 raise ValueError("Hop must have same length as model dimensions!") 200 for axis in range(self.nd): 201 shift = hop_vector[axis] 202 if self.pbc[axis]: 203 indices = np.roll(indices, shift=-shift, axis=axis) 204 else: 205 indices = np.add(indices, shift) 206 dimension = self.nd[axis] 207 indices[indices >= dimension] = None 208 return indices.tolist() 209 210 def _get_spin_matrix_indices(self, indices, spin_vector, structure=1): 211 """ 212 indices should be output of self._get_lattice_resolved_indices(...) or self._get_rolled_indices(...) 213 214 Examples 215 -------- 216 217 >>> A = Atom(fractional_coordinates=[0, 0]) 218 >>> B = Atom(fractional_coordinates=[0.5, 0]) 219 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 220 >>> lattice.cut(number_of_cells=3, direction=0) 221 >>> lattice.cut(number_of_cells=3, direction=1) 222 >>> lattice.n_spins = 2 223 224 >>> lattice.__set_indices__() 225 226 >>> indices = [[[12, 13, 14, 15], [16, 17, 18, 19], [20, 21, 22, 23]], [[24, 25, 26, 27], [28, 29, 30, 31], [32, 33, 34, 35]], [[0, 1, 2, 3], [4, 5, 6, 7], [8, 9, 10, 11]]] 227 >>> lattice._get_spin_matrix_indices(indices, structure=-1) 228 {(-1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]} 229 >>> lattice._get_spin_matrix_indices(indices, [0,1,0,0]) 230 {(1+0j): [[12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 0, 2, 4, 6, 8, 10, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 1, 3, 5, 7, 9, 11], [13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 33, 35, 1, 3, 5, 7, 9, 11, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 0, 2, 4, 6, 8, 10]]} 231 232 """ 233 234 indices = list(flatten(indices)) 235 236 if type(spin_vector)==list: 237 spin_vector = tensor(spin_vector, dtype=COMPLEX) 238 _pauli_vector = einsum('i,ijk->jk', spin_vector, pauli_vector) 239 structure_dict = {} 240 for i in range(2): 241 for j in range(2): 242 _structure = complex(_pauli_vector[i,j] * structure) 243 if abs(_structure)!=0: 244 _set_dict(structure_dict, structure, indices[i::2], indices[j::2]) 245 return structure_dict 246 247 def get_diagonal_indices(self, site, structure, pauli_vector=None): 248 """ 249 Examples 250 -------- 251 252 >>> A = Atom(fractional_coordinates=[0, 0]) 253 >>> B = Atom(fractional_coordinates=[0.5, 0]) 254 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 255 >>> lattice.cut(number_of_cells=3, direction=0) 256 >>> lattice.cut(number_of_cells=3, direction=1) 257 >>> lattice.n_spins = 2 258 259 >>> lattice.__set_indices__() 260 261 >>> lattice.get_diagonal_indices(site=A, structure=-1) 262 {(-1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]} 263 264 >>> lattice.get_diagonal_indices(site=A, structure=-1, pauli_vector=[0,0,0,1]) 265 {(1+0j): [[0, 4, 8, 12, 16, 20, 24, 28, 32, 1, 5, 9, 13, 17, 21, 25, 29, 33], [1, 5, 9, 13, 17, 21, 25, 29, 33, 0, 4, 8, 12, 16, 20, 24, 28, 32]]} 266 """ 267 268 indices = self._get_lattice_resolved_indices(site) 269 270 if type(pauli_vector)!=type(None): 271 structure = self._get_spin_matrix_indices(indices, pauli_vector, structure) 272 else: 273 indices = flatten(list(indices)) 274 structure = {complex(structure): [indices, indices]} 275 276 return structure 277 278 def get_offdiagonal_indices(self, site, neighbour, hop_vector=None, trs=True, structure=1, pauli_vector=None): 279 """ 280 Examples 281 -------- 282 283 >>> lattice.get_offdiagonal_indices(site=A, neighbour=B, hop_vector=[1,0], trs=True, structure=1j) 284 {1j: [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11]], -1j: [[14, 15, 18, 19, 22, 285 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]} 286 >>> lattice.get_offdiagonal_indices(site=A, neighbour=B, hop_vector=[1,0], trs=True, structure=1) 287 {(1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11], [14, 15, 18, 19, 22, 23, 288 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11, 0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 1 289 0, 11]]} 290 291 TODO add parameter to flip spin with Pauli vector 292 293 """ 294 if type(pauli_vector)!=type(None): 295 raise ValueError("Pauli vector not yet implemented!") 296 297 indices = self._get_lattice_resolved_indices(site) 298 alt_indices = self._get_lattice_resolved_indices(neighbour) 299 300 if type(hop_vector)!=type(None): 301 alt_indices = self._get_rolled_indices(alt_indices, hop_vector) 302 303 indices = flatten(list(indices)) 304 alt_indices = flatten(list(alt_indices)) 305 structure_dict = {} 306 if trs: 307 _indices = deepcopy(indices) 308 _alt_indices = deepcopy(alt_indices) 309 310 _set_dict(structure_dict, structure, indices, alt_indices) 311 312 if trs: 313 structure = scalar_conj(structure) 314 _set_dict(structure_dict, structure, _alt_indices, _indices) 315 316 return structure_dict 317 318 def _init_sparse_matrices(self): 319 320 self.__set_indices__() 321 322 self.anomalous = False 323 self.interacting = False 324 325 self.hamiltonian_sparse_matrix = {} 326 self.interaction_sparse_matrix = {} 327 self.hartree_sparse_matrix = {} 328 self.hartree_sparse_matrix = {} 329 self.fock_sparse_matrix = {} 330 self.gorkov_sparse_matrix = {} 331 self._sparse_matrices = [self.hamiltonian_sparse_matrix, self.interaction_sparse_matrix, self.hartree_sparse_matrix, self.fock_sparse_matrix, self.gorkov_sparse_matrix] 332 333 for parameter, list_of_parameters in self._diagonal.items(): 334 # Unpack dictionary 335 interaction = list_of_parameters["interaction"] 336 if interaction: 337 self.interacting = True 338 renormalise = list_of_parameters["renormalise"] 339 for parameters in list_of_parameters["parameters"]: 340 # Unpack dictionary 341 site = parameters["site"] 342 self.__dict__[parameter] = site.__dict__[parameter] 343 anomalous = parameters["anomalous"] 344 if anomalous: 345 self.anomalous = True 346 structure = parameters["structure"] 347 pauli_vector = parameters["pauli_vector"] 348 spin_orbital_matrix = parameters["spin_orbital_matrix"] 349 # Get {structure: [indices, alt indices]} 350 structure_dict = self.get_diagonal_indices(site=site, structure=structure, pauli_vector=pauli_vector) 351 if interaction: 352 _set_structure_dict(self.interaction_sparse_matrix, structure_dict, parameter, anomalous=False) 353 elif renormalise: 354 if anomalous: 355 _set_structure_dict(self.gorkov_sparse_matrix, structure_dict, parameter, anomalous) 356 else: 357 diagonals_dict, offdiagonals_dict = split_structure_dict(structure_dict) 358 _set_structure_dict(self.hartree_sparse_matrix, diagonals_dict, parameter, anomalous) 359 _set_structure_dict(self.fock_sparse_matrix, offdiagonals_dict, parameter, anomalous) 360 else: 361 _set_structure_dict(self.hamiltonian_sparse_matrix, structure_dict, parameter, anomalous) 362 363 for parameter, list_of_parameters in self._offdiagonal.items(): 364 # Unpack dictionary 365 interaction = list_of_parameters["interaction"] 366 if interaction: 367 self.interacting = True 368 renormalise = list_of_parameters["renormalise"] 369 for parameters in list_of_parameters["parameters"]: 370 # Unpack dictionary 371 site = parameters["site"] 372 self.__dict__[parameter] = site.__dict__[parameter] 373 neighbour = parameters["neighbour"] 374 hop_vector = parameters["vector"] 375 anomalous = parameters["anomalous"] 376 if anomalous: 377 self.anomalous = True 378 structure = parameters["structure"] 379 pauli_vector = parameters["pauli_vector"] 380 trs = parameters["trs"] 381 spin_orbital_matrix = parameters["spin_orbital_matrix"] 382 # Get {structure: [indices, alt indices]} 383 structure_dict = self.get_offdiagonal_indices(site=site, neighbour=neighbour, hop_vector=hop_vector, trs=trs, structure=structure, pauli_vector=pauli_vector) 384 if interaction: 385 _set_structure_dict(self.interaction_sparse_matrix, structure_dict, parameter, anomalous=False) 386 elif renormalise: 387 if anomalous: 388 _set_structure_dict(self.gorkov_sparse_matrix, structure_dict, parameter, anomalous) 389 else: 390 _set_structure_dict(self.hartree_sparse_matrix, structure_dict, parameter, anomalous) 391 _set_structure_dict(self.fock_sparse_matrix, structure_dict, parameter, anomalous) 392 raise ValueError("TODO: Split pauli_vector into diagonal and offdiagonal components") 393 else: 394 _set_structure_dict(self.hamiltonian_sparse_matrix, structure_dict, parameter, anomalous) 395 396 for _dict in self._sparse_matrices: 397 for key, items in _dict.items(): 398 value = self.__dict__[key] 399 for structure, map in items.items(): 400 map["update"] = value 401 items[structure]["indices"] = tensor(items[structure]["indices"]) 402 403 def _update_matrix(self): 404 405 if not hasattr(self, "_matrix"): 406 407 self._init_sparse_matrices() 408 409 if self.anomalous: 410 self.__dict__["_matrix"] = zeros([2*self.dof, 2*self.dof], dtype=COMPLEX) 411 else: 412 self.__dict__["_matrix"] = zeros([self.dof, self.dof], dtype=COMPLEX) 413 414 if self.interacting: 415 self.__dict__["interaction"] = zeros([self.dof, self.dof], dtype=COMPLEX) 416 417 for key, items in self.hamiltonian_sparse_matrix.items(): 418 value = self.__dict__[key] 419 for structure, map in items.items(): 420 [indices, alt_indices] = map["indices"] 421 update = map["update"] 422 anomalous = map["anomalous"] 423 if update != 0: 424 value = map["update"] 425 map["update"] = 0 426 _value = value * structure 427 if self.anomalous: 428 # Creates Hamiltonian of the following type: [[H0, Δ], [-H0, Δ*]] 429 if anomalous: 430 self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += _value 431 self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += scalar_conj(_value) 432 else: 433 self.__dict__["_matrix"][[indices, alt_indices]] += _value 434 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= _value 435 else: 436 self.__dict__["_matrix"][[indices, alt_indices]] += _value 437 438 for key, items in self.interaction_sparse_matrix.items(): 439 value = self.__dict__[key] 440 for structure, map in items.items(): 441 [indices, alt_indices] = map["indices"] 442 update = map["update"] 443 anomalous = map["anomalous"] 444 if update != 0: 445 value = map["update"] 446 map["update"] = 0 447 _value = value * structure 448 self.__dict__["interaction"][[indices, alt_indices]] += _value 449 450 for key, items in self.hartree_sparse_matrix.items(): 451 value = self.__dict__[key] 452 for structure, map in items.items(): 453 [indices, alt_indices] = map["indices"] 454 update = map["update"] 455 anomalous = map["anomalous"] 456 if update != 0: 457 value = map["update"] 458 map["update"] = 0 459 _value = value * structure 460 self.__dict__["_matrix"][[indices, alt_indices]] += _value 461 if self.anomalous: 462 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= _value 463 map["values"] = tensor([_value for _ in range(len(indices))]) 464 # elif "values" in map.keys(): 465 # self.__dict__["_matrix"][[indices, alt_indices]] += map["values"] 466 # if self.anomalous: 467 # self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= map["values"] 468 469 for key, items in self.fock_sparse_matrix.items(): 470 value = self.__dict__[key] 471 for structure, map in items.items(): 472 [indices, alt_indices] = map["indices"] 473 update = map["update"] 474 anomalous = map["anomalous"] 475 if update != 0: 476 value = map["update"] 477 map["update"] = 0 478 _value = value * structure 479 self.__dict__["_matrix"][[indices, alt_indices]] -= _value 480 if self.anomalous: 481 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] += _value 482 map["values"] = tensor([_value for _ in range(len(indices))]) 483 # elif "values" in map.keys(): 484 # self.__dict__["_matrix"][[indices, alt_indices]] -= map["values"] 485 # if self.anomalous: 486 # self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] += map["values"] 487 488 for key, items in self.gorkov_sparse_matrix.items(): 489 value = self.__dict__[key] 490 for structure, map in items.items(): 491 [indices, alt_indices] = map["indices"] 492 update = map["update"] 493 anomalous = map["anomalous"] 494 if update != 0: 495 value = map["update"] 496 map["update"] = 0 497 _value = value * structure 498 self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += _value 499 self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += scalar_conj(_value) 500 map["values"] = tensor([_value for _ in range(len(indices))]) 501 # elif "values" in map.keys(): 502 # self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += map["values"] 503 # self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += conj(map["values"]) 504 505 def extract_mean_fields(self): 506 507 # print("matrix:") 508 # print(self._matrix) 509 510 if not hasattr(self, "eigenvectors"): 511 raise Error("Run self.solve() before calling this function!") 512 for items in self.gorkov_sparse_matrix.values(): 513 for map in items.values(): 514 [indices, alt_indices] = map["indices"] 515 # TODO: add Fermi function 516 _alt_indices = alt_indices + self.dof 517 if "values" not in map: 518 map["values"] = 0 519 previous_values = deepcopy(map["values"]) 520 # map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, :self.dof], flip(self.eigenvectors[_alt_indices, self.dof:], [2])) 521 map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, :self.dof], self.eigenvectors[_alt_indices,: self.dof]) 522 values = map["values"] - previous_values 523 # print("##### gorkov #######") 524 # print(map["values"]) 525 # print("#################") 526 self.__dict__["_matrix"][[indices, _alt_indices]] += values 527 self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += conj(values) 528 529 for items in self.fock_sparse_matrix.values(): 530 for map in items.values(): 531 [indices, alt_indices] = map["indices"] 532 if "values" not in map: 533 map["values"] = 0 534 previous_values = deepcopy(map["values"]) 535 # TODO: add Fermi function 536 map["values"] = einsum('i,in,in->i', -self.interaction[indices,alt_indices], self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:])) 537 values = map["values"] - previous_values 538 self.__dict__["_matrix"][[indices, alt_indices]] -= values # Swap sign for anticommutativity!!! 539 if self.anomalous: 540 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] += values 541 542 for items in self.hartree_sparse_matrix.values(): 543 for map in items.values(): 544 [indices, alt_indices] = map["indices"] 545 # TODO: add Fermi function 546 if "values" not in map: 547 map["values"] = 0 548 previous_values = deepcopy(map["values"]) 549 map["values"] = einsum('ij,in,jn->i', self.interaction, self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:])) 550 # print("##### hartree #######") 551 # print(previous_values) 552 # print(map["values"]) 553 # print("#################") 554 # print(self._matrix) 555 values = map["values"] - previous_values 556 self.__dict__["_matrix"][[indices, alt_indices]] += values 557 if self.anomalous: 558 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= values 559 560 # print(self.hartree_sparse_matrix) 561 # print(self.fock_sparse_matrix) 562 # print(self.gorkov_sparse_matrix) 563 print(self._matrix) 564 565 # Normal: 566 # map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:])) 567 568 # def mean_field_step() 569 # for key, items in self.interaction_sparse_matrix.items(): 570 # value = self.__dict__[key] 571 # for structure, map in items.items(): 572 # [indices, alt_indices] = map["indices"] 573 # update = map["update"] 574 # anomalous = map["anomalous"] 575 # if update != 0: 576 # value = map["update"] 577 # map["update"] = 0 578 # _value = value * structure 579 # self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += _value 580 # self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += scalar_conj(_value) 581 # map["values"] = tensor([_value for _ in range(len(indices))]) 582 # else: 583 # self.__dict__["_matrix"][[indices, alt_indices + self.dof]] += map["values"] 584 # self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += conj(map["values"]) 585 586 ######################################################## 587 588 589 590 591 # def _get_indices(self, sites=None, hop=None, spin_flip=False, anomalous=False): 592 # _slice = [slice(None, None, None) for _ in range(self.nd)] 593 # cells = self.cells[*_slice] 594 595 # if type(hop)!=type(None): 596 # if len(hop)!=self.nd: 597 # raise ValueError("Hop must have same length as model dimensions!") 598 # for axis in range(self.nd): 599 # shift = hop[axis] 600 # cells = np.roll(cells, shift=-shift, axis=axis) 601 602 # cells = _Mapper(cells.tolist()) 603 604 # if type(sites)==type(None): 605 # sites = self.cell 606 607 # index = [0 for _ in range(self.nd)] 608 # index.append(0) 609 # test_cell = list(flatten(cells[*index])) 610 # sites_type = type(sites[0]) 611 # for site in sites: 612 # if type(site)!=sites_type: 613 # raise ValueError("All sites in sites must have same type!") 614 615 # def recursive(test_cell): 616 # test_cell = test_cell[0] 617 # structure_type = type(test_cell) 618 # if structure_type!=sites_type: 619 # _slice.append(slice(None, None, None)) 620 # test_cell = recursive(test_cell) 621 # return test_cell 622 623 # test_cell = recursive(test_cell) 624 625 # _slice.append(sites) 626 627 # if self.n_spins == 2: 628 # spin_flip = [0,1] 629 # if spin_flip: 630 # spin_flip = [1,0] 631 632 # _slice.append(spin_flip) 633 634 # cells = cells[*_slice] 635 636 # indices = cells.indices 637 # shape = np.shape(indices)[:-1] 638 # indices = flatten(cells.indices) 639 # if anomalous and self.anomalous: 640 # indices += self.n_leaves 641 642 # return list(indices), shape 643 644 def _slice_eigenvector(self, sites=None, hop=None, spin_flip=False, anomalous=False): 645 indices, shape = self._get_indices(sites=sites, hop=hop, spin_flip=spin_flip, anomalous=anomalous) 646 return self.eigenvectors[indices], shape 647 648 def density_matrix(self, temperature:float=0, sites=None, sites_to=None, hop=None, spin_flip=False, anomalous=False): 649 650 if type(sites)==type(None) and type(sites_to)!=type(None): 651 raise ValueError("If sites = None, then sites_to must also equal None!") 652 elif type(sites)!=type(None) and type(sites_to)==type(None): 653 sites_to = sites 654 elif type(sites)==type(None) and type(sites_to)==type(None): 655 sites = self.cell 656 sites_to = self.cell 657 elif len(sites)!=len(sites_to): 658 raise ValueError("sites and sites_to must have the same length!") 659 Ψ, shape = self._slice_eigenvector(sites=sites) 660 Ψc, shape = self._slice_eigenvector(sites=sites_to, hop=hop, spin_flip=spin_flip, anomalous=anomalous) 661 Ψc = conj(Ψc) 662 dm = multiply(Ψ, Ψc) 663 shape = shape + (Ψ.shape[-1],) 664 dm = reshape(dm, shape) 665 666 # if temperature==0 and not self.anomalous: 667 # return dm 668 # else: 669 return einsum('...n,n->...n', dm, self.fermi_distribution(temperature)) 670 671 def greens_function(self, energy:float or list, resolution:float, temperature:float=0, sites=None, sites_to=None, hop=None, spin_flip=False, anomalous=False, spin_polarised=False): 672 kwargs = {"temperature": temperature, "sites": sites, "sites_to": sites_to, "hop": hop, "spin_flip": spin_flip, "anomalous": anomalous} 673 dm = self.density_matrix(**kwargs) 674 try: 675 energy = energy[:, None] 676 poles = 1/(energy + 1.0j*resolution - self.eigenvalues) 677 ΨΨ = einsum('...n, en -> ...e', dm, poles) 678 except: 679 poles = 1/(energy + 1.0j*resolution - self.eigenvalues) 680 ΨΨ = einsum('...n, n -> ...', dm, poles) 681 if self.n_spins == 2 and not spin_polarised: 682 ΨΨ = sum(ΨΨ, axis=-1) 683 if (ΨΨ.shape)[-1] == 1: 684 ΨΨ = sum(ΨΨ, axis=-1) 685 return -1/pi*imag(ΨΨ)
def
fermi_distribution(self, temperature: float = 0):
102 def fermi_distribution(self, temperature:float=0): 103 if temperature == 0: 104 if self.anomalous: 105 distribution = ones(2*self.dof, dtype=COMPLEX) 106 else: 107 distribution = ones(self.dof, dtype=COMPLEX) 108 # Depopulate negative energy solutions 109 if self.anomalous: 110 ## Make use of BdG symmetry 111 distribution[..., :self.dof] = 0 112 else: 113 distribution[self.eigenvalues > 0] = 0 114 115 else: 116 distribution = 1/exp((self.eigenvalues/temperature)+1) 117 118 return distribution
def
get_diagonal_indices(self, site, structure, pauli_vector=None):
247 def get_diagonal_indices(self, site, structure, pauli_vector=None): 248 """ 249 Examples 250 -------- 251 252 >>> A = Atom(fractional_coordinates=[0, 0]) 253 >>> B = Atom(fractional_coordinates=[0.5, 0]) 254 >>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 255 >>> lattice.cut(number_of_cells=3, direction=0) 256 >>> lattice.cut(number_of_cells=3, direction=1) 257 >>> lattice.n_spins = 2 258 259 >>> lattice.__set_indices__() 260 261 >>> lattice.get_diagonal_indices(site=A, structure=-1) 262 {(-1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]} 263 264 >>> lattice.get_diagonal_indices(site=A, structure=-1, pauli_vector=[0,0,0,1]) 265 {(1+0j): [[0, 4, 8, 12, 16, 20, 24, 28, 32, 1, 5, 9, 13, 17, 21, 25, 29, 33], [1, 5, 9, 13, 17, 21, 25, 29, 33, 0, 4, 8, 12, 16, 20, 24, 28, 32]]} 266 """ 267 268 indices = self._get_lattice_resolved_indices(site) 269 270 if type(pauli_vector)!=type(None): 271 structure = self._get_spin_matrix_indices(indices, pauli_vector, structure) 272 else: 273 indices = flatten(list(indices)) 274 structure = {complex(structure): [indices, indices]} 275 276 return structure
Examples
>>> A = Atom(fractional_coordinates=[0, 0])
>>> B = Atom(fractional_coordinates=[0.5, 0])
>>> lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B])
>>> lattice.cut(number_of_cells=3, direction=0)
>>> lattice.cut(number_of_cells=3, direction=1)
>>> lattice.n_spins = 2
>>> lattice.__set_indices__()
>>> lattice.get_diagonal_indices(site=A, structure=-1)
{(-1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]}
>>> lattice.get_diagonal_indices(site=A, structure=-1, pauli_vector=[0,0,0,1])
{(1+0j): [[0, 4, 8, 12, 16, 20, 24, 28, 32, 1, 5, 9, 13, 17, 21, 25, 29, 33], [1, 5, 9, 13, 17, 21, 25, 29, 33, 0, 4, 8, 12, 16, 20, 24, 28, 32]]}
def
get_offdiagonal_indices( self, site, neighbour, hop_vector=None, trs=True, structure=1, pauli_vector=None):
278 def get_offdiagonal_indices(self, site, neighbour, hop_vector=None, trs=True, structure=1, pauli_vector=None): 279 """ 280 Examples 281 -------- 282 283 >>> lattice.get_offdiagonal_indices(site=A, neighbour=B, hop_vector=[1,0], trs=True, structure=1j) 284 {1j: [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11]], -1j: [[14, 15, 18, 19, 22, 285 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]} 286 >>> lattice.get_offdiagonal_indices(site=A, neighbour=B, hop_vector=[1,0], trs=True, structure=1) 287 {(1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11], [14, 15, 18, 19, 22, 23, 288 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11, 0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 1 289 0, 11]]} 290 291 TODO add parameter to flip spin with Pauli vector 292 293 """ 294 if type(pauli_vector)!=type(None): 295 raise ValueError("Pauli vector not yet implemented!") 296 297 indices = self._get_lattice_resolved_indices(site) 298 alt_indices = self._get_lattice_resolved_indices(neighbour) 299 300 if type(hop_vector)!=type(None): 301 alt_indices = self._get_rolled_indices(alt_indices, hop_vector) 302 303 indices = flatten(list(indices)) 304 alt_indices = flatten(list(alt_indices)) 305 structure_dict = {} 306 if trs: 307 _indices = deepcopy(indices) 308 _alt_indices = deepcopy(alt_indices) 309 310 _set_dict(structure_dict, structure, indices, alt_indices) 311 312 if trs: 313 structure = scalar_conj(structure) 314 _set_dict(structure_dict, structure, _alt_indices, _indices) 315 316 return structure_dict
Examples
>>> lattice.get_offdiagonal_indices(site=A, neighbour=B, hop_vector=[1,0], trs=True, structure=1j)
{1j: [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33], [14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11]], -1j: [[14, 15, 18, 19, 22,
23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11], [0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33]]}
>>> lattice.get_offdiagonal_indices(site=A, neighbour=B, hop_vector=[1,0], trs=True, structure=1)
{(1+0j): [[0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11], [14, 15, 18, 19, 22, 23,
26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 10, 11, 0, 1, 4, 5, 8, 9, 12, 13, 16, 17, 20, 21, 24, 25, 28, 29, 32, 33, 14, 15, 18, 19, 22, 23, 26, 27, 30, 31, 34, 35, 2, 3, 6, 7, 1
0, 11]]}
TODO add parameter to flip spin with Pauli vector
def
extract_mean_fields(self):
505 def extract_mean_fields(self): 506 507 # print("matrix:") 508 # print(self._matrix) 509 510 if not hasattr(self, "eigenvectors"): 511 raise Error("Run self.solve() before calling this function!") 512 for items in self.gorkov_sparse_matrix.values(): 513 for map in items.values(): 514 [indices, alt_indices] = map["indices"] 515 # TODO: add Fermi function 516 _alt_indices = alt_indices + self.dof 517 if "values" not in map: 518 map["values"] = 0 519 previous_values = deepcopy(map["values"]) 520 # map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, :self.dof], flip(self.eigenvectors[_alt_indices, self.dof:], [2])) 521 map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, :self.dof], self.eigenvectors[_alt_indices,: self.dof]) 522 values = map["values"] - previous_values 523 # print("##### gorkov #######") 524 # print(map["values"]) 525 # print("#################") 526 self.__dict__["_matrix"][[indices, _alt_indices]] += values 527 self.__dict__["_matrix"][[indices + self.dof, alt_indices]] += conj(values) 528 529 for items in self.fock_sparse_matrix.values(): 530 for map in items.values(): 531 [indices, alt_indices] = map["indices"] 532 if "values" not in map: 533 map["values"] = 0 534 previous_values = deepcopy(map["values"]) 535 # TODO: add Fermi function 536 map["values"] = einsum('i,in,in->i', -self.interaction[indices,alt_indices], self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:])) 537 values = map["values"] - previous_values 538 self.__dict__["_matrix"][[indices, alt_indices]] -= values # Swap sign for anticommutativity!!! 539 if self.anomalous: 540 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] += values 541 542 for items in self.hartree_sparse_matrix.values(): 543 for map in items.values(): 544 [indices, alt_indices] = map["indices"] 545 # TODO: add Fermi function 546 if "values" not in map: 547 map["values"] = 0 548 previous_values = deepcopy(map["values"]) 549 map["values"] = einsum('ij,in,jn->i', self.interaction, self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:])) 550 # print("##### hartree #######") 551 # print(previous_values) 552 # print(map["values"]) 553 # print("#################") 554 # print(self._matrix) 555 values = map["values"] - previous_values 556 self.__dict__["_matrix"][[indices, alt_indices]] += values 557 if self.anomalous: 558 self.__dict__["_matrix"][[indices + self.dof, alt_indices + self.dof]] -= values 559 560 # print(self.hartree_sparse_matrix) 561 # print(self.fock_sparse_matrix) 562 # print(self.gorkov_sparse_matrix) 563 print(self._matrix) 564 565 # Normal: 566 # map["values"] = einsum('i,in,in->i', self.interaction[indices,alt_indices], self.eigenvectors[indices, self.dof:], conj(self.eigenvectors[alt_indices, self.dof:]))
def
density_matrix( self, temperature: float = 0, sites=None, sites_to=None, hop=None, spin_flip=False, anomalous=False):
648 def density_matrix(self, temperature:float=0, sites=None, sites_to=None, hop=None, spin_flip=False, anomalous=False): 649 650 if type(sites)==type(None) and type(sites_to)!=type(None): 651 raise ValueError("If sites = None, then sites_to must also equal None!") 652 elif type(sites)!=type(None) and type(sites_to)==type(None): 653 sites_to = sites 654 elif type(sites)==type(None) and type(sites_to)==type(None): 655 sites = self.cell 656 sites_to = self.cell 657 elif len(sites)!=len(sites_to): 658 raise ValueError("sites and sites_to must have the same length!") 659 Ψ, shape = self._slice_eigenvector(sites=sites) 660 Ψc, shape = self._slice_eigenvector(sites=sites_to, hop=hop, spin_flip=spin_flip, anomalous=anomalous) 661 Ψc = conj(Ψc) 662 dm = multiply(Ψ, Ψc) 663 shape = shape + (Ψ.shape[-1],) 664 dm = reshape(dm, shape) 665 666 # if temperature==0 and not self.anomalous: 667 # return dm 668 # else: 669 return einsum('...n,n->...n', dm, self.fermi_distribution(temperature))
def
greens_function( self, energy: float, resolution: float, temperature: float = 0, sites=None, sites_to=None, hop=None, spin_flip=False, anomalous=False, spin_polarised=False):
671 def greens_function(self, energy:float or list, resolution:float, temperature:float=0, sites=None, sites_to=None, hop=None, spin_flip=False, anomalous=False, spin_polarised=False): 672 kwargs = {"temperature": temperature, "sites": sites, "sites_to": sites_to, "hop": hop, "spin_flip": spin_flip, "anomalous": anomalous} 673 dm = self.density_matrix(**kwargs) 674 try: 675 energy = energy[:, None] 676 poles = 1/(energy + 1.0j*resolution - self.eigenvalues) 677 ΨΨ = einsum('...n, en -> ...e', dm, poles) 678 except: 679 poles = 1/(energy + 1.0j*resolution - self.eigenvalues) 680 ΨΨ = einsum('...n, n -> ...', dm, poles) 681 if self.n_spins == 2 and not spin_polarised: 682 ΨΨ = sum(ΨΨ, axis=-1) 683 if (ΨΨ.shape)[-1] == 1: 684 ΨΨ = sum(ΨΨ, axis=-1) 685 return -1/pi*imag(ΨΨ)
Inherited Members
687class BdGPatch(StatisticalMechanics): 688 689 def __setattr__(self, name, value, update = True): 690 """Set attributes down the lineage""" 691 if name=="n_spins": 692 self._set_n_spins(value) 693 self.__dict__[f"{name}"] = value 694 else: 695 try: 696 previous_value = self.__dict__[f"{name}"] 697 for _dict in self._sparse_matrices: 698 items = _dict[name] 699 for structure in items.keys(): 700 items[structure]["update"] = value - previous_value 701 except: 702 pass 703 self.__dict__[f"{name}"] = value 704 [item.__setattr__(name, value, update = False) for item in self] 705 706 # @property 707 # def spins(self): 708 # return _Mapper(self) 709 710 def _set_n_spins(self, n_spins): 711 if n_spins<1 or n_spins>2: 712 raise ValueError("n_spins must be 1 or 2") 713 try: 714 if self.__dict__["_n_spins"]==2 and n_spins==1: 715 self.clear() 716 [leaf.parent.clear() for leaf in self.leaves] 717 except: 718 pass 719 if n_spins==2: 720 [leaf.extend([Spin('↑'), Spin('↓')]) for leaf in self.leaves] 721 722 def add(self, structure=1, pauli_vector:list=None, spin_orbital_matrix=None, anomalous:bool=False, interaction=False, renormalise=False, **kwarg): 723 """Add vertex term as a parameter described by a key-value pair""" 724 725 parameters = {"site": self, "anomalous": anomalous, "structure": structure, "pauli_vector": pauli_vector, "spin_orbital_matrix": spin_orbital_matrix} 726 727 if len(kwarg)==0: 728 raise ValueError("Requires a kwarg!") 729 730 if len(kwarg)>1: 731 raise ValueError("Accepts only one kwarg!") 732 733 key, value = list(kwarg.items())[0] 734 735 if type(value)!=float and type(value)!=int: 736 raise ValueError("Value must be real! (In order for the symmetries to apply correctly, put the imaginary number in structure and the magnitude in the value)") 737 738 if not key in self._diagonal: 739 self._diagonal[key] = {} 740 741 if not "parameters" in self._diagonal[key].keys(): 742 self._diagonal[key]["parameters"] = [] 743 744 self._diagonal[key]["parameters"].append(parameters) 745 self._diagonal[key]["interaction"] = interaction 746 self._diagonal[key]["renormalise"] = renormalise 747 748 self.__setattr__(key, value) 749 750 def hop(self, neighbour: Tree=None, structure=1, pauli_vector:list=None, spin_orbital_matrix=None, vector=None, trs:bool=True, anomalous:bool=False, interaction=False, renormalise=False, **kwarg): 751 """Add directed edge as a neighbour and a parameter (described by a 752 key-value pair) 753 """ 754 parameters = {"site": self, "neighbour": neighbour, "vector": vector, "trs": trs, "anomalous": anomalous, "structure": structure, "pauli_vector": pauli_vector, "spin_orbital_matrix": spin_orbital_matrix} 755 756 if len(kwarg)==0: 757 raise ValueError("Requires a kwarg!") 758 759 if len(kwarg)>1: 760 raise ValueError("Accepts only one kwarg!") 761 762 if type(neighbour)!=type(None): 763 if self.n_leaves!=neighbour.n_leaves: 764 raise ValueError("Self and neighbour have different number of leaves. Hopping is not well-defined!") 765 766 key, value = list(kwarg.items())[0] 767 768 if not key in self._offdiagonal: 769 self._offdiagonal[key] = {} 770 771 if not "parameters" in self._offdiagonal[key].keys(): 772 self._offdiagonal[key]["parameters"] = [] 773 774 self._offdiagonal[key]["parameters"].append(parameters) 775 self._offdiagonal[key]["interaction"] = interaction 776 self._offdiagonal[key]["renormalise"] = renormalise 777 778 self.__setattr__(key, value) 779 780 def is_anomalous(self): 781 if not hasattr(self, "anomalous"): 782 diagonals, offdiagonals = self._diagonal, self._offdiagonal 783 784 self.anomalous = False 785 for items in list(diagonals.values())+list(offdiagonals.values()): 786 for item in items["parameters"]: 787 if item["anomalous"]: 788 self.anomalous = True 789 break 790 791 return self.anomalous 792 793 # def has_interactions(self): 794 # if not hasattr(self, "interacting"): 795 # diagonals, offdiagonals = self._diagonal, self._offdiagonal 796 797 # self.interacting = False 798 # for items in list(diagonals.values())+list(offdiagonals.values()): 799 # if items["interaction"]: 800 # self.interacting = True 801 # break 802 # return self.interacting 803 804 def has_renormalisation_fields(self): 805 if not hasattr(self, "renormalise"): 806 diagonals, offdiagonals = self._diagonal, self._offdiagonal 807 808 self.renormalise = False 809 for items in list(diagonals.values())+list(offdiagonals.values()): 810 if items["renormalise"]: 811 self.renormalise = True 812 break 813 return self.renormalise 814 815 # def _get_spin_matrix_indices(self): 816 # if self.height == 1: 817 # _indices = self.indices 818 # _indices = [[x, y] for y in _indices for x in _indices] 819 # return _indices 820 # elif self.height>1: 821 # _indices = [] 822 # [_indices.extend(site._get_spin_matrix_indices()) for site in self] 823 # print(_indices) 824 # exit() 825 # return _indices 826 827 def read_value(self, list_of_parameters): 828 value = list_of_parameters["value"] 829 try: 830 value -= list_of_parameters["previous_value"] 831 except: 832 pass 833 return value 834 835 def _get_spin_indices(self, pauli_vector): 836 cache = {} 837 if type(pauli_vector)==list: 838 pauli_vector = tensor(pauli_vector, dtype=COMPLEX) 839 pauli_vector = einsum('i,ijk->jk', pauli_vector, pauli_vec) 840 if self.height>0: 841 spin_matrix_indices = self._get_spin_matrix_indices() 842 pauli_vector = pauli_vector.flatten() 843 for index in range(len(spin_matrix_indices)): 844 value = complex(pauli_vector[index%4]) 845 if abs(value)!=0: 846 _indices = spin_matrix_indices[index] 847 try: 848 cache[value][0].append(_indices[0]) 849 cache[value][1].append(_indices[1]) 850 except: 851 cache[value] = [[_indices[0]], [_indices[1]]] 852 return cache 853 854 def anomalous_cache(self, cache, anomalous): 855 if anomalous: 856 _cache = cache.clone().detach() 857 _cache[1,:] += self.n_leaves 858 _alt_cache = cache.clone().detach() 859 _alt_cache[0,:] += self.n_leaves 860 cache = [_cache, _alt_cache, anomalous] 861 else: 862 cache = [cache, cache + self.n_leaves, anomalous] 863 return cache 864 865 #def read_indices(self, list_of_parameters): 866 # try: 867 # return list_of_parameters["cache"] 868 # except: 869 # pass 870 871 # def set_cache(cache, structure, value, anomalous): 872 # try: 873 # cache[structure] = cat([cache[structure], value], axis=-1) 874 # except: 875 # cache[structure] = value 876 877 # if self.is_anomalous(): 878 # cache[structure] = self.anomalous_cache(cache[structure], anomalous) 879 880 # cache = {} 881 882 # for parameters in list_of_parameters["parameters"]: 883 # site = parameters["site"] 884 # structure = parameters["structure"] 885 # pauli_vector = parameters["pauli_vector"] 886 # spin_orbital_matrix = parameters["spin_orbital_matrix"] 887 # anomalous = parameters["anomalous"] 888 # try: 889 # vector = parameters["vector"] 890 # except: 891 # vector = None 892 # try: 893 # neighbour = parameters["neighbour"] 894 # except: 895 # neighbour = None 896 # try: 897 # trs = parameters["trs"] 898 # except: 899 # trs = False 900 901 # ################################################################ 902 # if type(vector)==type(None) and type(neighbour)==type(None): 903 # if type(spin_orbital_matrix)!=type(None): 904 # # indices = site._get_spin_orbital_indices(pauli_vector) 905 # raise ValueError("Not yet implemented!") 906 907 # elif type(pauli_vector)!=type(None): 908 # pauli_vector = tensor(pauli_vector, dtype=COMPLEX) 909 # pauli_vector *= structure 910 # for key, value in site._get_spin_indices(pauli_vector).items(): 911 # structure = key 912 # value = tensor(value) 913 # set_cache(cache, structure, value, anomalous) 914 # else: 915 # value = tensor([site.indices, site.indices]) 916 # set_cache(cache, structure, value, anomalous) 917 # ################################################################ 918 # else: 919 # if neighbour == None: 920 # neighbour = site 921 # vector = parameters["vector"] 922 923 # _slice = [slice(None, None, None) for _ in range(self.nd)] 924 # _neighbour_slice = deepcopy(_slice) 925 926 # if type(site)==Cell: 927 # cells = _Mapper([site]) 928 # neighbour_cells = _Mapper([neighbour]) 929 # elif type(site)!=Lattice: 930 # if site.depth>1: 931 # [_slice.append(slice(None,None,None)) for _ in range(site.height-1)] 932 # [_neighbour_slice.append(slice(None,None,None)) for _ in range(neighbour.height-1)] 933 # _slice.append(site) 934 # _neighbour_slice.append(neighbour) 935 # cells = self[*_slice] 936 # neighbour_cells = self[*_neighbour_slice] 937 # else: 938 # cells = self.cells 939 # neighbour_cells = self.cells 940 941 # if type(vector)!=type(None): 942 # _shape = self.cells.shape 943 # _neighbour_cells = np.empty(_shape, dtype=MatrixTree) 944 # nslice = [slice(None, None, None) for _ in range(self.nd)] 945 # _neighbour_cells[*nslice] = neighbour_cells 946 # for axis in range(self.nd): 947 # shift = vector[axis] 948 # _neighbour_cells = np.roll(_neighbour_cells, shift=shift, axis=axis) 949 # neighbour_cells = _neighbour_cells 950 # # else: 951 # # nslice = [slice(None, None, None) for _ in range(self.nd)] 952 # # cells = _cells[*nslice] 953 # # neighbour_cells = _neighbour_cells[*nslice] 954 955 # if type(pauli_vector)!=type(None): 956 957 # # Neighbour indices: 958 # _indices = [site._get_spin_indices(pauli_vector) for site in flatten(self.cells.tolist())] 959 # _neighbour_indices = [site._get_spin_indices(pauli_vector) for site in flatten(neighbour_cells)] 960 961 # keys = list(_indices[0].keys()) 962 963 # def func(x, y, key): 964 # try: 965 # x[key][0].extend(y[key][0]) 966 # x[key][1].extend(y[key][1]) 967 # except: 968 # x[key] = y[key] 969 970 # __indices = {} 971 # __neighbour_indices = {} 972 # [func(__indices, _index, key) for key in keys for _index in _indices] 973 # [func(__neighbour_indices, _index, key) for key in keys for _index in _neighbour_indices] 974 975 # indices = {} 976 # for key in __indices.keys(): 977 # indices[key] = [__indices[key][0], __neighbour_indices[key][1]] 978 979 # for key, value in indices.items(): 980 # structure = key 981 # value = tensor(value) 982 # set_cache(cache, structure, value, anomalous) 983 984 # else: 985 # try: 986 # indices = flatten(list(_Mapper(cells.tolist()).indices)) 987 # neighbour_indices = flatten(list(_Mapper(neighbour_cells.tolist()).indices)) 988 # except: 989 # indices = flatten(list(cells.indices)) 990 # neighbour_indices = flatten(list(neighbour.indices)) 991 # value = tensor([indices, neighbour_indices]) 992 # set_cache(cache, structure, value, anomalous) 993 # if trs: 994 # value = tensor([neighbour_indices, indices]) 995 # set_cache(cache, scalar_conj(structure), value, anomalous) 996 # # name = neighbour.__name__ 997 # # if type(pauli_vector)==type(None): 998 # # if site.depth == 0: 999 # # neighbour_indices = flatten([cell.indices for cell in cells]) 1000 # # elif site.depth == 1: 1001 # # neighbour_indices = flatten([cell[name].indices for cell in cells]) 1002 # # elif site.depth > 1: 1003 # # _slice = [slice(None, None, None) for _ in range(site.depth-1)] 1004 # # _slice.append(name) 1005 # # neighbour_indices = flatten([cell[*_slice].indices for cell in cells]) 1006 1007 # # if type(spin_orbital_matrix)!=type(None): 1008 # # # indices = site._get_spin_orbital_indices(pauli_vector) 1009 # # raise ValueError("Not yet implemented!") 1010 1011 # # elif type(pauli_vector)!=type(None): 1012 # # print(vector) 1013 # # exit() 1014 # # pauli_vector = tensor(pauli_vector, dtype=COMPLEX) 1015 # # pauli_vector *= structure 1016 # # for key, value in site._get_spin_indices(pauli_vector).items(): 1017 # # structure = key 1018 # # value = tensor(value) 1019 # # else: 1020 # # indices = site.indices 1021 # # print(cells) 1022 # # exit() 1023 # # exit() 1024 # # print(neighbour_indices) 1025 # # exit() 1026 # # print(self) 1027 # # print(site) 1028 # # exit() 1029 1030 # list_of_parameters["cache"] = cache 1031 1032 # return list_of_parameters["cache"] 1033 1034 @property 1035 def matrix(self): 1036 return self.get_matrix() 1037 1038 # @property 1039 def get_matrix(self): 1040 """ 1041 Construct the matrix representation of the tree 1042 1043 >>> A = MatrixTree() 1044 >>> B = MatrixTree() 1045 >>> matrix = MatrixTree(A, B) 1046 >>> matrix.add(μ=0.6) 1047 >>> A.hop(B, t=1) 1048 >>> matrix.matrix 1049 array([[0.6, 1. ], 1050 [0. , 0.6]]) 1051 1052 >>> A.t = 3 1053 >>> matrix.matrix 1054 array([[0.6, 3. ], 1055 [0. , 0.6]]) 1056 1057 >>> A.t *= 3 1058 >>> matrix.matrix 1059 array([[0.6, 9. ], 1060 [0. , 0.6]]) 1061 """ 1062 1063 self._update_matrix() 1064 1065 return self._matrix 1066 1067 def self_consistent_step(self): 1068 1069 matrix = self.get_matrix() 1070 1071 print(matrix) 1072 1073 exit()
def
add( self, structure=1, pauli_vector: list = None, spin_orbital_matrix=None, anomalous: bool = False, interaction=False, renormalise=False, **kwarg):
722 def add(self, structure=1, pauli_vector:list=None, spin_orbital_matrix=None, anomalous:bool=False, interaction=False, renormalise=False, **kwarg): 723 """Add vertex term as a parameter described by a key-value pair""" 724 725 parameters = {"site": self, "anomalous": anomalous, "structure": structure, "pauli_vector": pauli_vector, "spin_orbital_matrix": spin_orbital_matrix} 726 727 if len(kwarg)==0: 728 raise ValueError("Requires a kwarg!") 729 730 if len(kwarg)>1: 731 raise ValueError("Accepts only one kwarg!") 732 733 key, value = list(kwarg.items())[0] 734 735 if type(value)!=float and type(value)!=int: 736 raise ValueError("Value must be real! (In order for the symmetries to apply correctly, put the imaginary number in structure and the magnitude in the value)") 737 738 if not key in self._diagonal: 739 self._diagonal[key] = {} 740 741 if not "parameters" in self._diagonal[key].keys(): 742 self._diagonal[key]["parameters"] = [] 743 744 self._diagonal[key]["parameters"].append(parameters) 745 self._diagonal[key]["interaction"] = interaction 746 self._diagonal[key]["renormalise"] = renormalise 747 748 self.__setattr__(key, value)
Add vertex term as a parameter described by a key-value pair
def
hop( self, neighbour: tree.network.Tree = None, structure=1, pauli_vector: list = None, spin_orbital_matrix=None, vector=None, trs: bool = True, anomalous: bool = False, interaction=False, renormalise=False, **kwarg):
750 def hop(self, neighbour: Tree=None, structure=1, pauli_vector:list=None, spin_orbital_matrix=None, vector=None, trs:bool=True, anomalous:bool=False, interaction=False, renormalise=False, **kwarg): 751 """Add directed edge as a neighbour and a parameter (described by a 752 key-value pair) 753 """ 754 parameters = {"site": self, "neighbour": neighbour, "vector": vector, "trs": trs, "anomalous": anomalous, "structure": structure, "pauli_vector": pauli_vector, "spin_orbital_matrix": spin_orbital_matrix} 755 756 if len(kwarg)==0: 757 raise ValueError("Requires a kwarg!") 758 759 if len(kwarg)>1: 760 raise ValueError("Accepts only one kwarg!") 761 762 if type(neighbour)!=type(None): 763 if self.n_leaves!=neighbour.n_leaves: 764 raise ValueError("Self and neighbour have different number of leaves. Hopping is not well-defined!") 765 766 key, value = list(kwarg.items())[0] 767 768 if not key in self._offdiagonal: 769 self._offdiagonal[key] = {} 770 771 if not "parameters" in self._offdiagonal[key].keys(): 772 self._offdiagonal[key]["parameters"] = [] 773 774 self._offdiagonal[key]["parameters"].append(parameters) 775 self._offdiagonal[key]["interaction"] = interaction 776 self._offdiagonal[key]["renormalise"] = renormalise 777 778 self.__setattr__(key, value)
Add directed edge as a neighbour and a parameter (described by a key-value pair)
def
is_anomalous(self):
780 def is_anomalous(self): 781 if not hasattr(self, "anomalous"): 782 diagonals, offdiagonals = self._diagonal, self._offdiagonal 783 784 self.anomalous = False 785 for items in list(diagonals.values())+list(offdiagonals.values()): 786 for item in items["parameters"]: 787 if item["anomalous"]: 788 self.anomalous = True 789 break 790 791 return self.anomalous
def
has_renormalisation_fields(self):
804 def has_renormalisation_fields(self): 805 if not hasattr(self, "renormalise"): 806 diagonals, offdiagonals = self._diagonal, self._offdiagonal 807 808 self.renormalise = False 809 for items in list(diagonals.values())+list(offdiagonals.values()): 810 if items["renormalise"]: 811 self.renormalise = True 812 break 813 return self.renormalise
def
anomalous_cache(self, cache, anomalous):
854 def anomalous_cache(self, cache, anomalous): 855 if anomalous: 856 _cache = cache.clone().detach() 857 _cache[1,:] += self.n_leaves 858 _alt_cache = cache.clone().detach() 859 _alt_cache[0,:] += self.n_leaves 860 cache = [_cache, _alt_cache, anomalous] 861 else: 862 cache = [cache, cache + self.n_leaves, anomalous] 863 return cache
def
get_matrix(self):
1039 def get_matrix(self): 1040 """ 1041 Construct the matrix representation of the tree 1042 1043 >>> A = MatrixTree() 1044 >>> B = MatrixTree() 1045 >>> matrix = MatrixTree(A, B) 1046 >>> matrix.add(μ=0.6) 1047 >>> A.hop(B, t=1) 1048 >>> matrix.matrix 1049 array([[0.6, 1. ], 1050 [0. , 0.6]]) 1051 1052 >>> A.t = 3 1053 >>> matrix.matrix 1054 array([[0.6, 3. ], 1055 [0. , 0.6]]) 1056 1057 >>> A.t *= 3 1058 >>> matrix.matrix 1059 array([[0.6, 9. ], 1060 [0. , 0.6]]) 1061 """ 1062 1063 self._update_matrix() 1064 1065 return self._matrix
Construct the matrix representation of the tree
>>> A = MatrixTree()
>>> B = MatrixTree()
>>> matrix = MatrixTree(A, B)
>>> matrix.add(μ=0.6)
>>> A.hop(B, t=1)
>>> matrix.matrix
array([[0.6, 1. ],
[0. , 0.6]])
>>> A.t = 3
>>> matrix.matrix
array([[0.6, 3. ],
[0. , 0.6]])
>>> A.t *= 3
>>> matrix.matrix
array([[0.6, 9. ],
[0. , 0.6]])
1075class Spin(BdGPatch, MatrixTree): 1076 """ 1077 A class with a name: ↑ and ↓ spin, belonging 1078 as a child to an Orbital. 1079 1080 Attributes 1081 ---------- 1082 None 1083 1084 Methods 1085 ------- 1086 None 1087 1088 Examples 1089 -------- 1090 >>> spin = Spin('up') 1091 >>> print(spin) 1092 Spin(↑) 1093 1094 >>> spin = Spin('↓') 1095 >>> print(spin) 1096 Spin(↓) 1097 1098 >>> spin = Spin(1) 1099 >>> print(spin) 1100 Spin(↓) 1101 1102 Note 1103 ---- 1104 This class acts as a template for spin indices. 1105 """ 1106 def __init__(self, name='↑'): 1107 """Inherit list class method and add parent attribute to items""" 1108 super().__init__() 1109 if name in ['up', '↑', 0]: 1110 self.__name__ = '↑' 1111 elif name in ['dn', 'down', '↓', 1]: 1112 self.__name__ = '↓' 1113 else: 1114 raise ValueError('Varname must be up or dn!') 1115 self.__dict__["parent"] = None 1116 self.__dict__["_idx"] = None 1117 1118 def __str__(self): 1119 """Returns object class name""" 1120 return f'{self.__class__.__name__}({self.__name__})' 1121 1122 def __repr__(self): 1123 """Returns object class name""" 1124 return f'{self.__class__.__name__}({self.__name__})'
A class with a name: ↑ and ↓ spin, belonging as a child to an Orbital.
Attributes
None
Methods
None
Examples
>>> spin = Spin('up')
>>> print(spin)
Spin(↑)
>>> spin = Spin('↓')
>>> print(spin)
Spin(↓)
>>> spin = Spin(1)
>>> print(spin)
Spin(↓)
Note
This class acts as a template for spin indices.
Spin(name='↑')
1106 def __init__(self, name='↑'): 1107 """Inherit list class method and add parent attribute to items""" 1108 super().__init__() 1109 if name in ['up', '↑', 0]: 1110 self.__name__ = '↑' 1111 elif name in ['dn', 'down', '↓', 1]: 1112 self.__name__ = '↓' 1113 else: 1114 raise ValueError('Varname must be up or dn!') 1115 self.__dict__["parent"] = None 1116 self.__dict__["_idx"] = None
Inherit list class method and add parent attribute to items
Inherited Members
1126class Orbital(BdGPatch, MatrixTree, object): 1127 """ 1128 Orbital representation 1129 1130 Args 1131 ---- 1132 n (int): principal quantum number, determines the energy level and size of the orbital 1133 l (int): azimuthal quantum number, defines the shape of the orbital 1134 m (int): magnetic quantum number, defines the orientation of the orbital 1135 1136 Attributes 1137 ---------- 1138 1139 Methods 1140 ------- 1141 radial_function(self, r) 1142 angular_function(self, θ, φ) 1143 compute_wavefunction_cross_section(self, horizontal_interval, vertical_interval) 1144 compute_probability_density(self) 1145 plot_probability_density_cross_section(self, a0, dark_theme=False, colormap='rocket') 1146 TODO add 3D density plot 1147 1148 Examples 1149 -------- 1150 >>> orbital = Orbital(3,2,0) 1151 >>> print(orbital) 1152 Orbital(3d_{z^2}) 1153 1154 Note 1155 ---- 1156 Integrates over the unspecified axis 1157 1158 Reference 1159 --------- 1160 1. https://ssebastianmag.medium.com/computational-physics-with-python-hydrogen-wavefunctions-electron-density-plots-8fede44b7b12 1161 2. https://medium.com/@bldevries/a-symphony-of-spheres-animating-spherical-harmonics-in-blender-with-python-and-shape-keys-aa67b7ff3d93 1162 1163 """ 1164 def __init__(self, n:int=0, l:int=0, m:int=0): 1165 """Inherit list class method and add parent attribute to items""" 1166 1167 # Quantum numbers validation 1168 if not isinstance(n, int) or n < 1: 1169 raise ValueError('n should be an integer satisfying the condition: n >= 1') 1170 if not isinstance(l, int) or not (0 <= l < n): 1171 raise ValueError('l should be an integer satisfying the condition: 0 <= l < n') 1172 if not isinstance(m, int) or not (-l <= m <= l): 1173 raise ValueError('m should be an integer satisfying the condition: -l <= m <= l') 1174 1175 self.__dict__["quantum_numbers"] = [n,l,m] 1176 self.__dict__["n"] = n 1177 self.__dict__["l"] = l 1178 self.__dict__["m"] = m 1179 1180 super().__init__() 1181 self.__dict__["parent"] = None 1182 self.__dict__["_idx"] = None 1183 1184 def __str__(self): 1185 """Returns object class name""" 1186 return f'{self.__class__.__name__}({self.__name__})' 1187 1188 def __repr__(self): 1189 """Returns object class name""" 1190 return f'{self.__class__.__name__}({self.__name__})' 1191 1192 @property 1193 def __name__(self): 1194 n, l, m = self.n, self.l, self.m 1195 name = n 1196 if l==0: 1197 if m==0: 1198 name = 's' 1199 elif l==1: 1200 if m==0: 1201 name ='p_z' 1202 elif m==1: 1203 name ='p_x' 1204 elif m==-1: 1205 name ='p_y' 1206 elif l==2: 1207 if m==0: 1208 name ='d_{z^2}' 1209 elif m==1: 1210 name ='d_{xz}' 1211 elif m==-1: 1212 name ='d_{yz}' 1213 elif m==2: 1214 name ='d_{xy}' 1215 elif m==-2: 1216 name ='d_{x^2-z^2}' 1217 elif l==3: 1218 if m==0: 1219 name ='f_{z^3}' 1220 elif m==1: 1221 name ='f_{xz^2}' 1222 elif m==-1: 1223 name ='f_{yz^2}' 1224 elif m==2: 1225 name ='f_{xyz}' 1226 elif m==-2: 1227 name ='f_{z(x^2-z^2)}' 1228 elif m==3: 1229 name ='f_{x(x^2-3y^2)}' 1230 elif m==-3: 1231 name ='f_{y(3x^2-y^2)}' 1232 if l<=3: 1233 name=str(n)+name 1234 else: 1235 name = f'({n} {l} {m})' 1236 return name 1237 1238 def radial_function(self, r): 1239 """Compute the normalized radial part of the wavefunction using 1240 Laguerre polynomials and an exponential decay factor. 1241 Args: 1242 r (numpy.ndarray): radial coordinate 1243 Returns: 1244 numpy.ndarray: wavefunction radial component 1245 """ 1246 n = self.n 1247 l = self.l 1248 laguerre = sp.genlaguerre(n-l-1, 2*l+1) 1249 p = 2*r/n 1250 1251 constant_factor = np.sqrt( 1252 ((2/n)**3 * (sp.factorial(n-l-1))) / (2*n*(sp.factorial(n+l))) 1253 ) 1254 return constant_factor * np.exp(-p/2) * (p**l) * laguerre(p) 1255 1256 def compute_radial_density(self, rr): 1257 """ 1258 Compute the radial probability density for the 1259 list of radial points rr. 1260 """ 1261 area = 4*np.pi*np.square(rr) 1262 radius = np.square(self.radial_function(rr)) 1263 self.radial_density = radius * area 1264 return self.radial_density 1265 1266 def plot_radial_density(self, ax, rr): 1267 """ 1268 Compute the radial probability density for the 1269 list of radial points rr. 1270 """ 1271 if not hasattr(self, 'radial_density'): 1272 print('Ψ not computed. Computing...') 1273 self.compute_radial_density(rr) 1274 print('Ψ computed.') 1275 1276 ax.plot(rr, self.radial_density, label=self.__name__) 1277 return ax 1278 1279 def angular_function(self, θ, φ): 1280 """ Compute the normalized angular part of the wavefunction using 1281 Legendre polynomials and a phase-shifting exponential factor. 1282 Args: 1283 θ (numpy.ndarray): polar angle 1284 φ (int): azimuthal angle 1285 Returns: 1286 numpy.ndarray: wavefunction angular component 1287 """ 1288 1289 l = self.l 1290 m = self.m 1291 1292 legendre = sp.lpmv(m, l, np.cos(θ)) 1293 1294 constant_factor = ((-1)**m) * np.sqrt( 1295 ((2*l+1) * sp.factorial(l-np.abs(m))) / 1296 (4*np.pi*sp.factorial(l+np.abs(m))) 1297 ) 1298 return constant_factor * legendre * np.real(np.exp(1.j*m*φ)) 1299 1300 def compute_wavefunction(self, box): 1301 """ Compute the normalized wavefunction as a product 1302 of its radial and angular components. 1303 Args: 1304 interval (int, int, int): (start, end, number of points) 1305 Returns: 1306 numpy.ndarray: wavefunction 1307 """ 1308 1309 if not isinstance(box, BoxOfPoints): 1310 raise ValueError('Arguments should be type BoxOfPoints') 1311 1312 # x-y grid to represent electron spatial distribution 1313 x, y, z = box.meshgrid 1314 1315 # Use epsilon to avoid division by zero during angle calculations 1316 eps = np.finfo(float).eps 1317 1318 # Ψnlm(r,θ,φ) = Rnl(r).Ylm(θ,φ) 1319 r = np.sqrt((x**2 + y**2 + z**2)) 1320 θ = np.arccos(z / (r + eps)) 1321 φ = np.sign(y) * np.arccos(x / np.sqrt(x**2 + y**2 + eps)) 1322 self.Ψ = self.radial_function(r) * self.angular_function(θ,φ) 1323 return self.Ψ 1324 1325 def compute_probability_density(self): 1326 """ Compute the probability density of a given wavefunction. 1327 Args: 1328 psi (numpy.ndarray): wavefunction 1329 Returns: 1330 numpy.ndarray: wavefunction probability density 1331 """ 1332 self.density = np.abs(self.Ψ)**2 1333 return self.density 1334 1335 def plot_probability_density_cross_section(self, box:BoxOfPoints, crop=None, colormap='rocket'): 1336 """ Plot the probability density of the hydrogen 1337 atom's wavefunction for a given quantum state (n,l,m). 1338 Args: 1339 a0 (float): scale factor 1340 dark_theme (bool): If True, uses a dark background for the plot, defaults to False 1341 colormap (str): Seaborn plot colormap, defaults to 'rocket' 1342 """ 1343 1344 _labels = ['x','y','z'] 1345 labels = [] 1346 axes = [] 1347 for label, axis in zip(_labels, box.axes): 1348 if axis.number_of_points!=1: 1349 labels.append(label) 1350 axes.append(axis) 1351 1352 if len(axes)>2: 1353 raise ValueError("ERROR: BoxOfPoints must have one axis with only 1 point") 1354 elif len(axes)<2: 1355 raise ValueError("ERROR: BoxOfPoints has insufficiently many axes") 1356 1357 xlabel = '$' + labels[0] + '/' + axes[0].unit + '$' 1358 ylabel = '$' + labels[1] + '/' + axes[1].unit + '$' 1359 1360 n = self.n 1361 l = self.l 1362 m = self.m 1363 1364 # Colormap validation 1365 try: 1366 sns.color_palette(colormap) 1367 except ValueError: 1368 raise ValueError(f'{colormap} is not a recognized Seaborn colormap.') 1369 1370 # Configure plot aesthetics using matplotlib rcParams settings 1371 plt.rcParams['font.family'] = 'STIXGeneral' 1372 plt.rcParams['mathtext.fontset'] = 'stix' 1373 plt.rcParams['xtick.major.width'] = 4 1374 plt.rcParams['ytick.major.width'] = 4 1375 plt.rcParams['xtick.major.size'] = 15 1376 plt.rcParams['ytick.major.size'] = 15 1377 plt.rcParams['xtick.labelsize'] = 30 1378 plt.rcParams['ytick.labelsize'] = 30 1379 plt.rcParams['axes.linewidth'] = 4 1380 1381 fig, ax = plt.subplots(figsize=(16, 16.5)) 1382 plt.subplots_adjust(top=0.82) 1383 plt.subplots_adjust(right=0.905) 1384 plt.subplots_adjust(left=-0.1) 1385 1386 # Compute and visualize the wavefunction probability density 1387 if not hasattr(self, 'Ψ'): 1388 print('Ψ not computed. Computing...') 1389 self.compute_wavefunction(box) 1390 self.compute_probability_density() 1391 print('Ψ computed.') 1392 1393 prob_density = self.density.squeeze() 1394 1395 [x0,x1,y0,y1]=[axes[0].lower_bound, axes[0].upper_bound,axes[1].lower_bound, axes[1].upper_bound] 1396 1397 if crop!=None: 1398 idx,jdx = axes[0].get_cropped_index(*crop[0]) 1399 prob_density = prob_density[idx:jdx,:] 1400 1401 idy,jdy = axes[1].get_cropped_index(*crop[1]) 1402 prob_density = prob_density[:,idy:jdy] 1403 1404 xx = axes[0].get_cropped_values(*crop[0]) 1405 x0, x1 = np.min(xx), np.max(xx) 1406 yy = axes[1].get_cropped_values(*crop[1]) 1407 y0, y1 = np.min(yy), np.max(yy) 1408 1409 vmax = np.max(prob_density) 1410 # Here we transpose the array to align the calculated z-x plane with Matplotlib's y-x imshow display 1411 im = ax.imshow(prob_density, cmap=sns.color_palette(colormap, as_cmap=True), 1412 extent=[x0,x1,y0,y1], vmin=0, vmax=vmax) 1413 1414 cbar = plt.colorbar(im, fraction=0.046, pad=0.03) 1415 cbar.set_ticks([]) 1416 1417 plt.rcParams['text.color'] = '#000000' 1418 title_color = '#000000' 1419 ax.tick_params(axis='x', colors='#000000') 1420 ax.tick_params(axis='y', colors='#000000') 1421 ax.set_xlabel(xlabel, fontsize=36) 1422 ax.set_ylabel(ylabel, fontsize=36) 1423 1424 ax.set_title('Hydrogen Atom - Orbital Wavefunction', 1425 pad=130, fontsize=44, loc='left', color=title_color) 1426 ax.text(x0, 1.12*y0, ( 1427 r'$|ψ_{n \ell m}(r, θ, φ)|^{2} =' 1428 r' |R_{n\ell}(r) Y_{\ell}^{m}(θ, ψ)|^2$' 1429 ), fontsize=36) 1430 ax.text(0.95*x0, 0.9*y0, r'$({0}, {1}, {2})$'.format(n, l, m), color='#dfdfdf', fontsize=42) 1431 ax.text(1.24*x1, 0.322*y1, 'Probability density', rotation='vertical', fontsize=40) 1432 ax.text(1.08*x1, 1.08*y0, f'{vmax:.3}', fontsize=24) 1433 ax.text(1.12*x1, 1.09*y1, 0, fontsize=24) 1434 ax.text(1.24*x1, 0.82*y1, '−', fontsize=34) 1435 ax.text(1.24*x1, 0.78*y0, '+', fontsize=34, rotation='vertical') 1436 ax.invert_yaxis() 1437 1438 # Save and display the plot 1439 plt.savefig(f'({n},{l},{m})-{labels[0]}{labels[1]}.png', bbox_inches='tight')
Orbital representation
Args
n (int): principal quantum number, determines the energy level and size of the orbital l (int): azimuthal quantum number, defines the shape of the orbital m (int): magnetic quantum number, defines the orientation of the orbital
Attributes
Methods
radial_function(self, r) angular_function(self, θ, φ) compute_wavefunction_cross_section(self, horizontal_interval, vertical_interval) compute_probability_density(self) plot_probability_density_cross_section(self, a0, dark_theme=False, colormap='rocket') TODO add 3D density plot
Examples
>>> orbital = Orbital(3,2,0)
>>> print(orbital)
Orbital(3d_{z^2})
Note
Integrates over the unspecified axis
Reference
Orbital(n: int = 0, l: int = 0, m: int = 0)
1164 def __init__(self, n:int=0, l:int=0, m:int=0): 1165 """Inherit list class method and add parent attribute to items""" 1166 1167 # Quantum numbers validation 1168 if not isinstance(n, int) or n < 1: 1169 raise ValueError('n should be an integer satisfying the condition: n >= 1') 1170 if not isinstance(l, int) or not (0 <= l < n): 1171 raise ValueError('l should be an integer satisfying the condition: 0 <= l < n') 1172 if not isinstance(m, int) or not (-l <= m <= l): 1173 raise ValueError('m should be an integer satisfying the condition: -l <= m <= l') 1174 1175 self.__dict__["quantum_numbers"] = [n,l,m] 1176 self.__dict__["n"] = n 1177 self.__dict__["l"] = l 1178 self.__dict__["m"] = m 1179 1180 super().__init__() 1181 self.__dict__["parent"] = None 1182 self.__dict__["_idx"] = None
Inherit list class method and add parent attribute to items
def
radial_function(self, r):
1238 def radial_function(self, r): 1239 """Compute the normalized radial part of the wavefunction using 1240 Laguerre polynomials and an exponential decay factor. 1241 Args: 1242 r (numpy.ndarray): radial coordinate 1243 Returns: 1244 numpy.ndarray: wavefunction radial component 1245 """ 1246 n = self.n 1247 l = self.l 1248 laguerre = sp.genlaguerre(n-l-1, 2*l+1) 1249 p = 2*r/n 1250 1251 constant_factor = np.sqrt( 1252 ((2/n)**3 * (sp.factorial(n-l-1))) / (2*n*(sp.factorial(n+l))) 1253 ) 1254 return constant_factor * np.exp(-p/2) * (p**l) * laguerre(p)
Compute the normalized radial part of the wavefunction using Laguerre polynomials and an exponential decay factor. Args: r (numpy.ndarray): radial coordinate Returns: numpy.ndarray: wavefunction radial component
def
compute_radial_density(self, rr):
1256 def compute_radial_density(self, rr): 1257 """ 1258 Compute the radial probability density for the 1259 list of radial points rr. 1260 """ 1261 area = 4*np.pi*np.square(rr) 1262 radius = np.square(self.radial_function(rr)) 1263 self.radial_density = radius * area 1264 return self.radial_density
Compute the radial probability density for the list of radial points rr.
def
plot_radial_density(self, ax, rr):
1266 def plot_radial_density(self, ax, rr): 1267 """ 1268 Compute the radial probability density for the 1269 list of radial points rr. 1270 """ 1271 if not hasattr(self, 'radial_density'): 1272 print('Ψ not computed. Computing...') 1273 self.compute_radial_density(rr) 1274 print('Ψ computed.') 1275 1276 ax.plot(rr, self.radial_density, label=self.__name__) 1277 return ax
Compute the radial probability density for the list of radial points rr.
def
angular_function(self, θ, φ):
1279 def angular_function(self, θ, φ): 1280 """ Compute the normalized angular part of the wavefunction using 1281 Legendre polynomials and a phase-shifting exponential factor. 1282 Args: 1283 θ (numpy.ndarray): polar angle 1284 φ (int): azimuthal angle 1285 Returns: 1286 numpy.ndarray: wavefunction angular component 1287 """ 1288 1289 l = self.l 1290 m = self.m 1291 1292 legendre = sp.lpmv(m, l, np.cos(θ)) 1293 1294 constant_factor = ((-1)**m) * np.sqrt( 1295 ((2*l+1) * sp.factorial(l-np.abs(m))) / 1296 (4*np.pi*sp.factorial(l+np.abs(m))) 1297 ) 1298 return constant_factor * legendre * np.real(np.exp(1.j*m*φ))
Compute the normalized angular part of the wavefunction using Legendre polynomials and a phase-shifting exponential factor. Args: θ (numpy.ndarray): polar angle φ (int): azimuthal angle Returns: numpy.ndarray: wavefunction angular component
def
compute_wavefunction(self, box):
1300 def compute_wavefunction(self, box): 1301 """ Compute the normalized wavefunction as a product 1302 of its radial and angular components. 1303 Args: 1304 interval (int, int, int): (start, end, number of points) 1305 Returns: 1306 numpy.ndarray: wavefunction 1307 """ 1308 1309 if not isinstance(box, BoxOfPoints): 1310 raise ValueError('Arguments should be type BoxOfPoints') 1311 1312 # x-y grid to represent electron spatial distribution 1313 x, y, z = box.meshgrid 1314 1315 # Use epsilon to avoid division by zero during angle calculations 1316 eps = np.finfo(float).eps 1317 1318 # Ψnlm(r,θ,φ) = Rnl(r).Ylm(θ,φ) 1319 r = np.sqrt((x**2 + y**2 + z**2)) 1320 θ = np.arccos(z / (r + eps)) 1321 φ = np.sign(y) * np.arccos(x / np.sqrt(x**2 + y**2 + eps)) 1322 self.Ψ = self.radial_function(r) * self.angular_function(θ,φ) 1323 return self.Ψ
Compute the normalized wavefunction as a product of its radial and angular components. Args: interval (int, int, int): (start, end, number of points) Returns: numpy.ndarray: wavefunction
def
compute_probability_density(self):
1325 def compute_probability_density(self): 1326 """ Compute the probability density of a given wavefunction. 1327 Args: 1328 psi (numpy.ndarray): wavefunction 1329 Returns: 1330 numpy.ndarray: wavefunction probability density 1331 """ 1332 self.density = np.abs(self.Ψ)**2 1333 return self.density
Compute the probability density of a given wavefunction. Args: psi (numpy.ndarray): wavefunction Returns: numpy.ndarray: wavefunction probability density
def
plot_probability_density_cross_section(self, box: crystal.types.BoxOfPoints, crop=None, colormap='rocket'):
1335 def plot_probability_density_cross_section(self, box:BoxOfPoints, crop=None, colormap='rocket'): 1336 """ Plot the probability density of the hydrogen 1337 atom's wavefunction for a given quantum state (n,l,m). 1338 Args: 1339 a0 (float): scale factor 1340 dark_theme (bool): If True, uses a dark background for the plot, defaults to False 1341 colormap (str): Seaborn plot colormap, defaults to 'rocket' 1342 """ 1343 1344 _labels = ['x','y','z'] 1345 labels = [] 1346 axes = [] 1347 for label, axis in zip(_labels, box.axes): 1348 if axis.number_of_points!=1: 1349 labels.append(label) 1350 axes.append(axis) 1351 1352 if len(axes)>2: 1353 raise ValueError("ERROR: BoxOfPoints must have one axis with only 1 point") 1354 elif len(axes)<2: 1355 raise ValueError("ERROR: BoxOfPoints has insufficiently many axes") 1356 1357 xlabel = '$' + labels[0] + '/' + axes[0].unit + '$' 1358 ylabel = '$' + labels[1] + '/' + axes[1].unit + '$' 1359 1360 n = self.n 1361 l = self.l 1362 m = self.m 1363 1364 # Colormap validation 1365 try: 1366 sns.color_palette(colormap) 1367 except ValueError: 1368 raise ValueError(f'{colormap} is not a recognized Seaborn colormap.') 1369 1370 # Configure plot aesthetics using matplotlib rcParams settings 1371 plt.rcParams['font.family'] = 'STIXGeneral' 1372 plt.rcParams['mathtext.fontset'] = 'stix' 1373 plt.rcParams['xtick.major.width'] = 4 1374 plt.rcParams['ytick.major.width'] = 4 1375 plt.rcParams['xtick.major.size'] = 15 1376 plt.rcParams['ytick.major.size'] = 15 1377 plt.rcParams['xtick.labelsize'] = 30 1378 plt.rcParams['ytick.labelsize'] = 30 1379 plt.rcParams['axes.linewidth'] = 4 1380 1381 fig, ax = plt.subplots(figsize=(16, 16.5)) 1382 plt.subplots_adjust(top=0.82) 1383 plt.subplots_adjust(right=0.905) 1384 plt.subplots_adjust(left=-0.1) 1385 1386 # Compute and visualize the wavefunction probability density 1387 if not hasattr(self, 'Ψ'): 1388 print('Ψ not computed. Computing...') 1389 self.compute_wavefunction(box) 1390 self.compute_probability_density() 1391 print('Ψ computed.') 1392 1393 prob_density = self.density.squeeze() 1394 1395 [x0,x1,y0,y1]=[axes[0].lower_bound, axes[0].upper_bound,axes[1].lower_bound, axes[1].upper_bound] 1396 1397 if crop!=None: 1398 idx,jdx = axes[0].get_cropped_index(*crop[0]) 1399 prob_density = prob_density[idx:jdx,:] 1400 1401 idy,jdy = axes[1].get_cropped_index(*crop[1]) 1402 prob_density = prob_density[:,idy:jdy] 1403 1404 xx = axes[0].get_cropped_values(*crop[0]) 1405 x0, x1 = np.min(xx), np.max(xx) 1406 yy = axes[1].get_cropped_values(*crop[1]) 1407 y0, y1 = np.min(yy), np.max(yy) 1408 1409 vmax = np.max(prob_density) 1410 # Here we transpose the array to align the calculated z-x plane with Matplotlib's y-x imshow display 1411 im = ax.imshow(prob_density, cmap=sns.color_palette(colormap, as_cmap=True), 1412 extent=[x0,x1,y0,y1], vmin=0, vmax=vmax) 1413 1414 cbar = plt.colorbar(im, fraction=0.046, pad=0.03) 1415 cbar.set_ticks([]) 1416 1417 plt.rcParams['text.color'] = '#000000' 1418 title_color = '#000000' 1419 ax.tick_params(axis='x', colors='#000000') 1420 ax.tick_params(axis='y', colors='#000000') 1421 ax.set_xlabel(xlabel, fontsize=36) 1422 ax.set_ylabel(ylabel, fontsize=36) 1423 1424 ax.set_title('Hydrogen Atom - Orbital Wavefunction', 1425 pad=130, fontsize=44, loc='left', color=title_color) 1426 ax.text(x0, 1.12*y0, ( 1427 r'$|ψ_{n \ell m}(r, θ, φ)|^{2} =' 1428 r' |R_{n\ell}(r) Y_{\ell}^{m}(θ, ψ)|^2$' 1429 ), fontsize=36) 1430 ax.text(0.95*x0, 0.9*y0, r'$({0}, {1}, {2})$'.format(n, l, m), color='#dfdfdf', fontsize=42) 1431 ax.text(1.24*x1, 0.322*y1, 'Probability density', rotation='vertical', fontsize=40) 1432 ax.text(1.08*x1, 1.08*y0, f'{vmax:.3}', fontsize=24) 1433 ax.text(1.12*x1, 1.09*y1, 0, fontsize=24) 1434 ax.text(1.24*x1, 0.82*y1, '−', fontsize=34) 1435 ax.text(1.24*x1, 0.78*y0, '+', fontsize=34, rotation='vertical') 1436 ax.invert_yaxis() 1437 1438 # Save and display the plot 1439 plt.savefig(f'({n},{l},{m})-{labels[0]}{labels[1]}.png', bbox_inches='tight')
Plot the probability density of the hydrogen atom's wavefunction for a given quantum state (n,l,m). Args: a0 (float): scale factor dark_theme (bool): If True, uses a dark background for the plot, defaults to False colormap (str): Seaborn plot colormap, defaults to 'rocket'
Inherited Members
1441class Atom(BdGPatch, MatrixTree): 1442 """ 1443 Atomic representation in the real-space basis 1444 1445 Attributes 1446 ---------- 1447 1448 Methods 1449 ------- 1450 1451 Examples 1452 -------- 1453 >>> hydrogen = Atom([0.2,0.4], [Orbital(1,0,0)]) 1454 >>> print(hydrogen) 1455 Atom(hydrogen) 1456 1457 Note 1458 ---- 1459 """ 1460 def __init__(self, fractional_coordinates:list=[0,0,0], orbitals=[]): 1461 """Inherit list class method and add parent attribute to items""" 1462 self.__dict__["fractional_coordinates"] = fractional_coordinates 1463 self._name = varname() 1464 super().__init__(*orbitals) 1465 self.__dict__["parent"] = None 1466 self.__dict__["_idx"] = None 1467 for item in orbitals: 1468 item.__dict__["parent"] = self 1469 1470 @property 1471 def __name__(self): 1472 return self._name 1473 1474 def __str__(self): 1475 """Returns object class name""" 1476 return f'{self.__class__.__name__}({self.__name__})' 1477 1478 def __repr__(self): 1479 """Returns object class name""" 1480 return f'{self.__class__.__name__}({self.__name__})' 1481 1482 @property 1483 def orbitals(self): 1484 return _Mapper(self)
Atomic representation in the real-space basis
Attributes
Methods
Examples
>>> hydrogen = Atom([0.2,0.4], [Orbital(1,0,0)])
>>> print(hydrogen)
Atom(hydrogen)
Note
Atom(fractional_coordinates: list = [0, 0, 0], orbitals=[])
1460 def __init__(self, fractional_coordinates:list=[0,0,0], orbitals=[]): 1461 """Inherit list class method and add parent attribute to items""" 1462 self.__dict__["fractional_coordinates"] = fractional_coordinates 1463 self._name = varname() 1464 super().__init__(*orbitals) 1465 self.__dict__["parent"] = None 1466 self.__dict__["_idx"] = None 1467 for item in orbitals: 1468 item.__dict__["parent"] = self
Inherit list class method and add parent attribute to items
Inherited Members
1486class Molecule(BdGPatch, MatrixTree): 1487 """ 1488 Molecule in the real-space basis 1489 1490 Attributes 1491 ---------- 1492 1493 Methods 1494 ------- 1495 1496 Examples 1497 -------- 1498 >>> hydrogen_A = Atom(fractional_coordinates=[0.25, 0.33], orbitals=[Orbital(1, 0, 0)]) 1499 >>> hydrogen_B = Atom(fractional_coordinates=[0.75, 0.33], orbitals=[Orbital(1, 0, 0)]) 1500 >>> oxygen = Atom(fractional_coordinates=[0.5,0.66], orbitals=[Orbital(2, 1, 0), Orbital(2, 1, 1), Orbital(2, 1, -1)]) 1501 >>> water = Molecule(fractional_coordinates=[0.5,0.5], atoms=[hydrogen_A, hydrogen_B, oxygen]) 1502 >>> water.n_spins = 2 1503 >>> water['oxygen',:,:] 1504 [Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑)] 1505 1506 >>> oxygen['2p_z'].n_spins = 1 1507 >>> water['oxygen',:,:] 1508 [Spin(↓), Spin(↑), Spin(↓), Spin(↑)] 1509 1510 >>> water[:,:,:] 1511 [Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑)] 1512 1513 Note 1514 ---- 1515 """ 1516 def __init__(self, fractional_coordinates:list=[0,0,0], atoms=[]): 1517 """Inherit list class method and add parent attribute to items""" 1518 self.fractional_coordinates = fractional_coordinates 1519 self._name = varname() 1520 # name = f'{fractional_coordinates}' 1521 super().__init__(*atoms) 1522 self.__dict__["parent"] = None 1523 self.__dict__["_idx"] = None 1524 for item in atoms: 1525 item.__dict__["parent"] = self 1526 1527 @property 1528 def __name__(self): 1529 return self._name 1530 1531 def __str__(self): 1532 """Returns object class name""" 1533 return f'{self.__class__.__name__}({self.__name__})' 1534 1535 def __repr__(self): 1536 """Returns object class name""" 1537 return f'{self.__class__.__name__}({self.__name__})' 1538 1539 @property 1540 def atoms(self): 1541 return _Mapper(self)
Molecule in the real-space basis
Attributes
Methods
Examples
>>> hydrogen_A = Atom(fractional_coordinates=[0.25, 0.33], orbitals=[Orbital(1, 0, 0)])
>>> hydrogen_B = Atom(fractional_coordinates=[0.75, 0.33], orbitals=[Orbital(1, 0, 0)])
>>> oxygen = Atom(fractional_coordinates=[0.5,0.66], orbitals=[Orbital(2, 1, 0), Orbital(2, 1, 1), Orbital(2, 1, -1)])
>>> water = Molecule(fractional_coordinates=[0.5,0.5], atoms=[hydrogen_A, hydrogen_B, oxygen])
>>> water.n_spins = 2
>>> water['oxygen',:,:]
[Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑)]
>>> oxygen['2p_z'].n_spins = 1
>>> water['oxygen',:,:]
[Spin(↓), Spin(↑), Spin(↓), Spin(↑)]
>>> water[:,:,:]
[Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑), Spin(↓), Spin(↑)]
Note
Molecule(fractional_coordinates: list = [0, 0, 0], atoms=[])
1516 def __init__(self, fractional_coordinates:list=[0,0,0], atoms=[]): 1517 """Inherit list class method and add parent attribute to items""" 1518 self.fractional_coordinates = fractional_coordinates 1519 self._name = varname() 1520 # name = f'{fractional_coordinates}' 1521 super().__init__(*atoms) 1522 self.__dict__["parent"] = None 1523 self.__dict__["_idx"] = None 1524 for item in atoms: 1525 item.__dict__["parent"] = self
Inherit list class method and add parent attribute to items
Inherited Members
1557class Cell(BdGPatch, MatrixTree): 1558 """ 1559 Cell in the real-space basis 1560 1561 Examples 1562 -------- 1563 1564 >>> A = Atom(fractional_coordinates=[0, 0]) 1565 >>> B = Atom(fractional_coordinates=[0, 0.5]) 1566 >>> cell = Cell(fractional_coordinates=[0,0], sites=[A, B]) 1567 >>> cell 1568 Cell([0, 0]) 1569 1570 """ 1571 def __init__(self, fractional_coordinates, sites, basis = None): 1572 """Inherit list class method and add parent attribute to items""" 1573 self.fractional_coordinates = list(fractional_coordinates) 1574 if type(basis)!=type(None): 1575 self.__dict__["cartesian_coordinates"] = list(np.matmul(fractional_coordinates, basis)) 1576 super().__init__(*sites) 1577 for site in sites: 1578 site.__dict__["parent"] = self 1579 self.__dict__[site.__name__] = site 1580 1581 @property 1582 def __name__(self): 1583 return str(self.fractional_coordinates) 1584 1585 def __repr__(self): 1586 """Returns object class name""" 1587 return f'{self.__class__.__name__}({self.__name__})' 1588 1589 def __str__(self): 1590 """Returns object class name""" 1591 return f'{self.__class__.__name__}({self.__name__})'
Cell in the real-space basis
Examples
>>> A = Atom(fractional_coordinates=[0, 0])
>>> B = Atom(fractional_coordinates=[0, 0.5])
>>> cell = Cell(fractional_coordinates=[0,0], sites=[A, B])
>>> cell
Cell([0, 0])
Cell(fractional_coordinates, sites, basis=None)
1571 def __init__(self, fractional_coordinates, sites, basis = None): 1572 """Inherit list class method and add parent attribute to items""" 1573 self.fractional_coordinates = list(fractional_coordinates) 1574 if type(basis)!=type(None): 1575 self.__dict__["cartesian_coordinates"] = list(np.matmul(fractional_coordinates, basis)) 1576 super().__init__(*sites) 1577 for site in sites: 1578 site.__dict__["parent"] = self 1579 self.__dict__[site.__name__] = site
Inherit list class method and add parent attribute to items
Inherited Members
1593class Lattice(BdGPatch, MatrixTree): 1594 """ 1595 Lattice 1596 1597 Examples 1598 -------- 1599 >>> A = Atom(fractional_coordinates=[0, 0]) 1600 >>> B = Atom(fractional_coordinates=[0.5, 0.5]) 1601 >>> my_lattice = Lattice(basis=[[1,0.5],[0,1]], cell=[A, B]) 1602 >>> sublattice_half = my_lattice.cut(number_of_cells=3, direction=0) 1603 >>> sublattice_full = my_lattice.cut(number_of_cells=2, direction=1) 1604 1605 Add lattice parameter 1606 >>> my_lattice.t = 1 1607 >>> my_lattice[0,0].t 1608 1 1609 1610 Add lattice parameter to particular cell 1611 >>> my_lattice[1,0].t = 2 1612 >>> my_lattice[0].t 1613 [1, 1] 1614 >>> my_lattice.t 1615 1 1616 >>> my_lattice[1,0].A 1617 Atom(A) 1618 >>> my_lattice[0].A 1619 [Atom(A), Atom(A)] 1620 >>> my_lattice.A 1621 Atom(A) 1622 1623 # Cutout a sub-region 1624 # >>> square = my_lattice.cutout(boundary=[[-1, -1], [-1, 1], [1, 1], [1, -1]]) 1625 # >>> my_lattice.boundaries["square"] 1626 # {'boundary': [[-1, -1], [-1, 1], [1, 1], [1, -1]], 'cells': [Cell([0, 0]), Cell([0, 1]), Cell([1, 0]), Cell([1, 1])]} 1627 1628 Note 1629 ---- 1630 """ 1631 def __init__(self, basis:list=[[1,0,0],[0,1,0],[0,0,1]], cell=[]): 1632 """Inherit list class method and add parent attribute to items""" 1633 self.__dict__["_name"] = varname() 1634 self.__dict__["basis"] = np.array(basis) 1635 self.__dict__["inv_basis"] = np.linalg.inv(basis) 1636 self.__dict__["cell"] = cell 1637 self.__dict__["nd"] = len(self.basis) 1638 self.__dict__["n_cells"] = [1 for _ in range(self.nd)] 1639 self.__dict__["pbc"] = [True for _ in range(self.nd)] 1640 self.__dict__["boundaries"] = {} 1641 self.__dict__["subregions"] = {} 1642 self.__dict__["n_spins"] = 1 1643 self._initialise_cells() 1644 super().__init__() 1645 1646 def _initialise_cells(self): 1647 """ 1648 Initialises self.cells as an empty array --called by self.__init__() 1649 """ 1650 coord = [0 for _ in range(self.nd)] 1651 sites = deepcopy(self.cell) 1652 cell = Cell(fractional_coordinates = coord, sites = sites, basis = self.basis) 1653 for parent_site in self.cell: 1654 parent_site.clear() 1655 for child_site in cell: 1656 name = parent_site.__name__ 1657 if name == child_site.__name__: 1658 self.__dict__[name] = parent_site 1659 parent_site.extend(child_site) 1660 child_site.__dict__["parent"] = parent_site 1661 array = np.empty(self.n_cells, dtype = Cell) 1662 array[*coord] = cell 1663 self.__dict__["cells"] = array 1664 1665 # def __getattr__(self, name: str): 1666 # try: 1667 # return self.__dict__[name] 1668 # except: 1669 # return _Mapper([getattr(item, name) for item in _Mapper(self.cells.tolist())]) 1670 1671 # def __setattr__(self, name, value, update = True): 1672 # """ 1673 # Set attributes down the lineage, and add attribute to list of 1674 # attributes to be updates 1675 # """ 1676 # if update: 1677 # self.__cache__(name, value) 1678 1679 # self.__dict__[f"{name}"] = value 1680 1681 # [item.__setattr__(name, value, update = False) for item in self] 1682 1683 # def __setattr__(self, name: str, value): 1684 # if name in self.__dict__.keys(): 1685 # self.__dict__[name] = value 1686 # else: 1687 # [setattr(item, name, value) for item in _Mapper(self.cells.tolist())] 1688 1689 @property 1690 def __name__(self): 1691 return self._name 1692 1693 def __str__(self): 1694 """Returns object class name""" 1695 return f'{self.__name__}({[list(x) for x in self.basis]})' 1696 1697 def __repr__(self): 1698 """Returns object class name""" 1699 return f'{self.__class__.__name__}({[list(x) for x in self.basis]})' 1700 1701 def cut(self, number_of_cells: int, direction: int, glue_edges: bool = True): 1702 """ 1703 Returns a deepcopy of self 1704 """ 1705 1706 if not glue_edges: 1707 raise ValueError("TODO: Open boundary conditions not yet implemented!") 1708 1709 _slice = [(slice(0, self.n_cells[i])) for i in range(self.nd)] 1710 self.n_cells[direction] = number_of_cells 1711 cells = np.empty(self.n_cells, dtype=Cell) 1712 cells[*_slice] = self.cells 1713 self.cells = cells 1714 indices = np.moveaxis(np.indices(self.n_cells), 0, -1) 1715 indices = indices[cells == None] 1716 for index in indices: 1717 cell = Cell(fractional_coordinates=index, sites=deepcopy(self.cell), basis=self.basis) 1718 self.cells[*index] = cell 1719 1720 for parent_site in self.cell: 1721 parent_site.clear() 1722 for child_site in cell: 1723 name = parent_site.__name__ 1724 if name == child_site.__name__: 1725 self.__dict__[name] = parent_site 1726 parent_site.extend(child_site) 1727 child_site.__dict__["parent"] = parent_site 1728 cell = deepcopy(self.cell) 1729 cells = [Cell(fractional_coordinates=index, sites=[deepcopy(site) for site in self.cell], basis=self.basis) for index in indices] 1730 indices = np.array(indices).T 1731 self.cells[*indices] = cells 1732 1733 def mapper(parent_site, cell): 1734 parent_site.clear() 1735 for child_site in cell: 1736 name = parent_site.__name__ 1737 if name == child_site.__name__: 1738 self.__dict__[name] = parent_site 1739 parent_site.extend(child_site) 1740 child_site.__dict__["parent"] = parent_site 1741 1742 [[mapper(site, cell) for site in self.cell] for cell in cells] 1743 1744 self.pbc[direction] = glue_edges 1745 1746 self.init_cells() 1747 1748 try: 1749 name = varname() 1750 sublattice = Lattice(basis=copy(self.basis), cell=copy(self.cell)) 1751 sublattice._name = name 1752 sublattice.n_cells = copy(self.n_cells) 1753 sublattice.pbc = copy(self.pbc) 1754 sublattice.cells = deepcopy(self.cells) 1755 sublattice.extend(flatten(sublattice.cells)) 1756 return sublattice 1757 except: 1758 pass 1759 1760 @property 1761 def n_total_cells(self): 1762 return np.prod(self.n_cells) 1763 1764 def cutout(self, boundary, cut_from_infinite_plane = True): 1765 """ 1766 Cutout cells within the boundary. 1767 1768 The boundary is a list of points in the 1769 Cartesian. 1770 """ 1771 try: 1772 name = varname() 1773 except: 1774 name = len(self.boundaries) 1775 1776 if cut_from_infinite_plane: 1777 inv_boundary = np.matmul(boundary, self.inv_basis) 1778 polytope = Polytope(points=inv_boundary) 1779 inside, outside = polytope.cut() 1780 corners = polytope.frac_corners 1781 n_cells = corners[:,1] - corners[:,0] 1782 for d, n in enumerate(n_cells): 1783 self.cut(number_of_cells=int(n+1), direction=d) 1784 1785 polytope = Polytope(points=boundary) 1786 1787 # clip 1788 # ax = polytope.plot() 1789 # cartesian_coordinates = flatten(self.cartesian_coordinates, self.nd-1) 1790 # cartesian_coordinates = clean(cartesian_coordinates) 1791 # cartesian_coordinates = np.transpose(cartesian_coordinates) 1792 # # plt.scatter(*cartesian_coordinates) 1793 1794 # inside = np.matmul(inside, self.basis) 1795 # outside = np.matmul(outside, self.basis) 1796 # ax.scatter(*inside.T) 1797 # ax.scatter(*outside.T, c="r") 1798 # lim=[-5,20] 1799 # ax.set(xlim=lim, ylim=lim) 1800 # plt.show() 1801 # exit() 1802 1803 inside = inside.astype(int) 1804 outside = outside.astype(int) 1805 1806 cells = [cell.cartesian_coordinates for cell in self] 1807 cells = np.array(cells, dtype=float) 1808 1809 cartesian = np.matmul(inside, self.basis) 1810 for cartesian, cell in zip(cartesian, inside): 1811 self.cells[*cell].__dict__["fractional_coordinates"] = cell 1812 self.cells[*cell].__dict__["cartesian_coordinates"] = cartesian 1813 cells = [self.cells[*cell] for cell in inside] 1814 else: 1815 polytope = Polytope(points=boundary, basis=self.basis) 1816 cells = flatten_depth(self.cartesian_coordinates, self.nd-1) 1817 cells = [list(cell) for cell in cells] 1818 cells = np.array(cells, dtype=float) 1819 inside, outside = polytope.cut(cells) 1820 inside = np.matmul(inside, self.inv_basis).astype(int) 1821 outside = np.matmul(outside, self.inv_basis).astype(int) 1822 cells = [self.cells[*cell] for cell in inside] 1823 1824 for cell in outside: 1825 try: 1826 self.cells[*cell] = None 1827 except: 1828 pass 1829 1830 self.init_cells() 1831 1832 self.boundaries[name] = polytope 1833 self.subregions[name] = cells 1834 1835 return self.boundaries[name] 1836 1837 def init_cells(self): 1838 self.clear() 1839 self.extend(flatten(self.cells)) 1840 1841 def __set_indices__(self): 1842 """ 1843 """ 1844 self.initialised = True 1845 def f(l, i): l.__dict__["index"] = i 1846 [f(leaf, i) for i, leaf in enumerate(self.leaves)] 1847 for site in self.cell: 1848 _slice = [slice(None, None, None) for _ in range(self.nd+site.depth)] 1849 _slice.append(site.__name__) 1850 site.clear() 1851 site.extend(flatten(self[*_slice])) 1852 1853 def plt_boundary(self, ax=None): 1854 for boundary in self.boundaries: 1855 ax = boundary.plot(ax) 1856 return ax 1857 1858 def plt_lattice(self, ax=None): 1859 ax.scatter(*self.cartesian().T, c='b') 1860 return ax 1861 1862 def __getitem__(self,keys): 1863 cells = self.cells 1864 _keys = [] 1865 if type(keys)==int: 1866 return _Mapper(self.cells[keys]) 1867 if len(keys)>self.nd: 1868 _keys = keys[self.nd:] 1869 keys = keys[:self.nd] 1870 1871 cells = cells[keys] 1872 if type(cells)==np.ndarray: 1873 cells = _Mapper(cells.tolist()) 1874 1875 _slice = [] 1876 for key in keys: 1877 if key==slice(None,None,None): 1878 _slice.append(slice(None,None,None)) 1879 1880 if len(_keys)>0: 1881 cells = _Mapper(_Mapper(cells)[*_slice,*_keys]) 1882 1883 return cells 1884 1885 def plot(self, ax=None): 1886 if ax==None: 1887 ax = plt.axes() 1888 plt.tight_layout() 1889 cartesian_coordinates = flatten_depth(self.cartesian_coordinates, self.nd-1) 1890 cartesian_coordinates = clean(cartesian_coordinates) 1891 cartesian_coordinates = np.transpose(cartesian_coordinates) 1892 plt.scatter(*cartesian_coordinates)
Lattice
Examples
>>> A = Atom(fractional_coordinates=[0, 0])
>>> B = Atom(fractional_coordinates=[0.5, 0.5])
>>> my_lattice = Lattice(basis=[[1,0.5],[0,1]], cell=[A, B])
>>> sublattice_half = my_lattice.cut(number_of_cells=3, direction=0)
>>> sublattice_full = my_lattice.cut(number_of_cells=2, direction=1)
Add lattice parameter
>>> my_lattice.t = 1
>>> my_lattice[0,0].t
1
Add lattice parameter to particular cell
>>> my_lattice[1,0].t = 2
>>> my_lattice[0].t
[1, 1]
>>> my_lattice.t
1
>>> my_lattice[1,0].A
Atom(A)
>>> my_lattice[0].A
[Atom(A), Atom(A)]
>>> my_lattice.A
Atom(A)
Cutout a sub-region
>>> square = my_lattice.cutout(boundary=[[-1, -1], [-1, 1], [1, 1], [1, -1]])
>>> my_lattice.boundaries["square"]
{'boundary': [[-1, -1], [-1, 1], [1, 1], [1, -1]], 'cells': [Cell([0, 0]), Cell([0, 1]), Cell([1, 0]), Cell([1, 1])]}
Note
Lattice(basis: list = [[1, 0, 0], [0, 1, 0], [0, 0, 1]], cell=[])
1631 def __init__(self, basis:list=[[1,0,0],[0,1,0],[0,0,1]], cell=[]): 1632 """Inherit list class method and add parent attribute to items""" 1633 self.__dict__["_name"] = varname() 1634 self.__dict__["basis"] = np.array(basis) 1635 self.__dict__["inv_basis"] = np.linalg.inv(basis) 1636 self.__dict__["cell"] = cell 1637 self.__dict__["nd"] = len(self.basis) 1638 self.__dict__["n_cells"] = [1 for _ in range(self.nd)] 1639 self.__dict__["pbc"] = [True for _ in range(self.nd)] 1640 self.__dict__["boundaries"] = {} 1641 self.__dict__["subregions"] = {} 1642 self.__dict__["n_spins"] = 1 1643 self._initialise_cells() 1644 super().__init__()
Inherit list class method and add parent attribute to items
def
cut(self, number_of_cells: int, direction: int, glue_edges: bool = True):
1701 def cut(self, number_of_cells: int, direction: int, glue_edges: bool = True): 1702 """ 1703 Returns a deepcopy of self 1704 """ 1705 1706 if not glue_edges: 1707 raise ValueError("TODO: Open boundary conditions not yet implemented!") 1708 1709 _slice = [(slice(0, self.n_cells[i])) for i in range(self.nd)] 1710 self.n_cells[direction] = number_of_cells 1711 cells = np.empty(self.n_cells, dtype=Cell) 1712 cells[*_slice] = self.cells 1713 self.cells = cells 1714 indices = np.moveaxis(np.indices(self.n_cells), 0, -1) 1715 indices = indices[cells == None] 1716 for index in indices: 1717 cell = Cell(fractional_coordinates=index, sites=deepcopy(self.cell), basis=self.basis) 1718 self.cells[*index] = cell 1719 1720 for parent_site in self.cell: 1721 parent_site.clear() 1722 for child_site in cell: 1723 name = parent_site.__name__ 1724 if name == child_site.__name__: 1725 self.__dict__[name] = parent_site 1726 parent_site.extend(child_site) 1727 child_site.__dict__["parent"] = parent_site 1728 cell = deepcopy(self.cell) 1729 cells = [Cell(fractional_coordinates=index, sites=[deepcopy(site) for site in self.cell], basis=self.basis) for index in indices] 1730 indices = np.array(indices).T 1731 self.cells[*indices] = cells 1732 1733 def mapper(parent_site, cell): 1734 parent_site.clear() 1735 for child_site in cell: 1736 name = parent_site.__name__ 1737 if name == child_site.__name__: 1738 self.__dict__[name] = parent_site 1739 parent_site.extend(child_site) 1740 child_site.__dict__["parent"] = parent_site 1741 1742 [[mapper(site, cell) for site in self.cell] for cell in cells] 1743 1744 self.pbc[direction] = glue_edges 1745 1746 self.init_cells() 1747 1748 try: 1749 name = varname() 1750 sublattice = Lattice(basis=copy(self.basis), cell=copy(self.cell)) 1751 sublattice._name = name 1752 sublattice.n_cells = copy(self.n_cells) 1753 sublattice.pbc = copy(self.pbc) 1754 sublattice.cells = deepcopy(self.cells) 1755 sublattice.extend(flatten(sublattice.cells)) 1756 return sublattice 1757 except: 1758 pass
Returns a deepcopy of self
def
cutout(self, boundary, cut_from_infinite_plane=True):
1764 def cutout(self, boundary, cut_from_infinite_plane = True): 1765 """ 1766 Cutout cells within the boundary. 1767 1768 The boundary is a list of points in the 1769 Cartesian. 1770 """ 1771 try: 1772 name = varname() 1773 except: 1774 name = len(self.boundaries) 1775 1776 if cut_from_infinite_plane: 1777 inv_boundary = np.matmul(boundary, self.inv_basis) 1778 polytope = Polytope(points=inv_boundary) 1779 inside, outside = polytope.cut() 1780 corners = polytope.frac_corners 1781 n_cells = corners[:,1] - corners[:,0] 1782 for d, n in enumerate(n_cells): 1783 self.cut(number_of_cells=int(n+1), direction=d) 1784 1785 polytope = Polytope(points=boundary) 1786 1787 # clip 1788 # ax = polytope.plot() 1789 # cartesian_coordinates = flatten(self.cartesian_coordinates, self.nd-1) 1790 # cartesian_coordinates = clean(cartesian_coordinates) 1791 # cartesian_coordinates = np.transpose(cartesian_coordinates) 1792 # # plt.scatter(*cartesian_coordinates) 1793 1794 # inside = np.matmul(inside, self.basis) 1795 # outside = np.matmul(outside, self.basis) 1796 # ax.scatter(*inside.T) 1797 # ax.scatter(*outside.T, c="r") 1798 # lim=[-5,20] 1799 # ax.set(xlim=lim, ylim=lim) 1800 # plt.show() 1801 # exit() 1802 1803 inside = inside.astype(int) 1804 outside = outside.astype(int) 1805 1806 cells = [cell.cartesian_coordinates for cell in self] 1807 cells = np.array(cells, dtype=float) 1808 1809 cartesian = np.matmul(inside, self.basis) 1810 for cartesian, cell in zip(cartesian, inside): 1811 self.cells[*cell].__dict__["fractional_coordinates"] = cell 1812 self.cells[*cell].__dict__["cartesian_coordinates"] = cartesian 1813 cells = [self.cells[*cell] for cell in inside] 1814 else: 1815 polytope = Polytope(points=boundary, basis=self.basis) 1816 cells = flatten_depth(self.cartesian_coordinates, self.nd-1) 1817 cells = [list(cell) for cell in cells] 1818 cells = np.array(cells, dtype=float) 1819 inside, outside = polytope.cut(cells) 1820 inside = np.matmul(inside, self.inv_basis).astype(int) 1821 outside = np.matmul(outside, self.inv_basis).astype(int) 1822 cells = [self.cells[*cell] for cell in inside] 1823 1824 for cell in outside: 1825 try: 1826 self.cells[*cell] = None 1827 except: 1828 pass 1829 1830 self.init_cells() 1831 1832 self.boundaries[name] = polytope 1833 self.subregions[name] = cells 1834 1835 return self.boundaries[name]
Cutout cells within the boundary.
The boundary is a list of points in the Cartesian.
def
plot(self, ax=None):
1885 def plot(self, ax=None): 1886 if ax==None: 1887 ax = plt.axes() 1888 plt.tight_layout() 1889 cartesian_coordinates = flatten_depth(self.cartesian_coordinates, self.nd-1) 1890 cartesian_coordinates = clean(cartesian_coordinates) 1891 cartesian_coordinates = np.transpose(cartesian_coordinates) 1892 plt.scatter(*cartesian_coordinates)
Inherited Members
def
flatten_depth(items, depth: int):
def
flatten( l, ltypes=(<class 'list'>, <class 'tuple'>, <class 'tree.network._Mapper'>)):
1902def flatten(l, ltypes=(list, tuple, _Mapper)): 1903 if type(l)==np.ndarray: 1904 l=l.tolist() 1905 ltype = type(l) 1906 l = list(l) 1907 i = 0 1908 while i < len(l): 1909 while type(l[i]) in ltypes: 1910 if not l[i]: 1911 l.pop(i) 1912 i -= 1 1913 break 1914 else: 1915 l[i:i + 1] = l[i] 1916 i += 1 1917 return ltype(l)
def
clean(items):
def
friedel_oscillations():
1922def friedel_oscillations(): 1923 torch.set_default_dtype(torch.float32) 1924 1925 A = Atom(fractional_coordinates=[0, 0]) 1926 B = Atom(fractional_coordinates=[0.5, 0.5]) 1927 # A.hop(B, t=1) 1928 # A.hop(B, vector=[1,0], t=1) 1929 # B.hop(B, vector=[1,0], t=1) 1930 # A.hop(A, vector=[0,1], t=1) 1931 # B.hop(B, vector=[0,1], t=1) 1932 # A.add(μ=-3.95) 1933 lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A]) 1934 # lattice = Lattice(basis=[[1, 0],[0, 1]]) 1935 # square = lattice.cutout(boundary=[[-1, -1], [-1, 1], [1, 1], [1, -1]]) 1936 lattice.hop(vector=[1,0], t=1) 1937 lattice.hop(vector=[0,1], t=1) 1938 sublattice_half = lattice.cut(number_of_cells=11, direction=0) 1939 sublattice_half = lattice.cut(number_of_cells=11, direction=1) 1940 # lattice.n_spins = 2 1941 lattice.add(μ=3.95, structure=-1) 1942 # lattice.add(μ=-3.95, structure=-1, pauli_vector=[0,1,1,0]) 1943 # lattice.add(μ=-3.95, structure=1, pauli_vector=[0,0,0,1]) 1944 # lattice[1,0].add(V=-1., structure=-1, pauli_vector=[0,0,0,1]) 1945 lattice[5,5].add(V=1., structure=-1) 1946 # lattice[1,0].hop(lattice[0,0], t=1) 1947 # lattice.hop(vector=[1,0], t=1, pauli_vector=[0,0,0,1]) 1948 # lattice.hop(vector=[1,0], t=1, pauli_vector=[1,0,0,0]) 1949 # exit() 1950 # lattice.add(Δ=5j, anomalous=True, pauli_vector=[0,0,1,0]) 1951 # sublattice_half = my_lattice.cut(number_of_cells=3, direction=0) 1952 # sublattice_full = my_lattice.cut(number_of_cells=3, direction=1) 1953 # print(my_lattice.cells) 1954 # exit() 1955 # square = my_lattice.cutout(boundary=[[-10, -10], [-10, 10], [10, 10], [10, -10]]) 1956 # print(my_lattice[:,:,:,:]) 1957 1958 # print(my_lattice.indices) 1959 # exit() 1960 1961 # matrix = lattice.get_matrix() 1962 1963 1964 # exit() 1965 # matrix = np.real(matrix) 1966 # exit() 1967 # c = plt.imshow(matrix) 1968 # plt.colorbar(c) 1969 # plt.show() 1970 # exit() 1971 1972 w, v = lattice.solve() 1973 1974 # exit() 1975 # print(w.dtype) 1976 # print(v) 1977 # print(w) 1978 # print(np.shape(v)) 1979 # exit() 1980 # dm = my_lattice.density_matrix() 1981 # print(np.shape(dm)) 1982 # fermi = my_lattice.fermi_distribution(temperature=10) 1983 # fermi = my_lattice.thermal_density_matrix(temperature=0) 1984 # print(t1) 1985 # print(t2) 1986 # greens = my_lattice.greens_function(energy=0, resolution=0.1, sites=[A, B], sites_to=[B, A], hop=[1,0], spin_flip=True, anomalous = True) 1987 # greens = my_lattice.greens_function(energy=arange(-4,4,0.1), resolution=0.1, sites=[A, B], sites_to=[B, A], hop=[1,0], spin_flip=True, anomalous = True) 1988 # greens = my_lattice.greens_function(energy=0, resolution=0.1, sites=[A, B], sites_to=[B, A], hop=[1,0], spin_flip=True, anomalous = True, spin_polarised=True) 1989 greens = lattice.greens_function(energy=0, resolution=0.1) 1990 # energy = arange(-4,4,0.1) 1991 # greens = my_lattice.greens_function(energy=energy, resolution=0.1) 1992 # greens = sum(greens, axis=0) 1993 # plt.scatter(energy, greens) 1994 # plt.show() 1995 # exit() 1996 # greens = greens.reshape(11,11) 1997 c = plt.imshow(greens) 1998 plt.colorbar(c) 1999 plt.show() 2000 2001 exit() 2002 2003 # square = my_lattice.cutout(boundary=[[0, 0], [0, 10], [10, 10], [10, 0]]) 2004 2005 def plot(): 2006 ax = plt.axes() 2007 plt.tight_layout() 2008 my_lattice.plot(ax) 2009 # square.plot(ax) 2010 lim=[-1,3] 2011 ax.set(xlim=lim, ylim=lim) 2012 plt.show()
def
mean_field_theory():
2014def mean_field_theory(): 2015 2016 A = Atom(fractional_coordinates=[0, 0]) 2017 lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A]) 2018 lattice.n_spins = 2 2019 # lattice.hop(vector=[1,0], t=1, structure=-1) 2020 # lattice.hop(vector=[0,1], t=1, structure=-1) 2021 lattice.cut(number_of_cells=1, direction=0) 2022 lattice.cut(number_of_cells=1, direction=1) 2023 # lattice.add(μ=-1.7, structure=-1) 2024 # lattice.add(U=2.8, structure=-1, interaction=True) 2025 # lattice.add(φ=2.40, structure=-1, renormalise=True) 2026 lattice.add(Δ=1.00, anomalous=True, pauli_vector=[1,0,0,0], renormalise=True) 2027 2028 # lattice[5,5].add(V=1., structure=-1) 2029 2030 print(lattice._diagonal) 2031 print(lattice._offdiagonal) 2032 # print(lattice._hartree) 2033 # print(lattice._fock) 2034 print(lattice._gorkov) 2035 # print(lattice._interaction) 2036 exit() 2037 2038 print(lattice.get_matrix()) 2039 exit() 2040 2041 w, v = lattice.solve() 2042 exit() 2043 lattice.extract_fields(w, v) 2044 # lattice.self_consistent_step() 2045 2046 exit() 2047 2048 lattice.add(ρ=3.15, structure=-1, renormalise=True) 2049 lattice.add(φ=2.15, structure=-1, renormalise=True) 2050 lattice.add(Δ=3, anomalous=True, pauli_vector=[0,1j,0,0], renormalise=True) 2051 2052 matrix = lattice.matrix 2053 2054 print(matrix) 2055 2056 2057 exit() 2058 2059 # lattice.add(Δ=5j, anomalous=True, pauli_vector=[0,1,0,0]) 2060 2061 lattice.self_consistent_step() 2062 2063 # lattice.scmft(friction=friction,max_iterations=max_iterations,convergence_factor=convergence_factor) 2064 2065 2066 exit() 2067 2068 w, v = lattice.solve() 2069 2070 print(w)
def
test_indices():
2074def test_indices(): 2075 2076 A = Atom(fractional_coordinates=[0, 0]) 2077 B = Atom(fractional_coordinates=[0.5, 0]) 2078 lattice = Lattice(basis=[[1, 0],[0, 1]], cell=[A, B]) 2079 lattice.cut(number_of_cells=3, direction=0) 2080 lattice.cut(number_of_cells=3, direction=1) 2081 lattice.n_spins = 2 2082 lattice.hop(t=1, vector=[1,0], structure=-1) 2083 # lattice.hop(t=1, vector=[0,1], structure=-1) 2084 # lattice._initialise_cells 2085 2086 lattice.__set_indices__() 2087 2088 # >>> lattice.get_indices(site=A, anomalous=True, structure=-1, pauli_vector=[0,1,0,0]) 2089 # {(1+0j): [37, 45, 53, 61, 69, 40, 48, 56, 64]} 2090 exit() 2091 2092 print(lattice.get_diagonal_indices(site=A, anomalous=False, structure=-1, pauli_vector=[0,1,0,0])) 2093 exit() 2094 2095 print(lattice.get_offdiagonal_indices(site, neighbour, vector, trs, anomalous, structure, pauli_vector)) 2096 exit() 2097 2098 exit() 2099 2100 2101 # lattice.add(U=2.8, structure=-1, interaction=True) 2102 # lattice.add(φ=2.40, structure=-1, renormalise=True) 2103 # lattice.add(Δ=1.00, anomalous=True, pauli_vector=[1,0,0,0], renormalise=True) 2104 2105 # lattice[5,5].add(V=1., structure=-1)
lattice =
Lattice([[np.int64(1), np.int64(0)], [np.int64(0), np.int64(1)]])