Added an example of Ce HF calculation

python/converters/plovasp/examples/ce/.gitignore

@@ -0,0 +1,9 @@
+*
+!.gitignore
+!INCAR
+!POSCAR
+!POTCAR
+!KPOINTS
+!plo.cfg
+!*.sh
+!*.py

python/converters/plovasp/examples/ce/INCAR

@@ -0,0 +1,49 @@
+ADDGRID = TRUE
+AMIX = 0.1
+AMIN = 0.01
+BMIX = 0.0001
+#IMIX = 1
+#MIXFIRST = True
+ALGO = NORMAL
+EDIFF = 0.00001
+EDIFFG = -0.002
+EMAX = 13.0
+EMIN = -2.0
+ENCUT = 600
+ENAUG = 1100
+IBRION = -1
+ICHARG = 5
+ISIF = 4
+ISMEAR = 2
+ISPIN = 2
+ISYM = -1
+LASPH = TRUE
+LCHARG = TRUE
+LDAU = True
+LDAUJ = 0
+LDAUL = 3
+LDAUTYPE = 2
+LDAUU = 0.0
+LDAUPRINT = 2
+LMAXMIX = 6
+LORBIT = 13
+LREAL = False
+ROPT = 1e-3 1e-3 1e-3
+LWAVE = TRUE
+MAGMOM = 1.1
+#NBANDS = 16
+NEDOS = 501
+NELM = 100
+NELMDL = -7
+NELMIN = 1
+NPAR = 4
+NSIM = 1
+NSW = 0
+POTIM = 0.4
+PREC = Accurate
+SIGMA = 0.1
+SMASS = 0.5
+SYMPREC = 1.0e-6
+
+# Projectors
+LOCPROJ = 1 : f : Pr 1

python/converters/plovasp/examples/ce/KPOINTS

@@ -0,0 +1,4 @@
+Automatic kpoint scheme
+0
+Gamma
+4 4 4

python/converters/plovasp/examples/ce/POSCAR

@@ -0,0 +1,9 @@
+gamma-Ce
+   1.0
+     2.5790000000000000    2.5790000000000000    0.0000000000000000
+     0.0000000000000000    2.5790000000000000    2.5790000000000000
+     2.5790000000000000    0.0000000000000000    2.5790000000000000
+   Ce
+    1
+Direct
+   0.0  0.0  0.0

python/converters/plovasp/examples/ce/U_matrix.py

@@ -0,0 +1,252 @@
+from math import sqrt
+from scipy.misc import factorial as fact
+from itertools import product
+import numpy as np
+
+# The interaction matrix in desired basis
+# U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4)
+def U_matrix(l, radial_integrals=None, U_int=None, J_hund=None, basis="spherical", T=None):
+    """Calculate U matrix being given either radial_integrals or U_int and J_hund.
+       l = angular momentum of shell being treated (l=2 for d shell, l=3 for f shell)
+       radial_integrals = [F0,F2,F4,..] (default None)
+       U_int, J_hund = values to use to compute radial_integrals (default None),
+       basis = "spherical", "cubic", or "other",
+       T = transformation matrix from spherical to desired basis, if basis='other' (default None)"""
+
+    # Check all necessary information is present and consistent
+    if radial_integrals is None and (U_int is None and J_hund is None):
+        raise ValueError("U_matrix: provide either the radial_integrals or U_int and J_hund.")
+    if radial_integrals is None and (U_int is not None and J_hund is not None):
+        radial_integrals = U_J_to_radial_integrals(l, U_int, J_hund)
+    if radial_integrals is not None and (U_int is not None and J_hund is not None):
+        if len(radial_integrals)-1 != l:
+            raise ValueError("U_matrix: inconsistency in l and number of radial_integrals provided.")
+        if (radial_integrals - U_J_to_radial_integrals(l, U_int, J_hund)).any() != 0.0:
+            print "Warning: U_matrix: radial_integrals provided do not match U_int and J_hund. Using radial_integrals to calculate spherical U_matrix."
+
+    # Full interaction matrix
+    # Basis of spherical harmonics Y_{-2}, Y_{-1}, Y_{0}, Y_{1}, Y_{2}
+    # U^{spherical}_{m1 m2 m3 m4} = \sum_{k=0}^{2l} F_k angular_matrix_element(l, k, m1, m2, m3, m4)
+    U_matrix = np.zeros((2*l+1,2*l+1,2*l+1,2*l+1),dtype=float)
+
+    m_range = range(-l,l+1)
+    for n, F in enumerate(radial_integrals):
+        k = 2*n
+        for m1, m2, m3, m4 in product(m_range,m_range,m_range,m_range):
+            U_matrix[m1+l,m2+l,m3+l,m4+l] += F * angular_matrix_element(l,k,m1,m2,m3,m4)
+
+    # Transform from spherical basis if needed
+    if basis == "cubic": T = spherical_to_cubic(l)
+    if basis == "other" and T is None:
+        raise ValueError("U_matrix: provide T for other bases.")
+    if T is not None: U_matrix = transform_U_matrix(U_matrix, T)
+
+    return U_matrix
+
+# Convert full 4-index U matrix to 2-index density-density form
+def reduce_4index_to_2index(U_4index):
+    """Reduces the four-index matrix to two-index matrices."""
+
+    size = len(U_4index) # 2l+1
+    U  = np.zeros((size,size),dtype=float)      # matrix for same spin
+    Uprime = np.zeros((size,size),dtype=float)  # matrix for opposite spin
+
+    m_range = range(size)
+    for m,mp in product(m_range,m_range):
+        U[m,mp]  = U_4index[m,mp,m,mp].real - U_4index[m,mp,mp,m].real
+        Uprime[m,mp] = U_4index[m,mp,m,mp].real
+
+    return U, Uprime
+
+# Construct the 2-index matrices for the density-density form
+def U_matrix_kanamori(n_orb, U_int, J_hund):
+    """Calculate the Kanamori U and Uprime matrices."""
+
+    U  = np.zeros((n_orb,n_orb),dtype=float)      # matrix for same spin
+    Uprime = np.zeros((n_orb,n_orb),dtype=float)  # matrix for opposite spin
+
+    m_range = range(n_orb)
+    for m,mp in product(m_range,m_range):
+        if m == mp:
+            Uprime[m,mp] = U_int
+        else:
+            U[m,mp]  = U_int - 3.0*J_hund
+            Uprime[m,mp] = U_int - 2.0*J_hund
+
+    return U, Uprime
+
+# Get t2g or eg components
+def t2g_submatrix(U, convention=''):
+    """Return only the t2g part of the full d-manifold two- or four-index U matrix."""
+    if convention == 'wien2k': 
+        return subarray(U, len(U.shape)*[(2,3,4)])
+    else:
+        return subarray(U, len(U.shape)*[(0,1,3)])
+
+def eg_submatrix(U, convention=''):
+    """Return only the eg part of the full d-manifold two- or four-index U matrix."""
+    if convention == 'wien2k': 
+        return subarray(U, len(U.shape)*[(0,1)])
+    else:
+        return subarray(U, len(U.shape)*[(2,4)])
+
+# Transform the interaction matrix into another basis
+def transform_U_matrix(U_matrix, T):
+    """Transform the interaction matrix into another basis by applying matrix T."""
+    return np.einsum("ij,kl,jlmo,mn,op",np.conj(T),np.conj(T),U_matrix,np.transpose(T),np.transpose(T))
+
+# Rotation matrices: complex harmonics to cubic harmonics
+# Complex harmonics basis: ..., Y_{-2}, Y_{-1}, Y_{0}, Y_{1}, Y_{2}, ...
+def spherical_to_cubic(l, convention=''):
+    """Returns the spherical harmonics to cubic harmonics rotation matrix."""
+    size = 2*l+1
+    T = np.zeros((size,size),dtype=complex)
+    if convention != 'wien2k' and l == 2: 
+        raise ValueError("spherical_to_cubic: standard convention not implemented for l=1")
+    if l == 0:
+        cubic_names = ("s")
+    elif l == 1:
+        cubic_names = ("x","y","z")
+        T[0,0] = 1.0/sqrt(2);   T[0,2] = -1.0/sqrt(2)
+        T[1,0] = 1j/sqrt(2);    T[1,2] = 1j/sqrt(2)
+        T[2,1] = 1.0
+    elif l == 2:
+        if convention == 'wien2k':
+            cubic_names = ("z^2","x^2-y^2","xy","yz","xz")
+            T[0,2] = 1.0
+            T[1,0] = 1.0/sqrt(2);   T[1,4] = 1.0/sqrt(2)
+            T[2,0] = -1j/sqrt(2);   T[2,4] = 1j/sqrt(2)
+            T[3,1] = 1j/sqrt(2);    T[3,3] = -1j/sqrt(2)
+            T[4,1] = 1.0/sqrt(2);   T[4,3] = 1.0/sqrt(2)
+        else:
+            cubic_names = ("xy","yz","z^2","xz","x^2-y^2")
+            T[0,0] = 1j/sqrt(2);    T[0,4] = -1j/sqrt(2)
+            T[1,1] = 1j/sqrt(2);    T[1,3] = 1j/sqrt(2)
+            T[2,2] = 1.0
+            T[3,1] = 1.0/sqrt(2);   T[3,3] = -1.0/sqrt(2)
+            T[4,0] = 1.0/sqrt(2);   T[4,4] = 1.0/sqrt(2)
+    elif l == 3:
+        if convention == 'wien2k':
+            cubic_names = ("x(x^2-3y^2)","z(x^2-y^2)","xz^2","z^3","yz^2","xyz","y(3x^2-y^2)")
+            T[0,0] = 1.0/sqrt(2);    T[0,6] = -1.0/sqrt(2)
+            T[1,1] = 1.0/sqrt(2);    T[1,5] = 1.0/sqrt(2)
+            T[2,2] = 1.0/sqrt(2);    T[2,4] = -1.0/sqrt(2)
+            T[3,5] = 1.0
+            T[4,2] = 1j/sqrt(2);   T[4,4] = 1j/sqrt(2)
+            T[5,1] = 1j/sqrt(2);   T[5,5] = -1j/sqrt(2)
+            T[6,0] = 1j/sqrt(2);   T[6,6] = 1j/sqrt(2)
+        else:
+            cubic_names = ("y(3x^2-y^2)","xyz","yz^2","z^3","xz^2","z(x^2-y^2)","x(x^2-3y^2)")
+# Verified with the matrix generated in VASP
+            T[0,0] = 1j/sqrt(2);   T[0,6] = 1j/sqrt(2)
+            T[1,1] = 1j/sqrt(2);   T[1,5] = -1j/sqrt(2)
+            T[2,2] = 1j/sqrt(2);   T[2,4] = 1j/sqrt(2)
+            T[3,3] = 1.0
+            T[4,2] = 1.0/sqrt(2);    T[4,4] = -1.0/sqrt(2)
+            T[5,1] = 1.0/sqrt(2);    T[5,5] = 1.0/sqrt(2)
+            T[6,0] = 1.0/sqrt(2);    T[6,6] = -1.0/sqrt(2)
+    else: raise ValueError("spherical_to_cubic: implemented only for l=0,1,2,3")
+
+    return T
+
+# Names of cubic harmonics
+def cubic_names(l):
+    if l == 0 or l == 's':
+        return ("s")
+    elif l == 1 or l == 'p':
+        return ("x","y","z")
+    elif l == 2 or l == 'd':
+        return ("xy","yz","z^2","xz","x^2-y^2")
+    elif l == 't2g':
+        return ("xy","yz","xz")
+    elif l == 'eg':
+        return ("z^2","x^2-y^2")
+    elif l == 3 or l == 'f':
+        return ("x(x^2-3y^2)","z(x^2-y^2)","xz^2","z^3","yz^2","xyz","y(3x^2-y^2)")
+    else: raise ValueError("cubic_names: implemented only for l=0,1,2,3")
+
+# Convert U,J -> radial integrals F_k
+def U_J_to_radial_integrals(l, U_int, J_hund):
+    """Determines the radial integrals F_k from U_int and J_hund."""
+
+    F = np.zeros((l+1),dtype=float)
+    if l == 2:
+        F[0] = U_int
+        F[1] = J_hund * 14.0 / (1.0 + 0.63)
+        F[2] = 0.630 * F[1]
+    elif l == 3:
+        F[0] = U_int
+        F[1] = 6435.0 * J_hund / (286.0 + 195.0 * 451.0 / 675.0 + 250.0 * 1001.0 / 2025.0)
+        F[2] = 451.0 * F[1] / 675.0
+        F[3] = 1001.0 * F[1] / 2025.0
+    else: raise ValueError("U_J_to_radial_integrals: implemented only for l=2,3")
+
+    return F
+
+# Convert radial integrals F_k -> U,J
+def radial_integrals_to_U_J(l, F):
+    """Determines U_int and J_hund from the radial integrals."""
+
+    if l == 2:
+        U_int = F[0]
+        J_hund = F[1] * (1.0 + 0.63) / 14.0
+    elif l == 3:
+        U_int = F[0]
+        J_hund = F[1] * (286.0 + 195.0 * 451.0 / 675.0 + 250.0 * 1001.0 / 2025.0) / 6435.0
+    else: raise ValueError("radial_integrals_to_U_J: implemented only for l=2,3")
+
+    return U_int,J_hund
+
+# Angular matrix elements of particle-particle interaction
+# (2l+1)^2 ((l 0) (k 0) (l 0))^2 \sum_{q=-k}^{k} (-1)^{m1+m2+q} ((l -m1) (k q) (l m3)) ((l -m2) (k -q) (l m4))
+def angular_matrix_element(l, k, m1, m2, m3, m4):
+    result = 0
+    for q in range(-k,k+1):
+        result += three_j_symbol((l,-m1),(k,q),(l,m3))*three_j_symbol((l,-m2),(k,-q),(l,m4))*(-1.0 if (m1+q+m2) % 2 else 1.0)
+    result *= (2*l+1)**2 * (three_j_symbol((l,0),(k,0),(l,0))**2)
+    return result
+
+# Wigner 3-j symbols
+# ((j1 m1) (j2 m2) (j3 m3))
+def three_j_symbol(jm1, jm2, jm3):
+    j1, m1 = jm1
+    j2, m2 = jm2
+    j3, m3 = jm3
+
+    if (m1+m2+m3 != 0 or
+        m1 < -j1 or m1 > j1 or
+        m2 < -j2 or m2 > j2 or
+        m3 < -j3 or m3 > j3 or
+        j3 > j1 + j2 or
+        j3 < abs(j1-j2)):
+        return .0
+
+    result = -1.0 if (j1-j2-m3) % 2 else 1.0
+    result *= sqrt(fact(j1+j2-j3)*fact(j1-j2+j3)*fact(-j1+j2+j3)/fact(j1+j2+j3+1))
+    result *= sqrt(fact(j1-m1)*fact(j1+m1)*fact(j2-m2)*fact(j2+m2)*fact(j3-m3)*fact(j3+m3))
+
+    t_min = max(j2-j3-m1,j1-j3+m2,0)
+    t_max = min(j1-m1,j2+m2,j1+j2-j3)
+
+    t_sum = 0
+    for t in range(t_min,t_max+1):
+        t_sum += (-1.0 if t % 2 else 1.0)/(fact(t)*fact(j3-j2+m1+t)*fact(j3-j1-m2+t)*fact(j1+j2-j3-t)*fact(j1-m1-t)*fact(j2+m2-t))
+
+    result *= t_sum
+    return result
+
+# Clebsch-Gordan coefficients
+# < j1 m1 j2 m2 | j3 m3 > = (-1)^{j1-j2+m3} \sqrt{2j3+1} ((j1 m1) (j2 m2) (j3 -m3))
+def clebsch_gordan(jm1, jm2, jm3):
+    norm = sqrt(2*jm3[0]+1)*(-1 if jm1[0]-jm2[0]+jm3[1] % 2 else 1)
+    return norm*three_j_symbol(jm1,jm2,(jm3[0],-jm3[1]))
+
+# Create subarray containing columns in idxlist
+# e.g. idxlist = [(0),(2,3),(0,1,2,3)] gives
+#  column 0 for 1st dim,
+#  columns 2 and 3 for 2nd dim,
+#  columns 0,1,2 and 3 for 3rd dim.
+def subarray(a,idxlist,n=None) :
+    if n is None: n = len(a.shape)-1
+    sa = a[tuple(slice(x) for x in a.shape[:n]) + (idxlist[n],)]
+    return subarray(sa,idxlist, n-1) if n > 0 else sa

python/converters/plovasp/examples/ce/hf_solver.py

@@ -0,0 +1,335 @@
+
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Ferrero, O. Parcollet
+#
+# TRIQS is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# TRIQS. If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+from types import *
+#from pytriqs.applications.dft.U_matrix import *
+from U_matrix import *
+from pytriqs.gf.local import *
+#from hubbard_I import gf_hi_fullu, sigma_atomic_fullu
+import pytriqs.utility.mpi as mpi
+from itertools import izip
+import numpy
+import copy
+
+class Solver:
+    """
+       Hartree-Fock Solver
+    """
+
+    # initialisation:
+    def __init__(self, beta, l, n_msb=1025, use_spin_orbit=False, Nmoments=5, dudarev=False):
+
+        self.name = "Hartree-Fock"
+        self.beta = beta
+        self.l = l
+        self.Nmsb = n_msb
+        self.UseSpinOrbit = use_spin_orbit
+        self.Converged = False
+        self.Nspin = 2
+        self.Nmoments=Nmoments
+
+        self.dudarev = dudarev
+
+        assert use_spin_orbit == False, "Spin-orbit is not implemented"
+
+        self.Nlm = 2*l+1
+        if (use_spin_orbit):
+            # no blocks!
+            self.gf_struct = [ ('ud', range(2*self.Nlm)) ]
+        else:
+            # up/down blocks:
+            self.gf_struct = [ ('up', range(self.Nlm)), ('down', range(self.Nlm)) ]
+
+        # construct Greens functions:
+        self.a_list = [a for a,al in self.gf_struct]
+        glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb) for a,al in self.gf_struct]
+        self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
+        self.G_Old = self.G.copy()
+        self.G0 = self.G.copy()
+        self.Sigma = self.G.copy()
+        self.Sigma_Old = self.G.copy()
+
+
+
+    def solve(self, U_int, J_hund, T=None, verbosity=0, Iteration_Number=1, Test_Convergence=0.0001):
+        """Calculation of the impurity Greens function using Hubbard-I"""
+
+        if self.Converged :
+            mpi.report("Solver %(name)s has already converged: SKIPPING"%self.__dict__)
+            return
+
+        if mpi.is_master_node():
+            self.verbosity = verbosity
+        else:
+            self.verbosity = 0
+
+        #self.Nmoments = 5
+
+        ur,ujmn,umn=self.__set_umatrix(U=U_int,J=J_hund,T=T)
+
+
+        M = [x for x in self.G.mesh]
+        self.zmsb = numpy.array([x for x in M],numpy.complex_)
+
+#        # for the tails:
+#        tailtempl={}
+#        for sig,g in self.G:
+#            tailtempl[sig] = copy.deepcopy(g.tail)
+#            for i in range(9): tailtempl[sig][i] *= 0.0
+
+#        self.__save_eal('eal.dat',Iteration_Number)
+
+#        mpi.report( "Starting Fortran solver %(name)s"%self.__dict__)
+
+        self.Sigma_Old <<= self.Sigma
+        self.G_Old <<= self.G
+
+#        # call the fortran solver:
+#        temp = 1.0/self.beta
+#        gf,tail,self.atocc,self.atmag = gf_hi_fullu(e0f=self.ealmat, ur=ur, umn=umn, ujmn=ujmn,
+#                                                    zmsb=self.zmsb, nmom=self.Nmoments, ns=self.Nspin, temp=temp, verbosity = self.verbosity)
+
+        #self.sig = sigma_atomic_fullu(gf=self.gf,e0f=self.eal,zmsb=self.zmsb,ns=self.Nspin,nlm=self.Nlm)
+        def print_matrix(m):
+            for row in m:
+                print ''.join(map("{0:12.7f}".format, row))
+
+
+# Hartree-Fock solver
+        self.Sigma.zero()
+        dm = self.G.density()
+        if mpi.is_master_node():
+#            print
+#            print "  Reduced U-matrix:"
+#            print "  U:"
+#            print_matrix(ujmn)
+#            print "  Up:"
+#            print_matrix(umn)
+#
+##            sig_test = {bl: numpy.zeros((self.Nlm, self.Nlm)) for bl in dm}
+#            sig_test = {}
+#            sig_test['up'] = numpy.dot(umn, dm['up'].real) + numpy.dot(ujmn, dm['down'].real)
+#            sig_test['down'] = numpy.dot(umn, dm['down'].real) + numpy.dot(ujmn, dm['up'].real)
+#            print "  Sigma test:"
+#            print_matrix(sig_test['up'])
+            print
+            print "  Density matrix (up):"
+            print_matrix(dm['up'])
+            print
+            print "  Density matrix (down):"
+            print_matrix(dm['down'])
+
+        if self.dudarev:
+            Ueff = U_int - J_hund
+            corr_energy = 0.0
+            dft_dc = 0.0
+            for bl1 in dm:
+# (U - J) * (1/2 - n)
+                self.Sigma[bl1] << Ueff * (0.5 * numpy.identity(self.Nlm) - dm[bl1])
+# 1/2 (U - J) * \sum_{\sig} [\sum_{m} n_{m,m \sig} - \sum_{m1,m2} n_{m1,m2 \sig} n_{m2,m1 \sig}]
+                corr_energy += 0.5 * Ueff * (dm[bl1].trace() - (dm[bl1]*dm[bl1].conj()).sum()).real
+# V[n] * n^{\dagger}
+                dft_dc += (self.Sigma[bl1](0) * dm[bl1].conj()).sum().real
+
+        else:
+# !!!!!
+# !!!!! Mind the order of indices in the 4-index matrix!
+# !!!!!
+            for il1 in xrange(self.Nlm):
+                for il2 in xrange(self.Nlm):
+                    for il3 in xrange(self.Nlm):
+                        for il4 in xrange(self.Nlm):
+                            for bl1 in dm:
+                                for bl2 in dm:
+                                    self.Sigma[bl1][il1, il2] += ur[il1, il3, il2, il4] * dm[bl2][il3, il4]
+                                    if bl1 == bl2:
+                                        self.Sigma[bl1][il1, il2] -= ur[il1, il3, il4, il2] * dm[bl1][il3, il4]
+
+        if mpi.is_master_node() and self.verbosity > 0:
+            print
+            print "  Sigma (up):"
+            print_matrix(self.Sigma['up'](0).real)
+            print
+            print "  Sigma (down):"
+            print_matrix(self.Sigma['down'](0).real)
+#        if (self.verbosity==0):
+#            # No fortran output, so give basic results here
+#            mpi.report("Atomic occupancy in Hubbard I Solver  : %s"%self.atocc)
+#            mpi.report("Atomic magn. mom. in Hubbard I Solver : %s"%self.atmag)
+
+        # transfer the data to the GF class:
+        if (self.UseSpinOrbit):
+            nlmtot = self.Nlm*2         # only one block in this case!
+        else:
+            nlmtot = self.Nlm
+
+#        M={}
+#        isp=-1
+#        for a,al in self.gf_struct:
+#            isp+=1
+#            M[a] = numpy.array(gf[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot,:]).transpose(2,0,1).copy()
+#            for i in range(min(self.Nmoments,8)):
+#                tailtempl[a][i+1] = tail[i][isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot]
+#
+#        #glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb, data =M[a], tail =self.tailtempl[a])
+#        #                   for a,al in self.gf_struct]
+#        glist = lambda : [ GfImFreq(indices = al, beta = self.beta, n_points = self.Nmsb) for a,al in self.gf_struct]
+#        self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
+#        
+#        self.__copy_Gf(self.G,M,tailtempl)
+#
+#        # Self energy:
+#        self.G0 <<= iOmega_n
+#
+#        M = [ self.ealmat[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot] for isp in range((2*self.Nlm)/nlmtot) ]
+#        self.G0 -= M
+#        self.Sigma <<= self.G0 - inverse(self.G)
+#
+#        # invert G0
+#        self.G0.invert()
+
+        def test_distance(G1,G2, dist) :
+            def f(G1,G2) :
+                #print abs(G1.data - G2.data)
+                dS = max(abs(G1.data - G2.data).flatten())
+                aS = max(abs(G1.data).flatten())
+                if mpi.is_master_node():
+                    print "  Distances:", dS, " vs ", aS * dist
+                return dS <= aS*dist
+            return reduce(lambda x,y : x and y, [f(g1,g2) for (i1,g1),(i2,g2) in izip(G1,G2)])
+
+        mpi.report("\nChecking Sigma for convergence...\nUsing tolerance %s"%Test_Convergence)
+        self.Converged = test_distance(self.Sigma,self.Sigma_Old,Test_Convergence)
+
+        if self.Converged :
+            mpi.report("Solver HAS CONVERGED")
+        else :
+            mpi.report("Solver has not yet converged")
+
+        return corr_energy, dft_dc
+
+
+    def GF_realomega(self, ommin, ommax, N_om, U_int, J_hund, T=None, verbosity=0, broadening=0.01):
+        """Calculates the GF and spectral function on the real axis."""
+
+        delta_om = (ommax-ommin)/(1.0*(N_om-1))
+
+        omega = numpy.zeros([N_om],numpy.complex_)
+
+        ur,umn,ujmn=self.__set_umatrix(U=U_int,J=J_hund,T=T)
+
+        for i in range(N_om):
+            omega[i] = ommin + delta_om * i + 1j * broadening
+
+        tailtempl={}
+        for sig,g in self.G:
+            tailtempl[sig] = copy.deepcopy(g.tail)
+            for i in range(9): tailtempl[sig][i] *= 0.0
+
+        temp = 1.0/self.beta
+        gf,tail,self.atocc,self.atmag = gf_hi_fullu(e0f=self.ealmat, ur=ur, umn=umn, ujmn=ujmn,
+                                                    zmsb=omega, nmom=self.Nmoments, ns=self.Nspin, temp=temp, verbosity = verbosity)
+
+        # transfer the data to the GF class:
+        if (self.UseSpinOrbit):
+            nlmtot = self.Nlm*2         # only one block in this case!
+        else:
+            nlmtot = self.Nlm
+
+        M={}
+        isp=-1
+        for a,al in self.gf_struct:
+            isp+=1
+            M[a] = numpy.array(gf[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot,:]).transpose(2,0,1).copy()
+            for i in range(min(self.Nmoments,8)):
+                tailtempl[a][i+1] = tail[i][isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot]
+
+        #glist = lambda : [ GfReFreq(indices = al, window = (ommin, ommax), n_points = N_om, data = M[a], tail = self.tailtempl[a])
+        #                   for a,al in self.gf_struct]       # Indices for the upfolded G
+        glist = lambda : [ GfReFreq(indices = al, window = (ommin, ommax), n_points = N_om) for a,al in self.gf_struct]
+        self.G = BlockGf(name_list = self.a_list, block_list = glist(),make_copies=False)
+
+        self.__copy_Gf(self.G,M,tailtempl)
+
+        # Self energy:
+        self.G0 = self.G.copy()
+        self.Sigma = self.G.copy()
+        self.G0 <<= Omega + 1j*broadening
+
+        M = [ self.ealmat[isp*nlmtot:(isp+1)*nlmtot,isp*nlmtot:(isp+1)*nlmtot] for isp in range((2*self.Nlm)/nlmtot) ]
+        self.G0 -= M
+        self.Sigma <<= self.G0 - inverse(self.G)
+        self.Sigma.note='ReFreq'          # This is important for the put_Sigma routine!!!
+
+        #sigmamat = sigma_atomic_fullu(gf=gf,e0f=self.ealmat,zmsb=omega,nlm=self.Nlm,ns=self.Nspin)
+
+        #return omega,gf,sigmamat
+
+
+    def __save_eal(self,Filename,it):
+        f=open(Filename,'a')
+        f.write('\neff. atomic levels, Iteration %s\n'%it)
+        for i in range(self.Nlm*self.Nspin):
+            for j in range(self.Nlm*self.Nspin):
+                f.write("%10.6f %10.6f   "%(self.ealmat[i,j].real,self.ealmat[i,j].imag))
+            f.write("\n")
+        f.close()
+
+    def __copy_Gf(self,G,data,tail):
+        """ Copies data and tail to Gf object GF """
+        for s,g in G:
+            g.data[:,:,:]=data[s][:,:,:]
+            for imom in range(1,min(self.Nmoments,8)):
+                g.tail.data[1+imom,:,:]=tail[s][imom]
+
+
+    def set_atomic_levels(self,eal):
+        """ Helps to set correctly the variables for the atomic levels from a dictionary."""
+
+        assert (type(eal)==DictType), "Give a dictionary to set_atomic_levels!"
+
+        cnt = 0
+        self.ealmat[:,:] *= 0.0
+
+        for ind in eal:
+            self.Eff_Atomic_Levels[ind] = copy.deepcopy(eal[ind])
+
+            if self.UseSpinOrbit:
+                for ii in range(self.Nlm*2):
+                    for jj in range(self.Nlm*2):
+                        self.ealmat[ii,jj] = self.Eff_Atomic_Levels[ind][ii,jj]
+            else:
+                for ii in range(self.Nlm):
+                    for jj in range(self.Nlm):
+                        self.ealmat[cnt*self.Nlm + ii,cnt*self.Nlm + jj] = self.Eff_Atomic_Levels[ind][ii,jj]
+
+            cnt += 1
+
+
+    def __set_umatrix(self,U,J,T=None):
+        # U matrix:
+        # l = (Nlm-1)/2
+        # If T is specified, it is used to transform the Basis set
+        Umat = U_matrix(l=self.l, U_int=U, J_hund=J, basis='cubic', T=T)
+        U, Up = reduce_4index_to_2index(Umat)
+
+        return Umat, U, Up

python/converters/plovasp/examples/ce/plo.cfg

@@ -0,0 +1,15 @@
+[General]
+#EFERMI = 
+#BASENAME = plo_vasp
+#DOSMESH = -5.0 10.0 4001
+
+[Group 1]
+SHELLS = 1
+NORMALIZE = False
+EWINDOW = -41.45  41.8
+
+[Shell 1]
+LSHELL = 3
+IONS = 1
+
+

python/converters/plovasp/examples/ce/run_build.sh

@@ -0,0 +1,2 @@
+#!/bin/bash
+PYTHONPATH=$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src/python:$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src/c:$PYTHONPATH pytriqs $@

python/converters/plovasp/examples/ce/run_plovasp.sh

@@ -0,0 +1,2 @@
+PLOVASP_PATH=$HOME/Codes/triqs/triqs1.2/applications/dft_tools/dev_src
+PYTHONPATH=$PLOVASP_PATH/python:$PLOVASP_PATH/c:$PYTHONPATH pytriqs $PLOVASP_PATH/python/converters/plovasp/main.py $@

python/converters/plovasp/examples/ce/sc_dmft.py

@@ -0,0 +1,164 @@
+
+import os
+import errno
+import signal
+import sys
+import time
+import pytriqs.utility.mpi as mpi
+import converters.plovasp.main as plovasp
+from test_ham_hf import dmft_cycle
+
+debug = True
+#
+# Helper functions
+#
+def sigint_handler(signal, frame):
+    raise SystemExit(1)
+
+def is_vasp_lock_present():
+    return os.path.isfile('./vasp.lock')
+
+def is_vasp_running(vasp_pid):
+    """
+    Tests if VASP initial process is still alive.
+    """
+    pid_exists = False
+    if mpi.is_master_node():
+        try:
+            os.kill(vasp_pid, 0)
+        except OSError, e:
+            pid_exists = e.errno == errno.EPERM
+        else:
+            pid_exists = True
+
+    pid_exists = mpi.bcast(pid_exists)
+    return pid_exists
+
+def get_dft_energy():
+    """
+    Reads energy from the last line of OSZICAR.
+    """
+    with open('OSZICAR', 'r') as f:
+        nextline = f.readline()
+        while nextline.strip():
+            line = nextline
+            nextline = f.readline()
+#            print "OSZICAR: ", line[:-1]
+
+    try:
+        dft_energy = float(line.split()[2])
+    except ValueError:
+        print "Cannot read energy from OSZICAR, setting it to zero"
+        dft_energy = 0.0
+
+    return dft_energy
+
+class bcolors:
+    MAGENTA = '\033[95m'
+    BLUE = '\033[94m'
+    GREEN = '\033[92m'
+    YELLOW = '\033[93m'
+    RED = '\033[91m'
+    ENDC = '\033[0m'
+#
+# Main self-consistent cycle
+#
+def run_all(vasp_pid):
+    """
+    """
+    mpi.report("  Waiting for VASP lock to appear...")
+    while not is_vasp_lock_present():
+        time.sleep(1)
+
+    vasp_running = True
+
+    while vasp_running:
+        if debug: print bcolors.RED + "rank %s"%(mpi.rank) + bcolors.ENDC
+        mpi.report("  Waiting for VASP lock to disappear...")
+        mpi.barrier()
+        while is_vasp_lock_present():
+            time.sleep(1)
+#            if debug: print bcolors.YELLOW + " waiting: rank %s"%(mpi.rank) + bcolors.ENDC
+            if not is_vasp_running(vasp_pid):
+                mpi.report("  VASP stopped")
+                vasp_running = False
+                break
+
+        if debug: print bcolors.MAGENTA + "rank %s"%(mpi.rank) + bcolors.ENDC
+        err = 0
+        exc = None
+        try:
+            if debug: print bcolors.BLUE + "plovasp: rank %s"%(mpi.rank) + bcolors.ENDC
+            if mpi.is_master_node():
+                plovasp.generate_and_output_as_text('plo.cfg', vasp_dir='./')
+# Read energy from OSZICAR
+                dft_energy = get_dft_energy()
+        except Exception, exc:
+            err = 1
+
+        err = mpi.bcast(err)
+        if err:
+            if mpi.is_master_node():
+                raise exc
+            else:
+                raise SystemExit(1)
+
+        mpi.barrier()
+
+        try:
+            if debug: print bcolors.GREEN + "rank %s"%(mpi.rank) + bcolors.ENDC
+            corr_energy, dft_dc = dmft_cycle()
+        except:
+            if mpi.is_master_node():
+                print "  master forwarding the exception..."
+                raise
+            else:
+                print "  rank %i exiting..."%(mpi.rank)
+                raise SystemExit(1)
+        mpi.barrier()
+
+        if mpi.is_master_node():
+            total_energy = dft_energy + corr_energy - dft_dc
+            print
+            print "="*80
+            print "  Total energy: ", total_energy
+            print "  DFT energy: ", dft_energy
+            print "  Corr. energy: ", corr_energy
+            print "  DFT DC: ", dft_dc
+            print "="*80
+            print
+
+        if mpi.is_master_node() and vasp_running:
+            open('./vasp.lock', 'a').close()
+            
+    if mpi.is_master_node():
+        total_energy = dft_energy + corr_energy - dft_dc
+        with open('TOTENERGY', 'w') as f:
+            f.write("  Total energy: %s\n"%(total_energy))
+            f.write("  DFT energy: %s\n"%(dft_energy))
+            f.write("  Corr. energy: %s\n"%(corr_energy))
+            f.write("  DFT DC: %s\n"%(dft_dc))
+            f.write("  Energy correction: %s\n"%(corr_energy - dft_dc))
+
+    mpi.report("***Done")
+
+
+if __name__ == '__main__':
+    try:
+        vasp_pid = int(sys.argv[1])
+    except (ValueError, KeyError):
+        if mpi.is_master_node():
+            print "VASP process pid must be provided as the first argument"
+        raise
+#    if len(sys.argv) > 1:
+#        vasp_path = sys.argv[1]
+#    else:
+#        try:
+#            vasp_path = os.environ['VASP_DIR']
+#        except KeyError:
+#            if mpi.is_master_node():
+#                print "Path to VASP must be specified either as an argument of in VASP_DIR"
+#            raise
+    signal.signal(signal.SIGINT, sigint_handler)
+
+    run_all(vasp_pid)

python/converters/plovasp/examples/ce/sc_dmft.sh

@@ -0,0 +1,9 @@
+#/bin/bash
+
+NPROC=8
+
+rm -f vasp.lock
+stdbuf -o 0 mpirun -np $NPROC ~/Codes/vasp/build_20151202/vasp.5.4.2/build/std/vasp &
+
+mpirun -np $NPROC ./run_build.sh sc_dmft.py $(jobs -p) || kill %1
+

python/converters/plovasp/examples/ce/subm.sh

@@ -0,0 +1,28 @@
+#!/bin/bash -l
+#
+#SBATCH --job-name="JOBNAME"
+#SBATCH --time=36:00:00
+##SBATCH --nodes=8
+#SBATCH --ntasks=8
+#SBATCH --partition=qmac
+##SBATCH --mem-per-cpu=1024
+#SBATCH --output=job.%j.o
+#SBATCH --error=job.%j.e
+
+#======START=====
+# source $HOME/.profile
+# export KMP_AFFINITY=none
+
+# module load slurm
+source /cm/shared/apps/intel/composer_xe/composer_xe_2015.0.090/mkl/bin/mklvars.sh intel64
+# module switch PrgEnv-gnu PrgEnv-intel
+# module load intel/14.0.2.144
+#module load intel/mkl/64/11.2/2015.0.090
+
+export OMP_NUM_THREADS=1
+
+VASP_DIR=$HOME/Codes/vasp/build_20151202/vasp.5.4.2/build/std
+#mpirun -np $SLURM_NTASKS $VASP_DIR/vasp
+./sc_dmft.sh
+
+

python/converters/plovasp/examples/ce/test_ham_hf.py

@@ -0,0 +1,444 @@
+#from pytriqs.applications.dft.sumk_dft import *
+from sumk_dft import *
+#from pytriqs.applications.dft.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():
+                ar = HDFArchive(self.hdf_file,'r')
+                fermi_weights = ar['dft_misc_input']['dft_fermi_weights']
+                band_window = ar['dft_misc_input']['band_window']
+                del ar
+            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():
+            ar = HDFArchive(self.hdf_file,'r')
+            fermi_weights = ar['dft_misc_input']['dft_fermi_weights']
+            band_window = ar['dft_misc_input']['band_window']
+            del ar
+        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():
+        ar = HDFArchive(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=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()):
+            ar = HDFArchive(HDFfilename,'a')
+            S.Sigma <<= ar['SigmaF']
+            del ar
+            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)):
+                    ar = HDFArchive(HDFfilename,'r')
+                    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']
+                    del ar
+                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()):
+                ar = HDFArchive(HDFfilename,'a')
+                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
+                del ar
+
+            #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()):
+        ar = HDFArchive(HDFfilename)
+        itn = ar['iterations']
+        correnerg = ar['correnerg%s'%itn]
+        DCenerg = ar['DCenerg%s'%itn]
+        del ar
+        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()