HubbardHamiltonian

class hubbard.HubbardHamiltonian(TBHam, n=0, U=None, Uij=None, q=(0.0, 0.0), nkpt=[1, 1, 1], kT=1e-05, units='eV')

Bases: object

A class to create a Self Consistent field (SCF) object related to the mean-field Hubbard (MFH) model

The hubbard.HubbardHamiltonian class opens the possibility to include electron correlations in the tight-binding Hamiltonian by solving self-consistently the mean-field Hubbard Hamiltonian

It enables the convergence of several tight-binding described systems towards a user-defined tolerance criterion

It takes an input tight-binding Hamiltonian and updates the corresponding matrix elements according to the MFH model

Parameters

• TBHam (sisl.physics.Hamiltonian) – An unpolarized or spin-polarized tight-binding Hamiltonian

• n (numpy.ndarray, optional) – initial spin-densities vectors. The shape of n must be (spin deg. of freedom, no. of sites)

• U (float or np.ndarray, optional) – intra-orbital Coulomb repulsion parameter. If an array is passed, the length has to be equal to the number of orbitals

• Uij (np.ndarray, optional) – inter-orbital Coulomb repulsion parameter. The size of the matrix has to be equal to the number of orbitals

• q (array_like, optional) – One or two values specifying the total charge associated to each spin component. The array should contain as many values as the dimension of the problem. I.e., if the Hamiltonian is unpolarized q should be a single value

• nkpt (array_like or sisl.physics.BrillouinZone, optional) – Number of k-points along (a1, a2, a3) for Monkhorst-Pack BZ sampling

• kT (float, optional) – Temperature of the system in units of the Boltzmann constant

• units (str, optional) – (energy) units of U, kT and the TB Hamiltonian parameters (namely the hopping term and the onsite energies). If one uses another units then it should be specified here. electron volts are used by default

Note

The implementation to solve the extended Hubbard model (with inter-atomic interactions) has not been tested

Attributes

shape

The shape of the Hamiltonian matrix

Methods

DOS(egrid[, eta, spin, distribution])	Obtains the density of states (DOS) of the system with a distribution function
PDOS(egrid[, eta, spin, distribution])	Obtains the projected density of states (PDOS) of the system with a distribution function
__init__(TBHam[, n, U, Uij, q, nkpt, kT, units])	Initialize HubbardHamiltonian
calc_orbital_charge_overlaps([k, spin])	Obtain orbital (eigenstate) charge overlaps as \(\int dr |\psi_{\sigma\alpha}|^{4}\)
converge(calc_n_method[, tol, mixer, steps, ...])	Iterate Hamiltonian towards a specified tolerance criterion
copy()	Return a copy of the HubbardHamiltonian object
eigenstate(k[, spin])	Solve the eigenvalue problem at k and return it as a sisl.physics.electron.EigenstateElectron object containing all eigenstates
eigh([k, eigvals_only, spin])	Diagonalize Hamiltonian using the eigh routine
fermi_level([q, distribution])	Find the fermi level for a certain charge q at a certain kT
find_midgap()	Find the midgap for the system taking into account the up and dn different spectrums
get_Zak_phase([func, axis, nk, sub, eigvals])	Computes the intercell Zak phase for 1D (periodic) systems using sisl.physics.electron.berry_phase
get_bond_order([format, midgap])	Compute Huckel bond order
get_hash()
iterate(calc_n_method[, q, mixer])	Common method to iterate in a SCF loop that corresponds to the mean-field Hubbard approximation
normalize_charge()	Ensure the total up/down charge in pi-network equals Nup/Ndn
polarize_sublattices()	Quick way to polarize the lattices without checking that consequtive atoms actually belong to different sublattices
random_density()	Initialize spin polarization with random density
read_density(fn[, mode, group])	Read density from binary file
remove(atoms[, q])	Remove a subset of this sparse matrix by only retaining the atoms corresponding to atom
repeat(reps, axis)	Repeat the HubbardHamiltonian object along a specified axis to obtain a larger one
set_kmesh([nkpt])	Set the k-mesh for the HubbardHamiltonian
set_polarization(up[, dn])	Maximize spin polarization on specific atomic sites Optionally, sites with down-polarization can be specified
shift(E)	Shift the electronic structure by a constant energy (in-place operation)
spin_contamination([ret_exact])	Obtains the spin contamination for the MFH Hamiltonian following `Ref.
sub(atoms[, q])	Return a new HubbardHamiltonian object of a subset of selected atoms
tile(reps, axis)	Tile the HubbardHamiltonian object along a specified axis to obtain a larger one
update_hamiltonian()	Update spin Hamiltonian according to the extended Hubbard model, see for instance PRL 106, 236805 (2011) It updates the diagonal elements for each spin Hamiltonian following the mean field approximation:
write_density(fn[, mode, group])	Write density in a binary file
write_initspin(fn[, ext_geom, spinfix, ...])	Write spin polarization to SIESTA fdf-block This function only makes sense for spin-polarized calculations

DOS(egrid, eta=0.001, spin=[0, 1], distribution='Lorentzian')

Obtains the density of states (DOS) of the system with a distribution function

Parameters

• egrid (array_like) – Energy grid at which the DOS will be calculated.

• eta (float, optional) – Smearing parameter

• spin (int, optional) – If spin=0(1) it calculates the DOS for up (down) electrons in the system. If spin is not specified it returns DOS_up + DOS_dn.

• distribution (str or sisl.physics.distribution, optional) – distribution for the convolution, defaults to Lorentzian

See also

sisl.physics.get_distribution

sisl method to create distribution function

sisl.physics.electron.DOS

sisl method to obtain DOS

Returns

DOS – density of states at the given energies for the selected spin

Return type

numpy.ndarray

PDOS(egrid, eta=0.001, spin=[0, 1], distribution='Lorentzian')

Obtains the projected density of states (PDOS) of the system with a distribution function

Parameters

• egrid (array_like) – Energy grid at which the DOS will be calculated.

• eta (float, optional) – Smearing parameter

• spin (int, optional) – If spin=0(1) it calculates the DOS for up (down) electrons in the system. If spin is not specified it returns DOS_up + DOS_dn.

• distribution (str or sisl.distribution, optional) – distribution for the convolution, defaults to Lorentzian

See also

sisl.get_distribution

sisl method to create distribution function

sisl.physics.electron.PDOS

sisl method to obtain PDOS

Returns

PDOS – projected density of states at the given energies for the selected spin

Return type

numpy.ndarray

calc_orbital_charge_overlaps(k=[0, 0, 0], spin=0)

Obtain orbital (eigenstate) charge overlaps as \(\int dr |\psi_{\sigma\alpha}|^{4}\)

Where \(\sigma\) is the spin index and \(\alpha\) is the eigenstate index

Parameters

• k (array_like, optional) – k-point at which the eigenstate is going to be obtained

• spin (int, optional) – spin index

Returns

• ev (numpy.ndarray) – eigenvalues

• L (numpy.ndarray) – orbital charge overlaps

converge(calc_n_method, tol=1e-06, mixer=None, steps=100, max_iter=None, fn=None, print_info=False, func_args={})

Iterate Hamiltonian towards a specified tolerance criterion

This method calls iterate as many times as it needs until it reaches the specified tolerance

Parameters

• calc_n_method (callable) – method to obtain the spin-densities it must return the corresponding spin-densities (n) and the total energy (Etot)

• tol (float, optional) – tolerance criterion

• mixer (Mixer, optional) – sisl.mixing.Mixer instance, defaults to sisl.mixing.DIISMixer(0.7, history=7)

• steps (int, optional) – the code will print some relevant information (if print_info is set to True) about the convergence process when the number of completed iterations reaches a multiple of the specified steps. It also will store the densities in a binary file if fn is passed

• max_iter (int, optional) – maximum number of iterations before stopping

• fn (str, optional) – optionally, one can save the spin-densities during the calculation (when the number of completed iterations reaches the specified steps), by giving the name of the full name of the binary file

• func_args (dictionary, optional) – function arguments to pass to calc_n_method

See also

iterate

hubbard.calc_n

method to obtain n and Etot for tight-binding Hamiltonians with finite or periodic boundary conditions at a certain kT

hubbard.NEGF

class that contains the routines to obtain n and Etot for tight-binding Hamiltonians with open boundary conditions

sisl.mixing.AdaptiveDIISMixer

for adaptative DIIS (Pulay) mixing scheme

sisl.mixing.LinearMixer

for linear mixing scheme

Returns

difference between the ith and the (i-1)th iteration densities

Return type

dn

copy()

Return a copy of the HubbardHamiltonian object

eigenstate(k, spin=0)

Solve the eigenvalue problem at k and return it as a sisl.physics.electron.EigenstateElectron object containing all eigenstates

Parameters

• k (array_like) – k-point at which the eigenstate is going to be obtained

• spin (int, optional) – for spin=0(1) it solves the eigenvalue problem for the spin up (down) Hamiltonian

Returns

object

Return type

sisl.physics.electron.EigenstateElectron object

eigh(k=[0, 0, 0], eigvals_only=True, spin=0)

Diagonalize Hamiltonian using the eigh routine

Parameters

• k (array_like, optional) – k-point at which the eigenvalue want to be obtained

• eigvals_only (bool, optional) – if True only eigenvalues are returned, otherwise it also returns the eigenvectors

• spin (int, optional) – for spin = 0 (1) it solves the eigenvalue problem for the spin up (down) Hamiltonian

See also

sisl.physics.SparseOrbitalBZ.eigh

sisl class

Returns

eigenvalues

Return type

numpy.ndarray

fermi_level(q=[None, None], distribution='fermi_dirac')

Find the fermi level for a certain charge q at a certain kT

Parameters

• q (array_like, optional) – charge per spin channel. First index for spin up, second index for dn If the Hamiltonian is unpolarized q should have only one component otherwise it will take the first one

• distribution (str or sisl.distribution, optional) – distribution function

See also

sisl.physics.Hamiltonian.fermi_level

sisl class function

Returns

Ef – Fermi-level for each spin channel

Return type

numpy.array

find_midgap()

Find the midgap for the system taking into account the up and dn different spectrums

This method makes sense for insulators (where there is a bandgap)

Returns

midgap

Return type

float

get_Zak_phase(func=None, axis=0, nk=51, sub='filled', eigvals=False)

Computes the intercell Zak phase for 1D (periodic) systems using sisl.physics.electron.berry_phase

Parameters

• func (callable, optional) – function that creates a list of parametrized k-points to generate a new sisl.physics.BrillouinZone object parametrized in N separations

• axis (int, optional) – index for the periodic direction

• nk (int, optional) – number of k-points generated using the parameterization

• sub (int, optional) – number of bands that will be summed to obtain the Zak phase

• eigvals (bool, optional) – return the eigenvalues of the product of the overlap matrices, see sisl.physics.electron.berry_phase

Notes

If no func is passed it assumes the periodicity along the x-axis If no sub is passed it sums up to the last occuppied band (included)

Returns

Zak – Intercell Zak phase for the 1D system

Return type

float

get_bond_order(format='csr', midgap=0.0)

Compute Huckel bond order

Parameters

• format ({'csr', 'array', 'dense', 'coo', ...}) – the returned format of the matrix, defaulting to the scipy.sparse.csr_matrix, however if one always requires operations on dense matrices, one can always return in numpy.ndarray (‘array’) or numpy.matrix (‘dense’).

• midgap (float, optional) – energy value that separates filled states (lower energy) from empty states (higher energy)

Return type

the Huckel bond-order matrix object

get_hash()

iterate(calc_n_method, q=None, mixer=None, **kwargs)

Common method to iterate in a SCF loop that corresponds to the mean-field Hubbard approximation

The only thing that may change is the way in which the spin-densities (n) and total energy (Etot) are obtained where one needs to use the correct calc_n_method for the particular system.

Parameters

• calc_n_method (callable) – method to obtain the spin-densities it must return the corresponding spin-densities (n) and the total energy (Etot)

• q (array_like, optional) – total charge separated in spin-channels, q=[q_up, q_dn]

• mixer (Mixer, optional) – sisl.mixing.Mixer instance for the SCF loop, defaults to linear mixing with weight 0.7

See also

update_hamiltonian

hubbard.calc_n

method to obtain n and Etot for tight-binding Hamiltonians with finite or periodic boundary conditions at a certain kT

hubbard.NEGF.calc_n_open

method to obtain n and Etot for tight-binding Hamiltonians with open boundary conditions

sisl.mixing.AdaptiveDIISMixer

for adaptative DIIS (Pulay) mixing scheme

sisl.mixing.LinearMixer

for linear mixing scheme

Returns

dn – difference between the ith and the (i-1)th iteration densities

Return type

array_like

normalize_charge()

Ensure the total up/down charge in pi-network equals Nup/Ndn

polarize_sublattices()

Quick way to polarize the lattices without checking that consequtive atoms actually belong to different sublattices

random_density()

Initialize spin polarization with random density

read_density(fn, mode='r', group=None)

Read density from binary file

Parameters

• fn (str) – name of the file that is going to read from

• mode (str, optional) – mode in which the file is opened

• group (str, optional) – netCDF4 group

remove(atoms, q=(0, 0))

Remove a subset of this sparse matrix by only retaining the atoms corresponding to atom

Parameters

• atoms (array_like of int) – atomic ndices of removed atoms

• q (array_like, optional) – Two values specifying up, down electron occupations for the remaining subset of atoms after removal

See also

sisl.physics.Hamiltonian.remove, sisl.Geometry.remove

repeat(reps, axis)

Repeat the HubbardHamiltonian object along a specified axis to obtain a larger one

Parameters

• reps (int) – number of repetitions

• axis (int) – direction of tiling, 0, 1, 2 according to the cell-direction

See also

sisl.Geometry.repeat

sisl class method

Return type

A new larger hubbard.HubbardHamiltonian object

Note

We need to update this function so Uij also gets repeated like the geometry

set_kmesh(nkpt=[1, 1, 1])

Set the k-mesh for the HubbardHamiltonian

Parameters

nkpt (array_like or sisl.physics.BrillouinZone, optional) – k-mesh to be associated with the hubbard.HubbardHamiltonian instance

set_polarization(up, dn=())

Maximize spin polarization on specific atomic sites Optionally, sites with down-polarization can be specified

Parameters

• up (array_like) – atomic sites where the spin-up density is going to be maximized

• dn (array_like, optional) – atomic sites where the spin-down density is going to be maximized

property shape

The shape of the Hamiltonian matrix

shift(E)

Shift the electronic structure by a constant energy (in-place operation)

Parameters

E (float or (2,)) – the energy (in eV) to shift the electronic structure, if two values are passed the two spin-components get shifted individually

See also

sisl.physics.Hamiltonian.shift

spin_contamination(ret_exact=False)

Obtains the spin contamination for the MFH Hamiltonian following Ref. Chemical Physics Letters. 183 (5): 423–431.

\[\langle S^{2} \rangle_{MFH} = \langle S^{2} \rangle_{exact} + N_{\beta} - \sum_{ij}^{occ} |\langle \psi^{\alpha}_{i}|\psi^{\beta}_{j}\rangle|^{2}\]

Where the exact spin squared expectation value is obtained as

\[\langle S^{2} \rangle_{exact}=(\frac{ N_{\alpha}-N_{\beta}}{2})(\frac{N_{\alpha}-N_{\beta} }{2} + 1)\]

Notes

The current implementation works for non-periodic systems only.

Parameters

ret_exact (bool, optional) – If true this method will return also the exact spin squared expectation value

See also

sisl.physics.electron.spin_squared

sisl class function

Returns

• S_MFH (float) – expectation value for the MFH Hamiltonian

• S (float) – exact expectation value

sub(atoms, q=(0, 0))

Return a new HubbardHamiltonian object of a subset of selected atoms

Parameters

• atoms (int or array_like) – indices/boolean of all atoms to be kept

• q (array_like, optional) – Two values specifying up, down electron occupations for the subset of atoms

tile(reps, axis)

Tile the HubbardHamiltonian object along a specified axis to obtain a larger one

Parameters

• reps (int) – number of tiles (repetitions)

• axis (int) – direction of tiling, 0, 1, 2 according to the cell-direction

See also

sisl.Geometry.tile

sisl class method

Return type

A new larger hubbard.HubbardHamiltonian object

Note

We need to update this function so Uij also gets tiled like the geometry

update_hamiltonian()

Update spin Hamiltonian according to the extended Hubbard model, see for instance PRL 106, 236805 (2011) It updates the diagonal elements for each spin Hamiltonian following the mean field approximation:

\[\begin{split}H &= -\sum_{ij\sigma}t_{ij\sigma}c^{\dagger}_{i\sigma}c_{j\sigma} + \sum_{i}U_in_{i\uparrow}n_{i\downarrow} + \frac{1}{2}\sum_{i\neq j\sigma\sigma^\prime}U_{ij}n_{i\sigma}n_{j\sigma^\prime} \approx\\ & -\sum_{ij\sigma}t_{ij\sigma}c^{\dagger}_{i\sigma}c_{j\sigma} + \sum_{i\sigma} U_i \left\langle n_{i\sigma}\right\rangle n_{i\bar{\sigma}} + \frac{1}{2}\sum_{i\neq j\sigma}\left(U_{ij} + U_{ji}\right)\left(\langle n_{i\uparrow}\rangle + \langle n_{i\downarrow}\rangle\right)n_{j\sigma} + C\end{split}\]

The constant term \(C = -\sum_i U_i \langle n_{i\uparrow}\rangle\langle n_{i\downarrow}\rangle - \frac{1}{2}\sum_{i\neq j}U_{ij}\left(\langle n_{i\uparrow}\rangle+\langle n_{i\downarrow}\rangle\right)\left(\langle n_{j\uparrow}\rangle + \langle n_{j\downarrow}\rangle\right)\) will be added to the Hamiltonian in the iterate method, where the total energy is calculated

write_density(fn, mode='a', group=None)

Write density in a binary file

Parameters

• fn (str) – name of the file in which the densities are going to be stored

• mode (str, optional) – mode in which the file is opened

• group (str, optional) – netCDF4 group

write_initspin(fn, ext_geom=None, spinfix=True, mode='a', eps=0.1)

Write spin polarization to SIESTA fdf-block This function only makes sense for spin-polarized calculations

Parameters

• fn (str) – name of the fdf-file

• ext_geom (sisl.Geometry, optional) – an “external” geometry that contains the sp2-sites included in the simulation

• spinfix (bool, optional) – specifies if the Spin.Fix and Spin.Total lines are written to the fdf

• mode (str, optional) – mode in which the file is opened

• eps (float, optional) – atoms within this distance will be considered equivalent in case ext_geom != geometry (see sisl.Geometry.overlap)