Hubbard-I solver on SrVO3

[Tracy-Lee]

Dear all

I would like to use the Hubbard-I solver on SrVO3. I tried using Ce-gamma.py script. After the calculation of the chemical potential pytriqs complains in the assignment : S.G <<= SK.extract_Gloc()[0]

Can someone help me to build a script suitable for the calculation I want to make?

Thank you in advance.

[leopo]

Because SrVO3 is a metal, the Hubbard-I solver is not generally suitable for it.

Could you show the error message you got?

[Tracy-Lee]

The error message I got is:

Total Density = 0.999999;Chemical_Potential = 0.075647 Traceback (most recent call last): File "srvo3_04.py", line 76, in S.G <<= SK.extract_Gloc()[0] File "/opt/TRIQS/lib/python2.7/dist-packages/pytriqs/Base/GF_Local/GF.py", line 248, in ilshift for (i,g) in self : g.copyFrom(A[i]) RuntimeError: Triqs runtime error at /home/albert/TRIQS/TRIQS-git/triqs/gf_local/GF_Bloc_Base.hpp : 106

Trace is :

/opt/TRIQS/lib/libtriqs.so void check_have_same_structurestd::complex<double, std::complex >(GF_Bloc_Basestd::complex const&, GF_Bloc_Basestd::complex const&, bool, bool) 0x27e [0x2b20dba0ce4e] /opt/TRIQS/lib/libtriqs.so GF_Bloc_Basestd::complex::operator=(GF_Bloc_Basestd::complex const&) 0x20 [0x2b20dba15c10] /opt/TRIQS/lib/python2.7/dist-packages/pytriqs/Base/GF_Local/_pytriqs_GF.so boost::python::objects::caller_py_function_impl<boost::python::detail::caller<void (GF_BlocBase<std::complex >::)(GF_Bloc_Basestd::complex const&), boost::python::default_call_policies, boost::mpl::vector3<void, GF_Bloc_Base<std::complex >&, GF_Bloc_Basestd::complex const&> > >::operator()(object, object) 0xcf [0x2b20daa050df] /opt/TRIQS/lib/libboost_for_triqs.so boost::python::objects::function::call(object, object) const 0x1af [0x2b20db24c96f] /opt/TRIQS/lib/libboost_for_triqs.so 0xf1bc8 [0x2b20db24cbc8] /opt/TRIQS/lib/libboost_for_triqs.so boost::python::detail::exception_handler::operator()(boost::function0 const&) const 0x43 [0x2b20db237d23] /opt/TRIQS/lib/python2.7/dist-packages/pytriqs/boost/mpi.so boost::detail::function::function_obj_invoker2<boost::_bi::bind_t<bool, boost::python::detail::translate_exception<boost::mpi::python::object_without_skeleton, boost::mpi::python::translate_exception >, boost::_bi::list3boost::arg<1, boost::arg<2>, boost::_bi::valueboost::mpi::python::translate_exception > >, bool, boost::python::detail::exception_handler const&, boost::function0 const&>::invoke(boost::detail::function::function_buffer&, boost::python::detail::exception_handler const&, boost::function0 const&) 0x1c [0x2b20e0377f2c] /opt/TRIQS/lib/libboost_for_triqs.so boost::python::detail::exception_handler::operator()(boost::function0 const&) const 0x26 [0x2b20db237d06] /opt/TRIQS/lib/python2.7/dist-packages/pytriqs/boost/mpi.so boost::detail::function::function_obj_invoker2<boost::_bi::bind_t<bool, boost::python::detail::translate_exception<boost::mpi::exception, boost::mpi::python::translate_exception >, boost::_bi::list3boost::arg<1, boost::arg<2>, boost::_bi::valueboost::mpi::python::translate_exception > >, bool, boost::python::detail::exception_handler const&, boost::function0 const&>::invoke(boost::detail::function::function_buffer&, boost::python::detail::exception_handler const&, boost::function0 const&) 0x1c [0x2b20e036ef7c] /opt/TRIQS/lib/libboost_for_triqs.so boost::python::detail::exception_handler::operator()(boost::function0 const&) const 0x26 [0x2b20db237d06] /opt/TRIQS/lib/python2.7/dist-packages/pytriqs/Base/GF_Local/_pytriqs_GF.so boost::detail::function::function_obj_invoker2<boost::_bi::bind_t<bool, boost::python::detail::translateexception<std::string ()(std::string const&)>, boost::_bi::list3boost::arg<1, boost::arg<2>, boost::_bi::value<void (*)(std::string const&)> > >, bool, boost::python::detail::exception_handler const&, boost::function0 const&>::invoke(boost::detail::function::function_buffer&, boost::python::detail::exception_handler const&, boost::function0 const&) 0x13 [0x2b20daa01a23] /opt/TRIQS/lib/libboost_for_triqs.so boost::python::handle_exception_impl(boost::function0) 0x29 [0x2b20db237b09] /opt/TRIQS/lib/libboost_for_triqs.so 0xefe64 [0x2b20db24ae64] /usr/bin/python PyObject_Call 0x36 [0x4e9f36] /usr/bin/python PyEval_EvalFrameEx 0x86a [0x49846a] /usr/bin/python PyEval_EvalCodeEx 0x1a0 [0x49f1c0] /usr/bin/python [0x4a8960] /usr/bin/python PyObject_Call 0x36 [0x4e9f36] /usr/bin/python [0x4ec11a] /usr/bin/python PyObject_Call 0x36 [0x4e9f36] /usr/bin/python [0x479f97] /usr/bin/python PyEval_EvalFrameEx 0x644f [0x49e04f] /usr/bin/python PyEval_EvalCodeEx 0x1a0 [0x49f1c0] /usr/bin/python PyRun_FileExFlags 0xe1 [0x4a9081] /usr/bin/python PyRun_SimpleFileExFlags 0x1d1 [0x4a9311] /usr/bin/python Py_Main 0x55d [0x4aa8bd] /lib/x86_64-linux-gnu/libc.so.6 __libc_start_main 0xed [0x2b20d6d9376d] /usr/bin/python [0x41b9b1]

The first dimension of the two Green functions is not the same

[aichhorn]

There might be an issue with the Block structure of the Greens functions. The Hubbard I Solver does not compute an optimal structure of G, but assumes just G['up'] and G['down'], with indices 0,1,2 for the t2g orbitals of SrVO3. When initialising the SumK_LDA, did you use the UseLDABlocs = True option? Then the thing could crash. Set UseLDABlocs = False instead. But as Leonid said, using Hubbard I for SrVO3 is not at all recommended, results will be unphysical.

[Tracy-Lee]

This is the script I used:

from pytriqs.Wien2k.SumK_LDA import from pytriqs.Wien2k.SumK_LDA_Wien2k_input import from pytriqs.Solvers.HubbardI.Solver_HubbardI import Solver_HubbardI

LDAFilename = 'srvo3_04' Beta = 40 Uint = 2.70 JHund = 0.70 Loops = 3 # Number of DMFT sc-loops Mix = 0.7 # 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)

HDFfilename = LDAFilename+'.h5'

Convert DMFT input:

Can be commented after the first run

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

check if there are previous runs:

previous_runs = 0 previous_present = False

if MPI.IS_MASTER_NODE(): ar = HDF_Archive(HDFfilename,'a') if 'iterations' in ar: previous_present = True previous_runs = ar['iterations'] else: previous_runs = 0 previous_present = False del ar

MPI.barrier() previous_runs = MPI.bcast(previous_runs) previous_present = MPI.bcast(previous_present)

Init the SumK class

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

Norb = SK.corr_shells[0][3] l = SK.corr_shells[0][2]

Init the Solver:

S = Solver_HubbardI(Beta = Beta, Uint = Uint, JHund = JHund, l = l, Verbosity=2) S.Nmoments=10

if (previous_present):

load previous data:

MPI.report("Using stored data for initialisation")
if (MPI.IS_MASTER_NODE()):
    ar = HDF_Archive(HDFfilename,'a')
    S.Sigma <<= ar['SigmaF']
    del ar
S.Sigma = MPI.bcast(S.Sigma)
SK.load()

DMFT loop:

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

    itn = Iteration_Number + previous_runs

    # put Sigma into the SumK class:
    SK.put_Sigma(Sigmaimp = [ S.Sigma ])

    # Compute the SumK, possibly fixing mu by dichotomy
    if SK.Density_Required and (Iteration_Number > 0):
        Chemical_potential = SK.find_mu( precision = 0.000001 )
    else:
        MPI.report("No adjustment of chemical potential\nTotal density  = %.3f"%SK.total_density(mu=Chemical_potential))

    # Density:
    S.G <<= SK.extract_Gloc()[0]
    MPI.report("Total charge of Gloc : %.6f"%S.G.total_density())
    dm = S.G.density()

    if ((Iteration_Number==1)and(previous_present==False)):
 SK.SetDoubleCounting( dm, U_interact = Uint, J_Hund = JHund, orb = 0, useDCformula = DC_type)

    # set atomic levels:
    eal = SK.eff_atomic_levels()[0]
    S.set_atomic_levels( eal = eal )

    # update hdf5
    if (MPI.IS_MASTER_NODE()):
        ar = HDF_Archive(HDFfilename,'a')
        ar['Chemical_Potential%s'%itn] = Chemical_potential
        del ar

    # solve it:
    S.Solve()

    if (MPI.IS_MASTER_NODE()):
        ar = HDF_Archive(HDFfilename)
        ar['iterations'] = itn

    # Now mix Sigma and G:
    if ((itn>1)or(previous_present)):
        if (MPI.IS_MASTER_NODE()):
            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)

    if (MPI.IS_MASTER_NODE()):
        ar['SigmaF'] = S.Sigma
        ar['GF'] = S.G

    # after the Solver has finished, set new double counting: 
    dm = S.G.density()
    SK.SetDoubleCounting( dm, U_interact = Uint, J_Hund = JHund, orb = 0, useDCformula = DC_type )
    # correlation energy calculations:
    correnerg = 0.5 * (S.G * S.Sigma).total_density()
    MPI.report("Corr. energy = %s"%correnerg)
    if (MPI.IS_MASTER_NODE()):
        ar['correnerg%s'%itn] = correnerg
        ar['DCenerg%s'%itn] = SK.DCenerg
        del ar

    #Save stuff:
    SK.save()
    if (MPI.IS_MASTER_NODE()):
        print 'DC after solver: ',SK.dc_imp[SK.invshellmap[0]]

    # do some analysis:
    MPI.report("Orbital densities of impurity Green function:")
    dm1 = S.G.density()
    for s in dm1:
        MPI.report("Block %s: "%s)
        for ii in range(len(dm1[s])):
            str = ''
            for jj in range(len(dm1[s])):
                if (dm1[s][ii,jj].real>0):
                    str += "   %.4f"%(dm1[s][ii,jj].real)
                else:
                    str += "  %.4f"%(dm1[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.find_mu( precision = 0.000001 ) dN,d = SK.calc_DensityCorrection(Filename = LDAFilename+'.qdmft')

MPI.report("Trace of Density Matrix: %s"%d)

correlation energy:

if (MPI.IS_MASTER_NODE()): ar = HDF_Archive(HDFfilename) itn = ar['iterations'] correnerg = ar['correnerg%s'%itn] DCenerg = ar['DCenerg%s'%itn] del ar correnerg -= DCenerg[0] f=open(LDAFilename+'.qdmft','a') f.write("%.16f\n"%correnerg) f.close()

[leopo]

This script should be fine. However, you should change SrVO3.indmftpr there you will need to include the full V 3d shell as correlated (instead of including only t2g), hence, you substitute the old description of V

cubic 1 1 2 0 ! l included for each sort 0 0 2 0 ! l included for each sort 01 ! t2g Wannier functions will be constructed 0 ! no spin-orbit

with the following one

cubic 1 1 2 0 ! l included for each sort 0 0 0 0 ! l included for each sort 0 ! no spin-orbit

and rerun dmftproj

[Tracy-Lee]

Dear Leonid

Thank you for your help. If I use the full d shell pytriqs doesn’t give the error. Can you please explain why it is not possible to use only the t2g manifold as for the QMC-solver?

[leopo]

Generally, treating only one particular representation (e.g. t2g) instead of the full shell is not particularly useful in the case of Hub-I. In contrast to CT-QMC, the solver is very fast, and one can easily apply it to the full d-shell with almost no additional computational cost. For this reason this has not been implemented in Hub-I, though it is not difficult to do.

[Tracy-Lee]

Thank you Leonid.