Molecular Dynamics Simulation of an Elastic Material

First created: June 2002 Last updated: July 13, 2002
Molecular Dynamics Simulation of an Elastic Material

A crystalline face centered cubic solid consisting of point masses interacting through the Lennard Jones potential and one other model potential is characterized as an approximation for an isotropic linear elastic solid, and the linear elastic properties are found. Failure characteristics are also obtained. This model is then used to study vibrational modes of objects, as well as the motion of bouncing objects.

Molecular dynamics simulation of polycrystalline materials [A. M. Krivstov and M. Wiercigroch 2001].

The figure on the right shows a ball of radius 13.36 consisting of 10023 atoms. Each atom has a radius of 0.561231. The atoms are arranged in the "face centered cubic" structure, and are thus packed as tightly as they could be.

   To study the behavior of a solid, liquid or gas, a computer can be used to calculate the motions of all the individual molecules as they evolve in time. This approach is called "molecular dynamics" simulation.
   The motions of the individual atoms depend on the forces that they exert on each other. This depends on the potential energy, which depends on the distance between two atoms.
   A very commonly used approximate model potential energy formula is the Lennard Jones or "6-12" potential which has the following formula:

U(r) = 4ε( (σ/r)¹² - (σ/r)⁶ )
where ε and σ are constants; ε is in joules and σ is in metres. U is the potential energy in Joules and r is the center to center distance between two atoms.
   It should be mentioned that the "6-12" model is not a truly realistic model of interatomic forces. The Lennard Jones potential energy could only be a plausible model for a solid made of atoms of a noble gas, like Helium, Neon or Argon, which bond primarily through Van der Waals forces rather than ionic or covalent bonds. Normally, a solid cannot be exactly modelled in terms of two-body potentials of atoms. The details of the quantum mechanical behavior of the electrons in the solid need to be taken into account.
   The point of introducing the molecular dynamics model is that an elastic material can be simulated in a convenient and intuitively straightforward way.
   In the calculations introduced here, two different model potential energies are used. The first is the Lennard Jones potential. I set σ = 1 and choose ε = 0.0174989 so that dF/dr=1 at the equilibrium separation. The equilibrium separation distance, d_cc, is 1.1224621 times σ. The second model used is a quadratic potential

U(r) = ((r-d_cc)²)/2
   This model thus has the same equilibrium distance and the same dF/dr as the Lennard Jones model. In addition, only interactions between nearest neighbour atoms are considered. Consequently, for very small distortions of the lattice, the two models are the same.
   A solid object (such as a ball or dice) is modelled as an array of these atoms. Initially the atoms are arranged in a close packed structure (face-centered cubic). The density of atoms is exactly 1.000 atoms per cubic metre. (This is a bit of a surprise, since the distance between the centers of two atoms is 1.1224621) I assume that each atom has a mass of 1.00 kg.
   The two models have the following mechanical properties in common:

Crystal structure: Point particles could be assembled into a crystal in a number of ways, but I always use the face centered cubic structure. To obtain this structure, start with a cube of dimensions 1.58740 x 1.58740 x 1.58740. The volume of the cube is exactly 4. Place an atom at each of the eight corners of the cube. Only one eighth of each atom is inside the cube. Then place an atom at each of the six centers of the faces of the cube. One half of each atom is inside the cube. The distance across a face diagonal is equal to four times the atom radius. This is 1.41421 x 1.58740 = 2.244992. Thus, the radius of each atom is 0.561231. This makes the center to center distance 1.122461. Since a total of four atoms are inside the cube, and the volume of the cube is 4, the density of the atoms is exactly 1.
Density is 1.0 atom per cubic metre. It is also one kilogram per cubic metre.
Young's Modulus (also known as "elastic modulus" and "tensile modulus", symbol Y (other literature often uses E)) varies with the orientation of the lattice with respect to the axis along which strain is being applied. On average, Y is 1.215. It varies from 1.172 to 1.261. These figures were obtained from molecular dynamics computer simulations using the program md4ym.cpp which puts a cubic sample of material under gradually increasing uniform plane strain and measures the resulting stress. Young's modulus is stress/strain.
Bulk Modulus (symbol B (other literature often uses K)) is 0.8399473 (exactly (1/3)(2^(4/3)) ) I numerically calculated bulk modulus as a double check using md2bm.cpp, and found it was close to this value.
   The exact derivation of B goes as follows. Start with a sphere of radius R₀ and volume V₀ at equilibrium. The initial pressure is zero. Then compress the sphere from volume V₀ to V₁. During the compression, P(V) = A (V₀-V), where A is some constant. Work, W, is obtained by integrating P dV from V₀ to V₁. This is W = integral from V₀ to V₁ of A(V₀-V) dV, which is 0.5 A dV², where dV=V₀-V₁. The definition of bulk modulus B is B = dP/(dV/V) = V(dP/dV) = A V. But A = 2 W (dV)^-2 and so B = 2 W V (dV)^-2. All of this energy W goes into the potential energy of atom to atom bonds. Originally all bonds are of length d_cc. After compression all bonds are identically compressed to length d_cc - dx. The force is F = (dF/dr) dx. The total work done during the compression is 0.5 (dF/dr) dx². There are N atoms in the sphere and each atom has 12 neighbours. But each bond must be counted once, not twice. So the total energy stored is U = (12 N/2) 0.5 (dF/dr) dx². Inserting the expression for U into the expression for B, we get B = 2 (12 N/2) 0.5 (dF/dr) (V) (dx/dV)². Furthermore, if dx is small then 3 dx/d_cc = dV/V, and so dx/dV = d_cc/(3 V). This yields, B = 2 (12 N/2) 0.5 (dF/dr) (V) (d_cc)² (1/9) (1/V²). Also note that N = V ρ / m, so that

B =

2 ρ
----
3m

(dF/dr) (d_cc)²

This gives B = 0.8399472.
Poisson ratio (ν) can be obtained from ν=(3B-Y)/(6B). Since Y varies with orientation, Poisson's ratio also varies. Its limiting values are 0.250 to 0.267. The value of ν based on the average value of Y (1.215) is ν=0.2589.
Shear modulus (symbol G (others often use μ)) can be obtained from other elastic constants through the relation G = 0.5*Y/(1+ν). However, since both Y and ν vary together, it is more appropriate to use an equation depending on Y alone:
G = 3 B Y / (9B-Y)
This gives G ranging from 0.462 to 0.504 with the average of Y giving G = 0.4826.
The program md2sm.cpp was used to simulate an object placed under shear strain with upper and lower surfaces held fixed. It is important that the arrangement be such that the strain field is approximately constant, which requires that the object have one short dimension. Simulations of a objects of dimensions up to 39 by 39 by 6 (9100 atoms) gave values for shear modulus in reasonable agreement.
1st Lamé constant (symbol λ) Using B = λ + (2/3)G, λ = B - (2/3)G, so λ = 0.5182
Transverse speed of sound (symbol c_s) The speed of transverse waves is sqrt(G/ρ) = 0.6947 m/s.
Longitudinal speed of sound (symbol c_l) The speed of longitudinal (i.e. compressional) waves is sqrt((λ+2G)/ρ) = 1.2179 m/s.
These "constants" vary, since this solid is not an isotropic elastic material. It is a cubic crystal. Therefore, the relationship between stress and strain cannot be reduced to two elastic constants. However, it seems to be a useful approximation to view this solid as being isotropic, and to make use of it as a model of isotropic materials.
For a cubic crystal, linear elasticity theory gives three elastic constants, C₁₁,C₁₂ and C₄₄ [N. W. Ashcroft and N. D. Mermin, "Solid State Physics" Holt, Rinehart & Winston 1976, Table 22.1, page 445]. If we were dealing with a simple cubic lattice, then starting at any lattice point, there are six nearest neighbours, and therefore six directions in which the speed of sound must be equivalent. But dealing with a face centered cubic lattice there is even more symmetry, since from any lattice point there are twelve nearest neighbours, and consequently twelve directions along which the speed of sound must be the same.
(A concise and very convenient review of mechanics of materials can be found at this website.)

Linear Elastic Properties of Molecular Dynamics Solid:
Property	Symbol	Average value	Unit	Min. value	Max. value	Relation:
Mass of one atom	m_atom	1.0000000	kg	(exact)
Density	ρ	1.0000000	kg/m³	(exact)
Atomic force constant	dF/dr	1.0000000	N/m	(exact)
Atom center-center distance	d_cc	1.1224621	m	(exact)		d_cc=(sqrt(2)m_atom/ρ)^1/3
Bulk Modulus	B	0.8399473	Pa	(exact)		2ρ(d_cc)²(dF/dr)/(3m_atom)
Young's Modulus	Y	1.215	Pa	1.172	1.261	(use md4ym.cpp)
Poisson's ratio	ν	0.2589	- - -	0.250	0.267	ν=(3B-Y)/(6B)
Shear Modulus	G	0.4826	Pa	0.4623	0.5045	G=3BY/(9B-Y)
1st Lamé constant	λ	0.5182	Pa	0.5036	0.5317	λ = B - (2/3)G
Shear wave speed	C_s	0.6947	m/s	0.6799	0.7103	C_s=sqrt(G/ρ)
Longitudinal wave speed	C_l	1.2179	m/s	1.2068	1.2299	C_l=sqrt((λ+2G)/ρ)
C_s/C_l		0.5704	- - -	0.5634	0.5775

Material Failure:

Failure under plane stress is depicted (left) for a cube of 1729 atoms. The top and bottom surfaces (darker atoms) are slowly pulled apart. Strain increases by 0.0001 per iteration, and the time step per iteration is 0.2 seconds. Every 200th iteration is used for this animation. The Lennard Jones potential was used.

Failure under shear stress is depicted (right) for a cube of 1727 atoms. The top and bottom surfaces (darker atoms) are forced sideways. Strain increases by 0.0001 per iteration, and the time step per iteration is 0.2 seconds. Every 200th iteration is used for the animation. The Lennard Jones potential was used.

The same numerical computer experiment that allows Young's Modulus to be determined can also be used to obtain the failure characteristics of the material. The Lennard Jones model can be made to fail by applying excessive pulling strain (as opposed to compressive strain). It fails after it has elongated by about 17%. The quadratic potential model does not fail in the same way when it is pulled. However for the purposes of studying bouncing it is more relevant to look at compressive strain. The quadratic model undergoes failure when the strain passes 20%. The Lennard Jones model fails at a lower strain.
The computer program md2sm.cpp was used to simulate failure under excessive shear strain. For the Lennard Jones model, failure was found for shear strains varying from 0.125 to 0.16. The failure strain varied from 0.028 to 0.038. The quadratic model does not exhibit any noticeable failure under large shears, but simply stretches more and more.

Sphere vibrational frequencies:

The breathing mode (shown at right) for a sphere of radius 6.5, or 1159 atoms. It has a frequency of 0.0751 Hz. The simulation time step is 0.50 s. The animation time step is 1.0 s. (fig1j.gif)

The axial stretch mode (shown at left) for a sphere of radius 6.5 atoms. It has a frequency of 0.041 Hz. The simulation time step is 0.50 s. The animation time step is 1.0 s. (fig1m.gif)

The wet dog mode (shown at right) for a radius of 6.5 (1159 atoms). It has a frequency of 0.0364 Hz. The simulation time step is 0.50 s. The animation time step is 1.0 s. This is the lowest frequency mode of vibration. It is reminiscent of the manner in which a dog shakes its body to dry itself off. (fig1i.gif)

   The computer program md2mod.cpp was used to find vibrational frequencies of an elastic sphere with Young's Modulus 1.215, Poisson ratio 0.2589 and radius R. There are different modes of vibration. For the mode with isotropic dilation, (i.e. the "breathing mode") for a sphere of radius 14, the frequency, f is 0.0353 Hz. For spheres of radius from 11 to 14.7 (13725 atoms) it was found that f = 0.495/r. Surprisingly, this is not the lowest frequency vibrational mode. To find lower frequency modes, a sphere of radius 13 (9195 atoms) was randomly excited. The first four frequencies were found by Fourier analyzing the motions of two randomly chosen atoms. Two atoms were used in case a node happened to fall at or near one of the atoms. The first four were 0.0191 Hz, 0.0198 Hz, 0.0212 Hz and 0.0222 Hz. To clarify what kind of motion these frequencies correspond to, a sphere was "kicked" straight up from its bottom. This gave the following frequencies: 0.0221 Hz, 0.0251 Hz, 0.0279 Hz and 0.0327 Hz. After that, a sphere was kicked up and to the side from its bottom. This gave 0.0198 Hz, 0.0208 Hz 0.0212 Hz and 0.0221 Hz. This still did not produce the lowest frequency mode. To produce the 0.0191 Hz mode, a sphere was given an axial twist, or "wet dog" motion. This did produce the 0.0191 Hz peak. This was repeated for a sphere of radius 17.7 (13275 atoms) and f was 0.0169 Hz.
To systematically find all of the low frequency vibrational modes for an elastic sphere, I used md2mod3.cpp in which a sphere of up to 12000 atoms was given random initial velocities, and the motion of two atoms was then Fourier analyzed. This is done several times in case the atoms fall at a node. The quadratic potential was used.

Cylinder vibrational frequencies:
   Looking at a cylinder of radius R and height = 2*R*flatness, I found that the lowest frequency mode has a locus of nodes which forms a "+" sign. I call this the "saddle" mode. The C++ program is md2mod2.cpp. For example, I simulated a cylinder of radius 19.1 and flatness 0.3 (13107 atoms) and found a lowest frequency of 0.0667 Hz. I checked for convergence as the number of atoms was increased by varying r while keeping "flatness" constant. The frequency is inversely proportional to the radius in this case. The radius could be decreased to 5 before significant variation was noticed.
   It is of interest to see how the lowest frequency mode varies as flatness decreases. Doing this, I obtain the general formula for the lowest frequency:

   f lowest = (0.30/R) (flatness^0.77)

For example, consider a coin made of steel so that sqrt(Y/ρ) is approximately 5200 m/s. Assume r is 0.008 m and flatness is 0.07. In that case the lowest frequency of vibration will be 25??? kHz. This vibrational mode should explain the high pitched ringing of a coin bouncing on a hard surface. The lowest frequencies should be most relevant to energy loss of a tossed coin on a hard surface. The frequency of the "cup" mode is about 60% higher. There are other modes in between, and certainly many more modes of higher frequency. In some cases I looked for these modes by applying random velocities in random directions to randomly chosen atoms. I then fourier transformed the velocity as a function of time of two randomly chosen atoms. This usually reveals the lowest frequency mode, but in some cases repeated attempts are necessary. It was found that the lowest frequency mode always corresponded to the "saddle" mode. It was then more convenient to begin by creating a "saddle" distortion of the object, increasing the rate at which a clear frequency peak appears.
Another example is a polystyrene poker chip. Sqrt(Y/ρ) for polystyrene is on the order of 1500 m/s. Assume r = 0.015 m and flatness = 0.05. The frequency would be roughly 3000??? Hz. The audible sound when poker chips are rubbed together should be a mixture dominated by the lower frequency modes.
Yet another example is a ceramic (Corning Corelle) dinner plate. I can only guess that sqrt(Y/ρ) is 2000 m/s. The radius r is 0.10 m. The flatness is 0.02 or so. The frequency should be about 300??? Hz. It is interesting to take a real dinner plate and strike it in different places while holding it in different ways and hearing the different audible frequencies.

Cube vibrational frequencies:
   The lowest vibrational frequency of a cube corresponds to the "wet dog" mode. The frequency for this mode is:

   f lowest = 0.286 / d

where d is the side length of the cube. The dimensions of the cube are d by d by d. The next mode has a frequency 1.49 times higher. The "breathing" mode has a frequency 2.44 times higher. In this case the largest cube simulated had side length 23.6 which means its dimensions were 23.6 by 23.6 by 23.6, and it involved 13154 atoms, using md2mod2.cpp. Thus, the "wet dog" mode is the dominant distortion resulting from a long duration collision, such as a bounce on the floor. I sometimes used 6 point random initial velocities to ensure that the lowest mode was not missed.

Bouncing Spheres:

The figure on the right was generated by the C++ program md2graf3.cpp. It shows a ball made of 727 atoms and radius 5.5 falling with initial velocity 0.05 m/s. The ball rebounds with a velocity of 0.0476 m/s. Some of the energy can be seen to go into rotational kinetic energy. The simulation time step is 0.50 s, but each animation time step is 5 times longer, or 2.50 s. The quadratic potential was used. The floor is infinitely hard.

The figure on the left was generated by the same C++ program, and is also a ball of 727 atoms. This time it falls at 0.10 m/s. The time step is 0.10 seconds, and the animation time step is 2.50 s. The ball rebounds upwards with a speed of 0.0359 m/s. This simulation is 52 frames long. The Lennard Jones potential was used. (fig1f.gif)

The picture on the right shows a ball dropped at 0.20 m/s, with a simulation time step of 0.05 s and animation time step of 2.50 s. The animation is 75 frames long. The Lennard Jones potential is used. (fig1g.gif)

On left is an image of a 727 atom ball dropped at 0.40 m/s, also with an animation time step of 2.50 s. The animation is 51 frames long. The Lennard Jones potential is used. The simulation time step was 0.02 s. (fig1h.gif)

References:

A. M. Krivstov and M. Wiercigroch, "Molecular Dynamics Simulation of Mechanical Properties for Polycrystal Materials" Mater. Phys. Mech. 3(2001) 45-51 link to pdf ::

M.P.Allen and A.K. Tildesley, "Computer Simulation of Liquids" Clarendon Press, Oxford, 1987

J. M. Haile, "Molecular Dynamics Simulation: Elementary Methods" Wiley 1992

M. P. Allen, A. Warren and M. R. Wilson, Phys. Rev. E 57, p.5585 (1998).

B. L. Holian and P. S. Lomdahl, Science 280, p. 2085 (1998)

D. R. Collins, W. Smith, N. M. Harrison and T. R. Forrester, J. Mat. Chem. 7 p.2543 (1997)

J. Schiotz, T. Vegge, F. D. Di Tolla and K. W. Jacobsen, Physical Review B 60 p.11971 (1999)

B. L. Holian et al. Phys. Rev. A 43 p.2655 (1991)

VED Research Group home page link

Daniel Murray
Associate Professor
Math, Stats and Physics Unit
University of British Columbia Okanagan Kelowna, BC, Canada
Email: daniel DOT murray AT ubc DOT ca

Appendix: Material Properties
Information on material properties can be found at Matweb.

Material	Density	Young's Modulus	sqrt(Y/ρ)	ν	structure
	(g/cm³)	(GPa)	(m/s)
Aluminum	2.7	73	5200	0.33	fcc
Acrylic	1.17-1.20	3.2	1650
Brass	8.4	130	3930
Celluloid	1.4	1.3-1.5
Cellulose acetate	1.22-1.34	1.5-2.2	--
Copper	8.93			0.31	fcc
Diamond	3.516			0.1-0.29	diam.
Glass	2.60	46.2	4220
Gold	19.28			0.42	fcc
Iron	7.87			0.291	bcc
Lead	11.34			0.42	fcc
LDPE				0.38
nylon (polyamide)	1.1	1.8	--
nylon 6	1.1	3.0	1650	0.35
polystyrene	1.04 -1.09	1.5-3.5	--	0.33
polystyrene, impact				0.33-0.37
rubber (nat. vulc.)	0.95	0.00276	54
Silver	10.50			0.37	fcc
Steel	7.8	210	5190	0.29
Teflon				0.46

Poisson ratios are from www.matweb.com
Structures are from C. Kittel "Introduction to Solid State Physics", Wiley 1971

Hosted by www.Geocities.ws