Autodocs Operator

template<class T = danceq::internal::BasisU1<danceq::internal::ContainerTable<danceq::internal::State<128, 2, uint64_t>>>, class ScalarType = double, OperatorType OT = Hamiltonian> class Operator

The Operator class creates an operator based on the Basis class.

Template parameters

class T: It can be either BasisU1 or State class. If T refers to the BasisU1, the operator conserves the total particle number determined by input. However, if T refers to the State class the operator does not conserve the total particle number. There are two different Constructors made for both cases:
// Definition of the Operator class with number conservation using Operator_U1 = danceq::internal::Operator<Basis,std::complex<double>>; Operator_U1 H(L,n); // Definition of the Operator class without number conservation using Operator_NoU1 = danceq::inernal::Operator<State,std::complex<double>>; Operator_NoU1 H_full(L);
Default: danceq:: BasisU1 <danceq:: ContainerTable <danceq:: State <128,2,uint64_t>>>
class ScalarType: The primitive data type that is used. The class allows complex types, e.g., std::complex<double>.

Default: double
class OT: Defines the OperatorType to distinguish between Hamiltonian and Lindbladian systems.

Default: Hamiltonian

Things to know

Terms can be added with add_operator(…) or with add_jump_operator(…). Please follow the link(…) for more details.
New local operators that are not included as default operators can be added with define_operator(…).
The underlying State class can be accessed via Operator::state_class.
The class has full control over the local product space basis using the following functions: get_index(…), get_state(…), and increment_state(…). They are all defined regardless of the underlying basis T and OperatorType OT.
The Lindbladian OperatorType acts on the space of vectorized density matrices. The trace over such a density matrix can be evaluates using a Hamiltonian operator acting on the Hilbert space with trace_over_density_matrix(rho_vectorized).
The debugging level can be adjusted from 0 (no debugging) to 10 (maximal debugging):
#define dbug_level 10; // 10 maximal debugging, 0 is no debugging

Matrix Formats

create_DenseMatrix(): Dense matrix using two-dimensional std::vector
create_EigenDenseMatrix(): Dense matrix using Eigen
create_SparseMatrix(): Sparse matrix using the SparseMatrix class
create_PetscSparseMatrix(): Sparse matrix using Petsc (MPI)
create_PetscShellMatrix(): Matrix-free shell using Petsc (MPI)

Warning

The system size L has to be smaller or equal to the maximal number of sites MaxSites of the underlying State class. For Lindbladian systems, MaxSites has to be at least twice as big!

Warning

Internally, we inverted (\(i\rightarrow L-1-i\)) the sites on which an operator acts. This is due to the fact that the tensor product basis (defined on \((Q-1)^L\)) is revered compared to integer counting. This means that the tensor product \(S^x\otimes S^x\) has a basis structure (\(\vert 0;0\rangle\), \(\vert 0;1\rangle\), \(\vert 1;0\rangle\), \(\vert 1;1\rangle\)) which is inverted to the way how an integer is assigned.

Public Types

using state_class = typename T::state_class

State class.

State class can be retrieved from the Operator class:

// Definition of classes
using State = danceq::internal::State<256,4>;
using Operator = danceq::internal::Operator<State>;

// This equivalent to: State state;
Operator::state_class state;  

using scalartype = ScalarType

ScalarType.

The ScalarType can be accessed via:

typename Operator<T,ScalarType,OT>::scalartype

using scalartype_real = decltype(std::real(ScalarType()))

Real part of ScalarType.

The real ScalarType can be accessed via:

Operator<T,ScalarType>::scalartype_real

If ScalarType is std::complex<double> scalartype_real refers to double.

Public Functions

Operator(const uint64_t L_, const uint64_t n_)

Constructor with T as a Basis and given system size and particle number.

Initiates an Operator with an instance of the BasisU1 with system size L_ and particle number n_ using the Basis constructor without defining the number of subsystems. Template T has to refer to the Basis_class otherwise the constructor can not be used.

It calls add_default_operators() to add the default operators to available operators. It sets delete_basisptr to true such that the created basisptr is deleted.

Parameters:

L_ – System size
n_ – Particle number

Operator(const T &basis)

Constructor with T as a Basis.

Initiates an Operator from a given Basis_class. Template T has to refer to the Basis_class otherwise the constructor can not be used.

It calls add_default_operators() to add the default operators to available operators. It sets delete_basisptr to false such that the created basisptr is not deleted.

Parameters:: basis – Reference to the basis

Operator(const uint64_t L_)

Constructor with T as State class.

Initiates an Operator for a given length L_. Template T has to refer to State class otherwise the constructor can not be used.

It calls add_default_operators() to add the default operators to available operators.

Parameters:: L_ – System size

~Operator()

Deconstructor.

Deletes the basisptr if delete_basisptr is true.

std::vector<std::string> get_available_operators_list(void) const

Returns available local operators.

This are the keys of available_operators.

Returns:: list of operators

int32_t info(void) const

Prints all information of the class.

Prints various information such as dimension, the available_operators, and all input terms from add_operator(…).

Returns:: error_code

uint64_t get_dim(void) const

Returns dimension.

Returns:: dimension

int32_t add_default_operators(void)

Adds some local default operators.

The function adds the following local operator to available_operators: {"Sx","Sy","Sz","S+","S-","n","R","Rdag","P"+std::to_string(k)} for \(k=0,\dots,\textbf{Q}-1`\).

For spin operators, \(m,m^\prime =-S,-S+1,\dots,S\) with \(S = (\textbf{Q}-1)/2\):

\[\begin{split}\langle m^\prime \vert S^x \vert m \rangle &= \left( \delta_{m^\prime,m+1} + \delta_{m^\prime,m+1}\right)\frac{1}{2}\sqrt{S(S+1)-m^\prime m}\\ \langle m^\prime \vert S^y \vert m \rangle &= \left( \delta_{m^\prime,m+1} - \delta_{m^\prime,m+1}\right)\frac{1}{2i}\sqrt{S(S+1)-m^\prime m}\\ \langle m^\prime \vert S^z \vert m \rangle &= \delta_{m^\prime,m}m\\ \langle m^\prime \vert S^+ \vert m \rangle &= \delta_{m^\prime,m+1}\sqrt{S(S+1)-m^\prime m}\\ \langle m^\prime \vert S^- \vert m \rangle &= \delta_{m^\prime+1,m}\sqrt{S(S+1)-m^\prime m}\end{split}\]

\(\vert m\rangle\) refer to qudit states: \(\vert 0\rangle, \dots,\vert \textbf{Q}-1\rangle\).

We further define more general qudit operators:

\[\begin{split}\langle i \vert n \vert j \rangle &= \delta_{i,j} i\\ \langle i \vert R \vert j \rangle &= \delta_{i,j+1}\\ \langle i \vert R^\dagger \vert j \rangle &= \delta_{i+1,j}\\ \langle i \vert P_k \vert j \rangle &= \delta_{i,k}\delta_{k,j}\end{split}\]

Returns:: error_code

bool is_hermitian(void) const

Checks if the operator is hermitian (or symmetric if ScalarType is real).

\[M_{i,j} = \text{conj}\left(M_{j,i}\right) \, ,\,\forall i,j\]

Returns:: true if matrix is hermitian otherwise false

bool is_diagonal(void) const

Checks if the matrix is diagonal.

\[M_{i,j} = 0 \, ,\,\forall i\neq j\]

Returns:: true if matrix is diagonal otherwise false

std::map<std::string, danceq::internal::SparseMatrix<ScalarType>> get_available_operators(void) const

Returns the available operators.

Returns available_operators that act on a single site. The key value of the dictionary refers to its name, e.g., S+.

Returns:: available_operators

ScalarType trace_over_density_matrix(const Vector<ScalarType> &rho) const

Computes the trace over this operator times a (vectorized) density matrix.

Calling this function is only possible if this Operator is of OperatorType Hamiltonian. The vectorized density matrix needs to have the squared dimension of this operator dimension. Let \(O_{ij}\) be the matrix elements of this operator; the expectation value with respect to the density matrix \(\rho_{ij}\) is

\[\langle O\rangle_\rho = \text{Tr} \left( O\rho \right) =\sum_k \sum_l O_{kl}\rho_{lk} \,.\]

The function is compatible with MPI, and openMP. It further exploits the sparsity of the \(O\).

Parameters:: rho – vectorized density matrix
Returns:: expectation value

ScalarType trace_over_density_matrix(const Vec &rho) const

Computes the trace over this operator times a (vectorized) density matrix in Petsc Vec format.

This function is only available if Petsc has been included. It is equaivelnt to trace_over_density_matrix(rho).

Parameters:: rho – vectorized density matrix as Vec
Returns:: expectation value

int32_t set_Ntensor_max(const uint64_t Ntensor_max_input)

Sets Ntensor_max.

Parameters:: Ntensor_max_input – New value
Returns:: error_code

int32_t define_operator(std::string name, std::array<std::array<ScalarType, state_class::template_struct::template_Q>, state_class::template_struct::template_Q> matrix)

Adds an operator to the available operators.

Adds an operator to available_operators. If the name already exists, it will be replaced.

\[\langle i \vert O \vert j \rangle = \delta_{i,j+1}\]

std::array< std::array<double, Q>, Q> matrix;
for(uint64_t i_for = 0UL; i_for < Q; i_for++){
    for(uint64_t j_for = 0UL; j_for < Q; j_for++){
        if(i_for != j_for + 1){
            matrix[i_for][j_for] = 0.; // two-dimensional array has to be set to zer
        } else {
            matrix[i_for][j_for] = 1.;
        }
    }
}
H.define_operator("O",matrix);

Parameters:

name – Name of the operator
matrix – Interaction matrix of size Q times Q in form of an two-dimensional array

Returns:

error_code

int32_t add_operator(const ScalarType coef, const std::vector<uint64_t> sites_to_act, const std::vector<std::string> operator_string, const OperatorType OT_term = Hamiltonian, const bool print_warning = true)

Adds a new term to the operator.

coef is an overall coefficient, sites_to_act contains the sites on which the operator_string acts. OT_term is an OperatorType and can be either danceq::Hamiltonian or danceq::Lindbladian (to define a jump operator).

The operators listed in operator_string have to me included in available_operators.

Sites have to be within 0 and L-1.

sites_to_act and operator_string have to be of the same size.

If multiple operators act on the same site, the function internally defines their product and adds it to available_operators in the given order.

The following example adds a Hamiltonian term, see OperatorType, \(S^z_3S^z_4\) acting on sites 3 and 4:
H.add_operator(1., {3,4}, {"Sz","Sz"}, danceq::Hamiltonian);
Jump operators can be added with:
L.add_operator(1., {3}, {"Sz"}, danceq::Lindbladian);

All added operators are stored in operators_list. Internally, the function calls merge_terms(…) which sets up the sparse_matrices_for_computation to ensure a fast computation. Ntensor_max defines the maximal cluster size to generate a tensor matrix of size \(\textbf{Q}^{\textbf{Ntensor_max}}\). See here for more information.

Parameters:

coef – Overall coefficient
sites_to_act – Sites on which each operators act
operator_string – Operators that are applied to each site
OT_term – OperatorType (Hamiltonian/Lindbladian) (Default: Hamiltonian)
print_warning – Prints warning if true (Default: true)

Returns:

error_code

int32_t add_jump_operator(const ScalarType coef, const std::vector<uint64_t> sites_to_act, const std::vector<std::string> operators_to_apply, const bool print_warning = true)

Adds a new jump operator.

This function is a simple wrapper of add_operator(…) that fixes OperatorType OT_term=Lindbladian.

Parameters:

coef – Overall coefficient
sites_to_act – Sites on which each operators act
operators_to_apply – Operators that are applied to each site
print_warning – Prints warning if true (Default: true)

Returns:

error_code

int32_t add_diagonal(const ScalarType val)

Adds a diagonal shift.

Adds a constant to the diagonal of the operator:

\[O = O + \textbf{val} \cdot \mathbb{1}\otimes \dots \otimes \mathbb{1}\]

The diagonal is initially set to zero.

Parameters:: val – Adding val a diagonal contribution
Returns:: error_code

int32_t merge_terms(const bool transpose = false)

Identifies non-local subclusters to merge all operator strings that are fully supported on such a clusters.

The function groups sites into clusters to improve the performance and creates the std::tuple in sparse_matrices_for_computation. Each term that is listed in operator_list will be assigned to a cluster in order to optimize the computation.

An example: In an one-dimensional spin chain with periodic boundary condition and nearest-neighbor terms of size \(L=30\) we can define three clusters containing sites \(C_0=\{0,\dots,10\}\), \(C_1=\{10,\dots,20\}\), and \(C_2=\{0,20,\dots,29\}\). Hence, we only define three sparse matrices of size \(2^{11}\). For example, the terms \(S^x_3S^x_4\) and \(S^y_7S^y_8\) are assigned to the first cluster \(C_0\), and \(S^x_{16}S^x_{17}\) is assigned to the second cluster \(C_1\).

The maximal cluster size can is given by Ntensor_max and cane be set by calling set_Ntensor_max(…). Its default value is chosen such that the tensor matrix does not exceeds the dimension of 2028. However, if terms acting on more than Ntensor_max sites (e.g., long loop updates), Ntensor_max is increased. This function is called by add_operator(…).

The transposed problem will be encoded by setting transpose to true.

Parameters:: transpose – Sets the transposed problem
Returns:: error_code

uint64_t get_index(const state_class &state) const

Returns the index of a state.

If T refers to a BasisU1 the index is obtained using basisptr and get_index(…). If T refers to State Class the index is obtained using get_index_NoU1(0UL ,L).

Parameters:: state – State of interests
Returns:: Index

state_class get_state(const uint64_t index) const

Returns the state for a given index.

If T refers to a BasisU1 the index is obtained using basisptr and set_state(…). If T refers to State Class the index is obtained using set_state_from_index_NoU1(index, 0UL, L).

Parameters:: index – Index of interests
Returns:: state

int32_t increment_state(state_class &state, bool *overflow = nullptr) const

Increments the state.

If T refers to a BasisU1 the index is obtained using basisptr and increment(state). If T refers to State Class the index is obtained using increment_NoU1(0UL,L)(index, 0UL, L). If the pointer overflow is set this results in setting it to true otherwise it is false.

Parameters:

state – State to be incremented
overflow – Pointer for overflow detection

Returns:

error_code

std::vector<std::pair<ScalarType, uint64_t>> apply(const state_class &state, const bool unique_output = false) const

Applies the operator to a single state.

Creates all elements with coefficients that are obtained when the operator is applied to the input State state. The first entry in the std::pair refers to the coefficients and the second entry refers to the index within the used basis. The state can be retrieved using get_state(index), see the example below. The returned elements are ordered and double entries as well as zeros are removed.

auto elements = H.apply(psi);
std::cout << "H|psi> = ";
for (auto & element : Hpsi){
    State state__ = H.get_state(element.second); // retrieving state from basis index
    std::cout << "+ " << Helement.first << " " << state__.get_string(0UL, L) << " (Basis index: " << element.second << ")" << std::endl;
}
std::cout << std::endl;

If unique_output is true, double elements of the return std::vector are merged.

Parameters:

state – State of interests
unique_output – Boolean to control the unique output

Returns:

All created coefficients and the corresponding index

EigenMatrixType create_EigenDenseMatrix(void) const

Returns a dense matix with Eigen.

Creates and fills a dense matrix with Eigen:

auto matrix = H.create_EigenDenseMatrix();

Returns:: matrix

int32_t create_PetscSparseMatrix_reference(Mat &H, const bool show_info = false)

Returns a Petsc Mat as a MATMPIAIJ type.

Creates and fills a sparse matrix over several MPI processes using in a MATMPIAIJ format from Petsc. This functions can be called from create_PetscSparseMatrix(…).

Parameters:

H – Reference to the Petsc Mat for output
show_info – Prints information like memory usage if true

Returns:

error_code

Mat create_PetscSparseMatrix(const bool show_info = false)

Returns a Petsc Mat as a MATMPIAIJ type.

Creates and fills a sparse matrix over several MPI processes using in a MATMPIAIJ format from Petsc. This functions is a wrapper for create_PetscSparseMatrix_reference(…).

auto matrix = H.create_PetscSparseMatrix();

Parameters:: show_info – Prints information like memory usage if true
Returns:: matrix

int32_t create_PetscShellMatrix_reference(Mat &H, const bool show_info = false, const double overhead_in_GB_per_core = 2.) const

Returns a Petsc Mat as a MATSHELL type for matrix-free multiplication.

Creates a shell matrix over several MPI processes using in a MATSHELL format from Petsc. It is possible to control the maximal memory overhead per process during the multiplication using overhead_in_GB_per_core. More detailed information can be found here. This functions is called by create_PetscShellMatrix(…).

Parameters:

H – Reference to the Petsc Mat for output
show_info – Prints information like communication pattern used by MPI if true (Default: false)
overhead_in_GB_per_core – Maximal memory per core in GB used during the multiplication (Default: 2.)

Returns:

error_code

Mat create_PetscShellMatrix(const bool show_info = false, const double overhead_in_GB_per_core = 2.) const

Returns a Petsc Mat as a MATSHELL type for matrix-free multiplication.

auto matrix = H.create_PetscShellMatrix();

Parameters:

show_info – Prints information like communication pattern used by MPI if true (Default: false)
overhead_in_GB_per_core – Maximal memory per core in GB used during the multiplication (Default: 2.)

Returns:

matrix

ShellMatrix<Operator<T, ScalarType, OT>> create_ShellMatrix(const bool show_info = false, const double overhead_in_GB_per_core = 2.) const

Returns a ShellMatrix for matrix-free multiplication.

Creates a ShellMatrix that allows a matrix-free matrix-vector multiplication. The ShellMatrix works together with the Vector class. An instance of the Vector class with the correct memory layout can be retrieved using the function create_Vector(). Importantly, the ShellMatrix (and Vector class) is compatible with MPI and openMP. The syntax is equivalent in all cases. Additional memory for the communication between processes is requiered whe MPI is enabled. This can be limited by overhead_in_GB_per_core.

auto matrix = H.create_ShellMatrix();

Parameters:

show_info – Prints information like communication pattern used by MPI if true (Default: true)
overhead_in_GB_per_core – Maximal memory per core in GB used during the multiplication (Default: 2.)

Returns:

matrix

SparseMatrix<ScalarType> create_SparseMatrix(void)

Returns the matrix in a sparse format.

The operator is returned as a SparseMatrix. The function is compatible with openMP. We do not provide a MPI version as this should be done using Petsc via create_PetscSparseMatrix(…).

auto matrix = H.create_SparseMatrix();

Returns:: matrix

std::vector<std::vector<ScalarType>> create_DenseMatrix(void) const

Returns the operator in a dense format.

The operator is returned as a dense matrix in the form of a two-dimensional std::vector of size dim times dim. The function is compatible with openMP.

auto matrix = H.create_DenseMatrix();

Returns:: error_code

Private Functions

int32_t merge_given_terms(std::vector<uint64_t> &cluster, std::vector<std::pair<uint64_t, danceq::internal::ResolvedAction>> &list_of_terms, const bool transpose)

Merges given terms that are all defined on the given cluster.

This function is only called internally by merge_terms(…). It merges the terms given in list_of_terms which refers to operator_list into a single entry in sparse_matrices_for_computation that is used in apply(…). All terms are fully supported on the given cluster. The second argument in list_of_terms refers to the precise action on the wavefunction or density matrix as specified in ResolvedAction.

The transposed problem will be encoded by setting transpose to true.

Parameters:

cluster – Sites on which the matrix is defined
list_of_terms – Terms that are included
transpose – Boolean to generate the transposed problem

Returns:

error_code

Private Members

const T *basisptr

Pointer to a template T instance.

If template T refers to BasisU1, it points to an instance that is used to create the matrix. If template T refers to the State class, the pointer is not necessary and set to nullptr. It is deleted by the deconstructor if delete_basisptr is true.

uint64_t Ntensor_max

Defines the maximal number of sites of a non-local subcluster which contains all operators that are fully supported on it.

The idea is to group the system into non-local subclusters and to incorporate all operators strings that are fully supported on such a cluster. The maximal size of these clusters (if there is no operator string acting on more sites) is given by Ntensor_max. More information can be found here. The default value is chosen such that the dimension of the underlying sparse matrix does not exceed 2048. It can be set by calling set_Ntensor_max(…)

const bool delete_basisptr

Boolean to delete the basisptr.

The deconstructor deletes the basisptr if it is true.

const uint64_t L: System size.

const uint64_t n: Particle number sector. This is only set if T refers to the BasisU1.

uint64_t dim: Total dimension.

state_class mask_for_Lindbladian: A mask that is used for the bit arithmetic if terms act on density matrices.

uint64_t dim_for_Lindbladian: The dimension of the actual Hilbert space which is the square root of the total dimension in the case of a Lindbladian.

uint64_t Nbody_max: Maximal number of sites that participate in a tensor product.

std::vector<state_class> state_list

A sorted list of all states.

The vector stores a list of \(\textbf{Q}^{\textbf{Nbody_max}}\) sorted states that are used to apply a tensor product acting on Nbody_max sites (Q = 2, Nbody_max = 3):

\(\vert\Psi_0 \rangle = \vert 0;0;0\rangle\)	\(\vert\Psi_1 \rangle = \vert 0;0;1\rangle\)	\(\vert\Psi_2 \rangle = \vert 0;1;0\rangle\)	\(\vert\Psi_3 \rangle = \vert 0;1;1\rangle\)
\(\vert\Psi_4 \rangle = \vert 1;0;0\rangle\)	\(\vert\Psi_5 \rangle = \vert 1;0;1\rangle\)	\(\vert\Psi_6 \rangle = \vert 1;1;0\rangle\)	\(\vert\Psi_7 \rangle = \vert 1;1;1\rangle\)

uint64_t max_element: Maximal number of off-diagonal elements of all basis state.

ScalarType diagonal_const: Diagonal shift.

std::map<std::string, danceq::internal::SparseMatrix<ScalarType>> available_operators

Available local operators.

Dictionary maps names to local operators acting on single site, e.g.:

\[\langle m^\prime \vert S^x \vert m \rangle = \left( \delta_{m^\prime,m+1} + \delta_{m^\prime,m+1}\right)\frac{1}{2}\sqrt{S(S+1)-m^\prime m}\]

Default operators can be found here. New local operators can be added with define_operator(…).

std::vector<std::tuple<ScalarType, std::vector<uint64_t>, std::vector<std::string>, OperatorType>> operator_list

Operator list for bookkeeping.

This is only used for bookkeeping and not for computation. The first element of the tuple contains the overall coefficient. The second element stores the sites in form of std::vector on which each term acts. The third element, std::vector of std::string, is of the same size as the second element and contains the operators that are listed in available_operators. The last element refers to the OperatorType.

For example, the operator string \(3\cdot S^x_5S^x_6\) as a term of the Hamiltonian is stored as:

std::tuple<ScalarType,std::vector<uint64_t>, std::vector<std::string>, danceq::OperatorType> = std::make_tuple(3.,{5,6},{"Sx","Sx"},Hamiltonian)

An example for a jump operator is:

std::tuple<ScalarType,std::vector<uint64_t>, std::vector<std::string>, danceq::OperatorType> = std::make_tuple(.1,{0},{"Sx"},Lindbladian)

std::vector<std::tuple<std::vector<uint64_t>, std::vector<uint64_t>, std::vector<std::pair<state_class, ScalarType>>, state_class>> sparse_matrices_for_computation

Tensor matrices for computation.

The tuples are used to compute the action of the Operator to a State in the apply(…) function. The first element of each std::tuple contains the cluster on which each non-local tensor product acts. Similar to the SparseMatrix, the second element contains the range and the third element contains the data with the correct mask used to immetetly generate the new State. The last element is the mask that cuts out states defined by the cluster.

It is set up by merge_terms(…) which is called within add_operator(…).

bool bool_is_hermitian: Boolean to check if the input is hermitian.

bool bool_is_diagonal: Boolean to check if the input is diagonal.