NESSi
v1.0.2
The NonEquilibrium Systems Simulation Library
|
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). |
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:
Green's functions \(C(t,t')\) can be matrix-valued. We use the eigen3
library for handling matrices.
The following shortcuts are defined:
(For a complete documentation, see eigen3
)
Listed for reference only. Explanation in the sections below.
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 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):
ntau
,nt
, j=0,...,i ,nt
, i=0,...,j ,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:
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):
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.
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 . |
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
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() |
C.size1()>1
, M
must be a complex eigen matrix (cdmatrix
for double precision)C.size1()==1
, M
can be a scalar (cdouble
) or a matrixM
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() |
M
is a matrix, it is resized to a square matrix of dimension C.size1()
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.
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\).Ccc
is omitted, hermitian symmetry \(C=C^\ddagger\) is assumed.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.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).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.C.size1()>1
, M
must be a complex eigen matrix (cdmatrix
for double precision)C.size1()==1
, M
can be a scalar (cdouble
) or a matrixM
is a matrix, it is resized to a square matrix of dimension C.size1()
.The member functions print_to_file
and read_to_file
of a herm_matrix
implement text-file access, as apparent by example:
The file format is rather explicit (and storage intensive):
# nt ntau size1 sig
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}
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}
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.hdf5
format below.The hdf5
format, as well as some scripts to interpret the hdf5
files are discussed in a separate section hdf5 usage below.
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:
tstp=-1
, the timestep contains the Matsubara Green's function:ntau
,tstp>=0
, it contains the following real-time components:j=0,...,tstp
,j=0,...,tstp
,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 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.
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 . |
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.
simple print of all elements in case of a scalar Green's function, at tstp=-1 (Matsubara):
Remark: Some un-supported routines exist in which the dummy argument tstp
can be omitted.
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>
tstp
(integer): The timestep to be accessedtstp==-1
this copies the data from B
to A
as:i=0,...,ntau
tstp>=0
this copies the data from B to A as:j=0,...,tstp
j=0,...,tstp
j=0,...,ntau
A
is of type herm_matrix_timestep
, tstp==A.tstp()
is required, same for B
A
is of type herm_matrix
, tstp<=A.nt()
is required, same for B
A.ntau()==B.ntau()
and A.size1()==B.size1()
is requiredThe following code copies the first 5 timesteps and the Matsubara timestep from a Green's function B
to a Green's function A
:
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 accessedi1
, i2
, j1
, j2
(integer): Orbital matrix indiceststp==-1
this copies the data from B
to A
as:i=0,...,ntau
tstp>=0
this copies the data from B to A as:i=0,...,tstp
i=0,...,tstp
i=0,...,ntau
A
is of type herm_matrix_timestep
, tstp==A.tstp()
is required, same for B
A
is of type herm_matrix
, tstp<=A.nt()
is required, same for B
A.ntau()==B.ntau()
is required0 <= 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.
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.
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 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
).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 . |
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.
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 |
class | cntr::function<T> |
A contour function is a matrix or scalar-valued function \(f(t)\) which depends on physical time only:
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 nt
(integer): number of discretization points on the real time axis. 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
.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 |
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 |
f.size1()>1
, M
must be a square matrix (cdmatrix
for double
precision)f.size1()==1
, M
can be also a scalar (cdouble
) or a square matrixM
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()
M
is a scalar, only the (0,0) entry of \(f(t)\) is read.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.
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>
A
is of type herm_matrix_timestep
, then A.tstp()==tstp
is required (same for B
) A
is of type herm_matrix
, then A.nt()>=tstp
is required (same for B
)B.ntau()==A.ntau()
and B.size1()==A.size1()
required.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}
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
)A
is herm_matrix_timestep
, then tstp==A.tstp()
is required A
is herm_matrix
, then tstp<=A.nt()
is requiredtstp<=f.nt()
and f.size1()==A.size1()
is required.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.
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.
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):
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
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:
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
.
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 |
GG
can be herm_matrix_timestep<T>
or herm_matrix<T>
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:
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:
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:
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 |
The basic building blocks of Feynman diagrams are particle-hole bubbles and particle-particle bubbles, such as:
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
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\).(c1,c2)
, (a1,a2)
, or (b1,b2)
are omitted, they are assumed to be (0,0)
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
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.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}
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.
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 branchbeta
is the length of the imaginary-time contour (inverse temperature)X=A,B
, Xcc
is the hermitian conjugate of X
. If the argument is omitted, hermitian symmetry Xcc=X
is assumed.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
). SolveOrder=5
(=MAX_SOLVE_ORDER
) if there is no good reason against it.X.ntau()>SolveOrder
required for all C
, A
, Acc
, B
, Bcc
tstp>-1
: X.nt()>=SolveOrder
required for all C
, A
, Acc
, B
, Bcc
, f
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
:
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:
tstp>SolveOrder
or tstp==-1
:t<=tstp
of the input A
, Acc
, B
, Bcc
, and on f(t)
for t <=tstp
. This is called a causal time-dependence.0<tstp<=SolveOrder
:-1<=t<=SolveOrder
of the input A
, Acc
, B
, Bcc
, and f
.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]:
A
and Acc
for C
must be determined on all timesteps 0<=tstp<=SolveOrder
before C
is computed on timesteps 0<=tstp<=SolveOrder
.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.
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 |
A
, Acc
, B
, Bcc
as for cntr::convolution_timestep
aboveA.size1()==1
, M
can be a scalar std::complex<T> or a complex eigenmatrix cdmatrix
A.size1()>1
, M
must be complex eigenmatrix cdmatrix
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.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}
cntr::correlation_energy
, which implements the above integral directly.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 hermitianH
is a cntr::function<T>
, the function \(\epsilon\) in \eqref{dyson}.SolveOrder
\(\in\) 1,...,MAX_SOLVE_ORDER
, the order of accuracy for the solution. 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 branchbeta
is the length of the imaginary-time contour (inverse temperature)G.size1()==Sigma.size1()==H.size1()
, G.ntau()==Sigma.ntau()
, G.ntau() > SolveOrder
dyson_start
: G.nt(),Sigma.nt() >= SolveOrder
dyson_timestep
: G.nt(),Sigma.nt() >= tstp
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.
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.
Solution of \eqref{dyson} on all times for given self-energy
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.
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 |
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 branchbeta
is the length of the imaginary-time contour (inverse temperature)G
, F
, Fcc
, and Q
must have the same size1
and ntau
vie2_start
: G
, F
, Fcc
, and Q
must have nt>= SolveOrder
vie2_timestep
: G
, F
, Fcc
, and Q
must have nt >= tstp
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.
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.
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\).
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.
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 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 . |
A
and B
is cntr::herm_matrix<T>
or cntr::herm_matrix_timestep<T>
, the return type is T
(template parameter)A
and B
must have equal size1
and ntau
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
.
A typical application is to check for convergence in self-consistent simulations:
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 |
A
is cntr::herm_matrix<T>
or cntr::function<T>
A
is herm_matrix
, all entries \(A(t,t')\) on timestep tstp+1
are written. A
is function
, the value A(tstp+1)
is written.tstp>=ExtrapolationOrder
and A.nt()>=tstp+1
requiredExtrapolationOrder
must be between 1
and MAX_SOLVE_ORDER
(=5
).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. |
A.nt()>=0
requiredContinuity 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):
j=0,...,ntau
U
, which changes discontinuously for an interaction quench at time 0.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. |
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=1Acc
is the conjugate function to A
(similar as in Convolution). If A
is hermitian, simply provide Acc=A
.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\).For more information see [6], Chapter 8.
Along with the libcntr
library, we also provide an input file parser, which can be used in custom programs by including
Input files allow for passing variables to NESSi
-based programs. We use the following format:
Note the double underscores in front of the variable name. Let us consider an example input file input.txt
,
which is parsed by the code snipped
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
with nt+2
entries, we can read the electric field by referring to field.txt
in the input file input.txt
:
In the C++ program, the field can be read as follows:
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
Reading this field in a program
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.
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
:
cntr::herm_matrix
and cntr::herm_matrix_timestep
allow sending/receiving the data at one timestepcntr::distributed_timestep_array
can be used to handle all-to-all communication of Green's function on a given timestep.cntr::distributed_timestep_array
is derived from a class distributed_array
, which implements an array of data with distributed ownership.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 |
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.A
and B
must have equal size1
and ntau
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.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)
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> |
cntr::distributed_timestep_array<T>
contains an array T_k
of n
herm_matrix_timestep
s (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.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. |
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)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:
distributed_timestep_array
specifies which k
is owned by which rank.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
).Gk[k]
, and then copy the data of a given timestep tstp
to the TARRAY
(blue arrows in figure)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. class | cntr::distributed_array<T> |
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.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 |
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:
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.
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
.
In C++ this takes the form,
If the file data.h5
has been written previously with write_to_hdf5
, one can read it with the member function read_from_hdf5
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,
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,
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
As a part of the NESSi
package, we provide python tools, which include:
libcntr
hdf5
format via the h5py
python package.The python tools can be found in libcntr/python
(for python2) or libcntr/python3
(for python3 compatitability).
An input file for a custom program (see details in Creating and passing input files) can be generated using ReadCNTR.py
python 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.The python module unrolls the triangular storage of the ret
and les
components making it simple to work with time slices.
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
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
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.