notebook.community

Edit and run

Stochastic Series Expansion: Example!

We can also write a very basic SSE code for the spin-1/2 Heisenberg model in 1 or 2 dimensions. Naively, we might think this Hamiltonian has a sign problem, because:

\begin{align*} \hat{H} =& \sum_{\langle i,j \rangle } \hat{S}_i \cdot\hat{S}_j \\ =& \sum_{\langle i,j \rangle } \hat{S}^z_i \cdot \hat{S}^z_j + \frac{1}{2}\left( \hat{S}^+_i \hat{S}^-_j + \hat{S}^-_i \hat{S}^+_j \right) \end{align*}

We can see that the first term (the diagonal one) might generate negative amplitudes, but this is easy to fix! Just add a factor of 1/4 to the entire Hamiltonian. Now all matrix elements are positive, and life is good. We can continue, and write some code to insert diagonal operators.





In [1]:



    


import random
import math


random.seed()
sites = 16 ** 2
#number of time slices in the simulation in total - can be adjusted
L = sites / 2 
bondList    = [ 0 ] * (2 * sites + 1) #2D model
startConfig = [ 0 ] * sites
beta = 16.0 

#bondList[0] corresponds to the identity
for i in range(1,sites + 1):
    bondList[ i ]         = ( i - 1, (i)%16 + 16*(( i - 1)/16) )
    bondList[ i + sites ] = ( i - 1, (i - 1 + 16)%sites )
    if( (i - 1)%2 ):
        startConfig[ i - 1 ] = True
    else:
        startConfig[ i - 1 ] = False
operators = [ 0 ] * L
sweeps = 100
n = 0

The list operators will contain information about what operator (if any) we have inserted at some imaginary time slice. Since we can only insert one operator at each slice, we only need one list. We can uniquely encode which operator we insert using two numbers, a and b:

operator = b * 2 + a - 1

So that a = operator % 2 + 1 and b = operator/2. When a = b = 0, this is the identity. Otherwise, b encodes which bond the operator acts on, and a encodes whether it's diagonal or not. The sweeps variable tells us how many times we will traverse the operator list, doing updates, until we are satisfied we have equilibrated.

Diagonal Updates

During a diagonal update, we traverse the operator list and look for opportunities to insert or remove a diagonal operator. This is the only time we use the Metropolis updating scheme. If we encounter an off-diagonal operator, we simply update our stored spin state and move on to the next time slice. Similarly, if we reject a move (to either remove or add an operator), we just keep going to the next time slice. Remember, in order to insert a diagonal operator, the spins must be antiparallel! These diagonal operators change the weight of a path in the path integral, but do not change the path itself. DO NOT use SSE if the diagonal operators are the most important term!





In [2]:



    


def diagonalUpdate( startConfig, bondList, operator, n, beta ):
    time = 0
    while time < len(operator):
        #we should keep operator long enough that many of its elements are the identity
        if len(operator) < (4.0/3.0)*n:
            operator.extend( [0] * int(math.ceil(n/3.0)) )
            
        if( ( operator[ time ] > 0 ) and ( ( operator[ time ]%2 + 1 ) == 1 ) ): 
            #we have a diagonal operator at this slice 
            prob = ( 2.0 * float( len(operator) - n + 1 ) ) / ( float(len(bondList) - 1) * beta) 
            select = random.random()
            if ( (prob < 1.0 and select < prob) or prob > 1.0 ):
                operator[ time ] = 0 
                n               -= 1 
                
        elif( operator[ time ] == 0 ): #we have an identity at this slice 

            #have to propagate the spins to this point
            currentConfig = [ 0 ] * len( startConfig )
            currentConfig[ : ] = startConfig[ : ]
            for past in range(time):
                a = operator[ past ] % 2 + 1
                b = operator[ past ] / 2
                if( a == 2 ): #off diagonal
                    currentConfig[ bondList[b][0] ] = not (currentConfig[ bondList[b][0] ]) 
                    currentConfig[ bondList[b][1] ] = not (currentConfig[ bondList[b][1] ])
                    
            #now we can see if a diagonal update is possible
            testBond = random.randrange( 1,len(bondList) )
            if( currentConfig[ bondList[ testBond ][ 0 ] ] ^ currentConfig[ bondList[ testBond ][ 1 ] ] ):
                #spins must be antiparallel
                prob   = beta * float(len(bondList) - 1) / (2.0 * float( len(operator) - n ) )
                select = random.random()
                if ( (prob < 1.0 and select < prob) or prob > 1.0 ):
                    operator[ time ] = 2 * testBond
                    n               += 1 
        time += 1
    return n

These diagonal updates control the total number of operators in the list, and therefore also the energy. Here is an example of an operator list before and after an diagonal update:





In [3]:



    


from IPython.display import SVG
SVG("diagonalupdate.svg")



    


















    
Out[3]:







  
  
    
  
  
    
      
        image/svg+xml

To get off-diagonal operators in the list, we will need a loop update. In principle we could use some other type of off-diagonal update, but in practice these do not work very well (very unlikely to insert operators). The first part of performing the update is constructing the linked vertex list.

Linked Vertex Lists

We can represent the list of operators as a set of vertices. Each operator in this model has 4 vertices, but if we write down 3- or 4-site operators, we can have more vertices. Each operator has 4 "legs" which we'll connect to each other to form a closed loop. Then, we will flip all the operators on this loop. Since the total number of operators remains unchanged, we haven't changed the energy at all and this is ok to do. First, we have to construct the list of linked vertices. We'll use the number -1 for vertices that aren't connected to anything (they are the identity). The list should be bidirectional - it shouldn't matter what direction we traverse it in - so we will also require that vertexList[ vertexList[v] ] = v.





In [4]:



    


def constructVertexList( operators, bondList, sites ): 
    vertexList  = [ -1 ]*(4 * len(operators))
    lastVertex  = [ -1 ] * sites #which was the last vertex we saw at this site
    firstVertex = [ -1 ] * sites #which was the first vertex we saw at this site
    for time,op in enumerate(operators):
        if (op > 0):
            v0    = 4 * time #0-th leg of the operator
            b     = op/2
            sitei = bondList[ b ][ 0 ] 
            sitej = bondList[ b ][ 1 ] 
            v1    = lastVertex[ sitei ] 
            v2    = lastVertex[ sitej ]
            if( v1 > -1 ):
                vertexList[ v1 ] = v0
                vertexList[ v0 ] = v1
            else:
                firstVertex[ sitei ] = v0
            if( v2 > -1 ):
                vertexList[ v2 ] = v0 + 1
                vertexList[ v0 + 1 ] = v2
            else:
                firstVertex[ sitej ] = v0 + 1
            lastVertex[ sitei ] = v0 + 2
            lastVertex[ sitej ] = v0 + 3
    #we must connect vertices across the periodic imaginary time boundary
    for site in range(sites):
        f = firstVertex[ site ]
        l = lastVertex[ site ]
        if( f != -1 ):
            vertexList[ f ] = l 
            vertexList[ l ] = f 
    return vertexList,firstVertex

We need to carry firstVertex along because we can change our starting spin state as well (it will be a better guess for the groundstate/thermal average). Now we can perform the full loop update!

Below, you can see an example picture of an operator list, complete with a(p), b(p), operator(p), and the vertex list that would be generated. 1s have been added to the p, v, and vertexList to make it a little more obvious what's going on. If this isn't helpful, just subtract a 1 from all those variables. All vertex legs are labelled to help identify how things are connected.





In [5]:



    


from IPython.display import Image
i = Image("vertexlistsmall.png")
i



    


















    
Out[5]:







  
  
    
      
        image/svg+xml
        
        
      
    
  
  
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    p123456789101112
    a(p)102210012021
    b(p)204630026047
    op(p)409136004130914
    l=0[ 1] 31[ 5]  0[ 9] 43[13] 35[17]  4[21]  0[25]  0[29]  3[33] 15[37]  0[41] 20[45] 36
    l=1[ 2] 32[ 6]  0[10] 44[14] 47[18] 11[22]  0[26]  0[30] 19[34] 16[38]  0[42] 12[46] 48
    l=2[ 3] 29[ 7]  0[11] 48[15] 33[19] 30[23]  0[27]  0[31]  1[35] 13[39]  0[43]  9[47] 14
    l=3[ 4] 17[ 8]  0[12] 42[16] 34[20] 41[24]  0[28]  0[32]  2[36] 45[40]  0[44] 10[48] 46
    [v+1], vertexList[v]+1:
    1
    2
    3
    4
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    29
    30
    31
    32
    33
    34
    35
    36
    41
    42
    43
    44
    45
    46
    47
    48

Loop Update

We use -2 to indicate that we traversed and flipped this vertex, and -1 to show that we traversed it. We ignore all vertices with vertexList entries less than 0, since we can only traverse a loop once.





In [ ]:



    


def loopUpdate( operators, bondList, startConfig, sites ):

    legs = [1,0,3,2] #which leg is next to which
    Lists = constructVertexList( operators, bondList, sites )
    vertexList = Lists[ 0 ]
    firstVertex = Lists[ 1 ]
    while( (vertexList.count( -1 ) + vertexList.count( -2 ) ) < len(vertexList) ):
        #There's a 50% chance we flip the loop
        flip = (random.random() > 0.5)
        loopStart = 0
        for v, vv in enumerate( vertexList ):
            if( vv > -1 ):
                loopStart = v
                break #pick first valid vertex to start tracing the loop
        v = loopStart
        while( v < len(vertexList) ): 
            time = v/4
            a = operators[time]%2 + 1
            b = operators[time]/2
            vleg = 4*time + legs[ v%4 ] 
            vnext = vertexList[ vleg ]
            if flip: 
                #turn offdiagonal operators into diagonal
                if a == 2:
                    operators[ time ] = 2 * b 
                #turn diagonal operators into offdiagonal
                if a == 1:
                    operators[ time ] = 2 * b + 1 
                vertexList[ v ] = -2
                vertexList[ vleg ] = -2
            else: 
                vertexList[ v ] = -1
                vertexList[ vleg ] = -1
            
            v = vnext 
            if( v == loopStart ):
                break
    #update stored state
    #if a site was never touched by a loop, we can flip its spin
    #with probability 1/2
    #we also have to make sure things are still continuous across
    #the imaginary time boundary
    for site in range(sites):
        v = firstVertex[ site ]
        if( v == -1 and (random.random() > 0.5 )):
            startConfig[ site ] = not( startConfig[ site ] ) 
        elif( vertexList[ v ] == -2 ):
            startConfig[ site ] = not( startConfig[ site ] )

Now we can put everything together, sweep, and see what happens!





In [ ]:



    


import numpy
import pylab

ntrack = numpy.zeros( (sweeps), int )
Ltrack = numpy.zeros( (sweeps), int )
energy = numpy.zeros( (sweeps), float )
nmax = 0
for i in range(sweeps):

    n = diagonalUpdate( startConfig, bondList, operators, n, beta )
    loopUpdate( operators, bondList, startConfig, sites )
    ntrack[ i ] = n
    nmax        = max(nmax, n)
    Ltrack[ i ] = len(operators)
    energy[ i ] = numpy.average( ntrack[0:i] )/beta
    print i, n, len(operators)



    











In [ ]:



    


fig = pylab.figure(figsize=(15,15))
axes = fig.add_subplot(1,1,1)

axes.plot(range(sweeps), ntrack, label='n over MC sweeps')
axes.plot(range(sweeps), [nmax]*sweeps, label='Maximum n over MC sweeps')
axes.plot(range(sweeps), Ltrack, label='L over MC sweeps')
axes.set_title('Expansion Order As Monte Carlo Progresses')
axes.set_xlabel('Monte Carlo Sweeps')
axes.set_ylabel('Expansion Order')
axes.legend()
pylab.show()

Content source: kshyatt/Introduction-To-QMC

Similar notebooks:

notebook.community | gallery | about