From 3c500960dbd72b9a06009e0e92c9db9b515e6ecf Mon Sep 17 00:00:00 2001 From: antelmor Date: Thu, 14 Dec 2023 20:23:08 -0500 Subject: [PATCH 1/9] Merge process for noncollinear systems --- TB2J/io_merge.py | 484 ++++++++++-------------------------------- scripts/TB2J_merge.py | 11 +- 2 files changed, 121 insertions(+), 374 deletions(-) diff --git a/TB2J/io_merge.py b/TB2J/io_merge.py index 13fb3ef..f5e5e9b 100644 --- a/TB2J/io_merge.py +++ b/TB2J/io_merge.py @@ -1,435 +1,175 @@ import os import copy import numpy as np -from scipy.spatial.transform import Rotation +from itertools import product from TB2J.io_exchange import SpinIO -from TB2J.tensor_rotate import remove_components -from TB2J.Jtensor import DMI_to_Jtensor, Jtensor_to_DMI -from TB2J.tensor_rotate import Rzx, Rzy, Rzz, Ryz, Rxz +uz = np.array([[0.0, 0.0, 1.0]]) -def test_rotation_matrix(): - x = [1, 0, 0] - y = [0, 1, 0] - z = [0, 0, 1] - print(Rxz.apply(x)) - print(Ryz.apply(y)) - print(Rzx.apply(z)) - print(Rzy.apply(z)) +def get_rotation_matrix(magmoms): + dim = magmoms.shape[0] + v = magmoms / np.linalg.norm(magmoms, axis=-1).reshape(-1, 1) + n = v[:, [1, 0, 2]] + n[:, 0] *= -1 + n[:, -1] *= 0 + n /= np.linalg.norm(n, axis=-1).reshape(dim, 1) + z = np.repeat(uz, dim, axis=0) + A = np.stack([z, np.cross(n, z), n], axis=1) + B = np.stack([v, np.cross(n, v), n], axis=1) + R = np.einsum('nki,nkj->nij', A, B) -def recover_DMI_from_rotated_structure(Ddict, rotation): - """ - Recover the DMI vector from the rotated structure. - D: the dictionary of DMI vector in the rotated structure. - rotation: the rotation operator from the original structure to the rotated structure. - """ - for key, val in Ddict.items(): - Ddict[key] = rotation.apply(val, inverse=True) - return Ddict + Rnan = np.isnan(R) + if Rnan.any(): + nanidx = np.where(Rnan)[0] + R[nanidx] = I + R[nanidx, 2] = v[nanidx] + return R -def recover_Jani_fom_rotated_structure(Janidict, rotation): - """ - Recover the Jani tensor from the rotated structure. - Janidict: the dictionary of Jani tensor in the rotated structure. - rotation: the from the original structure to the rotated structure. - """ - R = rotation.as_matrix() - RT = R.T - for key, Jani in Janidict.items(): - # Note: E=Si J Sj , Si'= Si RT, Sj' = R Sj, - # Si' J' Sj' = Si RT R J RT R Sj => J' = R J RT - # But here we are doing the opposite rotation back to - # the original axis, so we replace R with RT. - Janidict[key] = RT @ Jani @ R - return Janidict +def transform_Jani(Jani, Ri, Rj): + new_Jani = Ri @ Jani @ Rj.T + new_Jani[:, 2] *= 0 + new_Jani[2, :] *= 0 + new_Jani = Ri.T @ new_Jani @ Rj -def recover_spinat_from_rotated_structure(spinat, rotation): - """ - Recover the spinat from the rotated structure. - spinat: the spinat in the rotated structure. - rotation: the rotation operator from the original structure to the rotated structure. - """ - for i, spin in enumerate(spinat): - spinat[i] = rotation.apply(spin, inverse=True) - return spinat + new_Jani = Rj @ new_Jani @ Ri.T + new_Jani[:, 2] *= 0 + new_Jani[2, :] *= 0 + new_Jani = Rj.T @ new_Jani @ Ri + return new_Jani.round(8) -# test_rotation_matrix() +def transform_DMI(DMI, Ri, Rj): -# R_xyz = [Rxz.as_matrix(), Ryz.as_matrix(), np.eye(3, dtype=float)] + new_DMI = Ri @ DMI + new_DMI[2] *= 0 + new_DMI = Ri.T @ new_DMI + new_DMI = Rj @ new_DMI + new_DMI[2] *= 0 + new_DMI = Rj.T @ new_DMI -def rot_merge_DMI(Dx, Dy, Dz): - Dx_z = Rzx.apply(Dx) - Dy_z = Rzy.apply(Dy) - D = ( - np.array([0.0, 0.5, 0.5]) * Dx_z - + np.array([0.5, 0.0, 0.5]) * Dy_z - + np.array([0.5, 0.5, 0.0]) * Dz - ) - return D + return new_DMI.round(8) +class SpinIO_merge(SpinIO): + def __init__(self, *args, **kwargs): + super(SpinIO_merge, self).__init__(*args, **kwargs) + self.projv = None + self._set_ind_atoms() -def rot_merge_DMI2(Dx, Dy, Dz): - Dx_z = Rzx.apply(Dx) - Dy_z = Rzy.apply(Dy) - D = ( - np.array([1, 0, 0]) * Dx_z - + np.array([0, 1, 0]) * Dy_z - + np.array([0, 0, 1]) * Dz - ) - return D + def _set_ind_atoms(self): + if not hasattr(self, 'nspin'): + self.nspin = len([i for i in self.index_spin if i >= 0]) + if not self.ind_atoms: + self.ind_atoms = dict([(i, i) for i in range(self.nspin)]) -def merge_DMI(Dx, Dy, Dz): - D = ( - np.array([0.0, 0.5, 0.5]) * Dx - + np.array([0.5, 0.0, 0.5]) * Dy - + np.array([0.5, 0.5, 0.0]) * Dz - ) - return D - - -def merge_DMI2(Dx, Dy, Dz): - D = np.array([1, 0, 0]) * Dx + np.array([0, 1, 0]) * Dy + np.array([0, 0, 1]) * Dz - return D - - -def swap_direction(m, idirections): - """ - swap two directions of a tensor m. - idirections: the index of the two directions. - """ - idirections = list(idirections) - inv = [idirections[1], idirections[0]] - n = np.copy(m) - n[:, idirections] = n[:, inv] - n[idirections, :] = n[inv, :] - return n + def _set_rotation_matrix(self, reference_cell=None): + if reference_cell is None: + self.sR = np.eye(3) + else: + self.sR = np.linalg.solve(reference_cell, self.atoms.cell.array) -def test_swap(): - m = np.reshape(np.arange(9), (3, 3)) - print(m) - print(swap_direction(m, [0, 1])) + self._set_ind_atoms() + spinat = self.spinat + ind_atoms = self.ind_atoms + + magmoms = spinat[list(ind_atoms.values())] + magmoms = magmoms @ self.sR + self.R = get_rotation_matrix(magmoms) -def merge_Jani(Janix, Janiy, Janiz): - # This is wrong. - # Jani = ( - # np.array([[0, 0, 0], [0, 1, 1], [0, 1, 1]]) * Janix - # + np.array([[1, 0, 1], [0, 0, 0], [1, 0, 1]]) * Janiy - # + np.array([[1, 1, 0], [1, 1, 0], [0, 0, 0]]) * Janiz - # ) / 2.0 - wx = np.array([[0, 0, 0], [0, 1, 1], [0, 1, 1]]) - wy = np.array([[1, 0, 1], [0, 0, 0], [1, 0, 1]]) - wz = np.array([[1, 1, 0], [1, 1, 0], [0, 0, 0]]) - Jani = (wx * Janix + wy * Janiy + wz * Janiz) / (wx + wy + wz) - return Jani + @classmethod + def load_pickle(cls, path='TB2J_results', fname='TB2J.pickle', reference_cell=None): + obj = super(SpinIO_merge, cls).load_pickle(path=path, fname=fname) + obj._set_rotation_matrix(reference_cell=reference_cell) + return obj -def read_pickle(path): - p1 = os.path.join(path, "TB2J_results", "TB2J.pickle") - p2 = os.path.join(path, "TB2J.pickle") +def read_pickle(path, reference_cell=None): + p1 = os.path.join(path, 'TB2J_results', 'TB2J.pickle') + p2 = os.path.join(path, 'TB2J.pickle') if os.path.exists(p1) and os.path.exists(p2): print(f" WARNING!: Both file {p1} and {p2} exist. Use default {p1}.") if os.path.exists(p1): - ret = SpinIO.load_pickle(os.path.join(path, "TB2J_results")) + ret = SpinIO_merge.load_pickle(os.path.join(path, 'TB2J_results'), reference_cell=reference_cell) elif os.path.exists(p2): - ret = SpinIO.load_pickle(path) + ret = SpinIO_merge.load_pickle(path, reference_cell=reference_cell) else: raise FileNotFoundError(f"Cannot find either file {p1} or {p2}") return ret - -class Merger2: - def __init__(self, paths, method): +class Merger(): + def __init__(self, *paths, main_path=None, method='structure'): + if method not in ['structure', 'spin']: + raise ValueError(f"Unrecognized method '{method}'. Available options are: 'structure' or 'spin'.") self.method = method - if method.lower() == "spin": - self.load_with_rotated_spin(paths) - elif method.lower() == "structure": - self.load_with_rotated_structure(paths) - else: - raise ValueError("method should be either 'spin' or 'structure'") - - def load_with_rotated_structure(self, paths): - """ - Merge TB2J results from multiple calculations. - :param paths: a list of paths to the TB2J results. - :param method: 'structure' or 'spin' - """ - self.paths = paths - if len(self.paths) != 3: - raise ValueError( - "The number of paths should be 3, with structure rotated from z to x, y, z" - ) - for i, path in enumerate(self.paths): - read_pickle(path) - self.indata = [read_pickle(path) for path in paths] - - self.dat = copy.deepcopy(self.indata[-1]) - # self.dat.description += ( - # "Merged from TB2J results in paths: \n " + "\n ".join(paths) + "\n" - # ) - Rotations = [Rzx, Rzy, Rzz] - for dat, rotation in zip(self.indata, Rotations): - dat.spinat = recover_spinat_from_rotated_structure(dat.spinat, rotation) - dat.dmi_ddict = recover_DMI_from_rotated_structure(dat.dmi_ddict, rotation) - dat.Jani_dict = recover_Jani_fom_rotated_structure(dat.Jani_dict, rotation) - - def load_with_rotated_spin(self, paths): - """ - Merge TB2J results from multiple calculations. - :param paths: a list of paths to the TB2J results. - :param method: 'structure' or 'spin' - """ - self.paths = paths - self.indata = [read_pickle(path) for path in paths] - self.dat = copy.deepcopy(self.indata[-1]) - # self.dat.description += ( - # "Merged from TB2J results in paths: \n " + "\n ".join(paths) + "\n" - # ) - - def merge_Jani(self): - """ - Merge the anisotropic exchange tensor. - """ - Jani_dict = {} - for key, Jani in self.dat.Jani_dict.items(): - R, i, j = key - weights = np.zeros((3, 3), dtype=float) - Jani_sum = np.zeros((3, 3), dtype=float) - for dat in self.indata: - Si = dat.get_spin_ispin(i) - Sj = dat.get_spin_ispin(j) - # print(f"{Si=}, {Sj=}") - Jani = dat.get_Jani(i, j, R, default=np.zeros((3, 3), dtype=float)) - Jani_removed, w = remove_components( - Jani, - Si, - Sj, - remove_indices=[[0, 2], [1, 2], [2, 2], [2, 1], [2, 0]], - ) - w = Jani_removed / Jani - Jani_sum += Jani * w # Jani_removed - print(f"{Jani* w=}") - weights += w - print(f"{weights=}") - if np.any(weights == 0): - raise RuntimeError( - "The data set to be merged does not give a complete anisotropic J tensor, please add more data" - ) - Jani_dict[key] = Jani_sum / weights - self.dat.Jani_dict = Jani_dict - - def merge_DMI(self): - """ - merge the DMI vector - """ - DMI = {} - for key, D in self.dat.dmi_ddict.items(): - R, i, j = key - weights = np.zeros((3, 3), dtype=float) - Dtensor_sum = np.zeros((3, 3), dtype=float) - for dat in self.indata: - Si = dat.get_spin_ispin(i) - Sj = dat.get_spin_ispin(j) - D = dat.get_DMI(i, j, R, default=np.zeros((3,), dtype=float)) - Dtensor = DMI_to_Jtensor(D) - Dtensor_removed, w = remove_components( - Dtensor, Si, Sj, remove_indices=[[0, 1], [1, 0]] - ) - Dtensor_sum += Dtensor * w # Dtensor_removed - weights += w - if np.any(weights == 0): - raise RuntimeError( - "The data set to be merged does not give a complete DMI vector, please add more data" - ) - DMI[key] = Jtensor_to_DMI(Dtensor_sum / weights) - self.dat.dmi_ddict = DMI - - def merge_Jiso(self): - """ - merge the isotropic exchange - """ - Jiso = {} - for key, J in self.dat.exchange_Jdict.items(): - R, i, j = key - weights = 0.0 - Jiso_sum = 0.0 - for dat in self.indata: - Si = dat.get_spin_ispin(i) - Sj = dat.get_spin_ispin(j) - J = dat.get_Jiso(i, j, R, default=0.0) - Jiso_sum += J # *np.eye(3, dtype=float) - weights += 1.0 - if np.any(weights == 0): - raise RuntimeError( - "The data set to be merged does not give a complete isotropic exchange, please add more data" - ) - Jiso[key] = Jiso_sum / weights - self.dat.exchange_Jdict = Jiso - - def write(self, path="TB2J_results"): - """ - Write the merged TB2J results to a folder. - :param path: the path to the folder to write the results. - """ - self.dat.description += ( - "Merged from TB2J results in paths: \n " + "\n ".join(self.paths) + "\n" - ) - if self.method == "spin": - self.dat.description += ( - ", which are from DFT data with various spin orientations. \n" - ) - elif self.method == "structure": - self.dat.description += ", which are from DFT data with structure with z axis rotated to x, y, z\n" - self.dat.description += "\n" - self.dat.write_all(path=path) + if main_path is None: + self.main_dat = read_pickle(paths[-1]) + else: + self.main_dat = read_pickle(main_path) -class Merger: - def __init__(self, path_x, path_y, path_z, method="structure"): - assert method in ["structure", "spin"] - self.dat_x = read_pickle(path_x) - self.dat_y = read_pickle(path_y) - self.dat_z = read_pickle(path_z) - self.dat = copy.copy(self.dat_z) - self.paths = [path_x, path_y, path_z] - self.method = method + self.dat = [read_pickle(path, reference_cell=self.main_dat.atoms.cell.array) for path in paths] + if main_path is not None: + self.dat.append(self.main_dat) def merge_Jani(self): Jani_dict = {} - Janixdict = self.dat_x.Jani_dict - Janiydict = self.dat_y.Jani_dict - Janizdict = self.dat_z.Jani_dict - for key, Janiz in Janizdict.items(): + for key in self.main_dat.Jani_dict.keys(): try: - R, i, j = key - keyx = R - keyy = R - Janix = Janixdict[(tuple(keyx), i, j)] - Janiy = Janiydict[(tuple(keyy), i, j)] + _, i, j = key + Jani_list = [] + for obj in self.dat: + Jani = obj.sR @ obj.Jani_dict[key] @ obj.sR.T + Jani_list.append(transform_Jani(Jani, obj.R[i], obj.R[j])) + Jani_list = np.stack(Jani_list) except KeyError as err: raise KeyError( "Can not find key: %s, Please make sure the three calculations use the same k-mesh and same Rcut." - % err - ) - if self.method == "spin": - Jani_dict[key] = merge_Jani(Janix, Janiy, Janiz) - else: - Jani_dict[key] = merge_Jani( - swap_direction(Janix, (0, 2)), swap_direction(Janiy, (1, 2)), Janiz - ) - self.dat.Jani_dict = Jani_dict + % err) + Jani_dict[key] = np.nan_to_num(np.mean(Jani_list, axis=0, where=Jani_list != 0.0)) + self.main_dat.Jani_dict = Jani_dict def merge_Jiso(self): - Jdict = {} - Jxdict = self.dat_x.exchange_Jdict - Jydict = self.dat_y.exchange_Jdict - Jzdict = self.dat_z.exchange_Jdict - for key, J in Jzdict.items(): + Jdict={} + for key in self.main_dat.exchange_Jdict.keys(): try: - Jx = Jxdict[key] - Jy = Jydict[key] - Jz = Jzdict[key] + J = np.mean([obj.exchange_Jdict[key] for obj in self.dat]) except KeyError as err: raise KeyError( - "Can not find key: %s, Please make sure the three calculations use the same k-mesh and same Rcut." - % err - ) - Jdict[key] = (Jx + Jy + Jz) / 3.0 - self.dat.exchange_Jdict = Jdict - - def merge_DMI(self): - dmi_ddict = {} - if self.dat_x.has_dmi and self.dat_y.has_dmi and self.dat_z.has_dmi: - Dxdict = self.dat_x.dmi_ddict - Dydict = self.dat_y.dmi_ddict - Dzdict = self.dat_z.dmi_ddict - for key, Dz in Dzdict.items(): - try: - R, i, j = key - keyx = R - keyy = R - Dx = Dxdict[(tuple(keyx), i, j)] - Dy = Dydict[(tuple(keyy), i, j)] - except KeyError as err: - raise KeyError( "Can not find key: %s, Please make sure the three calculations use the same k-mesh and same Rcut." - % err - ) - if self.method == "structure": - dmi_ddict[key] = rot_merge_DMI(Dx, Dy, Dz) - else: - dmi_ddict[key] = merge_DMI(Dx, Dy, Dz) - self.dat.dmi_ddict = dmi_ddict + % err) + Jdict[key] = J + self.main_dat.exchange_Jdict = Jdict + def merge_DMI(self): dmi_ddict = {} - try: - Dxdict = self.dat_x.debug_dict["DMI2"] - Dydict = self.dat_y.debug_dict["DMI2"] - Dzdict = self.dat_z.debug_dict["DMI2"] - for key, Dz in Dzdict.items(): + if all(obj.has_dmi for obj in self.dat): + for key in self.main_dat.dmi_ddict.keys(): try: - R, i, j = key - keyx = R - keyy = R - Dx = Dxdict[(tuple(keyx), i, j)] - Dy = Dydict[(tuple(keyy), i, j)] + _, i, j = key + DMI_list = [] + for obj in self.dat: + DMI = obj.sR @ obj.dmi_ddict[key] + DMI_list.append(transform_DMI(DMI, obj.R[i], obj.R[j])) + DMI_list = np.stack(DMI_list) except KeyError as err: raise KeyError( "Can not find key: %s, Please make sure the three calculations use the same k-mesh and same Rcut." - % err - ) - if self.method == "structure": - dmi_ddict[key] = rot_merge_DMI2(Dx, Dy, Dz) - elif self.method == "spin": - dmi_ddict[key] = merge_DMI2(Dx, Dy, Dz) - - self.dat.debug_dict["DMI2"] = dmi_ddict - except: - pass - - def write(self, path="TB2J_results"): - self.dat.description += ( - "Merged from TB2J results in paths: \n " + "\n ".join(self.paths) + "\n" - ) - if self.method == "spin": - self.dat.description += ( - ", which are from DFT data with spin along x, y, z orientation\n" - ) - elif self.method == "structure": - self.dat.description += ", which are from DFT data with structure with z axis rotated to x, y, z\n" - self.dat.description += "\n" - self.dat.write_all(path=path) - - -def merge(path_x, path_y, path_z, method, save=True, path="TB2J_results"): - m = Merger(path_x, path_y, path_z, method) - m.merge_Jiso() - m.merge_DMI() - m.merge_Jani() - if save: - m.write(path=path) - return m.dat - + % err) + DMI = np.nan_to_num(np.mean(DMI_list, axis=0, where=DMI_list != 0.0)) + dmi_ddict[key] = DMI + self.main_dat.dmi_ddict = dmi_ddict -def merge2(paths, method, save=True, path="TB2J_results"): - """ - Merge TB2J results from multiple calculations. - :param paths: a list of paths to the TB2J results. - :param method: 'structure' or 'spin' - :param save: whether to save the merged results. - :param path: the path to the folder to write the results. - """ - m = Merger2(paths, method) +def merge(*paths, main_path=None, method='structure', save=True, write_path='TB2J_results'): + m = Merger(*paths, main_path=main_path, method=method) m.merge_Jiso() m.merge_DMI() m.merge_Jani() if save: - m.write(path=path) + m.main_dat.write_all(path=write_path) return m.dat diff --git a/scripts/TB2J_merge.py b/scripts/TB2J_merge.py index d5ed8e9..38384fb 100755 --- a/scripts/TB2J_merge.py +++ b/scripts/TB2J_merge.py @@ -2,7 +2,7 @@ import argparse import os import sys -from TB2J.io_merge import merge, merge2 +from TB2J.io_merge import merge def main(): @@ -28,11 +28,18 @@ def main(): type=str, default="TB2J_results", ) + parser.add_argument( + "--main_path", + help="The path containning the reference structure.", + type=str, + default=None + ) args = parser.parse_args() # merge(*(args.directories), args.type.strip().lower(), path=args.output_path) # merge(*(args.directories), method=args.type.strip().lower(), path=args.output_path) - merge2(args.directories, args.type.strip().lower(), path=args.output_path) + #merge2(args.directories, args.type.strip().lower(), path=args.output_path) + merge(*args.directories, main_path=args.main_path, write_path=args.output_path) main() From d6a8bcd750fa41f5f7f0fa47c6fb1053495ab081 Mon Sep 17 00:00:00 2001 From: antelmor Date: Tue, 30 Jan 2024 12:55:38 -0500 Subject: [PATCH 2/9] Added definition of identity matrix --- TB2J/io_merge.py | 1 + 1 file changed, 1 insertion(+) diff --git a/TB2J/io_merge.py b/TB2J/io_merge.py index f5e5e9b..6de6256 100644 --- a/TB2J/io_merge.py +++ b/TB2J/io_merge.py @@ -4,6 +4,7 @@ from itertools import product from TB2J.io_exchange import SpinIO +I = np.eye(3) uz = np.array([[0.0, 0.0, 1.0]]) def get_rotation_matrix(magmoms): From 00e3d7db6d9d2d8d08a32603d4723d078e8ee9cf Mon Sep 17 00:00:00 2001 From: antelmor Date: Tue, 30 Jan 2024 14:28:51 -0500 Subject: [PATCH 3/9] Removed itertools line --- TB2J/io_merge.py | 1 - 1 file changed, 1 deletion(-) diff --git a/TB2J/io_merge.py b/TB2J/io_merge.py index 6de6256..c521353 100644 --- a/TB2J/io_merge.py +++ b/TB2J/io_merge.py @@ -1,7 +1,6 @@ import os import copy import numpy as np -from itertools import product from TB2J.io_exchange import SpinIO I = np.eye(3) From 7594cec2b1b1bad7e772ebb100e1fb32a8942fd9 Mon Sep 17 00:00:00 2001 From: antelmor Date: Tue, 19 Mar 2024 16:16:21 -0400 Subject: [PATCH 4/9] Non-collinear merging algorithm --- TB2J/io_merge.py | 175 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 175 insertions(+) create mode 100644 TB2J/io_merge.py diff --git a/TB2J/io_merge.py b/TB2J/io_merge.py new file mode 100644 index 0000000..498d080 --- /dev/null +++ b/TB2J/io_merge.py @@ -0,0 +1,175 @@ +import os +import copy +import numpy as np +from itertools import combinations_with_replacement, product +from TB2J.io_exchange import SpinIO + +u0 = np.zeros(3) +uy = np.array([0.0, 1.0, 0.0]) +uz = np.array([0.0, 0.0, 1.0]) + +def get_Jani_coefficients(a, R=np.eye(3)): + + if len(a) == 1: + u = a + v = a + else: + u = a[[0, 0, 1]] + v = a[[0, 1, 1]] + + ur = u @ R.T + vr = v @ R.T + coefficients = np.hstack([ur*vr, np.roll(ur, -1, axis=-1)*vr + np.roll(vr, -1, axis=-1)*ur]) + + return coefficients, u, v + +def get_projections(a, b, tol=1e-2): + + projections = np.empty((2, 3)) + if np.linalg.matrix_rank([a, b], tol=tol) == 1: + if np.linalg.matrix_rank([a, uy], tol=tol) == 1: + projections[0] = np.cross(a, uz) + else: + projections[0] = np.cross(a, uy) + projections[1] = np.cross(a, projections[0]) + projections /= np.linalg.norm(projections, axis=-1).reshape(-1, 1) + else: + projections[0] = np.cross(a, b) + projections[0] /= np.linalg.norm(projections[0]) + projections[1] = u0 + + return projections + +class SpinIO_merge(SpinIO): + def __init__(self, *args, **kwargs): + super(SpinIO_merge, self).__init__(*args, **kwargs) + self.projv = None + + def _set_projection_vectors(self): + + spinat = self.spinat + idx = [self.ind_atoms[i] for i in self.index_spin if i >= 0] + projv = {} + for i, j in combinations_with_replacement(range(self.nspin), 2): + a, b = spinat[idx][[i, j]] + projv[i, j] = get_projections(a, b) + projv[j, i] = projv[i, j] + + self.projv = projv + + @classmethod + def load_pickle(cls, path='TB2J_results', fname='TB2J.pickle'): + obj = super(SpinIO_merge, cls).load_pickle(path=path, fname=fname) + obj._set_projection_vectors() + + return obj + +def read_pickle(path): + p1 = os.path.join(path, 'TB2J_results', 'TB2J.pickle') + p2 = os.path.join(path, 'TB2J.pickle') + if os.path.exists(p1) and os.path.exists(p2): + print(f" WARNING!: Both file {p1} and {p2} exist. Use default {p1}.") + if os.path.exists(p1): + ret = SpinIO_merge.load_pickle(os.path.join(path, 'TB2J_results')) + elif os.path.exists(p2): + ret = SpinIO_merge.load_pickle(path) + else: + raise FileNotFoundError(f"Cannot find either file {p1} or {p2}") + return ret + +class Merger(): + def __init__(self, *paths, main_path=None): + self.dat = [read_pickle(path) for path in paths] + + if main_path is None: + self.main_dat = copy.deepcopy(self.dat[-1]) + else: + self.main_dat = read_pickle(main_path) + self.dat.append(copy.deepcopy(self.main_dat)) + + self._set_projv() + + def _set_projv(self): + + cell = self.main_dat.atoms.cell.array + rotated_cells = np.stack( + [obj.atoms.cell.array for obj in self.dat], axis=0 + ) + R = np.linalg.solve(cell, rotated_cells) + indices = range(len(self.dat)) + + proju = {}; projv = {}; coeff_matrix = {}; projectors = {}; + for key in self.main_dat.projv.keys(): + vectors = [obj.projv[key] for obj in self.dat] + coefficients, u, v = zip(*[get_Jani_coefficients(vectors[i], R=R[i]) for i in indices]) + projectors[key] = np.vstack([u[i] @ R[i].T for i in indices]) + coeff_matrix[key] = np.vstack(coefficients) + proju[key] = np.stack(u) + projv[key] = np.stack(v) + + self.proju = proju + self.projv = projv + self.coeff_matrix = coeff_matrix + self.projectors = projectors + + def merge_Jani(self): + Jani_dict = {} + proju = self.proju; projv = self.projv; coeff_matrix = self.coeff_matrix; + for key in self.main_dat.Jani_dict.keys(): + try: + R, i, j = key + u = proju[i, j] + v = projv[i, j] + Jani = np.stack([sio.Jani_dict[key] for sio in self.dat]) + projections = np.einsum('nmi,nij,nmj->nm', u, Jani, v).flatten() + except KeyError as err: + raise KeyError( + "Can not find key: %s, Please make sure the three calculations use the same k-mesh and same Rcut." + % err) + newJani = np.linalg.lstsq(coeff_matrix[i, j], projections, rcond=4e-1)[0] + Jani_dict[key] = np.array([ + [newJani[0], newJani[3], newJani[5]], + [newJani[3], newJani[1], newJani[4]], + [newJani[5], newJani[4], newJani[2]] + ]) + self.main_dat.Jani_dict = Jani_dict + + def merge_Jiso(self): + Jdict={} + for key in self.main_dat.exchange_Jdict.keys(): + try: + J = np.mean([obj.exchange_Jdict[key] for obj in self.dat]) + except KeyError as err: + raise KeyError( + "Can not find key: %s, Please make sure the three calculations use the same k-mesh and same Rcut." + % err) + Jdict[key] = J + self.main_dat.exchange_Jdict = Jdict + + + def merge_DMI(self): + dmi_ddict = {} + if all(obj.has_dmi for obj in self.dat): + projectors = self.projectors; proju = self.proju; + for key in self.main_dat.dmi_ddict.keys(): + try: + R, i, j = key + u = proju[i, j] + DMI = np.stack([sio.dmi_ddict[key] for sio in self.dat]) + projections = np.einsum('nmi,ni->nm', u, DMI).flatten() + except KeyError as err: + raise KeyError( + "Can not find key: %s, Please make sure the three calculations use the same k-mesh and same Rcut." + % err) + newDMI = np.linalg.lstsq(projectors[i, j], projections, rcond=4e-1)[0] + dmi_ddict[key] = newDMI + self.main_dat.dmi_ddict = dmi_ddict + +def merge(*paths, main_path=None, save=True, write_path='TB2J_results'): + m = Merger(*paths, main_path=main_path) + m.merge_Jiso() + m.merge_DMI() + m.merge_Jani() + if save: + m.main_dat.write_all(path=write_path) + return m.dat From 2ee97604500cfbd71384d585659f1c3a691d2b7d Mon Sep 17 00:00:00 2001 From: antelmor Date: Wed, 10 Apr 2024 19:44:49 -0400 Subject: [PATCH 5/9] Warning message for underdetermined equations --- TB2J/io_merge.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/TB2J/io_merge.py b/TB2J/io_merge.py index 498d080..c028b1c 100644 --- a/TB2J/io_merge.py +++ b/TB2J/io_merge.py @@ -1,5 +1,6 @@ import os import copy +import warnings import numpy as np from itertools import combinations_with_replacement, product from TB2J.io_exchange import SpinIO @@ -106,6 +107,14 @@ def _set_projv(self): coeff_matrix[key] = np.vstack(coefficients) proju[key] = np.stack(u) projv[key] = np.stack(v) + if np.linalg.matrix_rank(coeff_matrix[key], tol=1e-2) < 6: + warnings.warn(''' + WARNING: The matrix of equations to reconstruct the exchange tensors is + close to being singular. This happens when the magnetic moments between + different configurations are cloes to being parallel. You need to consider + more rotated spin configurations, otherwise the results might have a large + error.''' + ) self.proju = proju self.projv = projv @@ -126,7 +135,7 @@ def merge_Jani(self): raise KeyError( "Can not find key: %s, Please make sure the three calculations use the same k-mesh and same Rcut." % err) - newJani = np.linalg.lstsq(coeff_matrix[i, j], projections, rcond=4e-1)[0] + newJani = np.linalg.lstsq(coeff_matrix[i, j], projections, rcond=1e-2)[0] Jani_dict[key] = np.array([ [newJani[0], newJani[3], newJani[5]], [newJani[3], newJani[1], newJani[4]], From ba77918212c12d964218795ac06d3ac426c5450c Mon Sep 17 00:00:00 2001 From: antelmor Date: Thu, 11 Apr 2024 12:13:54 -0400 Subject: [PATCH 6/9] Create rotated structures for noncollinear systems --- TB2J/rotate_atoms.py | 42 ++++++++++++++++++++++-------------------- scripts/TB2J_rotate.py | 7 ++++++- 2 files changed, 28 insertions(+), 21 deletions(-) diff --git a/TB2J/rotate_atoms.py b/TB2J/rotate_atoms.py index ba8466f..1c63285 100755 --- a/TB2J/rotate_atoms.py +++ b/TB2J/rotate_atoms.py @@ -5,19 +5,24 @@ from TB2J.tensor_rotate import Rxx, Rxy, Rxz, Ryx, Ryy, Ryz, Rzx, Rzy, Rzz -def rotate_atom_xyz(atoms): +def rotate_atom_xyz(atoms, noncollinear=False): """ - given a atoms, return: - - 'z'->'x' - - 'z'->'y' - - 'z'->'z' + given a atoms, return rotated atoms: + atoms_1, ..., atoms_n, + where we considered n diffeerent roation axes. + + When noncollinear == True, more rotated structures + will be generated. """ - atoms_x = copy.deepcopy(atoms) - atoms_x.rotate(90, "y", rotate_cell=True) - atoms_y = copy.deepcopy(atoms) - atoms_y.rotate(90, "x", rotate_cell=True) - atoms_z = atoms - return atoms_x, atoms_y, atoms_z + + rotation_axes = [(1, 0, 0), (0, 1, 0)] + if noncollinear: + rotation_axes += [(1, 1, 0), (1, 0, 1), (0, 1, 1)] + + for axis in rotation_axes: + rotated_atoms = copy.deepcopy(atoms) + rotated_atoms.rotate(90, axis, rotate_cell=True) + yield rotated_atoms def rotate_atom_spin_one_rotation(atoms, Rotation): @@ -96,18 +101,15 @@ def check_ftype(ftype): print("=" * 40) -def rotate_xyz(fname, ftype="xyz"): +def rotate_xyz(fname, ftype="xyz", noncollinear=False): check_ftype(ftype) atoms = read(fname) atoms.set_pbc(True) - atoms_x, atoms_y, atoms_z = rotate_atom_xyz(atoms) + rotated = rotate_atom_xyz(atoms, noncollinear=noncollinear) - fname_x = f"atoms_x.{ftype}" - fname_y = f"atoms_y.{ftype}" - fname_z = f"atoms_z.{ftype}" + for i, rotated_atoms in enumerate(rotated): + write(f"atoms_{i+1}.{ftype}", rotated_atoms) + write(f"atoms_0.{ftype}", atoms) - write(fname_x, atoms_x) - write(fname_y, atoms_y) - write(fname_z, atoms_z) - print(f"The output has been written to {fname_x}, {fname_y}, {fname_z}") + print(f"The output has been written to the atoms_i.{ftype} files. atoms_0.{ftype} contains the reference structure.") diff --git a/scripts/TB2J_rotate.py b/scripts/TB2J_rotate.py index a135c2b..f2162e7 100755 --- a/scripts/TB2J_rotate.py +++ b/scripts/TB2J_rotate.py @@ -14,9 +14,14 @@ def main(): default="xyz", type=str, ) + parser.add_argument( + "--noncollinear", + action="store_true", + help="If present, six different configurations will be generated. These are required for non-collinear systems." + ) args = parser.parse_args() - rotate_xyz(args.fname, ftype=args.ftype) + rotate_xyz(args.fname, ftype=args.ftype, noncollinear=args.noncollinear) if __name__ == "__main__": From 5fb95543fde6b2e2241c1a61e3300aa493e7a430 Mon Sep 17 00:00:00 2001 From: antelmor Date: Fri, 12 Apr 2024 12:46:01 -0400 Subject: [PATCH 7/9] Script to rotate the spins given a density matrix --- TB2J/rotate_siestaDM.py | 36 ++++++++++++++++++++++++++++++++++++ scripts/TB2J_rotateDM.py | 21 +++++++++++++++++++++ setup.py | 1 + 3 files changed, 58 insertions(+) create mode 100644 TB2J/rotate_siestaDM.py create mode 100644 scripts/TB2J_rotateDM.py diff --git a/TB2J/rotate_siestaDM.py b/TB2J/rotate_siestaDM.py new file mode 100644 index 0000000..f0a353a --- /dev/null +++ b/TB2J/rotate_siestaDM.py @@ -0,0 +1,36 @@ +import sisl + +def rotate_siesta_DM(DM, noncollinear=False): + + angles_list = [ [0.0, 90.0, 0.0], [0.0, 90.0, 90.0] ] + if noncollinear: + angles_list += [[0.0, 45.0, 0.0], [0.0, 90.0, 45.0], [0.0, 45.0, 90.0]] + + for angles in angles_list: + yield DM.spin_rotate(angles) + +def read_label(fdf_fname): + + label = 'siesta' + with open(fdf_fname, 'r') as File: + for line in File: + corrected_line = line.lower().replace('.', '').replace('-', '') + if 'systemlabel' in corrected_line: + label = line.split()[1] + break + + return label + +def rotate_DM(fdf_fname, noncollinear=False): + + fdf = sisl.get_sile(fdf_fname) + DM = fdf.read_density_matrix() + label = read_label(fdf_fname) + + rotated = rotate_siesta_DM(DM, noncollinear=noncollinear) + + for i, rotated_DM in enumerate(rotated): + rotated_DM.write(f"{label}_{i+1}.DM") + DM.write(f"{label}_0.DM") + + print(f"The output has been written to the {label}_i.DM files. {label}_0.DM contains the reference density matrix.") diff --git a/scripts/TB2J_rotateDM.py b/scripts/TB2J_rotateDM.py new file mode 100644 index 0000000..63d215a --- /dev/null +++ b/scripts/TB2J_rotateDM.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python3 +import argparse +from TB2J.rotate_siestaDM import rotate_DM + +def main(): + parser = argparse.ArgumentParser(description="") + parser.add_argument( + "--fdf_fname", help="Name of the *.fdf siesta file." + ) + parser.add_argument( + "--noncollinear", + action="store_true", + help="If present, six different configurations will be generated. These are required for non-collinear systems." + ) + + args = parser.parse_args() + rotate_DM(args.fdf_fname, noncollinear=args.noncollinear) + + +if __name__ == "__main__": + main() diff --git a/setup.py b/setup.py index 874d7b5..f2e1646 100644 --- a/setup.py +++ b/setup.py @@ -19,6 +19,7 @@ "scripts/siesta2J.py", "scripts/abacus2J.py", "scripts/TB2J_rotate.py", + "scripts/TB2J_rotateDM.py", "scripts/TB2J_merge.py", "scripts/TB2J_magnon.py", "scripts/TB2J_magnon_dos.py", From 59b7b1780e66ad19849bca8982bd4046982bfa25 Mon Sep 17 00:00:00 2001 From: antelmor Date: Fri, 12 Apr 2024 20:23:03 -0400 Subject: [PATCH 8/9] Added documentation for merging process --- docs/src/rotate_and_merge.rst | 30 +++++++++++++++++++----------- 1 file changed, 19 insertions(+), 11 deletions(-) diff --git a/docs/src/rotate_and_merge.rst b/docs/src/rotate_and_merge.rst index 0235bf8..a8604d3 100644 --- a/docs/src/rotate_and_merge.rst +++ b/docs/src/rotate_and_merge.rst @@ -3,32 +3,40 @@ Averaging multiple parameters =============================== -As discussed in the previous section, the :math:`z` component of the DMI, and the :math:`xz`, :math:`yz`, :math:`zz`, :math:`zx`, :math:`zy`, components of the anisotropic exchanges are non-physical, and an :math:`xyz` average is needed to get the full set of magnetic interaction parameters if the spins are all along :math:`z`. In this case, scripts to rotate the structure and merge the results are provided, they are named TB2J\_rotate.py and TB2J\_merge.py. The TB2J\_rotate.py reads the structure file and generates three files containing the :math:`z\rightarrow x`, :math:`z\rightarrow y`,and the non-rotated structures. The output files are named atoms\_x, atoms\_y, atoms\_z. A large number of output file formats is supported thanks to the ASE library and the format of the output structure files is provided using the --format parameter. An example for using the rotate file is: +When the spins of sites :math:`i` and :math:`j` are along the directions :math:`\hat{\mathbff{m}}_i` and :math:`\hat{\mathbf{m}}_j`, respectively, the components of :math:`\mathbf{J}^{ani}_{ij}` and :math:`\mathbf{D}_{ij}` along those directions will be unphysical. In other words, if :math:`\hat{\mathbf{u}}` is a unit vector orthogonal to both :math:`\hat{\mathbff{m}}_i` and :math:`\hat{\mathbf{m}}_j`, we can only obtain the projections :math:`\hat{\mathbf{u}}^T \mathbf{J}^{ani}_{ij} \hat{\mathbf{u}}` and :math:`\hat{\mathbf{u}}^T \mathbf{D}_{ij} \hat{\mathbf{u}}`. Notice that for collinear systems, there will be two orthonormal vectors :math:`\hat{\mathbf{u}}` and :math:`\hat{\mathbf{v}}` that are also orthogonal to :math:`\hat{\mathbff{m}}_i` and :math:`\hat{\mathbf{m}}_j`. + +The projection for :math:`\mathbf{J}^{ani}_{ij}` can be written as + +:math:`\hat{\mathbf{u}}^T \mathbf{J}^{ani}_{ij} \hat{\mathbf{u}} = \hat{J}_{ij}^{xx} u_x^2 + \hat{J}_{ij}^{yy} u_y^2 + \hat{J}_{ij}^{zz} u_z^2 + 2\hat{J}_{ij}^{xy} u_x u_y + 2\hat{J}_{ij}^{yz} u_y u_z + 2\hat{J}_{ij}^{zx} u_z u_x,` + +where we considered :math:`\mathbf{J}^{ani}_{ij}` to be symmetric. This equation gives us a way of reconstructing :math:`\mathbf{J}^{ani}_{ij}` by performing TB2J calculations on rotated spin configurations. If we perform six calculations such that :math:`\hat{\mathbf{u}}` lies along six different directions, we obtain six linear equations that can be solved for the six independent components of :math:`\mathbf{J}^{ani}_{ij}`. We can also reconstruct the :math:`\mathbf{D}_{ij}` tensor in a similar way. Moreover, if the system is collinear then only three different calculations are needed. + +To account for this, TB2J provides scripts to rotate the structure and merge the results; they are named TB2J\_rotate.py and TB2J\_merge.py. The TB2J\_rotate.py reads the structue file and generates three(six) files containing the rotated structures whenever the system is collinear (non-collinear). The --noncollinear parameters is used to specify wheter the system is noncollinear. The output files are named atoms\_i (i = 0, ..., 5), where atoms\_0 contains the unrotated structure. A large number of file formats is supported thanks to the ASE library and the output structure files format is provided through the --format parameter. An example for using the rotate file with a collinear system is: .. code-block:: shell TB2J_rotate.py BiFeO3.vasp --ftype vasp +If te system is noncollinear, then we run the following instead: -.. note:: - - Some file format does not include the cell orientation, e.g. the cif file. They should not be used as the output format. - +.. code-block:: shell -The user has to perform DFT single point energy calculations for these three structures in different directories, keeping the spins along the $z$ direction, and run TB2J on each of them. After producing the TB2J results for the three rotated structures, we can merge the DMI results with the following command by providing the paths to the TB2J results of the three cases: + TB2J_rotate.py BiFeO3.vasp --ftype vasp --noncollinear -:: +The user has to perform DFT single point energy calculations for the generated structures in different directories, keeping the spins along the $z$ direction, and run TB2J on each of them. After producing the TB2J results for the rotated structures, we can merge the DMI results with the following command by providing the paths to the TB2J results of the three cases: - TB2J_merge.py BiFeO3_x BiFeO3_y BiFeO3_z --type structure +:: + TB2J_merge.py BiFeO3_1 BiFeO3_2 BiFeO3_0 -Note that the whole structure are rotated w.r.t. the laboratory axis but not to the cell axis. Therefore, the k-points should not be changed in both the DFT calculation and the TB2J calculation. +Here the last directory will be taken as the reference structure. Note that the whole structure are rotated w.r.t. the laboratory axis but not to the cell axis. Therefore, the k-points should not be changed in both the DFT calculation and the TB2J calculation. A new TB2J\_results directory is then made which contains the merged final results. -Another method is to do the DFT calculation with spins along the :math:`x`, :math:`y` and :math:`z` axis, respectively, and then merge the result with: +Another method is to do the DFT calculation with spins rotated globally. That is they are rotated with respect to an axis, but their relative orientations remain the same. This can be specified in the initial magnetic moments from a DFT calculation. For calculations done with SIESTA, there is a script that rotates the density matrix file along different directions. We can then use these density matrix files to run single point calculations to obtain the required rotated magnetic configurations. An example is: :: - TB2J_merge.py BiFeO3_x BiFeO3_y BiFeO3_z --type spin + TB2J_rotateDM.py --fdf_fname /location/of/the/siesta/*.fdf/file +As in the previous case, we can use the --noncollinear parameter to generate more configurations. The merging process is performed in the same way. From 730c3e45a385c9f888fe3d5609b011c778a77f1d Mon Sep 17 00:00:00 2001 From: antelmor Date: Fri, 12 Apr 2024 20:27:32 -0400 Subject: [PATCH 9/9] Fixed merge documentation --- docs/src/rotate_and_merge.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/rotate_and_merge.rst b/docs/src/rotate_and_merge.rst index a8604d3..0660033 100644 --- a/docs/src/rotate_and_merge.rst +++ b/docs/src/rotate_and_merge.rst @@ -3,7 +3,7 @@ Averaging multiple parameters =============================== -When the spins of sites :math:`i` and :math:`j` are along the directions :math:`\hat{\mathbff{m}}_i` and :math:`\hat{\mathbf{m}}_j`, respectively, the components of :math:`\mathbf{J}^{ani}_{ij}` and :math:`\mathbf{D}_{ij}` along those directions will be unphysical. In other words, if :math:`\hat{\mathbf{u}}` is a unit vector orthogonal to both :math:`\hat{\mathbff{m}}_i` and :math:`\hat{\mathbf{m}}_j`, we can only obtain the projections :math:`\hat{\mathbf{u}}^T \mathbf{J}^{ani}_{ij} \hat{\mathbf{u}}` and :math:`\hat{\mathbf{u}}^T \mathbf{D}_{ij} \hat{\mathbf{u}}`. Notice that for collinear systems, there will be two orthonormal vectors :math:`\hat{\mathbf{u}}` and :math:`\hat{\mathbf{v}}` that are also orthogonal to :math:`\hat{\mathbff{m}}_i` and :math:`\hat{\mathbf{m}}_j`. +When the spins of sites :math:`i` and :math:`j` are along the directions :math:`\hat{\mathbf{m}}_i` and :math:`\hat{\mathbf{m}}_j`, respectively, the components of :math:`\mathbf{J}^{ani}_{ij}` and :math:`\mathbf{D}_{ij}` along those directions will be unphysical. In other words, if :math:`\hat{\mathbf{u}}` is a unit vector orthogonal to both :math:`\hat{\mathbf{m}}_i` and :math:`\hat{\mathbf{m}}_j`, we can only obtain the projections :math:`\hat{\mathbf{u}}^T \mathbf{J}^{ani}_{ij} \hat{\mathbf{u}}` and :math:`\hat{\mathbf{u}}^T \mathbf{D}_{ij} \hat{\mathbf{u}}`. Notice that for collinear systems, there will be two orthonormal vectors :math:`\hat{\mathbf{u}}` and :math:`\hat{\mathbf{v}}` that are also orthogonal to :math:`\hat{\mathbf{m}}_i` and :math:`\hat{\mathbf{m}}_j`. The projection for :math:`\mathbf{J}^{ani}_{ij}` can be written as