diff --git a/report/figures/fig_1.pdf b/report/figures/fig_1.pdf deleted file mode 100644 index 629164b..0000000 Binary files a/report/figures/fig_1.pdf and /dev/null differ diff --git a/report/figures/lattice1_1.png b/report/figures/lattice1_1.png new file mode 100644 index 0000000..6d8f309 Binary files /dev/null and b/report/figures/lattice1_1.png differ diff --git a/report/figures/lattice1_2.png b/report/figures/lattice1_2.png new file mode 100644 index 0000000..302a5a2 Binary files /dev/null and b/report/figures/lattice1_2.png differ diff --git a/report/figures/lattice1_3.png b/report/figures/lattice1_3.png new file mode 100644 index 0000000..5ce6f52 Binary files /dev/null and b/report/figures/lattice1_3.png differ diff --git a/report/figures/m_v_T.png b/report/figures/m_v_T.png new file mode 100644 index 0000000..216da2e Binary files /dev/null and b/report/figures/m_v_T.png differ diff --git a/report/figures/m_v_t_beta_in.png b/report/figures/m_v_t_beta_in.png new file mode 100644 index 0000000..fa638f4 Binary files /dev/null and b/report/figures/m_v_t_beta_in.png differ diff --git a/report/figures/m_v_t_beta_out.png b/report/figures/m_v_t_beta_out.png new file mode 100644 index 0000000..ad578d3 Binary files /dev/null and b/report/figures/m_v_t_beta_out.png differ diff --git a/report/figures/m_v_t_per.png b/report/figures/m_v_t_per.png new file mode 100644 index 0000000..fad0c01 Binary files /dev/null and b/report/figures/m_v_t_per.png differ diff --git a/report/ising_report.tex b/report/ising_report.tex new file mode 100644 index 0000000..d67c921 --- /dev/null +++ b/report/ising_report.tex @@ -0,0 +1,150 @@ +\documentclass[]{article} + +\usepackage{amsmath} % AMS math package +\usepackage{amssymb} % AMS symbol package +\usepackage{bm} % bold math +\usepackage{graphicx} % Include figure files +\usepackage{dcolumn} % Align table columns on decimal point +\usepackage{multirow} % Multirow/column tables +\usepackage{hyperref} % Hyperlinks +\usepackage{caption} +\usepackage{subcaption} + +\begin{document} + +\title{Ising Model} +\author{Justin Browne} +\date{} +\maketitle + +\begin{abstract} +This work covers the development and application of a monte carlo program to simulate the Ising Model. The monte carlo method used is the Metropolis algorithm. The program is used to study how a lattice of spins reacts to different external fields, initial lattice configurations, temperatures, and boundary conditons. +\end{abstract} + +\section{Introduction} +The Ising Model is a model used to describe magnetism in a lattice in which each lattice site can have a spin $\sigma_{i} = \pm 1$. +The interaction is described by the Hamiltonian +\[ +H = - \sum_{} J_{ij} \sigma_i \sigma_j -\mu \sum_{i} h_i \sigma_i ~, +\] +where $J_{ij}$ is the interaction between lattice sites $i$ and $j$, $\mu$ is the magnetic moment, $h_{i}$ is the magnetic field at lattice site $i$, and $\sum_{i~j}$ indicates a sum among only nearest neighbors. + +This model is a very important model in statistical mechanics. It is a simple model, and was the first with a phase transition to be solved exactly \cite{mccoy2012}. + +\section{Implementation} +The implementation of the Ising model in this work has incorporated a few simplifications. +First, the field is assumed to be equal everywhere ($h_{i} = h$). +Second, the interaction between spins are assumed to be the same ($J_{ij} = J$). This means the Hamiltonian is +\[ +H = - J \sum_{} \sigma_i \sigma_j -\mu h \sum_{i} \sigma_i ~. +\] + +The system is simulated using the Metropolis algorithm\cite{metropolis1953} as follows: +\begin{enumerate} + \item Choose an initial lattice configuration. + \item Choose a random point, $P$, on the lattice. + \item If flipping the spin at $P$ is energetically favorable, flip the spin. + If it is not, flip the spin with a probability of $e^{-\beta \Delta \! H}$, where $\Delta \! H$ is the difference between the Hamiltonian before and after the proposed spin flip. + \[ + \sigma_{P,n+1} = + \begin{cases} + -\sigma_{P,n} & \text{if } \Delta \! H < 0 \\ + -\sigma_{P,n} & \text{if } \Delta \! H > 0 \wedge r < e^{- \beta \Delta \! H} \\ + \phantom{-}\sigma_{P,n} & \text{otherwise} + \end{cases} + \] + where $\beta = \frac{1}{k_{B} T}$ is the inverse temperature, and $r$ is a random number from a uniform distriution spanning $[0,1)$. + \item Repeat from step 2. +\end{enumerate} + +All simulations in this work are performed using a $50 \times 50$ grid, with $J=1$, unless otherwise specified. + +\section{System Convergence} +The system converged at rates that varied depending on many variables, such as the initial configuration, the external field, and the temperature. + +The temperature affects the rate of convergence because the temperature determines the probability of a spin flip, which would not be energetically favorable, occuring. This means the system is more likely to escape from a local minimum. Figure~\ref{fig:conv_beta} shows how convergence varies with $\beta$. + +\begin{figure}[htb] + \centering + \begin{subfigure}[b]{.9\textwidth} + \includegraphics[width=\textwidth]{figures/m_v_t_beta_out.png} + \end{subfigure}\\ + \begin{subfigure}[b]{.9\textwidth} + \includegraphics[width=\textwidth]{figures/m_v_t_beta_in.png} + \end{subfigure} + \caption{The magnetization versus iteration number for various temperatures, zoomed out (top) and in (bottom). The higher temperature (lower $\beta$) runs converge faster. For these simulaltions, the field interaction was $\mu h = +1$, and all of the spins were initially aligned down.} + \label{fig:conv_beta} +\end{figure} + +Whether the boundaries are periodic or not affects the convergence as well. +Since the sites on the edge of the lattice interact with fewer sites, their spins are more likely to flip. +This means they act as nucleation sites. +Figure~\ref{fig:nucleation} shows the lattice of a non-periodic lattice at various points through the simulation, and figure~\ref{fig:conv_per} shows how the convergence depends on the boundary conditions. + +\begin{figure}[htb] + \centering + \begin{subfigure}[b]{.3\textwidth} + \includegraphics[width=\textwidth]{figures/lattice1_1.png} + \end{subfigure} + \begin{subfigure}[b]{.3\textwidth} + \includegraphics[width=\textwidth]{figures/lattice1_2.png} + \end{subfigure} + \begin{subfigure}[b]{.3\textwidth} + \includegraphics[width=\textwidth]{figures/lattice1_3.png} + \end{subfigure} + \caption{The lattice of a simulation with $\beta = 1$, $\mu h = +1$, at the beginning of (left), midway through (center), and at the end of (right) of the simulation. The lattice boundaries are non-periodic. Spin flips tend to propagate from the edges toward the center.} + \label{fig:nucleation} +\end{figure} + +\begin{figure}[htb] + \centering + \includegraphics[width=\textwidth]{figures/m_v_t_per.png} + \caption{The magnetization versus iteration number for periodic and non-periodic boundaries. The simulation with non-periodic boundaries more quickly escaped the local minimum.} + \label{fig:conv_per} +\end{figure} + +\section{Temperature Dependence} +Varying the temperature of the system shows some interesting properties. +At low temperatures, the system prefers a configuration in which all the spins are aligned. +At high temperatures, the system tends to prefer a randomized configuration. +The temperature around which the system undergoes this change is called the Curie Temperature. +The system behavior around the Curie Temperature has been solved analytically \cite{osanger1944,montroll1963}. +The solution is of the form +\[ + M = (1-\sinh^{4}(2 \beta J))^{1/8} ~. +\] +This behavior can be seen in figure~\ref{fig:critical}. + +\begin{figure}[htb] + \centering + \includegraphics[width=\textwidth]{figures/m_v_T.png} + \caption{The magnetization versus temperature. The simulation results and the analytical solution are show. For these simulations, $\mu h = 0$, the spins were initially aligned up, and the boundaries were periodic.} + \label{fig:critical} +\end{figure} + +Above the Curie temperature, all conditions result in a random configuration. +However, the system experiences bifurcation below the Curie Temperature. +Depending on the initial configuration and the external field, the system will settle to different stable configurations. +With an external field, the spins tend to align with the field. +However, without an external field, the spins tend to form domains of the same spin. +This means that low temperature systems with no external field that are initialized randomly act randomly and chaotically. +Different initial conditions and different random seeding greatly affect the outcome of the simulation. + +\section{Further Study} +There are more aspects of the Ising model that could be explored. +This work only focused on systems with $J=1$, corresponds to ferromagnetic systems. +Studying antiferromagnetic systems could prove to be very interesting. +The scope of the program could be expanded to include non-uniform magnetic fields, non-uniform or long-range spin-spin interactions. + +\begin{thebibliography}{} + \bibitem{mccoy2012} + McCoy, B.~M., \& Maillard, J.\ 2012, arXiv:1203.1456 + \bibitem{metropolis1953} + Metropolis, N., Rosenbluth, A.~W., Rosenbluth, M.~N., Teller, A.~H., \& Teller, E.\ 1953, Journal of Chemical Physics, 21, 1087 + \bibitem{osanger1944} + Onsager, Lars (1944), "Crystal statistics. I. A two-dimensional model with an order-disorder transition", Phys. Rev., Series II 65: 117–149 + \bibitem {montroll1963} + Montroll, Elliott W. and Potts, Renfrey B. and Ward, John C., Journal of Mathematical Physics, 4, 308-322 (1963) +\end{thebibliography} + +\end{document} diff --git a/report/report.tex b/report/report.tex deleted file mode 100644 index 72d1253..0000000 --- a/report/report.tex +++ /dev/null @@ -1,156 +0,0 @@ -\documentclass[]{article} - -\usepackage{amsmath} % AMS math package -\usepackage{amssymb} % AMS symbol package -\usepackage{bm} % bold math -\usepackage{graphicx} % Include figure files -\usepackage{dcolumn} % Align table columns on decimal point -\usepackage{multirow} % Multirow/column tables -\usepackage{hyperref} % Hyperlinks - -\begin{document} - -\title{Manuscript Title:\\with Forced Linebreak}% Force line breaks with \\ -\author{Ann Author} -\date{\today}% It is always \today, today, but you can specify other dates manually -\maketitle - -\begin{abstract} -An article usually includes an abstract, a concise summary of the work -covered at length in the main body of the article. -\end{abstract} - - -\section{First-level heading} %Title for the section -\label{sec:level1} %Label for the section, to be used for referencing in other parts of the document - -This sample document demonstrates some common uses of \LaTeX\ in documents you'll write for ICCP. Further information can be found online at the \href{http://en.wikibooks.org/wiki/LaTeX}{\LaTeX\ wikibook} or in the distribution documentation. - -When we refer to commands in this example file, they always have their required formal arguments in normal \TeX{} format. In this format, \verb+#1+, \verb+#2+, etc. stand for required author-supplied arguments to commands. For example, in \verb+\section{#1}+ the \verb+#1+ stands for the title text of the author's section heading, and in \verb+\title{#1}+ the \verb+#1+ stands for the title text of the paper. - -Line breaks in section headings at all levels can be introduced using \textbackslash\textbackslash. A blank input line tells \TeX\ that the paragraph has ended. - -\subsection{\label{sec:level2} Second-level heading: Formatting} - -This file may be formatted in either the \texttt{preprint} or \texttt{reprint} style. \texttt{reprint} format mimics final journal output. The paper size may be specified with the option \texttt{letter}. These options are inserted in the square brackets inside the \texttt{\textbackslash documentclass[]\{article\}} command. - -\section{Math and Equations} -Inline math may be typeset using the \verb+$+ delimiters. Bold math symbols may be achieved using the \verb+bm+ package and the \verb+\bm{#1}+ command it supplies. For instance, a bold $\alpha$ can be typeset as \verb+$\bm{\alpha}$+ giving $\bm{\alpha}$. Fraktur and Blackboard (or open face or double struck) characters should be typeset using the \verb+\mathfrak{#1}+ and \verb+\mathbb{#1}+ commands respectively. Both are supplied by the \texttt{amssymb} package. For -example, \verb+$\mathbb{R}$+ gives $\mathbb{R}$ and \verb+$\mathfrak{G}$+ gives $\mathfrak{G}$ - -In \LaTeX\ there are many different ways to display equations, and a few preferred ways are noted below. Displayed math will center by default. Use the class option \verb+fleqn+ to flush equations left. - -Below we have numbered single-line equations; this is generally the most common type of equation in \LaTeX: -\begin{equation} -\label{eq:one} %Label for the equation, to be used for referencing in other parts of the document - i \hbar \frac{\partial \Psi}{\partial t} = -\frac{\hbar^2}{2m}\nabla^2 \psi + U(\mathbf{x}) \psi -\end{equation} - -When the \verb+\label{#1}+ command is used [cf. input for Eq.~(\ref{eq:one})], the equation can be referred to in text without knowing the equation number that \TeX\ will assign to it. Just use \verb+\ref{#1}+, where \verb+#1+ is the same name that was used in the \verb+\label{#1}+ command. - -Unnumbered single-line equations can be typeset using the \verb+equation*+ environment: -\begin{equation*} - \hat{f}(\xi) = \int_{-\infty}^\infty f(x) e^{-2\pi i x \xi} \, \mathrm{d}x -\end{equation*} - -\subsection{Multiline equations} - -You'll generally want to use the \verb+align+ environment for multiline equations. Use the \verb+\nonumber+ command at the end of each line (again denoted with \textbackslash\textbackslash) to avoid assigning a number, or \verb+align*+ to disable numbering entirely: -\begin{align} - \sin(2\theta) &= 2 \sin(\theta) \cos(\theta) \nonumber \\ - &= \frac{2\tan(\theta)}{1 + \tan^2(\theta)} -\end{align} -Note how alignment occurs at the (unprinted) \verb+&+ character. - -Do not use \verb+\label{#1}+ on a line of a multiline equation if \verb+\nonumber+ is also used on that line. Incorrect cross-referencing will result. Notice the use \verb+\text{#1}+ for using a Roman font within a math environment. - -\section{Cross-referencing} -\LaTeX\ will automatically number such things as sections, footnotes, equations, figure captions, and table captions for you. In order to reference them in text, use the \verb+\label{#1}+ and \verb+\ref{#1}+ commands\footnote{Check out the \texttt{cleveref} (one r) package, too!}. To reference a particular page, use the \verb+\pageref{#1}+ command. - -The \verb+\label{#1}+ should appear within the section heading (or immediately after), within the footnote text, within the equation, or within the table or figure caption. The \verb+\ref{#1}+ command is used in text at the point where the reference is to be displayed. Some examples: Section~\ref{sec:level1} on page~\pageref{sec:level1}, Table~\ref{tab:table1}, and Fig.~\ref{fig:epsart}. -\begin{figure}[b] - \centering - \includegraphics{figures/fig_1}% Imports a figure - does not automatically scale - \caption{\label{fig:epsart} A figure caption. The figure captions are automatically numbered.} -\end{figure} - -\section{Floats: Figures, Tables, etc.} -Figures (images) and tables are usually allowed to ``float'', which means that their placement is determined by \LaTeX, while the document is being typeset. - -Use the \texttt{figure} environment for a figure, the \texttt{table} environment for a table. In each case, use the \verb+\caption+ command within to give the text of the figure or table caption along with the \verb+\label+ command to provide a key for referring to this figure or table. Insert an image using either the \texttt{graphics} or \texttt{graphix} packages, which define the \verb+\includegraphics{#1}+ command. (The two packages differ in respect of the optional arguments used to specify the orientation, scaling, and translation of the image.) To create an alignment, use the \texttt{tabular} environment. - -\begin{table} - \centering - \begin{tabular}{ |l|l|l| } - \hline - \multicolumn{3}{ |c| }{Team sheet} \\ - \hline - Goalkeeper & GK & Paul Robinson \\ \hline - \multirow{4}{*}{Defenders} & LB & Lucus Radebe \\ - & DC & Michael Duburry \\ - & DC & Dominic Matteo \\ - & RB & Didier Domi \\ \hline - \multirow{3}{*}{Midfielders} & MC & David Batty \\ - & MC & Eirik Bakke \\ - & MC & Jody Morris \\ \hline - Forward & FW & Jamie McMaster \\ \hline - \multirow{2}{*}{Strikers} & ST & Alan Smith \\ - & ST & Mark Viduka \\ - \hline - \end{tabular} - \caption{\label{tab:table1}A table, demonstrating the use of the \texttt{multirow} package for spanning rows and columns.} -\end{table} - -The best place for the \texttt{figure} or \texttt{table} environment is immediately following its first reference in text; this sample document illustrates this practice for Fig.~\ref{fig:epsart}, which shows a figure that is small enough to fit in a single column. In exceptional cases, you will need to move the float earlier in the document: \LaTeX's float placement algorithms need to know about a full-page-width float sooner. - -The content of a table is typically a \texttt{tabular} environment, giving rows of type in aligned columns. Column entries separated by \verb+&+'s, and -each row ends with \textbackslash\textbackslash. The required argument for the \texttt{tabular} environment specifies how data are aligned in the columns. -For instance, entries may be centered, left-justified, right-justified, aligned on a decimal point. - -Extra column-spacing may be be specified as well, although \LaTeX sets this spacing so that the columns fill the width of the table. Horizontal rules are typeset using the \verb+\hline+ command. Rows with that columns span multiple columns can be typeset using the \verb+\multicolumn{#1}{#2}{#3}+ command (for example, see the first row of Table~\ref{tab:table1}). - -\appendix - -\section{Appendixes} - -To start the appendixes, use the \verb+\appendix+ command. This signals that all following section commands refer to appendixes instead of regular sections. Therefore, the \verb+\appendix+ command should be used only once---to setup the section commands to act as appendixes. Thereafter normal section commands are used. The heading for a section can be left empty. For example, -\begin{verbatim} -\appendix -\section{} -\end{verbatim} -will produce an appendix heading that says ``APPENDIX A'' and -\begin{verbatim} -\appendix -\section{Background} -\end{verbatim} -will produce an appendix heading that says ``APPENDIX A: BACKGROUND'' (note that the colon is set automatically). - -If there is only one appendix, then the letter ``A'' should not appear. This is suppressed by using the star version of the appendix command (\verb+\appendix*+ in the place of \verb+\appendix+). - -\section{A little more on appendixes} - -Observe that this appendix was started by using -\begin{verbatim} -\section{A little more on appendixes} -\end{verbatim} - -Note the equation number in an appendix: -\begin{equation} - E^2=p^2c^2 + m^2c^4. -\end{equation} - -\subsection{\label{app:subsec}A subsection in an appendix} - -You can use a subsection or subsubsection in an appendix. Note the numbering: we are now in Appendix~\ref{app:subsec}. - -Note the equation numbers in this appendix, produced with the subequations environment: -\begin{subequations} -\begin{align} - E^2 &= m^2c^4 \quad \text{(rest)}, \label{appa} \\ - E^2 &= m^2c^4 + p^2 c^2 \quad \text{(relativistic)}, \label{appb} \\ - E^2 &\approx p^2 c^2 \quad \text{(ultrarelativistic)} \label{appc} -\end{align} -\end{subequations} -They turn out to be Eqs.~(\ref{appa}), (\ref{appb}), and (\ref{appc}). - -\end{document} diff --git a/src/Makefile b/src/Makefile index eb24188..0269d58 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,19 +1,25 @@ -FC = gfortran -FFLAGS = -Wall -Wextra -march=native -O3 -fimplicit-none +FC = gfortran +FFLAGS = -Wall -Wextra -march=native -O3 -fimplicit-none #FFLAGS += -pedantic -fbounds-check -fmax-errors=1 -g -FFLAGS += $(shell pkg-config --cflags plplotd-f95) -LDFLAGS = $(shell pkg-config --libs plplotd-f95) -LIBS = +FFLAGS += $(shell pkg-config --cflags plplotd-f95) +LDFLAGS2 = $(shell pkg-config --libs plplotd-f95) +LIBS = COMPILE = $(FC) $(FFLAGS) LINK = $(FC) $(LDFLAGS) -TARGET = ising.exe # Name of final executable to produce -OBJS = plot.o ising.o # List of object dependencies +TARGET1 = ising.exe # Name of final executable to produce +OBJS1 = rand_tools.o ising.o # List of object dependencies -$(TARGET): $(OBJS) +TARGET2 = post.exe # Name of final executable to produce +OBJS2 = plot.o post.o # List of object dependencies + +$(TARGET1): $(OBJS1) $(LINK) -o $@ $^ $(LIBS) +$(TARGET2): $(OBJS2) + $(LINK) $(LDFLAGS2) -o $@ $^ $(LIBS) + %.o:%.f95 $(COMPILE) -c $< diff --git a/src/TODO b/src/TODO new file mode 100644 index 0000000..2943b80 --- /dev/null +++ b/src/TODO @@ -0,0 +1,4 @@ +1. Look at making reading in options file more robust +1. Look at modularizing reading in option files (variable scope might be a problem) +1. Add more options for post-processing +1. plplot stuff is kinda kludged together, tidy it up diff --git a/src/ising.f95 b/src/ising.f95 new file mode 100644 index 0000000..5925f81 --- /dev/null +++ b/src/ising.f95 @@ -0,0 +1,203 @@ +program ising + use rand_tools + implicit none + + integer,allocatable :: sigma(:,:) + integer :: nx,ny + integer :: nsteps + real(8) :: inter + real(8) :: field + real(8) :: beta + integer :: iostatus + integer :: uopts,uout + character(64) :: a,b + character(64) :: foptsname,foutname + integer :: initmode + integer :: numout + logical :: periodic + + call init_random_seed() + numout = 0 + + uopts = 50 + uout = 52 + + ! Default options + nx = 10 + ny = 10 + nsteps = 0 + inter = 1 + field = 0 + beta = 1E0 + foutname = 'out' + initmode = 1 + periodic = .false. + + if(iargc()>0) then + call getarg(1,foptsname) + open(unit=uopts,file=foptsname,action='read') + ! Read in options + ! Format is: + ! + ! Input that doesn't make sense is ignored + ! TODO: Work on handling weird input files + do + ! Read in two strings + read(uopts,*,IOSTAT=iostatus) a,b + ! EOF ends input + if(iostatus<0) then + exit + ! Parse input + else + if(a == 'nx') read(b,*,IOSTAT=iostatus) nx + if(a == 'ny') read(b,*,IOSTAT=iostatus) ny + if(a == 'nsteps') read(b,*,IOSTAT=iostatus) nsteps + if(a == 'inter') read(b,*,IOSTAT=iostatus) inter + if(a == 'field') read(b,*,IOSTAT=iostatus) field + if(a == 'beta') read(b,*,IOSTAT=iostatus) beta + if(a == 'initmode') read(b,*,IOSTAT=iostatus) initmode + if(a == 'foutname') read(b,*,IOSTAT=iostatus) foutname + if(a == 'periodic') read(b,*,IOSTAT=iostatus) periodic + if(iostatus/=0) cycle + end if + end do + close(unit=uopts) + end if + + ! Output input + write(*,*) '# Simulation Parameters' + write(*,*) '# nx=', nx + write(*,*) '# ny=', ny + write(*,*) '# nsteps=', nsteps + write(*,*) '# inter=', inter + write(*,*) '# field=', field + write(*,*) '# beta=', beta + write(*,*) '# initmode=', initmode + write(*,*) '# foutname=', foutname + write(*,*) '# periodic=', periodic + + allocate(sigma(nx,ny)) + open(unit=uout,file=foutname,status='replace',form='unformatted',access='stream',action='write') + write(unit=uout) nx,ny,numout + + call initialize(initmode) + write(unit=uout) sigma + + call montecarlo() + + rewind(unit=uout) + write(unit=uout) nx,ny,numout + + close(unit=uout) + deallocate (sigma) + +contains + ! Initialize sigma + ! initmode - specifies the way in which we will initialize the array + ! 1: every spin is randomly assigned + ! 2: every spin is up + ! 3: every spin is down + ! default: 1 + subroutine initialize(initmode) + integer, intent(in) :: initmode + real, allocatable :: tmparr(:,:) + + if(initmode==2) then ! up + sigma = 1 + else if(initmode==3) then ! down + sigma = -1 + else ! random (default) + allocate(tmparr(nx,ny)) + call random_number(tmparr) + sigma = floor(tmparr*2)*2-1 + deallocate(tmparr) + end if + end subroutine initialize + + ! Flip a random spin + ! If that was energetically favorable, keep it + ! If not, keep it depending on temperature and energy change + ! Repeat for nsteps + subroutine montecarlo() + integer :: i,j,k + integer,allocatable :: newsigma(:,:) + real :: tmp + real(8) :: Ediff + + + ! Repeat until a set number of states have been tried (or criteria met) + allocate(newsigma(nx,ny)) + do k=1,nsteps + ! Determine new state by flipping a random spin + call random_number(tmp) + i = ceiling(tmp*nx) + call random_number(tmp) + j = ceiling(tmp*ny) + + ! Determine whether to accept or reject new state + ! If the energy is less, stay + ! Else, stay with probability of exp(-beta*(H_1-H_0)) + newsigma = sigma + newsigma(i,j) = -sigma(i,j) + Ediff = neighb_contrib(i,j,newsigma) + field_contrib(i,j,newsigma) & + - neighb_contrib(i,j,sigma) - field_contrib(i,j,sigma) + call random_number(tmp) + if(tmp < exp(-beta*Ediff)) then + sigma = newsigma + numout=numout+1 + write(unit=uout) k,i,j + end if + end do + deallocate(newsigma) + + end subroutine montecarlo + + function neighb_contrib(i,j,s) result(contrib) + integer,intent(in) :: i,j + integer,intent(in) :: s(:,:) + real(8) :: contrib + + contrib = 0 + + ! Interaction between nearest neighbors + if(i /= 1) contrib = contrib - inter*s(i,j)*s(i-1,j) + if(j /= 1) contrib = contrib - inter*s(i,j)*s(i,j-1) + if(i /= nx) contrib = contrib - inter*s(i,j)*s(i+1,j) + if(j /= ny) contrib = contrib - inter*s(i,j)*s(i,j+1) + + ! Loop around the edges if using periodic boundaries + if(periodic) then + if(i == 1) contrib = contrib - inter*s(i,j)*s(nx,j) + if(j == 1) contrib = contrib - inter*s(i,j)*s(i,ny) + if(i == nx) contrib = contrib - inter*s(i,j)*s(1,j) + if(j == ny) contrib = contrib - inter*s(i,j)*s(i,1) + end if + + end function neighb_contrib + + function field_contrib(i,j,s) result(contrib) + integer,intent(in) :: i,j + integer,intent(in) :: s(:,:) + real(8) :: contrib + + contrib = 0 + + ! Interaction with external field + contrib = contrib - field*s(i,j) + end function field_contrib + + function hamiltonian() result(energy) + integer :: i,j + real(8) :: energy + + energy = 0 + + do j=1,ny + do i=1,nx + ! neighb_contrib is halved to prevent double counting + energy = energy + 0.5*neighb_contrib(i,j,sigma) + energy = energy + field_contrib(i,j,sigma) + end do + end do + end function hamiltonian +end program ising diff --git a/src/plot.f95 b/src/plot.f95 index 1fa8c67..6ab5dfb 100644 --- a/src/plot.f95 +++ b/src/plot.f95 @@ -6,32 +6,41 @@ module plot public plot_init, plot_close, plot_lattice contains + subroutine plot_init(sizex,sizey) + integer, intent(in) :: sizex,sizey + character(len=30) :: geomstr + real(8) :: asp_ratio - subroutine plot_init(latticeSize) - integer, intent(in) :: latticeSize - call plsdev("xcairo") - call plinit() + write (geomstr,"(i0,'x',i0)") sizex,sizey + asp_ratio = dble(sizey)/dble(sizex) !You can find default colors at !http://plplot.sourceforge.net/docbook-manual/plplot-html-5.9.9/plcol0.html - !call plscol0(0, 255, 255, 255) ! white - !call plscol0(1, 255, 0, 0) ! red - !call plscol0(2, 255, 77, 0) ! orange - !call plscol0(3, 255, 255, 0) ! yellow - !call plscol0(4, 0, 255, 0) ! green - !call plscol0(5, 0, 0, 255) ! blue - !call plscol0(6, 0, 255, 255) ! cyan - !call plscol0(7, 255, 0, 255) ! magenta - !call plscol0(8, 128, 128, 128) ! gray - !call plscol0(9, 0, 0, 0) ! black + call plscol0(0, 255, 255, 255) ! white + call plscol0(1, 255, 0, 0) ! red + call plscol0(2, 255, 77, 0) ! orange + call plscol0(3, 255, 255, 0) ! yellow + call plscol0(4, 0, 255, 0) ! green + call plscol0(5, 0, 0, 255) ! blue + call plscol0(6, 0, 255, 255) ! cyan + call plscol0(7, 255, 0, 255) ! magenta + call plscol0(8, 128, 128, 128) ! grey + call plscol0(9, 0, 0, 0) ! black - call plenv(0d0, latticeSize + 1d0, 0d0, latticeSize + 1d0, 0, 0) + call plsetopt("geometry",geomstr) ! Change? + call plsdev("png") + call plsfam(1,1,100000) + call plsfnam("%n.png") + call plinit() + call pladv(0) + call plvpas(0d0, 1d0, 0d0, 1d0, asp_ratio) + call plwind(5d-1, sizex+5d-1, 5d-1, sizey+5d-1) end subroutine plot_init subroutine plot_close() call plspause(.false.) - call plend() + !call plend() end subroutine plot_close subroutine plot_lattice(lattice) @@ -49,11 +58,11 @@ subroutine plot_lattice(lattice) do j = 1, size(lattice, 2) x = i; y = j if (lattice(i, j) .eq. 1) then - call plcol0(1) !default red - call plpoin([x], [y], 30) !30 denotes up-arrow glyph + call plcol0(8) ! grey + call plot_square(x,y) else ! Comment out the else clauses to speed drawing - call plcol0(11) !default cyan - call plpoin([x], [y], 31) !31 denotes down-arrow glyph + call plcol0(9) ! black + call plot_square(x,y) end if end do end do @@ -61,4 +70,13 @@ subroutine plot_lattice(lattice) call plflush() end subroutine plot_lattice + subroutine plot_square(i,j) + real(8), intent(in) :: i,j + real(8) :: x(4),y(4) + + x = (/i-.5,i-.5,i+.5,i+.5/) + y = (/j-.5,j+.5,j+.5,j-.5/) + + call plfill(x,y) + end subroutine end module plot diff --git a/src/post.f95 b/src/post.f95 new file mode 100644 index 0000000..b2d88f0 --- /dev/null +++ b/src/post.f95 @@ -0,0 +1,106 @@ +program post + use plot + + integer :: i + integer :: nx,ny,nsteps + integer,allocatable :: sigma(:,:) + integer :: uopts,uin,umag + character(64) :: a,b + character(64) :: foptsname,finname,magoutname + integer :: iostatus + real(8),allocatable :: m(:) + integer :: step,x,y + + uopts = 50 + uin = 51 + umag = 52 + + ! Default options + finname = 'out' + magoutname = 'mag.out' + + if(iargc()>0) then + call getarg(1,foptsname) + open(unit=uopts,file=foptsname,action='read') + ! Read in options + ! Format is: + ! + ! Input that doesn't make sense is ignored + ! TODO: Work on handling weird input files + do + ! Read in two strings + read(uopts,*,IOSTAT=iostatus) a,b + ! EOF ends input + if(iostatus<0) then + exit + ! Parse input + else + if(a == 'finname') read(b,*,IOSTAT=iostatus) finname + if(a == 'magoutname') read(b,*,IOSTAT=iostatus) magoutname + if(iostatus/=0) cycle + end if + end do + close(unit=uopts) + end if + + open(unit=uin,file=finname,form='unformatted',access='stream',action='read') + open(unit=umag,file=magoutname,action='write') + read(unit=uin) nx,ny,nsteps + + allocate(sigma(nx,ny)) + allocate(m(nsteps+1)) + read(unit=uin) sigma + call plot_init(nx,ny) + + call plbop() + call plot_lattice(sigma) + call pleop() + + i=1 + m(i) = magnetization() + write(umag,*) 0,m(i) + do + i = i+1 + read(unit=uin,iostat=iostatus) step,x,y + + if(iostatus>0) then + write(0,*) 'There was a problem reading in the file' + call exit(iostatus) + else if(iostatus<0) then + exit + else + sigma(x,y) = -sigma(x,y) + m(i) = magnetization() + write(umag,*) step,m(i) + +! call plbop() +! call plot_lattice(sigma) +! call pleop() + end if + end do + + call plbop() + call plot_lattice(sigma) + call pleop() + + call plot_close() + deallocate(m) + deallocate(sigma) + close(unit=umag) + close(unit=uin) + +contains + function magnetization() result(mag) + integer :: i,j + integer :: m + real(8) :: mag + + m = 0 + do j=1,ny + do i=1,nx + m=m+sigma(i,j) + end do + end do + mag = dble(m)/(nx*ny) + end function magnetization +end program post diff --git a/src/rand_tools.f95 b/src/rand_tools.f95 new file mode 100644 index 0000000..db7ad64 --- /dev/null +++ b/src/rand_tools.f95 @@ -0,0 +1,47 @@ +module rand_tools +contains + subroutine init_random_seed() + implicit none + integer, allocatable :: seed(:) + integer :: i, n, un, istat, dt(8), pid, t(2), s + integer(8) :: count, tms + + call random_seed(size = n) + allocate(seed(n)) + open(newunit=un, file="/dev/urandom", access="stream",& + form="unformatted", action="read", status="old", & + iostat=istat) + if (istat == 0) then + read(un) seed + close(un) + else + call system_clock(count) + if (count /= 0) then + t = transfer(count, t) + else + call date_and_time(values=dt) + tms = (dt(1) - 1970)*365_8 * 24 * 60 * 60 * 1000 & + + dt(2) * 31_8 * 24 * 60 * 60 * 1000 & + + dt(3) * 24 * 60 * 60 * 60 * 1000 & + + dt(5) * 60 * 60 * 1000 & + + dt(6) * 60 * 1000 + dt(7) * 1000 & + + dt(8) + t = transfer(tms, t) + end if + s = ieor(t(1), t(2)) + pid = getpid() + 1099279 ! Add a prime + s = ieor(s, pid) + if (n >= 3) then + seed(1) = t(1) + 36269 + seed(2) = t(2) + 72551 + seed(3) = pid + if (n > 3) then + seed(4:) = s + 37 * (/ (i, i = 0, n - 4) /) + end if + else + seed = s + 37 * (/ (i, i = 0, n - 1 ) /) + end if + end if + call random_seed(put=seed) + end subroutine init_random_seed +end module rand_tools