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
tntMpsMps.c
1 /*
2 Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
3 $LastChangedDate: 2016-07-29 10:32:20 +0100 (Fri, 29 Jul 2016) $
4 (c) University of Oxford 2014
5 */
6 
12 /* Include the header for the TNT MPS library */
13 #include "tntMpsInternal.h"
14 
15 
38 void tntMpsMpsConnect(tntNetwork mpsA,
41  tntNetwork mpsB)
42 {
43 
44  tntNode A, An; /* The current node and the next node of the unflipped MPS */
45  tntNode B, Bn; /* The current node and the next node of the MPS to flip */
46  tntNode Bf, Bfn; /* Flipped copy of the current node and the next node in mpsB */
47 
48  /* Before performing the flip operation, check that the MPS networks have the same number of sites */
49  if (tntMpsLength(mpsA) != tntMpsLength(mpsB)) {
50  tntErrorPrint("The MPS networks cannot be connected as they do not have the same number of sites"); /* NO_COVERAGE */
51  } /* NO_COVERAGE */
52 
53  /* Find the first node in wfa */
54  A = tntNodeFindFirst(mpsA);
55 
56  /* Find the first node in wfb */
57  B = tntNodeFindFirst(mpsB);
58 
59  /* Make a new flipped node
60  - First copy the node in wfb, specifying that the complex conjugate should be taken.
61  - Then change the id of a downwards facing physical leg to an upwards facing physical leg
62  */
63  Bf = tntNodeCopy(B, 1);
64  tntNodeMapLegs(Bf, "D=U");
65 
66  /* Join the flipped and unflipped node on the physical leg */
67  tntNodeJoin(A,"D",Bf,"U");
68 
69 
70  /* Move through the network adding a flipped node to each physical leg of the original wave function */
71  while (A != tntNodeFindLast(mpsA)) {
72 
73  /* identify the next node along in wfa and wfb */
74  An = tntNodeFindConn(A, "R");
75  Bn = tntNodeFindConn(B, "R");
76 
77  /* Make a new flipped node */
78  Bfn = tntNodeCopy(Bn, 1);
79  tntNodeMapLegs(Bfn, "D=U");
80 
81  /* Join the original flipped node and the new flipped node */
82  tntNodeJoin(Bf, "R", Bfn, "L");
83 
84  /* Join the new flipped node with its unflipped node on their physical legs */
85  tntNodeJoin(An,"D",Bfn,"U");
86 
87  /* Move along one site */
88  A = An;
89  B = Bn;
90  Bf = Bfn;
91  }
92 
93  return;
94 
95 }
96 
149 #ifdef MAKINGPUBLICDOCS
151 #else
152 tntComplex tntMpsMpsContract_(
153 #endif
154  tntNetwork sandwich,
155  int site_left,
156  int site_right,
157  tntNode beta,
158  tntNode gamma)
160 {
161 
162  /* Networks and nodes during the contraction */
163  tntNetwork mps_disc; /* The part of the network that will be equated to the identity and discarded during the contraction */
164  tntNode A, Af; /* The current MPS node, and flipped counterpart */
165  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. */
166  tntNode alpha; /* The large node that is a result of contracting the MPS */
167 
168  /* Counters and site information used during the contraction */
169  int sitecount = 0; /* the current site count */
170  int L; /* The number of MPS nodes in the network */
171 
172  /* Final results */
173  tntComplex contracted; /* The contracted network */
174 
175  /* Determine the length of the MPS. */
176  L = (int) tntMpsLength(sandwich);
177 
178  /* Deal with site_left and site_right */
179  if ((-1 == site_left) || (-1 == site_right)) {
180  /* set site_left and site_right to the beginning and end of the network if a full contraction is indicated */
181  site_left = 0;
182  site_right = L -1;
183  } else {
184  /* otherwise check that the values given are sensible, and modify network appropriately */
185 
186  /* Check site_left is not more than site_right */
187  if (site_left > site_right) tntErrorPrint("site_left (%d) is greater than site_right (%d)",site_left, site_right); /* NO_COVERAGE */
188 
189  /* Check site_right is not too big */
190  if (site_right >= L) tntErrorPrint("site_right (%d) is more than or equal to the length of the MPS (%d)", site_right, L); /* NO_COVERAGE */
191 
192  /* If the leftmost site required is more than zero, discard left parts of the MPS network */
193  if (site_left > 0) {
194 
195  /* First move to the correct site - the last of the sites to be discarded */
196  A = tntNodeFindFirst(sandwich);
197  for (sitecount = 1; sitecount < site_left; sitecount++) A = tntNodeFindConn(A, "R");
198  Af = tntNodeFindConn(A, "D");
199 
200  /* A and Af are the nodes that need to be split that will be the last nodes of the left network (to be discarded). */
201  nodes_L[0] = A;
202  nodes_L[1] = Af;
203 
204  /* 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 */
205  nodes_R[0] = A = tntNodeFindConn(A, "R");
206 
207  /* Find the corresponding flipped node. This will also belong to the right network */
208  nodes_R[1] = Af = tntNodeFindConn(A, "D");
209 
210  /* Split the network into two */
211  mps_disc = tntNetworkSplit(sandwich, /* network to split */
212  nodes_R[0], "s", tntNodeFindLast(sandwich), "R", /* info for first network (right part, keeping this) */
213  tntNodeFindFirst(sandwich), "s", nodes_L[0], "e", /* info for second network (left part, discarding this this) */
214  2, nodes_R, 2, nodes_L); /* info about nodes on split */
215 
216  /* Discard the right network */
217  tntNetworkFree(&mps_disc);
218 
219  if (NULL == beta) {
220  /* Join the left internal legs of the flipped and unflipped MPS node on the leftmost site */
221  tntNodeJoin(A,"L",Af,"L");
222  } else {
223  /* Join the left internal legs to the beta node */
224  beta = tntNodeCopy(beta);
225  tntNodeJoin(A,"L",beta,"R");
226  tntNodeJoin(Af,"L",beta,"S");
227  }
228  }
229 
230  /* If the right site required is less than the last site, discard right parts of the MPS Network */
231  if (site_right < L-1) {
232 
233  /* First move to the correct site - the first of the sites to be discarded */
234  A = tntNodeFindLast(sandwich);
235  for (sitecount = L-2; sitecount > site_right; sitecount--) A = tntNodeFindConn(A, "L");
236  Af = tntNodeFindConn(A, "D");
237 
238  /* A and Af are the nodes that need to be split that will be part of the right network (to be discarded) */
239  nodes_R[0] = A;
240  nodes_R[1] = Af;
241 
242  /* Find the node before A. This will now belong to the left network (to be kept). */
243  nodes_L[0] = A = tntNodeFindConn(A, "L");
244 
245  /* Find the corresponding flipped node. This will also belong to the right network */
246  nodes_L[1] = Af = tntNodeFindConn(A, "D");
247 
248  /* 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). */
249  mps_disc = tntNetworkSplit(sandwich, /* network to split */
250  tntNodeFindFirst(sandwich), "s", A, "e", /* info for first network (left part, keeping this) */
251  nodes_R[0], "s", tntNodeFindLast(sandwich), "e", /* info for second network (right part, discarding this this) */
252  2, nodes_L, 2, nodes_R); /* info about nodes on split */
253 
254  /* Discard the right network */
255  tntNetworkFree(&mps_disc);
256 
257  if (NULL == gamma) {
258  /* Join the right internal legs of the flipped and unflipped MPS node on the rightmost site */
259  tntNodeJoin(A,"R",Af,"R");
260  } else {
261  gamma = tntNodeCopy(gamma);
262  tntNodeJoin(A,"R",gamma,"L");
263  tntNodeJoin(Af,"R",gamma,"M");
264  }
265  }
266 
267  }
268 
269  /* Now that the unrequired parts of the network have been discarded, zip throught the network from left to right contracting nodes. */
270  sitecount = site_left;
271 
272  /* identify the first node */
273  A = tntNodeFindFirst(sandwich);
274  Af = tntNodeFindConn(A, "D");
275 
276  /* Form the first alpha tensor, using the precontracted node if it exists */
277  if (site_left > 0 && beta != NULL) {
278  beta = tntNodeContract(beta, A);
279  alpha = tntNodeContract(beta, Af, NULL, "LR=MS");
280  } else {
281 
282  alpha = tntNodeContract(A, Af, NULL, "LR=MS");
283  }
284 
285  /* Go through remaining sites, contracting them with alpha one by one */
286  for (sitecount = site_left+1; sitecount <= site_right; sitecount++) {
287 
288  /* Find the next nodes along */
289  A = tntNodeFindConn(alpha, "R");
290  Af = tntNodeFindConn(A, "D");
291 
292  /* Contract alpha with them */
293  alpha = tntNodeContract(alpha,A);
294  alpha = tntNodeContract(alpha,Af, NULL, "LR=MS");
295  }
296 
297  /* If a precontracted node has been inserted for the right part of the network, contract it with the alpha tensor */
298  if (site_right < L-1 && NULL != gamma) alpha = tntNodeContract(alpha, gamma, NULL, NULL);
299 
300  /* Get the one remaining element */
301  contracted = tntNodeGetFirstEl(alpha);
302 
303  /* Free the network */
304  tntNetworkFree(&sandwich);
305 
306  /* Return the result of the contraction */
307  return contracted;
308 
309 }
310 
329 void tntMpsSelfInit(tntNetwork mps,
330  tntNodeArray *betas,
331  tntNodeArray *gammas)
332 {
333  tntNetwork sandwich; /* The full MPS-MPS network */
334  tntNode A, Af; /* The current node A of the MPS and the flipped (hermitian conjugate) of the current node Af */
335  tntNode beta, gamma; /* The node that is the results of all previously contracted nodes from the right and left */
336  tntNode eye; /* Identity node having one element to act as terminating node */
337  tntComplexArray eyevals;/* Values for the identity node to act as termninating node */
338  int L, j; /* the length of the MPS network and variable used for looping */
339 
340 
341  /* Create the required identity node: it will have the same form as a beta or gamma node, but all legs will have dimension 1, and the value of it will be 1. */
342  eyevals = tntComplexArrayAlloc(1);
343  eyevals.vals[0].re = 1.0;
344  eye = tntNodeCreate(&eyevals, "LMRS");
345  /* Free the array values */
346  tntComplexArrayFree(&eyevals);
347 
348  /* Find the length of the network */
349  L = (int) tntMpsLength(mps);
350 
351  /* ------ First start contracting from the left to the right to form the nodes required for beta ---- */
352 
353  if (NULL != betas) {
354  *betas = tntNodeArrayAlloc(L);
355 
356  /* Create the full network */
357  sandwich = tntNetworkCopy(mps);
358  tntMpsMpsConnect(sandwich,mps);
359 
360  /* Find the MPS node, operator and flipped node corresponding to site 0 */
361  A = tntNodeFindFirst(sandwich);
362  Af = tntNodeFindConn(A, "D");
363 
364  /* Let beta be equal to the identity node */
365  beta = tntNodeCopy(eye);
366 
367  /* ----------------- SETTING QN INFO HERE ---------------- */
368  /* Put quantum number info (all zeros) on the legs of the node */
369  if (tntSymmTypeGet()) {
370  /* Set the quantum numbers for the identity node. */
371  tntNodeSetQN(beta,"L",NULL,TNT_QN_IN);
372  tntNodeSetQN(beta,"M",NULL,TNT_QN_IN);
373  tntNodeSetQN(beta,"R",NULL,TNT_QN_IN);
374  tntNodeSetQN(beta,"S",NULL,TNT_QN_IN);
375  }
376  /* ----------------- END OF SETTING QN INFO ----------------- */
377 
378  /* Insert a copy of the identity node at the start of the network and join to the other nodes */
379  tntNodeInsertAtStart(beta,"L","R",sandwich);
380  tntNodeJoin(beta,"S",Af,"L");
381 
382  /* Let beta be the first element in the array */
383  betas->vals[0] = tntNodeCopy(beta);
384 
385  /* Now loop through the network from left to right performing contractions */
386  for (j = 1; j < L; j++) {
387 
388  /* Contract beta with the MPS node, operator, and flipped MPS node corresponding to the next site */
389  beta = tntNodeListContract("LMRS", beta, A, Af);
390 
391  /* Copy the resulting node to the array */
392  betas->vals[j] = tntNodeCopy(beta);
393 
394  /* Find the three nodes that belong to the next site along */
395  A = tntNodeFindConn(beta, "R");
396  Af = tntNodeFindConn(beta, "S");
397 
398  }
399 
400  /* Free the network */
401  tntNetworkFree(&sandwich);
402  }
403 
404 
405  /* ------ Now contract from right to left to form the nodes required for gamma ---- */
406 
407  if (NULL != gammas) {
408  *gammas = tntNodeArrayAlloc(L);
409 
410  /* Create the full network */
411  sandwich = tntNetworkCopy(mps);
412  tntMpsMpsConnect(sandwich,mps);
413 
414  /* Find the MPS node, operator and flipped MPS node corresponding to site L-1 */
415  A = tntNodeFindLast(sandwich);
416  Af = tntNodeFindConn(A, "D");
417 
418  /* Let gamma be equal to the identity node */
419  gamma = tntNodeCopy(eye);
420 
421  /* ----------------- SETTING QN INFO HERE ---------------- */
422  /* Put quantum number info (all zeros) on the legs of the node if they are present on the MPO and MPS */
423  if (tntSymmTypeGet()) {
424 
425  tntIntArray legqn; /* Contains the quantum number for the relevant leg */
426 
427  /* quantum numbers for the first two legs come from the last MPS */
428  legqn = tntNodeGetQN(A,"R");
429  if (legqn.sz > 0) {
430  tntNodeSetQN(gamma,"R",&legqn,tntNodeGetLegDir(A,"R"));
431  tntNodeSetQN(gamma,"L",&legqn,-1*tntNodeGetLegDir(A,"R"));
432  tntIntArrayFree(&legqn);
433  }
434 
435  /* quantum numbers for the last two legs come from the last flipped MPS */
436  legqn = tntNodeGetQN(Af,"R");
437  if (legqn.sz > 0) {
438  tntNodeSetQN(gamma,"S",&legqn,tntNodeGetLegDir(Af,"R"));
439  tntNodeSetQN(gamma,"M",&legqn,-1*tntNodeGetLegDir(Af,"R"));
440  tntIntArrayFree(&legqn);
441  }
442  }
443  /* ----------------- END OF SETTING QN INFO ----------------- */
444 
445  /* Insert a copy of the identity node at the end of the network */
446  tntNodeInsertAtEnd(gamma,"L","R",sandwich); /* insert between terminating node and first MPS node */
447  tntNodeJoin(gamma,"M",Af,"R"); /* Join it to the flipped MPS node */
448 
449  /* Let gamma be the last element in the array */
450  gammas->vals[L-1] = tntNodeCopy(gamma);
451 
452  /* Now loop through the network from right to left performing contractions */
453  for (j = L-2; j >= 0; j--) {
454 
455  /* Contract beta with the MPS node, operator, and flipped MPS node corresponding to the next site left */
456  gamma = tntNodeListContract("RSLM", gamma, A, Af);
457 
458  /* Copy the resulting node to the array */
459  gammas->vals[j] = tntNodeCopy(gamma);
460 
461  /* Find the next three nodes belonging to the site to the left */
462  A = tntNodeFindConn(gamma, "L");
463  Af = tntNodeFindConn(gamma, "M");
464  }
465 
466  /* Free unrequired objects */
467  tntNetworkFree(&sandwich);
468  }
469 
470  tntNodeFree(&eye);
471 }
472 
490 tntComplex tntMpsMpsProduct(tntNetwork mpsA,
491  tntNetwork mpsB)
492 {
493 
494  tntNetwork sandwich; /* The network that will be contracted to find the product */
495 
496  /* `wfa` will be deleted during the network contraction, so use a copy so that the original wave function is untouched. */
497  sandwich = tntNetworkCopy(mpsA);
498 
499  /* First flip the network, and connect to mpsB (this function leave mpsB unchanged and connects a copy of mpsB) */
500  tntMpsMpsConnect(sandwich, mpsB);
501 
502  /* Now contract the nodes, performing a full contraction */
503  return tntMpsMpsContract(sandwich,-1,-1);
504 
505 }
506 
507 
544 double tntMpsSelfProduct(tntNetwork mps,
545  int orth_centre)
546 {
547 
548  /* Networks and nodes during the contraction */
549  tntNetwork mpsmps; /* The MPS-MPS network that will be contracted to find the expectation value */
550  tntComplex norm;/* The expectation value: real part is returned by the function. */
551  int L; /* Number of sites in the MPS */
552 
553  /* Find length of the MPS */
554  L = (int) tntMpsLength(mps);
555 
556  /* Check orth_centre is not too big */
557  if (orth_centre >= L) tntErrorPrint("The site for orth_centre (%d) is more than or equal to the length of the MPS (%d)", orth_centre, L); /* NO_COVERAGE */
558 
559  /* The wave function will be deleted during the network contraction, so use a copy so that the original wave function is untouched. */
560  mpsmps = tntNetworkCopy(mps);
561 
562  /* First flip the network, and connect to itself */
563  tntMpsMpsConnect(mpsmps, mps);
564 
565  /* Now contract the network - this destroys the mpsmps network */
566  norm = tntMpsMpsContract(mpsmps, orth_centre, orth_centre);
567 
568  /* return the real part of the norm (no need to check whether there is an imaginary part as inner product will always be real */
569  return (norm.re);
570 }
571 
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
void tntNodeSetQN(tntNode A, tntLegLabel legA, tntIntArray *qvals, int legdir)
Definition: tntNodeQN.c:37
tntComplex tntMpsMpsProduct(tntNetwork mpsA, tntNetwork mpsB)
Definition: tntMpsMps.c:490
void tntNetworkFree(tntNetwork *nwp)
Definition: tntNetwork.c:118
void tntNodeJoin(tntNode A, tntLegLabel legA, tntNode B, tntLegLabel legB)
Definition: tntNodeConn.c:52
tntNode tntNodeFindLast(tntNetwork nw)
Definition: tntNodeInfo.c:68
tntNodeArray tntNodeArrayAlloc(unsigned num_elems)
Definition: tntArray.c:428
void tntNodeMapLegs(tntNode A, tntLegLabel legmap)
Definition: tntNodeConn.c:125
void tntNodeInsertAtEnd(tntNode I, tntLegLabel legIlast, tntLegLabel legInwend, tntNetwork nw)
Definition: tntNodeConn.c:252
tntComplex tntMpsMpsContract(tntNetwork sandwich, int site_left, int site_right, tntNode beta, tntNode gamma)
Definition: tntMpsMps.c:150
void tntMpsMpsConnect(tntNetwork mpsA, tntNetwork mpsB)
Definition: tntMpsMps.c:38
tntNetwork tntNetworkSplit(tntNetwork nw, tntNode tnFirst_1, tntLegLabel leg_start_1, tntNode tnLast_1, tntLegLabel leg_end_1, tntNode tnFirst_2, tntLegLabel leg_start_2, tntNode tnLast_2, tntLegLabel leg_end_2, unsigned numSplit_1, tntNode *tnSplit_1, unsigned numSplit_2, tntNode *tnSplit_2)
Definition: tntNetwork.c:203
double tntMpsSelfProduct(tntNetwork mps, int orth_centre)
Definition: tntMpsMps.c:544
void tntNodeFree(tntNode *A)
Definition: tntNodeUtil.c:275
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:392
unsigned tntMpsLength(tntNetwork wf)
Definition: tntMpsUtil.c:17
tntIntArray tntNodeGetQN(tntNode A, tntLegLabel legA)
Definition: tntNodeQN.c:155
Definition: tnt.h:65
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray.c:620
tntNode tntNodeFindConn(tntNode A, tntLegLabel legA)
Definition: tntNodeInfo.c:29
tntNode tntNodeContract(tntNode A, tntNode B, tntLegLabel legMapAC, tntLegLabel legMapBC)
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
tntNetwork tntNetworkCopy(tntNetwork nw)
Definition: tntNetwork.c:56
void tntMpsSelfInit(tntNetwork mps, tntNodeArray *betas, tntNodeArray *gammas)
Definition: tntMpsMps.c:329
int tntSymmTypeGet(void)
Definition: tntSys.c:498
int tntNodeGetLegDir(tntNode A, tntLegLabel legA)
Definition: tntNodeInfo.c:295
tntComplex tntNodeGetFirstEl(tntNode A)
Definition: tntNodeInfo.c:86