diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml index 00b5d04..3662572 100644 --- a/.github/workflows/python-app.yml +++ b/.github/workflows/python-app.yml @@ -16,9 +16,9 @@ jobs: build: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Set up Python 3.10 - uses: actions/setup-python@v3 + uses: actions/setup-python@v5 with: python-version: "3.10" - name: Install dependencies diff --git a/README.md b/README.md index e1c4da4..e3b7068 100644 --- a/README.md +++ b/README.md @@ -7,6 +7,8 @@ TB2J is a open source python package for calculating the magnetic interaction parameters in Heisenberg models from DFT. It use the magnetic force theorem and take the local rigid spin rotation as a perturbation in the Green's function method. +The TB2J project is initialized in the PhyTheMa and Nanomat teams in the University of Liege. + The features include: - Calculates parameters in Heisenberg model, including isotropic exchange, anisotropic exchange, Dyzanoshinskii-Moriya interaction. - Can use the input from many DFT codes with Wannier90, e.g. Abinit, Quantum Espresso, Siesta, VASP, etc. diff --git a/TB2J/Jdownfolder.py b/TB2J/Jdownfolder.py index 55107e9..c77631e 100644 --- a/TB2J/Jdownfolder.py +++ b/TB2J/Jdownfolder.py @@ -17,9 +17,23 @@ def ind_to_indn(ind, n=3): return indn +class JR_model: + def __init__(self, JR, Rlist): + self.JR = JR + self.Rlist = Rlist + self.nR = len(Rlist) + + def get_Jq(self, q): + Jq = np.zeros(self.JR[0].shape, dtype=complex) + for iR, R in enumerate(self.Rlist): + phase = np.exp(2.0j * np.pi * np.dot(q, R)) + Jq += self.JR[iR] * phase + return Jq + + class JDownfolder: def __init__(self, JR, Rlist, iM, iL, qmesh, iso_only=False): - self.JR = JR + self.model = JR_model(JR, Rlist) self.Rlist = Rlist self.nR = len(Rlist) self.nM = len(iM) @@ -34,25 +48,18 @@ def __init__(self, JR, Rlist, iM, iL, qmesh, iso_only=False): self.nLn = self.nL * 3 self.iso_only = iso_only - def get_Jq(self, q): - Jq = np.zeros(self.JR[0].shape, dtype=complex) - for iR, R in enumerate(self.Rlist): - phase = np.exp(2.0j * np.pi * np.dot(q, R)) - Jq += self.JR[iR] * phase - return Jq - def get_JR(self): JR_downfolded = np.zeros((self.nR, self.nMn, self.nMn), dtype=float) Jq_downfolded = np.zeros((self.nqpt, self.nMn, self.nMn), dtype=complex) self.iMn = ind_to_indn(self.iM, n=3) self.iLn = ind_to_indn(self.iL, n=3) for iq, q in enumerate(self.qpts): - Jq = self.get_Jq(q) + Jq = self.model.get_Jq(q) Jq_downfolded[iq] = self.downfold_oneq(Jq) for iR, R in enumerate(self.Rlist): phase = np.exp(-2.0j * np.pi * np.dot(q, R)) JR_downfolded[iR] += np.real(Jq_downfolded[iq] * phase / self.nqpt) - return JR_downfolded + return JR_downfolded, self.Rlist def downfold_oneq(self, J): JMM = J[np.ix_(self.iMn, self.iMn)] @@ -63,17 +70,80 @@ def downfold_oneq(self, J): return Jn +class PWFDownfolder: + def __init__(self, JR, Rlist, iM, iL, qmesh, atoms=None, iso_only=False, **kwargs): + from lawaf.interfaces.magnon.magnon_downfolder import ( + MagnonWrapper, + MagnonDownfolder, + ) + + model = MagnonWrapper(JR, Rlist, atoms) + wann = MagnonDownfolder(model) + # Downfold the band structure. + index_basis = [] + for i in iM: + index_basis += list(range(i * 3, i * 3 + 3)) + params = dict( + method="projected", + # method="maxprojected", + kmesh=qmesh, + nwann=len(index_basis), + selected_basis=index_basis, + # anchors={(0, 0, 0): (-1, -2, -3, -4)}, + # anchors={(0, 0, 0): ()}, + # use_proj=True, + enhance_Amn=2.0, + ) + params.update(kwargs) + wann.set_parameters(**params) + print("begin downfold") + ewf = wann.downfold() + ewf.save_hr_pickle("downfolded_JR.pickle") + + # Plot the band structure. + wann.plot_band_fitting( + # kvectors=np.array([[0, 0, 0], [0.5, 0, 0], + # [0.5, 0.5, 0], [0, 0, 0], + # [.5, .5, .5]]), + # knames=['$\Gamma$', 'X', 'M', '$\Gamma$', 'R'], + cell=model.atoms.cell, + supercell_matrix=None, + npoints=100, + efermi=None, + erange=None, + fullband_color="blue", + downfolded_band_color="green", + marker="o", + ax=None, + savefig="downfold_band.png", + show=True, + ) + self.JR_downfolded = ewf.HwannR + self.Rlist = ewf.Rlist + + def get_JR(self): + return self.JR_downfolded, self.Rlist + + class JDownfolder_pickle: def __init__( - self, inpath, metals, ligands, outpath, qmesh=[7, 7, 7], iso_only=False + self, + inpath, + metals, + ligands, + outpath, + qmesh=[7, 7, 7], + iso_only=False, + method="pwf", + **kwargs ): self.exc = SpinIO.load_pickle(path=inpath, fname="TB2J.pickle") self.iso_only = (self.exc.dmi_ddict is None) or iso_only - self.metals = metals self.ligands = ligands self.outpath = outpath + self.method = method # read atomic structure self.atoms = self.exc.atoms @@ -83,7 +153,8 @@ def __init__( self.Rcut = None self._build_atom_index() self._prepare_distance() - self._downfold() + Jd, Rlist = self._downfold(**kwargs) + self._Jd_to_exchange(Jd, Rlist) def _build_atom_index(self): self.magnetic_elements = self.metals @@ -101,18 +172,33 @@ def _build_atom_index(self): self.nL = len(self.iL) self.nsite = self.nM + self.nL - def _downfold(self): + def _downfold(self, **kwargs): JR2 = self.exc.get_full_Jtensor_for_Rlist(asr=True) - d = JDownfolder( - JR2, - self.exc.Rlist, - iM=self.iM, - iL=self.iL, - qmesh=self.qmesh, - iso_only=self.iso_only, - ) - Jd = d.get_JR() + if self.method == "lowdin": + d = JDownfolder( + JR2, + self.exc.Rlist, + iM=self.iM, + iL=self.iL, + qmesh=self.qmesh, + iso_only=self.iso_only, + ) + Jd, Rlist = d.get_JR() + else: + d = PWFDownfolder( + JR2, + self.exc.Rlist, + iM=self.iM, + iL=self.iL, + qmesh=self.qmesh, + atoms=self.atoms, + iso_only=self.iso_only, + **kwargs + ) + Jd, Rlist = d.get_JR() + return Jd, Rlist + def _Jd_to_exchange(self, Jd, Rlist): self._prepare_distance() self._prepare_index_spin() self.Jdict = {} @@ -123,7 +209,7 @@ def _downfold(self): self.DMIdict = {} self.Janidict = {} - for iR, R in enumerate(d.Rlist): + for iR, R in enumerate(Rlist): for i, ispin in enumerate(self.index_spin): for j, jspin in enumerate(self.index_spin): if ispin >= 0 and jspin >= 0: diff --git a/TB2J/Jtensor.py b/TB2J/Jtensor.py index 617f613..375fb7e 100644 --- a/TB2J/Jtensor.py +++ b/TB2J/Jtensor.py @@ -79,7 +79,7 @@ def combine_J_tensor(Jiso=0.0, D=np.zeros(3), Jani=np.zeros((3, 3), dtype=float) :param Jani: 3x3 matrix anisotropic exchange :returns: A 3x3 matrix, the exchange paraemter in tensor form. """ - Jtensor = np.zeros((3, 3), dtype=float) + Jtensor = np.zeros((3, 3), dtype=complex) if Jiso is not None: Jtensor += np.eye(3, dtype=float) * Jiso if Jani is not None: diff --git a/TB2J/__init__.py b/TB2J/__init__.py index 1970223..085e838 100644 --- a/TB2J/__init__.py +++ b/TB2J/__init__.py @@ -1 +1 @@ -__version__ = "0.8.2.8" +__version__ = "0.9.0.1" diff --git a/TB2J/abacus/.gitignore b/TB2J/abacus/.gitignore index 75caa95..5efca39 100644 --- a/TB2J/abacus/.gitignore +++ b/TB2J/abacus/.gitignore @@ -1 +1,2 @@ TB2J_results/ +*.png diff --git a/TB2J/abacus/MAE.py b/TB2J/abacus/MAE.py new file mode 100644 index 0000000..28506ef --- /dev/null +++ b/TB2J/abacus/MAE.py @@ -0,0 +1,320 @@ +import numpy as np +from TB2J.abacus.abacus_wrapper import AbacusWrapper, AbacusParser +from TB2J.mathutils.rotate_spin import rotate_Matrix_from_z_to_axis +from TB2J.kpoints import monkhorst_pack +from TB2J.mathutils.fermi import fermi +from TB2J.mathutils.kR_convert import R_to_k +from scipy.linalg import eigh +from copy import deepcopy +from scipy.spatial.transform import Rotation +import matplotlib.pyplot as plt +from pathlib import Path +from TB2J.abacus.occupations import Occupations + +# TODO List: +# - [x] Add the class AbacusSplitSOCWrapper +# - [x] Add the function to rotate the XC part +# - [x] Compute the band energy at arbitrary angle + + +def get_occupation(evals, kweights, nel, width=0.1): + occ = Occupations(nel=nel, width=width, wk=kweights, nspin=2) + return occ.occupy(evals) + + +def get_density_matrix(evals=None, evecs=None, kweights=None, nel=None, width=0.1): + occ = get_occupation(evals, kweights, nel, width=width) + rho = np.einsum("kib, kb, kjb -> kij", evecs, occ, evecs.conj()) + return rho + + +def spherical_to_cartesian(theta, phi, normalize=True): + """ + Convert spherical coordinates to cartesian + """ + x = np.sin(theta) * np.cos(phi) + y = np.sin(theta) * np.sin(phi) + z = np.cos(theta) + vec = np.array([x, y, z]) + if normalize: + vec = vec / np.linalg.norm(vec) + return vec + + +class AbacusSplitSOCWrapper(AbacusWrapper): + """ + Abacus wrapper with Hamiltonian split to SOC and non-SOC parts + """ + + def __init__(self, *args, **kwargs): + HR_soc = kwargs.pop("HR_soc", None) + # nbasis = HR_soc.shape[1] + # kwargs["nbasis"] = nbasis + super().__init__(*args, **kwargs) + self._HR_copy = deepcopy(self._HR) + self.HR_soc = HR_soc + self.soc_lambda = 1.0 + self.nel = 16 + self.width = 0.1 + + @property + def HR(self): + return self._HR + self.HR_soc * self.soc_lambda + + def rotate_HR_xc(self, axis): + """ + Rotate SOC part of Hamiltonian + """ + for iR, R in enumerate(self.Rlist): + self._HR[iR] = rotate_Matrix_from_z_to_axis(self._HR_copy[iR], axis) + + def rotate_Hk_xc(self, axis): + """ + Rotate SOC part of Hamiltonian + """ + for ik in range(len(self._Hk)): + self._Hk[ik] = rotate_Matrix_from_z_to_axis(self._Hk_copy[ik], axis) + + def get_density_matrix(self, kpts, kweights=None): + rho = np.zeros((len(kpts), self.nbasis, self.nbasis), dtype=complex) + evals, evecs = self.solve_all(kpts) + occ = get_occupation(evals, kweights, self.nel, width=self.width) + rho = np.einsum( + "kib, kb, kjb -> kij", evecs, occ, evecs.conj() + ) # should multiply S to the the real DM. + return rho + + def rotate_DM(self, rho, axis): + """ + Rotate the density matrix + """ + for ik in range(len(rho)): + rho[ik] = rotate_Matrix_from_z_to_axis(rho[ik], axis) + return rho + + +class RotateHam: + def __init__(self, model, kmesh, gamma=True): + self.model = model + self.kpts = monkhorst_pack(kmesh, gamma_center=gamma) + self.kweights = np.ones(len(self.kpts), dtype=float) / len(self.kpts) + + def get_band_energy2(self): + for ik, kpt in enumerate(self.kpts): + Hk, Sk = self.model.gen_ham(kpt) + evals, evecs = eigh(Hk, Sk) + rho = np.einsum( + "ib, b, jb -> ij", + evecs, + fermi(evals, self.model.efermi, width=0.05), + evecs.conj(), + ) + eband1 = np.sum(evals * fermi(evals, self.model.efermi, width=0.05)) + eband2 = np.trace(Hk @ rho) + print(eband1, eband2) + + def get_band_energy(self, dm=False): + evals, evecs = self.model.solve_all(self.kpts) + occ = get_occupation( + evals, self.kweights, self.model.nel, width=self.model.width + ) + eband = np.sum(evals * occ * self.kweights[:, np.newaxis]) + # * fermi(evals, self.model.efermi, width=0.05) + if dm: + density_matrix = self.model.get_density_matrix(evecs) + return eband, density_matrix + else: + return eband + + def calc_ref(self): + # calculate the Hk_ref, Sk_ref, Hk_soc_ref, and rho_ref + self.Sk_ref = R_to_k(self.kpts, self.model.Rlist, self.model.SR) + self.Hk_xc_ref = R_to_k(self.kpts, self.model.Rlist, self.model._HR_copy) + self.Hk_soc_ref = R_to_k(self.kpts, self.model.Rlist, self.model.HR_soc) + self.rho_ref = np.zeros( + (len(self.kpts), self.model.nbasis, self.model.nbasis), dtype=complex + ) + + evals = np.zeros((len(self.kpts), self.model.nbasis), dtype=float) + evecs = np.zeros( + (len(self.kpts), self.model.nbasis, self.model.nbasis), dtype=complex + ) + + for ik, kpt in enumerate(self.kpts): + # evals, evecs = eigh(self.Hk_xc_ref[ik]+self.Hk_soc_ref[ik], self.Sk_ref[ik]) + evals[ik], evecs[ik] = eigh(self.Hk_xc_ref[ik], self.Sk_ref[ik]) + occ = get_occupation( + evals, self.kweights, self.model.nel, width=self.model.width + ) + # occ = fermi(evals, self.model.efermi, width=self.model.width) + self.rho_ref = np.einsum("kib, kb, kjb -> kij", evecs, occ, evecs.conj()) + + def get_band_energy_from_rho(self, axis): + """ + This is wrong!! Should use second order perturbation theory to get the band energy instead. + """ + eband = 0.0 + for ik, k in enumerate(self.kpts): + rho = rotate_Matrix_from_z_to_axis(self.rho_ref[ik], axis) + Hk_xc = rotate_Matrix_from_z_to_axis(self.Hk_xc_ref[ik], axis) + Hk_soc = self.Hk_soc_ref[ik] + Htot = Hk_xc + Hk_soc * self.model.soc_lambda + # Sk = self.Sk_ref[ik] + # evals, evecs = eigh(Htot, Sk) + # rho2= np.einsum("ib, b, jb -> ij", evecs, fermi(evals, self.model.efermi, width=0.05), evecs.conj()) + if ik == 0 and False: + pass + # print(f"{evecs[:4,0:4].real=}") + # print(f"{evals[:4]=}") + # print(f"{Hk_xc[:4,0:4].real=}") + # print(f"{Htot[:4,0:4].real=}") + # print(f"{Sk[:4,0:4].real=}") + # print(f"{rho[:4,0:4].real=}") + # print(f"{rho2[:4,0:4].real=}") + # eband1 = np.sum(evals * fermi(evals, self.model.efermi, width=0.05)) + # eband2 = np.trace(Htot @ rho2).real + # eband3 = np.trace(Htot @ rho).real + # print(eband1, eband2, eband3) + e_soc = np.trace(Hk_soc @ rho) * self.kweights[ik] * self.model.soc_lambda + eband += e_soc + return eband + + def get_band_energy_vs_angles( + self, + thetas, + psis, + ): + es = [] + # es2 = [] + # e,rho = self.model.get_band_energy(dm=True) + # self.calc_ref() + # thetas = np.linspace(*angle_range, npoints) + for i, theta, phi in enumerate(zip(thetas, psis)): + axis = spherical_to_cartesian(theta, phi) + self.model.rotate_HR_xc(axis) + # self.get_band_energy2() + e = self.get_band_energy() + es.append(e) + # es2.append(e2) + return es + + +def get_model_energy(model, kmesh, gamma=True): + ham = RotateHam(model, kmesh, gamma=gamma) + return ham.get_band_energy() + + +class AbacusSplitSOCParser: + """ + Abacus parser with Hamiltonian split to SOC and non-SOC parts + """ + + def __init__(self, outpath_nosoc=None, outpath_soc=None, binary=False): + self.outpath_nosoc = outpath_nosoc + self.outpath_soc = outpath_soc + self.binary = binary + self.parser_nosoc = AbacusParser(outpath=outpath_nosoc, binary=binary) + self.parser_soc = AbacusParser(outpath=outpath_soc, binary=binary) + spin1 = self.parser_nosoc.read_spin() + spin2 = self.parser_soc.read_spin() + if spin1 != "noncollinear" or spin2 != "noncollinear": + raise ValueError("Spin should be noncollinear") + + def parse(self): + nbasis, Rlist, HR, SR = self.parser_nosoc.Read_HSR_noncollinear() + nbasis2, Rlist2, HR2, SR2 = self.parser_soc.Read_HSR_noncollinear() + # print(HR[0]) + HR_soc = HR2 - HR + model = AbacusSplitSOCWrapper(HR, SR, Rlist, nbasis, nspin=2, HR_soc=HR_soc) + model.efermi = self.parser_soc.efermi + model.basis = self.parser_nosoc.basis + model.atoms = self.parser_nosoc.atoms + return model + + +def abacus_get_MAE( + path_nosoc, path_soc, kmesh, thetas, psis, gamma=True, outfile="MAE.txt" +): + """Get MAE from Abacus with magnetic force theorem. Two calculations are needed. First we do an calculation with SOC but the soc_lambda is set to 0. Save the density. The next calculatin we start with the density from the first calculation and set the SOC prefactor to 1. With the information from the two calcualtions, we can get the band energy with magnetic moments in the direction, specified in two list, thetas, and phis.""" + parser = AbacusSplitSOCParser( + outpath_nosoc=path_nosoc, outpath_soc=path_soc, binary=False + ) + model = parser.parse() + ham = RotateHam(model, kmesh, gamma=gamma) + es = [] + for theta, psi in zip(thetas, psis): + axis = spherical_to_cartesian(theta, psi) + model.rotate_HR_xc(axis) + e = ham.get_band_energy() + es.append(ham.get_band_energy()) + if outfile: + with open(outfile, "w") as f: + f.write("theta, psi, energy\n") + for theta, psi, e in zip(thetas, psis, es): + f.write(f"{theta}, {psi}, {e}\n") + return es + + +def test_AbacusSplitSOCWrapper(): + # path = Path("~/projects/2D_Fe").expanduser() + path = Path("~/projects/TB2Jflows/examples/2D_Fe/Fe_z").expanduser() + outpath_nosoc = f"{path}/soc0/OUT.ABACUS" + outpath_soc = f"{path}/soc1/OUT.ABACUS" + parser = AbacusSplitSOCParser( + outpath_nosoc=outpath_nosoc, outpath_soc=outpath_soc, binary=False + ) + model = parser.parse() + kmesh = [6, 6, 1] + + r = RotateHam(model, kmesh) + # thetas, es = r.get_band_energy_vs_theta(angle_range=(0, np.pi*2), rotation_axis="z", initial_direction=(1,0,0), npoints=21) + thetas, es, es2 = r.get_band_energy_vs_theta( + angle_range=(0, np.pi * 2), + rotation_axis="y", + initial_direction=(0, 0, 1), + npoints=11, + ) + # print the table of thetas and es, es2 + print("theta, e, e2") + for theta, e, e2 in zip(thetas, es, es2): + print(f"{theta=}, {e=}, {e2=}") + + plt.plot(thetas / np.pi, es - es[0], marker="o") + plt.plot(thetas / np.pi, es2 - es2[0], marker=".") + plt.savefig("E_along_z_x_z.png") + plt.show() + + +def abacus_get_MAE_cli(): + import argparse + + parser = argparse.ArgumentParser( + description="Get MAE from Abacus with magnetic force theorem. Two calculations are needed. First we do an calculation with SOC but the soc_lambda is set to 0. Save the density. The next calculatin we start with the density from the first calculation and set the SOC prefactor to 1. With the information from the two calcualtions, we can get the band energy with magnetic moments in the direction, specified in two list, thetas, and phis. " + ) + parser.add_argument("path_nosoc", type=str, help="Path to the calculation with ") + parser.add_argument("path_soc", type=str, help="Path to the SOC calculation") + parser.add_argument("thetas", type=float, nargs="+", help="Thetas") + parser.add_argument("psis", type=float, nargs="+", help="Phis") + parser.add_argument("kmesh", type=int, nargs=3, help="K-mesh") + parser.add_argument( + "--gamma", action="store_true", help="Use Gamma centered kpoints" + ) + parser.add_argument( + "--outfile", + type=str, + help="The angles and the energey will be saved in this file.", + ) + args = parser.parse_args() + abacus_get_MAE( + args.path_nosoc, + args.path_soc, + args.kmesh, + args.thetas, + args.psis, + gamma=args.gamma, + outfile=args.outfile, + ) + + +if __name__ == "__main__": + abacus_get_MAE_cli() diff --git a/TB2J/abacus/abacus_wrapper.py b/TB2J/abacus/abacus_wrapper.py index f7b0171..10dde17 100644 --- a/TB2J/abacus/abacus_wrapper.py +++ b/TB2J/abacus/abacus_wrapper.py @@ -16,10 +16,10 @@ class AbacusWrapper(AbstractTB): def __init__(self, HR, SR, Rlist, nbasis, nspin=1): - self.R2kfactor = -2j * np.pi + self.R2kfactor = 2j * np.pi self.is_orthogonal = False self._name = "ABACUS" - self.HR = HR + self._HR = HR self.SR = SR self.Rlist = Rlist self.nbasis = nbasis @@ -27,6 +27,14 @@ def __init__(self, HR, SR, Rlist, nbasis, nspin=1): self.norb = nbasis * nspin self._build_Rdict() + @property + def HR(self): + return self._HR + + @HR.setter + def set_HR(self, HR): + self._HR = HR + def _build_Rdict(self): if hasattr(self, "Rdict"): pass @@ -62,6 +70,8 @@ def gen_ham(self, k, convention=2): S = self.SR[iR] * phase # Sk += S + S.conjugate().T Sk += S + # Hk = (Hk + Hk.conj().T)/2 + # Sk = (Sk + Sk.conj().T)/2 elif convention == 1: # TODO: implement the first convention (the r convention) raise NotImplementedError("convention 1 is not implemented yet.") @@ -74,6 +84,14 @@ def solve(self, k, convention=2): Hk, Sk = self.gen_ham(k, convention=convention) return eigh(Hk, Sk) + def solve_all(self, kpts, convention=2): + nk = len(kpts) + evals = np.zeros((nk, self.nbasis), dtype=float) + evecs = np.zeros((nk, self.nbasis, self.nbasis), dtype=complex) + for ik, k in enumerate(kpts): + evals[ik], evecs[ik] = self.solve(k, convention=convention) + return evals, evecs + def HSE_k(self, kpt, convention=2): H, S = self.gen_ham(tuple(kpt), convention=convention) evals, evecs = eigh(H, S) diff --git a/TB2J/abacus/occupations.py b/TB2J/abacus/occupations.py new file mode 100644 index 0000000..6509bed --- /dev/null +++ b/TB2J/abacus/occupations.py @@ -0,0 +1,278 @@ +""" +This file is stolen from the hotbit programm, with some modification. +""" + +import numpy as np +from scipy.optimize import brentq +import sys + +from ase.dft.dos import DOS +from scipy import integrate + +# import numba + +# from numba import float64, int32 + +MAX_EXP_ARGUMENT = np.log(sys.float_info.max) + +# @numba.vectorize(nopython=True) +# def myfermi(e, mu, width, nspin): +# x = (e - mu) / width +# if x < -10: +# ret = 2.0 / nspin +# elif x > 10: +# ret = 0.0 +# else: +# ret = 2.0 / nspin / (math.exp(x) + 1) +# return ret + + +def myfermi(e, mu, width, nspin): + x = (e - mu) / width + return np.where(x < 10, 2.0 / (nspin * (np.exp(x) + 1.0)), 0.0) + + +class Occupations(object): + def __init__(self, nel, width, wk, nspin=1): + """ + Initialize parameters for occupations. + :param nel: Number of electrons + :param width: Fermi-broadening + :param wk: k-point weights. eg. If only gamma, [1.0] + :param nspin(optional): number of spin, if spin=1 multiplicity=2 else, multiplicity=1. + """ + self.nel = nel + self.width = width + self.wk = wk + self.nk = len(wk) + self.nspin = nspin + + def get_mu(self): + """Return the Fermi-level (or chemical potential).""" + return self.mu + + def fermi(self, mu): + """ + Occupy states with given chemical potential. + Occupations are 0...2; without k-point weights + """ + return myfermi(self.e, mu, self.width, self.nspin) + + def root_function(self, mu): + """This function is exactly zero when mu is right.""" + f = self.fermi(mu) + return np.einsum("i, ij->", self.wk, f) - self.nel + + def occupy(self, e, xtol=1e-11): + """ + Calculate occupation numbers with given Fermi-broadening. + + @param e: e[ind_k,ind_orb] energy of k-point, state a + Note added by hexu: With spin=2,e[k,a,sigma], it also work. only the *2 should be removed. + @param wk: wk[:] weights for k-points + @param width: The Fermi-broadening + + Returns: fermi[ind_k, ind_orb] + """ + self.e = e + eflat = e.flatten() + ind = np.argsort(eflat) + e_sorted = eflat[ind] + if self.nspin == 1: + m = 2 + elif self.nspin == 2: + m = 1 + n_sorted = (self.wk[:, None, None] * np.ones_like(e) * m).flatten()[ind] + + sum = n_sorted.cumsum() + if self.nel < sum[0]: + ifermi = 0 + elif self.nel > sum[-1]: + raise ("number of electrons larger than number of orbital*spin") + else: + ifermi = np.searchsorted(sum, self.nel) + try: + if ifermi == 0: + elo = e_sorted[0] + else: + elo = e_sorted[ifermi - 1] + if ifermi == len(e_sorted) - 1: + ehi = e_sorted[-1] + else: + ehi = e_sorted[ifermi + 1] + guess = e_sorted[ifermi] + dmu = np.max((self.width, guess - elo, ehi - guess)) + mu = brentq(self.root_function, guess - dmu, guess + dmu, xtol=xtol) + # mu = brent( + # self.root_function, + # brack=(guess - elo, guess, guess + dmu), + # tol=xtol) + except Exception as E: + # probably a bad guess + print("Error in finding Fermi level: ", E) + dmu = self.width + if self.nel < 1e-3: + mu = min(e_sorted) - dmu * 20 + elif self.nel - sum[-1] > -1e-3: + mu = max(e_sorted) + dmu * 20 + else: + # mu = brent( + # self.root_function, + # brack=(e_sorted[0] - dmu * 10, + # guess, + # e_sorted[-1] + dmu * 10), + # tol=xtol) + mu = brentq( + self.root_function, + e_sorted[0] - dmu * 20, + e_sorted[-1] + dmu * 20, + xtol=xtol, + ) + + if np.abs(self.root_function(mu)) > xtol * 1e4: + # raise RuntimeError( + # 'Fermi level could not be assigned reliably. Has the system fragmented?' + # ) + print( + "Fermi level could not be assigned reliably. Has the system fragmented?" + ) + + f = self.fermi(mu) + # rho=(self.eigenvecs*f).dot(self.eigenvecs.transpose()) + + self.mu, self.f = mu, f + return f + + def plot(self): + import pylab as pl + + for ik in range(self.nk): + pl.plot(self.e[ik, :], self.f[ik, :]) + pl.scatter(self.e[ik, :], self.f[ik, :]) + pl.title("occupations") + pl.xlabel("energy (Ha)") + pl.ylabel("occupation") + pl.show() + + +class GaussOccupations(Occupations): + def get_mu(self): + return self.mu + + def delta(self, energy): + """Return a delta-function centered at 'energy'.""" + x = -(((self.e - energy) / self.width) ** 2) + return np.exp(x) / (np.sqrt(np.pi) * self.width) + + def get_dos(self, npts=500): + eflat = self.e.flatten() + ind = np.argsort(eflat) + ##e_sorted = eflat[ind] + if self.nspin == 1: + m = 2 + elif self.nspin == 2: + m = 1 + # n_sorted = (self.wk * np.ones_like(self.e) * m).flatten()[ind] + dos = np.zeros(npts) + for w, e_n in zip(self.w_k, self.e_skn[0]): + for e in e_n: + dos += w * self.delta(e) + + def root_function(self, mu): + pass + + # @profile + def occupy(self, e, xtol=1e-8, guess=0.0): + self.e = e + dos = myDOS(kweights=self.wk, eigenvalues=e, width=self.width, npts=501) + edos = dos.get_energies() + d = dos.get_dos() + idos = integrate.cumtrapz(d, edos, initial=0) - self.nel + # f_idos = interpolate.interp1d(edos, idos) + # ret = optimize.fmin(f_idos, x0=edos[400], xtol=xtol, disp=True) + ifermi = np.searchsorted(idos, 0.0) + # self.mu = ret[0] + self.mu = edos[ifermi] + self.f = self.fermi(self.mu) + return self.f + + +class myDOS(DOS): + def __init__( + self, kweights, eigenvalues, nspin=1, width=0.1, window=None, npts=1001 + ): + """Electronic Density Of States object. + + calc: calculator object + Any ASE compliant calculator object. + width: float + Width of guassian smearing. Use width=0.0 for linear tetrahedron + interpolation. + window: tuple of two float + Use ``window=(emin, emax)``. If not specified, a window + big enough to hold all the eigenvalues will be used. + npts: int + Number of points. + + """ + self.npts = npts + self.width = width + # self.w_k = calc.get_k_point_weights() + self.w_k = kweights + self.nspins = nspin + # self.e_skn = np.array([[calc.get_eigenvalues(kpt=k, spin=s) + # for k in range(len(self.w_k))] + # for s in range(self.nspins)]) + # self.e_skn -= calc.get_fermi_level() + self.e_skn = np.array([eigenvalues.T]) # eigenvalues: iband, ikpt + + if window is None: + emin = None + emax = None + else: + emin, emax = window + + if emin is None: + emin = self.e_skn.min() - 10 * self.width + if emax is None: + emax = self.e_skn.max() + 10 * self.width + + self.energies = np.linspace(emin, emax, npts) + + # if width == 0.0: # To use tetrahedron method + # bzkpts = calc.get_bz_k_points() + # size, offset = get_monkhorst_pack_size_and_offset(bzkpts) + # bz2ibz = calc.get_bz_to_ibz_map() + # shape = (self.nspins,) + tuple(size) + (-1,) + # self.e_skn = self.e_skn[:, bz2ibz].reshape(shape) + # self.cell = calc.atoms.cell + + def get_idos(self): + e, d = self.get_dos() + return np.trapz(d, e) + + def delta(self, energy): + """Return a delta-function centered at 'energy'.""" + x = -(((self.energies - energy) / self.width) ** 2) + return np.exp(x) / (np.sqrt(np.pi) * self.width) + + def get_dos(self, spin=None): + """Get array of DOS values. + + The *spin* argument can be 0 or 1 (spin up or down) - if not + specified, the total DOS is returned. + """ + + if spin is None: + if self.nspins == 2: + # Spin-polarized calculation, but no spin specified - + # return the total DOS: + return self.get_dos(spin=0) + self.get_dos(spin=1) + else: + spin = 0 + + dos = np.zeros(self.npts) + for w, e_n in zip(self.w_k, self.e_skn[spin]): + for e in e_n: + dos += w * self.delta(e) + return dos diff --git a/TB2J/abacus/test_density_matrix.py b/TB2J/abacus/test_density_matrix.py new file mode 100644 index 0000000..589d901 --- /dev/null +++ b/TB2J/abacus/test_density_matrix.py @@ -0,0 +1,38 @@ +from scipy.linalg import eigh +import numpy as np + + +def gen_random_hermitean_matrix(n): + A = np.random.rand(n, n) + 1j * np.random.rand(n, n) + return A + A.conj().T + + +def gen_overlap_matrix(n): + A = np.random.rand(n, n) + 1j * np.random.rand(n, n) + return np.dot(A, A.conj().T) + + +def fermi_function(x, ef, beta): + return 1.0 / (np.exp(beta * (x - ef)) + 1) + + +def test(): + n = 10 + A = gen_random_hermitean_matrix(n) + S = gen_overlap_matrix(n) + beta = 0.1 + ef = 0 + + evals, evecs = eigh(A, S) + + etot = np.sum(evals * fermi_function(evals, ef, beta)) + + rho = np.einsum("ib,b,jb->ij", evecs, fermi_function(evals, ef, beta), evecs.conj()) + + etot2 = np.trace(np.dot(A, rho)) + + print(etot, etot2) + + +if __name__ == "__main__": + test() diff --git a/TB2J/exchange.py b/TB2J/exchange.py index 55ed7c4..fb60962 100644 --- a/TB2J/exchange.py +++ b/TB2J/exchange.py @@ -345,7 +345,7 @@ def set_tbmodels(self, tbmodels): self.norb = self.G.norb self.nbasis = self.G.nbasis # self.rho = np.zeros((self.nbasis, self.nbasis), dtype=complex) - self.rho = self.G.get_density_matrix().real + self.rho = self.G.get_density_matrix() self.A_ijR_list = defaultdict(lambda: []) self.A_ijR = defaultdict(lambda: np.zeros((4, 4), dtype=complex)) self.A_ijR_orb = dict() diff --git a/TB2J/green.py b/TB2J/green.py index d5d0884..7862311 100644 --- a/TB2J/green.py +++ b/TB2J/green.py @@ -7,6 +7,8 @@ from pathos.multiprocessing import ProcessPool import sys import pickle +import warnings +from TB2J.mathutils.fermi import fermi MAX_EXP_ARGUMENT = np.log(sys.float_info.max) @@ -26,19 +28,6 @@ def eigen_to_G(evals, evecs, efermi, energy): ) -def fermi(e, mu, width=0.01): - """ - the fermi function. - .. math:: - f=\\frac{1}{\exp((e-\mu)/width)+1} - - :param e,mu,width: e,\mu,width - """ - - x = (e - mu) / width - return np.where(x < MAX_EXP_ARGUMENT, 1 / (1.0 + np.exp(x)), 0.0) - - def find_energy_ingap(evals, rbound, gap=4.0): """ find a energy inside a gap below rbound (right bound), diff --git a/TB2J/mathutils/__init__.py b/TB2J/mathutils/__init__.py new file mode 100644 index 0000000..4d8f75e --- /dev/null +++ b/TB2J/mathutils/__init__.py @@ -0,0 +1 @@ +from .lowdin import Lowdin diff --git a/TB2J/mathutils/fermi.py b/TB2J/mathutils/fermi.py new file mode 100644 index 0000000..66d60ec --- /dev/null +++ b/TB2J/mathutils/fermi.py @@ -0,0 +1,22 @@ +import numpy as np +import warnings +import sys + +MAX_EXP_ARGUMENT = np.log(sys.float_info.max) + + +def fermi(e, mu, width=0.01): + """ + the fermi function. + .. math:: + f=\\frac{1}{\exp((e-\mu)/width)+1} + + :param e,mu,width: e,\mu,width + """ + x = (e - mu) / width + # disable overflow warning + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + ret = np.where(x < MAX_EXP_ARGUMENT, 1 / (1.0 + np.exp(x)), 0.0) + + return ret diff --git a/TB2J/mathutils/kR_convert.py b/TB2J/mathutils/kR_convert.py new file mode 100644 index 0000000..eaa2501 --- /dev/null +++ b/TB2J/mathutils/kR_convert.py @@ -0,0 +1,90 @@ +import numpy as np + + +def HR_to_k(HR, Rlist, kpts): + # Hk[k,:,:] = sum_R (H[R] exp(i2pi k.R)) + phase = np.exp(2.0j * np.pi * np.tensordot(kpts, Rlist, axes=([1], [1]))) + Hk = np.einsum("rlm, kr -> klm", HR, phase) + return Hk + + +def Hk_to_R(Hk, Rlist, kpts, kweights): + phase = np.exp(-2.0j * np.pi * np.tensordot(kpts, Rlist, axes=([1], [1]))) + HR = np.einsum("klm, kr, k->rlm", Hk, phase, kweights) + return HR + + +def k_to_R(kpts, Rlist, Mk, kweights=None): + """ + Transform k-space wavefunctions to real space. + params: + kpts: k-points + Rlist: list of R vectors + Mk: matrix of shape [nkpt, n1, n2] in k-space. + + return: + MR: matrix of shape [nR, n1, n2], the matrix in R-space. + + """ + nkpt, n1, n2 = Mk.shape + if kweights is None: + kweights = np.ones(nkpt, dtype=float) / nkpt + phase = np.exp(-2.0j * np.pi * np.tensordot(kpts, Rlist, axes=([1], [1]))) + MR = np.einsum("klm, kr, k -> rlm", Mk, phase, kweights) + return MR + + # nkpt, n1, n2 = Mk.shape + # nR = Rlist.shape[0] + # MR = np.zeros((nR, n1, n2), dtype=complex) + # if kweights is None: + # kweights = np.ones(nkpt, dtype=float)/nkpt + # for iR, R in enumerate(Rlist): + # for ik in range(nkpt): + # MR[iR] += Mk[ik] * np.exp(-2.0j*np.pi * np.dot(kpts[ik], R)) * kweights[ik] + # return MR + + +def R_to_k(kpts, Rlist, MR): + """ + Transform real-space wavefunctions to k-space. + params: + kpts: k-points + Rlist: list of R vectors + MR: matrix of shape [nR, n1, n2] in R-space. + + return: + Mk: matrix of shape [nkpt, n1, n2], the matrix in k-space. + + """ + phase = np.exp(2.0 * np.pi * 1j * np.tensordot(kpts, Rlist, axes=([1], [1]))) + Mk = np.einsum("rlm, kr -> klm", MR, phase) + + # nkpt, n1, n2 = Mk.shape + # nR = Rlist.shape[0] + # Mk = np.zeros((nkpt, n1, n2), dtype=complex) + # for iR, R in enumerate(Rlist): + # for ik in range(nkpt): + # Mk[ik] += MR[iR] * np.exp(2.0 * np.pi * 1j * np.dot(kpts[ik], R)) + return Mk + + +def R_to_onek(kpt, Rlist, MR): + """ + Transform real-space wavefunctions to k-space. + params: + kpt: k-point + Rlist: list of R vectors + MR: matrix of shape [nR, n1, n2] in R-space. + + return: + Mk: matrix of shape [n1, n2], the matrix in k-space. + + """ + phase = np.exp(2.0j * np.pi * np.dot(Rlist, kpt)) + Mk = np.einsum("rlm, r -> lm", MR, phase) + return Mk + # n1, n2 = MR.shape[1:] + # Mk = np.zeros((n1, n2), dtype=complex) + # for iR, R in enumerate(Rlist): + # Mk += MR[iR] * np.exp(2.0j*np.pi * np.dot(kpt, R)) + # return Mk diff --git a/TB2J/mathutils/lowdin.py b/TB2J/mathutils/lowdin.py new file mode 100644 index 0000000..bee9ccc --- /dev/null +++ b/TB2J/mathutils/lowdin.py @@ -0,0 +1,12 @@ +import numpy as np +from scipy.linalg import inv, eigh + + +def Lowdin(S): + """ + Calculate S^(-1/2). + Which is used in lowind's symmetric orthonormalization. + psi_prime = S^(-1/2) psi + """ + eigval, eigvec = eigh(S) + return eigvec @ np.diag(np.sqrt(1.0 / eigval)) @ (eigvec.T.conj()) diff --git a/TB2J/mathutils/rotate_spin.py b/TB2J/mathutils/rotate_spin.py new file mode 100644 index 0000000..af0c5b2 --- /dev/null +++ b/TB2J/mathutils/rotate_spin.py @@ -0,0 +1,35 @@ +import numpy as np +from TB2J.pauli import pauli_block_all, s0, s1, s2, s3, gather_pauli_blocks + + +def rotate_Matrix_from_z_to_axis(M, axis, normalize=True): + """ + Given a spinor matrix M, rotate it from z-axis to axis. + The spinor matrix M is a 2x2 matrix, which can be decomposed as I, x, y, z components using Pauli matrices. + """ + MI, Mx, My, Mz = pauli_block_all(M) + axis = axis / np.linalg.norm(axis) + # M_new = s0* MI + Mz * (axis[0] * s1 + axis[1] * s2 + axis[2] * s3) *2 + M_new = gather_pauli_blocks(MI, Mz * axis[0], Mz * axis[1], Mz * axis[2]) + return M_new + + +def test_rotate_Matrix_from_z_to_axis(): + M = np.array([[1.1, 0], [0, 0.9]]) + print(pauli_block_all(M)) + Mnew = rotate_Matrix_from_z_to_axis(M, [1, 1, 1]) + print(pauli_block_all(Mnew)) + print(Mnew) + + M = np.array( + [ + [-9.90532976e-06 + 0.0j, 0.00000000e00 + 0.0j], + [0.00000000e00 + 0.0j, -9.88431291e-06 + 0.0j], + ] + ) + print(M) + print(rotate_Matrix_from_z_to_axis(M, [0, 0, 1])) + + +if __name__ == "__main__": + test_rotate_Matrix_from_z_to_axis() diff --git a/TB2J/pauli.py b/TB2J/pauli.py index 0733919..ec34029 100644 --- a/TB2J/pauli.py +++ b/TB2J/pauli.py @@ -143,7 +143,24 @@ def pauli_block_all(M): return MI, Mx, My, Mz +def gather_pauli_blocks(MI, Mx, My, Mz): + """ + Gather the I, x, y, z component of a matrix. + """ + return np.kron(MI, s0) + np.kron(Mx, s1) + np.kron(My, s2) + np.kron(Mz, s3) + + +def test_gather_pauli_blocks(): + M = np.random.rand(4, 4) + MI, Mx, My, Mz = pauli_block_all(M) + M2 = gather_pauli_blocks(MI, Mx, My, Mz) + assert np.allclose(M, M2) + + def op_norm(M): + """ + Return the operator norm of a matrix. + """ return max(svd(M)[1]) diff --git a/TB2J/utils.py b/TB2J/utils.py index d9c74ae..1a7347b 100644 --- a/TB2J/utils.py +++ b/TB2J/utils.py @@ -87,6 +87,7 @@ def auto_assign_wannier_to_atom(positions, atoms, max_distance=0.1, half=False): """ pos = np.array(positions) atompos = atoms.get_scaled_positions(wrap=False) + cell = atoms.get_cell() ind_atoms = [] newpos = [] refpos = [] @@ -95,8 +96,9 @@ def auto_assign_wannier_to_atom(positions, atoms, max_distance=0.1, half=False): dp = p[None, :] - atompos # residual of d r = dp - np.round(dp) + r_cart = r @ cell # find the min of residual - normd = np.linalg.norm(r, axis=1) + normd = np.linalg.norm(r_cart, axis=1) iatom = np.argmin(normd) # ref+residual rmin = r[iatom] @@ -330,3 +332,82 @@ def simpson_nonuniform(x, f): result += f[N - 1] * (h[N - 1] ** 2 + 3 * h[N - 1] * h[N - 2]) / (6 * h[N - 2]) result -= f[N - 2] * h[N - 1] ** 3 / (6 * h[N - 2] * (h[N - 2] + h[N - 1])) return result + + +def simpson_nonuniform_weight(x): + """ + Simpson rule for irregularly spaced data. + x: list or np.array of floats + Sampling points for the function values + Returns + ------- + weight : list or np.array of floats + weight for the Simpson rule + For the function f(x), the integral is approximated as + $\int f(x) dx \approx \sum_i weight[i] * f(x[i])$ + """ + + weight = np.zeros_like(x) + N = len(x) - 1 + h = np.diff(x) + + for i in range(1, N, 2): + hph = h[i] + h[i - 1] + weight[i] += (h[i] ** 3 + h[i - 1] ** 3 + 3.0 * h[i] * h[i - 1] * hph) / ( + 6 * h[i] * h[i - 1] + ) + weight[i - 1] += ( + 2.0 * h[i - 1] ** 3 - h[i] ** 3 + 3.0 * h[i] * h[i - 1] ** 2 + ) / (6 * h[i - 1] * hph) + weight[i + 1] += ( + 2.0 * h[i] ** 3 - h[i - 1] ** 3 + 3.0 * h[i - 1] * h[i] ** 2 + ) / (6 * h[i] * hph) + + if (N + 1) % 2 == 0: + weight[N] += (2 * h[N - 1] ** 2 + 3.0 * h[N - 2] * h[N - 1]) / ( + 6 * (h[N - 2] + h[N - 1]) + ) + weight[N - 1] += (h[N - 1] ** 2 + 3 * h[N - 1] * h[N - 2]) / (6 * h[N - 2]) + weight[N - 2] -= h[N - 1] ** 3 / (6 * h[N - 2] * (h[N - 2] + h[N - 1])) + return weight + + +def trapz_nonuniform_weight(x): + """ + trapezoidal rule for irregularly spaced data. + x: list or np.array of floats + Sampling points for the function values + Returns + ------- + weight : list or np.array of floats + weight for the trapezoidal rule + For the function f(x), the integral is approximated as + $\int f(x) dx \approx \sum_i weight[i] * f(x[i])$ + """ + h = np.diff(x) + weight = np.zeros_like(x) + weight[0] = h[0] / 2.0 + weight[1:-1] = (h[1:] + h[:-1]) / 2.0 + weight[-1] = h[-1] / 2.0 + return weight + + +def test_simpson_nonuniform(): + x = np.array([0.0, 0.1, 0.3, 0.5, 0.8, 1.0]) + w = simpson_nonuniform_weight(x) + # assert np.allclose(w, [0.1, 0.4, 0.4, 0.4, 0.4, 0.1]) + assert np.allclose(simpson_nonuniform(x, x**8), 0.12714277533333335) + print("simpson_weight:", simpson_nonuniform_weight(x) @ x**8, 0.12714277533333335) + print("trapz_weight:", trapz_nonuniform_weight(x) @ x**8) + + x2 = np.linspace(0, 1, 500) + print(simpson_nonuniform_weight(x2) @ x2**8, 1 / 9.0) + print(simpson_nonuniform_weight(x2) @ x2**8) + print("simpson_weight:", simpson_nonuniform_weight(x2) @ x2**8) + print("trapz_weight:", trapz_nonuniform_weight(x2) @ x2**8) + + assert np.allclose(simpson_nonuniform(x, x**8), 1 / 9.0) + + +if __name__ == "__main__": + test_simpson_nonuniform() diff --git a/docs/conf.py b/docs/conf.py index 6c291bd..702c819 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -34,6 +34,7 @@ "sphinx_rtd_theme", "sphinx.ext.mathjax", "sphinx.ext.autodoc", + "sphinxemoji.sphinxemoji", ] # Add any paths that contain templates here, relative to this directory. diff --git a/docs/requirements.txt b/docs/requirements.txt index 906f386..376e690 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,3 +1,3 @@ recommonmark sphinx_rtd_theme - +sphinxemoji diff --git a/docs/src/Contributors.rst b/docs/src/Contributors.rst index ea653be..26dfadd 100644 --- a/docs/src/Contributors.rst +++ b/docs/src/Contributors.rst @@ -1,7 +1,7 @@ Contributors ============ -The main contributors to the TB2J package are: +The TB2J package is initially developped: * Xu He (mailhexu@gmail.com) @@ -11,6 +11,29 @@ The main contributors to the TB2J package are: * Eric Bousquet +from the Phythema and Nanomat groups at the University of Liège. +Part of the code is developped by Xu He during he was at ICN2 (Institut Catala de Nanociencia i Nanotecnologia) in Barcelona, Spain. + + +After the first public release, the following people has made significant contribution to the code: + +* Andres Tellez Mora (West Virginia University) + +* Aldo Romero (West Virginia University) + +* Andrey Rubakov (ICMol, University of Valencia) + +* Gan Jin (University of Science and Technology of China) + +* Zhenxiong Shen (University of Science and Technology of China) + + +We thank all the other people who have contributed to the code, including the ones who have reported bugs or typos, + suggested improvements, and asked questions. + + + + We welcome contributions. You can help us with: - extending the input format to other codes, e.g. first principles or tight binding code. - extending the output to other spin dynamics code. diff --git a/docs/src/ReleaseNotes.md b/docs/src/ReleaseNotes.md index 9c423c8..ff6c630 100644 --- a/docs/src/ReleaseNotes.md +++ b/docs/src/ReleaseNotes.md @@ -1,5 +1,9 @@ ## Release Notes ------------------------------------------------------------------------ +#### v0.9.0 March 22, 2024 +Improved merge method for anisotropic exchange and DMI. (thanks to Andres Tellez Mora!) + + #### v0.8.2 March 4, 2024 TB2J can now read the "tb.dat" file instead of the "hr.dat"+"centers.xyz" files. diff --git a/docs/src/abacus.md b/docs/src/abacus.md index 123a71c..08ce0a3 100644 --- a/docs/src/abacus.md +++ b/docs/src/abacus.md @@ -99,3 +99,5 @@ options: ``` + + diff --git a/docs/src/mae.md b/docs/src/mae.md new file mode 100644 index 0000000..9eeb102 --- /dev/null +++ b/docs/src/mae.md @@ -0,0 +1,127 @@ +### Computing Magnetocrystalline anisotropy energy (MAE) . + +:warning: This feature is currently under development and internal test. Do not use it for production yet. +:warning: This feature is only available with the ABACUS code. + + +To compute the magnetocrystalline anisotropy energy (MAE) of a magnetic system with the magnetic force theorem, two steps of DFT calculations are needed. + +- The first step is to do an collinear spin calculation. The density and the Hamiltonian is saved at this step. Note that the current implementation requires the SOC to be turned on in ABACUS, but setting the SOC strength to zero (soc_lambda=0). + +- The second step is to do a non-SCF non-collinear spin calculation with SOC turned on. The density is read from the previous step. In practice, one step of SCF calculation is done (as the current implementation does not write the Hamiltonian and the energy). The Hamiltonian should be saved in this step, too. + +Here is one example: +Step one: collinear spin calculation. Note that instead of using nspin=2, we use nspin=4, and lspinorb=1 to enable the SOC but set the soc\_lambda to 0.0 to turn off the SOC. This is to make the Hamiltonian saved in the spinor form, so it can be easily compared with the next step of a real calculation with SOC. + + +``` text +INPUT_PARAMETERS +# SCF calculation with SOC turned on, but soc_lambda=0. +calculation scf +nspin 4 +symmetry 0 +noncolin 1 +lspinorb 1 +ecutwfc 100 +scf_thr 1e-06 +init_chg atomic +out_mul 1 +out_chg 1 +out_dos 0 +out_band 0 +out_wfc_lcao 1 +out_mat_hs2 1 +ks_solver scalapack_gvx +scf_nmax 500 +out_bandgap 0 +basis_type lcao +gamma_only 0 +smearing_method gaussian +smearing_sigma 0.01 +mixing_type broyden +mixing_beta 0.5 +soc_lambda 0.0 +ntype 1 +dft_functional PBE +``` + + +- Step two: non-SCF non-collinear spin calculation. + +In this step, we need to start from the density saved in the previous step. So we can copy the output directory of the previous step to a the present directory. +To mimic the non-SCF calculation: + * The scf\_nmax should be set to 1 to end the scf calculation in one step. + * The "scf\_thr" should be set to a large value to make the calculation "converge" in one step. + * The mixing\_beta should be set to a small value to suppress the density mixing. + * Then we can set the "init\_chg" to "file" to read the density from the previous step. The lspionorb should be set to 1, and the soc_lambda should be set to a 1.0 to enable the SOC. +The "out\_mat\_hs2" should be set to 1 to save the Hamiltonian. + +:warning: Once ABACUS can output the Hamiltonian in the non-SCF calculation, we can use a "calculation=nscf" to do the non-SCF calculation. + + +```text +INPUT_PARAMETERS +# Non-SCF non-collinear spin calculation with SOC turned on. soc_lambda=1.0 +calculation scf +nspin 4 +symmetry 0 +noncolin 1 +lspinorb 1 +ecutwfc 100 +scf_thr 1000000.0 +init_chg file +out_mul 1 +out_chg 0 +out_mat_hs2 1 +ks_solver scalapack_gvx +scf_nmax 1 +basis_type lcao +smearing_method gaussian +smearing_sigma 0.01 +mixing_type broyden +mixing_beta 1e-06 +soc_lambda 1.0 +init_wfc file +ntype 1 +dft_functional PBE +``` + + +After the two steps of calculations, we can use the "abacus\_get\_MAE" function to compute the MAE. Essentially, the function reads the Hamiltonian from the two steps of calculations, and substract them to the the SOC part and the non-soc part. +Then the non-SOC part is rotated along the different directions, and the energy difference is computed. We can define the magnetic moment directions by list of thetas/psis in the spherical coordinates, and explore the energies. + +Here is an example of the usage of the function. + +```python + +import numpy as np +from TB2J.MAE import abacus_get_MAE + +def run(): + # theta, psi: along the xz plane, rotating from z to x. + thetas = np.linspace(0, 180, 19) * np.pi / 180 + psis = np.zeros(19) + abacus_get_MAE( + path_nosoc= 'soc0/OUT.ABACUS', + path_soc= 'soc1/OUT.ABACUS', + kmesh=[6,6,1], + gamma=True, + thetas=thetas, + psis=psis, + nel = 16, + ) + ) + +if __name__ == '__main__': + run() +``` + + +Here the soc0 and soc1 directories are where the two ABACUS calculations are done. We use a 6x6x1 Gamma-centered k-mesh for the integration (which is too small for practical usage!). And explore the energy with the magnetic moments in the x-z plane. +After running the python script aboves, a file named "MAE.txt" will be created, including the theta, psi, and the energies (in eV). + + + + + + diff --git a/docs/src/rotate_and_merge.rst b/docs/src/rotate_and_merge.rst index 0660033..0f2b11c 100644 --- a/docs/src/rotate_and_merge.rst +++ b/docs/src/rotate_and_merge.rst @@ -3,15 +3,19 @@ Averaging multiple parameters =============================== -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`. +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}}`. To obtain the other components, we need to rotate the spins or alternatively, rotate the structure while keeping the spin directions fixed. This method takes the approximation that the electronic strucuture is only slightly affected by the rotation of the spins, which will only lead to ne +gligible relative differences in the magnetic interaction parameters. This is a good approximation for systems with weak SOC (i.e. the SOC is much weaker than the exchange-correlation). But it can fail for systems with strong SOC. + +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 :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. +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. Note that when the system is only slightly noncollinear, e.g. AFM systems with weak-ferromagnetism due to spin canting, we can still treat it as collinear and three calculations is still enough. -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: +While rotating the spins can be done in the DFT calculation, the feature is sometimes not available for some DFT codes. In this case, an alternative is to rotate the structures while keeping the spins fixed. +To account for this, TB2J provides scripts to rotate the structure named TB2J\_rotate.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 whether 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 --ftype parameter. An example for using the rotate file with a collinear system is: .. code-block:: shell @@ -23,7 +27,9 @@ If te system is noncollinear, then we run the following instead: 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: +The user has to perform DFT single point energy calculations for the same structure with different spin orientations, or 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 results with the following command by providing the paths to the TB2J results of the three cases: :: @@ -31,8 +37,11 @@ The user has to perform DFT single point energy calculations for the generated s 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 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: :: @@ -40,3 +49,6 @@ Another method is to do the DFT calculation with spins rotated globally. That is 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. + + + diff --git a/docs/src/tutorial.rst b/docs/src/tutorial.rst index 087464d..f6cf847 100644 --- a/docs/src/tutorial.rst +++ b/docs/src/tutorial.rst @@ -8,6 +8,7 @@ Tutorial siesta.rst openmx.rst abacus.rst + mae.md parameters.rst rotate_and_merge.rst downfold.md diff --git a/scripts/TB2J_downfold.py b/scripts/TB2J_downfold.py index 2cf6114..160a5fb 100755 --- a/scripts/TB2J_downfold.py +++ b/scripts/TB2J_downfold.py @@ -58,6 +58,13 @@ def main(): default=False, ) + parser.add_argument( + "--method", + help="The method to downfold the exchange parameter. Options are Lowdin and PWF (projected Wannier function). ", + type=str, + default="Lowdin", + ) + args = parser.parse_args() if len(args.metals) == []: @@ -73,6 +80,7 @@ def main(): outpath=args.outpath, qmesh=args.qmesh, iso_only=args.iso_only, + method=args.method, ) diff --git a/setup.py b/setup.py index 2ab2127..16e417e 100644 --- a/setup.py +++ b/setup.py @@ -1,7 +1,7 @@ #!/usr/bin/env python from setuptools import setup, find_packages -__version__ = "0.8.2.8" +__version__ = "0.9.0.2" long_description = """TB2J is a Python package aimed to compute automatically the magnetic interactions (superexchange and Dzyaloshinskii-Moriya) between atoms of magnetic crystals from DFT Hamiltonian based on Wannier functions or Linear combination of atomic orbitals. It uses the Green's function method and take the local rigid spin rotation as a perturbation. The package can take the output from Wannier90, which is interfaced with many density functional theory codes or from codes based on localised orbitals. A minimal user input is needed, which allows for an easily integration into a high-throughput workflows. """ diff --git a/upload_to_pip.sh b/upload_to_pip.sh index e8422e3..6b1f1d5 100755 --- a/upload_to_pip.sh +++ b/upload_to_pip.sh @@ -1,4 +1,4 @@ #!/usr/bin/env bash rm ./dist/* -python3 setup.py sdist bdist_wheel -python3 -m twine upload --repository pypi dist/* --verbose +python3.11 setup.py sdist bdist_wheel +python3.11 -m twine upload --repository pypi dist/* --verbose