Tensor Network Theory Library  Beta release 1.2.1 A library of routines for performing TNT-based operations
tntMpsCreatePropagator.c
1 /*
2 Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
3 $LastChangedDate$
4 (c) University of Oxford 2014
5 */
6
14 /* Include internal header */
15 #include "tntMpsInternal.h"
16
75 tntNodeArray tntMpsCreatePropArray(unsigned L,
76  tntComplex h,
77  tntNodeArray *nnl,
78  tntNodeArray *nnr,
79  tntComplexArray *nnparam,
80  tntNodeArray *os,
81  tntComplexArray *osparam)
82 {
83
84  tntNodeArray Parr; /* Array of propagators to be returned by function */
85  tntComplex hs; /* Uniform scale factor multiplied by -i for real part, i for imaginary part */
86  unsigned i, j; /* Used for looping over terms and sites respectively */
87  unsigned varparam=1; /* Flag to determine whether or not the parameters vary accross the system */
88  tntComplexArray two_site_params; /* Array for parameters on each site */
89
90  /* Check that there are the same number of left and right nearest neighbour operators */
91  if ((NULL == nnl && NULL != nnr) || (NULL != nnl && NULL == nnr) || ((NULL != nnl && NULL != nnr) && nnl->sz != nnr->sz))
92  tntErrorPrint("Cannot create a propagator!|The number of left nearest-neighbour terms is not equal to the number of right nearest-neighbour terms!"); /* NO_COVERAGE */
93
94  /* 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 */
95  if (((nnl != NULL)&&(nnl->sz != 0))&&(NULL == nnparam || nnl->sz != nnparam->numrows)) {
96  tntErrorPrint("Cannot create a propagator!|The number of nearest neighbour parameters (%d) is not equal to the number of nearest neighbour terms (%d)!",nnparam->numrows,nnl->sz);
97  } /* NO_COVERAGE */
98
99  /* Check that there is a parameter for each of the onsite operators */
100  if (((os != NULL)&&(os->sz != 0)) && (NULL == osparam || os->sz != osparam->numrows)) {
101  tntErrorPrint("Cannot create a propagator!|The number of on-site parameters per site (%d) is not equal to the number of on-site terms (%d)!",osparam->numrows,os->sz);
102  } /* NO_COVERAGE */
103
104  /* Check that the time step is not non-negligible */
105  if (TNT_DEFAULT_TOL > fabs(h.re) && TNT_DEFAULT_TOL > fabs(h.im))
106  tntErrorPrint("Error! Cannot create spin half two-site gates as the time step size is less than the tolerance."); /* NO_COVERAGE */
107
108  /* Check that there is a non-zero number of operators */
109  if ((NULL == nnl || 0 == nnl->sz) && (NULL == os || 0 == os->sz))
110  tntErrorPrint("Cannot create a propogator, as there are no terms!"); /* NO_COVERAGE */
111
112  /* Turn off warning since we know we are going to strip QN */
114
115  /* Allocate the array to be returned */
116  Parr = tntNodeArrayAlloc(L-1);
117
118  /* First take real part of time step and multiply by -i and take the imaginary part of the timestep and multiply by i */
119  hs.re = -h.im;
120  hs.im = -h.re;
121
122  /* Determine whether or not the parameters vary */
123  if ((NULL == os || 0 == os->sz || 1 == osparam->numcols) && (NULL == nnl || 0 == nnl->sz || 1 == nnparam->numcols)) varparam = 0;
124
125  /* If there are no on-site terms, and the parameters are uniform, then all the nodes will be identical */
126  if ((NULL == os || 0 == os->sz) && 0 == varparam) {
127
128  two_site_params = tntComplexArrayAlloc(nnl->sz); /* Array for parameters on each site */
129  /* Set the parameters for that site */
130  for (i = 0; i < nnl->sz; i++) two_site_params.vals[i] = nnparam->vals[i];
131
132  Parr.vals[0] = tntMpsCreateTwoSiteOp(nnl,nnr,&two_site_params);
133
134  /* Take the matrix exponential to form the two-site propagator */
135  tntNodeExp(Parr.vals[0], hs, "DE", "UV");
136  for (j = 1; j < L-1; j++) Parr.vals[j] = tntNodeCopy(Parr.vals[0]);
137
138  /* Free the array for parameters on each site */
139  tntComplexArrayFree(&two_site_params);
140
141  /* Return the array of nodes */
142  return Parr;
143  }
144
145  /* If there are no on-site terms, and the parameters vary, then set a different node for each two site gate */
146  if ((NULL == os || 0 == os->sz) && 1 == varparam) {
147  two_site_params = tntComplexArrayAlloc(nnl->sz); /* Array for parameters on each site */
148  for (j = 0; j < L-1; j++) {
149  /* Set the parameters for that site */
150  for (i = 0; i < nnl->sz; i++) {
151  two_site_params.vals[i].re = nnparam->vals[i + j*(nnl->sz)].re;
152  two_site_params.vals[i].im = nnparam->vals[i + j*(nnl->sz)].im;
153  }
154  /* Create the node for that site */
155  Parr.vals[j] = tntMpsCreateTwoSiteOp(nnl,nnr,&two_site_params);
156  /* Take the matrix exponential to form the two-site propagator */
157  tntNodeExp(Parr.vals[j], hs, "DE", "UV");
158  }
159  /* Free the array for parameters on each site */
160  tntComplexArrayFree(&two_site_params);
161
162  /* Return the array of nodes */
163  return Parr;
164  }
165
166  /* ---------------------------------------------------------------- */
167  /* This part of the function is only reached if on-site terms exist */
168  /* ---------------------------------------------------------------- */
169
170  /* If the are on-site terms, then the identity node will need to be created */
171  tntNode eye = tntNodeCreateEyeOp(os->vals[0]); /* Identity node created using basis operator*/
172  tntNodeArray Lterms, Rterms; /* Arrays that contain both onsite and nearest neighbour left and right terms */
173
174  /* Create a node array that contains the nearest neighbour terms and the onsite terms */
175  Lterms = tntNodeArrayAlloc(nnl->sz + os->sz*2); /* twice the number of onsite terms, as require one term for .. */
176  Rterms = tntNodeArrayAlloc(nnl->sz + os->sz*2); /* ... operator on left site, and one term for operator on right site */
177
178  /* Create a complex array for the parameters */
179  two_site_params = tntComplexArrayAlloc(Lterms.sz); /* Array for parameters on each site */
180
181  /* First put the nearest neighbour terms in the array */
182  for (i = 0; i < nnl->sz; i++) {
183  Lterms.vals[i] = tntNodeCopy(nnl->vals[i]);
184  Rterms.vals[i] = tntNodeCopy(nnr->vals[i]);
185  }
186
187  /* Then put the on-site terms in the array */
188  for (i = 0; i < os->sz; i++) {
189  /* For the first term put the on-site operator on the left */
190  Lterms.vals[i*2 + nnl->sz] = tntNodeCopy(os->vals[i]);
191  Rterms.vals[i*2 + nnl->sz] = tntNodeCopy(eye);
192
193  /* For the first term put the on-site operator on the right */
194  Rterms.vals[i*2 + 1 + nnl->sz] = tntNodeCopy(os->vals[i]);
195  Lterms.vals[i*2 + 1 + nnl->sz] = tntNodeCopy(eye);
196  }
197
198  /* Now determine parameters for each two site gate, depending on whether they vary in space or not */
199  if ((0 == varparam) && (L > 2)) {
200  /* For the case where the operators are uniform, only the first and last nodes will be different, with
201  all the central nodes being the same */
202
203  /* Two site parameters are the same on every site */
204  for (i = 0; i < nnl->sz; i++) two_site_params.vals[i] = nnparam->vals[i];
205
206  /* For the first site, the left on-site term is 1*parameter, right on-site term is 0.5*parameter */
207  for (i = 0; i < os->sz; i++) {
208  two_site_params.vals[i*2 + nnl->sz].re = osparam->vals[i].re;
209  two_site_params.vals[i*2 + nnl->sz].im = osparam->vals[i].im;
210  two_site_params.vals[i*2 + 1 + nnl->sz].re = osparam->vals[i].re/2.0;
211  two_site_params.vals[i*2 + 1 + nnl->sz].im = osparam->vals[i].im/2.0;
212  }
213  /* Create the node and put into the first term in the array */
214  Parr.vals[0] = tntMpsCreateTwoSiteOp(&Lterms,&Rterms,&two_site_params);
215  /* Take the matrix exponential to form the two-site propagator */
216  tntNodeExp(Parr.vals[0], hs, "DE", "UV");
217
218  /* For the central sites, both the left and right on-site terms are 0.5*parameter */
219  for (i = 0; i < os->sz; i++) {
220  two_site_params.vals[i*2 + nnl->sz].re = osparam->vals[i].re/2.0;
221  two_site_params.vals[i*2 + nnl->sz].im = osparam->vals[i].im/2.0;
222  two_site_params.vals[i*2 + 1 + nnl->sz].re = osparam->vals[i].re/2.0;
223  two_site_params.vals[i*2 + 1 + nnl->sz].im = osparam->vals[i].im/2.0;
224  }
225  /* Create the node and put into the second term in the array */
226  Parr.vals[1] = tntMpsCreateTwoSiteOp(&Lterms,&Rterms,&two_site_params);
227  /* Take the matrix exponential to form the two-site propagator */
228  tntNodeExp(Parr.vals[1], hs, "DE", "UV");
229
230  /* Copy this node for all central terms */
231  for (j = 2; j<L-2; j++) {
232  Parr.vals[j] = tntNodeCopy(Parr.vals[1]);
233  }
234
235  /* For the last site, the left on-site term is 0.5*parameter, right on-site term is 1*parameter */
236  for (i = 0; i < os->sz; i++) {
237  two_site_params.vals[i*2 + nnl->sz].re = osparam->vals[i].re/2.0;
238  two_site_params.vals[i*2 + nnl->sz].im = osparam->vals[i].im/2.0;
239  two_site_params.vals[i*2 + 1 + nnl->sz].re = osparam->vals[i].re;
240  two_site_params.vals[i*2 + 1 + nnl->sz].im = osparam->vals[i].im;
241  }
242  /* Create the node and put into the first term in the array */
243  Parr.vals[L-2] = tntMpsCreateTwoSiteOp(&Lterms,&Rterms,&two_site_params);
244  /* Take the matrix exponential to form the two-site propagator */
245  tntNodeExp(Parr.vals[L-2], hs, "DE", "UV");
246  } else if ((1 == varparam) && (L > 2)) {
247
248  /* Otherwise every node in the array will be different */
249  double scL, scR; /* Amount to scale onsite parameter by - depends on site */
250  unsigned ind, indL, indR; /* index in array */
251
252  /* Set the terms for all sites */
253  for (j = 0; j<L-1; j++) {
254  /* Set the nearest neighbour terms for the central sites */
255  for (i = 0; i < nnl->sz; i++) {
256  ind = (1 == nnparam->numcols) ? i : i + j*nnl->sz;
257  two_site_params.vals[i] = nnparam->vals[ind];
258  }
259
260  /* Set the onsite terms for the central sites */
261  scL = (0==j) ? 1.0 : 2.0;
262  scR = (L-2 == j) ? 1.0 : 2.0;
263  for (i = 0; i < os->sz; i++) {
264  indL = (1 == osparam->numcols) ? i : i + j*os->sz;
265  indR = (1 == osparam->numcols) ? i : i + (j+1)*os->sz;
266  two_site_params.vals[i*2 + nnl->sz].re = osparam->vals[indL].re/scL;
267  two_site_params.vals[i*2 + nnl->sz].im = osparam->vals[indL].im/scL;
268  two_site_params.vals[i*2 + 1 + nnl->sz].re = osparam->vals[indR].re/scR;
269  two_site_params.vals[i*2 + 1 + nnl->sz].im = osparam->vals[indR].im/scR;
270  }
271
272  /* Create the node and put into the first term in the array */
273  Parr.vals[j] = tntMpsCreateTwoSiteOp(&Lterms,&Rterms,&two_site_params);
274  /* Take the matrix exponential to form the two-site propagator */
275  tntNodeExp(Parr.vals[j], hs, "DE", "UV");
276  }
277  } else {
278  /* Special behaviour if it is only a two site system */
279
280  /* nn terms cannot vary accross just two sites, so no need to check */
281  for (i = 0; i < nnl->sz; i++) two_site_params.vals[i] = nnparam->vals[i];
282
283  /* For onsite term in two site system use 1*parameter for both left and right site as there will be no double counting */
284  /* Take into account a possible variation in parameters by using the number of columns in osparam in the indexing */
285  for (i = 0; i < os->sz; i++) {
286  two_site_params.vals[i*2 + nnl->sz] = osparam->vals[i];
287  two_site_params.vals[i*2 + 1 + nnl->sz] = osparam->vals[i + (osparam->numcols-1)*os->sz];
288  }
289
290  /* Create the node and put into the first and only term in the array */
291  Parr.vals[0] = tntMpsCreateTwoSiteOp(&Lterms,&Rterms,&two_site_params);
292  /* Take the matrix exponential to form the two-site propagator */
293  tntNodeExp(Parr.vals[0], hs, "DE", "UV");
294  }
295
296  /* Free the arrays and nodes */
297  tntComplexArrayFree(&two_site_params);
298  tntNodeArrayFree(&Lterms);
299  tntNodeArrayFree(&Rterms);
300  tntNodeFree(&eye);
301
302  /* Turn the warning back on */
304
305  /* Return the array of nodes */
306  return Parr;
307 }
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
void tntNodeArrayFree(tntNodeArray *arr)
Definition: tntArray.c:668
void tntNodeExp(tntNode A, tntComplex c, tntLegLabel rowlegs, tntLegLabel collegs)
tntNode tntNodeCreateEyeOp(tntNode A)
Definition: tntNodeUtil.c:195
tntNodeArray tntNodeArrayAlloc(unsigned num_elems)
Definition: tntArray.c:428
void tntNodeFree(tntNode *A)
Definition: tntNodeUtil.c:275
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:392
Definition: tnt.h:65
tntNodeArray tntMpsCreatePropArray(unsigned L, tntComplex h, tntNodeArray *nnl, tntNodeArray *nnr, tntComplexArray *nnparam, tntNodeArray *os, tntComplexArray *osparam)
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray.c:650
double re
Definition: tnt.h:66
void tntSysQNClearWarnOn(void)
Definition: tntUtil.c:17
double im
Definition: tnt.h:67
void tntSysQNClearWarnOff(void)
Definition: tntUtil.c:35
tntNode tntMpsCreateTwoSiteOp(tntNodeArray *Lnodes, tntNodeArray *Rnodes, tntComplexArray *params)