SklogWiki - User contributions [en]

Grand canonical Monte Carlo

2012-11-22T15:29:04Z

134.109.16.22: /* Theoretical basis */ fixed typo and link

'''Grand-canonical ensemble Monte Carlo''' (GCEMC or GCMC) is a very versatile and powerful [[Monte Carlo]] technique that explicitly accounts for density fluctuations at fixed volume and [[temperature]]. This is achieved by means of trial insertion and deletion of molecules. Although this feature has made it the preferred choice for the study of [[Interface |interfacial]] phenomena, in the last decade
[[grand canonical ensemble | grand-canonical ensemble]] simulations have also found widespread applications in the study of bulk properties. Such applications had been hitherto limited by the very low particle insertion and deletion probabilities, but the development of the [[Configurational bias Monte Carlo |configurational bias]] grand canonical technique has very much improved the situation.

== Theoretical basis ==
In the grand canonical ensemble, one first chooses [[Random numbers |randomly]] whether
a trial particle insertion or deletion is attempted. If insertion is chosen,
a particle is placed with uniform probability density inside the system.
If deletion is chosen, then one deletes one out of <math>N</math> particles
randomly. The trial move is then accepted or rejected according to the
usual Monte Carlo lottery.
As usual, a trial move from an original state (<math>o</math>) to a new state (<math>n</math>) is accepted with probability

:<math> acc(o \rightarrow n) = min \left (1, q \right ) </math>

where <math>q</math> is given by:

:<math> q = \frac{ \alpha(n \rightarrow o)}{\alpha(o \rightarrow n)} \times \frac{f(n)}{f(o)} \qquad\qquad\text{(1)} </math>

Here, <math> \alpha(o \rightarrow n) </math> is the probability density of
attempting trial move from state <math>o</math> to state <math>n</math> (also known as underlying
probability), while <math>f(o)</math> is the probability density of state <math>o</math>.

As noted by Norman and Filinov <ref>G. E. Norman and V. S. Filinov "Investigations of phase transitions by a Monte-Carlo method", High Temperature '''7''' pp. 216-222 (1969)</ref>, evaluation of the proper acceptance rules requires very careful interpretation of the (classical) grand canonical probability density:

:<math> f_L({\mathbf r}_1,{\mathbf r}_2 ..., {\mathbf r}_N) \propto \frac{\Lambda^{-3N}}{N!} e^{\beta \mu N} e^{-\beta U_N}\qquad\qquad\text{(2)} </math>

where <math>N</math> is the total number of particles, <math>\mu</math> is the [[chemical potential]], <math>\beta := 1/k_B T</math> and <math>\Lambda</math> is the [[De Broglie thermal wavelength|de Broglie thermal wavelength]].

The sub-index ''L'' makes emphasis on a particular definition of [[microstate]] in which particles are assumed to be distinguishable. Thus, <math> f_L </math> should be interpreted as the probability density of a classical state with labelled particles, having labelled particle 1 in position <math>{\mathbf r}_1</math>, labelled particle 2 in position <math>{\mathbf r}_2</math> and so on. Since labelling of the particles has no physical significance whatsoever, there are <math>N!</math> identical states which result from permutation of the labels (this explains the <math>N!</math> term in the denominator). Hence, the probability of the significant microstate of undistinguishable particles, i.e., one with <math>N</math> particles at positions <math> {\mathbf r}_1</math>, <math> {\mathbf r}_2</math>, etc., irrespective of the labels, will be given by:

:<math> f_U( \{{\mathbf r}_1,{\mathbf r}_2 ..., {\mathbf r}_N \} ) \propto \sum_P
f({\mathbf r}_1,{\mathbf r}_2 ..., {\mathbf r}_N) = \Lambda^{-3N} e^{\beta \mu N} e^{-\beta U_N}</math>

where the <math>U</math> subindex stands for ''unlabelled'', the coordinates in parenthesis indicate that the positions are not attributed to any particular choice of labelling and the sum runs over all possible particle label permutations.

Upon trial insertion of an extra particle, one obtains:

:<math> \frac{f_U(N+1)}{f_U(N)} = \Lambda^{-3} e^{\beta \mu } e^{-\beta ( U_{N+1} - U_N )}\qquad\qquad\text{(3)} </math>

The probability density of attempting an insertion is <math> \alpha( N \rightarrow N+1 ) = \frac{1}{2} \frac{1}{V} </math>. The <math> 1/2 </math> factor accounts for the probability of attempting an insertion (from the choice of insertion or deletion). The <math> 1/V</math> factor results from placing the particle with uniform probability anywhere inside the simulation box. The reverse attempt (moving from state of <math>N+1</math> particles to the original <math>N</math> particle state) is chosen with probability:

:<math> \alpha( N+1 \rightarrow N ) = \frac{1}{2} \frac{1}{N+1} </math>

where the <math> 1/(N+1) </math> factor results from random removal of one among <math>N+1</math> particles. Therefore, the ratio of underlying probabilities is:

:<math> \frac{\alpha(n\rightarrow o)}{\alpha(o\rightarrow n)} =
\frac{\alpha( N+1 \rightarrow N )}{\alpha( N \rightarrow N+1)} = \frac{V}{N+1} \qquad\qquad\text{(4)}</math>

Substitution of Eq.(3) and Eq.(4) into Eq.(1) yields the
acceptance probability for attempted insertions:

:<math> acc(N \rightarrow N+1) = \frac{V \Lambda^{-3} }{N+1} e^{\beta \mu } e^{-\beta ( U_{N+1} - U_N )} </math>

For the inverse deletion process, similar arguments yield:

:<math> acc(N \rightarrow N-1) = \frac{N}{V \Lambda^{-3} } e^{-\beta \mu } e^{-\beta ( U_{N-1} - U_N )} </math>

Alternatively, one could derive the acceptance rules by considering the probability density of labelled states, Eq.(2) but taking into account that there are then <math>N+1</math> labelled microstates leading to the original <math>N</math> particle labelled state upon deletion (one for each possible label permutation of the deleted particle).

The same acceptance rules are obtained in reference books. Usually the problem of proper counting of states is circumvented by ignoring the labelling problem and assuming that the underlying probabilities for insertion and removal are equal.

== References ==
<references/>
'''Related reading'''
*[http://dx.doi.org/10.1080/00268977400102551 D. J. Adams "Chemical potential of hard-sphere fluids by Monte Carlo methods", Molecular Physics '''28''' pp. 1241-1252 (1974)]
*[http://dx.doi.org/10.1080/00268977500100221 D. J. Adams "Grand canonical ensemble Monte Carlo for a Lennard-Jones fluid", Molecular Physics '''29''' pp. 307-311 (1975)]
*[http://dx.doi.org/10.1063/1.2839302 Attila Malasics, Dirk Gillespie, and Dezső Boda "Simulating prescribed particle densities in the grand canonical ensemble using iterative algorithms", Journal of Chemical Physics '''128''' 124102 (2008)]
[[Category: Monte Carlo]]

Grand canonical Monte Carlo

2012-11-22T15:04:37Z

134.109.16.22: /* Theoretical basis */ typos corrected

'''Grand-canonical ensemble Monte Carlo''' (GCEMC or GCMC) is a very versatile and powerful [[Monte Carlo]] technique that explicitly accounts for density fluctuations at fixed volume and [[temperature]]. This is achieved by means of trial insertion and deletion of molecules. Although this feature has made it the preferred choice for the study of [[Interface |interfacial]] phenomena, in the last decade
[[grand canonical ensemble | grand-canonical ensemble]] simulations have also found widespread applications in the study of bulk properties. Such applications had been hitherto limited by the very low particle insertion and deletion probabilities, but the development of the [[Configurational bias Monte Carlo |configurational bias]] grand canonical technique has very much improved the situation.

== Theoretical basis ==
In the grand canonical ensemble, one first chooses [[Random numbers |randomly]] whether
a trial particle insertion or deletion is attempted. If insertion is chosen,
a particle is placed with uniform probability density inside the system.
If deletion is chosen, then one deletes one out of <math>N</math> particles
randomly. The trial move is then accepted or rejected according to the
usual Monte Carlo lottery.
As usual, a trial move from an original state (<math>o</math>) to a new state (<math>n</math>) is accepted with probability

:<math> acc(o \rightarrow n) = min \left (1, q \right ) </math>

where <math>q</math> is given by:

:<math> q = \frac{ \alpha(n \rightarrow o)}{\alpha(o \rightarrow n)} \times \frac{f(n)}{f(o)} \qquad\qquad\text{(1)} </math>

Here, <math> \alpha(o \rightarrow n) </math> is the probability density of
attempting trial move from state <math>o</math> to state <math>n</math> (also known as underlying
probability), while <math>f(o)</math> is the probability density of state <math>o</math>.

As noted by Norman and Filinov <ref>G. E. Norman and V. S. Filinov "Investigations of phase transitions by a Monte-Carlo method", High Temperature '''7''' pp. 216-222 (1969)</ref>, evaluation of the proper acceptance rules requires very careful interpretation of the (classical) grand canonical probability density:

:<math> f_L({\mathbf r}_1,{\mathbf r}_2 ..., {\mathbf r}_N) \propto \frac{\Lambda^{-3N}}{N!} e^{\beta \mu N} e^{-\beta U_N}\qquad\qquad\text{(2)} </math>

where <math>N</math> is the total number of particles, <math>\mu</math> is the [[chemical potential]], <math>\beta := 1/k_B T</math> and <math>\Lambda</math> is the de [[De Broglie thermal wavelength]].

The sub-index ''L'' makes emphasis on a particular definition of [[microstate]] in which particles are assumed to be distinguishable. Thus, <math> f_L </math> should be interpreted as the probability density of a classical state with labelled particles, having labelled particle 1 in position <math>{\mathbf r}_1</math>, labelled particle 2 in position <math>{\mathbf r}_2</math> and so on. Since labelling of the particles has no physical significance whatsoever, there are <math>N!</math> identical states which result from permutation of the labels (this explains the <math>N!</math> term in the denominator). Hence, the probability of the significant microstate of undistinguishable particles, i.e., one with <math>N</math> particles at positions <math> {\mathbf r}_1</math>, <math> {\mathbf r}_2</math>, etc., irrespective of the labels, will be given by:

:<math> f_U( \{{\mathbf r}_1,{\mathbf r}_2 ..., {\mathbf r}_N \} ) \propto \sum_P
f({\mathbf r}_1,{\mathbf r}_2 ..., {\mathbf r}_N) = \Lambda^{-3N} e^{\beta \mu N} e^{-\beta U_N}</math>

where the <math>U</math> subindex stands for ''unlabelled'', the coordinates in parenthesis indicate that the positions are not attributed to any particular choice of labelling and the sum runs over all possible particle label permutations.

Upon trial insertion of an extra particle, one obtains:

:<math> \frac{f_U(N+1)}{f_U(N)} = \Lambda^{-3} e^{\beta \mu } e^{-\beta ( U_{N+1} - U_N )}\qquad\qquad\text{(3)} </math>

The probability density of attempting an insertion is <math> \alpha( N \rightarrow N+1 ) = \frac{1}{2} \frac{1}{V} </math>. The <math> 1/2 </math> factor accounts for the probability of attempting an insertion (from the choice of insertion or deletion). The <math> 1/V</math> factor results from placing the particle with uniform probability anywhere inside the simulation box. The reverse attempt (moving from state of <math>N+1</math> particles to the original <math>N</math> particle state) is chosen with probability:

:<math> \alpha( N+1 \rightarrow N ) = \frac{1}{2} \frac{1}{N+1} </math>

where the <math> 1/(N+1) </math> factor results from random removal of one among <math>N+1</math> particles. Therefore, the ratio of underlying probabilities is:

:<math> \frac{\alpha(n\rightarrow o)}{\alpha(o\rightarrow n)} =
\frac{\alpha( N+1 \rightarrow N )}{\alpha( N \rightarrow N+1)} = \frac{V}{N+1} \qquad\qquad\text{(4)}</math>

Substitution of Eq.(3) and Eq.(4) into Eq.(1) yields the
acceptance probability for attempted insertions:

:<math> acc(N \rightarrow N+1) = \frac{V \Lambda^{-3} }{N+1} e^{\beta \mu } e^{-\beta ( U_{N+1} - U_N )} </math>

For the inverse deletion process, similar arguments yield:

:<math> acc(N \rightarrow N-1) = \frac{N}{V \Lambda^{-3} } e^{-\beta \mu } e^{-\beta ( U_{N-1} - U_N )} </math>

Alternatively, one could derive the acceptance rules by considering the probability density of labelled states, Eq.(2) but taking into account that there are then <math>N+1</math> labelled microstates leading to the original <math>N</math> particle labelled state upon deletion (one for each possible label permutation of the deleted particle).

The same acceptance rules are obtained in reference books. Usually the problem of proper counting of states is circumvented by ignoring the labelling problem and assuming that the underlying probabilities for insertion and removal are equal.

== References ==
<references/>
'''Related reading'''
*[http://dx.doi.org/10.1080/00268977400102551 D. J. Adams "Chemical potential of hard-sphere fluids by Monte Carlo methods", Molecular Physics '''28''' pp. 1241-1252 (1974)]
*[http://dx.doi.org/10.1080/00268977500100221 D. J. Adams "Grand canonical ensemble Monte Carlo for a Lennard-Jones fluid", Molecular Physics '''29''' pp. 307-311 (1975)]
*[http://dx.doi.org/10.1063/1.2839302 Attila Malasics, Dirk Gillespie, and Dezső Boda "Simulating prescribed particle densities in the grand canonical ensemble using iterative algorithms", Journal of Chemical Physics '''128''' 124102 (2008)]
[[Category: Monte Carlo]]

Metropolis Monte Carlo

2012-06-13T08:37:45Z

134.109.16.22: /* Main features */ grammar

The '''Metropolis Monte Carlo''' technique
<ref>[http://dx.doi.org/10.1063/1.1699114 Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller and Edward Teller, "Equation of State Calculations by Fast Computing Machines", Journal of Chemical Physics '''21''' pp.1087-1092 (1953)]</ref>
is a variant of the original [[Monte Carlo]] method proposed by [[Nicholas Metropolis]] and [[Stanislaw Ulam]] in 1949 <ref>[http://links.jstor.org/sici?sici=0162-1459%28194909%2944%3A247%3C335%3ATMCM%3E2.0.CO%3B2-3 Nicholas Metropolis and S. Ulam "The Monte Carlo Method", Journal of the American Statistical Association '''44''' pp. 335-341 (1949)]</ref>
== Main features ==
Metropolis Monte Carlo simulations can be carried out in different ensembles. For the case of one-component systems the usual ensembles are:
* [[Canonical ensemble]] (<math> NVT </math> )
* [[Isothermal-isobaric ensemble]] (<math> NpT </math>)
* [[Grand canonical ensemble]] (<math> \mu V T </math>)
In the case of mixtures, it is useful to consider the so-called [[Semi-grand ensembles]].
The purpose of these techniques is to sample representative configurations of the system at the corresponding
thermodynamic conditions. The sampling techniques make use of the so-called pseudo-[[Random numbers |random number]] generators.

== Configuration ==

A configuration is a microscopic realisation of the ''thermodynamic state'' of the system.
To define a configuration (denoted as <math> \left. X \right. </math> ) we usually require:
*The position coordinates of the particles
*Depending on the problem, other variables like volume, number of particles, etc.
The probability of a given configuration, denoted as <math> \Pi \left( X | k \right) </math>,
depends on the parameters <math> k </math> (e.g. [[temperature]], [[pressure]])

Example:

:<math> \Pi_{NVT}(X|T) \propto \exp \left[ - \frac{ U (X) }{k_B T} \right] </math>

In most of the cases <math> \Pi \left( X | k \right) </math> exhibits the following features:
* It is a function of many variables
* Only for a very small fraction of the configurational space the value of <math> \Pi \left( X | k \right) </math> is not negligible.
Due to these properties, Metropolis Monte Carlo requires the use of '''Importance Sampling''' techniques

== Importance sampling ==

Importance sampling is useful to evaluate average values given by:

: <math> \langle A(X|k) \rangle = \int dX \Pi(X|k) A(X) </math>

where:
* <math> \left. X \right. </math> represents a set of many variables,
* <math> \left. \Pi \right. </math> is a probability distribution function which depends on <math> X </math> and on the constraints (parameters) <math> k </math>
* <math> \left. A \right. </math> is an observable which depends on the <math> X </math>
Depending on the behavior of <math> \left. \Pi \right. </math> we can use to compute <math> \langle A(X|k) \rangle </math> different numerical methods:
* If <math> \left. \Pi \right. </math> is, roughly speaking, quite uniform: [[Monte Carlo Integration]] methods can be effective
* If <math> \left. \Pi \right. </math> has significant values only for a small part of the configurational space, Importance sampling could be the appropriate technique
====Outline of the Method====
* Random walk over <math> \left. X \right. </math>:

: <math> \left. X_{i+1}^{test} = X_{i} + \delta X \right. </math>

From the configuration at the i-th step one builds up a ''test'' configuration by slightly modifying some of the variables <math> X </math>
* The test configuration is accepted as the new (i+1)-th configuration with certain criteria (which depends basically on <math> \Pi </math>)
* If the test configuration is not accepted as the new configuration then: <math> \left. X_{i+1} = X_i \right. </math>
The procedure is based on the [[Markov chain]] formalism, and on the [[Perron-Frobenius theorem]].
The acceptance criteria must be chosen to guarantee that after a certain equilibration ''time'' a given configuration appears with probability given by <math> \Pi(X|k) </math>

== Temperature ==

The [[temperature]] is usually fixed in Metropolis Monte Carlo simulations, since in classical statistics the kinetic degrees of freedom (momenta) can be generally, integrated out. However, it is possible to design procedures to perform Metropolis Monte Carlo simulations in the [[Microcanonical ensemble| microcanonical ensemble]] (NVE).

See [[Monte Carlo in the microcanonical ensemble]]

== Boundary Conditions ==
The simulation of homogeneous systems is usually carried out using [[periodic boundary conditions]].

== Initial configuration ==

The usual choices for the initial configuration in fluid simulations are:

* an equilibrated configuration under similar conditions (for example see <ref>[http://dx.doi.org/10.1016/j.cpc.2005.02.006 Carl McBride, Carlos Vega and Eduardo Sanz "Non-Markovian melting: a novel procedure to generate initial liquid like phases for small molecules for use in computer simulation studies", Computer Physics Communications '''170''' pp. 137-143 (2005)]</ref>)

* an ordered lattice structure. For details concerning the construction of such structures see: [[Lattice Structures | lattice structures]].
== Advanced techniques ==
:''Main article: [[Monte Carlo]]''

== References ==
<references/>
'''Related reading'''
*[http://links.jstor.org/sici?sici=0162-1459%28194909%2944%3A247%3C335%3ATMCM%3E2.0.CO%3B2-3 Nicholas Metropolis and S. Ulam "The Monte Carlo Method", Journal of the American Statistical Association '''44''' pp. 335-341 (1949)]
*[http://library.lanl.gov/cgi-bin/getfile?00326886.pdf Herbert L. Anderson "Metropolis, Monte Carlo, and the MANIAC", Los Alamos Science '''14''' pp. 96-107 (1986)]
*[http://library.lanl.gov/cgi-bin/getfile?00326866.pdf N. Metropolis "The Beginnning of the Monte Carlo Method" Los Alamos Science '''15''' pp. 125-130 (1987)]
*[http://cise.aip.org/getpdf/servlet/GetPDFServlet?filetype=pdf&id=CSENFA000002000001000065000001&idtype=cvips&prog=normal Isabel Beichl and Francis Sullivan "The Metropolis Algorithm", Computing in Science & Engineering '''2''' Issue 1 pp. 65-69 (2000)]
*[http://dx.doi.org/10.1063/1.1632112 Marshall N. Rosenbluth "Genesis of the Monte Carlo Algorithm for Statistical Mechanics", AIP Conference Proceedings '''690''' pp. 22-30 (2003)]
*[http://dx.doi.org/10.1063/1.1632114 Marshall N. Rosenbluth "Proof of Validity of Monte Carlo Method for Canonical Averaging", AIP Conference Proceedings '''690''' pp. 32-38 (2003)]
*[http://dx.doi.org/10.1063/1.1632115 William W. Wood "A Brief History of the Use of the Metropolis Method at LANL in the 1950s", AIP Conference Proceedings '''690''' pp. 39-44 (2003)]
*[http://dx.doi.org/10.1063/1.1632124 David P. Landau "The Metropolis Monte Carlo Method in Statistical Physics", AIP Conference Proceedings '''690''' pp. 134-146 (2003)]

[[Category: Monte Carlo]]