Tensor Network Theory Library  Beta release 1.0
A library of routines for performing TNT-based operations
MPS Library

Modules

 Time evolution
 
 Variational Minimisation
 
 Deprecated
 
 Command line
 

Data Structures

struct  tntMpsExOp
 

Macros

#define TNT_MPS_L
 
#define TNT_MPS_R
 
#define TNT_MPS_D
 
#define TNT_MPS_U
 
#define TNT_MPS_SPIN
 
#define TNT_MPS_BOSON
 
#define TNT_MPS_TUL
 
#define TNT_MPS_TUR
 
#define TNT_MPS_TDL
 
#define TNT_MPS_TDR
 

Functions

void tntMpsCreateBosonOp (unsigned n_max, tntNode *b, tntNode *bdag, tntNode *n, tntNode *os_int, tntNode *eye)
 
tntNetwork tntMpsCreateConfig (unsigned L, const char *config)
 
tntNetwork tntMpsCreateMpo (unsigned L, tntNodeArray *nnl, tntNodeArray *nnr, tntComplexArray *nnparam, tntNodeArray *os, tntComplexArray *osparam)
 
tntNetwork tntMpsCreateRandom (unsigned L, unsigned D)
 
void tntMpsCreateSpinOp (unsigned TwoS, tntNode *Sx, tntNode *Sy, tntNode *Sz, tntNode *Sp, tntNode *Sm, tntNode *eye)
 
tntNetwork tntMpsCreateSymmRandom (unsigned L, int *qn)
 
tntNode tntMpsCreateTwoSiteOp (tntNodeArray *Lnodes, tntNodeArray *Rnodes, tntComplexArray *params)
 
void tntMpsExpecOutput (tntNetwork mps, tntMpsExOp *Op, int orthCentre, unsigned printOutput, unsigned saveOutput, char *saveprefix, unsigned counter)
 
void tntMpsMpoConnect (tntNetwork mps, tntNetwork mpo)
 
double tntMpsMpoContract (tntNetwork mpsmpo, int chi)
 
double tntMpsMpoProduct (tntNetwork mps, tntNetwork mpo, int chi)
 
tntNetwork tntMpsMpoMpsConnect (tntNetwork mps, tntNetwork mpo)
 
tntComplex tntMpsMpoMpsContract (tntNetwork sandwich)
 
tntComplex tntMpsMpoMpsProduct (tntNetwork mps, tntNetwork mpo)
 
void tntMpsMpsConnect (tntNetwork mpsA, tntNetwork mpsB)
 
tntComplex tntMpsMpsContract (tntNetwork sandwich, int site_left, int site_right)
 
tntComplex tntMpsMpsProduct (tntNetwork mpsA, tntNetwork mpsB)
 
double tntMpsSelfProduct (tntNetwork mps, int orth_centre)
 
void tntMpsOrth (tntNetwork mps, unsigned orth_centre)
 
void tntMpsPmpoProduct (tntNetwork mps, tntNodeArray *op, tntIntArray *sitenum)
 
tntComplex tntMpsPmpoMpsProduct (tntNetwork mps, tntNodeArray *op, tntIntArray *sitenum, int orth_centre)
 
double tntMpsTruncate (tntNetwork mps, unsigned orth_centre, int chi)
 
unsigned tntMpsLength (tntNetwork mps)
 
void tntMpsExOpFree (tntMpsExOp *ExOp)
 

Detailed Description

This library can be used in addition to the core library (libtnt.a), and contains network functions for matrix product states (MPS) in systems with open boundary conditions. The functions in this library are available by linking to the static library libtntMps.a.

These functions assume an MPS network of the following form:

mps.png

where the diamonds represent the network start and end, and the dashed lines represent singleton legs. The left internal leg is labelled with constant TNT_MPS_L (having value 1), the right internal leg is labelled with TNT_MPS_R (having value 2), and the physical leg is labelled with TNT_MPS_D (having value 3).

single_site_mps.png

The functions also currently accept three types of operators:

Matrix product operators (MPOs)

An MPO network is a site-wide operator, formed from nodes having two physical legs and two internal legs, as shown below.

mpo_op_sing.png

When joined along the internal legs to form a network, they can represent a site-wide operator that does not necessarily have to be able to be expressed as a product of single-site operators.

mpo_op_nw.png

Although these types of MPOs are not restricted to this, the library functions are currently able to build site-wide operators formed from a sum of nearest neighbour and on-site terms. The internal dimension (i.e. the size of the left and right legs) will be equal to the number of nearest-neighbour terms + 2.

To build an MPO, either use the function tntMpsProcessInitOptions() to create the system (including an MPO) from parameters specified on the command line, or prepare an MPO using single site operators (see below) within your function using tntMpsCreateMpo().

Single site operators or produt MPOs

These should have two physical legs, labelled as shown below.

single_site_op.png

They can be thought of as a subset of the MPOs described above, with an internal dimension of 1, so that the legs do not have to be specified.

To create the most common operators within your code use tntMpsCreateBosonOp() and tntMpsCreateSpinOp().

To load these types of operators from MATLAB, express the operator as a matrix, and then label the first index (rows) as 3, and the second index (columns) as 4. e.g. in your MATLAB script for creating the initialisation file:

sx = [0, 1;1, 0];
sx_id = [3, 4];
* 

and then use tntLoadNodes() to load these parameters from your .mat file as usual.

Two-site gates or propagators

These should have four physical legs, labelled as shown below.

two_site_gate.png

To load these types of operators from MATLAB, express each site operator as a matrix, use the kron() function to create two-site operators, and then reshape to a four dimensional tensor. If you use kron(id,op) to create the right-hand operator, and kron(op,id) to create the left-hand operator (where id is the identity matrix, and op is a matrix representing the single site operator), then the indices should be labelled 4, 3, 2, and 1. For example to create a two-site operator \(\hat{U}_{j,j+1}\) for evolution under the two site term \(\hat{h}_{j,j+1}= -\hat{\sigma}^z_{j}\hat{\sigma}^z_{j+1} + B_{\mathrm{left}}\hat{\sigma}^x_{j} + B_{\mathrm{right}}\hat{\sigma}^x_{j+1}\) use the following MATLAB code.

sx = [0, 1;1, 0]; sz = [1, 0;0, -1]; id = eye(2);
H = -kron(sz,id)*kron(id,sz) + B_left*kron(sx,id) + B_right*kron(id,sx);
U = expm(-1i*dt*H/2);
U = reshape(U,[2,2,2,2]);
U_id = [4, 3, 2, 1];

The above code creates an unitary operator suitable for using in TEBD to evolve the system under the transverse Ising Hamiltonian for half a time-step.

Both of these types of operators can be functional (i.e. be a function of other operators and parameters, see tntNode documentation for more information) or static. The code above creates a static node. If you wanted to create a functional version of U above, for example because you were varying B_left and B_right during your evolution (either in time or space) you could instead create the operator in your MATLAB script as follows:

Hidt = cell(4,1);
Hidt{1} = 'exp';
Hidt{2} = 1i*dt*kron(sz,id)*kron(id,sz)/2;
Hidt{3} = -1i*dt*kron(sx,id)/2;
Hidt{4} = -1i*dt*kron(id,sx)/2;

Hidt{2} = reshape(Hidt{2},[2,2,2,2]);
Hidt{3} = reshape(Hidt{3},[2,2,2,2]);
Hidt{4} = reshape(Hidt{4},[2,2,2,2]);

H_id = [4,3,2,1];

To use the operator, first use tntLoadNodes() to load these parameters from the .mat file (it is loaded in eactly the same way as static nodes). However before using it the parameters for each of the terms must be set e.g.

double J_curr, B_left_curr, B_right_curr;
tntLoadNodes(loadname,1,&U,"Hidt","H_id");
tntNodeSetRealParam(U,0,J_curr);
tntNodeSetRealParam(U,1,B_left_curr);
tntNodeSetRealParam(U,2,B_right_curr);

These parameters can then be changed depending on the position or the time step.

Macro Definition Documentation

#define TNT_MPS_BOSON

Code to denote the system is a boson system.

#define TNT_MPS_D

Label (value 3) of the downwards facing physical leg of the MPS node, MPO node, or one-site operator node.

#define TNT_MPS_L

Label (value 1) of the left facing internal leg of the MPS node or MPO node.

Examples:
dmrg.c.
#define TNT_MPS_R

Label (value 2) of the right facing internal leg of the MPS node or MPO node.

Examples:
dmrg.c.
#define TNT_MPS_SPIN

Code to denote the system is a spin system.

#define TNT_MPS_TDL

Label (value 3) of the top left physical leg of a two-site operator node.

#define TNT_MPS_TDR

Label (value 4) of the top left physical leg of a two-site operator node.

#define TNT_MPS_TUL

Label (value 1) of the top left physical leg of a two-site operator node.

#define TNT_MPS_TUR

Label (value 2) of the top left physical leg of a two-site operator node.

#define TNT_MPS_U

Label (value 4) of the upwards facing physical leg of the flipped MPS node, MPO node, or one-site operator node.

Function Documentation

void tntMpsCreateBosonOp ( unsigned  n_max,
tntNode b,
tntNode bdag,
tntNode n,
tntNode os_int,
tntNode eye 
)

Creates single site operators for bosonic systems, with two legs (i.e. no internal legs), numbered to match those expected by the MPS library functions:

single_site_op.png

Each leg will have dimension \(n_\mathrm{max}+1\), where \(n_\mathrm{max}\) is the maximum number of bosons on each site.

All nodes are created as static nodes.

Only pointers to the required nodes need be passed as arguments. If any of the nodes are not required, then pass NULL.

Returns
No return value.
Parameters
n_maxMaximum number of bosons allowed on each site
bNode for annihilation operator \(\hat{b}\)
bdagNode for creation operator \(\hat{b}^{\dagger}\)
nNode for \(n = \hat{b}^{\dagger}\hat{b}\)
os_intNode for onsite interaction term \(\hat{n}(\hat{n}-1) = \hat{b}^{\dagger}\hat{b}^{\dagger}\hat{b}\hat{b}\)
eyeNode for the identity matrix with dimension \(n_\mathrm{max}+1\)
35 {
36  tntComplexArray barr, eyearr; /* Ararys to create the operators: only b and identity are created from elements, and the remaining arrays created form these */
37  unsigned physdim = n_max+1; /* Physical dimension is 2 for a spin half system */
38  unsigned num_legs = 2, num_vals = physdim*physdim; /* Define array properties. Adding on an extra leg for contracting to make two site gates or to make MPO */
39  unsigned leg_ids[2] = {3,4}; /* Legs for the operators */
40  unsigned leg_dims[2] = {physdim,physdim}; /* the physical legs */
41  unsigned leg_map_t[5] = {0, 1, 2, 4, 3}; /* Map for taking transpose of operators */
42  int leg_map[5] = {0, 1, 2, 3, 4}; /* Map for contracting operators */
43  unsigned j; /* Row number */
44  tntNode tnA, tnB; /* Used for contracting nodes */
45  tntNode bl, bdagl; /* Local copy of nodes */
46  tntIntArray qninf; /* Quantum number information for the physical legs */
47 
48  /* Quantum number information to set for invariant operators only */
49  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
50  qninf = tntIntArrayAlloc(physdim);
51  for (j = 1; j < physdim; j++) {
52  qninf.vals[j] = j;
53  }
54  }
55 
56  /* Check whether each node pointer is NULL, before creating each node
57  (apart from b and bdag, which need to be created to create all the other nodes, but will be freed and set back to NULL at the end
58  of the routine if not required). */
59 
60  /* \b = [ 0 1 0 0 0 ...]
61  [ 0 0 sqrt(2) 0 0 ...]
62  [ 0 0 0 sqrt(3) 0 ...]
63  [ 0 0 0 0 2 ...]
64  [ 0 0 0 0 0 ...]
65  [ ..... ] */
66  /* Create empty array */
67  barr = tntComplexArrayAlloc(num_vals);
68 
69  /* Put elements in row j-1, col j */
70  for (j = 1; j < physdim; j++) {
71  barr.vals[j-1 + j*physdim].re = sqrt((double) j);
72  }
73 
74  /* Create node from elements */
75  bl = tntNodeCreate(&barr, num_legs, leg_ids, leg_dims);
76 
77  /* Free the array - it is no longer needed */
78  tntComplexArrayFree(&barr);
79 
80  /* take transpose of b to create bdag by swapping legs */
81  bdagl = tntNodeCopy(bl,0);
82  tntNodeMapLegs(bdagl,leg_map_t);
83 
84  /* n = bdag*b (U(1) invariant)
85  Make node through contracting annihilation and creation operators */
86  if (NULL != n) {
87  tnA = tntNodeCopy(bdagl,0);
88  tnB = tntNodeCopy(bl,0);
90 
91  *n = tntNodeContract(tnA, tnB, leg_map, leg_map);
92 
93  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
94  tntNodeSetQN(*n,TNT_MPS_D,&qninf,TNT_QN_IN);
96  }
97  }
98 
99  /* os_int = bdag*bdag*b*b (U(1) invariant)
100  Make node through contracting annihilation and creation operators */
101  if (NULL != os_int) {
102  /* First contract b with b */
103  tnA = tntNodeCopy(bl,0);
104  tnB = tntNodeCopy(bl,0);
105  tntNodeJoin(tnA,TNT_MPS_U,tnB,TNT_MPS_D);
106  tnB = tntNodeContract(tnA, tnB, leg_map, leg_map);
107 
108  /* Next contract bdag with bb - joining nodes such that bb is closest to MPS (i.e. the top node) */
109  tnA = tntNodeCopy(bdagl,0);
111  tnB = tntNodeContract(tnA, tnB, leg_map, leg_map);
112 
113  /* Next contract bdag with bdagbb - joining nodes such that bdagbb is closest to MPS (i.e. the top node) */
114  tnA = tntNodeCopy(bdagl,0);
115  tntNodeJoin(tnA,TNT_MPS_U,tnB,TNT_MPS_D); /* Make tnA the bottom node, tnB the top node (closet to MPS) */
116  *os_int = tntNodeContract(tnA, tnB, leg_map, leg_map);
117 
118  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
119  tntNodeSetQN(*os_int,TNT_MPS_D,&qninf,TNT_QN_IN);
120  tntNodeSetQN(*os_int,TNT_MPS_U,&qninf,TNT_QN_OUT);
121  }
122  }
123 
124  /* eye = [ 1 0]
125  [ 0 1] (U(1) invariant) */
126  if (NULL != eye) {
127  /* Allocate array */
128  eyearr = tntComplexArrayAlloc(num_vals);
129 
130  /* Put one along diagonal */
131  for (j = 0; j < physdim; j++) {
132  eyearr.vals[j*(physdim+1)].re = 1.0;
133  }
134  /* Create node from elements */
135  *eye = tntNodeCreate(&eyearr, num_legs, leg_ids, leg_dims);
136 
137  /* Free the array - it is no longer needed */
138  tntComplexArrayFree(&eyearr);
139 
140  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
141  tntNodeSetQN(*eye,TNT_MPS_D,&qninf,TNT_QN_IN);
142  tntNodeSetQN(*eye,TNT_MPS_U,&qninf,TNT_QN_OUT);
143  }
144  }
145 
146  /* Free the qn information if it was used */
147  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
148  tntIntArrayFree(&qninf);
149  }
150 
151  /* Free the nodes for b and bdag if not required, otherwise assign to pointer arguments */
152  if (NULL == b) {
153  tntNodeFree(&bl);
154  } else {
155  *b = bl;
156  }
157 
158  if (NULL == bdag) {
159  tntNodeFree(&bdagl);
160  } else {
161  *bdag = bdagl;
162  }
163 
164  return;
165 }
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:70
Simple 1D array for integer values.
Definition: tnt.h:148
tntComplex * vals
Pointer to the first element in the tntComplex array.
Definition: tnt.h:169
struct tnode * tntNode
Definition: tnt.h:118
#define TNT_SYMM_U1
Definition: tnt.h:65
tntNode tntNodeCreate(tntComplexArray *nodeVals, unsigned numLegs, unsigned legId[], unsigned legDim[])
Definition: tntNode_funcs.c:523
int tntSymmTypeGet(void)
Definition: tntMisc_funcs.c:699
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
#define TNT_QN_IN
Definition: tnt.h:71
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray_funcs.c:192
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
void tntNodeSetQN(tntNode tn, unsigned legid, tntIntArray *qvals, int legdir)
Definition: tntNode_funcs.c:1324
#define TNT_MPS_U
Definition: tntMps.h:71
Simple 1D array for complex values.
Definition: tnt.h:168
void tntNodeFree(tntNode *tn)
Definition: tntNode_funcs.c:420
double re
Definition: tnt.h:137
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:96
void tntNodeJoin(tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:462
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_QN_OUT
Definition: tnt.h:76
void tntNodeMapLegs(tntNode tn, unsigned legmap[])
Definition: tntNode_funcs.c:705
tntNetwork tntMpsCreateConfig ( unsigned  L,
const char *  config 
)

Creates an MPS network which represents a single configuration of a system. The the basis for the configuration is given by the basis operator set for the global system properties e.g. using the function tntSysBasisOpSet(), see below for more details. All internal legs will have dimension 1.

mps.png

The configuration should be expressed as a string of numbers from \(0\) to \(d-1\). If the length of the string is greater than the number of nodes in the network, then the extra characters (to the right) will be ignored. If the length of the string is smaller than the number of sites, then the rest of the sites will be filled with the first basis vector.

The MPS will have a norm of 1, and will have an orthogonality centre on every site.

The basis operator defines the properties of the physical legs, and should define the local physical eigenbasis for the system e.g. for a boson system this would be the number operator \(\hat{n}\), while for a spin system it would be the operator \(\hat{s}_z\). Operators having the correct format can be generated using the functions tntMpsCreateSpinOp() and tntMpsCreateBosonOp(), or can be loaded from an initialisation file. If a global physical symmetry has been set, then the quantum number information is also set on the physical and internal network legs. If loading your operators from an initialisation file, you must ensure that the upwards facing physical leg is labelled 4 (TNT_MPS_U), and the downwards facing physical leg is labelled 3 (TNT_MPS_D):

single_site_op.png
Returns
The network corresponding to the configuration matrix product state.
Parameters
LNumber of nodes in the MPS network
configString representing the configuration of the MPS
42 {
43 
45  tntNetwork mps; /* The mps network to create */
46  tntNode A; /* The current node */
47  tntComplexArray Aarr;
48  unsigned loop;
49  unsigned d; /* Dimension of the physical leg */
50  tntIntArray qnums_phys, qnums_int; /* Quantum numbers for the physical leg, and for the internal legs (1 internal leg only) */
51  unsigned config_length; /* Number of characters in the configuration */
52  unsigned bs; /* State for the current site */
53  unsigned numlegs = 3; /* Number of legs for a MPS */
54  unsigned legdim[3] = {1,1,1}; /* Leg dimensions for the random MPS: starting state is product state (physical leg is initialised to 1) */
55  unsigned legid[3] = {TNT_MPS_L,TNT_MPS_R,TNT_MPS_D}; /* Leg id's for the random MPS */
56  char err_msg[TNT_STRLEN]; /* character array to hold the error message when in needs to be built using local information */
57 
58  /* First create the empty network */
59  mps = tntNetworkCreate();
60 
61  /* Get the dimension of the physical leg */
62  if (NULL == tntSysBasisOpGet()) {
63  tntErrorPrint("Cannot create wave function using configuration as the basis operator for the system has not yet been set",TNT_MPS_ERRCODE); /* NO_COVERAGE */
64  } else { /* NO_COVERAGE */
65  legdim[2] = d = tntNodeGetLegDim(tntSysBasisOpGet(),TNT_MPS_D);
66  }
67 
68  if (tntSymmTypeGet()) {
69  /* Get the quantum numbers for the physical leg */
70  qnums_phys = tntNodeGetQN(tntSysBasisOpGet(), TNT_MPS_D);
71 
72  /* Allocate array for quantum numbers for internal leg. Initial value will be zero as required. */
73  qnums_int = tntIntArrayAlloc(1);
74  }
75 
76  /* Initialise array for holding node values */
77  Aarr = tntComplexArrayAlloc(d);
78 
79  /* Find the length of the provided configuration */
80  config_length = strlen(config);
81 
82  /* Keep inserting random nodes at the end. */
83  for (loop = 0; loop < L; loop++) {
84 
85  /* Find the non-zero entry using the character array, and convert the single character to an integer */
86  /* If loop is larger than character array, then fill with 0 instead */
87  if (loop < config_length) {
88  bs = (unsigned) config[loop] - '0';
89 
90  /* Check that the basis is valid */
91  if (bs > d-1) {
92  sprintf(err_msg,"Invalid configuration %s|The character %d is larger than the physical dimension (%d) allows.",config,bs,d); /* NO_COVERAGE */
93  tntErrorPrint(err_msg, TNT_MPS_ERRCODE); /* NO_COVERAGE */
94  }
95  } else {
96  bs = 0;
97  }
98 
99  /* Set the value in the array to form the node */
100  Aarr.vals[bs].re = 1.0;
101 
102  /* Create the node */
103  A = tntNodeCreate(&Aarr, numlegs, legid, legdim);
104 
105  /* Insert the node in the network */
107 
108  /* Set the array back to all zeros */
109  Aarr.vals[bs].re = 0.0;
110 
111  /* Now set the quantum number information if there are symmetries in the sysem: */
112  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
113  /* Physical leg is incoming leg and has same information as operator */
114  tntNodeSetQN(A,TNT_MPS_D,&qnums_phys,TNT_QN_IN);
115 
116  /* Left leg is incoming leg, and has value = to current total charge */
117  tntNodeSetQN(A,TNT_MPS_L,&qnums_int,TNT_QN_IN);
118 
119  /* Now increment current total charge of wave function to sum of incoming leg, and the basis state populated from physical leg */
120  qnums_int.vals[0] += qnums_phys.vals[bs];
121 
122  /* Right leg is outgoing leg, and has value to current total charge */
123  tntNodeSetQN(A,TNT_MPS_R,&qnums_int,TNT_QN_OUT);
124  }
125 
126  }
127 
128  /* Free the array for holding array values */
129  tntComplexArrayFree(&Aarr);
130 
131  /* Free the qn information if it was used */
132  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
133  tntIntArrayFree(&qnums_phys);
134  tntIntArrayFree(&qnums_int);
135  }
136 
137  /* return the completed matrix product state */
138  return mps;
139 }
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:70
Simple 1D array for integer values.
Definition: tnt.h:148
tntComplex * vals
Pointer to the first element in the tntComplex array.
Definition: tnt.h:169
struct tnode * tntNode
Definition: tnt.h:118
#define TNT_SYMM_U1
Definition: tnt.h:65
unsigned tntNodeGetLegDim(tntNode tn, unsigned legid)
Definition: tntNode_funcs.c:1069
tntNode tntNodeCreate(tntComplexArray *nodeVals, unsigned numLegs, unsigned legId[], unsigned legDim[])
Definition: tntNode_funcs.c:523
int tntSymmTypeGet(void)
Definition: tntMisc_funcs.c:699
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
void tntNodeInsertAtEnd(tntNode tnI, int legIlast, int legInwend, tntNetwork tnw)
Definition: tntNode_funcs.c:945
tntIntArray tntNodeGetQN(tntNode tn, unsigned legid)
Definition: tntNode_funcs.c:1442
#define TNT_MPS_L
Definition: tntMps.h:56
#define TNT_QN_IN
Definition: tnt.h:71
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray_funcs.c:192
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
void tntNodeSetQN(tntNode tn, unsigned legid, tntIntArray *qvals, int legdir)
Definition: tntNode_funcs.c:1324
struct tnetwork * tntNetwork
Definition: tnt.h:123
tntNetwork tntNetworkCreate(void)
Definition: tntNetwork_funcs.c:62
Simple 1D array for complex values.
Definition: tnt.h:168
#define TNT_STRLEN
Definition: tnt.h:59
double re
Definition: tnt.h:137
tntNode tntSysBasisOpGet(void)
Definition: tntMisc_funcs.c:647
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:96
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
#define TNT_QN_OUT
Definition: tnt.h:76
tntNetwork tntMpsCreateMpo ( unsigned  L,
tntNodeArray nnl,
tntNodeArray nnr,
tntComplexArray nnparam,
tntNodeArray os,
tntComplexArray osparam 
)

Creates an MPO network represeting a site-wide operator \(\hat{O}\) formed from a sum of nearest-neighbour and onsite terms.

\[ \hat{O} = \sum_{j=0}^{L-2}\sum_i^{n_n}\alpha_{i,j}\hat{o}^l_{i}\otimes\hat{o}^r_i + \sum_{j=0}^{L-1}\sum_i^{n_o}\beta_{i,j}\hat{o}^s_{i} \]

Nearest-neighbour operators \(\hat{o}^l_{i}\) and \(\hat{o}^r_i\) should be provided in arrays nnl and nnr respectively both having length \(n_n\). Onsite operators \(\hat{o}^s_{i}\) should be provided in array os having length \(n_o\). The operators should be single-site operators or product MPOs, i.e. no internal legs, and two physical legs with the legs labelled as follows:

single_site_op.png

All the operators should have the same dimension for the physical legs.

The parameters \(\alpha_{i,j}\) and \(\beta_{i,j}\) for the nearest neighbour are onsite terms are supplied in matrices nnparam and osparam respectively. The matrix must have a number of rows equal to the length \(n_n, n_o\) of its respective operators array, but there are two options for the number of columns:

  1. All the parameters are uniform accross the lattice. In this case the parameters array should have one column (which is the default behaviour if it is created with only one dimension specified). The parameter \(\alpha_{i,j}\) or \(\beta_{i,j}\) should be in position i,1 in the matrix for all sites.
  2. One or more of the parameters can vary accross the lattice. In this case the parameters matrix should have L-1 rows for nearest neighbour operators and L rows for onsite operators. The parameter \(\alpha_{i,j}\) or \(\beta_{i,j}\) for operator i and site j should be at position i,j in the matrix. Any uniform operators should have identical entries for all the sites.

A non-product site-wide MPO is then created which represents the sum of these operators, where now each operator in the network will have non-singleton dimension internal legs. The physical legs will have the same dimension as the original single-site operators, the internal legs will have a dimension equal to the number of nearest neighbour terms + 2. The legs are labelled as follows:

mpo_op_sing.png

They are connected to form the complete network:

mpo_op_nw.png
Returns
The network representing the matrix product operator.
Parameters
LLength of system.
nnlArray of nearest neighbour operators for left site. Send NULL if there are no nearest neighbour terms.
nnrArray of nearest neighbour operators for right site. Send NULL if there are no nearest neighbour terms.
nnparamArray of parameters for nearest neighbour operators. Send NULL if there are no nearest neighbour terms.
osArray of on-site operators. Send NULL if there are no on-site operators.
osparamParameters for the on-site operators. Send NULL if there are no on-site operators.
63 {
64  tntNetwork mponw; /* Network representing the MPO */
65  tntComplexArray mpogenvals; /* Values for one of the 'building blocks' of the MPO */
66  tntNode mpogen, mpobase, mpoc, eye, opss, opssterm, opnn; /* Nodes for building the MPO */
67  tntNode Lv, Rv; /* Left and right boundary vectors for the MPO */
68  unsigned numlegs = 2; /* Number of legs for mpo generator */
69  unsigned leg_id_mpo[2] = {TNT_MPS_L,TNT_MPS_R}; /* Id's of the MPO generator. One left leg, one right leg, and a leg to contract with the operator */
70  unsigned leg_dims_mpo[2] = {1,1}; /* Dimensions of the MPO generator (initialise to 1) */
71  unsigned leg_id_bv[2] = {TNT_MPS_L,TNT_MPS_R}; /* leg id for boundary vectors */
72  unsigned leg_dim_bvl[2] = {1,1}, leg_dim_bvr[2] = {1,1}; /* Dimensions of the boundary vectors (initialise to 1) */
73  unsigned Dmpo; /* Internal dimension of the MPO */
74  unsigned i, j; /* Used for looping over terms and sites respectively */
75  unsigned numnodes; /* The number of different nodes that need to be created */
76  tntComplex prm; /* The current paramter */
77  tntIntArray qnums_phys, qnums_int, qnums_op; /* Quantum numbers for the physical legs and the internal legs, and for the current operator */
78  unsigned symm_num_qn = tntSymmNumGet(); /* number of quantum numbers needed to specify a quantum number label */
79 
80  /* Check that there are the same number of left and right nearest neighbour operators */
81  if ((NULL == nnl && NULL != nnr) || (NULL != nnl && NULL == nnr) || ((NULL != nnl && NULL != nnr) && nnl->sz != nnr->sz)) {
82  tntErrorPrint("Cannot create a matrix product operator|The number of left nearest terms is not equal to the number of right nearest neighbour terms", TNT_MPS_ERRCODE); /* NO_COVERAGE */
83  } /* NO_COVERAGE */
84 
85  /* Check that there is a parameter for each of the left and right operators, or for each of the operators and each of the sites */
86  if (nnl != NULL) {
87  if (NULL == nnparam || nnl->sz != nnparam->numrows) {
88  tntErrorPrint("Cannot create a matrix product operator|The number of rows for the nearest neighbour parameters is not equal to the number of nearest neighbour terms",TNT_MPS_ERRCODE); /* NO_COVERAGE */
89  } /* NO_COVERAGE */
90  } /* NO_COVERAGE */
91 
92  /* Check that there is a parameter for each of the onsite operators */
93  if (os != NULL) {
94  if (NULL == osparam || os->sz != osparam->numrows) {
95  tntErrorPrint("Cannot create a matrix product operator|The number of rows for on-site parameters is not equal to the number of on-site terms",TNT_MPS_ERRCODE); /* NO_COVERAGE */
96  } /* NO_COVERAGE */
97  } /* NO_COVERAGE */
98 
99  /* Check that there is a non-zero number of operators */
100  if ((NULL == nnl || 0 == nnl->sz) && (NULL == os || 0 == os->sz)) {
101  tntErrorPrint("Cannot create a matrix product operator, as there are no terms",TNT_MPS_ERRCODE); /* NO_COVERAGE */
102  } /* NO_COVERAGE */
103 
104  /* Check that the length of the system is more than 1 */
105  if (L < 2) {
106  tntErrorPrint("Function to create MPO is only valid for a system length of two or more", TNT_MPS_ERRCODE); /* NO_COVERAGE */
107  } /* NO_COVERAGE */
108 
109  /* Calulcate how many nodes to create: one node if the operators are uniform, L nodes if any of the operators vary */
110  if ((NULL == os || 0 == os->sz || 1 == osparam->numcols) && (NULL == nnl || 0 == nnl->sz || 1 == nnparam->numcols)) {
111  numnodes = 1;
112  } else {
113  numnodes = L;
114  }
115 
116  /* Create the empty network that will contain all the MPOs and that will be returned by the function. */
117  mponw = tntNetworkCreate();
118 
119  /* Create the identity node basing it on the supplied single-site operators */
120  if (NULL != nnl && 0 != nnl->sz)
121  eye = tntNodeCreateEyeOp(nnl->vals[0]);
122  else
123  eye = tntNodeCreateEyeOp(os->vals[0]);
124 
125  /* The internal dimension of the MPO will be the number of nn terms + 2. */
126  Dmpo = (NULL == nnl) ? 2 : nnl->sz + 2;
127  leg_dims_mpo[0] = leg_dims_mpo[1] = Dmpo;
128  leg_dim_bvl[1] = leg_dim_bvr[0] = Dmpo;
129 
130  /* ----------------- SETTING QN INFO HERE ---------------- */
131  /* If symmetries are being used, allocate memory for the lists of quantum numbers for the internal legs, and assign the values using the supplied operators */
132  if (tntSymmTypeGet()) {
133 
134  /* Get the quantum numbers for the physical leg */
135  qnums_phys = tntNodeGetQN(tntSysBasisOpGet(), TNT_MPS_D);
136 
137  /* Allocate array of the correct size for the internal leg quantum numbers - all the values will be initialised to zero */
138  qnums_int = tntIntArrayAlloc(Dmpo*symm_num_qn);
139 
140  /* Loop through the nearest neighbour operators, finding the correct internal quantum number for each one */
141  if (nnl != NULL && 0 != nnl->sz) {
142  /* For each operator, find the quantum number that should be added to a new outgoing leg that will make the number of blocks = 1 */
143 
144  for (i = 0; i < nnl->sz; i++) {
145  /* First check whether the node is already covariant, if so no action required as required value for the internal leg is zero (initialised value) */
146  if (!tntNodeIsCovariant(nnr->vals[i])) {
147 
148  /* make a copy of the operator required */
149  /* Note: using right operator since the function tntNodeMakeCovariant() calculates QN for outgoing leg, and the right operator will be associated with the outgoing leg */
150  opnn = tntNodeCopy(nnr->vals[i],0);
151 
152  /* add a leg to it */
153  tntNodeAddLeg(opnn, 5);
154 
155  /* assign quantum numbers to the physical legs */
156  tntNodeSetQN(opnn,TNT_MPS_D,&qnums_phys,TNT_QN_IN);
157  tntNodeSetQN(opnn,TNT_MPS_U,&qnums_phys,TNT_QN_OUT);
158 
159  /* assign quantum number label to the remaining singleton leg, where the leg is chosen to make the node covariant */
161 
162  /* get the quantum number label */
163  qnums_op = tntNodeGetQN(opnn, 5);
164 
165  /* copy it to the relevant entry for the internal leg */
166  for (j = 0; j < symm_num_qn; j++) {
167  qnums_int.vals[(i+1)*symm_num_qn + j] = qnums_op.vals[j];
168  }
169 
170  /* free the array containing the quantum number label */
171  tntIntArrayFree(&qnums_op);
172 
173  /* free the copy of the operator */
174  tntNodeFree(&opnn);
175  }
176 
177  }
178  }
179 
180  }
181  /* ----------------- END OF SETTING QN INFO ----------------- */
182 
183  /* Create the generator for the MPO that doesn't depend on operators i.e. identity in first and last element */
184  mpogenvals = tntComplexArrayAlloc(Dmpo*Dmpo);
185  mpogenvals.vals[0].re = mpogenvals.vals[Dmpo*Dmpo - 1].re = 1.0;
186  mpogen = tntNodeCreate(&mpogenvals, 2, leg_id_mpo, leg_dims_mpo);
187 
188  /* Contract the gen with the identity operator to form the base mpo term */
189  /* Note that the contraction step will free both mpogen and eye */
190  mpobase = tntNodeContract(mpogen, eye, NULL, NULL);
191 
192  /* reset the numbers for the generator values */
193  mpogenvals.vals[0].re = mpogenvals.vals[Dmpo*Dmpo - 1].re = 0.0;
194 
195  /* Create 1 copy of the MPO base in the network if the operators are uniform, and L copies of the MPO base in the network
196  if any one of the operators vary through the system */
197  /* Insert first node in the network */
198  tntNodeInsertAtEnd(mpobase, TNT_MPS_L, TNT_MPS_R, mponw);
199 
200  /* Insert any remaining nodes in the network */
201  for (j = 1; j < numnodes; j++) {
202 
203  /* Make a copy of the node to insert */
204  mpoc = tntNodeCopy(mpobase,0);
205 
206  /* Insert the node in the network */
207  tntNodeInsertAtEnd(mpoc, TNT_MPS_L, TNT_MPS_R, mponw);
208 
209  }
210 
211  /* Now create the on-site terms: Create a different term for each site if the parameters vary. */
212  if (os != NULL && 0 != os->sz) {
213 
214  /* Populate the bottom-left element of the MPO generator */
215  mpogenvals.vals[Dmpo-1].re = 1.0;
216 
217  /* Create the node */
218  mpogen = tntNodeCreate(&mpogenvals, numlegs, leg_id_mpo, leg_dims_mpo);
219 
220  /* Reset the array for the MPO generator */
221  mpogenvals.vals[Dmpo-1].re = 0.0;
222 
223  for (j = 0; j < numnodes; j++) {
224 
225  /* Let the onsite term be equal to the first operator */
226  opss = tntNodeCopy(os->vals[0],0);
227 
228  /* Set the parameters depending on whether the on-site parameters also vary */
229  if (1 == osparam->numcols) {
230  prm.re = osparam->vals[0].re;
231  prm.im = osparam->vals[0].im;
232  } else {
233  prm.re = osparam->vals[j * os->sz].re;
234  prm.im = osparam->vals[j * os->sz].im;
235  }
236 
237  /* Scale the operator by the parameter */
238  tntNodeScaleComplex(opss, prm);
239 
240  /* Create each addition term, and add it to the total on-site operator */
241  for (i = 1; i < os->sz; i++) {
242  /* Copy operator */
243  opssterm = tntNodeCopy(os->vals[i],0);
244 
245  /* Set the parameters depending on whether the on-site parameters also vary */
246  if (1 == osparam->numcols) {
247  prm.re = osparam->vals[i].re;
248  prm.im = osparam->vals[i].im;
249  } else {
250  prm.re = osparam->vals[i + j * os->sz].re;
251  prm.im = osparam->vals[i + j * os->sz].im;
252  }
253 
254  /* Determine the correct parameter and scale */
255  tntNodeScaleComplex(opssterm, prm);
256 
257  /* Add to the total term */
258  tntNodeAdd(opss,opssterm);
259 
260  /* Free as it will no longer be required */
261  tntNodeFree(&opssterm);
262  }
263 
264  /* Create a copy of the MPO generator */
265  mpoc = tntNodeCopy(mpogen, 0);
266 
267  /* Contract the MPO generator with the on-site operator */
268  mpoc = tntNodeContract(mpoc, opss, NULL, NULL);
269 
270  /* Add the MPO term to the total MPO */
271  tntNodeAdd(mpobase, mpoc);
272 
273  /* free this generating node as it is no longer required */
274  tntNodeFree(&mpoc);
275 
276  /* Find the next mpobase term */
277  mpobase = tntNodeFindConn(mpobase, TNT_MPS_R);
278 
279  }
280 
281  /* Free the generating node */
282  tntNodeFree(&mpogen);
283 
284  }
285 
286  /* Now create the nearest neighbour terms */
287  if (nnl != NULL && 0 != nnl->sz) {
288  /* Create each term, using the MPO generator to insert each nn term to a different position in the MPO tensor. */
289  for (i = 0; i < nnl->sz; i++) {
290 
291  /* Make the MPO generator: position of left operator will be bottom row, and (i + 1)st column (if col numbering starts from zero) */
292  mpogenvals.vals[Dmpo-1 + (i+1)*Dmpo].re = 1.0;
293  mpogen = tntNodeCreate(&mpogenvals, 2, leg_id_mpo, leg_dims_mpo);
294 
295  /* Reset the array for the MPO generator */
296  mpogenvals.vals[Dmpo-1 + (i+1)*Dmpo].re = 0.0;
297 
298  /* Find the first MPO in the network */
299  mpobase = tntNodeFindFirst(mponw);
300 
301  /* Now loop through sites for the left operator of the pair.
302  Note there are L MPOs, but only L-1 parameters for the case of
303  varying parameters. Put the parameter on the L node, which although built for the last site, will be
304  destroyed by the boundary node. Therefore for the last node do not scale */
305 
306  for (j = 0; j < numnodes; j++) {
307  /* Copy the left operator and scale by parameter */
308  opnn = tntNodeCopy(nnl->vals[i],0);
309 
310  /* Determine the correct parameter */
311  if (j == L-1) {
312  prm.re = 1.0;
313  prm.im = 0.0;
314  } else {
315  /* Set the parameters depending on whether the nearest neighbour parameters also vary */
316  if (1 == nnparam->numcols) {
317  prm.re = nnparam->vals[i].re;
318  prm.im = nnparam->vals[i].im;
319  } else {
320  prm.re = nnparam->vals[i + j*nnl->sz].re;
321  prm.im = nnparam->vals[i + j*nnl->sz].im;
322  }
323  }
324 
325  tntNodeScaleComplex(opnn, prm);
326 
327  /* Make a copy of the MPO generator */
328  mpoc = tntNodeCopy(mpogen, 0);
329 
330  /* Contract the MPO generator and the single site operator */
331  mpoc = tntNodeContract(mpoc, opnn, NULL, NULL);
332 
333  /* Add the MPO term to the total MPO */
334  tntNodeAdd(mpobase, mpoc);
335 
336  /* free this generating node as it is no longer required */
337  tntNodeFree(&mpoc);
338 
339  /* Find the next MPO in the network */
340  mpobase = tntNodeFindConn(mpobase,TNT_MPS_R);
341 
342  }
343 
344  /* free this generating node as it is no longer required */
345  tntNodeFree(&mpogen);
346 
347  /* Make the MPO generator: position of right operator will be first column, and (i + 1)st row (if row numbering starts from zero) */
348  mpogenvals.vals[i+1].re = 1.0;
349  mpogen = tntNodeCreate(&mpogenvals, numlegs, leg_id_mpo, leg_dims_mpo);
350 
351  /* Reset the array for the MPO generator */
352  mpogenvals.vals[i+1].re = 0.0;
353 
354  /* Find the first MPO in the network */
355  mpobase = tntNodeFindFirst(mponw);
356 
357  /* Now loop through the sites for the right node of the pair. These
358  do not need to be scaled */
359  for (j = 0; j < numnodes; j++) {
360 
361  /* Copy the right operator (this does not need to be scaled by the parameter) */
362  opnn = tntNodeCopy(nnr->vals[i],0);
363 
364  /* Make a copy of the MPO generator */
365  mpoc = tntNodeCopy(mpogen, 0);
366 
367  /* Contract the MPO generator and the single site operator */
368  mpoc = tntNodeContract(mpoc, opnn, NULL, NULL);
369 
370  /* Add the MPO term to the total MPO */
371  tntNodeAdd(mpobase, mpoc);
372 
373  /* free this generating node as it is no longer required */
374  tntNodeFree(&mpoc);
375 
376  /* Find the next MPO in the network */
377  mpobase = tntNodeFindConn(mpobase,TNT_MPS_R);
378 
379  }
380 
381  /* free this generating node as it is no longer required */
382  tntNodeFree(&mpogen);
383  }
384  }
385 
386  /* ----------------- SETTING QN INFO HERE ---------------- */
387  /* Now put all the operators in the network in blocks form */
388  if (tntSymmTypeGet()) {
389 
390  mpoc = tntNodeFindFirst(mponw);
391 
392  /* Set the quantum numbers for the first MPO. */
393  tntNodeSetQN(mpoc,TNT_MPS_D,&qnums_phys,TNT_QN_IN);
394  tntNodeSetQN(mpoc,TNT_MPS_U,&qnums_phys,TNT_QN_OUT);
395  tntNodeSetQN(mpoc,TNT_MPS_L,&qnums_int,TNT_QN_IN);
396  tntNodeSetQN(mpoc,TNT_MPS_R,&qnums_int,TNT_QN_OUT);
397 
398  while (mpoc != tntNodeFindLast(mponw)) {
399 
400  /* Move on the the next node */
401  mpoc = tntNodeFindConn(mpoc,TNT_MPS_R);
402 
403  /* Set the quantum numbers for the MPO. */
404  tntNodeSetQN(mpoc,TNT_MPS_D,&qnums_phys,TNT_QN_IN);
405  tntNodeSetQN(mpoc,TNT_MPS_U,&qnums_phys,TNT_QN_OUT);
406  tntNodeSetQN(mpoc,TNT_MPS_L,&qnums_int,TNT_QN_IN);
407  tntNodeSetQN(mpoc,TNT_MPS_R,&qnums_int,TNT_QN_OUT);
408  }
409 
410 
411  /* free the quantum numbers */
412  tntIntArrayFree(&qnums_phys);
413  }
414  /* ----------------- END OF SETTING QN INFO ----------------- */
415 
416  /* Now if only one node has been created, create additional L-1 identical copies to put in the network */
417  if (1 == numnodes) {
418 
419  mpobase = tntNodeFindFirst(mponw);
420 
421  /* Keep inserting mpo nodes at the end. */
422  for (j = 1; j < L; j++) {
423 
424  /* Make a copy of the node to insert */
425  mpoc = tntNodeCopy(mpobase,0);
426 
427  /* Insert the node in the network */
428  tntNodeInsertAtEnd(mpoc, TNT_MPS_L, TNT_MPS_R, mponw);
429 
430  }
431  }
432 
433  /* Create the leftmost and rightmost MPO nodes by contracting with the boundary nodes */
434  /* Elements for the boundary nodes - reallocate array of values */
435  tntComplexArrayFree(&mpogenvals);
436  mpogenvals = tntComplexArrayAlloc(Dmpo);
437 
438  /* Left boundary vector: Want to pick out the final row, so choose row vector with last element 1 */
439  mpogenvals.vals[Dmpo-1].re = 1.0;
440  Lv = tntNodeCreate(&mpogenvals,numlegs,leg_id_bv,leg_dim_bvl);
441 
442  /* Right boundary vector: Want to pick out the first row, so choose column vector with the first element 1 */
443  mpogenvals.vals[Dmpo-1].re = 0.0;
444  mpogenvals.vals[0].re = 1.0;
445  Rv = tntNodeCreate(&mpogenvals,numlegs,leg_id_bv,leg_dim_bvr);
446 
447  /* ----------------- SETTING QN INFO HERE ---------------- */
448  /* Now put the boundary vectors in blocks form */
449  if (tntSymmTypeGet()) {
450 
451  mpoc = tntNodeFindFirst(mponw);
452 
453  /* Set the quantum numbers for the left boundary vector. */
454  tntNodeSetQN(Lv,TNT_MPS_L,&qnums_int,TNT_QN_IN);
455  tntNodeSetQN(Lv,TNT_MPS_R,&qnums_int,TNT_QN_OUT);
456 
457  /* Set the quantum numbers for the right boundary vector. */
458  tntNodeSetQN(Rv,TNT_MPS_L,&qnums_int,TNT_QN_IN);
459  tntNodeSetQN(Rv,TNT_MPS_R,&qnums_int,TNT_QN_OUT);
460 
461  /* Now quantum numbers for internal legs are finished with - free the array */
462  tntIntArrayFree(&qnums_int);
463 
464  }
465  /* ----------------- END OF SETTING QN INFO ----------------- */
466 
467  /* Insert the left and right boundary vectors in the network */
470 
471  /* Contract left and right boundary vectors their neighbouring MPOs */
472  mpoc = tntNodeFindConn(Lv, TNT_MPS_R);
473  Lv = tntNodeContract(Lv, mpoc, NULL, NULL);
474 
475  mpoc = tntNodeFindConn(Rv, TNT_MPS_L);
476  Rv = tntNodeContract(Rv, mpoc, NULL, NULL);
477 
478  /* Free the arrays and nodes that are no longer required */
479  tntComplexArrayFree(&mpogenvals);
480 
481  /* Return the MPO network */
482  return mponw;
483 }
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:70
Simple 1D array for integer values.
Definition: tnt.h:148
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
tntComplex * vals
Pointer to the first element in the tntComplex array.
Definition: tnt.h:169
struct tnode * tntNode
Definition: tnt.h:118
void tntNodeScaleComplex(tntNode tn, tntComplex val)
Definition: tntNode_funcs.c:646
tntNode tntNodeFindLast(tntNetwork tnw)
Definition: tntNode_funcs.c:685
tntNode tntNodeCreate(tntComplexArray *nodeVals, unsigned numLegs, unsigned legId[], unsigned legDim[])
Definition: tntNode_funcs.c:523
int tntSymmTypeGet(void)
Definition: tntMisc_funcs.c:699
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
void tntNodeAddLeg(tntNode tn, unsigned leglabel)
Definition: tntNode_funcs.c:401
void tntNodeInsertAtEnd(tntNode tnI, int legIlast, int legInwend, tntNetwork tnw)
Definition: tntNode_funcs.c:945
tntNode tntNodeCreateEyeOp(tntNode tn)
Definition: tntNode_funcs.c:270
tntIntArray tntNodeGetQN(tntNode tn, unsigned legid)
Definition: tntNode_funcs.c:1442
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
#define TNT_MPS_L
Definition: tntMps.h:56
unsigned tntSymmNumGet(void)
Definition: tntMisc_funcs.c:713
#define TNT_QN_IN
Definition: tnt.h:71
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray_funcs.c:192
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
void tntNodeInsertAtStart(tntNode tnI, int legInwstart, int legIfirst, tntNetwork tnw)
Definition: tntNode_funcs.c:905
unsigned tntNodeIsCovariant(tntNode tn)
Definition: tntNode_funcs.c:1256
void tntNodeSetQN(tntNode tn, unsigned legid, tntIntArray *qvals, int legdir)
Definition: tntNode_funcs.c:1324
struct tnetwork * tntNetwork
Definition: tnt.h:123
#define TNT_MPS_U
Definition: tntMps.h:71
unsigned numrows
The number of rows in the matrix.
Definition: tnt.h:171
tntNetwork tntNetworkCreate(void)
Definition: tntNetwork_funcs.c:62
Simple 1D array for complex values.
Definition: tnt.h:168
tntNode * vals
Pointer to the first element in the tntNode array.
Definition: tnt.h:179
void tntNodeFree(tntNode *tn)
Definition: tntNode_funcs.c:420
double re
Definition: tnt.h:137
tntNode tntSysBasisOpGet(void)
Definition: tntMisc_funcs.c:647
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:96
#define TNT_MPS_D
Definition: tntMps.h:66
void tntNodeMakeCovariantQN(tntNode tn)
Definition: tntNode_funcs.c:1465
#define TNT_MPS_R
Definition: tntMps.h:61
double im
Definition: tnt.h:138
#define TNT_QN_OUT
Definition: tnt.h:76
void tntNodeAdd(tntNode tnA, tntNode tnB)
Definition: tntNode_funcs.c:385
unsigned numcols
The number of columns in the matrix.
Definition: tnt.h:172
Use this type to define complex numbers. Any complex arithmetic needs to be done manually. For example to add two complex numbers A and B to give C write:
Definition: tnt.h:136
unsigned sz
The number of elements in the array.
Definition: tnt.h:180
tntNetwork tntMpsCreateRandom ( unsigned  L,
unsigned  D 
)

Creates an MPS network with random tensor elements and with maximum internal dimension \(D\) having open boundary conditions. The the basis for the configuration is given by the basis operator set for the global system properties e.g. using the function tntSysBasisOpSet(), see below for more details.

mps.png

The MPS will have a norm of 1, and will have an orthogonality centre on every site.

The basis operator defines the properties of the physical legs, and should define the local physical eigenbasis for the system e.g. for a boson system this would be the number operator \(\hat{n}\), while for a spin system it would be the operator \(\hat{s}_z\). Operators having the correct format can be generated using the functions tntMpsCreateSpinOp() and tntMpsCreateBosonOp(), or can be loaded from an initialisation file. If loading your operators from an initialisation file, you must ensure that the upwards facing physical leg is labelled 4 (TNT_MPS_U), and the downwards facing physical leg is labelled 3 (TNT_MPS_D):

single_site_op.png
Returns
The network corresponding to the configuration matrix product state.
Parameters
LNumber of nodes in the MPS network
DMaximum internal dimension for the state
36 {
37 
38  tntNetwork mps; /* The mps network to create */
39  tntNode A; /* The current node */
40  unsigned d; /* Dimension of the physical leg */
41  unsigned loop; /* used for looping */
42  unsigned numlegs = 3; /* Number of legs for a MPS */
43  unsigned legdim[3] = {D,D,1}; /* Leg dimensions for the random MPS (physical leg initialised to 1) */
44  unsigned legid[3] = {TNT_MPS_L,TNT_MPS_R,TNT_MPS_D}; /* Leg id's for the random MPS */
45 
46  /* First create the empty network */
47  mps = tntNetworkCreate();
48 
49  /* Get the dimension of the physical leg */
50  if (NULL == tntSysBasisOpGet()) {
51  tntErrorPrint("Cannot create random wave function as the basis operator for the system has not yet been set",TNT_MPS_ERRCODE); /* NO_COVERAGE */
52  } else { /* NO_COVERAGE */
53  legdim[2] = d = tntNodeGetLegDim(tntSysBasisOpGet(),TNT_MPS_D);
54  }
55 
56  /* First left leg will have dimension 1 */
57  legdim[0] = 1;
58 
59  /* Keep inserting random nodes at the end. */
60  for (loop = 0; loop < L; loop++) {
61 
62  /* First set the dimensions of the node: check whether the `exact' dimension is less than D (true at edges if d<D), and set appropriately */
63  if (loop <= L/2) {
64  if (pow(d,loop+1) > D || HUGE_VAL == pow(d,loop+1)) {
65  if (pow(d,L-loop-1) > D || HUGE_VAL == pow(d,L-loop-1))
66  legdim[1] = D;
67  else
68  legdim[1] = (unsigned) pow(d,L-loop-1); /* This line will be reached when there is an odd number of sites,
69  and D is between d^(L-1)/2 and d^(L+1)/2. */
70  }
71  else
72  legdim[1] = (unsigned) pow(d,loop+1);
73  } else {
74  if (pow(d,L-loop-1) > D || HUGE_VAL == pow(d,L-loop-1))
75  legdim[1] = D;
76  else
77  legdim[1] = (unsigned) pow(d,L-loop-1);
78  }
79 
80 
81  /* Create the node */
82  A = tntNodeCreate(NULL, numlegs, legid, legdim);
83 
84  if (loop) {
85  /* Take the SVD (if it is not the first node) to generate a random orthogonalised tensor */
86  /* Node will be created as a right orthonormalised node i.e. legs 2 and 3 assigned to U (which is node A), also connected to S and VT */
87  A = tntNodeSVD(A, (int *) legid+1, 2, TNT_MPS_L, TNT_MPS_R, NULL, -1, NULL);
88  }
89 
90  /* Copy the A node only into the network */
92 
93  /* Free the node group consisting of U (A), S and VT */
94  tntNodeGroupFree(&A);
95 
96  /* Next left leg has to have the same dimension as the previous right leg */
97  legdim[0] = legdim[1];
98  }
99 
100  return mps;
101 }
struct tnode * tntNode
Definition: tnt.h:118
unsigned tntNodeGetLegDim(tntNode tn, unsigned legid)
Definition: tntNode_funcs.c:1069
tntNode tntNodeCreate(tntComplexArray *nodeVals, unsigned numLegs, unsigned legId[], unsigned legDim[])
Definition: tntNode_funcs.c:523
void tntNodeInsertAtEnd(tntNode tnI, int legIlast, int legInwend, tntNetwork tnw)
Definition: tntNode_funcs.c:945
#define TNT_MPS_L
Definition: tntMps.h:56
tntNode tntNodeSVD(tntNode tn, int *startLegs, unsigned startNum, int legUS, int legSU, int *legMap, int connDim, double *err)
Definition: tntNode_funcs.c:64
void tntNodeGroupFree(tntNode *tn)
Definition: tntNode_funcs.c:1053
struct tnetwork * tntNetwork
Definition: tnt.h:123
tntNetwork tntNetworkCreate(void)
Definition: tntNetwork_funcs.c:62
tntNode tntSysBasisOpGet(void)
Definition: tntMisc_funcs.c:647
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
void tntMpsCreateSpinOp ( unsigned  TwoS,
tntNode Sx,
tntNode Sy,
tntNode Sz,
tntNode Sp,
tntNode Sm,
tntNode eye 
)

Creates single site operators for spin systems, with two legs (i.e. no internal legs), numbered to match those expected by the MPS library functions:

single_site_op.png

Each leg will have dimension \(2S+1\), where \(S\) is the spin of the particle on each site.

All nodes are created as static nodes.

Pointers to the required nodes should be passed as arguments. If any of the nodes are not required, then pass NULL.

Returns
No return value.
Parameters
TwoSTwice the spin \(S\) of the system.
SxNode for \(\hat{S}^x\)
SyNode for \(\hat{S}^y\)
SzNode for \(\hat{S}^z\)
SpNode for \(\hat{S}^+ = \hat{S}^x + \mathrm{i}\hat{S}^y\)
SmNode for \(\hat{S}^- = \hat{S}^x - \mathrm{i}\hat{S}^y\)
eyeNode for the identity matrix with dimension \(2S+1\)
36 {
37 
38  tntComplexArray Szarr, Sparr, Smarr, eyearr; /* Arrays to hold the tensor elements */
39  tntNode Spl = NULL, Sml = NULL; /* Local spin raising and lowering operators used to create Sx and Sy */
40  double m; /* z component of the spin */
41  double S = (double) TwoS*0.5; /* Spin of the system */
42  unsigned physdim = TwoS + 1; /* Physical dimension for the spin system */
43  unsigned num_legs = 2, num_vals = physdim*physdim; /* Define array properties. Adding on an extra leg for contracting to make two site gates or to make MPO */
44  unsigned leg_ids[2] = {3,4}; /* Legs for the operators */
45  unsigned leg_dims[2] = {physdim,physdim}; /* the physical legs */
46  unsigned j; /* index for arrays */
47  tntComplex ScaleSy = {0, 0.5}; /* tntComplex value to scale Sy by */
48  tntIntArray qninf; /* Quantum number information for the physical legs */
49 
50  /* Quantum number information to set for invariant operators only */
51  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
52  qninf = tntIntArrayAlloc(physdim);
53  for (j = 0; j < physdim; j++) {
54  qninf.vals[j] = -2*j-1+physdim;
55  }
56  }
57 
58  /* Check whether each node pointer is NULL, before creating each node */
59 
60  /* Sz = [ S 0 ... 0 0]
61  [ 0 S-1 ... 0 0]
62  [ : ]
63  [ 0 0 ... -S+1 0]
64  [ 0 0 ... 0 -S]
65  */
66  if (NULL != Sz) {
67 
68  /* Make array of values */
69  Szarr = tntComplexArrayAlloc(num_vals);
70  for (j = 0; j < physdim; j++) {
71  /* Calculate the value of the z componenet of the spin for each element */
72  m = S - j;
73  /* Put the value of the spin on the diagonal elements */
74  Szarr.vals[j*(physdim+1)].re = m;
75  }
76  *Sz = tntNodeCreate(&Szarr, num_legs, leg_ids, leg_dims);
77 
78  /* Free the array - it is no longer needed */
79  tntComplexArrayFree(&Szarr);
80 
81  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
82  tntNodeSetQN(*Sz,TNT_MPS_D,&qninf,TNT_QN_IN);
83  tntNodeSetQN(*Sz,TNT_MPS_U,&qninf,TNT_QN_OUT);
84  }
85  }
86 
87  /* Sp = [ 0 c_{S} 0 0 0 ...]
88  [ 0 0 c_{S-1} 0 0 ...]
89  [ 0 0 0 c_{-S+3} 0 ...]
90  [ 0 0 0 0 c_{-S+4} ...]
91  [ 0 0 0 0 0 ...]
92  [ ..... ]
93  with c_m = sqrt((S-m)(S-m+1))*/
94 
95  if ((NULL != Sp)||(NULL != Sx)||(NULL != Sy)) {
96 
97  /* Make array of values */
98  Sparr = tntComplexArrayAlloc(num_vals);
99  for (j = 1; j < physdim; j++) {
100  /* Calculate the value of the z componenet of the spin for each element */
101  m = S - j;
102  /* Put the value of coefficient for spin j on row j-1, column j */
103  Sparr.vals[(j-1) + j*physdim].re = sqrt((S-m)*(S+m+1.0));
104  }
105  Spl = tntNodeCreate(&Sparr, num_legs, leg_ids, leg_dims);
106 
107  /* Free the array - it is no longer needed */
108  tntComplexArrayFree(&Sparr);
109  }
110 
111 
112  if ((NULL != Sm)||(NULL != Sx)||(NULL != Sy)) {
113 
114  /* Make array of values */
115  Smarr = tntComplexArrayAlloc(num_vals);
116  for (j = 0; j < (physdim-1); j++) {
117  /* Calculate the value of the z componenet of the spin for each element */
118  m = S - j;
119  /* Put the value of coefficient for spin j on row j+1, column j */
120  Smarr.vals[(j+1) + j*physdim].re = sqrt((S+m)*(S-m+1.0));
121  }
122  Sml = tntNodeCreate(&Smarr, num_legs, leg_ids, leg_dims);
123 
124  /* Free the array - it is no longer needed */
125  tntComplexArrayFree(&Smarr);
126  }
127 
128  /* \sigma^x = 1/2*(Sp + Sm)*/
129  if (NULL != Sx) {
130  /* Make Sx a copy of the local version of Sp */
131  *Sx = tntNodeCopy(Spl,0);
132  /* Add the local copy of Sm */
133  tntNodeAdd(*Sx,Sml);
134  /* Scale the total node by 1/2 */
135  tntNodeScaleReal(*Sx,0.5);
136  }
137 
138  /* \sigma^y = -i/2*(Sp - Sm)*/
139  if (NULL != Sy) {
140  /* Make Sy a copy of the local version of Sp */
141  *Sy = tntNodeCopy(Spl,0);
142  /* Scale by -1 */
143  tntNodeScaleReal(*Sy,-1);
144  /* Add the local copy of Sm */
145  tntNodeAdd(*Sy,Sml);
146  /* Scale the total node by i/2 */
147  tntNodeScaleComplex(*Sy,ScaleSy);
148  }
149 
150  /* create a new version of sigma plus if it is required by user */
151  if (NULL != Sp) {
152  *Sp = tntNodeCopy(Spl,0);
153  }
154 
155  /* create a new version of sigma minus if it is required by user */
156  if (NULL != Sm) {
157  *Sm = tntNodeCopy(Sml,0);
158  }
159 
160  /* eye = [ 1 0 ....]
161  [ 0 1 ....]
162  .....
163  [ .... 0 1] */
164  if (NULL != eye) {
165  eyearr = tntComplexArrayAlloc(num_vals);
166  /* Put ones along diagonal */
167  for (j = 0; j < physdim; j++) {
168  eyearr.vals[j*(physdim+1)].re = 1.0;
169  }
170  *eye = tntNodeCreate(&eyearr, num_legs, leg_ids, leg_dims);
171 
172  /* Free the array - it is no longer needed */
173  tntComplexArrayFree(&eyearr);
174 
175  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
176  tntNodeSetQN(*eye,TNT_MPS_D,&qninf,TNT_QN_IN);
177  tntNodeSetQN(*eye,TNT_MPS_U,&qninf,TNT_QN_OUT);
178  }
179  }
180 
181  /* Free the qn information if it was used */
182  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
183  tntIntArrayFree(&qninf);
184  }
185 
186  /* Free local arrays used to build operators */
187  tntNodeFree(&Spl);
188  tntNodeFree(&Sml);
189 
190  return;
191 
192 }
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:70
Simple 1D array for integer values.
Definition: tnt.h:148
tntComplex * vals
Pointer to the first element in the tntComplex array.
Definition: tnt.h:169
struct tnode * tntNode
Definition: tnt.h:118
#define TNT_SYMM_U1
Definition: tnt.h:65
void tntNodeScaleComplex(tntNode tn, tntComplex val)
Definition: tntNode_funcs.c:646
tntNode tntNodeCreate(tntComplexArray *nodeVals, unsigned numLegs, unsigned legId[], unsigned legDim[])
Definition: tntNode_funcs.c:523
int tntSymmTypeGet(void)
Definition: tntMisc_funcs.c:699
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
#define TNT_QN_IN
Definition: tnt.h:71
void tntNodeScaleReal(tntNode tn, double val)
Definition: tntNode_funcs.c:623
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray_funcs.c:192
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
void tntNodeSetQN(tntNode tn, unsigned legid, tntIntArray *qvals, int legdir)
Definition: tntNode_funcs.c:1324
#define TNT_MPS_U
Definition: tntMps.h:71
Simple 1D array for complex values.
Definition: tnt.h:168
void tntNodeFree(tntNode *tn)
Definition: tntNode_funcs.c:420
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:96
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_QN_OUT
Definition: tnt.h:76
void tntNodeAdd(tntNode tnA, tntNode tnB)
Definition: tntNode_funcs.c:385
Use this type to define complex numbers. Any complex arithmetic needs to be done manually. For example to add two complex numbers A and B to give C write:
Definition: tnt.h:136
tntNetwork tntMpsCreateSymmRandom ( unsigned  L,
int *  qn 
)

Creates an MPS network with random tensor elements and with maximum internal dimension \(D\) having open boundary conditions, and with a given total quantum number. If there is no symmetry type set for the system, the quantum number is ignored, and a random MPS with no conservation is created, i.e. this function simply becomes equivalent to tntMpsCreateRandom(). Note: this function only currently works for the U(1) symmetry type. The the basis for the configuration is given by the basis operator set for the global system properties e.g. using the function tntSysBasisOpSet(), see below for more details.

mps.png

The MPS will have a norm of 1, and will have an orthogonality centre on every site.

The basis operator defines the properties of the physical legs, and should define the local physical eigenbasis for the system e.g. for a boson system this would be the number operator \(\hat{n}\), while for a spin system it would be the operator \(\hat{s}_z\). Operators having the correct format can be generated using the functions tntMpsCreateSpinOp() and tntMpsCreateBosonOp(), or can be loaded from an initialisation file. If loading your operators from an initialisation file, you must ensure that the upwards facing physical leg is labelled 4 (TNT_MPS_U), and the downwards facing physical leg is labelled 3 (TNT_MPS_D):

single_site_op.png
Returns
The network corresponding to the configuration matrix product state.
Parameters
LNumber of nodes in the MPS
qnArray containing the total quantum number label for the MPS. Size of array must be symm_num_qn.
38 {
39  tntNetwork mps; /* The mps network to create */
40  tntNode A; /* The current node */
41  unsigned d; /* Dimension of the physical leg */
42  unsigned loop; /* used for looping */
43  unsigned numlegs = 3; /* Number of legs for a MPS */
44  unsigned legdim[3] = {1,1,1}; /* Leg dimensions for the random MPS (physical leg initialised to 1) */
45  unsigned legid[3] = {TNT_MPS_L,TNT_MPS_R,TNT_MPS_D}; /* Leg id's for the random MPS */
46  double norm_sq; /* The norm of the MPS */
47  unsigned i, j, k, qind; /* used for looping and indexing*/
48  tntIntArray dqn; /* difference between succesive quantum numbers */
49  tntIntArray qmin, qmax; /* minimum and maximum quantum numbers for each quantum number label on the physical leg */
50  tntIntArray qminsum, qmaxsum; /* minimum and maximum possible quantum numbers for the internal leg, to allow the boundary qn's to be reached */
51  tntIntArray qnums_phys; /* Quantum numbers for the physical leg */
52  tntIntArray qnums_int[2]; /* Array of list of quantum numbers for left and right legs, where the parity of the site will be used to determine which of the quantum numbers is for the left and which is for the right site */
53  unsigned numqn, lblnum = tntSymmNumGet(); /* current number of qns, number of qns in a label */
54  unsigned qn_valid = 1, qn_appeared = 0; /* flag to say whether current qn label is allowed, and flag to say whether current qn has already appeared */
55  char err_msg[TNT_STRLEN]; /* Character array to hold the error message when it needs to be created using local information */
56  int numsites = L; /* number if sites as integer variable */
57 
58  /* Check there is a symmetry type set for the system */
59  if (!tntSymmTypeGet()) {
60  /* if there is no symmetry type set, simply create an unsymmetric MPS after printing warning */
61  tntWarningPrint("Trying to create random symmetric MPS but no symmetry type set|Unsymmetric random MPS with D=2 being created instead."); /* NO_COVERAGE */
62  mps = tntMpsCreateRandom(L,2); /* NO_COVERAGE */
63  return mps; /* NO_COVERAGE */
64  } /* NO_COVERAGE */
65 
66  /* Get the dimension of the physical leg */
67  if (NULL == tntSysBasisOpGet()) {
68  tntErrorPrint("Cannot create symmetric random wave function as the basis operator for the system has not yet been set",TNT_MPS_ERRCODE); /* NO_COVERAGE */
69  } else { /* NO_COVERAGE */
70  legdim[2] = d = tntNodeGetLegDim(tntSysBasisOpGet(),TNT_MPS_D);
71  }
72 
73  /* Get the quantum numbers for the physical leg */
74  qnums_phys = tntNodeGetQN(tntSysBasisOpGet(), TNT_MPS_D);
75 
76  /* allocate array for minimum and maximum qn for each label on physical leg */
77  qmin = tntIntArrayAlloc(lblnum);
78  qmax = tntIntArrayAlloc(lblnum);
79  dqn = tntIntArrayAlloc(lblnum);
80 
81  /* Get the minimum and maximum quantum number for each term in the label on the physical legs (not assuming that these are in order) */
82  for (i = 0; i < lblnum; i++) { /* loop through the qn's in a label */
83  /* set qmin and qmax equal to first entry initially */
84  qmin.vals[i] = qnums_phys.vals[i];
85  qmax.vals[i] = qnums_phys.vals[i];
86 
87  /* set difference to be zero initially */
88  dqn.vals[i] = 0;
89 
90  for (j = 1; j < d; j++) { /* loop through indices on physical leg */
91  /* index for current qn */
92  qind = i + j*lblnum;
93 
94  /* check if current index has a smaller qn, and reassign if it does */
95  if (qmin.vals[i] > qnums_phys.vals[qind])
96  qmin.vals[i] = qnums_phys.vals[qind];
97 
98  /* check if current index has a larger qn, and reassign if it does */
99  if (qmax.vals[i] < qnums_phys.vals[qind])
100  qmax.vals[i] = qnums_phys.vals[qind];
101 
102  /* if dqn not yet assigned assign difference between this qn and the first qn (NB, may be zero, but makes no difference if it is so no point checking first) */
103  /* note we are assuming that all differences are the same, or at least come in increasing order */
104  if (0 == dqn.vals[i])
105  dqn.vals[i] = abs(qnums_phys.vals[qind]-qnums_phys.vals[i]);
106 
107  } /* end of loop through indices on physical leg */
108 
109  } /* end of loop through qn's in a label */
110 
111  /* check that the total quantum number given actually makes sense */
112  for (i = 0; i < lblnum; i++) { /* loop through the quantum number labels */
113 
114  /* total qn must be more than or equal to (minimum qn on physical leg)*(number of physical legs) */
115  if (qn[i] < numsites*qmin.vals[i]) {
116  if (lblnum > 1) { /* NO_COVERAGE */
117  sprintf(err_msg, "Total quantum number %d (%d of %d qn labels) is invalid|Minimum allowable qn for this label is %d", qn[i], i+1, lblnum, L*qmin.vals[i]); /* NO_COVERAGE */
118  } else { /* NO_COVERAGE */
119  sprintf(err_msg, "Total quantum number %d is invalid - minimum allowable qn is %d", qn[i], L*qmin.vals[i]); /* NO_COVERAGE */
120  } /* NO_COVERAGE */
121  tntErrorPrint(err_msg, TNT_MPS_ERRCODE); /* NO_COVERAGE */
122  } /* NO_COVERAGE */
123 
124  /* total qn must be less than or equal to (maximum qn on physical leg)*(number of physical legs) */
125  if (qn[i] > numsites*qmax.vals[i]) {
126  if (lblnum > 1) { /* NO_COVERAGE */
127  sprintf(err_msg, "Total quantum number %d (%d of %d qn labels) is invalid|Maximum allowable qn for this label is %d", qn[i], i+1, lblnum, L*qmax.vals[i]); /* NO_COVERAGE */
128  } else { /* NO_COVERAGE */
129  sprintf(err_msg, "Total quantum number %d is invalid - maximum allowable qn is %d", qn[i], L*qmin.vals[i]); /* NO_COVERAGE */
130  } /* NO_COVERAGE */
131  tntErrorPrint(err_msg, TNT_MPS_ERRCODE); /* NO_COVERAGE */
132  } /* NO_COVERAGE */
133 
134  /* total qn must be a qn of one of the configurations written in terms of the basis operator.
135  Calculate this by taking the maximum total qn, the difference in qn when there is a lowering operation, and see if an integer number of flips will lead to this qn.
136  Alternatively if dqn[i] = 0 (i.e. all qns are the same for this label) then the total quantum number must be the value on the label * the number of sites */
137  if ((dqn.vals[i] !=0 && (numsites*qmax.vals[i]-qn[i])%dqn.vals[i] != 0) || (0 == dqn.vals[i] && qn[i] != numsites*qmax.vals[i])) {
138  if (lblnum > 1) { /* NO_COVERAGE */
139  sprintf(err_msg, "Total quantum number %d (%d of %d qn labels) is invalid|There are no configurations in a %d site system that have this qn", qn[i], i+1, lblnum, L); /* NO_COVERAGE */
140  } else { /* NO_COVERAGE */
141  sprintf(err_msg, "Total quantum number %d is invalid - there are no configurations in a %d site system that have this qn", qn[i], L); /* NO_COVERAGE */
142  } /* NO_COVERAGE */
143  tntErrorPrint(err_msg, TNT_MPS_ERRCODE); /* NO_COVERAGE */
144  } /* NO_COVERAGE */
145  }
146 
147  /* Now create the empty network */
148  mps = tntNetworkCreate();
149 
150  /* Allocate array for quantum numbers for the left internal leg. Initial value will be zero as required. */
151  qnums_int[0] = tntIntArrayAlloc(lblnum);
152 
153  /* Allocate arrays for the minimum and maximum allowable qn on each site to allow the total to be reached */
154  qminsum = tntIntArrayAlloc(lblnum);
155  qmaxsum = tntIntArrayAlloc(lblnum);
156 
157  /* First left leg will have dimension 1 */
158  legdim[0] = 1;
159 
160  /* Keep inserting random nodes at the end. */
161  for (loop = 0; loop < L; loop++) { /* loop over sites */
162 
163  /* --- Determine the quantum numbers for the right leg ---- */
164 
165  /* Allocate array for the possible quantum numbers for right internal leg. d*D is maximum number of possible labels for the first node */
166  qnums_int[(loop+1)%2] = tntIntArrayAlloc(d*lblnum*legdim[0]);
167 
168  /* Possible quantum numbers are given by sum of those on incoming physical leg and incoming left leg */
169  /* assign d*D possible qn's */
170 
171  for (i = 0; i < legdim[0]; i++) { /* loop through left internal leg indices */
172  for (j = 0; j < d; j++) { /* loop through physical leg indices */
173  for (k = 0; k < lblnum; k++) { /* loop through the qn's in a label */
174 
175  /* current qn is sum of qn label on physical leg and internal leg */
176  qnums_int[(loop+1)%2].vals[(j+d*i)*lblnum+k] = qnums_phys.vals[j*lblnum + k] + qnums_int[loop%2].vals[i*lblnum + k];
177 
178  } /* end of loop through qn's in a label */
179 
180  } /* end of loop through physical leg indices */
181 
182  } /* end of loop through left internal leg indices */
183 
184  /* Determine minimum and maximum allowable qn for the right internal leg to allow the final qn to be reached */
185  for (i = 0; i < lblnum; i++) { /* loop through the qn's in a label */
186 
187  /* quantum number must be more than or equal to qtotal_i - qmax*(number of sites to come) */
188  qminsum.vals[i] = qn[i] - qmax.vals[i]*(L-(loop+1));
189 
190  /* quantum number must be less than or qual to qtotal_i - qmin*(number of sites to come) */
191  qmaxsum.vals[i] = qn[i] - qmin.vals[i]*(L-(loop+1));
192 
193  } /* end of loop through qn's in a label */
194 
195  /* Now go through quantum numbers and remove any that are not allowed, as well as any duplicates */
196  numqn = 0; /* number of qns we are starting with: 0 assigned so far */
197  for (i = 0; i < d*legdim[0]; i++) { /* loop through new qns */
198 
199  /* First check whether the quantum number is within the range of allowed values */
200  for (j = 0; j < lblnum; j++) { /* loop through the qn's in a label */
201  /* reset flag */
202  qn_valid = 1;
203 
204  /* check the qn is valid */
205  if ((qnums_int[(loop+1)%2].vals[i*lblnum+j] < qminsum.vals[j]) || (qnums_int[(loop+1)%2].vals[i*lblnum+j] > qmaxsum.vals[j])) {
206  /* if it isn't, set the flag and break inner loop */
207  qn_valid = 0;
208  break;
209  }
210 
211  } /* end of loop through qn's in a label */
212 
213  /* Now check whether this quantum number has already appeared - loop through all qn currently assigned */
214  for (k = 0; k < numqn; k++) {
215  /* reset flag */
216  qn_appeared = 1;
217 
218  /* loop through every qn in the label, and see if it is the same as the current qn */
219  for (j = 0; j < lblnum; j++) {
220  if (qnums_int[(loop+1)%2].vals[i*lblnum+j] != qnums_int[(loop+1)%2].vals[k*lblnum+j]) {
221  qn_appeared = 0;
222  break;
223  }
224  }
225 
226  /* exit the loop if the qn has been found */
227  if (qn_appeared) {
228  break;
229  }
230 
231  }
232 
233  if (qn_valid&&(!qn_appeared)) {
234  /* if qn is valid and hasn't already appeared, increment the number of unique, allowed qns */
235  numqn++;
236 
237  /* Shuffle down entry if some values have already been removed */
238  if (numqn-1 != i) {
239  for (j = 0; j < lblnum; j++) {
240  qnums_int[(loop+1)%2].vals[(numqn-1)*lblnum+j] = qnums_int[(loop+1)%2].vals[i*lblnum+j];
241  }
242  }
243  }
244  }
245 
246  /* set the right internal leg dimension to be the number of possible qns */
247  legdim[1] = numqn;
248 
249  /* Create the node */
250  A = tntNodeCreate(NULL, numlegs, legid, legdim);
251 
252  /* --- Set qn information for the node --- */
253  /* Note: no warn version of setting QN info used since we know information is being lost */
254 
255  /* Physical leg is incoming leg and has same qn information as operator */
256  tntNodeSetQN_nowarn(A,TNT_MPS_D,&qnums_phys,TNT_QN_IN);
257 
258  /* Left leg is incoming leg, and has values = to right leg of previous tensor */
259  tntNodeSetQN_nowarn(A,TNT_MPS_L,qnums_int+loop%2,TNT_QN_IN);
260 
261  /* Now assign these values to the right outgoing leg */
262  tntNodeSetQN_nowarn(A,TNT_MPS_R,qnums_int+(loop+1)%2,TNT_QN_OUT);
263 
264  /* free left qns */
265  tntIntArrayFree(qnums_int+loop%2);
266 
267  /* Copy the A node only into the network */
269 
270  /* Free the node group consisting of U (A), S and VT */
271  tntNodeGroupFree(&A);
272 
273  /* Next left leg has to have the same dimension as the previous right leg */
274  legdim[0] = legdim[1];
275 
276  } /* end of loop over sites */
277 
278  /* Free the physical and right internal qns */
279  tntIntArrayFree(&qnums_phys);
280  tntIntArrayFree(qnums_int+L%2);
281 
282  /* Free arrays no longer required */
283  tntIntArrayFree(&qmin);
284  tntIntArrayFree(&qmax);
285  tntIntArrayFree(&dqn);
286  tntIntArrayFree(&qminsum);
287  tntIntArrayFree(&qmaxsum);
288 
289  /* orthonormalise the wave function */
290  /* orthogonalise */
291  tntMpsOrth(mps,0);
292  /* find the norm squared */
293  norm_sq = tntMpsSelfProduct(mps,0);
294  /* normalise by scaling the first node */
295  tntNodeScaleReal(tntNodeFindFirst(mps),1.0/sqrt(norm_sq));
296 
297  return mps;
298 
299 }
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:70
Simple 1D array for integer values.
Definition: tnt.h:148
struct tnode * tntNode
Definition: tnt.h:118
unsigned tntNodeGetLegDim(tntNode tn, unsigned legid)
Definition: tntNode_funcs.c:1069
tntNode tntNodeCreate(tntComplexArray *nodeVals, unsigned numLegs, unsigned legId[], unsigned legDim[])
Definition: tntNode_funcs.c:523
int tntSymmTypeGet(void)
Definition: tntMisc_funcs.c:699
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
void tntNodeInsertAtEnd(tntNode tnI, int legIlast, int legInwend, tntNetwork tnw)
Definition: tntNode_funcs.c:945
void tntNodeSetQN_nowarn(tntNode tn, unsigned legid, tntIntArray *qvals, int legdir)
Definition: tntNode_funcs.c:1392
tntIntArray tntNodeGetQN(tntNode tn, unsigned legid)
Definition: tntNode_funcs.c:1442
#define TNT_MPS_L
Definition: tntMps.h:56
unsigned tntSymmNumGet(void)
Definition: tntMisc_funcs.c:713
#define TNT_QN_IN
Definition: tnt.h:71
void tntNodeScaleReal(tntNode tn, double val)
Definition: tntNode_funcs.c:623
void tntNodeGroupFree(tntNode *tn)
Definition: tntNode_funcs.c:1053
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
double tntMpsSelfProduct(tntNetwork mps, int orth_centre)
Definition: tntMpsMps.c:355
struct tnetwork * tntNetwork
Definition: tnt.h:123
tntNetwork tntNetworkCreate(void)
Definition: tntNetwork_funcs.c:62
tntNetwork tntMpsCreateRandom(unsigned L, unsigned D)
Definition: tntMpsCreateRandom.c:34
void tntMpsOrth(tntNetwork mps, unsigned orth_centre)
Definition: tntMpsOrth.c:36
#define TNT_STRLEN
Definition: tnt.h:59
tntNode tntSysBasisOpGet(void)
Definition: tntMisc_funcs.c:647
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
#define TNT_QN_OUT
Definition: tnt.h:76
tntNode tntMpsCreateTwoSiteOp ( tntNodeArray Lnodes,
tntNodeArray Rnodes,
tntComplexArray params 
)

Creates a single two-site operator represented by a tntNode type, using single site operators. The two-site operator is given as the sum of the tensor product of a number \(n_n\) of pairs (given in the arrays Lnodes and Rnodes) of two one-site operators \(\hat{o}^l_i\) and \(\hat{o}^r_i\). Each pair of operators can be scaled by a parameter \(\alpha_i\) (given in the array params). In summary the returned node is given by \( \hat{T} = \sum_i^{n_n}\alpha_{i}\hat{o}^l_i\otimes\hat{o}^r_i \) or:

twositeop_func.png

The operators have four physical legs, each physical leg having the same dimension as the original physical legs, which are labelled as follows:

two_site_gate.png

All nodes will be static nodes.

If there is a global physical symmetry specified, then symmetry information is applied to the legs. The upwards facing legs are set to be outgoing legs, and the downwards facing legs are set to be incoming legs. Note that if the two-site operators are not invariant under the symmetry type chosen, then some information will be discarded. A warning will be printed (by the core level function) if this is the case.

Returns
The node representing the two site gate.
Parameters
LnodesArray of operators for the left legs
RnodesArray of operators for the right legs that should correspond to those for the left operators
paramsParameters for each two-site term
42 {
43 
44  tntNode Tterm, T; /* Nodes for building two-site operator */
45  tntNode tnA, tnB; /* Nodes used while contracting */
46  tntIntArray qnums_phys; /* Quantum numbers for the physical leg */
47  int idmapL[5] = {0,-1,-1,3,1}, idmapR[5] = {0,-1,-1,4,2}; /* Maps used while contracting to form two site operator for left and right nodes */
48  unsigned loop;
49 
50  /* First deal with the nearest neighbour terms */
51  for (loop = 0; loop < Lnodes->sz; loop++) {
52 
53  /* Copy over the nodes and add a leg to each one */
54  /* Left contribution */
55  tnA = tntNodeCopy(Lnodes->vals[loop],0);
56 
57  /* Right contribution */
58  tnB = tntNodeCopy(Rnodes->vals[loop],0);
59 
60  /* Contract them to form a two site operator */
61  Tterm = tntNodeContract(tnA, tnB, idmapL, idmapR);
62 
63  /* Scale by the parameter */
64  tntNodeScaleComplex(Tterm, params->vals[loop]);
65 
66  /* Add to the total Hamiltonian, unless it is the first term */
67  if (0 == loop) {
68  T = Tterm;
69  } else {
70  tntNodeAdd(T, Tterm);
71  /* Free the node for the current term */
72  tntNodeFree(&Tterm);
73  }
74  }
75 
76  /* Now put the unitary in the blocks form */
77  if (tntSymmTypeGet()) {
78  /* Get the quantum numbers for the physical leg */
79  qnums_phys = tntNodeGetQN(tntSysBasisOpGet(), TNT_MPS_D);
80 
81  /* Set the quantum numbers for the unitary - upwards facing legs (1 and 2) are outgoing, downwards facing legs (3 and 4) are incoming. */
82  tntNodeSetQN(T,1,&qnums_phys,TNT_QN_OUT);
83  tntNodeSetQN(T,2,&qnums_phys,TNT_QN_OUT);
84  tntNodeSetQN(T,3,&qnums_phys,TNT_QN_IN);
85  tntNodeSetQN(T,4,&qnums_phys,TNT_QN_IN);
86  /* Note - this may require some elements of the unitary to be discarded - no check at the moment for this */
87 
88  /* free the quantum numbers */
89  tntIntArrayFree(&qnums_phys);
90  }
91 
92  return T;
93 
94 }
Simple 1D array for integer values.
Definition: tnt.h:148
tntComplex * vals
Pointer to the first element in the tntComplex array.
Definition: tnt.h:169
struct tnode * tntNode
Definition: tnt.h:118
void tntNodeScaleComplex(tntNode tn, tntComplex val)
Definition: tntNode_funcs.c:646
int tntSymmTypeGet(void)
Definition: tntMisc_funcs.c:699
tntIntArray tntNodeGetQN(tntNode tn, unsigned legid)
Definition: tntNode_funcs.c:1442
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
#define TNT_QN_IN
Definition: tnt.h:71
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
void tntNodeSetQN(tntNode tn, unsigned legid, tntIntArray *qvals, int legdir)
Definition: tntNode_funcs.c:1324
tntNode * vals
Pointer to the first element in the tntNode array.
Definition: tnt.h:179
void tntNodeFree(tntNode *tn)
Definition: tntNode_funcs.c:420
tntNode tntSysBasisOpGet(void)
Definition: tntMisc_funcs.c:647
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_QN_OUT
Definition: tnt.h:76
void tntNodeAdd(tntNode tnA, tntNode tnB)
Definition: tntNode_funcs.c:385
unsigned sz
The number of elements in the array.
Definition: tnt.h:180
void tntMpsExOpFree ( tntMpsExOp ExOp)

Frees all the dynamically allocated memory associated with the structure for holding expectation value operators.

Returns
No return value.
Parameters
ExOpThe operator structure to free.
Examples:
dmrg.c, and tebd.c.
251 {
252 
253  /* Free all the node arrays */
254  tntNodeArrayFree(&ExOp->ExOpOs);
255  tntNodeArrayFree(&ExOp->ExOp2nn);
256  tntNodeArrayFree(&ExOp->ExOp2cs);
257  tntNodeArrayFree(&ExOp->ExOp2ap);
258 
259  /* Free all the string arrays */
260  tntStringArrayFree(&ExOp->LbOpOs);
261  tntStringArrayFree(&ExOp->LbOp2nn);
262  tntStringArrayFree(&ExOp->LbOp2cs);
263  tntStringArrayFree(&ExOp->LbOp2ap);
264 
265 }
tntStringArray LbOpOs
Array that contains strings for labelling the above operators.
Definition: tntMps.h:110
tntNodeArray ExOp2cs
Array that contains the centre-to-site expectation values to calculate. Will be a tntNodeArray of len...
Definition: tntMps.h:113
tntNodeArray ExOp2ap
Array that contains the all pairs of two-site expectation values to calculate. Will be a tntNodeArray...
Definition: tntMps.h:115
tntStringArray LbOp2cs
Array that contains strings for labelling the above operators.
Definition: tntMps.h:114
void tntStringArrayFree(tntStringArray *arr)
Definition: tntArray_funcs.c:233
void tntNodeArrayFree(tntNodeArray *arr)
Definition: tntArray_funcs.c:210
tntNodeArray ExOp2nn
Array that contains the nearest-neighbour expectation values to calculate.
Definition: tntMps.h:111
tntStringArray LbOp2ap
Array that contains strings for labelling the above operators.
Definition: tntMps.h:116
tntNodeArray ExOpOs
Array that contains the single site expectation values to calculate.
Definition: tntMps.h:109
tntStringArray LbOp2nn
Array that contains strings for labelling the above operators.
Definition: tntMps.h:112
void tntMpsExpecOutput ( tntNetwork  mps,
tntMpsExOp Op,
int  orthCentre,
unsigned  printOutput,
unsigned  saveOutput,
char *  saveprefix,
unsigned  counter 
)

Calculates various types of one-site and two-site expectation values, the operators for which are passed in the tntMpsExOp structure.

The function can output the expectation values to the screen as well as saving them to the named output file. Use the appropriate flags to specify whether results should be output to the screen or to an output file. The function will have no effect (and nothing will be calculated) if both the screen output flag is 0 and the printOutput flag is 0 - in this case a warning will be printed to screen but execution will continue.

Returns
No return value - the expectation values are output to the screen and/or to an output file.
See also
tntMpsProcessExpecOptions() for a function which converts command line arguments to the arrays of operators that are required by this function.
Parameters
mpsThe network representing the MPS. Unchanged by the function.
OpThe operators for calculating expectation values
orthCentreOrthogonality centre of the MPS.
printOutputSet to 1 to print expectation values to screen, 0 otherwise
saveOutputSet to 1 to save values to output file, 0 otherwise
saveprefixname for output file (without extension). Only used if saveOutput is 1.
countercan be used if saving multiple batches of expectation values e.g. for different timesteps. Added as a suffix to the variable name, if a non-zero value is given. Pass zero if a counter is not required.
Examples:
dmrg.c, and tebd.c.
37 {
38 
39  tntDoubleArray Exval; /* array for holding expectation values */
40  tntDoubleArray Exvalm; /* matrix for holding expectation values */
41  double normsq; /* Norm squared of the wave function. Always calculated */
42  unsigned loop; /* Used to loop through the different expectation value types */
43  tntNodeArray op; /* Node array for holding operators for expectation value currently being calculated */
44  tntIntArray sitenum; /* Array for holding site number for expectation value currently being calculated */
45  int L; /* Length of the MPS */
46  char count_label[TNT_STRLEN]; /* If the argument counter is not zero, then this variable will contain the current label appended but the count number */
47 
48 
49  if (!printOutput && !saveOutput) {
50  tntWarningPrint("No expectation values are being calculated as both output flags are set to off"); /* NO_COVERAGE */
51  return; /* NO_COVERAGE */
52  } /* NO_COVERAGE */
53 
54  /* Find the length of the MPS */
55  L = (int) tntMpsLength(mps);
56 
57  /* Calculate the norm squared */
58  normsq = tntMpsSelfProduct(mps,orthCentre);
59 
60  /* ----- Onsite expectation values ----- */
61 
62  if (Op->ExOpOs.sz > 0) {
63 
64  /* Allocate arrays of size one for onsite operators */
65  op = tntNodeArrayAlloc(1);
66  sitenum = tntIntArrayAlloc(1);
67 
68  /* Allocate memory for storing single-site expectation values */
69  Exval = tntDoubleArrayAlloc(L);
70 
71  if (printOutput) {
72  printf("\n~~~~~~~ Onsite expectation values ~~~~~~~~\n");
73  }
74 
75  }
76  /* Find the single site expectation values of the calculated ground state */
77  for (loop = 0; loop < Op->ExOpOs.sz; loop++) {
78 
79  /* Copy over the required operator */
80  op.vals[0] = tntNodeCopy(Op->ExOpOs.vals[loop],0);
81 
82  for (sitenum.vals[0] = 0; sitenum.vals[0] < L; (sitenum.vals[0])++) {
83  /* Find the expectation value on each site using a product MPO. */
84  Exval.vals[sitenum.vals[0]] = tntComplexToReal(tntMpsPmpoMpsProduct(mps, &op, &sitenum, orthCentre))/normsq;
85  }
86  if (printOutput) {
87  tntNamedDoubleArrayPrint(Op->LbOpOs.vals[loop],&Exval);
88  }
89 
90  /* Save this expectation value */
91  if (saveOutput) {
92  if (counter) {
93  sprintf(count_label, "%s_%d", Op->LbOpOs.vals[loop], counter);
94  tntSaveArrays(saveprefix,"", 0, 1, 0, &Exval, count_label);
95  } else {
96  tntSaveArrays(saveprefix,"", 0, 1, 0, &Exval, Op->LbOpOs.vals[loop]);
97  }
98 
99  }
100 
101  /* Free the copy of the operator */
102  tntNodeFree(&(op.vals[0]));
103  }
104  if (Op->ExOpOs.sz > 0) {
105 
106  /* free the arrays */
107  tntDoubleArrayFree(&Exval);
108  tntIntArrayFree(&sitenum);
109  tntNodeArrayFree(&op);
110 
111  if (printOutput) {
112  printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
113  }
114  }
115 
116 
117  /* ---------- Nearest-neighbour expectation values ----------- */
118 
119  if (Op->ExOp2nn.sz > 0) {
120 
121  /* Allocate arrays of size two for two-site operator site numbers */
122  sitenum = tntIntArrayAlloc(2);
123 
124  /* Allocate memory for storing single-site expectation values */
125  Exval = tntDoubleArrayAlloc(L-1);
126 
127  if (printOutput) {
128  printf("\n~~~~~~~ Nearest-neighbour expectation values ~~~~~~~~\n");
129  }
130 
131  }
132  for (loop = 0; loop < Op->ExOp2nn.sz; loop+=2) {
133 
134  /* Allocate arrays of size two for two-site operators */
135  op = tntNodeArrayAlloc(2);
136 
137  /* Copy over the required operator */
138  op.vals[0] = tntNodeCopy(Op->ExOp2nn.vals[loop],0);
139  op.vals[1] = tntNodeCopy(Op->ExOp2nn.vals[loop+1],0);
140 
141  for (sitenum.vals[0] = 0; sitenum.vals[0] < L-1; (sitenum.vals[0])++) {
142  /* assign value for second site */
143  sitenum.vals[1] = sitenum.vals[0]+1;
144  /* Find the expectation value on each site using a product MPO. */
145  Exval.vals[sitenum.vals[0]] = tntComplexToReal(tntMpsPmpoMpsProduct(mps, &op, &sitenum, orthCentre))/normsq;
146  }
147 
148  /* Free the operators */
149  tntNodeArrayFree(&op);
150 
151  /* Print the output */
152  if (printOutput) {
153  tntNamedDoubleArrayPrint(Op->LbOp2nn.vals[loop/2],&Exval);
154  }
155 
156  /* Save this expectation value */
157  if (saveOutput) {
158  if (counter) {
159  sprintf(count_label, "%s_%d", Op->LbOp2nn.vals[loop/2], counter);
160  tntSaveArrays(saveprefix,"", 0, 1, 0, &Exval, count_label);
161  } else {
162  tntSaveArrays(saveprefix,"", 0, 1, 0, &Exval, Op->LbOp2nn.vals[loop/2]);
163  }
164  }
165  }
166  if (Op->ExOp2nn.sz > 0) {
167 
168  /* free the arrays */
169  tntDoubleArrayFree(&Exval);
170  tntIntArrayFree(&sitenum);
171 
172  if (printOutput) {
173  printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
174  }
175  }
176 
177  /* ---------- Centre-site expectation values ----------- */
178 
179  if (Op->ExOp2cs.sz > 0) {
180 
181  /* Allocate arrays of size two for two-site operator site numbers */
182  sitenum = tntIntArrayAlloc(2);
183 
184  /* Allocate memory for storing single-site expectation values */
185  Exval = tntDoubleArrayAlloc(L);
186 
187  if (printOutput) {
188  printf("\n~~~~~~~ Centre-site expectation values ~~~~~~~~\n");
189  }
190 
191 
192 
193  }
194  for (loop = 0; loop < Op->ExOp2cs.sz; loop+=2) {
195 
196  /* Allocate first site that expectation value will be found on */
197  sitenum.vals[0] = L/2;
198 
199  /* Allocate arrays of size two for two-site operators */
200  op = tntNodeArrayAlloc(2);
201 
202  /* Copy over the required operator */
203  op.vals[0] = tntNodeCopy(Op->ExOp2cs.vals[loop],0);
204  op.vals[1] = tntNodeCopy(Op->ExOp2cs.vals[loop+1],0);
205 
206  for (sitenum.vals[1] = 0; sitenum.vals[1] < L; (sitenum.vals[1])++) {
207  /* Find the expectation value on each site using a product MPO. */
208  Exval.vals[sitenum.vals[1]] = tntComplexToReal(tntMpsPmpoMpsProduct(mps, &op, &sitenum, orthCentre))/normsq;
209  }
210 
211  /* Free the operators */
212  tntNodeArrayFree(&op);
213 
214  /* Print the output */
215  if (printOutput) {
216  tntNamedDoubleArrayPrint(Op->LbOp2cs.vals[loop/2],&Exval);
217  }
218 
219  /* Save this expectation value */
220  if (saveOutput) {
221  if (counter) {
222  sprintf(count_label, "%s_%d", Op->LbOp2cs.vals[loop/2], counter);
223  tntSaveArrays(saveprefix,"", 0, 1, 0, &Exval, count_label);
224  } else {
225  tntSaveArrays(saveprefix,"", 0, 1, 0, &Exval, Op->LbOp2cs.vals[loop/2]);
226  }
227  }
228  }
229  if (Op->ExOp2cs.sz > 0) {
230 
231  /* free the arrays */
232  tntDoubleArrayFree(&Exval);
233  tntIntArrayFree(&sitenum);
234 
235  if (printOutput) {
236  printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
237  }
238  }
239 
240 
241  /* ---------- All pairs of two-site expectation values ----------- */
242 
243  if (Op->ExOp2ap.sz > 0) {
244 
245  /* Allocate arrays of size two for sitenumbers for two-site operators */
246  sitenum = tntIntArrayAlloc(2);
247 
248  /* Allocate memory for storing single-site expectation values */
249  Exvalm = tntDoubleArrayAlloc(L,L);
250 
251  if (printOutput) {
252  printf("\n~~~~~~~~~ All pairs of two-site expectation values ~~~~~~~~~\n");
253  }
254 
255  }
256  for (loop = 0; loop < Op->ExOp2ap.sz; loop+=2) {
257  int ind; /* total matrix index */
258 
259  /* Allocate arrays of size two for two-site operators */
260  op = tntNodeArrayAlloc(2);
261 
262  /* Copy over the required operator */
263  op.vals[0] = tntNodeCopy(Op->ExOp2ap.vals[loop],0);
264  op.vals[1] = tntNodeCopy(Op->ExOp2ap.vals[loop+1],0);
265 
266  /* Find the expectation value on each pair of sites using a product MPO. */
267  for (sitenum.vals[0] = 0; sitenum.vals[0] < L; (sitenum.vals[0])++) {
268  for (sitenum.vals[1] = 0; sitenum.vals[1] < L; (sitenum.vals[1])++) {
269  ind = sitenum.vals[0] + L*sitenum.vals[1];
270  Exvalm.vals[ind] = tntComplexToReal(tntMpsPmpoMpsProduct(mps, &op, &sitenum, orthCentre))/normsq;
271  }
272  }
273 
274  /* Free the operators */
275  tntNodeArrayFree(&op);
276 
277  if (printOutput) {
278  tntNamedDoubleArrayPrint(Op->LbOp2ap.vals[loop/2],&Exvalm);
279  }
280 
281  /* Save this expectation value */
282  if (saveOutput) {
283  if (counter) {
284  sprintf(count_label, "%s_%d", Op->LbOp2ap.vals[loop/2], counter);
285  tntSaveArrays(saveprefix,"", 0, 1, 0, &Exvalm, count_label);
286  } else {
287  tntSaveArrays(saveprefix,"", 0, 1, 0, &Exvalm, Op->LbOp2ap.vals[loop/2]);
288  }
289  }
290 
291  }
292  if (Op->ExOp2ap.sz > 0) {
293 
294  /* free the arrays */
295  tntDoubleArrayFree(&Exvalm);
296  tntIntArrayFree(&sitenum);
297 
298  if (printOutput) {
299  printf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
300  }
301 
302  }
303 
304 }
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:70
tntStringArray LbOpOs
Array that contains strings for labelling the above operators.
Definition: tntMps.h:110
Simple 1D array for integer values.
Definition: tnt.h:148
void tntSaveArrays(const char *saveprefix, const char *saveinsert, unsigned num_int, unsigned num_dbl, unsigned num_cmp,...)
Definition: tntSave_funcs.c:133
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
tntNodeArray ExOp2cs
Array that contains the centre-to-site expectation values to calculate. Will be a tntNodeArray of len...
Definition: tntMps.h:113
double tntComplexToReal(tntComplex var)
Definition: tntMisc_funcs.c:1022
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
tntNodeArray ExOp2ap
Array that contains the all pairs of two-site expectation values to calculate. Will be a tntNodeArray...
Definition: tntMps.h:115
tntNodeArray tntNodeArrayAlloc(unsigned num_elems)
Definition: tntArray_funcs.c:110
tntStringArray LbOp2cs
Array that contains strings for labelling the above operators.
Definition: tntMps.h:114
Simple 1D array for double values.
Definition: tnt.h:158
double * vals
Pointer to the first element in the double array.
Definition: tnt.h:159
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
double tntMpsSelfProduct(tntNetwork mps, int orth_centre)
Definition: tntMpsMps.c:355
void tntDoubleArrayFree(tntDoubleArray *arr)
Definition: tntArray_funcs.c:177
Simple 1D array for tntNodes.
Definition: tnt.h:178
void tntNodeArrayFree(tntNodeArray *arr)
Definition: tntArray_funcs.c:210
tntNodeArray ExOp2nn
Array that contains the nearest-neighbour expectation values to calculate.
Definition: tntMps.h:111
tntDoubleArray tntDoubleArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:83
tntStringArray LbOp2ap
Array that contains strings for labelling the above operators.
Definition: tntMps.h:116
char ** vals
Pointer to the first element in the array of strings.
Definition: tnt.h:187
tntNode * vals
Pointer to the first element in the tntNode array.
Definition: tnt.h:179
tntNodeArray ExOpOs
Array that contains the single site expectation values to calculate.
Definition: tntMps.h:109
void tntNodeFree(tntNode *tn)
Definition: tntNode_funcs.c:420
#define TNT_STRLEN
Definition: tnt.h:59
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntStringArray LbOp2nn
Array that contains strings for labelling the above operators.
Definition: tntMps.h:112
tntComplex tntMpsPmpoMpsProduct(tntNetwork mps, tntNodeArray *op, tntIntArray *sitenum, int orth_centre)
Definition: tntMpsPmpoMps.c:56
unsigned sz
The number of elements in the array.
Definition: tnt.h:180
unsigned tntMpsLength ( tntNetwork  mps)

Counts the number of nodes in the MPS by moving along the internal legs. The result is not altered by any nodes connected to the physical legs (i.e. legs labelled TNT_MPS_D).

Returns
The number of nodes in the network.
Parameters
mpsThe network representing the MPS. Unchanged by the function.
209 {
210 
211  tntNode A; /* The current MPS node */
212  tntNode A_last; /* The last MPS node */
213  unsigned numsites = 0; /* The number of MPS nodes in the network */
214 
215  /* Loop through network determining the number of sites in the mps */
216  A = tntNodeFindFirst(mps);
217  A_last = tntNodeFindLast(mps);
218 
219  /* The loop will continue until the last site is reached */
220  while (1) {
221  /* Increment the number of sites in the network */
222  numsites++;
223 
224  /* If the last site has been reached, exit the loop, otherwise move to the next site */
225  if (A == A_last) {
226  break;
227  } else {
228  A = tntNodeFindConn(A, TNT_MPS_R);
229  }
230 
231  }
232 
233  /* Return the number of sites in the MPS network */
234  return numsites;
235 
236 }
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
struct tnode * tntNode
Definition: tnt.h:118
tntNode tntNodeFindLast(tntNetwork tnw)
Definition: tntNode_funcs.c:685
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
#define TNT_MPS_R
Definition: tntMps.h:61
void tntMpsMpoConnect ( tntNetwork  mps,
tntNetwork  mpo 
)

Creates a network that consists of an MPS connected to an MPO.

mps_mpo.png

As illustrated above, it removes the connection between starting and ending singleton legs of the MPO and the original network start and end after connecting them.

If the MPS represents a wave function \(|\psi\rangle\) and the MPO represents a site-wide operator \(\hat{O}\), the returned network represents \(\hat{O}|\psi\rangle \).

A copy of the MPO network is used, and it is connected to the original MPS network. Make a copy of your MPS network first if you require it to be unchanged.

Returns
No return value.
Parameters
mpsNetwork representing the MPS
mpoNetwork representing the MPO
32 {
33 
34  tntNode tnA; /* the current node A */
35  tntNode tno, tnoc, tno_prev = NULL; /* The node representing the operator tno on site j it's copy and j-1 of the current term in the MPO */
36  unsigned L, squeezeleg; /* Length of system, leg of MPO node to remove */
37  char err_msg[TNT_STRLEN]; /* buffer to hold error messages */
38 
39  /* First check that the MPS and the MPO have the same number of terms */
40  L = tntMpsLength(mps);
41 
42  if (tntMpsLength(mpo) != L) {
43  sprintf(err_msg, "Cannot connect the MPS to the MPO as there are not the same number of nodes in the two networks|" /* NO_COVERAGE */
44  "MPS contains %d nodes, MPO contains %d nodes",L,tntMpsLength(mpo)); /* NO_COVERAGE */
45  tntErrorPrint(err_msg, TNT_MPS_ERRCODE); /* NO_COVERAGE */
46  } /* NO_COVERAGE */
47 
48  /* Find the MPS node corresponding to site 0 */
49  tnA = tntNodeFindFirst(mps);
50 
51  /* Find the MPO node for site 0 */
52  tno = tntNodeFindFirst(mpo);
53 
54  /* Now loop through all sites, attaching a copy of the MPO term to the MPS term */
55  while (1) {
56 
57  /* Make a copy of the relevant term */
58  tnoc = tntNodeCopy(tno,0);
59 
60  /* Join a copy of the MPO to the MPS */
62 
63  /* Join to the previous MPO term, if it exists */
64  if (NULL != tno_prev) {
65  tntNodeJoin(tno_prev, TNT_MPS_R, tnoc, TNT_MPS_L);
66  }
67 
68  /* Find the next group of sites if they exist, or exit the loop if the last site has been reached */
69  if (tnA == tntNodeFindLast(mps)) {
70  break;
71  } else {
72  /* Find the next node along in the unflipped wave function */
73  tnA = tntNodeFindConn(tnA, TNT_MPS_R);
74 
75  /* Find the next node along in the MPO */
76  tno = tntNodeFindConn(tno, TNT_MPS_R);
77 
78  /* Set the current MPO term in the new network to the previous term */
79  tno_prev = tnoc;
80  }
81  }
82 
83  /* Find the first MPO term attached to the MPS, and remove the singleton starting leg */
84  /* first MPS node */
85  tnA = tntNodeFindFirst(mps);
86  /* MPO attached to first MPS node */
87  tno = tntNodeFindConn(tnA, TNT_MPS_D);
88  /* Leg to remove */
89  squeezeleg = TNT_MPS_L;
90  /* remove the leg */
91  tntNodeSqueeze(tno, 1, &squeezeleg);
92 
93  /* Find the first MPO term attached to the MPS, and remove the singleton starting leg */
94  /* first MPS node */
95  tnA = tntNodeFindLast(mps);
96  /* MPO attached to first MPS node */
97  tno = tntNodeFindConn(tnA, TNT_MPS_D);
98  /* Leg to remove */
99  squeezeleg = TNT_MPS_R;
100  /* remove the leg */
101  tntNodeSqueeze(tno, 1, &squeezeleg);
102 
103 }
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
struct tnode * tntNode
Definition: tnt.h:118
tntNode tntNodeFindLast(tntNetwork tnw)
Definition: tntNode_funcs.c:685
#define TNT_MPS_L
Definition: tntMps.h:56
#define TNT_MPS_U
Definition: tntMps.h:71
void tntNodeSqueeze(tntNode tn, unsigned squeezeNum, unsigned *squeezeLegs)
Definition: tntNode_funcs.c:1033
#define TNT_STRLEN
Definition: tnt.h:59
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
void tntNodeJoin(tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:462
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
double tntMpsMpoContract ( tntNetwork  mpsmpo,
int  chi 
)

Contracts an MPS-MPO network e.g. one generated using tntMpsMpoConnect() and at the same time zips through the network returning it to have a state with an orthogonality centre on the first site. The function also truncates during the zip procedure, setting the internal dimension to the parameter chi. If chi is set to -1 then no truncation is performed, other than discarding singular values that are smaller than the truncation tolerance, set through core library function tntSVDTruncTolSet().

The network passed as an argument must have the following form:

mps_mpo.png

First a left to right sweep is performed, contracting the MPO node with it's corresponding MPS node, and truncating to the truncation tolerance for singular values (use tntSVDTruncTolSet() to set this) rather than to chi. The network will now be in the MPS form. The routine then sweeps right to left, truncating to chi, using tntMpsTruncate().

Returns
Returns the sum of the truncation errors of all the SVDs.
Parameters
mpsmpoNetwork representing the MPS-MPO. Converted to an MPS network by the function.
chiMaximum internal dimension of the new MPS
124 {
125 
126  tntNode tnA, tnW; /* the current node A of the MPS and the current node W of the MPO */
127  tntNode tnC, tnCc; /* The node that is the result of contracting the MPS and its MPO node, and a copy of it taken to SVD */
128  tntNode tnI; /* The node that results from contracting the next MPS node along and it's MPO node, before contracting with tnC to the left to form a new C node */
129  tntNode tnU, tnUdag; /* First unitary that is a result of SVDing tnCc, and its conjugate */
130  unsigned L, j; /* the length of the MPS network, site counter */
131  double err_tot=0.0, err; /* total truncation error, and truncation error per SVD. */
132 
133  /* Maps */
134  int map_W[4] = {0,4,5,3}; /* Map used for tnW when contracting tnA to form tnC or tnI (legs of A keep their labels) */
135  int startlegs[2] = {1,3}; /* Legs assigned to U in the SVD */
136  unsigned map_Udag[4] = {0,2,1,3}; /* Map used to turn complex conjuate of U to hermitian conjugate of U following convention for direction of leg labels in MPS */
137 
138  /* Find the length of the network */
139  L = tntMpsLength(mpsmpo);
140 
141  /* ------ contracting from the left to the right, truncating to truncation tolerace, to result in a single MPS ---- */
142 
143  /* Find the MPS node corresponding to site 0 */
144  tnA = tntNodeFindFirst(mpsmpo);
145 
146  /* Find the MPO node for the current site */
147  tnW = tntNodeFindConn(tnA, TNT_MPS_D);
148 
149  /* Contract these two nodes to form tnC */
150  tnC = tntNodeContract(tnA, tnW, NULL, map_W);
151 
152  for (j = 0; j < L-1; j++) {
153 
154  /* Take a copy of tnC - this copy will not be connected to any other nodes in the network */
155  tnCc = tntNodeCopy(tnC,0);
156 
157  /* Take an SVD of tnCc, assigning the result back to the copy. Use -1 for chi to truncate to truncation tolerance only. */
158  tnCc = tntNodeSVD(tnCc, startlegs, 2, TNT_MPS_R, TNT_MPS_L, NULL, -1, &err);
159  err_tot += err;
160 
161  /* Make an unconnected copy of the first unitary - this will become the new MPS node to replace tnA when inserted back in the network */
162  tnU = tntNodeCopy(tnCc,0);
163 
164  /* the rest of the nodes are not needed -> discard them */
165  tntNodeGroupFree(&tnCc);
166 
167  /* Legs 1 and 3 of tnU define one dimension (lets call rows) of the Unitary matrix, leg 2 of tnU defines the other dimension (lets call columns) of the unitary matrix */
168 
169  /* Take a conjugate copy of tnU */
170  tnUdag = tntNodeCopy(tnU,1);
171 
172  /* Now map the legs of tnUdag so that leg2 -> leg 1 and leg 1 defines rows, and leg 1 -> leg 2 and leg 2 and 3 define columns.
173  * tnUdag is now hermitian conjuate of tnU, and when leg 2 of tnU contracts with leg 1 of tnUdag the result will be the identity */
174  tntNodeMapLegs(tnUdag, map_Udag);
175 
176  /* Now insert --tnU--tnUdag-- = the identity (when the physical legs are fused with the outer legs), back into the MPS network */
177  /* Take a different action depending on whether we are at the start of the network, or in the middle of the network */
178  if (0 == j) {
179  /* if we inserting are at the beginning of the network, we can put insert the nodes directly, without temporariyily (and illegally) connecting legs with non-matching dimensions */
180 
181  /* Insert the node at the start, which connects leg 2 of tnUdag to leg 1 of tnW */
182  tntNodeInsertAtStart(tnUdag, TNT_MPS_L, TNT_MPS_R, mpsmpo);
183 
184  /* Note that we defined legs 1 and 3 of tnC as rows, and legs 2 and 3 of tnUdag as columns - for a full insertion we also need to contract their physical legs */
185  tntNodeJoin(tnUdag, TNT_MPS_D, tnC, TNT_MPS_D);
186 
187  /* Now insert tnU at the start of the network. tnU becomes the new MPS node for this set (the previous MPS node is part of tnC) */
189  } else {
190  /* if we are inserting in the middle of the network, we can't simply insert the nodes in, as this would cause a mismatch between dimensions of connecting legs */
191  /* This operation is thus forbidden and would cause an error in the program */
192 
193  /* Instead, split the connections between tnC and the previous MPS node, and then join up nodes one by one */
194  /* first identify the previous node */
195  tnA = tntNodeFindConn(tnC, TNT_MPS_L);
196 
197  /* Split connections between tnA and tnC */
198  tntNodeSplit(tnA, tnC);
199 
200  /* connect tnA to tnU */
201  tntNodeJoin(tnA, TNT_MPS_R, tnU, TNT_MPS_L);
202 
203  /* connect tnU to tnUdag */
204  tntNodeJoin(tnU, TNT_MPS_R, tnUdag, TNT_MPS_L);
205 
206  /* Connect tnUdag to tnC, remembering they are joined on 2 legs */
207  tntNodeJoin(tnUdag, TNT_MPS_R, tnC, TNT_MPS_L);
208  tntNodeJoin(tnUdag, TNT_MPS_D, tnC, TNT_MPS_D);
209 
210  /* the network is now complete again */
211  }
212 
213  /* Now contract tnUdag and tnC, assigning the result back to tnC */
214  tnC = tntNodeContract(tnUdag,tnC,NULL,NULL);
215 
216  /* find the next MPS node along, and it's corresponding MPO node */
217  tnA = tntNodeFindConn(tnC, TNT_MPS_R);
218 
219  /* Find the MPO node for the current site */
220  tnW = tntNodeFindConn(tnA, TNT_MPS_D);
221 
222  /* Contract these two nodes to form tnI */
223  tnI = tntNodeContract(tnA, tnW, NULL, map_W);
224 
225  /* Contract these with tnC */
226  tnC = tntNodeContract(tnI,tnC,NULL,NULL);
227 
228  /* If tnA was the last node, tnC simply forms the final MPS node.
229  * Otherwise continue the loop to form further MPS nodes */
230  }
231 
232  /* Now sweep back right to left, this time truncating to the maximum internal dimension chi */
233  /* Note setting desired orthogonality centre to zero enures a full right to left sweep */
234  err = tntMpsTruncate(mpsmpo,0,chi);
235  err_tot += err;
236 
237  /* return the total truncation error */
238  return err_tot;
239 
240 }
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
struct tnode * tntNode
Definition: tnt.h:118
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
#define TNT_MPS_L
Definition: tntMps.h:56
tntNode tntNodeSVD(tntNode tn, int *startLegs, unsigned startNum, int legUS, int legSU, int *legMap, int connDim, double *err)
Definition: tntNode_funcs.c:64
void tntNodeGroupFree(tntNode *tn)
Definition: tntNode_funcs.c:1053
void tntNodeInsertAtStart(tntNode tnI, int legInwstart, int legIfirst, tntNetwork tnw)
Definition: tntNode_funcs.c:905
double tntMpsTruncate(tntNetwork mps, unsigned orth_centre, int chi)
Definition: tntMpsTruncate.c:44
void tntNodeSplit(tntNode tnA, tntNode tnB)
Definition: tntNode_funcs.c:747
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
void tntNodeJoin(tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:462
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
void tntNodeMapLegs(tntNode tn, unsigned legmap[])
Definition: tntNode_funcs.c:705
tntNetwork tntMpsMpoMpsConnect ( tntNetwork  mps,
tntNetwork  mpo 
)

Returns a network that consists of an MPS and its flipped counterpart with an MPO inserted between them i.e.

mps_mpo_mps.png

If the MPS represents a wave function \(|\psi\rangle\) and the MPO represents a site-wide operator \(\hat{O}\), the returned network represents \(\langle\psi|\hat{O}|\psi\rangle \).

As illustrated above, it removes the starting and ending singleton legs of the MPO after connecting them.

Copies of the original networks are used, and the original networks are unchanged.

Returns
The MPS-MPO-MPS network illustrated above.
Parameters
mpsNetwork representing the MPS
mpoNetwork representing the MPO
32 {
33 
34  tntNetwork sandwich; /* The network representing MPS-MPO-MPS */
35  tntNode tnA, tnAf; /* the current node A of the MPS and the flipped (hermitian conjugate) of the current node Af */
36  tntNode tno, tnoc, tno_prev = NULL; /* The node representing the operator tnh on site j it's copy and j-1 of the current term in the MPO */
37  unsigned L; /* Length of system */
38  char err_msg[TNT_STRLEN]; /* buffer to hold error messages */
39 
40  /* First check that the MPS and the MPO have the same number of terms */
41  L = tntMpsLength(mps);
42 
43  if (tntMpsLength(mpo) != L) {
44  sprintf(err_msg, "Cannot connect the MPS to the MPO as there are not the same number of nodes in the two networks|" /* NO_COVERAGE */
45  "MPS contains %d nodes, MPO contains %d nodes",L,tntMpsLength(mpo)); /* NO_COVERAGE */
46  tntErrorPrint(err_msg, TNT_MPS_ERRCODE); /* NO_COVERAGE */
47  } /* NO_COVERAGE */
48 
49  sandwich = tntNetworkCopy(mps); /* Set the network to be a copy of the original wave function */
50 
51  tntMpsMpsConnect(sandwich,sandwich); /* Connect the network to the complex conjugate of itself */
52 
53  /* Find the MPS node corresponding to site 0 */
54  tnA = tntNodeFindFirst(sandwich);
55 
56  /* Find the MPO node for site 0 */
57  tno = tntNodeFindFirst(mpo);
58 
59  /* Find flipped MPS node corresponding to site 0 */
60  tnAf = tntNodeFindConn(tnA, TNT_MPS_D);
61 
62  /* Now loop through all sites, insterting a copy of the MPO term between the wave function and its complex conjugate */
63  while (1) {
64 
65  /* Make a copy of the relevant term */
66  tnoc = tntNodeCopy(tno,0);
67 
68  /* Insert the copy of the current term of the MPO in the network */
70 
71  /* Join to the previous MPO term, if it exists */
72  if (NULL != tno_prev) {
73  tntNodeJoin(tno_prev, TNT_MPS_R, tnoc, TNT_MPS_L);
74  }
75 
76  /* Find the next group of sites if they exist, or exit the loop if the last site has been reached */
77  if (tnA == tntNodeFindLast(sandwich)) {
78  break;
79  } else {
80  /* Find the next node along in the unflipped wave function */
81  tnA = tntNodeFindConn(tnA, TNT_MPS_R);
82 
83  /* identify the flipped node tnf */
84  tnAf = tntNodeFindConn(tnA, TNT_MPS_D);
85 
86  /* Find the next node along in the MPO */
87  tno = tntNodeFindConn(tno, TNT_MPS_R);
88 
89  /* Set the current MPO term in the new network to the previous term */
90  tno_prev = tnoc;
91  }
92  }
93 
94  return sandwich;
95 }
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
struct tnode * tntNode
Definition: tnt.h:118
tntNode tntNodeFindLast(tntNetwork tnw)
Definition: tntNode_funcs.c:685
void tntMpsMpsConnect(tntNetwork mpsA, tntNetwork mpsB)
Definition: tntMpsMps.c:38
void tntNodeInsert(tntNode tnI, int legIA, int legIB, tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:879
#define TNT_MPS_L
Definition: tntMps.h:56
struct tnetwork * tntNetwork
Definition: tnt.h:123
#define TNT_MPS_U
Definition: tntMps.h:71
#define TNT_STRLEN
Definition: tnt.h:59
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
tntNetwork tntNetworkCopy(tntNetwork tnw)
Definition: tntNetwork_funcs.c:36
void tntNodeJoin(tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:462
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
tntComplex tntMpsMpoMpsContract ( tntNetwork  sandwich)

Contracts an MPS-MPO-MPS network sandwich e.g. one generated using tntMpsMpoMpsConnect() to return a scalar complex value.

The sandwich network passed as an argument must have the following form:

mps_mpo_mps.png

There is no requirement for the MPS networks forming the top and bottom of the network to be the same, or any requirement on the dimensions in the network. The network is destroyed by the function, so a copy should be made first if it will be required later.

Returns
A complex scalar value that is the result of the network contraction.
Parameters
sandwichNetwork representing the MPS-MPO-MPS sandwich. Destroyed by the function.
111 {
112  tntNode tnA, tnAf; /* the current node A of the MPS and the flipped (hermitian conjugate) of the current node Af */
113  tntNode tno; /* The node representing the operator tno on site j of the current term in the MPO */
114  tntNode beta; /* The node that is the results of all previously contracted nodes from the left */
115  tntComplex contracted; /* The scalar value of the contracted network */
116  unsigned L, j; /* the length of the MPS network, site counter */
117 
118  /* Maps for the leg id's on contraction */
119  int idmap_A[4] = {0,1,2,7}; /* Map used when contracting A with Heff */
120  int idmap_op[4] = {0,3,4,7}; /* Used for the operator when contracting with A or Heff */
121  int idmap_Afl[3] = {0,5,6}; /* Used for the flipped tensor contracting with the result of the A-operator contraction or Heff */
122 
123  /* Find the length of the network */
124  L = tntMpsLength(sandwich);
125 
126  /* ------ contracting from the left to the right to result in one node remaining (which should have size 1) ---- */
127 
128  /* Find the MPS node corresponding to site 0 */
129  tnA = tntNodeFindFirst(sandwich);
130 
131  /* Find the MPO node for site 0 */
132  tno = tntNodeFindConn(tnA, TNT_MPS_D);
133 
134  /* Find flipped MPS node corresponding to site 0 */
135  tnAf = tntNodeFindConn(tno, TNT_MPS_D);
136 
137  /* Contract these three nodes to find beta */
138  beta = tntNodeContract(tnA, tno, NULL,idmap_op);
139  beta = tntNodeContract(beta, tnAf, NULL, idmap_Afl);
140 
141 
142  /* Now loop through the network from left to right performing contractions */
143  for (j = 1; j < L; j++) {
144 
145  /* Find the next node along in the unflipped wave function */
146  tnA = tntNodeFindConn(beta, TNT_MPS_R);
147 
148  /* identify the operator for this term in the Hamiltonian tno */
149  tno = tntNodeFindConn(tnA, TNT_MPS_D);
150 
151  /* identify the flipped node tnf */
152  tnAf = tntNodeFindConn(tno, TNT_MPS_D);
153 
154  /* Contract beta with the MPS node */
155  beta = tntNodeContract(beta,tnA,NULL,idmap_A);
156 
157  /* Contract beta with the MPO term */
158  beta = tntNodeContract(beta,tno,NULL,idmap_op);
159 
160  /* Contract beta with the flipped MPS node (final contraction step) */
161  beta = tntNodeContract(beta,tnAf,NULL,idmap_Afl);
162 
163  }
164 
165  /* Get the value in the tensor */
166  contracted = tntNodeGetFirstEl(beta);
167 
168  /* Free the network */
169  tntNetworkFree(&sandwich);
170 
171  /* return the expectation value */
172  return contracted;
173 
174 }
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
struct tnode * tntNode
Definition: tnt.h:118
void tntNetworkFree(tntNetwork *tnwp)
Definition: tntNetwork_funcs.c:84
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
tntComplex tntNodeGetFirstEl(tntNode tn)
Definition: tntNode_funcs.c:770
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
Use this type to define complex numbers. Any complex arithmetic needs to be done manually. For example to add two complex numbers A and B to give C write:
Definition: tnt.h:136
tntComplex tntMpsMpoMpsProduct ( tntNetwork  mps,
tntNetwork  mpo 
)

Constructs an MPS-MPO-MPS a network that consists of an MPS and its flipped counterpart with an MPO inserted between them i.e.

mps_mpo_mps.png

and contracts the whole network to find the scalar value of this product. If the MPS represents a wave function \(|\psi\rangle\) and the MPO represents a site-wide operator \(\hat{O}\), the result of the contraction is \(\langle\psi|\hat{O}|\psi\rangle \). This result gives the expectation value of the operator \(\hat{O}\) if the wave function \(|\psi\rangle\) is normalised. If the wave function is not normalised, the use tntMpsMpsContract to find the norm squared (i.e. \(\langle\psi|\psi\rangle\)) and divide the result by that.

If the result of the contraction is real, the value itself is returned. If the result of the contraction is complex, then the absolute value is returned.

Copies of the original networks are used, and the original networks are unchanged.

Returns
The result of contracting the MPS-MPO-MPS network illustrated above.
Parameters
mpsNetwork representing the MPS. Unchanged by the function.
mpoNetwork representing the MPO. Unchanged by the function.
Examples:
dmrg.c.
197 {
198  tntNetwork sandwich; /* The full MPS-MPO-MPS network */
199 
200  /* Create the full network */
201  sandwich = tntMpsMpoMpsConnect(mps,mpo);
202 
203  /* Contract the full network */
204  return tntMpsMpoMpsContract(sandwich);
205 }
struct tnetwork * tntNetwork
Definition: tnt.h:123
tntNetwork tntMpsMpoMpsConnect(tntNetwork mps, tntNetwork mpo)
Definition: tntMpsMpoMps.c:30
tntComplex tntMpsMpoMpsContract(tntNetwork sandwich)
Definition: tntMpsMpoMps.c:110
double tntMpsMpoProduct ( tntNetwork  mps,
tntNetwork  mpo,
int  chi 
)

Constructs an MPS-MPO a network i.e.

mps_mpo.png

and contracts the whole network to return back to an MPS form with a maximum internal dimension \(\chi\) set by the input parameter chi. If the MPS represents a wave function \(|\psi\rangle\) and the MPO represents a site-wide operator \(\hat{O}\), the result of the contraction is \(\hat{O}|\psi\rangle \). This result gives the expectation value of the operator \(\hat{O}\) if the wave function \(|\psi\rangle\) is normalised.

If there is no truncation, this operation would result in the internal dimension of the mps increasing from \(D\) to \(kD\), where \(k\) is the internal dimension of the MPO. Setting chi to -1 will apply the MPO with no truncation other than truncating to the global truncation tolerance set by tntSVDTruncTolSet(). Setting chi to a positive value will apply further truncation to the value given.

First a left to right sweep is performed, contracting the MPO node with it's corresponding MPS node, and truncating to the truncation tolerance for singular values rather than to chi. The network will now be in the MPS form. The routine then sweeps right to left, truncating to chi, using tntMpsTruncate(). After the routine, the orthogonality centre of the MPS will be on the first site.

A copy of the MPO is used in the contraction, so the MPO passed as an argument is unchanged by the operation. The MPS passed as an argument is changed by the contraction by having the operator defined by the MPO applied to it, so make a copy of it first if the original MPS will be needed again.

Returns
The truncation error that results from truncating maximum internal dimension of the MPS to chi.
Parameters
mpsNetwork representing the MPS. Values and internal dimension changed by the function.
mpoNetwork representing the MPO. Unchanged by the function.
chiMaximum internal dimension of the new MPS
269 {
270 
271  /* Create the full network */
272  tntMpsMpoConnect(mps,mpo);
273 
274  /* Contract the full network and return the error */
275  return tntMpsMpoContract(mps,chi);
276 
277 }
void tntMpsMpoConnect(tntNetwork mps, tntNetwork mpo)
Definition: tntMpsMpo.c:30
double tntMpsMpoContract(tntNetwork mpsmpo, int chi)
Definition: tntMpsMpo.c:122
void tntMpsMpsConnect ( tntNetwork  mpsA,
tntNetwork  mpsB 
)

Flips an MPS mpsB and joins it to the original MPS mpsA along the physical legs. The flipping operation involves:

  • renumbering the downwards facing physical leg as an upwards facing physical leg
  • taking the complex conjugate of the tensor elements in each node
  • if symmetry information is defined, turning the physical leg from an incoming leg to an outgoing leg.

If the original MPSs mpsA, mpsB represent a wave functions \(|\psi_\mathrm{A}\rangle\), \(|\psi_\mathrm{B}\rangle\) then the new MPS network is \(\langle\psi_\mathrm{B}|\psi_\mathrm{A}\rangle\).

If any of the physical legs are already connected to something, this function will result in an error.

The function copies, node by node, each MPS node in mpsB and connects the flipped version to the network mpsA. Once the function is complete, the network mpsA will look like the following:

mpsA_mpsB.png

The network mpsB will be unchanged as a copy is used. Note that passing the same network for both arguments is permitted, although in this case mpsB is changed in the same way as mpsA since they are the same network.

Returns
No return value.
Parameters
mpsAThe MPS network to join the flipped network to. On entry: An MPS network with all physical legs unconnected. On exit: An MPS network with all physical legs connected to the flipped MPS network for mpsB.
mpsBThe MPS network to flip and connect to mpsA. Unchanged by the function.
42 {
43 
44  tntNode A_curr, A_nxt; /* The current node and the next node of the unflipped MPS */
45  tntNode B_curr, B_nxt; /* The current node and the next node of the MPS to flip */
46  tntNode Bf, Bf_nxt; /* Flipped copy of the current node and the next node in mpsB */
47  unsigned flipmap[4] = {0,1,2,4}; /* Map to change upwards facing physical leg to downwards facing physical leg */
48 
49  /* Before performing the flip operation, check that the MPS networks have the same number of sites */
50  if (tntMpsLength(mpsA) != tntMpsLength(mpsB)) {
51  tntErrorPrint("The MPS networks cannot be connected as they do not have the same number of sites",TNT_MPS_ERRCODE); /* NO_COVERAGE */
52  } /* NO_COVERAGE */
53 
54  /* Find the first node in wfa */
55  A_curr = tntNodeFindFirst(mpsA);
56 
57  /* Find the first node in wfb */
58  B_curr = tntNodeFindFirst(mpsB);
59 
60  /* Make a new flipped node
61  - First copy the node in wfb, specifying that the complex conjugate should be taken.
62  - Then change the id of a downwards facing physical leg to an upwards facing physical leg
63  */
64  Bf = tntNodeCopy(B_curr, 1);
65  tntNodeMapLegs(Bf, flipmap);
66 
67  /* Join the flipped and unflipped node on the physical leg */
68  tntNodeJoin(A_curr,TNT_MPS_D,Bf,TNT_MPS_U);
69 
70 
71  /* Move through the network adding a flipped node to each physical leg of the original wave function */
72  while (A_curr != tntNodeFindLast(mpsA)) {
73 
74  /* identify the next node along in wfa and wfb */
75  A_nxt = tntNodeFindConn(A_curr, TNT_MPS_R);
76  B_nxt = tntNodeFindConn(B_curr, TNT_MPS_R);
77 
78  /* Make a new flipped node */
79  Bf_nxt = tntNodeCopy(B_nxt, 1);
80  tntNodeMapLegs(Bf_nxt, flipmap);
81 
82  /* Join the original flipped node and the new flipped node */
83  tntNodeJoin(Bf, TNT_MPS_R, Bf_nxt, TNT_MPS_L);
84 
85  /* Join the new flipped node with its unflipped node on their physical legs */
86  tntNodeJoin(A_nxt,TNT_MPS_D,Bf_nxt,TNT_MPS_U);
87 
88  /* Move along one site */
89  A_curr = A_nxt;
90  B_curr = B_nxt;
91  Bf = Bf_nxt;
92  }
93 
94  return;
95 
96 }
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
struct tnode * tntNode
Definition: tnt.h:118
tntNode tntNodeFindLast(tntNetwork tnw)
Definition: tntNode_funcs.c:685
#define TNT_MPS_L
Definition: tntMps.h:56
#define TNT_MPS_U
Definition: tntMps.h:71
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
void tntNodeJoin(tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:462
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
void tntNodeMapLegs(tntNode tn, unsigned legmap[])
Definition: tntNode_funcs.c:705
tntComplex tntMpsMpsContract ( tntNetwork  sandwich,
int  site_left,
int  site_right 
)

Calculates the inner product of two MPSs joined to form a network as follows:

mpsA_mpsB.png

Note the MPSs forming the top and bottom parts of the sandwich do not have to be the same. If the original MPSs mpsA, mpsB represent a wave functions \(|\psi_\mathrm{A}\rangle\), \(|\psi_\mathrm{B}\rangle\) then this represents the inner product \(|\langle\psi_\mathrm{B}|\psi_\mathrm{A}\rangle|^2\).

This network can be created using function tntMpsMpsConnect(), where the MPSs forming the top or bottom of the network can be the same MPS of different MPSs.

The function can be made more efficient by make use of orthogonality centres of the MPS. All sites to the right of site_right should obey the right orthonormalisation condition:

right_orth.png

All sites to the left of site_left should obey the left orthonormalisation condition:

left_orth.png

This can then reduce the number of contractions required to find the expectation value. If either site_left or site_right is -1, then a full contraction is performed.

Note that after a left to right TEBD or DMRG sweep, the orthogonality centre will be the last site, and after a right to left sweep, the orthonormality centre will be the first site.

Returns
The complex value of the inner product of the MPS-MPS network, equivalent to \(|\langle\psi_\mathrm{B}|\psi_\mathrm{A}\rangle|^2\) if the MPSs represent wave functions.
Parameters
sandwichThe network representing the MPS-MPS sandwich. Destroyed by the function.
site_leftLeft-most site to keep. Any sites to the left of this are assumed to obey the left orthonormalisation condition. Use -1 for full contraction.
site_rightRight-most site to keep. Any sites to the left of this are assumed to obey the left orthonormalisation condition. Use -1 for full contraction.
133 {
134 
135  /* Maps for the leg id's on contraction */
136  int idmap_mps[8] = {0,1,2,5,4,8,6,7}; /* Used for the original mps when being multiplied by alpha */
137  int idmap_fl[4] = {0,3,4,-1}; /* Used for the flipped tensor when being multiplied by alpha or the original tensor*/
138 
139  /* Networks and nodes during the contraction */
140  tntNetwork mps_disc; /* The part of the network that will be equated to the identity and discarded during the contraction */
141  tntNode A, Af; /* The current MPS node, and flipped counterpart */
142  tntNode nodes_L[2], nodes_R[2]; /* Arrays for the nodes on the left and right side of the split to split the network in two. */
143  tntNode alpha; /* The large node that is a result of contracting the MPS */
144 
145  /* Counters and site information used during the contraction */
146  int sitecount = 0; /* the current site count */
147  int numsites; /* The number of MPS nodes in the network */
148  unsigned leg_s = TNT_MPS_L; /* Leg that points to the start of the reformed network */
149 
150  /* Final results */
151  char err_msg[TNT_STRLEN]; /* buffer for error messages */
152  tntComplex contracted; /* The contracted network */
153 
154  /* Determine the length of the MPS. */
155  numsites = (int) tntMpsLength(sandwich);
156 
157  /* Deal with site_left and site_right */
158  if ((-1 == site_left) || (-1 == site_right)) {
159  /* set site_left and site_right to the beginning and end of the network if a full contraction is indicated */
160  site_left = 0;
161  site_right = numsites -1;
162  } else {
163  /* otherwise check that the values given are sensible, and modify network appropriately */
164 
165  /* Check site_left is not more than site_right */
166  if (site_left > site_right) {
167  sprintf(err_msg,"The site for site_left (%d) is greater than the site for site_right (%d) which is not allowed",site_left, site_right); /* NO_COVERAGE */
168  tntErrorPrint(err_msg,TNT_MPS_ERRCODE); /* NO_COVERAGE */
169  } /* NO_COVERAGE */
170 
171  /* Check site_right is not too big */
172  if (site_right >= numsites) {
173  sprintf(err_msg,"The site for site_right (%d) is more than or equal to the length of the MPS (%d) which is not allowed", site_right, numsites); /* NO_COVERAGE */
174  tntErrorPrint(err_msg,TNT_MPS_ERRCODE); /* NO_COVERAGE */
175  } /* NO_COVERAGE */
176 
177  /* If the leftmost site required is more than zero, discard left parts of the MPS network */
178  if (site_left > 0) {
179 
180  /* First move to the correct site - the last of the sites to be discarded */
181  A = tntNodeFindFirst(sandwich);
182  for (sitecount = 1; sitecount < site_left; sitecount++) {
183  A = tntNodeFindConn(A, TNT_MPS_R);
184  }
185  Af = tntNodeFindConn(A, TNT_MPS_D);
186 
187  /* A and Af are the nodes that need to be split that will be the last nodes of the left network (to be discarded). */
188  nodes_L[0] = A;
189  nodes_L[1] = Af;
190 
191  /* find the node after tn. This will now belong to the right network, and will be the new first node in the network to be contracted */
192  nodes_R[0] = A = tntNodeFindConn(A, TNT_MPS_R);
193 
194  /* Find the corresponding flipped node. This will also belong to the right network */
195  nodes_R[1] = Af = tntNodeFindConn(A, TNT_MPS_D);
196 
197  /* Split the network into two: label the first network as the network to keep (right hand) and the second network the network to discard (left hand). */
198  mps_disc = tntNetworkSplit(sandwich, nodes_R[0], 6, tntNodeFindLast(sandwich), TNT_MPS_R, tntNodeFindFirst(sandwich), TNT_MPS_L, nodes_L[0], 4, 2, nodes_R, 2, nodes_L);
199 
200  /* Renumber the starting leg id */
201  leg_s = 6;
202 
203  /* Discard the right network */
204  tntNetworkFree(&mps_disc);
205 
206  /* Join the left internal legs of the flipped and unflipped MPS node on the leftmost site */
208  }
209 
210  /* If the right site required is less than the last site, discard right parts of the MPS Network */
211  if (site_right < numsites-1) {
212 
213  /* First move to the correct site - the first of the sites to be discarded */
214  A = tntNodeFindLast(sandwich);
215  for (sitecount = numsites-2; sitecount > site_right; sitecount--) {
216  A = tntNodeFindConn(A, TNT_MPS_L);
217  }
218  Af = tntNodeFindConn(A, TNT_MPS_D);
219 
220  /* A and Af are the nodes that need to be split that will be part of the right network (to be discarded) */
221  nodes_R[0] = A;
222  nodes_R[1] = Af;
223 
224  /* Find the node before A. This will now belong to the left network (to be kept). */
225  nodes_L[0] = A = tntNodeFindConn(A, TNT_MPS_L);
226 
227  /* Find the corresponding flipped node. This will also belong to the right network */
228  nodes_L[1] = Af = tntNodeFindConn(A, TNT_MPS_D);
229 
230  /* Split the network into two: label the first network as the network to keep (left hand) and the second network the network to discard (right hand). */
231  mps_disc = tntNetworkSplit(sandwich, tntNodeFindFirst(sandwich), leg_s, A, 7, nodes_R[0], 4, tntNodeFindLast(sandwich), TNT_MPS_R, 2, nodes_L, 2, nodes_R);
232 
233  /* Discard the right network */
234  tntNetworkFree(&mps_disc);
235 
236  /* Join the right internal legs of the flipped and unflipped MPS node on the rightmost site */
238 
239  }
240 
241  }
242 
243  /* Now that the unrequired parts of the network have been discarded, zip throught the network from left to right contracting nodes. */
244 
245  sitecount = site_left;
246 
247  /* identify the first node */
248  A = tntNodeFindFirst(sandwich);
249 
250  /* identify the corresponding flipped node Af */
251  Af = tntNodeFindConn(A, TNT_MPS_D);
252 
253  /* Perform the first contraction to form the alpha tensor */
254  alpha = tntNodeContract(A, Af, idmap_mps, idmap_fl);
255 
256  /* Go through remaining sites, contracting them with alpha one by one */
257  for (sitecount = site_left+1; sitecount <= site_right; sitecount++) {
258 
259  /* Find the next node along in the unflipped wave function */
260  A = tntNodeFindConn(alpha, TNT_MPS_R);
261 
262  /* identify the flipped node tnf */
263  Af = tntNodeFindConn(A, TNT_MPS_D);
264 
265  /* Contract alpha with the original wave function */
266  alpha = tntNodeContract(alpha,A,NULL,idmap_mps);
267 
268  /* Contract the alpha and the unflipped nodes */
269  alpha = tntNodeContract(alpha,Af,NULL,idmap_fl);
270  }
271 
272  /* Get the one remaining element */
273  contracted = tntNodeGetFirstEl(alpha);
274 
275  /* Free the network */
276  tntNetworkFree(&sandwich);
277 
278  /* Return the result of the contraction */
279  return contracted;
280 
281 }
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
struct tnode * tntNode
Definition: tnt.h:118
tntNode tntNodeFindLast(tntNetwork tnw)
Definition: tntNode_funcs.c:685
void tntNetworkFree(tntNetwork *tnwp)
Definition: tntNetwork_funcs.c:84
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
#define TNT_MPS_L
Definition: tntMps.h:56
struct tnetwork * tntNetwork
Definition: tnt.h:123
tntComplex tntNodeGetFirstEl(tntNode tn)
Definition: tntNode_funcs.c:770
#define TNT_STRLEN
Definition: tnt.h:59
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
void tntNodeJoin(tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:462
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
Use this type to define complex numbers. Any complex arithmetic needs to be done manually. For example to add two complex numbers A and B to give C write:
Definition: tnt.h:136
tntNetwork tntNetworkSplit(tntNetwork tnw, tntNode tnFirst_1, unsigned leg_start_1, tntNode tnLast_1, unsigned leg_end_1, tntNode tnFirst_2, unsigned leg_start_2, tntNode tnLast_2, unsigned leg_end_2, unsigned numSplit_1, tntNode *tnSplit_1, unsigned numSplit_2, tntNode *tnSplit_2)
Definition: tntNetwork_funcs.c:178
tntComplex tntMpsMpsProduct ( tntNetwork  mpsA,
tntNetwork  mpsB 
)

Calculates the inner product of two MPSs by first creating a 'flipped' copy joined along the physical legs, then contracting along the network.

mpsA_mpsB.png

If the original MPSs mpsA, mpsB represent a wave functions \(|\psi_\mathrm{A}\rangle\), \(|\psi_\mathrm{B}\rangle\) then this represents the overlap \(\langle\psi_\mathrm{B}|\psi_\mathrm{A}\rangle\).

If the result of the contraction is real, the value itself is returned. If the result of the contraction is complex, then the absolute value is returned.

See also
If the inner product of the MPS with itself is required, use tntMpsSelfProduct().
Returns
The inner product of the two MPSs, equivalent to \(\langle\psi_\mathrm{B}|\psi_\mathrm{A}\rangle\) if the MPSs represent wave functions, with the absolute value taken if the result is complex.
Parameters
mpsAThe network representing the first MPS. Unchanged by the function.
mpsBThe network representing the second MPS. Unchanged by the function.
Examples:
tebd.c.
303 {
304 
305  tntNetwork sandwich; /* The network that will be contracted to find the product */
306 
307  /* `wfa` will be deleted during the network contraction, so use a copy so that the original wave function is untouched. */
308  sandwich = tntNetworkCopy(mpsA);
309 
310  /* First flip the network, and connect to mpsB (this function leave mpsB unchanged and connects a copy of mpsB) */
311  tntMpsMpsConnect(sandwich, mpsB);
312 
313  /* Now contract the nodes, performing a full contraction */
314  return tntMpsMpsContract(sandwich,-1,-1);
315 
316 }
tntComplex tntMpsMpsContract(tntNetwork sandwich, int site_left, int site_right)
Definition: tntMpsMps.c:128
void tntMpsMpsConnect(tntNetwork mpsA, tntNetwork mpsB)
Definition: tntMpsMps.c:38
struct tnetwork * tntNetwork
Definition: tnt.h:123
tntNetwork tntNetworkCopy(tntNetwork tnw)
Definition: tntNetwork_funcs.c:36
void tntMpsOrth ( tntNetwork  mps,
unsigned  orth_centre 
)

Sweeps through the MPS performing a sequence of (untruncated) SVDs and contractions, to leave the every node apart from the orthogonality centre either left or right orthonormalised. After applying the routine, every node to the left of the orthogonality centre will obey the left orthonormalisation condition:

left_orth.png

and every node to the right of the orthogonality centre will obey the right orthonormalisation condition:

right_orth.png
Returns
No return value

The truncation tolerance is temporarily turned off before being returned to its initial value during this routine, so the truncation error will be exactly zero and is not returned. If truncation to the global truncation tolerance (set using tntSVDTruncTolSet()) is required, use tntMpsTruncate() with chi set to -1.

Parameters
mpsNetwork representing the MPS. Changed by function to be either left or right orthonormalised
orth_centreThe required orthogonality centre of the MPS.
38 {
39 
40  double err; /* truncation error from the orthonormalisation - not used as it will be exactly zero. */
41  double trunc_tol; /* current value of the truncation tolerance */
42 
43  /* Get the current truncation tolerance */
44  trunc_tol = tntSVDTruncTolGet();
45 
46  /* Turn the truncation tolerance off */
47  tntSVDTruncTolSet(-1.0);
48 
49  /* Call the routine to truncate, which will now simply orthogonalise since it will perform only untruncated SVDs */
50  err = tntMpsTruncate(mps, orth_centre, -1);
51 
52  /* we don't need the error. This line avoids compilation warnings */
53  err = (double) err;
54 
55  /* Set the truncation tolerance back to it's initial value. */
56  tntSVDTruncTolSet(trunc_tol);
57 
58  return;
59 
60 }
double tntSVDTruncTolGet(void)
Definition: tntMisc_funcs.c:587
double tntMpsTruncate(tntNetwork mps, unsigned orth_centre, int chi)
Definition: tntMpsTruncate.c:44
void tntSVDTruncTolSet(double truncTol)
Definition: tntMisc_funcs.c:569
tntComplex tntMpsPmpoMpsProduct ( tntNetwork  mps,
tntNodeArray op,
tntIntArray sitenum,
int  orth_centre 
)

Calculates the inner product of an MPS-MPS network with zero or more single-site operators op[i] representing a product MPO (PMPO) with the following form inserted between them.

single_site_op.png

This set of single site operators represent a site-wide MPO formed from a tensor product of operators on every site, where the operator on any site not given is the identity operator. Each of the specified operators is inserted on the relevant site, then contracted with the physical leg on the MPS to form a new MPS. The new MPS is then contracted with the original MPS which is flipped to form the bottom of the network. This function therefore represents an MPS-PMPO-MPS product.

mps_pmpo_mps.png

If multiple operators are given for a single site, then they will be applied to the wave function in the order they are given in the array i.e. in the figure above the operators given first will be on top.

The function can be made more efficient by specifying the orthogonality centre of the MPS. All sites to the right of the orthogonality centre should obey the right orthonormalisation condition:

right_orth.png

All sites to the left of the orthogonality centre should obey the left orthonormalisation condition:

left_orth.png

This means that, for example, if the orthogonality centre is site 2 the network contraction above can be simplified to:

expec_sing_part.png

This then reduces the number of contractions required to find the expectation value.

Note that after applying tntMpsPropST2scProduct(), or tntMpsVarMinMpo2sSweep() in direction TNT_MPS_L, the orthonormality centre will be on the first site. After applying tntMpsVarMinMpo2sSweep() in direction TNT_MPS_R, the orthonormality centre will be on the last site.

Returns
The product of the MPS-MPS network with the product MPO inserted between them. If the value is complex, then the absolute value is returned instead.
Parameters
mpsThe network representing the MPS. Unchanged by the function.
opArray giving the non-identity single-site operators which define the product MPO. Operators are assumed to have 2 legs, with upwards facing leg having number 4, and downwards facing leg having number 3. Unchanged by the function.
sitenumArray giving the site numbers for each of the operators given in op. Site numbers count from 0. They do not need to be in a particular order, but they should be in the same order as op. Unchanged by the function.
orth_centreThe orthogonality centre of the MPS wave function. Set to -1 if unknown.
63 {
64 
65  tntNetwork sandwich; /* The full MPS-MPO-MPS network */
66  int numsites; /* Number of sites in the MPS */
67  int site_left, site_right; /* Minimum and maximum site numbers for contracting the network */
68  unsigned loop; /* Used for looping over the elements in the arrays */
69 
70  char err_msg[TNT_STRLEN]; /* buffer for error message */
71 
72  /* Find length of the MPS */
73  numsites = (int) tntMpsLength(mps);
74 
75  /* Check orth_centre is not too big */
76  if (orth_centre >= numsites) {
77  sprintf(err_msg,"The site for orth_centre (%d) is more than or equal to the length of the MPS (%d) which is not allowed", orth_centre, numsites); /* NO_COVERAGE */
78  tntErrorPrint(err_msg,TNT_MPS_ERRCODE); /* NO_COVERAGE */
79  } /* NO_COVERAGE */
80 
81  /* the full network will be altered then deleted during the network contraction, so use a copy so that the original mps is untouched. */
82  sandwich = tntNetworkCopy(mps);
83 
84  /* Take the product of the MPS with the PMPO */
85  tntMpsPmpoProduct(sandwich,op,sitenum);
86 
87  /* Connect the altered MPS to a flipped copy of the original MPS */
88  tntMpsMpsConnect(sandwich,mps);
89 
90  /* Initialise values for finding the minimum and maximum site number in the sitenum array if orth_centre is not -1 */
91  /* (if it is -1 this is not necessary since a full contraction will be performed, therefore the position of the operators is irrelevant) */
92  if (orth_centre > 0) {
93  site_left = site_right = orth_centre;
94  } else {
95  site_left = -1;
96  site_right = -1;
97  }
98 
99  /* Check that the values given for the site numbers are sensible. At the same time, find the minimum and maximum site number that has an operator on */
100  for (loop = 0; loop < op->sz; loop++) {
101  if (sitenum->vals[loop] < 0) {
102  sprintf(err_msg,"The site (%d) for operator %d is less than zero", sitenum->vals[loop], loop); /* NO_COVERAGE */
103  tntErrorPrint(err_msg,TNT_MPS_ERRCODE); /* NO_COVERAGE */
104  } else if (sitenum->vals[loop] >= numsites) {
105  sprintf(err_msg,"The site (%d) for operator %d is more than or equal to the length of the MPS (%d)", sitenum->vals[loop], loop, numsites); /* NO_COVERAGE */
106  tntErrorPrint(err_msg,TNT_MPS_ERRCODE); /* NO_COVERAGE */
107  } else if (orth_centre > -1) {
108  /* if orth_centre > -1, then find the minimum and maximum site number so far encountered */
109  if (sitenum->vals[loop] < site_left) {
110  site_left = sitenum->vals[loop];
111  }
112  if (sitenum->vals[loop] > site_right) {
113  site_right = sitenum->vals[loop];
114  }
115  }
116  }
117 
118  /* Contract the full network */
119  return tntMpsMpsContract(sandwich,site_left,site_right);
120 }
tntComplex tntMpsMpsContract(tntNetwork sandwich, int site_left, int site_right)
Definition: tntMpsMps.c:128
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
void tntMpsMpsConnect(tntNetwork mpsA, tntNetwork mpsB)
Definition: tntMpsMps.c:38
void tntMpsPmpoProduct(tntNetwork mps, tntNodeArray *op, tntIntArray *sitenum)
Definition: tntMpsPmpo.c:35
struct tnetwork * tntNetwork
Definition: tnt.h:123
#define TNT_STRLEN
Definition: tnt.h:59
tntNetwork tntNetworkCopy(tntNetwork tnw)
Definition: tntNetwork_funcs.c:36
unsigned sz
The number of elements in the array.
Definition: tnt.h:180
void tntMpsPmpoProduct ( tntNetwork  mps,
tntNodeArray op,
tntIntArray sitenum 
)

Applies on more single-site operators op[i] to the physical legs of an MPS.

single_site_op.png

This set of single site operators represent a site-wide MPO formed from a tensor product of operators on every site, where the operator on any site not given is the identity operator. Each of the specified operators is inserted on the relevant site, then contracted with the physical leg on the MPS. This function therefore represents an MPS-PMPO product.

mps_pmpo.png

This results in another MPS, which replaces the MPS network passed as an argument. The physical legs of the MPS passed as an argument must be free, or an error will result.

Note that the single site operators applied can cause the dimension of the phsyical legs on the MPS to change, if the upwards and downwards facing legs have different dimensions.

Returns
No return value, but the input parameter mps is altered by having its product taken with the PMPO.
Parameters
mpsThe network representing the MPS. Remains an MPS network, but with altered nodes, after the function is applied.
opArray giving the non-identity single-site operators which define the product MPO. Operators are assumed to have 2 legs, with upwards facing leg having number 4, and downwards facing leg having number 3. Unchanged by the function.
sitenumArray giving the site numbers for each of the operators given in op. Site numbers count from 0. They do not need to be in a particular order, but they should be in the same order as op. Unchanged by the function.
42 {
43  /* Networks and nodes during the contraction */
44  tntNode A; /* The current MPS node */
45  tntNode opn; /* A copy of the operator to use for contracting */
46 
47  /* Counters and site information used during the contraction */
48  int sitecount; /* the current site count */
49  int numsites; /* The number of MPS nodes in the network */
50  unsigned opnum=0; /* Counter for the operators */
51  tntIntArray siterank; /* Rank for each site-number, to make then in ascending order */
52  tntIntArray sitenum_cp; /* Copy of the site numbers which will be used for sorting */
53  unsigned i, j; /* Used for looping */
54  int sitenummax, sitemax; /* Position of current maximum site, current maximum site */
55 
56  /* Other */
57  char err_msg[TNT_STRLEN]; /* buffer for error messages */
58 
59  /* Determine the length of the MPS. */
60  numsites = (int) tntMpsLength(mps);
61 
62  /* Check that operators if they have been provided */
63  if ((NULL == op) || (NULL == sitenum) || (0 == op->sz) || (0 == sitenum->sz)) {
64  tntWarningPrint("Leaving function tntMpsMpoProduct() with no action taken as no operators or sitenumbers have been provided"); /* NO_COVERAGE */
65  return; /* NO_COVERAGE */
66  } /* NO_COVERAGE */
67 
68  /* Allocate arrays for copies of the site numbers, and rank of site numbers */
69  sitenum_cp = tntIntArrayAlloc(op->sz);
70  siterank = tntIntArrayAlloc(op->sz);
71 
72  /* Check that the sites and orthogonality centre given are sensible, and at the same time copy over site numbers */
73  for (opnum = 0; opnum < op->sz; opnum++) {
74  if (sitenum->vals[opnum] >= numsites) {
75  sprintf(err_msg,"The site for the operator %d is %d, which is greater than the length of the MPS (%d).", opnum, sitenum->vals[opnum], numsites); /* NO_COVERAGE */
76  tntErrorPrint(err_msg, TNT_MPS_ERRCODE); /* NO_COVERAGE */
77  } /* NO_COVERAGE */
78 
79  if (sitenum->vals[opnum] < 0) {
80  sprintf(err_msg,"The site for the operator %d is %d, which is not allowed - site numbers should be zero or more", opnum, sitenum->vals[opnum]); /* NO_COVERAGE */
81  tntErrorPrint(err_msg, TNT_MPS_ERRCODE); /* NO_COVERAGE */
82  } /* NO_COVERAGE */
83 
84  sitenum_cp.vals[opnum] = sitenum->vals[opnum];
85  }
86 
87  /* First give each element of sitenum a rank, where siterank i gives the element of ith biggest site */
88  /* Loop through all legs from end to beginning. For each index bl, find the blth largest leg dimension */
89  for (i = opnum; i-->0; ) {
90  /* Reset marker for largest site number */
91  sitemax = -1;
92  for (j = opnum; j-- > 0; ) {
93  if (sitenum_cp.vals[j] > sitemax) {
94  sitemax = sitenum_cp.vals[j];
95  sitenummax = j;
96  }
97  }
98  /* Marks current largest site as being in jth position in array, and set its dimension to 0 to indicate that it no longer requires sorting */
99  sitenum_cp.vals[sitenummax] = -1;
100  siterank.vals[i] = sitenummax;
101  }
102 
103  /* Copy of site numbers no longer required */
104  tntIntArrayFree(&sitenum_cp);
105 
106  /* Reset counter for number of operators */
107  opnum = 0;
108 
109  /* Move through the network from left to right. On each site, check whether there is an operator to apply, and if contract it with the physical legs */
110 
111  /* identify the first node */
112  A = tntNodeFindFirst(mps);
113 
114  /* Go through remaining sites, contracting them with alpha one by one */
115  for (sitecount = 0; sitecount < numsites; sitecount++) {
116 
117  /* If there are operators on this site, join to the physical leg and contract with A */
118  while ((opnum < op->sz) && (sitenum->vals[siterank.vals[opnum]] == sitecount)) {
119  /* Make a copy of the operator */
120  opn = tntNodeCopy(op->vals[siterank.vals[opnum]],0);
121 
122  /* join the operator to the physical leg */
123  tntNodeJoin(opn, TNT_MPS_U, A, TNT_MPS_D);
124 
125  /* contract the first node with the operator on the physical leg */
126  A = tntNodeContract(A,opn,NULL,NULL);
127 
128  /* increment the counter for the operator number */
129  opnum++;
130  }
131 
132  /* Break if all the operators have been applied */
133  if (opnum == op->sz) {
134  break;
135  }
136 
137  /* Find the next node along in the MPS */
138  A = tntNodeFindConn(A, TNT_MPS_R);
139  }
140 
141  /* Free the site rank */
142  tntIntArrayFree(&siterank);
143 
144  return;
145 
146 }
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:70
Simple 1D array for integer values.
Definition: tnt.h:148
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
struct tnode * tntNode
Definition: tnt.h:118
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
#define TNT_MPS_U
Definition: tntMps.h:71
tntNode * vals
Pointer to the first element in the tntNode array.
Definition: tnt.h:179
#define TNT_STRLEN
Definition: tnt.h:59
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
void tntNodeJoin(tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:462
#define TNT_MPS_D
Definition: tntMps.h:66
#define TNT_MPS_R
Definition: tntMps.h:61
unsigned sz
The number of elements in the array.
Definition: tnt.h:150
unsigned sz
The number of elements in the array.
Definition: tnt.h:180
double tntMpsSelfProduct ( tntNetwork  mps,
int  orth_centre 
)

Calculates the self inner product of an MPS by first creating a 'flipped' copy joined along the physical legs, then contracting along the network.

mpsA_mpsA.png

If the MPS represents the wave function \(\psi\rangle\), then the product is \(\langle\psi|\psi\rangle\) i.e. the square of the 2-norm.

The function can be made more efficient by specifying the orthogonality centre of the MPS. All sites to the right of the orthogonality centre should obey the right orthonormalisation condition:

right_orth.png

All sites to the left of the orthogonality centre should obey the left orthonormalisation condition:

left_orth.png

This means that, for example, if the orthogonality centre is 0 the network contraction above can be simplified to:

norm_part.png

This then reduces the number of contractions required to find the product.

Note that after applying tntMpsPropST2scProduct(), or tntMpsVarMinMpo2sSweep() in direction TNT_MPS_L, the orthonormality centre will be on the first site. After applying tntMpsVarMinMpo2sSweep() in direction TNT_MPS_R, the orthonormality centre will be on the last site.

If the orthogonality centre is given as -1, then a full contraction is performed.

Returns
The self inner product of the MPS.
Parameters
mpsThe network representing the MPS. Unchanged by the function.
orth_centreThe orthogonality centre of the MPS wave function. Set to -1 if unknown.
357 {
358 
359  /* Networks and nodes during the contraction */
360  tntNetwork mpsmps; /* The MPS-MPS network that will be contracted to find the expectation value */
361  tntComplex norm;/* The expectation value: real part is returned by the function. */
362  int numsites; /* Number of sites in the MPS */
363 
364  char err_msg[TNT_STRLEN]; /* buffer for error message */
365 
366  /* Find length of the MPS */
367  numsites = (int) tntMpsLength(mps);
368 
369  /* Check orth_centre is not too big */
370  if (orth_centre >= numsites) {
371  sprintf(err_msg,"The site for orth_centre (%d) is more than or equal to the length of the MPS (%d) which is not allowed", orth_centre, numsites); /* NO_COVERAGE */
372  tntErrorPrint(err_msg,TNT_MPS_ERRCODE); /* NO_COVERAGE */
373  } /* NO_COVERAGE */
374 
375  /* The wave function will be deleted during the network contraction, so use a copy so that the original wave function is untouched. */
376  mpsmps = tntNetworkCopy(mps);
377 
378  /* First flip the network, and connect to itself */
379  tntMpsMpsConnect(mpsmps, mps);
380 
381  /* Now contract the network - this destroys the mpsmps network */
382  norm = tntMpsMpsContract(mpsmps, orth_centre, orth_centre);
383 
384  /* return the real part of the norm (no need to check whether there is an imaginary part as inner product will always be real */
385  return (norm.re);
386 }
tntComplex tntMpsMpsContract(tntNetwork sandwich, int site_left, int site_right)
Definition: tntMpsMps.c:128
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
void tntMpsMpsConnect(tntNetwork mpsA, tntNetwork mpsB)
Definition: tntMpsMps.c:38
struct tnetwork * tntNetwork
Definition: tnt.h:123
#define TNT_STRLEN
Definition: tnt.h:59
double re
Definition: tnt.h:137
tntNetwork tntNetworkCopy(tntNetwork tnw)
Definition: tntNetwork_funcs.c:36
Use this type to define complex numbers. Any complex arithmetic needs to be done manually. For example to add two complex numbers A and B to give C write:
Definition: tnt.h:136
double tntMpsTruncate ( tntNetwork  mps,
unsigned  orth_centre,
int  chi 
)

Sweeps through the MPS performing a sequence of truncated SVDs and contractions, to leave the every node apart from the orthogonality centre either left or right orthonormalised.

It does not assume anything about the initial orthonormality properties of the MPS. It simply starts from the left-most site performing truncated SVDs until the orthogonality centre is reached (obviously this has no effect if the orthogonality centre is the first site). It then starts from the right-most site performing truncated SVDs until the orthogonality centre is reached (obviously this has no effect if the orthogonality centre is the last site).

After applying the routine, every node to the left of the orthogonality centre will obey the left orthonormalisation condition:

left_orth.png

and every node to the right of the orthogonality centre will obey the right orthonormalisation condition:

right_orth.png

To internal dimension will be set by the minimum of the dimension chi given as an argument or the one that results from the global truncation tolerance (set using tntSVDTruncTolSet()). To truncate to the global truncation tolerance only, set the argument chi to -1.

Returns
The sum of the truncation errors of all the SVDs.
See also
Use tntMpsOrth() for a routine that orthogonalises the MPS without any truncation (i.e. with the truncation tolerance turned off).
Parameters
mpsNetwork representing the MPS. Changed by function to be either left or right orthonormalised and to have a maximum internal dimension of chi.
orth_centreThe required orthogonality centre of the MPS.
chiMaximum internal dimension of the new MPS. Use -1 to truncate to the truncation tolerance only.
47 {
48 
49  int start_legs[2] = {1,3}; /* Array for holding the types of the 'start legs' of the SVD */
50  tntNode A; /* the current node */
51  tntNode A_nxt; /* the next node along */
52  tntNode S, VT; /* nodes that result from taking SVD */
53  /* The maps for leg types: all leg types stay the same, so just needs to be consecutive */
54  unsigned numsites, j; /* Number of nodes and node counter in the MPS network */
55  char err_msg[TNT_STRLEN]; /* buffer for error messages */
56  double err, err_tot = 0.0; /* truncation error per SVD and total truncation error */
57 
58  /* Check the the number of nodes in the MPS network */
59  numsites = tntMpsLength(mps);
60 
61  /* Give an error if orth_centre is more than index of the last site */
62  if (orth_centre > numsites-1) {
63  sprintf(err_msg,"The site for the orthogonality centre (%d) is greater than the length of the MPS (%d)",orth_centre,numsites); /* NO_COVERAGE */
64  tntErrorPrint(err_msg,TNT_MPS_ERRCODE); /* NO_COVERAGE */
65  } /* NO_COVERAGE */
66 
67  /* Loop from first site to orthogonality centre, leaving the `twist' at the orthogonality centre */
68  A = tntNodeFindFirst(mps);
69  for (j = 0; j < orth_centre; j++) {
70  /* Choose take truncated SVD and add truncation error to total value. */
71  A = tntNodeSVD(A,start_legs,2,TNT_MPS_R,TNT_MPS_L,NULL,chi,&err);
72  err_tot += err;
73 
74  /* A is assigned to U so is now normalised. Contract S and VT with next node along. */
75 
76  /* Find S by finding node connected to the leg of the correct leg type, and assign to A */
78 
79  /* Find VT by finding node connected to the leg of the correct leg type, and assign to A_nxt */
80  VT = tntNodeFindConn(S,TNT_MPS_R);
81 
82  /* Contract S and VT and assign the result to A. */
83  A = tntNodeContract(S,VT,NULL,NULL);
84  /* A now points to the contracted node SVT */
85 
86  /* determine the next MPS node along */
87  A_nxt = tntNodeFindConn(A,TNT_MPS_R);
88 
89  /* Contract SVT with the next A matrix along and assign the result to A */
90  A = tntNodeContract(A,A_nxt,NULL,NULL);
91  }
92 
93 
94  /* Loop from the last site to the orthogality centre, leaving the `twist' at the orthogonality centre */
95  A = tntNodeFindLast(mps);
96  start_legs[0] = 2; /* Change the leg assigned to U in the SVD */
97 
98  for (j = numsites-1; j > orth_centre; j--) {
99  /* Choose take truncated SVD and add truncation error to total value. */
100  A = tntNodeSVD(A,start_legs,2,TNT_MPS_L,TNT_MPS_R,NULL,chi,&err);
101  err_tot += err;
102  /* A is assigned to U so is now normalised. Contract S and VT with next node along. */
103 
104  /* Find S by finding node connected to the leg of the correct leg type, and assign to S */
105  S = tntNodeFindConn(A,TNT_MPS_L);
106 
107  /* Find VT by finding node connected to the leg of the correct leg type, and assign to VT */
108  VT = tntNodeFindConn(S,TNT_MPS_L);
109 
110  /* Contract S and VT and assign the result to A. */
111  A = tntNodeContract(S,VT,NULL,NULL);
112  /* A now points to the contracted node SVT */
113 
114  /* determine the next MPS node along */
115  A_nxt = tntNodeFindConn(A,TNT_MPS_L);
116 
117  /* Contract SVT with the next A matrix along and assign the result to A */
118  A = tntNodeContract(A,A_nxt,NULL,NULL);
119  }
120 
121  return err_tot;
122 
123 }
tntNode tntNodeFindConn(tntNode tn, unsigned legA)
Definition: tntNode_funcs.c:438
unsigned tntMpsLength(tntNetwork mps)
Definition: tntMpsUtil.c:208
struct tnode * tntNode
Definition: tnt.h:118
tntNode tntNodeFindLast(tntNetwork tnw)
Definition: tntNode_funcs.c:685
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
#define TNT_MPS_L
Definition: tntMps.h:56
tntNode tntNodeSVD(tntNode tn, int *startLegs, unsigned startNum, int legUS, int legSU, int *legMap, int connDim, double *err)
Definition: tntNode_funcs.c:64
#define TNT_STRLEN
Definition: tnt.h:59
tntNode tntNodeFindFirst(tntNetwork tnw)
Definition: tntNode_funcs.c:664
#define TNT_MPS_R
Definition: tntMps.h:61