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
tntMpsCreateMpo.c
1 /*
2 Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
3 $LastChangedDate$
4 (c) University of Oxford 2014
5 */
6 /* Include the internal header for the TNT library */
7 #include "tntMpsInternal.h"
8 
51 #ifdef MAKINGPUBLICDOCS
52 tntNetwork tntMpsCreateMpo(unsigned L,
53 #else
54 tntNetwork tntMpsCreateMpo_(unsigned L,
55 #endif
56  tntNodeArray *nnl,
57  tntNodeArray *nnr,
58  tntComplexArray *nnparam,
59  tntNodeArray *os,
60  tntComplexArray *osparam,
61  unsigned ignoreQN)
62 {
63  tntNetwork mponw; /* Network representing the MPO */
64  tntComplexArray mpogenvals; /* Values for one of the 'building blocks' of the MPO */
65  tntNode mpogen, mpobase, mpoc, eye, opss, opssterm, opnn; /* Nodes for building the MPO */
66  tntNode Lv, Rv; /* Left and right boundary vectors for the MPO */
67  unsigned Dmpo; /* Internal dimension of the MPO */
68  unsigned i, j; /* Used for looping over terms and sites respectively */
69  unsigned numnodes; /* The number of different nodes that need to be created */
70  tntComplex prm; /* The current paramter */
71  tntIntArray qnums_phys, qnums_int, qnums_op; /* Quantum numbers for the physical legs and the internal legs, and for the current operator */
72 
73  /* Check that there are the same number of left and right nearest neighbour operators */
74  if ((NULL == nnl && NULL != nnr) || (NULL != nnl && NULL == nnr) || ((NULL != nnl && NULL != nnr) && nnl->sz != nnr->sz))
75  tntErrorPrint("Cannot create a matrix product operator|The number of left nearest terms is not equal to the number of right nearest neighbour terms"); /* NO_COVERAGE */
76 
77  /* 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 */
78  if ((nnl != NULL) && (NULL == nnparam || nnl->sz != nnparam->numrows))
79  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"); /* NO_COVERAGE */
80 
81  /* Check that there is a parameter for each of the onsite operators */
82  if ((os != NULL) && (NULL == osparam || os->sz != osparam->numrows))
83  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"); /* NO_COVERAGE */
84 
85  /* Check that there is a non-zero number of operators */
86  if ((NULL == nnl || 0 == nnl->sz) && (NULL == os || 0 == os->sz)) tntErrorPrint("Cannot create a matrix product operator, as there are no terms"); /* NO_COVERAGE */
87 
88  /* Check that the length of the system is more than 1 */
89  if (L < 2) tntErrorPrint("Function to create MPO is only valid for a system length of two or more"); /* NO_COVERAGE */
90 
91  /* Calulcate how many nodes to create: one node if the operators are uniform, L nodes if any of the operators vary */
92  if ((NULL == os || 0 == os->sz || 1 == osparam->numcols) && (NULL == nnl || 0 == nnl->sz || 1 == nnparam->numcols)) numnodes = 1;
93  else numnodes = L;
94 
95  /* Create the empty network that will contain all the MPOs and that will be returned by the function. */
96  mponw = tntNetworkCreate();
97 
98  /* Create the identity node basing it on the supplied single-site operators */
99  if (NULL != nnl && 0 != nnl->sz) eye = tntNodeCreateEyeOp(nnl->vals[0]);
100  else eye = tntNodeCreateEyeOp(os->vals[0]);
101 
102  /* The internal dimension of the MPO will be the number of nn terms + 2. */
103  Dmpo = (NULL == nnl) ? 2 : nnl->sz + 2;
104 
105  /* ----------------- SETTING QN INFO HERE ---------------- */
106  /* 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 */
107  if (tntSymmTypeGet() && (0 == ignoreQN)) {
108  /* Turn off warning since we know we are going to strip QN */
110 
111  /* Get the quantum numbers for the physical leg */
112  qnums_phys = tntNodeGetQN(tntSysBasisOpGet(), "D");
113 
114  /* Allocate array of the correct size for the internal leg quantum numbers - all the values will be initialised to zero */
115  qnums_int = tntIntArrayAlloc(Dmpo*tntSymmNumGet());
116 
117  /* Loop through the nearest neighbour operators, finding the correct internal quantum number for each one */
118  if (nnl != NULL) {
119  for (i = 0; i < nnl->sz; i++) {
120  /* Get quantum number for left operator since the function calculates QN for outgoing leg, and the left operator will be associated with the outgoing leg */
121  qnums_op = tntMpsOpGetQN(nnl->vals[i]);
122 
123  /* copy it to the relevant entry for the internal leg */
124  for (j = 0; j < tntSymmNumGet(); j++) qnums_int.vals[(i+1)*tntSymmNumGet() + j] = qnums_op.vals[j];
125 
126  /* free the array containing the quantum number label */
127  tntIntArrayFree(&qnums_op);
128  }
129  }
130 
131  /* Find the quantum number for the first index in the internal leg */
132  /* If there is an onsite term it will be equal to the quantum number for this */
133  /* Otherwise the resultant qn of the nearest neighbour operators needs to be found */
134  if (os != NULL && 0 != os->sz) {
135  qnums_op = tntMpsOpGetQN(os->vals[0]);
136  /* Copy the quantum numbers to the first position in the array */
137  for (j = 0; j < tntSymmNumGet(); j++) qnums_int.vals[j] = qnums_op.vals[j];
138 
139  tntIntArrayFree(&qnums_op);
140  } else {
141  tntIntArray qnums_op_second;
142 
143  qnums_op = tntMpsOpGetQN(nnl->vals[0]);
144  qnums_op_second = tntMpsOpGetQN(nnr->vals[0]);
145 
146  /* TODO: Need to make this more general when other symmetry types are added */
147  /* Add the quantum numbers together */
148  for (j = 0; j < tntSymmNumGet(); j++) qnums_int.vals[j] = qnums_op.vals[j] + qnums_op_second.vals[j];
149 
150  tntIntArrayFree(&qnums_op);
151  tntIntArrayFree(&qnums_op_second);
152  }
153 
154 
155  }
156  /* ----------------- END OF SETTING QN INFO ----------------- */
157 
158  /* Create the generator for the MPO that doesn't depend on operators i.e. identity in first and last element */
159  mpogenvals = tntComplexArrayAlloc(Dmpo*Dmpo);
160  mpogenvals.vals[0].re = mpogenvals.vals[Dmpo*Dmpo - 1].re = 1.0;
161  mpogen = tntNodeCreate(&mpogenvals, "LR", Dmpo, Dmpo);
162 
163  /* Contract the gen with the identity operator to form the base mpo term */
164  /* Note that the contraction step will free both mpogen and eye */
165  mpobase = tntNodeContract(mpogen, eye);
166 
167  /* reset the numbers for the generator values */
168  mpogenvals.vals[0].re = mpogenvals.vals[Dmpo*Dmpo - 1].re = 0.0;
169 
170  /* 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
171  if any one of the operators vary through the system */
172  /* Insert first node in the network */
173  tntNodeInsertAtEnd(mpobase, "L", "R", mponw);
174 
175  /* Insert any remaining nodes in the network */
176  for (j = 1; j < numnodes; j++) tntNodeInsertAtEnd(tntNodeCopy(mpobase), "L", "R", mponw);
177 
178  /* Now create the on-site terms: Create a different term for each site if the parameters vary. */
179  if (os != NULL && 0 != os->sz) {
180 
181  /* Populate the bottom-left element of the MPO generator and create node, then reset array */
182  mpogenvals.vals[Dmpo-1].re = 1.0;
183  mpogen = tntNodeCreate(&mpogenvals, "LR", Dmpo, Dmpo);
184  mpogenvals.vals[Dmpo-1].re = 0.0;
185 
186  for (j = 0; j < numnodes; j++) {
187 
188  /* Let the onsite term be equal to the first operator */
189  opss = tntNodeCopy(os->vals[0]);
190 
191  /* Set the parameters depending on whether the on-site parameters also vary */
192  if (1 == osparam->numcols) prm = osparam->vals[0];
193  else prm = osparam->vals[j * os->sz];
194 
195  /* Scale the operator by the parameter */
196  tntNodeScaleComplex(opss, prm);
197 
198  /* Create each addition term, and add it to the total on-site operator */
199  for (i = 1; i < os->sz; i++) {
200  /* Copy operator */
201  opssterm = tntNodeCopy(os->vals[i]);
202 
203  /* Set the parameters depending on whether the on-site parameters also vary */
204  if (1 == osparam->numcols) prm = osparam->vals[i];
205  else prm = osparam->vals[i + j * os->sz];
206 
207  /* Determine the correct parameter and scale */
208  tntNodeScaleComplex(opssterm, prm);
209 
210  /* Add to the total term */
211  tntNodeAdd(opss,opssterm);
212 
213  /* Free as it will no longer be required */
214  tntNodeFree(&opssterm);
215  }
216 
217  /* Create a copy of the MPO generator */
218  mpoc = tntNodeCopy(mpogen);
219 
220  /* Contract the MPO generator with the on-site operator */
221  mpoc = tntNodeContract(mpoc, opss);
222 
223  /* Add the MPO term to the total MPO */
224  tntNodeAdd(mpobase, mpoc);
225 
226  /* free this generating node as it is no longer required */
227  tntNodeFree(&mpoc);
228 
229  /* Find the next mpobase term */
230  mpobase = tntNodeFindConn(mpobase, "R");
231 
232  }
233  /* Free the generating node */
234  tntNodeFree(&mpogen);
235  }
236 
237  /* Now create the nearest neighbour terms */
238  if (nnl != NULL && 0 != nnl->sz) {
239  /* Create each term, using the MPO generator to insert each nn term to a different position in the MPO tensor. */
240  for (i = 0; i < nnl->sz; i++) {
241 
242  /* Make the MPO generator: position of left operator will be bottom row, and (i + 1)st column (if col numbering starts from zero) */
243  mpogenvals.vals[Dmpo-1 + (i+1)*Dmpo].re = 1.0;
244  mpogen = tntNodeCreate(&mpogenvals, "LR", Dmpo, Dmpo);
245 
246  /* Reset the array for the MPO generator */
247  mpogenvals.vals[Dmpo-1 + (i+1)*Dmpo].re = 0.0;
248 
249  /* Find the first MPO in the network */
250  mpobase = tntNodeFindFirst(mponw);
251 
252  /* Now loop through sites for the left operator of the pair.
253  Note there are L MPOs, but only L-1 parameters for the case of
254  varying parameters. Put the parameter on the L node, which although built for the last site, will be
255  destroyed by the boundary node. Therefore for the last node do not scale */
256 
257  for (j = 0; j < numnodes; j++) {
258  /* Copy the left operator and scale by parameter */
259  opnn = tntNodeCopy(nnl->vals[i]);
260 
261  /* Determine the correct parameter */
262  if (j == L-1) {
263  prm.re = 1.0;
264  prm.im = 0.0;
265  } else {
266  /* Set the parameters depending on whether the nearest neighbour parameters also vary */
267  if (1 == nnparam->numcols) prm = nnparam->vals[i];
268  else prm = nnparam->vals[i + j*nnl->sz];
269  }
270 
271  tntNodeScaleComplex(opnn, prm);
272 
273  /* Make a copy of the MPO generator */
274  mpoc = tntNodeCopy(mpogen);
275 
276  /* Contract the MPO generator and the single site operator */
277  mpoc = tntNodeContract(mpoc, opnn);
278 
279  /* Add the MPO term to the total MPO */
280  tntNodeAdd(mpobase, mpoc);
281 
282  /* free this generating node as it is no longer required */
283  tntNodeFree(&mpoc);
284 
285  /* Find the next MPO in the network */
286  mpobase = tntNodeFindConn(mpobase,"R");
287 
288  }
289 
290  /* free this generating node as it is no longer required */
291  tntNodeFree(&mpogen);
292 
293  /* Make the MPO generator: position of right operator will be first column, and (i + 1)st row (if row numbering starts from zero) */
294  mpogenvals.vals[i+1].re = 1.0;
295  mpogen = tntNodeCreate(&mpogenvals, "LR", Dmpo, Dmpo);
296 
297  /* Reset the array for the MPO generator */
298  mpogenvals.vals[i+1].re = 0.0;
299 
300  /* Find the first MPO in the network */
301  mpobase = tntNodeFindFirst(mponw);
302 
303  /* Now loop through the sites for the right node of the pair. These
304  do not need to be scaled */
305  for (j = 0; j < numnodes; j++) {
306 
307  /* Copy the right operator (this does not need to be scaled by the parameter) */
308  opnn = tntNodeCopy(nnr->vals[i]);
309 
310  /* Make a copy of the MPO generator */
311  mpoc = tntNodeCopy(mpogen);
312 
313  /* Contract the MPO generator and the single site operator */
314  mpoc = tntNodeContract(mpoc, opnn);
315 
316  /* Add the MPO term to the total MPO */
317  tntNodeAdd(mpobase, mpoc);
318 
319  /* free this generating node as it is no longer required */
320  tntNodeFree(&mpoc);
321 
322  /* Find the next MPO in the network */
323  mpobase = tntNodeFindConn(mpobase,"R");
324 
325  }
326 
327  /* free this generating node as it is no longer required */
328  tntNodeFree(&mpogen);
329  }
330  }
331 
332  /* ----------------- SETTING QN INFO HERE ---------------- */
333  /* Now put all the operators in the network in blocks form */
334  if (tntSymmTypeGet() && (0 == ignoreQN)) {
335 
336  mpoc = tntNodeFindFirst(mponw);
337 
338  /* Set the quantum numbers for the first MPO. */
339  tntNodeSetQN(mpoc,"D",&qnums_phys,TNT_QN_IN);
340  tntNodeSetQN(mpoc,"U",&qnums_phys,TNT_QN_OUT);
341  tntNodeSetQN(mpoc,"L",&qnums_int,TNT_QN_IN);
342  tntNodeSetQN(mpoc,"R",&qnums_int,TNT_QN_OUT);
343 
344  while (mpoc != tntNodeFindLast(mponw)) {
345 
346  /* Move on the the next node */
347  mpoc = tntNodeFindConn(mpoc,"R");
348 
349  /* Set the quantum numbers for the MPO. */
350  tntNodeSetQN(mpoc,"D",&qnums_phys,TNT_QN_IN);
351  tntNodeSetQN(mpoc,"U",&qnums_phys,TNT_QN_OUT);
352  tntNodeSetQN(mpoc,"L",&qnums_int,TNT_QN_IN);
353  tntNodeSetQN(mpoc,"R",&qnums_int,TNT_QN_OUT);
354  }
355 
356 
357  /* free the quantum numbers */
358  tntIntArrayFree(&qnums_phys);
359  }
360  /* ----------------- END OF SETTING QN INFO ----------------- */
361 
362  /* Now if only one node has been created, create additional L-1 identical copies to put in the network */
363  if (1 == numnodes) {
364  mpobase = tntNodeFindFirst(mponw);
365  /* Keep inserting copies of the mpo nodes at the end. */
366  for (j = 1; j < L; j++) tntNodeInsertAtEnd(tntNodeCopy(mpobase), "L", "R", mponw);
367  }
368 
369  /* Create the leftmost and rightmost MPO nodes by contracting with the boundary nodes */
370  /* Elements for the boundary nodes - reallocate array of values */
371  tntComplexArrayFree(&mpogenvals);
372  mpogenvals = tntComplexArrayAlloc(Dmpo);
373 
374  /* Left boundary vector: Want to pick out the final row, so choose row vector with last element 1 */
375  mpogenvals.vals[Dmpo-1].re = 1.0;
376  Lv = tntNodeCreate(&mpogenvals,"RL", Dmpo);
377 
378  /* Right boundary vector: Want to pick out the first row, so choose column vector with the first element 1 */
379  mpogenvals.vals[Dmpo-1].re = 0.0;
380  mpogenvals.vals[0].re = 1.0;
381  Rv = tntNodeCreate(&mpogenvals,"LR", Dmpo);
382 
383  /* ----------------- SETTING QN INFO HERE ---------------- */
384  /* Now put the boundary vectors in blocks form */
385  if (tntSymmTypeGet() && (0 == ignoreQN)) {
386 
387  mpoc = tntNodeFindFirst(mponw);
388 
389  /* Set the quantum numbers for the right boundary vector. */
390  tntNodeSetQN(Rv,"L",&qnums_int,TNT_QN_IN);
391  tntNodeSetQN(Rv,"R",&qnums_int,TNT_QN_OUT);
392 
393  /* Set the outgoing quantum numbers for the left boundary vector. */
394  tntNodeSetQN(Lv,"R",&qnums_int,TNT_QN_OUT);
395 
396  /* Incoming boundary vectors are always zero */
397  for (j = 0; j < tntSymmNumGet(); j++) qnums_int.vals[j] = 0;
398  tntNodeSetQN(Lv,"L",&qnums_int,TNT_QN_IN);
399 
400  /* Now quantum numbers for internal legs are finished with - free the array */
401  tntIntArrayFree(&qnums_int);
402 
403  /* Turn the warning back on */
405  }
406  /* ----------------- END OF SETTING QN INFO ----------------- */
407 
408  /* Insert the left and right boundary vectors in the network */
409  tntNodeInsertAtStart(Lv, "L", "R", mponw);
410  tntNodeInsertAtEnd(Rv, "L", "R", mponw);
411 
412  /* Contract left and right boundary vectors their neighbouring MPOs */
413  mpoc = tntNodeFindConn(Lv, "R");
414  Lv = tntNodeContract(Lv, mpoc);
415 
416  mpoc = tntNodeFindConn(Rv, "L");
417  Rv = tntNodeContract(Rv, mpoc);
418 
419  /* Free the arrays and nodes that are no longer required */
420  tntComplexArrayFree(&mpogenvals);
421 
422  /* Return the MPO network */
423  return mponw;
424 }
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
void tntNodeSetQN(tntNode A, tntLegLabel legA, tntIntArray *qvals, int legdir)
Definition: tntNodeQN.c:37
tntNetwork tntNetworkCreate(void)
Definition: tntNetwork.c:92
tntNode tntNodeFindLast(tntNetwork nw)
Definition: tntNodeInfo.c:68
unsigned tntSymmNumGet(void)
Definition: tntSys.c:512
tntNode tntNodeCreateEyeOp(tntNode A)
Definition: tntNodeUtil.c:195
void tntNodeInsertAtEnd(tntNode I, tntLegLabel legIlast, tntLegLabel legInwend, tntNetwork nw)
Definition: tntNodeConn.c:252
tntNetwork tntMpsCreateMpo(unsigned L, tntNodeArray *nnl, tntNodeArray *nnr, tntComplexArray *nnparam, tntNodeArray *os, tntComplexArray *osparam, unsigned ignoreQN)
void tntNodeScaleComplex(tntNode A, tntComplex val)
void tntNodeAdd(tntNode A, tntNode B)
void tntNodeFree(tntNode *A)
Definition: tntNodeUtil.c:275
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:392
tntIntArray tntNodeGetQN(tntNode A, tntLegLabel legA)
Definition: tntNodeQN.c:155
Definition: tnt.h:65
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray.c:620
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:29
tntIntArray tntMpsOpGetQN(tntNode op)
Definition: tntMpsUtil.c:51
tntNode tntNodeFindConn(tntNode A, tntLegLabel legA)
Definition: tntNodeInfo.c:29
tntNode tntNodeContract(tntNode A, tntNode B, tntLegLabel legMapAC, tntLegLabel legMapBC)
tntNode tntSysBasisOpGet(void)
Definition: tntSys.c:407
tntNode tntNodeCreate(tntComplexArray *nodeVals, tntLegLabel leglabels, unsigned dimleg1, unsigned dim...)
Definition: tntNodeUtil.c:354
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray.c:650
double re
Definition: tnt.h:66
tntNode tntNodeFindFirst(tntNetwork nw)
Definition: tntNodeInfo.c:50
void tntNodeInsertAtStart(tntNode I, tntLegLabel legInwstart, tntLegLabel legIfirst, tntNetwork nw)
Definition: tntNodeConn.c:213
int tntSymmTypeGet(void)
Definition: tntSys.c:498
void tntSysQNClearWarnOn(void)
Definition: tntUtil.c:17
double im
Definition: tnt.h:67
void tntSysQNClearWarnOff(void)
Definition: tntUtil.c:35