J. CHEM. SOC. FARADAY TRANS., 1994, 90(14), 2009-2013 Electronic States and the Metal-Insulator Transition in Caesium-Ammonia Solutions Zhihong Deng and Michael L. Klein* Department of Chemistry, University of Pennsylvania, Philadelphia, PA 19104-6323, USA Glenn J. Martyna Department of Chemistry, Indiana University, Bloomington, IN 474054001, USA The nature of the electronic states in caesium-ammonia solutions is examined from the insulating to the metallic regime using two different microscopic models (Z. Deng, G. J. Martyna and M. L. Klein, J. Chem. Phys., 1994, 100, 7590.). In the first model, the ammonia molecules are treated via a classical point-charge model and the cations as a positive neutralizing background. In the second model, the ammonia solvent is made fully polariz-able and the cations, here caesium, are explicitly included.The solvent and ions are treated classically and the electronic degrees of freedom are handled using the Car-Parrinello method and density functional theory. At 260 K, the models give the following picture of the electronic states as a function of caesium/electron concentra- tion: The dilute solution behaves like an electrolyte in which the electrons exist as polarons, on average spher- ical states localized in solvent cavities, far from the counter ions. At ca. 0.5 mol% metal (MPM), the solvated electrons are spin-paired and form large peanut-shaped species called bipolarons. The electrons still exist as bipolarons at ca.1 MPM but well separated from each other. The bipolarons exhibit no strong tendency to be oriented either parallel or perpendicular to each other.We have investigated the effect of system size on the structure and dynamics of the bipolarons at ca. 1 MPM and found it to be small. At higher concentration, ca. 2 MPM, the electrons exhibit a tendency to cluster and the electron density oscillates between localized and delocalized states. At much higher concentration, ca. 9 MPM, the solution behaves as a good liquid metal in which the electron density forms multi-tunnel-like extended states. The caesium cations are always solvated by the ammonia and are thereby isolated from close contact with the electron density. The role of the cations is assessed through a comparison of the results of the two models. At low metal concentration, the effect of the cations turns out to be rather small.However, the explicit inclusion of the ions is found to increase the metallic character of the solution at ca. 9 MPM. Our findings rationalize a large body of experimental data on this system. Metal-ammonia solutions are a classic chemical system that has attracted much interest and study.'-' The solutions exhibit a wide variety of phenomena as a function of metal concentration. However, for the most part, the effect of differ- ent metal counter ions is found to be rather small. At very low concentration, i.e. MPM, the solutions are char- acterized by isolated excess electrons called polarons. At higher concentrations, <2 MPM, the solutions are domi- nated by localized spin-paired species, referred to as bipol- arons.' The metal-insulator transition occurs at ca.4 MPM and, for concentrations greater than ca. 9 MPM, the solution behaves like a good liquid metal. Unlike other metal-ammonia solutions, caesium and ammonia are completely miscible and there is no two-phase coexistence region (miscibility gap). The solubility of caesium in ammonia is the largest among different metal solutes. Although the available experimental information is extensive, the microscopic under- standing of the electronic states of the system is mostly qual- itati~e.~-'~Little is known about the specific role played by the metal cations. The evolution of spin-pairing in dilute metal-ammonia volume of an isolated electron,' observations that suggest that the association is weak.In order to explain these results, Mott has proposed an electron distribution resembling that of a hydrogen molecule with a binding energy of a few tenths of an eV above that of two isolated electrons.' Previous simulation studies employing path integral Monte Carlo (PIMC) and Car-Parrinello local spin density functional (CP-LSDA) meth~ds'~" have confirmed Mott's inferences about the electronic states at ca. 1 MPM. In the singlet state, the electrons spin pair and form peanut-shaped cavities with peaks in the electron density about 6.5 8, apart.l4.' These spin-paired species are the so-called 'bipolarons'. A representative configuration of a singlet-state bipolaron in a periodic box of 256 ammonia molecules is shown in Plate 1.Despite this success, several issues remain to be explored. For example, how is the stability of the bipol- aron affected when the metal counter ions are explicitly included? Is there a finite size effect on the structure and dynamics of the bipolaron? To what extent are the bipol- arons associated at 1 MPM? Is there any difference in the properties of bipolarons at different metal concentrations? solutions has been extensively studied e~perimentally.',~,~,~ For metal concentrations above lop3 MPM, there is a marked drop in the spin susceptibility per electron. The con- ductivity per electron also falls off and reaches a minimum at ca. 0.1 MPM.' These observations have been interpreted to indicate that the electrons form associated spin-paired species such that at ca.0.1 MPM, spin-pairing is essentially complete with almost 80% of the electrons spin-paired.' However, the optical absorption spectrum of the solution is relatively unchanged from that at low metal concentration and the volume of the spin-paired species is approximately twice the In this paper, these issues are addressed with specific refer- ence to caesium-ammonia solutions. Two different models of caesium-ammonia solution are considered. l6 In the first model (I), the ammonia solvent is treated as a partially pol- arizable simple point-charge model and the cations as a uniform positive background.I6 In the second model (11), the ammonia solvent is made fully polarizable and the cations, here caesium, are explicitly included.The Car-Parrinello LSDA method is used to examine the electronic states in caesium-ammonia solutions at four different concentrations : 0.5, 1, 2 and 9 MPM. This detailed study, which complements our previous work using the first m~del,'~*'~*~~ will allow us to assess the role of the alkali-metal cations in determining the properties of the metal-ammonia solution. In addition, a possible system size effect on our previous low-concentration st~dies'~,'~is examined. To this end, a larger simulation cell (ca. 29 8, edge) containing 512 ammonia molecules plus two or four caesium atoms is employed to study the 0.5 and 1 MPM solutions. The structure, energetics and dynamics of the bipolarons obtained at these concentrations are com- pared to those obtained in earlier at 1 MPM using a smaller cell (ca.22 8, edge) containing 256 ammonia mol- ecules and two excess electrons. The optical conductivity and the qualitative nature of the electronic states for the two models are compared to each other and experiment. Antici- pating our results we will see that the most detailed model gives an excellent account of properties of the metallic solu- tion. Potential Models and Methodology In the two models of caesium-ammonia solution presented in this paper, the ammonia molecules are treated using the simple rigid point-charge ansatz.'6" This is a reasonable approximation as the ammonia molecule is known to retain its gas-phase geometry throughout the metal-ammonia phase diagram.The electronic degrees of freedom associated with the valence electrons are taken to interact with the ammonia and the alkali-metal cations via pseudo-potentials. Specifi- cally, the pseudo-potential used to describe electron-ammonia interactions contains the appropriate electrostatic interactions and both polarization and repulsion contribu- tions from the nitrogen.'6*20 For the e--Cs+ interaction, the s-wave pseudo-potential compilation of Bachelet is used as a local potential.'6.2 ' This simplification is justified because, in caesium-ammonia solutions, the caesium cations are solvated by ammonia molecules and isolated from close contact with the electron density.' LSDA is used to treat electronic exchange and correlation.' 6*22-24 The dipole polarizability of the ammonia nitrogen is treated in two ways.16 In the par- tially polarizable model, the electric field due to all the other ammonia molecules at a given nitrogen atom is assumed to vanish.The electronic contribution to the polarization energy is treated as first order in the polarizability using a pairwise additive scheme. In the fully polarizable model, both the ammonia-ammonia and ammonia-metal-ion contributions to the polarization energy are treated fully self-consistently and the electronic part is, again, treated using the pair approximation. In both cases, electron-electron terms in the polarization energy are neglected.The details of the poten- tials, as well as the treatment of the many-body polarization energy, are described in detail elsewhere.16 As in previous work, the Cs+-ammonia pair potential is taken to be of the following form2' The potential parameters were obtained from a fit to the position and value of the minimum total energy of the Cs+(NH3) complex determined by ab initio calculation^.^^^^^ The parameters used in eqn. (1) for both partially and fully polarizable ammonia can be found in Table 1. The Cs+-Cs+ interaction is described by the familiar Tosi-Fumi potential27 zi zj e2 vM+-M+(rij) = -+ B exp(-olrij) --c6--c8 (2)r6. r?.'ij V V J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 Table 1 Cs+-ammonia potential parameters model ion rM+&A V(r)/K C,,/Ehah2 C,/E,ag ro/aoa I Cs+ 3.0444 -8428 12764282.14 297.45 0.5 I1 Cs+ 3.0444 -8428 12683035.01 249.09 0.5 Table 2 Cs+-Cs+ potential parameters BIG u/a; C6/Eh C8/Eh 1923.9351 1.8765 158.845 1037.463 where the parameters are listed in Table 2.The polarizability of the caesium cations themselves is neglected. The effect of this approximation remains to be investigated. The electronic states of the system were determined using LSDA.22y23 In this scheme, the energy of an N-electron system is written as t2 r where n(r) = nl(r) + n-,(r), n,(r) = 1;:,I t+biU(r) and N,1' + N-, = N. The self-consistent equations for the Kohn- Sham (KS) orbitals are then =Ei $ia(r) (4) where the exchange correlation energy is defined as &(r) = drn(r)[e,(r) + ec(r)].The forms of the spin-polarized e,(r) and E,(r) are taken from Perdew and Z~nger.~~ The lowest spin state was studied in the present paper. The single-particle orbitals have been expanded in a plane- wave basis set using a cut-off of k,, = 124L for the small simulation box (L x 22 A), and k,,, = 18n/L for the large simulation box (L z29 A), where L is the edge of the simula- tion cell. Only the r point is used in the calculations because of the large size of the unit cell and the liquid nature of the system. The plane-wave cut-off was tested on a few represen- tative configurations and found to give energies to within a few per cent of the fully converged results. LSDA gives an expression for the ground-state energy.Thus, in all the calculations presented in this paper, the properties of the system on the ground-state electronic surface were determined. Clearly, this is an approximation in the metallic regime where the energy gap goes to zero. However, the present scheme has been shown to give good results in a variety of similar application^.^^.^^ Additional approximations are used in the evaluation of the optical con- ductivity. Here, we employ the Kubo-Greenwood (KG) rela-ti~n,~'*~l a(o) = -x IM[ XR)I'S[&P(R)-Ep) -ho] ) (5) where o is the spin state, f; is the occupation number of the spin state (a,i), R is the volume of the simulation cell, and J. CHEM. SOC. FARADAY TRANS., 1994, VOL.90 Plate 1 Electron density of a representative configuration of the singlet-state bipolaron at ca. 1 MPE taken from our previous CP-LSDA calculations using model I.14.15 The simulation cell (ca. 22 8, edge) consists of 256 ammonia molecules and two excess electrons. The outermost contour contains 95% of the electron density. Plate 2 Electron density of a representative configuration of a caesium-ammonia solution at ca. 0.5 MPM taken from the present Car- Parrinello LSDA calculations using model 11. The simulation cell (ca. 29 A edge) consists of 512 ammonia molecules and two caesium atoms. The Cs+ ions are shown with covalent radii (pink balls) and ammonia molecules with a ball-and-stick representation. The outermost contour contains 95% of the electron density.Deng et al. (Facing p. 2010) J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 Plate 3 Electron density of representative configurations of a caesium-ammonia solution at ca. 1 MPM taken from the present CP-LSDA calculations using model 11: in (a)the two bipolarons are parallel to each other while in (b)they are perpendicular. The simulation cell (ca. 29 A edge) consists of 512 ammonia molecules and four caesium atoms. The Csf ions are shown with covalent radii (pink balls) and ammonia molecules with a ball-and-stick representation. The outermost contour contains 95% of the electron density. J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 Plate 4 Electron density of representative configurations of a caesium-ammonia solution at ca.2 MPM taken from the present CP-LSDA calculations using model I1: (a) dimerized bipolaronic structure; (b)amoeba-like more extended structure. The simulation cell (ca. 22 A edge) consists of 256 ammonia molecules and four caesium atoms. The Cs+ ions are shown with covalent radii (pink balls) and ammonia molecules with a ball-and-stick representation. The outermost contour contains 95% of the electron density. J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 Plate 5 Electron density of representative configurations of caesium-ammonia solution at ca. 9 MPM taken from the present CP-LSDA calculations using model 11: (a) Cs+ ions are shown with covalent radii and ammonia molecules with a ball-and-stick representation. (b)As in (a)but with space-filling ammonia molecules.The system consists of 256 ammonia molecules and 24 caesium atoms. J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 201 1 Table 3 Simulation parametersa number of concentration /MPM number of valence electrons ammonia molecules box length /A 0.5 2 512 28.4 1 2 256 22.4 1 4 512 28.6 2 4 256 22.4 2 4 256 22.9 9 24 256 23.4 9 24 256 25.1 ' 1 au (mass) = 9.106 x kg. 1 au (time step) = 2.42 x lo-'' s. MT = (Y: I@IYq) is the momentum-operator matrix element between states (0,i) and (a,j). The latter is evaluated using the spin-up and spin-down single-particle KS states instead of the 'true' many-body eigenfunctions and eigen- values. Clearly, the use of the LSDA single-particle 'excited' states, represents an additional approximation as LSDA is a theory for the ground state.The average is over a canonical distribution of the nuclear positions (R)on the ground-state electronic surface. The structure and dynamics of the system were determined In ourusing the Car-Parrinello meth~d.~~,~~model of metal-ammonia solutions, the ammonia molecules and metal cations move on a Born-Oppenheimer (BO) energy surface generated by the electrons and the induced dipoles UdiRiI) = min (V;,C(R,), PJI + min {E,sC{Ri), $ill (6) Pi(a) Jli That is, the electrons must be in the instantaneous ground state and the induced dipoles must satisfy the minimization condition at every solvent configuration.16 The Car-Parrinello approach is applied to this system by treating the plane-wave expansion coefficients of the orbitals in the LSDA expression for the ground-state energy [cf.eqn. (3)] as dynamical variables and assigning them a small fictitious mass, me,and a very low temperature, K.32*33Similarly, the components of the induced dipoles on the nitrogen atoms are also treated as dynamical variables and assigned another small fictitious mass, mp, and a low temperature, Tp.34*35 The system is propagated according to a Hamiltonian consisting of the kinetic energy of the ions and ammonia molecules, the kinetic energy of the basis-set parameters, the kinetic energy of the induced dipoles, the pair additive portion of the ammonia-ammonia, ion-ion and ammonia-ion interactions, the electronic energy and the many-body polarization energy of the solvent.16 The fictitious masses, me and mp,have no physical meaning, but control the timescale of the motion of the electrons and the induced dipoles, respectively.In prac- tice, the fictitious masses of the basis-set parameters and the induced dipoles are adjusted for each concentration studied until an adiabatic separation of the dynamics of these fast variables and solvent is attained (see Table 3). Simulations were run for 30 ps for each state point studied using the velocity Verlet integration alg~rithm.~~ The large temperature differences between the slow ammonia molecules and metal ions (at 260 K) and fast but cold (T < 5 K) degrees of freedom associated with the induced dipoles and basis set parameters, are maintained using a modification of the Nose-Hoover canonical dynamics scheme, the MKT chain meth~d.'~*~~*~~Here, independent chains of thermostats are placed on each orbital, the induced dipole moments, the metal cations and the translational and rotational degrees of solvent molecules.16 The MKT chains were found to main- tain temperature control very well and, as a result, deviations model time step/au' mass, m,/aub mass, m,,/sub T, Tp K/K I1 8 64 600 0.02 5.0 I 16 256 600 0.02 5.0 I1 8 64 600 0.02 5.0 I 8 64 600 0.02 5.0 I1 8 64 600 0.02 5.0 I 4 16 600 0.02 5.0 I1 2 4 600 0.02 5.0 from the Born-Oppenheimer surface were observed to be <5% during a typical 30 ps trajectory.For higher concentrations (2 and 9 MPM), a smaller simu- lation cell (ca. 22 A edge) which consists of 256 ammonia molecules plus four and 24 caesium atoms, respectively, was used. Previous employed the simple model and this smaller simulation cell. The volume of the simulation cell is taken to be:' v,= VNHs + n, v,+ n, v, (7) where VNHJis the volume of 256 (512) ammonia molecules for the smaller (larger) simulation cell at a typical liquid density, pNHJ= 0.023 A-3, and the parameter V, = 65 A3 is the excess volume per electron and r/; = 125 A3 is the excess volume per Cs+ ion.' Results 0.5 and 1 MPM The caesium-ammonia solutions were examined at concen- trations corresponding to ca.0.5 and 1 MPM using the second model (fully polarizable ammonia and explicit caesium cations)16 and a large simulation box (ca. 29 A). At ca. 0.5 MPM, the valence electrons of the caesium atoms ionize and then associate to form a peanut-shaped bipolaron (see Plate 2). The distance between the average position of spin-up and spin-down density (see Fig. 1) is about 6.5 A tips Fig. 1 Distance between the average position of spin-up and spin- down density, as a function of time; for the singlet-state bipolaron in CQ. 0.5 MPM caesium-ammonia solution 2012 Table 4Average energes per bipolaron/E, 1 MPM I -0.19 0.18 0.06 -0.16 0.5 MPM I1 -0.20 0.19 0.06 -0.16 1 MPM I1 -0.19 0.18 0.06 -0.16 which is about the same as that of the 1 MPM singlet-state bipolaron previously studied.' 49L The bipolaron is observed to fluctuate and reorient.At times, the two electrons move closer than 6.5 8, with a concomitant shape change from peanut-like to spherical. However, when the electrons 'rebound' to the peanut shape, the bipolaron is oriented dif- ferently. The electrons at ca. 0.5 MPM never collapse to a single spherical cavity, but remain at least 3.5 8, apart (see Fig. 1). The caesium cations are solvated in the bulk ammonia solvent and are not in close contact with the elec- tron density. The average distance between the two caesium cations is about 12 A. The self-diffusion coefficient of caesium ions is about 0.2 A2 ps-l, less than half of that of the ammonia solvent, 0.5 A* ps-'.At ca. 1 MPM, the electron density is still localized. The electrons exist as separate bipolarons, see Plate 3. Each bipol- aron is very similar to those found using the simple model and a smaller simulation cell. 14*' The distance between the average position of spin-up and spin-down density in each bipolaron is about 6.5 A. The two bipolarons are well separated and constantly change their shape and orientation, moderated by the solvent fluctuation. The bipolarons sample all possible relative orientations, including parallel [Plate 3(a)] and perpendicular [Plate 3(b)] configurations. The caesium cations are, again, solvated in ammonia and isolated from close contact of the electron density. In order to quantify and compare the bipolarons obtained from different systems further, the average energy of the bipolaron and its component parts, V&3, V,,, V,, and kinetic energy (Ek)are presented in Table 4.The energetics of different bipolarons are almost identical.2and 9 MPM Higher-concentration caesium-ammonia solutions were also studied using the CP-LSDA method. At ca. 2 MPM, the electrons begin to associate and form clusters. Two types of cluster are observed. The first, shown in Plate qa),is a local- ized form, a dimer of bipolarons. This state occurs in 80% of the configurations. In Plate 4(b), the second cluster type which accounts for the remaining configurations is shown. Here, the electron density is delocalized in a tube-like struc- ture which spans the simulation cell.At ca. 9 MPM, the elec- tronic structure of the system is radically different (see Plate 5). The electrons are now delocalized in a multi-tunnel-like structure which spans the entire simulation cell. That is, the electron density remains excluded from the ammonia solvent, as at the low concentrations, and forms a bicontinuous con- struct. As tested by methods described elsewhere,' the7728 individual orbitals are also delocalized. No localized bipol- aronic structures are observed at this concentration. The single-particle electronic density of states, N(E), calcu-lated by averaging over KS eigenvalues [cf: eqn. (4)] of 40 configurations well separated in time, is presented as a func- tion of concentration in Fig.2. The gap narrows and closes as the metal concentration of the solution is increased. In Fig. 3, the optical conductivity is shown as a function of concentra- tion. At ca. 2 MPM, the extrapolated dc conductivity is very small (recall that the metal-insulator transition occurred at ca. 4 MPM). However, at ca. 9 MPM, a dc conductivity of J. CHEM. SOC. FARADAY TRANS., 1994, VOL. 90 LV I OEb!6'' IOi4' ' $!2' ' ' !I ' ' "012' 'Oi4"I EIeV Fig. 2 Single-particle electronic density of states, N(E), calculated by averaging over the Kohn-Sham eigenvalues of 40 configurations well separated in time using model 11. The Fermi energy defined as the eigen energy of the highest occupied orbital has been set to zero for each configuration.(a)2, (b)9 MPM. I-k 12345 E/eV EP Fig. 3 Optical conductivity of caesium-ammonia solutions at 260 K. (a)2, (b)9 MPM. (-) model 11, (---) model 1. about 2000 R-' cm-' is obtained by extrapolating our results to zero frequency, which is in reasonable agreement with the experimental value, 1500 R-' cm-'.' As studied above, at both 2 and 9 MPM, the Cs+ ions are, always, solvated by the ammonia and isolated from close contact with the electron density. The Csf-N radial distribu- tion function, which is presented in Fig. 4 (solid curve), indi- 4 6 8 10 12 r/A Fig. 4 Radial distribution functions for caesium-ammonia solution at ca. 9 MPM using model 11. (-) Cs+-N and (---) Cs+-Cs+ g(r) functions. J. CHEM. SOC.FARADAY TRANS., 1994, VOL. 90 cates that the first solvation shell for Cs' consists of about eight ammonia molecules. It was found that the local many- body polarization effects reduce the coordination number of Cs+ by about one ammonia molecule. The Cs+-Cs+ radial distribution function for caesium- ammonia solution at ca. 9 MPM is also shown in Fig. 4 (dashed curve). The first peak at about 4.5 A is a contact ion pair and the second peak at about 9 A is a solvent-separated ion pair. The area under the first peak gives, on average, one contact ion pair per Cs' ion, which suggests that the Cs' ions are dimerized. It remains to be ascertained whether or not this is a real effect or an artifact arising from neglect of the polarization of the Cs+ cation.Conclusion In summary, CP-LSDA calculations have been performed on caesium-ammonia solutions for a wide range of concentra- tions using two different microscopic models.16 In model I, the ammonia molecules are treated via a classical point- charge model and the cations as a positive neutralizing back- ground. In model 11, the ammonia solvent is made fully polarizable and the cations, here caesium, are explicitly included. The results were used to determine the role of the cations in the solutions and to assess finite size effects in our studies. The effect of the counter ions on the solutions was deter- mined by a comparison of the two models described above. At low concentrations (1 and 2 MPM), only minor differences were found.However, the effect of the counter ions was quite visible at high concentrations. Explicitly including the Cs + cations increases the metallic character of solutions at ca. 9 MPM. Specifically, the 'energy gap' decreases by a factor of seven and the extrapolated dc conductivity increases by about a factor of three. In fact, the optical conductivity calcu- lated for the system with the explicit ions is in reasonable accord with e~perirnent,'.~~ a finding which adds credibility to the present calculations. Finite size effects at low concentrations were examined using the most detailed model. The bipolarons at ca. 0.5 MPM are very similar to those at ca. 1 MPM. Essentially, no difference was found between studies at 1 MPM using a simulation cell with 256 ammonia molecules plus two excess electrons and a larger cell with 512 ammonia molecules and four caesium atoms.We would like to thank John Shelley for indispensable help in creating the three-dimensional graphics. Also, Z.D. thanks Cray for a Research Fellowship and G.J.M. thanks NSF for a Postdoctoral Research Associateship in Computational Science and Engineering (ASC 90-08812). The research described herein was supported by the National Science Foundation under CHE 92-24536. Some of the computations benefited from use of facilities provided by NSF/DMR 91- 20668. The bulk of the computing was performed at the Pitts- burgh Super Computing Center under grant MCA 933020. 201 3 References 1 J. C. Thompson, Electrons in Liquid Ammonia, Oxford Uni- versity Press, London, 1976.2 Physics and Chemistry of Electrons and Ions in Condensed Matter, ed. J. Acrivos, D. Reidel, New York, 1984. 3 Metal-Ammonia Solutions, ed. G. Lepoutre and M. J. Sienko, W. A. Benjamin, New York, 1963. 4 Electrons in Fluids, ed. J. Jortner and N. R. Kestner, Springer- Verlag, Berlin, 1973. 5 Colloque Weyl IV, J. Phys. Chem., 1975,79,2789. 6 Colloque Weyl V, J. Phys. Chem., 1980,84, 1065. 7 Colloque Weyl VI, J. Phys. Chem., 1984,88,3699. 8 Colloque Weyl VII, J. Phys. ZV,1991, 1. 9 N. F. Mott. Metal-Insulator Transition, Taylor & Francis, London, 2nd edn., 1990. 10 J. Jortner and M. H. Cohen, Phys. Rev. B, 1976,13, 1548. 11 K. Rapcewics and N. W. Ashcroft, Phys.Rev. B, 1991,44,4032. 12 M. Sprik, R. Impey and M. L. Klein, Phys. Rev. Lett., 1986, 56, 2326. 13 M. Marchi, M. Sprik and M. L. Klein, J. Phys. C, 1990,2, 5833. 14 Z. Deng, G. J. Martyna and M. L. Klein, Phys. Rev. Lett., 1992, 68,2496. 15 G. J. Martyna, Z. Deng and M. L. Klein, J. Chem. Phys., 1993, 98, 555. 16 Z. Deng, G. J. Martyna and M. L. Klein, J. Chem. Phys., 1994, loo. 17 Z. Deng, G. J. Martyna and M. L. Klein, Phys. Rev. Lett., 1993, 71, 267. 18 M. Sprik, R. Impey and M. L. Klein, J. Am. Chem. SOC., 1986, 56,2326. 19 F. Leclerq, P. Damay and P. Chieux, J. Phys. C, 1991,5, 150. 20 G.J. Martyna and M. L. Klein, J. Phys. Chem., 1991,95, 515. 21 G. Bachelet, D. Hamann and M. Schluter, Phys. Rev. B, 1982, 26,4 199. 22 Theory of the Inhomogeneous Electron Gas, ed.S. Lundqvist and N. H. March, Plenum Press, New York, 1983. 23 R. 0.Jones and 0.Gunnarsson, Rev. Mod. Phys., 1989,61,689. 24 J. P. Perdew and A. Zunger, Phys. Rev. B, 1981,23,5048. 25 S. F. Smith, J. Chandrasekhar and W. L. Jorgensen, J. Phys. Chem., 1982,86,3308. 26 J. Tse, personal communication, 1985. 27 M. J. L. Sangster and M. Dixon, Adv. Phys., 1976,25,247. 28 E. Fois, A. Selloni and M. Parrinello Phys. Rev. B, 1989, 39, 48 12. 29 R. Car G. Galli, R. Martin and M. Parrinello, Phys. Rev. Lett., 1989,63,988. 30 F. Wooten, Optical Properties of Solids, Academic Press, New York, 1972. 31 E. J. Mele, personal communication, 1992. 32 R. Car and M. Parrinello, Phys. Rev. Lett., 1985,55, 2471. 33 D. K. Remler and P. A. Madden, Mol. Phys., 1990,70,921. 34 M. Sprik and M. L. Klein, J. Chem. Phys., 1988,89,7556. 35 M. Sprik, J. Phys. Chem., 1991,95,2283. 36 W. C. Swope, H. C. Andersen, P. H. Berens and K. R. Wilson, J. Chem. Phys., 1982,76,637. 37 G. J. Martyna, M. Tuckerman and M. L. Klein, J. Chem. Phys., 1992,97,2635. 38 G. J. Martyna and M. L. Klein, Liouville Operator Analysis of the Car-Parrinello Method, preprint, 1994. 39 M. Schlauf and R. Winter, 2. Phys. Chem., 1988,156,207. Paper 4/00759J; Received 7th February, 1994