Yeah. That sounds like a really good deal. But I got a better one. How about... I give you the finger... and you give me my phone call?
QMBDQMBD Group Meeting 1/6/2020This is first introduction to using the toolkit for groundstates and time evolution of finite systems, using DMRG and TEBD, using the Bose-Hubbard model. The idea for this calculation comes from https://arxiv.org/abs/1310.0281. PreliminariesIn this tutorial, text that is in monospace font is text that will appear in the terminal. Commands that you need to type in will start with a $ sign. In some cases, the expected output will be written as well; this won't start with a $ sign. For example, $ echo hello hello If you are not familiar with the command line shell, then you should try this - type ' https://dev.to/awwsmm/101-bash-commands-and-tips-for-beginners-to-experts-30je#the-basics https://www.freecodecamp.org/news/basic-linux-commands-bash-tips-you-should-know/ Getting startedEveryone is sharing the same login directory, so firstly you should make a directory for your own use, perhaps using your name. mkdir <name> cd <name> Constructing the lattice fileTo get started with a calculation using the Matrix Product Toolkit, the first step is to construct what is known as a lattice file. This is a file that contains a description of the lattice (i.e. the degrees of freedom at each site), as well as associated operators, such as creation/annihilation operators, symmetry transformations, and short-cut operators that we can use to help construct a Hamiltonian. The lattice file is constructed by running a program to generate it (this is a short C++ program that is in the
$ bosehubbard-u1-trap Matrix Product Toolkit version (git) usage: bosehubbard-u1-trap [options] Allowed options: --help show this help message -N [ --NumBosons ] arg Maximum number of bosons per site [default 5] -t [ --trapwidth ] arg Trap width [default 50] -o [ --out ] arg output filename [required] Description: U(1) Bose-Hubbard model Author: IP McCulloch <ianmcc@physics.uq.edu.au> Operators: H_J - nearest-neighbor hopping H_U - on-site Coulomb repulsion N*(N-1)/2 H_V - nearest-neighbour Coulomb repulsion H_trap - harmonic trap, n(x)*(x-(L-1)/2)^2 Functions: (none) This shows that the only required option is the output filename. The other parameters we can leave as the defaults. The help message also shows some information about the Hamiltonian operators that will be constructed. In this case, there is a hopping term {$b^\dagger_i b_{i+1} + h.c.$}, local and nearest-neighbor Coulomb interactions, and a harmonic trap. Lets generate the lattice file.
$ bosehubbard-u1-trap -o lattice The output file is called $ ls lattice Groundstate wavefunctionNow we will obtain the groundstate wavefunction. The default trap width is 50 sites, so we need to construct a wavefunction with the same number of sites. For this tutorial, we will start with a random wavefunction, and use that as input for the DMRG algorithm to obtain the groundstate. To construct the random initial state, we use $ mp-random -l lattice -u 50 -q 20 -o psi The parameters here are Now we run the
Finally, the command we will run is $ mp-dmrg -w psi -H "lattice:H_J + 4*H_U + 0.05*H_trap" -m 100x20 Note: the quotes around the Hamiltonian are important, otherwise the shell will get confused about the spaces, * and other special characters in the Hamiltonian expression. This should only take a few seconds to run. There will be a lot of text printed on the terminal, but near the end it will show the final energy. Examining the groundstateThe
We can also use
Time evolutionThe next stage of the tutorial is to calculate the non-equilibrium real-time evolution if we quench the trap frequency. This induces breathing modes which can be probed using the 2nd moment of the centre of mass. There are several different algorithms for the time evolution of a tensor network, which basically amount to different ways to approximate the time evolution operator {$\exp[-iHt]$}. We will use the simplest algorithm, which is known as TEBD. This utilises the Lie-Trotter-Suzuki decomposition of the Hamiltonian, to decompose it into a sequence of 2-body unitary gates. This depends on the terms in the Hamiltonian being at most nearest-neighbour. The toolkit can decompose a Hamiltonian operator into a sum of 2-body gates automatically, in most cases. To use the
Lets copy the wavefunction to a new name, for consistency with the time-evolved wavefunctions. $ cp psi psi-t0.0 And construct the time evolution up to total time T=80, $ mp-tebd -w psi-t0.0 -H lattice:"H_J + 4*H_U + 0.04*H_trap" -t 0.2 -s 5 -n 400 -d 1e-8 -o psi-t Note that we changed the value of the trap strength from 0.05 to 0.04, which is a moderate strength quench that will excite several breathing modes.
This command will produce a sequence of wavefunctions named Properties of the time-evolutionThere are many things we could do to examine the time-evolved states. The width of the distribution is a good starting point. This is the operator $ for t in $(seq 0 80); do mp-expectation -r psi-t$t.0 lattice:"com(0)^2" ; done You could copy and paste this into your favourite graphing program, or save it to a file for later analysis. Note that we used the option The output should look something like: which shows that there are at least two oscillation frequencies that have been excited (this shape continues to at least time 120). AdditionalSome extra things that we can calculate:
There is a technical detail that we didn't mention yet - the lattice file is actually designed for an infinite system; the trapping potential is actually periodic and repeats every 50 sites. The $ mp-expectation psi lattice:"fsup(0,50, H_J + 4*H_U + 0.05*H_trap)" 9.306689016882382 Here we used the The above command will give the energy of the groundstate. We can also calculate the energy variance, {$\langle (H-E)^2 \rangle$}, via $ mp-expectation psi lattice:"(fsup(0,50, H_J + 4*H_U + 0.05*H_trap) - 9.306689016882382)^2" 8.738852557144128e-07 This shows that the energy variance is quite small, meaning that 100 basis states is enough to get a quite accurate wavefunction. Contrary to what you might expect, the error in the energy in MPS calculations is linear in the variance, so if we repeated the calculation with an increased basis size so that the variance decreased by some factor {$K$}, then the difference between the variational energy and the exact energy will decrease by the same factor (this is strictly true only for an infinite system -- for a finite system there is also a quadratic term). Using this, it is possible to do an extrapolation to zero variance and get an estimate of the exact groundstate energy.
We can also calculate the energy of the {$t=0$} with respect to the quenched Hamiltonian, $ mp-expectation psi lattice:"fsup(0,50, H_J + 4*H_U + 0.04*H_trap)" 5.217104818642688 In principle the energy should stay constant throughout the time evolution. Using TEBD, this will only be approximately true, because the wavefunction is projected to a smaller basis size after applying each 2-body gate. (There are more sophisticated time evolution algorithms such as TDVP that preserve unitarity exactly). If we calculate the energy of the {$t=80$} state, $ mp-expectation psi-t80.0 lattice:"fsup(0,50, H_J + 4*H_U + 0.04*H_trap)" 5.214969684653424 -4.206191323493645e-16i we see it has changed a bit, and there is also a numerically negligible imaginary part. So the time evolution isn't particularly accurate -- with some more computation time we could redo it with a smaller cutoff for the density matrix. (If we used a cutoff of {$10^{-10}$} then the energy is preserved to within {$10^{-4}$}, but the time evolution takes about 17 minutes rather than ~2 minutes.)
We can calculate the overlap of two different wavefunctions using $ mp-overlap psi-t1.0 psi-t0.0 #mp-overlap psi-t1.0 psi-t0.0 #Date: Mon, 01 Jun 2020 19:16:12 +1000 #size of wavefunction is 50 sites #real #imag #magnitude #argument(deg) 0.42809318292897 -0.80803786874454 0.91443368846269 -62.085649532926 This is the expectation value {$\langle \: \psi(t=1.0) \: \vert \: \psi(t=0) \: \rangle = \langle \: \psi \: \vert \: \exp[iH] \: \vert \: \psi \: \rangle$}, where {$H$} is the quenched Hamiltonian. {$\psi$} isn't an eigenstate of the quenched Hamiltonian so we cannot simplify this expression any further, but as the quench is fairly small we might expect that we could approximate the result as {$\exp i E$} = 0.48169 - 0.87634i. This expectation value should be time-translation-invariant, so we could also check eg {$\langle \psi(t=80.0) \: \vert \: \psi(t=79.0) \rangle$}. $ mp-overlap psi-t80.0 psi-t79.0 #mp-overlap psi-t80.0 psi-t79.0 #Date: Mon, 01 Jun 2020 19:23:23 +1000 #size of wavefunction is 50 sites #real #imag #magnitude #argument(deg) 0.42781229345049 -0.80962036510645 0.91570109425645 -62.147535109296 which is correct to just under 3 significant figures. (If we improve the accuracy of the time evolution by setting the cutoff to {$10^{-10}$}, the overlap improves accuracy to nearly 5 significant figures.} |