diff --git a/ising.f b/ising.f new file mode 100755 index 0000000..469d053 --- /dev/null +++ b/ising.f @@ -0,0 +1,173 @@ + PROGRAM ISING + ! 2D Monte Carlo simulation of Ising Model with Metropolis algorithm + ! Program written by Tridip Das (dastridi) + ! This program calculates magnetization with temperature variation + ! Susceptibility, Energy and Specific heat with temperature + ! Also finds the critical temperature at which magnetization is + ! lost. Scaled temperature is used, Boltzmann const, kB assumed as 1 + ! + IMPLICIT NONE + ! Variable Declaration + INTEGER :: I, J, M, N ! counters + INTEGER, ALLOCATABLE :: A(:,:),AOLD(:,:),ANEW(:,:) ! 2D array containing spins + INTEGER :: NROWS, NCOLS ! Defining # of rows and columns + ! + REAL :: TEMPERATURE, BETA, MAGNETIZATION, ENERGY + REAL :: MAG2_AVG, ENERGY2_AVG, SUSP, Cv + REAL :: DELTA_E, MAG_AVG, RAND_NUM, ENERGY_AVG + INTEGER :: TEMP_ITER, ITER_AVG, MC_ITER, TEMP_VAR + INTEGER :: NO_OF_MC_STEPS, NO_OF_ITER, TEMP_STEPS + REAL :: T_INITIAL, T_FINAL, TEMP_INTERVAL, START, FINISH + ! + ! + PRINT *, " MONTE CARLO SIMULATION OF 2D ISING MODEL " + PRINT *, " SIMULATION WITH METROPOLIS ALGORITHM " + PRINT *, " APPLYING PERIODIC BOUNDARY CONDITION " + PRINT *, "--------Code developed by Tridip Das------" + CALL CPU_TIME(START) + PRINT *, " Start time ", START + PRINT *, "------------------------------------------" + ! INPUT PARAMETERS FROM inp_file + ! OPEN(UNIT=11,FILE='inpfile',STATUS='OLD',ACTION='READ') + ! Uncomment below lines to get input from user + ! PRINT *, " PROVIDE THE NUMBER OF ROWS FOR ISING MODEL" + ! READ(*,*) NROWS + ! PRINT *, " PROVIDE THE NUMBER OF COLS FOR ISING MODEL" + ! READ(*,*) NCOLS + ! PRINT *, " NUMBER OF SIMULATION STEPS TO RUN FOR MC & +! & EQULIBRIATION " + ! READ(*,*) NO_OF_MC_STEPS + ! PRINT *, " NUMBER OF ITERATIONS TO PERFORM AT EACH TEMPERATURE & +! & STEP TO CALCULATE AVERAGE OVER MC PREDICTION" + ! READ(*,*) NO_OF_ITER + ! PRINT *, " PROVIDE INITIAL, FINAL TEMPERATURE AND INTERVAL" + ! PRINT *, " SEPARATED BY ENTER " + ! READ(*,*) T_INITIAL + ! READ(*,*) T_FINAL + ! READ(*,*) TEMP_INTERVAL + ! READ(11,*);READ(11,*) NO_OF_STEPS + ! READ(11,*);READ(11,*) NO_OF_ITER + ! READ(11,*);READ(11,*) T_INITIAL + ! READ(11,*);READ(11,*) T_FINAL + ! READ(11,*);READ(11,*) TEMP_INTERVAL + ! CLOSE(11) + ! + NROWS = 10 ; NCOLS = 10; + NO_OF_MC_STEPS = 1000000; NO_OF_ITER = 1 + T_INITIAL = 1.0; T_FINAL = 4.0; TEMP_INTERVAL = 0.1 + TEMP_STEPS = INT((T_FINAL - T_INITIAL)/TEMP_INTERVAL) + ! + ALLOCATE(A(NROWS+2,NCOLS+2)) + ALLOCATE(AOLD(NROWS+2,NCOLS+2)) + ALLOCATE(ANEW(NROWS+2,NCOLS+2)) + ! + ! Initialize Spin array + A(:,:) = 1 + ! Write in a OUTPUT file + ! + OPEN(UNIT=31,FILE='Spin_states.out',STATUS='REPLACE', & + & ACTION='WRITE') + ! + OPEN(UNIT=33,FILE='Magnetization.out',STATUS='REPLACE', & + & ACTION='WRITE') + ! WRITE(32,*) T_INITIAL + ! WRITE(32,*) TEMP_INTERVAL + ! WRITE(32,*) T_FINAL + ! + WRITE(33,*) "Temperature ", " Iteration ", " Magnetization" + ! + OPEN(UNIT=32,FILE='temperature_sweep.dat',STATUS='REPLACE', & + & ACTION='WRITE') + WRITE(32,*) " Temperature ", " Magnetization ", & + & " Susceptibility ", " Energy ", " Cv " + ! + MAGNETIZATION = 0.0; ENERGY = 0.0; + ! + TEMP_SCAN: DO TEMP_ITER = 0, TEMP_STEPS + ! + TEMPERATURE = T_INITIAL + TEMP_ITER * TEMP_INTERVAL + BETA = 1.0/TEMPERATURE + MAG_AVG = 0.0; MAG2_AVG = 0.0 + ENERGY_AVG = 0.0; ENERGY2_AVG = 0.0 + ! + TAKE_AVG: DO ITER_AVG = 1, NO_OF_ITER + AOLD = A + ! + MONTE_CARLO: DO MC_ITER = 1, NO_OF_MC_STEPS + ANEW = AOLD + ! Ramdomly change the spin + CALL RANDOM_NUMBER(RAND_NUM) + M = NINT((NROWS-1)*RAND_NUM+2) ! Choose a random row + ! + CALL RANDOM_NUMBER(RAND_NUM) + N = NINT((NCOLS-1)*RAND_NUM+2) ! Choose a random col + TEMP_VAR = -1*ANEW(M,N) + ANEW(M,N) = TEMP_VAR + ! Calculating change in energy due to new configuration + ! + DELTA_E = -2*ANEW(M,N)*(ANEW(M-1,N)+ANEW(M+1,N) & + & + ANEW(M,N-1)+ANEW(M,N+1)) + ! + CALL RANDOM_NUMBER(RAND_NUM) + ! + IF (EXP(-BETA*DELTA_E)> RAND_NUM) THEN + AOLD = ANEW + IF (M == 2) AOLD(NROWS+2,N) = TEMP_VAR + IF (M == NROWS+1) AOLD(1,N) = TEMP_VAR + IF (N == 2) AOLD(M,NCOLS+2) = TEMP_VAR + IF (N == NCOLS+1) AOLD(M,1) = TEMP_VAR + ELSE + CONTINUE + ENDIF + ! + MAGNETIZATION = ABS(SUM(AOLD(2:NROWS+1,2:NCOLS+1))/ & + & (NROWS*NCOLS*1.0)) + ! + IF ((MOD(NINT(TEMPERATURE*10),2) == 0) .AND. & + & MOD(MC_ITER,1000) == 0) THEN + WRITE(33,*) TEMPERATURE, MC_ITER, MAGNETIZATION + ENDIF + ! + ENDDO MONTE_CARLO + ! + MAG_AVG = MAG_AVG + MAGNETIZATION + MAG2_AVG = MAG2_AVG + MAGNETIZATION**2 + ! + DO I = 2, NROWS + 1 + DO J = 2, NCOLS + 1 + ENERGY = ENERGY -ANEW(I,J)*(ANEW(I-1,J)+ANEW(I+1,J) & + & +ANEW(I,J-1)+ANEW(I,J+1)) + ENDDO + ENDDO + ENERGY = ENERGY/(NROWS*NCOLS*2) + ENERGY_AVG = ENERGY_AVG + ENERGY + ENERGY2_AVG = ENERGY2_AVG + ENERGY**2 + ! + ENDDO TAKE_AVG + ! + MAG_AVG = MAG_AVG/NO_OF_ITER + MAG2_AVG = MAG2_AVG/NO_OF_ITER + ENERGY_AVG = ENERGY_AVG/NO_OF_ITER + ENERGY2_AVG = ENERGY2_AVG/NO_OF_ITER + ! + SUSP = BETA*(MAG2_AVG - MAG_AVG**2) + Cv = (BETA/TEMPERATURE)*(ENERGY2_AVG - ENERGY_AVG**2) + ! + WRITE(32,*) TEMPERATURE, MAG_AVG, SUSP, ENERGY, Cv + ! + ! WRITE(31,*) TEMPERATURE + ! WRITE(31,"(10I2)") AOLD(2:NROWS+1,2:NCOLS+1) + 1 + ! WRITE(31,*) "..................................." +! + ENDDO TEMP_SCAN + WRITE(31,*) AOLD(2:NROWS+1,2:NCOLS+1) + CLOSE(31) + CLOSE(32) + CLOSE(33) + print "(12I2)", AOLD + PRINT *, "........................................" + PRINT *, ENERGY/(NROWS*NCOLS*2) + CALL CPU_TIME(FINISH) + PRINT *, " End time ", FINISH + END PROGRAM ISING + diff --git a/magnetization_plot.m b/magnetization_plot.m new file mode 100755 index 0000000..e57821d --- /dev/null +++ b/magnetization_plot.m @@ -0,0 +1,30 @@ +% This program plots all the magnetization data for Ising model for +% iteration vs magnetization for each temperature +% BEFORE execution check ising.f for parameter setup in WRITE(33,*) file +clear all; +clc; +fileID = fopen('Magnetization.out'); +% +header = textscan(fileID,'%s',3); +% update data point count according to Magnetization.out +% Look into isisng.f for exact calculation +datain = textscan(fileID,'%f %d %f',16000); +fclose(fileID); +% +%plot_mag = zeros(1:100); +%X = zeros(1:100); +Magnetization = datain{3}; +for j = 1:16 + for i = 1:1000 % get this number from saved number of iteration + plot_mag(i,j)= Magnetization((j-1)*i+1); + X(i) = i*1000; + end +end +plot(X,plot_mag(:,1),X,plot_mag(:,6),X,plot_mag(:,8),X,plot_mag(:,11),X,plot_mag(:,16)) +grid on +xlabel('# of iterations') +ylabel('|magnetization|') +legend('T=1.0','T=2.0','T=2.4','T=3.0','T=4.0','Location','SouthWest') +set(gcf, 'PaperPosition', [0 0 5 5]); %Position plot at left hand corner with width 5 and height 5. +set(gcf, 'PaperSize', [5 5]); %Set the paper to have width 5 and height 5. +saveas(gcf, 'magnetization', 'pdf') %Save figure diff --git a/makeplots.R b/makeplots.R new file mode 100755 index 0000000..acb0a38 --- /dev/null +++ b/makeplots.R @@ -0,0 +1,39 @@ +wolff <- read.table("temperature_sweep.dat",header=T) +attach(wolff) #Attaches wolff to R search path; wolff$Value becomes Value + + +sus.max = which.max(Susceptibility) +cv.max = which.max(Cv) +print(paste("Susceptibility max occurs at:",Temperature[sus.max]) ,sep = " ") +print(paste("Cv max occurs at:",Temperature[cv.max]) ,sep = " ") + +par(mfrow = c(2,2), oma = c(0,0,2,0), mar = c(5,5,2,2) - .3) +#opens graphics device with 2x2 matrix for plots, 2 line outer margin + +plot(Temperature,Magnetization,type="p", + xlab="Temperature",ylab="") +abline(v = Temperature[sus.max],col = "red",lty = 3) +abline(v = Temperature[cv.max],col = "green",lty = 3) + +plot(Temperature,Susceptibility,type="p", + xlab="Temperature",ylab="Susceptibility") +abline(v = Temperature[sus.max],col = "red",lty = 3) +abline(v = Temperature[cv.max],col = "green",lty = 3) + +plot(Temperature,Energy,type="p", + xlab="Temperature",ylab="") +abline(v = Temperature[sus.max],col = "red",lty = 3) +abline(v = Temperature[cv.max],col = "green",lty = 3) + +plot(Temperature,Cv,type="p", + xlab="Temperature",ylab="") +abline(v = Temperature[sus.max],col = "red",lty = 3) +abline(v = Temperature[cv.max],col = "green",lty = 3) + + +mtext("Monte Carlo simulation of Ising model", line = .5, outer = T) + + +detach(wolff) + +dev.copy2eps(file="system_plots.eps") diff --git a/plotspin.m b/plotspin.m new file mode 100644 index 0000000..d67d0dc --- /dev/null +++ b/plotspin.m @@ -0,0 +1,13 @@ +function plotspin(i,j,k) + if k < 0 + xval = double(i) - 0.5; + yval = double(j) - 0.5; + text(xval,yval,'\downarrow','fontsize',15,'VerticalAlignment', ... + 'middle','HorizontalAlignment','center','Color','m') + else + xval = double(i) - 0.5; + yval = double(j) - 0.5; + text(xval,yval,'\uparrow','fontsize',15,'VerticalAlignment', ... + 'middle','HorizontalAlignment','center','Color','b') + end +end diff --git a/readme.txt b/readme.txt new file mode 100755 index 0000000..7fbb8ab --- /dev/null +++ b/readme.txt @@ -0,0 +1,19 @@ +To execute Ising model follow below instructions (download all the codes from github) +goto https://github.com/tridip66/ising or you found this 'text' here +1) Download ising.f from src folder +2) compile it in HPCC development node by: + f95 ising.f -o ising.exe +3) Execute it by ./ising.exe + +4) Download makeplots.R +Load R module to plot different parameters with temperature +5) module load R +6) R < makeplots.R [--save --no-save --vanilla] + +To visualize spin +7) Download spinviz.m and plotspin.m +8) module load Matlab +9) matlab -nodesktop < spinviz.m +10) you can download 'report' folder +11) Execute report by: pdflatex report.tex + diff --git a/report/figures/fig_1.pdf b/report/figures/fig_1.pdf old mode 100644 new mode 100755 index 629164b..a331634 Binary files a/report/figures/fig_1.pdf and b/report/figures/fig_1.pdf differ diff --git a/report/figures/fig_2.pdf b/report/figures/fig_2.pdf new file mode 100755 index 0000000..9360864 Binary files /dev/null and b/report/figures/fig_2.pdf differ diff --git a/report/figures/fig_3.pdf b/report/figures/fig_3.pdf new file mode 100755 index 0000000..2daafae Binary files /dev/null and b/report/figures/fig_3.pdf differ diff --git a/report/figures/fig_4.pdf b/report/figures/fig_4.pdf new file mode 100755 index 0000000..b0f3c37 Binary files /dev/null and b/report/figures/fig_4.pdf differ diff --git a/report/figures/fig_5.pdf b/report/figures/fig_5.pdf new file mode 100755 index 0000000..3e4e2cf Binary files /dev/null and b/report/figures/fig_5.pdf differ diff --git a/report/figures/fig_6.pdf b/report/figures/fig_6.pdf new file mode 100755 index 0000000..9706390 Binary files /dev/null and b/report/figures/fig_6.pdf differ diff --git a/report/figures/fig_7.pdf b/report/figures/fig_7.pdf new file mode 100755 index 0000000..b67f19e Binary files /dev/null and b/report/figures/fig_7.pdf differ diff --git a/report/report.tex b/report/report.tex old mode 100644 new mode 100755 index 72d1253..11b3416 --- a/report/report.tex +++ b/report/report.tex @@ -10,147 +10,111 @@ \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 +\title{Two-Dimensional Ising Model Simulation:\\with Metropolis Algorithm}% Force line breaks with \\ +\author{Tridip Das} +\date{February 18, 2015}% 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. +This study on Ising model investigates the variation of magnetization with temperature in a ferromagnetic material. +The model try to predict Curie Temperature for a ferromagnetic material with scaled properties. The Monte Carlo simulation with Metropolis algorithm has been +used to predict new spin states by randomly changing spin of a lattice and henceforth trying to minimize the overall energy. +This model also predicts Susceptibility and Specific heat variation with scaled temperature. \end{abstract} - -\section{First-level heading} %Title for the section +\section{Introduction} %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. +The Ising model is a simple representation of the ferromagnetic phase transition at Curie temperature. In this model at low temperature all the lattice site has same spin (say all are up). +The lattice sites preferably occupied by opposite (down) spins with increasing temperature and when it attains Curie temperature the spin distribution goes random in up or down. +The Hamiltonian for Ising model includes only neighbouring site interaction, the long range interaction are neglected. The spin directions may be "Up" (+1) or "Down" (-1). -\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} +\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 + H = - J \sum_{} S_i S_j \end{equation} +J is a coupling constant. For ferromagnetic model J $>$ 0. Here it is assumed to be 1. + +\section{Monte Carlo Simulation} +In the Ising model, each particle can occupy a site with up or down spin state. In Monte Carlo simulation one site is randomly choose, and the spin state of the particle is reversed. +The acceptance of the new state is determined with application of Metropolis algorithm. + +\subsection{Metropolis Algorithm} +The method is a modified Monte Carlo scheme. In this method instead of choosing configuration randomly, configurations are choose with a probability exp(-E/kT) and weight them evenly. +The change in energy is calculated as $\delta E$, if it is $<$ 0, change is accepted. Otherwise, the $\delta E$ with a weighting factor, $e^{-\beta \delta E}$ is compared with a random number $\in [0,1] $ and if the factor is greater, change is accepted. +\\ +\subsection{Calculating Observables} +Some quantitative information can be derived from spin array calculations. The magnetization is calculated from spin $\pm 1$ according to below formula: +\begin{equation} + m = \frac{1}{N}\sum_{i} \langle S_i \rangle +\end{equation} +\\ +The energy, $E$ is calculated as: +\begin{equation} + E = - J \sum_{} S_i S_j +\end{equation} +i,j are neighbouring lattice position of up, down, left and right.\\ +\\ +It is also possible to calculate heat capacity, $C_v$ at constant field according to below formula: +\begin{equation} + C_v = \frac{\beta}{T}[\langle E^2 \rangle - \langle E \rangle ^2] +\end{equation} +Where, $\beta = \frac{1}{k_B T}$ \\ -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*} +We can also calculate isothermal susceptibility, $\chi$ with below equation: +\begin{equation} + \chi = \beta[\langle M^2 \rangle - \langle M \rangle ^2] +\end{equation} -\subsection{Multiline equations} +\section{Results} +The Metropolis based Monte Carlo simulation of Ising model is implemented in Fortran. The code is available in github at \url{https://github.com/tridip66/ising/blob/master/src/ising.f}, also, other codes can be found in the parent directory, which are used +to generate the plots. +\begin{figure}[p] + \centering + \includegraphics[scale=0.9]{figures/fig_1}% Imports a figure - does not automatically scale + \caption{\label{fig:epsart} A magnetization vs Temperature plot for different number of iterations performed} +\end{figure} +\subsection{Parameter selection} +The number of iterations to be performed for the magnetization value to reach equilibrium is plotted in figure 1, a magnetization vs temperature curve. In the plot each line corresponds to different number of iterations. It is observed $10^4$ number of iterations is good enough for equilibrium but based on execution time $10^5$ iterations are performed (look into appendix A for details). +\\ +\begin{figure}[p] + % \centering + \includegraphics[scale=0.8]{figures/fig_4}% Imports a figure - does not automatically scale + \caption{\label{fig:epsart} Temperature vs diferent parameters plot for a 4x4 lattice} +\end{figure} -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. +The lattice site 10x10 is compared with 4x4 and it is observed that 10x10 provides a sharper change in magnetization (figure 3) at Curie temperature. Also, the specific heat and susceptibility change at Curie temperature is not sharp in 4x4 lattice as shown in figure 2. -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. +\subsection{Visualization of spin} +A Matlab code, \url{https://github.com/tridip66/ising/blob/master/spinviz.m} has been developed for visualization of spin states. The figure 4 shows a 10x10 lattice with randomly oriented spins at T=4. +\\ -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] +\begin{figure}[p] \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.} + \includegraphics[scale=0.8]{figures/fig_5}% Imports a figure - does not automatically scale + \caption{\label{fig:epsart} Temperature vs diferent parameters plot for a 10x10 lattice} \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} +\begin{figure}[p] \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} + \includegraphics[scale=0.5]{figures/fig_6}% Imports a figure - does not automatically scale + \caption{\label{fig:epsart} Visualization of spins in a 10x10 lattice at T=4} +\end{figure} -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} +\section*{References} +Charles Kittel and Herbert Kroemer. Thermal Physics. W. H. Freeman, January 1980. +\\ +\url{http://en.wikibooks.org/wiki/LaTeX/Mathematics} +\\ +\url{http://en.wikibooks.org/wiki/LaTeX/Importing_Graphics} +\\ \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}). +\section*{Appendix} +Figure 5 shows the iteration time increases linearly with number of iteration as expected. As $10^5$ steps takes around 2 minutes, which is reasonable time with accuracy, so it is selected for further calculations. +\begin{figure}[p] + \centering + \includegraphics[scale=0.6]{figures/fig_7}% Imports a figure - does not automatically scale + \caption{\label{fig:epsart} Iteration time vs number of iteration in 10x10 lattice} +\end{figure} \end{document} diff --git a/spinviz.m b/spinviz.m new file mode 100755 index 0000000..d03db18 --- /dev/null +++ b/spinviz.m @@ -0,0 +1,57 @@ +function tridip = spinviz +% This program visualizes the spin for Ising model for +% 2D lattice system at the end of the iteration +% BEFORE execution check ising.f for parameter setup in WRITE(31,*) file +clear all; +clc; +fileID = fopen('Spin_states.out'); +% +header = textscan(fileID,'%s %d %s %d',1); +NROWS = header{2}; NCOLS = header{4}; +lattice_sites = NROWS * NCOLS; +datain = textscan(fileID,'%d',lattice_sites); +spins = datain{1}; +fclose(fileID); +% +% Plot grid/box for spins to be displayed +hold on +plot([0 0],[0 NROWS],'k','LineWidth',2); +plot([0 NCOLS],[0 0],'k','LineWidth',2); +% +for i = 1:NROWS + plot([i i],[0 NROWS],'k','LineWidth',2); +end +% +for j = 1:NCOLS + plot([0 NROWS],[j j],'k','LineWidth',2); +end +% Plot spins inside box +for I = 1:NROWS + for J = 1:NCOLS + spin = spins((I-1)* NROWS + J); + plotspin(I,J,spin) + end +end + axis equal + axis off + set(gcf,'Color',[1,1,1]); +hold off +% +set(gcf, 'PaperPosition', [0 0 5 5]); %Position plot at left hand corner with width 5 and height 5. +set(gcf, 'PaperSize', [5 5]); %Set the paper to have width 5 and height 5. +saveas(gcf, 'spin_visualization', 'pdf') %Save figure +end +% +function plotspin(i,j,k) + if k < 0 + xval = double(i) - 0.5; + yval = double(j) - 0.5; + text(xval,yval,'\downarrow','fontsize',15,'VerticalAlignment', ... + 'middle','HorizontalAlignment','center','Color','m') + else + xval = double(i) - 0.5; + yval = double(j) - 0.5; + text(xval,yval,'\uparrow','fontsize',15,'VerticalAlignment', ... + 'middle','HorizontalAlignment','center','Color','b') + end +end \ No newline at end of file diff --git a/src/ising.f b/src/ising.f new file mode 100755 index 0000000..7a0a687 --- /dev/null +++ b/src/ising.f @@ -0,0 +1,192 @@ + PROGRAM ISING + ! 2D Monte Carlo simulation of Ising Model with Metropolis algorithm + ! Program written by Tridip Das (dastridi) + ! This program calculates magnetization with temperature variation + ! Susceptibility, Energy and Specific heat with temperature + ! Also finds the critical temperature at which magnetization is + ! lost. Scaled temperature is used, Boltzmann const, kB assumed as 1 + ! + IMPLICIT NONE + ! Variable Declaration + INTEGER :: I, J, M, N ! counters + INTEGER, ALLOCATABLE :: A(:,:),AOLD(:,:),ANEW(:,:) ! 2D array containing spins + INTEGER :: NROWS, NCOLS ! Defining # of rows and columns + ! + REAL :: TEMPERATURE, BETA, MAGNETIZATION, ENERGY + REAL :: MAG2_AVG, ENERGY2_AVG, SUSP, Cv + REAL :: DELTA_E, MAG_AVG, RAND_NUM, ENERGY_AVG + INTEGER :: TEMP_ITER, ITER_AVG, MC_ITER, TEMP_VAR + INTEGER :: NO_OF_MC_STEPS, NO_OF_ITER, TEMP_STEPS + REAL :: T_INITIAL, T_FINAL, TEMP_INTERVAL, START, FINISH + REAL :: SUMFOR_STDDEV, sTDdEV + REAL, ALLOCATABLE :: MAGNET_DATA(:) + ! + ! + PRINT *, " MONTE CARLO SIMULATION OF 2D ISING MODEL " + PRINT *, " SIMULATION WITH METROPOLIS ALGORITHM " + PRINT *, " APPLYING PERIODIC BOUNDARY CONDITION " + PRINT *, "--------Code developed by Tridip Das------" + CALL CPU_TIME(START) + PRINT *, " Start time ", START + PRINT *, "------------------------------------------" + ! INPUT PARAMETERS FROM inp_file + ! OPEN(UNIT=11,FILE='inpfile',STATUS='OLD',ACTION='READ') + ! Uncomment below lines to get input from user + ! PRINT *, " PROVIDE THE NUMBER OF ROWS FOR ISING MODEL" + ! READ(*,*) NROWS + ! PRINT *, " PROVIDE THE NUMBER OF COLS FOR ISING MODEL" + ! READ(*,*) NCOLS + ! PRINT *, " NUMBER OF SIMULATION STEPS TO RUN FOR MC & +! & EQULIBRIATION " + ! READ(*,*) NO_OF_MC_STEPS + ! PRINT *, " NUMBER OF ITERATIONS TO PERFORM AT EACH TEMPERATURE & +! & STEP TO CALCULATE AVERAGE OVER MC PREDICTION" + ! READ(*,*) NO_OF_ITER + ! PRINT *, " PROVIDE INITIAL, FINAL TEMPERATURE AND INTERVAL" + ! PRINT *, " SEPARATED BY ENTER " + ! READ(*,*) T_INITIAL + ! READ(*,*) T_FINAL + ! READ(*,*) TEMP_INTERVAL + ! READ(11,*);READ(11,*) NO_OF_STEPS + ! READ(11,*);READ(11,*) NO_OF_ITER + ! READ(11,*);READ(11,*) T_INITIAL + ! READ(11,*);READ(11,*) T_FINAL + ! READ(11,*);READ(11,*) TEMP_INTERVAL + ! CLOSE(11) + ! + NROWS = 10 ; NCOLS = 10; + NO_OF_MC_STEPS = 100000; NO_OF_ITER = 100 + T_INITIAL = 1.0; T_FINAL = 4.0; TEMP_INTERVAL = 0.1 + TEMP_STEPS = INT((T_FINAL - T_INITIAL)/TEMP_INTERVAL) + ! + ALLOCATE(A(NROWS+2,NCOLS+2)) + ALLOCATE(AOLD(NROWS+2,NCOLS+2)) + ALLOCATE(ANEW(NROWS+2,NCOLS+2)) + ALLOCATE(MAGNET_DATA(NO_OF_ITER)) + ! + ! Initialize Spin array + A(:,:) = 1 + ! Initialize Magnetization data storage + MAGNET_DATA(:) = 0 + ! Write in a OUTPUT file + ! Write all the final spin states for visualization + OPEN(UNIT=31,FILE='Spin_states.out',STATUS='REPLACE', & + & ACTION='WRITE') + WRITE(31,*) "ROWS ", NROWS, " COLS ", NCOLS + ! + ! Write magnetization data with monte carlo iteration to find the + ! equilibrium zone when magnetization value jumps around equilibrm + OPEN(UNIT=33,FILE='Magnetization.out',STATUS='REPLACE', & + & ACTION='WRITE') + WRITE(33,*) "Temperature ", " Iteration ", " Magnetization" + ! + ! Write all the temperature data with magnetization, + ! susceptibility, energy and others for the plot to find critical + ! temperature + OPEN(UNIT=32,FILE='temperature_sweep.dat',STATUS='REPLACE', & + & ACTION='WRITE') + WRITE(32,*) " Temperature ", " Magnetization ", & + & " Susceptibility ", " Energy ", " Cv ", " sTDdEV " + ! + MAGNETIZATION = 0.0; ENERGY = 0.0; + ! Below loop calculates magnetization at each temperature + TEMP_SCAN: DO TEMP_ITER = 0, TEMP_STEPS + ! + TEMPERATURE = T_INITIAL + TEMP_ITER * TEMP_INTERVAL + BETA = 1.0/TEMPERATURE + MAG_AVG = 0.0; MAG2_AVG = 0.0 + ENERGY_AVG = 0.0; ENERGY2_AVG = 0.0 + ! Below loop takes average on multiple Monte Carlo iteration + TAKE_AVG: DO ITER_AVG = 1, NO_OF_ITER + AOLD = A + ! Below loop performs Monte Carlo simulation by randomly changing + ! the spin states in the 2D lattice + MONTE_CARLO: DO MC_ITER = 1, NO_OF_MC_STEPS + ANEW = AOLD + ! Ramdomly change the spin + CALL RANDOM_NUMBER(RAND_NUM) + M = NINT((NROWS-1)*RAND_NUM+2) ! Choose a random row + ! + CALL RANDOM_NUMBER(RAND_NUM) + N = NINT((NCOLS-1)*RAND_NUM+2) ! Choose a random col + TEMP_VAR = -1*ANEW(M,N) + ANEW(M,N) = TEMP_VAR + ! Calculating change in energy due to new configuration + ! + DELTA_E = -2*ANEW(M,N)*(ANEW(M-1,N)+ANEW(M+1,N) & + & + ANEW(M,N-1)+ANEW(M,N+1)) + ! + CALL RANDOM_NUMBER(RAND_NUM) + ! + IF (EXP(-BETA*DELTA_E)> RAND_NUM) THEN + AOLD = ANEW + IF (M == 2) AOLD(NROWS+2,N) = TEMP_VAR + IF (M == NROWS+1) AOLD(1,N) = TEMP_VAR + IF (N == 2) AOLD(M,NCOLS+2) = TEMP_VAR + IF (N == NCOLS+1) AOLD(M,1) = TEMP_VAR + ELSE + CONTINUE + ENDIF + ! + MAGNETIZATION = ABS(SUM(AOLD(2:NROWS+1,2:NCOLS+1))/ & + & (NROWS*NCOLS*1.0)) + ! Uncomment the write(33,*) if want to print mag with # of iter + IF ((MOD(NINT(TEMPERATURE*10),2) == 0) .AND. & + & MOD(MC_ITER,1000) == 0) THEN + ! WRITE(33,*) TEMPERATURE, MC_ITER, MAGNETIZATION + ENDIF + ! + ENDDO MONTE_CARLO + ! + MAGNET_DATA(ITER_AVG) = MAGNETIZATION + MAG_AVG = MAG_AVG + MAGNETIZATION + MAG2_AVG = MAG2_AVG + MAGNETIZATION**2 + ! + DO I = 2, NROWS + 1 + DO J = 2, NCOLS + 1 + ENERGY = ENERGY -ANEW(I,J)*(ANEW(I-1,J)+ANEW(I+1,J) & + & +ANEW(I,J-1)+ANEW(I,J+1)) + ENDDO + ENDDO + ENERGY = ENERGY/(NROWS*NCOLS*2) + ENERGY_AVG = ENERGY_AVG + ENERGY + ENERGY2_AVG = ENERGY2_AVG + ENERGY**2 + ! + ENDDO TAKE_AVG + ! + MAG_AVG = MAG_AVG/NO_OF_ITER + MAG2_AVG = MAG2_AVG/NO_OF_ITER + ENERGY_AVG = ENERGY_AVG/NO_OF_ITER + ENERGY2_AVG = ENERGY2_AVG/NO_OF_ITER + ! + SUSP = BETA*(MAG2_AVG - MAG_AVG**2) + Cv = (BETA/TEMPERATURE)*(ENERGY2_AVG - ENERGY_AVG**2) + ! Calculate standard deviation on the magnetization data + SUMFOR_STDDEV = 0 + DO I = 1, NO_OF_ITER + SUMFOR_STDDEV = SUMFOR_STDDEV + & + & (MAGNET_DATA(I) - MAG_AVG)**2 + ENDDO + sTDdEV = SQRT(SUMFOR_STDDEV/NO_OF_ITER) + ! + ! Write temp, magnetization and other data to generate plot + ! in the temperature_sweep.dat + WRITE(32,*) TEMPERATURE, MAG_AVG, SUSP, ENERGY, Cv, sTDdEV + ! + ! WRITE(31,*) TEMPERATURE + ! WRITE(31,"(10I2)") AOLD(2:NROWS+1,2:NCOLS+1) + 1 + ! WRITE(31,*) "..................................." +! + ENDDO TEMP_SCAN + ! Write spin states to Spin_states.out for visualization of spins + WRITE(31,*) AOLD(2:NROWS+1,2:NCOLS+1) + CLOSE(31) + CLOSE(32) + CLOSE(33) + print "(12I2)", AOLD + PRINT *, "........................................" + ! PRINT *, ENERGY/(NROWS*NCOLS*2) + CALL CPU_TIME(FINISH) + PRINT *, " End time ", FINISH + END PROGRAM ISING +