diff --git a/doc/tutorials/images_scripts/NiO_local_lattice_GF.py.rst b/doc/tutorials/images_scripts/NiO_local_lattice_GF.py.rst new file mode 100644 index 00000000..f8352c2e --- /dev/null +++ b/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 diff --git a/doc/tutorials/images_scripts/maxent.py b/doc/tutorials/images_scripts/maxent.py new file mode 100644 index 00000000..8795990f --- /dev/null +++ b/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 diff --git a/doc/tutorials/images_scripts/maxent.py.rst b/doc/tutorials/images_scripts/maxent.py.rst new file mode 100644 index 00000000..7f553d11 --- /dev/null +++ b/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 diff --git a/doc/tutorials/images_scripts/nio.py b/doc/tutorials/images_scripts/nio.py index 1d238840..27f47f6a 100644 --- a/doc/tutorials/images_scripts/nio.py +++ b/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]) diff --git a/doc/tutorials/images_scripts/nio_Aw.png b/doc/tutorials/images_scripts/nio_Aw.png new file mode 100644 index 00000000..35b9b0c3 Binary files /dev/null and b/doc/tutorials/images_scripts/nio_Aw.png differ diff --git a/doc/tutorials/nio_csc.rst b/doc/tutorials/nio_csc.rst index a2d7e94d..1192a3bc 100644 --- a/doc/tutorials/nio_csc.rst +++ b/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) `_ 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) `_ 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 `, 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 `_ 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_it`, where `` 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`. + + +.. image:: images_scripts/nio_Aw.png + :width: 400 + :align: center diff --git a/python/converters/plovasp/plotools.py b/python/converters/plovasp/plotools.py index e1c42313..ff014eed 100644 --- a/python/converters/plovasp/plotools.py +++ b/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 diff --git a/python/converters/plovasp/proj_group.py b/python/converters/plovasp/proj_group.py index 767f90fe..8bb766ce 100644 --- a/python/converters/plovasp/proj_group.py +++ b/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 - 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 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) diff --git a/python/converters/vasp_converter.py b/python/converters/vasp_converter.py index 69a6202f..9993c315 100644 --- a/python/converters/vasp_converter.py +++ b/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] diff --git a/test/plovasp/proj_group/hk.out.h5 b/test/plovasp/proj_group/hk.out.h5 new file mode 100644 index 00000000..93698c08 Binary files /dev/null and b/test/plovasp/proj_group/hk.out.h5 differ diff --git a/test/plovasp/proj_group/test_one_site.py b/test/plovasp/proj_group/test_one_site.py index cdb4d33c..d4307891 100644 --- a/test/plovasp/proj_group/test_one_site.py +++ b/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) \ No newline at end of file diff --git a/test/plovasp/proj_group/test_one_site_compl.py b/test/plovasp/proj_group/test_one_site_compl.py new file mode 100644 index 00000000..1d9daa65 --- /dev/null +++ b/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) + + \ No newline at end of file