diff --git a/masci_tools/tools/greensf_calculations.py b/masci_tools/tools/greensf_calculations.py index 7e5d76a1c..ba3d6a692 100644 --- a/masci_tools/tools/greensf_calculations.py +++ b/masci_tools/tools/greensf_calculations.py @@ -11,13 +11,17 @@ 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 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 +from typing import Any, Callable +import h5py try: from typing import Literal @@ -31,12 +35,11 @@ # - custom decompositions (e.g. e2g/t2g) -def calculate_heisenberg_jij( - hdffileORgreensfunctions: FileLike | list[GreensFunction], - reference_atom: int, - onsite_delta: np.ndarray, - max_shells: int | 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 @@ -47,6 +50,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 """ @@ -56,6 +61,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 +81,24 @@ 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, 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('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 jij_constants['R'].append(dist) jij_constants['R_ij_x'].append(g1.atomDiff.tolist()[0]) @@ -87,8 +113,10 @@ def calculate_heisenberg_jij( def calculate_heisenberg_tensor(hdffileORgreensfunctions: FileLike | list[GreensFunction], reference_atom: int, - onsite_delta: np.ndarray, - max_shells: int | None = None) -> pd.DataFrame: + onsite_delta: np.ndarray | None = None, + max_shells: int | None = None, + 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 @@ -101,6 +129,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 """ @@ -110,18 +140,44 @@ 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: dist = round(dist, 12) + 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) - 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, 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('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) @@ -137,8 +193,9 @@ 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 + 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) @@ -252,3 +309,49 @@ 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, cutoff: float | None = None) -> 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')) + 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, 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[:meshpoints], + bxc_atomtype[lm, :meshpoints], + 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, 3 - lrep + m, 3 - lrep + mp] += gaunt_coeff * bxc_integrated + + bxc_mmp = bxc_mmp.real + if cutoff is not None: + bxc_mmp[np.abs(bxc_mmp) < cutoff] = 0 + return bxc_mmp diff --git a/masci_tools/tools/greensfunction.py b/masci_tools/tools/greensfunction.py index 9aeaf2f8b..3f2afa482 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) @@ -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 @@ -662,6 +663,36 @@ 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 + """ + + 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 + + for name in self._data.keys(): + 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 @@ -749,15 +780,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 +796,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 +832,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 +941,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 +952,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 +969,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 +981,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/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' ] 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): """