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
tntMpoExpecOutput.c
1 /*
2 Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
3 $LastChangedDate: 2013-10-09 10:02:48 +0100 (Wed, 09 Oct 2013) $
4 (c) University of Oxford 2014
5 */
6 
7 /* Include the header for the TNT MPS library */
8 #include "tntMpsInternal.h"
9 
10 /* Internal function to avoid repeated code */
11 void _tnt_MpoExpecOutputOneOp(tntComplexArray *Exval, int printOutput, int saveOutput, char *savename, char *label, unsigned counter);
12 
13 
32 void tntMpoExpecOutput(tntNetwork rho,
33  tntMpsExOp *Op,
34  unsigned printOutput,
35  unsigned saveOutput,
36  char *savename,
37  unsigned counter)
38 {
39 
40  tntComplexArray Exval; /* array for holding expectation values */
41  double norm; /* Norm of the density matrix. 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 
47  if (!printOutput && !saveOutput) {
48  tntWarningPrint("No expectation values are being calculated as both output flags are set to off"); /* NO_COVERAGE */
49  return; /* NO_COVERAGE */
50  } /* NO_COVERAGE */
51 
52  /* Find the length of the MPO */
53  L = (int) tntMpoLength(rho);
54 
55  /* Temporily turn off warnings, since we are going to be performing all operations on a copy of the wave function */
57 
58  /* Calculate the norm */
59  norm = tntComplexToReal(tntMpoTrace(rho));
60 
61  if (printOutput) {
62  tntPrintf("\n~~~~~~~ Norm is % #6.6g ~~~~~~~~\n\n",norm);
63  }
64  if (saveOutput) {
65  if (counter) {
66  tntDoubleParamsUpdate(savename, counter, norm);
67  } else {
68  tntDoubleParamsSave(savename, norm);
69  }
70 
71  }
72 
73  /* ----- Onsite expectation values ----- */
74 
75  if (Op->ExOpOs.sz > 0) {
76 
77  /* Allocate arrays of size one for onsite operators */
78  op = tntNodeArrayAlloc(1);
79  sitenum = tntIntArrayAlloc(1);
80 
81  if (printOutput) {
82  tntPrintf("\n~~~~~~~ Onsite expectation values ~~~~~~~~\n");
83  }
84 
85  }
86  /* Find the single site expectation values of the calculated ground state */
87  for (loop = 0; loop < Op->ExOpOs.sz; loop++) {
88 
89  /* Allocate memory for storing single-site expectation values */
90  Exval = tntComplexArrayAlloc(L);
91 
92  /* Copy over the required operator */
93  op.vals[0] = tntNodeCopy(Op->ExOpOs.vals[loop],0);
94 
95  for (sitenum.vals[0] = 0; sitenum.vals[0] < L; (sitenum.vals[0])++) {
96  /* Find the expectation value on each site using a product MPO. */
97  Exval.vals[sitenum.vals[0]] = tntMpoPmpoTrace(rho, &op, &sitenum);
98  Exval.vals[sitenum.vals[0]].re /= norm;
99  Exval.vals[sitenum.vals[0]].im /= norm;
100  }
101 
102  /* output and free values */
103  _tnt_MpoExpecOutputOneOp(&Exval, printOutput, saveOutput, savename, Op->LbOpOs.vals[loop], counter);
104 
105  /* Free the copy of the operator */
106  tntNodeFree(&(op.vals[0]));
107  }
108  if (Op->ExOpOs.sz > 0) {
109 
110  /* free the arrays */
111  tntIntArrayFree(&sitenum);
112  tntNodeArrayFree(&op);
113 
114  if (printOutput) {
115  tntPrintf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
116  }
117  }
118 
119 
120  /* ---------- Nearest-neighbour expectation values ----------- */
121 
122  if (Op->ExOp2nn.sz > 0) {
123 
124  /* Allocate arrays of size two for two-site operator site numbers */
125  sitenum = tntIntArrayAlloc(2);
126 
127  if (printOutput) {
128  tntPrintf("\n~~~~~~~ Nearest-neighbour expectation values ~~~~~~~~\n");
129  }
130 
131  }
132  for (loop = 0; loop < Op->ExOp2nn.sz; loop+=2) {
133 
134  /* Allocate memory for storing two-site expectation values */
135  Exval = tntComplexArrayAlloc(L-1);
136 
137  /* Allocate arrays of size two for two-site operators */
138  op = tntNodeArrayAlloc(2);
139 
140  /* Copy over the required operator */
141  op.vals[0] = tntNodeCopy(Op->ExOp2nn.vals[loop],0);
142  op.vals[1] = tntNodeCopy(Op->ExOp2nn.vals[loop+1],0);
143 
144  for (sitenum.vals[0] = 0; sitenum.vals[0] < L-1; (sitenum.vals[0])++) {
145  /* assign value for second site */
146  sitenum.vals[1] = sitenum.vals[0]+1;
147  /* Find the expectation value on each site using a product MPO. */
148  Exval.vals[sitenum.vals[0]] = tntMpoPmpoTrace(rho, &op, &sitenum);
149  Exval.vals[sitenum.vals[0]].re /= norm;
150  Exval.vals[sitenum.vals[0]].im /= norm;
151  }
152 
153  /* Free the operators */
154  tntNodeArrayFree(&op);
155 
156  /* output and free values */
157  _tnt_MpoExpecOutputOneOp(&Exval, printOutput, saveOutput, savename, Op->LbOp2nn.vals[loop/2], counter);
158 
159  }
160  if (Op->ExOp2nn.sz > 0) {
161 
162  /* free the arrays */
163  tntIntArrayFree(&sitenum);
164 
165  if (printOutput) {
166  tntPrintf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
167  }
168  }
169 
170  /* ---------- Centre-site expectation values ----------- */
171 
172  if (Op->ExOp2cs.sz > 0) {
173 
174  /* Allocate arrays of size two for two-site operator site numbers */
175  sitenum = tntIntArrayAlloc(2);
176 
177  if (printOutput) {
178  tntPrintf("\n~~~~~~~ Centre-site expectation values ~~~~~~~~\n");
179  }
180  }
181  for (loop = 0; loop < Op->ExOp2cs.sz; loop+=2) {
182 
183  /* Allocate memory for storing two-site expectation values */
184  Exval = tntComplexArrayAlloc(L);
185 
186  /* Allocate first site that expectation value will be found on */
187  sitenum.vals[0] = L/2;
188 
189  /* Allocate arrays of size two for two-site operators */
190  op = tntNodeArrayAlloc(2);
191 
192  /* Copy over the required operator */
193  op.vals[0] = tntNodeCopy(Op->ExOp2cs.vals[loop],0);
194  op.vals[1] = tntNodeCopy(Op->ExOp2cs.vals[loop+1],0);
195 
196  for (sitenum.vals[1] = 0; sitenum.vals[1] < L; (sitenum.vals[1])++) {
197  /* Find the expectation value on each site using a product MPO. */
198  Exval.vals[sitenum.vals[1]] = tntMpoPmpoTrace(rho, &op, &sitenum);
199  Exval.vals[sitenum.vals[1]].re /= norm;
200  Exval.vals[sitenum.vals[1]].im /= norm;
201  }
202 
203  /* Free the operators */
204  tntNodeArrayFree(&op);
205 
206  /* output and free values */
207  _tnt_MpoExpecOutputOneOp(&Exval, printOutput, saveOutput, savename, Op->LbOp2cs.vals[loop/2], counter);
208 
209  }
210  if (Op->ExOp2cs.sz > 0) {
211 
212  /* free the arrays */
213  tntIntArrayFree(&sitenum);
214 
215  if (printOutput) {
216  tntPrintf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
217  }
218  }
219 
220 
221  /* ---------- All pairs of two-site expectation values ----------- */
222 
223  if (Op->ExOp2ap.sz > 0) {
224 
225  /* Allocate arrays of size two for sitenumbers for two-site operators */
226  sitenum = tntIntArrayAlloc(2);
227 
228  if (printOutput) {
229  tntPrintf("\n~~~~~~~~~ All pairs of two-site expectation values ~~~~~~~~~\n");
230  }
231 
232  }
233  for (loop = 0; loop < Op->ExOp2ap.sz; loop+=2) {
234 
235  int ind; /* total matrix index */
236 
237  /* Allocate memory for storing two-site expectation values */
238  Exval = tntComplexArrayAlloc(L,L);
239 
240  /* Allocate arrays of size two for two-site operators */
241  op = tntNodeArrayAlloc(2);
242 
243  /* Copy over the required operator */
244  op.vals[0] = tntNodeCopy(Op->ExOp2ap.vals[loop],0);
245  op.vals[1] = tntNodeCopy(Op->ExOp2ap.vals[loop+1],0);
246 
247  /* Find the expectation value on each pair of sites using a product MPO. */
248  for (sitenum.vals[0] = 0; sitenum.vals[0] < L; (sitenum.vals[0])++) {
249  for (sitenum.vals[1] = 0; sitenum.vals[1] < L; (sitenum.vals[1])++) {
250  ind = sitenum.vals[0] + L*sitenum.vals[1];
251  Exval.vals[ind] = tntMpoPmpoTrace(rho, &op, &sitenum);
252  Exval.vals[ind].re /= norm;
253  Exval.vals[ind].im /= norm;
254  }
255  }
256 
257  /* Free the operators */
258  tntNodeArrayFree(&op);
259 
260  /* output and free values */
261  _tnt_MpoExpecOutputOneOp(&Exval, printOutput, saveOutput, savename, Op->LbOp2ap.vals[loop/2], counter);
262 
263  }
264  if (Op->ExOp2ap.sz > 0) {
265 
266  /* free the arrays */
267  tntIntArrayFree(&sitenum);
268 
269  if (printOutput) {
270  tntPrintf("\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n");
271  }
272 
273  }
274 
275  /* Turn the warning back on */
277 
278 }
279 
280 /* prints, saves and frees the complex array, converting to double if necessary, and taking account of flags and counters */
281 void _tnt_MpoExpecOutputOneOp(tntComplexArray *Exval, int printOutput, int saveOutput, char *savename, char *label, unsigned counter) {
282 
283  tntDoubleArray Exvalr; /* array for holding real part of expectation values */
284 
285  /* Now take different action depending on whether the expectation values are real or complex */
286  if (tntComplexArrayIsReal(Exval)) {
287  /* change complex array to real array */
288  Exvalr = tntComplexArrayToReal(Exval);
289  /* Print real output */
290  if (printOutput) {
291  tntNamedDoubleArrayPrint(label,&Exvalr);
292  }
293  /* Save this expectation value */
294  if (saveOutput) {
295  if (counter) tntDoubleArraysNamedUpdate(savename, counter,Exvalr, label);
296  else tntDoubleArraysNamedSave(savename, Exvalr, label);
297  }
298 
299  /* free real array */
300  tntDoubleArrayFree(&Exvalr);
301  } else {
302  /* Print complex output */
303  if (printOutput) {
304  tntNamedComplexArrayPrint(label,Exval);
305  }
306  /* Save complex output */
307  if (saveOutput) {
308  if (counter) tntComplexArraysNamedUpdate(savename, counter, *Exval, label);
309  else tntComplexArraysNamedSave(savename, *Exval, label);
310  }
311  /* Free complex array */
312  tntComplexArrayFree(Exval);
313  }
314 }
315 
double tntComplexToReal(tntComplex var)
Definition: tntUtil.c:108
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
#define tntMpsExOp
Definition: tntMps.h:75
void tntNodeArrayFree(tntNodeArray *arr)
Definition: tntArray.c:668
int tntComplexArrayIsReal(tntComplexArray *arr)
Definition: tntArray.c:570
void tntDoubleParamsSave(const char *fname, double varA, double varB, double varC, double varD, double varE)
Definition: tntSave.c:73
tntComplex tntMpoTrace(tntNetwork mpo)
Definition: tntMpoTrace.c:20
tntNodeArray tntNodeArrayAlloc(unsigned num_elems)
Definition: tntArray.c:428
tntDoubleArray tntComplexArrayToReal(tntComplexArray *arr)
Definition: tntArray.c:591
void tntDoubleArraysNamedUpdate(const char *fname, unsigned posn, tntDoubleArray arrA, const char *varnameA, tntDoubleArray arrB, const char *varnameB, tntDoubleArray arrC, const char *varnameC, tntDoubleArray arrD, const char *varnameD, tntDoubleArray arrE, const char *varnameE)
Definition: tntUpdate.c:302
void tntDoubleParamsUpdate(const char *fname, unsigned posn, double varA, double varB, double varC, double varD, double varE)
Definition: tntUpdate.c:88
void tntNodeFree(tntNode *A)
Definition: tntNodeUtil.c:275
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:392
void tntMpoExpecOutput(tntNetwork rho, tntMpsExOp *Op, unsigned printOutput, unsigned saveOutput, char *savename, unsigned counter)
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray.c:620
void tntDoubleArrayFree(tntDoubleArray *arr)
Definition: tntArray.c:635
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:29
void tntComplexArraysNamedSave(const char *fname, tntComplexArray arrA, const char *varnameA, tntComplexArray arrB, const char *varnameB, tntComplexArray arrC, const char *varnameC, tntComplexArray arrD, const char *varnameD, tntComplexArray arrE, const char *varnameE)
Definition: tntSave.c:327
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray.c:650
int tntPrintf(const char *format,...)
Definition: tntPrint.c:77
void tntSysQNClearWarnOn(void)
Definition: tntUtil.c:17
void tntSysQNClearWarnOff(void)
Definition: tntUtil.c:35
tntComplex tntMpoPmpoTrace(tntNetwork mpo, tntNodeArray *op, tntIntArray *sitenum)
Definition: tntMpoPmpo.c:35
void tntDoubleArraysNamedSave(const char *fname, tntDoubleArray arrA, const char *varnameA, tntDoubleArray arrB, const char *varnameB, tntDoubleArray arrC, const char *varnameC, tntDoubleArray arrD, const char *varnameD, tntDoubleArray arrE, const char *varnameE)
Definition: tntSave.c:271
void tntComplexArraysNamedUpdate(const char *fname, unsigned posn, tntComplexArray arrA, const char *varnameA, tntComplexArray arrB, const char *varnameB, tntComplexArray arrC, const char *varnameC, tntComplexArray arrD, const char *varnameD, tntComplexArray arrE, const char *varnameE)
Definition: tntUpdate.c:371