tntMpsProcessInitOptions.c

/*
 * Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch (Oxford) and Chris Goodyer (NAG)
 * $LastChangedDate$
 * (c) University of Oxford 2014
 */

/* Include the internal header for the TNT library */
#include "tntMpsInternal.h"

#include <getopt.h>
#include <ctype.h>

#define TNT_SH_NUM_TERMS 12
#define TNT_SH_NUM_NN_TERMS 9

#define TNT_BS_NUM_TERMS 7
#define TNT_BS_NUM_NN_TERMS 3

void tntMpsProcessSysOptions(int argc, 
                             char **argv, 
                             tntNetwork *mps, 
                             tntNetwork *mpo, 
                             tntNetwork *prop, 
                             tntMpsExOp *ExOp) 
{
    
    static struct option all_long_options[] = { /* Long options for all system */
    {"dt", required_argument, 0, 1},
    {"idt", required_argument, 0, 2},
    {"rand-state", required_argument, 0, 3},
    {"qnum-rand-state", required_argument, 0, 4},
    {"config-state", required_argument, 0, 5},
    {"qnum-config-state", required_argument, 0, 6},
    {"length", required_argument, 0, 7},
    {"system", required_argument, 0, 8},
    {"qnum-eye-state", no_argument, 0, 9},
    {"eye-state", no_argument, 0, 10},
    {0, 0, 0, 0}
    };
    static struct option sh_long_options[] = { /* Long options for a spin half system */
    {"Jxx", required_argument, 0, 1},
    {"Jyy", required_argument, 0, 2},
    {"Jzz", required_argument, 0, 3},
    {"Jxy", required_argument, 0, 4},
    {"Jyx", required_argument, 0, 5},
    {"Jyz", required_argument, 0, 6},
    {"Jzy", required_argument, 0, 7},
    {"Jzx", required_argument, 0, 8},
    {"Jxz", required_argument, 0, 9},
    {"Bx", required_argument, 0, 10},
    {"By", required_argument, 0, 11},
    {"Bz", required_argument, 0, 12},
    {"2S", required_argument, 0, 13},
    {0, 0, 0, 0}
    };
    static struct option bs_long_options[] = { /* Long options for a boson system */
    {"Jb", required_argument, 0, 1},
    {"Vb", required_argument, 0, 2},
    {"param-driv", required_argument, 0, 3},
    {"Ub", required_argument, 0, 4},
    {"mub", required_argument, 0, 5},
    {"E-harm", required_argument, 0, 6},
    {"coher-driv", required_argument, 0, 7},
    {"jc-harm", required_argument, 0, 8},
    {"n-max", required_argument, 0, 9},
    {0, 0, 0, 0}
    };
    int option_index;
    int c; /* Variable that contains the returned option code from getopt_long */
    char **argv_cp; /* Used to hold a copy of the arguments - this will be used by the function instead of the original arguments */
    unsigned numnn = 0; /* Number of nearest neighbour terms in the Hamiltonian */
    unsigned numos = 0; /* Number of onsite terms in the Hamiltonian */
    unsigned no_operators = 0; /* Marks whether either H or U should be calculated */
    double dt = 0.0, idt = 0.0; /* initial time step sizes */
    tntComplex dtc; /* The complex time step size */
    unsigned MPStype = 0; /* Use 1 for random, use 2 for config state, 0 if no MPS type has been assigned yet. */
    char config[TNT_STRLEN]; /* Code for the configuration */
    unsigned L=0; /* Length of system - set to zero initially to see if it has been set */
    unsigned D; /* Starting internal dimension of the random MPS */
    int Q; /* Quantum number for symmetric MPS */
    tntNodeArray nnL, nnR, os; /* Arrays for nearest neighbour and onsite operators */
    tntComplexArray nnparam, osparam; /* Arrays for the parameters for the nearest neighbour and onsite operators */
    unsigned loop, j;
    
    /* Set optind to 1 in case option getting routine has already been called previously */
    optind = 1;
    
    /* Set opterr to zero to prevent error message for unknown erguments, since they may be being treated by another process options function */
    opterr = 0;
    
    /* Make a copy of argv - getopt_long will permute the arguments, and they may be required by another processing function */
    argv_cp = (char **) malloc(argc*sizeof(char *));
    for (c = 0; c<argc; c++) {
        argv_cp[c] = argv[c];
    }
    
    /* 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 */
    if ((NULL == prop) && (NULL == mpo)) {
        no_operators = 1;
    }
    
    tntPrintf("----------------------------------------------------------\n");
    tntPrintf("Creating MPS system with the following parameters: \n");
    tntPrintf("----------------------------------------------------------\n");
    
    /* Loop over all input arguments */
    while ((c = getopt_long(argc, argv_cp, "", all_long_options, &option_index)) != -1) {
        switch (c) {
            case 1: /* This option relates to the real time step size */
                if (NULL == prop) {
                    tntWarningPrint("The input parameter dt is being ignored as U is NULL|");
                } else {
                    dt = atof(optarg);
                    tntPrintf("dt: <<<%g>>>\n",  dt);
                }
                break;
            case 2:  /* This option relates to the imaginary time step size */
                if (NULL == prop) {
                    tntWarningPrint("The input parameter idt is being ignored as U is NULL|");
                } else {
                    idt = atof(optarg);
                    tntPrintf("idt: <<<%g>>>\n",  idt);
                }
                break;
            case 3: /* This option relates to specifiying a random starting MPS */
                if (NULL == mps) {
                    tntWarningPrint("The input parameter rand-state is being ignored as argument mps is NULL|");
                } else if (MPStype) {
                    tntWarningPrint("The input parameter rand-state is being ignored as the MPS type has already been set|");
                } else {
                    MPStype = 1;
                    D = abs(atoi(optarg));
                    tntPrintf("Random MPS being creating with starting internal dimension: <<<%d>>>\n",D);
                }
                break;
            case 4: /* This option relates to specifying a given quantum number for the starting random MPS, and conserving this number */
                if (NULL == mps) {
                    tntWarningPrint("The input parameter qnum-eye-state is being ignored as argument mps is NULL|");
                } else if (MPStype) {
                    tntWarningPrint("The input parameter qnum-eye-state is being ignored as the MPS type has already been set|");
                } else {
                    MPStype = 3;
                    /* Extract total quantum number for the system */
                    Q = atoi(optarg);
                    /* Set the symmetry type for the system */
                    tntSymmTypeSet("U(1)", 1);
                    tntPrintf("Creating MPS with initial quantum number: <<<%d>>>\n",Q);
                }
                break;
            case 5: /* This option relates to specifying a configuration for the starting MPS */
                if (NULL == mps) {
                    tntWarningPrint("The input parameter config-state is being ignored as argument mps is NULL|");
                } else if (MPStype) {
                    tntWarningPrint("The input parameter config-state is being ignored as the MPS type has already been set|");
                } else {
                    MPStype = 2;
                    strcpy(config, optarg);
                    tntPrintf("Creating MPS with configuration: <<<%s>>>\n",config);
                }
                break;
            case 6: /* This option relates to specifying a configuration for the starting MPS, and conserving its quantum number */
                if (NULL == mps) {
                    tntWarningPrint("The input parameter qnum-state is being ignored as argument mps is NULL|");
                } else if (MPStype) {
                    tntWarningPrint("The input parameter qnum-state is being ignored as the MPS type has already been set|");
                } else {
                    MPStype = 2;
                    strcpy(config, optarg);
                    tntPrintf("Creating number conserving MPS with configuration: <<<%s>>>\n",config);
                    /* Set the symmetry type for the system */
                    tntSymmTypeSet("U(1)", 1);
                }
                break;
            case 7:  /* This option (required) relates to specifying the length of the starting MPS */
                L = abs(atoi(optarg));
                tntPrintf("Length: <<<%d>>>\n",L);
                break;
            case 8:  /* This option (required) relates to the type of system that will be created */
                if (tntSysTypeGet()) {
                    tntWarningPrint("The input parameter system is being ignored as the system type has already been set|");
                } else if (!strcmp(optarg,"spin")) {
                    tntPrintf("system: <<<spin>>>\n");
                    tntSysTypeSet(TNT_MPS_SPIN);
                } else if (!strcmp(optarg,"boson")) {
                    tntPrintf("system: <<<boson>>>\n");
                    tntSysTypeSet(TNT_MPS_BOSON);
                } else {
                    tntErrorPrint("System type %s is not known",optarg);
                }
                break;
            case 9:  /* This option (required) relates to the type of system that will be created */
                if (NULL == mps) {
                    tntWarningPrint("The input parameter qnum-eye-state is being ignored as argument mps is NULL|");
                } else if (MPStype) {
                    tntWarningPrint("The input parameter qnum-eye-state is being ignored as the start-state type has already been set|");
                } else {
                    MPStype = 4;
                    tntPrintf("Creating identity density matrix as start state, with quantum number conservation turned on\n");
                    tntSymmTypeSet("U(1)", 1);
                }
            case 10:  /* This option (required) relates to the type of system that will be created */
                if (NULL == mps) {
                    tntWarningPrint("The input parameter qnum-eye-state is being ignored as argument mps is NULL|");
                } else if (MPStype) {
                    tntWarningPrint("The input parameter qnum-eye-state is being ignored as the start-state type has already been set|");
                } else {
                    MPStype = 4;
                    tntPrintf("Creating identity density matrix as start state\n");
                }
            default:
                break;
        }
    }
    
    if (0 == L) {
        tntErrorPrint("The system cannot be initialised as the length has not been set"); /* NO_COVERAGE */
    } /* NO_COVERAGE */
    
    if ((0 == MPStype)&&(mps != NULL)) {
        tntErrorPrint("The wave function cannot be created as the starting state has not been set"); /* NO_COVERAGE */
    } /* NO_COVERAGE */
    
    if (TNT_MPS_SPIN == tntSysTypeGet()) { /* Do this bit if spin system chosen and there are operators */
        
        unsigned TwoS = 1; /* Spin for each particle in the system */
        tntComplex sh_params[TNT_SH_NUM_TERMS]; /* Multipliers for each of the terms in the Hamiltonian */
        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 */
        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 */
        
        /* Make a copy of argv - getopt_long will permute the arguments, and they may be required by another processing function */
        for (c = 0; c<argc; c++) {
            argv_cp[c] = argv[c];
        }
        
        /* Set optind to 1 in case option getting routine has already been called previously */
        optind = 1;
        
        tntPrintf("----------------------------------------------------------\n");
        tntPrintf("Creating spin system with the following parameters: \n");
        tntPrintf("----------------------------------------------------------\n");
        
        /* Loop over all input arguments */
        while ((c = getopt_long(argc, argv_cp, "", sh_long_options, &option_index)) != -1) {
            if (c <= TNT_SH_NUM_TERMS) { /* These options relate to terms in the Hamiltonian */
                /* Populate the array of parameters */
                sh_params[c-1] = tntAtoC(optarg);
                sh_param_present[c-1] = 1;
                if (fabs(sh_params[c-1].re) > DBL_EPSILON && fabs(sh_params[c-1].im) > DBL_EPSILON) {
                    tntPrintf("%s : <<<%g%+gi>>>\n", sh_param_strs[c-1], sh_params[c-1].re, sh_params[c-1].im);
                } else if (fabs(sh_params[c-1].re) > DBL_EPSILON) {
                    tntPrintf("%s : <<<%g>>>\n", sh_param_strs[c-1], sh_params[c-1].re);
                } else if (fabs(sh_params[c-1].im) > DBL_EPSILON) {
                    tntPrintf("%s : <<<%gi>>>\n", sh_param_strs[c-1], sh_params[c-1].im);
                } else {
                    tntPrintf("%s : <<<0>>>\n", sh_param_strs[c-1]);
                }
                if (c <= TNT_SH_NUM_NN_TERMS) {
                    numnn++;
                } else {
                    numos++;
                }
            } else if (13 == c) { /* Term for the spin for each particle */
                TwoS = abs(atoi(optarg));
                tntPrintf("2S : <<<%d>>>\n", TwoS);
            }
        }
        
        tntNode Sx, Sy, Sz, Sp, Sm, eye; /* Spin operators */
        
        /* Creating the Pauli operators */
        tntMpsCreateSpinOp(TwoS, &Sx, &Sy, &Sz, &Sp, &Sm, &eye);
        
        /* Set the operator for the basis */
        tntSysBasisOpSet(Sz);
        
        if (!no_operators) {
            
            /* First check whether the operators actually exist */
            if (0 == numnn && 0 == numos) {
                tntErrorPrint("Cannot create operators as all the parameters are zero"); /* NO_COVERAGE */
            }
            
            
            /* First step: if quantum number conservation is turned on, check for non-invariant operators */
            /* Only non-invariant operator combinations are: Jxx = Jyy, Jxy = -Jyx, Jzz anything, Bz anything */
            /* 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) */
            /* if non-allowed terms are present delete them and present warning to user (here) */
            if (TNT_SYMM_U1 == tntSymmTypeGet()) {
                /* Look for non-invariant parts, which cannot be re-expressed as an invariant + non-invariant part */
                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]) {
                    char full_msg[TNT_STRLEN];
                    strcpy(full_msg,"Your Hamiltonian is not invariant under the symmetry type chosen|");
                    for (loop = 5; loop < 11; loop++) {
                        if (sh_param_present[loop]) {
                            strcat(full_msg, sh_param_strs[loop]);
                            strcat(full_msg, " ");
                            sh_param_present[loop] = 0;
                            if (loop < TNT_SH_NUM_NN_TERMS) {
                                numnn--;
                            } else {
                                numos--;
                            }
                        }
                    }
                    strcat(full_msg, "being set to zero|");
                    tntWarningPrint(full_msg);
                }
            }
            
            /* Create arrays to hold the terms and parameters that are present */
            if (numnn) {
                nnL = tntNodeArrayAlloc(numnn);
                nnR = tntNodeArrayAlloc(numnn);
                nnparam = tntComplexArrayAlloc(numnn);
            } else {
                nnL.sz = 0;
                nnR.sz = 0;
                nnparam.numrows = nnparam.numcols = 0;
            }
            if (numos) {
                os = tntNodeArrayAlloc(numos);
                osparam = tntComplexArrayAlloc(numos);
            } else {
                os.sz = 0;
                osparam.numrows = osparam.numcols = 0;
            }
            
            /* reset the counters */
            numnn = numos = 0;
            
            /* Second step: if quantum number conservation is turned on, check for non-invariant operators Jxx, Jyy, Jxy, Jyx */
            /* These can be expressed as combinations of Sp and Sm (invariant) rather than Sx and Sy */
            /* Only non-invariant operator combinations are: Jxx = Jyy, Jxy = -Jyx. If different values are given present warning to user */
            if (TNT_SYMM_U1 == tntSymmTypeGet()) {
                /* Look for Jxx or Jyy term */
                if (sh_param_present[0] || sh_param_present[1]) {
                    if (fabs(sh_params[0].re-sh_params[1].re) + fabs(sh_params[0].im-sh_params[1].im) > DBL_EPSILON) {
                        tntWarningPrint("Your Hamiltonian is not invariant under the symmetry type chosen| Only Jxx=Jyy allowed - values being changed to average of Jxx and Jyy");
                    }
                    /* create the term 0.25*(Jxx+Jyy)(SpSm + SmSp) = JxxSxSx + JyySySy (+ any non-invariant parts being ignored) */
                    nnL.vals[numnn] = tntNodeCopy(Sp); /* Operator acting on left */
                    nnR.vals[numnn] = tntNodeCopy(Sm); /* Operator acting on right */
                    nnparam.vals[numnn].re = 0.25*(sh_params[0].re + sh_params[1].re); /* Copy over parameter */
                    nnparam.vals[numnn].im = 0.25*(sh_params[0].im + sh_params[1].im); /* Copy over parameter */
                    numnn++;
                    nnL.vals[numnn] = tntNodeCopy(Sm); /* Operator acting on left */
                    nnR.vals[numnn] = tntNodeCopy(Sp); /* Operator acting on right */
                    nnparam.vals[numnn].re = 0.25*(sh_params[0].re + sh_params[1].re); /* Copy over parameter */
                    nnparam.vals[numnn].im = 0.25*(sh_params[0].im + sh_params[1].im); /* Copy over parameter */
                    numnn++;
                    /* Since values have already been accounted for, set flag to zero so they are not included again */
                    sh_param_present[0] = sh_param_present[1] = 0;
                }
                /* Look for Jxy or Jyx term */
                if (sh_param_present[3] || sh_param_present[4]) {
                    if (fabs(sh_params[3].re + sh_params[4].re) + fabs(sh_params[3].im + sh_params[4].im) > DBL_EPSILON) {
                        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");
                    }
                    /* create the term 0.25i*(Jxy-Jyx)(SpSm - SmSp) = JxySxSy + JyxSySx (+ any non-invariant parts being ignored) */
                    nnL.vals[numnn] = tntNodeCopy(Sp); /* Operator acting on left */
                    nnR.vals[numnn] = tntNodeCopy(Sm); /* Operator acting on right */
                    nnparam.vals[numnn].im =  0.25*(sh_params[3].re - sh_params[4].re); /* Copy over parameter */
                    nnparam.vals[numnn].re = -0.25*(sh_params[3].im - sh_params[4].im); /* Copy over parameter */
                    numnn++;
                    nnL.vals[numnn] = tntNodeCopy(Sm); /* Operator acting on left */
                    nnR.vals[numnn] = tntNodeCopy(Sp); /* Operator acting on right */
                    nnparam.vals[numnn].im =  0.25*(sh_params[4].re - sh_params[3].re); /* Copy over parameter */
                    nnparam.vals[numnn].re = -0.25*(sh_params[4].im - sh_params[3].im); /* Copy over parameter */
                    numnn++;
                    /* Since values have already been accounted for, set flag to zero so they are not included again */
                    sh_param_present[3] = sh_param_present[4] = 0;
                }
            }
            
            
            /* Loop through parameters and assign them to the correct array */
            for (loop = 0; loop < TNT_SH_NUM_TERMS; loop++) {
                if (sh_param_present[loop]) {
                    switch (loop+1) {
                        case 1: /* Term Jxx */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 2: /* Term Jyy */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 3: /* Term Jzz */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 4: /* Term Jxy */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 5: /* Term Jyx */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 6: /* Term Jyz */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 7: /* Term Jzy */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sy); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 8: /* Term Jzx */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 9: /* Term Jxz */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(Sx); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(Sz); /* Operator acting on right */
                            nnparam.vals[numnn] = sh_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 10: /* Term Bx */
                            /* ---- Set the nodes to make the term ---- */
                            os.vals[numos] = tntNodeCopy(Sx); /* Operator acting on site */
                            osparam.vals[numos] = sh_params[loop]; /* Copy over parameter */
                            numos++;
                            break;
                        case 11: /* Term By */
                            /* ---- Set the nodes to make the term ---- */
                            os.vals[numos] = tntNodeCopy(Sy); /* Operator acting on site */
                            osparam.vals[numos] = sh_params[loop]; /* Copy over parameter */
                            numos++;
                            break;
                        case 12: /* Term Bz */
                            /* ---- Set the nodes to make the term ---- */
                            os.vals[numos] = tntNodeCopy(Sz); /* Operator acting on site */
                            osparam.vals[numos] = sh_params[loop]; /* Copy over parameter */
                            numos++;
                            break;
                        default:
                            break;
                    }
                }
            }
        }
        /* Free the spin operators */
        tntNodeFree(&Sx);
        tntNodeFree(&Sy);
        tntNodeFree(&Sz);
        tntNodeFree(&Sp);
        tntNodeFree(&Sm);
        tntNodeFree(&eye);
    }
    if (TNT_MPS_BOSON == tntSysTypeGet()) { /* Do this bit if boson system chosen and there are operators */
        
        double jc = (L-1.0)/2.0; /* Centre for the harmonic trap */
        unsigned n_max = 1; /* Maximum number of bosons on each site */
        tntComplex bs_params[TNT_BS_NUM_TERMS]; /* Multipliers for each of the terms in the Hamiltonian */
        unsigned bs_param_present[TNT_BS_NUM_TERMS] = {0,0,0,0,0,0,0}; /* Specifies whether each term is present */
        char bs_param_strs[TNT_BS_NUM_TERMS][12] = {"Jb","Vb","param-driv","Ub","mub","E-harm","coher-driv"}; /* String for each term */
        
        tntPrintf("----------------------------------------------------------\n");
        tntPrintf("Creating boson system with the following parameters: \n");
        tntPrintf("----------------------------------------------------------\n");
        
        
        /* Make a copy of argv - getopt_long will permute the arguments, and they may be required by another processing function */
        for (c = 0; c<argc; c++) {
            argv_cp[c] = argv[c];
        }
        
        /* Set optind to 1 in case option getting routine has already been called previously */
        optind = 1;
        
        /* Reset counters for number of operators */
        numnn = numos = 0;
        
        /* Loop over all input arguments */
        while ((c = getopt_long(argc, argv_cp, "", bs_long_options, &option_index)) != -1) {
            if (c <= TNT_BS_NUM_TERMS) { /* These options relate to terms in the Hamiltonian */
                /* Populate the array of parameters */
                bs_params[c-1] = tntAtoC(optarg);
                bs_param_present[c-1] = 1;
                if (fabs(bs_params[c-1].re) > DBL_EPSILON && fabs(bs_params[c-1].im) > DBL_EPSILON) {
                    tntPrintf("%s : <<<%g%+gi>>>\n", bs_param_strs[c-1], bs_params[c-1].re, bs_params[c-1].im);
                } else if (fabs(bs_params[c-1].re) > DBL_EPSILON) {
                    tntPrintf("%s : <<<%g>>>\n", bs_param_strs[c-1], bs_params[c-1].re);
                } else if (fabs(bs_params[c-1].im) > DBL_EPSILON) {
                    tntPrintf("%s : <<<%gi>>>\n", bs_param_strs[c-1], bs_params[c-1].im);
                } else {
                    tntPrintf("%s : <<<0>>>\n", bs_param_strs[c-1]);
                }
                if ((1 == c)||(3 == c)) { /* Two nearest neighbour terms for hopping or parametric driving */
                    numnn += 2;
                } else if (c <= TNT_BS_NUM_NN_TERMS) {
                    numnn++;
                } else if (7 == c) { /* Two onsite terms for coherent driving */
                    numos += 2;
                } else {
                    numos++;
                }
            } else if (8 == c) { /* Term for centre of harmonic potential */
                jc = atof(optarg);
                tntPrintf("jc-harm : <<<%g>>>\n", jc);
            } else if (9 == c) { /* Term for maximum number of atoms per site */
                n_max = abs(atoi(optarg));
                tntPrintf("n-max : <<<%d>>>\n", n_max);
            }
        }
        
        tntNode b_op, bdag_op, n_op, os_int, eye; /* Boson operators */
        
        /* Creating the boson operators */
        tntMpsCreateBosonOp(n_max, &b_op, &bdag_op, &n_op, &os_int, &eye);
        
        /* Set the operator for the basis */
        tntSysBasisOpSet(n_op);
        
        if (!no_operators) {
            
            /* First check whether the operators actually exist */
            if (0 == numnn && 0 == numos) {
                tntErrorPrint("Cannot create operators as all the parameters are zero"); /* NO_COVERAGE */
                exit(1); /* NO_COVERAGE */
            }
            
            /* First step: if quantum number conservation is turned on, check for non-invariant operators i.e. param-driv or coher-driv */
            if (TNT_SYMM_U1 == tntSymmTypeGet()) {
                /* Look for terms */
                if (bs_param_present[2] && bs_param_present[6]) {
                    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]);
                    bs_param_present[2] = bs_param_present[6] = 0;
                    numnn -= 2;
                    numos -= 2;
                } else if (bs_param_present[2]) {
                    tntWarningPrint("Your Hamiltonian is not invariant under the symmetry type chosen|%s being set to zero|", bs_param_strs[2]);
                    bs_param_present[2] = 0;
                    numnn -= 2;
                } else if (bs_param_present[6]) {
                    tntWarningPrint("Your Hamiltonian is not invariant under the symmetry type chosen|%s being set to zero|", bs_param_strs[6]);
                    bs_param_present[6] = 0;
                    numos -= 2;
                }
                
            }
            
            /* Create arrays to hold the terms and parameters that are present */
            if (numnn) {
                nnL = tntNodeArrayAlloc(numnn);
                nnR = tntNodeArrayAlloc(numnn);
                nnparam = tntComplexArrayAlloc(numnn);
            } else {
                nnL.sz = 0;
                nnR.sz = 0;
                nnparam.numrows = nnparam.numcols = 0;
            }
            if (numos) {
                os = tntNodeArrayAlloc(numos);
                if (bs_param_present[5]) /* If there is a harmonic trapping term, then there will be different onsite term for each site */
                    osparam = tntComplexArrayAlloc(numos,L);
                else
                    osparam = tntComplexArrayAlloc(numos);
            } else {
                os.sz = 0;
                osparam.numrows = osparam.numcols = 0;
            }
            
            /* reset the counters */
            numnn = numos = 0;
            
            /* Loop through parameters and assign them to the correct array */
            for (loop = 0; loop < TNT_BS_NUM_TERMS; loop++) {
                if (bs_param_present[loop]) {
                    switch (loop+1) {
                        case 1: /* Term Jb - two hopping terms required */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(b_op); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(bdag_op); /* Operator acting on right */
                            nnparam.vals[numnn].re = -bs_params[loop].re; /* Copy over parameter */
                            nnparam.vals[numnn].im = -bs_params[loop].im; /* Copy over parameter */
                            nnL.vals[numnn+1] = tntNodeCopy(bdag_op); /* Operator acting on left */
                            nnR.vals[numnn+1] = tntNodeCopy(b_op); /* Operator acting on right */
                            nnparam.vals[numnn+1].re = -bs_params[loop].re; /* Copy over cc parameter */
                            nnparam.vals[numnn+1].im = bs_params[loop].im; /* Copy over cc parameter */
                            numnn += 2;
                            break;
                        case 2: /* Term Vb: nearest neighbour interaction */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(n_op); /* Operator acting on left */
                            nnR.vals[numnn] = tntNodeCopy(n_op); /* Operator acting on right */
                            nnparam.vals[numnn] = bs_params[loop]; /* Copy over parameter */
                            numnn++;
                            break;
                        case 3: /* Term param-driv: parametric driving */
                            /* ---- Set the nodes to make the term ---- */
                            nnL.vals[numnn] = tntNodeCopy(b_op); /* Operator acting on left - set to b_op */
                            nnR.vals[numnn] = tntNodeCopy(b_op); /* Operator acting on right - set to b_op */
                            nnparam.vals[numnn] = bs_params[loop]; /* Copy over parameter */
                            
                            nnL.vals[numnn+1] = tntNodeCopy(bdag_op); /* Operator acting on right - set to bdag_op */
                            nnR.vals[numnn+1] = tntNodeCopy(bdag_op); /* Operator acting on left - set to bdag_op */
                            nnparam.vals[numnn+1].re = bs_params[loop].re; /* Copy over cc parameter */
                            nnparam.vals[numnn+1].im = -bs_params[loop].im; /* Copy over cc parameter */
                            
                            numnn += 2;
                            break;
                        case 4: /* Term Ub - onsite interaction */
                            /* ---- Set the nodes to make the term ---- */
                            os.vals[numos] = tntNodeCopy(os_int); /* Operator acting on site */
                            if (bs_param_present[5]) { /* If the harmonic trapping term is present, a term will be required for every site */
                                for (j = 0; j < L; j++) {
                                    osparam.vals[numos + j*os.sz].re = bs_params[loop].re/2.0; /* Copy over parameter */
                                    osparam.vals[numos + j*os.sz].im = bs_params[loop].im/2.0; /* Copy over parameter */
                                }
                            } else {
                                osparam.vals[numos].re = bs_params[loop].re/2.0; /* Copy over parameter */
                                osparam.vals[numos].im = bs_params[loop].im/2.0; /* Copy over parameter */
                            }
                            numos++;
                            break;
                        case 5: /* Term mub - chemical potential */
                            /* ---- Set the nodes to make the term ---- */
                            os.vals[numos] = tntNodeCopy(n_op); /* Operator acting on site */
                            if (bs_param_present[5]) { /* If the harmonic trapping term is present, a term will be required for every site */
                                for (j = 0; j < L; j++) {
                                    osparam.vals[numos + j*os.sz].re = -bs_params[loop].re; /* Copy over parameter */
                                    osparam.vals[numos + j*os.sz].im = -bs_params[loop].im; /* Copy over parameter */
                                }
                            } else {
                                osparam.vals[numos].re = -bs_params[loop].re; /* Copy over parameter */
                                osparam.vals[numos].im = -bs_params[loop].im; /* Copy over parameter */
                            }
                            numos++;
                            break;
                        case 6: /* Term Eharm: harmonic trapping term */
                            /* ---- Set the nodes to make the term ---- */
                            os.vals[numos] = tntNodeCopy(n_op); /* Operator acting on site */
                            for (j = 0; j < L; j++) {
                                /* Copy over parameter and covert to harmonic trap magnitude on that site*/
                                osparam.vals[numos + j*os.sz].re = 0.5*bs_params[loop].re*(j - jc)*(j - jc);
                                osparam.vals[numos + j*os.sz].im = 0.5*bs_params[loop].im*(j - jc)*(j - jc);
                            }
                            numos++;
                            break;
                        case 7: /* Term coher-driv:coherent driving term */
                            /* ---- Set the nodes to make the term ---- */
                            os.vals[numos] = tntNodeCopy(b_op); /* Operator acting on site */
                            os.vals[numos+1] = tntNodeCopy(bdag_op); /* Operator acting on site */
                            
                            if (bs_param_present[5]) { /* If the harmonic trapping term is present, a term will be required for every site */
                                for (j = 0; j < L; j++) {
                                    osparam.vals[numos + j*os.sz] = bs_params[loop]; /* Copy over parameter */
                                    osparam.vals[numos + 1 + j*os.sz].re =  bs_params[loop].re; /* Copy over cc parameter */
                                    osparam.vals[numos + 1 + j*os.sz].im = -bs_params[loop].im; /* Copy over cc parameter */
                                }
                            } else {
                                osparam.vals[numos] = bs_params[loop]; /* Copy over parameter */
                                osparam.vals[numos+1].re =  bs_params[loop].re; /* Copy over cc parameter */
                                osparam.vals[numos+1].im = -bs_params[loop].im; /* Copy over cc parameter */
                            }
                            numos += 2;
                            break;
                        default:
                            break;
                    }
                }
            }
        }
        /* Free the boson operators */
        tntNodeFree(&b_op);
        tntNodeFree(&bdag_op);
        tntNodeFree(&n_op);
        tntNodeFree(&os_int);
        tntNodeFree(&eye);
        
    } else if ((!no_operators) && (0 == tntSysTypeGet())) {
        tntErrorPrint("No operators can be created as the system type has not been set"); /* NO_COVERAGE */
    }
    
    
    if (NULL != mpo) {
        /* Call the function that makes MPOs */
        *mpo = tntMpsCreateMpo(L, &nnL, &nnR, &nnparam, &os, &osparam);
    }
    
    if (NULL != prop) {
        /* Contruct the time step */
        dtc.re = dt;
        dtc.im = idt;
        
        /* Call the function that makes the propagators */
        *prop = tntMpsCreatePropST2sc(L, dtc, &nnL, &nnR, &nnparam, &os, &osparam);
    }
    
    /* Create the starting wave function */
    if (MPStype) {
        if (1 == MPStype) {
            *mps = tntMpsCreateRandom(L,D);
        } else if (2 == MPStype) {
            *mps = tntMpsCreateConfig(L,config);
        } else if (3 == MPStype) {
            *mps = tntMpsCreateSymmRandom(L,&Q);
        } else if (4 == MPStype) {
            *mps = tntMpsCreateEyeMpo(L,NULL);
        }
    }
    
    if (!no_operators) {
        /* Free the operator arrays if they have been created */
        if (nnL.sz) {
            tntNodeArrayFree(&nnL);
            tntNodeArrayFree(&nnR);
            tntComplexArrayFree(&nnparam);
        }
        
        if (os.sz) {
            tntNodeArrayFree(&os);
            tntComplexArrayFree(&osparam);
        }
        
    }
    
    free(argv_cp);
    
    /* Now deal with expectation values - separate function used for this */
    tntMpsProcessExpecOptions(argc, argv, ExOp);
    
}