Skip to content

Ising module

animated_system

Support class for animations

Source code in bubblebox/ising.py 
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
class animated_system():
    """
    Support class for animations
    """
    def __init__(self, system = None, n_steps_per_vis = 5, interval = 1, phase_color = True):
        self.n_steps_per_vis = n_steps_per_vis
        self.system = system
        figsize = (6,6)



        plt.rcParams["figure.figsize"] = figsize

        self.phase_color = phase_color


        self.fig, self.ax = plt.subplots()    

        cc = np.random.uniform(0,1,(3, 6))
        cc[:,0] = np.array([0,0,0])

        cc = np.array([[0.        , 0.        , 0.        ],
                    [0.8834592 , 0.36962255, 0.21858202],
                    [0.64961546, 0.79727038, 0.55362479],
                    [0.22449319, 0.56457326, 0.60815318],
                    [0.75835695, 0.729311  , 0.54213821]]).T





        self.color = interp1d(np.linspace(-1,1, 5), cc)



        self.ani = FuncAnimation(self.fig, self.update, interval=interval, 
                                          init_func=self.setup_plot)#, blit=True,cache_frame_data=True)


    def update(self, j):
        for i in range(self.n_steps_per_vis):
            #self.system.lattice = np.random.randint(0,2,self.system.lattice.shape) #advance()
            self.system.advance()
        self.bubbles.set_color(self.color(self.system.Z).T)
        if self.phase_color:
            self.bubbles.set_color(self.color(self.system.identical_neighbors().ravel()).T)

        #self.infotext.set_text("%.2f" % self.system.energy())

        return self.bubbles,





    def setup_plot(self):


        #if len(self.system.lattice.shape)==2:
        x,y = np.array(np.meshgrid(np.arange(self.system.N), np.arange(self.system.N))).reshape(2, -1)

        s = 50000/self.system.N2
        self.bubbles = self.ax.scatter(x, y, c=self.color(self.system.Z.ravel()).T, s=s, marker = "8")
        self.ax.axis("off")
        self.infotext = self.ax.text(1,-2,"test", ha = "left", va = "center")
        plt.xlim(-1, self.system.N+1)
        plt.xlim(-1, self.system.N+1)

        return self.bubbles,

isingbox

Source code in bubblebox/ising.py 
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 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
116
117
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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
class isingbox():
    def __init__(self, N, kT, flipratio = 0.1):
        """
        A 2D (N by N) Ising model

        Arguments
        ---
        N          Lattice dimensions
        kT         Boltsmann factor times the temperature 
        flipratio  The ratio of spins to flip at each iteration

        Methods
        ---
        reset(kT)          reset the lattice to all spins aligned, and set kT
        energy_per_site()  compute the energy per site


        Example usage
        ---

        I = ising(10, kT = 1.0) # an Ising lattice with 10x10 spins, kT = 1.0

        I.advance() # do one Monte Carlo integration

        I.evolve(10) # do 10 Monte Carlo integrations

        energy = I.compute_energy() # compute energy

        mangetization = I.compute_magnetization() # compute magnetization

        energy, magnetization = I.sample(100, 1000) # collect 100 samples, separated by 1000 MC-iterations




        """

        self.kT = kT
        self.N = int(N)
        self.N2 = int(N**2)
        self.Z = np.ones(self.N**2, dtype = int)

        #perhaps not start out in this state?
        #self.Z[np.random.choice(self.N2, int(self.N2/2), replace = False)] = -1

        self.Z_ind = np.arange(self.N2)
        self.n = int(flipratio*self.N2)
        self.flipratio = flipratio

        self.env = self.energy_env_torus() #_torus(arange(self.N2))

        self.n_logdist = 1000000
        self.logdist = self.kT*np.log(np.random.uniform(0,1,self.n_logdist))/2.0


        self.t = np.zeros(10, dtype = float)

    def reset(self, kT):
        """
        Reset the system to a non-magnetized state ( |M| = 0 )

        and set kT to the specified value

        """
        self.kT = kT
        self.Z = np.ones(int(self.N**2), dtype = int)
        #self.Z[np.random.choice(int(self.N2), int(self.N2/2), replace = False)] = -1
        self.Z_ind = np.arange(self.N2)
        self.n = int(self.flipratio*self.N2)
        self.env = self.energy_env_torus()

    def energy_per_site(self,Z):
        """Compute sitewise energy (returns a lattice)"""
        return (Z[2:,1:-1]+Z[:-2,1:-1] +Z[1:-1,2:]+Z[1:-1,:-2]) #[1:-1,1:-1]



    def log_uniform(self, n):
        """A logarithmic random distribution"""
        return self.logdist[np.random.randint(0,self.n_logdist,n)]

    def energy_env_torus(self):
        """Compute energy per site for torus PBC"""
        z = self.Z.reshape((self.N,self.N))
        return (np.roll(z, 1,0) + np.roll(z,-1,0) + np.roll(z,1,1) + np.roll(z,-1,1)).ravel()

    def update_env(self, signs, z_ind):
        """update energy environment with a change in Z[z_ind]"""

        self.env[z_ind-1] += signs
        self.env[(z_ind+1)%self.N2] += signs
        self.env[z_ind-self.N] += signs
        self.env[(z_ind+self.N)%self.N2] += signs



    def advance(self):
        """Do one Monte Carlo step"""

        sel = np.random.choice(self.N2, self.n, replace = False) #bottleneck


        dz = self.Z[sel]
        denv = self.env[sel]

        dE = dz*denv

        flips = -dE>self.log_uniform(self.n)
        dz[flips] *= -1            
        self.Z[sel] = dz

        self.update_env(2*dz[flips], sel[flips]) #bottleneck



    def evolve(self, n_steps):
        """Do n_steps Monte Carlo steps"""

        for i in range(n_steps):
            self.advance()

    def compute_magnetization(self):
        """Compute magnetization"""
        return np.sum(self.Z)

    def compute_energy(self):
        """Compute total energy""" 
        return -np.sum(self.Z*self.env)/2.0


    def sample(self, n_samples, n_separations):
        """
        Collect n_samples separated by n_separations Monte Carlo steps

        Returns 
        energies_per_site, specific heat, absolute_magnetization_per_site
        """

        energy        = np.zeros(n_samples)
        magnetization = np.zeros(n_samples)

        energy[0] = self.compute_energy()
        magnetization[0] = self.compute_magnetization()
        for i in range(1, n_samples):
            self.evolve(n_separations)

            energy[i] = self.compute_energy()
            magnetization[i] = self.compute_magnetization()

        return np.mean(energy)/self.N2, np.var(energy)/(self.N2*self.kT**2), np.mean(magnetization)/self.N2

    def run(self, phase_color = False, n_steps_per_vis = 100):
        """Run simulation interactively"""
        self.run_system = animated_system(system = self, n_steps_per_vis=n_steps_per_vis, interval = 1, phase_color = phase_color)
        plt.show()

__init__(N, kT, flipratio=0.1)

A 2D (N by N) Ising model

Arguments

N Lattice dimensions kT Boltsmann factor times the temperature flipratio The ratio of spins to flip at each iteration

Methods

reset(kT) reset the lattice to all spins aligned, and set kT energy_per_site() compute the energy per site

Example usage

I = ising(10, kT = 1.0) # an Ising lattice with 10x10 spins, kT = 1.0

I.advance() # do one Monte Carlo integration

I.evolve(10) # do 10 Monte Carlo integrations

energy = I.compute_energy() # compute energy

mangetization = I.compute_magnetization() # compute magnetization

energy, magnetization = I.sample(100, 1000) # collect 100 samples, separated by 1000 MC-iterations

Source code in bubblebox/ising.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
def __init__(self, N, kT, flipratio = 0.1):
    """
    A 2D (N by N) Ising model

    Arguments
    ---
    N          Lattice dimensions
    kT         Boltsmann factor times the temperature 
    flipratio  The ratio of spins to flip at each iteration

    Methods
    ---
    reset(kT)          reset the lattice to all spins aligned, and set kT
    energy_per_site()  compute the energy per site


    Example usage
    ---

    I = ising(10, kT = 1.0) # an Ising lattice with 10x10 spins, kT = 1.0

    I.advance() # do one Monte Carlo integration

    I.evolve(10) # do 10 Monte Carlo integrations

    energy = I.compute_energy() # compute energy

    mangetization = I.compute_magnetization() # compute magnetization

    energy, magnetization = I.sample(100, 1000) # collect 100 samples, separated by 1000 MC-iterations




    """

    self.kT = kT
    self.N = int(N)
    self.N2 = int(N**2)
    self.Z = np.ones(self.N**2, dtype = int)

    #perhaps not start out in this state?
    #self.Z[np.random.choice(self.N2, int(self.N2/2), replace = False)] = -1

    self.Z_ind = np.arange(self.N2)
    self.n = int(flipratio*self.N2)
    self.flipratio = flipratio

    self.env = self.energy_env_torus() #_torus(arange(self.N2))

    self.n_logdist = 1000000
    self.logdist = self.kT*np.log(np.random.uniform(0,1,self.n_logdist))/2.0


    self.t = np.zeros(10, dtype = float)

advance()

Do one Monte Carlo step

Source code in bubblebox/ising.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def advance(self):
    """Do one Monte Carlo step"""

    sel = np.random.choice(self.N2, self.n, replace = False) #bottleneck


    dz = self.Z[sel]
    denv = self.env[sel]

    dE = dz*denv

    flips = -dE>self.log_uniform(self.n)
    dz[flips] *= -1            
    self.Z[sel] = dz

    self.update_env(2*dz[flips], sel[flips]) #bottleneck

compute_energy()

Compute total energy

Source code in bubblebox/ising.py 
131
132
133
def compute_energy(self):
    """Compute total energy""" 
    return -np.sum(self.Z*self.env)/2.0

compute_magnetization()

Compute magnetization

Source code in bubblebox/ising.py 
127
128
129
def compute_magnetization(self):
    """Compute magnetization"""
    return np.sum(self.Z)

energy_env_torus()

Compute energy per site for torus PBC

Source code in bubblebox/ising.py 
87
88
89
90
def energy_env_torus(self):
    """Compute energy per site for torus PBC"""
    z = self.Z.reshape((self.N,self.N))
    return (np.roll(z, 1,0) + np.roll(z,-1,0) + np.roll(z,1,1) + np.roll(z,-1,1)).ravel()

energy_per_site(Z)

Compute sitewise energy (returns a lattice)

Source code in bubblebox/ising.py 
77
78
79
def energy_per_site(self,Z):
    """Compute sitewise energy (returns a lattice)"""
    return (Z[2:,1:-1]+Z[:-2,1:-1] +Z[1:-1,2:]+Z[1:-1,:-2]) #[1:-1,1:-1]

evolve(n_steps)

Do n_steps Monte Carlo steps

Source code in bubblebox/ising.py 
121
122
123
124
125
def evolve(self, n_steps):
    """Do n_steps Monte Carlo steps"""

    for i in range(n_steps):
        self.advance()

log_uniform(n)

A logarithmic random distribution

Source code in bubblebox/ising.py 
83
84
85
def log_uniform(self, n):
    """A logarithmic random distribution"""
    return self.logdist[np.random.randint(0,self.n_logdist,n)]

reset(kT)

Reset the system to a non-magnetized state ( |M| = 0 )

and set kT to the specified value

Source code in bubblebox/ising.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
def reset(self, kT):
    """
    Reset the system to a non-magnetized state ( |M| = 0 )

    and set kT to the specified value

    """
    self.kT = kT
    self.Z = np.ones(int(self.N**2), dtype = int)
    #self.Z[np.random.choice(int(self.N2), int(self.N2/2), replace = False)] = -1
    self.Z_ind = np.arange(self.N2)
    self.n = int(self.flipratio*self.N2)
    self.env = self.energy_env_torus()

run(phase_color=False, n_steps_per_vis=100)

Run simulation interactively

Source code in bubblebox/ising.py 
157
158
159
160
def run(self, phase_color = False, n_steps_per_vis = 100):
    """Run simulation interactively"""
    self.run_system = animated_system(system = self, n_steps_per_vis=n_steps_per_vis, interval = 1, phase_color = phase_color)
    plt.show()

sample(n_samples, n_separations)

Collect n_samples separated by n_separations Monte Carlo steps

Returns energies_per_site, specific heat, absolute_magnetization_per_site

Source code in bubblebox/ising.py 
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
def sample(self, n_samples, n_separations):
    """
    Collect n_samples separated by n_separations Monte Carlo steps

    Returns 
    energies_per_site, specific heat, absolute_magnetization_per_site
    """

    energy        = np.zeros(n_samples)
    magnetization = np.zeros(n_samples)

    energy[0] = self.compute_energy()
    magnetization[0] = self.compute_magnetization()
    for i in range(1, n_samples):
        self.evolve(n_separations)

        energy[i] = self.compute_energy()
        magnetization[i] = self.compute_magnetization()

    return np.mean(energy)/self.N2, np.var(energy)/(self.N2*self.kT**2), np.mean(magnetization)/self.N2

update_env(signs, z_ind)

update energy environment with a change in Z[z_ind]

Source code in bubblebox/ising.py 
92
93
94
95
96
97
98
def update_env(self, signs, z_ind):
    """update energy environment with a change in Z[z_ind]"""

    self.env[z_ind-1] += signs
    self.env[(z_ind+1)%self.N2] += signs
    self.env[z_ind-self.N] += signs
    self.env[(z_ind+self.N)%self.N2] += signs