Tensor Network Theory Library  Beta release 1.0
A library of routines for performing TNT-based operations
Command line

Functions

void tntMpsProcessTebdCLOptions (int argc, char **argv, unsigned *numsteps, unsigned *tbigstep, int *update_chi, double *emax, unsigned *delta_chi, unsigned *max_chi)
 
void tntMpsProcessDmrgCLOptions (int argc, char **argv, unsigned *max_it, double *precision)
 
void tntMpsProcessExpecOptions (int argc, char **argv, tntMpsExOp *ExOp)
 
void tntMpsProcessSysOptions (int argc, char **argv, tntNetwork *mps, tntNetwork *mpo, tntNetwork *prop, tntMpsExOp *ExOp)
 

Detailed Description

These functions are useful for allowing various flags on the command line to be accepted and processed for setting system parameters useful for simulations.

Function Documentation

void tntMpsProcessDmrgCLOptions ( int  argc,
char **  argv,
unsigned *  max_it,
double *  precision 
)

Function that processes additional commonly required input command line arguments for the MPS DMRG algorithm. For general input command line arguments use tntProcessCLOptions() in addition to this, and for input options that initialise the system use tntMpsProcessInitOptions().

Command line options:

Long nameShort nameDescription
--help-dmrg <n> Print this help
--max-it <n> Maximum number of iterations. If not given default 50 is set.
--trunc-type <str> How to calculate the truncation error using the discarded singular values. See tntTruncType() for allowed options.
--precision <val> When energy difference between successive DMRG iterations is smaller than this value, the calculation ends. If not given default 1e-4 is used.

Note that, with the exception of --trunc-type which calls the function tntTruncType() to set the truncation type, this function just sets the required variables. The way these variables will be used is implementation dependent. If any of the variables are not required or should not be allowed to be set from the command line for your implementation, simply pass the NULL pointer in place of that argument.

Parameters
argcInput argument count.
argvInput arguments.
max_itMaximum number of iterations.
precisionMaximum difference in energy between iterations for calculation to complete
Examples:
dmrg.c.
237 {
238 
239  static struct option long_options[] = {
240  {"max-it", required_argument, 0, 1},
241  {"precision", required_argument, 0, 2},
242  {"trunc-type", required_argument, 0, 3},
243  {"help-dmrg", no_argument, 0, 4},
244  {0, 0, 0, 0}
245  };
246  int option_index;
247  int c;
248  char **argv_cp;
249 
250  /* Make a copy of argv - getopt_long will permute the arguments, and they may be required by another processing function */
251  argv_cp = (char **) malloc(argc*sizeof(char *));
252  for (c = 0; c<argc; c++) {
253  argv_cp[c] = argv[c];
254  }
255 
256  /* Set optind to 1 in case it has already been called previously */
257  optind = 1;
258 
259  /* Set opterr to zero to prevent error message for unknown erguments, since they may be being treated by another process options function */
260  opterr = 0;
261 
262  /* Set variables to a preset value, so can check if it is unset at the end */
263  if (max_it != NULL) *max_it = 0;
264  if (precision != NULL) *precision = -100.0;
265 
266  /* Loop over all input arguments */
267  while ((c = getopt_long (argc, argv_cp, "", long_options, &option_index)) != -1) {
268  switch (c) {
269  case 1 :
270  if (max_it != NULL) {
271  *max_it = atoi(optarg);
272  printf("Maximum number of iterations <<<%d>>>\n", *max_it);
273  }
274  break;
275  case 2 :
276  if (precision != NULL) {
277  *precision = atof(optarg);
278  printf("Maximum difference in energy between DMRG iterations <<<%7.2e>>>\n", *precision);
279  }
280  break;
281  case 3 :
282  tntTruncType(optarg);
283  printf("Truncation errors will be calculated using %s.\n", optarg);
284  break;
285  case 4 :
286  printf("Options for TNT DMRG routines\n");
287  printf("\n");
288  printf("%s <parameters>:\n", argv[0]);
289  printf("Input parameters:\n");
290  printf("\n");
291  printf("--help-dmrg | : Print this help.\n");
292  printf("--max-it <n> | : Maximum number of iterations. If not given default 50 is set. \n");
293  printf("--precision <val> | : When energy difference between successive DMRG iterations is smaller than this value, the calculation ends. \n");
294  printf(" If not given default 1e-4 is used. \n");
295  printf("--trunc-type <str> | : How to calculate the truncation error using the discarded singular values.\n");
296  exit(0);
297  default: /* '?' Do nothing since argument may be treated by another function */
298  break;
299  }
300 
301  }
302 
303 
304  /* Set unset variables to default values */
305  if ((max_it != NULL ) && (0 == *max_it)) {
306  *max_it = 50;
307  printf("max_it: Using default <<<%d>>>\n", *max_it);
308  }
309  if ((precision != NULL) && (*precision < -99.0)) {
310  *precision = 1e-4;
311  printf("precision: using default <<<%7.2e>>>\n", *precision);
312  }
313 
314  free(argv_cp);
315 
316  return;
317 
318 }
void tntTruncType(const char *funcname)
Definition: tntMisc_funcs.c:497
void tntMpsProcessExpecOptions ( int  argc,
char **  argv,
tntMpsExOp ExOp 
)

Function that processes input command line arguments to create the variables needed for calculating expecation values for spin-half or boson simulations using matrix product states (MPS). Using the options given on the command line, the function can create arrays containing operators for:

  1. On-site expectation values
  2. Nearest neighbour expectation values
  3. All pairs of sites for two-site operators
  4. Centre to all other sites for two site operators.

It stores the operators required in the tntMpsExOp structure.

The system type e.g. 'spin' or 'boson' must previously have been set before calling this function.

If the flag for a given expecation value is given, then the corresponding operators are placed in the structure. This structure can then be passed to a relevant tier 2 function e.g. tntMpsExpecOutput() to calculate expectation values when they are required.

Onsite expectation values

Setting these flags will indicate that the value \(\langle\psi|\hat{o}_j|\psi\rangle\) should be calculated for every site \(j\) i.e. \(L\) values will be calculated.

For bosonic systems the allowed flags are:

Long nameDescription
--Ex1N Calculate \(\langle\psi|\hat{n}|\psi\rangle\)
--Ex1Int Calculate \(\langle\psi|\hat{n}(\hat{n}-1)/2|\psi\rangle\)
--Ex1Nsq Calculate \(\langle\psi|\hat{n}^2|\psi\rangle\)
--Ex1Q1 Calculate \(\langle\psi|\hat{b} + \hat{b}^{\dagger}|\psi\rangle\)
--Ex1Q2 Calculate \(\langle\psi|\mathrm{i}(\hat{b} - \hat{b}^{\dagger})|\psi\rangle\)

For spin systems the allowed flags are:

Long nameDescription
--Ex1Sx Calculate \(\langle\psi|\hat{S}^{x}|\psi\rangle\)
--Ex1Sy Calculate \(\langle\psi|\hat{S}^{y}|\psi\rangle\)
--Ex1Sz Calculate \(\langle\psi|\hat{S}^{z}|\psi\rangle\)

Two-site expectation values

Setting these flags will indicate that the value \(\langle\psi|\hat{o}^{[1]}_{i}\hat{o}^{[2]}_{j}|\psi\rangle\) should be calculated for pairs of sites \(i\) and \(j\). There are three different choices for the pairs of sites that will be chosen, each with a two-character code set as the argument of the flag (where an '=' must be used) i.e. Ex2<opname>=**.

  • For nearest-neighbour expectation values it will be calculated for \(i=0\) to \(L-2\) and \(j=i+1\) i.e. \(L-1\) values will be calculated. Use the code nn to specify this e.g. --Ex2SzSz=nn.

  • For centre-site expectation values it will be calculated for \(i=i_c\) where \(i_c\) is the central site \(\mathrm{floor}(L/2)\) and \(j=0\) to \(L-1\) i.e. \(L\) values will be calculated. Use the code cs to specify this e.g. --Ex2SzSz=cs.

  • For all pairs of expectation values it will be calculated for both \(i=0\) to \(L-1\) and \(j=0\) to \(L-1\) i.e. \(L^2\) values will be calculated. Use the code ap to specify this e.g. --Ex2SzSz=ap.

If no two-character code is given, or the two-character code is not recognised, then nearest-neighbour pairs will be chosen.

For bosonic systems the allowed flags are:

Long nameDescription
--Ex2bdagb Calculate \(\langle\psi|\hat{b}^{\dagger}_{i}\hat{b}_{j}|\psi\rangle\).
--Ex2den Calculate \(\langle\psi|\hat{n}_{i}\hat{n}_{j}|\psi\rangle\).

For spin systems the allowed flags are:

Long nameDescription
--Ex2SxSx Calculate \(\langle\psi|\hat{S}^{x}_{i}\hat{S}^{x}_{j}|\psi\rangle\).
--Ex2SySy Calculate \(\langle\psi|\hat{S}^{y}_{i}\hat{S}^{y}_{j}|\psi\rangle\).
--Ex2SzSz Calculate \(\langle\psi|\hat{S}^{z}_{i}\hat{S}^{z}_{j}|\psi\rangle\).
--Ex2SxSy Calculate \(\langle\psi|\hat{S}^{x}_{i}\hat{S}^{y}_{j}|\psi\rangle\).
--Ex2SySx Calculate \(\langle\psi|\hat{S}^{y}_{i}\hat{S}^{x}_{j}|\psi\rangle\).
--Ex2SySz Calculate \(\langle\psi|\hat{S}^{y}_{i}\hat{S}^{z}_{j}|\psi\rangle\).
--Ex2SzSy Calculate \(\langle\psi|\hat{S}^{z}_{i}\hat{S}^{y}_{j}|\psi\rangle\).
--Ex2SzSx Calculate \(\langle\psi|\hat{S}^{z}_{i}\hat{S}^{x}_{j}|\psi\rangle\).
--Ex2SxSz Calculate \(\langle\psi|\hat{S}^{x}_{i}\hat{S}^{z}_{j}|\psi\rangle\).
--Ex2SpSm Calculate \(\langle\psi|\hat{S}^{+}_{i}\hat{S}^{-}_{j}|\psi\rangle\).
Returns
No return value.
Parameters
argcInput argument count.
argvInput arguments.
ExOpStructure that contains nodes that represent operators for finding expectation values. Unitialised on entry.
204 {
205 
206  int option_index;
207  int c; /* Variable that contains the returned option code from getopt_long */
208  char **argv_cp; /* Used to hold a copy of the arguments - this will be used by the function instead of the original arguments */
209  unsigned numos = 0; /* Number of onsite expectation values to calculate */
210  unsigned numts[TNT_NUM_TSTYPES] = {0,0,0}; /* Number of two-site expectation values of each type to calculate */
211  char ts_labels[TNT_NUM_TSTYPES][TNT_STRLEN] = {"nearest-neighbour","centre-site","all"};
212  char ts_codes[TNT_NUM_TSTYPES][TNT_TSTYPE_STRLEN] = {"nn","cs","ap"};
213  tntNodeArray *ts_narrays[TNT_NUM_TSTYPES]; /* Array of pointers to tntNodeArray arguments to be populated */
214  tntStringArray *ts_sarrays[TNT_NUM_TSTYPES]; /* Array of pointers to tntStringArray arguments to be populated */
215  int loop, tstype; /* loop counter for terms, type of two-site term */
216 
217  int numosterms; /* Number of onsite terms */
218  int numtsterms; /* Number of two-site terms */
219  struct option *os_options; /* Options available for onsite terms. Must all have no_argument set as argument type. */
220  struct option *ts_options; /* Options available for two-site terms. Must all have optional_argument set as argument type. */
221  tntStringArray osop_strs, tsop_strs; /* String for each one-site and two-site term */
222  tntNodeArray osop_ops, tsop_ops[2]; /* Operators for each one-site and each site of each two-site term */
223  tntIntArray osop_present, tsop_present; /* Specifies whether each one-site term or two-site term is present */
224 
225  printf("----------------------------------------------------------\n");
226  printf("The following expectation values will be calculated in the simulation: \n");
227 
228  /* Assign the expectation value operator arrays to the arrays for two-site expectation values */
229  ts_narrays[TNT_TSTYPE_NN] = &(ExOp->ExOp2nn);
230  ts_narrays[TNT_TSTYPE_CS] = &(ExOp->ExOp2cs);
231  ts_narrays[TNT_TSTYPE_AP] = &(ExOp->ExOp2ap);
232 
233  ts_sarrays[TNT_TSTYPE_NN] = &(ExOp->LbOp2nn);
234  ts_sarrays[TNT_TSTYPE_CS] = &(ExOp->LbOp2cs);
235  ts_sarrays[TNT_TSTYPE_AP] = &(ExOp->LbOp2ap);
236 
237  /* Set allowed options for each type of system */
238  if (TNT_MPS_SPIN == tntSysTypeGet()) { /* Do this bit if the system is a spin system - all bosonic command line arguments will be ignored. */
239 
240 
241  /* Assign the arrays that will be used for dealing with the command line options */
242 
243  /* Long options for spin onsite terms */
244  static struct option sp_os_options[] = {
245  {"Ex1Sx", no_argument, 0, TNT_SP_OS_SX},
246  {"Ex1Sy", no_argument, 0, TNT_SP_OS_SY},
247  {"Ex1Sz", no_argument, 0, TNT_SP_OS_SZ},
248  {0, 0, 0, 0}
249  };
250 
251  /* Long options for spin two-site terms */
252  static struct option sp_ts_options[] = {
253  {"Ex2SxSx", optional_argument, 0, TNT_SP_TS_SXSX},
254  {"Ex2SySy", optional_argument, 0, TNT_SP_TS_SYSY},
255  {"Ex2SzSz", optional_argument, 0, TNT_SP_TS_SZSZ},
256  {"Ex2SxSy", optional_argument, 0, TNT_SP_TS_SXSY},
257  {"Ex2SySz", optional_argument, 0, TNT_SP_TS_SYSZ},
258  {"Ex2SzSx", optional_argument, 0, TNT_SP_TS_SZSX},
259  {"Ex2SySx", optional_argument, 0, TNT_SP_TS_SYSX},
260  {"Ex2SzSy", optional_argument, 0, TNT_SP_TS_SZSY},
261  {"Ex2SxSz", optional_argument, 0, TNT_SP_TS_SXSZ},
262  {"Ex2SpSm", optional_argument, 0, TNT_SP_TS_SPSM},
263  {0, 0, 0, 0}
264  };
265 
266  /* Create all the operators needed for processing the command line options */
267  tntNode Sx, Sy, Sz, Sp, Sm; /* Spin operators */
268  unsigned TwoS = tntNodeGetLegDim(tntSysBasisOpGet(),TNT_MPS_D) - 1; /* Physical dimension of the spin system */
269 
270  /* Create the spin operators */
271  tntMpsCreateSpinOp(TwoS, &Sx, &Sy, &Sz, &Sp, &Sm, NULL);
272 
273  /* Set variables that are system dependent */
274  numosterms = TNT_SP_NUMOSTERMS; /* Number of onsite terms */
275  numtsterms = TNT_SP_NUMTSTERMS; /* Number of two-site terms */
276  os_options = sp_os_options; /* Options available for onsite terms. Must all have no_argument set as argument type. */
277  ts_options = sp_ts_options; /* Options available for two-site terms. Must all have optional_argument set as argument type. */
278 
279  /* String for each one-site term */
280  osop_strs = tntStringArrayAlloc(numosterms);
281  strcpy(osop_strs.vals[TNT_SP_OS_SX], "Sx");
282  strcpy(osop_strs.vals[TNT_SP_OS_SY], "Sy");
283  strcpy(osop_strs.vals[TNT_SP_OS_SZ], "Sz");
284 
285  /* Operators for each one-site term */
286  osop_ops = tntNodeArrayAlloc(numosterms);
287  osop_ops.vals[TNT_SP_OS_SX] = tntNodeCopy(Sx,0);
288  osop_ops.vals[TNT_SP_OS_SY] = tntNodeCopy(Sy,0);
289  osop_ops.vals[TNT_SP_OS_SZ] = tntNodeCopy(Sz,0);
290 
291  /* String for each two-site term */
292  tsop_strs = tntStringArrayAlloc(numtsterms);
293  strcpy(tsop_strs.vals[TNT_SP_TS_SXSX], "SxSx");
294  strcpy(tsop_strs.vals[TNT_SP_TS_SYSY], "SySy");
295  strcpy(tsop_strs.vals[TNT_SP_TS_SZSZ], "SzSz");
296  strcpy(tsop_strs.vals[TNT_SP_TS_SXSY], "SxSy");
297  strcpy(tsop_strs.vals[TNT_SP_TS_SYSZ], "SySz");
298  strcpy(tsop_strs.vals[TNT_SP_TS_SZSX], "SzSx");
299  strcpy(tsop_strs.vals[TNT_SP_TS_SYSX], "SySx");
300  strcpy(tsop_strs.vals[TNT_SP_TS_SZSY], "SzSy");
301  strcpy(tsop_strs.vals[TNT_SP_TS_SXSZ], "SxSz");
302  strcpy(tsop_strs.vals[TNT_SP_TS_SPSM], "SpSm");
303 
304  /* Operators for each site of two-site term */
305  tsop_ops[0] = tntNodeArrayAlloc(numtsterms); tsop_ops[1] = tntNodeArrayAlloc(numtsterms);
306  tsop_ops[0].vals[TNT_SP_TS_SXSX] = tntNodeCopy(Sx,0); tsop_ops[1].vals[TNT_SP_TS_SXSX] = tntNodeCopy(Sx,0);
307  tsop_ops[0].vals[TNT_SP_TS_SYSY] = tntNodeCopy(Sy,0); tsop_ops[1].vals[TNT_SP_TS_SYSY] = tntNodeCopy(Sy,0);
308  tsop_ops[0].vals[TNT_SP_TS_SZSZ] = tntNodeCopy(Sz,0); tsop_ops[1].vals[TNT_SP_TS_SZSZ] = tntNodeCopy(Sz,0);
309  tsop_ops[0].vals[TNT_SP_TS_SXSY] = tntNodeCopy(Sx,0); tsop_ops[1].vals[TNT_SP_TS_SXSY] = tntNodeCopy(Sy,0);
310  tsop_ops[0].vals[TNT_SP_TS_SYSZ] = tntNodeCopy(Sy,0); tsop_ops[1].vals[TNT_SP_TS_SYSZ] = tntNodeCopy(Sz,0);
311  tsop_ops[0].vals[TNT_SP_TS_SZSX] = tntNodeCopy(Sz,0); tsop_ops[1].vals[TNT_SP_TS_SZSX] = tntNodeCopy(Sx,0);
312  tsop_ops[0].vals[TNT_SP_TS_SYSX] = tntNodeCopy(Sy,0); tsop_ops[1].vals[TNT_SP_TS_SYSX] = tntNodeCopy(Sx,0);
313  tsop_ops[0].vals[TNT_SP_TS_SZSY] = tntNodeCopy(Sz,0); tsop_ops[1].vals[TNT_SP_TS_SZSY] = tntNodeCopy(Sy,0);
314  tsop_ops[0].vals[TNT_SP_TS_SXSZ] = tntNodeCopy(Sx,0); tsop_ops[1].vals[TNT_SP_TS_SXSZ] = tntNodeCopy(Sz,0);
315  tsop_ops[0].vals[TNT_SP_TS_SPSM] = tntNodeCopy(Sp,0); tsop_ops[1].vals[TNT_SP_TS_SPSM] = tntNodeCopy(Sm,0);
316 
317  /* Free the spin operators */
318  tntNodeFree(&Sx);
319  tntNodeFree(&Sy);
320  tntNodeFree(&Sz);
321  tntNodeFree(&Sp);
322  tntNodeFree(&Sm);
323 
324  } else if (TNT_MPS_BOSON == tntSysTypeGet()) { /* Do this bit if the system is a bosonic system - all spin command line arguments will be ignored. */
325 
326  /* Assign the arrays that will be used for dealing with the command line options */
327 
328  /* Long options for boson onsite terms */
329  static struct option bs_os_options[] = { /* Long options for bosonic onsite terms */
330  {"Ex1N", no_argument, 0, TNT_BS_OS_N},
331  {"Ex1Int", no_argument, 0, TNT_BS_OS_INT},
332  {"Ex1Nsq", no_argument, 0, TNT_BS_OS_NSQ},
333  {"Ex1Q1", no_argument, 0, TNT_BS_OS_Q1},
334  {"Ex1Q2", no_argument, 0, TNT_BS_OS_Q2},
335  {0, 0, 0, 0}
336  };
337 
338  /* Long options for boson two-site terms */
339  static struct option bs_ts_options[] = { /* Long options for nearest-neighbour bosonic two-site terms */
340  {"Ex2bdagb", optional_argument, 0, TNT_BS_TS_BDB},
341  {"Ex2denden", optional_argument, 0, TNT_BS_TS_DD},
342  {0, 0, 0, 0}
343  };
344 
345 
346  /* Create all the operators needed for processing the command line options */
347  tntNode b, bdag, n, os_int, nsq, ncp, Q1, Q2; /* Boson operators */
348  unsigned nmax = tntNodeGetLegDim(tntSysBasisOpGet(),TNT_MPS_D) - 1; /* Maximum number of bosons on each site */
349  tntComplex imag = {0,1.0}; /* i */
350 
351  /* Create the boson operators */
352  tntMpsCreateBosonOp(nmax, &b, &bdag, &n, &os_int, NULL);
353 
354  /* ---- Create the nodes to make the term n^2 ---- */
355  nsq = tntNodeCopy(n,0);
356  ncp = tntNodeCopy(n,0);
358  nsq = tntNodeContract(ncp,nsq,NULL,NULL);
359 
360  /* ---- Create the nodes to make the term b + bdag ---- */
361  Q1 = tntNodeCopy(b,0);
362  tntNodeAdd(Q1,bdag);
363 
364  /* ---- Create the nodes to make the term i(b - bdag) ---- */
365  Q2 = tntNodeCopy(bdag,0);
366  tntNodeScaleReal(Q2,-1.0);
367  tntNodeAdd(Q2,b);
368  tntNodeScaleComplex(Q2,imag);
369 
370  /* Set variables that are system dependent */
371  numosterms = TNT_BS_NUMOSTERMS; /* Number of onsite terms */
372  numtsterms = TNT_BS_NUMTSTERMS; /* Number of two-site terms */
373  os_options = bs_os_options; /* Options available for onsite terms. Must all have no_argument set as argument type. */
374  ts_options = bs_ts_options; /* Options available for two-site terms. Must all have optional_argument set as argument type. */
375 
376  /* String for each one-site term */
377  osop_strs = tntStringArrayAlloc(numosterms);
378  strcpy(osop_strs.vals[TNT_BS_OS_N], "density");
379  strcpy(osop_strs.vals[TNT_BS_OS_INT], "interaction");
380  strcpy(osop_strs.vals[TNT_BS_OS_NSQ], "n_squared");
381  strcpy(osop_strs.vals[TNT_BS_OS_Q1], "Q1");
382  strcpy(osop_strs.vals[TNT_BS_OS_Q2], "Q2");
383 
384  /* Operators for each one-site term */
385  osop_ops = tntNodeArrayAlloc(numosterms);
386  osop_ops.vals[TNT_BS_OS_N] = tntNodeCopy(n,0);
387  osop_ops.vals[TNT_BS_OS_INT] = tntNodeCopy(os_int,0);
388  osop_ops.vals[TNT_BS_OS_NSQ] = tntNodeCopy(nsq,0);
389  osop_ops.vals[TNT_BS_OS_Q1] = tntNodeCopy(Q1,0);
390  osop_ops.vals[TNT_BS_OS_Q2] = tntNodeCopy(Q2,0);
391 
392  /* String for each two-site term */
393  tsop_strs = tntStringArrayAlloc(numtsterms);
394  strcpy(tsop_strs.vals[TNT_BS_TS_BDB], "bdagb");
395  strcpy(tsop_strs.vals[TNT_BS_TS_DD], "denden");
396 
397  /* Operators for each site of two-site term */
398  tsop_ops[0] = tntNodeArrayAlloc(numtsterms); tsop_ops[1] = tntNodeArrayAlloc(numtsterms);
399  tsop_ops[0].vals[TNT_BS_TS_BDB] = tntNodeCopy(bdag,0); tsop_ops[1].vals[TNT_BS_TS_BDB] = tntNodeCopy(b,0);
400  tsop_ops[0].vals[TNT_BS_TS_DD] = tntNodeCopy(n,0); tsop_ops[1].vals[TNT_BS_TS_DD] = tntNodeCopy(n,0);
401 
402  /* Free the boson operators */
403  tntNodeFree(&b);
404  tntNodeFree(&bdag);
405  tntNodeFree(&n);
406  tntNodeFree(&os_int);
407  tntNodeFree(&nsq);
408  tntNodeFree(&Q1);
409  tntNodeFree(&Q2);
410 
411  } else if (0 == tntSysTypeGet()) {
412  printf("Error setting up operators for expectation values: System type has not been set yet."); /* NO_COVERAGE */
413  exit(1); /* NO_COVERAGE */
414  } else { /* NO_COVERAGE */
415  printf("Error setting up operators for expectation values: System type %d is not known.", tntSysTypeGet()); /* NO_COVERAGE */
416  exit(1); /* NO_COVERAGE */
417  } /* NO_COVERAGE */
418 
419  /* -------------------------------------------------------------- */
420  /* ---- Now general part that doesn't depend on system type ----- */
421  /* -------------------------------------------------------------- */
422 
423 
424  /* Set opterr to zero to prevent error message for unknown erguments, since they may be being treated by another process options function */
425  opterr = 0;
426 
427  /* Allocate memory for arrays that determine whether or not a given term is present */
428  osop_present = tntIntArrayAlloc(numosterms);
429  tsop_present = tntIntArrayAlloc(numtsterms);
430 
431  /* Initialise arrays for whether terms are present - if term not present value -1 is assigned */
432  for (loop = 0; loop < numosterms; loop++)
433  osop_present.vals[loop] = -1;
434  for (loop = 0; loop < numtsterms; loop++)
435  tsop_present.vals[loop] = -1;
436 
437  /* --- Loop over the input arguments for onsite terms first --- */
438 
439  /* Set optind to 1 as option getting routine has already been called previously */
440  optind = 1;
441 
442  /* Make a copy of argv - getopt_long will permute the arguments, and they will be required by another processing function */
443  argv_cp = (char **) malloc(argc*sizeof(char *));
444  for (c = 0; c < argc; c++) {
445  argv_cp[c] = argv[c];
446  }
447 
448  /* Determine which options are present */
449  c = getopt_long(argc, argv_cp, "", os_options, &option_index);
450  while ((c = getopt_long(argc, argv_cp, "", os_options, &option_index)) != -1) {
451  /* Populate the array of flags */
452  if ((c <= numosterms)&&(-1 == osop_present.vals[c])) {
453  osop_present.vals[c] = 1;
454  printf(" - %s \n", osop_strs.vals[c]);
455  numos++;
456  }
457  }
458 
459  /* Allocate arrays to hold the terms and parameters that are present */
460  if (numos) {
461  ExOp->ExOpOs = tntNodeArrayAlloc(numos);
462  ExOp->LbOpOs = tntStringArrayAlloc(numos);
463  } else {
464  ExOp->ExOpOs.sz = 0;
465  ExOp->LbOpOs.sz = 0;
466  }
467 
468  /* reset the counter */
469  numos = 0;
470 
471  /* Loop through parameters and assign them to the array */
472  for (loop = 0; loop < numosterms; loop++) { /* Loop over all possible terms */
473  if (1 == osop_present.vals[loop]) { /* Check if this term is present */
474  ExOp->ExOpOs.vals[numos] = tntNodeCopy(osop_ops.vals[loop],0); /* Onsite operator */
475  strcpy(ExOp->LbOpOs.vals[numos], osop_strs.vals[loop]); /* Label for on-site operator */
476  numos++;
477  } /* end of check if term is present */
478  } /* End of loop over all possible onsite terms */
479 
480  /* -------------- Now deal with two-site terms ------------------ */
481 
482  /* Make a copy of argv - getopt_long will permute the arguments, and they may be required later */
483  for (c = 0; c<argc; c++)
484  argv_cp[c] = argv[c];
485 
486  /* Set optind to 1 as option getting routine has already been called previously */
487  optind = 1;
488 
489  /* Loop over options for these set of two-site terms, making a note if they are present, and increasing the count */
490  while ((c = getopt_long(argc, argv_cp, "", ts_options, &option_index)) != -1) {
491  /* Determine the type of two-site term. Default is nearest neighbour */
492  tstype = TNT_TSTYPE_NN;
493  if (optarg) {
494  for (tstype = 0; tstype < TNT_NUM_TSTYPES; tstype++) {
495  if (0 == strcmp(optarg,ts_codes[tstype]))
496  break; /* break if string is the same as a code - ensuring tstype takes the value of the term present */
497  }
498  }
499  /* Populate the array of flags */
500  if ((c <= numtsterms)&&(tsop_present.vals[c] != 1)) {
501  /* if the end of the loop was reached and no match found, set to default type, with warning for user */
502  if (TNT_NUM_TSTYPES == tstype) {
503  printf("Warning: type %s not known for two-site operators. Defaulting to nearest neighbour for expectation of %s.\n", optarg, tsop_strs.vals[c]);
504  tstype = TNT_TSTYPE_NN;
505  }
506  /* set the flag to the type of two-site term */
507  tsop_present.vals[c] = tstype;
508  /* print information for user */
509  printf(" - %s for %s pairs\n", tsop_strs.vals[c], ts_labels[tstype]);
510  /* increment the relevant counter */
511  numts[tstype]++;
512  }
513  } /* end of loop over options */
514 
515  /* Now the number has been counted, loop over type of two-site terms */
516  for (tstype = 0; tstype < TNT_NUM_TSTYPES; tstype++) {
517 
518  /* Allocate arrays to hold the terms and parameters that are present */
519  if (numts[tstype]) {
520  *(ts_narrays[tstype]) = tntNodeArrayAlloc(numts[tstype]*2); /* array twice size of number of terms, as two operators for each term */
521  *(ts_sarrays[tstype]) = tntStringArrayAlloc(numts[tstype]);
522  } else {
523  ts_narrays[tstype]->sz = 0;
524  ts_sarrays[tstype]->sz = 0;
525  }
526 
527  /* reset counter */
528  numts[tstype] = 0;
529  } /* End of loop over type of two-site terms */
530 
531  /* Loop through parameters and assign them to the correct array */
532  for (loop = 0; loop < numtsterms; loop++) { /* Loop over all possible terms */
533  if (-1 != tsop_present.vals[loop]) { /* Check if this term is present */
534  /* set the type of two-site term */
535  tstype = tsop_present.vals[loop];
536 
537  /* ---- Set the nodes to make the term ---- */
538  ts_narrays[tstype]->vals[2*numts[tstype]] = tntNodeCopy(tsop_ops[0].vals[loop],0); /* Operator acting on left */
539  ts_narrays[tstype]->vals[2*numts[tstype] + 1] = tntNodeCopy(tsop_ops[1].vals[loop],0); /* Operator acting on right */
540  strcpy(ts_sarrays[tstype]->vals[numts[tstype]], tsop_strs.vals[loop]); /* Label for two-site operator */
541  strcat(ts_sarrays[tstype]->vals[numts[tstype]], "_"); /* Append underscore */
542  strcat(ts_sarrays[tstype]->vals[numts[tstype]], ts_codes[tstype]); /* Append code for type of two-site operator */
543  numts[tstype]++;
544  } /* End of check if terms are present */
545  } /* End of loop over all possible two-site terms */
546 
547  printf("----------------------------------------------------------\n");
548 
549  /* Free the arrays */
550  free(argv_cp);
551  tntStringArrayFree(&osop_strs);
552  tntStringArrayFree(&tsop_strs);
553  tntNodeArrayFree(&osop_ops);
554  tntNodeArrayFree(tsop_ops);
555  tntNodeArrayFree(tsop_ops+1);
556  tntIntArrayFree(&osop_present);
557  tntIntArrayFree(&tsop_present);
558 
559 }
void tntMpsCreateSpinOp(unsigned TwoS, tntNode *Sx, tntNode *Sy, tntNode *Sz, tntNode *Sp, tntNode *Sm, tntNode *eye)
Definition: tntMpsCreateSpinOp.c:29
tntIntArray tntIntArrayAlloc(unsigned numrows, unsigned numcols)
Definition: tntDummy_funcs.c:70
tntStringArray LbOpOs
Array that contains strings for labelling the above operators.
Definition: tntMps.h:110
Simple 1D array for integer values.
Definition: tnt.h:148
tntStringArray tntStringArrayAlloc(unsigned num_elems)
Definition: tntArray_funcs.c:136
tntNodeArray ExOp2cs
Array that contains the centre-to-site expectation values to calculate. Will be a tntNodeArray of len...
Definition: tntMps.h:113
struct tnode * tntNode
Definition: tnt.h:118
void tntNodeScaleComplex(tntNode tn, tntComplex val)
Definition: tntNode_funcs.c:646
int tntSysTypeGet(void)
Definition: tntMisc_funcs.c:618
unsigned tntNodeGetLegDim(tntNode tn, unsigned legid)
Definition: tntNode_funcs.c:1069
int * vals
Pointer to the first element in the integer array.
Definition: tnt.h:149
#define TNT_MPS_BOSON
Definition: tntMps.h:99
tntNodeArray ExOp2ap
Array that contains the all pairs of two-site expectation values to calculate. Will be a tntNodeArray...
Definition: tntMps.h:115
unsigned sz
The number of elements in the array.
Definition: tnt.h:188
tntNodeArray tntNodeArrayAlloc(unsigned num_elems)
Definition: tntArray_funcs.c:110
tntNode tntNodeContract(tntNode tnA, tntNode tnB, int *legMapAC, int *legMapBC)
Definition: tntDummy_funcs.c:34
void tntMpsCreateBosonOp(unsigned n_max, tntNode *b, tntNode *bdag, tntNode *n, tntNode *os_int, tntNode *eye)
Definition: tntMpsCreateBosonOp.c:29
void tntNodeScaleReal(tntNode tn, double val)
Definition: tntNode_funcs.c:623
tntStringArray LbOp2cs
Array that contains strings for labelling the above operators.
Definition: tntMps.h:114
#define TNT_MPS_SPIN
Definition: tntMps.h:95
void tntStringArrayFree(tntStringArray *arr)
Definition: tntArray_funcs.c:233
void tntIntArrayFree(tntIntArray *arr)
Definition: tntArray_funcs.c:162
Simple 1D array for tntNodes.
Definition: tnt.h:178
void tntNodeArrayFree(tntNodeArray *arr)
Definition: tntArray_funcs.c:210
#define TNT_MPS_U
Definition: tntMps.h:71
tntNodeArray ExOp2nn
Array that contains the nearest-neighbour expectation values to calculate.
Definition: tntMps.h:111
tntStringArray LbOp2ap
Array that contains strings for labelling the above operators.
Definition: tntMps.h:116
Simple 1D array for strings (maximum length TNT_STRLEN).
Definition: tnt.h:186
char ** vals
Pointer to the first element in the array of strings.
Definition: tnt.h:187
tntNode * vals
Pointer to the first element in the tntNode array.
Definition: tnt.h:179
tntNodeArray ExOpOs
Array that contains the single site expectation values to calculate.
Definition: tntMps.h:109
void tntNodeFree(tntNode *tn)
Definition: tntNode_funcs.c:420
#define TNT_STRLEN
Definition: tnt.h:59
tntNode tntSysBasisOpGet(void)
Definition: tntMisc_funcs.c:647
tntNode tntNodeCopy(tntNode tn, int conj)
Definition: tntNode_funcs.c:486
tntStringArray LbOp2nn
Array that contains strings for labelling the above operators.
Definition: tntMps.h:112
void tntNodeJoin(tntNode tnA, int legA, tntNode tnB, int legB)
Definition: tntNode_funcs.c:462
#define TNT_MPS_D
Definition: tntMps.h:66
void tntNodeAdd(tntNode tnA, tntNode tnB)
Definition: tntNode_funcs.c:385
Use this type to define complex numbers. Any complex arithmetic needs to be done manually. For example to add two complex numbers A and B to give C write:
Definition: tnt.h:136
unsigned sz
The number of elements in the array.
Definition: tnt.h:180
void tntMpsProcessSysOptions ( int  argc,
char **  argv,
tntNetwork mps,
tntNetwork mpo,
tntNetwork prop,
tntMpsExOp ExOp 
)

Function that processes input command line arguments to create the system variables needed for spin-half or boson simulations using matrix product states (MPS). Using the options given on the command line, the function can create:

  • A network representing the initial MPS wave function.
  • A network representing the matrix product operator for the system Hamiltonian.
  • An array of nodes U each of which give a two-site propagator representing evolution under the system Hamiltonian.

If any of these is not required, then set the relevant argument to NULL.

The system to build can either be a spin-half system or a boson system, each with their own terms. It is not possible to mix terms for the different systems.

Spin system

The spin system Hamiltonian will have the following form:

\[ H = \sum_{j=0}^{L-2} J_{xx}\hat{S}^{x}_{j}\hat{S}^{x}_{j+1} + J_{yy}\hat{S}^{y}_{j}\hat{S}^{y}_{j+1} + J_{zz}\hat{S}^{z}_{j}\hat{S}^{z}_{j+1}+\sum_{j=0}^{L-1} B_{x}\hat{S}^{x}_{j} + B_{y}\hat{S}^{y}_{j} + B_{z}\hat{S}^{z}_{j} \]

\[ +\sum_{j=0}^{L-2} J_{xy}\hat{S}^{x}_{j}\hat{S}^{y}_{j+1} + J_{yx}\hat{S}^{y}_{j}\hat{S}^{x}_{j+1} + J_{yz}\hat{S}^{y}_{j}\hat{S}^{z}_{j+1} + J_{zy}\hat{S}^{z}_{j}\hat{S}^{y}_{j+1} + J_{zx}\hat{S}^{z}_{j}\hat{S}^{x}_{j+1} + J_{xz}\hat{S}^{x}_{j}\hat{S}^{z}_{j+1}. \]

Where the \(\hat{S}_{j}^{x,y,z}\) are the spin operators for site \(j\) in units where \(\hbar = 1\).

Boson system

The boson system Hamiltonian will have the following form:

\[ H = \sum_{j=0}^{L-2}\left[-J (\hat{b}_{j}^{\dagger}\hat{b}_{j+1} + \mathrm{h.c.}) + V\hat{n}_{j}\hat{n}_{j+1} + \Delta(\hat{b}_{j}\hat{b}_{j+1}+\hat{b}_{j}^{\dagger}\hat{b}_{j+1}^{\dagger}) \right] \]

\[ + \sum_{j=0}^{L-1} \left[\frac{U}{2} \hat{n}_{j}(\hat{n}_{j}-1) - \mu \hat{n}_{j} + \frac{1}{2}\epsilon_{\mathrm{harm}}\hat{n}_{j} (j - j_c)^2 + \Omega(\hat{b}_{j}+\hat{b}_{j}^{\dagger}) \right], \]

Where \(\hat{b}_j\) is the bosonic annihilation operator for site \(j\), and \(\hat{n}_j = \hat{b}_j^{\dagger}\hat{b}_j\).

Operators

For both of the above Hamiltonians, the terms are only included if the relevant parameter is given in the command line options - note that only real values are currently accepted.

The Hamiltonian is built as a network representing the matrix product operator (MPO) using function tntMpsCreateMpo(), which contains the sum of all the terms, and has the following form:

mpo_h.png

The physical dimension (i.e. the dimension of leg 3 and leg 4) will be \(2S+1\) for the spin system, and \(n_{\mathrm{max}} + 1\) for the boson system, where \(\hat{n}_{\mathrm{max}}\) is the maximum number of bosons allowed on each site. The internal dimension (i.e. the dimension of leg 1 and leg 2) will be equal to 2 + the number of nearest neighbour terms.

If one or both of the parameters --dt= \(\delta t\) and --idt = \(\delta t_i\) are given as a command line option the network representing the evolution in real or imaginary time can be built. The time steps are assumed to be scaled by \(\hbar\), and the propagators are formed into a network formed of two-site propagators arranged in a left-to-right then right-to-left staircase i.e. forming a second order Suzuki-Trotter decomposition. Each two-site propagator evolves its respective pair of sites for half the time-step i.e. each two site propagator will have the form

\[ U_{j,j+1} = \mathrm{e}^{-\mathrm{i}h_{j,j+1} \delta t/2}\times\mathrm{e}^{-h_{j,j+1} \delta t_i/2}, \]

and when the entire sequence of left-to-right and right-to-left staircases are applied they result in a total propagator for a timestep:

\[ U_{\mathrm{total}} = \mathrm{e}^{-\mathrm{i}H \delta t}\times\mathrm{e}^{-H \delta t_i}, \]

The resulting network is suitable for applying using the MPS library function tntMpsPropST2scProduct(), which multplies an MPS by a staircase propagator i.e. it evolves the MPS one time-step. If neither of dt or idt is specified prop cannot be constructed, and if the argument prop is not set to NULL then an error will result.

Start state

There are four options for the start state, and network representing the starting wave function will be created if the argument mps is a pointer to a tntNetwork variable, which should be unitialised on entry to the function, and if one of these options is given. If a wave function is not required (because you are loading it from an initialisation file), set the argument of the function to NULL.

Random state

This creates a random state having internal dimension \(D\) - this is useful as a starting state when using DMRG or imaginary time evolution to find the ground state. The initial wave function will be orthonormalised with the orthogonality centre on the first site. This wave function is then in the correct format for an initial left to right sweep using the TEBD or DMRG algorithms.

Quantum number conserving random state

This creates a random start state with a given quantum number Q, and sets the system properties to conserve this quantum number during any subsequent simulations. If the sytem is a bosonic system, then the quantum number Q representes the total number of bosons in the system. If the system is a spin system the the quantum number Q represents \(2\sum\hat{S}^{z}_{j}\).

Note that keeping track of symmetry information to conserve the quantum number causes any part of the MPO, two-site gates, and wave function that do not obey the U(1) symmetry to be discarded, and only the symmetry-invariant parts of the tensors are kept. This greatly speeds up the calculations.

Configuration state

This creates a starting configuration e.g. for a spin system:

\[ |\psi\rangle = |\uparrow\rangle_0\otimes|\uparrow\rangle_1\otimes|\downarrow\rangle_2\otimes|\downarrow\rangle_3\otimes|\uparrow\rangle_4\otimes\downarrow\rangle_5 \]

or for a boson system:

\[ |\psi\rangle = |0\rangle_0\otimes|1\rangle_1\otimes|1\rangle_2\otimes|2\rangle_3\otimes|1\rangle_4\otimes|1\rangle_5\otimes|0\rangle_6 \]

which will have internal dimension 1. This is useful as the start state of real time evolution. To select the first option, use --rand-state <n>, where <n> is an integer that gives the internal dimension. To select the second option, use --config-state=<n> where <n> is a number representing the configuration. The number is given by writing the configuration as a string of digits, with the maximum digit being ( \(d-1\)) e.g. 00110 for the spin example above (0 signifies spin up, 1 signifies spin down) and 0112110 for the boson example above. If the number of digits in the configuration is less than the number of sites specified, then all sites to the right will be filled with 0.

Quantum number conserving configuration state

This creates a configuration state in exactly the same way as above, however quantum numbers are added to the internal indices, and a flat is set in the system configuration, to ensure that the quantum number of the initial state is conserved in all subsequent simulations. If the sytem is a bosonic system, then the total number of bosons is conserved. If the system is a spin system then the total spin \(2\sum\hat{S}^{z}_{j}\) is conserved.

Note that keeping track of symmetry information to conserve the quantum number causes any part of the MPO, two-site gates, and wave function that do not obey the U(1) symmetry to be discarded, and only the symmetry-invariant parts of the tensors are kept. This greatly speeds up the calculations.

List of command line options

Long nameDescription
--system [spin | boson] System type. Required if operators are to be built.
--dt <val> Size of real time step (real double). Required to build U if –idt is not given. Default is 0.
--idt <val> Size of imaginary time step (real double). Required to build U if –dt is not given. Default is 0.
--rand-state <n> Create a random start state with internal dimension D given by the integer argument.
--qnum-rand-state <n> Create a random start state with total quantum number Q given by the integer argument, and conserve the total quantum number in subsequent simulations.
--config-state <str>Create a starting state having a configuration given by the supplied string.
--qnum-config-state <n> Create a starting state having a configuration given by the supplied string, and conserve the total quantum number of this configuration in subsequent simulations.
--length <n> Number of sites in the MPS. Required.

If --system=spin is chosen, then the following options are available:

Long nameDescription
--2S <n> Twice the spin \(S\) of each particle in the system. Default is 1 i.e. a spin half system.
--Jxx <val> Multiplier (real double) for \(\hat{S}^{x}_{j}\hat{S}^{x}_{j+1}\). Default is 0.
--Jyy <val> Multiplier (real double) for \(\hat{S}^{y}_{j}\hat{S}^{y}_{j+1}\). Default is 0.
--Jzz <val> Multiplier (real double) for \(\hat{S}^{z}_{j}\hat{S}^{z}_{j+1}\). Default is 0.
--Jxy <val> Multiplier (real double) for \(\hat{S}^{x}_{j}\hat{S}^{y}_{j+1}\). Default is 0.
--Jyx <val> Multiplier (real double) for \(\hat{S}^{y}_{j}\hat{S}^{x}_{j+1}\). Default is 0.
--Jyz <val> Multiplier (real double) for \(\hat{S}^{y}_{j}\hat{S}^{z}_{j+1}\). Default is 0.
--Jzy <val> Multiplier (real double) for \(\hat{S}^{z}_{j}\hat{S}^{y}_{j+1}\). Default is 0.
--Jzx <val> Multiplier (real double) for \(\hat{S}^{z}_{j}\hat{S}^{x}_{j+1}\). Default is 0.
--Jxz <val> Multiplier (real double) for \(\hat{S}^{x}_{j}\hat{S}^{z}_{j+1}\). Default is 0.
--Bx <val> Multiplier (real double) for \(\hat{S}^{x}_{j}\). Default is 0.
--By <val> Multiplier (real double) for \(\hat{S}^{y}_{j}\). Default is 0.
--Bz <val> Multiplier (real double) for \(\hat{S}^{z}_{j}\). Default is 0.

If --system=boson is chosen, then the following options are available:

Long nameDescription
--n-max <n> Maximum number of bosons allowed per site. Default is 1.
--Jb <val> Hopping energy i.e. multiplier (real double) for \(-(\hat{b}_{j}^{\dagger}\hat{b}_{j+1}+\mathrm{h.c.})\). Default is 0.
--Ub <val> On-site interaction energy i.e. multiplier (real double) for \(\hat{n}_{j}(\hat{n}_{j}-1)/2\). Default is 0.
--Vb <val> Nearest neighbour interaction energy i.e. multiplier (real double) for \(\hat{n}_{j}\hat{n}_{j+1}\). Default is 0.
--mub <val> Chemical potential i.e. multiplier (real double) for \(-\hat{n}_{j}\). Default is 0.
--E-harm <val> Energy scale for harmonic potential i.e. multiplier (real double) for \(\hat{n}_{j}(j-j_c)^2/2\). Default is 0.
--jc-harm <val> Centre of harmonic potential \(j_c\). Only applied when --E-harm is set. Default is centre of the system.
--coher-driv <val> Coherent driving term i.e. multiplier (real double) for \(\hat{b}_{j}+\hat{b}_{j}^{\dagger}\). Default is 0.
--param-driv <val> Parametric driving term i.e. multiplier (real double) for \(\hat{b}_{j}\hat{b}_{j+1}+\hat{b}_{j}^{\dagger}\hat{b}_{j+1}^{\dagger}\). Default is 0.

This function can also be used to create a structure that contains operators for: 1) On-site expectation values 2) Nearest neighbour expectation values 3) All pairs of sites for two-site operators 4) Centre to all other sites for two site operators.

It does this by calling tntMpsProcessExpecOptions() - see the documentation for this function for more information.

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

Function that processes additional commonly required input command line arguments for the MPS TEBD algorithm. For general input command line arguments use tntProcessCLOptions() in addition to this, and for input options that initialise the system use tntMpsProcessInitOptions().

Command line options:

Long nameShort nameDescription
--help-tebd <n> Print this help
--timesteps <n> -t <n> Number of timesteps to be performed.
--tbigstep <n> -b <n>Number of timesteps to be performed between each calculation of expectation values. If not given, then expectation values are calculated at the end of the routine only.
--trunc-type <str> How to calculate the truncation error using the discarded singular values. See tntTruncType() for allowed options.
--update-chi-u Update chi with increasing truncation error. Turned off by default.
--error-tol <val>Truncation error tolerance for increasing chi. Default is 1e-6.
--delta-chi <n> Amount by which to increase chi. Default value 50.
--max-chi <n> Value at which to stop increasing chi, even if the error exceeds error-tol. Default value 500.

Note that, with the exception of --trunc-type which calls the function tntTruncType() to set the truncation type, this function just sets the required variables. The way these variables will be used is implementation dependent. If any of the variables are not required or should not be allowed to be set from the command line for your implementation, simply pass the NULL pointer in place of that argument

Parameters
argcInput argument count.
argvInput arguments.
numstepsNumber of timesteps to run for.
tbigstepNumber of timesteps between each calculation of expectation values.
update_chiFlag for allowing chi to change.
emaxTolerance for truncation error for changing chi.
delta_chiAmount by which to increase chi.
max_chiValue at which to stop increasing chi.
Examples:
tebd.c.
68 {
69 
70 
71  static struct option long_options[] = {
72  {"timesteps", required_argument, 0, 't'},
73  {"tbigstep", required_argument, 0, 'b'},
74  {"trunc-type", required_argument, 0, 1},
75  {"update-chi", no_argument, 0, 'u'},
76  {"error-tol", required_argument, 0, 2},
77  {"delta-chi", required_argument, 0, 3},
78  {"max-chi", required_argument, 0, 4},
79  {"help-tebd", no_argument, 0, 5},
80  {0, 0, 0, 0}
81  };
82  int option_index;
83  int c;
84  char **argv_cp;
85 
86  /* Make a copy of argv - getopt_long will permute the arguments, and they may be required by another processing function */
87  argv_cp = (char **) malloc(argc*sizeof(char *));
88  for (c = 0; c<argc; c++) {
89  argv_cp[c] = argv[c];
90  }
91 
92  /* Set optind to 1 in case it has already been called previously */
93  optind = 1;
94 
95  /* Set opterr to zero to prevent error message for unknown erguments, since they may be being treated by another process options function */
96  opterr = 0;
97 
98  /* Set variables to a preset value, so can check if it is unset at the end */
99  if (numsteps != NULL) *numsteps = 0;
100  if (update_chi != NULL) *update_chi = 0;
101  if (tbigstep != NULL) *tbigstep = 0;
102  if (emax != NULL) *emax = -100.0;
103  if (delta_chi != NULL) *delta_chi = 0;
104  if (max_chi != NULL) *max_chi = 0;
105 
106  /* Loop over all input arguments */
107  while ((c = getopt_long (argc, argv_cp, "ht:b:u", long_options, &option_index)) != -1) {
108  switch (c) {
109  case 't':
110  if (numsteps != NULL) {
111  *numsteps = atoi(optarg);
112  printf("numsteps : <<<%d>>>\n", *numsteps);
113  }
114  break;
115  case 'b' :
116  if (tbigstep != NULL) {
117  *tbigstep = atoi(optarg);
118  printf("tbigstep: <<<%d>>>\n", *tbigstep);
119  }
120  break;
121  case 1 :
122  tntTruncType(optarg);
123  printf("Truncation errors will be calculated using %s.\n", optarg);
124  break;
125  case 'u' :
126  if (update_chi != NULL) {
127  *update_chi = 1;
128  printf("Enabling automatic chi updating based on truncation error.\n");
129  }
130  break;
131  case 2 :
132  if (emax != NULL) {
133  *emax = atof(optarg);
134  printf("Truncation error tolerance per sweep for updating chi set to <<<%7.2e>>>\n", *emax);
135  }
136  break;
137  case 3 :
138  if (delta_chi != NULL) {
139  *delta_chi = atoi(optarg);
140  printf("Amount by which to change chi when updating chi set to <<<%d>>>\n", *delta_chi);
141  }
142  break;
143  case 4 :
144  if (max_chi != NULL) {
145  *max_chi = atoi(optarg);
146  printf("Maximum value of chi to use when updating chi set to <<<%d>>>\n", *max_chi);
147  }
148  break;
149  case 5 :
150  printf("Options for TNT TEBD routines\n");
151  printf("\n");
152  printf("%s <parameters>:\n", argv[0]);
153  printf("Input parameters:\n");
154  printf("\n");
155  printf("--help-tebd | : Print this help.\n");
156  printf("--timesteps <n> | -t <n> : Number of timesteps to be performed. If not given default 50 is set. \n");
157  printf("--tbigstep <n> | -b <n> : Number of timesteps to be performed between each calculation of expectation values. \n");
158  printf(" If not given, then expectation values are calculated at the end of the routine only.\n");
159  printf("--trunc-type <str> | : How to calculate the truncation error using the discarded singular values.\n");
160  printf("--update-chi | -u : Update chi with increasing truncation error. Turned off by default. \n");
161  printf("--error-tol <val> | : Truncation error tolerance for increasing chi. Default is 1e-6. \n");
162  printf("--delta-chi | : Amount by which to increase chi. Default value 50. \n");
163  printf("--max-chi | : Value at which to stop increasing chi, even if the error exceeds error-tol. Default value 500. \n");
164  exit(0);
165  default: /* '?' Do nothing since argument may be treated by another function */
166  break;
167  }
168 
169  }
170 
171 
172  /* Set unset variables to default values */
173  if ((numsteps != NULL ) && (0 == *numsteps)) {
174  *numsteps = 50;
175  printf("numsteps: Using default <<<%d>>>\n", *numsteps);
176  }
177  if ((tbigstep != NULL ) && (0 == *tbigstep)) {
178  *tbigstep = *numsteps;
179  printf("tbigstep: Using default <<<%d>>>\n", *tbigstep);
180  }
181  if ((update_chi != NULL) && (1 == *update_chi)) {
182  if ((emax != NULL) && (*emax < -99.0)) {
183  *emax = 1e-6;
184  printf("Truncation error tolerance per sweep for updating chi set to default <<<%7.2e>>>\n", *emax);
185  }
186  if ((delta_chi != NULL ) && (0 == *delta_chi)) {
187  *delta_chi = 50;
188  printf("Amount by which to change chi when updating chi set to default <<<%d>>>\n", *delta_chi);
189  }
190  if ((max_chi != NULL ) && (0 == *max_chi)) {
191  *max_chi = 500;
192  printf("Maximum value of chi to use when updating chi set to default <<<%d>>>\n", *max_chi);
193  }
194  }
195 
196  free(argv_cp);
197 
198  return;
199 }
void tntTruncType(const char *funcname)
Definition: tntMisc_funcs.c:497