tests for plovasp's H(k) and complement. Further work on NiO tutorial

doc/tutorials/images_scripts/NiO_local_lattice_GF.py.rst

@@ -0,0 +1,9 @@
+.. _NiO_local_lattice_GF.py:
+
+NiO_local_lattice_GF.py
+-----------
+
+Download :download:`NiO_local_lattice_GF.py <./NiO_local_lattice_GF.py>`.
+
+.. literalinclude:: NiO_local_lattice_GF.py
+   :language: python

doc/tutorials/images_scripts/maxent.py

@@ -0,0 +1,51 @@
+from pytriqs.gf import *
+from pytriqs.archive import *
+from triqs_maxent import *
+
+filename = 'vasp'
+
+ar = HDFArchive(filename+'.h5','a')
+if 'iteration_count' in ar['DMFT_results']:
+    iteration_offset = ar['DMFT_results']['iteration_count']+1
+    G_latt = ar['DMFT_results']['Iterations']['G_latt_orb_it'+str(iteration_offset-1)]
+
+
+tm = TauMaxEnt(cost_function='bryan', probability='normal')
+
+print(G_latt['up'][0,0])
+t2g_orbs = [0,1,3]
+eg_orbs = [2,4]
+op_orbs = [5,6,7]
+
+orbs = [t2g_orbs, eg_orbs, op_orbs]
+#orbs = [t2g_orbs]
+
+for orb in orbs:
+
+    print '\n'+str(orb[0])+'\n'
+
+    gf = 0*G_latt['up'][0,0]
+    for iO in orb:
+        gf = gf + G_latt['up'][iO,iO]
+    tm.set_G_iw(gf)
+    tm.omega =LinearOmegaMesh(omega_min=-20, omega_max=20, n_points=201)
+    tm.alpha_mesh = LogAlphaMesh(alpha_min=0.01, alpha_max=20000, n_points=60)
+
+    tm.set_error(1.e-3)
+    result=tm.run()
+    result.get_A_out('LineFitAnalyzer')
+    
+    if 'iteration_count' in ar['DMFT_results']:
+        iteration_offset = ar['DMFT_results']['iteration_count']+1
+        for oo in orb:
+            ar['DMFT_results']['Iterations']['G_latt_orb_w_o'+str(oo)+'_it'+str(iteration_offset-1)] = result.analyzer_results['LineFitAnalyzer']['A_out']
+        ar['DMFT_results']['Iterations']['w_it'+str(iteration_offset-1)] = result.omega
+
+
+    # you may be interested in the details of the line analyzer:
+    # from pytriqs.plot.mpl_interface import oplot
+    #plt.figure(2)
+    #result.analyzer_results['LineFitAnalyzer'].plot_linefit()
+    #plt.savefig('ana'+str(orb[0])+'.pdf',fmt='pdf')
+
+del ar

doc/tutorials/images_scripts/maxent.py.rst

@@ -0,0 +1,9 @@
+.. _maxent.py:
+
+maxent.py
+-----------
+
+Download :download:`maxent.py <./maxent.py>`.
+
+.. literalinclude:: maxent.py
+   :language: python

doc/tutorials/images_scripts/nio.py

@@ -3,11 +3,17 @@ import numpy as np
 import pytriqs.utility.mpi as mpi
 from pytriqs.archive import *
 from pytriqs.gf import *
+import sys, pytriqs.version as triqs_version
 from triqs_dft_tools.sumk_dft import *
 from triqs_dft_tools.sumk_dft_tools import *
 from pytriqs.operators.util.hamiltonians import *
 from pytriqs.operators.util.U_matrix import *
 from triqs_cthyb import *
+import triqs_cthyb.version as cthyb_version
+import triqs_dft_tools.version as dft_tools_version
+
+import warnings
+warnings.filterwarnings("ignore", category=FutureWarning)
 
 filename = 'vasp'
 
@@ -32,7 +38,6 @@ for i_sh in range(len(SK.deg_shells)):
 n_orb = SK.corr_shells[0]['dim']
 spin_names = ['up','down']
 orb_names = [i for i in range(0,n_orb)]
-orb_hyb = False
 
 #gf_struct = set_operator_structure(spin_names, orb_names, orb_hyb)
 gf_struct = SK.gf_struct_solver[0]
@@ -74,6 +79,7 @@ p["perform_tail_fit"] = True
 
 # Double Counting: 0 FLL, 1 Held, 2 AMF
 DC_type = 0
+DC_value = 59.0
 
 # Prepare hdf file and and check for previous iterations
 n_iterations = 10
@@ -83,13 +89,23 @@ if mpi.is_master_node():
     ar = HDFArchive(filename+'.h5','a')
     if not 'DMFT_results' in ar: ar.create_group('DMFT_results')
     if not 'Iterations' in ar['DMFT_results']: ar['DMFT_results'].create_group('Iterations')
+    if not 'DMFT_input' in ar: ar.create_group('DMFT_input')
+    if not 'Iterations' in ar['DMFT_input']: ar['DMFT_input'].create_group('Iterations')
+    if not 'code_versions' in ar['DMFT_input']: ar['DMFT_input'].create_group('code_versio\
+ns')
+    ar['DMFT_input']['code_versions']["triqs_version"] = triqs_version.version
+    ar['DMFT_input']['code_versions']["triqs_git"] = triqs_version.git_hash
+    ar['DMFT_input']['code_versions']["cthyb_version"] = cthyb_version.version
+    ar['DMFT_input']['code_versions']["cthyb_git"] = cthyb_version.cthyb_hash
+    ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
+    ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.dft_tools_hash
     if 'iteration_count' in ar['DMFT_results']: 
         iteration_offset = ar['DMFT_results']['iteration_count']+1
         S.Sigma_iw = ar['DMFT_results']['Iterations']['Sigma_it'+str(iteration_offset-1)]
         SK.dc_imp = ar['DMFT_results']['Iterations']['dc_imp'+str(iteration_offset-1)]
         SK.dc_energ = ar['DMFT_results']['Iterations']['dc_energ'+str(iteration_offset-1)]
         SK.chemical_potential = ar['DMFT_results']['Iterations']['chemical_potential'+str(iteration_offset-1)].real
-
+    ar['DMFT_input']["dmft_script_it"+str(iteration_offset)] = open(sys.argv[0]).read()
 iteration_offset = mpi.bcast(iteration_offset)
 S.Sigma_iw = mpi.bcast(S.Sigma_iw)
 SK.dc_imp = mpi.bcast(SK.dc_imp)
@@ -97,6 +113,7 @@ SK.dc_energ = mpi.bcast(SK.dc_energ)
 SK.chemical_potential = mpi.bcast(SK.chemical_potential)
 
 # Calc the first G0
+SK.symm_deg_gf(S.Sigma_iw,orb=0)
 SK.put_Sigma(Sigma_imp = [S.Sigma_iw])
 SK.calc_mu(precision=0.01)
 S.G_iw << SK.extract_G_loc()[0]
@@ -105,7 +122,7 @@ SK.symm_deg_gf(S.G_iw, orb=0)
 #Init the DC term and the self-energy if no previous iteration was found
 if iteration_offset == 0:
     dm = S.G_iw.density()
-    SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type)
+    SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value)
     S.Sigma_iw << SK.dc_imp[0]['up'][0,0]
 
 mpi.report('%s DMFT cycles requested. Starting with iteration %s.'%(n_iterations,iteration_offset))
@@ -120,12 +137,13 @@ for it in range(iteration_offset, iteration_offset + n_iterations):
     # Solve the impurity problem
     S.solve(h_int = H, **p)
     if mpi.is_master_node(): 
+        ar['DMFT_input']['Iterations']['solver_dict_it'+str(it)] = p
         ar['DMFT_results']['Iterations']['Gimp_it'+str(it)] = S.G_iw
         ar['DMFT_results']['Iterations']['Gtau_it'+str(it)] = S.G_tau
         ar['DMFT_results']['Iterations']['Sigma_uns_it'+str(it)] = S.Sigma_iw
     # Calculate double counting
     dm = S.G_iw.density()
-    SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type)
+    SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value)
     # Get new G
     SK.symm_deg_gf(S.Sigma_iw,orb=0)
     SK.put_Sigma(Sigma_imp=[S.Sigma_iw])

doc/tutorials/nio_csc.rst

@@ -3,7 +3,7 @@
 DFT and projections
 ==================================================
 
-We will perform charge self-consitent DFT+DMFT calcluations for the charge-transfer insulator NiO. We still start from scratch and provide all necessary input files to do the calculations.
+We will perform charge self-consitent DFT+DMFT calcluations for the charge-transfer insulator NiO. We start from scratch and provide all necessary input files to do the calculations: First for doing a single-shot calculation and then for charge-selfconsistency.
 
 VASP setup
 -------------------------------
@@ -31,7 +31,7 @@ We do this by invoking :program:`plovasp plo.cfg` which is configured by an inpu
 
 .. literalinclude:: images_scripts/plo.cfg
 
-Here, in `[General]' we set the basename and the grid for calculating the density of
+Here, in `[General]` we set the basename and the grid for calculating the density of
 states. In `[Group 1]` we define a group of two shells which are orthonormalized with
 respect to states in an energy window from `-9` to `2` for all ions simultanously
 (`NORMION = False`). We define the two shells, which correspond to the Ni d states
@@ -52,8 +52,23 @@ DMFT
 dmft script
 ------------------------------
 
-Since the python script for performing the dmft loop pretty much resembles that presented in the tutorial on :ref:`srvo3`, we will not go into detail here but simply provide the script :ref:`nio.py`. Following Kunes et al. in `PRB 75 165115 (2007) <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.75.165115>`_ we use :math:`U=8` and :math:`J=1`. Here, we use :math:`\beta=5` instead of :math:`\beta=10` to speed up the calculations.
+Since the python script for performing the dmft loop pretty much resembles that presented in the tutorial on :ref:`srvo3`, we will not go into detail here but simply provide the script :ref:`nio.py`. Following Kunes et al. in `PRB 75 165115 (2007) <https://journals.aps.org/prb/abstract/10.1103/PhysRevB.75.165115>`_ we use :math:`U=8` and :math:`J=1`. We se;ect :math:`\beta=5` instead of :math:`\beta=10` to ease the problem slightly. For simplicity we fix the double-counting potential to :math:`\mu_{DC}=59` eV by::
+
+  DC_value = 59.0
+  SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_value=DC_value)
+
+For sensible results run this script in parallel on at least 20 cores. As a quick check of the results, we can compare the orbital occupation from the paper cited above (:math:`n_{eg} = 0.54` and :math:`n_{t2g}=1.0`) and those from the cthyb output (check lines `Orbital up_0 density:` for a t2g  and `Orbital up_2 density:` for an eg orbital). They should coincide well.
+
 
 Local lattice Green's function for all projected orbitals
 ----------------------
-We calculate the local lattice Green's function - now also for the uncorrelated orbitals, i.e., the O p states. Therefor we use :download:`NiO_local_lattice_GF.py <images_scripts/NiO_local_lattice_GF.py`
+We calculate the local lattice Green's function - now also for the uncorrelated orbitals, i.e., the O p states, for what we use the script :ref:`NiO_local_lattice_GF.py`. The result is saved in the h5 file as `G_latt_orb_it<n_it>`, where `n_it>` is the number of the last DMFT iteration.
+
+Spectral function on real axis: MaxEnt
+----------------------
+To compare to results from literature we make use of the `maxent triqs application <https://triqs.github.io/maxent/master/>`_ and calculate the spectral function on real axis. Use this script to perform a crude but quick calculation: :ref:`maxent.py` using a linear real axis and a line-fit analyzer to determine the optimal :math:`\alpha`. The output is saved in the h5 file in `DMFT_results/Iterations/G_latt_orb_w_o<n_o>_it<n_it>`, where `<n_o>` is the number of the orbital and `n_it` is again the number of the last iteration. The real axis information is stored in `DMFT_results/Iterations/w_it<n_it>`.
+
+
+.. image:: images_scripts/nio_Aw.png
+    :width: 400
+    :align: center

python/converters/plovasp/plotools.py

@@ -128,10 +128,14 @@ def generate_plo(conf_pars, el_struct):
     for gr_par in conf_pars.groups:
         pgroup = ProjectorGroup(gr_par, pshells, eigvals)
         pgroup.orthogonalize()
-        if gr_par['complement']:
-            pgroup.complement(eigvals)
+        if pgroup.complement:
+            pgroup.calc_complement(eigvals)
         if conf_pars.general['hk']:
             pgroup.calc_hk(eigvals)
+            #testout = 'hk.out.h5'
+            #from pytriqs.archive import HDFArchive
+            #with HDFArchive(testout, 'w') as h5test:
+            #    h5test['hk'] = pgroup.hk
 # DEBUG output
         print "Density matrix:"
         nimp = 0.0

python/converters/plovasp/proj_group.py

@@ -62,11 +62,11 @@ class ProjectorGroup:
         self.ishells = gr_pars['shells']
         self.ortho = gr_pars['normalize']
         self.normion = gr_pars['normion']
+        self.complement = gr_pars['complement']
 
         self.shells = shells
 
 # Determine the minimum and maximum band numbers
-        ib_win, ib_min, ib_max = self.select_bands(eigvals)
         if 'bands' in gr_pars:
             nk, nband, ns_band = eigvals.shape
             ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32)
@@ -81,20 +81,15 @@ class ProjectorGroup:
         self.ib_min = ib_min
         self.ib_max = ib_max
         self.nb_max = ib_max - ib_min + 1
-        
-        print self.ib_win
-        print self.ib_min
-        print self.ib_max
-        print self.nb_max
+       
         
         
-        
-        if gr_pars['complement']:
+        if self.complement:
             n_bands = self.ib_win[:,:,1] - self.ib_win[:,:,0]+1
             n_orbs = sum([x.ndim for x in self.shells])
             assert np.all( n_bands == n_bands[0,0] ), "At each band the same number of bands has to be selected for calculating the complement (to end up with an equal number of orbitals at each k-point)."
             if n_orbs == n_bands[0,0]:
-                gr_pars['complement'] = False
+                self.complement = False
                 print "\nWARNING: The total number of orbitals in this group is  "
                 print "equal to the number of bands. Setting COMPLEMENT to FALSE!\n"
         
@@ -207,8 +202,6 @@ class ProjectorGroup:
 
         _, ns, nk, _, _ = self.shells[0].proj_win.shape
         
-        #print(p_mat.shape)
-        print block_maps, ndim
         self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128)
 # Note that 'ns' and 'nk' are the same for all shells
         for isp in xrange(ns):
@@ -228,7 +221,7 @@ class ProjectorGroup:
                         nlm = i2 - i1 + 1
                         shell = self.shells[ish]
                         p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
-                        
+                
                 self.hk[isp,ik,:,:] = np.dot(p_mat*eigvals[ik,bmin:bmax,isp],
                                         p_mat.transpose().conjugate())
 
@@ -238,7 +231,7 @@ class ProjectorGroup:
 # complement
 #
 ################################################################################
-    def complement(self,eigvals):
+    def calc_complement(self,eigvals):
         """
         Calculate the complement for a group of projectors.
             
@@ -282,19 +275,20 @@ class ProjectorGroup:
                         p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
                 orbs_done = 1*ndim
                 p_full[0,isp,ik,:ndim,:] = p_mat
-                while orbs_done < self.nb_max:                    
-                    proj_work = p_full[0,isp,ik,:,:] 
+                while orbs_done < self.nb_max:
 #We calculate the overlap of all bloch states: sum_l <n|l><l|m>
-                    overlap = np.dot(proj_work.transpose().conjugate(),proj_work)
+                    overlap = np.dot(p_full[0,isp,ik,:orbs_done,:].transpose().conjugate(),p_full[0,isp,ik,:orbs_done,:])
 # work is the projector onto the orthogonal complment <n| ( 1 - sum_l |l><l| ) |m>
                     work = np.eye(self.nb_max) - overlap
 # calculate the norm of the projected bloch function
-                    norm = np.sqrt(np.sum(work*work.transpose().conjugate(),axis=1))
+                    norm = np.sqrt(np.sum(work*work.transpose(),axis=1))
 # select the bloch function leading to the largest norm
                     max_ind = np.argmax(norm)
 # normalize and put it to the projectors
-                    p_full[0,isp,ik,orbs_done,:] = work[max_ind]/norm[max_ind]
+                    p_full[0,isp,ik,orbs_done,:] = work[:,max_ind].conjugate()/norm[max_ind]
+
                     orbs_done += 1
+
         sh_pars = {}
         sh_pars['lshell'] = -1
         sh_pars['ions'] = {'nion':1,'ion_list':[[1]]}
@@ -303,6 +297,7 @@ class ProjectorGroup:
         sh_pars['ib_min']  = bmin
         sh_pars['ib_max']  = bmax
         sh_pars['ib_win']  = self.ib_win
+        
         self.shells.append(ComplementShell(sh_pars,p_full[:,:,:,ndim:,:],False))
         self.ishells.append(self.ishells[-1]+1)
                                                 

python/converters/vasp_converter.py

@@ -379,7 +379,7 @@ class VaspConverter(ConverterTools):
             things_to_save = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required',
                           'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat',
                           'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping',
-                          'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr','proj_or_hk']
+                              'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr','proj_or_hk']
             if self.proj_or_hk == 'hk' or True:
                 things_to_save.append('proj_mat_csc')
             for it in things_to_save: ar[self.dft_subgrp][it] = locals()[it]

test/plovasp/proj_group/test_one_site.py

@@ -25,7 +25,8 @@ class TestProjectorGroup(mytest.MyTestCase):
 
     Scenarios:
     - **test** that orthogonalization is correct
-    - **test** that NORMION = True gives the same result
+    - **test** that NORMION = True gives the same results
+    - **test that HK = TRUE gives correct H(k)
     """
     def setUp(self):
         conf_file = _rpath + 'example.cfg'
@@ -88,5 +89,13 @@ class TestProjectorGroup(mytest.MyTestCase):
 #        self.assertFileEqual(testout, expected_file)
         expected_file = _rpath + 'projortho.out.h5'
         self.assertH5FileEqual(testout, expected_file)
-
-
+        
+    def test_hk(self):
+        self.proj_gr.orthogonalize()
+        self.proj_gr.calc_hk(self.eigvals)
+        
+        testout = _rpath + 'hk.test.h5'
+        with HDFArchive(testout, 'w') as h5test:
+            h5test['hk'] = self.proj_gr.hk
+        expected_file = _rpath + 'hk.out.h5'
+        self.assertH5FileEqual(testout, expected_file)

test/plovasp/proj_group/test_one_site_compl.py

@@ -0,0 +1,94 @@
+
+import os
+import rpath
+_rpath = os.path.dirname(rpath.__file__) + '/'
+
+import numpy as np
+from triqs_dft_tools.converters.plovasp.vaspio import VaspData
+from triqs_dft_tools.converters.plovasp.elstruct import ElectronicStructure
+from triqs_dft_tools.converters.plovasp.inpconf import ConfigParameters
+from triqs_dft_tools.converters.plovasp.proj_shell import ProjectorShell
+from triqs_dft_tools.converters.plovasp.proj_group import ProjectorGroup
+from pytriqs.archive import HDFArchive
+import mytest
+
+################################################################################
+#
+# TestProjectorGroup
+#
+################################################################################
+class TestProjectorGroupCompl(mytest.MyTestCase):
+    """
+    Class:
+
+    ProjectorGroupCompl(sh_pars, proj_raw)
+
+    Scenarios:
+    - **test** that unequal number of bands at different k-points gives error
+    - **test** that COMLEMENT=TRUE gives orthonormal projectors
+    """
+    def setUp(self):
+        conf_file = _rpath + 'example.cfg'
+        self.pars = ConfigParameters(conf_file)
+        self.pars.parse_input()
+        vasp_data = VaspData(_rpath + 'one_site/')
+        self.el_struct = ElectronicStructure(vasp_data)
+
+        efermi = self.el_struct.efermi
+        self.eigvals = self.el_struct.eigvals - efermi
+        
+        struct = self.el_struct.structure
+        kmesh = self.el_struct.kmesh
+        
+        self.proj_sh = ProjectorShell(self.pars.shells[0], vasp_data.plocar.plo, vasp_data.plocar.proj_params, kmesh, struct, 0)
+        
+
+    def test_num_bands(self):
+        self.pars.groups[0]['complement'] = True        
+        err_mess = "At each band the same number"
+        with self.assertRaisesRegexp(AssertionError, err_mess):
+            self.proj_gr = ProjectorGroup(self.pars.groups[0], [self.proj_sh], self.eigvals)        
+        
+    def test_compl(self):
+        self.pars.groups[0]['complement'] = True
+        self.pars.groups[0]['bands'] = [10, 25]
+        
+        self.proj_gr = ProjectorGroup(self.pars.groups[0], [self.proj_sh], self.eigvals)
+        
+        self.proj_gr.orthogonalize()
+        self.proj_gr.calc_complement(self.eigvals)
+        
+        temp = self.proj_gr.normion
+        self.proj_gr.normion = False
+        block_maps, ndim = self.proj_gr.get_block_matrix_map()
+        self.proj_gr.normion = temp
+
+        _, ns, nk, _, _ = self.proj_gr.shells[0].proj_win.shape
+        
+# Note that 'ns' and 'nk' are the same for all shells
+        for isp in xrange(ns):
+            for ik in xrange(nk):
+                print('ik',ik)
+                bmin = self.proj_gr.ib_win[ik, isp, 0]
+                bmax = self.proj_gr.ib_win[ik, isp, 1]+1
+                
+                nb = bmax - bmin 
+                p_mat = np.zeros((ndim, nb), dtype=np.complex128)
+                #print(bmin,bmax,nb)
+# Combine all projectors of the group to one block projector
+                for bl_map in block_maps:
+                    p_mat[:, :] = 0.0j  # !!! Clean-up from the last k-point and block!
+                    for ibl, block in enumerate(bl_map):
+                        i1, i2 = block['bmat_range']
+                        ish, ion = block['shell_ion']
+                        nlm = i2 - i1 + 1
+                        shell = self.proj_gr.shells[ish]
+                        p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
+                        
+                overlap_L = np.dot(p_mat.conjugate().transpose(),p_mat)
+                overlap_N = np.dot(p_mat,p_mat.conjugate().transpose())
+                
+                assert np.all(np.abs(np.eye(overlap_N.shape[0]) - overlap_N) < 1e-13)
+                assert np.all(np.abs(np.eye(overlap_L.shape[0]) - overlap_L) < 1e-13)
+
+