Tensor Network Theory Library  Beta release 1.2.1
A library of routines for performing TNT-based operations
 All Data Structures Functions Variables Groups Pages
tntMpsCreateSymmRandom.c
1 /*
2 Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
3 $LastChangedDate$
4 (c) University of Oxford 2014
5 */
6 
7 /* Include the header for the TNT library */
8 #include "tntMpsInternal.h"
9 
33 tntNetwork tntMpsCreateSymmRandom(unsigned L,
34  int *qn)
35 {
36  tntNetwork mps; /* The mps network to create */
37  tntNode A; /* The current node */
38  unsigned d; /* Dimension of the physical leg */
39  unsigned loop; /* used for looping */
40  unsigned legdim[2] = {1,1}; /* Leg dimensions for the internal legs */
41  unsigned i, j, k, qind; /* used for looping and indexing*/
42  tntIntArray dqn; /* difference between succesive quantum numbers */
43  tntIntArray qmin, qmax; /* minimum and maximum quantum numbers for each quantum number label on the physical leg */
44  tntIntArray qminsum, qmaxsum; /* minimum and maximum possible quantum numbers for the internal leg, to allow the boundary qn's to be reached */
45  tntIntArray qnums_phys; /* Quantum numbers for the physical leg */
46  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 */
47  unsigned numqn, lblnum = tntSymmNumGet(); /* current number of qns, number of qns in a label */
48  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 */
49  double trunc_tol; /* current value of the truncation tolerance */
50 
51  /* Get the current truncation tolerance */
52  trunc_tol = tntSVDTruncTolGet();
53 
54  /* Turn the truncation tolerance off so all SVDs are exact */
55  tntSVDTruncTolSet(-1.0);
56 
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"); /* NO_COVERAGE */
69  } else { /* NO_COVERAGE */
71  }
72 
73  /* Get the quantum numbers for the physical leg */
74  qnums_phys = tntNodeGetQN(tntSysBasisOpGet(), "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]) qmin.vals[i] = qnums_phys.vals[qind];
96 
97  /* check if current index has a larger qn, and reassign if it does */
98  if (qmax.vals[i] < qnums_phys.vals[qind]) qmax.vals[i] = qnums_phys.vals[qind];
99 
100  /* if dqn not yet assigned assign difference between this qn and the first different qn */
101  /* note we are assuming that all differences are the same, or at least come in increasing order */
102  if (0 == dqn.vals[i]) dqn.vals[i] = abs(qnums_phys.vals[qind]-qnums_phys.vals[i]);
103 
104  } /* end of loop through indices on physical leg */
105  } /* end of loop through qn's in a label */
106 
107  /* check that the total quantum number given actually makes sense */
108  for (i = 0; i < lblnum; i++) { /* loop through the quantum number labels */
109 
110  /* total qn must be more than or equal to (minimum qn on physical leg)*(number of physical legs) */
111  if (qn[i] < ((int) L)*qmin.vals[i]) {
112  if (lblnum > 1) { /* NO_COVERAGE */
113  tntErrorPrint("Total quantum number %d (%d of %d qn labels) is invalid|"
114  "Minimum allowable qn for this label is %d",
115  qn[i], i+1, lblnum, L*qmin.vals[i]); /* NO_COVERAGE */
116  } else { /* NO_COVERAGE */
117  tntErrorPrint("Total quantum number %d is invalid|"
118  "Minimum allowable qn is %d",
119  qn[i], L*qmin.vals[i]); /* NO_COVERAGE */
120  } /* NO_COVERAGE */
121  } /* NO_COVERAGE */
122 
123  /* total qn must be less than or equal to (maximum qn on physical leg)*(number of physical legs) */
124  if (qn[i] > ((int) L)*qmax.vals[i]) {
125  if (lblnum > 1) { /* NO_COVERAGE */
126  tntErrorPrint("Total quantum number %d (%d of %d qn labels) is invalid|"
127  "Maximum allowable qn for this label is %d",
128  qn[i], i+1, lblnum, L*qmax.vals[i]); /* NO_COVERAGE */
129  } else { /* NO_COVERAGE */
130  tntErrorPrint("Total quantum number %d is invalid|"
131  "Maximum allowable qn is %d",
132  qn[i], L*qmax.vals[i]); /* NO_COVERAGE */
133  } /* NO_COVERAGE */
134  } /* NO_COVERAGE */
135 
136  /* total qn must be a qn of one of the configurations written in terms of the basis operator.
137  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.
138  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 */
139  if ((dqn.vals[i] !=0 && (L*qmax.vals[i]-qn[i])%dqn.vals[i] != 0) || (0 == dqn.vals[i] && qn[i] != L*qmax.vals[i])) {
140  if (lblnum > 1) { /* NO_COVERAGE */
141  tntErrorPrint("Total quantum number %d (%d of %d qn labels) is invalid|"
142  "There are no configurations in a %d site system that have this qn",
143  qn[i], i+1, lblnum, L); /* NO_COVERAGE */
144  } else { /* NO_COVERAGE */
145  tntErrorPrint("Total quantum number %d is invalid"
146  "There are no configurations in a %d site system that have this qn",
147  qn[i], L); /* NO_COVERAGE */
148  } /* NO_COVERAGE */
149  } /* NO_COVERAGE */
150  }
151 
152  /* Now create the empty network */
153  mps = tntNetworkCreate();
154 
155  /* Allocate array for quantum numbers for the left internal leg. Initial value will be zero as required. */
156  qnums_int[0] = tntIntArrayAlloc(lblnum);
157 
158  /* Allocate arrays for the minimum and maximum allowable qn on each site to allow the total to be reached */
159  qminsum = tntIntArrayAlloc(lblnum);
160  qmaxsum = tntIntArrayAlloc(lblnum);
161 
162  /* First left leg will have dimension 1 */
163  legdim[0] = 1;
164 
165  /* Keep inserting random nodes at the end. */
166  for (loop = 0; loop < L; loop++) { /* loop over sites */
167 
168  /* --- Determine the quantum numbers for the right leg ---- */
169 
170  /* Allocate array for the possible quantum numbers for right internal leg. d*D is maximum number of possible labels for the first node */
171  qnums_int[(loop+1)%2] = tntIntArrayAlloc(d*lblnum*legdim[0]);
172 
173  /* Possible quantum numbers are given by sum of those on incoming physical leg and incoming left leg */
174  /* assign d*D possible qn's */
175 
176  for (i = 0; i < legdim[0]; i++) { /* loop through left internal leg indices */
177  for (j = 0; j < d; j++) { /* loop through physical leg indices */
178  for (k = 0; k < lblnum; k++) { /* loop through the qn's in a label */
179 
180  /* current qn is sum of qn label on physical leg and internal leg */
181  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];
182 
183  } /* end of loop through qn's in a label */
184 
185  } /* end of loop through physical leg indices */
186 
187  } /* end of loop through left internal leg indices */
188 
189  /* Determine minimum and maximum allowable qn for the right internal leg to allow the final qn to be reached */
190  for (i = 0; i < lblnum; i++) { /* loop through the qn's in a label */
191 
192  /* quantum number must be more than or equal to qtotal_i - qmax*(number of sites to come) */
193  qminsum.vals[i] = qn[i] - qmax.vals[i]*(L-(loop+1));
194 
195  /* quantum number must be less than or qual to qtotal_i - qmin*(number of sites to come) */
196  qmaxsum.vals[i] = qn[i] - qmin.vals[i]*(L-(loop+1));
197 
198  } /* end of loop through qn's in a label */
199 
200  /* Now go through quantum numbers and remove any that are not allowed, as well as any duplicates */
201  numqn = 0; /* number of qns we are starting with: 0 assigned so far */
202  for (i = 0; i < d*legdim[0]; i++) { /* loop through new qns */
203 
204  /* First check whether the quantum number is within the range of allowed values */
205  for (j = 0; j < lblnum; j++) { /* loop through the qn's in a label */
206  /* reset flag */
207  qn_valid = 1;
208 
209  /* check the qn is valid */
210  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])) {
211  /* if it isn't, set the flag and break inner loop */
212  qn_valid = 0;
213  break;
214  }
215 
216  } /* end of loop through qn's in a label */
217 
218  /* Now check whether this quantum number has already appeared - loop through all qn currently assigned */
219  for (k = 0; k < numqn; k++) {
220  /* reset flag */
221  qn_appeared = 1;
222 
223  /* loop through every qn in the label, and see if it is the same as the current qn */
224  for (j = 0; j < lblnum; j++) {
225  if (qnums_int[(loop+1)%2].vals[i*lblnum+j] != qnums_int[(loop+1)%2].vals[k*lblnum+j]) {
226  qn_appeared = 0;
227  break;
228  }
229  }
230 
231  /* exit the loop if the qn has been found */
232  if (qn_appeared) {
233  break;
234  }
235 
236  }
237 
238  if (qn_valid&&(!qn_appeared)) {
239  /* if qn is valid and hasn't already appeared, increment the number of unique, allowed qns */
240  numqn++;
241 
242  /* Shuffle down entry if some values have already been removed */
243  if (numqn-1 != i) {
244  for (j = 0; j < lblnum; j++) {
245  qnums_int[(loop+1)%2].vals[(numqn-1)*lblnum+j] = qnums_int[(loop+1)%2].vals[i*lblnum+j];
246  }
247  }
248  }
249  }
250 
251  /* set the right internal leg dimension to be the number of possible qns */
252  legdim[1] = numqn;
253 
254  /* Create the node */
255  A = tntNodeCreate(NULL, "LRD", legdim[0], legdim[1], d);
256 
257  /* --- Set qn information for the node --- */
258  /* Note: no warn version of setting QN info used since we know information is being lost */
259 
260  /* Physical leg is incoming leg and has same qn information as operator */
261  tntNodeSetQN_nowarn(A,"D",&qnums_phys,TNT_QN_IN);
262 
263  /* Left leg is incoming leg, and has values = to right leg of previous tensor */
264  tntNodeSetQN_nowarn(A,"L",qnums_int+loop%2,TNT_QN_IN);
265 
266  /* Now assign these values to the right outgoing leg */
267  tntNodeSetQN_nowarn(A,"R",qnums_int+(loop+1)%2,TNT_QN_OUT);
268 
269  /* free left qns */
270  tntIntArrayFree(qnums_int+loop%2);
271 
272  /* scale the node by its norm */
274 
275  /* Insert A into the network */
276  tntNodeInsertAtEnd(A, "L", "R", mps);
277 
278  /* Next left leg has to have the same dimension as the previous right leg */
279  legdim[0] = legdim[1];
280 
281  } /* end of loop over sites */
282 
283  /* Free the physical and right internal qns */
284  tntIntArrayFree(&qnums_phys);
285  tntIntArrayFree(qnums_int+L%2);
286 
287  /* Free arrays no longer required */
288  tntIntArrayFree(&qmin);
289  tntIntArrayFree(&qmax);
290  tntIntArrayFree(&dqn);
291  tntIntArrayFree(&qminsum);
292  tntIntArrayFree(&qmaxsum);
293 
294  /* orthonormalise the wave function */
295  tntMpsOrthNorm(mps,0);
296 
297  /* Set the truncation tolerance back to it's initial value. */
298  tntSVDTruncTolSet(trunc_tol);
299 
300  return mps;
301 
302 }
303 
tntNetwork tntNetworkCreate(void)
Definition: tntNetwork.c:92
unsigned tntSymmNumGet(void)
Definition: tntSys.c:512
void tntNodeInsertAtEnd(tntNode I, tntLegLabel legIlast, tntLegLabel legInwend, tntNetwork nw)
Definition: tntNodeConn.c:252
void tntSVDTruncTolSet(double truncTol)
Definition: tntSys.c:253
double tntNodeGetNorm(tntNode A)
Definition: tntNodeInfo.c:554
void tntMpsOrthNorm(tntNetwork mps, unsigned orth_centre)
Definition: tntMpsOrth.c:83
tntIntArray tntNodeGetQN(tntNode A, tntLegLabel legA)
Definition: tntNodeQN.c:155
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray.c:620
double tntSVDTruncTolGet(void)
Definition: tntSys.c:312
tntNetwork tntMpsCreateRandom(unsigned L, unsigned D)
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:29
tntNode tntSysBasisOpGet(void)
Definition: tntSys.c:407
tntNode tntNodeCreate(tntComplexArray *nodeVals, tntLegLabel leglabels, unsigned dimleg1, unsigned dim...)
Definition: tntNodeUtil.c:354
void tntNodeSetQN_nowarn(tntNode A, tntLegLabel legA, tntIntArray *qvals, int legdir)
Definition: tntNodeQN.c:105
void tntNodeScaleReal(tntNode A, double val)
int tntSymmTypeGet(void)
Definition: tntSys.c:498
unsigned tntNodeGetLegDim(tntNode A, tntLegLabel legA)
Definition: tntNodeInfo.c:246
tntNetwork tntMpsCreateSymmRandom(unsigned L, int *qn)