J. Phys. Chem. B 2008, 112, 1699-1705
1699
Elastic Properties of Soft Sphere Crystal from Monte Carlo Simulations Konstantin V. Tretiakov* and Krzysztof W. Wojciechowski Institute of Molecular Physics, Polish Academy of Sciences, Smoluchowskiego 17/19, 60-179 Poznan´ , Poland ReceiVed: August 29, 2007; In Final Form: NoVember 21, 2007
Elastic properties of faced centered cubic (fcc) crystals composed of soft spheres, interacting through potentials of the form u(r) ∼ r-n, have been investigated by Monte Carlo (MC) simulations. It is shown that both the softness parameter (n-1) and temperature strongly influence the elastic properties of the studied system. The simulations show explicitly that when T > 0 the elastic constants of the hard sphere crystal can be obtained by taking the limit n f ∞ of soft spheres. When T f 0 for any finite n, the elastic constants of the soft spheres tend to those of the static model. At all temperatures and softness parameters studied here, n, the Poisson’s ratio in [110]⊥ direction is negative.
I. Introduction Various physical systems and phenomena have been described using the inverse power potential1-36
u(r) )
(σr )
n
(1)
where r is the separation between two particles, σ is the particle diameter, sets the energy scale, and the exponent n is a parameter that determines the hardness of the potential (softness is ∼1/n). Changing the exponent n, one can describe different physical systems from the one component plasma15 (n ) 1), through liquid alkali metals31 (small n), melting of aluminum, molybdenum, and light actinides30 (n ) 6), molecular liquids and solids1,4,6,8,9,12,14,16,17,21-26,28 (n ) 6, 12), and granular materials, powders, and colloids27,33 (large n), to the hard spheres37-39 (n f ∞). The inverse power potential is particularly useful for perturbation theories of the equation of state of classical fluids and solids.3,8,40,41 Structural properties of various soft sphere systems (e.g., radial distribution functions for n ) 4, 6, 9, 12 soft sphere fluids) have been investigated in detail in refs 6 and 12. The equation of state for systems interacting via purely repulsive inverse power potential has been studied extensively by computer simulations both for finite values of n and in the limit n f ∞.1,7,15,27,37 Recently, virial coefficients up to the seventh one have been calculated for such systems.32 Influence of the softness of the pair potential on the melting transition and stability of different phases has been investigated in pioneering works of Hoover4,5 and others.14,16,34,36,42,43 It has been found that, depending on the softness, the system interacting via the n-inverse-power potential freezes either into the face-centered cubic (fcc) crystal structure or into the body-centered cubic (bcc) one. This softness-driven transition has been extensively studied by Hoover5 and others.16,42 Specifically, Agrawal and Kofke43 have shown that the fcc structure is the thermodynamically favored solid phase for soft spheres when n > 7; however, for n < 6, the system freezes into the bcc structure. Systems of particles interacting through inverse power potentials are also of increasing interest within the colloid and * Corresponding author. E-mail:
[email protected].
interface science community. One of the reasons is that particles of different softness can be used for various applications. Hence, not only phase diagrams and structures of soft sphere systems are of interest, but also a comprehensive description of the thermodynamic, mechanical, and transport properties such as that performed by Heyes and co-workers for r-n fluids.17,18,21-26,28 In particular, transport properties of soft sphere liquids were determined in refs 18 and 22-26. The thermal conductivity was also studied for soft sphere solids13 with n ) 12. One should stress that the self-similarity of the n-inversepower-potential, that is, the lack of any characteristic length scale in systems interacting through this potential, makes it an interesting subject in the context of classical elasticity theory, which also does not exhibit any characteristic length scale. Elastic moduli and the viscoelastic properties of r-n fluids have been calculated in refs 17, 21-23, 26, and 29. Elastic properties of the limiting case n f ∞ have been also extensively studied.39,44-48 Elastic properties of soft sphere solids have been investigated only very recently by the present authors,39 who compared the low-temperature limit with the high-temperature results, and by Branka and Heyes,49 who studied the system at two densities. Here, we extend the previous results by changing the temperature smoothly from small to large values. In the present work, we study the elastic properties of threedimensional soft spheres in the fcc phase. This work is a part of a larger project which concentrates on investigation of the Poisson’s ratio in various model systems. We expect that studies of simple and well-defined models will lead to a better understanding and description of a new class of materials50-56 which exhibit anomalous (negative) Poisson’s ratio.57 Such unusual materials are of interest both for fundamental research and for applications. Searching for mechanisms which can decrease the Poisson’s ratio of model systems can help in better understanding, modeling, and designing new materials with anomalous (negative) Poisson’s ratios. The aim of this paper is threefold. First, we investigate the influence of the exponent n on the elastic properties of soft sphere crystalline phases. Second, we present the simulation results of the soft spheres to show explicitly that the elastic properties of hard spheres correspond to n f ∞ and T > 0, whereas the static (zero temperature) systems are obtained in
10.1021/jp076929o CCC: $40.75 © 2008 American Chemical Society Published on Web 01/24/2008
1700 J. Phys. Chem. B, Vol. 112, No. 6, 2008
Tretiakov and Wojciechowski
Figure 1. Softness (1/n) dependence of (a) the bulk modulus, (b) elastic constant C/11, (c) the elastic constant C/12 and (d) the elastic constant C/44 at constant pressure p* ) 37.69.
Figure 2. Temperature dependence of the bulk modulus (a), elastic constant C/11 (b), elastic constant C/12 (c), and elastic constant C/44 (d), at pressure p* ) 37.69. The solid circles represent results of the hard sphere system obtained in ref 39.
the limit T f 0. Third, we determine explicitly the influence of the particle motions (i.e., positive temperature) on the Poisson’s
ratio of the hard sphere system. Since in the simulations described, in the present paper n g 12, all of the simulations
Elastic Properties of Soft Sphere Crystal
J. Phys. Chem. B, Vol. 112, No. 6, 2008 1701 the [110] direction then | is [001] direction and ⊥ is [1h10] direction. By using the following notation,
1 B ) (C11 + 2C12 + p) 3
(4)
1 µ1 ) (C11 - C12 - 2p) 2
(5)
µ2 ) C44 - p
(6)
then in the case of an fcc crystal under pressure, the Poisson’s ratio is equal to Figure 3. Pressure dependence of the Poisson’s ratio of hard sphere system along three symmetry directions in crystal. Open symbols correspond to results from ref 48.
were performed for the face-centered cubic (fcc) crystalline phase which is stable for such values of n.43 The structure of the paper is as follows. In section II, we outline the simulation details. In section III, we present and discuss the results of simulations. The results are compared with those obtained both for a static lattice model interacting by the n-inverse-power potential and for the hard sphere model. The last section (IV) contains a summary and conclusions. II. Simulation Details The Monte Carlo simulations were performed in the NPT ensemble by a method following the Parrinello-Rahman idea of averaging strain fluctuations,58 which has been further developed in refs 59 and 60. The version of this method applied here is based on refs 61-63. The cubic symmetry of the fcc structure implies that the strain expansion of the free enthalpy (Gibbs free energy, G) reads
1 G ) G(0) + B11Vp(εxx2 + εyy2 + εzz2) 2 + B12Vp(εxxεyy + εyyεzz + εxxεzz) + 2B44Vp(εxy2 + εyz2 + εxz2)
(2)
where Vp is the volume of the reference state which corresponds to the equilibrium state at p; εii are the components of strain tensor, and the elastic constants Bij and Cij fulfill the relations: 64,65
B11 ) C11 - p
B12 ) C12 + p
B44 ) C44 - p (3)
In present paper, we used the Voigt’s notation where the indices 1-6 are used for xx, yy, zz, xy, xz, yz, respectively. For an isotropic solid, the Poisson’s ratio is defined as the negative ratio of the transverse strain to the longitudinal strain when the stress parallel to the longitudinal direction is changed.57 In the case of an anisotropic solid, the Poisson’s ratio depends on both the direction along which the increase in stress is made and the direction along which the resulting change is measured. If the stress is increased along kˆ , then it is sufficient to know the Poisson’s ratio for two mutually orthogonal directions, both of which are perpendicular to kˆ . In this way, one can obtain the Poisson’s ratio along any other direction associated with this change using the standard rules for rotating tensors. All details are discussed in Wallace.64 The two mutually orthogonal directions will be denoted by | and ⊥. For example, if kˆ is in
ν[100] )
3B - 2µ1 6B + 2µ1
(7)
ν[111] )
3B - 2µ2 6B + 2µ2
(8)
ν[110]| ) ν[110]⊥ )
6Bµ2 - 4µ1µ2 3B(3µ1 + µ2) + 4µ1µ2 3B(3µ1 - µ2) - 4µ1µ2 3B(3µ1 + µ2) + 4µ1µ2
(9)
(10)
Because of the cubic symmetry of the system, ν does not depend on the transverse direction for the [100] and [111] longitudinal directions (i.e., ν(|) ) ν(⊥)). Monte Carlo simulations were performed for particles interacting via the inverse power potential (1) for several values where n ) 12, 16, 24, 48, 96, 192, 384, 768. All of the simulations have been done in the reduced units of the energy E* ) E/, the density F* ) Fσ3 ) Nσ3/V, and the temperature is T* ) kBT/. Thus, the dimensionless pressure is p* ) pσ3/, the dimensionless bulk modulus is B* ) Bσ3/, and the dimensionless elastic constants are C/ij ) Cijσ3/. A standard interaction cutoff of 2.5σ for n e 48 and 2.0σ for n > 48 was used. Two kinds of trial motions were used. The first, concerned changes of the sphere positions and its acceptance ratio was kept close to 30%. The second type of motion corresponded to changes of the components of the symmetric box matrix and was tried about N1/2 times less frequently than the sphere motions. These box motions determined the size and the shape of the box, and their acceptance ratio was close to 20%. We simulated systems of 256 particles with a periodic boundary condition for which the T ) 0 ground state configuration is the fcc lattice occupying a cubic box of side 4x2a0 (a0 is the nearest-neighbor distance). In ref 39, it has been shown that simulations of systems as small as N ) 256 give the elastic constants and the Poisson’s ratio which differ by only a few percent from the results obtained by extrapolation to the N f ∞ limit. Typical lengths of the runs were 5 × 106 trial steps per particle (Monte Carlo cycles) after an equilibration of 106 MC cycles. III. Results and Discussion A. Static Limit (T ) 0). The question concerning the influence of the particle motions (i.e., positive temperature) on the elastic properties of the soft sphere system can be answered by considering a static, that is, zero temperature, fcc
1702 J. Phys. Chem. B, Vol. 112, No. 6, 2008
Tretiakov and Wojciechowski
Figure 4. Poisson’s ratio versus the softness at p* ) 37.69: (a) in [100] direction, (b) in [111] direction, (c) in [110]| direction, and (d) in [110]⊥ direction.
Figure 5. Temperature dependence of the Poisson’s ratio at p* ) 37.69: (a) in [100] direction, (b) in [111] direction, (c) in [110]| direction, and (d) in [110]⊥ direction.
lattice whose nearest-neighboring sites (distanced by a) interact by the potential (eq 1). The pressure, bulk modulus, elastic
constants, and the Poisson’s ratio of such a lattice are the following:
Elastic Properties of Soft Sphere Crystal
pstatic ) 23/2 na-(n+3)
(3n + 1)p n ) ( + 1)p 2
J. Phys. Chem. B, Vol. 112, No. 6, 2008 1703
(11)
Bstatic )
static
(12)
Cstatic 11
static
(13)
1 static Cstatic 12 ) C11 2
(14)
1 ) Cstatic Cstatic 44 2 11
(15)
n+6 3n + 6
(16)
n+8 5n + 10
(17)
νstatic [100] ) νstatic [111] ) static ) ν[110]|
static ν[110]⊥ )
(n - 2)(n + 6) 2n2 - 5n - 18
n + 18 -2n2 + 5n + 18
(18) (19)
The Poisson’s ratios of the static model in the limit n f ∞ are equal ν[100] ) 1/3, ν[111] ) 1/5, ν[110]| ) 1/2, and ν[110]⊥ ) 0. It is worth noticing that the Poisson’s ratio in the [110]⊥ direction is negative for n > 4.5 as follows from eq 19; however, only for n > 7 can the thermodynamically stable fcc structure be observed.43 B. Bulk Modulus and the Elastic Constants. In Figure 1, we show the dependence of the dimensionless bulk modulus and the elastic constants on the softness parameter, 1/n, for different temperatures 0.001 e T* e 1 at the pressure p* ) 37.69. It is easy to see in Figure 1 that the results for T* )
0.001 are very close to the static case described by eqs 12-15. It can be also seen that in the limit n f ∞ the elastic constants reach the values corresponding to those of the hard sphere crystal at the given pressure and (a positive) temperature. The appropriate values are described in Figure 1 by arrows and symbols HS. One should stress that the n dependence of the C12 is nonmonotonic which has an influence on the Poisson’s ratio. Analyzing the softness dependence (Figure 1) of elastic constants, we can say that the system with a given n becomes harder when temperature is decreased. Temperature dependencies of the bulk modulus and the elastic constants are presented in Figure 2. One finds that, at high temperatures when 0.1 < T* e 1.0 and n f ∞, both the bulk modulus and the elastic constants of soft spheres tend toward results obtained for hard spheres at the same pressure. It is wellknown that the hard sphere model is an athermal system; therefore, such a comparison is made possible by introducing the following temperature for the hard sphere system
T* )
p* p/HS
(20)
where p/HS ) pσ3/kBT is dimensionless pressure of hard spheres and p* ) pσ3/ is dimensionless pressure of a soft sphere system. From the temperature dependencies of the elastic constants (Figure 2), we conclude that for some fixed temperature an increase of the exponent n in the potential makes the system more stiff from the point of view of its elastic properties. C. Poisson’s Ratio. At the beginning of the section, we present the Poisson’s ratio of hard spheres along the main symmetry axis in Figure 3. We find very good agreement of our results with literature data. It is worth noting that in the [110]⊥ direction the Poisson ratio is negative. Recently, the Poisson’s ratio of the soft spheres in the [100] direction has been calculated.66 Here, we present the analysis of the Poisson’s ratio of soft spheres in all crystallographic
Figure 6. Poisson’s ratio versus the inverse temperature-scaled pressure: (a) in [100] direction, (b) in [111] direction, (c) in [110]| direction, and (d) in [110]⊥ direction.
1704 J. Phys. Chem. B, Vol. 112, No. 6, 2008 directions. In Figure 4, the Poisson’s ratio of the soft sphere systems at the pressure p* ) 37.69 is plotted versus the softness parameter 1/n. We find that the qualitative behavior of the Poisson’s ratio in the [110]⊥ direction (Figure 4d) is completely different than in other directions. The [110]⊥ Poisson’s ratio grows with exponent n, whereas in other directions, we observe decreasing of the Poisson’s ratio with this parameter. One can see that the limiting values for n f ∞ tend to the Poisson’s ratio obtained for the hard sphere system at the same pressure which is marked by an arrow in Figure 3. In the case when T* f 0, the values of the Poisson’s ratio tend toward the values of the static model. It is also worth noting that for large n > 384 in the range 0.25 e T* e 0.75 the values of the Poisson’s ratio are lower than in the hard sphere limit. In the latter case, the Poisson’s ratio is lower than that for the static system. A similar effect has been observed for the soft disk system in two dimensions.67 The temperature dependence of the Poisson’s ratio is shown in Figure 5. It is easy to notice that an increase of the power n in the interaction potential leads to a decrease of the Poisson’s ratio for soft spheres in the whole temperature range considered. It can also be seen that by increasing the particle motions (i.e., increasing the temperature) one can decrease the Poisson’s ratio of soft spheres with respect to the static case (i.e., to the zero temperature limit) for n J 48 in [100], [111], and [110]| directions. In the [110]⊥ direction, we observe a different and more complicated dependence (Figure 5d). At all temperatures and for each studied n, the Poisson’s ratio in [110]⊥ is negative. For very soft spheres when n ) 12, 16, the Poisson’s ratio in the [110]⊥ direction decreases monotonically with increasing temperature. A particularly interesting conclusion which follows from analysis of Figure 5 is that when n f ∞ there is a “jump” of the Poisson’s ratios in directions [100], [111], [110]|, and [110]⊥ between its value 1/3, 1/5, 1/2, and 0 obtained for T ) 0 and the values 0.2, 0.044, 0.26, and -0.07 obtained for T > 0, respectively. Comparison of the elastic properties of soft spheres with the corresponding hard sphere limit can be made clearly using a common variable, such as temperature (cf. Figure 5). It is known, however, that the soft sphere system may also be represented in a scaled form, in which the temperature can be scaled out of the problem. Specifically, the temperature-scaled pressure can be expressed as p˜ ) p*/(T*)1+(3/n) ) (pσ3/)(/ kBT)1+(3/n).49 In Figure 6, we present the Poisson’s ratio as a function of the inverse temperature-scaled pressure. Comparing these alternative representations, one finds that the inverse temperature-scaled pressure dependencies (Figure 6) of the Poisson’s ratio look very similar to the temperature ones (Figure 5). IV. Summary and Conclusions Elastic properties of the soft sphere system were determined by analysis of the box matrix evolution in Monte Carlo simulations. The elastic constants and Poisson’s ratio were determined for a wide range of the softness parameter (1/n) of the soft sphere potential and for a wide range of temperatures. These data can be used to construct temperature-softness dependence of these properties. Simulations of the soft spheres indicate that at T ) 0.001 the difference between elastic constants for that model and for the static model is less than two percent when n g 12. On the other hand, simulations of the soft spheres at T* ) 1 show that
Tretiakov and Wojciechowski the elastic constants of hard spheres differ by less than two percent from those of the soft spheres system when n g 384. We have shown that particle motions decrease the Poisson’s ratio of the soft sphere system with respect to the static case (i.e., to the zero temperature case) for n J 48. We have also found that by increasing the exponent n in the interaction potential one can decrease the Poisson’s ratio in [100], [111], and [110]| directions at high pressures/low temperatures. In a separate work, we will show that in some systems combination of these two effects leads to negative Poisson’s ratios. Acknowledgment. K.V.T. thanks K. J. M. Bishop (Northwestern University) for kindly reading and commenting on the manuscript. Part of this work was supported by Grant N20207032/ 1512 (2007-2010) of the Polish Ministry of Science and Higher Education. Part of the calculations were performed at the Poznan´ Computing and Networking Center (PCSS). References and Notes (1) Hoover, W. G.; Ross, M.; Johnson, K. W.; Henderson, D.; Barker, J. A.; Brown, B. C. J. Chem. Phys. 1970, 52, 4931. (2) Hansen, J.-P. Phys. ReV. A 1970, 2, 221. (3) Andersen, H. C.; Weeks, J. D.; Chandler, D. Phys. ReV. A 1971, 4, 1597. (4) Hoover, W. G.; Gray, S. G.; Johnson, K. W. J. Chem. Phys. 1971, 55, 1128. (5) Hoover, W. G.; Young, D. A.; Grover, R. J. Chem. Phys. 1972, 56, 2207. (6) Hansen, J.-P.; Weis, J. J. Mol. Phys. 1972, 23, 853. (7) Cape, J. N.; Woodcock, L. V. Chem. Phys. Lett. 1978, 59, 271. (8) Ross, M. J. Chem. Phys. 1979, 71, 1567. (9) Cape, J. N.; Woodcock, L. V. J. Chem. Phys. 1980, 72, 976. (10) Broughton, J. Q.; Gilmer, G. H.; Weeks, J. D. J. Chem. Phys. 1981, 75, 5128. (11) Broughton, J. Q.; Gilmer, G. H.; Weeks, J. D. Phys. ReV. B 1982, 25, 4651. (12) Young, D. A.; Rogers, F. J. J. Chem. Phys. 1984, 81, 2789. (13) Ladd, A. C.; Moran, B.; Hoover, W. G. Phys. ReV. B 1986, 34, 5058. (14) Laird, B. B.; Haymet, A. D. J. J. Chem. Phys. 1989, 91, 3638. (15) Stringfellow, G. S.; DeWitt, H. E.; Slattery, W. L. Phys. ReV. A 1990, 41, 1105. (16) Laird, B. B.; Haymet, A. D. J. Mol. Phys. 1992, 75, 71. (17) Heyes, D. M.; Aston, P. J. J. Chem. Phys. 1994, 100, 2149. (18) Heyes, D. M. J. Phys.: Condens. Matter 1994, 6, 6409. (19) Rothen, F.; Pieran´ski, P. Phys. ReV. E 1996, 53, 2828. (20) Wojciechowski, K. W.; Klos, J. J. Phys. A: Math. Gen. 1996, 29, 3963. (21) Heyes, D. M. J. Chem. Phys. 1997, 107, 1963. (22) Heyes, D. M.; Powles, J. G. Mol. Phys. 1998, 95, 259. (23) Powles, J. G.; Heyes, D. M. Mol. Phys. 2000, 98, 917. (24) Heyes, D. M.; Powles, J. G. Mol. Phys. 2001, 99, 1077. (25) Heyes, D. M.; Powles, J. G.; Rickayzen, G. Mol. Phys. 2002, 100, 595. (26) Rickayzen, G.; Powles, J. G.; Heyes, D. M. J. Chem. Phys. 2003, 118, 11048. (27) Bran´ka, A. C.; Heyes, D. M. Mol. Phys. 2004, 102, 2049. (28) Bran´ka, A. C.; Heyes, D. M. Phys. ReV. E 2004, 69, 021202. (29) Bran´ka, A. C.; Heyes, D. M. Comput. Method. Sci. Technol. 2004, 10 (2), 127. (30) Ross, M.; Yang, L. H.; Boehler, R. Phys. ReV. B 2004, 70, 184112. (31) Ben-Amotz, D.; Stell, G. J. Chem. Phys. 2004, 120, 4844. (32) Wheatley, R. J. J. Phys. Chem. B 2005, 109, 7463. (33) Bran´ka, A. C.; Heyes, D. M. Mol. Phys. 2005, 103, 2359. (34) Cardenas, M.; Tosi, M. P. Phys. Lett. A 2005, 336, 423. (35) Hynninen, A.-P.; Dijkstra, M. Phys. ReV. Lett. 2005, 94, 138303. (36) Davidchach, R. L.; Laird, B. B. Phys. ReV. Lett. 2005, 94, 086102. (37) Hoover, W. G.; Ree, F. H. J. Chem. Phys. 1968, 49, 3609. (38) Frenkel, D.; Ladd, A. J. C. J. Chem. Phys. 1984, 18, 3188. (39) Tretiakov, K. V.; Wojciechowski, K. W. J. Chem. Phys. 2005, 123, 074509. (40) Kang, H. S.; Lee, C. S.; Ree, T. J. Chem. Phys. 1985, 82, 414. (41) Kang, H. S.; Ree, T.; Ree, F. H. J. Chem. Phys. 1986, 84, 4547. (42) Ogura, H.; Matsuda, H.; Ogawa, T.; Ogita, N.; Ueda, A. Prog. Theor. Phys. 1977, 58, 419. (43) Agrawal, R.; Kofke, D. A. Phys. ReV. Lett. 1995, 74, 122. (44) Hoover, W. G.; Salsburg, Z. W. J. Chem. Phys. 1971, 54, 2129.
Elastic Properties of Soft Sphere Crystal (45) Frenkel, D.; Ladd, A. J. C. Phys. ReV. Lett. 1987, 59, 1169. (46) Runge, K. J.; Chester, G. V. Phys. ReV. A 1987, 36, 4852. (47) Farago, O.; Kantor, Y. Phys. ReV. E 2000, 61, 2478. (48) Pronk, S.; Frenkel, D. Phys. ReV. Lett. 2003, 90, 255501. (49) Bran´ka, A. C.; Heyes, D. M. Mol. Simul. 2005, 13, 937. (50) Lakes, R. AdV. Mater. 1993, 5, 293. (51) Baughman, R. H.; Shacklette, J. M.; Zakhidov, A. A.; Stafstrom, S. Nature 1998, 392, 362. (52) Evans, K. E.; Alderson, A. AdV. Mater. 2000, 12, 617. (53) Evans, K. E.; Alderson, K. L. Eng. Sci. Ed. J. 2000, 4, 148. (54) Wojciechowski, K. W. In Properties and Aplications of Nanocrystalline Alloys from Amorphous Precursors; Idzikowski, B., Sˇ vec, P., Miglierini, M., Eds.; Kluwer Academic Publishers: Dordrecht, 2005; pp 241-252. (55) Konyok, D. A.; Wojciechowski, K. W.; Pleskachevskii, Y. M.; Shilko, S. V. Mekh. Kompoz. Material. Konstr. 2004, 10, 35, in Russian. (56) See, e.g., http://www.ifmpan.poznan.pl/zp10/auxet/ main.html.
J. Phys. Chem. B, Vol. 112, No. 6, 2008 1705 (57) Landau, L. D.; Lifshits, E. M.; Kosevich, A. M.; Pitaevskii, I. P. Theory of Elasticity; Pergamon Press: London, 1986. (58) Parrinello, M.; Rahman, A. J. Chem. Phys. 1982, 76, 2662. (59) Ray, J. R.; Rahman, A. J. Chem. Phys. 1984, 80, 4423. (60) Ray, J. R.; Rahman, A. J. Chem. Phys. 1985, 82, 4243. (61) Wojciechowski, K. W.; Tretiakov, K. V. Comput. Phys. Commun. 1999, 121-122, 528. (62) Wojciechowski, K. W. Comput. Method. Sci. Technol. 2002, 8 (2), 77. (63) Wojciechowski, K. W.; Tretiakov, K. V.; Bran´ka, A. C.; Kowalik, M. J. Chem. Phys. 2003, 119, 939. (64) Wallace, D. C. In Solid State Physics; Ehrenreich, H., Seitz, F., Turnbull, D., Eds.; Academic Press: New York, 1970; p 301. (65) Wojciechowski, K. W. Mol. Phys. Rep. 1995, 10, 129. (66) Tretiakov, K. V.; Wojciechowski, K. W. Mater. Sci. (Poland) 2007, 25, 541. (67) Tretiakov, K. V.; Wojciechowski, K. W. ReV. AdV. Mater. Sci. 2007, 14, 104.
