Keldysh Green functions¶
The Keldysh Green function (KGF) code allows for calculations of non-equilibrium transport calculations where electron exchange and correlation effect are threated using many body perturbation theory such as Hartree-Fock, second Born and the GW approximation. It is recommended that you go through the ASE/GPAW electron transport exercice to get familiar with the general transport setup and definitions used in ase and gpaw and the KGF code.
Download and Installation¶
The KGF code is currently being merged into the development version of GPAW and is expected to be part of the GPAW package in the near future. The latest revision can be obtained from svn:
$ svn checkout https://svn.fysik.dtu.dk/projects/KeldyshGF/trunk KeldyshGF
Installation is completed by adding the path of KeldyshGF to the PYTHONPATH environment variable.
Doing a KGF calculation¶
The KGF code can perform finite bias non-equilibrium calculation of a molecular junction using various electron exchange and correlation approximations. It is assumed that interactions are only included in a central region. The KGF code can handle both model Hamiltonians of the the Pariser-Parr-Pople (extended Hubbard) as well as abinitio calculated Hamiltonians.
A KGF calculatoin normally involves the following steps:
- Setting up the non-interacting lead and scattering Hamiltonian.
- Setting up a non-interacting GF
- Setting up various self-energies to handle Hartree, exchange and correlation
- Runnig the calculation
Example: Pariser-Parr-Model Hamiltonian¶
To do an electron transport calculation using a model Hamiltonian the parameters of both the non-interacting part as well as the interacting part of the Hamiltonian need to be explicitly specified. The non-interacting part h_ij describe kinetic energy and electron-electron interaction part in the PPP approximation is on the form V_ij = v_iijj, where the v_ijkl’s are two electron Coulomb integrals. To get started consider a simple four site interacting model. The four the x’s in the figure below represent the sites where electron-electron interactions are included. The o’s (dashes) represents non-interacting sites.:
Left lead Molecule Right Lead --------------------------------- o o o o x| x x x x |x o o o o --------------------------------- 0 1 2 3 4 5
The numbers refers to indix numbers in the Green functions - the Green function will be a 6x6 matrix where the subspace corresponding to the molecule will be the central 4x4 matrix. Leads are treated as simple nearest neighbour tight-binding chains with a principal layer size of one.
The following parameters will be used to simulate a physisorbed molecule:
||-20.0||intra lead hopping|
||-2.4||intra molecule hopping|
where V is the matrix:
V = [[ 0. 7.45 4.54 3.18 2.42 0. ] [ 7.45 11.26 7.45 4.54 3.18 2.42] [ 4.54 7.45 11.26 7.45 4.54 3.18] [ 3.18 4.54 7.45 11.26 7.45 4.54] [ 2.42 3.18 4.54 7.45 11.26 7.45] [ 0. 2.42 3.18 4.54 7.45 0. ]]
In Python code the input parameters can generated like this:
from __future__ import print_function import numpy as np pos = np.asarray([1.450 * a for a in range(6)]) n = len(pos) d2_ii = np.zeros((n, n)) for i1, i2 in np.ndindex(n, n): d2_ii[i1, i2] = ((pos[i1]-pos[i2])**2).sum() d_ii = np.sqrt(d2_ii) U = 11.26 t_m = -2.4 t_lm = -1.0 # Two particle interaction V = U / np.sqrt(1.0 + 0.6117*d2_ii) V[0,-1] = V[-1,0] = 0.0 V[0,0] = V[-1,-1] = 0.0 # One particle part mask_c = d_ii < 1.5 mask_c *= d_ii > 0.0 mask_c = mask_c.astype(int) nbf = len(d_ii) h = t_m * np.ones((nbf, nbf)) * mask_c h[0, 1] = h[1, 0] = t_lm h[-2,-1] = h[-1,-2] = t_lm # H_lead h1 = np.zeros((2,2), complex) t_l = -20.0 h1[0,1] = t_l h1[1,0] = t_l nbf = len(h) # H_scat H = np.zeros((nbf+2, nbf+2), complex) H[0,1] = H[1,0] = t_l H[-2,-1] = H[-1,-2] = t_l H[1:-1, 1:-1] = h # Hartree potential of the ions (Z=1) ion_shift = np.zeros((n, n)) for i in range(n): ion_shift[i, i] += -0.5 * V[i, i] for k in range(n): if k!=i: ion_shift[i, i] += -1.0 * V[i, k] H = H.astype(complex) V = V.astype(complex) ion_shift = ion_shift.astype(complex) if __name__=='__main__': print(np.around(H.real, 2)) print(np.around(V.real ,2)) print(np.around(ion_shift.diagonal().real, 2))
We begin by performing an equilibrium calculation (zero bias). An equilibrium involces setting the relevant Green’s functions and self-energies. All Green’s functions are represented on energy grid which should have a grid spacing fine enough to resolve all spectreal feautres. In practise this accomplished by choosing an energy grid spacing about half the size of the infinitesimal eta appearing in the Green’s functions (which is given a finite value in numerical calculations).
In Python code an equilibrium non-interacting calculatoins followed by a Hartree-Fock calculations and a GW calculation look like this:
from __future__ import print_function import pickle import numpy as np from gpaw.mpi import world from parms import V, H, h1, ion_shift from kgf.GreenFunctions import NonEqNonIntGreenFunction, NonEqIntGreenFunction from kgf.ScatteringHamiltonianMatrix import ScatteringHamiltonianMatrix from kgf.LeadHamiltonianMatrix import LeadHamiltonianMatrix from kgf.CurrentCalculator import CurrentCalculator from kgf.selfenergies import NonEqConstantSelfEnergy from kgf.selfenergies.GW2index import Hartree2index, Fock2index, GW2index hmat = ScatteringHamiltonianMatrix(leftprincipallayer=1, rightprincipallayer=1, hamiltonianmatrix=H) left = LeadHamiltonianMatrix(principallayer=1, hamiltonianmatrix=h1) right = LeadHamiltonianMatrix(principallayer=1, hamiltonianmatrix=h1) energies = np.arange(-96.0, 96.0, 0.02) de = energies - energies assert abs(energies).min() < 1.0e-8 fermi = 0.0 g0 = NonEqNonIntGreenFunction(hmat, left, right, E_Fermi=fermi, energies=energies) g0.SetInfinitesimal(2*de) Ntot = g0.GetTotalParticleNumber() if world.rank==0: print('Total electron number:', Ntot) se_null = NonEqConstantSelfEnergy(g0, np.zeros((6, 6), complex), 'null') se_ion = NonEqConstantSelfEnergy(g0, ion_shift, 'shift') se_hartree = Hartree2index(g0, V, initialhartree=np.zeros((6, 6), complex)) se_fock = Fock2index(g0, V) se_gw = GW2index(g0, V) gf = NonEqIntGreenFunction([se_null, se_ion, se_hartree, se_fock, se_gw]) orbitals = range(1, 5) pulay = (0.03, 0.25, 1) # Non interacting calculation se_ion.SetActive(False) se_hartree.SetActive(False) se_fock.SetActive(False) se_gw.SetActive(False) gf.WriteSpectralInfoToNetCDFFile('nonint.nc', diagonalize=True, orbitals=orbitals, spectral='individual') # HF calculation se_ion.SetActive(True) se_hartree.SetActive(True) se_fock.SetActive(True) gf.SelfConsistent(log='HF.log', pulay=pulay) gf.WriteSpectralInfoToNetCDFFile('HF.nc', diagonalize=True, orbitals=orbitals, spectral='individual') # GW calculation se_gw.SetActive(True) gf.SelfConsistent(log='GW.log', pulay=pulay) gf.WriteSpectralInfoToNetCDFFile('GW.nc', diagonalize=True, orbitals=orbitals, spectral='individual')
The self-consistent solution is obtained by mixing Green’s function
using a pulay mixing scheme, which is controllod by three parameters
(tol:float , mix:float : float, history:int). The self-consistent
cycles is initianted by the GF method
The self-consistent iterations can be monited if written to the
Saving calculated data to a NetCDFFile¶
The GF method
WriteSpectralInfoToNetCDFFile is used to save all
the calculated data such as spectral functions, transmission function etc.
to a NetCDFFile.