- Email: [email protected]
Molecular dynamics simulations of cluster collisions
COMPUTATIONAL MATERIALS SCIENCE ELSEVIER
Computational Materials Science 2 (1994) 585-588
Molecular dynamics simulations of cluster collisions G. Seifert Institut fi~r Theoretische Physik, Technische Universith't Dresden, D-01069 Dresden, Germany
J. Schulte National Superconductiong Cyclotron Laboratory, Michigan State University, East Lansing, Michigan 488241 USA
(Received 25 October 1993; accepted 20 December 1993)
Abstract
Molecular dynamics simulations are very powerful to get more insight into the dynamics of collision processes. Simulations of such processes have been successfully performed using empirical interaction potentials. More recently the combination of the Density-Functional-Theory within the Local Density Approximation (DFT-LDA) with Molecular Dynamics (MD) has been realized. A simplified L C A O - D F T - L D A scheme, which enables to consider the electronic states in the calculation of the forces on the atoms, was developed and applied for simulations of molecule-cluster and cluster-cluster collisions. In the present paper the authors demonstrate the extension of such calculations under the consideration of the statistics in cluster collision processes by using parallel computer techniques.
1. Introduction
Collisions between atomic clusters and atomic clusters with atoms or molecules are becoming a very attractive new field in cluster research [1]. The study of such collisions may be helpful for l e a r n i n g m o r e about the general dynamics of collisions of finite many particle systems. Obviously, similarities and also differences are expected in comparison to a t o m - a t o m , a t o m - m o l e cule, m o l e c u l e - m o l e c u l e as well as nuclear heavy ion collisions. Furthermore, collision processes between clusters may play also a crucial role in particle growth processes [2]. As to the theoretical examination of the collision dynamics of a t o m - m o l e c u l e reactions quasiclassical trajectory calculations are well estab-
lished since the pioneering investigations of Karplus and others in the sixtees (see e.g. [3]). However, even such classical trajectory calculations are up to now restricted to systems with a few atoms due to the complexity of the potential hypersurfaces. Traditional quantum chemical ab initio methods allow calculations for rather small systems and are up to now practically not applicable for studying collisions dynamics of more complex systems, due to their high computational effort. On the other hand empirical potential models or tight-binding methods give useful potentials for simulations of the dynamics near the equilibrium structures of molecules, clusters and solids. However, the application of empirical potentials or tight binding methods to processes far from equilibrium structures, e.g. in high energy
0927-0256/94/$07.00 © 1994 - Elsevier Science B.V. All rights reserved SSDI 0927-0256(94)00029-C
586
G, Seifert, J. Schulte / Computational Materials Science 2 (1994) 585-588
collisions, is rather questionable. Furthermore, empirical potential models cannot consider the interplay between electronic structure and the structural atomic configuration. Density-Functional Theory (DFT) within the Local Density Approximation (LDA) allows the calculation of the potential hypersurface of even large systems with a reasonable computational effort and satisfactory accuracy. In this way nonempirical (ab initio) MD simulations became realizable [4]. In order to be able to simulate collision processes of clusters (up to approximately 100 atoms) for simulation periods of several picoseconds, we developed a method based on LDA. The electronic wave functions (KohnSham orbitals) are written as linear combinations of atomic orbitals (LCAO). The effective one particle potential in the Kohn-Sham Hamiltonian is constructed from potentials of neutral atoms. Such treatment is physically comparable with the so called Harris-Functional [5]. Furthermore, the short range part of the interatomic potential energy in E [ R I . . . R N] is described by an empirical potential. Such treatment may be viewed as a "hybrid" between an ab initio MD based on DFT and the usage of purely empirical potentials. It has the advantage over the latter in overcoming the "transferability" problem, where parametrized potentials based on experimental data can be difficult to transfer from one system to another. On the other hand, it requires less computational effort than the ab initio molecular dynamics. Furthermore, besides neutral systems, collisions of charged clusters may be handled also. The method described here was already successfully applied for describing structural properties of phosphorus [6] clusters, amorphous carbon systems and hydrocarbons [12] using MD and the technique of simulated annealing. Also collisions between clusters have been simulated with this method (Na n + Na n, n = 8,9 [11] C60 + C60 [13]). The relatively small computational effort allows to consider statistics in the collision process, and in this way to calculate reaction cross sections and reaction rates of collisions of clusters by integration over impact parameters, averaging over different relative orientations of the collision
partners, consideration of internal energies and calculations at different collision energies. To realize the consideration of statistics in collision processes (calculation of several thousand MD trajectories) the program was implemented on a parallel computer (nCUBE2).
2. Method
The electronic wave functions (Kohn-Sham orbitals ~O(r)) are written as linear combinations of atomic orbitals 4~(r - R j), centered at the nuclei j (LCAO-ansatz):
4t(r) = ~_, C~c~,(r - Rj)
(1)
/x
The valence wave functions are considered in the LCAO ansatz. Each wave function is represented by a set of 12 Slater type functions. The effective one particle potential V~ff in the Kohn-Sham Hamiltonian /~ = ~"+ Vaf(r) is approximated as a sum of potentials of neutral atoms. Consistent with this approximation one has to neglect several contributions to the Hamiltonian matrix elements h ~ (see [7] and [8]). Furthermore, we use an empirical ansatz for the repulsive force between the atoms at small distances. Newton's equations of motion:
M k R = -OE/OR k
(2)
are integrated numerically. For a detailed description and applications of the method see Refs. [6,9]. The statistics can be considered by calculations of trajectories with different initial conditions. The calculations of trajectories with different impact parameters b and different relative orientations of the colliding particles. The reaction rates may then be calculated by integrating the reaction probability P(b) (number of the "reactive" trajectories divided by the total numbers of trajectories at a given b over the impact parameter b [10]:
k = ~
mp(b)2rrb db.
(3).
G. Seifert, J. Schulte / Computational Materials Science 2 (1994) 585-588 T b is the t e m p e r a t u r e corresponding to the relative c.o.m, collision energy, and /z is the reduced mass of the colliding system. The reaction cross section is then given by:
o"r = k / V ,
(4),
where ~ is the constant prefactor in Eq. (12). The integral given in Eq. (12) can be computed using standard Monte-Carlo integration techniques. This means, the initial impact p a r a m e t e r s are chosen from a normalized b 2 distribution,
f(b)
db=~l~b2m] k 0,
db'
if0
(5).
otherwise
where b m is a maximum impact parameter. The maximum impact p a r a m e t e r b m is determined by calculating the reaction onset impact p a r a m e t e r of the respective extreme collision orientations. To consider the orientational statistics, or in other words, to avoid collision orientation correlation, the colliding particles are randomly rotated through Euler's angles about their internal axes. Internal t e m p e r a t u r e or thermal collisions may be simulated by considering a Maxwell-Boltzmann velocity distribution (see e.g. [101). The realization of the collision statistics is ideally suited for parallel computing. The different trajectories or sets of trajectories can run parallel on different processors of a parallel computer. The program has been implemented on a multiple instruction multiple data ( M I M D ) parallel super computer (nCUBE2).
3. Applications We have presented a new scheme for handling molecular collision processes on the basis of an approximate D F T - L D A method combined with molecular dynamics (MD) and considering statistics. The method was applied to the reaction of carbon clusters with D 2. Especially the collisions of C 2 and C~- with D 2 w e r e simulated. The simulations could give new insights into the mechanism of this reaction. For neutral C~
587
we found as the preferred reaction channel the formation of an acetylenelike D - C 2 - D complex. Whereas for C~- the formation of C 2 - D dominates strongly over that of D - C 2 - D . This result shows clearly the importance of the dynamics in the understanding of such reactions, since in both cases (C 2 and C~) the acetylenelike ( D - C 2 - D ) complex is the energetically favoured product. The calculated reaction rate constant for C~-+ D 2 is in agreement with experimental finding [15]. For a more detailed discussion of these results see Ref. [14]. The advantage of the method is that it can be easily applied to much larger systems, and there are clear indications for new qualitative aspects in the dynamics of collisions for liarge systems with many internal degrees of freedom [1,2,11]. That is, the principle o f microscopic reversibility no longer holds [2]. Investigations on larger systems (C60 + C60,C60 + D 2) employing the presented molecular dynamics statistical ensemble m e t h o d are in progress, and will be published elsewhere.
Acknowledgement This work has been partly supported by N A T O under Grant No. C R G 930351, and by a computational grant by nCUBE.
References [1] R. Schmidt, G. Seifert and H.-O. Lutz, in: Proc. 88th WE-Heraeus-Seminar, Lecture Notes in Physics, eds. R. Schmidt, H.O. Lutz und R. Dreizler (Springer, Berlin, 1992) p. 128. [2] E. Blaisten-Barojas and R.M. Zachariah, Phys. Rev. B 45 (1992) 4403. [3] M. Karplus, R.N. Porter and R.D. Sharma, J. Chem. Phys. 48 (1965) 3259. [4] R. Car and M. Parrine|lo, Phys. Rev. Lett 55 (1985) 2471. [5] J. Harris: Phys. Rev. B 31 (1986) 1770. [6] G. Seifert, and R.O. Jones, Z. Phys. D 20 (1991) 77; J. Chem. Phys. 96 (1992) 7564. [7] G. Seifert, H. Eschrig and W. Bieger, Z. Phys. Chem. (Leipzig) 267 (1986) 529. [8] G. Seifert and H. Eschrig: Phys. Stat. Sol. B 127 (1985) 573.
588
G. Seifert, J. Schulte / Computational Materials Science 2 (1994) 585-588
[9] G. Seifert and R. Schmidt: New J. Chem. 16 (1992) 1145. [10] J. Schulte, R.R. Lucchese and W.H. Marlow: J. Chem. Phys. 99 (2) (1993) 1178. [11] R. Schmidt, G. Seifert and H.O. Lutz: Physics Letters A 158 (1991) 231. [12] P. Blaudeck, Th. Frauenheim, G. Seifert and E. Fromm: J. Phys. Condensed Matter 4 (1992) 6389.
[13] G. Seifert and R. Schmidt: Journal of Modern Physics B 96 (1992) 3845. [14] J. Schulte and G. Seifert, Chem. Phys. Lett. 221 (1994) 230. [15] S.W. McElvany, B.I. Dunlap and A. O'Keefe: J. Chem. Phys. 86 (1987) 715.
Computational Materials Science 2 (1994) 585-588
Molecular dynamics simulations of cluster collisions G. Seifert Institut fi~r Theoretische Physik, Technische Universith't Dresden, D-01069 Dresden, Germany
J. Schulte National Superconductiong Cyclotron Laboratory, Michigan State University, East Lansing, Michigan 488241 USA
(Received 25 October 1993; accepted 20 December 1993)
Abstract
Molecular dynamics simulations are very powerful to get more insight into the dynamics of collision processes. Simulations of such processes have been successfully performed using empirical interaction potentials. More recently the combination of the Density-Functional-Theory within the Local Density Approximation (DFT-LDA) with Molecular Dynamics (MD) has been realized. A simplified L C A O - D F T - L D A scheme, which enables to consider the electronic states in the calculation of the forces on the atoms, was developed and applied for simulations of molecule-cluster and cluster-cluster collisions. In the present paper the authors demonstrate the extension of such calculations under the consideration of the statistics in cluster collision processes by using parallel computer techniques.
1. Introduction
Collisions between atomic clusters and atomic clusters with atoms or molecules are becoming a very attractive new field in cluster research [1]. The study of such collisions may be helpful for l e a r n i n g m o r e about the general dynamics of collisions of finite many particle systems. Obviously, similarities and also differences are expected in comparison to a t o m - a t o m , a t o m - m o l e cule, m o l e c u l e - m o l e c u l e as well as nuclear heavy ion collisions. Furthermore, collision processes between clusters may play also a crucial role in particle growth processes [2]. As to the theoretical examination of the collision dynamics of a t o m - m o l e c u l e reactions quasiclassical trajectory calculations are well estab-
lished since the pioneering investigations of Karplus and others in the sixtees (see e.g. [3]). However, even such classical trajectory calculations are up to now restricted to systems with a few atoms due to the complexity of the potential hypersurfaces. Traditional quantum chemical ab initio methods allow calculations for rather small systems and are up to now practically not applicable for studying collisions dynamics of more complex systems, due to their high computational effort. On the other hand empirical potential models or tight-binding methods give useful potentials for simulations of the dynamics near the equilibrium structures of molecules, clusters and solids. However, the application of empirical potentials or tight binding methods to processes far from equilibrium structures, e.g. in high energy
0927-0256/94/$07.00 © 1994 - Elsevier Science B.V. All rights reserved SSDI 0927-0256(94)00029-C
586
G, Seifert, J. Schulte / Computational Materials Science 2 (1994) 585-588
collisions, is rather questionable. Furthermore, empirical potential models cannot consider the interplay between electronic structure and the structural atomic configuration. Density-Functional Theory (DFT) within the Local Density Approximation (LDA) allows the calculation of the potential hypersurface of even large systems with a reasonable computational effort and satisfactory accuracy. In this way nonempirical (ab initio) MD simulations became realizable [4]. In order to be able to simulate collision processes of clusters (up to approximately 100 atoms) for simulation periods of several picoseconds, we developed a method based on LDA. The electronic wave functions (KohnSham orbitals) are written as linear combinations of atomic orbitals (LCAO). The effective one particle potential in the Kohn-Sham Hamiltonian is constructed from potentials of neutral atoms. Such treatment is physically comparable with the so called Harris-Functional [5]. Furthermore, the short range part of the interatomic potential energy in E [ R I . . . R N] is described by an empirical potential. Such treatment may be viewed as a "hybrid" between an ab initio MD based on DFT and the usage of purely empirical potentials. It has the advantage over the latter in overcoming the "transferability" problem, where parametrized potentials based on experimental data can be difficult to transfer from one system to another. On the other hand, it requires less computational effort than the ab initio molecular dynamics. Furthermore, besides neutral systems, collisions of charged clusters may be handled also. The method described here was already successfully applied for describing structural properties of phosphorus [6] clusters, amorphous carbon systems and hydrocarbons [12] using MD and the technique of simulated annealing. Also collisions between clusters have been simulated with this method (Na n + Na n, n = 8,9 [11] C60 + C60 [13]). The relatively small computational effort allows to consider statistics in the collision process, and in this way to calculate reaction cross sections and reaction rates of collisions of clusters by integration over impact parameters, averaging over different relative orientations of the collision
partners, consideration of internal energies and calculations at different collision energies. To realize the consideration of statistics in collision processes (calculation of several thousand MD trajectories) the program was implemented on a parallel computer (nCUBE2).
2. Method
The electronic wave functions (Kohn-Sham orbitals ~O(r)) are written as linear combinations of atomic orbitals 4~(r - R j), centered at the nuclei j (LCAO-ansatz):
4t(r) = ~_, C~c~,(r - Rj)
(1)
/x
The valence wave functions are considered in the LCAO ansatz. Each wave function is represented by a set of 12 Slater type functions. The effective one particle potential V~ff in the Kohn-Sham Hamiltonian /~ = ~"+ Vaf(r) is approximated as a sum of potentials of neutral atoms. Consistent with this approximation one has to neglect several contributions to the Hamiltonian matrix elements h ~ (see [7] and [8]). Furthermore, we use an empirical ansatz for the repulsive force between the atoms at small distances. Newton's equations of motion:
M k R = -OE/OR k
(2)
are integrated numerically. For a detailed description and applications of the method see Refs. [6,9]. The statistics can be considered by calculations of trajectories with different initial conditions. The calculations of trajectories with different impact parameters b and different relative orientations of the colliding particles. The reaction rates may then be calculated by integrating the reaction probability P(b) (number of the "reactive" trajectories divided by the total numbers of trajectories at a given b over the impact parameter b [10]:
k = ~
mp(b)2rrb db.
(3).
G. Seifert, J. Schulte / Computational Materials Science 2 (1994) 585-588 T b is the t e m p e r a t u r e corresponding to the relative c.o.m, collision energy, and /z is the reduced mass of the colliding system. The reaction cross section is then given by:
o"r = k / V ,
(4),
where ~ is the constant prefactor in Eq. (12). The integral given in Eq. (12) can be computed using standard Monte-Carlo integration techniques. This means, the initial impact p a r a m e t e r s are chosen from a normalized b 2 distribution,
f(b)
db=~l~b2m] k 0,
db'
if0
(5).
otherwise
where b m is a maximum impact parameter. The maximum impact p a r a m e t e r b m is determined by calculating the reaction onset impact p a r a m e t e r of the respective extreme collision orientations. To consider the orientational statistics, or in other words, to avoid collision orientation correlation, the colliding particles are randomly rotated through Euler's angles about their internal axes. Internal t e m p e r a t u r e or thermal collisions may be simulated by considering a Maxwell-Boltzmann velocity distribution (see e.g. [101). The realization of the collision statistics is ideally suited for parallel computing. The different trajectories or sets of trajectories can run parallel on different processors of a parallel computer. The program has been implemented on a multiple instruction multiple data ( M I M D ) parallel super computer (nCUBE2).
3. Applications We have presented a new scheme for handling molecular collision processes on the basis of an approximate D F T - L D A method combined with molecular dynamics (MD) and considering statistics. The method was applied to the reaction of carbon clusters with D 2. Especially the collisions of C 2 and C~- with D 2 w e r e simulated. The simulations could give new insights into the mechanism of this reaction. For neutral C~
587
we found as the preferred reaction channel the formation of an acetylenelike D - C 2 - D complex. Whereas for C~- the formation of C 2 - D dominates strongly over that of D - C 2 - D . This result shows clearly the importance of the dynamics in the understanding of such reactions, since in both cases (C 2 and C~) the acetylenelike ( D - C 2 - D ) complex is the energetically favoured product. The calculated reaction rate constant for C~-+ D 2 is in agreement with experimental finding [15]. For a more detailed discussion of these results see Ref. [14]. The advantage of the method is that it can be easily applied to much larger systems, and there are clear indications for new qualitative aspects in the dynamics of collisions for liarge systems with many internal degrees of freedom [1,2,11]. That is, the principle o f microscopic reversibility no longer holds [2]. Investigations on larger systems (C60 + C60,C60 + D 2) employing the presented molecular dynamics statistical ensemble m e t h o d are in progress, and will be published elsewhere.
Acknowledgement This work has been partly supported by N A T O under Grant No. C R G 930351, and by a computational grant by nCUBE.
References [1] R. Schmidt, G. Seifert and H.-O. Lutz, in: Proc. 88th WE-Heraeus-Seminar, Lecture Notes in Physics, eds. R. Schmidt, H.O. Lutz und R. Dreizler (Springer, Berlin, 1992) p. 128. [2] E. Blaisten-Barojas and R.M. Zachariah, Phys. Rev. B 45 (1992) 4403. [3] M. Karplus, R.N. Porter and R.D. Sharma, J. Chem. Phys. 48 (1965) 3259. [4] R. Car and M. Parrine|lo, Phys. Rev. Lett 55 (1985) 2471. [5] J. Harris: Phys. Rev. B 31 (1986) 1770. [6] G. Seifert, and R.O. Jones, Z. Phys. D 20 (1991) 77; J. Chem. Phys. 96 (1992) 7564. [7] G. Seifert, H. Eschrig and W. Bieger, Z. Phys. Chem. (Leipzig) 267 (1986) 529. [8] G. Seifert and H. Eschrig: Phys. Stat. Sol. B 127 (1985) 573.
588
G. Seifert, J. Schulte / Computational Materials Science 2 (1994) 585-588
[9] G. Seifert and R. Schmidt: New J. Chem. 16 (1992) 1145. [10] J. Schulte, R.R. Lucchese and W.H. Marlow: J. Chem. Phys. 99 (2) (1993) 1178. [11] R. Schmidt, G. Seifert and H.O. Lutz: Physics Letters A 158 (1991) 231. [12] P. Blaudeck, Th. Frauenheim, G. Seifert and E. Fromm: J. Phys. Condensed Matter 4 (1992) 6389.
[13] G. Seifert and R. Schmidt: Journal of Modern Physics B 96 (1992) 3845. [14] J. Schulte and G. Seifert, Chem. Phys. Lett. 221 (1994) 230. [15] S.W. McElvany, B.I. Dunlap and A. O'Keefe: J. Chem. Phys. 86 (1987) 715.