Skip to content

Add support for La3Ni2O7 #46

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Julpe wants to merge 1 commit into
base: main
Choose a base branch
from
Draft

Add support for La3Ni2O7 #46

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions moldga/config.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -148,6 +148,8 @@ def __init__(self):
self.mu: float = 0.0
self.n: float = 0.0
self.n_bands: int = 1
self.nd_bands: int = 1
self.np_bands: int = 0
self.occ: np.ndarray = np.ndarray(0)
self.occ_k: np.ndarray = np.ndarray(0)
self.occ_dmft: np.ndarray = np.ndarray(0)
Expand Down
94 changes: 87 additions & 7 deletions moldga/dga_io.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,13 @@ def uniquify_path(path: str = None):
return path


def load_from_w2dyn_file_and_update_config() -> tuple[GreensFunction, SelfEnergy, LocalFourPoint, LocalFourPoint]:
def load_from_w2dyn_file_and_update_config(
combine_two_atoms_to_one_obj: bool = False,
) -> tuple[GreensFunction, SelfEnergy, LocalFourPoint, LocalFourPoint]:
"""
Loads data from the w2dyn file and updates the config file.
If combine_atoms_to_one_obj is True, we are doubling the orbital space and putting the data from the
inequivalent atoms into the diagonal blocks after we average over them.
"""
file = w2dyn_aux.W2dynFile(fname=str(os.path.join(config.dmft.input_path, config.dmft.fname_1p)))

Expand All @@ -51,36 +55,106 @@ def load_from_w2dyn_file_and_update_config() -> tuple[GreensFunction, SelfEnergy
config.lattice.interaction.vpp = file.get_vpp()

config.sys.mu = file.get_mu()
config.sys.n_bands = file.get_nd() + file.get_np()
config.sys.nd_bands = file.get_nd()
config.sys.np_bands = file.get_np()
config.sys.n_bands = config.sys.nd_bands + config.sys.np_bands
config.sys.n = file.get_totdens()
config.sys.occ_dmft = 2 * np.mean(file.get_rho1(), axis=(1, 3))

if combine_two_atoms_to_one_obj:
config.sys.n_bands *= 2
config.sys.nd_bands *= 2
config.sys.np_bands *= 2
config.sys.n *= 2

config.sys.occ_dmft = np.zeros((config.sys.n_bands, config.sys.n_bands))
rho_1_mean = 0.5 * (np.mean(file.get_rho1(), axis=(1, 3)) + np.mean(file.get_rho1(ineq=2), axis=(1, 3)))
config.sys.occ_dmft[0:2, 0:2] = 2 * rho_1_mean
config.sys.occ_dmft[2:4, 2:4] = 2 * rho_1_mean

if config.sys.n == 0:
config.sys.n = 2 * np.sum(config.sys.occ_dmft)

file2 = w2dyn_aux.W2dynG4iwFile(fname=str(os.path.join(config.dmft.input_path, config.dmft.fname_2p)))
g2_dens = LocalFourPoint(file2.read_g2_full_multiband(config.sys.n_bands, name="dens"), channel=SpinChannel.DENS)
g2_magn = LocalFourPoint(file2.read_g2_full_multiband(config.sys.n_bands, name="magn"), channel=SpinChannel.MAGN)
g2_dens = LocalFourPoint(
file2.read_g2_full_multiband(file.get_nd() + file.get_np(), name="dens"), channel=SpinChannel.DENS
)
g2_magn = LocalFourPoint(
file2.read_g2_full_multiband(file.get_nd() + file.get_np(), name="magn"), channel=SpinChannel.MAGN
)

if combine_two_atoms_to_one_obj:
g2_dens_2 = LocalFourPoint(
file2.read_g2_full_multiband(file.get_nd() + file.get_np(), ineq=2, name="dens"), channel=SpinChannel.DENS
)
g2_magn_2 = LocalFourPoint(
file2.read_g2_full_multiband(file.get_nd() + file.get_np(), ineq=2, name="magn"), channel=SpinChannel.MAGN
)

def construct_large_g2(g2_1: LocalFourPoint, g2_2: LocalFourPoint) -> LocalFourPoint:
g2_mean = 0.5 * (g2_1.mat + g2_2.mat)
del g2_2
g2_1.mat = np.zeros((4, 4, 4, 4, *g2_mean.shape[4:]))
g2_1.mat[0:2, 0:2, 0:2, 0:2] = g2_mean
g2_1.mat[2:4, 2:4, 2:4, 2:4] = g2_mean
del g2_mean
return g2_1

g2_dens = construct_large_g2(g2_dens, g2_dens_2)
g2_magn = construct_large_g2(g2_magn, g2_magn_2)

file2.close()

update_frequency_boxes(g2_dens.niw, g2_dens.niv)

def extend_orbital(arr: np.ndarray) -> np.ndarray:
return np.einsum("i...,ij->ij...", arr, np.eye(config.sys.n_bands))
return np.einsum("i...,ij->ij...", arr, np.eye(arr.shape[0]))

giw_spin_mean = np.mean(file.get_giw(), axis=1) # [band,spin,niv]
niv_dmft = giw_spin_mean.shape[-1] // 2
niv_cut = config.box.niw_core + config.box.niv_full + 10
giw_spin_mean = giw_spin_mean[..., niv_dmft - niv_cut : niv_dmft + niv_cut]
g_dmft = GreensFunction(extend_orbital(giw_spin_mean))

if combine_two_atoms_to_one_obj:
giw_spin_mean_2 = np.mean(file.get_giw(ineq=2), axis=1)[..., niv_dmft - niv_cut : niv_dmft + niv_cut]
giw_spin_mean = 0.5 * (giw_spin_mean + giw_spin_mean_2)
del giw_spin_mean_2
giw_spin_mean_large = np.zeros((2 * giw_spin_mean.shape[0], *giw_spin_mean.shape[1:]))
giw_spin_mean_large[0:2] = giw_spin_mean
giw_spin_mean_large[2:4] = giw_spin_mean
g_dmft = GreensFunction(extend_orbital(giw_spin_mean_large))
del giw_spin_mean_large

siw_spin_mean = np.mean(file.get_siw(), axis=1) # [band,spin,niv]
siw_spin_mean = extend_orbital(siw_spin_mean)[None, None, None, ...]
siw_dc_spin_mean = np.mean(file.get_dc(), axis=-1) # [band,spin]
siw_dc_spin_mean = extend_orbital(siw_dc_spin_mean)[None, None, None, ..., None]
siw_spin_mean = siw_spin_mean[..., niv_dmft - niv_cut : niv_dmft + niv_cut]
sigma_dmft = SelfEnergy(siw_spin_mean, estimate_niv_core=True) + siw_dc_spin_mean

if combine_two_atoms_to_one_obj:
siw_spin_mean_2 = np.mean(file.get_siw(ineq=2), axis=1)
siw_spin_mean_2 = extend_orbital(siw_spin_mean_2)[None, None, None, ...]
siw_dc_spin_mean_2 = np.mean(file.get_dc(ineq=2), axis=-1)
siw_dc_spin_mean_2 = extend_orbital(siw_dc_spin_mean_2)[None, None, None, ..., None]
siw_spin_mean_2 = siw_spin_mean_2[..., niv_dmft - niv_cut : niv_dmft + niv_cut]
siw_spin_mean = 0.5 * (siw_spin_mean + siw_spin_mean_2)
del siw_spin_mean_2
siw_dc_spin_mean = 0.5 * (siw_dc_spin_mean + siw_dc_spin_mean_2)
del siw_dc_spin_mean_2

siw_spin_mean_large = np.zeros(
(1, 1, 1, 2 * siw_spin_mean.shape[3], 2 * siw_spin_mean.shape[3], siw_spin_mean.shape[-1])
)
siw_spin_mean_large[:, :, :, 0:2, 0:2, ...] = siw_spin_mean
siw_spin_mean_large[:, :, :, 2:4, 2:4, ...] = siw_spin_mean
siw_dc_spin_mean_large = np.zeros_like(siw_spin_mean_large)
siw_dc_spin_mean_large[:, :, :, 0:2, 0:2, ...] = siw_dc_spin_mean
siw_dc_spin_mean_large[:, :, :, 2:4, 2:4, ...] = siw_dc_spin_mean
sigma_dmft = SelfEnergy(siw_spin_mean_large, estimate_niv_core=True) + siw_dc_spin_mean_large
del siw_spin_mean_large, siw_dc_spin_mean_large

del giw_spin_mean, siw_spin_mean, siw_dc_spin_mean

file.close()
Expand Down Expand Up @@ -198,11 +272,17 @@ def set_hamiltonian(er_type: str, er_input: str | list, int_type: str, int_input
if int_type == "one_band_from_dmft" or int_type == "" or int_type is None:
return ham.single_band_interaction(config.lattice.interaction.udd)
elif int_type == "kanamori_from_dmft":
return ham.kanamori_interaction(
config.sys.n_bands,
return ham.kanamori_interaction_dp(
config.sys.nd_bands,
config.sys.np_bands,
config.lattice.interaction.udd,
config.lattice.interaction.upp,
config.lattice.interaction.udp,
config.lattice.interaction.jdd,
config.lattice.interaction.jpp,
config.lattice.interaction.jdp,
config.lattice.interaction.vdd,
config.lattice.interaction.vpp,
)
elif int_type == "custom":
if not isinstance(int_input, str):
Expand Down
73 changes: 64 additions & 9 deletions moldga/hamiltonian.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,26 +107,81 @@ def interaction_orbital_diagonal(self, u: float, n_bands: int = 1) -> "Hamiltoni

return self._add_interaction_term(interaction_elements)

def kanamori_interaction(self, n_bands: int, udd: float, jdd: float, vdd: float = None) -> "Hamiltonian":
def kanamori_interaction_d(self, n_bands: int, udd: float, jdd: float, vdd: float = None) -> "Hamiltonian":
"""
Adds the Kanamori interaction terms to the Hamiltonian. The interaction terms are defined by the Hubbard 'udd' (U),
the exchange interaction 'jdd' (J) and the pair hopping 'vdd' (V or sometimes U'). 'vdd' is an optional parameter, if left empty,
it is set to U-2J.
Adds the Kanamori interaction terms ONLY for d orbitals to the Hamiltonian.
The interaction terms are defined by the Hubbard 'udd' (U),
the exchange interaction 'jdd' (J) and the pair hopping 'vdd' (V or sometimes U').
'vdd' is an optional parameter, if left empty, it is set to V=U-2J.
"""
return self.kanamori_interaction_dp(nd_bands=n_bands, udd=udd, jdd=jdd, vdd=vdd)

def kanamori_interaction_p(self, n_bands: int, upp: float, jpp: float, vpp: float = None) -> "Hamiltonian":
"""
Adds the Kanamori interaction terms ONLY for p orbitals to the Hamiltonian.
The interaction terms are defined by the Hubbard 'udd' (U),
the exchange interaction 'jpp' (J) and the pair hopping 'vpp' (V or sometimes U').
'vpp' is an optional parameter, if left empty, it is set to V=U-2J.
"""
return self.kanamori_interaction_dp(np_bands=n_bands, upp=upp, jpp=jpp, vpp=vpp)

def kanamori_interaction_dp(
self,
nd_bands: int = 0,
np_bands: int = 0,
udd: float = 0.0,
upp: float = 0.0,
udp: float = 0.0,
jdd: float = 0.0,
jpp: float = 0.0,
jdp: float = 0.0,
vdd: float = None,
vpp: float = None,
) -> "Hamiltonian":
"""
Adds the full Kanamori interaction terms for d and p orbitals to the Hamiltonian.
The interaction terms are defined by the local interaction Hubbard U,
the exchange interaction J and the pair hopping V or sometimes U'.
vdd (vpp) (vdp) are optional parameters, if left empty, they are set to V=U-2J.
"""

# Default V=U-2J
if vdd is None:
vdd = udd - 2 * jdd
if vpp is None:
vpp = upp - 2 * jpp

def orb_type(i: int) -> str:
return "d" if i < nd_bands else "p" # d bands come first

def get_params(o1: int, o2: int) -> tuple[float, float, float]:
# Return correct (U, J, V) depending on orbital types of o1, o2
t1, t2 = orb_type(o1), orb_type(o2)

if t1 == "d" and t2 == "d":
return udd, jdd, vdd
elif t1 == "p" and t2 == "p":
return upp, jpp, vpp
else:
return 0, jdp, udp # udp is used as "Vdp", since there is no Udp possible that fits U_llll

r_loc = [0, 0, 0]
n_tot = nd_bands + np_bands

interaction_elements = []
for a, b, c, d in it.product(range(n_bands), repeat=4):

for a, b, c, d in it.product(range(n_tot), repeat=4):
bands = [a + 1, b + 1, c + 1, d + 1]

# choose parameters based on (a,b) pair (equivalently (c,d))
u, j, v = get_params(a, b)

if a == b == c == d: # U_{llll}
interaction_elements.append(InteractionElement(r_loc, bands, udd))
elif (a == d and b == c) or (a == b and c == d): # U_{lmml} or U_{llmm}
interaction_elements.append(InteractionElement(r_loc, bands, jdd))
interaction_elements.append(InteractionElement(r_loc, bands, u))
elif (a == d and b == c) or (a == b and c == d): # U_{lmml}, U_{llmm}
interaction_elements.append(InteractionElement(r_loc, bands, j))
elif a == c and b == d: # U_{lmlm}
interaction_elements.append(InteractionElement(r_loc, bands, vdd))
interaction_elements.append(InteractionElement(r_loc, bands, v))

return self._add_interaction_term(interaction_elements)

Expand Down
Loading
Loading