Strain-Strain Fluctuations
5.3 Cubic Crystals
5.3.1 Elastic Constants of Hard-Sphere Crystals
The elastic properties of hard-sphere crystals have been analyzed in the following papers [Fre87, Lai92,Far00,Pro03,Tre05], which we use as reference for our calculations. The values are listed in appendixC2.
Strain Histograms
0.01 0.1 1 10 100
-0.1 -0.05 0 0.05 0.1
P( <ε11 - ε22>Lb )
<ε11 - ε22>Lb Lb/L = 0.16
Lb/L = 0.24 Lb/L = 0.32
Lb/L = 0.40 Lb/L = 0.48 Lb/L = 0.56
Lb/L = 0.64 Lb/L = 0.72
Figure 5.3: Representative strain histograms ofhε11−ε22iLb of a hard core system with fcc-lattice symmetry.
In figure5.3we plot the strain histograms ofhε11−ε22iLb for a hard-core crystal of fcc symmetry in a semi-log plot for a selection of relative system sizes Lb/L. Obviously, the histograms are all centered around zero, which indicates that the overall system is not strained. All the strain distributionshε11−ε22iLb perfectly match with the Gaussian distribution. This is also true for all the other strain combinations being considered. From fits of the Gaussian distribution to these strain histograms we can obtain the block size dependent strain fluctuationsSiL˜b. After additional
scaling these values withx≡Lb/Lwe obtain the scaling curves forx·SL˜ib as functions ofxwhich are described in the following.
Block Analysis of Cubic Blocks
The free energy in the harmonic form of a cubic crystal is given by equation (3.60). It contains seven linear combinations of Lagrangian strainse˜i, which are defined by equation (3.59). Conse-quently, we can calculate seven different strain-strain fluctuationshei˜ei˜ifor just three independent elastic constantsA= (C11+2C12)/3, B= (C11−C12)/3, andC44. Thus, some of the block size dependent strain correlations as a function of the relative system sizeLb/Lare expected to show an identical scaling behavior.
0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Lb/L*Si(Lb/L)
Lb/L S1
S2 S3 S4 S5 S6 S7
Figure 5.4: Block size dependent strain fluctuations of a cubic crystal with hard-core interaction which are evaluated for cubic blocks (N=4000,ρ=1.13).
In figure5.4 we plot the scaling behavior of all the seven strain correlationsLb/L·Si˜(Lb/L) ex-plicitly for a face-centered cubic crystal ofN=4000 particles and the particle densityρ=1.13.
These curves have been obtained after a simulation run of 106 MC steps. The strain correlations are evaluated for cubic blocks. Obviously, we get an identical scaling behavior for the three strain correlationsS2=S3=S4, which are connected to the elastic constantB, and the three correlations S5=S6=S7, which are connected to the elastic constantC44. This internal consistency check, which is embodied in the block analysis method, makes us confident, that we did not build in any inadvertent errors into our algorithm. The scaling behavior of all the curves in figure5.4 shows the expected limiting behavior,i.e.forLb/L→0 andLb/L→1 all curves approach zero.
It is known, that for cubic hard-core systems the elastic constantAhas the largest value, followed by the constantC44, and the constantB(cf. literature values listed in sectionC2). Therefore, in the limitLb/L→0 the slopes of strain correlationsLb/L·Si˜(Lb/L), which are connected to the inverse values of the elastic constants, are expected to reflect this order. From figure5.4it is obvious, that the slope of the linear extrapolation in the limitLb/L→0 is maximal for the three curvesS2,S3, S4, which determine the value ofB. The slope of the curve ofS1in the limitLb/L→0 is smaller than for the three curvesS5,S6,S7, which reflects that the value ofAis a little bit greater than the valueC44.
Block Analysis of Spherical Blocks
Instead of using cubic blocks we have also performed simulations with spherical blocks.4 Please note, that for spherical blocks we plot the strain correlations as a function of (Vb/V)1/3. We employ the following relations
Lb
L = Vb
V 1/3
= 43πR3b L3
!1/3
=π 6
1/3Db
L =0.805996Db
L , (5.2)
where Db denotes the diameter of the spherical block. The values ofDb belong to the interval Db∈[0,L], and thus the maximum value of(Vb/V)1/3 is 0.805996. The largest spherical block which fits into a cubic simulation box of size L contains only 52% of the total volume of the simulation box. To take the strain-strain correlations of all particles into account it is important to perform the average over many block locations. By default 100 different block positions are evaluated.
0 0.001 0.002 0.003 0.004 0.005
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
(Vb/V)1/3 *Sγγ(Vb)
(Vb/V)1/3 A = 92.3 +/- 1.0
B = 40.1 +/- 0.1 C44 = 82.7 +/- 0.7
Figure 5.5: Elastic constants for a cubic crystal with hard-core interaction evaluated for spherical blocks (N=4000,ρ =1.13). For the constantBthe evaluation of the data points is restricted the intervalx= (Vb/V)1/3= [0.1,0.7]. For the two other elastic constants the data points in the interval x= (Vb/V)1/3= [0.1,0.81] are evaluated. The given errors of the elastic constants are due to the errors of the fitting variableS∞˜i .
In figure 5.5 we show a representative plot of the scaling behavior of the block size dependent strain-strain fluctuations evaluated for a selection of spherical block sizes. The curves are quali-tative similar to the block size dependent strain correlations for cubic block, except for the curve x·SLBBb which is bend forx→xmax=0.805996.
4By default the strain-strain fluctuations are evaluated for cubic blocks, but on definition of the compiler variable SPHERICAL SUBBOXspherical blocks are used instead.
Elastic Constants from the Strain-Strain Fluctuations
From the block size dependency of the strain-strain fluctuationsSL˜jb we determine the values for the compliances S∞j˜ by fitting the finite-size scaling function to the simulation data points. For spherical blocks it reads (cf. equation (3.100) and section3.4.1)
x·SL˜ib(x) =S∞˜i x 1−e−αx(αx+1)
−γx4
. (5.3)
In the case of cubic blocks the explicit form is more complex, because we have to use Φ(α) of equation (3.98) instead of Φ(α) of equation (3.96). The three constantsS∞i˜,α ≡xL/ξ, and γ ≡Ψ(L/ξ) are used as fit parameters. The solid lines in figure 5.5 are the resulting fits of equation (5.3) to the data.
10 100 1000
0.1 1
P, BT [kBT/σ3 ]
(ρ/ρ0)-1
P P from [Far00, Sti67]
P from [Fre87]
BT BT from [Far00, Sti67]
BT from [Fre87]
10 100 1000
0.1 1
C12, C44 [kBT/σ3 ]
(ρ/ρ0)-1
C12 C12 from [Far00]
C12 from [Fre87]
C44 C44 from [Far00]
C44 from [Fre87]
Figure 5.6: Results for the elastic constants of a cubic crystal with hard-core interaction. The reference lattice is a perfect fcc lattice. Shown are the pressureP, the isothermal bulk modulus BT, and the constantsC12andC44. Also plotted are the literature values of the references [Fre87,Far00,Sti67].
In the next step the elastic constants are determined from the compliances S∞˜j by the relations (cf. equation (3.105))
βA= 1
f1SAA∞ , βB= 1
f2SBB∞ , and βC44= 1
2S44∞. (5.4)
We obtain very similar results for either cubic or spherical blocks. Figure 5.6 summarizes the simulation results and compares them to the literature values of Frenkel [Fre87,Pro03] and Kan-tor [Far00]. Our simulation results are marked by the colored symbols, and the black data points indicate the results of Frenkel [Fre87,Pro03]. The lines represent the expected values from the free-volume predictions of Stillinger and Salsburg [Sti67]. For the elastic constantsC12andC44we used the values for the constantsA1andA2as they were determined by Farago and Kantor [Far00].
We determined the pressure Pby using the expression of equation (3.125) for the calculation of the components of the stress tensorσ and the relationP=Tr(σ)/3.
Please note, that in order to obtain values, which are in good agreement with the literature values, we have to choose f1=f2=3 in the case of hard sphere systems of face-centered cubic symmetry.
These factors are determined heuristically. In the case of 2D systems the factors f1andf2are equal to two, which was explained by the argument, that the block analysis method accounts only for effects of the “inner” system, but not for the thermal bath. The heuristically determined factors f1 and f2are an indication, that the direct generalization of the block analysis method, as described in section3.4.1, has some deficiencies.
Kantor Method
For comparison we also determined the elastic constants of cubic hard-sphere crystals by the Kantor method (cf. section 3.4.3). In figure5.7 we show the result for a fcc-crystal consisting
0 0.002 0.004 0.006 0.008 0.01
εn /a 0
50 100 150 200
P, C αβ [k BT/ a3 ]
C11 = 114.3 C12 = 22.5 C44 = 73.0 P = 14.6
Figure 5.7: Elastic constants evaluated by using the Kantor method.
of N =4000 hard spheres and the volume fraction η =1.1. We plot the pressure P and the three elastic constantsC11,C12, andC44 for the volume fractionη=1.1 as a function ofεn (see the definition in equation (3.127)) which was normalized by the hard sphere radius a. The data points are obtained after averaging over 10 independent simulation runs of 2·106MC steps each.
The solid curves are the weighted least squares fits of the third order polynomials in εn to the simulation data. Extrapolation toεn=0 (particle contact value) returns the values of the pressure and the elastic constants. These constants have the following values in reduced units P=14.6, C11=114.33,C12=22.5, andC44=73.0. This is in good agreement with the literature values, as listed in AppendixC2. But from figure5.7we can directly conclude that the values for the elastic constantsC11andC12 have a high uncertainty, because there are large increases near the contact in the intervalεn= [0,0.004]a.
We conclude, that the Kantor method has a poor convergence and the computation effort is very high.
5.3.2 Elastic Constants of Lennard-Jones Crystals Squire Method
With our simulation algorithm we exactly reproduced the results of Cowley by the implementation of the fluctuation formula of Squire, Holt, and Hoover [Squ69] of equation (3.122).5 Already, simulations of very small systems with onlyN=256 Lennard-Jones particles had been sufficient to obtain the correct results. We made additional simulation runs withN=4000 particles, but we did not detect any size dependency. This shows, that the Squire method is a very efficient method to determine elastic constants in a simulation, but it implies a-priori the explicit knowledge of the differentiable particle pair interaction potential (cf. equations (3.121) and (3.122)). Thus, it is not applicable for singular potentials as for example the hard-sphere potential.
Block Analysis Method
Figure5.8 shows the scaling behavior of the block size dependent strain-strain fluctuations plus the fits of the finite-size scaling function of equation (5.3) for a face-centered cubic crystal with Lennard-Jones pair interaction.
In figure5.9we summarize the values of elastic constants obtained by our method and compare them with the results published by Cowley [Cow83]. To bring our results in accordance with the results of Cowley we have to use the factors f1=2 and f2=3 (cf. equation (5.4)). These heuristic factors f1and f2are different to those needed for hard-sphere systems.
5The executable programSEC3D, which calculates the elastic constants according to the fluctuation formula of Squire, Holt, and Hoover given in equation (3.122), can be generated by the commandmake SEC3Dand with the compiler variablesSQUIRE EC,LENNARD JONES,FCC 100, andDIMENSION=3being defined. Optionally, one can also define the compiler variableALL TENSOR ELEMENTS, where all 21 independent tensor elements ofCi jklare determined.
0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
(Vb/V)1/3 *Sγγ(Vb)
(Vb/V)1/3 A = 102.00
B = 29.20 C44 = 81.07
Figure 5.8: Elastic constants for a face-centered cubic crystal withN=4000 Lennard-Jones par-ticles at the reduced temperatureT =0.3.
0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
T [ε/kB] 0
100 200 300 400 500 600
C ik / Nk BT
C11 C12 C44
C11, Cowley C12, Cowley C44, Cowley
Figure 5.9: Comparison of the elastic constants for a Lennard-Jones crystal of a fcc-structure with the results of Cowley [Cow83].
Comment on the Particle Monte Carlo Moves
In order to check for possible systematic errors of our Metropolis Monte Carlo routine we com-pared the results for different realizations of particle selection. First, we sequentially moved each particle in all three spatial directions according to the particle number being given to it during the initialization of the crystal lattice. Later we selected from a uniform random distribution a particle number for a particle to be moved. In this case a MC step is defined as the MC moves of N particles of the crystal, which consists ofN particles.6 Consequently, in a single MC step not necessarily every particle is attempted to be moved. Additionally, we decoupled the MC moves in each coordinate direction. All three sequencing methods gave us identical results of the block size dependent strain-strain fluctuations.
System Size Dependency of the Block Size Dependent Strain Fluctuations
Longer simulation runs with up to 5·106 MC steps did result only in minor differences of the block size dependent strain fluctuations, which can be explained by better statistics.
0 0.0005 0.001 0.0015 0.002 0.0025 0.003
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Lb/L*SΑΑ(Lb)
Lb/L N=4000
N=6912 N=32000 N=62500 N=108000
Figure 5.10: Representative plot of the system size dependency ofLb/L·SLAAb for a cubic Lennard-Jones crystal atT =0.3.
The results of the block size dependent strain-strain fluctuations as a function of the relative system size Lb/Ldepend on the number of simulated particles. In figure5.10we show a representative plot for the correlationLb/L·SLAAb of a cubic Lennard-Jones crystal atT=0.3. But, this effect does not explain the discrepancies between our results and the values of the literature, which lead to the introduction of the heuristic factors f1and f2.
6The random selection of which particle is selected to move is used by default in our algorithm.
Elastic and Plastic Deformation of the Crystal
In section4.3.1we introduced (in the context of the method of Falk and Langer [Fal98] for local strain calculation) the variableD2(t), which provides a measure of the local instantaneous devi-ation at the position ri0(t) from affine deformation, i.e. from purely elastic deformation of the crystal. We calculated this variable in our simulations of a hard-core cubic crystal for a selection of particle densitiesρ. The resulting distributions ofD2(t)are shown in figure5.11.
0 0.5 1 1.5 2 2.5 3 3.5
0 0.1 0.2 0.3 0.4 0.5
number of counts [105 ]
< D2(t) >
ρ=1.12 ρ=1.13 ρ=1.15 ρ=1.17 ρ=1.25
Figure 5.11: Distribution of the variableD2(t)of equation (4.17) for a hard-core crystal with fcc symmetry. The results are obtained forN=4000 particles after a simulation run of 2·106MC steps.
Clearly, the distribution ofD2(t) has peaks at values greater than zero. With increasing particle densityρ,i.e.with increasing crystal stiffness, the distribution ofD2(t)is peaked at lower values and therefore is of thermal origin. During the simulation run we performed checks of the configu-ration symmetry for each Monte Carlo move. We rejected Monte Carlo moves of particles where a particle would leave the cell being defined by its 12 neighboring reference lattice atoms.7 Nev-ertheless, we find small spontaneous particle displacements which are not connected to an elastic deformation of the crystal. This can be deduced from the non-zero residual value ofD2(t).
The instantaneous displacement field of particle iis given by a superposition of an elastic (cf.
equation (3.1)) and a plastic (non-affine) part [Bar06,Lub02],i.e.
ui(R,t) =uelastici +uplastici . (5.5) Consequently, for the strains we have
ei(R,t) =eelastici +eplastici , (5.6) which are not correlated with each other,i.e.heelastici eplastici i=0. The whole crystal lattice responds elastically to applied uniform strains, which implies that RdRh(eplastici )2i=0. These conditions are expected to be sufficient to include the instantaneous non-affine displacements into our model.
7This configuration check is applied, when the compiler variableLATTICE STRUCTURE CHECKis defined.
Generally, for the determination of elastic constants in the vicinity of the melting transition, the non-affine displacements need to be subtracted to obtain accurate results.
All the results of the block analysis method, which we listed in this section, lead to the conclusion, that the assumption of isotropic strain-strain fluctuations, which we made in the analytic derivation of the finite-size scaling form, is not valid for three-dimensional systems.