Tensor Network Theory Library  Beta release 1.2.1 A library of routines for performing TNT-based operations
tntThermal_cl.c

Calculates the density operator for a thermal state $$\rho=\mathrm{e}^{-\beta H}$$ (where $$\beta$$ is the inverse temperature) by performing imaginary time evolution on an identity operator. The state is represented as an MPO.

This uses the relation $$\mathrm{e}^{-\beta H} = \left( \mathrm{e}^{-\Delta t H} \right)^M \mathbf{1} \left( \mathrm{e}^{-\Delta t H} \right)^M$$, with $$\Delta t = \frac{\beta}{2M}$$. See Matrix Product Density Operators: Simulation of Finite-Temperature and Dissipative Systems, F . Verstraete, J. J. GarcĂ­a-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 for further details.

The Hamiltonian for the thermal state is defined using command line parameters. Use --idt to set the imaginary time-step $$\Delta t$$ and -t or --timesteps to set $$M$$.

### Setting system parameters

The program uses the appropriate library functions for handling command line parameters to allow the system to be initialised.

Command line parameters allowed for this program are:

Long nameShort nameDescription
--chi <n>-c <n> The maximum internal dimension $$\chi$$ of nodes in the network. If not given it is set to -1 meaning no truncation will be performed.
--TOL <val>-T <val> Tolerance for zeroing matrix values in SVD. See tntSVDTolSet() for more information. If not given default of 1e-14 will be used.
--directory <dir> -d <dir> Output directory and prefix for the output files.
--timesteps <n> -t <n> Number of timesteps to be performed.
--tbigstep <n> -b <n>Number of timesteps to be performed between each calculation of expectation values. If not given, then expectation values are calculated at the end of the routine only.
--trunc-type <str> How to calculate the truncation error using the discarded singular values. See tntTruncType() for allowed options.
--dt <val> Size of real time step (real double). Required to build propagator if –idt is not given. Default is 0.
--idt <val> Size of imaginary time step (real double). Required to build propagator if –dt is not given. Default is 0.
--system [spin-half | boson] System type. Required if operators are to be built.
--qnum-rand-state <n> Create a random start state with total quantum number Q given by the integer argument, and conserve the total quantum number during evolution.
--eye-state <str>Start with an identity matrix.
--qnum-eye-state <n> Start with an identity matrix and keep track of global physical symmetries.
--length <n> Number of sites in the MPS. Required.

If --system=spin is chosen, then the following options are available:

Long nameDescription
--2S <n> Twice the spin $$S$$ of each particle in the system. Default is 1 i.e. a spin half system.
--Jxx <cval> Adds the term $$J_{xx}\hat{S}^{x}_{j}\hat{S}^{x}_{j+1}$$ to the Hamiltonian for a complex $$J_{xx} = \mathtt{<cval>}$$.
--Jyy <cval> Adds the term $$J_{yy}\hat{S}^{y}_{j}\hat{S}^{y}_{j+1}$$ to the Hamiltonian for a complex $$J_{yy} = \mathtt{<cval>}$$.
--Jzz <cval> Adds the term $$J_{zz}\hat{S}^{z}_{j}\hat{S}^{z}_{j+1}$$ to the Hamiltonian for a complex $$J_{zz} = \mathtt{<cval>}$$.
--Jxy <cval> Adds the term $$J_{xy}\hat{S}^{x}_{j}\hat{S}^{y}_{j+1}$$ to the Hamiltonian for a complex $$J_{xy} = \mathtt{<cval>}$$.
--Jyx <cval> Adds the term $$J_{yx}\hat{S}^{y}_{j}\hat{S}^{x}_{j+1}$$ to the Hamiltonian for a complex $$J_{yx} = \mathtt{<cval>}$$.
--Jyz <cval> Adds the term $$J_{yz}\hat{S}^{y}_{j}\hat{S}^{z}_{j+1}$$ to the Hamiltonian for a complex $$J_{yz} = \mathtt{<cval>}$$.
--Jzy <cval> Adds the term $$J_{zy}\hat{S}^{z}_{j}\hat{S}^{y}_{j+1}$$ to the Hamiltonian for a complex $$J_{zy} = \mathtt{<cval>}$$.
--Jzx <cval> Adds the term $$J_{zx}\hat{S}^{z}_{j}\hat{S}^{x}_{j+1}$$ to the Hamiltonian for a complex $$J_{zx} = \mathtt{<cval>}$$.
--Jxz <cval> Adds the term $$J_{xz}\hat{S}^{x}_{j}\hat{S}^{z}_{j+1}$$ to the Hamiltonian for a complex $$J_{xz} = \mathtt{<cval>}$$.
--Bx <cval> Adds the term $$B_x\hat{S}^{x}_{j}$$ to the Hamiltonian for a complex $$B_x = \mathtt{<cval>}$$.
--By <cval> Adds the term $$B_y\hat{S}^{y}_{j}$$ to the Hamiltonian for a complex $$B_y = \mathtt{<cval>}$$.
--Bz <cval> Adds the term $$B_z\hat{S}^{z}_{j}$$ to the Hamiltonian for a complex $$B_z = \mathtt{<cval>}$$.

If --system=boson is chosen, then the following options are available:

Long nameDescription
--n-max <n> Maximum number of bosons allowed per site. Default is 1.
--Jb <cval> Adds a hopping term to the Hamiltonian with the form $$-(J\hat{b}_{j}^{\dagger}\hat{b}_{j+1}+J^{\ast}\hat{b}_{j}\hat{b}_{j+1}^{\dagger})$$ for a complex $$J = \mathtt{<cval>}$$.
--Ub <cval> Adds an on-site interaction term to the Hamiltonian with the form $$\frac{U}{2}\hat{n}_{j}(\hat{n}_{j}-1)$$ for a complex $$U = \mathtt{<cval>}$$.
--Vb <cval> Adds a nearest-neighbour interaction term to the Hamiltonian with the form $$V\hat{n}_{j}\hat{n}_{j+1}$$ for a complex $$V = \mathtt{<cval>}$$.
--mub <cval> Adds a chemical potential term to the Hamiltonian with the form $$-\mu\hat{n}_{j}$$ for a complex $$\mu = \mathtt{<cval>}$$.
--E-harm <cval> Adds a harmonic potential term to the Hamiltonian with the form $$\epsilon_{\mathrm{harm}}\hat{n}_{j}(j-j_c)^2/2$$ for a complex $$E_{h} = \mathtt{<cval>}$$.
--jc-harm <val> Centre of harmonic potential $$j_c$$. Only applied when --E-harm is set and must be a real value. Default is centre of the system.
--coher-driv <cval> Adds a coherent driving term to the Hamiltonian with the form $$\Omega\hat{b}_{j}+\Omega^{\ast}\hat{b}_{j}^{\dagger}$$ for a complex $$\Omega = \mathtt{<cval>}$$.
--param-driv <cval> Adds a parametric driving term to the Hamiltonian with the form $$\Delta\hat{b}_{j}\hat{b}_{j+1}+\Delta^{\ast}\hat{b}_{j}^{\dagger}\hat{b}_{j+1}^{\dagger}$$ for a complex $$\Delta = \mathtt{<cval>}$$.

As described in the table some parameters can take complex values - these must be printed without spaces e.g. the following types of numbers would be recognised 2.1,-1e4,10i,3-2i,-4i+1e-5 but 3 + 3i would not.

These terms are processed to create a ::tntNetwork for a staircase propagator network representing a time-step of the system-wide evolution operator using the function tntMpsCreatePropST2sc(). It also creates a ::tntNetwork representing the starting MPS as well as setting the relevant simulation parameters to their default values or the values passed on the command line.

This program also processes input command line arguments to create the variables needed for calculating expectation values for spin or boson simulations. Using the options given on the command line, the function can create arrays containing operators for:

1. On-site expectation values
2. Nearest neighbour expectation values
3. All pairs of sites for two-site operators
4. Centre to all other sites for two site operators.

Setting the following flags will indicate that the value $$\langle|\psi|\hat{o}_j|\psi\rangle$$ should be calculated for every site $$j$$ i.e. $$L$$ values will be calculated.

For bosonic systems the allowed flags are:

Long nameDescription
--Ex1N Calculate $$\langle\psi|\hat{n}|\psi\rangle$$
--Ex1Int Calculate $$\langle\psi|\hat{n}(\hat{n}-1)/2|\psi\rangle$$
--Ex1Nsq Calculate $$\langle\psi|\hat{n}^2|\psi\rangle$$
--Ex1Q1 Calculate $$\langle\psi|\hat{b} + \hat{b}^{\dagger}|\psi\rangle$$
--Ex1Q2 Calculate $$\langle\psi|\mathrm(i)(\hat{b} - \hat{b}^{\dagger})|\psi\rangle$$

For spin systems the allowed flags are:

Long nameDescription
--Ex1Sx Calculate $$\langle\psi|\hat{S}^{x}|\psi\rangle$$
--Ex1Sy Calculate $$\langle\psi|\hat{S}^{y}|\psi\rangle$$
--Ex1Sz Calculate $$\langle\psi|\hat{S}^{z}|\psi\rangle$$

Setting the following flags will indicate that the value $$\langle|\psi|\hat{o}^{[1]}_{j_1}\hat{o}^{[2]}_{j_2}|\psi\rangle$$ should be calculated for pairs of sites $$j_1$$ and $$j_2$$. There are three different choices for the pairs of sites that will be chosen, each with a two-character code set as the argument of the flag (where an '=' must be used) i.e. Ex2<opname>=**.

• For nearest-neighbour expectation values it will be calculated for $$i=0$$ to $$L-2$$ and $$j=i+1$$ i.e. $$L-1$$ values will be calculated. Use the code nn to specify this e.g. --Ex2SzSz=nn.
• For centre-site expectation values it will be calculated for $$i=i_c$$ where $$i_c$$ is the central site $$\mathrm{floor}(L/2)$$ and $$j=0$$ to $$L-1$$ i.e. $$L$$ values will be calculated. Use the code cs to specify this e.g. --Ex2SzSz=cs.
• For all pairs of expectation values it will be calculated for both $$i=0$$ to $$L-1$$ and $$j=0$$ to $$L-1$$ i.e. $$L^2$$ values will be calculated. Use the code ap to specify this e.g. --Ex2SzSz=ap.

If no two-character code is given, or the two-character code is not recognised, then nearest-neighbour pairs will be chosen.

For bosonic systems the allowed flags are:

Long nameDescription
--Ex2bdagb Calculate $$\langle\psi|\hat{b}^{\dagger}_{j_1}\hat{b}_{j_2}\rangle$$.
--Ex2den Calculate $$\langle\psi|\hat{n}_{j_1}\hat{n}_{j_2}\rangle$$.

For spin systems the allowed flags are:

Long nameDescription
--Ex2SxSx Calculate $$\langle\psi|\hat{S}^{x}_{j_1}\hat{S}^{x}_{j_2}\rangle$$.
--Ex2SySy Calculate $$\langle\psi|\hat{S}^{y}_{j_1}\hat{S}^{y}_{j_2}\rangle$$.
--Ex2SzSz Calculate $$\langle\psi|\hat{S}^{z}_{j_1}\hat{S}^{z}_{j_2}\rangle$$.
--Ex2SxSy Calculate $$\langle\psi|\hat{S}^{x}_{j_1}\hat{S}^{y}_{j_2}\rangle$$.
--Ex2SySx Calculate $$\langle\psi|\hat{S}^{y}_{j_1}\hat{S}^{x}_{j_2}\rangle$$.
--Ex2SySz Calculate $$\langle\psi|\hat{S}^{y}_{j_1}\hat{S}^{z}_{j_2}\rangle$$.
--Ex2SzSy Calculate $$\langle\psi|\hat{S}^{z}_{j_1}\hat{S}^{y}_{j_2}\rangle$$.
--Ex2SzSx Calculate $$\langle\psi|\hat{S}^{z}_{j_1}\hat{S}^{x}_{j_2}\rangle$$.
--Ex2SxSz Calculate $$\langle\psi|\hat{S}^{x}_{j_1}\hat{S}^{z}_{j_2}\rangle$$.
--Ex2SpSm Calculate $$\langle\psi|\hat{S}^{+}_{j_1}\hat{S}^{-}_{j_2}\rangle$$.

### Evolving the state

The program uses the building block functions provided in the MPS library to perform the steps of the TEBD algorithm. The starting wave function is first copied to a separate variable. For each time-step, the staircase network representing a time-step of the propagator is multiplied with the MPS representing the current wave function using tntMpoPropST2scProduct().

The expectation values are calculated every 'big time step' and saved. The overlap with the starting wave function is also calculated.

### Running the program

To compile this main file change to the source directory then type  make tntThermal_cl . You may first have to amend makefile.inc to point to your copies of the required external linear algebra and MATLAB libraries.

To run the file you may also have to set up certain environment variables. To do this edit then run ./set_tnt_vars.sh in the scripts directory.

Now return to the root examples directory, and call ./bin/tntThermal_cl with the parameters of your choice.

/*
Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
Date: $LastChangedDate$
(c) University of Oxford 2015
*/
/* Include the header for the TNT library */
#include "tntMps.h"
int main(int argc, char **argv)
{
tntNetwork rho; /* network representing the MPO */
tntNetwork propB, propT, H; /* Propagator for the evolution, and the MPO Hamiltonian */
char saveprefix[TNT_STRLEN]; /* Path to output file */
int numsteps, tbigstep; /* Number of steps for run and between expectation values */
int tstep = 0, bigstep = 0; /* Counters for the above */
double err = 0.0, err_step = 0.0; /* The truncation error */
tntComplex E; /* The current energy */
int chi; /* maximum internal dimension for the time evolution */
tntMpsExOp ExOp; /* Structure for holding all the operators for calculating expectation values */
/* Initialize the TNT library */
/* Process the options supplied on the command line */
tntProcessCLOptions(argc, argv, NULL, NULL, saveprefix, &chi, NULL);
tntMpsProcessTebdCLOptions(argc, argv, &numsteps, &tbigstep, NULL, NULL, NULL, NULL);
tntMpsProcessSysOptions(argc, argv, &rho, &H, &propB, &ExOp);
propT = tntNetworkCopy(propB);
/* Print initial expectation values */
tntMpoExpecOutput(rho, &ExOp, 1, 1, saveprefix, 0);
E = tntMpoMpoTrace(rho,H);
tntComplexParamsUpdate(saveprefix,bigstep,E);
printf("--------------------------------------------------------------------------------\n");
printf(" Starting imaginary time evolution on identity operator to find thermal state\n");
printf("--------------------------------------------------------------------------------\n");
/* Run the simulation for all the time steps. */
for (tstep = 1; tstep <= numsteps; tstep++) {
/* Time step then renormalise */
err_step = tntMpoPropST2scProduct(rho, propB, propT, chi);
err += err_step;
/* Calculate expectation values every tbigstep */
if (0 == tstep%tbigstep || tstep == numsteps) {
printf("Completed %d out of %d steps. \n", tstep, numsteps);
bigstep++;
tntMpoExpecOutput(rho, &ExOp, 1, 1, saveprefix, bigstep);
E = tntMpoMpoTrace(rho,H);
tntComplexParamsUpdate(saveprefix,bigstep,E);
tntIntParamsUpdate(saveprefix, bigstep, tstep);
tntDoubleParamsUpdate(saveprefix, bigstep, err);
}
}
/* Free all the dynamically allocated nodes and associated dynamically allocated arrays. */
tntNetworkFree(&propT);
tntNetworkFree(&propB);
/* Finish with the TNT library */
return 0;
}