# 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.