Updated: Feb 15, 2003
Molecular Dynamic Simulation of an
Elastic Solid with Hexagonal Symmetry

A molecular dynamics simulation is developed for an elastic solid with general hexagonal symmetry. This allows the vibrational frequencies of solid objects made of a single crystal with hexagonal symmetry to be determined. The method is tested for the isotropic special case for exactly known frequencies for an isotropic sphere.

   A C++ computer program hc2mod.cpp is available to accurately compute the vibrational frequencies of an elastic sphere to within the order of 1%. A separate computer program hex2.cpp is used to determined the five force constants starting from the five elastic constants of the material. The lattice used is hexagonal with lattice constants a = 1 and c = 0.93. The ratio c/a must be chosen so as to avoid having some negative force constants. (For the case when Poisson ratio = 0.45, c was changed to 0.80 to avoid instability.) However, this c/a need not equal the c/a ratio of the actual crystal.
   Five different force constants are required in order to represent the general elasticity of a material with hexagonal symmetry. Such a material has five independent elastic constants, C11, C12, C13, C33 and C44. Other elastic constants are dependent on these through: C11=C22=C33, C12=C21, C13=C31=C23=C32, C44=C55, C66=(C11-C12)/2. Finally, constants Cij and Cji with i=1,2,3 and j=4,5,6 are all 0.
   In selecting force rules for the lattice, it is desirable that the connections be as short as possible in order to quickly achieve convergence to the continuum limit as the object increases in number of particles. It is not possible to represent general elastic properties in terms of two-particle force laws alone.
   The first three force constants are simple two-particle force laws: A spring with force constant ksp1 connects particles at distance a (assuming a is unequal to c). Each particle has six such neighbours. A spring with force constant ksp3 connects particles at distance c (assuming a is unequal to c). Each particle has two such neighbours. A spring with force constant ksp4 connects particles at distance sqrt(a2+c2). Each particle has twelve such neighbours.
   There is a three point force term involving a force constant k3pt. This force is calculated not for pairs of atoms, but for equilateral triangular groups of three atoms. Apart from atoms on the surface, each atom is involved in six of these clusters. Consider one such group with atoms numbered j = 0, 1, 2. 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 = (1/3) Σ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 = k3pt (RCMp - rjp) (D - a2)

   Since, for an undistorted lattice, each cluster is a triangle of side a, the distance from each atom to the center of mass is a/sqrt(3). D, the sum of the squared distances to the center of mass, is nominally 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 3-point force can equivalently be defined in terms of the potential energy,

U( r0, r1, r2) = (??coef??) k3pt (D - a2)2

   There is a six point force term involving a force constant k6pt. This force is calculated for groups of six atoms. Apart from atoms on the surface, each atom is involved in twelve of these clusters. Consider one such group with atoms numbered j = 0, 1, 2, ... 5. 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 = (1/6) Σ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 = k6pt (RCMp - rjp) (D - 6(a2/3+c2/4))

   For an undistorted lattice, the distance from each atom to the center of mass of the cluster is (a2/3+c2/4). D, the sum of the squared distances to the center of mass, is nominally 6 (a2/3+c2/4). 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 six-point force can equivalently be defined in terms of the potential energy,

U( r0, r1, r2, r3, r4, r5) = (??coef??) k6pt (D - 6 (a2/3+c2/4))2


Table I. Force constants for objects with
isotropic elasticity
(Using C++ program hex2.cpp.)
constant units ν=0.15 ν=0.25 ν=0.35 ν=0.45
a (metres) 1 1 1 1
c (metres) 0.93 0.93 0.93 0.8
C11 (GPa) 90 90 90 90
C44 (GPa) 37.06 30 - - 8.182
k6pt (???) -3.287e9 0 4.298e9 1.181e10
ksp1 (N/m) 5.659e10 4.581e10 3.171e10 9.211e9
k3pt (???) -3.909e9 0 5.112e9 -1.890e9
ksp3 (N/m) 4.970e10 3.548e10 1.690e10 1.807e10
ksp4 (N/m) 2.145e10 1.737e10 1.202e10 4.842e9

Figure 1. Fourier transforms of lattice simulations of isotropic elastic spheres with ν=0.35 are shown. The initial disturbance of the lattice is as shown at lower left. For each plot, the horizontal axis is η = ω R/CS. The vertical axis is 1/R, with infinite R at the top. R at the bottom is 3. The maximum radius and number of atoms, N, used is shown at the lower right. (hc2mod.cpp hcp35sp0.gif hcp35to2.gif...)

Figure 2. Fourier transforms of lattice simulations of isotropic elastic spheres with ν=0.45 are shown. The initial disturbance of the lattice is as shown at lower left. For each plot, the horizontal axis is η = ω R/CS. The vertical axis is 1/R, with infinite R at the top. R at the bottom is 3. The maximum radius and number of atoms, N, used is shown at the lower right. (hc2mod.cpp hcp45to2.gif...hcp45sp0.gif )

Figure 3. Fourier transforms of lattice simulations of isotropic elastic spheres with ν=0.15 are shown. The initial disturbance of the lattice is as shown at lower left. For each plot, the horizontal axis is η = ω R/CS. The vertical axis is 1/R, with infinite R at the top. R at the bottom is 3. The maximum radius and number of atoms, N, used is shown at the lower right. (hc2mod.cpp hcp15to2.gif...hcp15sp0.gif )

Table II: Comparison for ν = 0.35
Typel η
(simul.)
η
(exact)
% diff.
sph 0 5.80 5.790 +0.2%
tor 2 2.50 2.502 -0.1%
sph 2 2.65 2.652 -0.1%
sph 1 3.65 3.628 +0.6%
tor 3 3.95 3.865 +2.2%
sph 3 3.96 3.957 +0.1%

Table III: Comparison for ν = 0.45
Typel η
(simul.)
η
(exact)
% diff.
tor 2 2.52 2.502 +0.7%
sph 2 2.66 2.662 -0.1%
sph 1 no peak 3.798 N/A
tor 3 3.95 3.865 +2.2%
sph 3 3.98 3.989 -0.2%
sph 0 10.05 10.006 +0.4%

Table IV: Comparison for ν = 0.15
Typel η
(simul.)
η
(exact)
% diff.
tor 2 2.50 2.502 -0.1%
sph 2 2.65 2.625 +1.0%
sph 1 3.18 3.195 -0.5%
sph 0 3.68 3.680 0%
tor 3 3.93 3.865 +1.7%
sph 3 3.93 3.867 +1.6%

   The agreement between the molecular dynamics simulations and exactly known frequencies for spheres with isotropic elasticity is within 2%, and usually within 1%. Only in the l = 1 spheroidal case at ν = 0.45 did a peak fail to appear. The rate of convergence is variable, the general pattern being that convergence is best for the lowest frequency modes. The estimated frequencies are more often higher than the exact frequencies.

References:

D. B. Murray "Eight Point Force Molecular Dynamical Estimates of Vibrational Frequencies of an Isotropic Elastic Sphere" 2002 link to article

D. B. Murray "Breathing Mode Vibrations of an Isotropic Elastic Sphere Surrounded by a Fluid Medium" 2002 link to article

D. B. Murray "Vibrational Frequency of Silicon Nanoparticles" 2002 link to article

http://www.chem.vt.edu/chem-dept/tissue/nano/
"While the elastic modes of spherical nanoparticles can be calculated analytically, this is not the case for nonspherical nanoparticles." "Visscher et al. ...mode structure for ellipsoidal nanoparticles..." ::

W. M. Visscher, A. Migliori, T. M. Bell and R. A. Reinhart, "The normal modes of free vibration of inhomogenous and anisotropic elastic objects" J. Acoust. Soc. Am. (Journal of the Acoustic Society of America) Vol. 90 (1991) 2154-2161.

A. Stekel, J. L. Sarrao, T. M. Bell, M. Lei, R. G. Leisure, W. M. Visscher, A. Migliori, "Method of Identification of the Vibrational Modes of a Rectangular Parallelipiped" J. Acoust. Soc. Am. vol. 92 no. 2 pt. 1 (August 1992) pp. 663-668.

A. Migliori, J. L. Sarrao, W. M. Visscher, T. M. Bell, M. Lei, Z. Fisk and R. G. Leisure, "Resonant ultrasound spectroscopic techniques for measurement of the elastic moduli of solids" Physica B 183 1-24 (1993)

M. K�gl, L. Gaul A 3-D Boundary Element Method for Dynamic Analysis of Anisotropic Elastic Solids Computer Modeling in Engineering & Sciences, 1(4) (2000), pp. 27-43 abstract


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 here.

CdS properties:
http://www.efunda.com/materials/piezo/material_data/matdata_output.cfm?Material_ID=CdS
Piezo Data: CdS
Hexagonal
??? density 5684 kg/m^3. ???? seems wrong
Compliance SE 10-12 m2/N: (inverse of 6x6 stiffness C matrix ! )
20.69 -9.99 -5.81 0 0 0
-9.99 20.69 -5.81 0 0 0
-5.81 -5.81 16.97 0 0 0
0 0 0 66.49 0 0
0 0 0 0 66.49 0
0 0 0 0 0 61.36
http://www.issp.ac.ru/lpcbc/DANDP/cds_adv.html
Density 4825 g/cc (agrees) ???!!!!
Young's modulus 45 GPa
http://www.vidrine.com/iropmat3.htm
Gives density as 4.82 g/cc. (agrees)
??? CdSe density 5.67 g/cc. !!!???
http://ncsr.csci-va.com/materials/cdse.asp
CdSe (Wurtzite) 5.81 g/cc

Article on CdSe modelled using pair potential
Table III gives calc and experi. elastic constants for wurtzite, zinc blende and rocksalt forms
CdSe Experimental values (Wurtzite = hexagonal):
C11 = 74.6 GPa
C12 = 46.1 GPa
C13 = 39.4 GPa
C33 = 81.7 GPa
C44 = 13.0 GPa
C66 = 14.3 GPa
B = 53.4 GPa
from: Eran Rabani "An interatomic pair potential for cadmium selenide" Journal of Chemical Physics vol. 116, no. 1 (2002) pdf ::

Michael Conry "Notes on Wave propagation in anisotropic elastic solids"
CdS hexagonal
C11=90.7 GPa
C12=58.1 GPa
C33=93.8 GPa
C13=51.0 GPa
C44=15.04 GPa
Source: pdf ::

Has general axis formula for Young's modulus!
1/E = (C11+C12) +...
http://www.crystran.co.uk/optics.htm
Sapphire Al2O3
C11=496
C12=164
C13=115
C33=498
C44=148
B = 240 GPa
http://www.tydex.ru/materials/materials2/sapphire.html

W. Q. Chen, H. J. Ding, R. Q. Xu "Three-dimensional static analysis of multi-layered piezoelectric hollow spheres via the state space method" International Journal of Solids and Structures 38 (2001) 4921-4936 ::

Chen, W. Q. and Ding, H. J. (2001) Free vibrations of multi-layered spherically isotropic hollow spheres. International Journal of Mechanical Science, Vol. 43(3), 667-680. ::

H. Cohen, A. H. Shah and C. V. Ramakrishnan 1972 Acustica 26 329-333, "Free Vibrations of a Spherically Isotropic Hollow Sphere"

P. R. Heyliger and A. Jilani 1992 International Journal of Solids and Structures 29 2689-2708 "The Free Vibrations of Inhomogeneous Elastic Cylinders and Spheres"

Chen, W. Q., Cai, J. B., Ye, G. R. and Ding, H. J. (2000) "On eigenfrequencies of an anisotropic sphere". ASME Journal of Applied Mechanics, Vol. 67(2), 422-424.

Chen, W. Q. "Effect of radial inhomogeneity on natural frequencies of an anisotropic hollow sphere" Journal of Sound and Vibration, Vol. 226(4), 787-794.(1999) ::

Chen, W. Q. and Ding, H. J. (1997) "On free vibrations of an embedded anisotropic spherical shell" ASME Journal of Pressure Vessel Technology, Vol. 119(4), 481-487.

Ding, H. J. and Chen, W. Q. (1996) "Natural frequencies of an elastic spherically isotropic hollow sphere submerged in a compressible fluid medium". Journal of Sound and Vibration, Vol. 192(1), 173-198.

W. Q. Chen and L. Z. Wang, "Free Vibrations of Functionally Graded Piezoceramic Hollow Spheres with Radial Polarization" Journal of Sound and Vibration vol. 251, pages 103-114 (2002) ::

W. Q. Chen, "Vibration Theory of Non-Homogeneous Spherically Isotropic Piezoelastic Bodies" Journal of Sound and Vibration" vol 236 pages 833-860 (2000) ::

Palko JW, Kriven WM, Sinogeikin SV, Bass JD, Sayir A "Elastic constants of yttria (Y2O3) monocrystals to high temperatures"
JOURNAL OF APPLIED PHYSICS 89 (12): 7791-7796 JUN 15 2001
http://hercules.geology.uiuc.edu/~stas/Abstracts/JAP01_89_7791.htm
Room temp: C11=223.7(0.6) C12=112.4(1.1) C44=74.6(0.7) GPa K = 149.5(1.0) GS=66.3(0.8) [Voigt-Reuss-Hill average]

Formula for Voigt and Reuss for cubic crystal (in french; has error in C' definition)
http://www.ens-lyon.fr/LST/WEBHP/elasticite/node2.html
GV = (2C'+3C44)/5 = Voigt Bound
GR = 5C'C44/(2C44+3C') = Reuss Bound
GS=(GV+GR)/2 = Voigt-Reuss-Hill average
C'=(C11-C12)/2
Also, for a cubic crystal, the VRH bulk modulus is (C11+2C12)/3

"Gas-phase-condensed Y2O3:Eu3+ nanoparticles of < 10 nm contain multiple phases: monoclinic Y2O3:Eu3+, cubic Y2O3:Eu3+, monoclinic Eu2O3, and a disordered phase. Larger particle sizes predominantly form in the metastable monoclinic phase" "Annealing as-prepared 4-nm Y2O3:Eu3+ at 500-900 C produces cubic-phase Y2O3:Eu3+ and removes all other phases. " link

List of recent publications by Takagahara group (in Japanese)
http://www.hiei.kit.ac.jp/~takaghra/

Noninertial mechanism for electronic energy relaxation in nanocrystals Ho-Soon Yang,* Michael R. Geller, and W. M. Dennis Department of Physics and Astronomy, University of Georgia, Athens, Georgia 30602-2451 pdf ::

H. S. Yang, M. R. Geller Dennis, "New Mechanism for Electronic Energy Relaxation in nanocrystals" Aug 15 2000 preprint

H. S. Yang, M. R. Geller, W. M. Dennis, "Noninertial Mechanism for electronic energy relaxation in nanocrystals" Phys. Rev. B vol. 62 (2000) 9398. ::

T. Takagahara, J. Lumin. 70, 129-143 (1996)
"Electron-phonon interactions in semiconductor nanocrystals"

U. Banin, G. Cerullo, et al. "Quantum confinement and ultrafast dephasing dynamics in InP " Phys Rev. B vol 55, March 1997 ::

A. M. Alcalde, G. E. Marques, G. Weber. T. L. Reinecke "Electron acoustic phonon scattering rates in II-IV quantum dots: contribution of the macroscopic deformation potential" Solid State Communications 116 (2000) 247-252 pdf ::

Hosted by www.Geocities.ws

1