LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-22 20:34:38 +01:00
Code Issues Releases Wiki Activity
800ef7c0b4
dft_tools/python/converters/plovasp/examples/ce/test_ham_hf.py
Manuel a8c7569830 change del to with when reading hdf
2018-12-06 17:28:49 -05:00

438 lines
19 KiB
Python
Raw Blame History

#from triqs_dft_tools.sumk_dft import *
from sumk_dft import *
#from triqs_dft_tools.converters.wien2k_converter import *
from converters.vasp_converter import *
#from pytriqs.applications.impurity_solvers.hubbard_I.hubbard_solver import Solver
from hf_solver import Solver
import shutil
class TestSumkDFT(SumkDFT):
# calculate and save occupancy matrix in the Bloch basis for VASP charge denity recalculation
def calc_density_correction(self, filename='GAMMA', dm_type='wien2k'):
r"""
Calculates the charge density correction and stores it into a file.
The charge density correction is needed for charge-self-consistent DFT+DMFT calculations.
It represents a density matrix of the interacting system defined in Bloch basis
and it is calculated from the sum over Matsubara frequecies of the full GF,
..math:: N_{\nu\nu'}(k) = \sum_{i\omega_{n}} G_{\nu\nu'}(k, i\omega_{n})
The density matrix for every `k`-point is stored into a file.
Parameters
----------
filename : string
Name of the file to store the charge density correction.
Returns
-------
(deltaN, dens) : tuple
Returns a tuple containing the density matrix `deltaN` and
the corresponing total charge `dens`.
"""
assert type(filename) == StringType, "calc_density_correction: filename has to be a string!"
assert dm_type in ('vasp', 'wien2k'), "'type' must be either 'vasp' or 'wienk'"
ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO]
dens = {sp: 0.0 for sp in spn}
# Fetch Fermi weights and energy window band indices
if dm_type == 'vasp':
fermi_weights = 0
band_window = 0
if mpi.is_master_node():
with HDFArchive(self.hdf_file,'r') as ar:
fermi_weights = ar['dft_misc_input']['dft_fermi_weights']
band_window = ar['dft_misc_input']['band_window']
fermi_weights = mpi.bcast(fermi_weights)
band_window = mpi.bcast(band_window)
# Convert Fermi weights to a density matrix
dens_mat_dft = {}
for sp in spn:
dens_mat_dft[sp] = [fermi_weights[ik, ntoi[sp], :].astype(numpy.complex_) for ik in xrange(self.n_k)]
# Set up deltaN:
deltaN = {}
for sp in spn:
deltaN[sp] = [numpy.zeros([self.n_orbitals[ik,ntoi[sp]],self.n_orbitals[ik,ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray):
G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw")
for bname,gf in G_latt_iw:
deltaN[bname][ik][:, :] = G_latt_iw[bname].density()
dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density()
if dm_type == 'vasp':
# In 'vasp'-mode subtract the DFT density matrix
diag_inds = numpy.diag_indices(self.n_orbitals[ik, ntoi[bname]])
deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik]
dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real
# mpi reduce:
for bname in deltaN:
for ik in range(self.n_k):
deltaN[bname][ik] = mpi.all_reduce(mpi.world, deltaN[bname][ik], lambda x,y : x+y)
dens[bname] = mpi.all_reduce(mpi.world, dens[bname], lambda x,y : x+y)
mpi.barrier()
# now save to file:
if dm_type == 'wien2k':
if mpi.is_master_node():
if self.SP == 0:
f = open(filename,'w')
else:
f = open(filename+'up','w')
f1 = open(filename+'dn','w')
# write chemical potential (in Rydberg):
f.write("%.14f\n"%(self.chemical_potential/self.energy_unit))
if self.SP != 0: f1.write("%.14f\n"%(self.chemical_potential/self.energy_unit))
# write beta in rydberg-1
f.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit))
if self.SP != 0: f1.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit))
if self.SP == 0: # no spin-polarization
for ik in range(self.n_k):
f.write("%s\n"%self.n_orbitals[ik,0])
for inu in range(self.n_orbitals[ik,0]):
for imu in range(self.n_orbitals[ik,0]):
valre = (deltaN['up'][ik][inu,imu].real + deltaN['down'][ik][inu,imu].real) / 2.0
valim = (deltaN['up'][ik][inu,imu].imag + deltaN['down'][ik][inu,imu].imag) / 2.0
f.write("%.14f %.14f "%(valre,valim))
f.write("\n")
f.write("\n")
f.close()
elif self.SP == 1: # with spin-polarization
# dict of filename: (spin index, block_name)
if self.SO == 0: to_write = {f: (0, 'up'), f1: (1, 'down')}
if self.SO == 1: to_write = {f: (0, 'ud'), f1: (0, 'ud')}
for fout in to_write.iterkeys():
isp, sp = to_write[fout]
for ik in range(self.n_k):
fout.write("%s\n"%self.n_orbitals[ik,isp])
for inu in range(self.n_orbitals[ik,isp]):
for imu in range(self.n_orbitals[ik,isp]):
fout.write("%.14f %.14f "%(deltaN[sp][ik][inu,imu].real,deltaN[sp][ik][inu,imu].imag))
fout.write("\n")
fout.write("\n")
fout.close()
elif dm_type == 'vasp':
# assert self.SP == 0, "Spin-polarized density matrix is not implemented"
if mpi.is_master_node():
with open(filename, 'w') as f:
f.write(" %i -1 ! Number of k-points, default number of bands\n"%(self.n_k))
for sp in spn:
for ik in xrange(self.n_k):
ib1 = band_window[0][ik, 0]
ib2 = band_window[0][ik, 1]
f.write(" %i %i %i\n"%(ik + 1, ib1, ib2))
for inu in xrange(self.n_orbitals[ik, 0]):
for imu in xrange(self.n_orbitals[ik, 0]):
if self.SP == 0:
valre = (deltaN['up'][ik][inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0
valim = (deltaN['up'][ik][inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0
else:
valre = deltaN[sp][ik][inu, imu].real
valim = deltaN[sp][ik][inu, imu].imag
f.write(" %.14f %.14f"%(valre, valim))
f.write("\n")
else:
raise NotImplementedError("Unknown density matrix type: '%s'"%(dm_type))
return deltaN, dens
def calc_hamiltonian_correction(self, filename='GAMMA'):
r"""
Calculates the charge density correction and stores it into a file.
The charge density correction is needed for charge-self-consistent DFT+DMFT calculations.
It represents a density matrix of the interacting system defined in Bloch basis
and it is calculated from the sum over Matsubara frequecies of the full GF,
..math:: N_{\nu\nu'}(k) = \sum_{i\omega_{n}} G_{\nu\nu'}(k, i\omega_{n})
The density matrix for every `k`-point is stored into a file.
Parameters
----------
filename : string
Name of the file to store the charge density correction.
Returns
-------
(deltaN, dens) : tuple
Returns a tuple containing the density matrix `deltaN` and
the corresponing total charge `dens`.
"""
assert type(filename) == StringType, "calc_density_correction: filename has to be a string!"
ntoi = self.spin_names_to_ind[self.SO]
spn = self.spin_block_names[self.SO]
dens = {sp: 0.0 for sp in spn}
# Fetch Fermi weights and energy window band indices
fermi_weights = 0
band_window = 0
if mpi.is_master_node():
with HDFArchive(self.hdf_file,'r') as ar:
fermi_weights = ar['dft_misc_input']['dft_fermi_weights']
band_window = ar['dft_misc_input']['band_window']
fermi_weights = mpi.bcast(fermi_weights)
band_window = mpi.bcast(band_window)
# Set up deltaH:
deltaH = {}
for sp in spn:
deltaH[sp] = [numpy.zeros([self.n_orbitals[ik,ntoi[sp]],self.n_orbitals[ik,ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)]
ikarray = numpy.array(range(self.n_k))
for ik in mpi.slice_array(ikarray):
sigma_minus_dc = [s.copy() for s in self.Sigma_imp_iw]
sigma_minus_dc = self.add_dc('iw')
beta = self.Sigma_imp_iw[0].mesh.beta # override beta if Sigma_iw is present
n_iw = len(self.Sigma_imp_iw[0].mesh)
block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ]
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]
glist = lambda : [ GfImFreq(indices=inner,beta=beta,n_points=n_iw) for block,inner in gf_struct]
G_latt = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False)
G_latt.zero()
for icrsh in range(self.n_corr_shells):
for bname, gf in G_latt:
gf += self.upfold(ik,icrsh,bname,sigma_minus_dc[icrsh][bname],gf)
for sp in spn:
deltaH[sp][ik][:, :] = G_latt[sp](0) # Any Matsubara frequency will do
# G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw")
# for bname,gf in G_latt_iw:
# deltaN[bname][ik][:, :] = G_latt_iw[bname].density()
# dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density()
# if dm_type == 'vasp':
## In 'vasp'-mode subtract the DFT density matrix
# diag_inds = numpy.diag_indices(self.n_orbitals[ik, ntoi[bname]])
# deltaN[bname][ik][diag_inds] -= dens_mat_dft[bname][ik]
# dens[bname] -= self.bz_weights[ik] * dens_mat_dft[bname][ik].sum().real
# mpi reduce:
for bname in deltaH:
for ik in range(self.n_k):
deltaH[bname][ik] = mpi.all_reduce(mpi.world, deltaH[bname][ik], lambda x,y : x+y)
mpi.barrier()
# now save to file:
if mpi.is_master_node():
with open(filename, 'w') as f:
f.write("H %i -1 ! Number of k-points, default number of bands\n"%(self.n_k))
for sp in spn:
for ik in xrange(self.n_k):
ib1 = band_window[0][ik, 0]
ib2 = band_window[0][ik, 1]
f.write(" %i %i %i\n"%(ik + 1, ib1, ib2))
for inu in xrange(self.n_orbitals[ik, 0]):
for imu in xrange(self.n_orbitals[ik, 0]):
if self.SP == 0:
valre = (deltaH['up'][ik][inu, imu].real + deltaH['down'][ik][inu, imu].real) / 2.0
valim = (deltaH['up'][ik][inu, imu].imag + deltaH['down'][ik][inu, imu].imag) / 2.0
else:
valre = deltaH[sp][ik][inu, imu].real
valim = deltaH[sp][ik][inu, imu].imag
f.write(" %.14f %.14f"%(valre, valim))
f.write("\n")
return deltaH
def dmft_cycle():
lda_filename = 'vasp'
beta = 400
U_int = 4.00
J_hund = 0.70
Loops = 1 # Number of DMFT sc-loops
Mix = 1.0 # Mixing factor in QMC
DC_type = 0 # 0...FLL, 1...Held, 2... AMF, 3...Lichtenstein
DC_Mix = 1.0 # 1.0 ... all from imp; 0.0 ... all from Gloc
useBlocs = False # use bloc structure from LDA input
useMatrix = True # use the U matrix calculated from Slater coefficients instead of (U+2J, U, U-J)
chemical_potential_init=-0.0 # initial chemical potential
use_dudarev = True
HDFfilename = lda_filename+'.h5'
# Convert DMFT input:
# Can be commented after the first run
Converter = VaspConverter(filename=lda_filename)
Converter.convert_dft_input()
#check if there are previous runs:
previous_runs = 0
previous_present = False
if mpi.is_master_node():
with HDFArchive(HDFfilename,'a') as ar:
if 'iterations' in ar:
previous_present = True
previous_runs = ar['iterations']
else:
previous_runs = 0
previous_present = False
mpi.barrier()
previous_runs = mpi.bcast(previous_runs)
previous_present = mpi.bcast(previous_present)
# Init the SumK class
#SK=SumkDFT(hdf_file=lda_filename+'.h5',use_dft_blocks=False)
SK=TestSumkDFT(hdf_file=lda_filename+'.h5',use_dft_blocks=False)
Norb = SK.corr_shells[0]['dim']
l = SK.corr_shells[0]['l']
# Init the Hubbard-I solver:
S = Solver(beta = beta, l = l, dudarev=use_dudarev)
# DEBUG
#SK.put_Sigma(Sigma_imp=[S.Sigma])
#dH = SK.calc_hamiltonian_correction(filename='GAMMA')
# END DEBUG
chemical_potential=chemical_potential_init
# load previous data: old self-energy, chemical potential, DC correction
if (previous_present):
mpi.report("Using stored data for initialisation")
if (mpi.is_master_node()):
with HDFArchive(HDFfilename,'a') as ar:
S.Sigma <<= ar['SigmaF']
things_to_load=['chemical_potential','dc_imp']
old_data=SK.load(things_to_load)
chemical_potential=old_data[0]
SK.dc_imp=old_data[1]
S.Sigma = mpi.bcast(S.Sigma)
chemical_potential=mpi.bcast(chemical_potential)
SK.dc_imp=mpi.bcast(SK.dc_imp)
# DMFT loop:
for Iteration_Number in range(1,Loops+1):
itn = Iteration_Number + previous_runs
# put Sigma into the SumK class:
S.Sigma.zero() # !!!!
SK.put_Sigma(Sigma_imp = [ S.Sigma ])
# Compute the SumK, possibly fixing mu by dichotomy
if SK.density_required and (Iteration_Number > 1):
chemical_potential = SK.calc_mu( precision = 0.001 )
else:
mpi.report("No adjustment of chemical potential\nTotal density = %.3f"%SK.total_density(mu=chemical_potential))
# Density:
S.G <<= SK.extract_G_loc()[0]
mpi.report("Total charge of Gloc : %.6f"%S.G.total_density())
# calculated DC at the first run to have reasonable initial non-interacting atomic level positions
if ((Iteration_Number==1)and(previous_present==False)):
if use_dudarev:
dc_value_init = 0.0
else:
dc_value_init=U_int/2.0
dm=S.G.density()
SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type, use_dc_value=dc_value_init)
# calculate non-interacting atomic level positions:
# eal = SK.eff_atomic_levels()[0]
# S.set_atomic_levels( eal = eal )
# solve it:
corr_energy, dft_dc = S.solve(U_int = U_int, J_hund = J_hund, verbosity = 1)
# Now mix Sigma and G:
if ((itn>1)or(previous_present)):
if (mpi.is_master_node()and (Mix<1.0)):
with HDFArchive(HDFfilename,'r') as ar:
mpi.report("Mixing Sigma and G with factor %s"%Mix)
if ('SigmaF' in ar):
S.Sigma <<= Mix * S.Sigma + (1.0-Mix) * ar['SigmaF']
if ('GF' in ar):
S.G <<= Mix * S.G + (1.0-Mix) * ar['GF']
S.G = mpi.bcast(S.G)
S.Sigma = mpi.bcast(S.Sigma)
# after the Solver has finished, set new double counting:
# dm = S.G.density()
# if not use_dudarev:
# SK.calc_dc( dm, U_interact = U_int, J_hund = J_hund, orb = 0, use_dc_formula = DC_type )
# correlation energy calculations:
correnerg = 0.5 * (S.G * S.Sigma).total_density()
mpi.report("Corr. energy = %s"%correnerg)
# store the impurity self-energy, GF as well as correlation energy in h5
if (mpi.is_master_node()):
with HDFArchive(HDFfilename,'a') as ar:
ar['iterations'] = itn
ar['chemical_cotential%s'%itn] = chemical_potential
ar['SigmaF'] = S.Sigma
ar['GF'] = S.G
ar['correnerg%s'%itn] = correnerg
ar['DCenerg%s'%itn] = SK.dc_energ
#Save essential SumkDFT data:
things_to_save=['chemical_potential','dc_energ','dc_imp']
SK.save(things_to_save)
# print out occupancy matrix of Ce 4f
# mpi.report("Orbital densities of impurity Green function:")
# for s in dm:
# mpi.report("Block %s: "%s)
# for ii in range(len(dm[s])):
# str = ''
# for jj in range(len(dm[s])):
# if (dm[s][ii,jj].real>0):
# str += " %.4f"%(dm[s][ii,jj].real)
# else:
# str += " %.4f"%(dm[s][ii,jj].real)
# mpi.report(str)
mpi.report("Total charge of impurity problem : %.6f"%S.G.total_density())
# find exact chemical potential
#if (SK.density_required):
# SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
# SK.chemical_potential = SK.calc_mu( precision = 0.01 )
#dN, d = SK.calc_density_correction(filename='GAMMA', dm_type='vasp')
#mpi.report("Trace of Density Matrix: %s"%d)
mpi.report("Storing Hamiltonian correction GAMMA...")
SK.put_Sigma(Sigma_imp=[S.Sigma])
dH = SK.calc_hamiltonian_correction(filename='GAMMA')
# shutil.copyfile('GAMMA', 'it%i.GAMMA'%(itn))
# store correlation energy contribution to be read by Wien2ki and then included to DFT+DMFT total energy
if (mpi.is_master_node()):
with HDFArchive(HDFfilename) as ar:
itn = ar['iterations']
correnerg = ar['correnerg%s'%itn]
DCenerg = ar['DCenerg%s'%itn]
correnerg -= DCenerg[0]
f=open(lda_filename+'.qdmft','a')
f.write("%.16f\n"%correnerg)
f.close()
return corr_energy, dft_dc
if __name__ == '__main__':
dmft_cycle()
Reference in New Issue View Git Blame Copy Permalink