See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/232722754

MOLECULAR DYNAMICS SIMULATION OF DENSE GASES

Conference Paper · June 2011

CITATIONS READS

0 308

3 authors, including:

Nishanth Dongari Yong-hao Zhang

University of Strathclyde The University of Edinburgh
51 PUBLICATIONS 911 CITATIONS 207 PUBLICATIONS 6,269 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Nishanth Dongari on 01 June 2014.

The user has requested enhancement of the downloaded file.

 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 1

GASMEMS2011-20

MOLECULAR DYNAMICS SIMULATION OF DENSE

GASES
N. Dongari∗ , Y. H. Zhang and J. M. Reese
Department of Mechanical Engineering, University of Strathclyde, Glasgow G1 1XJ, UK.
∗
e-mail: nishanth.dongari@strath.ac.uk

Abstract
In this work, we investigate the dense gas behaviour by measuring microscopic gas properties, using
the molecular dynamics (MD) simulation method. The probability of finding a given number of
molecules simultaneously in a molecule’s intermolecular potential, the time that molecules spend
during the inter-molecular collision process and the ratio of potential to kinetic energy have been
measured as a function of density of the gas. Our numerical measurements cover a wide range of
fluid densities from S.T.P. conditions to compressible gas to a density hundred times as great, for
a cubic geometry with all periodic boundaries. The results show that the ideal gas equation of
state becomes ineffective and multiple molecule collisions increase when compared to the binary
collisions, with increase in density. This has direct implications for the limit of the dilute gas
approximation employed when deriving the Boltzmann equation. MD results also compared with
the hard-sphere based Chapman-Enskog equation for dense gases and significant deviations noticed
beyond moderately high density gases.

1 Introduction
Gas flows in narrow conduits have been receiving increased attention over the last decade with
the rapid development of microfabrication techniques. In the literature, many studies assume
that gas flows in microchannels are dynamically similar to rarefied gas flows in macro conduits
at high altitudes, if both flows have the same Knudsen number, defined as Kn = λ/h, with λ
and h representing the mean free path (MFP) and the characteristic length scale of the system,
respectively (Gad-el-hak 1999). However, this basic assumption for gas flows in nanochannels may
be questionable as discussed in the following. In a general system, two other characteristic lengths
should also be considered: the molecular diameter σ and the mean molecular spacing δ. The ratio
of these two lengths is a criterion to judge whether the gas is dilute or dense. For instance, it has
been shown that when σ/δ < 1/7, the gas can be regarded as sufficiently dilute and its dynamics
can be described by the Boltzmann equation while when σ/δ > 1/7, the gas is so dense that
imperfect gas behaviour cannot be ignored a more general kinetic equation is needed (Chapman
and Cowling 1970).
For rarefied gas flows at high altitudes, since the gas density is rather low, the mean free path
can increase up to being comparable to or even higher than the macroscopic length scale, and the
mean molecular spacing becomes much larger than the molecular diameter δ À σ. However, this is
not always true for gas flows in nanochannels, as the mean free path and mean molecular spacing
can still be very small, (l, σ) ∼ σ. This indicates that the gas in micro/nanochannles can be
rather dense and the treatment of rarefied gas dynamics may not be valid under such conditions.
Therefore, the study of dense gas dynamics is imperative and the limitations of classical kinetic
theory need to be explored.
In the rarefied gases, the transfer of molecular properties has been regarded as solely due to the
free motions of molecules between collisions; there is also, however, a transfer at collisions, over the
distance separating the centres of the two colliding molecules, during the brief time of encounter,
i.e. instantaneous transfer of energy and momentum. In the classical rigid spherical model, it is
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 2

assumed that collisions are instantaneous, and multiple encounters can be neglected. However, in
a gas at high density a molecule in the field of force of others during a large part of its motion, and
multiple encounters can be significant. In addition, molecules also spend considerable time during
the inter-molecular collisions and they may not be instantaneous any more. An adequate discussion
of dense gases for general molecular models must allow for multiple and repeated encounters, and
also for the finite number of molecules which at any instant are actually undergoing encounters.
In the case of hard-spheres, the Enskog equation (Chapman and Cowling 1970), with a simple but
essential modification, has been studied for moderately high density gases. However, a comparable
theory for realistic continuous potential functions is still lacking. Experimental investigations show
that significant deviations exist when compared with simulation results for hard-spheres (Egelstaff
et al. 1982). It seems reasonable to attribute this to differences in the dynamics of inter-molecular
collisions; however, one cannot be certain of this interpretation without full calculation using a
realistic potential.
There are great difficulties in investigating gases experimentally, although bulk properties can
be measured reliably. Interpretation of these integral flow parameters, however, is problematic
because there can be multiple causes of deviation of these quantities from those expected from
conventional macroscopic theory; assigning their relative contributions is difficult. Fitting of
macroscopic theories to experimental data with the aid of tuning parameters is unlikely to be
physically realistic.
In this situation, particle based simulation methods have come to the forefront, such as molecular
dynamics (MD) and direct simulation Monte Carlo (DSMC). To test this, fundamental properties
of the gas – ideal gas law, multiple molecular collisions relative to binary collisions, probability of
finding given number of molecules simultaneously exchanging momentum, time spent by a molecule
during the collision etc., need to be assessed over a range of gas densities. These properties can
be measured directly in molecular dynamics (MD) simulations of a moving gas in unconfined and
confined spaces (Rapaport 2004). MD is an appropriate method, compared to DSMC (Bird 1994),
to estimate these fundamental properties as it is the only deterministic approach, allowing for
realistic molecular behaviour, i.e. molecular attractions, repulsions, movements and scatterings.
The major objective of this paper is to assess the theoretical basis for the classical description
of dense gas dynamics with the aid of molecular dynamics simulation data. Such numerical
experiments are useful for assessing the range of validity of classical kinetic theories and their
application to model gas flows in miniature devices.

2 Methodology
2.1 Molecular dynamics modelling
We use the open source C++ CFD toolbox OpenFOAM (Open Field Operation and Manipu-
lation), which includes the MD solver mdF oam implemented by Macpherson et al. (2007) and
Macpherson and Reese (2008), although modified as we describe below. The major features of
the mdF oam solver comprise of the capability to perform both steady and transient MD simula-
tions for multi-species conditions, to model arbitrary 3D geometries using unstructured polyhedral
meshes, and unlimited parallel processing. Monatomic Lennard-Jones molecules are simulated,
such that the short-range repulsive/long-range attractive potential between two molecules sepa-
rated by a distance r is: ·³ ´ ¸
r −12 ³ r ´−6
VLJ (r) = 4² − , (1)
σ σ
where ² is related to the interaction strength of the molecules and σ is the hard-sphere diameter of
the gas molecules (see also section below). We investigate Argon gas, which has σ = 3.404 × 10−10
m, ² = 17.2857 × 10−22 J, and a molar mass m = 39.948 kg/kmol (Bird 1994). The molecules are
spatially distributed in a domain of interest with a random initial Gaussian velocity distribution
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 3

corresponding to a prescribed gas temperature (Kennard 1938). They are then allowed to relax
to their steady-state configuration before we take measurements.
All interactions between nonbonded particles should be calculated but, since the LJ potential
vanishes at larger r, only the interactions between molecules within a certain cutoff radius rc need
to be evaluated in the MD simulation. We choose a cut-off radius of rc = 3.0σ, as any larger
value has a negligible effect on our results but increases the computational cost of the simulations.
The equations of molecular motion are integrated using a leapfrog scheme with a step size of
10−4 τ , where τ = (mσ 2 /²) (∼ 10−12 s). Conventionally a time step of 10−3 τ is sufficient to achieve
accurate results for gases (Rapaport 2004), however we have implemented lower time step value
so that there is a negligible error in the tracking distance at the collision interface between two
molecules. The spatially-resolved parameters are written at intervals of 10τ , and the total run
time of the simulations to relax to steady-state is at least 3000τ .

2.2 Collision handling algorithm

It is well-known from the kinetic theory of gases that viscosity and thermal conductivity can be
interpreted in terms of the collisions of gas molecules, and of the free paths of the molecules be-
tween collisions. Hence, it is imperative to carry out collision measurements accurately. A collision
is defined by the magnitude of the distance R between a molecule i, and the closest molecule j,
being smaller than the target collision distance rcol . At each time-step a molecule can be at four
different parts of a collision:
(a). The molecule just collided between this time-step and the previous one. Conditions: if
(R < rcol , and collisionStatus = f alse), measure final set of properties and set collisionStatus
of molecule to true and timecol [id] = 0, where timecol stores the value of time spent by a molecule
id during the inter-molecular collision with an other molecule.
(b). The molecule is still in collision. Conditions: if (R < rcol , and collisionStatus = true),
timecol = +∆tm , where ∆tm is the MD time step.
(c). The molecule left the collision. Conditions: if (R ≥ rcol , and collisionStatus = true), set
collisionStatus to f alse.
(d). The molecule is in free flight. Conditions: if (R ≥ rcol , and collisionStatus = f alse),
measure and append properties.
R and collisionStatus are two properties defined in the molecule class. R is updated by the inter-
molecular force calculation part of the existing mdF oam code, while collisionStatus is updated
by this measurement model.
The way velocities are calculated in the Leap-Frog algorithm, might given an error if one is to
calculate the travelled distance of a molecule in the previous time-step from velocities. If measured
at the end of the time-step, (in which a second update of velocity takes place after the molecule
tracking step), measurement of displacement from the velocity of the molecules should be properly
given as:
∆r = mag(vi (t) − 0.5∆tm ai (t))∆tm , (2)
where vi (t) is the velocity of molecule i at the end of the time-step, ai (t) is its acceleration
calculated in a sub-step of the time-step.
The current method records the molecular free paths in every time step, and its major advantage
is that we can calculate the probability that molecules travel a distance without experiencing
a collision, i.e. the free path distribution function (Dongari et al. 2011). The average of these
free paths of all simulated molecules is the mean free path. Additionally, the collision handling
algorithm allows the evaluation of multiple collisions, i.e. collision events which involve three or
more molecules simultaneously and the time spent during the collision process. This means the
dilute gas assumption employed in deriving the Boltzmann equation can be explored.
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 4

6
10
4

1 2

0
80 800

0.1

MD Data
0.01
Theory

0.001
1 10 * 100

Figure 1: Variation of β as a function of normalized density and MD results comparison with the
Chapman-Enskog theory. β quantifies the extent of deviation from ideal gas equation of state.

2.3 The collision radius

We need to identify the occurrence of inter-molecular collisions in the MD simulations. While
the instance of a collision is unambiguous only if the molecular interaction is a step potential
(i.e. hard-sphere molecules), the time step required for simulations in that case would need to be
vanishingly small in order to capture the collision instant reliably . For the continuous and more
physically-realistic LJ potential that we use, we instead declare a collision to have occurred if two
molecules are closer together than a distance rcol . Molecules in the simulation start to exchange
momentum as soon as they are within the MD cut-off distance rc , and they start to repel each
other when they are closer than σ. Hence, the closeness parameter rcol must lie between 3.0σ and
σ.
An appropriate value for rcol can be evaluated through MD simulations √ to obtain the theoretical
2
mean free path for unconfined Argon gas — λk = m/(NAv πρσ 2) (Cercignani 2000), where
NAv is Avogadros number, 6.0221415 × 1023 and ρ is gas density. At STP, λAr is 6.4 × 10−8 m .
Using rcol = σ, the MD simulation retrieves this mean free path value within 1% (Dongari et
al. 2011). As a sensitivity check, if we use rcol = 2.0σ, the MD simulation result is 75% less than
the theoretical mean free path value. So the appropriate value for rcol is σ, and all the simulations
in this paper have been carried out using this value.

3 Results and Discussion

We simulate a cubic spatial geometry with periodic boundary conditions, i.e. if a molecule leaves
the domain, it is reintroduced on the opposite side of the domain, retaining the same velocity.
The side-length of this simulation domain is chosen to be the MFP of Argon at STP conditions,
and the domain is filled with 51505 molecules to achieve a gas density of 1.79 kg/m3 . Number of
molecules is kept constant for all simulations of various gas density conditions by appropriately
changing the side-length of the cube and the temperature is maintained at 298 K for all simulations
using the Brensden temperature controller (Borg et al. 2010). For the initial validation of our
MD results for dense gas conditions, we have compared the numerical measurements with the
Chapman-Enskog equation for dense gases. As the mean molecular spacing of a gas is reduced
and becomes comparable to its molecular diameter, the gas system is dense and the ideal gas
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 5



0.2

PE ⁄ KE

0.02

0.002
1 10 * 100

Figure 2: Variation of ratio of potential energy to kinetic energy as a function of normalized

density.

equation of state becomes insufficient. Alternatively, through the Chapman-Enskog procedure,

it can be demonstrated that the dense gas equation of state can be expressed as (Chapman and
Cowling 1970):
p = ρRT (1 + bρχ ), (3)
|{z}
β

with p, ρ, T and R representing the pressure, density, temperature and the gas constant, re-
spectively. χ is the pair correlation function and can be expressed as a polynomial of the fluid
density
5
χ = 1 + bρ + 0.2869(bρ)2 + 0.1103(bρ)3 + 0.0386(bρ)4 , (4)
8
with b = 2πσ 3 /(3m), where m is the mass of a molecule. MD results for β are compared with
the Eq. (4) in Fig. (1), as a function of normalized density ρ∗ , here density is normalized with
the atmospheric density value of Argon gas–1.79 kg/m3 . MD results are in good agreement with
the theoretical predictions up to ρ∗ ∼ 20 and further they qualitatively follow the theoretical
curve up to ρ∗ ∼ 100, with minor deviations. Beyond the value of ρ∗ ∼ 100, Chapman-Enskog
equation significantly deviate from the numerical measurements, see inset figure in Fig. (1). The
comparison is in agreement with the experimental findings that the hard-sphere based dense gas
equation is accurate only up to moderately dense gases (Egelstaff et al. 1982). In the present
work, we have investigated in detail to illustrate the limitations of Chapman-Enskog dense gas
theory by measuring microscopic properties, as discussed below.
The rigid hard-spherical model for gases has the mathematical advantage in the connection that
energy and momentum are merely represented in their kinetic forms. Hence, pre- and post-collision
processes between two molecules can be represented with their velocities, respectively. However,
when a realistic molecular interaction case such as LJ gas is considered, during the collision
process molecules can gain potential energy due to their inter-molecular interaction. The hard-
sphere model can be equivalent to the realistic molecular interaction case, only when the potential
energy is negligible when compared to the kinetic energy of molecules. This condition can hold
for dilute gases where the rate of collisions, i.e. intermolecular interactions are comparatively low
w.r.t molecular time scale (τ ). With increase in fluid density the inter-molecular interactions can
become quite frequent and the potential energy may significantly become comparable with the
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 6

1 1.E+0 0.8 1.E+0

*= 2.023
*= 0.921 1.E-2
1.E-2
0.8
1.E-4
0.6
Probability

Probability
1.E-4
1.E-6
0.6
1.E-8 1.E-6
0.4
0 1 2 3 4 5 0 1 2 3 4 5
0.4






0.2
0.2

0 0
0 1 2 3 4 5 0 1 2 3 4 5
Number of Molecules (a) Number of Molecules (b)

0.6 1.E+0 0.4

*= 4.975 MD Data
1.E-1

1.E-2 0.3
*= 9.976
Probability

Probability

0.4
1.E-3

1.E-4 0.2
0 1 2 3 4 5

0.2


0.1

0 0
0 1 2 3 4 5 0 1 2 3 4 5
Number of Molecules (c) Number of Molecules (d)

0.3


 0.2
*= 52.622



*= 18.876
0.15
Probability

Probability

0.2

0.1

0.1
0.05

0 0
0 1 2 3 4 5 0 1 2 3 4 5
Number of Molecules (e) Number of Molecules (f)

Figure 3: Probability of finding a given number of other molecules simultaneously in a molecule’s

intermolecular potential - LJ Argon gas, 298K. Probability variation is presented at different
densities and logarithmic plots are shown in the inset figures of (a), (b) and (c).
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 7

1
Finding zero molecules in rcut

0.8
Probability
0.6

0.4



0.2

0
0.9 9 * 90

Figure 4: Variation of probability of finding no molecule simultaneously in a molecule’s inter-

molecular potential as a function of normalized density.

1.E+0



Probability

1.E-3

1.E-6

Finding more then 6 molecules in rcut

1.E-9
0.9 9 * 90

Figure 5: Variation of probability of finding more than equivalent to 6 molecules simultaneously

in a molecule’s intermolecular potential, as a function of normalized density.
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 8

1.E+0 1.E+0 1.E+0

*= 4.795 *= 9.976
1.E-1 *= 0.921 1.E-1 1.E-1
Probability

Probability

Probability
1.E-2 1.E-2 1.E-2

1.E-3 1.E-3 1.E-3

1.E-4 1.E-4 1.E-4

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
(a) (b) (c)

1.E+0 1.E+0 1.E+0

*= 18.876 *= 52.623 *= 94.51
1.E-1 1.E-1 1.E-1
Probaiblity
Probability

Probability
1.E-2 1.E-2 1.E-2

1.E-3 1.E-3 1.E-3

1.E-4 1.E-4 1.E-4

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35
* (d) (e) (f)

Figure 6: Probability distribution of time spent by molecule in collision at various density condi-
tions - LJ Argon gas, 298K.


Percentage of Multi Collisions

 



  


1 10 * 100 1000

Figure 7: Variation of relative percentage of multiple particle collisions w.r.t. binary collisions, as
a function of normalized density.
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 9

kinetic energy.
Fig. (2) shows the ratio of the potential energy (PE) to the kinetic energy (KE) as a function
of normalized density. PE is negligible up to ρ∗ ∼ 10 and beyond this value it is comparable to
KE. Increment in PE suggests that the inter-molecular interactions are quite frequent and free
flight assumption of gas molecules may breakdown for moderately high dense gases. PE increases
linearly with increase in density and around ρ∗ ∼ 500 the ratio (PE/KE) become more than unity.
This can be attributed to the transition from gaseous phase to vapor phase and with further
increase in density up to the liquid state density value, the ratio asymptotically reaches a constant
value.
Figure (3) demonstrates the probability of finding a given number of other molecules simultane-
ously within a molecule’s intermolecular potential – rcut , at different gas densities. Number of
molecules within rcut of every molecule are measured at every time step and the probabilities are
obtained using 10 million samples. These measurements are imperative to asses the momentum
and energy transfer between molecules from the microscopic point of view. Fig. (3a) shows the
distribution at a density value close to the atmospheric condition and probability of finding zero
molecules is around 0.89 and in this case the chance of a certain molecule interaction with other
molecules is very low during its motion. This suggests that the transfer of molecular properties
can approximately be regarded as due solely to the free motions of molecules between collisions,
as employed in the classical kinetic theory, i.e. transport properties such as viscosity and thermal
conductivity can linearly be expressed in terms of mean free path, molecular velocity and density.
As the gas density increases, chances of finding multiple number of molecules within rcut is sig-
nificant and they start to dominate the probability of finding no zero molecules even at slightly
dense gas (ρ∗ ∼ 10, see Fig. (3d)). In consequence an additional mechanism for the transfer of
momentum and energy becomes important, which is negligible at ordinary densities.
Figures (4) and (5) show the probability of finding zero molecules and more than 6 molecules
within the rcut , respectively, as a function of normalised gas density. One can notice that finding
zero molecules is completely negligible beyond ρ∗ ∼ 50 and probability of finding more than 6
molecules is close to unity value from ρ∗ ∼ 90. Hence, beyond ρ∗ ∼ 90, finite number of molecules
which at any instant are actually undergoing through inter-molecular interactions and continuously
exchanging momentum and energy with neighboring molecules. At such conditions existing kinetic
theory definitions for molecular transfer properties, such as free path or mean free path may not
be valid, like for liquids.
Figure (6) shows the probability distribution for the time that molecules spend during the inter-
molecular collision process. For a certain molecule, the measurement of collision time starts when
the distance between the nearest molecule is less than rcol i.e. σ and the time ends when the
distance between the same molecule is more than rcol . At the standard density (Fig.6a), majority
of molecules spend lesser times (τ ∗ < 0.1) during the collision and the data has a very short
tail. As the density increases, the data starts to exhibit the long tail nature; the probability of
spending a long period in collision is relatively high and beyond ρ∗ ∼ 50 (Fig.6e), more than
90 % of molecules have longer collision periods (τ ∗ > 0.1). Long tail of the data is interpreted as
molecules forming temporary pairings, due to the closed packing of molecules at high densities.
Another of the assumptions underpinning the Boltzmann equation that of molecular chaos, where
the motion of individual particles is uncorrelated, something which is not applicable to molecules
which spend significant periods interacting.
Figure (7) depicts the relative percentage of multiple collisions w.r.t binary collisions, as a func-
tion of normalized density. Multiple collisions are evaluated by measuring number of molecules are
simultaneously within the rcol of a certain molecule. Employment of binary collisions assumption
for dilute gases is very imperative to derive the kinetic expression for mean free path and con-
sequently derivation of transport properties (Cercignani 2000). MD results convey that multiple
collisions are significant from ρ∗ ∼ 10 onwards.
 Proceedings of the 3rd GASMEMS Workshop-Bertinoro, June 9-11, 2011 10

4 Conclusion
View publication stats

The results of numerical molecular dynamics experiments support the argument that the molec-
ular behaviour of dense gases is better described using continuous LJ based potential functions
and the limitations of hard-sphere dense gas equation are explored. Potential energy is compara-
ble to kinetic energy at moderately high densities and explains the physical reasons behind the
limitations of ideal gas law, free flights of gas molecules and instantaneous transfer of energy and
momentum between molecules. MD measurements of time spent by molecules during their colli-
sions are distinguished by exhibiting a long tail, which has direct implications for the derivation
of Boltzmann equation based on dilute gas approximation.
The MD simulations in this paper are for a unbounded cubic geometry with periodic walls, and
for isothermal gases with no flow. It is important to explore any new description of the molecular
behaviour of dense gases through simulations of dense gases in complex nano scale geometries,
with specular-diffusive and explicit walls, and also for non-isothermal cases.

Acknowledgments
The research leading to these results has received funding from the European Community’s Seventh
Framework Programme FP7/2007-2013 under grant agreement ITN GASMEMS no 215504.

References
Bird, G. A., 1994, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford
University Press, New York.
Borg, M. K., Macpherson, G. B. and Reese, J. M., 2010, Controllers for imposing continuum-to-
molecular boundary conditions in arbitrary fluid flow geometries. Molecular Simulation, 36(10):745-
757.
Cercignani, C., 2000, Rarefied gas dynamics: from basic concepts to actual calculations. Cam-
bridge University Press, New York.
Chapman, S. and Cowling, T. G., 1970, The Mathematical Theory of Non-Uniform Gases,
Cambridge University Press, New York.
Dongari, N., Zhang, Y. H. and Reese, J. M., 2011, Molecular free path distribution in rarefied
gases. Journal of Physics D: Applied Physics, 44, 125502.
Egelstaff, P. A., Glaser, W., Litchinsky, D., Schneider, E. and Suck, J. B., 1982, Two-body
time correlations in (and the structure of) dense krypton gas. Physical Review A 27, 11061115.
Gad-el-Hak, M., 1999, The fluid mechanics of microdevices, J. Fluids Eng., 121, 5-33.
Kennard, E. H., 1938, Kinetic theory of gases with an introduction to statistical mechanics,
McGraw-Hill, New York.
Macpherson, G. B. and Reese, J.M., 2008, Molecular dynamics in arbitrary geometries: parallel
evaluation of pair forces. Molecular Simulation, 34(1), 97-115.
Macpherson, G. B., Borg, M. K. and Reese, J. M., 2007, Parallel generation of molecular
dynamics initial configurations in arbitrary geometries. Molecular Simulation, 33(15), 1199-1212.
OpenFOAM, the Open Source CFD Toolbox, http://www.opencfd.co.uk.
Rapaport, D. C., 2004, The art of molecular dynamics simulation, 2nd ed. Academic Press.