compMBP/1_classical_montecarlo_metropolis.py at main · kneureither/compMBP · GitHub

1
2
3
4
5
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
import numpy as np
import matplotlib.pyplot as plt
import time

from numba import jit


@jit(nopython=True)
def energy(system, i, j, L):
    """Energy function of spins connected to site (i, j)."""
    return -1. * system[i, j] * (system[np.mod(i - 1, L), j] + system[np.mod(i + 1, L), j] +
                                 system[i, np.mod(j - 1, L)] + system[i, np.mod(j + 1, L)])


@jit
def prepare_system(L):
    """Initialize the system."""
    system = 2 * (0.5 - np.random.randint(0, 2, size=(L, L)))
    return system


@jit(nopython=True)
def measure_energy(system):
    L = system.shape[0]
    E = 0
    for i in range(L):
        for j in range(L):
            E += energy(system, i, j, L) / 2.
    return E


@jit(nopython=True)
def metropolis_loop(system, T, N_sweeps, N_eq, N_flips):
    """ Main loop doing the Metropolis algorithm."""
    E = measure_energy(system)
    L = system.shape[0]
    E_list = []
    for step in range(N_sweeps + N_eq):
        i = np.random.randint(0, L)
        j = np.random.randint(0, L)

        dE = -2. * energy(system, i, j, L)
        if dE <= 0.:
            system[i, j] *= -1
            E += dE
        elif np.exp(-1. / T * dE) > np.random.rand():
            system[i, j] *= -1
            E += dE

        if step >= N_eq and np.mod(step, N_flips) == 0:
            # measurement
            E_list.append(E)
    return np.array(E_list),system


def magnetization():
    """ Scan through some temperatures """
    # Set parameters here
    L = 6  # Linear system size
    N_sweeps = np.arange(2000,12000, 1000)  # Number of steps for the measurements
    N_eq = 1000  # Number of equilibration steps before the measurements start
    N_flips = 10  # Number of steps between measurements
    N_bins = 10  # Number of bins use for the error analysis

    T_range = np.arange(1.5, 3.1, 0.5)

    system = prepare_system(L)
    for T in T_range:
        plt.figure()
        plt.title('Temperature: {}'.format(T))
        M_list = []
        for N_sweeps_i in N_sweeps:
            print(system)
            Es, updated_system = metropolis_loop(system, T, N_sweeps_i, N_eq, N_flips)

            M = sum(updated_system) / L ** 2
            M_list.append(M)
        print(M_list)
        plt.plot(N_sweeps, M_list, 'o')
    plt.show()

if __name__ == "__main__":
    """ Scan through some temperatures """
    # Set parameters here
    # L = np.arange(4, 10)  # Linear system size
    # N_sweeps = 5000  # Number of steps for the measurements
    # N_eq = 1000  # Number of equilibration steps before the measurements start
    # N_flips = 10  # Number of steps between measurements
    # N_bins = 10  # Number of bins use for the error analysis
    #
    # T_range = np.arange(1.5, 3.1, 0.1)
    #
    # for l in L:
    #     C_list = []
    #     E_list = []
    #     M_list = []
    #
    #     system = prepare_system(l)
    #     for T in T_range:
    #         C_list_bin = []
    #         for k in range(N_bins):
    #             Es,_ = metropolis_loop(system, T, N_sweeps, N_eq, N_flips)
    #
    #             mean_E = np.mean(Es)
    #             mean_E2 = np.mean(Es**2)
    #
    #
    #             C_list_bin.append(1. / T**2. / l**2. * (mean_E2 - mean_E**2))
    #
    #         M = sum(system) / l**2
    #         M_list.append(M)
    #         E_list.append(mean_E/(l**2))
    #         C_list.append([np.mean(C_list_bin), np.std(C_list_bin) / np.sqrt(N_bins)])
    #
    #         print(T, mean_E, C_list[-1])
    #
    #     # Plot the results
    #     C_list = np.array(C_list)
    #     plt.figure(1)
    #     plt.plot(T_range, E_list, 'o-', label='l: {}'.format(l))
    #
    #
    #     plt.figure(2)
    #     plt.errorbar(T_range, C_list[:, 0], C_list[:, 1], label='l: {},std: {}'.format(l, np.mean(C_list[:, 1])))
    #
    # Tc = 2. / np.log(1. + np.sqrt(2))
    # print(Tc)
    # plt.axvline(Tc, color='r', linestyle='--')
    # plt.xlabel('$T$')
    # plt.ylabel('$c$')
    # plt.legend()
    # plt.show()


    magnetization()