Skip to content

Commit

Permalink
Merge pull request #86 from issp-center-dev/dc_type
Browse files Browse the repository at this point in the history 
A new parameter [system] dc_type is introduced
  • Loading branch information
@shinaoka
shinaoka authored Nov 20, 2020
2 parents e982882 + c748310 commit 97ac5e8
Show file tree
Hide file tree
Showing 4 changed files with 116 additions and 78 deletions.
156 changes: 88 additions & 68 deletions python/dmft_core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -211,6 +211,64 @@ def solve_impurity_model(solver_name, solver_params, mpirun_command, basis_rot,
return r


def calc_dc(dc_type, u_mat, dens_mat, spin_block_names, use_spin_orbit):

# dim_tot is the dimension of spin x orbit for SO = 1 or that of orbital for SO=0
# dim_tot = self._dim_sh[ish]
# TODO: num_orb can be deleted, if inverse transformation of .tools.to_spin_full_U_matrix is implemented}
num_orb = int(u_mat.shape[0] / 2)

dc_imp_sh = {} # This will be returned
#
# Hartree-Fock contribution
#
if dc_type == "HF_DFT" or dc_type == "HF_imp":
if use_spin_orbit:
dim_tot = dens_mat["ud"].shape[0] # 2 * num_orb
dc_imp_sh["ud"] = numpy.zeros((dim_tot, dim_tot), numpy.complex_)
dc_imp_sh["ud"] += numpy.einsum("ijkl, jl->ik", u_mat, dens_mat["ud"], optimize=True) # Hartree
dc_imp_sh["ud"] -= numpy.einsum("ijkl, jk->il", u_mat, dens_mat["ud"], optimize=True) # Fock
else:
for sp1 in spin_block_names:
u_mat_reduce = u_mat[0:num_orb, 0:num_orb, 0:num_orb, 0:num_orb] # TODO: make a function
dens_mat_tot = sum(dens_mat.values()) # spin-sum of density matrix

dc_imp_sh[sp1] = numpy.zeros((num_orb, num_orb), numpy.complex_)
dc_imp_sh[sp1] += numpy.einsum("ijkl, jl->ik", u_mat_reduce, dens_mat_tot, optimize=True) # Hartree
dc_imp_sh[sp1] -= numpy.einsum("ijkl, jk->il", u_mat_reduce, dens_mat[sp1], optimize=True) # Fock
#
# Fully-Localized Limit (FLL)
#
elif dc_type == "FLL":
if use_spin_orbit:
raise NotImplementedError
else:
u_sum = numpy.einsum("ijij", u_mat[0:num_orb, 0:num_orb, 0:num_orb, 0:num_orb], optimize=True)
j_sum = numpy.einsum("ijji", u_mat[0:num_orb, 0:num_orb, 0:num_orb, 0:num_orb], optimize=True)

u_ave = u_sum / num_orb**2 # U
u_j_ave = (u_sum - j_sum) / (num_orb*(num_orb-1)) # U-J
j_ave = u_ave - u_j_ave # J

n_sp = {sp: numpy.sum(numpy.diag(dens_mat[sp])).real for sp in spin_block_names}
n_tot = sum(n_sp.values())

print("\n FLL formula")
print(" U_ave = {0:.3f}".format(u_ave))
print(" J_ave = {0:.3f}".format(j_ave))
print(" n_sp = {}".format(n_sp))
print(" n_tot = {}".format(n_tot))

dc_imp_sh = {}
for sp in spin_block_names:
dc = u_ave * (n_tot - 0.5) - j_ave * (n_sp[sp] - 0.5)
dc_imp_sh[sp] = numpy.identity(num_orb) * dc
else:
raise ValueError("Here should not be reached")

return dc_imp_sh


class DMFTCoreSolver(object):
def __init__(self, seedname, params, output_file='', output_group='dmft_out', read_only=False, restart=False):
"""
Expand Down Expand Up @@ -591,92 +649,49 @@ def set_dc_imp(self, dm_sh):
"""

dc_type = self._params['system']['dc_type']
print("\n set DC (dc_type = {})".format(dc_type))

def print_matrix(mat):
dim = mat.shape[0]
for i1 in range(dim):
print(" ", end="")
for i2 in range(dim):
print("{0:.3f} ".format(mat[i1, i2]), end="")
print("")

# Loop over inequivalent shells
_dc_imp = []
for ish in range(self._n_inequiv_shells):
u_mat = self._Umat[self._sk.inequiv_to_corr[ish]]
dens_mat = dm_sh[self._sk.inequiv_to_corr[ish]]

# dim_tot is the dimension of spin x orbit for SO = 1 or that of orbital for SO=0
dim_tot = self._dim_sh[ish]
num_orb = int(u_mat.shape[0] / 2)
# TODO: num_orb can be deleted, if inverse transformation of .tools.to_spin_full_U_matrix is implemented}

dens_mat = dm_sh[self._sk.inequiv_to_corr[ish]]

print("")
print(" DC for inequivalent shell {0}".format(ish))
# print U matrix, J matrix, density matrix
print("\n DC for inequivalent shell {0}".format(ish))
print("\n 2-index U:".format(ish))
for i1 in range(num_orb):
print(" ", end="")
for i2 in range(num_orb):
print("{0:.3f} ".format(u_mat[i1, i2, i1, i2]), end="")
print("")
print_matrix(numpy.einsum("ijij->ij", u_mat[0:num_orb, 0:num_orb, 0:num_orb, 0:num_orb]))
print("\n 2-index J:".format(ish))
for i1 in range(num_orb):
print(" ", end="")
for i2 in range(num_orb):
print("{0:.3f} ".format(u_mat[i1, i2, i2, i1]), end="")
print("")

print_matrix(numpy.einsum("ijji->ij", u_mat[0:num_orb, 0:num_orb, 0:num_orb, 0:num_orb]))
print("\n Local Density Matrix:".format(ish))
for sp1 in self._spin_block_names:
print(" Spin {0}".format(sp1))
for i1 in range(num_orb):
print(" ", end="")
for i2 in range(num_orb):
print("{0:.3f} ".format(dens_mat[sp1][i1, i2]), end="")
print("")
print_matrix(dens_mat[sp1][0:num_orb, 0:num_orb])

if self._use_spin_orbit:
dc_imp_sh = {}
dc_imp_sh["ud"] = numpy.zeros((dim_tot, dim_tot), numpy.complex_)
for s1, i1, s2, i2 in product(range(2), range(num_orb), range(2), range(num_orb)):
#
# Hartree
#
dc_imp_sh["ud"][i1 + s1 * num_orb, i2 + s1 * num_orb] += numpy.sum(
u_mat[i1, 0:num_orb, i2, 0:num_orb] * dens_mat["ud"][s2 * num_orb:s2 * num_orb + num_orb,
s2 * num_orb:s2 * num_orb + num_orb]
)
#
# Exchange
#
dc_imp_sh["ud"][i1 + s1 * num_orb, i2 + s2 * num_orb] -= numpy.sum(
u_mat[i1, 0:num_orb, 0:num_orb, i2]
* dens_mat["ud"][s2 * num_orb:s2 * num_orb + num_orb, s1 * num_orb:s1 * num_orb + num_orb]
)
_dc_imp.append(dc_imp_sh)
else:
dc_imp_sh = {}
for sp1 in self._spin_block_names:
dc_imp_sh[sp1] = numpy.zeros((num_orb, num_orb), numpy.complex_)
for i1, i2 in product(range(num_orb), repeat=2):
#
# Hartree
#
for sp2 in self._spin_block_names:
dc_imp_sh[sp1][i1, i2] += \
numpy.sum(u_mat[i1, 0:num_orb, i2, 0:num_orb] * dens_mat[sp2][:, :])
#
# Exchange
#
dc_imp_sh[sp1][i1, i2] += \
- numpy.sum(u_mat[i1, 0:num_orb, 0:num_orb, i2] * dens_mat[sp1][:, :])
_dc_imp.append(dc_imp_sh)
# calculated DC
_dc_imp.append(calc_dc(dc_type, u_mat, dens_mat, self.spin_block_names, self._use_spin_orbit))

# print DC self energy
print("\n DC Self Energy:")
for sp1 in self._spin_block_names:
print(" Spin {0}".format(sp1))
for i1 in range(dim_tot):
print(" ", end="")
for i2 in range(dim_tot):
print("{0:.3f} ".format(_dc_imp[ish][sp1][i1, i2]), end="")
print("")
print_matrix(_dc_imp[ish][sp1])
print("")

# Now copy dc_imp to all correlated shells
self._dc_imp = [_dc_imp[self._sk.corr_to_inequiv[icrsh]] for icrsh in range(self._n_corr_shells)]


# Now copy dc_imp to all correlated shells
self._dc_imp = [_dc_imp[self._sk.corr_to_inequiv[icrsh]] for icrsh in range(self._n_corr_shells)]

def do_steps(self, max_step):
"""
Expand All @@ -692,6 +707,7 @@ def do_steps(self, max_step):

previous_present = self._previous_runs > 0
with_dc = self._params['system']['with_dc']
dc_type = self._params['system']['dc_type']
sigma_mix = self._params['control']['sigma_mix'] # Mixing factor of Sigma after solution of the AIM
output_group = self._output_group

Expand Down Expand Up @@ -779,7 +795,6 @@ def quantities_to_check():
symmetrize_spin(new_Gimp_iw[ish])
symmetrize_spin(new_Sigma_iw[ish])


# Update Sigma_iw and Gimp_iw.
# Mix Sigma if requested.
if iteration_number > 1 or previous_present:
Expand All @@ -795,6 +810,11 @@ def quantities_to_check():
if len(symm_generators[ish]) > 0:
self._sh_quant[ish].Sigma_iw << symmetrize(self._sh_quant[ish].Sigma_iw, symm_generators[ish])

# update DC correction
if with_dc and dc_type == "HF_imp":
dm_imp = [new_Gimp_iw[ish].density() for ish in range(self._n_inequiv_shells)]
self.set_dc_imp(dm_imp)

self._quant_to_save_history.update(chemical_potential=self._chemical_potential,
dc_imp=self._dc_imp,
dc_energ=self._dc_energ)
Expand Down
33 changes: 24 additions & 9 deletions python/impurity_solvers/alps_cthyb_seg.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,6 +131,12 @@ def func(u):
print("", file=f)


def set_tail(g_block):
for bname, gf in g_block:
gf.tail.zero()
gf.tail[1] = numpy.identity(gf.N1)


class ALPSCTHYBSEGSolver(SolverBase):

def __init__(self, beta, gf_struct, u_mat, n_iw=1025):
Expand Down Expand Up @@ -267,15 +273,24 @@ def _read(key):

# (3) Copy results into
# self._Sigma_iw
swdata = numpy.zeros((2, self.n_orb, self.n_iw), dtype=complex)
with HDFArchive('sim.h5', 'r') as f:
for orb in range(self.n_orb):
for spin in range(2):
swdata_array = f["S_omega"][str(orb*2+spin)]["mean"]["value"]
assert swdata_array.dtype == numpy.complex
assert swdata_array.shape == (self.n_iw,)
swdata[spin, orb, :] = swdata_array
assign_from_numpy_array(self._Sigma_iw, swdata, self.block_names)
# self._Gimp_iw

def set_blockgf_from_h5(sigma, group):
swdata = numpy.zeros((2, self.n_orb, self.n_iw), dtype=complex)
with HDFArchive('sim.h5', 'r') as f:
for orb in range(self.n_orb):
for spin in range(2):
swdata_array = f[group][str(orb*2+spin)]["mean"]["value"]
assert swdata_array.dtype == numpy.complex
assert swdata_array.shape == (self.n_iw,)
swdata[spin, orb, :] = swdata_array
assign_from_numpy_array(sigma, swdata, self.block_names)

set_blockgf_from_h5(self._Sigma_iw, "S_omega")
set_blockgf_from_h5(self._Gimp_iw, "G_omega")

if triqs_major_version == 1:
set_tail(self._Gimp_iw)

# self.quant_to_save['nn_equal_time']
nn_equal_time = numpy.zeros((2*self.n_orb, 2*self.n_orb), dtype=float)
Expand Down
2 changes: 2 additions & 0 deletions python/impurity_solvers/null_solver.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,8 @@ def solve(self, rot, mpirun_command, params_kw):
# (3) Copy results into
# self._Sigma_iw
# self._Gimp_iw
# Only _Sigma_iw is used for updating Delta(iw), but _Gimp_iw is also required.
# _Gimp_iw used for, e.g., checking the self-consistency condition.

def name(self):
return "null"
3 changes: 2 additions & 1 deletion python/program_options.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,7 +76,8 @@ def create_parser(target_sections=None):
parser.add_option("system", "mu", float, 0.0, "Initial chemical potential.")
parser.add_option("system", "prec_mu", float, 0.0001,
"Threshold for calculating chemical potential with the bisection method.")
parser.add_option("system", "with_dc", bool, False, "Whether or not use double counting correction (See below)")
parser.add_option("system", "with_dc", bool, False, "Whether or not use double-counting correction (See below)")
parser.add_option("system", "dc_type", str, 'HF_DFT', "Chosen from 'HF_DFT' (default), 'HF_imp', 'FLL'")

# [impurity_solver]
parser.add_option("impurity_solver", "name", str, 'null',
Expand Down

0 comments on commit 97ac5e8

Please sign in to comment.