LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2025-05-01 04:34:57 +02:00
Code Issues Releases Wiki Activity

[doc] NiO VASP CSC tutorial updates for new interface

Browse Source
This commit is contained in:
the-hampel 2025-02-19 20:09:15 +01:00
parent bd039be7b2
commit 24d102465b
7 changed files with 156 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
doc/tutorials/nio_csc_vasp/INCAR
View File

@ -1,7 +1,7 @@
System = NiO
# for PLOT DOS
# ISMEAR = -5
KPAR = 5
# better convergence for small kpt grids
ISMEAR = 2
SIGMA = 0.05
@ -16,9 +16,6 @@ EMAX = 10
LMAXMIX = 6
# switch off all symmetries
ISYM = -1
# project to Ni d and O p states
LORBIT = 14
LOCPROJ = "1 : d : Pr

32
doc/tutorials/nio_csc_vasp/INCAR.CSC Normal file
View File

@ -0,0 +1,32 @@
System = NiO
KPAR = 5
# better convergence for small kpt grids
ISMEAR = 2
SIGMA = 0.05
# converge wave functions
EDIFF = 1.E-7
NELMIN = 35
# the energy window to optimize projector channels (absolute)
EMIN = -3
EMAX = 10
LMAXMIX = 6
# project to Ni d and O p states
LORBIT = 14
LOCPROJ = "1 : d : Pr
2 : p : Pr"
# CSC flags
ICHARG = 5
NELM = 1000
NELMIN = 1000
NELMDL = -1
IMIX=1
BMIX=0.5
AMIX=0.02
LSYNCH5=True

7
doc/tutorials/nio_csc_vasp/INCAR.CSC.rst Normal file
View File

@ -0,0 +1,7 @@
.. _INCAR.CSC:
INCAR.CSC
-----------
.. literalinclude:: INCAR.CSC
:language: bash

2
doc/tutorials/nio_csc_vasp/KPOINTS
View File

@ -1,4 +1,4 @@
Automatically generated mesh
0
Gamma
6 6 6
11 11 11

104
doc/tutorials/nio_csc_vasp/nio.py
View File

@ -1,7 +1,6 @@
from itertools import *
import numpy as np
import triqs.utility.mpi as mpi
from h5 import *
from h5 import HDFArchive
from triqs.gf import *
import sys, triqs.version as triqs_version
from triqs_dft_tools.sumk_dft import *
@ -18,14 +17,11 @@ warnings.filterwarnings("ignore", category=FutureWarning)
filename = 'vasp'
beta = 5.0
mesh = MeshImFreq(beta=beta, S='Fermion', n_iw=1000)
mesh = MeshImFreq(beta=beta, S='Fermion', n_iw=500)
SK = SumkDFT(hdf_file=filename+'.h5', use_dft_blocks=False, mesh=mesh)
Sigma = SK.block_structure.create_gf(mesh=mesh)
SK.put_Sigma([Sigma])
G = SK.extract_G_loc()
G = SK.extract_G_loc(transform_to_solver_blocks=False, with_Sigma=False)
SK.analyse_block_structure_from_gf(G, threshold=1e-3)
for i_sh in range(len(SK.deg_shells)):
num_block_deg_orbs = len(SK.deg_shells[i_sh])
@ -46,7 +42,7 @@ 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, n_iw=1000)
S = Solver(beta=beta, gf_struct=gf_struct, n_iw=500)
# Construct the Hamiltonian and save it in Hamiltonian_store.txt
H = Operator()
@ -70,13 +66,12 @@ mpi.report('Up Matrix set to:\n%s' % Upmat)
p = {}
p["max_time"] = -1
p["random_name"] = ""
p["random_seed"] = 123 * mpi.rank + 567
p["length_cycle"] = 100
p["n_warmup_cycles"] = 8000
p["n_cycles"] = 200000
p["length_cycle"] = 400
p["n_warmup_cycles"] = 3000
p["n_cycles"] = 20000
p["fit_max_moment"] = 4
p["fit_min_n"] = 30
p["fit_max_n"] = 50
p["fit_min_w"] = 20
p["fit_max_w"] = 30
p["perform_tail_fit"] = True
# Double Counting: 0 FLL, 1 Held, 2 AMF
@ -88,32 +83,33 @@ n_iterations = 10
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_versions')
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.triqs_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.triqs_dft_tools_hash
ar['DMFT_input']['sumk_block_structure'] = SK.block_structure
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()
with HDFArchive(filename+'.h5', 'a') as ar:
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_versions')
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.triqs_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.triqs_dft_tools_hash
ar['DMFT_input']['sumk_block_structure'] = SK.block_structure
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)
@ -146,10 +142,11 @@ 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
with HDFArchive(filename+'.h5', 'a') as ar:
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)
@ -165,13 +162,14 @@ for it in range(iteration_offset, iteration_offset + n_iterations):
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
with HDFArchive(filename+'.h5', 'a') as ar:
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
mpi.report('-------------')
if mpi.is_master_node():
del ar

104
doc/tutorials/nio_csc_vasp/nio_csc.py
View File

@ -1,7 +1,6 @@
from itertools import *
import numpy as np
import triqs.utility.mpi as mpi
from h5 import *
from h5 import HDFArchive
from triqs.gf import *
import sys, triqs.version as triqs_version
from triqs_dft_tools.sumk_dft import *
@ -25,13 +24,11 @@ def dmft_cycle():
Converter.convert_dft_input()
beta = 5.0
mesh = MeshImFreq(beta=beta, S='Fermion', n_iw=1000)
mesh = MeshImFreq(beta=beta, S='Fermion', n_iw=500)
SK = SumkDFT(hdf_file=filename+'.h5', use_dft_blocks=False, mesh=mesh)
Sigma = SK.block_structure.create_gf(mesh=mesh)
SK.put_Sigma([Sigma])
G = SK.extract_G_loc()
G = SK.extract_G_loc(transform_to_solver_blocks=False, with_Sigma=False)
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])
@ -52,7 +49,7 @@ def dmft_cycle():
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, n_iw=1000)
S = Solver(beta=beta, gf_struct=gf_struct, n_iw=500)
# Construct the Hamiltonian and save it in Hamiltonian_store.txt
H = Operator()
@ -75,13 +72,12 @@ def dmft_cycle():
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["length_cycle"] = 400
p["n_warmup_cycles"] = 3000
p["n_cycles"] = 20000
p["fit_max_moment"] = 4
p["fit_min_n"] = 30
p["fit_max_n"] = 50
p["fit_min_w"] = 20
p["fit_max_w"] = 30
p["perform_tail_fit"] = True
# Double Counting: 0 FLL, 1 Held, 2 AMF
@ -93,32 +89,32 @@ def dmft_cycle():
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_versions')
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.triqs_cthyb_hash
ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
ar['DMFT_input']['code_versions']["dft_tools_git"] = dft_tools_version.triqs_dft_tools_hash
ar['DMFT_input']['sumk_block_structure'] = SK.block_structure
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()
with HDFArchive(filename+'.h5', 'a') as ar:
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_versions')
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.triqs_cthyb_hash
ar['DMFT_input']['code_versions']["dft_tools_version"] = dft_tools_version.version
ar['DMFT_input']['code_versions']["dft_tools_git"] = dft_tools_version.triqs_dft_tools_hash
ar['DMFT_input']['sumk_block_structure'] = SK.block_structure
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)
@ -151,10 +147,11 @@ def dmft_cycle():
# 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
with HDFArchive(filename+'.h5', 'a') as ar:
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,
@ -171,13 +168,14 @@ def dmft_cycle():
mpi.report('Total charge of Gloc : %.6f' % S.G_iw.total_density().real)
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
with HDFArchive(filename+'.h5', 'a') as ar:
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...')
@ -196,10 +194,8 @@ def dmft_cycle():
SK.calc_dc(dm, U_interact=U, J_hund=J, orb=0, use_dc_formula=DC_type, use_dc_value=DC_value)
if mpi.is_master_node():
ar['DMFT_results']['Iterations']['corr_energy_it'+str(it)] = correnerg
ar['DMFT_results']['Iterations']['dc_energy_it'+str(it)] = SK.dc_energ[0]
if mpi.is_master_node():
del ar
with HDFArchive(filename+'.h5', 'a') as ar:
ar['DMFT_results']['Iterations']['corr_energy_it'+str(it)] = correnerg
ar['DMFT_results']['Iterations']['dc_energy_it'+str(it)] = SK.dc_energ[0]
return correnerg, SK

22
doc/tutorials/nio_csc_vasp/nio_csc.rst
View File

@ -5,6 +5,8 @@ DFT and projections
We will perform DFT+DMFT calculations 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-self consistency).
**Note: This example works with VASP version 6.5.0 or newer with hdf5 support enabled**
VASP setup
-------------------------------
We start by running a simple VASP calculation to converge the charge density initially.
@ -17,8 +19,7 @@ for our many-body calculation.
.. literalinclude:: INCAR
`LORBIT = 14` takes care of optimizing the projectors in the energy window defined
by `EMIN` and `EMAX`. We switch off all symmetries with `ISYM=-1` since symmetries
are not implemented in the later DMFT scripts. Finally, we select the relevant orbitals
by `EMIN` and `EMAX`. Finally, we select the relevant orbitals
for atom 1 (Ni, d-orbitals) and 2 (O, p-orbitals) by the two `LOCPROJ` lines.
For details refer to the VASP wiki on the `LOCPROJ <https://cms.mpi.univi
e.ac.at/wiki/index.php/LOCPROJ>`_ flag. The projectors are stored in the file `LOCPROJ`.
@ -26,8 +27,9 @@ e.ac.at/wiki/index.php/LOCPROJ>`_ flag. The projectors are stored in the file `L
PLOVASP
------------------------------
Next, we post-process the projectors, which VASP stored in the file `LOCPROJ`.
We do this by invoking :program:`plovasp plo.cfg` which is configured by an input file, e.g., named :ref:`plo.cfg`.
Next, we post-process the projectors, which VASP stored in the file `vaspout.h5/results/locproj`.
You can also take a look at the text file `LOCPROJ` which holds the equivalent information.
By invoking :program:`plovasp plo.cfg` we run the converter, which is configured by an input file, e.g., named :ref:`plo.cfg`.
.. literalinclude:: plo.cfg
@ -43,7 +45,7 @@ optional but later used in the post-processing.
Converting to hdf5 file
-------------------------------
We run the whole conversion to a dft_tools readable h5 archive by running the converter script provided :program:`python converter.py`
We run the whole conversion to a dft_tools readable h5 archive by running the converter script provided :program:`python converter.py` . This actually includes the `plovasp` set in the first line:
.. literalinclude:: converter.py
@ -80,24 +82,26 @@ 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::
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 (see also `INCAR.CSC`)::
ICHARG = 5
NELM = 1000
NELMIN = 1000
NELMDL = -1
IMIX=1
BMIX=0.5
AMIX=0.02
LSYNCH5=True
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`. We change the mixing here to stabilize the updating, which can be problem for charge ordered systems. Vasp 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.
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` or `vaspgamma.h5` if VASP is compiled with hdf5 support. We change the mixing here to stabilize the updating, which can be problem for charge ordered systems. Vasp 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 straightforwardly (for a working example see :ref:`nio_csc.py`). The most important new lines are::
We take the respective converged DFT and DMFT calculations from before as a starting point. I.e., we copy the `WAVECAR` and `nio.h5` together with the other :program:`VASP` input files (copy INCAR.CSC here) 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 straightforwardly (for a working example see :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-Galitski formula. We also slightly increase the tolerance for the detection of blocks since the DFT calculation now includes some QMC noise.
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` :file:`vaspgamma.h5`. The correlation energy is calculated via Migdal-Galitski formula. We also slightly increase the tolerance for the detection of blocks since the starting point now includes some QMC noise.
To help convergence, we keep the density (i.e., the GAMMA file) fixed for a few DFT iterations. We do so since VASP uses an iterative diagonalization. Within these steps we still need to update the projectors and recalculate the GAMMA file to not shuffle eigenvalues around by accident.