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

Example of using the MPS library for time evolution using parameters specified on the command line, and calculation of expectation values and wave function overlap. The program does not require an initialisation file.

The start state is evolved using the TEBD algorithm, with the site-wide propagator decomposed using a second-order Suzuki-Trotter decomposition into a sequence of two-site propagators. The internal dimension can grow up to a maximum $$\chi$$ during the simulation.

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.
--rand-state <n> Create a random start state with internal dimension D given by the integer argument.
--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.
--config-state <str>Create a starting state having a configuration given by the supplied string.
--qnum-config-state <n> Create a starting state having a configuration given by the supplied string, and conserve the total quantum number of this configuration during evolution.
--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 tntMpsPropST2scProduct().

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 tntEvolve_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/tntEvolve_cl with the parameters of your choice. For example:

./bin/tntEvolve_cl -d output/xxz --system spin --length 20 --qnum-config-state 1111111111 --Jxx 1 --Jyy 1 --Jzz 0.5 -c 20 --Ex1Sz -t 1000 -b 50 --dt 0.01 to evolve a system with half spins up and half spins down under the XXZ Hamiltonian. Try changing the ZZ coupling and seeing how this affects the overlap between the start and finish state.

./bin/tntEvolve_cl -d output/bh_alt --system boson --length 31 --qnum-config-state 0101010101010101010101010101010 --n-max 1 --Ub 5 --Jb -1 --E-harm 0.02 --Ex1N -t 400 -b 10 --dt 0.01 -c 200 to evolve a system with alternating filled and empty sites under the Bose-Hubbard Hamiltonian with exactly 8 bosons. This is a smaller scale of the types of calculations performed in Trotzky, S., Chen, Y.-A., Flesch, A., McCulloch, I. P., Schollwöck, U., Eisert, J., & Bloch, I. (2011). Probing the relaxation towards equilibrium in an isolated strongly correlated 1D Bose gas. Nature Physics, 8(4), 8.

/*
Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
Date: $LastChangedDate$
(c) University of Oxford 2015
*/
#include "tntMps.h"
int main(int argc, char **argv)
{
tntNetwork wf, wfstart, prop; /* networks representing the MPS wave function and propagator */
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 */
int chi; /* Maximum internal dimension for the MPS */
double overlap = 1.0; /* Overlap - starts at zero */
tntMpsExOp ExOp; /* Defines all the operators for calculating expectation values */
/* Initialize the TNT library */
/* Process the command line line options */
tntProcessCLOptions(argc, argv, NULL, NULL, saveprefix, &chi, NULL);
tntMpsProcessTebdCLOptions(argc, argv, &numsteps, &tbigstep, NULL, NULL, NULL, NULL);
tntMpsProcessSysOptions(argc, argv, &wf, NULL, &prop, &ExOp);
/* Initial expectation values */
tntMpsExpecOutput(wf, &ExOp, 0, 1, 1, saveprefix, 0);
/* Make a copy of the initial wavefunction and save initial overlap */
wfstart = tntNetworkCopy(wf);
tntDoubleParamsSave(saveprefix,overlap);
printf("--------------------------------------------------------------------------------\n");
printf(" Starting time evolution from command line parameters \n");
printf("--------------------------------------------------------------------------------\n");
/* Run the simulation for all the time steps. */
for (tstep = 1; tstep <= numsteps; tstep++) {
err_step = tntMpsPropST2scProduct(wf, prop, 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++;
tntMpsExpecOutput(wf, &ExOp, 0, 1, 1, saveprefix, bigstep);
overlap = tntComplexToReal(tntMpsMpsProduct(wf,wfstart));
tntIntParamsUpdate(saveprefix, bigstep, tstep);
tntDoubleParamsUpdate(saveprefix, bigstep, err, overlap);
}
}
printf("The truncation error is %g.\n", err);
/* Free all the dynamically allocated nodes and associated dynamically allocated arrays. */
tntNetworkFree(&wfstart);
/* Finish with the TNT library */
return 0;
}