U on oxygen in SrVO3

[Tracy-Lee]

Dear developers

I'm a Wien2k user. I would like to see the effect of an U on the oxigen in SrVO3. I used the input I have attached below. The calculation runs fine but I need to know whether my input gives what I need or not.

Thanks in advance

---

from pytriqs.Wien2k.SumK_LDA import from pytriqs.Wien2k.SumK_LDA_Wien2k_input import from pytriqs.Wien2k.Solver_MultiBand import from pytriqs.Base.GF_Local import from pytriqs.Base.Archive.HDF_Archive import *

LDAFilename='srvo3' U1 = 2.7 U2 = 0.5 J1 = 0.65 J2 = 0.05 Beta = 40 Loops = 10 # Number of DMFT sc-loops Mix = 1.0 # Mixing factor of Sigma after solution of the AIM DeltaMix = 1.0 # Mixing factor of Delta as input for the AIM DC_type = 1 # DC type: 0 FLL, 1 Held, 2 AMF useBlocs = True # use bloc structure from LDA input useMatrix = False # True: Slater parameters, False: Kanamori parameters U+2J, U, U-J use_spinflip = False # use the full rotational invariant interaction? prec_mu = 0.0001 QMCcycles = 20000
Length_cycle = 200 Warming_iterations = 2000

Converter = SumK_LDA_Wien2k_input(Filename=LDAFilename,repacking=False) Converter.convert_DMFT_input() MPI.barrier()

previous_runs = 0 previous_present = False if MPI.IS_MASTER_NODE(): ar = HDF_Archive(LDAFilename+'.h5','a') if 'iterations' in ar: previous_present = True previous_runs = ar['iterations'] del ar previous_runs = MPI.bcast(previous_runs) previous_present = MPI.bcast(previous_present)

if previous runs are present, no need for recalculating the bloc structure:

calc_blocs = useBlocs and (not previous_present)

SK=SumK_LDA(HDFfile=LDAFilename+'.h5',UseLDABlocs=calc_blocs)

Norb1 = SK.corr_shells[0][3] Norb2 = SK.corr_shells[1][3] l1 = SK.corr_shells[0][2] l2 = SK.corr_shells[1][2] print Norb1, Norb2 print l1, l2 deg_orbs1=SK.deg_shells[0] deg_orbs2=SK.deg_shells[1] print deg_orbs1 print deg_orbs2

print SK.GFStruct_Solver[0] print SK.GFStruct_Solver[1]

print SK.map[0] print SK.map[1]

S1=Solver_MultiBand(Beta=Beta,U_interact=U1,J_Hund=J1,Norb=Norb1,useMatrix=useMatrix, T=SK.T[0] ,GFStruct=SK.GFStruct_Solver[0],map=SK.map[1], l=l1, deg_orbs=SK.deg_shells[0], use_spinflip=use_spinflip)

S2=Solver_MultiBand(Beta=Beta,U_interact=U2,J_Hund=J2,Norb=Norb2,useMatrix=useMatrix, T=SK.T[0] ,GFStruct=SK.GFStruct_Solver[1],map=SK.map[1], l=l2, deg_orbs=SK.deg_shells[1], use_spinflip=use_spinflip)

S1.N_Cycles = QMCcycles S1.Length_Cycle = Length_cycle S1.N_Warmup_Cycles = Warming_iterations

S2.N_Cycles = QMCcycles S2.Length_Cycle = Length_cycle S2.N_Warmup_Cycles = Warming_iterations

if (previous_present): if (MPI.IS_MASTER_NODE()): ar = HDF_Archive(LDAFilename+'.h5','a') S1.Sigma <<= ar['SigmaF1'] S2.Sigma <<= ar['SigmaF2'] del ar S1.Sigma = MPI.bcast(S1.Sigma) S2.Sigma = MPI.bcast(S2.Sigma)

for IterationNumber in range(1,Loops+1) :

  SK.symm_deg_GF(S1.Sigma,orb=0)                                 # symmetrise Sigma
  SK.symm_deg_GF(S2.Sigma,orb=1)                                 # symmetrise Sigma
  SK.put_Sigma(Sigmaimp = [ S1.Sigma,S2.Sigma ])                 # put Sigma into the SumK class:

  Chemical_potential = SK.find_mu( precision = prec_mu )         # find the chemical potential
  S1.G <<= SK.extract_Gloc()[0]                                 # calculation of the local Green function
  S2.G <<= SK.extract_Gloc()[1]                                 # calculation of the local Green function
  MPI.report("Total charge of G1loc : %.6f"%S1.G.total_density())
  MPI.report("Total charge of G2loc : %.6f"%S2.G.total_density())

  if ((IterationNumber==1)and(previous_present==False)):
      # Init the DC term and the real part of Sigma, if no previous run was found:
      dm1 = S1.G.density()
      SK.SetDoubleCounting( dm1, U_interact = U1, J_Hund = J1, orb = 0, useDCformula = DC_type)
      S1.Sigma <<= GF_Initializers.Const(SK.dc_imp[0]['up'][0,0])

      dm2 = S2.G.density()
      SK.SetDoubleCounting( dm2, U_interact = U2, J_Hund = J2, orb = 1, useDCformula = DC_type)
      S2.Sigma <<= GF_Initializers.Const(SK.dc_imp[1]['up'][0,0])

  # now calculate new G0:
  if (MPI.IS_MASTER_NODE()):
      # We can do a mixing of Delta in order to stabilize the DMFT iterations:
      S1.G0 <<= S1.Sigma + inverse(S1.G)
      S2.G0 <<= S2.Sigma + inverse(S2.G)
      ar = HDF_Archive(LDAFilename+'.h5','a')
      if ((IterationNumber>1) or (previous_present)):
          MPI.report("Mixing input Delta with factor %s"%DeltaMix)
          Delta1 = DeltaMix * S1.G0.Delta() + (1.0-DeltaMix) * ar['DeltaF']
          Delta2 = DeltaMix * S2.G0.Delta() + (1.0-DeltaMix) * ar['DeltaF']
          S1.G0 <<= S1.G0 + S1.G0.Delta() - Delta1
          S2.G0 <<= S2.G0 + S2.G0.Delta() - Delta2

      ar['DeltaF'] = S1.G0.Delta() + S1.G0.Delta()
      S1.G0 <<= inverse(S1.G0)
      S2.G0 <<= inverse(S2.G0)
      del ar
      #### S.G0 <<= inverse(S.Sigma+inverse(S.G))

  S1.G0 = MPI.bcast(S1.G0)
  S2.G0 = MPI.bcast(S2.G0)

  # Solve the impurity problem:
  S1.Solve()
  S2.Solve()

  # solution done, do the post-processing:
  MPI.report("Total charge of impurity1 problem : %.6f"%S1.G.total_density())
  MPI.report("Total charge of impurity2 problem : %.6f"%S2.G.total_density())

  # Now mix Sigma and G with factor Mix, if wanted:
  if ((IterationNumber>1) or (previous_present)):
      if (MPI.IS_MASTER_NODE()):
          ar = HDF_Archive(LDAFilename+'.h5','a')
          MPI.report("Mixing Sigma and G with factor %s"%Mix)
          S1.Sigma <<= Mix * S1.Sigma + (1.0-Mix) * ar['SigmaF1']
          S2.Sigma <<= Mix * S2.Sigma + (1.0-Mix) * ar['SigmaF2']
          S1.G <<= Mix * S1.G + (1.0-Mix) * ar['GF1']
          S2.G <<= Mix * S2.G + (1.0-Mix) * ar['GF2']
          del ar
      S1.G = MPI.bcast(S1.G)
      S2.G = MPI.bcast(S2.G)
      S1.Sigma = MPI.bcast(S1.Sigma)
      S2.Sigma = MPI.bcast(S2.Sigma)

  # Write the final Sigma and G to the hdf5 archive:
  if (MPI.IS_MASTER_NODE()):
      ar = HDF_Archive(LDAFilename+'.h5','a')
      ar['iterations'] = previous_runs + IterationNumber
      ar['SigmaF1'] = S1.Sigma
      ar['SigmaF2'] = S2.Sigma
      ar['GF1'] = S1.G
      ar['GF2'] = S1.G
      del ar

  # Now set new double counting:
  dm1 = S1.G.density()
  dm2 = S1.G.density()
  SK.SetDoubleCounting( dm1, U_interact = U1, J_Hund = J1, orb = 0, useDCformula = DC_type)
  SK.SetDoubleCounting( dm2, U_interact = U2, J_Hund = J2, orb = 1, useDCformula = DC_type)

  #Save stuff:
  SK.save()

import os LDAFilename = os.getcwd().rpartition('/')[2]

SK.symm_deg_GF(S1.Sigma,orb=0) SK.symm_deg_GF(S2.Sigma,orb=1) S1.G <<= inverse(S1.G0) - S1.Sigma S2.G <<= inverse(S2.G0) - S2.Sigma S1.G.invert() S2.G.invert()

find exact chemical potential

SK.put_Sigma(Sigmaimp = [ S1.Sigma,S2.Sigma ]) SK.Chemical_potential = SK.find_mu( precision = 0.000001 ) dN,d = SK.calc_DensityCorrection(Filename = LDAFilename+'.qdmft') SK.save()

correnerg1 = 0.5 * (S1.G * S1.Sigma).total_density() correnerg2 = 3 * 0.5 * (S2.G * S2.Sigma).total_density()

correnerg1 -= SK.DCenerg[0] correnerg2 -= SK.DCenerg[1]

if (MPI.IS_MASTER_NODE()): f=open(LDAFilename+'.qdmft','a') f.write("%.16f\n"%correnerg1) f.write("%.16f\n"%correnerg2) f.close()

[aichhorn]

Before looking at this input, the first thing you have to do is to construct Wannier functions also for oxygen with dmftproj. Following the instructions in the manual this should be no problem, and I guess you did this already. Otherwise the script looks fine, except the last lines. You have to write the sum of the correlation energies at the end of the .qdmft file, not two separate lines. And the factor 3 should also be in front of the DC energy. But it you are not preliminarily interested in total energies, you can forget about these details. Another bug: There is S1 instead of S2 two times when storing the Greens functions, and when setting the double counting. There might be more subtle details like this (I did not try the script myself)...

[parcollet]

Inactive for 3 months. Closing.