Keldysh Green functions¶
The Keldysh Green function (KGF) code allows for calculations of nonequilibrium transport calculations where electron exchange and correlation effect are threated using many body perturbation theory such as HartreeFock, 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 nonequilibrium 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 PariserParrPople (extended Hubbard) as well as abinitio calculated Hamiltonians.
A KGF calculatoin normally involves the following steps:
Setting up the noninteracting lead and scattering Hamiltonian.
Setting up a noninteracting GF
Setting up various selfenergies to handle Hartree, exchange and correlation
Runnig the calculation
Example: PariserParrModel Hamiltonian¶
To do an electron transport calculation using a model Hamiltonian the parameters of both the noninteracting part as well as the interacting part of the Hamiltonian need to be explicitly specified. The noninteracting part h_ij describe kinetic energy and electronelectron 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 electronelectron interactions are included. The o’s (dashes) represents noninteracting 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 tightbinding chains with a principal layer size of one.
The following parameters will be used to simulate a physisorbed molecule:
parameter 
value 
description 


20.0 
intra lead hopping 

1.0 
leadmolecule hopping 

2.4 
intra molecule hopping 

electronelectron interaction 
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 selfenergies. 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 noninteracting calculatoins followed by a HartreeFock calculations and a GW calculation look like this:
from __future__ import print_function
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.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[1]  energies[0]
assert abs(energies).min() < 1.0e8
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')
Selfconsistency¶
The selfconsistent 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 selfconsistent
cycles is initianted by the GF method selfconsistent
.
The selfconsistent iterations can be monited if written to the
logfiles files.
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.