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
tntMpsProcessInitOptions.c
1 /*
2  * Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch (Oxford) and Chris Goodyer (NAG)
3  * $LastChangedDate$
4  * (c) University of Oxford 2014
5  */
6 
12 /* Include the internal header for the TNT library */
13 #include "tntMpsInternal.h"
14 
15 #include <getopt.h>
16 #include <ctype.h>
17 
18 #define TNT_SH_NUM_TERMS 12
19 #define TNT_SH_NUM_NN_TERMS 9
20 
21 #define TNT_BS_NUM_TERMS 7
22 #define TNT_BS_NUM_NN_TERMS 3
23 
210 void tntMpsProcessSysOptions(int argc,
211  char **argv,
212  tntNetwork *mps,
213  tntNetwork *mpo,
214  tntNetwork *prop,
215  tntMpsExOp *ExOp)
216 {
217 
218  static struct option all_long_options[] = { /* Long options for all system */
219  {"dt", required_argument, 0, 1},
220  {"idt", required_argument, 0, 2},
221  {"rand-state", required_argument, 0, 3},
222  {"qnum-rand-state", required_argument, 0, 4},
223  {"config-state", required_argument, 0, 5},
224  {"qnum-config-state", required_argument, 0, 6},
225  {"length", required_argument, 0, 7},
226  {"system", required_argument, 0, 8},
227  {"qnum-eye-state", no_argument, 0, 9},
228  {"eye-state", no_argument, 0, 10},
229  {0, 0, 0, 0}
230  };
231  static struct option sh_long_options[] = { /* Long options for a spin half system */
232  {"Jxx", required_argument, 0, 1},
233  {"Jyy", required_argument, 0, 2},
234  {"Jzz", required_argument, 0, 3},
235  {"Jxy", required_argument, 0, 4},
236  {"Jyx", required_argument, 0, 5},
237  {"Jyz", required_argument, 0, 6},
238  {"Jzy", required_argument, 0, 7},
239  {"Jzx", required_argument, 0, 8},
240  {"Jxz", required_argument, 0, 9},
241  {"Bx", required_argument, 0, 10},
242  {"By", required_argument, 0, 11},
243  {"Bz", required_argument, 0, 12},
244  {"2S", required_argument, 0, 13},
245  {0, 0, 0, 0}
246  };
247  static struct option bs_long_options[] = { /* Long options for a boson system */
248  {"Jb", required_argument, 0, 1},
249  {"Vb", required_argument, 0, 2},
250  {"param-driv", required_argument, 0, 3},
251  {"Ub", required_argument, 0, 4},
252  {"mub", required_argument, 0, 5},
253  {"E-harm", required_argument, 0, 6},
254  {"coher-driv", required_argument, 0, 7},
255  {"jc-harm", required_argument, 0, 8},
256  {"n-max", required_argument, 0, 9},
257  {0, 0, 0, 0}
258  };
259  int option_index;
260  int c; /* Variable that contains the returned option code from getopt_long */
261  char **argv_cp; /* Used to hold a copy of the arguments - this will be used by the function instead of the original arguments */
262  unsigned numnn = 0; /* Number of nearest neighbour terms in the Hamiltonian */
263  unsigned numos = 0; /* Number of onsite terms in the Hamiltonian */
264  unsigned no_operators = 0; /* Marks whether either H or U should be calculated */
265  double dt = 0.0, idt = 0.0; /* initial time step sizes */
266  tntComplex dtc; /* The complex time step size */
267  unsigned MPStype = 0; /* Use 1 for random, use 2 for config state, 0 if no MPS type has been assigned yet. */
268  char config[TNT_STRLEN]; /* Code for the configuration */
269  unsigned L=0; /* Length of system - set to zero initially to see if it has been set */
270  unsigned D; /* Starting internal dimension of the random MPS */
271  int Q; /* Quantum number for symmetric MPS */
272  tntNodeArray nnL, nnR, os; /* Arrays for nearest neighbour and onsite operators */
273  tntComplexArray nnparam, osparam; /* Arrays for the parameters for the nearest neighbour and onsite operators */
274  unsigned loop, j;
275 
276  /* Set optind to 1 in case option getting routine has already been called previously */
277  optind = 1;
278 
279  /* Set opterr to zero to prevent error message for unknown erguments, since they may be being treated by another process options function */
280  opterr = 0;
281 
282  /* Make a copy of argv - getopt_long will permute the arguments, and they may be required by another processing function */
283  argv_cp = (char **) malloc(argc*sizeof(char *));
284  for (c = 0; c<argc; c++) {
285  argv_cp[c] = argv[c];
286  }
287 
288  /* Check whether either of H or U needs to be built. If not, then all the parameters will be ignored. If it does, then create all the spin half operators */
289  if ((NULL == prop) && (NULL == mpo)) {
290  no_operators = 1;
291  }
292 
293  tntPrintf("----------------------------------------------------------\n");
294  tntPrintf("Creating MPS system with the following parameters: \n");
295  tntPrintf("----------------------------------------------------------\n");
296 
297  /* Loop over all input arguments */
298  while ((c = getopt_long(argc, argv_cp, "", all_long_options, &option_index)) != -1) {
299  switch (c) {
300  case 1: /* This option relates to the real time step size */
301  if (NULL == prop) {
302  tntWarningPrint("The input parameter dt is being ignored as U is NULL|");
303  } else {
304  dt = atof(optarg);
305  tntPrintf("dt: <<<%g>>>\n", dt);
306  }
307  break;
308  case 2: /* This option relates to the imaginary time step size */
309  if (NULL == prop) {
310  tntWarningPrint("The input parameter idt is being ignored as U is NULL|");
311  } else {
312  idt = atof(optarg);
313  tntPrintf("idt: <<<%g>>>\n", idt);
314  }
315  break;
316  case 3: /* This option relates to specifiying a random starting MPS */
317  if (NULL == mps) {
318  tntWarningPrint("The input parameter rand-state is being ignored as argument mps is NULL|");
319  } else if (MPStype) {
320  tntWarningPrint("The input parameter rand-state is being ignored as the MPS type has already been set|");
321  } else {
322  MPStype = 1;
323  D = abs(atoi(optarg));
324  tntPrintf("Random MPS being creating with starting internal dimension: <<<%d>>>\n",D);
325  }
326  break;
327  case 4: /* This option relates to specifying a given quantum number for the starting random MPS, and conserving this number */
328  if (NULL == mps) {
329  tntWarningPrint("The input parameter qnum-eye-state is being ignored as argument mps is NULL|");
330  } else if (MPStype) {
331  tntWarningPrint("The input parameter qnum-eye-state is being ignored as the MPS type has already been set|");
332  } else {
333  MPStype = 3;
334  /* Extract total quantum number for the system */
335  Q = atoi(optarg);
336  /* Set the symmetry type for the system */
337  tntSymmTypeSet("U(1)", 1);
338  tntPrintf("Creating MPS with initial quantum number: <<<%d>>>\n",Q);
339  }
340  break;
341  case 5: /* This option relates to specifying a configuration for the starting MPS */
342  if (NULL == mps) {
343  tntWarningPrint("The input parameter config-state is being ignored as argument mps is NULL|");
344  } else if (MPStype) {
345  tntWarningPrint("The input parameter config-state is being ignored as the MPS type has already been set|");
346  } else {
347  MPStype = 2;
348  strcpy(config, optarg);
349  tntPrintf("Creating MPS with configuration: <<<%s>>>\n",config);
350  }
351  break;
352  case 6: /* This option relates to specifying a configuration for the starting MPS, and conserving its quantum number */
353  if (NULL == mps) {
354  tntWarningPrint("The input parameter qnum-state is being ignored as argument mps is NULL|");
355  } else if (MPStype) {
356  tntWarningPrint("The input parameter qnum-state is being ignored as the MPS type has already been set|");
357  } else {
358  MPStype = 2;
359  strcpy(config, optarg);
360  tntPrintf("Creating number conserving MPS with configuration: <<<%s>>>\n",config);
361  /* Set the symmetry type for the system */
362  tntSymmTypeSet("U(1)", 1);
363  }
364  break;
365  case 7: /* This option (required) relates to specifying the length of the starting MPS */
366  L = abs(atoi(optarg));
367  tntPrintf("Length: <<<%d>>>\n",L);
368  break;
369  case 8: /* This option (required) relates to the type of system that will be created */
370  if (tntSysTypeGet()) {
371  tntWarningPrint("The input parameter system is being ignored as the system type has already been set|");
372  } else if (!strcmp(optarg,"spin")) {
373  tntPrintf("system: <<<spin>>>\n");
374  tntSysTypeSet(TNT_MPS_SPIN);
375  } else if (!strcmp(optarg,"boson")) {
376  tntPrintf("system: <<<boson>>>\n");
377  tntSysTypeSet(TNT_MPS_BOSON);
378  } else {
379  tntErrorPrint("System type %s is not known",optarg);
380  }
381  break;
382  case 9: /* This option (required) relates to the type of system that will be created */
383  if (NULL == mps) {
384  tntWarningPrint("The input parameter qnum-eye-state is being ignored as argument mps is NULL|");
385  } else if (MPStype) {
386  tntWarningPrint("The input parameter qnum-eye-state is being ignored as the start-state type has already been set|");
387  } else {
388  MPStype = 4;
389  tntPrintf("Creating identity density matrix as start state, with quantum number conservation turned on\n");
390  tntSymmTypeSet("U(1)", 1);
391  }
392  case 10: /* This option (required) relates to the type of system that will be created */
393  if (NULL == mps) {
394  tntWarningPrint("The input parameter qnum-eye-state is being ignored as argument mps is NULL|");
395  } else if (MPStype) {
396  tntWarningPrint("The input parameter qnum-eye-state is being ignored as the start-state type has already been set|");
397  } else {
398  MPStype = 4;
399  tntPrintf("Creating identity density matrix as start state\n");
400  }
401  default:
402  break;
403  }
404  }
405 
406  if (0 == L) {
407  tntErrorPrint("The system cannot be initialised as the length has not been set"); /* NO_COVERAGE */
408  } /* NO_COVERAGE */
409 
410  if ((0 == MPStype)&&(mps != NULL)) {
411  tntErrorPrint("The wave function cannot be created as the starting state has not been set"); /* NO_COVERAGE */
412  } /* NO_COVERAGE */
413 
414  if (TNT_MPS_SPIN == tntSysTypeGet()) { /* Do this bit if spin system chosen and there are operators */
415 
416  unsigned TwoS = 1; /* Spin for each particle in the system */
417  tntComplex sh_params[TNT_SH_NUM_TERMS]; /* Multipliers for each of the terms in the Hamiltonian */
418  unsigned sh_param_present[TNT_SH_NUM_TERMS] = {0,0,0,0,0,0,0,0,0,0,0,0}; /* Specifies whether each term is present */
419  char sh_param_strs[TNT_SH_NUM_TERMS][4] = {"Jxx","Jyy","Jzz","Jxy","Jyx","Jyz","Jzy","Jzx","Jxz","Bx","By","Bz"}; /* String for each term */
420 
421  /* Make a copy of argv - getopt_long will permute the arguments, and they may be required by another processing function */
422  for (c = 0; c<argc; c++) {
423  argv_cp[c] = argv[c];
424  }
425 
426  /* Set optind to 1 in case option getting routine has already been called previously */
427  optind = 1;
428 
429  tntPrintf("----------------------------------------------------------\n");
430  tntPrintf("Creating spin system with the following parameters: \n");
431  tntPrintf("----------------------------------------------------------\n");
432 
433  /* Loop over all input arguments */
434  while ((c = getopt_long(argc, argv_cp, "", sh_long_options, &option_index)) != -1) {
435  if (c <= TNT_SH_NUM_TERMS) { /* These options relate to terms in the Hamiltonian */
436  /* Populate the array of parameters */
437  sh_params[c-1] = tntAtoC(optarg);
438  sh_param_present[c-1] = 1;
439  if (fabs(sh_params[c-1].re) > DBL_EPSILON && fabs(sh_params[c-1].im) > DBL_EPSILON) {
440  tntPrintf("%s : <<<%g%+gi>>>\n", sh_param_strs[c-1], sh_params[c-1].re, sh_params[c-1].im);
441  } else if (fabs(sh_params[c-1].re) > DBL_EPSILON) {
442  tntPrintf("%s : <<<%g>>>\n", sh_param_strs[c-1], sh_params[c-1].re);
443  } else if (fabs(sh_params[c-1].im) > DBL_EPSILON) {
444  tntPrintf("%s : <<<%gi>>>\n", sh_param_strs[c-1], sh_params[c-1].im);
445  } else {
446  tntPrintf("%s : <<<0>>>\n", sh_param_strs[c-1]);
447  }
448  if (c <= TNT_SH_NUM_NN_TERMS) {
449  numnn++;
450  } else {
451  numos++;
452  }
453  } else if (13 == c) { /* Term for the spin for each particle */
454  TwoS = abs(atoi(optarg));
455  tntPrintf("2S : <<<%d>>>\n", TwoS);
456  }
457  }
458 
459  tntNode Sx, Sy, Sz, Sp, Sm, eye; /* Spin operators */
460 
461  /* Creating the Pauli operators */
462  tntMpsCreateSpinOp(TwoS, &Sx, &Sy, &Sz, &Sp, &Sm, &eye);
463 
464  /* Set the operator for the basis */
465  tntSysBasisOpSet(Sz);
466 
467  if (!no_operators) {
468 
469  /* First check whether the operators actually exist */
470  if (0 == numnn && 0 == numos) {
471  tntErrorPrint("Cannot create operators as all the parameters are zero"); /* NO_COVERAGE */
472  }
473 
474 
475  /* First step: if quantum number conservation is turned on, check for non-invariant operators */
476  /* Only non-invariant operator combinations are: Jxx = Jyy, Jxy = -Jyx, Jzz anything, Bz anything */
477  /* if Jxx, Jyy, Jxy or Jyx are present convert to being expressed as combinations of Sp and Sm (invariant) rather than Sx and Sy (done below) */
478  /* if non-allowed terms are present delete them and present warning to user (here) */
479  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
480  /* Look for non-invariant parts, which cannot be re-expressed as an invariant + non-invariant part */
481  if (sh_param_present[5] || sh_param_present[6] || sh_param_present[7] || sh_param_present[8] || sh_param_present[9] || sh_param_present[10]) {
482  char full_msg[TNT_STRLEN];
483  strcpy(full_msg,"Your Hamiltonian is not invariant under the symmetry type chosen|");
484  for (loop = 5; loop < 11; loop++) {
485  if (sh_param_present[loop]) {
486  strcat(full_msg, sh_param_strs[loop]);
487  strcat(full_msg, " ");
488  sh_param_present[loop] = 0;
489  if (loop < TNT_SH_NUM_NN_TERMS) {
490  numnn--;
491  } else {
492  numos--;
493  }
494  }
495  }
496  strcat(full_msg, "being set to zero|");
497  tntWarningPrint(full_msg);
498  }
499  }
500 
501  /* Create arrays to hold the terms and parameters that are present */
502  if (numnn) {
503  nnL = tntNodeArrayAlloc(numnn);
504  nnR = tntNodeArrayAlloc(numnn);
505  nnparam = tntComplexArrayAlloc(numnn);
506  } else {
507  nnL.sz = 0;
508  nnR.sz = 0;
509  nnparam.numrows = nnparam.numcols = 0;
510  }
511  if (numos) {
512  os = tntNodeArrayAlloc(numos);
513  osparam = tntComplexArrayAlloc(numos);
514  } else {
515  os.sz = 0;
516  osparam.numrows = osparam.numcols = 0;
517  }
518 
519  /* reset the counters */
520  numnn = numos = 0;
521 
522  /* Second step: if quantum number conservation is turned on, check for non-invariant operators Jxx, Jyy, Jxy, Jyx */
523  /* These can be expressed as combinations of Sp and Sm (invariant) rather than Sx and Sy */
524  /* Only non-invariant operator combinations are: Jxx = Jyy, Jxy = -Jyx. If different values are given present warning to user */
525  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
526  /* Look for Jxx or Jyy term */
527  if (sh_param_present[0] || sh_param_present[1]) {
528  if (fabs(sh_params[0].re-sh_params[1].re) + fabs(sh_params[0].im-sh_params[1].im) > DBL_EPSILON) {
529  tntWarningPrint("Your Hamiltonian is not invariant under the symmetry type chosen| Only Jxx=Jyy allowed - values being changed to average of Jxx and Jyy");
530  }
531  /* create the term 0.25*(Jxx+Jyy)(SpSm + SmSp) = JxxSxSx + JyySySy (+ any non-invariant parts being ignored) */
532  nnL.vals[numnn] = tntNodeCopy(Sp); /* Operator acting on left */
533  nnR.vals[numnn] = tntNodeCopy(Sm); /* Operator acting on right */
534  nnparam.vals[numnn].re = 0.25*(sh_params[0].re + sh_params[1].re); /* Copy over parameter */
535  nnparam.vals[numnn].im = 0.25*(sh_params[0].im + sh_params[1].im); /* Copy over parameter */
536  numnn++;
537  nnL.vals[numnn] = tntNodeCopy(Sm); /* Operator acting on left */
538  nnR.vals[numnn] = tntNodeCopy(Sp); /* Operator acting on right */
539  nnparam.vals[numnn].re = 0.25*(sh_params[0].re + sh_params[1].re); /* Copy over parameter */
540  nnparam.vals[numnn].im = 0.25*(sh_params[0].im + sh_params[1].im); /* Copy over parameter */
541  numnn++;
542  /* Since values have already been accounted for, set flag to zero so they are not included again */
543  sh_param_present[0] = sh_param_present[1] = 0;
544  }
545  /* Look for Jxy or Jyx term */
546  if (sh_param_present[3] || sh_param_present[4]) {
547  if (fabs(sh_params[3].re + sh_params[4].re) + fabs(sh_params[3].im + sh_params[4].im) > DBL_EPSILON) {
548  tntWarningPrint("Your Hamiltonian is not invariant under the symmetry type chosen|Only Jxy=-Jyx allowed - values being changed to Jxy_new = -Jyx_new = (Jxy - Jyx)/2");
549  }
550  /* create the term 0.25i*(Jxy-Jyx)(SpSm - SmSp) = JxySxSy + JyxSySx (+ any non-invariant parts being ignored) */
551  nnL.vals[numnn] = tntNodeCopy(Sp); /* Operator acting on left */
552  nnR.vals[numnn] = tntNodeCopy(Sm); /* Operator acting on right */
553  nnparam.vals[numnn].im = 0.25*(sh_params[3].re - sh_params[4].re); /* Copy over parameter */
554  nnparam.vals[numnn].re = -0.25*(sh_params[3].im - sh_params[4].im); /* Copy over parameter */
555  numnn++;
556  nnL.vals[numnn] = tntNodeCopy(Sm); /* Operator acting on left */
557  nnR.vals[numnn] = tntNodeCopy(Sp); /* Operator acting on right */
558  nnparam.vals[numnn].im = 0.25*(sh_params[4].re - sh_params[3].re); /* Copy over parameter */
559  nnparam.vals[numnn].re = -0.25*(sh_params[4].im - sh_params[3].im); /* Copy over parameter */
560  numnn++;
561  /* Since values have already been accounted for, set flag to zero so they are not included again */
562  sh_param_present[3] = sh_param_present[4] = 0;
563  }
564  }
565 
566 
567  /* Loop through parameters and assign them to the correct array */
568  for (loop = 0; loop < TNT_SH_NUM_TERMS; loop++) {
569  if (sh_param_present[loop]) {
570  switch (loop+1) {
571  case 1: /* Term Jxx */
572  /* ---- Set the nodes to make the term ---- */
573  nnL.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on left */
574  nnR.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on right */
575  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
576  numnn++;
577  break;
578  case 2: /* Term Jyy */
579  /* ---- Set the nodes to make the term ---- */
580  nnL.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on left */
581  nnR.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on right */
582  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
583  numnn++;
584  break;
585  case 3: /* Term Jzz */
586  /* ---- Set the nodes to make the term ---- */
587  nnL.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on left */
588  nnR.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on right */
589  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
590  numnn++;
591  break;
592  case 4: /* Term Jxy */
593  /* ---- Set the nodes to make the term ---- */
594  nnL.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on left */
595  nnR.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on right */
596  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
597  numnn++;
598  break;
599  case 5: /* Term Jyx */
600  /* ---- Set the nodes to make the term ---- */
601  nnL.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on left */
602  nnR.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on right */
603  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
604  numnn++;
605  break;
606  case 6: /* Term Jyz */
607  /* ---- Set the nodes to make the term ---- */
608  nnL.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on left */
609  nnR.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on right */
610  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
611  numnn++;
612  break;
613  case 7: /* Term Jzy */
614  /* ---- Set the nodes to make the term ---- */
615  nnL.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on left */
616  nnR.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on right */
617  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
618  numnn++;
619  break;
620  case 8: /* Term Jzx */
621  /* ---- Set the nodes to make the term ---- */
622  nnL.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on left */
623  nnR.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on right */
624  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
625  numnn++;
626  break;
627  case 9: /* Term Jxz */
628  /* ---- Set the nodes to make the term ---- */
629  nnL.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on left */
630  nnR.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on right */
631  nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
632  numnn++;
633  break;
634  case 10: /* Term Bx */
635  /* ---- Set the nodes to make the term ---- */
636  os.vals[numos] = tntNodeCopy(Sx); /* Operator acting on site */
637  osparam.vals[numos] = sh_params[loop]; /* Copy over parameter */
638  numos++;
639  break;
640  case 11: /* Term By */
641  /* ---- Set the nodes to make the term ---- */
642  os.vals[numos] = tntNodeCopy(Sy); /* Operator acting on site */
643  osparam.vals[numos] = sh_params[loop]; /* Copy over parameter */
644  numos++;
645  break;
646  case 12: /* Term Bz */
647  /* ---- Set the nodes to make the term ---- */
648  os.vals[numos] = tntNodeCopy(Sz); /* Operator acting on site */
649  osparam.vals[numos] = sh_params[loop]; /* Copy over parameter */
650  numos++;
651  break;
652  default:
653  break;
654  }
655  }
656  }
657  }
658  /* Free the spin operators */
659  tntNodeFree(&Sx);
660  tntNodeFree(&Sy);
661  tntNodeFree(&Sz);
662  tntNodeFree(&Sp);
663  tntNodeFree(&Sm);
664  tntNodeFree(&eye);
665  }
666  if (TNT_MPS_BOSON == tntSysTypeGet()) { /* Do this bit if boson system chosen and there are operators */
667 
668  double jc = (L-1.0)/2.0; /* Centre for the harmonic trap */
669  unsigned n_max = 1; /* Maximum number of bosons on each site */
670  tntComplex bs_params[TNT_BS_NUM_TERMS]; /* Multipliers for each of the terms in the Hamiltonian */
671  unsigned bs_param_present[TNT_BS_NUM_TERMS] = {0,0,0,0,0,0,0}; /* Specifies whether each term is present */
672  char bs_param_strs[TNT_BS_NUM_TERMS][12] = {"Jb","Vb","param-driv","Ub","mub","E-harm","coher-driv"}; /* String for each term */
673 
674  tntPrintf("----------------------------------------------------------\n");
675  tntPrintf("Creating boson system with the following parameters: \n");
676  tntPrintf("----------------------------------------------------------\n");
677 
678 
679  /* Make a copy of argv - getopt_long will permute the arguments, and they may be required by another processing function */
680  for (c = 0; c<argc; c++) {
681  argv_cp[c] = argv[c];
682  }
683 
684  /* Set optind to 1 in case option getting routine has already been called previously */
685  optind = 1;
686 
687  /* Reset counters for number of operators */
688  numnn = numos = 0;
689 
690  /* Loop over all input arguments */
691  while ((c = getopt_long(argc, argv_cp, "", bs_long_options, &option_index)) != -1) {
692  if (c <= TNT_BS_NUM_TERMS) { /* These options relate to terms in the Hamiltonian */
693  /* Populate the array of parameters */
694  bs_params[c-1] = tntAtoC(optarg);
695  bs_param_present[c-1] = 1;
696  if (fabs(bs_params[c-1].re) > DBL_EPSILON && fabs(bs_params[c-1].im) > DBL_EPSILON) {
697  tntPrintf("%s : <<<%g%+gi>>>\n", bs_param_strs[c-1], bs_params[c-1].re, bs_params[c-1].im);
698  } else if (fabs(bs_params[c-1].re) > DBL_EPSILON) {
699  tntPrintf("%s : <<<%g>>>\n", bs_param_strs[c-1], bs_params[c-1].re);
700  } else if (fabs(bs_params[c-1].im) > DBL_EPSILON) {
701  tntPrintf("%s : <<<%gi>>>\n", bs_param_strs[c-1], bs_params[c-1].im);
702  } else {
703  tntPrintf("%s : <<<0>>>\n", bs_param_strs[c-1]);
704  }
705  if ((1 == c)||(3 == c)) { /* Two nearest neighbour terms for hopping or parametric driving */
706  numnn += 2;
707  } else if (c <= TNT_BS_NUM_NN_TERMS) {
708  numnn++;
709  } else if (7 == c) { /* Two onsite terms for coherent driving */
710  numos += 2;
711  } else {
712  numos++;
713  }
714  } else if (8 == c) { /* Term for centre of harmonic potential */
715  jc = atof(optarg);
716  tntPrintf("jc-harm : <<<%g>>>\n", jc);
717  } else if (9 == c) { /* Term for maximum number of atoms per site */
718  n_max = abs(atoi(optarg));
719  tntPrintf("n-max : <<<%d>>>\n", n_max);
720  }
721  }
722 
723  tntNode b_op, bdag_op, n_op, os_int, eye; /* Boson operators */
724 
725  /* Creating the boson operators */
726  tntMpsCreateBosonOp(n_max, &b_op, &bdag_op, &n_op, &os_int, &eye);
727 
728  /* Set the operator for the basis */
729  tntSysBasisOpSet(n_op);
730 
731  if (!no_operators) {
732 
733  /* First check whether the operators actually exist */
734  if (0 == numnn && 0 == numos) {
735  tntErrorPrint("Cannot create operators as all the parameters are zero"); /* NO_COVERAGE */
736  exit(1); /* NO_COVERAGE */
737  }
738 
739  /* First step: if quantum number conservation is turned on, check for non-invariant operators i.e. param-driv or coher-driv */
740  if (TNT_SYMM_U1 == tntSymmTypeGet()) {
741  /* Look for terms */
742  if (bs_param_present[2] && bs_param_present[6]) {
743  tntWarningPrint("Your Hamiltonian is not invariant under the symmetry type chosen|%s and %s being set to zero|", bs_param_strs[2], bs_param_strs[6]);
744  bs_param_present[2] = bs_param_present[6] = 0;
745  numnn -= 2;
746  numos -= 2;
747  } else if (bs_param_present[2]) {
748  tntWarningPrint("Your Hamiltonian is not invariant under the symmetry type chosen|%s being set to zero|", bs_param_strs[2]);
749  bs_param_present[2] = 0;
750  numnn -= 2;
751  } else if (bs_param_present[6]) {
752  tntWarningPrint("Your Hamiltonian is not invariant under the symmetry type chosen|%s being set to zero|", bs_param_strs[6]);
753  bs_param_present[6] = 0;
754  numos -= 2;
755  }
756 
757  }
758 
759  /* Create arrays to hold the terms and parameters that are present */
760  if (numnn) {
761  nnL = tntNodeArrayAlloc(numnn);
762  nnR = tntNodeArrayAlloc(numnn);
763  nnparam = tntComplexArrayAlloc(numnn);
764  } else {
765  nnL.sz = 0;
766  nnR.sz = 0;
767  nnparam.numrows = nnparam.numcols = 0;
768  }
769  if (numos) {
770  os = tntNodeArrayAlloc(numos);
771  if (bs_param_present[5]) /* If there is a harmonic trapping term, then there will be different onsite term for each site */
772  osparam = tntComplexArrayAlloc(numos,L);
773  else
774  osparam = tntComplexArrayAlloc(numos);
775  } else {
776  os.sz = 0;
777  osparam.numrows = osparam.numcols = 0;
778  }
779 
780  /* reset the counters */
781  numnn = numos = 0;
782 
783  /* Loop through parameters and assign them to the correct array */
784  for (loop = 0; loop < TNT_BS_NUM_TERMS; loop++) {
785  if (bs_param_present[loop]) {
786  switch (loop+1) {
787  case 1: /* Term Jb - two hopping terms required */
788  /* ---- Set the nodes to make the term ---- */
789  nnL.vals[numnn] = tntNodeCopy(b_op); /* Operator acting on left */
790  nnR.vals[numnn] = tntNodeCopy(bdag_op); /* Operator acting on right */
791  nnparam.vals[numnn].re = -bs_params[loop].re; /* Copy over parameter */
792  nnparam.vals[numnn].im = -bs_params[loop].im; /* Copy over parameter */
793  nnL.vals[numnn+1] = tntNodeCopy(bdag_op); /* Operator acting on left */
794  nnR.vals[numnn+1] = tntNodeCopy(b_op); /* Operator acting on right */
795  nnparam.vals[numnn+1].re = -bs_params[loop].re; /* Copy over cc parameter */
796  nnparam.vals[numnn+1].im = bs_params[loop].im; /* Copy over cc parameter */
797  numnn += 2;
798  break;
799  case 2: /* Term Vb: nearest neighbour interaction */
800  /* ---- Set the nodes to make the term ---- */
801  nnL.vals[numnn] = tntNodeCopy(n_op); /* Operator acting on left */
802  nnR.vals[numnn] = tntNodeCopy(n_op); /* Operator acting on right */
803  nnparam.vals[numnn] = bs_params[loop]; /* Copy over parameter */
804  numnn++;
805  break;
806  case 3: /* Term param-driv: parametric driving */
807  /* ---- Set the nodes to make the term ---- */
808  nnL.vals[numnn] = tntNodeCopy(b_op); /* Operator acting on left - set to b_op */
809  nnR.vals[numnn] = tntNodeCopy(b_op); /* Operator acting on right - set to b_op */
810  nnparam.vals[numnn] = bs_params[loop]; /* Copy over parameter */
811 
812  nnL.vals[numnn+1] = tntNodeCopy(bdag_op); /* Operator acting on right - set to bdag_op */
813  nnR.vals[numnn+1] = tntNodeCopy(bdag_op); /* Operator acting on left - set to bdag_op */
814  nnparam.vals[numnn+1].re = bs_params[loop].re; /* Copy over cc parameter */
815  nnparam.vals[numnn+1].im = -bs_params[loop].im; /* Copy over cc parameter */
816 
817  numnn += 2;
818  break;
819  case 4: /* Term Ub - onsite interaction */
820  /* ---- Set the nodes to make the term ---- */
821  os.vals[numos] = tntNodeCopy(os_int); /* Operator acting on site */
822  if (bs_param_present[5]) { /* If the harmonic trapping term is present, a term will be required for every site */
823  for (j = 0; j < L; j++) {
824  osparam.vals[numos + j*os.sz].re = bs_params[loop].re/2.0; /* Copy over parameter */
825  osparam.vals[numos + j*os.sz].im = bs_params[loop].im/2.0; /* Copy over parameter */
826  }
827  } else {
828  osparam.vals[numos].re = bs_params[loop].re/2.0; /* Copy over parameter */
829  osparam.vals[numos].im = bs_params[loop].im/2.0; /* Copy over parameter */
830  }
831  numos++;
832  break;
833  case 5: /* Term mub - chemical potential */
834  /* ---- Set the nodes to make the term ---- */
835  os.vals[numos] = tntNodeCopy(n_op); /* Operator acting on site */
836  if (bs_param_present[5]) { /* If the harmonic trapping term is present, a term will be required for every site */
837  for (j = 0; j < L; j++) {
838  osparam.vals[numos + j*os.sz].re = -bs_params[loop].re; /* Copy over parameter */
839  osparam.vals[numos + j*os.sz].im = -bs_params[loop].im; /* Copy over parameter */
840  }
841  } else {
842  osparam.vals[numos].re = -bs_params[loop].re; /* Copy over parameter */
843  osparam.vals[numos].im = -bs_params[loop].im; /* Copy over parameter */
844  }
845  numos++;
846  break;
847  case 6: /* Term Eharm: harmonic trapping term */
848  /* ---- Set the nodes to make the term ---- */
849  os.vals[numos] = tntNodeCopy(n_op); /* Operator acting on site */
850  for (j = 0; j < L; j++) {
851  /* Copy over parameter and covert to harmonic trap magnitude on that site*/
852  osparam.vals[numos + j*os.sz].re = 0.5*bs_params[loop].re*(j - jc)*(j - jc);
853  osparam.vals[numos + j*os.sz].im = 0.5*bs_params[loop].im*(j - jc)*(j - jc);
854  }
855  numos++;
856  break;
857  case 7: /* Term coher-driv:coherent driving term */
858  /* ---- Set the nodes to make the term ---- */
859  os.vals[numos] = tntNodeCopy(b_op); /* Operator acting on site */
860  os.vals[numos+1] = tntNodeCopy(bdag_op); /* Operator acting on site */
861 
862  if (bs_param_present[5]) { /* If the harmonic trapping term is present, a term will be required for every site */
863  for (j = 0; j < L; j++) {
864  osparam.vals[numos + j*os.sz] = bs_params[loop]; /* Copy over parameter */
865  osparam.vals[numos + 1 + j*os.sz].re = bs_params[loop].re; /* Copy over cc parameter */
866  osparam.vals[numos + 1 + j*os.sz].im = -bs_params[loop].im; /* Copy over cc parameter */
867  }
868  } else {
869  osparam.vals[numos] = bs_params[loop]; /* Copy over parameter */
870  osparam.vals[numos+1].re = bs_params[loop].re; /* Copy over cc parameter */
871  osparam.vals[numos+1].im = -bs_params[loop].im; /* Copy over cc parameter */
872  }
873  numos += 2;
874  break;
875  default:
876  break;
877  }
878  }
879  }
880  }
881  /* Free the boson operators */
882  tntNodeFree(&b_op);
883  tntNodeFree(&bdag_op);
884  tntNodeFree(&n_op);
885  tntNodeFree(&os_int);
886  tntNodeFree(&eye);
887 
888  } else if ((!no_operators) && (0 == tntSysTypeGet())) {
889  tntErrorPrint("No operators can be created as the system type has not been set"); /* NO_COVERAGE */
890  }
891 
892 
893  if (NULL != mpo) {
894  /* Call the function that makes MPOs */
895  *mpo = tntMpsCreateMpo(L, &nnL, &nnR, &nnparam, &os, &osparam);
896  }
897 
898  if (NULL != prop) {
899  /* Contruct the time step */
900  dtc.re = dt;
901  dtc.im = idt;
902 
903  /* Call the function that makes the propagators */
904  *prop = tntMpsCreatePropST2sc(L, dtc, &nnL, &nnR, &nnparam, &os, &osparam);
905  }
906 
907  /* Create the starting wave function */
908  if (MPStype) {
909  if (1 == MPStype) {
910  *mps = tntMpsCreateRandom(L,D);
911  } else if (2 == MPStype) {
912  *mps = tntMpsCreateConfig(L,config);
913  } else if (3 == MPStype) {
914  *mps = tntMpsCreateSymmRandom(L,&Q);
915  } else if (4 == MPStype) {
916  *mps = tntMpsCreateEyeMpo(L,NULL);
917  }
918  }
919 
920  if (!no_operators) {
921  /* Free the operator arrays if they have been created */
922  if (nnL.sz) {
923  tntNodeArrayFree(&nnL);
924  tntNodeArrayFree(&nnR);
925  tntComplexArrayFree(&nnparam);
926  }
927 
928  if (os.sz) {
929  tntNodeArrayFree(&os);
930  tntComplexArrayFree(&osparam);
931  }
932 
933  }
934 
935  free(argv_cp);
936 
937  /* Now deal with expectation values - separate function used for this */
938  tntMpsProcessExpecOptions(argc, argv, ExOp);
939 
940 }
941 
942 
tntNode tntNodeCopy(tntNode A)
Definition: tntNodeUtil.c:304
#define tntMpsExOp
Definition: tntMps.h:75
void tntNodeArrayFree(tntNodeArray *arr)
Definition: tntArray.c:668
void tntMpsProcessSysOptions(int argc, char **argv, tntNetwork *mps, tntNetwork *mpo, tntNetwork *prop, tntMpsExOp *ExOp)
int tntSysTypeGet(void)
Definition: tntSys.c:376
void tntMpsCreateSpinOp(unsigned TwoS, tntNode *Sx, tntNode *Sy, tntNode *Sz, tntNode *Sp, tntNode *Sm, tntNode *eye)
void tntMpsProcessExpecOptions(int argc, char **argv, tntMpsExOp *Op)
tntNodeArray tntNodeArrayAlloc(unsigned num_elems)
Definition: tntArray.c:428
void tntSysTypeSet(int sysnum)
Definition: tntSys.c:361
tntNetwork tntMpsCreateMpo(unsigned L, tntNodeArray *nnl, tntNodeArray *nnr, tntComplexArray *nnparam, tntNodeArray *os, tntComplexArray *osparam, unsigned ignoreQN)
void tntMpsCreateBosonOp(unsigned Nmax, tntNode *b, tntNode *bd, tntNode *n, tntNode *os_int, tntNode *eye)
void tntNodeFree(tntNode *A)
Definition: tntNodeUtil.c:275
tntComplexArray tntComplexArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntArray.c:392
tntNetwork tntMpsCreateConfig(unsigned L, const char *config)
Definition: tnt.h:65
tntNetwork tntMpsCreateRandom(unsigned L, unsigned D)
void tntSymmTypeSet(const char *symmTypeStr, unsigned num_qn)
Definition: tntSys.c:439
tntComplex tntAtoC(const char *str)
Definition: tntUtil.c:123
void tntSysBasisOpSet(tntNode basisOp)
Definition: tntSys.c:389
void tntComplexArrayFree(tntComplexArray *arr)
Definition: tntArray.c:650
double re
Definition: tnt.h:66
int tntSymmTypeGet(void)
Definition: tntSys.c:498
int tntPrintf(const char *format,...)
Definition: tntPrint.c:77
double im
Definition: tnt.h:67
tntNetwork tntMpsCreateEyeMpo(unsigned L, tntNode basisOp)
tntNetwork tntMpsCreatePropST2sc(unsigned L, tntComplex dtc, tntNodeArray *nnL, tntNodeArray *nnR, tntComplexArray *nnparam, tntNodeArray *os, tntComplexArray *osparam)
tntNetwork tntMpsCreateSymmRandom(unsigned L, int *qn)