# Example: Static Bose Hubbard ModelΒΆ

Next we handle the more difficult task of exploring the two parameter phase space of the Bose-Hubbard Hamiltonian. In one spatial dimension the Bose-Hubbard Hamiltonian takes the form shown in Eq. (1). The first term, whose magnitude is controlled by , is known as the tunneling term of the model, a large promotes particles tunneling between nearest neighbor sites. The second term, whose magnitude is controlled by , describes repulsive on site interaction, finally the last term in the sum accounts for the energy associated with the chemical potential of the system.

(1)

As usual we begin my importing necessary libraries. There are a couple new
libraries in the lines below. From `matplotlib`

we import `cm`

, this
allows us to use color maps in our plots. This of course requires the
installation of `matplotlib`

, see Sec. Dependencies. We also
import the time library, we use this to time our simulation, allowing us
to get instant and exact feedback about how long our simulations take to run.

1 2 3 4 5 6 | ```
import MPSPyLib as mps
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from time import clock
import sys
``` |

Towards the end of the script we define a function for plotting the results
of our simulations, The function is useful because its sets the limits of
plots using the data provided to it, but not that defining it requires
`matplotlib`

to be installed, along with *LaTeX*, and additional fonts as
described in the dependencies section, see Sec. Dependencies.

118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 | ```
def plotIt(jvalues, muvalues, dependentvalues):
"""
Scatter plot for the phase diagram of the Bose Hubbard model.
**Arguments**
jvalues : floats
x-values of points
muvalues : floats
y-values of points
dependentvalues : floats
value for coloring points
"""
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.scatter(jvalues, muvalues, c=dependentvalues, cmap=cm.jet)
plt.xlim((np.min(jvalues), np.max(jvalues)))
plt.ylim((np.min(muvalues), np.max(muvalues)))
plt.xlabel(r"\textbf{tunneling} " r"$t/U$", fontsize=16)
plt.ylabel(r"\textbf{chemical potential} " r"$\mu/U$", fontsize=16)
cbar = plt.colorbar()
cbar.set_label(r"\textbf{Quantum Depletion}", fontsize=16)
plt.show()
``` |

The first time we run the file, we set the post process flag to `False`

via the command line argument `--PostProcess=False`

which is passed to
the main function:

9 | ```
def main(PostProcess=False):
``` |

Inside the main function, we set an initial time, this is stored in seconds.

21 | ```
t0 = clock()
``` |

We build the appropriate operators and construct the interaction operator
using the built in linear algebra of python. Note that computationally it
is not necessary to compute the term involving , since in the
canonical ensemble this only introduces an overall shift in the energy of
the system, the term
is computed by subtracting the ground state energy `Output['energy']`

across simulations with different numbers of particles.

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 | ```
# Build operators
Operators = mps.BuildBoseOperators(6)
Operators['interaction'] = 0.5 * (np.dot(Operators['nbtotal'],
Operators['nbtotal'])
- Operators['nbtotal'])
# Define Hamiltonian MPO
H = mps.MPO(Operators)
H.AddMPOTerm('bond', ['bdagger','b'], hparam='t', weight=-1.0)
H.AddMPOTerm('site', 'interaction', hparam='U', weight=1.0)
# ground state observables
myObservables = mps.Observables(Operators)
# Site terms
myObservables.AddObservable('site', 'nbtotal', 'n')
# correlation functions
myObservables.AddObservable('corr', ['nbtotal', 'nbtotal'], 'nn')
myObservables.AddObservable('corr', ['bdagger', 'b'], 'spdm')
``` |

We let OSMPS set the maximum bond dimension using default criteria,
see Sec. Variational ground state search or `convergence.MPSConvParam`

for details.

41 | ```
myConv = mps.MPSConvParam(max_num_sweeps=7)
``` |

We set `U=1`

, and generate lists for parameters we will change across
simulations. The first is the list of values our tunneling parameter will
range over, `tlist`

, we will sample 20 points from the interval
. The system size is set to `L=10`

and the number
of particles is set to range from `1`

to `11`

, with 11 sampling points.

43 44 45 46 47 | ```
U = 1.0
tlist = np.linspace(0,0.4,20)
parameters = []
L = 10
Nlist = np.linspace(1, 11, 11)
``` |

Using these lists of numbers we generate all of the necessary parameter lists for our simulation.

49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 | ```
for N in Nlist:
for t in tlist:
parameters.append({
'simtype' : 'Finite',
# Directories
'job_ID' : 'Bose_Hubbard_statics',
'unique_ID' : 't_' + str(t) + 'N_' + str(N),
'Write_Directory' : 'TMP/',
'Output_Directory' : 'OUTPUTS/',
# System size and Hamiltonian parameters
'L' : L,
't' : t,
'U' : U,
# Specification of symmetries and good quantum numbers
'Abelian_generators' : ['nbtotal'],
# Define Filling
'Abelian_quantum_numbers' : [N],
# Convergence parameters
'verbose' : 1,
'logfile' : True,
'MPSObservables' : myObservables,
'MPSConvergenceParameters' : myConv
``` |

If `PostProcess=False`

we run the simulations, check the time in seconds,
and subtract our original time to measure the amount of time `mps.runMPS`

took to run.

72 73 74 75 76 77 78 79 | ```
# Write Fortran-readable main files
MainFiles = mps.WriteFiles(parameters, Operators, H,
PostProcess=PostProcess)
# Run the simulations and quit if we are not just Post
if(not PostProcess):
mps.runMPS(MainFiles, RunDir='./')
print('Elapsed time', clock() - t0)
``` |

Once we have run the simulation once already, we can run the python file
again with `PostProcess=True`

and enter the following `if`

statement.
We initialize the list of `numpy`

arrays `alldata`

, and use
`mps.ReadStaticObservables(parameters)`

to get all of the output data from
our simulations and store it in `Outputs`

. We enter a `for`

loop which
loops through the values in tlist, each time we pass through the loop we use
`mps.GetObservables(Outputs,'t',t)`

to grab all the output data which has
`'t':t`

, we store these dictionaries in `Outputs2`

. We then make
`len(Nlist)`

number of copies of this particular `t`

value, one for each
of the several number of simulations done at this values of `'t'`

. Within
this `for`

loop we also initialize three empty lists `mulist`

,
`energylist`

, and `depletionlist`

. In the `for`

loop on line 89, we ask
for the single particle density matrix associated with each output dictionary
in `Outputs2`

. We use the built in linear algebra provided by numpy to
compute the quantum depletion of the state [PO56], by computing
its eigenspectrum, this is stored in `depletion`

. Next we append this to
`depletionlist`

, ask for the energy of the state, store it in `energy`

,
and append it to `energylist`

. Finally we ask python to print the status
of `'converged'`

.

On line 105 we set the first value of `mulist`

to the energy of the 1
particle simulation, the essential assumption here is that the zero particle
state has zero energy. With `energylist`

in hand we can compute
for the rest of our simulations, again, at a particular value of .
Finally `mulist`

, `tinternal`

, and `depletion`

are appended to respective
numpy arrays in `alldata`

. This process in repeated for each value of
. We end by plotting the results of these simulations on line 114.
So that if line 114 is to run without error, you must install `matplotlib`

,
*LaTeX*, and additional fonts as described in the dependencies section, see
Sec. Dependencies.

81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 | ```
# Postprocessing and plotting
# ---------------------------
alldata = [np.array([]), np.array([]), np.array([])]
Outputs = mps.ReadStaticObservables(parameters)
for t in tlist:
mulist = []
energylist = []
depletionlist = []
Outputs2 = mps.GetObservables(Outputs, 't', t)
tinternal = t * np.ones(len(Nlist))
for Output in Outputs2:
spdm = Output['spdm']
spdmeigs, U = np.linalg.eigh(spdm)
maxeig = np.max(spdmeigs)
eigsum = np.sum(spdmeigs)
depletionlist.append(1 - (maxeig / eigsum))
energylist.append(Output['energy'])
print(Output['t'])
print('Simulation Converged?:', Output['converged'])
mulist.append(energylist[0])
for i in range(1, len(energylist), 1):
mu = energylist[i] - energylist[i - 1]
mulist.append(mu)
alldata[0] = np.append(alldata[0], tinternal)
alldata[1] = np.append(alldata[1], mulist)
alldata[2] = np.append(alldata[2], depletionlist)
plotIt(alldata[0], alldata[1], alldata[2])
``` |

The plot should then look like

Alternatively, we can export data in neatly formatted columns by replacing line 111 with the following, which follows the same syntax described in section Example: Static Ising Model.

```
outfilename = 'BoseStatics.dat'
outfile = open(outfilename, 'w')
for i in range(len(alldata[0])):
outfile.write('%30.15E'%(alldata[0][i]) + '%30.15E'%(alldata[1][i]) +
'%30.15E'%(alldata[2][i]) + '\n')
outfile.close()
```

Now we can continue in the next example with a real time evolution.