- Email: [email protected]
Parallelization in classical molecular dynamics simulation and applications
Computational Materials Science 37 (2006) 146–151 www.elsevier.com/locate/commatsci
Parallelization in classical molecular dynamics simulation and applications S.L. Chaplot Solid State Physics Division, Bhabha Atomic Research Centre, Mumbai 400 085, India
Abstract Over the last four decades, molecular dynamics simulation has been a well-established technique to study a wide variety of physical phenomena. The technique basically involves solving the equations of motion of a system of particles (e.g., atoms or even stars interacting via a potential), and following their trajectories. At a microscopic level, MD simulations have been extensively used in the study of the structures and dynamics, phase transitions and thermodynamic properties of solids and liquids. Simulations of thermal and highpressure processes (often beyond laboratory conditions), shock-stress propagation and defect dynamics have provided unique insights into material properties. This review paper will briefly discuss the parallelization of our in-house code and its uses at the BARC parallel computer. Ó 2006 Elsevier B.V. All rights reserved. Keywords: Molecular dynamics simulation; Thermal properties; High-pressure processes, Shock propagation
1. Introduction Molecular dynamics simulation (MDS) is the name of a technique in which the trajectories of a system of interacting particles are calculated by numerically solving their classical equations of motion [1a]. A very wide range of applications, e.g., the study of phase transitions and dynamics of microscopic defects in solids, has been developed over the past few decades. In fact, this rapid development has kept pace with the growing availability of computing power. An essential requirement for solving the classical equations of motion is the ability to calculate the forces on the particles. There are many ways to do this. The forces may be obtained from an interparticle potential; or from the force constants in a coupled system. Alternatively, the forces in an atomic system could be evaluated from the first principles quantum mechanical calculation of the electronic energies [2]. Accordingly, the MDS may be called
E-mail address: [email protected] 0927-0256/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.12.013
the classical MDS or the ab initio MDS. The two are complementary while the former is useful for extensive studies and for relatively more complex systems. Starting with the early simulations by Alder and Wainright [1a] using hard-core potentials and Rahman [1a] using Lennard–Jones potentials, there has been a great deal of progress in terms of various applications in the last four decades. However, the progress was relatively slow in the sixties and seventies as extensive computing facilities were not widely available. In late seventies we studied the defect vibrational modes in rare gas solids and identified them as local modes or the resonance modes [1b]. A major milestone was the development of the constant-pressure technique by Parrinello and Rahman [1c], which facilitated studies of phase transitions involving changes in the crystallographic unit cells. Simulations relevant to materials science involve study of a variety of physical systems, such as the bulk solids and liquids, clusters and nano-particles, surfaces and interfaces, and solid lattices with point and extended defects, etc. The simulations may involve the characteristic parameters of temperature (T), volume (V), total energy (E), external
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
stress (S), total number of particles (N), etc. The quantities like T, S, V, and E, are easily calculated at any instant of time using the atomic positions and velocities. In turn, any of these quantities (T, S, E, V) can be controlled by suitably constraining the atomic trajectories in the simulations [1c]. Simulations commonly employed are those in which various combinations, e.g., (NVT), (NPT), (NVE), or (NE) are kept fixed. It is now convenient to perform simulations with N about 10 000 to more than a million atoms. Bulk properties are generated by assuming periodic boundary conditions. However, the time scale of the microscopic MDS remains limited to about a few nano-seconds, i.e., a few million time steps of a few femto seconds each. This time sale is sufficient for a large variety of dynamical phenomena at the atomic level including critical propagation of extended defects and many of the phase transitions. It is easy to make a contact between the MDS and the various scattering experiments, such as the Raman and neutron scattering experiments from various materials. Such experiments essentially measure the correlations among the positions of atoms at the same time or at different times. These correlations can be directly calculated from the simulated atomic trajectories and compared with experiments. The simulations are validated when one obtains a good agreement between the simulated results and the experiments. Then the simulations can be used to visualize and study the microscopic details that are not easily accessible from the experimental observations. For example, one could study the details of the collective phenomena associated with phase transitions, defect dynamics, plastic deformation or damage. The simulations may also be extended to various sample environments, which may not have been achieved in the experiments (with due caution about the applicability of the simulations). The simulations can be only as good as the interatomic potential employed. Usually one identifies the nature of the chemical bonds and uses suitable functional forms with adjustable parameters, which may be derived from the available experimental data or results from ab initio calculations. It is advisable to test the validity of the potential for the conditions, such as the interparticle distances and the coordinations that are likely to be encountered in the simulations. We have developed empirical interatomic potentials in variety of complex ionic and molecular solids [3], tested their lattice dynamical predictions against neutron scattering and other experimental data, and used them in molecular dynamics simulations at high temperature and pressure. At Trombay, we have developed the necessary softwares for lattice dynamics (DISPR) [4] and for molecular dynamics simulations (MOLDY), which may be applicable to systems of arbitrary size and symmetry. In particular, these softwares are also applicable to systems containing molecular ions for which the long-ranged Coulomb interaction is handled via the Ewald technique. Our in-house code MOLDY has been parallelized and installed on the BARC-parallel computer.
147
In this review paper we discuss a few examples of MDS largely from our own work. First we illustrate our simulation studies of the crystal to amorphous transition due to static pressure as well as due to shock propagation in alpha-quartz. Simulations of a propagating shock wave [5] using realistic Coulomb and short-ranged interatomic potentials have perhaps been reported for the first time by us [6]. We then discuss our simulations of fullerene solid where the orientational disorder in the structure has been of special interest. NMR and neutron quasi-elastic scattering experiments indicated dynamical reorientations of individual molecules. However, X-ray and neutron diffuse scattering experiments revealed that the orientations of neighbouring molecules are correlated. Specifically, neutron diffuse quasi-elastic scattering showed the presence of correlations extending to at least second neighbours, and its analysis by molecular dynamics simulations revealed rich information on the orientational correlations. Certain models of the intermolecular potential involving electrostatic interaction between bond charges are able to account for the orientational distribution function and the correlations. 2. Parallelization in classical molecular dynamics simulation As noted above, the simulation requires the calculation of the force on each atom in any given configuration. This calculation involves forces between all pairs of interacting atoms. For a macro cell or cluster of about 10 000 atoms or more, there would be millions of such pairs depending on the range of interatomic interaction. Most of the computation time is used in the force calculation, which is repetitive at each time step of the simulation. Therefore it is most effective to parallelize the force calculation. Essentially we parallelize the calculation over the millions of interactions by simply distributing them over the available processors. This is called force-loop parallelization and this delivers nearly perfect parallelization efficiency up to a very large number of parallel processors. Our in-house code MOLDY has been parallelized and installed on the BARC-parallel computer that involves a cluster of intel-P4 processors. The number of processors is an input to the code. In the above-mentioned scheme of force-loop parallelization, the communication is only required between the master processor and the slaves, and not among the slaves. However, we have implemented the communications over a tree structure; i.e., the master communicates to two slaves that further communicate to another two and so on up to the input number of slaves. The optimum number of slaves depends on the number of atoms in the simulation. Nearly linear gain of the CPU time may be achieved over several tens of slaves. 3. Static pressure and shock-induced amorphization Amorphous or glassy materials are usually formed by rapid quenching of a melt. However, in case of a large
148
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
number of compounds, especially those composed of polyhedral networks such as quartz, application of pressure often leads to a crystal to amorphous transition. Typically, static or shock pressures of about 2–20 GPa are involved. Molecular dynamics simulations are useful in understanding, and even predicting, such transformations [7]. Earlier we carried out MD simulations under static pressures on quartz, which clearly showed the first-order crystal to amorphous transitions [7a]. While in the parent crystalline phases the silicon atoms are known to have a four-fold tetrahedral coordination of oxygen atoms, the MD simulations revealed that in the transformed amorphous phase there is a mixture of four, five and six-fold coordinated atoms at high pressures. The average coordination number further increases with increasing pressure and reduces with decrease of pressure. The average coordination number in quartz glass reduced from about 5.3 to about 4.5 on release of pressure which is close but different from the quartz glass prepared from melt. The MD simulations also showed that the transition pressure decreased with increase of temperature. As in case of usual MD simulations, the changes in the structural variables with pressure including the equation of state are also obtained. These simulated values in the parent crystalline phase revealed that the bond angles and the bond lengths reach certain limiting values before the event of transformation begins. We have also carried out a MD simulation of shock pressure induced amorphization of alpha-quartz [6]. A very long periodic macro cell of 100–200 nm length is used along the shock direction that contains about 15 000 atoms. These large scale MD simulation have been performed on BARC-parallel computer ANUPAM. We have monitored the shock pressure front, temperature and the average coordination of the silicon atoms along the shock-axis as a function of time after the shock initiation. These simulations reveal a very interesting result, namely, that the shock-induced transformations may occur in a short time scale of a few ps, which is, however, two orders greater than the shortest vibrational period of about 30 fs. This result is found consistent with experimental observations [9]. Fig. 1 gives the shock front profiles at a few pressures. The shock front is expected [8f] to have one or more steps depending on whether the pressure is below the Hugoniot elastic limit, below the phase transition pressure, just above it or much above it. Split shock wave profiles have been observed in experiments on a-quartz [8e]. Fig. 1 clearly shows such a behavior. The shock front is found to be rather sharp with a width of a few nm only. The calculated equation of state for the shock along the c-axis is given in Fig. 2, which includes a comparison with shock data [8]. The agreement is generally very good. The simulated Hugoniot elastic limit is about 16 GPa for shock along the c-axis, which is in good agreement with the shock experimental values of about 14 GPa [8]. The transformation is simulated at 21 GPa and 500 K. While earlier experiments [8a,8b] reported the transition at 14.5 GPa and
Fig. 1. The shock front profiles after 5 ps of initiating the shock along the a-quartz c-axis at a few pressures depicting the multiple split shock (see text). The values of the shock stress are indicated.
Fig. 2. The simulated equation of state and temperature for the shock stress along the a-quartz c-axis. The experimental data are shown by squares [8a] and circles [8b].
475 K, more recently [8e] a transition pressure of 23 GPa has been clearly observed. Thus our simulated transition pressure is in good agreement with experiments. The shock temperature remains at about 350 K for the pressures below the transition pressure of about 21 GPa, but increases substantially at high pressures reaching about 2500 K at 60 GPa. Snapshots of atomic positions under the shock at a few pressures are shown in Fig. 3. While the plot for 8 GPa indicates the a-quartz structure, the plot for 52 GPa indicates amorphization. After the phase transformation above 21 GPa the structure is found to be amorphous with several silicon atoms having acquired a coordination of 5 or 6 oxygen atoms with an increase in the Si–O bond lengths and very broad bond-angle distributions. These results have some similarity with the microscopic arrangement obtained under static pressure [5a] although the distributions appear to be broader under shock stress.
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
Fig. 3. A snap shot of a layer of the structure showing some of the SiOx (x = 4–6) polyhedra (c-axis horizontal, a-axis vertical) under the shock stress along the a-quartz c-axis of 8 and 52 GPa (upper and lower, respectively).
This evidence from the MD simulations of the shock experiments is very important since an in situ identification of the shocked sample phase in experiments has not been yet obtained. In simulations, the recovered phase after the shock unloading from the transformed phase is also retained to be amorphous with an average silicon coordination of about 4.5 only and a density of about 25% more than the a-quartz density at ambient pressure. These results are consistent with experiments [8d] and recent NMR data on shock-recovered samples [10]. A signature of the important difference between the shock and static pressure, at the same value of the compression, is revealed from Raman experiments [11]. The Raman line of the A1 (466 cm 1) phonon of alpha-quartz has been measured both with static hydrostatic pressure and under shock loading along the z-axis. For the same density compression, the line shifts under shock loading were found to be considerably greater than that under static loading. This A1 phonon is associated with the symmetric bending of the Si–O–Si bond, which is believed to be the primary mechanism of compression in alpha-quartz. We have simulated the variation in the A1 Raman phonon both under static and shock loading. The simulation results are in fair agreement with the experiments. 4. Fullerene solid The fullerene is a nearly spherical molecule made of sixty carbon atoms. Its surface contains 30 C@C double
149
bonds, 60 C–C single bonds, 12 pentagons and 20 hexagons. In the crystal, the molecular centres arrange in a fcc lattice (see [12] for a review). At T above 260 K, the molecules rotate nearly freely and no librations are observed. Below 260 K, in the simple cubic lattice, the molecules stay in two possible orientations of nearly equal energies, so their relative population varies from 60:40 at 250 K to 85:15 at T below 90 K. In the majority (minority) state, a double-bond centre of one molecule faces opposite a pentagon (hexagon) of a neighbouring molecule. The intermolecular interaction in fullerene is largely of the van der Waals type since the nearest carbon atoms of the neighbouring molecules are separated by about the sum of their van der Waals radii. Using previously known vdW potential, either from the organic molecules or from graphite, it is possible to nearly predict the lattice constant, bulk modulus and the phonon dispersion relation of the rigid molecular translations. However, both the orientational structure and the librational frequencies are not explained by the vdW potential. Various empirical models [12] were developed with considerable success that involved electrostatic charges. We developed a split-bond-charge model [13] in which a small negative charge ( 0.27 e) was placed at two split positions radially inwards and outwards a double-bond centre, and a compensating positive charge (0.27 e) placed on the carbon atom. Among all the available models, this model produced the best agreement with the experimental phonon dispersion relation of the external modes and the associated specific heat at low temperatures [13]. The short-range correlations in a disordered structure may be conveniently studied by the wave-vector dependence of the diffuse scattering of X-rays or neutrons from a single crystal. Above 260 K in the disordered phase, the diffuse scattering experiments [14] indicated significant short-range order of the orientations extending to several neighbours beyond 2 nm. The order was found to be much more pronounced in neutron quasi-elastic scattering than in X-ray scattering [15], which indicated its dynamical nature. We have carried out MD simulation on 256 molecular periodic cells, to understand the observed pattern of the diffuse scattering and thereby reveal the nature of the short-range order [13,14]. Simulation using the splitbond-charge model produces the pattern of the atomic density on the molecular surface in good agreement with that derived from the diffraction experiments [16] (Fig. 4). So also the calculated diffuse scattering from the MD simulation compares very well with the observed diffuse scattering. The simulated orientational correlations indicate preference for the double bonds on one molecule to face opposite the pentagons of the neighbouring molecules (Fig. 5). In contrast, the MD simulation using the vdW model produces an atomic density pattern in clear disagreement with experiment (Fig. 4), and predicts preference for the double bond to hexagon correlations (Fig. 5). In this way, these MD simulations allowed to visualize the com-
150
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
5. Summary and outlook The results presented in this review illustrate the power of MD simulations in revealing the microscopic details of the structures and the dynamics in complex materials. Some more examples of our work on high-temperature superconducting ceramic oxides and minerals of geophysical interest may be found in Ref. [3]. MDS have also been employed in more complex solids such as proteins [17], and in studying large-scale material properties such as crack propagation [18]. As more computing power becomes available, one hopes to take up more challenging simulations of physical phenomena of higher complexity, of masoscopic dimensions and sub-micro second time scales relevant to materials science. References
Fig. 4. Atomic number density on the spherical surface of fullerene, shown as excess over the average density (upper and lower figures). From MD simulations [14] using the split-bond-charge model and the van der Waals potential model, respectively. The brighter and darker areas show excess and deficit densities, respectively.
Fig. 5. The pair correlation functions as calculated from the MD simulations [14] using (a) the split-bond-charge model and (b) the van der Waals potential model. The excess from the average pair correlation function is shown for the pairs of sites on the surface of neighbouring molecules, which reveals the dominant correlation at the nearest average contact separation of about 3 A. DB-double bond, PG-pentagon, and HG-hexagon.
plex details of the microscopic structure and the dynamics, and also helped identify suitable interaction potential.
[1] (a) B.J. Alder, T.E. Wainright, J. Chem. Phys. 31 (1959) 459; A. Rahman, Phys. Rev. A 136 (1964) 105; S.L. Chaplot, Curr. Sci. 55 (1986) 949; S.L. Chaplot, Bull. Mater. Sc. 22 (1999) 279; M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon, Oxford, 1987; J.P. Hansen, M.L. Klein, Phys. Rev. B 13 (1976) 878; S.L. Chaplot, Indian J. Pure Appl. Phys. 32 (1994) 560; (b) K.R. Rao, S.L. Chaplot, in: K.R. Rao (Ed.), Current Trends in Lattice Dynamics, Indian Physics Association, Bombay, 1979, p. 589; (c) M. Parrinello, A. Rahman, Phys. Rev. Lett. 45 (1980) 1196; R. Martonak, A. Lailo, M. Parrinello, Phys. Rev. Lett. 90 (2003) 075505. [2] R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471; R.M. Wentzcovitch, J.L. Martins, G.D. Price, Phys. Rev. Lett. 70 (1993) 3947. [3] S.L. Chaplot, Phys. Rev. B 37 (1988) 7435; S.L. Chaplot, Phase Transitions 19 (1989) 49; S.L. Chaplot, Phys. Rev. B 42 (1990) 2149; S.L. Chaplot, Phys. Rev. B 45 (1992) 4885; S.L. Chaplot, K.R. Rao, Phys. Rev. B 35 (1987) 9771; K.R. Rao et al., Science 236 (1987) 64; K.R. Rao et al., Phys. Chem. Miner. 16 (1988) 83; S.L. Chaplot, N. Choudhury, K.R. Rao, Am. Mineral. 83 (1998) 937; S.L. Chaplot, N. Choudhury, Am. Mineral. 86 (2001) 752–761; R. Mittal et al., Phys. Rev. Lett. 86 (2001) 4692–4695; M.N. Rao et al., Solid State Commun. 121 (2002) 333–338; A. Sen, S.L. Chaplot, R. Mittal, Phys. Rev. B 68 (2003) 134105; P. Goel, N. Choudhury, S.L. Chaplot, Phys. Rev. B 70 (2004) 174307. [4] S.L. Chaplot, Report B.A.R.C.-972, 1978; S.L. Chaplot, unpublished, 1992. [5] (a) A.B. Belonoshko, Science 275 (1997) 955; (b) D.H. Robertson et al., in: High-Pressure Shock Compression of Solids III, Springer-Verlag, New York, 1998, p. 37; (c)B.H. Holian, in: S.C. Schmidt, D.P. Dandekar, J.W. Forbes (Eds.), Shock Compression of Condensed Matter-1997, AIP Conf. Proc. 429, 1998, p. 173; V.V. Zhakhovskii et al., Phys. Rev. Lett. 83 (1999) 1175. [6] S.L. Chaplot, S.K. Sikka, Phys. Rev. B 61 (2000) 11205. [7] (a) S.L. Chaplot, S.K. Sikka, Phys. Rev. B 47 (1993) 5710; (b) M.S. Somayazulu et al., J. Phys.: Condens. Matter 5 (1993) 6345; (c) J.S. Tse, D.D. Klug, Phys. Rev. Lett. 67 (1991) 3559; (d) J. Badro et al., Phys. Rev. Lett. 76 (1996) 772. [8] (a) J. Wackerle, J. Appl. Phys. 33 (1962) 922; (b) G.R. Fowles, J. Geophys. Res. 72 (1967) 5729; (c) J.W. Swegle, J. Appl. Phys. 68 (1990) 1563;
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
[9] [10] [11] [12] [13]
(d) L.C. Chhabildas, in: Y.M. Gupta (Ed.), Shock Waves in Condensed Matter, Plenum, 1986; (e) Y.N. Zhugin, in: W.A. Trzeciakowski (Ed.), High Pressure Science and Technology, World Scientific Pub. Co., Singapore, 1996, p. 953; (f) G.E. Duvall, R.A. Graham, Rev. Mod. Phys. 49 (1977) 523. M.D. Knudson, Y.M. Gupta, A.B. Kunz, Phys. Rev. B 59 (1999) 11704. P.S. Fiske et al., Am. Mineral. 83 (1998) 1285. R.L. Gustavsen, Y.M. Gupta, J. Appl. Phys. 75 (1994) 2837, and ref. therein. J.D. Axe et al., Solid State Phys. 48 (1994) 149. L. Pintschovius, S.L. Chaplot, Z. Phys. B 98 (1995) 527; S.L. Chaplot, L. Pintschovius, Fullerene Sci. Technol. 3 (1995) 707.
151
[14] L. Pintschovius et al., Phys. Rev. Lett. 75 (1995) 2843; L. Pintschovius et al., Physica B 219 & 220 (1996) 148; S.L. Chaplot, L. Pintschovius, Int. J. Mod. Phys. B 13 (1999) 217; S.L. Chaplot, P.S. Schiebel, L. Pintschovius, Fullerene Sci. Technol. 9 (2001) 363–374. [15] R. Moret et al., J. Phys. I 3 (1993) 1085. [16] P.C. Chow et al., Phys. Rev. Lett. 69 (1992) 2943; P. Schiebel et al., Acta Cryst. A 52 (1996) 176; W.I.F. David et al., Proc. R. Soc. (London) A 442 (1993) 129. [17] R. Elber, M. Karplus, Science 235 (1987) 318. [18] F.F. Abraham, H. Gao, Phys. Rev. Lett. 84 (2000) 3113.
Parallelization in classical molecular dynamics simulation and applications S.L. Chaplot Solid State Physics Division, Bhabha Atomic Research Centre, Mumbai 400 085, India
Abstract Over the last four decades, molecular dynamics simulation has been a well-established technique to study a wide variety of physical phenomena. The technique basically involves solving the equations of motion of a system of particles (e.g., atoms or even stars interacting via a potential), and following their trajectories. At a microscopic level, MD simulations have been extensively used in the study of the structures and dynamics, phase transitions and thermodynamic properties of solids and liquids. Simulations of thermal and highpressure processes (often beyond laboratory conditions), shock-stress propagation and defect dynamics have provided unique insights into material properties. This review paper will briefly discuss the parallelization of our in-house code and its uses at the BARC parallel computer. Ó 2006 Elsevier B.V. All rights reserved. Keywords: Molecular dynamics simulation; Thermal properties; High-pressure processes, Shock propagation
1. Introduction Molecular dynamics simulation (MDS) is the name of a technique in which the trajectories of a system of interacting particles are calculated by numerically solving their classical equations of motion [1a]. A very wide range of applications, e.g., the study of phase transitions and dynamics of microscopic defects in solids, has been developed over the past few decades. In fact, this rapid development has kept pace with the growing availability of computing power. An essential requirement for solving the classical equations of motion is the ability to calculate the forces on the particles. There are many ways to do this. The forces may be obtained from an interparticle potential; or from the force constants in a coupled system. Alternatively, the forces in an atomic system could be evaluated from the first principles quantum mechanical calculation of the electronic energies [2]. Accordingly, the MDS may be called
E-mail address: [email protected] 0927-0256/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2005.12.013
the classical MDS or the ab initio MDS. The two are complementary while the former is useful for extensive studies and for relatively more complex systems. Starting with the early simulations by Alder and Wainright [1a] using hard-core potentials and Rahman [1a] using Lennard–Jones potentials, there has been a great deal of progress in terms of various applications in the last four decades. However, the progress was relatively slow in the sixties and seventies as extensive computing facilities were not widely available. In late seventies we studied the defect vibrational modes in rare gas solids and identified them as local modes or the resonance modes [1b]. A major milestone was the development of the constant-pressure technique by Parrinello and Rahman [1c], which facilitated studies of phase transitions involving changes in the crystallographic unit cells. Simulations relevant to materials science involve study of a variety of physical systems, such as the bulk solids and liquids, clusters and nano-particles, surfaces and interfaces, and solid lattices with point and extended defects, etc. The simulations may involve the characteristic parameters of temperature (T), volume (V), total energy (E), external
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
stress (S), total number of particles (N), etc. The quantities like T, S, V, and E, are easily calculated at any instant of time using the atomic positions and velocities. In turn, any of these quantities (T, S, E, V) can be controlled by suitably constraining the atomic trajectories in the simulations [1c]. Simulations commonly employed are those in which various combinations, e.g., (NVT), (NPT), (NVE), or (NE) are kept fixed. It is now convenient to perform simulations with N about 10 000 to more than a million atoms. Bulk properties are generated by assuming periodic boundary conditions. However, the time scale of the microscopic MDS remains limited to about a few nano-seconds, i.e., a few million time steps of a few femto seconds each. This time sale is sufficient for a large variety of dynamical phenomena at the atomic level including critical propagation of extended defects and many of the phase transitions. It is easy to make a contact between the MDS and the various scattering experiments, such as the Raman and neutron scattering experiments from various materials. Such experiments essentially measure the correlations among the positions of atoms at the same time or at different times. These correlations can be directly calculated from the simulated atomic trajectories and compared with experiments. The simulations are validated when one obtains a good agreement between the simulated results and the experiments. Then the simulations can be used to visualize and study the microscopic details that are not easily accessible from the experimental observations. For example, one could study the details of the collective phenomena associated with phase transitions, defect dynamics, plastic deformation or damage. The simulations may also be extended to various sample environments, which may not have been achieved in the experiments (with due caution about the applicability of the simulations). The simulations can be only as good as the interatomic potential employed. Usually one identifies the nature of the chemical bonds and uses suitable functional forms with adjustable parameters, which may be derived from the available experimental data or results from ab initio calculations. It is advisable to test the validity of the potential for the conditions, such as the interparticle distances and the coordinations that are likely to be encountered in the simulations. We have developed empirical interatomic potentials in variety of complex ionic and molecular solids [3], tested their lattice dynamical predictions against neutron scattering and other experimental data, and used them in molecular dynamics simulations at high temperature and pressure. At Trombay, we have developed the necessary softwares for lattice dynamics (DISPR) [4] and for molecular dynamics simulations (MOLDY), which may be applicable to systems of arbitrary size and symmetry. In particular, these softwares are also applicable to systems containing molecular ions for which the long-ranged Coulomb interaction is handled via the Ewald technique. Our in-house code MOLDY has been parallelized and installed on the BARC-parallel computer.
147
In this review paper we discuss a few examples of MDS largely from our own work. First we illustrate our simulation studies of the crystal to amorphous transition due to static pressure as well as due to shock propagation in alpha-quartz. Simulations of a propagating shock wave [5] using realistic Coulomb and short-ranged interatomic potentials have perhaps been reported for the first time by us [6]. We then discuss our simulations of fullerene solid where the orientational disorder in the structure has been of special interest. NMR and neutron quasi-elastic scattering experiments indicated dynamical reorientations of individual molecules. However, X-ray and neutron diffuse scattering experiments revealed that the orientations of neighbouring molecules are correlated. Specifically, neutron diffuse quasi-elastic scattering showed the presence of correlations extending to at least second neighbours, and its analysis by molecular dynamics simulations revealed rich information on the orientational correlations. Certain models of the intermolecular potential involving electrostatic interaction between bond charges are able to account for the orientational distribution function and the correlations. 2. Parallelization in classical molecular dynamics simulation As noted above, the simulation requires the calculation of the force on each atom in any given configuration. This calculation involves forces between all pairs of interacting atoms. For a macro cell or cluster of about 10 000 atoms or more, there would be millions of such pairs depending on the range of interatomic interaction. Most of the computation time is used in the force calculation, which is repetitive at each time step of the simulation. Therefore it is most effective to parallelize the force calculation. Essentially we parallelize the calculation over the millions of interactions by simply distributing them over the available processors. This is called force-loop parallelization and this delivers nearly perfect parallelization efficiency up to a very large number of parallel processors. Our in-house code MOLDY has been parallelized and installed on the BARC-parallel computer that involves a cluster of intel-P4 processors. The number of processors is an input to the code. In the above-mentioned scheme of force-loop parallelization, the communication is only required between the master processor and the slaves, and not among the slaves. However, we have implemented the communications over a tree structure; i.e., the master communicates to two slaves that further communicate to another two and so on up to the input number of slaves. The optimum number of slaves depends on the number of atoms in the simulation. Nearly linear gain of the CPU time may be achieved over several tens of slaves. 3. Static pressure and shock-induced amorphization Amorphous or glassy materials are usually formed by rapid quenching of a melt. However, in case of a large
148
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
number of compounds, especially those composed of polyhedral networks such as quartz, application of pressure often leads to a crystal to amorphous transition. Typically, static or shock pressures of about 2–20 GPa are involved. Molecular dynamics simulations are useful in understanding, and even predicting, such transformations [7]. Earlier we carried out MD simulations under static pressures on quartz, which clearly showed the first-order crystal to amorphous transitions [7a]. While in the parent crystalline phases the silicon atoms are known to have a four-fold tetrahedral coordination of oxygen atoms, the MD simulations revealed that in the transformed amorphous phase there is a mixture of four, five and six-fold coordinated atoms at high pressures. The average coordination number further increases with increasing pressure and reduces with decrease of pressure. The average coordination number in quartz glass reduced from about 5.3 to about 4.5 on release of pressure which is close but different from the quartz glass prepared from melt. The MD simulations also showed that the transition pressure decreased with increase of temperature. As in case of usual MD simulations, the changes in the structural variables with pressure including the equation of state are also obtained. These simulated values in the parent crystalline phase revealed that the bond angles and the bond lengths reach certain limiting values before the event of transformation begins. We have also carried out a MD simulation of shock pressure induced amorphization of alpha-quartz [6]. A very long periodic macro cell of 100–200 nm length is used along the shock direction that contains about 15 000 atoms. These large scale MD simulation have been performed on BARC-parallel computer ANUPAM. We have monitored the shock pressure front, temperature and the average coordination of the silicon atoms along the shock-axis as a function of time after the shock initiation. These simulations reveal a very interesting result, namely, that the shock-induced transformations may occur in a short time scale of a few ps, which is, however, two orders greater than the shortest vibrational period of about 30 fs. This result is found consistent with experimental observations [9]. Fig. 1 gives the shock front profiles at a few pressures. The shock front is expected [8f] to have one or more steps depending on whether the pressure is below the Hugoniot elastic limit, below the phase transition pressure, just above it or much above it. Split shock wave profiles have been observed in experiments on a-quartz [8e]. Fig. 1 clearly shows such a behavior. The shock front is found to be rather sharp with a width of a few nm only. The calculated equation of state for the shock along the c-axis is given in Fig. 2, which includes a comparison with shock data [8]. The agreement is generally very good. The simulated Hugoniot elastic limit is about 16 GPa for shock along the c-axis, which is in good agreement with the shock experimental values of about 14 GPa [8]. The transformation is simulated at 21 GPa and 500 K. While earlier experiments [8a,8b] reported the transition at 14.5 GPa and
Fig. 1. The shock front profiles after 5 ps of initiating the shock along the a-quartz c-axis at a few pressures depicting the multiple split shock (see text). The values of the shock stress are indicated.
Fig. 2. The simulated equation of state and temperature for the shock stress along the a-quartz c-axis. The experimental data are shown by squares [8a] and circles [8b].
475 K, more recently [8e] a transition pressure of 23 GPa has been clearly observed. Thus our simulated transition pressure is in good agreement with experiments. The shock temperature remains at about 350 K for the pressures below the transition pressure of about 21 GPa, but increases substantially at high pressures reaching about 2500 K at 60 GPa. Snapshots of atomic positions under the shock at a few pressures are shown in Fig. 3. While the plot for 8 GPa indicates the a-quartz structure, the plot for 52 GPa indicates amorphization. After the phase transformation above 21 GPa the structure is found to be amorphous with several silicon atoms having acquired a coordination of 5 or 6 oxygen atoms with an increase in the Si–O bond lengths and very broad bond-angle distributions. These results have some similarity with the microscopic arrangement obtained under static pressure [5a] although the distributions appear to be broader under shock stress.
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
Fig. 3. A snap shot of a layer of the structure showing some of the SiOx (x = 4–6) polyhedra (c-axis horizontal, a-axis vertical) under the shock stress along the a-quartz c-axis of 8 and 52 GPa (upper and lower, respectively).
This evidence from the MD simulations of the shock experiments is very important since an in situ identification of the shocked sample phase in experiments has not been yet obtained. In simulations, the recovered phase after the shock unloading from the transformed phase is also retained to be amorphous with an average silicon coordination of about 4.5 only and a density of about 25% more than the a-quartz density at ambient pressure. These results are consistent with experiments [8d] and recent NMR data on shock-recovered samples [10]. A signature of the important difference between the shock and static pressure, at the same value of the compression, is revealed from Raman experiments [11]. The Raman line of the A1 (466 cm 1) phonon of alpha-quartz has been measured both with static hydrostatic pressure and under shock loading along the z-axis. For the same density compression, the line shifts under shock loading were found to be considerably greater than that under static loading. This A1 phonon is associated with the symmetric bending of the Si–O–Si bond, which is believed to be the primary mechanism of compression in alpha-quartz. We have simulated the variation in the A1 Raman phonon both under static and shock loading. The simulation results are in fair agreement with the experiments. 4. Fullerene solid The fullerene is a nearly spherical molecule made of sixty carbon atoms. Its surface contains 30 C@C double
149
bonds, 60 C–C single bonds, 12 pentagons and 20 hexagons. In the crystal, the molecular centres arrange in a fcc lattice (see [12] for a review). At T above 260 K, the molecules rotate nearly freely and no librations are observed. Below 260 K, in the simple cubic lattice, the molecules stay in two possible orientations of nearly equal energies, so their relative population varies from 60:40 at 250 K to 85:15 at T below 90 K. In the majority (minority) state, a double-bond centre of one molecule faces opposite a pentagon (hexagon) of a neighbouring molecule. The intermolecular interaction in fullerene is largely of the van der Waals type since the nearest carbon atoms of the neighbouring molecules are separated by about the sum of their van der Waals radii. Using previously known vdW potential, either from the organic molecules or from graphite, it is possible to nearly predict the lattice constant, bulk modulus and the phonon dispersion relation of the rigid molecular translations. However, both the orientational structure and the librational frequencies are not explained by the vdW potential. Various empirical models [12] were developed with considerable success that involved electrostatic charges. We developed a split-bond-charge model [13] in which a small negative charge ( 0.27 e) was placed at two split positions radially inwards and outwards a double-bond centre, and a compensating positive charge (0.27 e) placed on the carbon atom. Among all the available models, this model produced the best agreement with the experimental phonon dispersion relation of the external modes and the associated specific heat at low temperatures [13]. The short-range correlations in a disordered structure may be conveniently studied by the wave-vector dependence of the diffuse scattering of X-rays or neutrons from a single crystal. Above 260 K in the disordered phase, the diffuse scattering experiments [14] indicated significant short-range order of the orientations extending to several neighbours beyond 2 nm. The order was found to be much more pronounced in neutron quasi-elastic scattering than in X-ray scattering [15], which indicated its dynamical nature. We have carried out MD simulation on 256 molecular periodic cells, to understand the observed pattern of the diffuse scattering and thereby reveal the nature of the short-range order [13,14]. Simulation using the splitbond-charge model produces the pattern of the atomic density on the molecular surface in good agreement with that derived from the diffraction experiments [16] (Fig. 4). So also the calculated diffuse scattering from the MD simulation compares very well with the observed diffuse scattering. The simulated orientational correlations indicate preference for the double bonds on one molecule to face opposite the pentagons of the neighbouring molecules (Fig. 5). In contrast, the MD simulation using the vdW model produces an atomic density pattern in clear disagreement with experiment (Fig. 4), and predicts preference for the double bond to hexagon correlations (Fig. 5). In this way, these MD simulations allowed to visualize the com-
150
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
5. Summary and outlook The results presented in this review illustrate the power of MD simulations in revealing the microscopic details of the structures and the dynamics in complex materials. Some more examples of our work on high-temperature superconducting ceramic oxides and minerals of geophysical interest may be found in Ref. [3]. MDS have also been employed in more complex solids such as proteins [17], and in studying large-scale material properties such as crack propagation [18]. As more computing power becomes available, one hopes to take up more challenging simulations of physical phenomena of higher complexity, of masoscopic dimensions and sub-micro second time scales relevant to materials science. References
Fig. 4. Atomic number density on the spherical surface of fullerene, shown as excess over the average density (upper and lower figures). From MD simulations [14] using the split-bond-charge model and the van der Waals potential model, respectively. The brighter and darker areas show excess and deficit densities, respectively.
Fig. 5. The pair correlation functions as calculated from the MD simulations [14] using (a) the split-bond-charge model and (b) the van der Waals potential model. The excess from the average pair correlation function is shown for the pairs of sites on the surface of neighbouring molecules, which reveals the dominant correlation at the nearest average contact separation of about 3 A. DB-double bond, PG-pentagon, and HG-hexagon.
plex details of the microscopic structure and the dynamics, and also helped identify suitable interaction potential.
[1] (a) B.J. Alder, T.E. Wainright, J. Chem. Phys. 31 (1959) 459; A. Rahman, Phys. Rev. A 136 (1964) 105; S.L. Chaplot, Curr. Sci. 55 (1986) 949; S.L. Chaplot, Bull. Mater. Sc. 22 (1999) 279; M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon, Oxford, 1987; J.P. Hansen, M.L. Klein, Phys. Rev. B 13 (1976) 878; S.L. Chaplot, Indian J. Pure Appl. Phys. 32 (1994) 560; (b) K.R. Rao, S.L. Chaplot, in: K.R. Rao (Ed.), Current Trends in Lattice Dynamics, Indian Physics Association, Bombay, 1979, p. 589; (c) M. Parrinello, A. Rahman, Phys. Rev. Lett. 45 (1980) 1196; R. Martonak, A. Lailo, M. Parrinello, Phys. Rev. Lett. 90 (2003) 075505. [2] R. Car, M. Parrinello, Phys. Rev. Lett. 55 (1985) 2471; R.M. Wentzcovitch, J.L. Martins, G.D. Price, Phys. Rev. Lett. 70 (1993) 3947. [3] S.L. Chaplot, Phys. Rev. B 37 (1988) 7435; S.L. Chaplot, Phase Transitions 19 (1989) 49; S.L. Chaplot, Phys. Rev. B 42 (1990) 2149; S.L. Chaplot, Phys. Rev. B 45 (1992) 4885; S.L. Chaplot, K.R. Rao, Phys. Rev. B 35 (1987) 9771; K.R. Rao et al., Science 236 (1987) 64; K.R. Rao et al., Phys. Chem. Miner. 16 (1988) 83; S.L. Chaplot, N. Choudhury, K.R. Rao, Am. Mineral. 83 (1998) 937; S.L. Chaplot, N. Choudhury, Am. Mineral. 86 (2001) 752–761; R. Mittal et al., Phys. Rev. Lett. 86 (2001) 4692–4695; M.N. Rao et al., Solid State Commun. 121 (2002) 333–338; A. Sen, S.L. Chaplot, R. Mittal, Phys. Rev. B 68 (2003) 134105; P. Goel, N. Choudhury, S.L. Chaplot, Phys. Rev. B 70 (2004) 174307. [4] S.L. Chaplot, Report B.A.R.C.-972, 1978; S.L. Chaplot, unpublished, 1992. [5] (a) A.B. Belonoshko, Science 275 (1997) 955; (b) D.H. Robertson et al., in: High-Pressure Shock Compression of Solids III, Springer-Verlag, New York, 1998, p. 37; (c)B.H. Holian, in: S.C. Schmidt, D.P. Dandekar, J.W. Forbes (Eds.), Shock Compression of Condensed Matter-1997, AIP Conf. Proc. 429, 1998, p. 173; V.V. Zhakhovskii et al., Phys. Rev. Lett. 83 (1999) 1175. [6] S.L. Chaplot, S.K. Sikka, Phys. Rev. B 61 (2000) 11205. [7] (a) S.L. Chaplot, S.K. Sikka, Phys. Rev. B 47 (1993) 5710; (b) M.S. Somayazulu et al., J. Phys.: Condens. Matter 5 (1993) 6345; (c) J.S. Tse, D.D. Klug, Phys. Rev. Lett. 67 (1991) 3559; (d) J. Badro et al., Phys. Rev. Lett. 76 (1996) 772. [8] (a) J. Wackerle, J. Appl. Phys. 33 (1962) 922; (b) G.R. Fowles, J. Geophys. Res. 72 (1967) 5729; (c) J.W. Swegle, J. Appl. Phys. 68 (1990) 1563;
S.L. Chaplot / Computational Materials Science 37 (2006) 146–151
[9] [10] [11] [12] [13]
(d) L.C. Chhabildas, in: Y.M. Gupta (Ed.), Shock Waves in Condensed Matter, Plenum, 1986; (e) Y.N. Zhugin, in: W.A. Trzeciakowski (Ed.), High Pressure Science and Technology, World Scientific Pub. Co., Singapore, 1996, p. 953; (f) G.E. Duvall, R.A. Graham, Rev. Mod. Phys. 49 (1977) 523. M.D. Knudson, Y.M. Gupta, A.B. Kunz, Phys. Rev. B 59 (1999) 11704. P.S. Fiske et al., Am. Mineral. 83 (1998) 1285. R.L. Gustavsen, Y.M. Gupta, J. Appl. Phys. 75 (1994) 2837, and ref. therein. J.D. Axe et al., Solid State Phys. 48 (1994) 149. L. Pintschovius, S.L. Chaplot, Z. Phys. B 98 (1995) 527; S.L. Chaplot, L. Pintschovius, Fullerene Sci. Technol. 3 (1995) 707.
151
[14] L. Pintschovius et al., Phys. Rev. Lett. 75 (1995) 2843; L. Pintschovius et al., Physica B 219 & 220 (1996) 148; S.L. Chaplot, L. Pintschovius, Int. J. Mod. Phys. B 13 (1999) 217; S.L. Chaplot, P.S. Schiebel, L. Pintschovius, Fullerene Sci. Technol. 9 (2001) 363–374. [15] R. Moret et al., J. Phys. I 3 (1993) 1085. [16] P.C. Chow et al., Phys. Rev. Lett. 69 (1992) 2943; P. Schiebel et al., Acta Cryst. A 52 (1996) 176; W.I.F. David et al., Proc. R. Soc. (London) A 442 (1993) 129. [17] R. Elber, M. Karplus, Science 235 (1987) 318. [18] F.F. Abraham, H. Gao, Phys. Rev. Lett. 84 (2000) 3113.