A lot happened
This commit is contained in:
@@ -0,0 +1 @@
|
||||
We have shown, that our implementation of the twodimensional Ising Model without the presence of a magentic field behaves, as expected from theoretical predictions. Especially showing relative errors in the order of \qtyrange{e-3}{e-4}{\percent} in quantities that are linear in the spin quantum numbers. Furthermore the statistical phenomena of phase change at the critical temperature has been proven. In future work it may be possible to improve the accuracy of the phase change behavior given more computing resources. We additionaly propose a parallelization approach for calculating the system evolution for multiple systems concurrently. We obtain a parallel efficiency of $\mathcal{E} = \num{.61}$ with the potential of further speed up under ideal conditions.
|
||||
@@ -0,0 +1,5 @@
|
||||
The Ising model is a mathematical representation of ferromagnetism in statistical mechanics. This model consists of discrete variables, known as spins, which can take on one of two values, typically represented as +1 (up) or -1 (down). These spins are arranged on a lattice and interact with their nearest neighbors. The model is used to study phase transitions, such as the transition from a magnetized to a non-magnetized state, and has applications in physics, biology and the social sciences.
|
||||
|
||||
It allows us to measure the statistical properties of the system, such as magnetization, susceptibility and specific heat, as functions of temperature. In order to efficiently sample the Ising model, we approximate the system state using Markov Chain Monte Carlo (MCMC) methods with rejection sampling. These are a class of algorithms that sample complex probability distributions by constructing a Markov chain with the desired distribution as its equilibrium distribution. In the context of the Ising model, these methods enable us to generate spin configurations representative of the system's behavior at a given temperature. Parallel computation can significantly speed up the sampling process, enabling the concurrent evaluation of multiple systems. This is particularly useful for large numbers of computations under different conditions, such as at various temperatures.
|
||||
|
||||
In \autoref{sec:methods}, the methods employed to implement MCMC sampling of the Ising model are delineated, including the specific algorithms employed and the parallelization techniques utilized. In \autoref{sec:results_and_discussion} the results of the aforementioned simulations are presented, including measurements of key statistical properties of the Ising model. The implications of these findings are discussed as well. Finally, in \autoref{sec:conclusion} we offer a summary of our work and put forward several suggestions for future research directions.
|
||||
@@ -0,0 +1,97 @@
|
||||
\subsection{Ising Model}
|
||||
In a two-dimensional system of size $L \times L$, each lattice site $i$ contains a spin variable $s_i$ that can take on values of either +1 (up) or -1 (down). The energy of a given configuration of spins is described by the sum over nearest-neighbor interactions:
|
||||
\begin{equation}
|
||||
E = -J \sum_{\langle i,j \rangle} s_i s_j.
|
||||
\end{equation}
|
||||
|
||||
The coupling constant, $J$, is the parameter that quantifies the strength of the interaction between neighboring spins. The sum is taken over all pairs of nearest neighbors, defined as $\langle i,j \rangle$. In order to facilitate the process of summation without the occurrence of double counting, the total energy can be expressed in the following manner:
|
||||
\begin{equation} \label{eq:total_energy}
|
||||
E = -\frac{J}{2} \sum_{i} \sum_{j \in \text{NN}(i)} s_i s_j,
|
||||
\end{equation}
|
||||
where $\text{NN}(i)$ represents the set of nearest neighbors of site $i$. Depending on the configuration of the neighboring spins, the flipping of a single spin, $s_k$, will result in a change in the energy of the system by an amount
|
||||
\begin{equation}
|
||||
\Delta E = \Delta s_k \cdot J \sum_{j \in \text{NN}(k)} s_j,
|
||||
\end{equation}
|
||||
where $\Delta s_k = s_k^{\text{new}} - s_k^{\text{old}} \in \{-2, 2\}$ is the change in the spin variable at site $k$. To substantiate this claim, it is essential to recognize that the energy contribution from site $k$ prior to the flip is given by:
|
||||
\begin{equation}
|
||||
E_{\text{old}} = -J s_k^{\text{old}} \sum_{j \in \text{NN}(k)} s_j.
|
||||
\end{equation}
|
||||
The discrepancy in the factor 2 in comparison to \autoref{eq:total_energy} emerges due to the fact that $s_k$ is a constituent of the sum over nearest neighbors for each of its neighbors, resulting in double counting. Subsequent to the flipping of the spin at site $k$, the contribution to the energy from that site is updated, and after the flip, it becomes
|
||||
\begin{equation}
|
||||
E_{\text{new}} = -J s_k^{\text{new}} \sum_{j \in \text{NN}(k)} s_j.
|
||||
\end{equation}
|
||||
The alteration in energy resulting from the flip is then
|
||||
\begin{equation}
|
||||
\Delta E = E_{\text{new}} - E_{\text{old}} = -J (s_k^{\text{new}} - s_k^{\text{old}}) \sum_{j \in \text{NN}(k)} s_j.
|
||||
\end{equation}
|
||||
This expression demonstrates how the energy change depends exclusively on the local configuration of spins surrounding site $k$.
|
||||
|
||||
The total magnetization of the system $M$ is derived from the summation of all spin variables
|
||||
\begin{equation}
|
||||
M = \sum_{i} s_i.
|
||||
\end{equation}
|
||||
Both extensive quantities, defined as those that scale with the system size, can be normalized per spin to obtain intensive quantities
|
||||
\begin{equation}
|
||||
\epsilon = \frac{E}{N}, \quad m = \frac{M}{N},
|
||||
\end{equation}
|
||||
where $N = L^2$ is the total number of spins in the lattice.
|
||||
|
||||
The Ising model displays a specific heat capacity, $C_V$, and magnetic susceptibility, $\chi$, both of which are defined as follows, once more normalized per spin:
|
||||
\begin{equation} \label{eq:heat_capacity}
|
||||
\frac{C_V}{N} = \frac{N \left(\expect{\epsilon^2} - \expect{\epsilon}^2\right)}{k_B T^2},
|
||||
\end{equation}
|
||||
\begin{equation} \label{eq:susceptibility}
|
||||
\frac{\chi}{N} = \frac{N \left(\expect{m^2} - \expect{m}^2\right)}{k_B T},
|
||||
\end{equation}
|
||||
where $k_B$ is the Boltzmann constant and $T$ is the temperature of the system. The notation $\expect{\cdot}$ denotes the expectation value, which can be estimated through averaging over sampled configurations. The average over a sufficiently large number of Monte Carlo samples will converge to the true expectation value. Given that expectation values are computed from sampled configurations, it is imperative to ensure that the system's state space is sampled correctly. In the Ising model, it is assumed that the configurations follow the Boltzmann distribution. This distribution dictates the probability of a configuration with energy $E$ at temperature $T$ as follows:
|
||||
\begin{equation}
|
||||
P(E) = \frac{e^{-E/(k_B T)}}{Z},
|
||||
\end{equation}
|
||||
with the partition function $Z$ defined as
|
||||
\begin{equation}
|
||||
Z = \sum_{\text{all states}} e^{-E/(k_B T)}.
|
||||
\end{equation}
|
||||
|
||||
\subsubsection{Critical Temperature}
|
||||
In two dimensions, the Ising model undergoes a continuous phase transition at the critical temperature $T_C$ separating an ordered low-temperature phase from a disordered high-temperature phase. In the context of the model, spontaneous magnetization is exhibited when the temperature falls below the critical temperature, $T_C$. This phenomenon occurs because, even in the absence of an external magnetic field, the magnetization $m$ remains non-zero. As the temperature approaches $T_C$ from below, the magnetization decreases continuously, eventually vanishing according to a characteristic power law. Concurrently, fluctuations in energy and magnetization exhibit a marked increase, as evidenced by the divergent behavior of both the heat capacity \autoref{eq:heat_capacity}, and the susceptibility \autoref{eq:susceptibility}. Precisely at $T_C$, the system becomes scale-invariant, and configurations manifest correlations that extend across the entire lattice. In numerical simulations, such phenomena manifest as strong critical slowing down, necessitating sufficiently long sampling times to resolve the sharply peaked response functions. The exact solution for the square lattice with nearest-neighbor interactions yields the following result:
|
||||
$k_B T_C / J = 2 / \ln(1 + \sqrt{2}) \approx \num{2.269}$, providing a precise reference point for validating numerical results.
|
||||
|
||||
\subsubsection{Ergodic Property}
|
||||
The implementation of Monte Carlo methods within the Ising model is contingent upon the assumption of ergodicity, which stipulates that the Markov chain generated by successive spin updates must possess the capacity to attain all physically relevant configurations within a sufficient time frame. In practice, this implies that single-spin flip dynamics must form a connected set over the configuration space, thereby avoiding traps in metastable regions. The Metropolis–Hastings update rule, which is employed for sampling purposes, ensures detailed balance by ensuring that the stationary distribution of the Markov chain corresponds to the Boltzmann distribution. Ergodicity guarantees the convergence of long-time averages of observables, such as $\epsilon$ and $m$, to ensemble averages. As the critical point is approached, ergodicity remains valid; however, the time required to decorrelate successive configurations grows significantly due to the emergence of large correlated domains. This phenomenon, often referred to as critical slowing down, necessitates meticulous care in determining autocorrelation times to obtain statistically independent samples and reliable estimates of thermodynamic quantities.
|
||||
When an ergodic Markov chain is considered, the equivalence between ensemble averages and long-time averages is guaranteed in the limit of a sufficient number of Monte Carlo steps. For any observable $A$, the ensemble expectation value $\expect{A}$ obtained from the Boltzmann distribution coincides with the time-averaged estimate $\hat{A}$ extracted from the simulation. It is imperative to note that the aforementioned equality is only valid once the chain has thoroughly explored configuration space. This underscores the significance of ergodicity and adequate sampling length as essential prerequisites for reliable numerical measurements. However, it does allow for the utilization of an adequately sampled time average in place of the expectation value. The time-averaged estimate is defined as follows:
|
||||
\begin{equation}
|
||||
\hat A(t) = \frac{1}{t} \int_0^t \mathrm{d} t' A(t').
|
||||
\end{equation}
|
||||
|
||||
|
||||
\subsection{Markov Chain Monte Carlo}
|
||||
To ensure the proper sampling of the Boltzmann distribution and to facilitate efficient exploration of the configuration space, we employ Markov Chain Monte Carlo (MCMC) methods. MCMC methods generate a sequence of configurations (states) where each configuration depends only on the previous one, forming a Markov chain. The fundamental objective is to establish a transition probability between states that satisfies detailed balance, thereby ensuring the system's convergence to the desired equilibrium distribution. In each Monte Carlo cycle, the objective is to flip each spin in the lattice precisely once. The acceptance of a proposed spin flip is determined using the Metropolis-Hastings algorithm, which involves calculating the energy change, denoted by $\Delta E$, associated with the flip. In the event that $\Delta E \leq 0$ the flip is accepted without condition, as it results in a lower energy state. In the event that $\Delta E > 0$, the flip is accepted with a probability that is determined by the Boltzmann factors
|
||||
\begin{equation} \label{eq:acceptance_prob}
|
||||
P_{\text{accept}} = e^{-\Delta E/(k_B T)}.
|
||||
\end{equation}
|
||||
This acceptance criterion ensures that the system can explore higher energy states, allowing it to escape local minima and sample the configuration space more effectively. The acceptance of the flip is contingent upon the random number $r$ being drawn uniformly from the interval $[0, 1)$. More specifically, the flip is accepted if $r < P_{accept}$.
|
||||
|
||||
\subsection{Parallelization}
|
||||
Parallelization of MCMC methods is challenging due to the sequential nature of Markov chains. However, it is possible to exploit the independence of multiple simulations at different temperatures or initial conditions to achieve parallelism. The implementation of multiple independent Markov chains operating concurrently, with each chain simulating the Ising model at an distinct temperature, has been demonstrated to result in a substantial acceleration of the overall sampling process. Each chain functions independently, facilitating straightforward parallelization across multiple processors. This approach is particularly beneficial when studying temperature-dependent properties of the Ising model, as it enables simultaneous exploration of the system's behavior across a range of temperatures.
|
||||
|
||||
\subsection{Implementation Details}
|
||||
|
||||
The simulations were executed using a custom C++ implementation of the two-dimensional Ising model with periodic boundary conditions. Each simulation instance constructs an $L\times L$ lattice with spins $s_{ij} \in \{\pm 1\}$. The lattice is initialized in one of two ways: either in a fully ordered configuration or with randomly assigned spins drawn from a uniform distribution. The generation of random numbers is achieved via a Mersenne Twister engine (\texttt{std::mt19937}), with each run being seeded individually. To ensure that periodic boundary conditions are met, the neighboring indices are computed using modulo arithmetic
|
||||
\begin{equation}
|
||||
\begin{aligned}
|
||||
\text{NN}(i,j) = \{&((i+1) \mod L, j), \\&((i-1 + L) \mod L, j), \\&(i, (j+1) \mod L), \\&(i, (j-1 + L) \mod L)\}.
|
||||
\end{aligned}
|
||||
\end{equation}
|
||||
|
||||
The spin dynamics follow the standard single-spin Metropolis algorithm. A Monte Carlo step involves the selection of a lattice site uniformly at random, the computation of the local energy change $\Delta E$ associated with the flipping of that spin, and the acceptance of the flip with probability $P_\text{accept}$. In order to minimize computational expenditure, all conceivable Boltzmann factors (\cref{eq:acceptance_prob}) associated with the four-neighbor interaction are initially calculated and subsequently stored in a lookup table.
|
||||
|
||||
During the simulation, the instantaneous energy and absolute magnetization are evaluated after each attempted update. The values and their squares are then accumulated to compute ensemble averages and variances. Subsequently, the code derives the heat capacity and magnetic susceptibility from these moments, as outlined in \cref{eq:heat_capacity,eq:susceptibility}. Optional time-resolved data-energy, magnetization, and other intermediate quantities, are stored at each Monte Carlo cycle when evolution recording is enabled.
|
||||
|
||||
Furthermore, to support large parameter sweeps, the code incorporates a scheduling layer that manages multiple independent simulations. The configuration of each run is delineated by a specific lattice size, a designated temperature, the number of Monte Carlo cycles employed, the initialization mode, and the RNG seed. A master seed generator ensures statistical independence across runs. The scheduler distributes the runs across a number of CPU threads specified by the user, executing them concurrently. Upon completion of a given run, an output file is generated, containing the averaged observables. When the relevant option is enabled, an additional file is also generated, which contains the full time evolution of the system.
|
||||
|
||||
|
||||
\subsubsection{Tools and Usage of AI}
|
||||
For concurrent execution, the implementation utilizes the C++ library \texttt{OpenMP}. The executables were compiled using \texttt{g++ (GCC) 15.2.1 20251022 (Red Hat 15.2.1-3)} with optimization flag \texttt{-O3} enabled. All executions were carried out on an AMD Ryzen 7 PRO 4750U Processor, equipped with 16 GB of RAM, operating on the Fedora Linux 42 distribution. A large impact of IO operations on performance is not anticipated, given the utilization of an NVMe SSD for all read and write operations.
|
||||
|
||||
For the purposes of analysis, the \texttt{Python 3.13} programming language was utilized in conjunction with the libraries \texttt{NumPy}, \texttt{Pandas}, \texttt{SciPy}, and \texttt{Matplotlib}. In the course of composing this report, the large language model \texttt{ChatGPT} was employed to proofread various sections and to propose enhancements in terms of phrasing and grammar. All figures were created manually, without the use of AI assistance. During the implementation of the code, the large language model provided by \texttt{GitHub Copilot} was used to improve coding speed, but not for generating large code blocks. A comprehensive review and testing of the code was conducted to ensure its correctness.
|
||||
|
||||
@@ -0,0 +1,100 @@
|
||||
\subsection{Validation against Analytical Results}
|
||||
In \cref{app:analytical} of the appendix, a summary of the expected behavior of a correctly sampled lattice of $L=2$ is provided, in terms of its expectation values. Subsequent to the generation of said results, a comparison is formulated against a set of 10 simulations, which were executed for a duration of \num{10000} Monte Carlo cycles at a temperature of $k_B T = \num{1.0}$. To assess the uncertainty associated with each ensemble mean, the standard deviation of the ensemble values is employed. The uncertainty in the relative error is derived through the application of Gaussian error propagation to the ensemble variances. The results are presented in \cref{tab:comparison} and demonstrate a high degree of congruence with the anticipated analytical expectation values for the energy and magnetization of the system. Both quantities manifest relative deviations in the order of \qtyrange{e-3}{e-4}{\percent}. The comparison indicates discrepancies of \qty{3.4959}{\percent} for the specific heat capacity of the system, accompanied by a substantial uncertainty range of $>\qty{30}{\percent}$ . In a similar vein, the relative error of the susceptibility demonstrates a large deviation of $\qty{4.2046+-35.9493}{\percent}$. The significant uncertainty hinders the rise of any tension between the theoretical prediction and the simulated quantity. Consequently, the observed alignment is deemed to be highly satisfactory in accordance with theoretical expectations, thereby validating the implementation of the model.
|
||||
|
||||
A comparatively extended simulation time frame was selected for the lattice size. However, a convergence analysis reveals that a minimum of \num{1000} steps are necessary to achieve results that are equivalent in quality. In order to execute this analysis, the deviation between time averages $\hat \epsilon, |\hat m|$ and theoretical predictions $\expect{\epsilon}, \expect{|m|}$ as a function of the Monte Carlo cycle for each of the simulations. The results are displayed in \cref{fig:convergence_2x2}. Following a brief period of equilibration, the presence of artifacts becomes evident within the range of \numrange{e-2}{e-3}, with an upper boundary delineating the maximum deviation. This observation signifies the stability of the system within its equilibrium state.
|
||||
|
||||
\begin{figure}
|
||||
\centering
|
||||
\includegraphics[width=\columnwidth]{../../src/project4/python/convergence_analysis.pdf}
|
||||
\caption{Evolution of simulated 2x2 lattice structures over the duration of \num{10000} Monte Carlo cycles. Comparison of ensemble averages $\hat \epsilon, |\hat m|$ and theoretical predictions $\expect{\epsilon}, \expect{|m|}$. We use the values from \cref{tab:analytical_expectation_values} as the prediction values.}
|
||||
\label{fig:convergence_2x2}
|
||||
\end{figure}
|
||||
|
||||
|
||||
\begin{table}
|
||||
\caption{Comparison of ensemble properties with results obtained from analytical calculations in \cref{tab:analytical_expectation_values}. Ensemble properties sampled from 10 simulations with $L=2$, \num{10000} Monte Carlo cycles and random initialization.}
|
||||
\label{tab:comparison}
|
||||
\begin{tabular}{ccc}
|
||||
\toprule
|
||||
Quantity & Ensemble Average & Relative Error\\
|
||||
\midrule
|
||||
$\expect{\epsilon}$ & \num{-1.9958+-0.0014} & \qty{0.0071+-0.0692}{\percent} \\
|
||||
$\expect{|m|}$ & \num{0.9987+-.0005} & \qty{0.0004+-0.0463}{\percent} \\
|
||||
$\expect{\frac{C_V}{N}}$ & \num{0.033+-.011} & \qty{3.4959+-34.3238}{\percent} \\
|
||||
$\expect{\frac{\chi}{N}}$ & \num{0.0038+-.0014} & \qty{4.2046+-35.9493}{\percent} \\
|
||||
\bottomrule
|
||||
\end{tabular}
|
||||
\end{table}
|
||||
|
||||
\subsection{Equilibration of the System} \label{subsec:equil}
|
||||
In order to provide further validation of the working principle of our implementation, it is possible to evaluate the equilibration/thermalization process. The system evolution is recorded by initializing the system from either a fully ordered state, where all the spins are parallel, or a fully disordered state, where the spin directions are chosen from an entirely uncorrelated uniform distribution. At a constant temperature, a specific degree of ordering is anticipated to arise. Below the critical temperature of the system, the energetic advantage of ordering is sufficient to cause nearly complete sorting of the system. At temperatures above the critical temperature, a disruption in the microscopic ordering of spins becomes possible, resulting in the loss of magnetization at the macroscopic level of the system. The assessment of disruptions to the system's optimal order necessitates the utilization of the normalized energy of the system. An optimal ordering of the system will result in minimal energy, while any disruption to this ordered structure will result in an increase in energy. A steady normalized energy indicates successful equilibration, i.e., the energy is in a steady state relative to the system's temperature. In order to illustrate the discrepancies occurring on very small time scales in comparison to the behavior observed on the long time scale, the following quantities are graphed: $\epsilon$ is the current normalized energy of the system at a specific time step and the averaged normalized energy $\hat\epsilon$:
|
||||
\begin{equation}
|
||||
\hat \epsilon(t) = \frac{1}{t} \int_0^t \mathrm{d} t' \epsilon(t').
|
||||
\end{equation}
|
||||
|
||||
For a system of $L=20$, a total of \num{2000} Monte Carlo cycles are simulated and the normalized energies $\epsilon, \hat \epsilon$ are recorded. This procedure is repeated for two temperatures: $k_BT = \num{1.0}$ and $k_B T = \num{2.4}$ and the two extremes of ordering, i.e., fully ordered and fully randomized initial configurations.
|
||||
|
||||
\begin{figure}
|
||||
\centering
|
||||
\includegraphics[width=\columnwidth]{../../src/project4/python/energy_convergence.pdf}
|
||||
\caption{Evolution of the normalized current energy $\epsilon$ (slightly colored line) and the normalized time-averaged energy $\hat{\epsilon}$ (fully colored lines) for systems of lattice size $L=20$. }
|
||||
\label{fig:energy_evolution}
|
||||
\end{figure}
|
||||
|
||||
The results of this study are shown in \cref{fig:energy_evolution}. The observed behavior is consistent with the anticipated outcomes. Specifically, it has been determined that the system approaches a state of equilibrium energy for both temperatures after approximately \num{1000} Markov Chain Monte Carlo cycles. When the fully unordered state is employed as the initial condition, a higher-than-equilibrium energy state is achieved, accompanied by a rapid decay towards the thermalized state. The temperature of $k_B T = 1$ results in complete ordering of the system, irrespective of the initial configuration. For lower temperatures, fewer stochastic deviations from equilibrium are observed, indicating that the system remains predominantly ordered even on short time scales. For the higher temperature point that was sampled, there are noticeable deviations from the time-averaged energy $\hat \epsilon$. On the long-time-scale average, the system nevertheless approaches an equilibrium state with respect to the energy. In summary, the model's behavior aligns with the theoretical framework.
|
||||
|
||||
\subsection{Phase Transition Behavior}
|
||||
In order to observe the transition between the two phases of the two-dimensional Ising model, a simulation was conducted in which randomly initialized systems were observed at different temperature points for \num{50000} MCMC cycles. A range of temperatures, from $k_B T = \numrange{1.0}{2.4}$ is scanned at an offset between simulations of $k_B \Delta T = \num{1e-3}$. In order to accelerate the simulation, no output is recorded during the simulation, and the systems parameters are extracted solely at the end of the simulation. This procedure is carried out for four lattice sizes, namely for $L \in \{40, 60, 80, 100\}$. The resulting system parameters are then plotted as a function of the set temperature and the lattice size. A plot containing the observables can be found in \cref{fig:observables}. The average energy of the system increases linearly with the temperature, irrespective of the system's size. For temperatures above the critical temperature, the system exhibits minimal spread, with larger systems demonstrating increased energies. In the context of absolute magnetization, a decline is observed above the critical temperature threshold. This phenomenon can be attributed to the disruption of the global magnetization, which results in the manifestation of averaging effects within the lattice. The magnitude of the drop is more pronounced with increasing lattice sizes, as the averaging effect is statistically more significant. For temperatures below the critical temperature, a decrease in average absolute magnetization with an increase in temperature is observed, as the added dynamic in the system reduces perfect ordering. The specific heat capacity and the susceptibility demonstrate a diverging tendency around the critical temperature. A comparative analysis of the temperature-dependent specific heat capacity reveals a uniform trend across all systems, with the exception of a shift in the point of divergence. This variation in the critical temperature is attributed to the influence of the lattice size on the system's behavior. As the lattice size increases, the degree of divergence concomitantly rises.
|
||||
|
||||
\begin{figure}
|
||||
\centering
|
||||
\includegraphics[width=\columnwidth]{../../src/project4/python/observable_vs_temperature.pdf}
|
||||
\caption{Calculated normalized energy $\expect{\epsilon}$, magnetization $\expect{|m|}$, specific heat capacity $\expect{\frac{C_V}{N}}$ and susceptibility $\expect{\frac{\chi}{N}}$ as a function of system temperature for different lattice sizes. Systems were evolved for \num{50e3} MCMC cycles.}
|
||||
\label{fig:observables}
|
||||
\end{figure}
|
||||
|
||||
\subsubsection{Limit of Infinite Lattice Size}
|
||||
The diverging behavior is used in a second step to approximate the point of the critical temperature for each lattice size. For both the specific heat capacity and the susceptibility data a fit to the function
|
||||
\begin{equation} \label{eq:fitfunc}
|
||||
f(T) = A e^{-\frac{(T-\mu)^2}{2\sigma^2}} + C
|
||||
\end{equation}
|
||||
is performed. As $\mu$ denotes the maximum point of the Gaussian function, it is employed as an estimator for the critical temperature. The results of all fits are summarized in \cref{fig:fit}. The fit is performed using the function \texttt{optimize.curve\_fit} from the \texttt{scipy} package, with the minimization of squared residuals. Subsequently, the values obtained for the critical temperature in dependence of the lattice size are employed to approximate the critical temperature of an infinitely sized lattice. To elaborate, the critical temperature is plotted against the lattice size, and a secondary least squares minimization is performed on the relation.
|
||||
\begin{equation}
|
||||
T_C(L) = T_C(\infty) + a L^{-1}.
|
||||
\end{equation}
|
||||
It is evident that two distinct values are obtained for the critical temperature of infinite lattices. Specifically, the first value is determined to be \num{2.284} and is derived from the data points extracted from the susceptibility. The second value, \num{2.271}, is derived from the measurements on the specific heat capacity. The fit to obtain those results is shown in \cref{fig:critical_temp_fit}. The data points obtained via the specific heat capacity demonstrate a high degree of agreement with the anticipated result, $k_B T_C = 2.269$. The discrepancy between the measurements obtained from the susceptibility suggests the necessity for a more sophisticated function to ascertain the critical temperatures from the simulations, i.e., improving \cref{eq:fitfunc}. In order to achieve more precise measurements, it is necessary to increase the number of samples taken in the region of the critical temperature. Due to the scarcity of computing resources, this step has not been implemented in the present study.
|
||||
|
||||
\begin{figure}
|
||||
\centering
|
||||
\includegraphics[width=\columnwidth]{../../src/project4/python/phase_transition_analysis.pdf}
|
||||
\caption{Fit to the specific heat capacity $\expect{\frac{C_V}{N}}$ and susceptibility $\expect{\frac{\chi}{N}}$ as a function of system temperature for different lattice sizes. Fit to \cref{eq:fitfunc}.}
|
||||
\label{fig:fit}
|
||||
\end{figure}
|
||||
|
||||
\begin{figure}
|
||||
\centering
|
||||
\includegraphics[width=\columnwidth]{../../src/project4/python/critical_temperature_extrapolation.pdf}
|
||||
\caption{Fit of the values obtained from \cref{fig:fit} to the relation of the system size dependent critical temperature.}
|
||||
\label{fig:critical_temp_fit}
|
||||
\end{figure}
|
||||
|
||||
\subsection{Performance Analysis}
|
||||
The performance of the simulation was evaluated through benchmarking in both a single-threaded and a 12-thread multithreaded configuration using the system time utility. For each of the executions, a total of 500 identical systems of lattice size $L = \num{100}$ were simulated for \num{1000} Monte Carlo cycles. Large numbers of identical systems are deliberately selected to minimize the impact of a subset of threads that complete execution before others. This process enhances the overall run time without exhausting all available resources. In order to facilitate the uninterrupted execution of background processes, a configuration of 12 threads was selected on a system with a total of 16 threads. The single-threaded execution required a wall time of \qty{217.634}{\s}, while the multithreaded run completed in \qty{29.580}{\s}. This corresponds to a speedup of
|
||||
\begin{equation}
|
||||
\mathcal S = \frac{217.634}{29.580} \approx 7.36 ,
|
||||
\end{equation}
|
||||
which is substantial, though below the ideal linear scaling expected for perfectly parallel work on 12 threads. The resulting parallel efficiency is therefore $\mathcal E = \mathcal S/12 \approx 0.61$.
|
||||
|
||||
The reported CPU times provide additional insight into the system's performance. The single-threaded version accumulated \qty{215.592}{\s} of user time, consistent with full utilization of one core. In contrast, the multithreaded run consumed \qty{339.713}{\s} of user time over a wall interval of \qty{29.580}{\s}. This corresponds to an average of roughly $11.5$ concurrently active hardware threads. This finding suggests that the parallel implementation utilizes the majority of available cores during the majority of the runtime.
|
||||
|
||||
One potential explanation for the observed limited parallel efficiency is overhead due to concurrent processing, including task delegation, I/O, and an augmented thermal load on the CPU, which can result in thermal throttling. Given that the CPU utilized for the evaluation is part of a notebook computer, the final point is particularly salient, given the prevalence of inadequate cooling capacity in compact designs. Subsequent studies are required to ascertain the impact of all factors on the reduced parallel efficiency.
|
||||
|
||||
\subsection{Statistical Properties of the Ising Model}
|
||||
The evaluation of the probability distribution of the normalized energy is achieved through a comprehensive analysis of the distribution of epsilon for systems of $L=20$ at the two temperatures previously specified in \cref{subsec:equil}. For each set of initial conditions, 25 identical simulations were executed, each with a different seed, and the resulting trajectories were collected. The results are then displayed in a normalized histogram, using the counts of each unique energy value present in the dataset. The approximate probability density function is displayed in \cref{fig:pdf}.
|
||||
|
||||
\begin{figure}
|
||||
\centering
|
||||
\includegraphics[width=\columnwidth]{../../src/project4/python/energy_histograms.pdf}
|
||||
\caption{Normalized counts per discrete energy value of simulations using $L=20$. Total of 50 simulation of \num{2000} cycles.}
|
||||
\label{fig:pdf}
|
||||
\end{figure}
|
||||
Binary file not shown.
@@ -0,0 +1,187 @@
|
||||
% USEFUL LINKS:
|
||||
% -------------
|
||||
%
|
||||
% - UiO LaTeX guides: https://www.mn.uio.no/ifi/tjenester/it/hjelp/latex/
|
||||
% - Mathematics: https://en.wikibooks.org/wiki/LaTeX/Mathematics
|
||||
% - Physics: https://ctan.uib.no/macros/latex/contrib/physics/physics.pdf
|
||||
% - Basics of Tikz: https://en.wikibooks.org/wiki/LaTeX/PGF/Tikz
|
||||
% - All the colors! https://en.wikibooks.org/wiki/LaTeX/Colors
|
||||
% - How to make tables: https://en.wikibooks.org/wiki/LaTeX/Tables
|
||||
% - Code listing styles: https://en.wikibooks.org/wiki/LaTeX/Source_Code_Listings
|
||||
% - \includegraphics https://en.wikibooks.org/wiki/LaTeX/Importing_Graphics
|
||||
% - Learn more about figures: https://en.wikibooks.org/wiki/LaTeX/Floats,_Figures_and_Captions
|
||||
% - Automagic bibliography: https://en.wikibooks.org/wiki/LaTeX/Bibliography_Management (this one is kinda difficult the first time)
|
||||
%
|
||||
% (This document is of class "revtex4-1", the REVTeX Guide explains how the class works)
|
||||
% REVTeX Guide: http://www.physics.csbsju.edu/370/papers/Journal_Style_Manuals/auguide4-1.pdf
|
||||
%
|
||||
%
|
||||
% COMPILING THE .pdf FILE IN THE LINUX TERMINAL
|
||||
% ---------------------------------------------
|
||||
%
|
||||
% [terminal]$ pdflatex report_example.tex
|
||||
%
|
||||
% Run the command twice, always.
|
||||
%
|
||||
% When using references, footnotes, etc. you should run the following chain of commands:
|
||||
%
|
||||
% [terminal]$ pdflatex report_example.tex
|
||||
% [terminal]$ bibtex report_example
|
||||
% [terminal]$ pdflatex report_example.tex
|
||||
% [terminal]$ pdflatex report_example.tex
|
||||
%
|
||||
% This series of commands can of course be gathered into a single-line command:
|
||||
% [terminal]$ pdflatex report_example.tex && bibtex report_example.aux && pdflatex report_example.tex && pdflatex report_example.tex
|
||||
%
|
||||
% ----------------------------------------------------
|
||||
|
||||
|
||||
|
||||
% \documentclass[english,notitlepage,reprint,nofootinbib]{revtex4-2} % defines the basic parameters of the document
|
||||
\documentclass[english,notitlepage,reprint,nofootinbib]{revtex4-2} % defines the basic parameters of the document
|
||||
% If you want a single-column, remove "reprint"
|
||||
\usepackage{silence}
|
||||
\WarningFilter{revtex4-2}{Repair the float}
|
||||
|
||||
% Allows special characters (including æøå)
|
||||
\usepackage[utf8]{inputenc}
|
||||
\usepackage[english]{babel}
|
||||
|
||||
% Note that you may need to download some of these packages manually, it depends on your setup.
|
||||
% It may be usefult to download TeXMaker, because it includes a large library of the most common packages.
|
||||
|
||||
\usepackage{amssymb} % mathematical symbols (physics imports amsmath)
|
||||
\usepackage{amsmath}
|
||||
\usepackage{graphicx} % include graphics such as plots
|
||||
\usepackage{xcolor} % set colors
|
||||
\usepackage{hyperref} % automagic cross-referencing
|
||||
\usepackage{listings} % display code
|
||||
\usepackage{subfigure} % imports a lot of cool and useful figure commands
|
||||
\usepackage{float}
|
||||
%\usepackage[section]{placeins}
|
||||
\usepackage{algorithm}
|
||||
\usepackage{booktabs}
|
||||
\usepackage[noend]{algpseudocode}
|
||||
\usepackage{subfigure}
|
||||
\usepackage{tikz}
|
||||
\usepackage[noabbrev,nameinlink]{cleveref}
|
||||
\usepackage[separate-uncertainty=true]{siunitx}
|
||||
% Define amu = Dalton
|
||||
\DeclareSIUnit\amu{u}
|
||||
|
||||
\usepackage{todonotes}
|
||||
\usetikzlibrary{quantikz}
|
||||
% defines the color of hyperref objects
|
||||
% Blending two colors: blue!80!black = 80% blue and 20% black
|
||||
\hypersetup{ % this is just my personal choice, feel free to change things
|
||||
colorlinks,
|
||||
linkcolor={red!50!black},
|
||||
citecolor={blue!50!black},
|
||||
urlcolor={blue!80!black}}
|
||||
|
||||
|
||||
% ===========================================
|
||||
|
||||
\newcommand{\expect}[1]{\left\langle#1\right\rangle}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\title{Efficient Monte-Carlo Sampling of Ising Models} % self-explanatory
|
||||
\author{Lars Bogner} % self-explanatory
|
||||
\date{\today} % self-explanatory
|
||||
\affiliation{{\url{https://github.uio.no/larsbog/FYS4150}}} % ignore this, but keep it.
|
||||
|
||||
%This is how we create an abstract section.
|
||||
\begin{abstract}
|
||||
This work presents a numerical study of the two-dimensional Ising model using a custom implementation of Markov Chain Monte Carlo sampling with single-spin Metropolis updates. The model is simulated on square lattices with periodic boundary conditions, and the statistical properties of the system, i.e., energy, magnetization, heat capacity, and susceptibility, are extracted across a range of temperatures. Multiple independent Markov chains are executed in parallel to accelerate parameter sweeps, enabling efficient exploration of temperature-dependent behavior. Validation against analytical results for the $L=2$ lattice demonstrates excellent agreement for energy and magnetization, confirming the correctness of the implementation. Detailed equilibration studies highlight the role of ergodicity and illustrate the contrasting dynamical behavior below and above the critical temperature. The simulations reproduce the characteristic signatures of the continuous phase transition at $T_C$, including the emergence of large correlated domains and the divergence of response functions. The methods and performance characteristics of the implementation provide a robust foundation for further investigations of critical phenomena in lattice spin systems.
|
||||
\end{abstract}
|
||||
\maketitle
|
||||
\tableofcontents
|
||||
|
||||
|
||||
% ===========================================
|
||||
\section{Introduction}
|
||||
\input{chapters/introduction.tex}
|
||||
|
||||
% ===========================================
|
||||
\section{Methods}\label{sec:methods}
|
||||
%
|
||||
\input{chapters/methods.tex}
|
||||
% ===========================================
|
||||
\section{Results and discussion}\label{sec:results_and_discussion}
|
||||
%
|
||||
\input{chapters/results.tex}
|
||||
|
||||
% ===========================================
|
||||
\section{Conclusion}\label{sec:conclusion}
|
||||
|
||||
\input{chapters/conclusion.tex}
|
||||
% ===========================================
|
||||
\appendix
|
||||
\section{Analytical expectation values for the 2x2 lattice} \label{app:analytical}
|
||||
|
||||
\begin{table}
|
||||
\caption{Possible configurations of a 2x2 lattice with periodic boundary conditions. Assuming $|s_i| = 1$ and properties of each configuration.}
|
||||
\label{tab:2x2_configurations}
|
||||
\begin{tabular}{cccc}
|
||||
\toprule
|
||||
$\sum_i \delta_{s_i,1}$ & Degeneracy & Energy $E$ & Magnetization $M$ \\
|
||||
\midrule
|
||||
0 & 1 & -8J & -4 \\
|
||||
1 & 4 & 0 & -2 \\
|
||||
2 & 4 & 0 & 0 \\
|
||||
2 & 2 & 8J & 0 \\
|
||||
3 & 4 & 0 & 2 \\
|
||||
4 & 1 & -8J & 4 \\
|
||||
\bottomrule
|
||||
\end{tabular}
|
||||
\end{table}
|
||||
|
||||
In \cref{tab:2x2_configurations}, a comprehensive list of all possible configurations of a $2\times 2$ lattice with periodic boundary conditions is provided, along with their respective degeneracies, energies, and magnetizations. It is possible to compute the partition function, $Z$, by utilizing the aforementioned information
|
||||
\begin{equation}
|
||||
Z = 2 e^{8J\beta} + 2 e^{-8J\beta} + 12,
|
||||
\end{equation}
|
||||
where $\beta = 1/(k_B T)$. From the partition function, we can derive the expectation values for the energy per spin state as
|
||||
\begin{equation}
|
||||
\begin{aligned}
|
||||
\expect{\epsilon} &= -\frac{J}{ZN} \sum_{i} E_i \cdot \text{deg}_i \cdot e^{-\beta E_i}\\
|
||||
&= \frac{6}{Z} \left(e^{-8\beta} - e^{8 \beta}\right),
|
||||
\end{aligned}
|
||||
\end{equation}
|
||||
with $N=4$. In general an expectation value of the ensemble of a general quantity $A$ can be expressed as
|
||||
\begin{equation}
|
||||
\expect{A} = \frac{1}{ZN} \sum_i A_i \cdot \text{deg}_i \cdot e^{-\beta E_i}.
|
||||
\end{equation}
|
||||
The following quantities are furthermore obtained analytically:
|
||||
\begin{align}
|
||||
\expect{\epsilon^2} &= \frac{72}{Z} \left(e^{-8\beta} + e^{8 \beta}\right)\\
|
||||
\expect{|m|} &= \frac{2}{Z} \left(2+e^{8\beta}\right)\\
|
||||
\expect{m^2} &= \frac{8}{Z} \left(1 + e^{8\beta}\right).
|
||||
\end{align}
|
||||
Using those expression it is also straightforward to obtain $\expect{\frac{C_V}{N}}$ and $\expect{\frac{\chi}{N}}$.
|
||||
|
||||
For the numerical evaluation at $k_BT = J = 1$ of all expectation values, \texttt{Python} is used. The expectation values shown in \cref{tab:analytical_expectation_values} are obtained from the numerical analysis.
|
||||
|
||||
\begin{table}
|
||||
\centering
|
||||
\caption{Numerical expectation values for the intensive thermodynamic properties of the 2x2 lattice.}
|
||||
\label{tab:analytical_expectation_values}
|
||||
\begin{tabular}{cc}
|
||||
\toprule
|
||||
Quantity & Numerical Expectation Value \\
|
||||
\midrule
|
||||
$\expect{\epsilon}$ & \num{-1.996} \\
|
||||
$\expect{|m|}$ & \num{0.999} \\
|
||||
$\expect{\frac{C_V}{N}}$ & \num{11.984} \\
|
||||
$\expect{\frac{\chi}{N}}$ & \num{2.996} \\ \bottomrule
|
||||
\end{tabular}
|
||||
\end{table}
|
||||
|
||||
% ===========================================
|
||||
\onecolumngrid
|
||||
% \bibliographystyle{apalike}
|
||||
\bibliographystyle{unsrt}
|
||||
\bibliography{ref}
|
||||
|
||||
|
||||
\end{document}
|
||||
@@ -0,0 +1,4 @@
|
||||
|
||||
real 0m29.580s
|
||||
user 5m39.713s
|
||||
sys 0m0.236s
|
||||
@@ -1,7 +1,6 @@
|
||||
#include "scheduling.hpp"
|
||||
|
||||
int main() {
|
||||
// Just a dummy run
|
||||
Scheduler scheduler(12, "results", true);
|
||||
scheduler.setMasterSeed(2025);
|
||||
scheduler.addRun(100, 2.3, 1000, true, 500);
|
||||
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -0,0 +1,4 @@
|
||||
|
||||
real 3m37.634s
|
||||
user 3m35.592s
|
||||
sys 0m0.207s
|
||||
@@ -3,6 +3,7 @@
|
||||
#include <cmath>
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
#include "ising.hpp"
|
||||
|
||||
|
||||
@@ -161,14 +162,19 @@ void IsingModel::saveResults(const std::string& filename) const {
|
||||
throw std::runtime_error("Could not open file for writing: " + filename);
|
||||
}
|
||||
|
||||
outFile << "Lattice Size: " << latticeSize << "\n";
|
||||
outFile << "Temperature: " << temp << "\n";
|
||||
outFile << "Number of Steps: " << n_steps << "\n";
|
||||
outFile << "Average Energy: " << calculateEnergy() << "\n";
|
||||
outFile << "Average Magnetization: " << calculateMagnetization() << "\n";
|
||||
outFile << "Heat Capacity: " << calculateHeatCapacity() << "\n";
|
||||
outFile << "Susceptibility: " << calculateSusceptibility() << "\n";
|
||||
double average_energy = calculateEnergy();
|
||||
double average_magnetization = calculateMagnetization();
|
||||
double heat_capacity = calculateHeatCapacity();
|
||||
double susceptibility = calculateSusceptibility();
|
||||
|
||||
outFile << std::fixed << std::setprecision(15)
|
||||
<< "Lattice Size: " << latticeSize << "\n"
|
||||
<< "Temperature: " << temp << "\n"
|
||||
<< "Number of Steps: " << n_steps << "\n"
|
||||
<< "Average Energy: " << average_energy << "\n"
|
||||
<< "Average Magnetization: " << average_magnetization << "\n"
|
||||
<< "Heat Capacity: " << heat_capacity << "\n"
|
||||
<< "Susceptibility: " << susceptibility << "\n";
|
||||
outFile.close();
|
||||
}
|
||||
|
||||
@@ -188,7 +194,7 @@ void IsingModel::saveEvolution(const std::string& filename) const {
|
||||
|
||||
outFile << "Step,CurrentEnergy,CurrentMagnetization,AvgEnergy,AvgMagnetization\n";
|
||||
for (const auto& record : evolutionData) {
|
||||
outFile << record[0] << "," << record[1] << "," << record[2] << "," << record[3] << "," << record[4] << "\n";
|
||||
outFile << std::fixed << std::setprecision(15) << record[0] << "," << record[1] << "," << record[2] << "," << record[3] << "," << record[4] << "\n";
|
||||
}
|
||||
|
||||
outFile.close();
|
||||
|
||||
@@ -2,6 +2,7 @@
|
||||
#include "omp.h"
|
||||
#include <iostream>
|
||||
#include <filesystem>
|
||||
#include <iomanip>
|
||||
|
||||
Scheduler::Scheduler(int numThreads, std::string outPutDir, bool recordEvolution) : numThreads(numThreads), outPutDir(outPutDir), recordEvolution(recordEvolution) {
|
||||
// Create output directory if it doesn't exist
|
||||
@@ -33,7 +34,8 @@ void Scheduler::executeRun(const RunConfig& config) {
|
||||
// Output results (could be saved to file or processed further)
|
||||
#pragma omp critical
|
||||
{
|
||||
std::cout << "Run (Size: " << config.size
|
||||
std::cout << std::fixed << std::setprecision(8)
|
||||
<< "Run (Size: " << config.size
|
||||
<< ", Temp: " << config.temperature
|
||||
<< ", Cycles: " << config.cycles
|
||||
<< ", Seed: " << config.seed << ") -> "
|
||||
|
||||
Reference in New Issue
Block a user