Skip to content

Commit 8624918

Browse files
simongravellesimongravelle
authored
Merge pull request #13 from mdcourse/improved-project
Implement more detail to the project
2 parents 8c79d2d + c37451a commit 8624918

File tree

1 file changed

+123
-69
lines changed

1 file changed

+123
-69
lines changed

Diff for: ‎docs/source/projects/project1.rst

+123-69
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
.. _project1-label:
22

3-
A Monte Carlo study of argon
3+
A Monte Carlo Study of Argon
44
============================
55

6-
Using Monte Carlo simulations, the equation of state of 3D fluid of argon is simulated
7-
for a large range of density.
6+
*In short --* Using Monte Carlo simulations, the equation of state (EoS) for a
7+
3D fluid of argon is simulated across a wide range of densities.
88

99
.. figure:: project1/avatar-dm.webp
1010
:alt: The fluid made of argon atoms and simulated using monte carlo and python.
@@ -18,58 +18,19 @@ for a large range of density.
1818
:align: right
1919
:class: only-light
2020

21-
In 1957, Wood and Parker reported the first numerical study of a 3D fluid that was
22-
using an attractive Lennard-Jones potential :cite:`woodMonteCarloEquation1957` (before
23-
them, Metropolis et al. used repulsive hard-sphere
24-
potentials only :cite:`metropolis1953equation, metropolis1953simulated`).
25-
In their study, Wood and Parker used a Monte Carlo algorithm to predict the
26-
Equation of State (EoS) of neutral particles interacting with Lennard-Jones potential,
27-
whose parameters were chosen to match those of argon gas. Their results show
28-
good agreement with experimental measurements of Michels on argon :cite:`michels1949isotherms`.
21+
In 1957, four years after the initial implementation of a Monte Carlo simulation by
22+
Metropolis et al. :cite:`metropolis1953equation, metropolis1953simulated`, Wood and
23+
Parker successfully implemented a 3D Monte Carlo simulation of a fluid using a full
24+
Lennard-Jones potential :cite:`woodMonteCarloEquation1957`. In their study, Wood and
25+
Parker used a Monte Carlo algorithm to predict the Equation of State (EoS) of neutral
26+
particles whose parameters were chosen to match those of argon gas. Their results
27+
showed good agreement with experimental measurements by Michels on argon :cite:`michels1949isotherms`.
28+
Here, we take advantage of our code to reproduce the results of Wood and Parker for
29+
varying density. To follow this project, only Monte Carlo moves are needed. All the
30+
chapters up to the :ref:`chapter7-label` chapter must have been completed.
2931

30-
Here, we take advantage of our code and try to reproduce the results of Wood and Parker
31-
for varying density. To follow this project, only Monte Carlo moves are needed,
32-
and all the chapters up to :ref:`chapter7-label` must have been followed.
33-
34-
Parameters choice
35-
-----------------
36-
37-
Some parameters are taken *exactly* the same as those from Wood and Parker:
38-
39-
- :math:`\sigma = 3.405~\text{Å}` (calculated as :math:`\sigma = r^* 2^{-1/6}`
40-
where :math:`r^* = 3.822~\text{Å}`),
41-
- :math:`\epsilon = 0.238~\text{kcal/mol}`,
42-
- :math:`m = 39.948~\text{g/mol}`,
43-
- :math:`T = 328.15~\text{K}` (or :math:`T = 55~^\circ\text{C}`).
44-
45-
In addition, some parameters that are not specified by Wood and Parker are
46-
freely chosen, but it should not impact the result too much, except maybe
47-
for the cut-off:
48-
49-
- :math:`d_\text{mc} = \sigma/5`,
50-
- :math:`r_\text{c} = 2.5 \sigma`.
51-
52-
Finally, since an average modern laptop is much faster than Wood and Parker's
53-
IBM 701 calculators (if you are curious, here is the |IBM-wiki| page of the IBM 701),
54-
let us choose a number of atoms that is slightly larger
55-
than theirs, i.e. :math:`N_\text{atom} = 200`.
56-
57-
.. |IBM-wiki| raw:: html
58-
59-
<a href="https://en.wikipedia.org/wiki/IBM_701" target="_blank">wikipedia</a>
60-
61-
In order to cover the same density range as Wood and Parker, let us vary the volume
62-
of the box, :math:`V`, so that :math:`V/V^* \in [0.75, 7.6]`, where
63-
:math:`V^* = 2^{-1/2} N_\text{A} r^{*3} = 23.77 \text{centimeter}^3/\text{mole}`
64-
is the molar volume.
65-
66-
Equation of State
67-
-----------------
68-
69-
The Equation of State (EoS) is a fundamental relationship in thermodynamics and
70-
statistical mechanics that describes how the state variables of a system, such as
71-
pressure :math:`p`, volume :math:`V`, and temperature :math:`T`, are interrelated.
72-
Here, let us extract the pressure of the fluid for different density values.
32+
Prepare the Python file
33+
-----------------------
7334

7435
In a Python script, let us start by importing the *constants* module of *SciPy*, and
7536
the UnitRegistry of *Pint*.
@@ -81,17 +42,17 @@ the UnitRegistry of *Pint*.
8142
ureg = UnitRegistry()
8243
ureg = UnitRegistry(autoconvert_offset_to_baseunit = True)
8344
84-
Let up also import *NumPy*, *sys*, and *multiprocessing* to launch multiple
85-
simulations in parallel.
45+
Let up also import *NumPy*, *sys*, as well as *multiprocessing* to launch
46+
multiple simulations in parallel:
8647

8748
.. code-block:: python
8849
8950
import sys
90-
import multiprocessing
51+
import multiprocessing --
9152
import numpy as np
9253
9354
Provide the full path to the code. If the code was written in the same folder
94-
as the current Python script, then *path_to_code = "./"*:
55+
as the current Python script, then simply write:
9556

9657
.. code-block:: python
9758
@@ -114,28 +75,121 @@ the right units to these variables using the UnitRegistry:
11475
Na = cst.Avogadro/ureg.mole # avogadro
11576
R = kB*Na # gas constant
11677
117-
Then, let us write a function called *launch_MC_code* that will be used to call
118-
our Monte Carlo script with a chosen value of :math:`\tau = v / v^*`. As a first
119-
step, all the required variables are being defined, and the *box_size* is calculated
120-
according to the chosen value of :math:`\tau`:
78+
Matching the Parameters by Wood and Parker
79+
------------------------------------------
80+
81+
The interaction parameters must be taken *exactly* as those from the original
82+
publication by Wood and Parker :cite:`woodMonteCarloEquation1957`. Wood and Parker
83+
provide the value for the position of the energy minimum, :math:`r^* = 3.822~\text{Å}`,
84+
which can be related to :math:`\sigma` as:
85+
86+
.. math::
87+
88+
\sigma = r^* \times 2^{-1/6} = 3.405~\text{Å}
89+
90+
Wood and Parker also specify the interaction energy for the Lennard-Jones potential,
91+
:math:`\epsilon / k_\text{B} = 119.76\,\text{K}`, from which one can calculate
92+
:math:`\epsilon = 0.238~\text{kcal/mol}` using the Boltzmann constant
93+
:math:`k_\text{B} = 1.987 \times 10^{-3}\,\text{kcal/(mol·K)}`. From the reduced
94+
temperature, :math:`k_\text{B} T/ \epsilon = 2.74`, one can calculate the temperature
95+
:math:`T = 328.15~\text{K}` (or :math:`T = 55~^\circ\text{C}`).
96+
97+
Within the Python script, write:
98+
99+
.. code-block:: python
100+
101+
epsilon = (119.76*ureg.kelvin*cst.Boltzmann*cst.N_A).to(ureg.kcal/ureg.mol)
102+
r_star = 3.822*ureg.angstrom
103+
sigma = r_star / 2**(1/6)
104+
m_argon = 39.948*ureg.gram/ureg.mol
105+
T = (55 * ureg.degC).to(ureg.degK)
106+
107+
Here, the mass of argon atoms is used. Note that the mass parameter does not
108+
matter for Monte Carlo simulations. Since an average modern laptop is much
109+
faster than Wood and Parker's IBM 701 calculators, let us use a number of
110+
particles that is larger than their own numbers of 32 or 108 particles:
111+
112+
.. code-block:: python
113+
114+
N_atom = 200
115+
116+
If you are curious, here is the |IBM-wiki| page of the IBM 701.
117+
118+
.. |IBM-wiki| raw:: html
119+
120+
<a href="https://en.wikipedia.org/wiki/IBM_701" target="_blank">wikipedia</a>
121+
122+
There are also some parameters that are not explicited by Wood and Parker that
123+
we have to choose freely, such as the maximum displacement for the Monte Carlo
124+
move, and the cut-off for the Lennard-Jones interaction:
125+
126+
.. code-block:: python
127+
128+
cut_off = sigma*2.5 # angstrom
129+
displace_mc = sigma/5 # angstrom
130+
131+
The choice of :math:`d_\text{mc}` does not impact the result, but only impacts
132+
the efficiency of the simulation. If :math:`d_\text{mc}` is too large, most move
133+
will be rejected. If :math:`d_\text{mc}` is too small, it will take a very large
134+
number of steps for the particles to explore the space. On the other hand, the
135+
choice of cutoff may impact the result. The choice of :math:`r_c = 2.5 \sigma`
136+
is relatively standard and should do the job.
137+
138+
Set the Simulation Density
139+
--------------------------
140+
141+
To cover the same density range as Wood and Parker, we will vary the volume of
142+
the box, :math:`V`, so that :math:`V/V^* \in [0.75, 7.6]`, where
143+
:math:`V^* = 2^{-1/2} N_\text{A} r^{*3} = 23.77 \, \text{cm}^3/\text{mol}` is the
144+
molar volume.
145+
146+
In the Python script, add the following:
147+
148+
.. code-block:: python
149+
150+
volume_star = r_star**3 * Na * 2**(-0.5)
151+
volume = N_atom*volume_star*tau/Na
152+
L = volume**(1/3)
153+
154+
The tau parameter (:math:`\tau = V/V^*`) will be used as the control parameter.
155+
156+
Equation of State
157+
-----------------
158+
159+
The Equation of State (EoS) is a fundamental relationship in thermodynamics and
160+
statistical mechanics that describes how the state variables of a system, such as
161+
pressure :math:`p`, volume :math:`V`, and temperature :math:`T`, are interrelated.
162+
Here, let us extract the pressure of the fluid for different density values.
163+
164+
To easily launch multiple simulation in parallel, let us create a
165+
function called *launch_MC_code* that will be used to call
166+
our Monte Carlo script with a chosen value of :math:`\tau = v / v^*`.
167+
168+
Create a new function, so that the current script ressemble the following:
121169

122170
.. code-block:: python
123171
124172
def launch_MC_code(tau):
125173
126-
epsilon = (119.76*ureg.kelvin*kB*Na).to(ureg.kcal/ureg.mol) # kcal/mol
127-
r_star = 3.822*ureg.angstrom # angstrom
128-
sigma = r_star / 2**(1/6) # angstrom
129-
N_atom = 200 # no units
174+
epsilon = (119.76*ureg.kelvin*cst.Boltzmann*cst.N_A).to(ureg.kcal/ureg.mol)
175+
r_star = 3.822*ureg.angstrom
176+
sigma = r_star / 2**(1/6)
130177
m_argon = 39.948*ureg.gram/ureg.mol
131-
T = (55 * ureg.degC).to(ureg.degK) # 55°C
132-
volume_star = r_star**3 * Na * 2**(-0.5)
133-
cut_off = sigma*2.5
178+
T = (55 * ureg.degC).to(ureg.degK)
179+
180+
N_atom = 200
181+
182+
cut_off = sigma*2.5 # angstrom
134183
displace_mc = sigma/5 # angstrom
184+
185+
volume_star = r_star**3 * Na * 2**(-0.5)
135186
volume = N_atom*volume_star*tau/Na
136187
L = volume**(1/3)
188+
137189
folder = "outputs_tau"+str(tau)+"/"
138190
191+
A folder name that will be different for every value of :math:`\tau` was also added.
192+
139193
Then, let us call the *MinimizeEnergy* class to create and pre-equilibrate the
140194
system. The use of MinimizeEnergy will help us avoid starting the Monte Carlo
141195
simulation with too much overlap between the atoms:

0 commit comments

Comments
 (0)