Tensor Network Theory Library  Beta release 1.2.1 A library of routines for performing TNT-based operations
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