tntEvolve_if.c

Example of using the MPS library evolving a one-dimensional system, where the start state can either be a ground state or a product MPS. The product state, Hamiltonian that defines the ground state, and the Hamiltonian for time evolution are given in an initialisation file.

The ground state is found using the DMRG algorithm - see tntMpsVarMinMpo2sStep() for more details.

The time evolution is carried out using the TEBD algorithm - see tntMpsPropST2scProduct() for more details.

Running the program

To compile this main file change to the source directory then type make tntEvolve_if . You may first have to amend makefile.inc to point to your copies of the required external linear algebra and MATLAB libraries.

To run the file you may also have to set up certain environment variables. To do this edit then run ./set_tnt_vars.sh in the scripts directory.

You then need to create an initialisation file. See the examples ending in .m given in the 'inititialisation' directory and listed here, and the page initfiles. You would then call the function with arguments for the initialisation file and where to save the data e.g.

./bin/tntEvolve_if -d output/savename.mat -i initialisation/inname.mat

Program input

There are some pre-made initialisation files you can load in directly, all found in the initialisation directory and listed in the examples in this documentation (files ending in .m).

Program output

The program will output the expectation values specified in the initialisation file.

In addition the program will always output:

timestep

An array of timesteps for which expectation values were taken. A timestep of -1 corresponds to the ground state.

t

An array of times at which expectation values were taken (i.e. timestep \(\times\) dt). A time of -1 corresponds to the ground state.

normsq

The squared norm of the state for each time-step.

err

The truncation error for each time-step.

tntSystem

Global system parameters used for the calculation.

tntLibVersion

Version of the library used.

tntLibType

An integer code which describes which external libraries (e.g. ARPACK, NetCDF) where compiled in with this version of the library.

/*
Authors: Sarah Al-Assam, Stephen Clark and Dieter Jaksch
Date:    $LastChangedDate$
(c) University of Oxford 2016
*/

#include "tntMps.h"

unsigned checkParams(tntComplexArray *nntsc, tntComplexArray *ostsc, int tstep);
void setParams(tntComplexArray *nnparam, tntComplexArray *nntsc, tntComplexArray *nnpsc, 
                  tntComplexArray *osparam, tntComplexArray *ostsc, tntComplexArray *ospsc, int tstep);

int main(int argc, char **argv)
{

    char initfile[TNT_STRLEN];          /* Path to initialisation data */
    char savename[TNT_STRLEN];          /* Path to output file */
    
    tntExOp ExOp;                       /* Defines all the operators for calculating expectation values */
    tntNetwork wf_gs, wf_start, wf_ev;  /* The MPS wave function for the ground state, start state and evolved state */
    tntNetwork H, mod, prop;            /* The MPO Hamiltonian, MPO modifier and the staircase propagator */
    int chi, L;                         /* Maximum internal dimension for the MPS, length of the system */
    int mod_renorm;                     /* Flag to renormalise the state when modifier applied */
    int calc_ol_gs, calc_ol_is;         /* Flags for whether or not to calculate overlaps with ground state or initial state */
    int calc_energy;                    /* Flag to specify whether energy of evolved state should be calculated */
    int U1symm_dyn;                     /* Flag to specify whether symmetry information is turned on for dynamic calculation */
    
    tntNodeArray nnl, nnr, os;          /* operators for building the hamiltonian */
    tntComplexArray nnparam, osparam;   /* parameters for building the hamiltonian */
    tntComplexArray nntsc, ostsc;       /* parameters for scaling the parameters of the time-varying hamiltonian */
    tntComplexArray nnpsc, ospsc;       /* scaled parameters of the time-varying hamiltonian */
    tntIntArray sitenum;                /* sitenumbers for the product MPO */
    double err, precision, E, Eprev;    /* The truncation error, precision for the ground state, energy  */
    tntIntArray qn_tot;                 /* Quantum numbers for the ground state */
    int i, i_max;                       /* Iteration counter, Maximum number of iterations */
    tntNodeArray HeffL, HeffR;          /* Precontracted nodes for the right and left sides of the network */
    
    tntIntArray modifypmpo;             /* Flag for type of each modifier to load */
    char osnm[15], ospnm[15], snnm[15]; /* Variable names for the modifiers to load */
    
    
    int numsteps, tbigstep;             /* Number of steps for run and between expectation values */
    int tstep = -1, bigstep = 0;        /* Counters for the above */ 
    tntComplex dt, t={-1,0};            /* Time step size, the current time */
    int changed_params;                 /* Flag for whether parameters in Hamiltonian need to be updated */
    double ol;                          /* output from the overlap calculations */
    
    /* Initialize the TNT library */
    tntInitialize();
    
    /* Process the command line options */
    tntProcessCLOptions(argc, argv, initfile, NULL, savename, NULL, NULL);
        
    /* Load general simulation parameters  */
    tntIntParamsLoad(initfile, chi, L);
    
    /* Load the operators required for expectation values */
    tntExOpsLoad(initfile, ExOp);

    /*************************************************************************************************/ 
    /*                                     Calculating ground state                                  */
    /*************************************************************************************************/ 
    
    /* Load precision with which to find the ground state */
    tntDoubleParamsLoad(initfile, precision);
    
    /* only perform calculation if the precision is > 0 */
    if (precision > 0.0) {
        
        /* Create the MPO Hamiltonian */
        tntNodeArraysNamedLoad(initfile,nnl,"nnlg",nnr,"nnrg",os,"osg");
        tntComplexArraysNamedLoad(initfile,nnparam,"nnparamg",osparam,"osparamg");
        H = tntMpsCreateMpo(L,&nnl,&nnr,&nnparam,&os,&osparam);
        tntNodeArrayFree(&nnl);
        tntNodeArrayFree(&nnr);
        tntNodeArrayFree(&os);
        tntComplexArrayFree(&nnparam);
        tntComplexArrayFree(&osparam);
        
        /* Load the starting wave function - if it's length is zero make a random wave function */
        tntNetworksLoad(initfile, wf_gs);
        if (0 == tntMpsLength(wf_gs)) {
            tntNetworkFree(&wf_gs);
            tntIntArraysLoad(initfile, qn_tot);
            /* Set a random number seed based on the current time */
            tntSysRandSeedSet(time(NULL));
            wf_gs = tntMpsCreateSymmRandom(L, qn_tot.vals);
        } else {
            tntMpsOrthNorm(wf_gs,0);
        }
        
        /* Initialise the nodes that will be needed for the DMRG sweep */
        tntMpsVarMinMpoInit(wf_gs, H, &HeffL, &HeffR);
        
        /* Determine the energy of the start state (the norm of the state is assumed to be 1) */
        E = tntComplexToReal(tntMpsMpoMpsProduct(wf_gs, H));
        tntDoubleParamsNamedSave(savename,E,"E_it");

        /* Perform minimization sweeps (first load maximum number of iterations) */
        tntIntParamsLoad(initfile, i_max);
        for (i = 0; i < i_max; i++) {
            Eprev = E;
            E = tntMpsVarMinMpo2sStep(wf_gs, chi, H, &HeffL, &HeffR, &err);

            tntDoubleParamsNamedUpdate(savename,i+1,E,"E_it");
            tntPrintf("Calculating ground state, Delta E_%d is %3.3e\n", i, Eprev-E);
            if ((Eprev - E < precision) && (Eprev - E > 0.0)) break;
        }
        
        tntNodeArrayFree(&HeffL);
        tntNodeArrayFree(&HeffR);
        tntNetworkFree(&H);
        
        tntMpsExpecOutput(wf_gs, &ExOp, 0, 1, 1, savename, bigstep);
        tntDoubleParamsUpdate(savename, bigstep, err, E); 
    } 
    tntIntParamsUpdate(savename, bigstep, tstep);
    tntComplexParamsUpdate(savename, bigstep, t);
    
    /*************************************************************************************************/ 
    /*                            Create start state for time evolution                              */
    /*************************************************************************************************/  
    
    /* determine whether to continue conserving quantum number */
    tntIntParamsLoad(initfile, U1symm_dyn);
    if (!U1symm_dyn) tntSymmTypeUnset();
    
    /* Load flags for modifiers to apply */
    tntIntArraysLoad(initfile, modifypmpo); 
    
    /* load the starting wave function */
    tntNetworksLoad(initfile, wf_start);
    
    /* if the length is zero, then start from the ground state instead (if it exists) */
    /* if no state exists, start from a vacuum */
    if (0 == tntMpsLength(wf_start)) {
        tntNetworkFree(&wf_start);
        if (precision > 0.0) wf_start = tntNetworkCopy(wf_gs);
        else wf_start = tntMpsCreateConfig(L,"");
    }
    
    err = 0;
    
    /* Create and apply each MPO one by one */
    for (i = 0; i < modifypmpo.sz; i++) {
        sprintf(osnm,"osm_%d",i);
        tntNodeArraysNamedLoad(initfile,os,osnm);
        
        /* create a product or sum MPO depending on the flag */
        if (modifypmpo.vals[i]) {
            sprintf(snnm,"snm_%d",i);
            tntIntArraysNamedLoad(initfile,sitenum,snnm);
            mod = tntMpsCreateProductMpo(L,&os,&sitenum);
            tntIntArrayFree(&sitenum);
        } else {
            sprintf(ospnm,"osparamm_%d",i);
            tntComplexArraysNamedLoad(initfile,osparam,ospnm);
            mod = tntMpsCreateMpo(L,NULL,NULL,NULL,&os,&osparam);
            tntComplexArrayFree(&osparam);
        }
        tntNodeArrayFree(&os);
        
        /* Multiply the MPO with the start state */
        tntPrintf("Applying operator O_%d to the starting base state\n",i+1);
        err += tntMpsMpoProduct(wf_start,mod,chi);
    }
    
    /* Renormalise the wave function if appropriate */
    tntIntParamsLoad(initfile, mod_renorm);
    if (mod_renorm) tntMpsOrthNorm(wf_start, 0);
    else            tntMpsOrth(wf_start, 0);
    
    /*************************************************************************************************/
    /*                                      Perform time evolution                                   */
    /*************************************************************************************************/  
    err = 0; tstep = 0;
    /* Load number of time steps and flags  */
    tntIntParamsLoad(initfile, numsteps, tbigstep, calc_ol_gs, calc_ol_is, calc_energy);
    if (numsteps > 0) {
        
        /* Load all the parameters needed for time evolution */
        tntNodeArraysNamedLoad(initfile,nnl,"nnlt",nnr,"nnrt",os,"ost");
        tntComplexArraysNamedLoad(initfile,nnparam,"nnparamt",osparam,"osparamt");
        tntComplexArraysNamedLoad(initfile,nnpsc,"nnparamt",ospsc,"osparamt");
        tntComplexArraysNamedLoad(initfile,nntsc,"nnparamt_time",ostsc,"osparamt_time");
        tntComplexParamsLoad(initfile, dt);
        
        /* Make a copy of the start state to evolve */
        wf_ev = tntNetworkCopy(wf_start);
        
        /* Free the other wave functions if not required */
        if (!calc_ol_gs && precision > 0) tntNetworkFree(&wf_gs);
        if (!calc_ol_is) tntNetworkFree(&wf_start);
        
        tntPrintf("Starting time evolution\n");

        /* Run the simulation for all the time steps. */
        err = 0.0;
        for (tstep = 0; tstep <= numsteps; tstep++) {
            
            /* Check whether the propagator needs to be redefined after a change in time varying parameter */
            if (checkParams(&nntsc, &ostsc, tstep)) {
                setParams(&nnparam, &nntsc, &nnpsc, &osparam, &ostsc, &ospsc, tstep);
                prop = tntMpsCreatePropST2sc(L, dt, &nnl, &nnr, &nnpsc, &os, &ospsc);
                changed_params = 1;
            }
            
            if (tstep) err += tntMpsPropST2scProduct(wf_ev, prop, chi);
            
            /* Calculate expectation values every tbigstep */
            if (0 == tstep%tbigstep || tstep == numsteps) {
                tntPrintf("Calculating time evolution, %2.0f%% complete\n",(100.0*tstep)/(1.0*numsteps));

                bigstep++;
                
                tntMpsExpecOutput(wf_ev, &ExOp, 0, 0, 1, savename, bigstep);
                if (calc_ol_gs) {
                    ol = tntComplexToReal(tntMpsMpsProduct(wf_gs,wf_ev));
                    tntDoubleParamsNamedUpdate(savename, bigstep, ol, "ol_gs");
                }
                if (calc_ol_is) {
                    ol = tntComplexToReal(tntMpsMpsProduct(wf_start,wf_ev));
                    tntDoubleParamsNamedUpdate(savename, bigstep, ol, "ol_is");
                }
                if (calc_energy) {
                    if (changed_params) { /* check whether Hamiltonian has changed since last big step */
                        if (0 != tstep) tntNetworkFree(&H);
                        H = tntMpsCreateMpo(L, &nnl, &nnr, &nnpsc, &os, &ospsc);
                        changed_params = 0; /* reset flag */
                    }
                    E = tntComplexToReal(tntMpsMpoMpsProduct(wf_ev,H));
                    tntDoubleParamsUpdate(savename, bigstep, E);   
                }
                tntIntParamsUpdate(savename, bigstep, tstep);
                tntDoubleParamsUpdate(savename, bigstep, err);    
                t.re = tstep*dt.re;
                t.im = tstep*dt.im;
                tntComplexParamsUpdate(savename, bigstep, t);
            }
        }

        tntNetworkFree(&prop);
        if (calc_energy) tntNetworkFree(&H);
        tntNodeArrayFree(&nnl);
        tntNodeArrayFree(&nnr);
        tntNodeArrayFree(&os);
        tntComplexArrayFree(&nnparam);
        tntComplexArrayFree(&osparam);
        tntComplexArrayFree(&nntsc);
        tntComplexArrayFree(&ostsc);
        tntComplexArrayFree(&nnpsc);
        tntComplexArrayFree(&ospsc);
        tntNetworkFree(&wf_ev);
        if (calc_ol_gs && precision > 0) tntNetworkFree(&wf_gs);
        if (calc_ol_is) tntNetworkFree(&wf_start);
    }

    tntMpsExOpFree(&ExOp);
    
    /* Finish with the TNT library */
    tntFinalize();

    return 0;
}

/* Checks whether parameters have changed since the last time step */
unsigned checkParams(tntComplexArray *nntsc, /* Temporal scaling parameters per time-step for nn terms */
                     tntComplexArray *ostsc, /* Temporal scaling parameters per time-step for os terms */
                     int tstep)              /* Current time step */
{
    
    unsigned i;             /* Loop counter */
    double x1, x2, y1, y2;  /* Real and imaginary parts of the previous and current parameter */
    
    if (0 == tstep) return 1;
    
    if (nntsc->sz) {
        for (i = 0; i < nntsc->numrows; i++) {
            x1 = nntsc->vals[i + nntsc->numrows*(tstep-1)].re;
            x2 = nntsc->vals[i + nntsc->numrows*tstep].re;
            y1 = nntsc->vals[i + nntsc->numrows*(tstep-1)].im;
            y2 = nntsc->vals[i + nntsc->numrows*tstep].im;
            
            if (x1 != x2 || y1 != y2) return 1;
        }
    }
    if (ostsc->sz) {
        for (i = 0; i < ostsc->numrows; i++) {
            x1 = ostsc->vals[i + ostsc->numrows*(tstep-1)].re;
            x2 = ostsc->vals[i + ostsc->numrows*tstep].re;
            y1 = ostsc->vals[i + ostsc->numrows*(tstep-1)].im;
            y2 = ostsc->vals[i + ostsc->numrows*tstep].im;
            
            if (x1 != x2 || y1 != y2) return 1;
        }
    }
    
    return 0;
}

/* Scales the spatially varying hamiltonian parameters by the time varying hamiltonian parameters for the current time step */
void setParams(tntComplexArray *nnparam,  /* Unscaled spatially varying terms for nn terms */
               tntComplexArray *nntsc,    /* Temporal scaling parameters per time-step for nn terms */
               tntComplexArray *nnpsc,    /* Scaled spatially varying terms for nn terms */
               tntComplexArray *osparam,  /* Unscaled spatially varying terms for nn terms */
               tntComplexArray *ostsc,    /* Temporal scaling parameters per time-step for os terms */
               tntComplexArray *ospsc,    /* Scaled spatially varying terms for nn terms */
               int tstep)                 /* Current time step */
{
    unsigned i, j, M;  /* Row and column counters */
    double x1, x2, y1, y2;  /* Real and imaginary parts of the scaling value and unscaled parameter */
    
    if (nntsc->sz && nnparam->sz) {
        M = nnparam->numrows;
        for (i = 0; i < M; i++) { /* Row counter is for terms */
            x1 =  nntsc->vals[i + M*tstep].re;
            y1 =  nntsc->vals[i + M*tstep].im;
            for (j = 0; j < nnparam->numcols; j++) { /* column counter is for sites */
                x2 = nnparam->vals[i + M*j].re;
                y2 = nnparam->vals[i + M*j].im;
                
                nnpsc->vals[i + M*j].re = x1*x2 - y1*y2;
                nnpsc->vals[i + M*j].im = x1*y2 + y1*x2;
            }
        }
    }
    if (ostsc->sz && osparam->sz) {
        M = osparam->numrows;
        for (i = 0; i < M; i++) { /* Row counter is for terms */
            x1 =  ostsc->vals[i + M*tstep].re;
            y1 =  ostsc->vals[i + M*tstep].im;
            for (j = 0; j < osparam->numcols; j++) { /* column counter is for sites */
                x2 = osparam->vals[i + M*j].re;
                y2 = osparam->vals[i + M*j].im;
                
                ospsc->vals[i + M*j].re = x1*x2 - y1*y2;
                ospsc->vals[i + M*j].im = x1*y2 + y1*x2;
            }
        }
    }
}