uncompleted work on csc NiO tutorial

doc/tutorials/images_scripts/NiO_local_lattice_GF.py

@@ -11,7 +11,7 @@ from triqs_cthyb import *
 import warnings
 warnings.filterwarnings("ignore", category=FutureWarning)
 
-filename = 'vasp'
+filename = 'nio'
 SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
 
 beta = 5.0

doc/tutorials/images_scripts/maxent.py

@@ -2,7 +2,7 @@ from pytriqs.gf import *
 from pytriqs.archive import *
 from triqs_maxent import *
 
-filename = 'vasp'
+filename = 'nio'
 
 ar = HDFArchive(filename+'.h5','a')
 if 'iteration_count' in ar['DMFT_results']:

doc/tutorials/images_scripts/nio.py

@@ -15,7 +15,7 @@ import triqs_dft_tools.version as dft_tools_version
 import warnings
 warnings.filterwarnings("ignore", category=FutureWarning)
 
-filename = 'vasp'
+filename = 'nio'
 
 SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
 

doc/tutorials/images_scripts/nio_csc.py

@@ -0,0 +1,201 @@
+from itertools import *
+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
+from triqs_dft_tools.converters.vasp_converter import *
+
+
+import warnings
+warnings.filterwarnings("ignore", category=FutureWarning)
+
+
+def dmft_cycle():
+    filename = 'nio'
+    
+    Converter = VaspConverter(filename=filename)
+    Converter.convert_dft_input()
+    
+    SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
+    
+    beta = 5.0 
+     
+    Sigma = SK.block_structure.create_gf(beta=beta)
+    SK.put_Sigma([Sigma])
+    G = SK.extract_G_loc()
+    SK.analyse_block_structure_from_gf(G, threshold = 1e-2)
+    for i_sh in range(len(SK.deg_shells)):
+        num_block_deg_orbs = len(SK.deg_shells[i_sh])
+        mpi.report('found {0:d} blocks of degenerate orbitals in shell {1:d}'.format(num_block_deg_orbs, i_sh))
+        for iblock in range(num_block_deg_orbs):
+            mpi.report('block {0:d} consists of orbitals:'.format(iblock))
+            for keys in SK.deg_shells[i_sh][iblock].keys():
+                mpi.report('  '+keys)
+    
+    # Setup CTQMC Solver
+    
+    n_orb = SK.corr_shells[0]['dim']
+    spin_names = ['up','down']
+    orb_names = [i for i in range(0,n_orb)]
+    
+    #gf_struct = set_operator_structure(spin_names, orb_names, orb_hyb)
+    gf_struct = SK.gf_struct_solver[0]
+    mpi.report('Sumk to Solver: %s'%SK.sumk_to_solver)
+    mpi.report('GF struct sumk: %s'%SK.gf_struct_sumk)
+    mpi.report('GF struct solver: %s'%SK.gf_struct_solver)
+    
+    S = Solver(beta=beta, gf_struct=gf_struct)
+    
+    # Construct the Hamiltonian and save it in Hamiltonian_store.txt
+    H = Operator() 
+    U = 8.0
+    J = 1.0
+    
+    
+    U_sph = U_matrix(l=2, U_int=U, J_hund=J)
+    U_cubic = transform_U_matrix(U_sph, spherical_to_cubic(l=2, convention=''))
+    Umat, Upmat = reduce_4index_to_2index(U_cubic)
+    
+    H = h_int_density(spin_names, orb_names, map_operator_structure=SK.sumk_to_solver[0], U=Umat, Uprime=Upmat)
+    
+    # Print some information on the master node
+    mpi.report('Greens function structure is: %s '%gf_struct)
+    mpi.report('U Matrix set to:\n%s'%Umat)
+    mpi.report('Up Matrix set to:\n%s'%Upmat)
+    
+    # Parameters for the CTQMC Solver
+    p = {}
+    p["max_time"] = -1
+    p["random_name"] = ""
+    p["random_seed"] = 123 * mpi.rank + 567
+    p["length_cycle"] = 100
+    p["n_warmup_cycles"] = 2000
+    p["n_cycles"] = 20000
+    p["fit_max_moment"] = 4
+    p["fit_min_n"] = 30
+    p["fit_max_n"] = 50
+    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 = 1
+    
+    iteration_offset = 0
+    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)
+    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]
+    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,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))
+    
+    
+    
+    # The infamous DMFT self consistency cycle
+    for it in range(iteration_offset, iteration_offset + n_iterations):
+        mpi.report('Doing iteration: %s'%it)
+        
+        # Get G0
+        S.G0_iw << inverse(S.Sigma_iw + inverse(S.G_iw))
+        # 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,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])
+        SK.calc_mu(precision=0.01)
+        S.G_iw << SK.extract_G_loc()[0]
+        
+        # print densities
+        for sig,gf in S.G_iw:
+            mpi.report("Orbital %s density: %.6f"%(sig,dm[sig][0,0]))
+        mpi.report('Total charge of Gloc : %.6f'%S.G_iw.total_density())
+    
+        if mpi.is_master_node(): 
+            ar['DMFT_results']['iteration_count'] = it
+            ar['DMFT_results']['Iterations']['Sigma_it'+str(it)] = S.Sigma_iw
+            ar['DMFT_results']['Iterations']['Gloc_it'+str(it)] = S.G_iw
+            ar['DMFT_results']['Iterations']['G0loc_it'+str(it)] = S.G0_iw
+            ar['DMFT_results']['Iterations']['dc_imp'+str(it)] = SK.dc_imp
+            ar['DMFT_results']['Iterations']['dc_energ'+str(it)] = SK.dc_energ
+            ar['DMFT_results']['Iterations']['chemical_potential'+str(it)] = SK.chemical_potential
+    
+    
+    
+    
+    if mpi.is_master_node():
+        print 'calculating mu...'
+    SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
+    
+    if mpi.is_master_node():
+        print 'calculating GAMMA'
+    SK.calc_density_correction(dm_type='vasp')
+    
+    if mpi.is_master_node():
+        print 'calculating energy corrections'
+    
+    correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density()
+    
+    dm = S.G_iw.density() # compute the density matrix of the impurity problem
+    SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type,use_dc_value=DC_value)
+    dc_energ = SK.dc_energ[0]
+    
+    if mpi.is_master_node(): 
+        ar['DMFT_results']['Iterations']['corr_energy_it'+str(it)] = correnerg
+        ar['DMFT_results']['Iterations']['dc_energy_it'+str(it)] = dc_energ
+    
+    if mpi.is_master_node(): del ar
+        
+    return correnerg, dc_energ

doc/tutorials/images_scripts/nio_csc.py.rst

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

doc/tutorials/images_scripts/plo.cfg

@@ -7,6 +7,7 @@ SHELLS = 1 2
 EWINDOW = -9 2
 NORMION = False
 NORMALIZE = True
+BANDS = 2 10
 
 [Shell 1] # Ni d shell
 LSHELL = 2

doc/tutorials/nio_csc.rst

@@ -3,7 +3,8 @@
 DFT and projections
 ==================================================
 
-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.
+We will perform 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
 -------------------------------
@@ -62,7 +63,7 @@ For sensible results run this script in parallel on at least 20 cores. As a quic
 
 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, 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.
+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
 ----------------------
@@ -72,3 +73,27 @@ To compare to results from literature we make use of the `maxent triqs applicati
 .. image:: images_scripts/nio_Aw.png
     :width: 400
     :align: center
+
+.. Charge self-consistent DMFT
+.. ==================================================
+
+
+.. In this part we will perform charge self-consistent DMFT calculations. To do so we have to adapt the VASP `INCAR` such that :program:`VASP` reads the updated charge density after each step. We add the lines
+
+..   ICHARG = 5
+..   NELM = 1000
+..   NELMIN = 1000
+
+.. which makes VASP wait after each step of its iterative diagonalization until the file vasp.lock is created. It then reads the update of the charge density in the file `GAMMA`. It is terminated by an external script after a desired amount of steps, such that we deactivate all automatic stoping criterion by setting the number of steps to a very high number.
+
+.. We take the respective converged DFT and DMFT calculations from before as a starting point. I.e., we copy the `CHGCAR` and `nio.h5` together with the other :program:`VASP` input files and :file:`plo.cfg` in a new directory. We use a script called :program:`vasp_dmft` to invoke :program:`VASP` in the background and start the DMFT calculation together with :program:`plovasp` and the converter. This script assumes that the dmft sript contains a function `dmft_cycle()` and also the conversion from text files to the h5 file. Additionally we have to add a few lines to calculate the density correction and calculate the correlation energy. We adapt the script straightforardly and call it :ref:`nio_csc.py`. The most important new lines are
+
+..  SK.chemical_potential = SK.calc_mu( precision = 0.000001 )
+..  SK.calc_density_correction(dm_type='vasp')
+..  correnerg = 0.5 * (S.G_iw * S.Sigma_iw).total_density()
+
+.. where the chemical potential is determined to a greater precision than before, the correction to the dft density matrix is calculated and output to the file :file:`GAMMA`. The correlation energy is calculated via Migdal-Galitzki formula. We also slightly increase the tolerance for the detection of blocks since the DFT calculation now includes some QMC noise.
+
+.. We can start the whole machinery by excectuing
+
+..  vasp_dmft -n <n_procs> -i <n_iters> -p <vasp_exec>  nio_csc.py