# Example: Dynamics in the Bose-Hubbard ModelΒΆ

We want to simulate the Bose-Hubbard model, and look at some basic dynamics. So first we write our Hamiltonian that defines this model.

(1)

Equation (1) is the general expression for the Hamiltonian of the Bose-Hubbard model. For this exercise we work in the canonical ensemble (fixed particle number), and we ignore the chemical potential and write the Bose-Hubbard model as Eq. (2).

(2)

Here we have our creation and annihilation operators, , and our number operator, , as usual. Additionally, our second term is scaled by the parameter that describes site interactions, .

First we must include the necessary python libraries. The two that are
needed are the MPSPyLib, which includes how python processes the MPS matrix
inputs and format the user's Hamiltonian into the proper format for the
Fortran compiler to run. The post-process flag is an argument in the main
function and can be passed as command line argument as `--PostProcess=F`

or `--PostProcess=F`

. We have imported
two different libraries, and use the handles mps and np to refer to functions
within each library. Thus any time a function is preceded with one of these
handles, (i.e. `Ops.BuildBoseOperators()`

)

1 2 3 4 5 | ```
import MPSPyLib as mps
import numpy as np
from sys import version_info
import matplotlib.pyplot as plt
import sys
``` |

Next we need to create which operators we want to observe in our Hamiltonian.
So for the Bose-Hubbard model, we want bosonic operators, so we call
`Ops.BuildBoseOperators()`

, with a maximum of 5 particles allowed on
a given lattice site. This is the term inside the parentheses. We then want
build a custom operator: we are choosing to look at our interaction term,
which will be defined as
. The
following lines correspond to creating this operator, and setting up our
system.

20 21 22 23 24 | ```
# Build operators
Operators = mps.BuildBoseOperators(5)
Operators['interaction'] = 0.5 * (np.dot(Operators['nbtotal'],
Operators['nbtotal'])
- Operators['nbtotal'])
``` |

Then we build our Hamiltonian. We need to add two terms (or two sums) to our
Hamiltonian. So we first define a variable that will act as our Hamiltonian,
which in this case will be `H`

.

26 27 28 29 | ```
# 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)
``` |

Note that we have two different types of MPO terms present in our Hamiltonian
(here we have `'bond'`

and `'site'`

terms). This is because the summations
for each term in the Hamiltonian are different. The bond summation considers
nearest-neighbor interactions, and has different site indices,
, for each operator. MPS takes each operator fed to
the function (in this case the `'bdagger'`

, `'b'`

expression) and assigns
indices in order of the list. The other type of summation is simply over all
of the sites. This type of MPO term expects only one operator (or a single
expression of operators).

Next we define what observables we want to measure in our ground state. We can look at the number operator at each site

31 32 33 34 | ```
# Ground state observables
myObservables = mps.Observables(Operators)
# Site terms
myObservables.AddObservable('site', 'nbtotal', 'n')
``` |

Then we want to define our correlation operators, when looks at neighboring sites. We define two observables, and . This is done by the following commands:

35 36 37 | ```
# Correlation functions
myObservables.AddObservable('corr', ['nbtotal','nbtotal'], 'nn')
myObservables.AddObservable('corr', ['bdagger','b'], 'spdm')
``` |

Next we set the convergence parameters for our simulation. This is important to set the convergence for finding our ground state to a reasonable tolerance to get an accurate ground state, but should be optimized along with simulation run times.

39 40 41 42 43 44 45 | ```
# Convergence parameters
myConv = mps.MPSConvParam(max_bond_dimension=20, max_num_sweeps=6)
myConv.AddModifiedConvergenceParameters(0, ['max_bond_dimension',
'local_tol'], [50, 1E-14])
myKrylovConv = mps.KrylovConvParam(max_num_lanczos_iter=20,
lanczos_tol=1E-6)
``` |

We have modified our convergence parameters to check that our max bond dimension in our system is no larger than 50 and our local tolerance is . We set our Krylov-Convergence parameters to do a maximum of 20 Lanczos Iterations with a tolerance of for each calculation. It should be noted that computation time is largely impacted by these convergence parameters. If tolerances are set too small for a given system, computation time can be greatly increased.

In this simulation we also define a second set of observables to show that we can look at different observables for different portions of our simulation. Here we want to look at the number operator for each site as our system is evolving. So we define a new observable dictionary:

47 48 49 50 | ```
# Dynamics time evolution observables
dynObservables = mps.Observables(Operators)
# Site terms
dynObservables.AddObservable('site', 'nbtotal', 'n')
``` |

The next step is to define values for our Hamiltonian we created, assign what location and name our simulation output will be run with, and define what type of simulation we want to run. We first define our parameters and in our Hamiltonian. We also define the number of lattice sites.

52 53 54 55 56 57 | ```
# Specify constants and parameter lists
U = 10.0
t = 1.0
staticsParameters = []
L = 6
tlist = [5.0, 20.0]
``` |

Then we want to run a simulation over two different variations, where the value of U is changed at different rates. This will define our quenching function which will allow us to ramp our function value . This is done through:

60 61 62 63 64 65 | ```
# Quench function ramping down
def Ufuncdown(t, tau=tau):
return 10.0 + 2.0 * (1.0 - 10.0) * t / tau
# Quench function ramping back up
def Ufuncup(t, tau=tau):
return 1.0 + 2.0 * (10.0 - 1.0) * (t - 0.5 * tau) / tau
``` |

We then define our quenches, which will change the value of one of our parameters, in this case .

71 72 73 74 75 | ```
Quenches = mps.QuenchList(H)
Quenches.AddQuench(['U'], 0.5 * tau, min(0.5 * tau / 100.0, 0.1),
[Ufuncdown], ConvergenceParameters=myKrylovConv)
Quenches.AddQuench(['U'], 0.5 * tau, min(0.5 * tau / 100.0, 0.1),
[Ufuncup], ConvergenceParameters=myKrylovConv)
``` |

We note that `Dynamics.AddQuench()`

takes the Hamiltonian, the
variable in our Hamiltonian which we want to ramp, the range of values for
which we will apply this function, assigns a function which calculates the
new value, and the set of convergence parameters that we want to abide by.

Then we create the dictionary that will be used for our simulation. This is done through the following lines of code:

73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 | ```
staticsParameters.append({
'simtype' : 'Finite',
# Directories
'job_ID' : 'Bose_Hubbard_mod',
'unique_ID' : 'tau_'+str(tau),
'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'],
'Abelian_quantum_numbers' : [L],
# Convergence parameters
'verbose' : 1,
'logfile' : True,
'MPSObservables' : myObservables,
'MPSConvergenceParameters' : myConv,
'Quenches' : Quenches,
'DynamicsObservables' : dynObservables
``` |

Here we have three different sections of assignments. We have our simulation
references, `'job\_ID'`

which is what each simulation file will be named,
with a unique ending given by the value assigned to the key `'unique\_ID'`

.
Then we assign where these files should be saved. We want our compilation
files to be saved in the local folder `TMP`

and all of our outputs to the
local folder `OUTPUTS`

.

Next we define our system. We need to call every unassigned variable in our
Hamiltonian, and also specify our system size. We do this through the key
assignments for `'L'`

, our system size, `'t'`

and `'U'`

the variables
discussed in our Hamiltonian in Eq. (2).

Next we have the section which we can define any symmetries or quantum
numbers for our system. These are set by the keys `'Abelian_generators'`

and `'Abelian_quantum_numbers'`

.

And finally we describe which convergence parameters we want to use for each
section of our analysis, and which observables we want to measure for each
portion of our analysis. This is done through the lines *89-92*.

Next we have the primary function calls that convert this python dictionary
to the fortran files that will be used to actually run the MPS algorithms.
This is done by calling the dictionary we just defined, and passing the
operators, the Hamiltonian, and the option to post process data for our
simulation. Then we call the function `tools.runMPS()`

to run the
simulation.

95 96 97 98 99 100 101 102 | ```
# Write Fortran-readable main files
MainFiles = mps.WriteFiles(staticsParameters, Operators, H,
PostProcess=PostProcess)
# Run the simulations and quit if not just post processing
if(not PostProcess):
mps.runMPS(MainFiles, RunDir='./')
``` |

Then to check that this simulation works correctly, we can use the following code write our time data for version of our simulation. Here we have run two simulations with and . This can be done through the following commands:

100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 | ```
# Postprocessing and plotting
# ---------------------------
Outputs=mps.ReadStaticObservables(staticsParameters)
DynOutputs=mps.ReadDynamicObservables(staticsParameters)
ii = 0
le_list = [[], []]
for t in tlist:
myfile='outmod'+str(t)+'.dat'
hfile=open(myfile,'w')
for p in DynOutputs[ii]:
hfile.write('%30.15E'%(p['time']) + '%30.15E'%(p['U']) +
'%30.15E'%(p['Loschmidt_Echo'].real) +
'%30.15E'%(p['Loschmidt_Echo'].imag) +
'%30.15E'%(p['bond_dimension']) + '\n')
le_list[ii].append(abs(p['Loschmidt_Echo'])**2)
hfile.close()
``` |

In the last part we collected already the information of the Losschmidt echo over time. The remaining lines plot the overlap with the initial state over time for the two different quench times:

126 127 128 129 130 131 132 133 134 135 136 137 138 139 | ```
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.plot(np.linspace(0, 1, len(le_list[0])), le_list[0], 'r-',
label="tq=5")
plt.plot(np.linspace(0, 1, len(le_list[1])), le_list[1], 'g-',
label="tq=20")
plt.xlabel(r"\textbf{Time in percent} " r"$t / t_q$", fontsize=16)
plt.ylabel(r"\textbf{Loschmidt Echo} "
r"$| \langle \psi(t) | \psi(0) \rangle |^2$",
fontsize=16)
plt.legend(loc="lower left")
plt.show()
``` |