Skip to content

WIP: Improvements for the DBSE tutorial #41

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
hmenke wants to merge 2 commits into
base: unstable
Choose a base branch
from
Draft

WIP: Improvements for the DBSE tutorial #41

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7887c66
Improvements for the DBSE tutorial
hmenke Sep 19, 2023
f1e5219
Find chemical potential via lattice Dyson equation
hmenke Sep 25, 2023
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
116 changes: 81 additions & 35 deletions doc/user_guide/dmft_susceptibility_dbse/common.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ def setup_dmft_calculation(p):
p.iter = 0

# -- Block structure of GF

p.num_orbitals = 3
p.spin_names = ('up','do')
p.orb_names = list(range(p.num_orbitals))
Expand All @@ -67,27 +67,38 @@ def setup_dmft_calculation(p):
p.fundamental_operators = [
c_operator(f'{s}_{o}', 0) for s, o in
itertools.product(p.spin_names, p.orb_names)]

# -- Local Hamiltonian

p.J = p.U * p.J_over_U
KanMat1, KanMat2 = U_matrix_kanamori(p.num_orbitals, p.U, p.J)
p.solve.h_int = h_int_kanamori(p.spin_names, p.num_orbitals, KanMat1, KanMat2, p.J, False)
KanMat1, KanMat2 = U_matrix_kanamori(
n_orb=p.num_orbitals,
U_int=p.U,
J_hund=p.J,
)
p.solve.h_int = h_int_kanamori(
spin_names=p.spin_names,
n_orb=p.num_orbitals,
U=KanMat1,
Uprime=KanMat2,
J_hund=p.J,
off_diag=False,
)

# -- Sr2RuO4 Wannier90 model from GPAW

from tight_binding_model import tight_binding_model
H = tight_binding_model()

p.kmesh = H.get_kmesh(n_k = (p.n_k, p.n_k, p.n_k))
p.e_k = H.fourier(p.kmesh)

# -- Initial zero guess for the self-energy
p.sigma_w = Gf(mesh=MeshImFreq(p.init.beta, 'Fermion', p.init.n_iw), target_shape=[6, 6])

p.sigma_w = Gf(mesh=MeshImFreq(p.init.beta, 'Fermion', p.init.n_iw), target_shape=[2 * p.num_orbitals, 2 * p.num_orbitals])
p.sigma_w.zero()
p.sigma_w << +p.mu * np.eye(2*p.num_orbitals)

return p


Expand All @@ -101,7 +112,7 @@ def target_function(mu, p0, ps):

if len(ps) != 0: p0 = copy.deepcopy(ps[-1])
p0.mu = mu

ps += solve_self_consistent_dmft(p0, verbose=False)

if filename:
Expand All @@ -112,7 +123,7 @@ def target_function(mu, p0, ps):
with HDFArchive(tmp_filename, 'w') as a:
#a['ps'] = ParameterCollections(ps)
a['ps'] = ps

p = ps[-1]
mpi.report(f'--> solve_self_consistent_dmft_fix_N: target_function ' + \
f'iter = {p.iter} mu = {p.mu}')
Expand All @@ -138,6 +149,43 @@ def target_function(mu, p0, ps):
return ps


def set_chemical_potential(p0):
step = 0
p = copy.deepcopy(p0)
def f(mu):
nonlocal p, step
step += 1
p.g_w = lattice_dyson_g_w(mu, p.e_k, p.sigma_w)
p.N = p.g_w.total_density().real
mpi.report(f"--> Step: {step:3d} || mu = {mu} || Density: {p.N}")
return p.N - p.N_target

from scipy.optimize import toms748
p.mu = toms748(f, p.mu_min, p.mu_max, xtol=1e-5, maxiter=100)
mpi.report(f"--> Result: || mu = {p.mu} || Density: {p.N}")
return p


def solve_self_consistent_dmft_fix_mu(p, verbose=False):

ps = []
for dmft_iter in range(p.sc_iter_max):
mpi.report(f'--> DMFT Iteration: {p.iter}')
p = set_chemical_potential(p)
p = dmft_self_consistent_step(p, verbose=verbose)
ps.append(p)
mpi.report(f'--> DMFT Convergence: dG = {p.dG:2.2E}')
mpi.report(f'--> Density: N = {p.N:2.2E}')
mpi.report(f'--> Chempot: mu = {p.mu:2.2E}')
mpi.report(f'--> Densmat: rho_up = \n{p.rho.real[:p.num_orbitals,:p.num_orbitals]}')
mpi.report(f'--> Densmat: rho_do = \n{p.rho.real[p.num_orbitals:,p.num_orbitals:]}')
if p.dG < p.G_tol: break

if p.dG > p.G_tol: mpi.report('--> Warning: DMFT Not converged!')
else: mpi.report(f'--> DMFT Converged: dG = {p.dG:2.2E}')
return ps


def solve_self_consistent_dmft(p, verbose=False):

ps = []
Expand All @@ -148,8 +196,8 @@ def solve_self_consistent_dmft(p, verbose=False):
mpi.report(f'--> DMFT Convergence: dG = {p.dG:2.2E}')
mpi.report(f'--> Density: N = {p.N:2.2E}')
mpi.report(f'--> Chempot: mu = {p.mu:2.2E}')
mpi.report(f'--> Densmat: rho_up = \n{p.rho.real[:3,:3]}')
mpi.report(f'--> Densmat: rho_do = \n{p.rho.real[3:,3:]}')
mpi.report(f'--> Densmat: rho_up = \n{p.rho.real[:p.num_orbitals,:p.num_orbitals]}')
mpi.report(f'--> Densmat: rho_do = \n{p.rho.real[p.num_orbitals:,p.num_orbitals:]}')
if p.dG < p.G_tol: break

if p.dG > p.G_tol: mpi.report('--> Warning: DMFT Not converged!')
Expand All @@ -163,18 +211,18 @@ def dmft_self_consistent_step(p, verbose=False):
p.iter += 1

p.g_w = lattice_dyson_g_w(p.mu, p.e_k, p.sigma_w)

p.rho = p.g_w.density()
p.N = np.sum(np.diagonal(p.rho)).real

# -- Solve Impurity problem

p.g0_w = p.g_w.copy()
p.g0_w << inverse(inverse(p.g_w) + p.sigma_w)

cthyb = triqs_cthyb.Solver(**p.init.dict())

BlockGf_from_Gf_matrix_valued(cthyb.G0_iw, p.g0_w)
BlockGf_from_Gf_matrix_valued(cthyb.G0_iw, p.g0_w, p)

solve = copy.deepcopy(p.solve)
solve.n_cycles = int(np.round(solve.n_cycles / mpi.size))
Expand All @@ -184,7 +232,7 @@ def dmft_self_consistent_step(p, verbose=False):
p.G0_w = cthyb.G0_iw.copy()
p.G_tau_raw = cthyb.G_tau
p.Delta_infty = cthyb.Delta_infty

G_tau = fit_dlr(p.G_tau_raw, p)

p.dG = np.max(np.abs(BlockGf_data(p.G_tau - G_tau))) \
Expand All @@ -195,11 +243,11 @@ def dmft_self_consistent_step(p, verbose=False):
p.G_w = p.G0_w.copy()
p.G_w << Fourier(p.G_tau)

p.Sigma_w = p.G0_w.copy()
p.Sigma_w = p.G0_w.copy()
p.Sigma_w << inverse(p.G0_w) - inverse(p.G_w)
Gf_matrix_valued_from_BlockGf(p.sigma_w, p.Sigma_w)

Gf_matrix_valued_from_BlockGf(p.sigma_w, p.Sigma_w, p)

return p


Expand All @@ -209,7 +257,7 @@ def fit_dlr(G_tau, p):

G_tau = G_tau.copy()
cmesh = MeshDLR(p.init.beta, 'Fermion', w_max=p.w_max, eps=p.eps)
print(cmesh)
mpi.report(cmesh)

block_mat = np.diag([1, 1, 2, 1, 1, 2])
sym = BlockSymmetrizer(len(cmesh), block_mat)
Expand All @@ -223,13 +271,13 @@ def fit_dlr(G_tau, p):
verbose=mpi.is_master_node(),
)

G_tau_mat = Gf(mesh=G_tau.mesh, target_shape=[6, 6])
Gf_matrix_valued_from_BlockGf(G_tau_mat, G_tau)
G_tau_mat = Gf(mesh=G_tau.mesh, target_shape=[2 * p.num_orbitals, 2 * p.num_orbitals])
Gf_matrix_valued_from_BlockGf(G_tau_mat, G_tau, p)

H_delta = Operator()
for i in range(6):
for i in range(2 * p.num_orbitals):
op = p.fundamental_operators[i]
H_delta += p.Delta_infty[i][0,0].real * dagger(op) * op
H_delta += p.Delta_infty[i][0,0].real * dagger(op) * op

H = p.solve.h_int + H_delta

Expand All @@ -238,25 +286,23 @@ def fit_dlr(G_tau, p):
from triqs.gf.gf_factories import make_gf_imtime
G_mat_fit = make_gf_imtime(G_c_mat, len(G_tau_mat.mesh))

BlockGf_from_Gf_matrix_valued(G_tau, G_mat_fit)
BlockGf_from_Gf_matrix_valued(G_tau, G_mat_fit, p)
return G_tau


def get_matrix_block_pairs():
blocks = [ f'{s}_{o}' for s, o in itertools.product( ['up', 'do'], range(3) ) ]
def get_matrix_block_pairs(p):
blocks = [ f'{s}_{o}' for s, o in itertools.product(p.spin_names, p.orb_names) ]
matrix_block_pairs = [ (idx, bidx) for idx, bidx in enumerate(blocks) ]
return matrix_block_pairs


def Gf_matrix_valued_from_BlockGf(G_mat, G_block):
for idx, bidx in get_matrix_block_pairs():
def Gf_matrix_valued_from_BlockGf(G_mat, G_block, p):
for idx, bidx in get_matrix_block_pairs(p):
G_mat[idx, idx] << G_block[bidx][0, 0]
return G_mat


def BlockGf_from_Gf_matrix_valued(G_block, G_mat):
for idx, bidx in get_matrix_block_pairs():
def BlockGf_from_Gf_matrix_valued(G_block, G_mat, p):
for idx, bidx in get_matrix_block_pairs(p):
G_block[bidx][0,0] << G_mat[idx, idx]
return G_block