LCPQ
/
dft_tools
mirror of https://github.com/triqs/dft_tools
3
0
Fork 0
Code Issues Releases Wiki Activity

[feat] add option to use DLR mesh in Sumk

This commit is contained in:
Alexander Hampel 2023-12-14 11:10:59 -05:00
parent 4951b15fa4
commit 11abd8dc06
3 changed files with 70 additions and 15 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
python/triqs_dft_tools/sumk_dft.py
View File

@ -58,7 +58,7 @@ class SumkDFT(object):
The value of magnetic field to add to the DFT Hamiltonian.
The contribution -h_field*sigma is added to diagonal elements of the Hamiltonian.
It cannot be used with the spin-orbit coupling on; namely h_field is set to 0 if self.SO=True.
mesh: MeshImFreq or MeshReFreq, optional. Frequency mesh of Sigma.
mesh: MeshImFreq, MeshDLRImFreq, or MeshReFreq, optional. Frequency mesh of Sigma.
beta : real, optional
Inverse temperature. Used to construct imaginary frequency if mesh is not given.
n_iw : integer, optional
@ -115,6 +115,9 @@ class SumkDFT(object):
self.mesh_values = np.linspace(self.mesh(self.mesh.first_index()),
self.mesh(self.mesh.last_index()),
len(self.mesh))
elif isinstance(mesh, MeshDLRImFreq):
self.mesh = mesh
self.mesh_values = np.array([iwn.value for iwn in mesh.values()])
elif isinstance(mesh, MeshReFreq):
self.mesh = mesh
self.mesh_values = np.linspace(self.mesh.w_min, self.mesh.w_max, len(self.mesh))
@ -563,6 +566,8 @@ class SumkDFT(object):
mesh = Sigma_imp[0].mesh
if isinstance(mesh, MeshImFreq):
mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
elif isinstance(mesh, MeshDLRImFreq):
mesh_values = np.array([iwn.value for iwn in mesh.values()])
else:
mesh_values = np.linspace(mesh.w_min, mesh.w_max, len(mesh))
else:
@ -570,9 +575,11 @@ class SumkDFT(object):
mesh_values = self.mesh_values
elif not mesh is None:
assert isinstance(mesh, MeshReFreq) or isinstance(mesh, MeshImFreq), "mesh must be a triqs MeshReFreq or MeshImFreq"
assert isinstance(mesh, (MeshReFreq, MeshDLRImFreq, MeshImFreq)), "mesh must be a triqs MeshReFreq or MeshImFreq"
if isinstance(mesh, MeshImFreq):
mesh_values = np.linspace(mesh(mesh.first_index()), mesh(mesh.last_index()), len(mesh))
elif isinstance(mesh, MeshDLRImFreq):
mesh_values = np.array([iwn.value for iwn in mesh.values()])
else:
mesh_values = np.linspace(mesh.w_min, mesh.w_max, len(mesh))
else:
@ -586,12 +593,8 @@ class SumkDFT(object):
gf_struct = [(spn[isp], block_structure[isp])
for isp in range(self.n_spin_blocks[self.SO])]
block_ind_list = [block for block, inner in gf_struct]
if isinstance(mesh, MeshImFreq):
glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)])
for block, inner in gf_struct]
else:
glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)])
for block, inner in gf_struct]
glist = lambda: [Gf(mesh=mesh, target_shape=[len(inner),len(inner)])
for block, inner in gf_struct]
G_latt = BlockGf(name_list=block_ind_list,
block_list=glist(), make_copies=False)
G_latt.zero()
@ -603,7 +606,7 @@ class SumkDFT(object):
for ibl, (block, gf) in enumerate(G_latt):
ind = ntoi[spn[ibl]]
n_orb = self.n_orbitals[ik, ind]
if isinstance(mesh, MeshImFreq):
if isinstance(mesh, (MeshImFreq, MeshDLRImFreq)):
gf.data[:, :, :] = (idmat[ibl] * (mesh_values[:, None, None] + mu + self.h_field*(1-2*ibl))
- self.hopping[ik, ind, 0:n_orb, 0:n_orb])
else:
@ -647,7 +650,10 @@ class SumkDFT(object):
assert len(Sigma_imp) == self.n_corr_shells,\
"put_Sigma: give exactly one Sigma for each corr. shell!"
if isinstance(self.mesh, MeshImFreq) and all(isinstance(gf.mesh, MeshImFreq) and isinstance(gf, Gf) and gf.mesh == self.mesh for bname, gf in Sigma_imp[0]):
if (isinstance(self.mesh, (MeshImFreq, MeshDLRImFreq)) and
all(isinstance(gf.mesh, (MeshImFreq, MeshDLRImFreq)) and
isinstance(gf, Gf) and
gf.mesh == self.mesh for bname, gf in Sigma_imp[0])):
# Imaginary frequency Sigma:
self.Sigma_imp = [self.block_structure.create_gf(ish=icrsh, mesh=Sigma_imp[icrsh].mesh, space='sumk')
for icrsh in range(self.n_corr_shells)]

8
test/python/calc_mu.py
View File

@ -31,7 +31,8 @@ from triqs_dft_tools.sumk_dft import *
class test_solver(unittest.TestCase):
def setUp(self):
self.iw_mesh = MeshImFreq(beta=40, S='Fermion', n_iw=300)
self.iw_mesh = MeshImFreq(beta=40, statistic='Fermion', n_iw=300)
self.dlr_mesh = MeshDLRImFreq(beta=40, statistic='Fermion', w_max=10, eps=1e-10)
self.w_mesh = MeshReFreq(n_w=1001, window=(-3,3))
# magic reference value for the Wien2k SVO t2g example
self.ref_mu = 0.281
@ -42,6 +43,11 @@ class test_solver(unittest.TestCase):
mu = sumk.calc_mu(method='dichotomy', precision=0.001, delta=0.1)
self.assertTrue(abs(self.ref_mu - mu) < 0.01)
def test_dichotomy_dlr(self):
sumk = SumkDFT('SrVO3.ref.h5', mesh=self.dlr_mesh)
mu = sumk.calc_mu(method='dichotomy', precision=0.001, delta=0.1)
self.assertTrue(abs(self.ref_mu - mu) < 0.01)
def test_dichotomy_real(self):
sumk = SumkDFT('SrVO3.ref.h5', mesh=self.w_mesh)
mu = sumk.calc_mu(method='dichotomy', precision=0.001, delta=0.1, broadening = 0.01, beta=1000)

51
test/python/srvo3_Gloc.py
View File

@ -20,18 +20,25 @@
################################################################################
from h5 import *
from triqs.gf import *
from triqs.gf import MeshImFreq, MeshDLRImFreq
from triqs_dft_tools.sumk_dft import *
from triqs_dft_tools.converters.wien2k import *
from triqs.operators.util import set_operator_structure
from triqs.utility.comparison_tests import *
from triqs.utility.h5diff import h5diff
import time
# Basic input parameters
beta = 40
n_iw = 1025
# Init the SumK class
SK=SumkDFT(hdf_file='SrVO3.ref.h5',use_dft_blocks=True)
# classic full Matsubara mesh
mpi.report(f"{'#'*12}\nregular Matsubara mesh test\n")
# Init the SumK class (reference data with n_iw=1025)
iw_mesh = MeshImFreq(n_iw=n_iw,beta=beta, statistic='Fermion')
SK=SumkDFT(hdf_file='SrVO3.ref.h5',mesh=iw_mesh,use_dft_blocks=True)
num_orbitals = SK.corr_shells[0]['dim']
l = SK.corr_shells[0]['l']
@ -39,15 +46,51 @@ spin_names = ['down','up']
orb_hybridized = False
gf_struct = set_operator_structure(spin_names,num_orbitals,orb_hybridized)
glist = [ GfImFreq(target_shape=(bl_size,bl_size),beta=beta) for bl, bl_size in gf_struct]
glist = [ Gf(target_shape=(bl_size,bl_size),mesh=iw_mesh) for bl, bl_size in gf_struct]
Sigma_iw = BlockGf(name_list = [bl for bl, bl_size in gf_struct], block_list = glist, make_copies = False)
SK.set_Sigma([Sigma_iw])
if mpi.is_master_node():
start_time = time.time()
Gloc = SK.extract_G_loc()
if mpi.is_master_node():
mpi.report(f'extract_G_loc time: {(time.time()-start_time)*1000:.1f} msec')
if mpi.is_master_node():
with HDFArchive('srvo3_Gloc.out.h5','w') as ar:
ar['Gloc'] = Gloc[0]
if mpi.is_master_node():
h5diff("srvo3_Gloc.out.h5","srvo3_Gloc.ref.h5")
mpi.report(f"{'#'*12}\n")
# DLR Matsubara mesh
mpi.report(f"{'#'*12}\nDLR Matsubara mesh test\n")
dlr_mesh = MeshDLRImFreq(beta=beta, statistic='Fermion', w_max=10, eps=1e-10)
SK=SumkDFT(hdf_file='SrVO3.ref.h5',mesh=dlr_mesh,use_dft_blocks=True)
glist_dlr = [ Gf(target_shape=(bl_size,bl_size),mesh=dlr_mesh) for bl, bl_size in gf_struct]
Sigma_dlr = BlockGf(name_list = [bl for bl, bl_size in gf_struct], block_list = glist_dlr, make_copies = False)
SK.set_Sigma([Sigma_dlr])
if mpi.is_master_node():
start_time = time.time()
Gloc_dlr_iw = SK.extract_G_loc()
if mpi.is_master_node():
mpi.report(f'extract_G_loc time: {(time.time()-start_time)*1000:.1f} msec')
with HDFArchive('srvo3_Gloc.out.h5','a') as ar:
ar['Gloc_dlr'] = make_gf_imfreq(make_gf_dlr(Gloc_dlr_iw[0]),n_iw=n_iw)
# get full Giw and compare
Gloc_iw_full = make_gf_imfreq(make_gf_dlr(Gloc_dlr_iw[0]),n_iw=n_iw)
assert_block_gfs_are_close(Gloc[0], Gloc_iw_full)
mpi.report(f"{'#'*12}\n")