\documentclass[12pt]{article}
\usepackage{eamc-template}
\usepackage[round,sort]{natbib}
\usepackage[brazil]{babel}
% \usepackage[latin1]{inputenc}
\usepackage[utf8]{inputenc}
\usepackage{graphicx,url}
\usepackage{float}
\usepackage{tabularx}
\usepackage{booktabs}
\usepackage{blindtext}
%necessario para formulas matematicas
\usepackage{amsmath}
%necessario para alguns simbolos matematicos, como \quad
\usepackage{amssymb}
%habilita modo codigo dentro do texto
\usepackage{verbatim}
%%%%%%% formatacao para escrita de pseudocodigo
\usepackage{algorithm}
\usepackage[noend]{algpseudocode}
\usepackage{pifont}
\makeatletter
\def\BState{\State\hskip-\ALG@thistlm}
\makeatother
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{hyperref} % necessario para usar \url{} e \href{}
%para nao ficar o retangulo em volta dos links, apenas muda a cor dos caracteres
\hypersetup{ colorlinks,
linkcolor=blue,
filecolor=blue,
urlcolor=blue,
citecolor=blue }
%package para fazer figuras a e b automaticamente
\usepackage{subfigure}
%\usepackage[pdftex]{graphicx}
\sloppy
\title{Scientific computing in Statistical Mechanics: time decay of orientation order in 2D hard disk system}
\author{A.G. França\inst{1}, R.S. Grisotto\inst{2} B.M. Rocha\inst{3}, F.A.A.M.N. Soares\inst{3}, H.A. Fernandes\inst{1}, P.F. Gomes\inst{1}}
\address{Instituto de Ciência Exatas, Universidade Federal de Goiás\\
Jataí, GO -- Brasil
\nextinstitute
Instituto de Computação -- Universidade Estadual de Campinas\\
Campinas, SP -- Brasil
\nextinstitute
Instituto de Informática -- Universidade Federal de Goiás\\
Goiânia, GO -- Brasil
\email{andreygfranca@gmail.com, paulofisicajatai@gmail.com }
% \email{ andreygfranca, @gmail.com,rafaelgrisotto@gmail.com, rocha3runo@gmail.com, fabrizzio.soares@gmail.com, ha.fernandes@gmail.com,paulofisicajatai@gmail.com}
}
\begin{document}
\maketitle
\begin{abstract}
In this work we successfully applied diverse computational techniques to calculate important quantities in one problem in Statistical Mechanics: the hard disk system in two dimensions. We calculated the global and local orientation order and the time decay constant of the global orientation correlation function. We also computed the Voronoi construction to visuallize the spatial distribution of the local orientation order.
\end{abstract}
\begin{resumo}
Neste trabalho apresentamos diversas técnicas computacionais no cálculo de grandezas de interesse em um sistema de Mecânica Estatística: discos rígidos em duas dimensões. Calculamos a ordem de orientação local e global e a constante de decaimento da função de correlação da orientação global. Fizemos também a diagramação de Voronoi para visualização da distribuição espacial da orientação local.
\end{resumo}
\section{Introduction}
Solidification and vaporization of water are two commom examples of phase transitions in physics. The temperature and pressure defines which phase the water is. Another parameter called order parameter is used to verify which phase the system is. In the case of the water it is the molecule organization: in the solid phase the molecules form an ordered hexagonal state, in the liquid state the molecules are close to each other but with no organization and in the gas state there is no correlation between the molecules. One interesting system to study phase transitions is the hard disks in two dimensions [\cite{Isobe2015}], composed of $N$ disks of radius $R$ confined in a 2D planar rectangular box of sides $L_x$ and $L_y$ (see Fig. \ref{disk3}). Each disk $k$ has its center at position $(x_k,y_k)$. The superposition between two disks is forbidden so the minimal distance between them is $2R$. One parameter is the local orientation of each disk $k$ defined as [\cite{Deutschlander2014,Engel2013}]:
\begin{equation}
\psi_k(x_k,y_k) = \dfrac{1}{N_k} \sum_{l=1}^6 e^{6i \theta_{kl}} , \label{defiorientalocal}
\end{equation}
where $i = \sqrt{-1}$ and the sum is taken over the six nearest-neighbors. The angle $\theta_{kl}$ is between a line passing through the center of the two discs and the $x$ axis (see Fig. \ref{angulo_vizinho}). The $\psi_k$ is defined so that $\psi_k = 0$ in the liquid phase and $\psi_k = 1$ in the solid one. The spatial distribution of the parameter through the volume is also important. The order parameter can have different ranges defined by its dependency on the distance $r$ [\cite{Deutschlander2014}]: exponential decay $e^{-r/\xi}$ means short range, algebraic one $r^{-\xi}$ means quasi-long range and a constant variation means a long range. So the liquid state in the 2D hard disk system has short range and the solid state has long range for $\psi_k$. However, a third phase was theoretical and experimentally observed: the so called hexatic phase [\cite{Strandburg1988}] (see Fig. \ref{transicoes}) and it has quasi-long range orientational order [\cite{Bernard2011}]. The hexatic phase is an exclusive feature of 2D systems [\cite{Isobe2015}].
\subsection{Global orientation order}
The global orientation is the absolute value of the average of the local one:
\begin{equation}
\Psi = \left\vert \dfrac{1}{N} \sum_{k=1}^N \psi_k \right\vert. \label{orientacaoglobal}
\end{equation}
The phase of the system depends on the density $\eta = (N\pi R^2)/(L_x L_y)$ (see Fig. \ref{transicoes}). The solid phase has an hexagonal symmetry (see Fig. \ref{triangular_N1024_eta0p5}) and the liquid phase has random positions for the disks (see Fig. \ref{final_q_N1024_eta0p5_Q200}). We study the time behavior of the orientation order correlation function:
\begin{equation}
C(\delta) = \sum_{t=0}^Q \Phi(t) \Phi (t+\delta) , \label{correcdelta}
\end{equation}
which should be normalized so $C(0)=1$. The normalized orientation order is $\Phi(t) = \Psi(t) - \Psi_{\infty}$ where $\Psi_{\infty} = \lim_{t \rightarrow \infty} \Psi ( t)$\footnote{In practical terms, $\Psi_{\infty} = \Psi (t=Q)$ where $Q$ is the number of evaluated Monte Carlo iterations.}. The correlation function has an exponential decay [\cite{Bernard2009}] of the form $C(\delta) \approx e^{-t/\tau}$ where the decay time $\tau$ is a function of the number of disks, $\tau = \tau (N)$. This function depends of the phase:
\begin{equation}
\tau (N) = \left\{
\begin{array}{rl}
N \log N, \text{ liquid phase} \\
\\
cN^{\alpha}, \text{ hexatic phase} \label{conshexatica} \\
\\
N, \text{ solid phase}
\end{array}
\right.
\end{equation}
In this work we study the time decay $\tau$ of the global orientation in the liquid-hexatic transition as function of $N$ in the 2D hard disk system. We use different computation technique to evaluate the time evolution and the orientation order. The goal is to verify the exponential decay of the orientation order correlation function and then obtain the behavior of the time constant as function of the number of disks $N$. To our knowledge, and complete study of the $\tau$ as function of $N$ in the hexatic-liquid coexistence phase ($0.70 < \eta < 0.716$) has not been done.
\begin{figure}[!h]
\centering
\subfigure{a)
\includegraphics[width= 2.5 in]{disks3b.pdf}
\label{disk3}
}
\subfigure{b)
\includegraphics[width=1.5 in]{angulo_vizinho.png}
\label{angulo_vizinho}
}
\caption{\subref{disk3} Illustration of the 2D hard disk system. \subref{angulo_vizinho} Illustration of the angle $\theta_{jk}$ between two discs. This Fig. was taken from Ref. [\cite{Bernard2011thesis}].}
\end{figure}
\begin{figure}[!h]
\centering
\subfigure{a)
\includegraphics[width=3.4 in]{transicoes.pdf}
\label{transicoes}
}
\subfigure{b)
\includegraphics[width=1.0 in]{move_ECMC.png}
\label{move_ECMC}
}
\caption{\subref{transicoes} Identification of the phase as function of the density $\eta = (N\pi R^2)/(L_x L_y)$. For $\eta <0.700$ the system is liquid (L), for $0.700 < \eta < 0.716$ it is in the coexistence between the liquid and hexatic (L-H) phases and for $0.716 < \eta \lesssim 0.720$ the pase is hexatic (H). For $\eta \gtrsim 0.72$ the phase is solid (S). \subref{move_ECMC} Illustration of the movement in the Event-chain algorithm. Collective move with the total displacement $\ell$. The Figures were taken from Ref. [\cite{Bernard2011thesis}].}
\end{figure}
\subsection{Simulation Method}
We used two methods for the time evolution of the system. The first one is the Event Chain Monte Carlo method [\cite{Michel2014}] where many disks move in sequence in each iteration. Each disk moves until it reaches another disk (when their distance $=2R$). Then the first disk stops and the second one starts to move until it reaches the third disk. The process repeats until the sum of the distances of all disks are equal to an input value $\ell$ [\cite{Bernard2009}]. The method is efficient because many disks move in one iteration. The other method is the parallel Markov Chain Monte Carlo one, calculated with the Hard Particle Monte Carlo module (HPMC) [\cite{Anderson2016}] of HOOMD-BLUE package\footnote{\href{http://glotzerlab.engin.umich.edu/hoomd-blue}{glotzerlab.engin.umich.edu/hoomd-blue}.} [\cite{Anderson2008,Glaser2015}]. This implementation has been used in many type of problems including disks [\cite{Anderson2017}].
To calculate $\psi_k$ we need the distances between the disk $k$ and all other disks. Using brute force one requires $N-1$ calculations. Then, to calculate the average $\Psi$, another $N$ calculations will be done with total of $(N-1)N \simeq N^2$ calculations when $N$ is large. The distances are also required in the time evolution of the system. We implemented this algorithm in C language. However, for large $N$ we used the method List of Cells [\cite{Frenkel02}] which is faster. In this one a grid of $N$ cells is created so there is one disk per cell. When we need to get the neighbor of a disk, we just look into the neighbor cells. The scaling of this method is $N$. We implemented this list of cells method in C$++$.
\section{Results}
Our C$++$ implementation is illustrated in the Algorithm \ref{pseudocodigo}. The distance calculations are performed in the functions \verb|newL| and \verb|abs_Psi|, where the \verb|listcell| structure is used to identify the neighbor disk. The matrix \verb|LxLy| keeps the positions $x(t),y(t)$ of all disks in the iteration $t$\footnote{Although the concepts of the method and the definition of $\Psi$ are not too complicated, the code itself presents many numerical challenges [\cite{Bernard2011thesis}].}. In the Fig. \ref{diagrama_celulas} there is a simple example of the list with 9 cells. To get the neighbor of a disk, just check the neighbor cells (each cell has pointer pointing to them and a linked list to keep the id of the disks). The Fig. \ref{time_cells_brute} shows a comparison of execution time between brute force and list of cells methods. For small $N$, brute force is better because there is a cost to create the cell grid. However, for $N=128^2$ the list of cells is already one order of magnitude better, and this improvement increases with $N$.
%
%\begin{algorithm}
%\caption{Time evolution with Event Chain Monte Carlo}\label{evolucalslks}
%\begin{algorithmic}[1]
%\Procedure{MCCE}{}
%\State \textbf{Input}: $N$, $\ell$, $\eta$, $S$, $Q$, $L_x$
%\State listcell = \emph{create\_list}($N$)
%\State \textbf{for} {k = 1,$S$}
%\State \qquad LxLy = \emph{initial}($N$,$\eta$, $L_x$)
%\State \qquad \emph{add\_disks}($N$,listcell)
%\State \qquad seed = k
%\State \qquad \textbf{for} {t = 1,$Q$}
%\State \qquad \qquad LxLy = \emph{newL}(LxLy,$N$,$\eta$,$\ell$,seed,listcell)
%\State \qquad \qquad Psi(t) = \emph{abs\_Psi}($\eta$,LxLy,listcell)
%\State \qquad \textbf{end};
%\State \textbf{end};
%\State file = \textit{time\_Psi.csv}
%\State \textbf{write}(t,Psi,file)
%\EndProcedure
%\end{algorithmic}
%\label{pseudocodigo}
%\end{algorithm}
\begin{algorithm}
\caption{Time evolution with Event Chain Monte Carlo}\label{evolucalslks}
\begin{algorithmic}[1]
\Procedure{MCCE}{}
\State \textbf{Input}: $N, \ell, \eta, S, Q, L_x$
\State listcell $\gets$ \Call{create\_list}{$N$}
\For{k $\in S$}
\State LxLy $\gets$ \Call{initial}{$N$, $\eta$, $L_x$}
\State \Call{add\_disks}{$N$, listcell}
\State seed $\gets k$
\For{t $\in Q$}
\State LxLy $\gets$ \Call{newL}{LxLy, $N, \eta, \ell$, seed, listcell}
\State $Psi(t) \gets$ \Call{abs\_Psi}{$\eta$, LxLy, listcell}
\EndFor
\EndFor
\State file $\gets$ \textit{time\_Psi.csv}
\State \textbf{write}(t,$Psi$,file)
\EndProcedure
\end{algorithmic}
\label{pseudocodigo}
\end{algorithm}
\begin{figure}[!h]
\centering
\subfigure{a)
\includegraphics[width=2.8 in]{diagrama_celulas.pdf}
\label{diagrama_celulas}
?}
\subfigure{b)
\includegraphics[width=2.6 in]{time_cells_brute.pdf}
\label{time_cells_brute}
}
\caption{\subref{diagrama_celulas} System with 9 cells, each one has pointer pointing to the neighbors and a linked list. \subref{time_cells_brute} Comparison between the two methods to calculate the distances between the disks: brute force in blue and list of cells in red.}
\end{figure}
In total, the code has six input parameters: (1) number of disks $N$, (2) the density $\eta$, (3) the number of iterations $Q$, (4) the side $L_x$ of the box, (5) the number $S$ of samples in order to take the average (and decrease the error) and (6) the total distance $\ell$. We used $S=$ 5, $L_x=$ 100 and $\ell=\sqrt{N}$ and $Q=$ 30 000 in all simulations. All the graphics presented here were done in Python (Matplotlib package). The way to achieve the liquid state configuration shown in Fig. \ref{final_q_N1024_eta0p5_Q200} was to evaluate the system with $Q=200$ iterations and $\eta = 0.5$ using our C$++$ implementation and an initial hexagonal configuration (Fig. \ref{triangular_N1024_eta0p5}).
\begin{figure}[!h]
\centering
\subfigure{a)
\includegraphics[width=2.3 in]{triangular_N1024_eta0p5.pdf}
\label{triangular_N1024_eta0p5}
}
\subfigure{b)
\includegraphics[width=2.3 in]{final_q_N1024_eta0p5_Q200.pdf}
\label{final_q_N1024_eta0p5_Q200}
}
\caption{Simulation with $N=32^2$ and $\eta=0.5$ using our C$++$ implementation. \subref{triangular_N1024_eta0p5} Initial configuration with triangular symmetry. \subref{final_q_N1024_eta0p5_Q200} Configuration with $t=200$.}
\end{figure}
The correlation functions $C(\delta)$ were calculated for $N=64^2$, $128^2$, $256^2$ and $512^2$. Due to time limitations we were not able to perform simulations with a larger number of $S$. The results shown in the Fig. \ref{figCorr_list} were noisy. To find the decay constant it is needed more samples. Another result is shown in the Fig. \ref{voronoi_128}, which is a kind of graphic used to visually identify the phases [\cite{Bernard2011,Engel2013}]. First a Voronoi construction is created where each cell contains only one disk. Then, the local orientation order $\psi_k$ is calculated for each disk and a color coding is created. The disk in each cell are removed and each cell are painted with the color of its disk orientation order [\cite{Bernard2011thesis}]. If $\psi_k$ is alligned with $\Psi$ we have $\cos \theta_k=1$ and if they are perpendicular $\cos \theta_k=0$, using the associated vectors with $\psi_k$ and $\Psi$. Note that in Fig. \ref{voronoi_128} there are $N=128^2$ disks. The number of sides of each cell is the number of its neighbors cells. As the solid phase of the hard disk system has hexagonal symmetry, the average number of neighbors cell is 6, which explains why most of the cells are hexagons. As the density is $\eta = 0.7$, the configuration of Fig. \ref{voronoi_128} is liquid.
The functions of HOOMD-BLUE package are organized in modules, which should be installed as regular Python packages (within Conda environment). The results for $C(\delta)$ are shown in the Fig. \ref{figuraCorr} for $\eta = 0.705$ and $N=64^2$, $128^2$, $256^2$ and $512^2$. To calculate the decay time $\tau$ we adjust an exponential function $f(\delta) = A_0+A_1e^{-\delta / \tau}$ to the data. The blue circles in Fig. \ref{Corr_fit_eta705} are the obtained values of $\tau$. The second step is to determine $\alpha$ from Eq. \ref{conshexatica}. So we take the natural logarithm $\ln \tau = B + \alpha \ln N$ (with $B=\ln c$) and make a linear fit for $\ln \tau$ vs. $\ln N$. The slope will be the constant $\alpha$. Indeed, they follow a linear behavior as the red solid line indicates in Fig. \ref{Corr_fit_eta705}\footnote{The code that generated the results presented here are available at github: \href{https://github.com/AndreyGFranca/mcec/tree/master/EAMC_2018}{github.com/AndreyGFranca/mcec}.}.
\begin{figure}[!h]
\centering
\subfigure{a)
\includegraphics[width=2.5 in]{figCorr_list.pdf}
\label{figCorr_list}
}
\subfigure{b)
\includegraphics[width=2.8 in]{voronoi_128.pdf}
\label{voronoi_128}
}
\subfigure{c)
\includegraphics[width=2.5 in]{figuraCorr.pdf}
\label{figuraCorr}
}
\subfigure{d)
\includegraphics[width=2.5 in]{Corr_fit_eta705.pdf}
\label{Corr_fit_eta705}
}
\caption{Time behavior of correlation function $C(\delta)$ (Eq. \ref{correcdelta}). \subref{figCorr_list} Using our C$++$ implementation. \subref{voronoi_128} Voronoi construction where each cell contains only one disk. \subref{figuraCorr} Using HOOMD. \subref{Corr_fit_eta705} Linear fit of the decay $\tau$ vs. the number of disks $N$ using the results from HOOMD case.}
\end{figure}
\section{Conclusions}
We compute the correlation function of the orientation order in the 2D hard disk system and verified its exponential behavior. We then obtain the expected behavior of the decay constant $\tau$ as function of the number of disks $N$ for $\eta = 0.705$. We successfully applied various computation techiniques to this Statistical mechanical problem, showing the great benefit which one can have from the interplay between Statistical Mechanics and computational modelling.
\section{Aknowledgements}
This work was funded in part by CNPq (A. G. França and R.S. Grisotto). We thank W. Krauth for the suggestions.
\section{References}
%\bibliographystyle{./bib/plainnatEN}
\bibliographystyle{plainnat}
%\bibliographystyle{./bib/plainnatBR}
\bibliography{referencias_EAMC}
%\bibliography{./bib/eamc-template}
\end{document}
