Data Analysis

4.3.1 Numerical Computation of the Lagrangian Strains

Numerical calculation of the Lagrangian strains has been implemented by various methods into the simulation procedure. These are:

• Fourier transform method: The displacement field is Fourier transformed and the La-grangian strains are calculated in thek-space according to equation (3.25). Reverse Fourier transformation gives the Lagrangian strains in the position space.

• Octahedron method:

Figure 4.1: The nearest neighbor cluster of an fcc-crystal, which is given by a particle and its twelve near-est neighbors. One of these tetrahedrons is marked by the particles in red. The figure is taken from [Rei06f].

The displacement derivatives of a particle with indexα =1,...,N can be approximately computed using the displacements of the twelve nearest neighbor particles with indexβ = 1,...,12 from the system of linear equations (i=1,2,3)

u^β_i (t)−u^α_i (t)≈u^α_{i j}(t)·(R^β_j −R^α_j). (4.13) Remind that u^α_{i j}(t) =∂u^α_i /∂R^α_j are the displacement gradients. These can be computed from the analysis of the nearest neighbor cluster of a particle with indexα. Such a nearest neighbor cluster is illustrated in figure4.1for a fcc-crystal. It allows for the construction of eight different tetrahedrons containing the central particle with indexα. For each tetrahe-dron we can write down the equation

det







∆X ∆Y ∆Z ∆u_i(t)

∆X1 ∆Y1 ∆Z1 ∆ui,1(t)

∆X₂ ∆Y₂ ∆Z₂ ∆u_i,2(t)

∆X₃ ∆Y₃ ∆Z₃ ∆u_i,3(t)





=0 (4.14)

where the particleα has the index 0, the three other particles of the tetrahedron are labeled by the indicesn=1,2,3, andR= (X,Y,Z)used in∆Xn≡Xn−X0. Solving equation (4.14) for the unknown ∆u_i(t)≡u_i−u_i,0 and dividing by ∆X ≡X−X₀, ∆Y ≡Y−Y₀ or ∆Z≡ Z−Z₀leads to nine expressions for the nine unknown displacement gradient tensor elements ui j(t)which can be computed numerically.⁸ The Lagrangian strainsε_{i j}^α of the particleα at position R^α are obtained from the symmetric combinations of the displacement gradient tensor elements averaged over all eight possible tetrahedrons,i.e.ε_{i j}^α=¹₂

u^α_ij+u^α_ji .

8These elementsu^α_{i j}(t)are explicitly written in terms of the known distances and displacements of the three other

• Singular value decomposition method: In 3D the system of linear equations of (4.13) consists of twelve equations for nine unknown derivatives. Such an overdetermined system of linear equations can be approximately solved in the least square sense by the method of Singular Value Decomposition(SVD).⁹This method has the advantage that is can also be applied to non-perfect reference lattices where the “Octahedron method” does not work.

• Method of Falk and Langer: This method introduced by Falk and Langer [Fal98] searches for the best mapping of the nearest neighbor configuration of a given particle in the actual configuration onto the corresponding nearest neighbor configuration within the reference lattice. It computes instantaneous local displacement gradientsu_{i j} by minimizing the devi-ations from uniform deformdevi-ations.

To make this clearer, let me introduce the following definitions

∆rⁿ_i(t) ≡ rⁿ_i(t)−r⁰_i(t) (4.15)

∆Rⁿ_i ≡ Rⁿ_i −R⁰_i. (4.16)

The expression (4.15) defines the displacements within the actual configuration of the par-ticle with indexnrelative to the central one having the indexn=0. The indexnruns over all particles within the interaction range of the central particle located atr_i⁰(t). The second definition (4.16) is the analog of (4.15) within the reference lattice. The total mean-square difference between the actual relative separation∆rⁿ_i(t)and the relative separation∆Rⁿ_i that the particles would have if they were in a region of uniform strain is

D²(t) =^NNeighbor

∑

n=1

∑

d i=1

"

∆rⁿ_i(t)−

∑

^d

j=1(δi j+ui j(t))∆Rⁿ_i

#₂

, (4.17)

From definition (3.2) is is clear thatD²(t)is zero for homogeneous (affine) deformation as being demanded.

particles of the tetrahedron from the central particle of the nearest neighbor cluster:

u^α_i1(t) = (∆Y₂∆Z₃−∆Y₃∆Z₂)∆u_i,1+ (∆Y₃∆Z₁−∆Y₁∆Z₃)∆u_i,2+ (∆Y₁∆Z₂−∆Y₂∆Z₁)∆u_i,3 (∆Y₂∆Z₃−∆Y₃∆Z₂)∆X₁+ (∆Y₃∆Z₁−∆Y₁∆Z₃)∆X₂+ (∆Y₁∆Z₂−∆Y₂∆Z₁)∆X₃ u^α_i2(t) = (∆X₃∆Z₂−∆X₂∆Z₃)∆u_i,1+ (∆X₁∆Z₃−∆X₃∆Z₁)∆u_i,2+ (∆X₂∆Z₁−∆X₁∆Z₂)∆u_i,3

(∆Y₂∆Z₃−∆Y₃∆Z₂)∆X₁+ (∆Y₃∆Z₁−∆Y₁∆Z₃)∆X₂+ (∆Y₁∆Z₂−∆Y₂∆Z₁)∆X₃ u^α_i3(t) = (∆X₂∆Y₃−∆X₃∆Y₂)∆u_i,1+ (∆X₃∆Y₁−∆X₁∆Y₃)∆u_i,2+ (∆X₁∆Y₂−∆X₂∆Y₁)∆u_i,3

(∆Y₂∆Z₃−∆Y₃∆Z₂)∆X₁+ (∆Y₃∆Z₁−∆Y₁∆Z₃)∆X₂+ (∆Y₁∆Z₂−∆Y₂∆Z₁)∆X₃

9The SVD is based on the theorem of linear algebra which states: AnyM×NmatrixAwhose number of rowsMis greater or equal to its number of columnsNhas a singular value decomposition into

A_M×N=U_M×N·S_N×N·V^T_N×N.

HereU_M×Ndenotes a column-orthogonal matrix,S_N×Na diagonal matrix with non-negative singular valuesσ_i= S_ii, andV^T_N×Nis the transpose of an orthogonal matrix. Applying this theorem to the linear system

A_M×N·x_N=b_N results in the solution

x=V·S⁻¹·U·b which minimizes||A·x=b||₂in the least square sense. [Pre92]

We used the functionsint gsl linalg SV decomp(...) andint gsl linalg SV solve(...) for numerical com-putation of the SVD which are part of the GNU Scientific Library [Gal06].

The displacement gradient tensor ui j which minimizes the expression (4.17) can be found by computing

X_{i j}(t)≡^NNeighbor

∑

n=1 ∆rⁿ_i(t)∆Rⁿ_j and Y_{i j}≡^NNeighbor

∑

n=1 ∆Rⁿ_i∆Rⁿ_j (4.18) which are connected to the displacement gradient tensor by¹⁰

u_{i j}(t) =

∑

^d

k=1

X_ik(t)Y_jk⁻¹−δ_{i j}. (4.19) The instantaneous local Lagrangian strains areε_{i j}(t) =¹₂(u_ij(t) +u_ji(t))according to their definition by equation (3.8). The minimum value ofD²(t)is a measure for the local devi-ation at positionr⁰_i(t) from the affine (homogeneous) deformation, i.e. a measure for the plastic (non-local) deformation of the crystal.

4.3.2 Simulation Details of the Block Analysis Method

The block analysis described in subsection3.4.1is performed every 10 Monte Carlo steps¹¹. The essential steps of the simulation algorithm for the block analysis procedure according to Sen-guptaet al.[Sen00b] are:

(a) Calculation of the elastic strains from the sampled configurations for a given reference lattice using one of the methods described in subsection4.3.1.

(b) This is followed by theblock analysis,i.e.a coarse-graining procedure of averaging fluctu-ations of strains over larger and larger blocks of lengthL_b≤Lwithin the simulation box of sizeL. The block sizes are chosen as¹²

L_b= L M_b = n

50L, withn=1,...,50. (4.20)

10Auxiliary calculation: Finding the minimum ofD²(t)as defined in equation (4.17):

0 =^! ∂D²(t)

11This can be modified by the variableDiffMCSin the parameter fileparameter.ini.

12The size difference of the block lengthsL_b of two consecutive block sizes can be adjusted by the parameter Subbox SizeDifference, which has the default value 0.02L_b/L.

Each block of sizeL_bis placed 100 times at a random position within the simulation box.

The strainsε^L_{i j}^b averaged over all block sizes and over all 100 blocks of each block size are accumulated together with the fluctuationsε^L_{i j}^bε^L_kl^brequired for the symmetry of the reference lattice.

(c) At the end of the simulation run wecalculate the ensemble average(average over all MC steps of the production run) of the block averaged quantitiesε^L_{i j}^bε^L_kl^b to obtain the compli-ance tensorS_{i jkl}^Lb =D

ε^L_{i j}^bε^L_kl^bE

. Additionally, we calculate the superpositionsS^L_j˜^b ≡D e^L_j˜^be^L_j˜^bE (cf. subsection3.4.1) needed for the given symmetry of the reference lattice.

(d) Subsequently, in an evaluation procedure weconvert the data for the strain fluctuations into elastic constantsby fitting the scaling form (3.100) to the data and using the connection of equation (3.105) the compliances and the elastic constants.

4.3.3 Order Parameters

The following order parameters are added in the simulation algorithm. They enable us to check the instantaneous local and global lattice structure.

Orientational Order Parameter in 2D

Hexagonal order of the system can be checked by the local orientational order parameter [Nel02]

ψ6,i= 1 N_b

Nb

j=1

∑

e^{6iθi j} (4.21)

which returns a measure of the orientational order based upon the distribution of anglesθ_{i j} (mea-sured with respect to a fixed axis) of the lines joining a particleiwith its surroundingNbneighbors.

The value ofψ_6,iis invariant of system translation and rotation. In the case of a perfect hexagonal lattice the distribution of angles θ_{i j} is expected to show six peaks separated by an angle of 60^◦ each. The above defined order parameterψ6,i is constructed in such a way that it has the values one for the local hexagonal order and zero for no local orientational order. The corresponding global order parameter can be obtained by averaging over allNparticles,

Ψ6= 1

N

∑

N j=1ψ6,i

. (4.22)

Bond Orientational Order Parameter in 3D

Steinhardt et al.[Ste83] first proposed so-calledbond orientational order parametersto charac-terize geometries of atomic clusters or crystal lattices in 3D. They associate spherical harmonics Y_lm(θ,φ) with every bond joining a particle with its nearest neighbors. A bond are lines which result from the connection of a particular particle with its nearest neighbors. For a set of bonds a rotationally invariant local order parameterql(r^α)of the particle at positionr^αcan be constructed

from

q_l(r^α)≡



 4π 2l+1

m=l

∑

m=−l

1 N_b

Nb

i=1

∑

Y_lm²(θ_i,φ_i)

2



1/2

, (4.23)

whereN_bis the number of bonds taken into account, the numberlis an even integer¹³, andθ_iand φ_iare the polar and azimuthal angles of thei-th bond, respectively. For the localq_l(r^α)the bonds to the 12 closely packed nearest neighbor particles are taken into account. Alternatively, one can also take all particles within a cutoff radiusr_cinto account to avoid the numerical determination of the nearest neighbors. A global parameter Q_l is defined by averaging over all local order parametersql:

Q_l= 1 N

∑

N

α=1q_l(r^α). (4.24)

Some values of Q_l for perfect crystal lattice structures are listed in table4.1. We implemented the calculation of local and global bond-orientational order parameters Q_l to check for possible defects of the crystal lattice structure.

cluster Q₄ Q₆

fcc 0.15887 0.57452

hcp 0.08089 0.48476

bcc 0.51059

simple cubic 0.35355

random order 0.0 0.0

Table 4.1:Q₄andQ₆for various perfect crystal lattices.

13For evenlthe spherical harmonics are invariant under inversion. Thus, we need not to assign a direction with a particular bond.

Im Dokument Computer-Simulationen zu Strukturen und Phasenumwandlungen in Modell-Kolloiden (Seite 63-69)