/
02-wolff.c
119 lines (104 loc) · 3.38 KB
/
02-wolff.c
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
/*****************************************************************************
*
* Tutorial: How to ALPSize your applications
*
* Copyright (C) 2005-2010 by Synge Todo <wistaria@comp-phys.org>
*
* Permission is hereby granted, free of charge, to any person or organization
* obtaining a copy of the software and accompanying documentation covered by
* this license (the "Software") to use, reproduce, display, distribute,
* execute, and transmit the Software, and to prepare derivative works of the
* Software, and to permit third-parties to whom the Software is furnished to
* do so, all subject to the following:
*
* The copyright notices in the Software and this entire statement, including
* the above license grant, this restriction and the following disclaimer,
* must be included in all copies of the Software, in whole or in part, and
* all derivative works of the Software, unless such copies or derivative
* works are solely in the form of machine-executable object code generated by
* a source language processor.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
* SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
* FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
*
*****************************************************************************/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define L 32
#define N (L*L)
#define T 2.2
#define MCSTEP (1 << 15)
#define MCTHRM (MCSTEP >> 3)
#define SEED 93812
double random_01() {
return (double)(rand()) / RAND_MAX;
}
int main() {
int x, y, k, s, mcs;
/* setting up square lattice */
int nn[N][4];
for (y = 0; y < L; ++y)
for (x = 0; x < L; ++x) {
nn[x+L*y][0] = ((x+L-1)%L) + L*y;
nn[x+L*y][1] = ((x+1)%L) + L*y;
nn[x+L*y][2] = x + L*((y+L-1)%L);
nn[x+L*y][3] = x + L*((y+1)%L);
}
/* random number generator */
srand(SEED);
/* spin configuration */
int spin[N];
for (s = 0; s < N; ++s) spin[s] = 1;
int sz = N;
/* stack for uninspected sites */
int stck[N];
int is = 0;
/* connecting probability */
double pc = 1 - exp(-2./T);
/* measurement */
double m = 0;
double m2 = 0;
double m4 = 0;
/* timer */
clock_t tm = clock();
for (mcs = 0; mcs < MCSTEP + MCTHRM; ++mcs) {
s = random_01() * N;
int so = spin[s];
spin[s] = -so;
stck[0] = s;
is = 1;
int cs = 0;
while (is) {
++cs;
int sc = stck[--is];
for (k = 0; k < 4; ++k) {
int sn = nn[sc][k];
if (spin[sn] == so && random_01() < pc) {
stck[is++] = sn;
spin[sn] = -so;
}
}
}
sz -= 2 * so * cs;
if (mcs >= MCTHRM) {
double dsz = sz / (double)N;
m += dsz;
m2 += dsz * dsz;
m4 += dsz * dsz * dsz * dsz;
}
}
/* output results */
printf("Magnetization = %f\n", m / MCSTEP);
printf("Magnetization^2 = %f\n", m2 / MCSTEP);
printf("Magnetization^4 = %f\n", m4 / MCSTEP);
printf("Binder Ratio of Magnetization = %f\n", m2 * m2 / m4 / MCSTEP);
fprintf(stderr, "Elapsed time = %f sec\n", difftime(clock(), tm) / CLOCKS_PER_SEC);
return 0;
}