diff --git a/doc/tutorials/images_scripts/NiO_local_lattice_GF.py b/doc/tutorials/images_scripts/NiO_local_lattice_GF.py index 02dbe051..21cf2eb6 100644 --- a/doc/tutorials/images_scripts/NiO_local_lattice_GF.py +++ b/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 diff --git a/doc/tutorials/images_scripts/maxent.py b/doc/tutorials/images_scripts/maxent.py index 8795990f..841d8067 100644 --- a/doc/tutorials/images_scripts/maxent.py +++ b/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']: diff --git a/doc/tutorials/images_scripts/nio.py b/doc/tutorials/images_scripts/nio.py index 27f47f6a..0195e797 100644 --- a/doc/tutorials/images_scripts/nio.py +++ b/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) diff --git a/doc/tutorials/images_scripts/nio_Aw.png b/doc/tutorials/images_scripts/nio_Aw.png index 35b9b0c3..45b62fd0 100644 Binary files a/doc/tutorials/images_scripts/nio_Aw.png and b/doc/tutorials/images_scripts/nio_Aw.png differ diff --git a/doc/tutorials/images_scripts/nio_csc.py b/doc/tutorials/images_scripts/nio_csc.py new file mode 100644 index 00000000..993293da --- /dev/null +++ b/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 diff --git a/doc/tutorials/images_scripts/nio_csc.py.rst b/doc/tutorials/images_scripts/nio_csc.py.rst new file mode 100644 index 00000000..cd339bfe --- /dev/null +++ b/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 diff --git a/doc/tutorials/images_scripts/plo.cfg b/doc/tutorials/images_scripts/plo.cfg index 91c72333..c370bf17 100644 --- a/doc/tutorials/images_scripts/plo.cfg +++ b/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 diff --git a/doc/tutorials/nio_csc.rst b/doc/tutorials/nio_csc.rst index 1192a3bc..23f60f43 100644 --- a/doc/tutorials/nio_csc.rst +++ b/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`, 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`, where `` 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 -i -p nio_csc.py