Tensor Network Theory Library  Beta release 1.2.1
A library of routines for performing TNT-based operations
 All Data Structures Functions Variables Groups Pages
tntGS_cl.c

Example of using the MPS library for finding the ground state of a 1D spin or boson Hamiltonian using parameters specified on the command line. The program does not require an initialisation file.

The ground state is found using two-site updates, and the internal dimension can grow up to a maximum \(\chi\) during the simulation. The precision (i.e. energy difference between sweeps to indicate ground state has been found) is set to a default of 1e-4, but be changed from the command line.

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 the 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.
--max-it <n> Maximum number of iterations. If not given default 50 is set.
--precision <val> When energy difference between successive DMRG iterations is smaller than this value, the calculation ends. If not given default 1e-4 is used.
--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.
--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.
--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 an MPO representing the system Hamiltonian. 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.

It stores the operators required in the tntMpsExOp structure. If the flag for a given expecation value is given, then the corresponding operators are placed in the structure. This structure is then passed to tntMpsExpecOutput() to calculate expectation values when they are required.

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>=**.

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\).

Finding the ground state

The program uses the building block functions provided in the MPS library to perform the steps of the DMRG algorithm. Function tntMpsVarMinMpoInit() sets up the 'effective Hamiltonians' for each pair of sites using the starting random wave function and the MPO representing the Hamiltonian.

Minimisation sweeps are then performed until the difference in the energy between sweeps is less than the precision. Once the energy has converged, expectation values are calculated, and the data is saved to a MATLAB output file.

Running the program

To compile this main file change to the source directory then type make tntGS_cl . You may first have to amend makefile.inc to point to your copies of the 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/tntGS_cl with the parameters of your choice. For example:

./bin/tntGS_cl -d output/bh_harm --system boson --length 41 --qnum-rand-state 45 --n-max 3 --Ub 4 --Jb -1 --E-harm 0.25 -c 20 --Ex1N --Ex2bdagb=ap -T -1 to find the ground state of the Bose-Hubbard model with exactly 10 bosons in a harmonic trapping potential. Try increasing the interaction and seeing the difference in the density distribution and the single particle density matrix.

bhmodeldensity.png
/*
Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
Date: $LastChangedDate$
(c) University of Oxford 2015
*/
/* Include the header for the TNT library, including MPS functions */
#include "tntMps.h"
int main(int argc, char **argv)
{
tntNetwork wf, H; /* networks representing the MPS wave function and the MPO Hamiltonian */
tntNodeArray HeffL, HeffR; /* Precontracted nodes for the right and left sides of the network */
char saveprefix[TNT_STRLEN]; /* Path to output file */
int i, i_max; /* Iteration counter, Maximum number of iterations */
double prec; /* Precision for the calculation */
double E, Eprev; /* Energy for this iteration */
double err = 0.0; /* The truncation error */
int chi = -1; /* Maximum internal dimension for the MPS */
tntMpsExOp ExOp; /* Defines all the operators for calculating expectation values */
/* Initialize the TNT library */
/* Set a random number seed based on the current time */
tntSysRandSeedSet(time(NULL));
/* Create the system paramaters from command line arguments */
tntProcessCLOptions(argc, argv, NULL, NULL, saveprefix, &chi, NULL);
tntMpsProcessDmrgCLOptions(argc, argv, &i_max, &prec);
tntMpsProcessSysOptions(argc, argv, &wf, &H, NULL, &ExOp);
/* Initialise the nodes that will be needed for the DMRG sweep */
tntMpsMpoMpsInit(wf, H, &HeffL, &HeffR);
/* Determine the energy of the start state (the norm of the state is assumed to be 1) */
tntDoubleParamsSave(saveprefix,E);
printf("--------------------------------------------------------------------------------\n");
printf(" Starting minimisation sweeps to find ground state from command line parameters \n");
printf("--------------------------------------------------------------------------------\n");
printf("Starting energy is %g. \n", E);
/* Perform minimization sweeps */
for (i = 0; i < i_max; i++) {
Eprev = E;
E = tntMpsVarMinMpo2sStep(wf, chi, H, &HeffL, &HeffR, &err);
tntDoubleParamsUpdate(saveprefix,i+1,E);
printf("Energy is %g. Difference is %g\n", E, Eprev-E);
if ((Eprev - E < prec) && (Eprev - E > 0.0)) break;
}
tntMpsExpecOutput(wf, &ExOp, 0, 1, 1, saveprefix, 0);
/* Free all the dynamically allocated nodes and associated dynamically allocated arrays. */
return 0;
}