THE TWO-DIMENSIONAL ISING MODEL OF FERROMAGNETISM USING THE METROPOLIS ALGORITHM
Abstract:
The properties of a two-dimensional lattice of spins were investigated using the Ising Model. A Monte Carlo method based on the Metropolis algorithm was implemented. This algorithm generates a sequence of samples whose statistical distribution approaches the Boltzmann distribution of the magnetic system as more iterations are executed. Under no applied field, the magneti sation, the internal energy and the heat capac ity were found to vary significantly within the vicinity of a critical temperature, at which a continuous phase transition occurs. In particular, plotting the magnetisation against temperature allowed to estimate the critical temperature asTc = 2.272 0.004 and the critical exponent as β = 0 .109 0.018. The theoretic al values predicted by Onsager, Tc = 2 /(1 + ln 2) and β = 1 /8, are both found within erro r range. When a non-zero field was applied, the system was found to have a non-zero magnetisation above the critical temperature and a single equilibrium state below it.
±
1
±
√
INTRODUCTION
low Tc ) according to the following equation:
The magnetic properties of ferromagnets can be investigated using a simplified model of spin-spin interactions known as the Ising model. In particular, the two-dimensional Ising model consists of a square lattice in which each site has either spin up (+1) or down (-1). Within such mode l the energy of a N x N lattice can be written as: N2
E=
−J ∑ si s j − µH ∑ si
(1)
i=1
Spins on adjacent sites, called nearest neighbours, in-
M=
− 1
sinh(ln(1 +
√
Tc 2) ) T
1/ 8
− 4
(2)
Within the vicinity of the critical temperat ure, this expression can be Taylor expanded, in which case the magnetization is found to vary as M ∝ (Tc T )β , where β = 1/8. An analysis of the computational aspects of the problem is disclosed in the next section . Section 3 includes a description of the implementation of the necessary algorithm to solve the problem. In subsection 3.1, the performance of the program is discussed, including steps taken to improve its efficiency. Section 4 presents the quantitative and qualitative results
−
teract a pairwise manner with strength J, known as the in exchange energy coefficient . In addition, the externally applied magnetic field H couples linearly with the magnetisation. µ is the magnetic moment of each spin.
of the problem a discussion of their nificance. Theand overall conclusions are physical present edsigin section 5.
In 1944, Onsager published the analytic solution to the two-dimensional Ising model [1]. The critical temperature in the thermo dynamic limit (i.e. an infinitely large lattice, N ∞) is predicted to be Tc = 2 /(1 + ln 2). Below this temperature, the material behaves like a ferromagnet, as the spins tend to align in the most stable configuration even in the absence of a field. Above it, the spins tend to be randomly oriented unless the field is non-zero, as expected for a paramagnet. Predicted by Onsager and rigourously proved by Yang in 1952 [2], the magnetisation is expected to depend on the temperature (be-
2
√
→
ANALYSIS
The statistical mechanics of complex physical systems usually poses problems that are very difficult to solve by analyt ical methods. Numerical simulation techniques are therefore indispensable tools to the study of such problems, which range from phase transitions in binary alloys to neural networks. The Ising model of ferromagnetism exhibits this feature clearly. For example, if we considered a small 10 x 10 lattice of spins, the total number of microstates would be of 2 100 1.2 x 10 30 . Summing
3
over such a large number of microstates to perform the ensemble averages involved in the calculations of the thermodynamic variables of the system would be very inefficient. Instead of choosing configurations randomly and then weighting them with the appropriate Boltzmann distribution, a more efficient approach to this problem involves choosing the spin configurations with a probability ∝ e−E /kB T (where E is the energy of the configuration and T is the temperature set by the reservoir) and then weighti ng them evenly. This is the underlying idea behind the Metropolis algorit hm [3]. The Metropolis algorithm consists of the following steps:
IMPLEMENTATION
The program listing can be found in the Appendix. A class, Matrix, which can be found in file Matrix.h, was created to represent two-dimensional arrays in a simple way. There is only one function de-
The N 2 operations described above may be thought of as a single time step. In order to reach a quasi equilibrium state, the simulation must be run for several time steps. At this point the thermodynamic and magnetic properties of the system may be measured by repeating several time steps and then averaging over the relevant quantities stored at each time instant. This time average is equivalent to an ensemble aver-
fined withincol). the class, which allows to call a matrix as lattice(row, random_generator generates a random number between 0 and 1 based on gsl_rng_uniform. This function is called by initialize_lattice_random to generate a random set of spins. Alternatively, the lattice may be initialised with all spins up using the function initialize_lattice_uniform. The function energy returns the energy of the spin configuration using equation 1. This function computes the sum over the nearest neighbours of all lattice sites by calling the function nearest_neighbours, which implements periodic boun dary conditions. J and µ are defined in this function and are both set to unity for convenience. The Boltzmann constant, kB , is also defined as a constant and set to unity. The function magnetisation returns the magnetization per unit spin. Given that µ is set to 1, this function simply sums all spins and divides by N 2 . Metropolis executes steps 3-4 described in section 2. The row and column of the lattice site considered are passed as parameter s. This function creates a new lattice that differs from the previous one by a flipped spin. The flipping is performed by the function flip_spin. The previous lattice configuration is saved so that the energy difference between the two can be calculated. To perform steps 4.(b ) and 4.(c) the Boltzmann factor is compared to a random number generated by random_generator. Step 5 can be performed in two disti nct ways : all lattice sites may be visited systematically line af-
age. In order to ensure that each lattice site had four nearest neighbours, periodic boundary conditions were used. Hence, when the selecte d site was along the border of the lattice, its neighbours were taken from the opposite side of the lattice. This effectively reduced the geometry to a torus, thus avoiding undesirable effects at the edges. Even though the two-dimensional Ising model of ferromagnetism was solved analytically by Onsager, this numerical simulation technique is in fact of great value, as it may be generalized to higher dimensions, for which no analytical solution has been found, and to more complex interactions (e.g. next-nearestneighbour interactions). Nevertheless, the quantitative results obtained through this method were compared to the analytical results derived by Onsager.
ter2 line ( single_Monte_Carlo_step_regular) or N lattice sites may be picked randomly ( single_ Monte_Carlo_step_random). In the latter case, the row and column of the site considered correspond to random_ the integral part of N times the output of generator. In the main, the size of the lattice, the initial temperature and the applied magnetic field are initialized. The program has two modes: mode_1 simply takes M single time steps, saves the magnetization per spin at each instant and prints the final lattice configuration using Gnuplot. mode_2 covers a range of temperatures and saves the following physical variables at equilibrium for each temperature: magnetization per site, energy per site, their standard deviations and heat capacity. The former mode is useful to underst and what the equilibrium states at different temperatures
1. Create a N x N lattice 2. Select a lattice site 3. Determine the energy change ∆E involved in flipping the spin 4. Decide whether to flip the spin or leave it unchanged as follows: (a) If ∆E < 0 flip the spin (b) If ∆E > 0 and exp( ∆E /kB T ) > r , where r is a random number generated in the interval [0,1], flip the spin
−
(c) Otherwise, leave the spin unchanged 5. Repeat steps 2-4 for N 2 lattice sites
2
are, whereas the latter is useful to plot the relevant physical quantities against temperature. The heat capacity is estimated using two methods: calculating the discrete derivative C = dE /dT and using the Fluctuation-Dissipation Theorem, according to which:
σE2 C = kB T 2
val allowed to significantly reduce the run time, especially for large N. M was also the number of time steps taken to compute the ensemble average of the relevant quantities. This also justifies the need to have a greater M closer to the critical temperature, as the fluctuations become significant in this region. The CPU time required to perform 1000 single time steps (STS) at a fixed temperature for a 20 x 20 lattice was 30 s, whereas for a 30 x 30 lattice the same procedure took 2 min 40 s. Performing 100 STS for T [1.5, 2.0] [3.0, 3.5], with δT = 0.1, and 1000 STS for T [2.0, 3.0], with δ T = 0.05, for a 10 x 10 lattice took 1 min 30 s. The same result for N = 20 and N = 30 required 24 min and 1h 58 min of CPU time, respectively.
(3)
where σ E is the standard deviation of energy. In the former case, the energy and the temperature of the previous iteration must be saved in order to calculate the appropriate ratio of differences. Close to the critical temperature, the fluctuations are so large and the energy gap required to move from one equilibrium state (all spins down) to another (all spins up) is so low that the average magnetization per site tends to switch from 1 to -1 and vice-versa. This effect tends to artificially decrease the value of the magnetization just below the critical temperature. To avoid this, the absolute value of the magnetization was printed in mode 2.
∈
4 4.1
3.1
∈
RESULTS AND DISCUSSION No Applied Field (H = 0)
Performance By visualizing the lattice configuration and the magnetisation against time, a continuous phase transition was observed to occur near T = 2 .30 [ J /kB ]. Above this critical temperature Tc , the equilibrium state was a randomly aligned lattice (figure 1) with average magnetization zero (figure 2). Below this temperature, two different equilibrium states were observed: one with all spins down and another with all spins up. The system would attain one of the two states depending on the initial configuration (figure 3).
The two initial lattice configurations were used at different temperatures. The uniform lattic e configuration reached quasi-equilibrium more rapidly than the random one below the critical temperature, whereas above it the random configuration was only slightly faster to reach near-equilibrium. Given that the initial temperature in mode 2 was always set below the critical temperature, the lattice was always initialized with all spins up. For the following temperatures the initial configuration was the equilibrium configuration of the previous temperature, therefore the initialization of the lattice was only relevant for the starting temperature. In addition, there were two alternative ways of picking the spin to flip . The two functions, single_Monte_Carlo_step_regularand single_ Monte_Carlo_step_random, were tested for different temperatures, but no significant differences in efficiency were note d. As a result, the spins were flipped sequentially using single_Monte_Carlo_ step_regular. Moreover, since the plots against temperature have their most distinctive features within the vicinity of the critical temperature, the temperature increment was set to 0.05 for the interval T [2.0, 3.0] and 0.1 elsewhere. The system also requires more time to attain near-equilibrium close to the critical temperature, so the number of single time steps M taken to reach equilibrium was of 1000 for T [2.0, 3.0] and 100 elsewhere. Having 1000 time steps only in this inter-
∈
Figure 1: Spin configuration of a 20 x 20 latti ce at temperature T = 4 .0 [ J /kB ] after 100 single time steps . Black squares represent spin ’down’ whereas white squares cor-
∈
respondabove to spinTc’up’. As expected, the spins are random ly aligned .
3
sation vs. temperature graph approaches zero magnetization more smoothly than the theoretical prediction (figure 4), so the critical temperature could only be estimated as Tc = 2 .3 0.1[J /kB ]. Onsager’s critical temperature is within this interval.
Magnetization vs. Time (N = 20, T = 4.0) 1
±
] [µ 0.5 et si re p 0 n o tia itz e n g a-0.5 M
-1 0
1
Magnetization vs. Temperature (N = 30) Data Onsager Analytic Solution
0.9
10
20
30
40
50
60
70
80
90
] 0.8 µ [ e its 0.7 er 0.6 p n 0.5 o it az it 0.4 e n 0.3 ag M0.2
100
Number of Time Steps
Figure 2: Magnetization versus time for a 20 x 20 lattice at temperature T = 4.0 [J /kB ]. Since the temperature is above Tc , the average magnetization near equilibrium is zero.
0.1 0 1 .4 1.6 1.8
2.2
2.4 2.6 2.8
3
Temperature [J/kB]
3.2
3.4 3.6
Figure 4: Magnetisation versus tempe rature for a 30 x 30 lattice (red) and Onsager’s analytic solution (green, equation 2). As expected, the magnetisati on decreases abruptl y close to Tc . However, the numerical estimate diverges from the theoretical prediction just above the critical temperature, as the magnetisa tion goes more smooth ly to zero. This is due to the finite size of the system.
Magnetization vs. Time (N = 20, T = 1.5) 1.5
2
Uniform Random
1
] [µ e ist 0.5 re p n o it 0 az tie n-0.5 g a M
Energy vs. Temperature (N = 30) -0.6
-1
-0.8 -1.5 0
10
20
30
40
50
60
70
80
90
100 ] -1 J [ e it s-1.2 r e p y-1.4 g r e n E -1.6
Number of Time Steps
Figure 3: Magnetization time for a 20 x 20 lattice at temperature T = 1 .5 [ J /kversus B ] for two different initial configurations. The two initial configurations gave rise to two different equilibrium states, one with all spins up and another with all spins down. Due to the absence of an externally applied field, the system is invariant under the flip of all spins, and therefore these two equilibrium states are energetically equivalent and equally likely to be observed.
-1.8 -2 1 .4 1.6 1.8
By plotting the magnetisation and the energy against temperature a phase transition could be identified, as the magnetisation decreases sharply to zero close to the critical temperature and the energy per site increases in the same region. The same feature s were observed for different values of N. However, none of the two graphs allowed to accurately estimate the critical temperature by visual inspection. Since the lattices considered were of finite size, the magneti-
2
2.2
2.4 2.6 2.8
3
Temperature [J/kB]
3.2
3.4 3.6
Figure 5: Energy per site versus temperature for a 30 x 30 lattice. As expected, the energy increase s across the phase transition. The variation of the energy is not as steep as for the magnetisation, therefore a better estimate of the critical temperature is found from the M vs. T graph.
The heat capacity C was measured via two different methods. The graph C vs. T obtained by calculating the discrete derivative dE/dT was far less
4
clear than the one obtained from the FluctuationDissipation Theorem due to the random errors in the simulation (figure 6). Nevertheless, the best graph does allow to make a better estimate of Tc , since the heat capacity peaks at this temper ature: for N = 30, Tc = 2 .3 0.05 (figure 7). The heat capacity is expected to diverge according to Onsager’s model, but such divergence is not observed in lattices of finite size.
Onsager’s prediction that the magnetisation varies as M ∝ (Tc T )β was used to make an accurate estimate of both the critical temperature Tc and the universal exponent β for different values of N (10, 15, 20, 25, 30). The function f (T ) = C1 (Tc T )β was fit to the data points from the numerical simulation using the Gnuplot command fit . However, given that this expression is only valid for temperatures just below the critical temperature, only some data points were selected. Figure 8 clearly show s that the data points at the top section of the graph with negative curvature do follow the analyti c solution very closel y. Those data points were the ones considered when fitting to the model. In addition, the uncertainties in each of the measurements were taken to be the standard deviation of fluctuations in the magnetisation.
−
−
±
2500
Heat Capacity vs. Temperature (N = 30)
2000 ] B k [
1500 ty i c a p a C1000 t a e H
Magnetization vs. Temperature with Fitting (N = 30) 1.2
Complete Data Fitting Data f(x)
500
] 1 [µ 0 1 .4 1.6 1.8
2
2.2 2.4 2.6 2.8
3
Temperature [J/kB]
site 0.8 er p n o 0.6 tia z tie n 0.4 g a M
3.2 3.4 3.6
Figure 6: Heat Capacity versus Temperature for a 30 x 30 lattice using the Discrete Deriviative dE/dT. Although the general shape of the curve is noticeable, the noise due to the random errors in the simulation is significant.
0.2 0 1 .4
1.6
1.8
2
2.2
2.4
2.6
2.8
3
3.2
Temperature [J/kB]
3.4
3.6
Heat Capacity vs. Temperature (N = 30) 2200
Figure 8: Magnetisation versus tempe rature for a 30 x 30 lattice. The complete set of data points, the selected data (and respective error bars) and the fitting model are presented in red, blue and magenta, respectively. For N = 30, Tc = 2.277 0.018[J /kB ] and β = 0 .09 0.02.
2000 1800 ]1600 B [k
1400
y itc1200 a p a1000 C t a 800 e H 600
±
Table 1: Summary of the estimated values of the critical temperature Tc and universal coefficient β as parameters of the model
400 200 0 1 .4 1.6 1.8
±
Table 1 summarizes the estimates of Tc and β for the different values of N considered:
2
2.2 2.4 2.6 2.8
3
Temperature [J/kB]
3.2 3.4 3.6
N 10 15 20 25 30
Figure 7: Heat Capacity versus Temperature (red) for a 30 x 30 lattice using the Fluctuation-Dissipation Theorem. An spline (blue) was added to visualize the peaked curve more clearly and hence to estimate the critical temperature more accurately. In addition, a vertical line (green) was plotted to identify the data point of maximum heat capacity. From
Tc [J /kB ] 2.35 0.03 2.30 0.02 2.299 0.019 2.286 0.004 2.277 0.018
± ± ± ± ±
β 0 .7 0.40 0 .17 0 .11 0 .09
± 0 .3 ± 0.15 ± 0.06 ± 0.06 ± 0.02
Onsager’s prediction of the critical temperature, Tc = 2/(1 + ln 2), only applies in the limit N ∞.
Tc = 2.3 this graph, 0.05, which is vs. more accurate than the previous estimate made from the M T graph.
±
√
5
→
To extrapolate the numerical predictions to N ∞, ln(Tc ) was plotted against ln (N ): the slope of the linear regression was found to be close to -2. Hence, plotting Tc against 1/N 2 should have an intersect with the vertical axis corresponding to Tc as derived by Onsager (figure 9). The critical temperature was estimated as Tc = 2.272 0.004 [ J /kB ].
→
ity class. A universality class is characte rized by the dimensionality of space (in this case, d = 2, since the lattice is two-dimensional) and the number of components of the order parameter (in this case, n = 1, since the spin is a scalar). Regardless of the size N, all spin lattices belong to this d = 2, n = 1 universality class, and therefore the predicted value of β from the numerical simulati on should correspond to the theoretical value of 1/8 for all N.
± Critical Temperature versus 1/N
2.38
2
However, the two estimates of β for N = 10, 15 are clearly far from the expected value. Figure 10 shows M vs. T (with the fitting model) for N = 10. Clearly, selecting the data points that should be used to fit the model is not as clear as for N = 30, for example. Since that selection was greatly influenced by the previous knowledge about the critical temperature (namely, for example, from C vs. T) and given the large correlation of 0.445 between Tc and β, the two estimates of β for N = 10, 15 were neglected. The average of the three values for N = 20, 25, 30 gives a final estimate of β = 0 .109 0.018. The theo retical value is within error range.
2.36 2.34 2.32 c
T
2.3 2.28 2.26 2.24 0
0.002
0.004
0.006
0.008
0.01
1/N2
±
Figure 9: Tc versus 1/N 2 and respective linear regression. The intersect was estimated as Tc = 2.272 0.004 [J /kB ]
±
4.2 Magnetization versus Temperature with Fitting (N = 10)
Complete Data Model Fitting Data
1
] [µ te si 0.5 er p n o tia 0 z
When a non-zero magnetic field was applied, a nonzero magnetisation was found above the critical temperature. The magnetisation has the same sign as the applied field, that is, if the sign of the field is switched, the sign of the magnetisation changes as well. In addition, the magnetisation tends to increase with the applied field. Figure 11 illustrates these observations.
eti n g a -0.5 M
-1 1 .4 1.6 1.8
Non-Zero Applied Field (H = 0)
2
2.2 2.4 2.6 2.8
3
Temperature [J/kB]
3.2 3.4 3.6
Below the critical temperature, there is now a single stable equilibrium state. By considering the same initial configurations as in figure 4, for an applied field H = 1.0 [J /µ] both lattices reach the same equilibrium state, contrary to previous observations. The application of the magnetic field explicitly breaks the symmetry of the system under the switch of all spins. Hence, there is now a preferred direction, so one of the uniform states is raised in energy while the other is lowered.
Figure 10: Magnetisation versus temperature for a 10 x 10 lattice. The complete set of data points, the selected data and respective error bars and the fitting model are presented in red, blue and magenta, respecti vely. Clearly the graph is not as clear as for N = 30 (figure 8), for example, which justifies the difficulty in selecting the appropriate points to fit the model.
As the name suggests, the universal exponent β should be the same for all states of the same universal-
6
Magnetization versus Time: Varying Field
1.5
H H H H
= = = =
1.0 0.5 0.1 0.5
REFERENCES
1
] [µ tei 0.5 s re p n 0 io atz it e-0.5 n g a M
[1] L. Onsager,
dimensional model with an order-disorder transition,” Phys. Rev., vol. 65, pp. 117–149, Feb 1944. [2] C. N. Yang, “The spontaneous magnetization of a two-dimensional ising model,” Phys. Rev. , vol. 85, pp. 808–816, Mar 1952. [3] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller, “Equation of State Calculations by Fast Computing Machines,” Journal of Chemical Physics, vol. 21, pp. 1087–1092, June 1953.
-1 -1.5
0
100
200
300
400
500
600
700
800
900
1000
Number of Single Time Steps
Figure 11: Magnetisation versus time for varying fields at temperature T = 3.0 [ J /kB ]. Red, green, blue and magenta data points correspond to H = 1.0, 0.5, 0.1, -0.5 [ J /µ], respectively. It is possible to conclude that increasing the applied field increases the magnetisation, that is, spins become more aligned. Moreover, switching the sign of the applied field also switches the magentisation at equilibrium.
5
“Crystal statistics. i. a two-
CONCLUSIONS
A Monte Carlo method using the Metropolis algorithm was developed to study the prope rties of a two-dimensional lattice of spins. A continuous phase transition was identified at a critical temperature. Magnetisation, energy and heat capacity were observed to have the same temperature dependence for different values of N. Below this temperature, spins tend to align, whereas above it spins are randomly oriented at equilibrium. The criti cal temp erature was found to decrease with 1 /N 2 and was estimated as Tc = 2 .272 0.004 in the thermodynamic limit N ∞. The universal exponent β was determined as β = 0.109 0.018. When a magnetic field was applied, only one stable equilibrium state was found below the critical tempe rature. Above this temperature, the magnetisation was determined to be non-zero and dependent on the magnitude and sign of the field.
→
± ±
7
6
Appendix: Program Script
6.1
Matrix.h
// Class of matrices class Matrix { private: std::vector m_data; int m_rows; int m_cols; public: Matrix( int n_ro ws, int n_cols, double init = 0.0 ) : m_data( n_rows * n_cols, init ), m_rows( n_rows ), m_cols( n_cols ) { } ˜Matrix(){} double& operator() ( int i, int j ) { return m_dat a[i * m_cols + j]; } };
6.2
ising_bm442.cc
// Computing Project, Part II Physics // CRSID bm442 // Project A: The Ising Model of a Ferromagnet // // // // // // // // // // // // // // // //
The program asks to inse rt the mode: 1 / 2. 1
Takes M single tim e steps, sav es magnetization at e ach ti me i nstant i n fi le " M_vs_t.txt" an d automatically shows lattice after the M steps by calling Gnuplot.
2
Covers a range of tem peratures and saves the magnetisation per spin, the energy per spin, their standard deviations and the heat capacity in file "output.txt"
The initi al temperature and field , the size of the lattice and the number of single steps for mode 1 are set in the main function.
#include
8
#include #include #include #include #include #include #include #include #include
"gsl/gsl_rng.h" "Matrix.h"
#define k_B 1.0
// Boltzmann Constant
// Function that calculates the product of spins for nearest neighbours of one site double nearest_neighbours( Matrix& lattice, int row, int col, int N) { double nn_sum = 0.0; int a, b, c, d; // Periodic Boundary Conditions if ( col - 1 < 0) { a = N - 1; } else { a = col - 1; } if ( col + 1 > N - 1) { b = 0; } else { b = col + 1; } if ( row - 1 < 0) { c = N - 1; } else { c = row - 1; } if ( row + 1 > N - 1) { d = 0; } else
9
{ d = row + 1; } // 4 nearest neighbours nn_sum nn_sum nn_sum nn_sum
= = = =
nn_sum nn_sum nn_sum nn_sum
+ + + +
((double)lattice( ((double)lattice( ((double)lattice( ((double)lattice(
row, row, row, row,
col col col col
)) )) )) ))
* * * *
((double)lattice( ((double)lattice( ((double)lattice( ((double)lattice(
row, a row, b c, col d, col
)); )); )); ));
return nn_sum; }
// Function that calculates the magnetization per site double magnetisation( Matrix& lattice, int N) { double M = 0.0; for ( int ro w = 0; row < N; ro w++ ) {
for ( int co l = 0; col < N; co l++ ) { M = M + ((double)lattice( row, col )); }
} // Dividing by Nˆ2 to get ener gy per latt ice site return M / (double)pow( N, 2 ); }
// Function that calculates the energy of a spin configuration double energ y( Matrix& lattice, int N, double H ) { double double double double
J miu nn_sum spin_sum
= 1.0; = 1.0; = 0.0; = 0.0;
// Nearest Neighbours’ Coupling Constant // Magnetic Moment of Electron // Sum over nearest neighbours (initialization) // Sum of spins (initialization)
for ( int ro w = 0; row < N; ro w++ ) { for ( int co l = 0; col < N; co l++ ) { nn_sum = nn_sum + nearest_neighbours( lattice, row, col, N ); spin_sum = spin_sum + ((double)lattice( row, col )); } }
10
// Avoiding double counting of nearest neighbours interactions nn_sum = nn_sum / 2.0; return (-J * nn_sum - H * miu * spin_sum); }
// Function that returns a random numbe r between 0 and 1 using the GSL function gsl_rng. double random_generator( unsigned long int seed = 0) { static gsl_r ng* rng = 0; if ( rng == 0 ) rng = gsl_rng_alloc( gsl_rng_default ); if ( seed != 0 ) gsl_rng_set( rng, seed ); return gsl_rng_uniform( rng ); }
// Function that flips spin void flip_spin( Matrix& { if (lattice( row, col { lattice( row, col ) } else { lattice( row, col ) } }
lattice, int row, int col ) ) == 1.0) = -1.0;
= 1.0;
// Metropolis Algorithm void Metropolis( Matrix& lattice, int row, int col, int N, double T, double H ) { // Creating new micr ostate with one spin flipp ed at ( row, col ) Matrix new_microstate( N, N ); for (int i = 0; i < N; i++ ) { for (int j = 0; j < N; j++ ) { new_microstate( i, j ) = lattice( i, j ); } }
11
// All spins equal to
previous configuration so far. Flip one spin:
flip_spin( new_microstate, row, col );
// Change in energy for transition double delta _E = energy( new_microstate, N, H ) - energy( lattice, N, H ); // Random number between 0 and 1 double p = random_generator( ); if ( delta_E < 0.0 ) { flip_spin( lattice, row, col ); } else if ( exp(-delta_E/(k_B * T)) > p ) { flip_spin( lattice, row, col ); } } // Function that executes Single Time Step in a regular sequence void single_Monte_Carlo_step_regular( Matrix& lattice, int N, double T, double H ) { // Going through all latt ice site s to decide whet her or not to flip spin for ( int ro w = 0; row < N; ro w++ ) { for ( int co l = 0; col < N; co l++ ) { Metropolis( lattice, row, col, N, T, H ); } } }
// Function that executes Single Time Step in a random sequence void single_Monte_Carlo_step_random( Matrix& lattice, int N, double T, double H ) { // Going through Nˆ2 latt ice site s to decide whet her or not to flip spin for ( int i = 0; i < N * N; i++ ) { int row = (int) (random_generator( ) * (double) N); int col = (int) (random_generator( ) * (double) N);
Metropolis( lattice, row, col, N, T, H );
12
} }
// Function that sets initial random set of spins void initialize_lattice_random( Matrix& lattice, int N ) { for ( int ro w = 0; row < N; ro w++ ) { for ( int co l = 0; col < N; co l++ ) { double q = random_generator( ); if ( q < 0.5 ) { lattice( row, col ) = 1.0; } else { lattice( row, col ) = -1.0; }
}
} }
// Function that sets initial uniform set of spins void initialize_lattice_uniform( Matrix& lattice, int N ) { for ( int ro w = 0; row < N; ro w++ ) { for ( int co l = 0; col < N; co l++ ) { lattice( row, col ) = 1.0; } } } void mode_1( Matrix& lattice, int N, int M, double T_0, double H_0 ) { double T = T_0; double H = H_0; std::ofstream f("init"); f f f f f
<< "t = 1 \n"; << "show grid\n"; << "set term x11 persist \n"; << "set palette defined (-1 \"black\", 1 \"white\")" << "\n"; << "se t xrange [1: " << N << "] \n";
f << "se t yrange [1: " << N << "] \n";
13
f << "se t xtics 0,1," << N << "\n "; f << "se t ytics 0,1," << N << "\n "; f << "plot ’lattice_data.txt’ matri x every ::1::(" << N << "+1) with image" << "\n"; f.close(); std::ofstream g("lattice_data.txt"); std::ofstream h("M_vs_t.txt"); for ( int i = 0; i < M; i++ ) { single_Monte_Carlo_step_regular( lattice, N, T, H ); h << i << " " << magn etisation( lattice, N ) << "\n"; if( i == M - 1 ) { // To ensure that we start counting spins from row 1 rather than 0 // I add these values, which will not be seen because the range // goes from 1 to N, thereby neglecting this 0th row for( int j = 0; j < N + 1; j++ ) { g << "1 "; } g << "\n" ; for ( int ro w = 0; row < N; ro w++ ) { // To ensure that we start counting spins from column 1, as above g << "1 "; for ( int co l = 0; col < N; co l++ ) { g << lat tice( row, col ) << " "; } g << "\n" ; } } } g.close(); h.close(); system("gnuplot init"); } void mode_2( Matrix& lattice, int N, int M, double T_0, double H_0 ) { // Initialization of Dummy Variables for calculations
14
double double double double double double double double double double double double double double double double double
magnetisation_sum = 0.0; energy_sum = 0.0; squared_magnetisation_sum = 0.0; squared_energy_sum = 0.0; mean_magnetisation = 0.0; mean_energy = 0.0; mean_energy_per_site = 0.0; mean_squared_magnetisation = 0.0 ; mean_squared_energy = 0.0; squared_std_magnetisation = 0.0; std_magnetisation = 0.0; squared_std_energy = 0.0; std_energy = 0.0; heat_capacity_1 = 0.0; heat_capacity_2 = 0.0; mean_energy_previous = 0.0; T_previous = 0.0;
// Initialization of Temperature and Magnetic Field double T = T_0; double H = H_0; std::ofstream j("output"); while ( T < 3.5 ) { // More Monte Carlo Steps required to attain equilibrium close to transition if ( T == T_0 ) { M = 100; } else { M = 1000; } // Reaching equilibrium for ( int i = 0; i < M; i++ ) { single_Monte_Carlo_step_regular( lattice, N, T, H ); } // Performing calculations of properties of ensemble for ( int j = 0; j < M; j++ ) { single_Monte_Carlo_step_regular( lattice, N, T, H );
magnetisation_sum = magnetisation_sum + magnetisation( lattice, N );
15
squared_magnetisation_sum = squared_magnetisation_sum + pow( magnetisation( lattice, N ), 2 ); energy_sum = energy_sum + energy( lattice, N, H ); squared_energy_sum = squared_energy_sum + pow( energy( lattice, N, H ), 2 ); } // Dividing by number of iterations M to get mean quantities mean_magnetisation = magnetisation_sum / (double) M; mean_energy = energy_sum / (double) M; mean_squared_magnetisation = squared_magnetisation_sum / (double) M; mean_squared_energy = squared_energy_sum / (double) M;
mean_energy_per_site = mean_energy / (double) (N * N);
squared_std_magnetisation = fabs(mean_squared_magnetisation - pow( mean_magnetisation, 2 )); std_magnetisation = sqrt(squared_std_magnetisation); squared_std_energy = fabs(mean_squared_energy - pow( mean_energy, 2 )); std_energy = sqrt(squared_std_energy);
// Heat Capacity = dE/dT, but we require two values of T, so it has to wait for 2nd loop if( T > T_0 ) { heat_capacity_1 = (mean_energy - mean_energy_previous) / (T - T_previous); } // Heat Capacity from Fluctuation-Dissipation Theorem heat_capacity_2 = squared_std_energy / pow( T, 2 );
// Printing to file j << T << " " << fabs(mean_magnetisation) << " " << mean_energy_per_site << " "; j << std_magnetisation << " " << std_energy << " "; j << heat_capacity_1 << " " << heat_capacity_2 << "\n";
// Previous Values of T and are requ ired to compute Heat Capacity = dE/dT mean_energy_previous = mean_energy;
16
T_previous = T; // Reinitialization of Summing Variables energy_sum = 0.0; magnetisation_sum = 0.0; squared_energy_sum = 0.0; squared_magnetisation_sum = 0.0;
// Plots have most of their distinctive features close to critical temperature // so more data points are required in its vici nity if ( T < 2.0 || T > 3.0 ) { T = T + 0.1; } else {
T = T + 0.05;
} } } int main() { int int Matrix double double int
N = 10; M = 1000; lattice( N, N ); T_0 = 1.5; H_0 = 0.0; mode;
std::cout << "Insert mode: "; std::cin >> mode; // Initializing lattice initialize_lattice_uniform( lattice, N ); switch (mode) { case 1: mode_1( lattice, N, M, T_0, H_0 ); break; case 2: mode_2( lattice, N, M, T_0, H_0 ); break; default:
17
// Size of Spin Lattice // Number of Single Monte Carlo Steps // Square Lattice // Initial Temperature // Initial External Magnetic Field // Selects mode of program
std::cout << "No valid mode selected\n"; break; } return 0; }
18
