# This is an example script for calculating the time-evolution of an anti-neel state (|down,up,down,up,...>) under XXZ-Heisenberg Hamiltonian using W-II stepper for initial growth of bond dimension and then switching over to a tdvp 2-site stepper
# Input parameter are 
#  $1: path to executable directory
#  $2: path to model directory (use path to your models-git which is at the same location as the executables) containing besides others a model-file which declares a Heisenberg Hamiltonian
#  $3: system size, default value is 40
#  $4: value of anisotropy parameter J, default is J=1.0
#  $5: value of sz-interaction Jz, default is Jz=1.0
#  $6: value of longitudinal magn. field hz, default is hz=0.0

# at first initialize default values
L=40
J=1.0
Jz=1.0
hz=0.0

# read input and update model parameters if passed
case $# in
	3)
		L=$3
		;;
	4)
		L=$3
		J=$4
		;;
	5)
		L=$3
		J=$4
		Jz=$5
		;;
        6)
                L=$3
                J=$4
                Jz=$5
		hz=$6
                ;;
	0|1)
		echo "invalid number of input parameter, at least path to executable directory and model-file must be specified"
		exit 1
		;;
	*)
		echo "max. number of interpreted input parameters is 6, I'm going to ignore the remaining" $(( $#-6 )) "parameters"
esac

EXEC_DIR=$1
PRM_DIR=$2

# force to use only one thread
export OMP_NUM_THREADS=1
export MKL_NUM_THREADS=1

# do some clean-ups
rm -rf log/*
rm -rf teo/*
rm lattice.*
rm *.symmps

# now create two lattices with real as well as complex fundamental datatype
# short option -F defines fundamental datatype of lattice
# short option -b defines basis system used for this lattice (default is custom, containing all basis sets with only one conserved quantum number)
# short option -g defines symmetry generator for this lattice (default is Sz but can be set to any set of local operators, consistent with the number of conserved quantities supprted by the chosen basis type. Note, that in order to get information about the common symmetry generators for each basis type by running sym-create-lattice with short-options -b BASISTYPE -G where BASISTYPE is a valid local basis)
# short option -L specifies the numebr of lattice sites
# short option -d specifies the dimension of the local operators generating the basis at each lattice site (default is 2, but can be set to larger values which then are combined properly with the chosen local basis type, for instance S=1 systems can be realized by setting local dim to 3)
# short option -o defines output file in which lattice is stored
$EXEC_DIR/link/sym-create-lattice -F d -b custom -g Sz -L $L -d 2 -o lattice.real
$EXEC_DIR/link/sym-create-lattice -F c -b custom -g Sz -L $L -d 2 -o lattice.complex

# we proceed by constructing the Heisenberg hamiltonian into the complex lattice (we use the complex latticee because the time-evolution is complex for obvious reasons)
# short option -i defines lattice input file specifying the Hilbertspace on which the MPO is constructed
# short option -r defines relative path of the file containing the FSM with respect to the value set in prm_dir
# short option -m defines model filename cpontaining the desired FSM (default value is models.prm so collecting FSMs in files with these names you can drop this option)
# short option -M defines FSM identifier in model file from which MPO is constructed
# short option -p sets parameter for FSM used to create MPO
# short option -o defines output filename for a lattice containing the created MPO
# short option -n defines identifier of created MPO under which it is stored in output lattice-file
$EXEC_DIR/link/sym-create-mpo --prm_dir $PRM_DIR -i lattice.complex -r Teq0/spin/1D/ -m models.prm -M heisenberg -p "J=${J}|Jz=${Jz}|hz=${hz}" -o lattice.complex -n H

# now create anti-neel state, we do this as MPS with real fundamental datatype to demonstrate how such MPS are automatically converted into complex MPS when they are loaded into an executable with complex lattice
# short option -i defines lattice input file containing Hilbert space and local basis
# short option -t defines type of created mps, choice 'anti-neel' means that sites are initialized alternatively in eigenstate of the local symmetry generators with minimum/maximum eigenvalue, i.e., s_i=-0.5/+0.5 in our case
# short option -o defines dest state filename
$EXEC_DIR/link/sym-create-mps -i lattice.real -t anti-neel -o anti-neel.symmps

# to test if we everything worked so far we measure the Sz expectation value of the fresh state
# short option -r is a boolean flag which, if set, resets the generated datafile erasing all content which may previously have been written into
# short option -n is a boolean flag which, if set, normalizes the expectation value(s) calculated with the states norm/overlap sqrt(<bra|ket>)
# short option -i defines lattice input file containing Hilbert space and local basis
# short option -k defines ket state of scalar-product (if bra is not specified it is set to the value of ket)
# short option -o defines a set of local operators for which the scalar product is to be evaluated
# short option -t defines a tag which is used to label the obtained expectation values in the output file (we set it to 0.000000 for reasons which become clear later)
$EXEC_DIR/link/sym-scalar-product -rni lattice.real -o 'Sz' -k anti-neel.symmps -t 0.000000 -d teo
# Note, that now there is a directory teo which contains a folder local_operators in which two files: Sz_real, Sz_imag are placed containing the real and imaginary part of the evaluated scalar product. For each local operator specifed via the -o short option a new file is created, if it does not already exists. If the file already exists the behaviour of the output is controlled by the -r flag. If set then the file's content will be removed, otherwise the data will be appended where in the first column the passed value for the tag is written and the second column is reserved for the norm/overlap of the input bra/ket-states

# Now comes the tricky part, we want to a time-evolution of this state which is a product state. Therefore, we can not begin the time-evolution using single-site TDVP-solver. However, as the problem has only nearest neighbor interactions we could in principle use 2-site TDVP solver which treats nearest neighbor interactions exactly and therefore does not pickup a projection error. This changes drastically if we have longer ranged interactions. To demonstrate a way around this problem we start the time-evolution employing an MPO stepper to grow the bond-dimension before switching to 2-site TDVP

# start with MPO time-evolution using W-II MPO stepper in a second order decomposition scheme
# short option -i defines lattice input file containing Hilbert space and local basis
# short option -w defines permitted maximum bond dimension of MPS-site tensors during time-evolution
# short option -d defines truncated weight threshold of MPS-site tensors
# short option -k defines data-type of time-argument \tau: c(omplex) means real-time evolution \tau= \delta t, d(ouble) means imaginary time-evolution \tau= \mathrm i \delta t
# short option -t defines initial time of time evolution
# short option -D defines time-step size |\delta t|
# short option -T defines final time of time evolution
# short option -I defines initial state for time evolution
# short option -O defines MPO used to construct time-stepper MPO loaded from input lattice
# short option -n defines max. number of sweeps during variational application of MPO-stepper
# short option -e defines escape precision for variational application of MPO-stepper
# short option -o defines checkpoint directory in which calculated time steps (MPS) are written
# short option -X is a flag which, if set, stores the constructed time stepper MPO(s) into the lattice, can be usefull if the lattice may be re-used later to save stepper construction time
# short option -V defines output verbosity (values > 7 also yield debug outputs) which we decrease to reduce runtime-information output
$EXEC_DIR/link/sym-time-evol-WII -i lattice.complex -w 100 -d 1e-10 -k c -t 0.0 -D 0.05 -T 1.0 -I anti-neel.symmps -O H -n 10 -e 1e-10 -o teo -XV 1
# If you check the output you might notice, that the norm of the time evolved states change during time-evolution. This is an artefact of the time-stepper being non-unitary and can be compensated for by passing the flag -z which forces renormalization after evolving the state by a time step \delta T

# Finally we switch to the 2-site TDVP solver, note that the arguments are merely the same except for those specifying the local krylov solver. Note, that we can choose a larger time-step \delta t since 2-site tdvp has no Trotter error
# short option -i defines lattice input file containing Hilbert space and local basis
# short option -w defines permitted maximum bond dimension of MPS-site tensors during time-evolution
# short option -d defines truncated weight threshold of MPS-site tensors
# short option -k defines data-type of time-argument \tau: c(omplex) means real-time evolution \tau= \delta t, d(ouble) means imaginary time-evolution \tau= \mathrm i \delta t
# short option -t defines initial time of time evolution
# short option -D defines time-step size |\delta t|
# short option -T defines final time of time evolution
# short option -I defines initial state for time evolution, here we pass the last checkpoint from the prev. MPO time-evolution
# short option -O defines MPO used to perform time evolution loaded from input lattice
# short option -s defines the TDVP stepper-type, can be choosen as single site (1site) or two site (2site) stepper, note that 1site is by a factor of the local dimensionoi faster but unable to increase bond dimension of MPS and may have larger projection error if, in particular if MPS bond-dimension is small
# short option -n defines minimal number of Krylov vectors used to construct Krylov subspace when solving the local Schrödinger equation
# short option -N defines maximal number of Krylov vectors used to construct Krylov subspace when solving the local Schrödinger equation
# short option -E defines escape precision for the local Krylov solver obtained by calculating the error described by Hochbruck and Lubich
# short option -V defines output verbosity (values > 7 also yield debug outputs) which we increase it a bit to obtain more runtime-information output
$EXEC_DIR/link/sym-time-evol-tdvp -i lattice.complex -w 100 -d 1e-10 -k c -t 1.0 -D 0.1 -T 20.0 -I teo/state.1.000000 -O H -s 2site -n 3 -N 10 -E 1e-08 -o teo -V 9

# Now the time evolution with TDVP 2-site solver is running. While waiting for it to return you can launch the script meas_teo.sh in the examples directory of your symmps folder.
# This script waits for the next checkpoint to be finished and performs a scalar-product to measure the Sz-expectation values after each time-step.
# Since we have already reset the data-files above by measuring the initial anit-neel state passing the short option -r, the measure script will only append new expectation values to the existing file yielding a matrix of datasets which can be plotted, for instance, in gnuplot via the command: "p 'teo/local_operators/Sz_real' matrix u 1:(0.1*$2):w ev ::2 w image".
# You obtain an image of the total time evolution of the expectation value <ψ|S^z|ψ> where |ψ> is the initial anti-neel state