Ising2D/ising.cpp at main · vataspro/Ising2D · 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
#include "ising.h"


Ising2D::Ising2D(int size, double beta, double external_field) :
                 lattice(size, std::vector<int>(size, 0)),
                 N(size),
                 b(beta),
                 H(external_field) {

    initialiseLattice();
    for (int k=0; k<THERMAL_ITERS; k++) {
        sweep();
    }
}


void Ising2D::initialiseLattice() {
        for (int i=0; i<N; i++) {
        ranlxs_fill_ising_vector(lattice[i]);
    }

}


void Ising2D::sweep() {

    for (int k=0; k<SWEEPS_PER_TRAJ; k++) {

        for (int i=0; i<N; i++) {
            for (int j=0; j<N; j++) {
                double energy = lattice[i][j] *
                                    ( lattice[i][(j+1)%N]     +
                                        lattice[i][(N+j-1)%N] +
                                        lattice[(i+1)%N][j]   +
                                        lattice[(N+i-1)%N][j] );

                energy += H * lattice[i][j];

                if ((energy <= 0) || (ranlxs_rand() < exp(-2. * b  * energy)))
                    lattice[i][j] *= -1;
            }
        }
    }
}


double Ising2D::getMagnetisation() const {

    double s = 0;
    for (int i=0; i<N; i++) {
        for (int j=0; j<N; j++) {
            s += lattice[i][j];
        }
    }

    return s / (N*N);
}


double Ising2D::getMagnetisationSquared() const {

    double s = 0;
    for (int i=0; i<N; i++) {
        for (int j=0; j<N; j++) {
            s += lattice[i][j]*lattice[i][j];
        }
    }

    return s / (N*N);
}


int Ising2D::computeEnergy() const {
    int energy = 0;
    int Lx = lattice.size();          // number of rows
    int Ly = lattice[0].size();       // number of columns

    for (int i = 0; i < Lx; ++i) {
        for (int j = 0; j < Ly; ++j) {
            int s = lattice[i][j];
            // Periodic boundary conditions
            int right = lattice[i][(j + 1) % Ly];
            int down = lattice[(i + 1) % Lx][j];

            energy -= s * (right + down);
        }
    }

    return energy;
}


void Ising2D::print() const {

    for (int i=0; i<N; i++) {
        for (int j=0; j<N; j++) {
            std::cout << std::setw(2) << std::setfill(' ') <<  lattice[i][j] << " ";
        }
        std::cout << "\n";
    }
}


void Ising2D::save(std::string filename) const {

    std::ofstream out(filename);
    if (!out) {
        std::cerr << "error: failed to open file.\n";
        exit(1);
    }

    for (const auto row : lattice) {
        for (const auto elem : row) {
            out << elem << " ";
        }
        out << "\n";
    }

    out.close();
}