3
0
mirror of https://github.com/triqs/dft_tools synced 2024-06-29 08:24:54 +02:00

work in NiO tutorial

This commit is contained in:
Malte Schüler 2019-07-08 09:55:22 +02:00
parent 1eecf80b66
commit a882ffa575
22 changed files with 445 additions and 15 deletions

View File

@ -18,6 +18,6 @@ add_custom_target(doc_sphinx ALL DEPENDS ${sphinx_top} ${CMAKE_CURRENT_BINARY_DI
# ---------------------------------
install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/html/ COMPONENT documentation DESTINATION share/doc/triqs_dft_tools
FILES_MATCHING
REGEX "\\.(html|pdf|png|gif|jpg|js|xsl|css|py|txt|inv|bib)$"
REGEX "\\.(html|pdf|png|gif|jpg|js|xsl|css|py|txt|inv|bib|cfg)$"
PATTERN "_*"
)

View File

@ -5,9 +5,7 @@ Interface with VASP
===================
.. warning::
The VASP interface is in the alpha-version and the VASP part of it is not
yet publicly released. The documentation may, thus, be subject to changes
before the final release.
The VASP interface is in the alpha-version. The documentation may, thus, be subject to changes before the final release.
*Limitations of the alpha-version:*
@ -26,7 +24,8 @@ in the :ref:`PLOVasp User's Guide <plovasp>`. Here, a quick-start guide is prese
The VASP interface relies on new options introduced since version
5.4.x. In particular, a new INCAR-option `LOCPROJ`
and new `LORBIT` modes 13 and 14 have been added.
and new `LORBIT` modes 13 and 14 have been added as well as the new ICHARG
mode 5 for charge self-consistent calculations
Option `LOCPROJ` selects a set of localized projectors that will
be written to file `LOCPROJ` after a successful VASP run.
@ -41,7 +40,7 @@ with the indices corresponding to the site position in the POSCAR file;
`<projector type>` chooses a particular type of the local basis function.
The recommended projector type is `Pr 2`. The formalism for this type
of projectors is presented in
`M. Schüler et al. 2018 J. Phys.: Condens. Matter 30 475901 <https://doi.org/10.1088/1361-648X/aae80a>`_.
`M. Schüler et al. 2018 J. Phys.: Condens. Matter 30 475901 <https://doi.org/10.1088/1361-648X/aae80a>`_. For details on `LOCPROJ` also have a look in the `VASP wiki <https://cms.mpi.univie.ac.at/wiki/index.php/LOCPROJ>`_
The allowed labels of the local states defined in terms of cubic
harmonics are:

View File

@ -55,21 +55,25 @@ A PLOVasp input file can contain three types of sections:
Section [General]
"""""""""""""""""
The entire section is optional and it contains three parameters:
The entire section is optional and it contains four parameters:
* **BASENAME** (string): provides a base name for output files.
Default filenames are :file:`vasp.*`.
* **DOSMESH** ([float float] integer): if this parameter is given,
the projected density of states for each projected orbital will be
evaluated and stored to files :file:`pdos_<n>.dat`, where `n` is the
orbital index. The energy
mesh is defined by three numbers: `EMIN` `EMAX` `NPOINTS`. The first two
evaluated and stored to files :file:`pdos_<s>_<n>.dat`, where `s` is the
shell index and `n` the ion index. The energy mesh is defined by three
numbers: `EMIN` `EMAX` `NPOINTS`. The first two
can be omitted in which case they are taken to be equal to the projector
energy window. **Important note**: at the moment this option works
only if the tetrahedron integration method (`ISMEAR = -4` or `-5`)
is used in VASP to produce `LOCPROJ`.
* **EFERMI** (float): provides the Fermi level. This value overrides
the one extracted from VASP output files.
* **HK** (True/False): If True, the projectors are applied the the Kohn-Sham
eigenvalues which results in a Hamitlonian H(k) in orbital basis. The H(k)
is written for each group to a file :file:`Basename.hk<Ng>`. It is recommended
to also set `COMPLEMENT = True` (see below). Default is False.
There are no required parameters in this section.
@ -94,6 +98,8 @@ given by `LSHELL` must be present in the LOCPROJ file.
There are additional optional parameters that allow one to transform
the local states:
* **CORR** (True/False): Determines if shell is correlated or not. At least one
shell has to be correlated. Default is True.
* **TRANSFORM** (matrix): local transformation matrix applied to all states
in the projector shell. The matrix is defined by a (multiline) block
of floats, with each line corresponding to a row. The number of columns
@ -131,6 +137,13 @@ Optional group parameters:
condition will be enforced on each site separately but the Wannier functions
on different sites will not be orthogonal. If `NORMION = False`, the Wannier functions
on different sites included in the group will be orthogonal to each other.
* **BANDS** (int int): the energy window specified by two ints: band index of
lowest band and band index of highest band. Using this overrides the selection
in `EWINDOW`.
* **COMPLEMENT** (True/False). If True, the orthogonal complement is calculated
resulting in unitary (quadratic) projectors, i.e., the same number of orbitals
as bands. It is required to have an equal number of bands in the energy window
at each k-point. Default is False.
.. _transformation_file:

View File

@ -23,3 +23,12 @@ Full charge self consistency with Wien2k: :math:`\gamma`-Ce
tutorials/ce-gamma-fscs_wien2k
A full example with VASP: NiO
-----------------------------------------------------------
.. toctree::
:maxdepth: 2
tutorials/nio_csc

View File

@ -0,0 +1,18 @@
System = NiO
ISMEAR = -5
# the energy window to optimize projector channels
EMIN = -3
EMAX = 7
NBANDS = 16
LMAXMIX = 6
# switch off all symmetries
ISYM = -1
# project to Ni d and O p states
LORBIT = 14
LOCPROJ = 1 : d : Pr 1
LOCPROJ = 2 : p : Pr 1

View File

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

View File

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

View File

@ -0,0 +1,7 @@
.. _KPOINTS:
KPOINTS
-----------
.. literalinclude:: KPOINTS
:language: bash

View File

@ -0,0 +1,90 @@
from itertools import *
import numpy as np
import pytriqs.utility.mpi as mpi
from pytriqs.archive import *
from pytriqs.gf import *
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 warnings
warnings.filterwarnings("ignore", category=FutureWarning)
filename = 'vasp'
SK = SumkDFT(hdf_file = filename+'.h5', use_dft_blocks = False)
beta = 5.0
# We analyze the block structure of the Hamiltonian
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-3)
# Setup CTQMC Solver
n_orb = SK.corr_shells[0]['dim']
spin_names = ['up','down']
orb_names = [i for i in range(0,n_orb)]
# Print some information on the master node
iteration_offset = 0
Sigma_iw = 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 'iteration_count' in ar['DMFT_results']:
iteration_offset = ar['DMFT_results']['iteration_count']+1
print('offset',iteration_offset)
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
iteration_offset = mpi.bcast(iteration_offset)
Sigma_iw = mpi.bcast(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)
SK.put_Sigma(Sigma_imp = [Sigma_iw])
ikarray = numpy.array(range(SK.n_k))
# set up the orbitally resolved local lattice greens function:
n_orbs = SK.proj_mat_csc.shape[2]
spn = SK.spin_block_names[SK.SO]
mesh = Sigma_iw.mesh
block_structure = [range(n_orbs) for sp in spn]
gf_struct = [(spn[isp], block_structure[isp])
for isp in range(SK.n_spin_blocks[SK.SO])]
block_ind_list = [block for block, inner in gf_struct]
glist = lambda: [GfImFreq(indices=inner, mesh=mesh)
for block, inner in gf_struct]
G_latt_orb = BlockGf(name_list=block_ind_list,
block_list=glist(), make_copies=False)
G_latt_orb.zero()
for ik in mpi.slice_array(ikarray):
G_latt_KS = SK.lattice_gf(ik=ik, beta=beta)
G_latt_KS *= SK.bz_weights[ik]
for bname, gf in G_latt_orb:
add_g_ik = gf.copy()
add_g_ik.zero()
add_g_ik << SK.downfold(ik, 0, bname, G_latt_KS[bname], gf, shells='csc', ir=None)
gf << gf + add_g_ik
G_latt_orb << mpi.all_reduce(
mpi.world, G_latt_orb, lambda x, y: x + y)
mpi.barrier()
if mpi.is_master_node():
ar['DMFT_results']['Iterations']['G_latt_orb_it'+str(iteration_offset-1)] = G_latt_orb
if mpi.is_master_node(): del ar

View File

@ -0,0 +1,10 @@
NiO
4.17 #taken from 9x9x9 with sigma=0.2 ismear=2
+0.0000000000 +0.5000000000 +0.5000000000
+0.5000000000 +0.0000000000 +0.5000000000
+0.5000000000 +0.5000000000 +0.0000000000
Ni O
1 1
Direct
+0.0000000000 +0.0000000000 +0.0000000000
+0.5000000000 +0.5000000000 +0.5000000000

View File

@ -0,0 +1,7 @@
.. _POSCAR:
POSCAR
-----------
.. literalinclude:: POSCAR
:language: bash

View File

@ -0,0 +1,3 @@
from triqs_dft_tools.converters.vasp_converter import *
Converter = VaspConverter(filename = 'nio')
Converter.convert_dft_input()

View File

@ -0,0 +1,10 @@
.. _converter.py:
converter.py
-------------------
Download :download:`converter.py <./converter.py>`.
.. literalinclude:: converter.py
:language: python

View File

@ -0,0 +1,149 @@
from itertools import *
import numpy as np
import pytriqs.utility.mpi as mpi
from pytriqs.archive import *
from pytriqs.gf import *
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 *
filename = 'vasp'
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-3)
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)]
orb_hyb = False
#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"] = 20000
p["n_cycles"] = 200000
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
# Prepare hdf file and and check for previous iterations
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 '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
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.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)
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_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)
# 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(): del ar

View File

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

View File

@ -0,0 +1,27 @@
[General]
BASENAME = nio
DOSMESH = -21 55 400
[Group 1]
SHELLS = 1 2
EWINDOW = -9 2
NORMION = False
NORMALIZE = True
[Shell 1] # Ni d shell
LSHELL = 2
IONS = 1
CORR = True
TRANSFORM = 1.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 0.0
0.0 0.0 1.0 0.0 0.0
0.0 0.0 0.0 1.0 0.0
0.0 0.0 0.0 0.0 1.0
[Shell 2] # O p shell
LSHELL = 1
IONS = 2
CORR = False
TRANSFORM = 1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0

View File

@ -0,0 +1,9 @@
.. _plo.cfg:
plo.cfg
-----------
Download :download:`plo.cfg<./plo.cfg>`.
.. literalinclude:: plo.cfg
:language: bash

59
doc/tutorials/nio_csc.rst Normal file
View File

@ -0,0 +1,59 @@
.. _nio_csc:
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.
VASP setup
-------------------------------
We start by running a simple VASP calculation to converge the charge density initially.
Use the :ref:`INCAR`, :ref:`POSCAR`, and :ref:`KPOINTS` provided and use your
own :file:`POTCAR` file.
Let us take a look in the :file:`INCAR`, where we have to specify local orbitals as basis
for our many-body calculation.
.. literalinclude:: images_scripts/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
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`.
plovasp
------------------------------
Next, we postprocess 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`.
.. literalinclude:: images_scripts/plo.cfg
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
and the O p states. Only the Ni shell is treated as correlated (`CORR = True`), i.e.,
is supplemented with a Coulomb interaction later in the DMFT calculation.
Converting to hdf5 file
-------------------------------
We gather the output generated by :program:`plovasp` into a hdf5 archive which :program:`dft_tools` is able to read. We do this by running :program:`python converter.py` on the script :ref:`converter.py`:
.. literalinclude:: images_scripts/converter.py
Now we are all set to perform a dmft calculation.
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.
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`

View File

@ -381,7 +381,7 @@ class VaspConverter(ConverterTools):
'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']
if self.proj_or_hk == '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]

View File

@ -97,13 +97,13 @@ class SumkDFT(object):
# Read input from HDF:
things_to_read = ['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',
'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'proj_mat_csc', 'bz_weights', 'hopping',
'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr','proj_or_hk']
self.subgroup_present, self.value_read = self.read_input_from_hdf(
subgrp=self.dft_data, things_to_read=things_to_read)
if self.proj_or_hk == 'hk':
self.subgroup_present, self.value_read = self.read_input_from_hdf(
subgrp=self.dft_data, things_to_read=['proj_mat_csc'])
# if self.proj_or_hk == 'hk':
# self.subgroup_present, self.value_read = self.read_input_from_hdf(
# subgrp=self.dft_data, things_to_read=['proj_mat_csc'])
if self.symm_op:
self.symmcorr = Symmetry(hdf_file, subgroup=self.symmcorr_data)

Binary file not shown.