From 14e892bbeace387bf2298d582417f29e266a3391 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Tue, 13 Sep 2022 14:32:46 +0200 Subject: [PATCH 01/32] Add support to `get_coefficient` for getting the radial dependence of the green's function --- masci_tools/tools/greensfunction.py | 66 ++++++++++++++++++++++++----- tests/tools/test_greensfunction.py | 10 +++++ 2 files changed, 66 insertions(+), 10 deletions(-) diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 9aeaf2f8b..59bd8f937 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -749,15 +749,15 @@ def get_coefficient(self, name: CoefficientName, spin: int | None = None, radial if name.startswith('ulo'): r = self.radial_functions['ulo'] elif name.startswith('u'): - r = self.radial_functions['u'][self.atomType - 1] + r = self.radial_functions['u'][self.atomType - 1, :, self.l] else: - r = self.radial_functions['d'][self.atomType - 1] + r = self.radial_functions['udot'][self.atomType - 1, :, self.l] if name.endswith('ulo'): rp = self.radial_functions['ulop'] elif name.endswith('u'): - rp = self.radial_functions['u'][self.atomTypep - 1] + rp = self.radial_functions['u'][self.atomTypep - 1, :, self.lp] else: - rp = self.radial_functions['d'][self.atomTypep - 1] + rp = self.radial_functions['udot'][self.atomTypep - 1, :, self.lp] if not self.onsite: #The radial functions are stored as r*u(r), meaning #when multiplied they produce the right factor for @@ -765,6 +765,23 @@ def get_coefficient(self, name: CoefficientName, spin: int | None = None, radial #independent so we need to multiply each radial function by r r *= self.radial_functions['rmsh'][self.atomType - 1] rp *= self.radial_functions['rmsh'][self.atomTypep - 1] + if name == 'uloulo': + radial_part = np.einsum('lxsr,pyst->lpxyrt', r, rp) + else: + radial_part = np.einsum('...xsr,...yst->...xyrt', r, rp) + elif name == 'uloulo': + radial_part = np.einsum('lxsr,pysr->lpxyr', r, rp) + else: + radial_part = np.einsum('...xsr,...ysr->...xyr', r, rp) + + if spin is not None: + spin1, spin2 = self.to_spin_indices(spin) + if self.onsite: + radial_part = radial_part[..., spin1, spin2, :] + else: + radial_part = radial_part[..., spin1, spin2, :, :] + radial_part = radial_part.T + else: if spin is not None: spin1, spin2 = self.to_spin_indices(spin) @@ -784,7 +801,18 @@ def get_coefficient(self, name: CoefficientName, spin: int | None = None, radial return np.swapaxes(data, -2, -1) if radial: - raise NotImplementedError() + if spin is not None: + if name == 'uloulo': + return np.einsum('...ij,r...ij->r...ij', data, radial_part) + if 'lo' in name: + return np.einsum('...i,r...i->r...i', data, radial_part) + return np.einsum('...,r...->r...', data, radial_part) + + if name == 'uloulo': + return np.einsum('...ijklm,r...ijkl->r...ijklm', data, radial_part) + if 'lo' in name: + return np.einsum('...ijkl,r...ijk->r...ijkl', data, radial_part) + return np.einsum('...ijl,r...ij->r...ijl', data, radial_part) if spin is not None: if name == 'uloulo': @@ -882,7 +910,8 @@ def energy_dependence(self, mp: int | None = None, spin: int | None = None, imag: bool = True, - both_contours: bool = False) -> np.ndarray: + both_contours: bool = False, + radial: bool = False) -> np.ndarray: r""" Select data with energy dependence @@ -892,6 +921,7 @@ def energy_dependence(self, :param both_contours: bool id True the data is not added for both energy contours :param imag: bool if True and both_contours is False the imaginary part:math:`\frac{1}{2i}\left[G\left(z\right)-G\left(z^\ast\right)\right]` is returned otherwise the real part :math:`\frac{1}{2}\left[G\left(z\right)+G\left(z^\ast\right)\right]` + :param radial: bool if True the green's function will be returned with radial dependence :returns: numpy array with the selected data """ @@ -908,12 +938,10 @@ def energy_dependence(self, else: mp_index = slice(self.lmax - self.l, self.lmax + self.lp + 1, 1) - kwargs: dict[str, Any] = { - 'spin': spin, - } + kwargs: dict[str, Any] = {'spin': spin, 'radial': radial} if self.sphavg: gf = self.get_coefficient('sphavg', **kwargs)[:, m_index, mp_index, ...] - else: + elif not radial: gf = self.get_coefficient('uu', **kwargs)[:,m_index,mp_index,...] \ + self.get_coefficient('ud', **kwargs)[:,m_index,mp_index,...] \ + self.get_coefficient('du', **kwargs)[:,m_index,mp_index,...] \ @@ -922,6 +950,24 @@ def energy_dependence(self, + np.sum(self.get_coefficient('ulou', **kwargs)[:,m_index,mp_index,...], axis=-1) \ + np.sum(self.get_coefficient('dulo', **kwargs)[:,m_index,mp_index,...], axis=-1) \ + np.sum(self.get_coefficient('uloulo', **kwargs)[:,m_index,mp_index,...], axis=(-1,-2)) + elif self.onsite: + gf = self.get_coefficient('uu', **kwargs)[:,:,m_index,mp_index,...] \ + + self.get_coefficient('ud', **kwargs)[:,:,m_index,mp_index,...] \ + + self.get_coefficient('du', **kwargs)[:,:,m_index,mp_index,...] \ + + self.get_coefficient('dd', **kwargs)[:,:,m_index,mp_index,...] \ + + np.sum(self.get_coefficient('uulo', **kwargs)[:,:,m_index,mp_index,...], axis=-1) \ + + np.sum(self.get_coefficient('ulou', **kwargs)[:,:,m_index,mp_index,...], axis=-1) \ + + np.sum(self.get_coefficient('dulo', **kwargs)[:,:,m_index,mp_index,...], axis=-1) \ + + np.sum(self.get_coefficient('uloulo', **kwargs)[:,:,m_index,mp_index,...], axis=(-1,-2)) + else: + gf = self.get_coefficient('uu', **kwargs)[:,:,:,m_index,mp_index,...] \ + + self.get_coefficient('ud', **kwargs)[:,:,:,m_index,mp_index,...] \ + + self.get_coefficient('du', **kwargs)[:,:,:,m_index,mp_index,...] \ + + self.get_coefficient('dd', **kwargs)[:,:,:,m_index,mp_index,...] \ + + np.sum(self.get_coefficient('uulo', **kwargs)[:,:,:,m_index,mp_index,...], axis=-1) \ + + np.sum(self.get_coefficient('ulou', **kwargs)[:,:,:,m_index,mp_index,...], axis=-1) \ + + np.sum(self.get_coefficient('dulo', **kwargs)[:,:,:,m_index,mp_index,...], axis=-1) \ + + np.sum(self.get_coefficient('uloulo', **kwargs)[:,:,:,m_index,mp_index,...], axis=(-1,-2)) if both_contours: return gf diff --git a/tests/tools/test_greensfunction.py b/tests/tools/test_greensfunction.py index 664bf0c04..15a1948fc 100644 --- a/tests/tools/test_greensfunction.py +++ b/tests/tools/test_greensfunction.py @@ -102,6 +102,11 @@ def test_greensfunction_radial(test_file): assert gf.trace_energy_dependence(spin=1).shape == (128,) assert gf.trace_energy_dependence(spin=1).dtype == float + assert isinstance(gf.energy_dependence(spin=1, radial=True), np.ndarray) + assert not np.isnan(gf.energy_dependence(spin=1, radial=True)).any() + assert gf.energy_dependence(spin=1, radial=True).shape == (757, 128, 5, 5) #(jr,nz,2*l+1,2*l+1) + assert gf.energy_dependence(spin=1, radial=True).dtype == float + def test_list_elements(test_file): """ @@ -196,6 +201,11 @@ def test_greensfunction_radial_complete_spin(test_file): assert gf.trace_energy_dependence().shape == (128, 2, 2) assert gf.trace_energy_dependence().dtype == float + assert isinstance(gf.energy_dependence(radial=True), np.ndarray) + assert not np.isnan(gf.energy_dependence(radial=True)).any() + assert gf.energy_dependence(radial=True).shape == (757, 128, 5, 5, 2, 2) #(jr,nz,2*l+1,2*l+1, spin1, spin2) + assert gf.energy_dependence(radial=True).dtype == float + def test_greensfunction_sphavg_full_spin_matrix(test_file): """ From cb5972d1b0aa15df4cec4a303f9b72e17ae61d07 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Fri, 16 Sep 2022 11:09:09 +0200 Subject: [PATCH 02/32] Allow onsite_delta to be a matrix in the spherical harmonics space Also added IO for constructing this full matrix from entry in greensf.hdf file --- masci_tools/tools/greensf_calculations.py | 93 ++++++++++++++++++++--- masci_tools/tools/greensfunction.py | 24 +++--- pyproject.toml | 3 +- 3 files changed, 98 insertions(+), 22 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 7e5d76a1c..00a34e0e3 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -12,12 +12,16 @@ from .greensfunction import GreensFunction, intersite_shells, intersite_shells_from_file from masci_tools.io.common_functions import get_pauli_matrix +from masci_tools.util.constants import HTR_TO_EV import numpy as np import pandas as pd from scipy import constants +from scipy.interpolate import interp1d +from sympy.physics.wigner import gaunt from collections import defaultdict from typing import Any +import h5py try: from typing import Literal @@ -34,7 +38,7 @@ def calculate_heisenberg_jij( hdffileORgreensfunctions: FileLike | list[GreensFunction], reference_atom: int, - onsite_delta: np.ndarray, + onsite_delta: np.ndarray | None = None, max_shells: int | None = None, ) -> pd.DataFrame: r""" @@ -56,6 +60,14 @@ def calculate_heisenberg_jij( else: shells = intersite_shells_from_file(hdffileORgreensfunctions, reference_atom, max_shells=max_shells) + if onsite_delta is None: + if isinstance(hdffileORgreensfunctions, list): + raise ValueError('For calculating the Bxc from the greensf.hdf file ' + 'The file has to be provided') + onsite_delta = calculate_bxc_mmp_matrix(hdffileORgreensfunctions) + + delta_is_matrix = len(onsite_delta.shape) > 2 + jij_constants: dict[str, list[Any]] = defaultdict(list) for dist, g1, g2 in shells: @@ -68,11 +80,17 @@ def calculate_heisenberg_jij( gij = g1.energy_dependence(both_contours=True, spin=1) gji = g2.energy_dependence(both_contours=True, spin=2) - delta_square = onsite_delta[g1.atomType - 1, g1.l] * onsite_delta[g1.atomTypep - 1, g1.l] weights = np.array([g1.weights, -g1.weights.conj()]).T - integral = np.einsum('zm,zijm,zjim->', weights, gij, gji) - jij = 0.5 * 1 / (8.0 * np.pi * 1j) * delta_square * integral + if delta_is_matrix: + delta_i = onsite_delta[g1.atomType - 1, g1.l] + delta_j = onsite_delta[g1.atomTypep - 1, g1.l] + integral = np.einsum('zm,no,zopm,pq,zqnm->', weights, delta_i, gij, delta_j, gji) + jij = 0.5 * 1 / (8.0 * np.pi * 1j) * integral + else: + delta_square = onsite_delta[g1.atomType - 1, g1.l] * onsite_delta[g1.atomTypep - 1, g1.l] + integral = np.einsum('zm,zijm,zjim->', weights, gij, gji) + jij = 0.5 * 1 / (8.0 * np.pi * 1j) * delta_square * integral jij_constants['R'].append(dist) jij_constants['R_ij_x'].append(g1.atomDiff.tolist()[0]) @@ -87,7 +105,7 @@ def calculate_heisenberg_jij( def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[GreensFunction], reference_atom: int, - onsite_delta: np.ndarray, + onsite_delta: np.ndarray | None = None, max_shells: int | None = None) -> pd.DataFrame: r""" Calculate the Heisenberg exchange tensor :math:`\mathbf{J}` from Green's functions using the formula @@ -110,6 +128,14 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens else: shells = intersite_shells_from_file(hdffileORgreensfunctions, reference_atom, max_shells=max_shells) + if onsite_delta is None: + if isinstance(hdffileORgreensfunctions, list): + raise ValueError('For calculating the Bxc from the greensf.hdf file ' + 'The file has to be provided') + onsite_delta = calculate_bxc_mmp_matrix(hdffileORgreensfunctions) + + delta_is_matrix = len(onsite_delta.shape) > 2 + jij_tensor: dict[str, list[Any]] = defaultdict(list) for dist, g1, g2 in shells: @@ -121,7 +147,12 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens gij = g1.energy_dependence(both_contours=True) gji = g2.energy_dependence(both_contours=True) - delta_square = onsite_delta[g1.atomType - 1, g1.l] * onsite_delta[g1.atomTypep - 1, g1.l] + if delta_is_matrix: + delta_i = onsite_delta[g1.atomType - 1, g1.l] + delta_j = onsite_delta[g1.atomTypep - 1, g1.l] + else: + delta_square = onsite_delta[g1.atomType - 1, g1.l] * onsite_delta[g1.atomTypep - 1, g1.l] + weights = np.array([g1.weights, -g1.weights.conj()]).T jij_tensor['R'].append(dist) @@ -136,9 +167,13 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens sigmai = get_pauli_matrix(sigmai_str) #type: ignore[arg-type] sigmaj = get_pauli_matrix(sigmaj_str) #type: ignore[arg-type] - - integral = np.einsum('zm,ab,zijbcm,cd,zjidam->', weights, sigmai, gij, sigmaj, gji) - jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * delta_square * integral + if delta_is_matrix: + integral = np.einsum('zm,no,ab,zopbcm,pq,cd,zqndam->', weights, delta_i, sigmai, gij, delta_j, + sigmaj, gji) + jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * integral + else: + integral = np.einsum('zm,ab,zijbcm,cd,zjidam->', weights, sigmai, gij, sigmaj, gji) + jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * delta_square * integral jij_tensor[f'J_{sigmai_str}{sigmaj_str}'].append(jij.real * 1000) #Convert to meV return pd.DataFrame.from_dict(jij_tensor) @@ -252,3 +287,43 @@ def calculate_hybridization(greensfunction: GreensFunction) -> np.ndarray: gz = greensfunction.trace_energy_dependence() delta = 1 / (2 * greensfunction.l + 1) * np.linalg.inv(gz) return delta.real + + +def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000) -> np.ndarray: + """ + Calculate the Bxc potential in the basis of products of spherical harmonics + (Same basis as greens functions) + + :param file: greensf.hdf file (atleast version 9) + :param radial_mesh_points: how many points are used for radial integration + """ + with h5py.File(file, 'r') as hdffile: + if hdffile.get('meta').attrs['version'] < 9: + raise ValueError('The greensf.hdf file must be atleast of version 9' + 'to extract the Bxc from it') + #Read in the Bxc and radial meshes + bxc = np.array(hdffile.get('bxc/data')) + bxc = bxc[..., 0] + 1j * bxc[..., 1] + bxc *= HTR_TO_EV + log_rmesh = np.array(hdffile.get('RadialFunctions/rmsh')) + + #Now calculate all the different atomtype and orbitals that could be needed + #Since this is relativaley cheap we can just do it for everything + + bxc_mmp = np.zeros((bxc.shape[-1], 4, 7, 7)) + + for atomtype, (bxc_atomtype, logmesh_atomtype) in enumerate(zip(bxc, log_rmesh)): + rmesh = np.arange(0, max(logmesh_atomtype), max(logmesh_atomtype) / radial_mesh_points) + for lrep in range(4): + for lpot in range(2 * lrep + 1): + for mpot in range(-lpot, lpot + 1): + lm = lpot * (lpot + 1) + mpot + bxc_interpolated = interp1d(logmesh_atomtype, bxc_atomtype[lm, :], fill_value='extrapolate') + bxc_integrated = np.trapz(bxc_interpolated(rmesh), rmesh) + + for m in range(2 * lrep + 1): + for mp in range(2 * lrep + 1): + gaunt_coeff = (-1)**mp * gaunt(lrep, lpot, lrep, -m + lrep, mpot, mp - lrep) + bxc_mmp[atomtype, lrep, m, mp] += gaunt_coeff * bxc_integrated + + return bxc_mmp diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 59bd8f937..02d04f0f7 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -414,20 +414,20 @@ def _read_element_header(hdffile: h5py.File, index: int) -> GreensfElement: """ group_name = _get_greensf_group_name(hdffile) - element = hdffile.get(f'/{group_name}/element-{index}') - - l = element.attrs['l'][0] - lp = element.attrs['lp'][0] - atomType = element.attrs['atomType'][0] - atomTypep = element.attrs['atomTypep'][0] - sphavg = element.attrs['l_sphavg'][0] == 1 - onsite = element.attrs['l_onsite'][0] == 1 - contour = element.attrs['iContour'][0] - kresolved = element.attrs.get('l_kresolved', [0])[0] == 1 - atomDiff = np.array(element.attrs['atomDiff']) + element = hdffile[f'/{group_name}/element-{index}'].attrs + + l = element['l'][0] + lp = element['lp'][0] + atomType = element['atomType'][0] + atomTypep = element['atomTypep'][0] + sphavg = element['l_sphavg'][0] == 1 + onsite = element['l_onsite'][0] == 1 + contour = element['iContour'][0] + kresolved = element.get('l_kresolved', [0])[0] == 1 + atomDiff = np.array(element['atomDiff']) atomDiff[abs(atomDiff) < 1e-12] = 0.0 atomDiff *= BOHR_A - nLO = element.attrs['numLOs'][0] + nLO = element['numLOs'][0] return GreensfElement(l, lp, atomType, atomTypep, sphavg, onsite, kresolved, contour, nLO, atomDiff) diff --git a/pyproject.toml b/pyproject.toml index d4d2c2796..e1ea8b793 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,6 +25,7 @@ requires-python = ">=3.7" dependencies = [ 'numpy', 'scipy', + 'sympy', 'matplotlib', 'h5py', 'pandas', @@ -122,7 +123,7 @@ module = [ 'lxml.builder', 'IPython.*', 'pandas', - 'scipy', + 'scipy.*', 'ase.*', 'spglib' ] From df1808bca6404473fc288fdbfa2a9740d5cbe72d Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Fri, 16 Sep 2022 11:26:04 +0200 Subject: [PATCH 03/32] Fix array dimensions and padding --- masci_tools/tools/greensf_calculations.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 00a34e0e3..19facd09d 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -83,8 +83,8 @@ def calculate_heisenberg_jij( weights = np.array([g1.weights, -g1.weights.conj()]).T if delta_is_matrix: - delta_i = onsite_delta[g1.atomType - 1, g1.l] - delta_j = onsite_delta[g1.atomTypep - 1, g1.l] + delta_i = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] integral = np.einsum('zm,no,zopm,pq,zqnm->', weights, delta_i, gij, delta_j, gji) jij = 0.5 * 1 / (8.0 * np.pi * 1j) * integral else: @@ -148,8 +148,8 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens gji = g2.energy_dependence(both_contours=True) if delta_is_matrix: - delta_i = onsite_delta[g1.atomType - 1, g1.l] - delta_j = onsite_delta[g1.atomTypep - 1, g1.l] + delta_i = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] else: delta_square = onsite_delta[g1.atomType - 1, g1.l] * onsite_delta[g1.atomTypep - 1, g1.l] @@ -310,7 +310,7 @@ def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000) -> #Now calculate all the different atomtype and orbitals that could be needed #Since this is relativaley cheap we can just do it for everything - bxc_mmp = np.zeros((bxc.shape[-1], 4, 7, 7)) + bxc_mmp = np.zeros((bxc.shape[0], 4, 7, 7)) for atomtype, (bxc_atomtype, logmesh_atomtype) in enumerate(zip(bxc, log_rmesh)): rmesh = np.arange(0, max(logmesh_atomtype), max(logmesh_atomtype) / radial_mesh_points) @@ -324,6 +324,6 @@ def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000) -> for m in range(2 * lrep + 1): for mp in range(2 * lrep + 1): gaunt_coeff = (-1)**mp * gaunt(lrep, lpot, lrep, -m + lrep, mpot, mp - lrep) - bxc_mmp[atomtype, lrep, m, mp] += gaunt_coeff * bxc_integrated + bxc_mmp[atomtype, lrep, 3 - lrep + m, 3 - lrep + mp] += gaunt_coeff * bxc_integrated return bxc_mmp From ade24c8d654d6d59aeafa90297186e609795fdb5 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Wed, 21 Sep 2022 13:52:32 +0200 Subject: [PATCH 04/32] Optimize einsum operations in jij calculation to avoid duplicate calculations --- masci_tools/tools/greensf_calculations.py | 71 +++++++++++++++++------ 1 file changed, 53 insertions(+), 18 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 19facd09d..fc8122351 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -11,7 +11,7 @@ from masci_tools.util.typing import FileLike from .greensfunction import GreensFunction, intersite_shells, intersite_shells_from_file -from masci_tools.io.common_functions import get_pauli_matrix +from masci_tools.io.common_functions import get_pauli_matrix, get_spin_rotation from masci_tools.util.constants import HTR_TO_EV import numpy as np @@ -83,14 +83,32 @@ def calculate_heisenberg_jij( weights = np.array([g1.weights, -g1.weights.conj()]).T if delta_is_matrix: - delta_i = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - integral = np.einsum('zm,no,zopm,pq,zqnm->', weights, delta_i, gij, delta_j, gji) - jij = 0.5 * 1 / (8.0 * np.pi * 1j) * integral + delta_i = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) + delta_j = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) + + delta_i[..., 0, 0] = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j[..., 0, 0] = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j[..., 1, 1] = delta_j[..., 0, 0] + delta_i[..., 1, 1] = delta_i[..., 0, 0] + + #Rotate into global spin frame + alpha, alphap = g1._angle_alpha #pylint: disable=protected-access + beta, betap = g1._angle_beta #pylint: disable=protected-access + + rot_spin = get_spin_rotation(-alpha, -beta) + rotp_spin = get_spin_rotation(-alphap, -betap) + + delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) + delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) + + gdeltaij = np.einsum('zij,jk...->zik...', gij, delta_i) + gdeltaji = np.einsum('zij,jk...->zik...', gji, delta_j) else: - delta_square = onsite_delta[g1.atomType - 1, g1.l] * onsite_delta[g1.atomTypep - 1, g1.l] - integral = np.einsum('zm,zijm,zjim->', weights, gij, gji) - jij = 0.5 * 1 / (8.0 * np.pi * 1j) * delta_square * integral + gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] + gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] + + integral = np.einsum('zm,zijm,zjim->', weights, gdeltaij, gdeltaji) + jij = 0.5 * 1 / (8.0 * np.pi * 1j) * integral jij_constants['R'].append(dist) jij_constants['R_ij_x'].append(g1.atomDiff.tolist()[0]) @@ -148,10 +166,30 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens gji = g2.energy_dependence(both_contours=True) if delta_is_matrix: - delta_i = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_i = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) + delta_j = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) + + delta_i[..., 0, 0] = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j[..., 0, 0] = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j[..., 1, 1] = delta_j[..., 0, 0] + delta_i[..., 1, 1] = delta_i[..., 0, 0] + + #Rotate into global spin frame + alpha, alphap = g1._angle_alpha #pylint: disable=protected-access + beta, betap = g1._angle_beta #pylint: disable=protected-access + + rot_spin = get_spin_rotation(-alpha, -beta) + rotp_spin = get_spin_rotation(-alphap, -betap) + + delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) + delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) + + gdeltaij = np.einsum('zij,jk...->zik...', gij, delta_i) + gdeltaji = np.einsum('zij,jk...->zik...', gji, delta_j) + else: - delta_square = onsite_delta[g1.atomType - 1, g1.l] * onsite_delta[g1.atomTypep - 1, g1.l] + gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] + gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] weights = np.array([g1.weights, -g1.weights.conj()]).T @@ -167,13 +205,10 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens sigmai = get_pauli_matrix(sigmai_str) #type: ignore[arg-type] sigmaj = get_pauli_matrix(sigmaj_str) #type: ignore[arg-type] - if delta_is_matrix: - integral = np.einsum('zm,no,ab,zopbcm,pq,cd,zqndam->', weights, delta_i, sigmai, gij, delta_j, - sigmaj, gji) - jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * integral - else: - integral = np.einsum('zm,ab,zijbcm,cd,zjidam->', weights, sigmai, gij, sigmaj, gji) - jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * delta_square * integral + + integral = np.einsum('zm,ab,zijbcm,cd,zjidam->', weights, sigmai, gdeltaij, sigmaj, gdeltaji) + jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * integral + jij_tensor[f'J_{sigmai_str}{sigmaj_str}'].append(jij.real * 1000) #Convert to meV return pd.DataFrame.from_dict(jij_tensor) From d576f465e0d765a5fb62de077e4142a81b497e22 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Wed, 21 Sep 2022 14:22:19 +0200 Subject: [PATCH 05/32] fix for last commit --- masci_tools/tools/greensf_calculations.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index fc8122351..fc8a714cd 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -101,8 +101,8 @@ def calculate_heisenberg_jij( delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - gdeltaij = np.einsum('zij,jk...->zik...', gij, delta_i) - gdeltaji = np.einsum('zij,jk...->zik...', gji, delta_j) + gdeltaij = np.einsum('ij...,zjk...->zik...', delta_i,gij) + gdeltaji = np.einsum('ij...,zjk...->zik...', delta_j,gji) else: gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] @@ -184,9 +184,8 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - gdeltaij = np.einsum('zij,jk...->zik...', gij, delta_i) - gdeltaji = np.einsum('zij,jk...->zik...', gji, delta_j) - + gdeltaij = np.einsum('ij...,zjk...->zik...', delta_i,gij) + gdeltaji = np.einsum('ij...,zjk...->zik...', delta_j,gji) else: gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] From f44ada7d0bcf6e3fc4fc7a666f824ab63d188876 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Wed, 21 Sep 2022 14:27:27 +0200 Subject: [PATCH 06/32] Add missing rotation --- masci_tools/tools/greensf_calculations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index fc8a714cd..c943953c1 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -101,8 +101,8 @@ def calculate_heisenberg_jij( delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - gdeltaij = np.einsum('ij...,zjk...->zik...', delta_i,gij) - gdeltaji = np.einsum('ij...,zjk...->zik...', delta_j,gji) + gdeltaij = np.einsum('ijab,zjkbc->zikac', delta_i,gij) + gdeltaji = np.einsum('ijab,zjkbc->zikac', delta_j,gji) else: gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] @@ -184,8 +184,8 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - gdeltaij = np.einsum('ij...,zjk...->zik...', delta_i,gij) - gdeltaji = np.einsum('ij...,zjk...->zik...', delta_j,gji) + gdeltaij = np.einsum('ijab,zjkbc->zikac', delta_i,gij) + gdeltaji = np.einsum('ijab,zjkbc->zikac', delta_j,gji) else: gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] From f8354ab7145919ef5956522255e7c3b151b708ff Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Wed, 21 Sep 2022 14:31:45 +0200 Subject: [PATCH 07/32] And another fix --- masci_tools/tools/greensf_calculations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index c943953c1..856a66dba 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -101,8 +101,8 @@ def calculate_heisenberg_jij( delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - gdeltaij = np.einsum('ijab,zjkbc->zikac', delta_i,gij) - gdeltaji = np.einsum('ijab,zjkbc->zikac', delta_j,gji) + gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i,gij) + gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j,gji) else: gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] @@ -184,8 +184,8 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - gdeltaij = np.einsum('ijab,zjkbc->zikac', delta_i,gij) - gdeltaji = np.einsum('ijab,zjkbc->zikac', delta_j,gji) + gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i,gij) + gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j,gji) else: gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] From a520bdc6ed09f9d9f6db2f6f8f5724230c4d379b Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 23 Sep 2022 13:03:02 +0000 Subject: [PATCH 08/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- masci_tools/tools/greensf_calculations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 856a66dba..be3216172 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -101,8 +101,8 @@ def calculate_heisenberg_jij( delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i,gij) - gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j,gji) + gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) + gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) else: gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] @@ -184,8 +184,8 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i,gij) - gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j,gji) + gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) + gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) else: gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] From f3451399fc1e225ac8a34cab8eedf6ad97a75869 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Fri, 23 Sep 2022 15:35:28 +0200 Subject: [PATCH 09/32] fix last commits --- masci_tools/tools/greensf_calculations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index be3216172..01f08773b 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -104,8 +104,8 @@ def calculate_heisenberg_jij( gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) else: - gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] - gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] + gdeltaij = onsite_delta[g1.atomType - 1, g1.l] * gij + gdeltaji = onsite_delta[g1.atomTypep - 1, g1.l] * gji integral = np.einsum('zm,zijm,zjim->', weights, gdeltaij, gdeltaji) jij = 0.5 * 1 / (8.0 * np.pi * 1j) * integral @@ -187,8 +187,8 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) else: - gdeltaij *= onsite_delta[g1.atomType - 1, g1.l] - gdeltaji *= onsite_delta[g1.atomTypep - 1, g1.l] + gdeltaij = onsite_delta[g1.atomType - 1, g1.l] * gij + gdeltaji = onsite_delta[g1.atomTypep - 1, g1.l] * gji weights = np.array([g1.weights, -g1.weights.conj()]).T From c7d407b47c5501a45f45fa8f0287c3b0ebd7e5f4 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 08:33:54 +0200 Subject: [PATCH 10/32] Allow passing transformation functions to jij calculation Useful for calculating jij decompositions (e.g. eg/t2g) --- masci_tools/tools/greensf_calculations.py | 26 +++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 01f08773b..89f5cbc8b 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -20,7 +20,7 @@ from scipy.interpolate import interp1d from sympy.physics.wigner import gaunt from collections import defaultdict -from typing import Any +from typing import Any, Callable import h5py try: @@ -40,6 +40,7 @@ def calculate_heisenberg_jij( reference_atom: int, onsite_delta: np.ndarray | None = None, max_shells: int | None = None, + transform_func: Callable | None = None ) -> pd.DataFrame: r""" Calculate the Heisenberg exchange constants form Green's functions using the formula @@ -51,6 +52,8 @@ def calculate_heisenberg_jij( :param reference_atom: integer index of the atom to calculate the Jijs from :param onsite_delta: List of floats containing the onsite exchange splitting for each atom type and l-channel :param max_shells: optional int, if given only the first max_shells shells are constructed + :param transform_func: Callable, if given the green's functions (and bxc if it's a matrix) are passed to it + and are expected to be returned in a compatible shape :returns: pandas DataFrame containing all the Jij constants """ @@ -80,6 +83,10 @@ def calculate_heisenberg_jij( gij = g1.energy_dependence(both_contours=True, spin=1) gji = g2.energy_dependence(both_contours=True, spin=2) + if transform_func is not None: + gij = transform_func(gij) + gji = transform_func(gji) + weights = np.array([g1.weights, -g1.weights.conj()]).T if delta_is_matrix: @@ -101,6 +108,10 @@ def calculate_heisenberg_jij( delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) + if transform_func is not None: + delta_i = transform_func(delta_i) + delta_i = transform_func(delta_i) + gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) else: @@ -124,7 +135,8 @@ def calculate_heisenberg_jij( def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[GreensFunction], reference_atom: int, onsite_delta: np.ndarray | None = None, - max_shells: int | None = None) -> pd.DataFrame: + max_shells: int | None = None, + transform_func: Callable | None = None) -> pd.DataFrame: r""" Calculate the Heisenberg exchange tensor :math:`\mathbf{J}` from Green's functions using the formula @@ -137,6 +149,8 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens :param reference_atom: integer index of the atom to calculate the Jijs from :param onsite_delta: List of floats containing the onsite exchange splitting for each atom type and l-channel :param max_shells: optional int, if given only the first max_shells shells are constructed + :param transform_func: Callable, if given the green's functions (and bxc if it's a matrix) are passed to it + and are expected to be returned in a compatible shape :returns: pandas DataFrame containing all the J_xx, J_xy, etc. constants """ @@ -165,6 +179,10 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens gij = g1.energy_dependence(both_contours=True) gji = g2.energy_dependence(both_contours=True) + if transform_func is not None: + gij = transform_func(gij) + gji = transform_func(gji) + if delta_is_matrix: delta_i = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) delta_j = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) @@ -184,6 +202,10 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) + if transform_func is not None: + delta_i = transform_func(delta_i) + delta_i = transform_func(delta_i) + gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) else: From fa9c1a74053cb74e1c18a40d83038485d9613533 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 26 Sep 2022 06:35:12 +0000 Subject: [PATCH 11/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- masci_tools/tools/greensf_calculations.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 89f5cbc8b..4c19cb8c4 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -35,13 +35,11 @@ # - custom decompositions (e.g. e2g/t2g) -def calculate_heisenberg_jij( - hdffileORgreensfunctions: FileLike | list[GreensFunction], - reference_atom: int, - onsite_delta: np.ndarray | None = None, - max_shells: int | None = None, - transform_func: Callable | None = None -) -> pd.DataFrame: +def calculate_heisenberg_jij(hdffileORgreensfunctions: FileLike | list[GreensFunction], + reference_atom: int, + onsite_delta: np.ndarray | None = None, + max_shells: int | None = None, + transform_func: Callable | None = None) -> pd.DataFrame: r""" Calculate the Heisenberg exchange constants form Green's functions using the formula From 432610aae773f65e8bca7be31485b85417a65360 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 08:45:46 +0200 Subject: [PATCH 12/32] fix typo --- masci_tools/tools/greensf_calculations.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 4c19cb8c4..b58609765 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -108,7 +108,7 @@ def calculate_heisenberg_jij(hdffileORgreensfunctions: FileLike | list[GreensFun if transform_func is not None: delta_i = transform_func(delta_i) - delta_i = transform_func(delta_i) + delta_j = transform_func(delta_j) gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) @@ -202,7 +202,7 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens if transform_func is not None: delta_i = transform_func(delta_i) - delta_i = transform_func(delta_i) + delta_j = transform_func(delta_j) gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) From 368636ba933abf8c245ab3ce4d9eed9ddfc5a877 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 08:49:40 +0200 Subject: [PATCH 13/32] fix dtype --- masci_tools/tools/greensf_calculations.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index b58609765..9a6b1cc65 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -364,7 +364,7 @@ def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000) -> #Now calculate all the different atomtype and orbitals that could be needed #Since this is relativaley cheap we can just do it for everything - bxc_mmp = np.zeros((bxc.shape[0], 4, 7, 7)) + bxc_mmp = np.zeros((bxc.shape[0], 4, 7, 7), dtype=complex) for atomtype, (bxc_atomtype, logmesh_atomtype) in enumerate(zip(bxc, log_rmesh)): rmesh = np.arange(0, max(logmesh_atomtype), max(logmesh_atomtype) / radial_mesh_points) @@ -380,4 +380,4 @@ def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000) -> gaunt_coeff = (-1)**mp * gaunt(lrep, lpot, lrep, -m + lrep, mpot, mp - lrep) bxc_mmp[atomtype, lrep, 3 - lrep + m, 3 - lrep + mp] += gaunt_coeff * bxc_integrated - return bxc_mmp + return bxc_mmp.real From 829b5ebcb202006ff264e99a545d3b55172279ed Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 09:43:49 +0200 Subject: [PATCH 14/32] Fix rotation of onsite exchange splitting --- masci_tools/tools/greensf_calculations.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 9a6b1cc65..b82c36186 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -182,14 +182,17 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens gji = transform_func(gji) if delta_is_matrix: - delta_i = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) - delta_j = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) + delta_i = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2, 2), dtype=complex) + delta_j = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2, 2), dtype=complex) - delta_i[..., 0, 0] = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j[..., 0, 0] = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_i[..., 0, 0] = -onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j[..., 0, 0] = -onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] delta_j[..., 1, 1] = delta_j[..., 0, 0] delta_i[..., 1, 1] = delta_i[..., 0, 0] + delta_i[..., 1] = delta_i[..., 0].conj() + delta_j[..., 1] = delta_j[..., 0].conj() + #Rotate into global spin frame alpha, alphap = g1._angle_alpha #pylint: disable=protected-access beta, betap = g1._angle_beta #pylint: disable=protected-access @@ -197,15 +200,15 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens rot_spin = get_spin_rotation(-alpha, -beta) rotp_spin = get_spin_rotation(-alphap, -betap) - delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) - delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) + delta_i = np.einsum('ij,xyjk...,km->xyim...', rot_spin, delta_i, rot_spin.T.conj()) + delta_j = np.einsum('ij,xyjk...,km->xyim...', rotp_spin, delta_j, rotp_spin.T.conj()) if transform_func is not None: delta_i = transform_func(delta_i) delta_j = transform_func(delta_j) - gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) - gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) + gdeltaij = np.einsum('ijabm,zjkbcm->zikacm', delta_i, gij) + gdeltaji = np.einsum('ijabm,zjkbcm->zikacm', delta_j, gji) else: gdeltaij = onsite_delta[g1.atomType - 1, g1.l] * gij gdeltaji = onsite_delta[g1.atomTypep - 1, g1.l] * gji From 31c640dc07ce0bc19c05841fb616b24cc08a117d Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 09:44:18 +0200 Subject: [PATCH 15/32] Precalculate energy integrals in tensor calculation --- masci_tools/tools/greensf_calculations.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index b82c36186..ee32b145e 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -222,13 +222,15 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens jij_tensor['Atom i'].append(f"{g1.extras['atom_label']}({g1.extras['element']})") jij_tensor['Atom j'].append(f"{g1.extras['atom_labelp']}({g1.extras['elementp']})") + gdeltaij = np.einsum('zm,zijbcm->', weights, gdeltaij) + gdeltaji = np.einsum('zm,zijbcm->', weights, gdeltaji) for sigmai_str in ('x', 'y', 'z'): for sigmaj_str in ('x', 'y', 'z'): sigmai = get_pauli_matrix(sigmai_str) #type: ignore[arg-type] sigmaj = get_pauli_matrix(sigmaj_str) #type: ignore[arg-type] - integral = np.einsum('zm,ab,zijbcm,cd,zjidam->', weights, sigmai, gdeltaij, sigmaj, gdeltaji) + integral = np.einsum('ab,ijbc,cd,jida->', weights, sigmai, gdeltaij, sigmaj, gdeltaji) jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * integral jij_tensor[f'J_{sigmai_str}{sigmaj_str}'].append(jij.real * 1000) #Convert to meV From ded5f5334defbba05ca3d1204dbaf52148ea68a9 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 09:46:11 +0200 Subject: [PATCH 16/32] fix last commits --- masci_tools/tools/greensf_calculations.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index ee32b145e..f6c77ae13 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -185,10 +185,10 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens delta_i = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2, 2), dtype=complex) delta_j = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2, 2), dtype=complex) - delta_i[..., 0, 0] = -onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j[..., 0, 0] = -onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j[..., 1, 1] = delta_j[..., 0, 0] - delta_i[..., 1, 1] = delta_i[..., 0, 0] + delta_i[..., 0, 0, 0] = -onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j[..., 0, 0, 0] = -onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j[..., 1, 1, 0] = delta_j[..., 0, 0, 0] + delta_i[..., 1, 1, 0] = delta_i[..., 0, 0, 0] delta_i[..., 1] = delta_i[..., 0].conj() delta_j[..., 1] = delta_j[..., 0].conj() From 5f2550d735dafcd79d064b4bc68a2ff49ad00a6f Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 09:47:55 +0200 Subject: [PATCH 17/32] fix einsum --- masci_tools/tools/greensf_calculations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index f6c77ae13..8ed087a6c 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -230,7 +230,7 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens sigmai = get_pauli_matrix(sigmai_str) #type: ignore[arg-type] sigmaj = get_pauli_matrix(sigmaj_str) #type: ignore[arg-type] - integral = np.einsum('ab,ijbc,cd,jida->', weights, sigmai, gdeltaij, sigmaj, gdeltaji) + integral = np.einsum('ab,ijbc,cd,jida->', sigmai, gdeltaij, sigmaj, gdeltaji) jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * integral jij_tensor[f'J_{sigmai_str}{sigmaj_str}'].append(jij.real * 1000) #Convert to meV From 58005c80c75473d9ed394afbb6607ab859ab2bf3 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 09:50:09 +0200 Subject: [PATCH 18/32] revert --- masci_tools/tools/greensf_calculations.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 8ed087a6c..956c5fc35 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -222,15 +222,13 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens jij_tensor['Atom i'].append(f"{g1.extras['atom_label']}({g1.extras['element']})") jij_tensor['Atom j'].append(f"{g1.extras['atom_labelp']}({g1.extras['elementp']})") - gdeltaij = np.einsum('zm,zijbcm->', weights, gdeltaij) - gdeltaji = np.einsum('zm,zijbcm->', weights, gdeltaji) for sigmai_str in ('x', 'y', 'z'): for sigmaj_str in ('x', 'y', 'z'): sigmai = get_pauli_matrix(sigmai_str) #type: ignore[arg-type] sigmaj = get_pauli_matrix(sigmaj_str) #type: ignore[arg-type] - integral = np.einsum('ab,ijbc,cd,jida->', sigmai, gdeltaij, sigmaj, gdeltaji) + integral = np.einsum('zm,ab,zijbcm,cd,zjidam->', weights, sigmai, gdeltaij, sigmaj, gdeltaji) jij = 1 / 4 * 1 / (8.0 * np.pi * 1j) * integral jij_tensor[f'J_{sigmai_str}{sigmaj_str}'].append(jij.real * 1000) #Convert to meV From 4dd3be73ca46965ac8ad01a79be19bdd7e9e7a4e Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 10:30:43 +0200 Subject: [PATCH 19/32] Add option to average diagonal --- masci_tools/tools/greensf_calculations.py | 9 ++++++++- masci_tools/tools/greensfunction.py | 12 ++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 956c5fc35..cad3654a5 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -134,7 +134,8 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens reference_atom: int, onsite_delta: np.ndarray | None = None, max_shells: int | None = None, - transform_func: Callable | None = None) -> pd.DataFrame: + transform_func: Callable | None = None, + average_diagonal: bool = False) -> pd.DataFrame: r""" Calculate the Heisenberg exchange tensor :math:`\mathbf{J}` from Green's functions using the formula @@ -172,6 +173,12 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens dist = round(dist, 12) + if average_diagonal: + g1.to_local_frame() + g2.to_local_frame() + g1.average_spindiagonal() + g2.average_spindiagonal() + g1.to_global_frame() g2.to_global_frame() gij = g1.energy_dependence(both_contours=True) diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 02d04f0f7..1db3629cb 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -662,6 +662,18 @@ def _set_spin_matrix(self, name: CoefficientName, spin_matrix: np.ndarray) -> No data[:, :, :, 2, ...] = spin_matrix[:, :, :, 1, 0, ...] data[:, :, :, 3, ...] = spin_matrix[:, :, :, 0, 1, ...] + def average_spindiagonal(self) -> None: + """ + Average out ferromagnetic contributions on the local spin-diagonals + """ + + for name in self._data.keys(): + data = self._get_spin_matrix(name) + for m in range(2 * self.l + 1): + data[:, m, m, 0, 0] = (data[:, m, m, 0, 0] + data[:, m, m, 1, 1]) / 2 + data[:, m, m, 1, 1] = data[:, m, m, 0, 0] + self._set_spin_matrix(name, data) + def to_global_frame(self) -> None: """ Rotate the Green's function into the global real space and spin space frame From e66d97aaf2a1d9bed91eb4c100653e584a8e9e6a Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 11:05:19 +0200 Subject: [PATCH 20/32] fix --- masci_tools/tools/greensfunction.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 1db3629cb..c0394e17e 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -669,7 +669,7 @@ def average_spindiagonal(self) -> None: for name in self._data.keys(): data = self._get_spin_matrix(name) - for m in range(2 * self.l + 1): + for m in range(self.lmax): data[:, m, m, 0, 0] = (data[:, m, m, 0, 0] + data[:, m, m, 1, 1]) / 2 data[:, m, m, 1, 1] = data[:, m, m, 0, 0] self._set_spin_matrix(name, data) From b804f46ea2d09f6641d6e04ee7bd0980ebea5198 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 11:09:23 +0200 Subject: [PATCH 21/32] add possibility of cutoff --- masci_tools/tools/greensf_calculations.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index cad3654a5..a74498294 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -353,7 +353,7 @@ def calculate_hybridization(greensfunction: GreensFunction) -> np.ndarray: return delta.real -def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000) -> np.ndarray: +def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000, cutoff: float | None=None) -> np.ndarray: """ Calculate the Bxc potential in the basis of products of spherical harmonics (Same basis as greens functions) @@ -390,4 +390,7 @@ def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000) -> gaunt_coeff = (-1)**mp * gaunt(lrep, lpot, lrep, -m + lrep, mpot, mp - lrep) bxc_mmp[atomtype, lrep, 3 - lrep + m, 3 - lrep + mp] += gaunt_coeff * bxc_integrated - return bxc_mmp.real + bxc_mmp = bxc_mmp.real + if cutoff is not None: + bxc_mmp[np.abs(bxc_mmp) Date: Mon, 26 Sep 2022 09:10:34 +0000 Subject: [PATCH 22/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- masci_tools/tools/greensf_calculations.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index a74498294..3d0628041 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -353,7 +353,7 @@ def calculate_hybridization(greensfunction: GreensFunction) -> np.ndarray: return delta.real -def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000, cutoff: float | None=None) -> np.ndarray: +def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000, cutoff: float | None = None) -> np.ndarray: """ Calculate the Bxc potential in the basis of products of spherical harmonics (Same basis as greens functions) @@ -392,5 +392,5 @@ def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000, cut bxc_mmp = bxc_mmp.real if cutoff is not None: - bxc_mmp[np.abs(bxc_mmp) Date: Mon, 26 Sep 2022 11:57:56 +0200 Subject: [PATCH 23/32] fix --- masci_tools/tools/greensf_calculations.py | 2 -- masci_tools/tools/greensfunction.py | 18 ++++++++++++++---- 2 files changed, 14 insertions(+), 6 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 3d0628041..1f1e84715 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -174,8 +174,6 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens dist = round(dist, 12) if average_diagonal: - g1.to_local_frame() - g2.to_local_frame() g1.average_spindiagonal() g2.average_spindiagonal() diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index c0394e17e..5f37c5d1a 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -667,12 +667,22 @@ def average_spindiagonal(self) -> None: Average out ferromagnetic contributions on the local spin-diagonals """ + alpha, alphap = self._angle_alpha + beta, betap = self._angle_beta + if not self._local_real_frame: + rot_real_space = get_wigner_matrix(self.l, alpha, beta) + rotp_real_space = get_wigner_matrix(self.lp, alphap, betap) + + for name, data in self._data.items(): + data = np.einsum('ij,xjk...,km->xim...', rot_real_space.T.conj(), data, rotp_real_space) + self._data[name] = data + self._local_real_frame = True + for name in self._data.keys(): - data = self._get_spin_matrix(name) for m in range(self.lmax): - data[:, m, m, 0, 0] = (data[:, m, m, 0, 0] + data[:, m, m, 1, 1]) / 2 - data[:, m, m, 1, 1] = data[:, m, m, 0, 0] - self._set_spin_matrix(name, data) + self._data[name][:, m, m, 0] = (self._data[name][:, m, m, 0] + self._data[name][:, m, m, min(1,self.nspins-1)]) / 2 + if self.nspins > 1: + self._data[name][:, m, m, 1] = self._data[name][:, m, m, 0] def to_global_frame(self) -> None: """ From da4a582c4ba43fa596fd479202e8b2e53d34e334 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 26 Sep 2022 09:59:15 +0000 Subject: [PATCH 24/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- masci_tools/tools/greensfunction.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 5f37c5d1a..3deced157 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -680,7 +680,8 @@ def average_spindiagonal(self) -> None: for name in self._data.keys(): for m in range(self.lmax): - self._data[name][:, m, m, 0] = (self._data[name][:, m, m, 0] + self._data[name][:, m, m, min(1,self.nspins-1)]) / 2 + self._data[name][:, m, m, 0] = (self._data[name][:, m, m, 0] + + self._data[name][:, m, m, min(1, self.nspins - 1)]) / 2 if self.nspins > 1: self._data[name][:, m, m, 1] = self._data[name][:, m, m, 0] From b4895916c5ae92aa311f86300c850253b4a58336 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 13:48:26 +0200 Subject: [PATCH 25/32] try --- masci_tools/tools/greensf_calculations.py | 5 +++-- masci_tools/tools/greensfunction.py | 20 ++++++++++---------- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 1f1e84715..a498e9ead 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -173,12 +173,13 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens dist = round(dist, 12) + g1.to_global_frame() + g2.to_global_frame() + if average_diagonal: g1.average_spindiagonal() g2.average_spindiagonal() - g1.to_global_frame() - g2.to_global_frame() gij = g1.energy_dependence(both_contours=True) gji = g2.energy_dependence(both_contours=True) diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 3deced157..146077ee5 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -667,16 +667,16 @@ def average_spindiagonal(self) -> None: Average out ferromagnetic contributions on the local spin-diagonals """ - alpha, alphap = self._angle_alpha - beta, betap = self._angle_beta - if not self._local_real_frame: - rot_real_space = get_wigner_matrix(self.l, alpha, beta) - rotp_real_space = get_wigner_matrix(self.lp, alphap, betap) - - for name, data in self._data.items(): - data = np.einsum('ij,xjk...,km->xim...', rot_real_space.T.conj(), data, rotp_real_space) - self._data[name] = data - self._local_real_frame = True + # alpha, alphap = self._angle_alpha + # beta, betap = self._angle_beta + # if not self._local_real_frame: + # rot_real_space = get_wigner_matrix(self.l, alpha, beta) + # rotp_real_space = get_wigner_matrix(self.lp, alphap, betap) + + # for name, data in self._data.items(): + # data = np.einsum('ij,xjk...,km->xim...', rot_real_space.T.conj(), data, rotp_real_space) + # self._data[name] = data + # self._local_real_frame = True for name in self._data.keys(): for m in range(self.lmax): From 62a541db74337c4ea2e0f78da96570ee1966f0f4 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Mon, 26 Sep 2022 14:16:15 +0200 Subject: [PATCH 26/32] correct isotropic calculation --- masci_tools/tools/greensf_calculations.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index a498e9ead..15262c421 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -88,13 +88,8 @@ def calculate_heisenberg_jij(hdffileORgreensfunctions: FileLike | list[GreensFun weights = np.array([g1.weights, -g1.weights.conj()]).T if delta_is_matrix: - delta_i = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) - delta_j = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2), dtype=complex) - - delta_i[..., 0, 0] = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j[..., 0, 0] = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j[..., 1, 1] = delta_j[..., 0, 0] - delta_i[..., 1, 1] = delta_i[..., 0, 0] + delta_i = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] #Rotate into global spin frame alpha, alphap = g1._angle_alpha #pylint: disable=protected-access @@ -110,8 +105,8 @@ def calculate_heisenberg_jij(hdffileORgreensfunctions: FileLike | list[GreensFun delta_i = transform_func(delta_i) delta_j = transform_func(delta_j) - gdeltaij = np.einsum('ijab,zjkbc...->zikac...', delta_i, gij) - gdeltaji = np.einsum('ijab,zjkbc...->zikac...', delta_j, gji) + gdeltaij = np.einsum('ij,zjk...->zik...', delta_i, gij) + gdeltaji = np.einsum('ij,zjk...->zik...', delta_j, gji) else: gdeltaij = onsite_delta[g1.atomType - 1, g1.l] * gij gdeltaji = onsite_delta[g1.atomTypep - 1, g1.l] * gji From 3e3b99358456abee42e5c6cef9052da646487039 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Tue, 27 Sep 2022 09:19:28 +0200 Subject: [PATCH 27/32] more fixes --- masci_tools/tools/greensf_calculations.py | 10 -------- masci_tools/tools/greensfunction.py | 28 +++++++++++++++-------- 2 files changed, 18 insertions(+), 20 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 15262c421..5597cf0a4 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -91,16 +91,6 @@ def calculate_heisenberg_jij(hdffileORgreensfunctions: FileLike | list[GreensFun delta_i = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - #Rotate into global spin frame - alpha, alphap = g1._angle_alpha #pylint: disable=protected-access - beta, betap = g1._angle_beta #pylint: disable=protected-access - - rot_spin = get_spin_rotation(-alpha, -beta) - rotp_spin = get_spin_rotation(-alphap, -betap) - - delta_i = np.einsum('ij,xyjk,km->xyim', rot_spin, delta_i, rot_spin.T.conj()) - delta_j = np.einsum('ij,xyjk,km->xyim', rotp_spin, delta_j, rotp_spin.T.conj()) - if transform_func is not None: delta_i = transform_func(delta_i) delta_j = transform_func(delta_j) diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 146077ee5..d4c5113ac 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -667,24 +667,32 @@ def average_spindiagonal(self) -> None: Average out ferromagnetic contributions on the local spin-diagonals """ - # alpha, alphap = self._angle_alpha - # beta, betap = self._angle_beta - # if not self._local_real_frame: - # rot_real_space = get_wigner_matrix(self.l, alpha, beta) - # rotp_real_space = get_wigner_matrix(self.lp, alphap, betap) + alpha, alphap = self._angle_alpha + beta, betap = self._angle_beta + if not self._local_real_frame: + rot_real_space = get_wigner_matrix(self.l, alpha, beta) + rotp_real_space = get_wigner_matrix(self.lp, alphap, betap) - # for name, data in self._data.items(): - # data = np.einsum('ij,xjk...,km->xim...', rot_real_space.T.conj(), data, rotp_real_space) - # self._data[name] = data - # self._local_real_frame = True + for name, data in self._data.items(): + data = np.einsum('ij,xjk...,km->xim...', rot_real_space.T.conj(), data, rotp_real_space) + self._data[name] = data for name in self._data.keys(): - for m in range(self.lmax): + for m in range(2*self.lmax+1): self._data[name][:, m, m, 0] = (self._data[name][:, m, m, 0] + self._data[name][:, m, m, min(1, self.nspins - 1)]) / 2 if self.nspins > 1: self._data[name][:, m, m, 1] = self._data[name][:, m, m, 0] + if not self._local_real_frame: + rot_real_space = get_wigner_matrix(self.l, alpha, beta, inverse=True) + rotp_real_space = get_wigner_matrix(self.lp, alphap, betap, inverse=True) + + for name, data in self._data.items(): + data = np.einsum('ij,xjk...,km->xim...', rot_real_space.T.conj(), data, rotp_real_space) + self._data[name] = data + + def to_global_frame(self) -> None: """ Rotate the Green's function into the global real space and spin space frame From 36fc0788cb3186a75ccfd83c6c634be6fc8d3d00 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 27 Sep 2022 07:20:52 +0000 Subject: [PATCH 28/32] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- masci_tools/tools/greensfunction.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index d4c5113ac..0588352e4 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -678,7 +678,7 @@ def average_spindiagonal(self) -> None: self._data[name] = data for name in self._data.keys(): - for m in range(2*self.lmax+1): + for m in range(2 * self.lmax + 1): self._data[name][:, m, m, 0] = (self._data[name][:, m, m, 0] + self._data[name][:, m, m, min(1, self.nspins - 1)]) / 2 if self.nspins > 1: @@ -692,7 +692,6 @@ def average_spindiagonal(self) -> None: data = np.einsum('ij,xjk...,km->xim...', rot_real_space.T.conj(), data, rotp_real_space) self._data[name] = data - def to_global_frame(self) -> None: """ Rotate the Green's function into the global real space and spin space frame From c0eab84f3c44c6ae86fd7845b171923791afb6d8 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Wed, 26 Apr 2023 09:46:50 +0200 Subject: [PATCH 29/32] Revert changes for onsite delta in tensor calculation --- masci_tools/tools/greensf_calculations.py | 56 ++++++----------------- 1 file changed, 15 insertions(+), 41 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 5597cf0a4..bdfe2d673 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -81,26 +81,22 @@ def calculate_heisenberg_jij(hdffileORgreensfunctions: FileLike | list[GreensFun gij = g1.energy_dependence(both_contours=True, spin=1) gji = g2.energy_dependence(both_contours=True, spin=2) - if transform_func is not None: - gij = transform_func(gij) - gji = transform_func(gji) - weights = np.array([g1.weights, -g1.weights.conj()]).T if delta_is_matrix: delta_i = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - if transform_func is not None: - delta_i = transform_func(delta_i) - delta_j = transform_func(delta_j) - gdeltaij = np.einsum('ij,zjk...->zik...', delta_i, gij) gdeltaji = np.einsum('ij,zjk...->zik...', delta_j, gji) else: gdeltaij = onsite_delta[g1.atomType - 1, g1.l] * gij gdeltaji = onsite_delta[g1.atomTypep - 1, g1.l] * gji + if transform_func is not None: + gdeltaij = transform_func(gdeltaij) + gdeltaji = transform_func(gdeltaji) + integral = np.einsum('zm,zijm,zjim->', weights, gdeltaij, gdeltaji) jij = 0.5 * 1 / (8.0 * np.pi * 1j) * integral @@ -158,52 +154,30 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens dist = round(dist, 12) - g1.to_global_frame() - g2.to_global_frame() - if average_diagonal: g1.average_spindiagonal() g2.average_spindiagonal() + g1.to_global_frame() + g2.to_global_frame() + gij = g1.energy_dependence(both_contours=True) gji = g2.energy_dependence(both_contours=True) - if transform_func is not None: - gij = transform_func(gij) - gji = transform_func(gji) - if delta_is_matrix: - delta_i = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2, 2), dtype=complex) - delta_j = np.zeros((2 * g1.l + 1, 2 * g1.l + 1, 2, 2, 2), dtype=complex) - - delta_i[..., 0, 0, 0] = -onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j[..., 0, 0, 0] = -onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - delta_j[..., 1, 1, 0] = delta_j[..., 0, 0, 0] - delta_i[..., 1, 1, 0] = delta_i[..., 0, 0, 0] - - delta_i[..., 1] = delta_i[..., 0].conj() - delta_j[..., 1] = delta_j[..., 0].conj() - - #Rotate into global spin frame - alpha, alphap = g1._angle_alpha #pylint: disable=protected-access - beta, betap = g1._angle_beta #pylint: disable=protected-access - - rot_spin = get_spin_rotation(-alpha, -beta) - rotp_spin = get_spin_rotation(-alphap, -betap) - - delta_i = np.einsum('ij,xyjk...,km->xyim...', rot_spin, delta_i, rot_spin.T.conj()) - delta_j = np.einsum('ij,xyjk...,km->xyim...', rotp_spin, delta_j, rotp_spin.T.conj()) - - if transform_func is not None: - delta_i = transform_func(delta_i) - delta_j = transform_func(delta_j) + delta_i = onsite_delta[g1.atomType - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] + delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] - gdeltaij = np.einsum('ijabm,zjkbcm->zikacm', delta_i, gij) - gdeltaji = np.einsum('ijabm,zjkbcm->zikacm', delta_j, gji) + gdeltaij = np.einsum('ij,zjk...->zik...', delta_i, gij) + gdeltaji = np.einsum('ij,zjk...->zik:::', delta_j, gji) else: gdeltaij = onsite_delta[g1.atomType - 1, g1.l] * gij gdeltaji = onsite_delta[g1.atomTypep - 1, g1.l] * gji + if transform_func is not None: + gdeltaij = transform_func(gdeltaij) + gdeltaji = transform_func(gdeltaji) + weights = np.array([g1.weights, -g1.weights.conj()]).T jij_tensor['R'].append(dist) From 78e12474b72c3027d1086b6e4b0d116411739e80 Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Wed, 26 Apr 2023 09:47:27 +0200 Subject: [PATCH 30/32] fix bxc calculation by including mesh dimensions --- masci_tools/tools/greensf_calculations.py | 5 +++-- masci_tools/tools/greensfunction.py | 1 + 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index bdfe2d673..b881100f3 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -328,19 +328,20 @@ def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000, cut bxc = bxc[..., 0] + 1j * bxc[..., 1] bxc *= HTR_TO_EV log_rmesh = np.array(hdffile.get('RadialFunctions/rmsh')) + rmesh_dimensions = np.array(hdffile.get('RadialFunctions/jri')) #Now calculate all the different atomtype and orbitals that could be needed #Since this is relativaley cheap we can just do it for everything bxc_mmp = np.zeros((bxc.shape[0], 4, 7, 7), dtype=complex) - for atomtype, (bxc_atomtype, logmesh_atomtype) in enumerate(zip(bxc, log_rmesh)): + for atomtype, (bxc_atomtype, logmesh_atomtype, meshpoints) in enumerate(zip(bxc, log_rmesh, rmesh_dimensions)): rmesh = np.arange(0, max(logmesh_atomtype), max(logmesh_atomtype) / radial_mesh_points) for lrep in range(4): for lpot in range(2 * lrep + 1): for mpot in range(-lpot, lpot + 1): lm = lpot * (lpot + 1) + mpot - bxc_interpolated = interp1d(logmesh_atomtype, bxc_atomtype[lm, :], fill_value='extrapolate') + bxc_interpolated = interp1d(logmesh_atomtype[:meshpoints], bxc_atomtype[lm, :meshpoints], fill_value='extrapolate') bxc_integrated = np.trapz(bxc_interpolated(rmesh), rmesh) for m in range(2 * lrep + 1): diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 0588352e4..3f2afa482 100644 --- a/masci_tools/tools/greensfunction.py +++ b/masci_tools/tools/greensfunction.py @@ -474,6 +474,7 @@ def _read_gf_element(file: Any, index: int) -> tuple[GreensfElement, dict[str, A recipe = _get_radial_recipe(group_name, index, gf_element.contour, nLO=gf_element.nLO, version=version) data, attributes = h5reader.read(recipe=recipe) + attributes['file_version'] = version return gf_element, data, attributes From 3121140f898f8c77c7924055a4eb5590973b27bb Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Wed, 26 Apr 2023 09:50:36 +0200 Subject: [PATCH 31/32] formatting --- masci_tools/tools/greensf_calculations.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index b881100f3..3497accc9 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -341,7 +341,9 @@ def calculate_bxc_mmp_matrix(file: FileLike, radial_mesh_points: int = 4000, cut for lpot in range(2 * lrep + 1): for mpot in range(-lpot, lpot + 1): lm = lpot * (lpot + 1) + mpot - bxc_interpolated = interp1d(logmesh_atomtype[:meshpoints], bxc_atomtype[lm, :meshpoints], fill_value='extrapolate') + bxc_interpolated = interp1d(logmesh_atomtype[:meshpoints], + bxc_atomtype[lm, :meshpoints], + fill_value='extrapolate') bxc_integrated = np.trapz(bxc_interpolated(rmesh), rmesh) for m in range(2 * lrep + 1): From 23f9191edf33f88ba50644d227500c2ed7d848af Mon Sep 17 00:00:00 2001 From: janssenhenning Date: Tue, 16 May 2023 12:30:43 +0200 Subject: [PATCH 32/32] Fix einsum --- masci_tools/tools/greensf_calculations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 3497accc9..ba3d6a692 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -169,7 +169,7 @@ def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[Greens delta_j = onsite_delta[g1.atomTypep - 1, g1.l, 3 - g1.l:4 + g1.l, 3 - g1.l:4 + g1.l] gdeltaij = np.einsum('ij,zjk...->zik...', delta_i, gij) - gdeltaji = np.einsum('ij,zjk...->zik:::', delta_j, gji) + gdeltaji = np.einsum('ij,zjk...->zik...', delta_j, gji) else: gdeltaij = onsite_delta[g1.atomType - 1, g1.l] * gij gdeltaji = onsite_delta[g1.atomTypep - 1, g1.l] * gji