mirror of
https://github.com/triqs/dft_tools
synced 2024-07-17 16:33:46 +02:00
ad8c4e75fe
Adding Elk-TRIQS interface (first iteration) This interface reads in Elk's ground-state files / projectors generated by a specific Elk interface code version (https://github.com/AlynJ/Elk_interface-TRIQS). The interface can perform charge-self consistent DFT+DMFT calculations using the aforementioned Elk code version, including spin orbit-coupling. Hence, this is the first open source interfaced DFT code to triqs with FCSC support. The commit includes detailed documentation and tutorials on how to use this interface. Moreover, further new post-processing routines are added for Fermi surface plots and spectral functions (A(w)) from the elk inputs. The interface was tested by A. James and A. Hampel. However, this is the first iteration of the interface and should be used with care. Please report all bugs. List of changes: --------------- - sumk.py: added cacluation of charge density correction for elk (dm_type='elk'). - sumk_dft_tools.py: added new post-processing functions to calculate the Fermi surface and A(w) from the Elk inputs. - documentation and tutorial files amended for this interface. - added python tests for the Elk converter.
166 lines
6.2 KiB
Python
166 lines
6.2 KiB
Python
import triqs.utility.mpi as mpi
|
|
from triqs.operators.util import *
|
|
from h5 import HDFArchive
|
|
from triqs_cthyb import *
|
|
from triqs.gf import *
|
|
from triqs_dft_tools.sumk_dft import *
|
|
from triqs_dft_tools.converters.elk import *
|
|
|
|
dft_filename='SrVO3'
|
|
beta = 40
|
|
loops = 2 # Number of DMFT sc-loops
|
|
sigma_mix = 1.0 # Mixing factor of Sigma after solution of the AIM
|
|
use_blocks = True # use bloc structure from DFT input
|
|
prec_mu = 0.0001
|
|
h_field = 0.0
|
|
|
|
|
|
# Solver parameters
|
|
p = {}
|
|
p["max_time"] = -1
|
|
p["length_cycle"] = 200
|
|
p["n_warmup_cycles"] = 10000
|
|
p["n_cycles"] = 1000000
|
|
p["perform_tail_fit"] = True
|
|
p["fit_max_moment"] = 4
|
|
p["fit_min_n"] = 30
|
|
p["fit_max_n"] = 60
|
|
|
|
# If conversion step was not done, we could do it here. Uncomment the lines it you want to do this.
|
|
Converter = ElkConverter(filename=dft_filename, repacking=True)
|
|
Converter.convert_dft_input()
|
|
mpi.barrier()
|
|
|
|
previous_runs = 0
|
|
previous_present = False
|
|
if mpi.is_master_node():
|
|
with HDFArchive(dft_filename+'.h5','a') as f:
|
|
if 'dmft_output' in f:
|
|
ar = f['dmft_output']
|
|
if 'iterations' in ar:
|
|
previous_present = True
|
|
previous_runs = ar['iterations']
|
|
else:
|
|
f.create_group('dmft_output')
|
|
previous_runs = mpi.bcast(previous_runs)
|
|
previous_present = mpi.bcast(previous_present)
|
|
|
|
SK=SumkDFT(hdf_file=dft_filename+'.h5',use_dft_blocks=use_blocks,h_field=h_field)
|
|
|
|
n_orb = SK.corr_shells[0]['dim']
|
|
l = SK.corr_shells[0]['l']
|
|
spin_names = ["up","down"]
|
|
orb_names = [i for i in range(n_orb)]
|
|
|
|
## KANAMORI DENSITY-DENSITY (for full Kanamori use h_int_kanamori)
|
|
# Define interaction paramters, DC and Hamiltonian
|
|
U = 4.0
|
|
J = 0.65
|
|
dc_type = 1 # DC type: 0 FLL, 1 Held, 2 AMF
|
|
# Construct U matrix for density-density calculations
|
|
Umat, Upmat = U_matrix_kanamori(n_orb=n_orb, U_int=U, J_hund=J)
|
|
# Construct density-density Hamiltonian
|
|
h_int = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat)
|
|
|
|
## SLATER HAMILTONIAN
|
|
## Define interaction paramters, DC and Hamiltonian
|
|
#U = 9.6
|
|
#J = 0.8
|
|
#dc_type = 0 # DC type: 0 FLL, 1 Held, 2 AMF
|
|
## Construct Slater U matrix
|
|
#U_sph = U_matrix(l=2, U_int=U, J_hund=J)
|
|
#U_cubic = transform_U_matrix(U_sph, spherical_to_cubic(l=2, convention='wien2k'))
|
|
#Umat = t2g_submatrix(U_cubic, convention='wien2k')
|
|
## Construct Slater Hamiltonian
|
|
#h_int = h_int_slater(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U_matrix=Umat)
|
|
|
|
# Use GF structure determined by DFT blocks
|
|
gf_struct = [(block, indices) for block, indices in SK.gf_struct_solver[0].items()]
|
|
|
|
# Construct Solver
|
|
S = Solver(beta=beta, gf_struct=gf_struct)
|
|
|
|
if previous_present:
|
|
chemical_potential = 0
|
|
dc_imp = 0
|
|
dc_energ = 0
|
|
if mpi.is_master_node():
|
|
with HDFArchive(dft_filename+'.h5','r') as ar:
|
|
S.Sigma_iw << ar['dmft_output']['Sigma_iw']
|
|
chemical_potential,dc_imp,dc_energ = SK.load(['chemical_potential','dc_imp','dc_energ'])
|
|
S.Sigma_iw << mpi.bcast(S.Sigma_iw)
|
|
chemical_potential = mpi.bcast(chemical_potential)
|
|
dc_imp = mpi.bcast(dc_imp)
|
|
dc_energ = mpi.bcast(dc_energ)
|
|
SK.set_mu(chemical_potential)
|
|
SK.set_dc(dc_imp,dc_energ)
|
|
|
|
for iteration_number in range(1,loops+1):
|
|
if mpi.is_master_node(): print("Iteration = ", iteration_number)
|
|
|
|
SK.symm_deg_gf(S.Sigma_iw,ish=0) # symmetrise Sigma
|
|
SK.set_Sigma([ S.Sigma_iw ]) # set Sigma into the SumK class
|
|
chemical_potential = SK.calc_mu( precision = prec_mu ) # find the chemical potential for given density
|
|
S.G_iw << SK.extract_G_loc()[0] # calc the local Green function
|
|
mpi.report("Total charge of Gloc : %.6f"%S.G_iw.total_density())
|
|
|
|
# Init the DC term and the real part of Sigma, if no previous runs found:
|
|
if (iteration_number==1 and previous_present==False):
|
|
dm = S.G_iw.density()
|
|
SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
|
|
S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
|
|
|
|
# Calculate new G0_iw to input into the solver:
|
|
S.G0_iw << inverse(S.Sigma_iw + inverse(S.G_iw))
|
|
|
|
# Solve the impurity problem:
|
|
S.solve(h_int=h_int, **p)
|
|
|
|
# Solved. Now do post-processing:
|
|
mpi.report("Total charge of impurity problem : %.6f"%S.G_iw.total_density())
|
|
|
|
# Now mix Sigma and G with factor sigma_mix, if wanted:
|
|
if (iteration_number>1 or previous_present):
|
|
if mpi.is_master_node():
|
|
with HDFArchive(dft_filename+'.h5','r') as ar:
|
|
mpi.report("Mixing Sigma and G with factor %s"%sigma_mix)
|
|
S.Sigma_iw << sigma_mix * S.Sigma_iw + (1.0-sigma_mix) * ar['dmft_output']['Sigma_iw']
|
|
S.G_iw << sigma_mix * S.G_iw + (1.0-sigma_mix) * ar['dmft_output']['G_iw']
|
|
|
|
S.G_iw << mpi.bcast(S.G_iw)
|
|
S.Sigma_iw << mpi.bcast(S.Sigma_iw)
|
|
|
|
# Write the final Sigma and G to the hdf5 archive:
|
|
if mpi.is_master_node():
|
|
with HDFArchive(dft_filename+'.h5','a') as ar:
|
|
ar['dmft_output']['iterations'] = iteration_number + previous_runs
|
|
ar['dmft_output']['G_tau'] = S.G_tau
|
|
ar['dmft_output']['G_iw'] = S.G_iw
|
|
ar['dmft_output']['Sigma_iw'] = S.Sigma_iw
|
|
ar['dmft_output']['G0-%s'%(iteration_number)] = S.G0_iw
|
|
ar['dmft_output']['G-%s'%(iteration_number)] = S.G_iw
|
|
ar['dmft_output']['Sigma-%s'%(iteration_number)] = S.Sigma_iw
|
|
|
|
# Set the new double counting:
|
|
dm = S.G_iw.density() # compute the density matrix of the impurity problem
|
|
SK.calc_dc(dm, U_interact = U, J_hund = J, orb = 0, use_dc_formula = dc_type)
|
|
|
|
# Save stuff into the user_data group of hdf5 archive in case of rerun:
|
|
SK.save(['chemical_potential','dc_imp','dc_energ'])
|
|
|
|
#For FCSC calculations
|
|
#output the density matrix for Elk interface
|
|
dN, d = SK.calc_density_correction(dm_type='elk')
|
|
#correlation energy via the Migdal formula
|
|
correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density()
|
|
#subtract the double counting energy
|
|
correnerg -= SK.dc_energ[0]
|
|
#convert to Hartree
|
|
correnerg = correnerg/SK.energy_unit
|
|
#save the correction to energy
|
|
if (mpi.is_master_node()):
|
|
f=open('DMATDMFT.OUT','a')
|
|
f.write("%.16f\n"%correnerg)
|
|
f.close()
|
|
|