diff --git a/setup.py b/setup.py index b55efe9..62b34df 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ setuptools.setup( name="tb-mean-field-hubbard", - version="1.3.1", + version="1.4.0", author="Kristjan Eimre", author_email="kristjaneimre@gmail.com", description="Package to run tight-binding mean field hubbard calculations", diff --git a/tb_mean_field_hubbard/__init__.py b/tb_mean_field_hubbard/__init__.py index 95e1917..b169645 100644 --- a/tb_mean_field_hubbard/__init__.py +++ b/tb_mean_field_hubbard/__init__.py @@ -1,3 +1,3 @@ from .mfh import MeanFieldHubbardModel -__version__ = "1.3.1" +__version__ = "1.4.0" diff --git a/tb_mean_field_hubbard/mfh.py b/tb_mean_field_hubbard/mfh.py index 077fce0..9828188 100644 --- a/tb_mean_field_hubbard/mfh.py +++ b/tb_mean_field_hubbard/mfh.py @@ -20,12 +20,18 @@ class MeanFieldHubbardModel: - def __init__(self, ase_geom, t_list=[2.7], charge=0, multiplicity=1, bond_cutoff='auto'): + def __init__(self, ase_geom, t_list=[2.7], charge=0, multiplicity='auto', bond_cutoff='auto'): self.t_list = t_list - self.multiplicity = multiplicity self.charge = charge + if multiplicity == 'auto': + self.multiplicity = 0 + self.relax_multiplicity = True + else: + self.multiplicity = multiplicity + self.relax_multiplicity = False + self.ase_geom = ase_geom self.num_atoms = len(ase_geom) @@ -37,7 +43,7 @@ def __init__(self, ase_geom, t_list=[2.7], charge=0, multiplicity=1, bond_cutoff if bond_cutoff == 'auto': cc_len = utils.get_cc_bond_len(self.ase_geom) # graphene 2nd nn distance is 1.73*cc, so use halfway there as the cutoff - bond_cutoff = 1.37*cc_len + bond_cutoff = 1.37 * cc_len self.neighbor_list = ase.neighborlist.neighbor_list('ij', self.ase_geom, bond_cutoff) @@ -76,19 +82,21 @@ def _load_spin_guess(self, ase_geom, flip_alpha_majority=True): def _set_up_tb_model(self): - # check parameter validity: - self.num_el = self.num_atoms - self.charge - if self.multiplicity % 2 == self.num_el % 2: - raise Exception("ERROR: Charge & multiplicity combination not allowed!") + if not self.relax_multiplicity: + + if self.multiplicity % 2 == self.num_el % 2: + raise Exception("ERROR: Charge & multiplicity combination not allowed!") - # determine spin populations + # determine spin populations + self.num_spin_el = [ + (self.num_el + (self.multiplicity - 1)) // 2, + (self.num_el - (self.multiplicity - 1)) // 2, + ] - self.num_spin_el = [ - (self.num_el + (self.multiplicity - 1)) // 2, - (self.num_el - (self.multiplicity - 1)) // 2, - ] + else: + self.num_spin_el = [0, 0] lat = [[1.0, 0.0], [0.0, 1.0]] @@ -128,8 +136,9 @@ def visualize_spin_guess(self): def print_parameters(self): print("Total number of electrons: %d" % self.num_el) - print("α electrons: %d" % self.num_spin_el[0]) - print("β electrons: %d" % self.num_spin_el[1]) + if not self.relax_multiplicity: + print("α electrons: %d" % self.num_spin_el[0]) + print("β electrons: %d" % self.num_spin_el[1]) ### ------------------------------------------------------------------- ### MFH routines @@ -141,11 +150,26 @@ def new_spin_res_dens(self, evals_list, evecs_list): new_spin_resolved_dens = (1.0 - mix_coef) * np.copy(self.spin_resolved_dens) - for i_spin, n_el in enumerate(self.num_spin_el): - - for i_mo in range(n_el): - - new_spin_resolved_dens[:, i_spin] += mix_coef * np.abs(evecs_list[i_spin][i_mo])**2 + if not self.relax_multiplicity: + + for i_spin, n_el in enumerate(self.num_spin_el): + for i_mo in range(n_el): + new_spin_resolved_dens[:, i_spin] += mix_coef * np.abs(evecs_list[i_spin][i_mo])**2 + + else: + # in case of relaxed multiplicity, fill the lowest orbitals, without considering spin + i_a = 0 + i_b = 0 + for _i_el in range(self.num_el): + en_a = evals_list[0][i_a] + en_b = evals_list[1][i_b] + if en_a < en_b: + new_spin_resolved_dens[:, 0] += mix_coef * np.abs(evecs_list[0][i_a])**2 + i_a += 1 + else: + new_spin_resolved_dens[:, 1] += mix_coef * np.abs(evecs_list[1][i_b])**2 + i_b += 1 + self.num_spin_el = [i_a, i_b] return new_spin_resolved_dens @@ -211,6 +235,9 @@ def run_mfh(self, u, print_iter=False, plot=False, energy_tol=1e-6, mag_tol=1e-4 self.evals = evals self.evecs = evecs + if self.relax_multiplicity: + self.multiplicity = max(self.num_spin_el) - min(self.num_spin_el) + 1 + # post-process # constant energy term in the MFH [ U*sum() ] @@ -349,8 +376,9 @@ def plot_sts_map(self, def report(self, num_orb=2, sts_h=3.5, sts_broad=0.05): - print("abs. magnetization: %14.6f" % self.abs_mag) - print("energy: %14.6f" % self.energy) + print(f"multiplicity: {self.multiplicity:12d}") + print(f"abs. magnetization: {self.abs_mag:12.4f}") + print(f"energy: {self.energy: 12.4f}") print("---") print("spin density:") plt.figure(figsize=self.figure_size) @@ -362,9 +390,9 @@ def report(self, num_orb=2, sts_h=3.5, sts_broad=0.05): self.plot_evals(self.evals) gap_a, gap_b, gap = self.gaps(self.evals) - print("gap alpha: %.6f" % gap_a) - print("gap beta: %.6f" % gap_b) - print("gap eff.: %.6f" % gap) + print(f"gap alpha: {gap_a:.4f}") + print(f"gap beta: {gap_b:.4f}") + print(f"gap eff.: {gap:.4f}") print("---") print("frontier orbitals:")