comparison

COMPARISON OF THE LATTICE-DYNAMICSTHERMODYNAMIC PROPERTIES

Several paths are available for calculating macroscopic thermodynamic properties for a system with given interparticle forces. Most paths are approximate. The lattice dynamics approximationl) treats correctly all terms in the energy which are quadratic in the particle displacements. The cell-model approximation2) in which a single particle moves in the field of its fixed neighbors modifies the quadratic terms but includes an estimate of anharmonic' corrections. More sophisticated theories3) treat anharmonic perturbations analytically (the actual calculations require fast computers) but seem complicated enough to attract few follo\vers.
The approximate methods have the important advantage of being quick and inexpensive to calculate. The computer experiments giving exact thermodynamic properties, either by following the motion of the particles or by sampling the configuration space 4), are relatively expensive because they require so much computer time, particularly if high precision is necessary. To halve the statistical uncertainty in computer-experiment results requires quadrupling the computer time used up. The precision also depends on the sensitivity of the property measured to fluctuations in pressure and energy. Because successively higher derivatives of the free energy involve higher moments of the pressure-tensor component distributions and energy distributions, the time necessary to characterize derivatives increases rapidly with derivative order. Second-order elastic constants and the specific heat involve second moments; third-order elastic constants involve third moments; and so on. So far only first-and second-derivative quantities have been examined. Thus, if the approximations should prove to come close to exact results, they would give us a useful shortcut to accurate thermodynamic properties. It was in the hope of establishing their usefulness that we undertook these calculations.
In this paper we compare the results of solid-phase Monte Carlo experiments on 108 particles interacting with the Lennard-Jones and exponential-six potentials with the predictions of 108-particle lattice dynamics and the cell model. In addition to the energy and pressure, we compare the secondderivative quantities: specific heat, Griineisen y, and elastic constants, with approximate predictions. Although either computer method, molecular dynamics or Monte Carlo, can be extended to quantum systems by using the \Vigner-Kirkwood Planck's constant expansion of the free energy5), we have made our comparisons using classical mechanics. To find out how important the small size and classical nature of our systems are, we use the latticedynamics method to investigate the number dependence of all of the thermodynamic properties and quantum corrections to the elastic constants. Quan-t tum corrections to other thermodynamic properties have been calculated elsewhere6).
In section 2 we describe the lattice-dynamics calculations. The method has been in use for over 50 years, although many so-called lattice-dynamics calculations of elastic constants have actually been calculations of the elastic response of a static lattice. We take lattice vibrations explicitly into account. The first correct harmonic calculation using normal-mode vibrations was announced by Feldman7). He obtained expressions for the second-order elastic constants which involved the vibration frequencies and their strain derivatives. The frequencies were then obtained by the usual method of diagonalizing the dynamical matrix. The frequency-shift derivatives were calculated by means of perturbation theory. We calculate the work of deforming the crystal by the alternative procedure of computing the free energy numerically for several slightly different values of the strain and then fitting the results to a strain polynomial. Besides avoiding the tedious algebra of Feldman's analytic approach, the numerical method is more readily generalized to higher-order elastic constants.
In section 3 we describe the cell-model calculations. The cell model, although actually only appropriate for solids, was first used in an attempt to describe gases and liquids 8). In the cell model the effect of heating the crystal lattice is approximated by a "one-particle" model in which a single particle moves in the field of its fixed neighbors. Neglecting interparticle correlations by approximating an N-body problem by a one-body problem for classical systems, most reasonable at low temperatures and is exact only in the static-lattice limit.
The advantage of the cell model over the harmonic approximation lies in its ability to estimate anharmonic contributions from potential-energy terms beyond the quadratic ones. The cell model has often been used to calculate energy, pressure, and specific heat 9). Our elastic-constant calculations are a new use of this model. Henkel10) has studied a similar model, a quantum cell model in which the potential was expanded in powers of displacement and the quartic contributions were treated by perturbation theory. If the perturbations were ignored, Henkel's work would reduce to the usual harmonic Einstein model.
In section 4 we compare the tabulated results from both approximations and consider the dependence of the results on number of particles. Vie also discuss some interesting cancellations found in the course of the Monte Carlo calculations. In section 5 we assess the importance of quantum effects on the elastic constants.
2. Lattice-dynamics calculations.
The lattice dynamics calculations are based on the approximation of truncating a Taylor expansion of the lattice potential energy after the quadratic terms in the particle displacements. Sometimes this is called the "quasi-harmonic" approximation. The coefficients in the Taylor series expansion are calculated from the assumed force law, and the lattice sites are chosen to match the structure of the lattice being described. The expansion is made about a configuration in which each atom is fixed at its average position in a perfect crystal with fixed center of mass. Thus the truncated potential depends upon the size and shape of the assumed static lattice configuration. By changing to normal-mode coordinates the quasi-harmonic Hamiltonian can be rewritten as a sum of 3N 3 independent harmonic oscillator Hamiltonians. The partition function of an oscillator, either quantum or classical, is known11) so that the quasi-harmonic thermodynamic properties of the system can be calculated.
The partition function, Z = exp(-A/kT), where A is the Helmholtz energy and kT is Boltzmann's constant times the absolute temperature, can be written as a product of single oscillator partition functions:

In the classical limit only the first term in the expansion, kT/hv, is kept. from the quadratic Hamiltonian; the center of mass contribution is Zcm. For accurate work on small crystals Zcm has to be included if comparisons are made with Monte Carlo calculations in which the center of mass is allowed to move.
Thermodynamic properties can all be derived from the partition function. Temperature derivatives can be evaluated explicitly to compute the energy and specific heat:

"Strain" derivatives are harder to evaluate. The strains are defined in terms of three vectors colinear with the edges common to one corner of a parallelepiped produced by deforming a cube of crystal with sidelength a. If aI, a2, and a3 are the vectors, then

are the six independent strains. The dynamical matrix then gives a complicated implicit relation for the vibration frequencies as functions of the strains. Because there is no convenient expression for v(1]) it is easiest to proceed numerically. Both the average pressure tensor component
where fi and ri are the Cartesian coordinates of rand r, and 1]1 is the only nonzero strain. In the strained configuration the Hamiltonian is again expanded, the quadratic terms kept, and the result diagonalized, giving a new set of frequencies and the free energy A (1]1). The diagonalization can be visualized in terms of plane-wave solutions of the classical equations of motion. The waves propagate through the crystal with wavelengths and directions imposed by the shape of the crystal and described by wave vectors chosen from a convenient Brillouin zone12). If y is a wave vector in the unstrained crystal then the corresponding y in the strained crystal is:

where Yi and Yi are the Cartesian coordinates of y and y. Applying the Born-von Karman normal-mode analysis for several values of 1]1, the first four or five coefficients in the expahsion,

can be determined. By choosing values of 1J1 separated by 0.0001, both the pressure and Cil were determined with four-figure accuracy in this way. By considering two simultaneous strains along the one axis and the two axis, and analyzing the resulting free energy changes as a series in 1]1 and

The generalization of this technique to calculation of higher-order elastic constants or to mixed strain-temperature derivatives is straightforward. For crystals of lower symmetry one needs to calculate more elastic constants and hence one considers more different combinations of strain. Other strains would also be needed for cubic crystals to determine the six nonzero thirdorder constants or the eleven nonzero fourth-order constants and the mixed strain-temperature derivatives.The adiabatic elastic constants

where 5 is the entropy of the crystal. Because 5(1]1) and (oT/o5)n can be calculated exactly for a quasi-harmonic crystal, the correction term can be evaluated numerically.In order to make a comparison with the 108-particle Monte Carlo calculations13, 14) we have calculated the thermodynamic properties for a 108particle system with the same periodic boundaries and Hamiltonian as those used in the Monte Carlo work. Using the nearest-image convention, each particle in the crystal interacts with 107 neighbors according to the LennardJones 6-12 potential

These potentials have both been used principally to describe rare gases and have shown themselves to fit these rea(materials well. Although we picked these potentials because of their value in describing rare gases we expect

that our general conclusions in comparing approximate calculations with exact computer experiments will be valid for potentials describing interactions in salts or metals as well.Table I gives both the Monte Carlo and the lattice-dynamics results. The three temperatures span the range from about 0.48TtriPIO to O.95TtriPle and the densities correspond closely to zero pressure. Because the static-lattice contributions to the thermodynamic properties present no theoretical problems (they are correctly calculated by any theory) we have tabulated separately the thermal contributions to the thermodynamic properties. The data in the table show that the lattice-dynamics elastic constants are quite· close to the Monte Carlo values at three different temperatures and for both potentials tested.