diff --git a/report/.DS_Store b/report/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/report/.DS_Store differ diff --git a/report/figures/.DS_Store b/report/figures/.DS_Store new file mode 100644 index 0000000..5008ddf Binary files /dev/null and b/report/figures/.DS_Store differ diff --git a/report/figures/GammaValues.jpg b/report/figures/GammaValues.jpg new file mode 100644 index 0000000..4e5af4d Binary files /dev/null and b/report/figures/GammaValues.jpg differ diff --git a/report/figures/HeatCapacity.jpg b/report/figures/HeatCapacity.jpg new file mode 100644 index 0000000..307c87f Binary files /dev/null and b/report/figures/HeatCapacity.jpg differ diff --git a/report/figures/InternalEnergy.jpg b/report/figures/InternalEnergy.jpg new file mode 100644 index 0000000..b00990b Binary files /dev/null and b/report/figures/InternalEnergy.jpg differ diff --git a/report/figures/Magnetization.jpg b/report/figures/Magnetization.jpg new file mode 100644 index 0000000..1e5701f Binary files /dev/null and b/report/figures/Magnetization.jpg differ diff --git a/report/figures/Susceptibility.jpg b/report/figures/Susceptibility.jpg new file mode 100644 index 0000000..be8ba6b Binary files /dev/null and b/report/figures/Susceptibility.jpg differ diff --git a/report/figures/UTexasHeatCapacity.png b/report/figures/UTexasHeatCapacity.png new file mode 100644 index 0000000..d62e7bc Binary files /dev/null and b/report/figures/UTexasHeatCapacity.png differ diff --git a/report/figures/UTexasMagnetization.jpg b/report/figures/UTexasMagnetization.jpg new file mode 100644 index 0000000..084c732 Binary files /dev/null and b/report/figures/UTexasMagnetization.jpg differ diff --git a/report/monte-carlo-simulation.tex b/report/monte-carlo-simulation.tex new file mode 100644 index 0000000..c4bd2b8 --- /dev/null +++ b/report/monte-carlo-simulation.tex @@ -0,0 +1,209 @@ +\documentclass[10pt,letterpaper]{article} + +\usepackage[english]{babel} +\usepackage[utf8x]{inputenc} +\usepackage{amsmath} +\usepackage{graphicx} +\usepackage[colorinlistoftodos]{todonotes} + +\usepackage{cite} +\usepackage{graphicx} +\usepackage{float} + +\title{Monte Carlo Simulation of Ising Model Lattice} +\author{Stephen Mee} + +\begin{document} +\maketitle + +\begin{abstract} + +Using the Wolff algorithm for a range of temperature values, the critical temperature and other macroscopic and microscopic values are found for a two-dimensional lattice of electron spin states. + +\end{abstract} + +\section{Introduction} + +A two-dimensional lattice of electron spin states is initialized. For each temperature value, a random site is chosen and added to the cluster specified in the Wolff algorithm. Using periodic boundary conditions, the 4 nearest neighbors of the chosen site are added to a "queue" if they pass an exponential check. The sites that are in the queue are subsequently added to the cluster and the algorithm to add their nearest neighbors is run again. Sites that have been added to the cluster immediately have the sign of their spin flipped. + +After a predetermined number of clusters are built on the lattice, a range of values are calculated and plotted vs. temperature. These values include: magnetization, critical temperature, magnetic susceptibility, internal energy, heat capacity, and critical exponents of the magnetization near the critical temperature. + +\section{Wolff Algorithm} + +The Wolff algorithm relies on the idea of flipping entire "clusters" of spin states at once, as opposed to the Metropolis algorithm which only flips individual states. The condition to add sites to the cluster is as follows: + +\begin{equation} +R < 1 - \exp\frac{-2J}{k_{b}T} +\end{equation} + +Where $R$ is a number picked from a random distribution ranging from 0 to 1, $J$ is the coupling constant, $k_{B}$ is the Boltzmann constant, and $T$ is the temperature. In the code, both $J$ and $k_{B}$ have arbitrarily been set to one. + +The size of the cluster is a function of temperature; for low temperatures we expect the cluster to encompass the entire lattice while for high temperatures we expect the clusters to be much smaller. Because of this temperature dependence, the number of clusters being built on a single lattice also needs be a function of temperature to ensure accurate results. Each cluster will sample less at high temperatures, so more clusters will be needed to cover the entire lattice. A temperature-dependent cluster number has been built into the code to compensate for this. + + + +\section{Code Results} + +All results were found with a lattice of N = 22500 that was initialized with all spin-up sites. The temperature ranged from T = 0.01K - 4.0K in iterations of 0.01K. From 0.01K - 2.0K there are 10 clusters built per lattice, 2.01K - 2.5K there are $2*\sqrt[]{N}$, and from 2.51K - 4.0K there are $2N$. The cluster process is repeated 20 times for any given temperature to get an average magnetization and calculate the standard deviation. + +\subsection{Magnetization} + +The magnetization is found by: + +\begin{equation} +M = \frac{1}{N} \sum_{i}\sum_{j}{S_{i,j}} +\end{equation} + +Where $S_{i,j}$ is the spin state of the $[i,j]$ index of the electron lattice. The magnetization is expected to stay near one for low temperatures then plummet to zero near the "critical temperature", following a roughly S-shaped curve. + +\begin{figure}[H] +\centering +\includegraphics[width=120mm]{Magnetization.jpg} +\caption{Magnetization versus temperature for N = 22500} +\end{figure} + +\pagebreak + +\subsection{Magnetic Susceptibility} + +The susceptibility is defined as: + +\begin{equation} +\chi = \frac{\partial{M}}{\partial{T}} +\end{equation} + +Assuming many temperature iterations, this can approximated as: + +\begin{equation} +\chi = \frac{\Delta{M}}{\Delta{T}} = \frac{M_{n+1} - M_{n}}{T_{n+1} - T_{n}} +\end{equation} + +Based off the plot for magnetization one would expect the susceptibility to look like a sharp Gaussian in the negative direction, with the peak located at the point where the magnetization was most rapidly changing. + +\begin{figure}[H] +\centering +\includegraphics[width=120mm]{Susceptibility.jpg} +\caption{Magnetic Susceptibility vs Temperature for N = 22500} +\end{figure} + +\subsection{Internal Energy} + +The internal energy is defined as: + +\begin{equation} +E = - J \sum_{i,NN}S_{i}S_{NN} +\end{equation} + +Where $J$ is the coupling constant, $S_{i}$ is a single spin, and $S_{NN}$ is representative of the spins of the nearest neighbors of $S_{i}$. + +When summing over all indexes the algorithm must be careful to not double-count the same interaction. This was bypassed by only considering a single horizontal and vertical neighbor for any given particle. Assuming periodic boundary conditions, summing over all $S_{i}$ using this method manages to include all nearest neighbor interactions. + +\begin{figure}[H] +\centering +\includegraphics[width=120mm]{InternalEnergy.jpg} +\caption{Internal Energy vs Temperature for N = 22500} +\end{figure} + +\subsection{Heat Capacity} + +The heat capacity is defined as: + +\begin{equation} +C = \frac{\partial{E}}{\partial{T}} +\end{equation} + +Similarly to the susceptibility, this can be rewritten as: + +\begin{equation} +C = \frac{\Delta{E}}{\Delta{T}} = \frac{E_{n+1} - E_{n}}{T_{n+1} - T{n}} +\end{equation} + +Once again, the heat capacity would be expected to be a Gaussian in the positive direction with the peak at the moment of the most rapid change. + +\begin{figure}[H] +\centering +\includegraphics[width=120mm]{HeatCapacity.jpg} +\caption{Heat Capacity versus Temperature for N = 22500} +\end{figure} + +\subsection{Critical Temperature} + +The critical temperature can be defined as both the point at which the magnetization and the internal energy are most rapidly changing. Equivalently, we are looking for the temperature where the concavity of the magnetization and internal energy switches sign. This is written mathematically as: + +\begin{equation} +0 = \frac{\partial^{2}M}{\partial{T}^{2}} +\end{equation} + +\begin{equation} +0 = \frac{\partial^{2}E}{\partial{T}^{2}} +\end{equation} + +By solving for T in both cases, we can conclude: + +\begin{equation} +T_{c} \approx 2.21K +\end{equation} + +This is supported by the simulation. Without averaging the magnetization over many runs, 2.2K is approximately when the calculated values began to fluctuate wildly. The statistical noise after the simulation passes the critical temperature is still readily apparent in Figure 4. + +\subsection{Critical Exponents} + +It was decided that the calculation of the critical exponent $\gamma$ for the magnetization would yield the cleanest results. Around the critical temperature, the magnetization obeys the following power law: + +\begin{equation} +M = |T - T_{c}|^{\gamma} +\end{equation} + +Which can be rewritten as: + +\begin{equation} +\gamma = \frac{\ln(M)}{\ln(|T - T_{c}|)} = \ln(M - |T - T_{c}|) +\end{equation} + +The calculation of gamma requires a lot of precise data in the region surrounding the critical temperature. Because of this, it was decided the original run would not provide good enough statistics. A second simulation run was preformed with N = 22500, $T_{init}$ = 2.18K, $T_{end}$ = 2.24K, and $T_{iter}$ = 0.002K to provide better data. The gamma values calculated for each temperature iteration are as follows: + +\begin{figure}[H] +\centering +\includegraphics[width=120mm]{GammaValues.jpg} +\caption{A range of calculated gamma values around the critical temperature} +\end{figure} + +After the simulation passes the critical temperature the gamma values begin to fluctuate, so only temperature values up to 2.224K are plotted above. By taking the average of all the above values above we can estimate: + +\begin{equation} +\gamma = 0.119 \approx 0.125 +\end{equation} + +Which is the accepted value of the critical index. + +\section{Conclusion} + +The Wolff algorithm provides more concise and accurate data for low temperatures at the expense of accuracy at high temperatures. For the same effective input the Wolff algorithm also runs faster than the Metropolis algorithm, but suffers from being a much more specific model that cannot apply to many other cases. + +At high temperatures the clusters become very small (single sites), and the Wolff Algorithm devolves into the Metropolis algorithm. Because of this more clusters must be run at high temperatures for accurate. + +Overall this simulation produces results consistent with the Metropolis algorithm simulation and the theory behind the Ising model. The Ising model predicts a critical temperature $T_{C}$ around 2.27K, very close to the calculated 2.21K. Below are plots of magnetization and heat capacity calculated by Nikos Drakos from the University of Leeds and Ross Moore from Macquarie University. Both plots can be seen to match the results presented in this report. + +\begin{figure}[H] +\centering +\includegraphics[width=50mm]{UTexasMagnetization.jpg} +\caption{Drakos, Nikos, and Ross Moore. Figure 115. Digital image. University of Texas, 29 Mar. 2006. Web. 31 Mar. 2015. } +\end{figure} + +\begin{figure}[H] +\centering +\includegraphics[width=50mm]{UTexasHeatCapacity.png} +\caption{Drakos, Nikos, and Ross Moore. Figure 116. Digital image. University of Texas, 29 Mar. 2006. Web. 31 Mar. 2015. .} +\end{figure} + + + + + + + + + + + +\end{document} \ No newline at end of file diff --git a/src/.DS_Store b/src/.DS_Store new file mode 100644 index 0000000..a411092 Binary files /dev/null and b/src/.DS_Store differ diff --git a/src/a.out b/src/a.out new file mode 100755 index 0000000..6cb9209 Binary files /dev/null and b/src/a.out differ diff --git a/src/a.out.dSYM/Contents/Info.plist b/src/a.out.dSYM/Contents/Info.plist new file mode 100644 index 0000000..3679a65 --- /dev/null +++ b/src/a.out.dSYM/Contents/Info.plist @@ -0,0 +1,20 @@ + + + + + CFBundleDevelopmentRegion + English + CFBundleIdentifier + com.apple.xcode.dsym.a.out + CFBundleInfoDictionaryVersion + 6.0 + CFBundlePackageType + dSYM + CFBundleSignature + ???? + CFBundleShortVersionString + 1.0 + CFBundleVersion + 1 + + diff --git a/src/a.out.dSYM/Contents/Resources/DWARF/a.out b/src/a.out.dSYM/Contents/Resources/DWARF/a.out new file mode 100644 index 0000000..64b056a Binary files /dev/null and b/src/a.out.dSYM/Contents/Resources/DWARF/a.out differ diff --git a/src/main.cpp b/src/main.cpp new file mode 100644 index 0000000..e29531e --- /dev/null +++ b/src/main.cpp @@ -0,0 +1,386 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using std::cout; +using std::endl; +using std::vector; +using std::ofstream; +using std::queue; +using std::pair; +using std::make_pair; + +void addNearestNeighbors(vector< vector > &lattice, queue < pair > &index_queue, int x, int y, int dimension, double J, double kb, double T, std::default_random_engine generator){ + + std::uniform_real_distribution distribution(0.0,1.0); + generator.seed(rand()); + + int left_coor = x - 1; + int right_coor = x + 1; + int up_coor = y + 1; + int down_coor = y - 1; + + if (left_coor == -1){ + left_coor = dimension - 1; + } + if (right_coor == dimension){ + right_coor = 0; + } + if (up_coor == dimension){ + up_coor = 0; + } + if (down_coor == -1){ + down_coor = dimension - 1; + } + + int current = lattice[x][y]; + + int up_neighbor = lattice[x][up_coor]; + int down_neighbor = lattice[x][down_coor]; + int right_neighbor = lattice[right_coor][y]; + int left_neighbor = lattice[left_coor][y]; + + double R = distribution(generator); + if (R < (1.0 - exp(-2.0*J/kb/T)) && (current != left_neighbor)){ + + index_queue.push(make_pair(left_coor, y)); + lattice[left_coor][y] *= -1; + + } + + R = distribution(generator); + if (R < (1.0 - exp(-2.0*J/kb/T)) && (current != right_neighbor)){ + + index_queue.push(make_pair(right_coor,y)); + lattice[right_coor][y] *= -1; + + } + + R = distribution(generator); + if (R < (1.0 - exp(-2.0*J/kb/T)) && (current != up_neighbor)){ + + index_queue.push(make_pair(x, up_coor)); + lattice[x][up_coor] *= -1; + + } + + R = distribution(generator); + if (R < (1.0 - exp(-2.0*J/kb/T)) && (current != down_neighbor)){ + + index_queue.push(make_pair(x, down_coor)); + lattice[x][down_coor] *= -1; + + } + +} + +double calculateInternalEnergy(vector< vector > lattice, int dimension, double J){ + + double energy = 0; + + for (int i = 0; i < dimension; i++){ + for (int j = 0; j < dimension; j++){ + + int up_coor = j + 1; + int right_coor = i + 1; + + if (up_coor == dimension){ + up_coor = 0; + } + + if (right_coor == dimension){ + right_coor = 0; + } + + energy += lattice[i][up_coor] * lattice[right_coor][j]; + + } + } + + return -J * energy; + +} + +vector calculateHeatCapacity(vector interalEnergyList, vector temperatureList){ + + vector heatCapacityList; + int limit = interalEnergyList.size(); + limit = limit - 1; + + for (int i = 0; i < limit; i++){ + + double deltaE = interalEnergyList[i+1] - interalEnergyList[i]; + double deltaT = temperatureList[i+1] - temperatureList[i]; + + heatCapacityList.push_back(fabs(deltaE / deltaT)); + + } + + return heatCapacityList; + +} + +vector calculateSusceptibility(vector magnetizationList, vector temperatureList){ + + vector susceptibility; + int limit = magnetizationList.size(); + limit = limit - 1; + + for (int i = 0; i < limit; i++){ + + double deltaMag = magnetizationList[i+1] - magnetizationList[i]; + double deltaTemp = temperatureList[i+1] - temperatureList[i]; + + susceptibility.push_back(-1.0 * fabs(deltaMag / deltaTemp)); + } + + return susceptibility; +} + +double calculateCriticalTemperature(vector susceptibility, vector temperatureList){ + + double maxS = 0; + + double criticalTemp; + + int sLimit = susceptibility.size(); + + for (int i = 0; i < sLimit; i++){ + + if (fabs(susceptibility[i]) > maxS){ + + maxS = fabs(susceptibility[i]); + criticalTemp = temperatureList[i]; + + } + + } + + return criticalTemp; + +} + +void buildCluster(vector< vector > &lattice, int dimension, int x, int y, double J, double kb, double T, std::default_random_engine generator){ + + // cout << "In build cluster" << endl; + + std::uniform_real_distribution distribution(0.0,1.0); + generator.seed(rand()); + + queue < pair > index_queue; + lattice[x][y] *= -1; + + addNearestNeighbors(lattice, index_queue, x, y, dimension, J, kb, T, generator); + + while (!index_queue.empty()){ + + pair current_index = index_queue.front(); + index_queue.pop(); + int new_x = current_index.first; + int new_y = current_index.second; + + addNearestNeighbors(lattice, index_queue, new_x, new_y, dimension, J, kb, T, generator); + + } + + /*std::ostream_iterator out_it(std::cerr, " "); + for(auto row: lattice) { + std::copy(row.begin(), row.end(), out_it); + std::cerr << std::endl; + } + cout << endl;*/ + +} + +double findAverage(vector< vector > lattice, int dimension){ + + double sum = 0.0; + + for (int i = 0; i < dimension; i++){ + + for (int j = 0; j < dimension; j++){ + + sum += lattice[i][j]; + } + } + + return fabs(sum / (dimension * dimension)); +} + +double findStdDev(vector stdDevCalcList, int dimension){ + + double summation = 0.0; + double z = 0; + + for (int i = 0; i < dimension; i++){ + + z += stdDevCalcList[i]; + } + + double avg = z / dimension; + + for (int i = 0; i < dimension; i++){ + + summation += pow((avg - stdDevCalcList[i]), 2); + } + + double variance = summation / (dimension * dimension); + return sqrt(variance); +} + +int main () { + + srand(time(NULL)); + + int dimension = 100; + double start_T = 0.01; // Kelvin + double end_T = 4.0; // Kelvin + double iter_T = 0.01; // Kelvin + const double kb = 1.0; + const double J = 1.0; + + std::uniform_real_distribution distribution(0.0,1.0); + std::default_random_engine generator; + generator.seed(rand()); + + vector magnetizationList; + vector temperatureList; + vector stdDevCalcList; + vector stdDevList; + vector interalEnergyList; + + /*std::ostream_iterator out_it(std::cerr, " "); + for(auto row: lattice) { + std::copy(row.begin(), row.end(), out_it); + std::cerr << std::endl; + }*/ + + double magnetization; + + int clusterN = 10; + int stdDevRuns = 20; + + // Iterate through every temperature value + for (double T = start_T; T <= end_T; T += iter_T){ + + magnetization = 0; + vector< vector > lattice(dimension, vector (dimension,1)); + + // Uncomment below to randomize initial lattice + + /*for (int i = 0; i < dimension; i++){ + for (int j = 0; j < dimension; j++){ + + double rand_float = distribution(generator); + + if (rand_float < 0.5){ + + lattice[i][j] = -1; + + } else{ + + lattice[i][j] = 1; + + } + } + }*/ + + if (T > 2.0 && T <= 2.5){ + clusterN = 2 * dimension; + } + if (T > 2.5){ + clusterN = 2 * dimension * dimension; + } + + double totMag = 0.0; + + for (int k = 0; k < stdDevRuns; k++){ + + for (int i = 0; i <= clusterN; i++){ + + int i_rand = rand() % dimension; + int j_rand = rand() % dimension; + //int i_rand = 0; + //int j_rand = 0; + buildCluster(lattice, dimension, i_rand, j_rand, J, kb, T, generator); + + + } // End of Metropolis Test loop + + magnetization = findAverage(lattice, dimension); + totMag += magnetization; + stdDevCalcList.push_back(magnetization); + + } + + double avgMag = totMag / stdDevRuns; + magnetizationList.push_back(avgMag); + temperatureList.push_back(T); + + double stdDev = findStdDev(stdDevCalcList, stdDevCalcList.size()); + stdDevList.push_back(stdDev); + + double internalEnergy = calculateInternalEnergy(lattice, dimension, J); + interalEnergyList.push_back(internalEnergy); + + cout << "For temperature " << T << ", magnetization: " << avgMag << ", Std Dev: " << stdDev << endl; + + // output << T << "\t" << "\t" << avgMag << "\t" << "\t" << stdDev << endl; + + + } // End of temperature iteration loop + + ofstream output; + output.open("Output.txt"); + output << "Temperature" << "\t" << "Magnetization" << "\t" << "Standard Deviation" << "\t" << "Internal Energy" << "\t" << "Susceptibility" << "\t" << "Heat Capacity" << endl; + + vector susceptibility = calculateSusceptibility(magnetizationList, temperatureList); + vector heatCapacityList = calculateHeatCapacity(interalEnergyList, temperatureList); + + int arraySize = magnetizationList.size(); + + for (int i = 0; i < arraySize; i++){ + + if (i == arraySize - 1){ + + output << temperatureList[i] << "\t" << "\t" << magnetizationList[i] << "\t" << "\t" << stdDevList[i] << "\t" << "\t" << interalEnergyList[i] << endl; + + }else{ + + output << temperatureList[i] << "\t" << "\t" << magnetizationList[i] << "\t" << "\t" << stdDevList[i] << "\t" << "\t" << interalEnergyList[i] << "\t" << "\t" << susceptibility[i] << "\t" << "\t" << heatCapacityList[i] << endl; + + } + } + + double criticalTemp = calculateCriticalTemperature(susceptibility, temperatureList); + + // Gamma calculation loop below + + double gamma; + int gammaLimit = magnetizationList.size(); + double Tc = 2.21; + for (int i = 0; i < gammaLimit; i++){ + + double tempDiff = fabs(temperatureList[i] - Tc); + gamma = log(magnetizationList[i]) / log(tempDiff); + + cout << gamma << endl; + + } + + cout << "Critical Temperature: " << criticalTemp << endl; + + output.close(); + cout << "Completed." << endl; + cout << "Output.txt file created" << endl; + return 0; + +} \ No newline at end of file