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 t, is known as the tunneling term of the model, a large t promotes particles tunneling between nearest neighbor sites. The second term, whose magnitude is controlled by U, 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)\hat{H} = -t \sum_{i=1}^{L-1} (\hat{\mathrm{b}}^{\dagger}_i
 \hat{\mathrm{b}}_{i+1}
 + \hat{\mathrm{b}}_i \hat{\mathrm{b}}^{\dagger}_{i+1})
 + \frac{1}{2} U \sum_{i=1}^{L} \hat{\mathrm{n}}_i
 (\hat{\mathrm{n}}_i-\hat{I})
 - \mu\sum_{i=1}^{L} \hat{\mathrm{n}}_i\,,

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 \mu, since in the canonical ensemble this only introduces an overall shift in the energy of the system, the \frac{\mu}{U} 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 0 \leq t \leq 0.4. 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 \mu for the rest of our simulations, again, at a particular value of t. Finally mulist, tinternal, and depletion are appended to respective numpy arrays in alldata. This process in repeated for each value of t. 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

../_images/Tut2_BoseHubbard_Depletion.jpg

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.