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 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:

parameter value description
t_ll -20.0 intra lead hopping
t_lm -1.0 lead-molecule hopping
t_mm -2.4 intra molecule hopping
V   electron-electron 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 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,

left = LeadHamiltonianMatrix(principallayer=1,
right = LeadHamiltonianMatrix(principallayer=1,

energies = np.arange(-96.0, 96.0, 0.02)
de = energies[1] - energies[0]
assert abs(energies).min() < 1.0e-8
fermi = 0.0
g0 = NonEqNonIntGreenFunction(hmat, left, right, E_Fermi=fermi,
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
# HF calculation
gf.SelfConsistent(log='HF.log', pulay=pulay)

# GW calculation
gf.SelfConsistent(log='GW.log', pulay=pulay)


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 selfconsistent. The self-consistent 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.