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)\hat{H} = -t \sum_{\langle i, j \rangle}
 \left[\hat{b}_i^{\dagger} \hat{b}_j + \mathrm{h.c.} \right]
 + U \sum_{i} \frac{\hat{n}_i \left(\hat{n}_i - 1 \right)}{2}
 - \mu \sum \limits_i \hat{n}_i

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)\hat{H} = -t \sum_{\langle i,j\rangle}
 \left[\hat{b}_i^{\dagger} \hat{b}_j + \mathrm{h.c.} \right]
 + U \sum_{i} \frac{\hat{n}_i \left(\hat{n}_i - 1 \right)}{2}

Here we have our creation and annihilation operators, \hat{b}_i^\dagger \text{ and } \hat{b}_j, and our number operator, \hat{n}_i = \hat{b}_i^\dagger \hat{b}_i, as usual. Additionally, our second term is scaled by the parameter that describes site interactions, U.

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 \frac{1}{2} \left( \hat{n} \cdot \hat{n}- \hat{n} \right). 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, i \text{ and } j, 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, \hat{n}_{i} \hat{n}_{j} and \hat{b}^{\dagger}_{i} \hat{b}_{j}. 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 10^{-6}. We set our Krylov-Convergence parameters to do a maximum of 20 Lanczos Iterations with a tolerance of 10^{-6} 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 U and t 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 U. 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 U.

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 t=5.0 and t=10.0. 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()