Updated: January 21, 2008
Eight Point Force Molecular Dynamical Estimates of
Vibrational Frequencies of an Isotropic Elastic Sphere
Vibrational frequencies of a isotropic elastic sphere are calculated, based on molecular dynamics simulations of a simple cubic lattice with first and second neighbour coupling as well as an eight point force term. Good agreement is found with exactly known frequencies for Poisson ratio ν = 0.15, 0.35 and 0.45.

Introduction
   Molecular dynamics simulations can be used to determine the physical properties of materials [Lee pdf ::] [Iitaka and Ebisuzaki pdf ::]. Simulations of a face centered cubic lattice of atoms interacting through a nearest neighbour potential (linear springs) yields reasonable agreement with exact results for the vibrational modes of a sphere [Murray 2002] even though the lattice does not have isotropic elastic properties. Simulation of a simple cubic lattice interacting through first and second neighbour spring potentials furnishes isotropic elastic constants, although only for Poisson ratio ν = 0.25 [Murray 2002]. Simulations that include only two-body forces are not able to describe the range of material properties [Lee pdf ::]. The present work incorporates an eight point force term that opens up the range of Poisson ratios that can be simulated.

Elasticity Theory
   An isotropic elastic material has two independent elastic constants, which are any two of the Young's modulus (Y), Bulk modulus (B), Shear modulus (μ) and Poisson ratio (ν).
   For an isotropic elastic material (or for a crystalline material with cubic symmetry) the components of the stress and strain tensors are related by the linear equation:

σxx  C11 C12 C12   exx
σyy  C12 C11 C12   eyy
σzz= C12 C12 C11   ezz
σyz  C44   eyz
σxz  C44   exz
σxy  C44   exy

where for an isotropic material
C12 = C11 - 2 C44
so that just two elastic constants, C11 and C44 are sufficient to specific the behavior of an isotropic material.
   These may be used to obtain bulk modulus

B = (C11 + 2 C12) / 3

and shear modulus
μ = (C11-C12)/2

Note that the more familiar relation μ = C44 applies for an isotropic material.
   Plane wave vibrations travelling in an elastic solid can be longitudinally (L) or transversely (T) polarized. Longitudinal waves are associated with compression of the material. Transverse waves create shear strains while keeping the density constant. In an isotropic elastic material with density ρ, the speed of longitudinal waves is CL = sqrt(C11/ρ) and transverse ("shear") waves move at speed CS = sqrt(C44/ρ).

Crystal Lattice Simulation
   The frequencies of the vibrational modes of a crystal sphere are determined by constructing a crystal lattice from point masses connected by springs. The resulting structure is given an initial velocity or initial strain, and the subsequent motion is simulated. Motions of individual atoms are Fourier analyzed to yield the spectrum of vibrational frequencies.
   The mass of the sphere is m. The mass of one atom in the sphere is matom. N atoms are used. The radius of the sphere is R. The atoms are arranged on a simple cubic lattice with side a. The density of the material is ρ = matoma-3. Run time limitations restrict N to 200000 or less, but accurate results are obtainable with just a few thousand.
   Each atom has six nearest neighbours, each a distance a away. A spring of force constant ksp1 connects nearest neighbour atoms. A spring of force constant ksp2 connects each atom to its 12 second neighbours, each at distance sqrt(2) a. In addition, the alternative of a spring of force constant ksp3 connecting each atom to its 8 third neighbours each at distance sqrt(3) a will be considered below.
   There is also an eight point force term involving a force constant k8pt. This force is calculated not for pairs of atoms, but for cubic groups of eight atoms. Apart from atoms on the surface, each atom is involved in eight of these clusters. Consider one such group with atoms numbered j = 0, 1, 2, ... 7. Let p = 0, 1, 2 label the x, y and z axes. The coordinates of the atoms are rjp. The center of mass of the cluster is at RCM where

RCMp = 0.125 Σj rjp

Define the quantity D by

D = Σj,p (RCMp - rjp)2

The p-component of the force on atom j is given by

Fjp = (k8pt/(a2)) (RCMp - rjp) (D - 6 a2)

   Since, for an undistorted lattice, each cluster is a cube of dimensions a × a × a, the distance from each atom to the center of mass is sqrt(3)a/2. D, the sum of the squared distances to the center of mass, is nominally 6 a2. So there is no force when the lattice is undistorted. From its construction, this force respects the translational and rotational symmetry of space just as the two-body forces do. It also satisfies the requirement that the net force exerted by a cluster on itself is zero. Furthermore, since the forces emanate directly outward from the center of mass, the net torque exerted on the cluster by itself is also zero. This 8-point force can equivalently be defined in terms of the potential energy,

U( r0, r1, r2, r3, r4, r5, r6, r7) = (k8pt/(4 a2)) (D - 6 a2)2

The equivalence of the force formula to this potential energy was numerically verified by the short C++ computer program same8ptf.cpp.

Bulk Elastic Constants:
   It is necessary to relate the spring constants to the bulk elastic constants, C11, C12 and C44. This was previously done for ksp1 and ksp2 [Murray 2002]. First, recall that the bulk modulus

B = (C11+2C12)/3

and from considering an infinitessimal isotropic dilatation of the material by δ = dV/V, and the potential energy stored, U = (1/2)V B δ2,

a (C11+2 C12) = ksp1 + 4 ksp2 + 24 k8pt

   Next, considering a longitudinal wave moving in the [100] direction with speed sqrt(C11/ρ),

a C11 = ksp1 + 2 ksp2 + 8 k8pt

A short C++ computer program 8ptc11.cpp was used to numerically reconfirm the coefficients of ksp1 and ksp1 and determine the coefficient of k8pt.
   Finally, considering a transverse wave moving in the [100] direction moving at speed sqrt(C44/ρ),

a C44 = ksp2
where once again 8ptc11.cpp was used.
   Combining these three results, the necessary values of the three force constants are obtained:

ksp1 = a (C11 - (C12 + C44))
ksp2 = a C44
k8pt = (1/8) a (C12 - C44)

Comparison With Exact Analytic Solution:
   The classic problem of the vibrational modes of an isotropic elastic sphere has been extensively studied. [Lamb 1882] [Love 1944] [Sato and Usami 1962] [Auld 1973] [Nishiguchi and Sakuma 1981] [Bastrukov 1994] [Graff 1994] [Bastrukov et al. 1998] [McDaniel and Holt 2000 pdf ::] [Ye 2000 pdf] [Patton and Geller 2001 pdf ::]. This problem finds diverse application in geophysics [geophysics introduction], seismology [simple introduction to earthquake waves], solid stars [Bastrukov 1996], silicon nanoparticles [Patton and Geller 2001 pdf ::] and aqueous foam spheres [McDaniel and Holt 2000 pdf ::], and an accurate list of the exact frequencies has recently been reported [Murray 2002].
   Consider a freely suspended isotropic elastic sphere not subject to body forces. The sphere has radius R, uniform density ρ, shear modulus μ and Poisson ratio ν. These parameters fully describe the problem. If the sphere is disturbed, it will begin to vibrate. A given eigenfunction has frequency ω (in radians per second). The speed of transverse waves (shear waves) is

Cs = sqrt(μ/ρ)

   It is convenient to express normal mode frequencies in terms of the dimensionless variable

η = ω R / Cs

   Modes are of two types: spheroidal and torsional. In spheroidal modes, the velocity field is of the form

u(r,θ,φ,t) exp(-i η Cs t / R ) = Al ∇ [ j(l, kl r ) P(l,cos(θ)) ]
+ Bl ∇× ∇× [ r j(l, ks r ) P(l,cos(θ)) ]

where r is the position vector, j(l,x) are spherical Bessel functions of the first kind and P(l,x) are Legendre polynomials. The ratio Bl/Al must be chosen so as to satisfy the boundary conditions [Murray 2000].
   The expression for the displacement field u for torsional modes is
u(x,y,z,t) exp(-i η Cs t / R ) = Cl ∇× [ (x,y,z) j(l, ks r ) P(l,cos(θ)) ]
    In the calculation, distance units are normally chosen so that the center to center distance between two atoms is 1.
   If the velocities of the atoms became too large, the forces would become nonlinear. Since we are interested in linear elastic behavior, the velocities of the atoms are kept small enough that the forces are linear. Normally, in these simulation, the peak speed of an atom is on the order of 1x10-5 Cs. If the speed is made too small, there are numerical precision problems.
   The sharpness of the Fourier spectrum peaks depends on the amount of time over which the integrations are made.

   Results of these simulations for the case ν = 0.25 have been previously presented [Murray 2002].

Figure 1. Dimensionless frequency parameter η for vibrational modes of an elastic sphere with Poisson ratio ν = 0.15 are shown. The vertical axis is the reciprocal of the sphere's radius. The upper edge of the chart (above the black band) corresponds to infinite R. The bottom edge of the chart corresponds to R=3. The maximum R and maximum number of atoms used are plotted at the lower right. (cc3mod.cpp p15tor2.gif, p15sph2.gif, p15sph1.gif, p15sph0.gif, ... p15tor4.gif )

Figure 2. Dimensionless frequency parameter η for vibrational modes of an elastic sphere with Poisson ratio ν = 0.35 are shown. The vertical axis is the reciprocal of the sphere's radius. The upper edge of the chart (above the black band) corresponds to infinite R. The bottom edge of the chart corresponds to R=3. The maximum R and maximum number of atoms used are plotted at the lower right. (cc3mod.cpp p35tor2.gif, p35sph2.gif, p35sph1.gif, p35sph0.gif, ... p35tor4.gif )

Figure 3. Dimensionless frequency parameter η for vibrational modes of an elastic sphere with Poisson ratio ν = 0.45 are shown. The vertical axis is the reciprocal of the sphere's radius. The upper edge of the chart (above the black band) corresponds to infinite R. The bottom edge of the chart corresponds to R=3. The maximum R and maximum number of atoms used are plotted at the lower right. (cc3mod.cpp p45tor2.gif, p45sph2.gif, p45sph1.gif, p45sph0.gif, ... p45tor4.gif )

Table I: Comparison for ν = 0.15
Typel η
(simul.)		 η
(exact)		 % diff.
sph 0 3.70 3.680 +0.5%
sph 1 3.23 3.195 +1.1%
sph 2 2.65 2.625 +1.0%
sph 2 4.57 4.566 +0.1%
sph 3 3.91 3.867 +1.1%
sph 4 4.98 4.925 1.1%
tor 2 2.65 2.625 +1.0%
tor 3 3.91 3.867 +1.1%
tor 4 5.15 5.095 +1.1%

Table II: Comparison for ν = 0.35
Typel η
(simul.)		 η
(exact)		 % diff.
sph 0 5.83 5.790 +0.7%
sph 1 3.63 3.628 +0.1%
sph 1 7.33 7.288 +0.6%
sph 1 9.05 8.887 +1.8%
sph 2 2.66 2.652 +0.3%
sph 2 5.10 5.140 -0.8%
sph 3 4.05 3.957 +2.4%
sph 4 5.10 5.080 +0.4%
tor 2 2.52 2.502 +0.7%
tor 3 3.83 3.865 -0.9%
tor 4 5.08 5.095 -0.3%

Table III: Comparison for ν = 0.45
Typel η
(simul.)		 η
(exact)		 % diff.
sph 0 10.12 10.006 +1.1%
sph 1 - - - 3.798 - - -
sph 2 2.68 2.662 +0.7%
sph 3 4.02 3.989 +0.8%
tor 2 2.50 2.502 -0.1%
tor 3 3.85 3.865 -0.4%

   The values of η from the simulations as shown in Table I, II and III were obtained by simple extrapolation of Figures 1, 2 and 3). The results in Table I, II and III show that the lattice calculation is able to accurately calculate the vibrational frequencies of an isotropic elastic sphere. This provides an independent check on the analysis as well as on the correctness of the computer program. The rapidity of the convergence as N increases demonstrates that reasonably accurate frequencies could have been obtained even if the simulation been limited to a much smaller number of atoms.

Alternative Approach: Third Neighbour Coupling
   Another way of setting up the calculation is to omit the second neighbour couplings but instead introduce a force constant ksp3 between each atom and its 8 third neighbours. Making use of results in md11.htm,
ksp1 = a (C11 - C12)
ksp3 = (3a
---
 4		)   C44
k8pt = ( a
---
 8		) (C12 - C44)

   The advantage of this approach is that ksp1 remains positive as long as C12 < C11. This is true for all materials that I am aware of. (If ksp1 is negative the lattice becomes unstable, making simulation impossible.) This allows materials with Zener anisotropy (A = 2C44/(C11-C12)) above 2 (e.g. InAs, InP, GaN (ZB) and crystalline Ag, Fe, Ni, Cu and Au) to be handled.

References:

H. Lamb, "On the vibrations of an elastic sphere" Proc. London Math. Soc. 13, 189 (1881-1882).
A. E. H. Love, A Treatise on the Mathematical Theory of Elasticity, (Dover, New York, 1944).
Y. Sato and T. Usami, Geophys. Mag. 31, 15 (1962)
C. Kittel, "Introduction to Solid State Physics" Wiley 4th Ed. 1971, see table 4 on page 38
B. A. Auld, Acoustic Fields and Waves in Solids, (John Wiley & Sons, New York, 1973).
N. Nishiguchi and T. Sakuma, Sol. Stat. Comm. 38, 1073 (1981).
Sergey I Bastrukov, Phys. Rev. E 49, 3166 (1994).
F. K. Graff, Wave Motion in Elastic Solid, (Ohio University Press, Ohio, 1994).
S. Bastrukov, Phys. Rev. E (1996)
Bastrukov et al., Physica A 250, 435 (1998).
Zhen Ye, "On the Low Frequency Elastic Response of a Spherical Particle" Chinese Journal of Physics Vol. 38, pages 103-110, (2000). pdf
J. Gregory McDaniel and R. Glynn Holt "Measurement of Aqueous foam Rheology by acoustic levitation" Phys. Rev. E Rapid Communications volume 61 R2204 (2000). pdf ::
Kelly Patton & Michael Geller "Phonon Spectrum in a Nanoparticle mechanically coupled to a substrate" Journal of Luminescence 94-95 (2001) 747-570 pdf ::
Michael R. Geller et al. "Theory of electron-phonon dynamics in insulating nanoparticles" Preprint 11 July 2001 t / d where vt is the bulk transverse sound velocity." > pdf ::
Kelly R. Patton, Michael R. Geller, "Phonon Spectrum in a Nanoparticle Mechanically Coupled to a Substrate" (Preprint Tues, June 26, 2001) abstract, preprint pdf ::
G. Cerullo, S. De Silvestri and U. Banin "Size-dependent dynamics of coherent acoustic phonons in nanocrystal quantum dots" Phys. Rev. B volume 60 (July 15, 1999) pdf ::
H. Portales, L. Saviot, E. Duval, M. Gaudry, E. Cottancin, M. Pellarin, J. Lermé and M. Broyer "Resonant Raman Scattering by Quadrupolar Vibrations of Ni-Ag Core-shell nanoparticles" (Preprint Mar 22 2002) pdf ::
Daniel T. Simon and M. Geller "Electron Phonon Dynamics in an Ensemble of Nearly Isolated Nanoparticles" (Preprint June 8, 2002) To appear in Phys. Rev. B volume 64 PRB preprint :: preprint pdf ::
Daniel B. Murray, "Molecular Dynamics Simulation of an Elastic Material" 2002 (link to article)
Daniel B. Murray, "Vibrational Frequencies of an Elastic Sphere" 2002 (link to article)

Byeong-Joo Lee "A Semi-Empirical Atomistic Approach In Materials Science and Engineering" pdf ::

also
http://www.postech.ac.kr/~calphad/files/kissme.pdf

http://www.ees4.lanl.gov/nonlinear/iwnem/abstracts.html TJ Ulrich, KR McCall, PA Johnson, TW Darling, A Migliori "Application of Resonant Ultrasound Spectroscopy (RUS) to Determine the Elastic Properties of Rock Samples " Summer Workshop 1999
T. Iitaka and T. Ebisuzaki "First principles calculation of elastic properties of solid argon at high pressures" pdf ::

Albert Brown "..Ultrasonics.." ::
Daniel Murray
Associate Professor
Math, Stats & Physics Unit
University of British Columbia - Okanagan
Kelowna, BC, Canada
daniel "dot" murray "at" ubc "dot" ca

For a list of related articles click on the link.
Hosted by www.Geocities.ws

1