NESSi  v1.0.2
The NonEquilibrium Systems Simulation Library
Manual

Preliminaries

Namespaces

cntr All relevant classes and functions of libcntr are defined here
fourier Fourier transformation routines. Used almost exclusively internally.

Needed for developers only:

integration Weights for quadrature and differentiation routines.
linalg Linear algebra routines (out-dated).

Scalar and matrix types

Complex numbers:

The numerical routines in NESSi are built on standard complex-valued algebra.

For the most commonly-used double precision numbers, the following shortcuts are defined:

#include <complex>
#define cdouble std::complex<double>
#define II cdouble(0.0,1.0)
Matrices and linear algebra:

Green's functions \(C(t,t')\) can be matrix-valued. We use the eigen3 library for handling matrices.

The following shortcuts are defined:

#include <eigen3/Eigen/Dense>
// integer precision variable-size vector
#define ivector VectorXi
// double precision variable-size vector
#define dvector VectorXd
// complex-valued double precision variable-size vector
#define cdvector VectorXcd
// integer precision variable-size matrix
#define imatrix MatrixXi
// double precision variable-size matrix
#define dmatrix MatrixXd
// complex-valued double precision variable-size matrix
#define cdmatrix MatrixXcd
Example:
cdmatrix mm(2,2); // a complex-valued 2x2 matrix
cdmatrix mm2;
// set mm to sigma_x
mm.setZero();
mm(0,1)=1;
mm(1,0)=1;
mm2=mm*mm; // mm^2
// etc ...

(For a complete documentation, see eigen3)

Some global #defines

Listed for reference only. Explanation in the sections below.

#define CNTR_PI 3.14159265358979323846
#define CNTR_MAT_FOURIER 0
#define CNTR_MAT_CG 1
#define CNTR_MAT_FIXPOINT 2
#define FERMION -1
#define BOSON 1
#define GREEN cntr::herm_matrix<double>
#define GREEN_TSTP cntr::herm_matrix_timestep<double>
#define CFUNC cntr::herm_function<double>

↑ to the top of the page ↑

Green functions

Overview

class cntr::herm_matrix<T>

This class contains the data structures for representing two-time contour functions \( C(t,t') \) on the equidistantly discretized time-contour \(\mathcal{C}\) with points \(\{i\Delta t: i=0,1,2,...,{\tt nt}\}\) along the real time branch, and \(\{i\Delta\tau,i=0,1,...,{\tt ntau}\}\) along the imaginary branch. The Green's functions are defined by the following parameters:

  • T (template parameter): Precision, usually set to double; we use the global definition
    #define GREEN cntr::herm_matrix<double>
  • nt and ntau (integer): number of distretization points on the real and imaginary time axis.
  • size1 (integer): orbital dimension. Each element \(C(t,t')\) is a square matrix of dimension sixe1 \(\times\)size1.
  • sig (integer): Takes the values FERMION (defined as -1) or BOSON (+1).

herm_matrix stores the following Keldysh components on the equidistantly discretized Keldysh contour with time discretization \(\Delta t\) and \(\Delta \tau\) on the real and imaginary time, respectively (The real-time domain of a herm_matrix object is illustrated in the Figure):

If you don't see it, something is missing
  • Matsubara component \( C^\mathrm{M}(i\Delta \tau) \) for i=0,...,ntau,
  • retarded component \( C^\mathrm{R}(i\Delta t,j\Delta t) \) for i=0,...,nt, j=0,...,i ,
  • lesser component \( C^<(i\Delta t,j\Delta t) \) for j=0,...,nt, i=0,...,j ,
  • left-mixing component \( C^\rceil(i\Delta t,j\Delta \tau) \) for i=0,...,nt, j=0,...,ntau.

If nt = -1, only the Matsubara component is stored.

Size arguments of herm_matrix are returned by the following read-only member functions:

GREEN G;
cout << "The object G of type herm_matrix has the following dimensions:" << endl;
cout << "nt= " << G.nt() << endl;
cout << "ntau= " << G.ntau() << endl;
cout << "size1= " << G.size1() << endl;
cout << "sig= " << G.sig() << endl;
Hermitian symmetry and conjugate:

In the following, we call the components mentioned above the hermitian domain of a two-time contour function. The Hermitian conjugate \(C^\ddagger\) of a Green's function \(C\) is defined ( \(\xi=\)+1 and -1 for BOSON and FERMION, respectively):

  • \( [C^\ddagger]^\mathrm{R}(t,t^\prime) = \left( C^\mathrm{A}(t^\prime,t)\right)^\dagger \ , \)
  • \( [C^\ddagger]^\gtrless(t,t^\prime) = -\left( C^\gtrless(t^\prime,t) \right)^\dagger \ , \)
  • \( [C^\ddagger]^{\rceil}(t,\tau) = - \xi \left( C^\lceil(\beta-\tau,t) \right)^\dagger \ , \)
  • \( [C^\ddagger]^{\mathrm{M}}(\tau) =\left( C^{\mathrm{M}}(\tau) \right)^\dagger \ . \)

Apparently, \(C=(C^\ddagger)^\ddagger\).

Hence, if a Green's function has hermitian symmetry ( \(C=C^\ddagger\)),it is sufficient to use one object of type herm_matrix and store the hermitian domain. If a Green's function is not hermitian, one must store the hermitian domain of \( C\) and of \( C^\ddagger\) in two separate herm_matrix variables.

For a summary of all member functions, see Summary: member functions of herm_matrix and herm_matrix_timestep[_view]. The following paragraphs give detailed explanations of some functionalities.

Constructors

herm_matrix<T>() Default constructor, does not allocate memory, sets nt=-2
herm_matrix<T>(int nt,int ntau,int size1,int sig) Allocate memory, and sets all entries to 0. It requires nt>=-1, ntau>0, size1>0, and sig \(\in\{\)FERMION,BOSON \(\}\).
herm_matrix<T>(int nt,int ntau) Equivalent to herm_matrix(int nt,int ntau,int size1,int sig) with size1=1 and sig=FERMION.

Accessing individual elements

The following routines allow to read/write the elements of a contour function \(C(t,t')\) stored as herm_matrix at individual time arguments from/to another variable M. The latter can be either a scalar of type std::complex<T>, or a complex square matrix defined in Eigen (see Scalar and matrix types).

Member functions of herm_matrix:

Member functions of herm_matrix allow to read/write elements of \(C\) in the hermitian domain. To read elements of \(C\) outside the hermitian domain, one can use the non-member element access functions defined later.

The following member functions set components of a contour function \(C\) in the hermitian domain from M:

C.set_les(i,j,M) \( C^<(i\Delta t,j\Delta t)\) is set to M required: 0 <= j <= C.nt(), 0<=i<=j
C.set_ret(i,j,M) \( C^R(j\Delta t,i\Delta t)\) is set to M required: 0 <= i <= C.nt(), 0<=j<=i
C.set_tv(i,j,M) \( C^{\rceil}(i\Delta t,j\Delta \tau)\) is set to M required: 0 <= i<= C.nt(), 0 <= j<= C.ntau()
C.set_mat(i,M) \( C^M(i\Delta \tau)\) is set to M required: 0 <= i <= C.ntau()
  • If C.size1()>1, M must be a complex eigen matrix (cdmatrix for double precision)
  • If C.size1()==1, M can be a scalar (cdouble) or a matrix
  • If M is a Matrix, it must be a square matrix of dimension size1

The following member functions read components of a contour function \(C\) in the hermitian domain to M:

C.get_les(i,j,M) M is set to \( C^<(i\Delta t,j\Delta t)\) required: 0 <= j <= C.nt(), 0<=i<=j
C.get_ret(i,j,M) M is set to \( C^R(i\Delta t,j\Delta t)\) required: 0 <= i <= C.nt(), 0<=j<=i
C.get_tv(i,j,M) M is set to \( C^{\rceil}(i\Delta t,j\Delta \tau)\) required: 0 <= i<= C.nt(), 0 <= j<= C.ntau()
C.get_mat(i,M) M is set to \( C^M(i\Delta \tau)\) required: 0 <= i <= C.ntau()
  • If M is a matrix, it is resized to a square matrix of dimension C.size1()
  • If M is a scalar, only the (0,0) entry of \(C(t,t')\) is read.

Remark: Member functions set_les and set_ret can return also values outside the hermitian domain, but this functionality is no longer supported. Similarly, there are non-supported member functions like get_gtr(...) which return other Keldysh components. Use instead the general access functions described in the next paragraph.

General element access:

The following functions read components of a contour function \(C\) to M:

cntr::get_les(i,j,M,C,Ccc) M is set to \( C^<(i\Delta t,j\Delta t)\) required: 0 <= i,j <= C.nt()
cntr::get_gtr(i,j,M,C,Ccc) M is set to \( C^>(i\Delta t,j\Delta t)\) required: 0 <= i,j <= C.nt()
cntr::get_ret(i,j,M,C,Ccc) M is set to \( C^>(i\Delta t,j\Delta t)-C^<(i\Delta t,j\Delta t)\) required: 0 <= i,j <= C.nt()
cntr::get_tv(i,j,M,C,Ccc) M is set to \( C^{\rceil}(i\Delta t,j\Delta \tau)\) required: 0 <= i<= C.nt(), 0 <= j<= C.ntau()
cntr::get_vt(i,j,M,C,Ccc) M is set to \( C^{\lceil}(i\Delta \tau,j\Delta t)\) required: 0 <= i<= C.ntau(), 0 <= j<= C.nt()
cntr::get_mat(i,M,C,Ccc) M is set to \( C^M(i\Delta \tau)\) required: 0 <= i <= C.ntau()
  • C and Ccc are (size-matched) arguments of type herm_matrix which store the hermitian domain of \(C\) and its hermitian conjugate \(C^\ddagger\).
  • If the argument Ccc is omitted, hermitian symmetry \(C=C^\ddagger\) is assumed.
  • If M is a matrix, it is resized to a square matrix of dimension C.size1().
  • M may always be a scalar; in this case only the (0,0) entry of \(C(t,t')\) is read.
Note
All set_XXX and get_XXX element access functions are not optimised for massive read/write operations. (In particular, set_XXX and get_XXX are not used inside the numerically costly tasks like the solution of integral equations). Mostly, data are anyway read and written timestep-wise (see Timeslices).

Density matrix

The equal-time elements of Green's functions contain expectation values of physical observables. If \(C_{ab}(t,t') = -i \langle T_\mathcal{C} \hat A_a(t)\hat B_b(t')\rangle\) is the contour ordered expectation value of two operators \(\hat A\) and \(\hat B\) (such as creation and annihilation operators with orbital indices a and b) then the time-dependent expectation value \(\rho_{ab}(t) = \langle B_{b} A_{a} \rangle_t \) is given by ( \(\xi\) is the BOSON/FERMION sign)

\begin{align} &\rho_{ab}(t) = i \xi C^<_{a,b}(t,t) \,\,\,\text{for general real times} \\ &\rho_{ab} = - C^M_{a,b}(\beta) \,\,\,\text{for the equilibrium initial state} \end{align}

We therefore define the Density matrix of a two-time Green's function on the discretized contour at timestep tstp as

\begin{align} \text{DensityMatrix}[C]({\tt tstp}) & = - C^M({\tt ntau}\,\Delta\tau) \,\,\,\text{for the equilibrium initial state}\,({\tt tstp}=-1) \\ & = i \xi C^<_{a,b}({\tt tstp}\Delta t,{\tt tstp}\Delta t) \,\,\,\text{for}\,{\tt tstp}>=0. \end{align}

The following member function of a cntr::herm_matrix<T> object C writes the Density matrix to an object M:

C.density_matrix(int tstp,M) M is set to \(\text{DensityMatrix}[C]({\tt tstp})\)
  • tstp <= C.nt() required.
  • If C.size1()>1, M must be a complex eigen matrix (cdmatrix for double precision)
  • If C.size1()==1, M can be a scalar (cdouble) or a matrix
  • If M is a matrix, it is resized to a square matrix of dimension C.size1().

File I/O

Reading/writing to text files

The member functions print_to_file and read_to_file of a herm_matrix implement text-file access, as apparent by example:

int nt=10;
int ntau=10;
int size1=2;
int sig=FERMION;
GREEN A(nt,ntau,size1,sig);
// ... set elements of A ...
//WRITE:
// create a file filename.txt and store the data of A:
A.print_to_file("filename.txt");
// READ
GREEN B;
// if the file "filename.txt" has been written previously with print_to_file ,
// the parameters (nt,ntau,size1,sig) and the data of B are modified
// according to the information in the file:
B.read_from_file("filename.txt");
// now B is resized to B.nt()=10, B.ntau()=10,B.size1()=2,B.sig()=FERMION, and the data of B match the data of A.

The file format is rather explicit (and storage intensive):

  • First line: # nt ntau size1 sig
  • The following lines list the Matsubara component in the format:
    mat: i Re_C^mat(i)_{0,0} Im_C^mat(i)_{0,0} ... Re_C^mat(i)_{size1,size1} Im_C^mat(i)_{size1,size1}
  • If nt>-1, the following lines list the other components in the format:
    XXX: i j Re_C^XXX(i,j)_{0,0} Im_C^XXX(i,j)_{0,0} ... Re_C^XXX(i,j)_{size1,size1} Im_C^XXX(i,j)_{size1,size1}
    where XXX is ret for \(C^R\), les for \(C^<\), tv for \(C^{\rceil}\), and i,j loop through the corresponding time arguments of the hermitian domain.
Note
Writing to text files is supposed as a quick way to generate human-readible data. For compressed storage, one should use the hdf5 format below.
Reading/writing to hdf5 files

The hdf5 format, as well as some scripts to interpret the hdf5 files are discussed in a separate section hdf5 usage below.

↑ to the top of the page ↑

Timeslices

Overview

class cntr::herm_matrix_timestep<T>

For a Green's function \(C\) of type cntr::herm_matrix, we define the timestep tstp, with tstp \(\in\{\)-1,0,1,...,nt \(\}\), as the following data set:

  • If tstp=-1, the timestep contains the Matsubara Green's function:
    • Matsubara component \( C^\mathrm{M}(j\Delta \tau) \) for j=0,...,ntau,
  • If tstp>=0, it contains the following real-time components:
    • \( C^\mathrm{R}({\tt tstp}\Delta t,j \Delta t) \) for j=0,...,tstp,
    • \( C^<(j \Delta t,{\tt tstp} \Delta t) \)for j=0,...,tstp,
    • \( C^\rceil({\tt tstp}\Delta t,j\Delta \tau) \) for j=0,...,ntau.

The data of the herm_matrix Green's function \(C\) are therefore precisely a union of the timesteps -1,...,nt. The class cntr::herm_matrix_timestep<T> is the container to store these data. It is characterized by the following parameters:

  • T (template parameter): Precision, usually set to double; we use the definition
    #define GREEN_TSTP cntr::herm_matrix_timestep<double>
  • tstp (integer): The timestep variable, tstp>=-1.
  • ntau (integer): number of discretization points on the imaginary time axis, `ntau>=0
  • size1 (integer): orbital dimension. Each element \(C(t,t')\) is a square matrix of dimension sixe1 \(\times\)size1.
  • sig (FERMION or BOSON).

For a summary of all member functions, see Summary: member functions of herm_matrix and herm_matrix_timestep[_view]. The Following paragraphs give detailed explanations of some functionalities.

Constructors

herm_matrix_timestep<T>() Default constructor, does no allocate memory and sets nt=-2.
herm_matrix_timestep<T>(int tstp,int ntau,int size1,int sig) Allocate memory, set all entries to 0. It requires tstp>=-1, ntau>0, size1>0, and sig=FERMION or sig=BOSON.

Accessing individual matrix elements.

Element access, as well as read-out of the density matrix follows the same syntax as for the Green's function herm_matrix<T>. See Sections Density matrix and Accessing individual elements.

Example 1:
int ntau=100;
int size=2;
int tstp=10;
cdmatrix M(size1,size); // an Eigen matrix
GREEN_TSTP tA(tstp,ntau,size,FERMION); // allocate a timestep
// ... set entries of M ...
tA.set_les(8,tstp,M); // ok: A^<(8,10) set to M
tA.set_les(tstp,8,M); // error: time arguments do not match the domain of tA
cntr::get_les(tstp,8,M,tA,tA); // ok: M set to A^<(8,10), assuming that A is hermitian
Example 2:

simple print of all elements in case of a scalar Green's function, at tstp=-1 (Matsubara):

int ntau=100;
int size=1;
tstp=-1;
GREEN_TSTP tA(tstp,ntau,size,FERMION);
// ... do something with tA ...
for(int j=0;j<=ntau;j++){
tA.get_mat(j,x);
cout << "Amat(j) at j=" << j << " : " << x << endl;
}

Remark: Some un-supported routines exist in which the dummy argument tstp can be omitted.

Timestep-wise data exchange between Green's functions

In typical applications, data between different Green's functions are exchanged at once for an entire timestep. We provide two main member functions which allow the exchange of data of one entire timestep between variables of type cntr::herm_matrix and cntr::herm_matrix_timestep:

A.set_timestep(tstp, B) Copies the data of timestep tstp from B to A
  • A and B are either cntr::herm_matrix<T> or cntr::herm_matrix_timestep<T>
    Types can be mixed to enable data exchange between the types in both directions
  • tstp (integer): The timestep to be accessed
  • For tstp==-1 this copies the data from B to A as:
    • \(A^M(i\Delta \tau) = B^M(i\Delta \tau)\) for i=0,...,ntau
  • For tstp>=0 this copies the data from B to A as:
    • \(A^R({\tt tstp}\Delta t,j \Delta t) = B^R({\tt tstp}\Delta t,j \Delta t)\) for j=0,...,tstp
    • \(A^<(j \Delta t,{\tt tstp}\Delta t) = B^<(j \Delta t, {\tt tstp}\Delta t)\) for j=0,...,tstp
    • \(A^{\rceil}({\tt tstp}\Delta t,j\Delta\tau) = B^{\rceil}({\tt tstp}\Delta t,j\Delta\tau)\) for j=0,...,ntau
  • Size consistency is required:
    • If A is of type herm_matrix_timestep, tstp==A.tstp() is required, same for B
    • If A is of type herm_matrix, tstp<=A.nt() is required, same for B
    • A.ntau()==B.ntau() and A.size1()==B.size1() is required
Example:

The following code copies the first 5 timesteps and the Matsubara timestep from a Green's function B to a Green's function A:

int nt=10;
int ntau=10
int size1=2;
GREEN A(nt,ntau,size1,FERMION);
GREEN B(5,ntau,size1,FERMION);
// ... do something with B ...
for(int tstp=-1;tstp<=5;tstp++) A.set_timestep(tstp,B);


A.set_matrixelement(tstp, i1,i2,B,j1,j2) Read out only certain(orbital) matrix elements at timestep tstp from B and write them to A
  • A and B are either cntr::herm_matrix or cntr::herm_matrix_timestep; The two types can be mixed.
  • tstp (integer): The timestep to be accessed
  • i1, i2, j1, j2 (integer): Orbital matrix indices
  • For tstp==-1 this copies the data from B to A as:
    • \(A^M(i\Delta \tau)_{i1,i2} = B^M(i\Delta \tau)_{j1,j2}\) for i=0,...,ntau
  • For tstp>=0 this copies the data from B to A as:
    • \(A^R({\tt tstp}\Delta t,i \Delta t)_{i1,i2} = B^R({\tt tstp}\Delta t,i \Delta t)_{j1,j2}\) for i=0,...,tstp
    • \(A^<(i \Delta t,{\tt tstp}\Delta t)_{i1,i2} = B^<(i \Delta t,{\tt tstp}\Delta t)_{j1,j2}\) for i=0,...,tstp
    • \(A^{\rceil}({\tt tstp}\Delta t,i\Delta\tau)_{i1,i2} = B^{\rceil}({\tt tstp}\Delta t,i\Delta\tau)_{j1,j2}\) for i=0,...,ntau
  • size consistency is required:
    • If A is of type herm_matrix_timestep, tstp==A.tstp() is required, same for B
    • If A is of type herm_matrix, tstp<=A.nt() is required, same for B
    • A.ntau()==B.ntau() is required
    • 0 <= i1,i2 <= A.size1(), 0 <= j1,j2 <= B.size1()

set_matrixelement is useful in the context of constructing Feynman diagrams with multiorbital Green's functions.

Example:

The following (little useful) code permutes the (0,1) and (1,0) matrix elements of timestep 7 of a herm_matrix Green's function, using the herm_matrix_timestep as a temporary variable.

int nt=10;
int ntau=10
int size1=2;
int sig=FERMION;
GREEN G(nt,ntau,size1,sig);
// ... do something with G
{
int tstp=7;
GREEN tG(tstp,ntau,1,sig); // allocate a timestep of orbital dimension 1 as temp. storage
tG.set_matrixelement(tstp,0,0,G,0,1); // save (0,1) matrix element
G.set_matrixelement(tstp,0,1,G,1,0); // copy (1,0) matrix element to (0,1)
G.set_matrixelement(tstp,1,0,tG,0,0); // copy (0,0) of tG to (1,0) matrix-element of G
}

Map

Overview

class cntr::herm_matrix_timestep_view<T>

The purpose of this class is to create a map to an pre-existing object of type cntr::herm_matrix or cntr::herm_matrix_timestep. One can then use this class without making a physical copy of the original data. The structure is inherited from the existing object. The usage of this class is mainly reserved for active developers as it is employed to enhance a performance of low-lying routines and MPI communications. It is characterized by the following parameters:

  • T (template parameter): Precision, usually set to double; we use the definition
    #define GREEN_TSTP_VIEW cntr::herm_matrix_timestep_view<double>
  • tstp (integer): The timestep variable, tstp>=-1.
  • ntau (integer): number of discretization points on the imaginary time axis, `ntau>=0
  • size1 (integer): orbital dimension. Each element \(C(t,t')\) is a square matrix of dimension sixe1 \(\times\)size1.
  • sig (FERMION or BOSON).

Constructors

herm_matrix_timestep_view<T>() Default constructor, does no allocate memory and sets nt=-2.
herm_matrix_timestep_view<T>(int tstp,int ntau,int size1, int size2, int sig) Creates a map with all entries to 0. It requires tstp>=-1, ntau>0, size1>0, size2>0, and sig=FERMION or sig=BOSON.
herm_matrix_timestep_view(int tstp, herm_matrix<T> &g) Creates a map to the predefined herm_matrix object g at a given time tstp. Similar constructor exist for herm_matrix_timestep.

Accessing and manipulation

In general, all functions decribed for herm_matrix_timestep argument in Overview can be replaced by a herm_matrix_timestep_view argument following the same syntax.

↑ to the top of the page ↑

Summary: member functions of herm_matrix and herm_matrix_timestep[_view]

The following table gives a brief summary of member functions of a herm_matrix and herm_matix_timesstep object A, with reference to the documentation section:

Constructors and file i/o of herm_matrix
herm_matrix<T>() Default constructor Constructors
herm_matrix<T>(int nt,int ntau,int size1,int sig) Constructor Constructors
void print_to_file(const char* filename) create a textfile filename and write human readable data (large files!) File I/O
void read_from_file(const char* filename) read a textfile filename previously created with print_to_file, and resize/set A to content of file File I/O
void write_to_hdf5(const char *filename, const char *groupname) write a data group groupname into a hdf5 file filename hdf5 usage
void read_from_hdf5(const char *filename, const char *groupname) read a group groupname from hdf5 file filename and resize/set A to content of file hdf5 usage
Constructors of herm_matrix_timestep
herm_matrix_timestep() Default constructor Constructors
herm_matrix_timestep(int nt,int ntau,int size1,int sig) Constructor Constructors
Constructors of herm_matrix_timestep_view
herm_matrix_timestep_view() Default constructor Constructors
herm_matrix_timestep_view(int nt,int ntau,int size1,int size2,int sig) Constructor Constructors
herm_matrix_timestep_view(int tstp, GG &g) Constructor, GG = herm_matrix<T> or herm_matrix_timestep<T> Constructors
Member functions of herm_matrix, herm_matrix_timestep, herm_matrix_timestep_view with the same syntax
void set_[XXX](time_arguments,matrix &M) Set component XXX=mat,ret,les,tv at given time arguments of the hermitian domain from a matrix or scalar M Accessing individual elements
void get_[XXX](time_arguments,matrix &M) read component XXX=mat,ret,les,tv at given time arguments of the hermitian domain to a matrix or scalar M Accessing individual elements
void density_matrix(int tstp,matrix &M) read density matrix at time tstp to a matrix or scalar M Density matrix
void set_timestep(int tstp,GG &B) set data at timestep tstp to data of B, where B is herm_matrix<T> or herm_matrix_timestep<T> Timestep-wise data exchange between Green's functions
void set_matrixelement(int tstp,int i1,int i2,GG &B,int j1,int j2) set matrixelement (i1,i2) at timestep tstp to matrixelemnt (ja,j2) of B, where B is herm_matrix<T> or herm_matrix_timestep<T> Timestep-wise data exchange between Green's functions
void set_timestep_zero(int tstp) Set elements at timestep tstp to 0 Simple algebra
void smul(int tstp, std::complex<T> x) \(A(t,t')\) set to \( xA(t,t') \) on timestep tstp of A Simple algebra
void incr_timestep(int tstp,GG &B, std::complex<T> x) \(A(t,t')\) set to \( A(t,t')+xB(t,t') \) on timestep tstp of A, where B is herm_matrix<T> or herm_matrix_timestep<T> Simple algebra
void left_multiply(int tstp,cntr::function<T> &f, T x) \( A(t,t')\) set to \( x f(t) A(t,t')\) at timestep tstp of A Multiplication with one-time contour functions
void right_multiply(int tstp,cntr::function<T> &f, T x) \( A(t,t')\) set to \( x A(t,t') f(t')\) at timestep tstp of A Multiplication with one-time contour functions
void left_multiply_hermconj(int tstp,cntr::function<T> &f, T x) \( A(t,t')\) set to \( x f(t)^\dagger A(t,t')\) at timestep tstp of A Multiplication with one-time contour functions
void right_multiply_hermconj(int tstp,cntr::function<T> &f, T x) \( A(t,t')\) set to \( x A(t,t') f(t')^\dagger\) at timestep tstp of A Multiplication with one-time contour functions
void Bcast_timestep(int tstp, int root) Broadcast the data of \(A(t,t')\) on timestep tstp from rank root to all other ranks. Communicating timesteps
void Reduce_timestep(int tstp, int root) In place reduction, replacing \(A(t,t')\) on rank root by \( \sum_{\text{rank }\,j} A_{\text{at rank}\,j}(t,t')\) for one timestep tstp. Communicating timesteps
void Send_timestep(int tstp, int dest, int tag) Send the data of \(A(t,t')\) on timestep tstp to rank dest with a message tage tag. Communicating timesteps
void Recv_timestep(int tstp, int root, int tag) Receive a message with tag tag from rang root which contains the data of \(A(t,t')\) on timestep Communicating timesteps

↑ to the top of the page ↑

Contour Functions

Overview

class cntr::function<T>

A contour function is a matrix or scalar-valued function \(f(t)\) which depends on physical time only:

  • On the imaginary time branch, it takes one given constant value \( f_{-1} \)
  • On the real-time branches, its value is the same for the same time \(t\) on the upper and lower contour branch.

A contour function is typically used to store time-dependent parameters of a system, or a time-dependent Hamiltonian.

The class cntr::function<T> is the container to store these data. It is characetrized by the following parameters:

  • T (template parameter): Precision, usually set to double; we use the definition
    #define CFUNC cntr::herm_function<double>
  • nt (integer): number of discretization points on the real time axis.
    For nt=-1, only the (constant) value on the imaginary time axis is stored.
  • size1 (integer): orbital dimension. Each element \(f(t)\) is a square matrix of dimension sixe1 \(\times\)size1.

Constructors

function<T>() Default constructor, does not allocate memory and sets nt=-2.
function<T>(int nt,int size1) Allocate memory, set all entries to 0. It requires tstp>=-1, size1>0

Accessing individual matrix elements.

The following routines allow to read/write the elements of a contour function \(f(t)\) stored as cntr::function at individual time arguments from/to another variable M. The latter can be either a scalar of type std::complex<T>, or a complex square matrix defined in Eigen (see Scalar and matrix types)

The following member functions of cntr::function<T> set components of a contour function \(f(t)\) from M:

f.set_value(j,M) For j=-1 \( f_{-1}\) is set to M
For 0<=j<=f.nt(), \( f(j\Delta t)\) is set to M
  • If f.size1()>1, M must be a square matrix (cdmatrix for double precision)
  • If f.size1()==1, M can be also a scalar (cdouble) or a square matrix
  • If M is a matrix, it must be a square matrix of dimension f.size1()

The following member functions of cntr::function<T> read components of a contour function \(f\) to M:

f.get_value(j,M) For j=-1, M is set to \( f_{-1}\)
For 0<=j<=f.nt(), M is set to \( f(j\Delta t)\)

If M is a matrix, it is resized to a square matrix of dimension f.size1()

  • If M is a scalar, only the (0,0) entry of \(f(t)\) is read.

↑ to the top of the page ↑

Simple algebra

Simple algebra operations include scalar multiplication, incrementation of a Green's function with another Green's function, and multiplication of Green's function with one-time functions. These operations are performed usually on one timestep of a Greens function. They are implemented as member functions of the cntr::herm_matrix and cntr::herm_matrix_timestep classes, with the same syntax for both objects.

Scalar multiplication and incrementation

A.smul(int tstp, std::complex<T> x) \(A(t,t')\) set to \( xA(t,t') \) on timestep tstp of A
A.set_timestep_zero(int tstp) \(A(t,t')\) set to \(0 \) on timestep tstp of A
A.incr_timestep(int tstp,GG &B, std::complex<T> x) \(A(t,t')\) set to \( A(t,t')+xB(t,t') \) on timestep tstp of A
  • A and B can be of type herm_matrix<T> and herm_matrix_timestep<T>
  • If A is of type herm_matrix_timestep, then A.tstp()==tstp is required (same for B)
    If A is of type herm_matrix, then A.nt()>=tstp is required (same for B)
  • size consistency B.ntau()==A.ntau() and B.size1()==A.size1() required.
Example:

The following operation considers a set of Green's functions \( G_k\) indexed by momentum \(k\), and computes the Fourier transform at a given j=2:

\begin{equation} G_{j}(t,t') = \frac{1}{N_k}\sum_{k=0}^{N_k-1} e^{i2\pi k j /N_k} G_k(t,t').\end{equation}

int nt=10;
int ntau=100;
int size1=1;
int nk=10;
GREEN_TSTP tGav;
std::vector<GREEN> Gk(nk); // a vector of nk=10 Green's functions
for(int k=0;k<nk;k++) Gk[k]=GREEN(nt,ntau,size1,FERMION); // allocation
// ...
// ... do something with Gk
// ...
// do Fourier transform on timestep -1 (Matsubara)
// and store result in timestep variable tGav:
{
double nkinv=1.0/nk;
int j=2;
tstp=-1;
tGav=GREEN_TSTP(tstp,ntau,size1,FERMION); // properly resize tGav
tGav.set_timestep_zero(tstp); // in principle not needed, as tGav is initialized with 0 already
for(int k=0;k<nk;k++){
cdouble weight=exp(II*2.0*M_PI*j*k*nkinv);
tGav.incr_timestep(tstp,Gk[k],weight);
}
tGav.smul(tstp,nkinv);
}
// the syntax is analogous if tGav is replaced by a herm_matrix type Green's function

Multiplication with one-time contour functions

A.left_multiply(int tstp,cntr::function<T> &f, T x) \( A(t,t')\) set to \( x f(t) A(t,t')\) at timestep tstp of A
A.right_multiply(int tstp,cntr::function<T> &f, T x) \( A(t,t')\) set to \( x A(t,t') f(t')\) at timestep tstp of A
A.left_multiply_hermconj(int tstp,cntr::function<T> &f, T x) \( A(t,t')\) set to \( x f(t)^\dagger A(t,t')\) at timestep tstp of A
A.right_multiply_hermconj(int tstp,cntr::function<T> &f, T x) \( A(t,t')\) set to \( x A(t,t') f(t')^\dagger\) at timestep tstp of A
  • A is of type cntr::herm_matrix<T> or cntr::herm_matrix_timestep<T>
  • x is a scalar of type T, can be omitted (default x=1)
  • The following size consistency is required:
    • If A is herm_matrix_timestep, then tstp==A.tstp() is required
      If A is herm_matrix, then tstp<=A.nt() is required
    • tstp<=f.nt() and f.size1()==A.size1() is required.
  • The last two routines multiply wth the adjoint matrix \(f(t)^\dagger\) of \(f(t)\)
Example:

The following example computes a retarded interaction \( W(t,t') = U(t) \chi(t,t') U(t')\) out of a given susceptibility \(\chi\), where \( U(t)\) is a time-dependent bare interaction matrix element.

int nt=10;
int ntau=100;
int size1=1;
GREEN W(nt,ntau,size1,BOSON);
GREEN chi(nt,ntau,size1,BOSON);
CFUNC U(nt,size1);
// ...
// ... do something to set U and chi
// ...
// compute W on all timesteps
for(int tstp=-1;tstp<=nt;tstp++){
W.set_timestep(tstp,chi);
W.left_multiply(tstp,U);
W.right_multiply(tstp,U);
}

↑ to the top of the page ↑

Free Green's functions

We provide a number of tools for calculating GFs for non-interacting systems, either from a single-particle Hamiltonian or from a given spectral function.

Equilibrium Green's functions for a given density of states

In many applications (for representing a fermionic bath, for example), one needs to construct GFs with a given density of states (DoS) \(A(\omega)\):

\begin{align} G(t,t') = \int d \omega \,A(\omega) g_{\omega}(t,t') \ , \end{align}

where \(g_\omega(t,t')\) is the single-orbital noninteracting GF with energy \(\omega\). The latter is defined as ( \(F_\xi(x)=1/(e^{\beta x}-\xi )\) is the Fermi/Bose distribution):

  • \(g_\omega^M (\tau) = - e^{-\tau(\omega-\mu)} F_\xi(-(\omega-\mu))\)
  • \(g_\omega^R (t,t') = -i e^{-i(t-t')\omega}\)
  • \(g_\omega^< (t,t') = -i\xi e^{-i(t-t')\omega} F_\xi(\omega-\mu)\)
  • \(g_\omega^{\rceil} (t,\tau) = -i\xi e^{-it\omega} e^{\tau(\omega-\mu)} F_\xi(\omega-\mu)\)

This task is accomplished by cntr::green_equilibrium:

cntr::green_equilibrium(herm_matrix<T> &G, dos_function &dos, double beta, double h, double mu=0.0, int limit = 100, int nn = 20) Computes the free GF with given density of states dos. output: object G of type cntr::herm_matrix

The user provides the inverse temperature beta, the time step h and the chemical potential mu (set to 0.0 by default). Parameters limit and nn determine the adapative sampling method implemented in fourier. Use the default settings unless there is a special to reason not to. dos is an object representing a class DoS. This can represent an arbitrary density of state. As a useful example to show the minimal call properties, we provide the Bethe DoS \(A(\omega) = V^2/(2\pi) \sqrt{ 4 V^2 - \omega^2} \), represented by the class

class bethedos {
public:
double hi_; // Higher edge of the density of states
double lo_; // Lower edge of the density of states
double V_; // 4V corresponds to the bandwidth
bethedos() {
V_ = 1;
lo_ = -2;
hi_ = 2;
}
double operator()(double x) {
double arg = 4.0 * V_ * V_ - x * x;
double num = V_ * V_ * 3.14159265358979323846 * 2;
return (arg < 0 ? 0.0 : sqrt(arg) / num);
}
};

The member variables lo_, hi_ and the function double operator()(double x) are the basic building blocks for any DoS class. For constructing a GF for a Bethe DoS with band width of 1 centered around \(\omega=0\), we would use the following code:

// define parameters, nt, ntau, beta, ...
// define herm_matrix G
dos.V_ = 0.25;
dos.lo_ = -0.5;
dos.hi_ = 0.5;
cntr::green_equilibrium(G, dos, beta, h);

We also provide the function cntr::green_equilibrium_bethe, which constructs the GF for the Bethe DoS with V_=1.0 without invoking the dos class directly.

Another example of a useful DoS is a smoothened box-like shape \(A(\omega) = C (f_\nu(\omega-b) - f_\nu(\omega-a))\), where \(f_\nu(\omega)\) is a Fermi function with effective temperature \(\nu\), while \([a, b]\) defines the interval for DoS. Further details can be found in the documentation of cntr::smooth_box.

Green's function for a constant or time-dependent Hamiltonian

Another often encountered example is the GF of a single-particle Hamiltonian \(\epsilon(t)\), defined by

\begin{align} \label{tttsfwg} \big[ i \partial_t + \mu - \epsilon(t) \big] G(t,t')= \delta_{\mathcal{C}}(t,t') \ . \end{align}

We provide routines for both constant and time-dependent Hamiltonians:

cntr::green_from_H(herm_matrix<T> &G,T mu,cdmatrix &eps,T beta,T h) Computes the free GF for constant Hamiltonian eps. output: object G of type cntr::herm_matrix
cntr::green_from_H(int tstp, GG &G,T mu,cdmatrix &eps,T beta,T h) Computes the free GF for constant Hamiltonian eps at time step tstp. output: object G
cntr::green_from_H(herm_matrix<T> &G,T mu,cntr::function<T> &eps,T beta,T h,bool fixHam,int SolveOrder,int cf_order) Computes the free GF for the time-dependent Hamiltonian eps. output: object G of type cntr::herm_matrix
cntr::green_from_H(herm_matrix<T> &G,T mu,cntr::function<T> &eps,T beta,T h,int SolveOrder,int cf_order) Computes the free GF for the time-dependent Hamiltonian eps. output: object G of type cntr::herm_matrix
cntr::green_from_H(int tstp, GG &G,T mu,cntr::function<T> &eps,T beta,T h,bool fixHam,int SolveOrder,int cf_order)

Computes the free GF for the time-dependent Hamiltonian eps at time step tstp. output: object G

  • In the routines with a \(tstp\)-Argument, GG can be herm_matrix_timestep<T> or herm_matrix<T>
  • For a time-dependent Hamiltonian, eps stores the Hamiltonian at all times. The GF is obtained by numerically solving the equation of motion \eqref{tttsfwg}. The numerical parameters (cf_order,SolveOrder,fixHam) are explained in the second Example below.

Example 1: The GF for a time-independent Hamiltonian can be constructed as follows:

int nt=100;
int ntau=100;
int size1=2;
double mu=0.0;
double h=0.05;
double beta=10.0;
GREEN G(nt,ntau,size1,FERMION);
cdmatrix eps0(size1,size1);
// define the Hamiltonian
eps0(0,0) = -1.0;
eps0(0,1) = 0.5 + 0.5*II;
eps0(1,0) = std::conj(eps(0,1));
eps0(1,1) = 1.0;
// construct GF
cntr:::green_from_H(G, mu, eps0, beta, h);

Example 2: For time-dependent Hamiltonians, we have implemented the commutator-free decomposition of the time-evolution operator from ref. [1] in second (cf_order=2) and fourth (cf_order=4) order. By default, the fourth-order scheme is used. The method requires the evaluation of \(\epsilon(t)\) at times between the interval \([(n-1) h, n h]\); thus, we use polynomial interpolation of order SolveOrder. By default, the maximum order is used. The following examples shows how to use cntr::green_from_H for this case:

// define nt, ntau, h, beta, mu ...
int size1=2;
cdmatrix eps0(size1,size1);
// defines eps0 as in the above example
cdmatrix eps1(size1,size1);
eps1(0,1) = 1.0;
eps1(1,0) = 1.0;
GREEN G(nt,ntau,size1,FERMION);
CFUNC eps(nt, size1);
eps.set_value(-1, eps0);
cdmatrix eps_at_t(size1,size1);
for(int tstp=0; tstp<=nt; tstp++){
eps_at_t = eps0 + sin(tstp*h)*eps1;
eps.set_value(tstp, eps_at_t);
}
cntr::green_from_H(G, mu, eps, beta, h, true);

A typical example where the Hamiltonian is not yet known at a certain time step is a time-dependent mean-field calculation, where the density matrix at \(t_n = nh \) is needed to construct \(\epsilon(n h)\). For this purpose, we have implemented a polynomial extrapolation as a predictor for \(\epsilon(n h)\), triggered by the argument fixHam=False. In constrast, fixHam=True indicates that \(\epsilon(n h)\) is already known and an interpolation in \([(n-1) h, n h]\) can be employed. The following code snippet exemplifies a corresponding predictor-corrector scheme:

// define nt, ntau, h, beta, mu ...
// define SolveOrder
GREEN G(nt,ntau,size1,FERMION);
CFUNC eps(nt, size1); c
// determine eps, G for time steps tstp=-1,...,SolveOrder
for(int tstp=SolveOrder+1; tstp<=nt; tstp++){
// predictor
cntr::green_from_H(tstp, G, mu, eps, beta, h, false, SolveOrder);
// update density from G at tstp, update eps at tstp
// corrector
cntr::green_from_H(tstp, G, mu, eps, beta, h, true, SolveOrder);
}

Bosonic Green's functions

All the above functions work for fermionic or bosonic GFs. However, for problems with electron-phonon coupling (or similar effects), we provide functions for calculating the free phonon propagator \(D_0(t,t') = -i\langle T_\mathcal{C} \hat{X}(t)\hat{X}(t')\rangle \), with \(\hat{X} = (\hat{a} + \hat{a}^\dagger)/\sqrt{2}\):

cntr::green_single_pole_XX_timestep(int tstp, herm_matrix_timestep<T> &D0, double w, double beta, double h) Computes the free phonon-type propagator with frequency w at time step tstp. output: object D0 of type cntr::herm_matrix_timestep
cntr::green_single_pole_XX_timestep(int tstp, herm_matrix<T> &D0, double w, double beta, double h) Computes the free phonon-type propagator with frequency w at time step tstp. output: object D0 of type cntr::herm_matrix
cntr::green_single_pole_XX_timestep(herm_matrix<T> &D0, double w, double beta, double h) Computes the free phonon-type propagator with frequency w for all time steps. output: object D0 of type cntr::herm_matrix

↑ to the top of the page ↑

Diagram utilities

The basic building blocks of Feynman diagrams are particle-hole bubbles and particle-particle bubbles, such as:

Particle hole and particle particle diagrams

The figure shows a typical particle-hole (left) and particle-particle (right) bubble, where \(V^c_{ab}\) is an interaction vertex depending on three orbital indices, and \(A\) and \(B\) are Green's functions. The diagrams above correspond to the following algebraic expressions:

\begin{equation}\label{diagram1} \text{particle-hole bubble:}\,\,\,C_{c,c'}(t,t') = \sum_{a,a'} \sum_{b,b'} V^c_{a,b} \Big[ iA_{a,a'}(t,t') B_{b',b}(t',t) \Big] V^{c'}_{a',b'}, \end{equation}

\begin{equation}\label{diagram2} \text{particle-particle bubble:}\,\,\,C_{c,c'}(t,t') = \sum_{a,a'} \sum_{b,b'} V^c_{a,b} \Big[ i A_{a,a'}(t,t') B_{b,b'}(t,t') \Big] V^{c'}_{a',b'}. \end{equation}

Note on hermitian conjugation: The complex \(i\) factor has been inserted for convenience. With this definition, the hermitian conjugate \(C^\ddagger\) of \(C\) is obtained simply by replacing A and B by their hermitian conjugates, and \( V^c_{a,b}\) by its element-wise complex conjugate. In particular, if \(V\) is real and A and B are hermitian symmetric, then also C is hermitian symmetric (as a matrix-valued GF; its individual components are not hermitian, but \( (C_{1,2})^\ddagger = C_{2,1}\)).

The basic building blocks are therefore products of the type \(C(t,t')=iA_{a,a'}(t,t') B_{b',b}(t',t)\) and \(C(t,t')=iA_{a,a'}(t,t') B_{b,b'}(t,t')\), where \(t,t'\) are arbitrary times on the contour. To evaluate those products, one uses Langreth rules. For example

\begin{equation}\label{langreth-product-1} C(t,t')=iA(t,t')B(t',t) \,\,\,\Rightarrow\,\,\, C^<(t,t')=iA^<(t,t')B^>(t',t)\end{equation}

Detailed expressions are given in Ref. [6], Chapter 3.

We provide two basic functions that allow to compute particle-particle and particle-hole Bubbles on a given timestep:

cntr::Bubble1(int tstp, GG &C,int c1,int c2,GG &A,GG &Acc,int a1,int a2,GG &B,GG &Bcc,int b1,int b2) \(C_{c1,c2}(t,t')=i A_{a1,a2}(t,t') B_{b2,b1}(t',t) \) is calculated on timestep tstp of C for given two-time objects A,B with matrix indices a1, a2 and b1, b2, respectively. output: object C with matrix indices c1, c2
cntr::Bubble2(int tstp, GG &C,int c1,int c2,GG &A,GG &Acc,int a1,int a2,GG &B,GG &Bcc,int b1,int b2) \(C_{c1,c2}(t,t')=i A_{a1,a2}(t,t') B_{b1,b2}(t,t') \) is calculated in on timestep tstp of C for given two-time objects A,B with matrix indices a1, a2 and b1, b2, respectively. output: object C with matrix indices c1, c2
  • C, A, Acc, B,Bcc can be of type cntr::herm_matrix or cntr::herm_matrix_timestep
  • For X=A or B: Xcc contains the hermitian conjugate \(X^\ddagger\) of \(X\). If the argument Xcc is omitted, it is assumed that X is hermitian, \(X=X^\ddagger\).
  • If the index pairs (c1,c2), (a1,a2), or (b1,b2) are omitted, they are assumed to be (0,0)
  • The following size requirements apply:
    • For all X = C, A, Acc, B,Bcc: X.tstp()==tstp required if X is of type herm_matrix_timestep, and X.nt()>=tstp required if X is of type ‘herm_matrix’,
    • A.size1()>=a1,a2 required, analogous for C, A, Acc, B,Bcc
    • X.ntau() must be equal for all X = C, A, Acc, B,Bcc
  • From Eq. \eqref{langreth-product-1} one can see that in general A and B must be known outside the hermitian domain in order to compute C in the hermitian domain, which is why the hermitian conjugates \(A^\ddagger\) and \(B^\ddagger\) must be provided as arguments.
Example 1:

The following example calculates the second-order self-energy for a general time-dependent interaction \(U(t)\) out of a scalar local Green's function, \(\Sigma(t,t') = U(t) G(t,t') G(t',t) U(t') G(t,t')\). The diagram is factorised in a particle-hole bubble \(\chi\) and a particle-particle bubble:

\begin{equation} \Sigma(t,t') = -i W(t,t') G(t,t'), \,\,\,\, W(t,t')=U(t) \chi(t,t') U(t'), \,\,\,\,\, \chi(t,t') = iG(t,t') G(t',t). \end{equation}

int nt=10;
int ntau=100;
int size1=1;
GREEN G(nt,ntau,size1,FERMION);
GREEN Sigma(nt,ntau,size1,FERMION);
CFUNC U(nt,size1);
// ...
// ... do something to set G and U
// ...
// compute Sigma on all timesteps:
for(int tstp=-1;tstp<=nt;tstp++){
GREEN_TSTP tW(tstp,ntau,size1,BOSON); // used as temporary variable; Note that a bubble of two Fermion GF is bosonic
cntr::Bubble1(tstp,W,G,G); // W(t,t') is set to chi=ii*G(t,t')G(t',t);
// Note: Because G is scalar and have hermitian symmetry, many arguments can be replaced by default. The above extends to
// cntr::Bubble1(tstp,W,0,0,G,G,0,0,G,G,0,0,);
W.left_multiply(tstp,U);
W.right_multiply(tstp,U); // now W is set to W(t,t')=U(t)chi(t,t')U(t')
cntr::Bubble2(tstp,Sigma,G,W); // Sigma(t,t') is set to Sigma=ii*G(t,t')W(t',t);
// Note that by construction, W is hermitian, see above
Sigma.smul(tstp,-1.0); // the final -1 sign.
}
Example 2:

Evaluation of the diagrams \eqref{diagram1} and \eqref{diagram2}, where A and B are two-dimensional, C is one-dimensional, and V is a corresponding contraction. The contraction over internal indices is done by an explicit four-dimensional summation.

int nt=10;
int ntau=100;
int size1=2;
int sizec=1;
GREEN A(nt,ntau,size1,FERMION);
GREEN B(nt,ntau,size1,FERMION);
GREEN CPH(nt,ntau,sizec,BOSON);
GREEN CPP(nt,ntau,sizec,BOSON);
// ...
// ...
// ...
// compute particle-hole bubble (CPH) and particle-particle bubble (CPP) on all timesteps:
for(int tstp=-1;tstp<=nt;tstp++){
GREEN_TSTP ctmp(tstp,ntau,1,BOSON); // temporary
CPH.set_timestep_zero(tstp);
CPP.set_timestep_zero(tstp);
for(int a1=0;a1<size1;a1++){
for(int a2=0;a2<size1;a2++){
for(int b1=0;b1<size1;b1++){
for(int b2=0;b2<size1;b2++){
double v1=get_vertex(a1,b1); // supposed to return V_{a1,b1}
double v2=get_vertex(a2,b2); // supposed to return V_{a2,b2}
cntr::Bubble1(ctmp,0,0,A,A,a1,a2,B,B,b1,b2); // Note the order of the arguments b1,b2!
CPH.incr_timestep(tstp,ctmp,v1*v2);
cntr::Bubble2(ctmp,0,0,A,A,a1,a2,B,B,b1,b2);
CPP.incr_timestep(tstp,ctmp,v1*v2);
}
}
}
}
}

↑ to the top of the page ↑

Convolution

Full convolution

The convolution of two Green's functions A and B is defined by the integral

\begin{equation} \label{convolution001} C(t,t^\prime) \equiv [A\ast B]C(t,t^\prime) = \int_{\mathcal{C}} d \bar{t}\, A(t,\bar{t}) B(\bar{t},t^\prime), \end{equation}

where the product of \(A(t,\bar{t})\) and \(B(\bar{t},t^\prime)\) is understood as a matrix product in orbital indices. The convolution is one of the most basic operation when dealing with nonequilibrium Green's functions. The evaluation of the integral uses Langreth rules, to derive separate integrals for each Keldysh component of C in terms of the Keldysh components od A and B. Precise expressions are given in Ref. [6], Chapter 3.

We provide a number of routines to calculate the generalized convolution, where a one-time function f(t) is sandwiched between A and B:

\begin{equation} \label{convolution002} C(t,t^\prime) = \int_{\mathcal{C}} d \bar{t}\, A(t,\bar{t}) f(\bar{t}) B(\bar{t},t^\prime). \end{equation}

This integral is computed by the routine:

cntr::convolution(GG &C, GG &A, GG &Acc,cntr::function<T> &f, GG &B, GG &Bcc T beta, T dt, int SolveOrder) Obtain \( C(t,t^\prime) = \int_{\mathcal{C}} d \bar{t}\, A(t,\bar{t}) f(\bar{t}) B(\bar{t},t^\prime)\) in the full two-time plane with given two-time objects A and B, and a function f. output: object C
cntr::convolution_timestep(int tstp, GG &C, GG &A, GG &Acc,cntr::function<T> &f, GG &B, GG &Bcc T beta, T dt, int SolveOrder) Obtain \( C(t,t^\prime) = \int_{\mathcal{C}} d \bar{t}\, A(t,\bar{t}) f(\bar{t}) B(\bar{t},t^\prime)\) on timestep tstp of C with given two-time objects A and B, and a function f. output: object C
cntr::convolution_timestep_omp(int omp_num_threads,int tstp, GG &C, GG &A, GG &Acc,cntr::function<T> &f, GG &B, GG &Bcc, T beta, T dt, int SolveOrder) Same as cntr::convolution_timestep, using shared memory parallelization. output: object C
  • C, A, Acc, B, Bcc are of type cntr::herm_matrix
  • f is of type cntr::function. If the argument is omitted, \(f(t)=1\) is assumed.
  • dt is the time-discretisation step \(\Delta t\) on the real-time branch
  • beta is the length of the imaginary-time contour (inverse temperature)
  • For X=A,B, Xcc is the hermitian conjugate of X. If the argument is omitted, hermitian symmetry Xcc=X is assumed.
  • The arguments must satisfy the following size consistency:
    • X.nt()>=tstp for all C, A, Acc, B, Bcc, f
    • X.ntau() equal for all C, A, Acc, B, Bcc
    • X.size1() equal for all C, A, Acc, B, Bcc, f
  • 0 <= SolveOrder <= MAX_SOLVE_ORDER is the integration order (currently, MAX_SOLVE_ORDER=5).
    Use SolveOrder=5 (=MAX_SOLVE_ORDER) if there is no good reason against it.
  • To be able to perform the integration up to that order, the additional size requirements hold:
    • X.ntau()>SolveOrder required for all C, A, Acc, B, Bcc
    • If tstp>-1: X.nt()>=SolveOrder required for all C, A, Acc, B, Bcc, f
  • The second variant convolution_timestep_omp uses shared memory openMP paralelization to distribute the evaluation of \(C(t,t')\) for different time arguments on the timestep tstp over omp_num_threads tasks.

The numerical implementation is detailed in Ref. [6], Chapter 11.

Note on hermitian conjugate: For the hermitian conjugate one finds

\begin{equation} \label{convcc} [C^\ddagger](t,t^\prime) = \int_{\mathcal{C}} d \bar{t}\, B^\ddagger(t,\bar{t}) f(\bar{t})^\dagger A^\ddagger(\bar{t},t^\prime).\end{equation}

In general, therefore, \(C\) is not hermitian. cntr::convolution_timestep computes C only on the hermitian domain. To restore the full function C, one needs \(C^\ddagger\) on hermitian domain. Following Eq.~ \eqref{convcc}, \(C^\ddagger\) is computed bt a second call to cntr::convolution_timestep, storing the result to a separate herm_matrix variable Ccc:

cntr::convolution_timestep(tstp,C,A,Acc,f,B,Bcc,beta,dt,SolveOrder);
cntr::convolution_timestep(tstp,Ccc,Acc,A,fcc,Bcc,B,beta,dt,SolveOrder); // fcc is the function fcc(t)=f(t)^+

Remark: The library provides a number of un-supported routines which calculate only certain Keldysh components for the convolution, as well as variants with different argument lists.

Note on causal time-dependence and time-stepping:

The exact convolution is always causal. The numerical implementation is exactly causal for timesteps tstp>SolveOrder, while the numerical error on timesteps 0<=tstp<=SolveOrder can depend on the input (A, Acc, B, Bcc, f) on all timesteps -1<=tstp<=SolveOrder. More precisely, the result for C on timestep tstp depends on the following input data:

  • If tstp>SolveOrder or tstp==-1:
    Results depends only on timesteps t<=tstp of the input A, Acc, B, Bcc, and on f(t) for t <=tstp. This is called a causal time-dependence.
  • If 0<tstp<=SolveOrder:
    Results depends on timesteps -1<=t<=SolveOrder of the input A, Acc, B, Bcc, and f.
Example:

The example computes the convolution

\begin{equation} C = G_0 \ast \Sigma \ast G_0 \end{equation}

which is part of the dyson series. The convolution is calculated in two steps:

\begin{equation} A = G_0 \ast \Sigma, \,\,\,\, C = A \ast G_0. \end{equation}

Because A is not hermitian, three calls to cntr::convolution_timestep are needed, which calculate both \(A\) and \(A^\ddagger\) on the [hermitian domain]:

int nt=10;
int ntau=20;
int size1=2;
int SolveOrder=5;
double dt=0.01; // problem dependent
double beta=10.0;
GREEN G0(nt,ntau,size1,FERMION);
GREEN Sigma(nt,ntau,size1,FERMION);
GREEN C(nt,ntau,size1,FERMION);
// the temparary A and its herm. conjugate Acc
GREEN A(nt,ntau,size1,FERMION);
GREEN Acc(nt,ntau,size1,FERMION);
// do anything to fill G0 and Sigma
// calculate A and Acc
for(int tstp=-1;tstp<=nt;tstp++){
// because G0 and Sigma are hermitian, some arguments are omitted
cntr::convolution_timestep(tstp,A,G0,Sigma,beta,dt,SolveOrder);
cntr::convolution_timestep(tstp,Acc,Sigma,G0,beta,dt,SolveOrder);
}
// now compute C
for(int tstp=-1;tstp<=nt;tstp++){
cntr::convolution_timestep(tstp,C,A,Acc,G0,G0,beta,dt,SolveOrder);
}
Note
As mentioned in the comment on causal time-dependency, the input A and Acc for C must be determined on all timesteps 0<=tstp<=SolveOrder before C is computed on timesteps 0<=tstp<=SolveOrder.
// ...
int SolveOrder=5;
// INCORRECT:
for(int tstp=-1;tstp<=nt;tstp++){
cntr::convolution_timestep(tstp,A,G0,Sigma,beta,dt,SolveOrder);
cntr::convolution_timestep(tstp,Acc,Sigma,G0,beta,dt,SolveOrder);
// does not work for tstp=0...4, because for tstp=0...5 routine needs input A and Acc on all times 0...5
cntr::convolution_timestep(tstp,C,A,Acc,G0,G0,beta,dt,SolveOrder);
}
// Alternative (CORRECT):
// convolution on matsubara:
tstp=-1;
cntr::convolution_timestep(tstp,A,G0,Sigma,beta,dt,SolveOrder);
cntr::convolution_timestep(tstp,Acc,Sigma,G0,beta,dt,SolveOrder);
cntr::convolution_timestep(tstp,C,A,Acc,G0,G0,beta,dt,SolveOrder);
// convolution on warm-up steps:
for(int tstp=0;tstp<=SolveOrder;tstp++){
cntr::convolution_timestep(tstp,A,G0,Sigma,beta,dt,SolveOrder);
cntr::convolution_timestep(tstp,Acc,Sigma,G0,beta,dt,SolveOrder);
}
for(int tstp=0;tstp<=SolveOrder;tstp++){
cntr::convolution_timestep(tstp,C,A,Acc,G0,G0,beta,dt,SolveOrder);
}
// convolution on all other steps:
for(int tstp=SolveOrder+1;tstp<=nt;tstp++){
cntr::convolution_timestep(tstp,A,G0,Sigma,beta,dt,SolveOrder);
cntr::convolution_timestep(tstp,Acc,Sigma,G0,beta,dt,SolveOrder);
cntr::convolution_timestep(tstp,C,A,Acc,G0,G0,beta,dt,SolveOrder);
}
// of course, it is also correct to first compute A and Acc on the whole time-domain, as in the code above,
// but in real-life applications with self-consistent Green's function simulations, G0 and Sigma
// are typically not known a priory, but are itself only calculated timestep after timestep.
OpenMP parallelization

The variant convolution_timestep_omp uses shared memory openMP paralelization to distribute the evaluation of \(C(t,t')\) for different time arguments of the timestep tstp over omp_num_threads tasks. The implementation is free of any race condition if omp_num_threads is set to the number of threads in the current team.

Convolution density matrix

In typical applications one is often only interested in the density matrix of the convolution integral \eqref{convolution001} or \eqref{convolution002} (see section Density matrix for the definition of the density matrix ). The following routine may be of use:

cntr::convolution_density_matrix(int tstp, Type &M, GG &C, GG &A, GG &Acc,cntr::function<T> &f, GG &B, GG &Bcc,T beta, T dt,int SolveOrder) Obtain density matrix for Green function \( C(t,t^\prime) = \int_{\mathcal{C}} d \bar{t}\, A(t,\bar{t}) f(\bar{t}) B(\bar{t},t^\prime)\) with given two-time objects A,B at time tstp, and write result to M
  • same comments on input A, Acc, B, Bcc as for cntr::convolution_timestep above
  • If A.size1()==1, M can be a scalar std::complex<T> or a complex eigenmatrix cdmatrix
    If A.size1()>1, M must be complex eigenmatrix cdmatrix
  • If M is a complex eigenmatrix, it is automatically resized to size1 \(\times\)size1
  • f is of type cntr::function. If the argument is omitted, \(f(t)=1\) is assumed.
Example:

Calculate the correlation energy for a two-particle interaction from the self-energy:

\begin{equation} E_{int} = \frac{1}{2} \text{tr} \Big[ \text{DensityMatrix}\big(\Sigma \ast G \big)(t)\Big] \end{equation}

int nt=10;
int ntau=20;
int size1=2;
int SolveOrder=5;
double dt=0.01; // time-discretization
double beta=10.0;
GREEN G(nt,ntau,size1,FERMION);
GREEN Sigma(nt,ntau,size1,FERMION);
// calculate G and Sigma on all timesteps
for(int tstp=-1;tstp<=nt;tstp++){
cdmtarix mm;
cntr::convolution_density_matrix(tstp,mm,Sigma,G,beta,dt,SolveOrder);
cout << "Eint at t= " << tstp << " : " << 0.5*mm.trace().real() << endl;
}
Note
The user can also use the function cntr::correlation_energy, which implements the above integral directly.

↑ to the top of the page ↑

Dyson equation

The Dyson equation for the Green's function \(G(t,t^\prime)\) can be written as

\begin{equation} \label{dyson} i\partial_t G(t,t^\prime) + \mu G(t,t^\prime) - \epsilon(t) G(t,t^\prime) - \int_\mathcal{C} d\bar t\, \Sigma(t,\bar t) G(\bar t,t^\prime) = \delta_{\mathcal{C}}(t,t^\prime). \end{equation}

This equation is to be solved for \(G(t,t^\prime)\) for given input \(\epsilon(t)\) and \(\Sigma(t,t^\prime)\), and the KMS boundary conditions. All quantities \(\Sigma(t,t^\prime)\), \(G(t,t^\prime)\), and \(\epsilon(t)\) are square matrices. The equation is an integro-differential form of the Dyson series

\begin{equation} G=G_0+G_0\ast \Sigma\ast G, \end{equation}

where the free Green's function \(G_0\) is determined by the differential equation \( i\partial_t G(t,t^\prime) + \mu G(t,t^\prime) - \epsilon(t) G(t,t^\prime) = \delta_{\mathcal{C}}(t,t^\prime)\) (cf Sec. Free Green's functions). For details further , see [6], Chapter 12.

It is assumed that \(\Sigma=\Sigma^\ddagger\) is hermitian, and \(\epsilon(t) =\epsilon(t)^\dagger\), which implies that also the solution \(G\) possesses hermitian symmetry. Because of the hermitian symmetry, \(G\) can also be determined from the equivalent conjugate equation

\begin{equation} \label{dyson_cc} -i\partial_{t^\prime} G(t,t^\prime) + \mu G(t,t^\prime)- G(t,t^\prime) \epsilon(t^\prime)- \int_\mathcal{C} d\bar t \,G(t,\bar t) \Sigma(\bar t,t^\prime) = \delta_{\mathcal{C}}(t,t^\prime). \end{equation}

Because of its causal nature, the Dyson equation can be solved in a time-stepping manner (see explanation at input/output relation below). The following routines are used to solve the Dyson equation:

void cntr::dyson(GG &G, T mu, cntr::function<T> &H, GG &Sigma, T beta, T dt, int SolveOrder, MAT_METHOD) Solve \eqref{dyson} for G in the full two-time plane with given Sigma, mu, H.
void cntr::dyson_mat(GG &G, T mu, cntr::function<T> &H, GG &Sigma, T beta, int SolveOrder, MAT_METHOD) Solve \eqref{dyson} for G on timestep tstp=-1 with given Sigma, mu, H.
void cntr::dyson_start(herm_matrix<T> &G, T mu, cntr::function<T> &H,GG &Sigma, T beta, T dt, int SolveOrder) Solve \eqref{dyson} for G of type cntr::herm_matrix on all timesteps 0 <= tstp <= SolveOrder with given Sigma, mu, H
void cntr::dyson_timestep(int tstp,GG &G, T mu, cntr::function<T> &H,GG &Sigma, T beta, T dt, int SolveOrder) Solve \eqref{dyson} for G on a timesteps tstp with tstp > SolveOrder with given Sigma, mu, H
void cntr::dyson_timestep_omp(int omp_num_threads, int tstp,GG &G, T mu, cntr::function<T> &H,GG &Sigma, T beta, T dt, int SolveOrder) Same as cntr::dyson_timestep, using shared memory parallelization. output: object G
  • G and Sigma are Green's functions of type cntr::herm_matrix<T>, assumed to be hermitian
  • H is a cntr::function<T>, the function \(\epsilon\) in \eqref{dyson}.
  • SolveOrder \(\in\) 1,...,MAX_SOLVE_ORDER, the order of accuracy for the solution.
    Use SolveOrder=5 (=MAX_SOLVE_ORDER) if there is no good reason against it.
  • dt is the time-discretitation step \(\Delta t\) on the real-time branch
  • beta is the length of the imaginary-time contour (inverse temperature)
  • Size requirements:
    • G.size1()==Sigma.size1()==H.size1(), G.ntau()==Sigma.ntau(), G.ntau() > SolveOrder
    • For dyson_start: G.nt(),Sigma.nt() >= SolveOrder
    • For dyson_timestep: G.nt(),Sigma.nt() >= tstp
  • Special numerical parameters for dyson_mat: MAT_METHOD = CNTR_MAT_FOURIER or CNTR_MAT_FIXPOINT. Different methods to solve the equation on the Matsubara contour (see [6]). CNTR_MAT_FIXPOINT is default and should be used.

The implementation is detailed in [6], Chapter 12.

Input/Output relation and time-stepping
  • dyson_mat: Sigma is read on timestep tstp=-1, H is read on times tstp=-1
    G is written on timestep tstp=-1
  • dyson_start: Sigma is read on timestep tstp=-1,...,SolveOrder, H is read on times tstp=-1,...,SolveOrder, G is read on timestep tstp=-1
    G is written on timestep tstp=0,...,SolveOrder
  • dyson_timestep: Sigma is read on timestep t=-1,...,tstp, H is read on times t=-1,...,tstp, G is read on timestep t=-1,...tstp-1
    G is written on timestep t=tstp

Because of this causal structure, the dyson equation can be solved by time-stepping: First G is determined on timestep tstp=-1 using dyson_mat. The result enters the determination of G on timesteps 0,...,SolveOrder, done with dyson_start. The equation is solved successively for timesteps tstp=SolveOrder+1,SolveOrder+2,..., where the result at timestep tstp depends on all previous timesteps.

Example:

Solution of \eqref{dyson} on all times for given self-energy

int nt=10;
int ntau=20;
int size1=2;
int SolveOrder=5;
double dt=0.01; // time-discretiation
double beta=10.0; // inverse temperature
double mu=1.433; // chemical potential
GREEN G(nt,ntau,size1,FERMION);
GREEN Sigma(nt,ntau,size1,FERMION);
CFUNC H(nt,size1);
//...
// ... do something to set H and Sigma on timestep -1
// solve for G on timestep -1:
cntr::dyson_mat(G,mu,H,Sigma,beta,SolveOrder,CNTR_MAT_FIXPOINT);
// ... do something to set H and Sigma on timesteps 0 ... SolveOrder=5
// solve for G on timesteps 0 ... SolveOrder=5
cntr::dyson_start(G,mu,H,Sigma,beta,dt,SolveOrder);
//all other timesteps:
for(int tstp=SolveOrder+1;tstp<=nt;tstp++){
// ... do something to set H and Sigma on timestep tstp
cntr::dyson_timestep(tstp,G,mu,H,Sigma,beta,dt,SolveOrder);
}
OpenMP parallelization

The variant dyson_timestep_omp uses shared memory openMP parallelization to distribute the evaluation of \(G(t,t')\) for different time arguments of the timestep tstp over omp_num_threads tasks. Note that dyson_timestep_omp and dyson_timestep follow a slightly different implementation (see [6], Section 12.4), thus the result differs by the numerical error.

The openMP parallelization is handled internally in dyson_timestep_omp. A race condition can not occur if omp_num_threads is set to the number of openMP threads in the current team.

↑ to the top of the page ↑

VIE2

Non-singular VIE2

The second important equation is an integral equation of the form

\begin{equation} \label{vie2} G(t,t^\prime) + \int_\mathcal{C} d\bar t\, F(t,\bar t) G(\bar t,t^\prime) = Q(t,t^\prime) \quad \Leftrightarrow\quad (1+F)*G=Q, \end{equation}

This linear equation is to be solved for \(G(t,t^\prime)\) for a given input kernel \(F(t,t^\prime)\), its hermitian conjugate \(F^\ddagger(t,t^\prime)\), and a source term \(Q(t,t^\prime)\). A typical physical application of \eqref{vie2} is the summation of a random phase approximation (RPA) series for a susceptibility \(\chi\),

\begin{equation} \label{eq:rpa_vie2} \chi = \chi_0 + \chi_0\ast V \ast \chi_0 + \chi_0\ast V \ast \chi_0\ast V \ast \chi_0 + \cdots = \chi_0 + \chi_0\ast V \ast \chi. \end{equation}

In the solution of the equation, We assume that both \(Q\) and \(G\) are hermitian. In general, the hermitian symmetry would not hold for an arbitrary input \(F\) and \(Q\). However, it does hold when \(F\) and \(Q\) satisfy the relation \(F\ast Q=Q\ast F^\ddagger\) and \(Q=Q^\ddagger \), which is the case for the typical applications such as the RPA series. In this case, there is an equivalent conjugate equation

\begin{equation} \label{vie2_cc} G(t,t^\prime) + \int_\mathcal{C} d\bar t \,G(t,\bar t) F^\ddagger(\bar t,t^\prime) = Q(t,t^\prime) \quad \Leftrightarrow\quad G*(1+F^\ddagger)=Q. \end{equation}

In the implementation, the solution of Eq \eqref{vie2} is reduced to a Volterra Integral Equation of 2nd kind,which is why the corresponding routines are called vie2.

Because of its causal nature, the vie2 equation can be solved in a time-stepping manner (see explanation of the input/output relation below). The following routines are used to solve the vie2 equation:

void cntr::vie2_mat(GG &G, GG &F, GG &Fcc, GG &Q, T beta, MAT_METHOD, int SolveOrder) Solve \eqref{vie2} for G on timestep tstp=-1 with given two-time objects F and Q
void cntr::vie2_start(GG &G, GG &F, GG &Fcc, GG &Q,T beta,T dt,int SolveOrder) Solve \eqref{vie2} for G on all timesteps 0 <= tstp <= SolveOrder with given two-time objects F and Q
void cntr::vie2_timestep(int tstp,GG &G, GG &F, GG &Fcc, GG &Q,T beta,T dt,int SolveOrder) Solve \eqref{vie2} for G on a timesteps tstp with tstp > SolveOrder with given two-time objects F and Q
void cntr::vie2_timestep_omp(int omp_num_threads,int tstp,GG &G, GG &F, GG &Fcc, GG &Q,T beta,T dt,int SolveOrder) Same as cntr::vie2_timestep, using shared memory parallelization. output: object G
  • The type GG of G, F, Fcc, and Q is cntr::herm_matrix<T> F and Fcc store the hermitian domain of \(F\) and \(F^\ddagger\), respectively.
  • SolveOrder \(\in\) 1,...,MAX_SOLVE_ORDER, the order of accuracy for the solution. Use SolveOrder=5 (=MAX_SOLVE_ORDER) if there is no good reason against it.
  • dt is the time-discretisation step \(\Delta t\) on the real-time branch
  • beta is the length of the imaginary-time contour (inverse temperature)
  • Size requirements:
    • G, F, Fcc, and Q must have the same size1 and ntau
    • For vie2_start: G, F, Fcc, and Q must have nt>= SolveOrder
    • For vie2_timestep: G, F, Fcc, and Q must have nt >= tstp
  • Special numerical parameters for vie2_mat: MAT_METHOD = CNTR_MAT_FOURIER or CNTR_MAT_FIXPOINT. Different methods to solve the equation on the Matsubara contour (see [6]). CNTR_MAT_FIXPOINT is default and should be used.

The implementation is detailed in [6], Chapter 9.

Input/Output relation and time-stepping
  • vie2_mat: F and Q are read on timestep tstp=-1
    G is written on timestep tstp=-1
  • vie2_start: F and Q are read on timestep tstp=-1,...,SolveOrder, G is read on timestep tstp=-1
    G is written on timestep tstp=0,...,SolveOrder
  • vie2_timestep: F and Q are read on timestep t=-1,...,tstp, G is read on timestep t=-1,...tstp-1
    G is written on timestep t=tstp

Because of this causal structure, the vie2 equation can be solved by time-steppig, analogous to dyson: First G is determined on timestep tstp=-1 using vie2_mat. The result enters the determination of G on timesteps 0,...,SolveOrder with vie2_start. The equation is solved successively for timesteps tstp=SolveOrder+1,SolveOrder+2,..., where the result at timestep tstp depends on all previous timesteps.

Example:

Solution of the RPA series \( \chi = \chi_0 + \chi_0\ast W\ast \chi\) for \(\chi\), where \(\chi_0\) is a known (hermitian) two-time function, and \(W\) is a known one-time function (also hermitian, \(W(t)=W(t)^\dagger\)). The equation is cast in the form \eqref{vie2} with \( F (t,t')= -\chi_0(t,t') W (t')\), \( F^\ddagger (t,t')= -W (t)\chi_0(t,t')\), and \(Q=\chi_0\).

int tstp;
int nt=10;
int ntau=20;
int size1=2;
int SolveOrder=5;
double dt=0.01; // time-discretiation
double beta=10.0; // inverse temperature
GREEN chi(nt,ntau,size1,BOSON);
GREEN chi0(nt,ntau,size1,BOSON);
CFUNC W(nt,size1);
GREEN F(nt,ntau,size1,BOSON);
GREEN Fcc(nt,ntau,size1,BOSON);
//
// ... do something to determine chi0 and W on timestep -1
// determine F, Fcc, and solve for chi on timestep -1:
tstp=-1;
F.set_timestep(tstp,chi0);
F.right_multiply(tstp,W,-1.0);
Fcc.set_timestep(tstp,chi0);
Fcc.left_multiply(tstp,W,-1.0);
cntr::vie2_mat(chi,F,Fcc,chi0,beta,SolveOrder,CNTR_MAT_FIXPOINT);
// .... do something to determine chi0 and W on timesteps 0 ... SolveOrder=5
// get F, Fcc, and solve for chi on timesteps 0 ... SolveOrder=5:
for(int tstp=0;tstp<=SolveOrder;tstp++){
F.set_timestep(tstp,chi0);
F.right_multiply(tstp,W,-1.0);
Fcc.set_timestep(tstp,chi0);
Fcc.left_multiply(tstp,W,-1.0);
}
cntr::vie2_start(chi,F,Fcc,chi0,beta,dt,SolveOrder);
//all other timesteps:
for(int tstp=SolveOrder+1;tstp<=nt;tstp++){
// .... do something to determine chi0 and W on timesteps tstp
F.set_timestep(tstp,chi0);
F.right_multiply(tstp,W,-1.0);
Fcc.set_timestep(tstp,chi0);
Fcc.left_multiply(tstp,W,-1.0);
cntr::vie2_timestep(tstp,chi,F,Fcc,chi0,beta,dt,SolveOrder);
}
OpenMP parallelization

The variant vie2_timestep_omp uses shared memory openMP paralelization to distribute the evaluation of \(G(t,t')\) for different time arguments of the timestep tstp over omp_num_threads tasks. Note that vie2_timestep_omp and vie2_timestep follow a slightly different implementation (see [6], Section 13.4), thus the result differs by the numerical error. Similar to dyson_timestep_omp, the implementation prevents any race conditions if omp_num_threads is set to the number of threads in the current team.

Singular VIE2

The singular vie2 is a generalized solver of vie2, where propagators have an instantaneous term, for example

\begin{equation} \label{Gsin} G(t,t^\prime) = \delta_{\mathcal{C}}(t,t^\prime) G^\delta(t) + G^R(t,t^\prime), \end{equation}

where \(G^\delta\) is the instantaneous part and \(G^R(t,t^\prime)\) is the retarded part. The generalization of Eq. \eqref{vie2} to the following example leads to

\begin{equation} \label{vie2sin} (1+F^\delta + F^R)*(G^\delta + G^R)=Q^\delta + Q^R \end{equation}

First, for a given timestep \(t\), we solve this equation for an instantaneous part leading to

\begin{equation} G^\delta(t)= (1+ F^\delta(t))^{-1} Q^\delta(t) . \end{equation}

The retarded part of the vie2 can be rewritten in a form where a routine for the non-singular vie2 can be employed, see VIE2

\begin{equation} (1+[F^\delta + F^R])*G^R= Q^R - F^R * G^\delta, \end{equation}

where \(G^\delta\) is known from the solution of the instantaneous part.

A typical physical usage of the singular version of vie2 is an evaluation of the retarded interaction \(t\) within the GW approximation, see Lattice GW for an example

\begin{equation} [1-V(t) \chi_0(t,t^\prime)]* W(t,t^\prime)= V(t) \delta(t,t^\prime), \end{equation}

where \(V\) is an instantaneous interaction, \(\chi_0\) is a bare bubble and \(W\) is a retarded interaction entering the GW self-energy.

The following routines are used to solve Eq. \eqref{vie2sin}:

void cntr::vie2_timestep_sin(int tstp,GG &G, cntr::function<T> &Gsin, GG &F, GG &Fcc, cntr::function<T> &Fsin,GG &Q, cntr::function<T> &Qsin,T beta,T dt,int SolveOrder)` Solve \eqref{vie2sin} for G for a given timesteps tstp
void cntr::vie2_timestep_sin_omp(int tstp,GG &G, cntr::function<T> &Gsin, GG &F, GG &Fcc, cntr::function<T> &Fsin,GG &Q, cntr::function<T> &Qsin,T beta,T dt,int SolveOrder) Same as cntr::vie2_timestep_sin, using shared memory parallelization.

These routines use the same input as cntr::vie2_timestep, see Non-singular VIE2, with an addition of a singular part of the propagator Gsin presented by cntr::function<T>.

↑ to the top of the page ↑

Utilities

Comparing Green's functions

To compare the data of two Greens functions on a given timestep, one can use

T cntr::distance_norm2(int tstp,GType &A,GType &B) Returns a difference measure \(\Delta[A,B]_{\tt tstp}\) (see below) between A and B on timestep tstp.
  • The type GType of A and B is cntr::herm_matrix<T> or cntr::herm_matrix_timestep<T>, the return type is T (template parameter)
  • Size requirements:
    • A and B must have equal size1 and ntau
    • for X=A,B: If X is herm_matrix, then X.nt()>=tstp is required, if X is herm_matrix_timestep, then X.tstp()==tstp is reqiured.

The difference is defined as the L \(_2\)-norm difference \(|| M ||\) for \(M=A(t,t')-B(t,t')\) of the individual elements, summed over all time-arguments of the timestep:

\begin{align} \label{eq:distnorm} \text{For tstp}=-1:\,\,\, \Delta[A,B] =& \sum_{i=0}^{\tt ntau} ||A^\mathrm{M}(i\Delta\tau)-B^M(i\Delta\tau)|| \\ \text{For tstp}>=0:\,\,\, \Delta[A,B] =& \sum_{i=0}^{\tt tstp} \big( ||A^\mathrm{R}({\tt tstp}\,\Delta t,i\Delta t)-B^\mathrm{R}({\tt tstp}\,\Delta t,i\Delta t)|| + ||A^<(i\Delta t,{\tt tstp}\,\Delta t)-B^<(i\Delta t,{\tt tstp}\,\Delta t)|| \big) \nonumber \\ &+ \sum_{i=0}^{\tt ntau} ||A^{\rceil}({\tt tstp}\,\Delta t,i\Delta\tau)-B^{\rceil}({\tt tstp}\,\Delta t,i\Delta\tau)|| \end{align}

\(|| M ||\) is the standard matrix L \(_2\) norm.

It is also possible to compare the individual components entering Eq. \eqref{eq:distnorm} separately. For \({\tt tstp}>=0\), one obtains the differences of the retarded/lesser/left-mixing components by cntr::distance_norm2_ret/cntr::distance_norm2_les/cntr::distance_norm2_tv. The interface is identical to cntr::distance_norm2.

Example:

A typical application is to check for convergence in self-consistent simulations:

// int nt= ...
// int ntau= ...
// int sig=...
GREEN G(nt,ntau,size1,sig);
// int tstp= ...
// some iterative procedure to determine G on timestep tstp:
{
GREEN_TSTP tG(tstp,ntau,size1,sig); //temporary variable
double convergence_error;
int iter_max=100;
for(int iter=0;iter<=iter_max;iter++){
tG.set_timestep(tstp,G); // store values of G before iteration
// ... some code to update G on timestep tstp ...
convergence_error = cntr::distance_norm2(tstp,G,tG);
if( convergence_error < some_sufficiently_small_number) break;
}
if(iter>iter_max){
cout << "no convergence!" << endl;
}
}

Timestep extrapolation

Extrapolation:
void cntr::extrapolate_timestep(int tstp,Gtype &A,int ExtrapolationOrder) Extrapolate from timesteps t=tstp,tstp-1,...,tstp-ExtrapolationOrder to timestep tstp+1, using Polynomial extrapolation. output: extrapolated object A
  • The type GType of A is cntr::herm_matrix<T> or cntr::function<T>
  • If A is herm_matrix, all entries \(A(t,t')\) on timestep tstp+1 are written.
    If A is function, the value A(tstp+1) is written.
  • Size requirements: tstp>=ExtrapolationOrder and A.nt()>=tstp+1 required
  • The ExtrapolationOrder must be between 1 and MAX_SOLVE_ORDER (=5).
Continuation from Imaginary time to Real time 0:
void cntr::set_t0_from_mat(herm_matrix<T> &A) Sets the values at timestep tstp=0 from tstp=-1 assuming that \(A(t,t')\) is continuous.
  • Size requirements: A.nt()>=0 required

Continuity of \(A(t,t')\) implies that the point 0 on the lower real-time branch and on the imaginary time branch of the contour are identical. This implies timestep 0 can be set from timestep -1 using the relation ( \(\xi\) is the FERMION or BOSON sign):

  • \( G^{\rceil}(0,j\Delta \tau) = \xi G^M(({\tt ntau} - j)\Delta \tau)\) for j=0,...,ntau
  • \( G^{<}(0,0) = G^{\rceil}(0,0)\)
  • \( G^{R}(0,0) = i G^M(0) - G^{<}(0,0)\)
Note
A two-time functions need not be continuous if it explicitly depends on a parameter which changes discontinuously from the imaginary time to real-time contour. For example, the self-energy depends on the interaction U, which changes discontinuously for an interaction quench at time 0.

Differentiation

In order to perform differentiation on the real part of the contour, one can use:

void deriv1_timestep(int tstp, Gtype &dA, Gtype &A, Gtype &Acc, T beta, T h, int SolverOrder) calculates the left derivative dA at time step tstp, i.e. \(dA(t,t') = id/dt A(t,t')\), when A and its complex conjugate Acc are given.
void deriv2_timestep(int tstp, Gtype &dA, Gtype &A, Gtype &Acc, T beta, T h, int SolverOrder, ) calculates the right derivative dA at time step tstp, i.e. \(dA(t,t') = -id/dt' A(t,t')\), when A and its complex conjugate Acc are given.
  • The type Gtype of A, Acc, and dA is cntr::herm_matrix<T>. The return type is T (template parameter). Currently, the routines are implemented only for objects of size1=1
  • Acc is the conjugate function to A (similar as in Convolution). If A is hermitian, simply provide Acc=A.
  • For the computation of timestep tstp, \(A(t,t')\) is addressed at \(t,t' \leq \max(n,k)\), where \(k\) is the solver order. Note that tstp=0..k can be computed only if A is given for \(t,t'\leq k\).
  • the possible contribution proportional to \(\delta(t,t’)\), arrisung from the step discontinuity of a GFs at equal times, is subtracted in these routines.

For more information see [6], Chapter 8.

↑ to the top of the page ↑

Creating and passing input files

Along with the libcntr library, we also provide an input file parser, which can be used in custom programs by including

#include "cntr/utils/read_inputfile.hpp"

Input files allow for passing variables to NESSi-based programs. We use the following format:

__var=value

Note the double underscores in front of the variable name. Let us consider an example input file input.txt,

__nt=1000
__ntau=400
__h=0.01
__beta=10.0

which is parsed by the code snipped

int nt,ntau;
double h,beta;
find_param("input.txt", "__nt=", nt)
find_param("input.txt", "__ntau=", ntau)
find_param("input.txt", "__h=", h)
find_param("input.txt", "__beta=", beta)

Furthermore, the input file parser can read vectors from file. A typical application is a time-dependent parameter such as an external electric field \(E(t)\). We assume the field to be represented on the same grid as in the calculation: \(E_j = E(j h)\), \( j=0,\dots,n_t\). We also include the value in thermal equilibrium \(E_{-1}\) (hence, \( E(t)\) is effectively represented by a cntr::function). Preparing an input file field.txt like

0.0
0.0
0.1
0.2
...

with nt+2 entries, we can read the electric field by referring to field.txt in the input file input.txt:

__Efield=--field.txt

In the C++ program, the field can be read as follows:

std::vector<double> Efield;
read_tvector("input.txt", "__Efield=", Efield, nt);

Resizing to the corresponding size is part of the function read_tvector.

Multidimensional data can also be read. Extending the above example to a two-dimensional electric field, the input file field.txt would look like

0.0 0.0
0.0 0.0
0.1 0.1
0.2 0.2
...

Reading this field in a program

int dim=2;
std::vector<dvector> Efield;
read_tvector("input.txt", "__Efield=", Efield, dim, nt);

To ease the usage of input files, we provide the python module ReadCNTR.py (available for both python2 and python3), which allows to create input files from python dictionaries. For more details see Python tools.

↑ to the top of the page ↑

MPI parallelization

In large-scale applications, one often encounters problems that can be parallelized by distributing the calculation of different Green's functions across different computing ranks. A typical example are lattice simulations, where the solution of the Dyson equation for the Green's function \(G_k\) on each point \(k\) on a discrete momentum grid can be performed in parallel. The library provides some features for distributed-memory parallelization via MPI:

Communicating timesteps

Member functions of A, where A is of type cntr::herm_matrix and cntr::herm_matrix_timestep:

void Bcast_timestep(int tstp, int root) Broadcast the data of \(A(t,t')\) on timestep tstp from rank root to all other ranks.
void Reduce_timestep(int tstp, int root) In place reduction, replacing \(A(t,t')\) on rank root by \( \sum_{\text{rank }\,j} A_{\text{at rank}\,j}(t,t')\) for one timestep tstp.
void Send_timestep(int tstp, int dest, int tag) Send the data of \(A(t,t')\) on timestep tstp to rank dest with a message tage tag.
void Recv_timestep(int tstp, int root, int tag) Receive a message with tag tag from rang root which contains the data of \(A(t,t')\) on timestep tstp.

A simple reduction operation:

cntr::Reduce_timestep(int tstp,int root,GG &A,GG &B) \(A(t,t')\) at root is set to \( \sum_{\text{rank }\,j} B_{\text{at rank}\,j}(t,t')\) for all time arguments on the timestep tstp
  • The type GG of A and B is cntr::herm_matrix<T> or cntr::herm_matrix_timestep<T>
  • B must exist on each rank, A on the root rank.
  • Size requirements:
    • A and B must have equal size1 and ntau
    • for X=A,B: If X is herm_matrix, then X.nt()>=tstp is required, if X is herm_matrix_timestep, then X.tstp()==tstp is reqiured.
Example:

In the following example, each rank owns the Green's function \(G_k\) for one momentum \(k\) out of a given discrete grid of \(N_k\) momenta, which are equidistantly spaced over the first Brillouin zone \([0,2\pi)\) of a one-dimensional lattice. Rank 0 owns a local Self-energy, which is send to all ranks to solve the Dyson equation

\begin{equation} i\partial_t G_k(t,t^\prime) + \mu G_k(t,t^\prime) - \epsilon_k(t) G_k(t,t^\prime) - \int_\mathcal{C} d\bar t\, \Sigma_{loc}(t,\bar t) G_k(\bar t,t^\prime) = \delta_{\mathcal{C}}(t,t^\prime). \end{equation}

The data are then collected to compute a local Green's function \(G_{loc} = \sum_{k} w_k G_k\) on rank 0, with some given weights w_k (which are just w_k=1/N_k in the example). (In real DMFT applications, e.g., \(\Sigma_{loc}\) is the re-calculated from \(G_{loc}\), and the process is iterated to convergence)

// MPI is initialized as:
MPI::Init(argc,argv);
int nranks=MPI::COMM_WORLD.Get_size();
int tid=MPI::COMM_WORLD.Get_rank();
int root=0;
int nk=ranks; // Here : # of momentum points = # of MPI ranks
double local_k = (2.0*M_PI/nk)*tid; // momentum kept at the local rank
double wk=1.0/nk;
// allocate Gk and epsilonk with same size nt,ntau,size1,sig on all ranks
GREEN Gk(nt,ntau,size1,sig);
CFUNC epsilonk(nk,size1);
GREEN Sigma,Gloc;
// the variable Sigma and Gloc is allocated only on root:
if(tid==root){
Sigma=GREEN(nt,ntau,size1,sig);
Gloc=GREEN(nt,ntau,size1,sig);
}
// ... do something to initialize function epsilonk et every rank, depending on local_k
// typical simulation loop at timestep tstp (tstp>SolveOrder):
// send Sigma to all ranks:
Sigma.Bcast_timestep(tstp,root);
// Solve dyson equation separately on each rank:
cntr::dyson_timestep(tstp,Gk,mu,H,Sigma,SolveOrder,beta,dt);
// compute the local Green's function at the root, using a temporary variable:
GREEN_TSTP tGk(tstp,ntau,size1,sig);
tGk.set_timestep(tstp,Gk);
tGk.smul(tstp,wk);
cntr::Reduce_timestep(tstp,root,Gloc,tGk);

Distributed timestep array

While the previous example was based on simple root to all reduction and broadcasting operations, calculating lattice susceptibilities and self-energies often requires an all-to-all communication of a set of Green's functions at one timestep. We provide a class cntr::distributed_timestep_array which is customized for this application:

class cntr::distributed_timestep_array<T>
  • The cntr::distributed_timestep_array<T> contains an array T_k of n herm_matrix_timesteps (each with the same tstp, ntau, size1 and sig parameters, see Timeslices), which are indexed by an index k \(\in\{\)0,...,n-1 \(\}\). k can be, e.g., a momentum label for lattice simulations.
  • The full array is accessible at each MPI rank, but each timestep T_k is owned by precisely one rank (a vector tid_map[k] stores the id of the rank which owns timestep k). Typically, each rank first performs local operations on the timesteps T_k which it owns. The all timesteps are broadcasted from the owning rank to all other ranks, so that finally each timestep T_k is known to all ranks.

Important member functions of cntr::distributed_timestep_array<T>:

distributed_timestep_array(int n,int nt,int ntau,int size,int sig,bool mpi) Constructor. ntau, size and sig are parameters of the individual timesteps. They cannot be changed later. The tstp parameter of the timesteps can later be adjusted (see reset_tstp) up to a maximal value of nt. mpi should be true (otherwise each rank owns each timestep, and mpi operations are ignored). The constructor automatically determines the ownership map tid_map (as evenly as possible). Different distributed_timestep_array objects with the same parameters have the same tid_map.
int nt() returns nt. Analogous for size(), ntau(), tstp(), sig(), n()
void reset_tstp(int tstp) The tstp parameter of all timesteps is reset to tstp. tstp<=nt required.
int tid(), int ntasks() returns rank id (tid) of the local process and number of MPI ranks (ntasks)
void clear(void) All data set to 0.
bool rank_owns(int k) return true if local rank owns timestep with index k, return false otherwise. ‘0 <= k < n’ required'
cntr::herm_matrix_timestep_view<T> &G(int k) Return a handle to the timestep of index k. See Note below.
void mpi_bcast_block(int k) Broadcast timestep \(T_k\) from its owner to all other ranks.
void mpi_bcast_all(void) Broadcast all timesteps from their owner to all other ranks.
std::vector<int> data().tid_map() A way to return the vector tid_map, where tid_map[k] for 0 <= k < n is the rank which owns timestep T_k.
data() Returns a handle to an underlying cntr::distributed_array on which the distributed_timestep_array is built. Not needed in standard applications.
Note
In the implementation, the distributed_timestep_array is not an array of cntr::herm_matrix_timestep<T> variables, but a sufficiently large contiguous data block and an array of shallow cntr::herm_matrix_timestep_view<T> objects. This allows for an easier adjustment of the size of the timesteps. The herm_matrix_timestep_view<T> object has the same functionality as the herm_matrix_timestep (see Overview)
Example:

The usage of the distributed_timestep_array can be understood best by means of an example. In the example below we consider a vector (std::vector<GREEN> Gk) of momentum-dependent Green's functions \(G_k\) on a grid of nk=6 momenta (k=0,...,5). The memory-intensive Green's functions should be distributed on ntasks=3 MPI ranks, such that each \(G_k\) is available only on one rank. One can use a distributed_timestep_array object TARRAY to make a given timestep of all Gk available to each rank, as illustrated in the Figure:

If you don't see it, something is missing.
  • The distributed_timestep_array specifies which k is owned by which rank.
  • The vector Gk is available at each rank, but at a given rank memory for Gk[k] is only allocated if k is owned by that rank (i.e., if TARRAY.rank_owns(k)==true).
  • Each rank can perform local operations on its own Green's functions Gk[k], and then copy the data of a given timestep tstp to the TARRAY (blue arrows in figure)
  • The global MPI operation TARRAY.mpi_bcast_all(); then makes all timesteps visible to all ranks. One can then perform any operation, such as averaging the Green's function on one timestep (in the example below) or computing a momentum dependent self-energy or polarization diagram at one timestep which requires an internal momentum summation.
    // MPI initialized with ntasks=6 ranks, local rank has rank-ID tid
    int nk=6;
    // ... set nt,ntau,size1,sig ...
    GREEN Gloc(nt,ntau,size1,sig);
    std::vector<GREEN> Gk(nk);
    cntr::distributed_timestep_array<double> TARRAY(nk,nt,ntau,size1,FERMION,true);
    for(int k=0;k<nk;k++){
    if(TARRAY.rank_owns(k)){
    cout << "rank " << tid << " owns k= " << k << endl;
    Gk[k]=GREEN(nt,ntau,size1,sig); // allocate memory for full Green's function Gk only on ranks which own k
    }else{
    cout << "rank " << tid << " does not own k= " << k << endl;
    }
    }
    // typical simulations at a given timestep:
    TARRAY.reset_tstp(tstp);
    for(int k=0;k<nk;k++){
    if(TARRAY.rank_owns(k)){
    // ... rank owns kdo some heavy numerics on Gk[k]
    // then copy timestep to TARRAY:
    TARRAY.G(k).set_timestep(tstp,Gk[k]); // read in data from Gk[k] to the array
    }
    }
    TARRAY.mpi_bcast_all();
    // now on all ranks and for all k TARRAY.G(k) contains the data of Gk[k]
    // this can be used, e.g., to calculate a local Green's function:
    Gloc.set_timestep_zero(tstp);
    for(int k=0;k<nk;k++) Gloc.incr_timestep(tstp,gk_timesteps.G(k),1.0/nk);
    // or do more complicated stuff, such as calculating a self-energy diagram with internal k summations

Distributed array

class cntr::distributed_array<T>
  • The cntr::distributed_array<T> contains an array of n equally sized data blocks, each containing blocksize data of type T, which are indexed by an index k \(\in\{\)0,...,n-1 \(\}\). All data are stored contiguous in memory.
  • The full data set is accessible at each MPI rank, but each block is owned by precisely one rank (a vector tid_map[k] stores the id of the rank which owns timestep k). MPI routines are used send the blocks between ranks.

Important member functions of cntr::distributed_array<T>:

distributed_array(int n,int maxlen,bool mpi) Constructor. maxlen sets the maximum block size. mpi should be true (otherwise each rank owns each timestep, and mpi operations are ignored).
void reset_blocksize(int blocksize) The number of elements in each block is set to blocksize. blocksize<=maxlen required.
int n() returns n. Analogous for maxlen(), blocksize()
int tid(), int ntasks() returns rank id (tid) of the local process and number of MPI ranks (ntasks)
void clear(void) All data set to 0.
bool rank_owns(int k) return true if local rank owns timestep with index k, return false otherwise. ‘0 <= k < n’ required'
std::vector<int> tid_map() Return the vector tid_map, where tid_map[k] for 0 <= k < n is the rank which owns block k.
T* block(int k) Return a pointer to the first element of block k. 0 <= k < n required.
int numblock_rank(void) Return the number of blocks owned by the local rank.
void mpi_bcast_block(int k) Broadcast data of block k from its owner to all other ranks.
void mpi_bcast_all(void) Broadcast the data all block from their owner to all other ranks.
void mpi_send_block(int k,int dest) Send Block k from its owner to rank dest.
void mpi_gather(int dest) MPI Gather all blocks from their owner at rank dest

↑ to the top of the page ↑

hdf5 usage

Overview

libcntr uses the hdf5 to store the basic data types for contour functions to disk. hdf5 is an open source library and file format for numerical data which is widely used in the field of scientific computing. The format has two building blocks:

  • data sets : general multi-dimensional arrays of a single type
  • groups : containers which can hold data sets and other groups.

Hence, by nesting groups, it is possible to store arbitrarily complicated structured data, and to create a file-system-like hierarchy where groups can be indexed using standard POSIX format, e.g. /path/to/data.

The libcntr library comes with helper functions to store the basic contour response function data types in hdf5 with a predefined structure of groups and data sets, defined in the header cntr/hdf5/hdf5_interface.hpp. For example a herm_matrix response function is stored as a group with a data set for each contour component mat ( \(g^M(\tau)\)), ret ( \(g^R(t, t')\) ), les ( \(g^<(t, t')\) ), and tv ( \(g^\rceil(t, \tau)\) ), respectively, see Section Green functions. The retarded and lesser components are stored in upper and lower triangular contiguous time order respectively. In the libcntr hdf5 format each component is stored as a rank 3 array where the first index is time, imaginary time, or triangular contiguous two-time, and the remaining two indices are orbital indices.

Reading/writing to hdf5 files

To store a contour GF of type cntr::herm_matrix, one writes its components into a group of a hdf5 file using the member function write_to_hdf5.

Example

In C++ this takes the form,

#include <cntr/cntr.hpp>
..
// Create a contour Green's function
int nt = 200, ntau = 400, norb = 1;
GREEN A(nt, ntau, norb, FERMION);
// Open HDF5 file and write components of the Green's function A into a group g.
std::string filename = "data.h5";
A.write_to_hdf5(filename.c_str(), "g");

If the file data.h5 has been written previously with write_to_hdf5, one can read it with the member function read_from_hdf5

// Open HDF5 file and read group g. The result is saved into the Green's function B
GREEN B;
B.read_from_hdf5(filename.c_str(), "g");

the parameters (nt,ntau,size1,sig) and the data of B are modified according to the information in the file (Similar to reading/writing to text files discussed in File I/O).

To understand the structure of the resulting hdf5 file data.h5 we inspect it with the h5ls command line program that can be used to list all groups and data sets in a hdf5 file,

$ h5ls -r data.h5
...
/g Group
/g/element_size Dataset {1}
/g/les Dataset {20301, 1, 1}
/g/mat Dataset {401, 1, 1}
/g/nt Dataset {1}
/g/ntau Dataset {1}
/g/ret Dataset {20301, 1, 1}
/g/sig Dataset {1}
/g/size1 Dataset {1}
/g/size2 Dataset {1}
/g/tv Dataset {80601, 1, 1}

where we see that apart from the contour components the Green's function group g contains additional information about the dimensions and the Fermi/Bose statistics (sig \( = \mp 1\)), for details see the API documentation of cntr::herm_matrix and Section Green functions.

To understand the dimensions of the contour components we can look at the number of imaginary time steps ntau and number of real time steps nt using the h5dump command line utility,

$ h5dump -d /g/ntau data.h5
HDF5 "data.h5" {
DATASET "/g/ntau" {
DATATYPE H5T_STD_I32LE
DATASPACE SIMPLE { ( 1 ) / ( 1 ) }
DATA {
(0): 400
}
}
}
$ h5dump -d /g/nt data.h5
HDF5 "data.h5" {
DATASET "/g/nt" {
DATATYPE H5T_STD_I32LE
DATASPACE SIMPLE { ( 1 ) / ( 1 ) }
DATA {
(0): 200
}
}
}

which shows that the dimensions are \(n_\tau = 400\) and \(n_t=200\). The size of the /g/mat component reveals that this corresponds to \(n_\tau + 1 = 401\) imaginary time points. The mixed /g/tv component has a slow time index and a fast imaginary time index and is of size \((n_t + 1)(n_\tau + 1) = 80601\) while the two time triangular storage of the /g/ret and /g/les components contains \((n_t + 1)(n_t + 2)/2 = 20301\) elements.

To simplify postprocessing of contour GFs NESSi also provides the python module ReadCNTRhdf5.py for reading the hdf5 format (using the python modules numpy and h5py) producing python objects with the contour components as members. For details see reading Green's functions from hdf5

↑ to the top of the page ↑

Python tools

As a part of the NESSi package, we provide python tools, which include:

  • scripts for pre-processing to assist the use of programs based on libcntr
  • scripts for reading and post-processing Green's functions from hdf5 format via the h5py python package.

The python tools can be found in libcntr/python (for python2) or libcntr/python3 (for python3 compatitability).

creation of an input file

An input file for a custom program (see details in Creating and passing input files) can be generated using ReadCNTR.py python module:

Example
from ReadCNTR import write_input_file
inp = {
'nt': 1000, # number of points on the real axis
'ntau': 500, # number of points on the Matsubara (imaginary) axis
'h': 0.01, # timestep interval
'beta': 10.0, # inverse temperature
'Efield': '--field.txt' # electrical field readed from a file
}
write_input_file(inp, 'input.txt') # creates an input file
Note
To make the python module available to a python script, the module ReadCNTR.py needs to be in the python include path. For python2, this is achieved by setting export PYTHONPATH=<path>/nessi/libcntr/python, while export PYTHONPATH=<path>/nessi/libcntr/python3 makes the python3 verion accessible.

reading Green's functions from hdf5

The python module unrolls the triangular storage of the ret and les components making it simple to work with time slices.

Example 1

To store the imaginary part of the retarded Green's function \(\textrm{Im}[G_{i,j}^\mathrm{R}(t, t^\prime=5 \text{h})]\) at orbital indices \(i=j=0\) as a function of \(t\) into a text file we may use the commands

import h5py
from ReadCNTRhdf5 import read_group
# read the components of the contour function from hdf5 file
with h5py.File('data.h5', 'r') as fd:
g = read_group(fd).g
# choose a timeslice for the retarded Green's function
indxt=5
G_ret=g.ret[:,indxt,0,0].imag
with open("data_output.txt",'w') as output: # create a text file
output.write("%s \t" % G_ret) # write data into it
Example 2

The same way one can save the lesser component of the Green's function \(\textrm{Im}[G^<(t=5 \text{h}, t^\prime)]\) as a function of \(t^\prime\) into a text file

import h5py
from ReadCNTRhdf5 import read_group
# read the components of the contour function from hdf5 file
with h5py.File('data.h5', 'r') as fd:
g = read_group(fd).g
# choose a timeslice for the lesser Green's function
indxt=5
G_les=g.les[indxt,:,0,0].imag
with open("data_output.txt",'w') as output: # create a text file
output.write("%s \t" % G_les) # write data into it

postprocessing analysis

The ReadCNTRhdf5 python module contains also function for reading parameters and observables from an output file. The following minimal example script is used to interpret the output files in the examples provided in section Example programs.

Example
import h5py
from ReadCNTRhdf5 import read_imp_h5file
# Plot
ax0=plt.subplot(1,1,1)
imp_filename = '{}/data.h5'.format(output_file)
data = read_imp_h5file(imp_filename,['obs','parm'])
Nt=data.parm.Nt # read out the number of timesteps on the real axis from the output file
dt=data.parm.h # read out the timestep interval from the output file
tpts = np.linspace(0.0,Nt*dt,Nt+1) # define the real time axis
ekin=np.real(data.obs.Ekin.data[1:,0,0]-data.obs.Ekin.data[0,0,0]) # read out the value of an observable (Ekin) from the output file and calculate the difference between the equilibrium (Ekin.data[0,0,0]) and nonequlibrium values (Ekin.data[1:,0,0])
plt.plot(tpts,ekin,label="Ekin") # plot data
ax0.set_xlabel(r't ') # set x label
ax0.set_ylabel(r"$\mathrm{E}_\mathrm{kin}") # set y label
ax0.legend(loc="best",bbox_to_anchor=(0.45,0.3),frameon=False,fancybox=False,handlelength=1)
plt.savefig('Ekin.pdf')
plt.show()

↑ to the top of the page ↑