arXiv:1504.07665v1 [physics.comp-ph] 28 Apr 2015

Elastic properties of mono– and polydisperse two–dimensional
crystals of hard–core repulsive Yukawa particles
Jakub W. Narojczyk1 , Paweł M. Pigłowski1 , Krzysztof W. Wojciechowski2 and
Konstantin V. Tretiakov∗1
1

Institute of Molecular Physics, Polish Academy of Sciences, Smoluchowskiego 17, 60-179
Poznań, Poland
2
President St. Wojciechowski PWSZ in Kalisz, Nowy Swiat 4, 62-800 Kalisz, Poland
March 1, 2018

Abstract
Monte Carlo simulations of mono– and polydisperse two–dimensional crystals are reported.
The particles in the studied system, interacting through hard–core repulsive Yukawa potential,
form a solid phase of hexagonal lattice. The elastic properties of crystalline Yukawa systems are
determined in the N pT ensemble with variable shape of the periodic box. Effects of the Debye
screening length (κ−1 ), contact value of the potential (ǫ), and the size polydispersity of particles
on elastic properties of the system are studied. The simulations show that the polydispersity
of particles strongly influences the elastic properties of the studied system, especially on the
shear modulus. It is also found that the elastic moduli increase with density and their growth
rate depends on the screening length. Shorter screening length leads to faster increase of elastic
moduli with density and decrease of the Poisson’s ratio. In contrast to its three-dimensional
version, the studied system is non-auxetic, i.e. shows positive Poisson’s ratio.

1

Introduction

Colloidal suspensions are systems consisting of small particles of sizes from tens to hundreds nanometers dispersed in a solvent [1]. The particles can form colloidal crystals not only in three dimensions
(3D) [2] but also in two dimensions (2D) where they create a hexagonal lattice [3]. When particles are
endowed with electrical charge and placed in a dispersing medium, they exhibit exponential decay of
electric potential due to screening eﬀect from ions in the solvent [4]. This decay is described by the
screening length κ−1 (known also as the Debye length) which depends on temperature, permittivity,
and concentration of ions in the solvent. Interactions between particles in such systems are decribed
by Hard–Core Repulsive Yukawa Potential (HCRYP) [5, 6], both in two [7] and three [5, 8] dimensions. The potential consists of two parts: the hard core (reﬂecting the inﬂuence of ﬁnite sizes of
particles) and the repulsive part describing spherically symmetric Coulomb screening eﬀect. HCRYP
systems can behave (depending on κ−1 ) as hard particle systems (perfect screening) or as Wigner
crystal (no screening) [9].
More than ten years ago the liquid–solid phase diagrams for two-dimensional (2D) [7] and threedimensional (3D) [4, 5, 10, 11] HCRYP systems have been studied. The research showed that in
∗

kvt@ifmpan.poznan.pl

1

186

0.5

63

190

185

188

62.5

ν

Y∗

184

µ∗

B∗

0.496
62

183

186
184

0.492
61.5

182

182

0.488

61

181
0

0.005

0.01

0.015

0.02

1/N

0

0.005

0.01

0.015

180

0

0.02

0.005 0.01 0.015 0.02

1/N

1/N

0

0.005

0.01

0.015

0.02

1/N

Figure 1: Size dependence of the elastic moduli B (a), µ (b), Poisson’s ratio (c), and Young’s
modulus (d) for the system consisting of N = 56, 224, 504, 896, 1400 particles, at p∗ = 30, βǫ = 20
and κσ = 10. The lines represent linear ﬁts to the obtained data.
2D HCRYP particles can crystallize into a hexagonal lattice, whereas the in 3D they form f.c.c.
or b.c.c. structures (depending on external conditions). Various thermodynamic properties and
physical phenomena have been studied for the Yukawa potential, both in 2D [7, 9, 12–17] and 3D [4–
6, 10, 11, 18–27]. Recently the elastic properties of three dimensional monodisperse Yukawa system
were studied and it was shown that the system is partially auxetic [25].
The term auxetic [28] refers to the class of materials and models that exhibit unusual elastic
properties. They expand (shrink) transversally when stretched (compressed) longitudinally. Their
Poisson’s ratio [29] (in contrast to common materials) is negative. Due to their exceptional properties,
auxetic materials and models gather ever growing interest in recent years [30].
It is worth noting that usually colloidal systems are not composed of identical particles. Thus
it is meaningful to consider the enhancement of the model by introduction of the particle size dispersion [21, 22]. The study shows that the dispersion in sizes of particles can signiﬁcantly impact
the thermodynamic [31, 32] (especially elastic [33]) properties of the system. The study of polydisperse hard disks system not only indicated a strong dependence of elastic constants on particles
size polydispersity, but also revealed surprising behaviour of the Poisson’s ratio in the close–packing
limit [34]. In this context, the subject of this work are both mono– and polydisperse systems.
This work encompasses the results of study performed on two dimensional hexagonal systems
composed of Yukawa particles. The elastic properties of mono– and polydisperse two–dimensional
crystals of hard–core repulsive Yukawa particles are described. Beside studies of the inﬂuence of
screening length, the temperature and the particles size polydispersity on elastic properties, the
authors posed additional question: whether this system will also exhibit auxetic properties. It is
worth to add that the elastic properties of such a system have never been investigated.
The structure of this manuscript is the following. In section 2 basic formulas and deﬁnitions of
physical quantities necessary to describe the elasticity of the studied system, as well as the description
of models and methods are provided. In section 3 the simulation results and their discussion are
presented. The last section (4) provides the summary of the work.

2
2.1

Preliminaries
Elastic properties

of the system under non–zero external pressure p which is subjected to reversible deformation can
be described by expanding the free enthalpy (Gibbs free energy) with respect to the Lagrange strain
tensor components. The change of the Gibbs free energy corresponding the elastic deformation of
hexagonal lattice under external pressure p expressed with respect to the macroscopic elastic moduli,
2

1200

400

0.7

300

0.6

160

1000

80

600

ν

µ∗

800

B∗

p∗

120

200

0.5

400

0
0.4

0.4

100

40

βǫ = 20
βǫ = 40
βǫ = 80

0.6

0.8
∗

1

βǫ = 20
βǫ = 40
βǫ = 80

200
0
0.4

1.2

0.6

0.8
∗

1

0
0.4

1.2

βǫ = 20
βǫ = 40
βǫ = 80

0.6

ρ

ρ

0.8
∗

1

0.3
0.4

1.2

βǫ = 20
βǫ = 40
βǫ = 80

0.6

0.8
∗

1

1.2

ρ

ρ

Figure 2: Thermodynamic characteristics of the monodisprese hard–core repulsive Yukawa system
at κσ = 10. The pressure p∗ (a); bulk modulus B ∗ (b); shear modulus µ∗ (c) and the Poisson’s ratio
ν(d) are shown as functions of density for three selected values of temperatures. The lines are drawn
to guide the eye.
800

2000

1

160
κσ = 5
κσ = 10
κσ = 20

1500

80

600

κσ = 5
κσ = 10
κσ = 20

0.8

1000

ν

0.6

µ∗

B∗

p∗

120

κσ = 5
κσ = 10
κσ = 20

400

0.4

500

40
0
0.4

0.6

0.8
∗

1

1.2

200

0
0.4

0.6

0.8
∗

1

0
0.4

1.2

ρ

ρ

0.2

κσ = 5
κσ = 10
κσ = 20

0
0.6

0.8
∗

ρ

1

1.2

0.6

0.8

ρ∗

1

Figure 3: Thermodynamic characteristic of the monodisperse hard–core repulsive Yukawa system at
βǫ = 20 as functions of density: (a) pressure p∗ , (b) the bulk modulus B ∗ , (c) the shear modulus µ∗ ,
and (d) the Poisson’s ratio. The lines are drawn to guide the eye.
is of the form:
1
B(εxx + εyy )2
2
1
µ[4ε2xy + (εxx − εyy )2 ] ,
(1)
+
2
where ε is the strain tensor; Vp is the reference volume (at equilibrium) of the system under pressure
p; B is the bulk modulus, and µ is the shear modulus [35]. It follows from eq. (1) that hexagonal
lattice is elastically isotropic [29].
∆G/Vp =

2.2

The Poisson’s ratio

is the negative ratio of the relative changes of transverse to longitudinal dimensions when only the
longitudinal stress tensor component is inﬁnitesimally changed. For hexagonal structure in two
dimensions under pressure p one can express this quantity in the form [35]:
ν=

B−µ
.
B+µ

(2)

It is worth adding that for 2D systems the Poisson’s ratio can change in the range of −1 ≤ ν ≤ 1 [35].

2.3

Particle size dispersion

was introduced randomly according to the normal distribution with the ﬁrst two central moments
equal to σ and δ, respectively. The standard deviation of the particle size distribution divided by
3

the mean diameter of the particle, called the polydispersity parameter [32]:
q
1
δ=
hσ 2 i − hσi2 ,
hσi

(3)

was used to control the degree of disorder in the system. It should be noted that, since the system
consists of ﬁnite number of particles, it is important to control the values of hσi and δ as well
as possible. The simplest way to meet this condition is to repeatedly generate diﬀerent sets of
random numbers until the one is found whose properties fall, within a speciﬁed accuracy, near the
desired values [36]. In this work, however, a diﬀerent approach was used. After generating a set
of n Gaussian pseudo random numbers (using method described in [37]), the corresponding average
σ (n) and standard deviation δ(n) are used to convert this set to the standard normal distribution
Z = (X − σ (n) )/δ(n) . Next one can convert obtained values such that they exactly match any desired
distribution X ′ = (Z − σ)/δ, with the machine accuracy. This approach ensures the compatibility
of distributions for diﬀerent samples with very high accuracy.

2.4

The model description

starts with introducing reduced thermodynamic quantities used in the simulations:
reduced value
inverse temperature
number density
pressure
bulk modulus
shear modulus

symbol
βǫ
ρ∗
p∗
B∗
µ∗

therm. equiv.
ǫ/kB T
ρσ 2
βP σ 2
Bσ 2 /kB T
µσ 2 /kB T

where T is the temperature and kB is the Boltzman constant. The studied system consists of N
particles interacting via HCRYP. For the monodispersed systems, the potential of interaction is of
the form:
(
∞,
r < σ,
βu(r) =
(4)
exp[−κ(r−σ)]
βǫ
, r ≥ σ,
r/σ
where r is the distance between the centers of particles, σ is the diameter of the hard core, ǫ is
the value of pair interaction of particles at contact (i.e. at r = σ) and κ−1 is the Debye screening
length. In the case of systems with non–zero size polydispersity of particles a generalization of the
equation (4) proposed by Colombo and Dijkstra [22] was used:
(
∞,
r < σ ij ,
(5)
βu(r) =
σi σj
βǫ σr exp[−κ(r − σ ij )], r ≥ σ ij ,

2.5

Computer simulations

have been performed in the isothermal–isobaric ensemble (N pT ) by Monte Carlo method. In order to
calculate the elastic moduli, the approach based on the averaging of strain ﬂuctuations proposed by
Parrinello–Rahman [38] and further developed by Ray and Rahman [39, 40], was used. The details
of the method were described in ref. [41]. It is worth noting that in the case of a potential with
long–range interactions the minimum image method (MIM) [42] can be used. Our studies showed
that to obtain accurate results for following parameters (κ ≤ 4, β ≥ 80, rcut ≤ 2.5) of studied
potential the MIM or another method, which takes into account the long–range interactions, must
be used.
4

The typical studied systems were composed of N = 224 particles under periodic boundary conditions. Samples for simulation were chosen
√ such that the edge ratio was as close to 1 as possible.
The particles were placed in the 14a × 8 3a box, where a is the distance between the nearest
neighbors in crystal. The term nearest neighbor refers to particles that share a common edge of
their Dirichlet polygons. Results from simulations were averaged over 107 Monte Carlo cycles after prior equilibration of 106 cycles. To check the dependence of the elastic moduli on the size of
studied systems, the system of N = 56, 224, 504, 896, 1400 particles have been studied, for which
5 · 106 , 107 , 3 · 107 , 6 · 107 , 108 MC cycles were performed. The interaction potential implemented
for simulations of polydisperse systems is described by the equation (5). For monodisperse systems
the latter reduces to the form of eq. (4). Elastic moduli for each phase point were determined by
averaging over 10 independent simulations, each of which was performed for a diﬀerent structure
generated from a proper distribution.

Results and discussion

p∗

120
90

150

κσ = 5
δ
δ
δ
δ
δ

= 0.00
= 0.02
= 0.04
= 0.06
= 0.08

120
90

150

κσ = 10
δ
δ
δ
δ
δ

= 0.00
= 0.02
= 0.04
= 0.06
= 0.08

120

p∗

150

p∗

3

90

60

60

60

30

30

30

0

0
0.6

0.8

∗

1

κσ = 20
δ
δ
δ
δ
δ

= 0.00
= 0.02
= 0.04
= 0.06
= 0.08

0
0.6

0.8

ρ

∗

1

0.6

0.8

ρ∗

ρ

1

Figure 4: Equation of state of the hard–core repulsive Yukawa system at βǫ = 20 for κσ = 5 (a),
κσ = 10 (b), and κσ = 20 (c). The lines are drawn to guide the eye.

X = µ∗ B ∗
= 0.02
= 0.04
= 0.06
= 0.08

κσ = 5

1.4
1.2

1

δ
δ
δ
δ

X = µ∗ B ∗
= 0.02
= 0.04
= 0.06
= 0.08

κσ = 10

1.4
1.2

X δ /X 0

X δ /X 0

1.2

δ
δ
δ
δ

X δ /X 0

1.4

1
0.8

0.8

0.6

0.6

0.6

0.8

ρ∗

1

0.6

0.8

ρ∗

1

X = µ∗ B ∗
= 0.02
= 0.04
= 0.06
= 0.08

κσ = 20

1

0.8

0.6

δ
δ
δ
δ

0.8

ρ∗

1

Figure 5: Relative change of elastic moduli due to the non–zero size polydispersity for κσ = 5 (a),
κσ = 10 (b) and κσ = 20 (c). X stands for B or µ. The lines are drawn to guide the eye.
Based on the equation of state in Ref. [7], values of the potential’s parameters corresponding to
solid phases were chosen to allow for comparison of simulation results with experimental data, e.g.
βε ≈ 20 corresponds to colloidal systems in a low dielectric solvent at room temperature [22]. The
discussion of elastic properties of given systems that follows, is presented in the form of the three
5

1.4

1.2

1.2

µδ /µ0

B δ /B 0

1.4

1
p = 150

0.8

30

0

0.02

0.04

0.06

30

1
0.8

κσ = 5
κσ = 10
κσ = 20

0.6

p = 150

κσ = 5
κσ = 10
κσ = 20

0.6
0

0.08

δ

0.02

0.04

0.06

0.08

δ

Figure 6: Relative change of (a) the bulk modulus and (b) the shear modulus due to the non–zero
size polydispersity, studied at βǫ = 20, for two diﬀerent values of pressure. The lines are drawn to
guide the eye.
cases: the constant screening length (κσ)−1 (sec. 3.1); the constant temperature (βǫ)−1 (sec. 3.2);
and the inﬂuence of particle size polydispersity (sec. 3.3) on the mentioned above properties.
In Figure 1 are shown dependencies of the elastic moduli and Poisson’s ratio with respect to the
size of the system. It may be noted that the simulations performed for N = 224 yield results that
diﬀer only a few percent from the values extrapolated to the thermodynamic limit (N → ∞). One
can add that the agreement is quite good for all the studied system sizes. In view of the above, all
further discussion applies to systems containing N = 224 particles.

3.1

The case of constant screening length

is illustrated in Figure 2. The elastic moduli as functions of density for (κσ)−1 = 10−1 are presented
there. One can observe a monotonic increase of both elastic moduli and pressure with increasing
density. A decrease of the temperature results in an increase of pressure and elastic moduli at entire
range of studied densities. Furthermore, it causes a shift of the entire isotherm describing the lower
temperature in the direction of high densities.
The density dependence of Poisson’s ratio is presented in Fig. 2d. It is immediately apparent
that for selected κ values of Poisson’s ratio for diﬀerent temperatures are very close one to another
at high densities. This suggests that the inﬂuence of particles’ hard cores becomes dominant at
these densities. At low densities an increase of temperature causes an increase of the Poisson’s ratio,
however the dependence on temperature is rather weak.

3.2

The case of constant temperature

was investigated for (βǫ)−1 = 20−1 . Figure 3 illustrates the density dependence of pressure (a)
as well as the bulk modulus (b), the shear modulus (c) and the Poisson’s ratio (d). Similarly to
previous case, an increase of the pressure and the elastic moduli with the density can be observed.
However, the growth rate of the elastic moduli is diﬀerent and it clearly depends on the screening
length. An increase of the elastic moduli with increasing density is slower for long–range potential
((κσ)−1 = 5−1 ) than for interactions with shorter screening length ((κσ)−1 = 20−1 ).
Figure 3d illustrates the dependence of Poisson’s ratio on density for systems with diﬀerent
screening lengths. It is noticeable that for a given κ the Poisson’s ratio weakly depends on the density
of the system. However, the screening length itself essentially inﬂuences the value of Poisson’s ratio.
This leads to conclusion that one can control the Poisson’s ratio by modifying the Debye screening

6

1

1
κσ = 5

1
κσ = 10

κσ = 20

0.6

0.6

0.6

0.4

ν

0.8

ν

0.8

ν

0.8

0.4
δ
δ
δ
δ
δ

0.2
0
0.5

0.6

0.7

0.8
∗

= 0.00
= 0.02
= 0.04
= 0.06
= 0.08

0.9

0.4
δ
δ
δ
δ
δ

0.2

1

0
0.6

0.7

0.8

ρ

∗

= 0.00
= 0.02
= 0.04
= 0.06
= 0.08

0.9

ρ

1

δ
δ
δ
δ
δ

0.2
0
0.7

0.8

0.9

ρ∗

= 0.00
= 0.02
= 0.04
= 0.06
= 0.08

1

Figure 7: Density dependence of the Poisson’s ratio of the polydisperse hard–core repulsive Yukawa
system at βǫ = 20, for three values of the screening length κσ = 5 (a), κσ = 10 (b), and κσ = 20 (c).
length (e.g. by changing the properties of the electrolyte).

3.3

Influence of size polydispersity

of the particles’ hard cores on elastic properties is discussed in this section. One can start from the
examination of the equation of state presented in the Figure 4. As in the previous case, presented
data correspond to systems at reduced temperature (βǫ)−1 = 20−1 for three values of screening
length (κσ)−1 : 5−1 , 10−1 , and 20−1 . The choice of these values is dictated by two factors. Firstly,
studies of Poisson’s ratio in case of monodisperse system (see sec. 3.1) showed that its temperature
dependence is rather week while the screening length has a signiﬁcant impact on its value. Secondly,
the present study is mainly focused on looking for mechanisms aﬀecting Poisson’s ratio. Thus the
inﬂuence of the screening length on the elastic properties and Poisson’s ratio in systems with non-zero
polydispersity of particles sizes was considered at a selected (constant) temperature. Systems with
four diﬀerent values of polydispersity parameter δ have been investigated. Its inﬂuence manifests
as the change of number density (i.e. the ratio of the number of particles per 2D ’volume’) for a
given pressure (see Fig. 4). An increase of the size polydispersity causes a decrease of the number
density at a given pressure. The decrease, however, depends on the screening length. For long–range
interactions (Fig. 4a) the impact of size dispersion is signiﬁcantly smaller and noticeable only at high
pressures. Once the screening length is getting smaller (increasing κ), rising δ enforces the system
to expand even at lower pressures.
Figure 5 illustrates changes of elastic moduli caused by the size polydispersity relative to monodisperse system for diﬀerent screening lengths. The relations are presented as functions of density. One
can notice that the inﬂuence of the size polydispersity is very diﬀerent for both the moduli. In the
case of the bulk modulus, a signiﬁcant increase can be observed only for high densities and long
interactions (Fig 5a). For lower densities the inﬂuence of the size polydispersity on the bulk modulus
is negligible. Furthermore, when the eﬀect of particles’ hard cores becomes dominant (increasing
κ), one can observe only minor impact of the non–zero polydispersity on B in the whole range of
studied densities (Fig. 5c). This eﬀect is further illustrated in the Figure 6a, where the relative bulk
modulus with respect the polydispersity parameter is presented for three screening lengths at high
and low pressures. It can be seen there that only for long–range interactions and high pressure one
can observe a signiﬁcant inﬂuence of polydispersity on the bulk modulus, and this impact quickly
decays with the decrease of the screening length.
The eﬀect of the size polydispersity on the shear modulus µ is stronger than on the bulk modulus. In Figure 5a one can see that an increase of δ decreases the shear modulus. From the other
7

hand shortening of the screening length increases an inﬂuence of the polydispersity on the shear
modulus (Fig. 5 and 6). Thus, one can conclude that for a system with given size polydispersity
of particles one can additionally control its (the system) susceptibility to the shear deformations by
changing the Debye screening length (e.g. by modifying the properties of the electrolyte). Unlike
the bulk modulus, pressure has a little impact on a decrease of µ with the increasing parameter of
polydispersity (Fig. 6b). Only in the case of shortest screening lengths (κσ = 20) the diﬀerence
between the considered properties of systems under low and high pressures is noticeable. For smaller
values of κ only an increase of the polydispersity parameter is responsible for a decrease of the shear
modulus. It is worth noting, that most signiﬁcant decrease of the shear modulus is observed for the
shortest screening length.
Figure 7 shows the density dependence of the Poisson’s ratio for diﬀerent values of polydispersity
parameter and various screening lengths. It can be seen that for long-range interactions the size
polydispersity of particles aﬀects weakly the Poisson’s ratio and it weakly depends on the density.
With shortening the screening length, the impact of polydispersity becomes more noticeable, getting
very signiﬁcant role for the short-range interactions.

4

Summary and conclusions

The inﬂuence of the temperature, screening lengths and size polydispersity on elastic properties of
the two-dimensional hexagonal structure of hard–core, repulsive Yukawa particles was investigated by
Monte Carlo simulations in the N pT ensemble with variable shape of the periodic box. The obtained
results show that inﬂuence of the screening length on the elastic properties of the studied systems
is signiﬁcant. It was also shown that size polydispersity impacts the elastic properties, especially
the shear modulus, which decreases with decreasing the size polydispersity parameter. However, an
increase of particle size dispersion is reﬂected in the increase of Poisson’s ratio. It was shown that
the 2D Yukawa model is non–auxetic but the value of the Poisson’s ratio decreases with decreasing
screening length. This mechanism can be useful in designing systems with low Poisson’s ratio.

Acknowledgement
We are grateful to Dr. Mikołaj Kowalik (Pennsylvania State University) for discussions. This work
was supported partially by the Polish National Science Center grants DEC-2012/05/B/ST3/03255.
The calculations were partially performed at Poznań Supercomputing and Networking Center (PCSS).

References
[1] A. P. Gast and W. B. Russel, Phys. Today 51, 24 (1998).
[2] W. B. Russel, D. A. Saville, and W. R. Schowalter, Colloidal Dispersions (Cambridge
University Press, Cambridge, 1991).
[3] P. Pieranski, Phys. Rev. Lett. 45, 569 (1980).
[4] S. Auer and D. Frenkel, J. Phys.: Condens. Matter 14, 7667–7680 (2002).
[5] F. E. Azhar, M. Baus, and J. P. Ryckaert, J. Chem. Phys 112, 5121–5126 (2000).
[6] A. P. Hynninen and M. Dijkstra, J. Phys.: Condens. Matter 15, 3557–3567 (2003).

8

[7] R. Asgari, B. Davoudi, and B. Tanatar, Phys. Rev. E 64, 41406 (2001).
[8] A. P. Hynninen and M. Dijkstra, Phys. Rev. E 68, 021407 (2003).
[9] L. Assoud, R. Messina, and H. Lowen, J. Chem. Phys. 129, 164511 (2008).
[10] B. Davoudi, M. Kohandel, M. Mohammadi, and B. Tanatar, Phys. Rev. E 62, 6977
(2000).
[11] N. M. Dixit and C. F. Zukoski, J. Phys.: Condens. Matter 15, 1531–1552 (2003).
[12] T. Glanz and H. Lowen, J. Phys.: Condens. Matter 24, 464114 (2012).
[13] P. Hartmann, G. J. Kalman, Z. Donko, and K. Kutasi, Phys. Rev. E 72, 26409 (2005).
[14] Y. V. Khrustalyov and O. S. Vaulina, Phys. Rev. E 85, 46405 (2012).
[15] O. S. Vaulina and X. G. Koss, Phys. Lett. A 378, 719–722 (2014).
[16] O. S. Vaulina and X. G. Koss, Phys. Lett. A 378, 3475–3479 (2014).
[17] A. Radzvilavicius, Phys. Rev. E 85, 51111 (2012).
[18] M. Heinen, P. Holmqvist, A. J. Banchio, and G. Nagele, J. Chem. Phys 134, 44532
(2011).
[19] S. Zhou and X. Zhang, J. Phys. Chem. B 107, 5294–5299 (2003).
[20] K. Krerner, M. O. Robbins, and G. S. Grest, Phys. Rev. Lett. 57, 2694–2697 (1986).
[21] M. N. van der Linden, A. van Blaaderen, and M. Dijkstra, J. Chem. Phys. 138, 114903
(2013).
[22] J. Colombo and M. Dijkstra, J. Chem. Phys 134, 154504 (2011).
[23] J. Gapinski, G. Nagele, and A. Patkowski, J. Chem. Phys 136, 24507 (2012).
[24] J. Gapinski, G. Nagele, and A. Patkowski, J. Chem. Phys 141, 124505 (2014).
[25] K. V. Tretiakov and K. W. Wojciechowski, Phys. Status Solidi B 251, 383–387 (2014).
[26] G. J. Pauschenwein and G. Kahl, J. Phys.: Condens. Matter 21, 474202 (2009).
[27] E. J. Meijer and F. E. Azhar, J. Chem. Phys 106, 4678 (2007).
[28] K. E. Evans, Endeavour, New Series 15, 170 (1991).
[29] L. D. Landau and E. M. Lifshits, Theory of Elasticity, 3rd Edition (Pergamon Press, Oxford, 1993).
[30] K. L. Alderson1, A. Alderson, J. N. Grima, and K. W. Wojciechowski, Phys. Status
Solidi B 251, 263–266 (2014).
[31] S. E. Phan, W. B. Russel, Z. Cheng, J. Zhu, P. M. Chaikin, J. H. Dunsmuir, and R. H.
Ottewill, Phys. Rev. E 54, 6633–6645 (1996).
[32] S. E. Phan, W. B. Russel, J. X. Zhu, and P. M. Chaikin, J. Chem. Phys. 108, 9789 (1998).
9

[33] K. W. Wojciechowski and J. Narojczyk, Rev. Adv. Mater. Sci. 12(2), 120–126 (2006).
[34] K. V. Tretiakov and K. W. Wojciechowski, J. Chem. Phys. 136, 204506 (2012).
[35] K. W. Wojciechowski, Phys. Lett. A 137, 60–64 (1989).
[36] K. V. Tretiakov and K. W. Wojciechowski, Phys. Status Solidi B 250(10), 2020–2029
(2013).
[37] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1987).
[38] M. Parrinello and A. Rahman, J. Chem. Phys. 76, 2662 (1982).
[39] J. R. Ray and A. Rahman, J. Chem. Phys. 80, 4423 (1984).
[40] J. R. Ray and A. Rahman, J. Chem. Phys. 82, 4243 (1985).
[41] K. W. Wojciechowski, K. V. Tretiakov, and M. Kowalik, Phys. Rev. E 67, 036121
(2003).
[42] K. V. Tretiakov and K. W. Wojciechowski, Comput. Phys. Commun. 189, 77–83 (2015).

10

