Dynamics of the Carbon Nuclei in C60 Studied by Feynman Path

J. Phys. Chem. 1995,99, 12401-12408

12401

Dynamics of the Carbon Nuclei in c60 Studied by Feynman Path-Integral Quantum Monte Carlo Simulations Michael C. Bohm*" Fritz-Haber-Institut der Max-Planck-Gesellschaji, Faradayweg 4-6, 0-14195 Berlin, Gemtany

Rafael Ramirez" Instituto de Ciencia de Materiales, Consejo Superior de Investigaciones Cientijicas, E-28006, Madrid, Spain Received: January 26, 1995; In Final Form: June 9, 1995@

The dynamics of the carbon nuclei in the CW molecule have been studied by Feynman path-integral (PI) quantum Monte Carlo (QMC) simulations and within a harmonic oscillator approximation. The following finite-temperature properties have been calculated: kinetic, potential, and total energy, radial and angular distribution functions (rdf, adf) as well as parameters describing the degree of atomic deflocalization. D, measures the overall spatial atomic fluctuations, D, the associated quantum, and D,the classical (=thermal) contributions. PI simulations render possible the derivation of the quantum delocalization DS. At T = 0 K this quantity measures the amplitude of the zero-point vibrations which correspond to the spatial uncertainty of the atoms according to the Heisenberg uncertainty principle. Comparison of results derived by quantum simulations and classical ones demonstrates the importance of quantum effects in the whole temperature range below room temperature. The possible influence of quantum fluctuations on microscopic electronic structure properties of superconducting fullerides is briefly mentioned. PI simulations as well as simple analytical results show that the spatial uncertainty of the C atoms in CW is comparable to the difference between the lengths of the short and long CC bonds in the icosahedral molecule. Consequences of such large zero-point motions over a large temperature interval for the microscopic interpretation of diffraction data are concisely summarized. Reliable agreement between the PI-derived rdf and the experimental curve based on neutron diffraction data is observed. The interaction between the C atoms has been modeled by the potential of Tersoff which has been developed to reproduce several properties of fullerenes such as binding energy and geometrical parameters.

I. Introduction The molecular and solid state properties of fullerenes and fullerides have absorbed the interest of many chemical and physical laboratories. The fascinating material properties of this new modification of carbon have been summarized in detailed review articles covering both experimental and theoretical topics.'-3 On one side, the majority of the chemical work has been dedicated to the characterization of molecular properties. The analysis of any collective phenomena encountered in fullerides such as low-temperature superconductivity in alkaliand alkaline-doped c 6 0 materials has been left aside. On the other side, it has been this low-temperature superconductivity in the first place which called forth comprehensive work in physical laboratories. Fulleride superconductors differ from many conventional ones by the lighter mass of the corresponding atoms. Note that superconductors usually are formed by atoms from the lower part of the periodic system of the elements. Conceming the atomic masses encountered in superconducting materials, similarities between fullerides and organic charge transfer (CT) salts which have been characterized in the past two decades become e ~ i d e n t . ~ ~ ~ Different pairing mechanisms have been suggested for these nonconventional superconductors to account for the differences between organic materials on one side and conventional superconductors on the other. For fullerides theoretical pairing

' Permanent address:

Institut fur Physikalische Chemie, Physikalische Chemie 111, Technische Hochschule Darmstadt, Petersenstrasse 20, D-64287 Darmstadt, Germany. Abstract published in Advance ACS Absrracfs, July 15, 1995.

*

0022-365419512099-12401$09.00/0

models have been discussed which cover the way from the standard electron-phonon (el-ph) coupling of the BardeenCooper-Schrieffer (BCS) type via a superposition of el-ph and electrodhole pairing up to bare charge carrier It is common among the different phenomenological models that the dynamics of the atomic nuclei have been ignored. As a result of the light carbon mass in fulleride superconductors, as well as in other organic representatives, it can be expected that quantum effects such as zero-point motions may be of large importance. In an improved theory of superconductivity the possible influence of the zero-point vibrations on microscopic (electronic) parameters of the system (Le., electronic density of states (DOS)at the Fermi energy ( E F ) , Fermi surface (FS) properties, el-ph coupling constants, etc.) should be considered. Although a number of simulation studies of the molecular dynamics (MD) and Monte Carlo (MC) type on the finitetemperature properties of fullerides can be found in the literature, quantum effects associated with nuclear degrees-of-freedom have been neglected completely. For previous classical simulation work see refs 12-17. An exception in this context is an investigation by Kohanoff et al.I* where zero-point motion effects of Cm have been studied by the Car-Parrinello method.Ig The generated trajectory with this degree of sophistication is a classical one. Quantum effects have been modeled in that contribution within the harmonic approximation. It is the purpose of the present investigation to study finitetemperature (T 0 K) properties of the C ~ molecule O which are related to the dynamics of the carbon nuclei by a quantum approach. We believe that the theoretical treatment of a single molecule may be an important precursor to improving our

*

0 1995 American Chemical Society

12402 J. Phys. Chem., Vol. 99, No. 33, 1995 understanding of possible collective solid state properties of fullerides. Detailed work has demonstrated that mainly intramolecular modes are decisive for the generation of many collective properties in this class of compounds.z0 We have employed the Feynman path-integral (PI) quantum Monte Carlo (QMC) approach as a theoretical too1.2’,2z This method has become a powerful simulation technique to evaluate statistical mechanical properties of quantum systems. The range of applicability of PI techniques has been described in several review^.^^-^^ In this work we adopt the PI formalism to calculate the total energy of the Cm molecule, the radial and angular distribution functions (rdf, adf), and the deAocalization parameters of the atomic nuclei which are related to zero-point motions. The PI results will be correlated both with available experimental data such as the rdf measured by neutron scattering and with theoretical data from a simpler analytic model based on the harmonic oscillator (ho) description. In the PI formulation employed it is straightforward to compare results of quantum simulations with those of classical ones. Hence, it becomes possible to discriminate between quantum degreesof-freedom and classical (Le., thermal) ones contributing to the dynamics of the c 6 0 nuclei. In previous studies we have used PI techniques to investigate the quantum delocalization in crystalline silicon27and the quantum behavior of hydrogen or muonium in siliconz8and to construct the phase diagram of the one-dimensional (1D) Hubbard chain (=electronic degrees-offreedom).29 Bosonic degrees-of-freedom have been studied in refs 27 and 28, while a fermionic problem has been considered in ref 29. For a preliminary PI investigation of Cm see ref 30. Apart from the actual interest in the material properties of fullerides we additionally want to emphasize two points of more general interest. The fiist one concems the dynamics of carbon atoms in n systems with bond length altemation. Cm is a nonaltemant molecule3’ which belongs to the icosahedral point group I h . The carbon atoms form 12 pentagons and 20 hexagons with 2 topologically different carbon-carbon bonds. Thirty short bonds are shared between neighboring hexagons [(6-6) = C=C] and 60 long bonds are formed by pentagon-hexagon contacts [(6--5) = C-C]. The mean values of the two bond lengths measured by different methods (e.g., X-ray, NMR, and neutron and electron diffraction) amount to 140.2 and 144.6 pm.32 This bond length altemation of 4.4 pm will be of importance in the later discussion. The value is characteristic for many altemant and nonaltemant hydrocarbon n systems. In section I11 we discuss the spatial fluctuations of the C atoms of Cm around the respective mean values. As will be demonstrated in detail, the spatial distribution functions confined to the two bonds show remarkable overlap. The second issue of major importance is the quantitative understanding of atomic fluctuations in crystals which is essential for the proper interpretation of diffraction data. These diffraction techniques are very successful in the determination of averaged structures in the crystallographic unit cell. On the other hand it is very difficult to derive trustworthy information concerning the microscopic origin of spatial atomic fluctuations by these methods. The fluctuations are caused by dynamic effects as well as static disorder. Path-integral simulations of the present type are very useful for discriminating between both effects. The organization of the present contribution is as follows. In section I1 the theoretical background of the PI approach is summarized concisely. In the next section the computational results are discussed and compared with experimental data (rdf of c60). The article ends with a final resume in section IV.

Bohm and Rm’rez 11. Theoretical Background

Since several general reviews on the theoretical basis of the Feynman PI QMC approach can be found in the literature, we will only summarize the main principle^.*^-^^ A key element of the PI decomposition described below is the isomorphism between a single quantum particle (here, each C nucleus of the Cm molecule) and a cyclic chain of N classical particles. For an ensemble of many quantum particles we thus have one separate cyclic chain for each. This isomorphism allows the use of conventional classical simulation methods, e.g., MD or MC to calculate the thermal expectation values of the original quantum system. The associated Hamiltonian reads

H=

-z[-[;+;+$)] h2 a2 a2 p=l

2mp axp

ay,

+ V(R)

(1)

az,

where V(R) is the potential acting on the quantum particles (here, atomic nuclei). mp is the atomic mass of the particle p . R symbolizes a vector in 3P dimensional space where P stands for the number of atoms. The components of R are the Cartesian coordinates rpof the atoms; see eq 2.

R = (rl, r2, ..., r,)

(2)

In the present work the interaction between the C atoms V(R) is described by the empirical Tersoff potential,33which has a simpler analytical structure and is better documented in the literature than other model potentials for carbon polytypes. It works best for 4-fold coordinated carbon. The accuracy of the Tersoff potential is at least comparable to the parametrizations of Brenner34or C h e l i k o w ~ k ywhich ~ ~ have been developed for fullerenes. In eq 3 we define the partition function of the quantum system Z by the trace of the statistical density matrix @(R,R’;/3).This quantity can be considered as the propagator for the evolution of the system in imaginary time; see below

Z = Tr[exp(-/3H)] = Tr[@(R,R’$)]

(3)

/3 is given by ( k ~ T ) - lwith the Boltzmann constant k

~ As . is well known, the quantity $3 h. has the character of an imaginary time in quantum simulations. The PI evaluation of Z requires the discretization of the statistical density matrix along cyclic paths which are fragmented into N (time-) slices. This decomposition reads

Note the cyclic structure of the above expression. Starting from a spatial distribution of atoms described by the vector R1, the first vector moves the ensemble to R2, the second to R3, etc. The Nth factor finally retums the ensemble to Rl. Equation 4 becomes exact in the limit N 00 where the partition function can be expressed in form of a so-called shorttime approximation (iPhlN is the size of each time-slice). For N the “true” many-body potential V(R) does not have sufficient time to feel itself. We have

-

-

Dynamics of Carbon Nuclei in Cm

J. Phys. Chem., Vol. 99, No. 33, 1995 12403

In eq 5 and in the following equations we have omitted the index p in the atomic masses because there is only one atomic species in the fullerene molecule we have studied. In analogy to eq 3 the 3 P dimensional vectors Rj+l and Rj are defined as P

on the harmonic oscillator model. The overall atomic fluctuations Dt, including quantum and classical (=thermal) degreesof-freedom (D, and DJ,are defined by the root-mean-square deviation of the atomic coordinate rp around its respective equilibrium value (rep).

p= I

The 3D vector rp, in (6) labels the position of the carbon atom p and the index j refers to the coordinate along the path. Definition 5 for Z has the same analytical form as that of a classical distribution function for a cyclic chain of N images (=beads) coupled through harmonic forces with an effective force constant k = Nm/P2h2. k has its origin in the kinetic energy operator of the Hamiltonian. The potential acting on each bead of the chain corresponds to 1IN of the “true” manybody potential V(R) occuring in eq 1. As can be seen, quantum effects facilitate the penetration of the nuclei into regions of high potential. This “tunneling” of the quantum particles is made possible by the renormalization of V(R) in eq 5. It is obvious that the result of a classical approach can be obtained by setting N = 1. In this case any quantum delocalization is prevented. The index j in eq 5 has the following meaning. Within each of the P chains the nucleus at time-slice j is “harmonically coupled“ to i t s j 1 and j - 1 images. The coupling between chains via the potential V(R) requires the same imaginary time argument j , a condition which guarantees an instantaneous interaction. In eq 5 distinguishable particles have been assumed, Le., the many-particle calculations are of the bosonic type. PI simulations of many-electron systems (=fermionic degrees-offreedom) with their characteristic sign problems are still difficult to implement.36 For a previous theoretical approach to treat the sign problem in PI simulations of monocyclic networks (Le., the 1D Hubbard chain with cyclic boundary conditions) see ref 37. The thermal expectation values of t h e c a molecule have been derived by using the classical Metropolis Monte Carlo sampling technique.38 The partition function Z has been sampled according to eq 5 with V(Rj) described by the above Tersoff potential. In section I11 we discuss the radial and angular distribution functions of Cm. The rdf g(r) is given by

+

P

g(r) = &(R;P)&R

-

Ir, - rql)drl...drp

Pf4

(7)

In eq 10 we have employed the following abbreviations: (i) angular brackets represent averaged values at any given temperature; (ii) index c in the definition of rcPsymbolizes that we have evaluated the center of gravity of the path coordinate for nucleus p in the Cm molecule. The two quantities occuring in the square-root of eq 10 are defined as N

j= I

and N

In eq 13 D,is related to D, and Dc:

D,=

,ID: + 0,’

(13)

In the harmonic oscillator approximation simple analytic formulas have been derived for 0:”and D p ; see eqs 14 and 15. This approach leads to a simple access to the quantum delocalization

0:“=

Jx

&)

2mw coth( 2kBT

Equation 14, for example, has been used in the Car-Parrinello study of ref 18. In the next section we will analyze the “experimentally” deduced delocalization parameters DP,Dr, and DY of the Cm molecule as derived in the harmonic oscillator approximation. The “experimental” origin is due to the use of the measured phonon DOS eeXp(w)in eqs 14 and 15.39 The analytically accessible quantum delocalization DP in the ho approximation has a simple pendant in the PI formalism. The corresponding quantity DY is expressed in eq 16:

A similar expression holds for the adf. To check the accuracy of the energy derived by the PI implementation, we compare the theoretical result with the experimental value accessible via the measured phonon DOS distribution ~ ~ ~ ~For( Cwa this ) . ~ ~ quantity has been measured by neutron scattering. w symbolizes the vibrational frequency. Within the harmonic approximation In connection with eq 10 it is of some interest to mention that (ho) E is related to eexp(w)via eq 8. the first terms in the square-roots of eqs 10 and 16 are identical. The definition of the two delocalization parameters, however, E= dw hw[n(w) 0.5]gexp(w) (8) differs in the second term. Df has a well-defined physical meaning in the PI approach.25 It measures the root-mean-square n(w) is the average number of excited phonons of frequency radius of the cyclic path (i.e., the radius of gyration) associated w . In the Bose-Einstein model n(w) is given by with the quantum particle. Physically, this parameter has its origin in the mean extension of the springs in the cyclic chain n ( o ) = exp - - 1 (9) in thermal equilibrium. As a net result, the chain as a whole spreads out. And this effect is measured by eq 16. It should be pointed out that D, can be set equal to D : only in the Below, we introduce quantum and classical (=thermal) atomic harmonic oscillator approximation: D , = DF = D5;’. deAocalization parameters. Alternatively, they are accessible At T = through PI simulations or a simplified analytical approach based 0 K this quantity gives the amplitude of the zero-point motion.

hw

+

[

(3 I-1

12404 J. Phys. Chem., Vol. 99, No. 33, 1995

Of course, it also corresponds to the spatial uncertainty in the Heisenberg uncertainty principle. In order to evaluate numerically stable PI results, the following computational conditions have been employed in the quantum simulations. At each temperature studied we have generated 6 x lo4 quantum paths per carbon atom for the calculation of ensemble-averaged thermal properties. For the system equilibration lo4 quantum paths have been employed. The number of imaginary time-slices has been fixed by the condition NT = 3000 K. In the temperature interval studied this criterion leads to N numbers from N 5 4 to N = 60. With these computational boundaries a convergence in the total energy with an error less than 3% at 50 K is guaranteed. Other finitetemperature properties are feasible with comparable accuracy. All calculations have been performed on workstations of the DEC ALPHA type. At the lowest temperature studied in the quantum simulations (T = 50 K) one PI calculation took more than 40 h of CPU time. The whole theoretical investigation required several days of CPU time. 111. Results and Discussion

The empirical Tersoff potential conserves the icosahedral Ih symmetry of the Cm molecule. It leads to a bond length of 145 pm for the 30 (6-6)C=C “double” bonds and to 150 pm for the 60 (6-5)C-C “single” bonds. So as far as the bond length altemation is concerned, experiment32(see section I) and computational simulation coincide almost quantitatively. It is in the absolute bond length values where the Tersoff potential slightly overestimates the two CC bond lengths by roughly 5 pm. Other model potentials in the l i t e r a t ~ r eyield ~ ~ , about ~ ~ the same overestimation. The calculated elongation in the bond lengths causes a theoretical C a diameter of 737 pm which exceeds the experimental value by 31 pm. With the Tersoff potential we calculate a local breathing mode of 691 cm-I, in which each C atom moves radially from the molecular center, and the so-called pentagonal pinch modes of 1303 and 1325 cm-I. In the latter local modes two different bond lengths are modulated causing a decrease or increase in the area of the c 6 0 pentagons. Comparison with the fully symmetric ag normal modes which develop from the above local modes shows the capability, but also the limitations, of the model potential. The all-in-phase breathing mode is found at 496 cm-I, while the two pinch modes occur at 1458 and 1470 cm-’ (solid state). For previous experimental and theoretical studies on the vibrational properties of Cm we refer to refs 18, 20, 41, and 42. In Figure 1 we have portrayed the temperature dependence of the kinetic, potential, and total energy of the C a molecule in a T window between 0 and 800 K. The three energy elements have been calculated by PI quantum simulations and in the classical limit. As already mentioned, the latter is restored by setting N = 1 in the computational approach. We have adopted the potential energy of the minimum configuration as intemal zero. In an absolute energy scale the Tersoff potential yields -6.73 eV per carbon atom. In the diagram it is seen that the difference between quantum and classical simulations increases with decreasing temperature. Even at room temperature (RT) both quantities differ by roughly a factor of 2. Already this energy difference indicates that quantum effects such as zeropoint vibrations are of particular importance in fullerened fullerides within a rather large temperature interval. The classical simulations reproduce the well-known linear T dependence for all quantities displayed in Figure l. Quantum simulations, on the other hand, yield almost constant values between 0 and 200 K. This behavior demonstrates that only

Bohm and Ram’rez

0.3

i

oe

n

E 0

Y

s

I

,

1

I

1

kinetic

-

0.2

(d

/s

Q

W

0

200

400

800

600

Temperature (K) Figure 1. Temperature dependence of the kinetic, potential, and total energies of the C ~ Omolecule according to quantum and classical simulations. All values in the figure are relative to the potential energy for the minimum energy configuration.

1.061

0

I

I

I

1.04

x

.d Y

3 x

1.02

?ic Q

a ..-( Y

1.

E Q +-)

0

a 0.98I 0

I

I

200

I

1

400

I

1

600

800

Temperature (K) Figure 2. Virial coefficient f calculated for the C ~ molecule O as a function of temperature according to quantum and classical simulations. Additionally, we show f h o = 1 for the harmonic oscillator.

zero-point vibrations are operative in this temperature range. The zero-point energy derived by the model potential of ref 33 amounts to roughly 0.184 eV per atom. This value is approximately 0.02 eV smaller than the experimental zero-point energy which has been derived from the measured phonon DOS under the adoption of the harmonic appro xi ma ti or^.^^ The 9% difference between the calculated and experimentally extrapolated zero-point energy is almost temperature independent. The discrepancy between both quantities may be caused by (i) shortcomings in the potential employed, (ii) experimental uncertainties in the phonon DOS, (iii) limitations in the applicability of the harmonic approximation, or (iv) computational boundaries employed in the PI simulations (Le., number of time-slices). In Figure 2 we have plotted the vinal coefficientf= (V)/(K) as a function of temperature with (V) and ( K ) representing the

Dynamics of Carbon Nuclei in Cm

J. Phys. Chem., Vol. 99,No. 33, 1995 12405

T =300 K classical

v1

Y

0

T=50K classical

n v) .d Y

c

1

1 -

1

0:.

: ’:

d

0.3

0.2

:

: :

0.1

n-

e

: : : :

200

250

130

140

150

160

170

Distance (pm) : : A .

1

T=50K

150

0.4

W

m

100

c 1

2

T= 5K classical

1

. I

300

Distance (pm) Figure 3. Radial distribution function rdf for the C nuclei at distances between 100 and 300 pm. The three upper panels correspond to classical simulations while the bottom diagram is based on a quantum simulation.

thermal averages for the potential and kinetic energy, respectively. In analogy to Figure 1 again we have compared the results of quantum (f,) and classical (f,) simulations. The oscillations in the quantum and classical virial coefficient in the figure are within the error bars of the QMC simulations; they do not represent any physics. For convenience we have added the virial coefficient of the harmonic oscillatorp”. Of course, p” = 1 does hold whatever the temperature is. fc approaches the ho value at T = 0 K. It is enhanced with increasing temperature while f, is a decreasing function of T. Extrapolation offq andfc to higher temperatures shows that both curves converge toward the same limit. Obviously, we find the trivial coincidence offc andp” at zero temperature. In both simulations (v> exceeds ( K ) in the T interval studied. The behavior of the classical curve is easily understandable. With increasing temperature anharmonic effects related to the internuclear potential become more important as the nuclei vibrate with larger amplitudes around their mean positions. The quantum result is more subtle. At high temperatures it converges correctly toward the classical limit. The comparatively large value of the virial coefficientf, at low temperatures is a typical many-particle effect. In our theoretical work we have identified this fact by performing a quantum simulation of only one carbon nucleus, leaving the other nuclei fixed at their equilibrium positions. In this case the fq values for the one-particle problem do not differ sizably from thefc parameters. Therefore, we assume that the large value of f4 a t l o w temperatures has its origin in many-particle effects. Radial distribution functions for CC contacts up to 300 pm as derived in the classical limit (rdf,) and for the quantum description (rdf,) are displayed in Figure 3. In the diagram the nearest-neighbor (1,2) and next-nearest-neighbor (1,3) separations occur around 150 and 250 pm, respectively. As already expected on the basis of the computational results of Figure 1, T dependent rdfs are only realized in the classical simulations. Here, the peaks associated with the short and long CC bonds

Figure 4. Radial distribution function rdf for the C nuclei at distances confined to nearest-neighbor contacts. The sum curve (solid line) has been decomposed into two distributions associated with the short (6-6) C=C bonds (broken curve) and the long (6-5) C-C bonds (dotted curve). The theoretical curve has been calculated via a PI quantum simulation at T = 50 K.

are well resolved in the low-temperature region ( T = 5 and 50 K). With increasing T the width of the distributions is broadened. At RT the overlap of the peaks associated with the short and long CC bonds is even large enough to cause a collapse into one single envelope. Conceming the quantum simulation, sizeable zero-point fluctuations give rise to strongly overlapping distributions for the two nearest-neighbor distances around 150 pm as well as for the two next-nearest-neighbors at distances of about 250 pm for low temperatures. We point out that the quantum rdf is T independent in the interval below 200 K. In the diagram we have selected T = 50 K for the representation of the quantum rdf. In Figure 4 we have displayed the quantum rdf (=rdf,) for nearest-neighbor pairs in an enlarged mesh. Here, we have decomposed the sum curve into two partial rdf curves associated with the short (6-6) and long (6-5) bonds. The strong overlap of the two distribution functions is evident. The mean half-width of the two quantum distributions exceeds 5 pm. Deviations from an idealized Gaussian behavior can be detected in the figure. The fact that the left and right wings of the maximum are not exactly identical can be traced back to deviations from the harmonic oscillator behavior. In Figure 5 we compare the theoretical rdf, curve for an isolated Cm molecule to the experimentally constructed radial distribution function obtained in the solid state by neutron diffraction measurement^.^^ For previous theoretical investigations of the rdf of Cm see refs 13, 18, and 44. The reliability of such a comparison has to be justified. Comparison of the solid state c 6 0 bond lengths (X-ray data) and gas phase electron diffraction values has demonstrated that the lengths of the CC bonds are not changed in the transition from the gas phase to the solid.32 Of course, such a “length conservation” is expected in molecular crystals of the van der Waals type with weak intermolecular interactions. Our previous theoretical analysis o f the Cm molecule and the corresponding solid has shown that the electronic structure of the Cm unit is more or less identical in both casesg Therefore, the correlation in Figure 5 seems to be justified. In the diagram we have compared the experimental radial distribution function and PI-based rdf, curve of Cm. The distribution obtained in the PI simulations has been convoluted to the maximum wave vector Qmaxof 0.45 pm-’ accessible in the neutron diffraction e~periment.4~ Comparison with Figure 3 shows that the fine structure of the second peak around 250

Bohm and Ram’rez

12406 J. Phys. Chem., Vol. 99,No. 33, 1995

0.6

1

1

1

1

1

----- experimental

1

I

2 b

I

I

1

I

1

1

4

T =300 K classical

- 1 rn

I\

T=50K

I

T = 5K classical

100

150

200

250

300

Distance (pm) Figure 5. Comparison of the rdf obtained in the present PI approach (solid curve) with the experimental distribution (neutron diffraction) at distances between 100 and 300 pm. The theoretical and the experimental rdfs have been convoluted to the experimental Qmax parameter

T=50K quantum

0

-0.7

.................................

-0.5

-0.3

-0.1

cos( angle)

E 0

Figure 7. Angular distribution function of the Cm molecule as obtained by classical and quantum simulations. The conditions employed (temperature, quantum and classical setups) correspond to those in Figures 3 and 6.

Y

cd

\

4 C 0 P rcl

0

b s E

z

130

140

150

160

170

Distance (pm) Figure 6. Number of bonds per carbon atom of the Cm molecule as a function of the internuclear separation in the interval between 130 and 170 pm. The sequence of the classical and quantum simulations corresponds to the one of Figure 3. The total integrated number of bonds has been symbolized by a solid curve. The broken curve refers to the distribution function for the short bond per atom, while the dotted curve symbolizes the two long bonds per C atom.

pm (quantum simulation, T = 50 K) is suppressed due to the Qmaxboundary employed. The present PI simulation yields a rdf of good accuracy. The maximum of the experimental and theoretical peaks differ as a result of the 5 pm overestimation of the CC bond lengths by the Tersoff potential. Of course, this deviation becomes larger with larger interatomic distances. Note that the calculated diameter of the C a molecule exceeds the experimental value by roughly 31 pm. Despite the limitations of the simple potential Figure 5 demonstrates the capability of PI simulations to derive T t 0 K properties of quantum systems. Figure 6 completes the information from Figure 3. In the diagram we have displayed the integrated number of CC bonds at a given distance for any C atom of C a in an interval between 130 and 170 pm. Each C nucleus has three nearest neighbors

(solid curve). The distribution function confined to the two long bonds is symbolized by dotted lines in the schematic representation, while the integrated function for the short “double” bond is represented by broken lines. In analogy to Figure 3 the three upper diagrams refer to classical simulations. We have adopted the same T values as chosen in Figure 3. In the bottom panel the integrated number of bonds as derived in the quantum description is portrayed (T = 50 K). The quantum fluctuations of the C nuclei cause the distribution functions for the (6-6)C=C and (6-5)C-C bonds to be strongly r dependent in the same interval even at low temperature. Figure 6 is another clear representation of the strong nuclear dynamics of the C atoms of c60 as a result of zero-point vibrations. Note that the width of the quantum distribution observed at 50 K even exceeds the width of the classical simulation at room temperature; compare this with Figure 3. The angular distribution function adf (Figure 7) shows large broadening due to quantum fluctuations as well. They cause a small overlap between both distribution functions. The two peaks in the diagram correspond to the bond angles within the six- and five-membered rings. In the classical simulations we observed a remarkable T dependence which is analogous to the one realized in Figures 3 and 6. With the aid of Figure 8 we finally analyze the atomic localization parameters introduced in the foregoing section. The topmost full curve gives the overall atomic fluctuations in the harmonic oscillator approximation 0: according to eq 14. The dotted curve refers to D,””, Le., the thermal spatial uncertainty in the ho limit derived via eq 15. Once 0: and DP are known, the harmonic oscillator Dt can be calculated via eq 13. In the diagram Dt is symbolized by the broken curve. We point out that all curves have been calculated by adopting the experimental phonon DOS.39 At room temperature roughly comparable contributions to the net atomic fluctuations from

Dynamics of Carbon Nuclei in Cm lor

I

I

I

,

J. Phys. Chem., Vol. 99, No. 33, 1995 12407 1

1

,

,

C

T

0

200

400

600

800

Temperature (K) Figure 8. Delocalization parameters of the carbon nuclei of the Cm molecule as a function of temperature. The solid and dotted curves correspond to the parameters 0:”and DF in the harmonic oscillator

approximation; they have been derived from the experimental phonon spectrum. The quantum delocalization in the ho description DF is symbolized by the broken curve. The open circles are PI results for the quantum delocalization DY. The open squares have been derived from the latter PI curve by means of a projection onto the directions of the principal axes. The topmost data points correspond to the spatial fluctuations into the out-of-plane direction (Le,, radial breathing mode), while the remaining two sets of data points correspond to in-plane fluctuations (Le., pentagonal pinch modes). quantum and thermal degrees-of-freedom occur. With decreasing T the quantum fluctuations are enhanced at the expense of the thermal fluctuations. It is of some interest to mention that in ref 18 a ho approximation for Dt = 0:”has been discussed which is based on density functional calculations. The zeropoint amplitude of about 7 pm derived in that work is comparable to the present “experimental” ho value. Let us consider the PI-based Dp’curve in Figure 8 (open circles). The deviations between D! and DP are rather small in the T interval above 150 K. With decreasing T the DE;’distribution is shifted above the DP curve. All quantities of Figure 8 introduced up to now referred to the 3D vector r,. The overall Dg” fluctuations have been decomposed into three components associated with the principal axes of the c60 molecule. One component refers to the radial breathing mode (upper collection of squares in Figure 8), while two refer to in-plane pinch modes (open squares in the diagram). The difference between the two types of zero-point amplitudes is caused by different vibrational frequencies for the respective modes. The PI results demonstrate that the quantum delocalization of the C nuclei in Cm exceeds 7 pm even at low temperatures. The amplitude of the zero-point motions exceeds the static (=equilibrium) difference in the lengths of the short and long CC bonds which amounts to 4.4pm in this carbon system. The overall atomic fluctuations Dp” shown in Figure 8, and obviously also Or’, are roughly constant in the T interval between 0 and 200 K. In this T window almost constant dynamical properties of the C nuclei are conserved. The relative contributions of thermal and quantum degrees-of-freedom, however, depend on the temperature. It is important to emphasize that a classical treatment fails to reproduce the dynamics of the C atoms in the fullerenes/fullerides even at room temperature. It is our opinion that zero-point amplitudes as derived for the Cm molecule may be characteristic for other hydro/carbon

n systems as well. The precise knowledge of these parameters is of great importance for a correct interpretation of crystallographic data. On the basis of diffraction experiments, it is almost impossible to explain the microscopic origin of measured atomic fluctuations. They can be caused by vibrational motions or by static disorder. Predominance of the second effect leads to fluctuations that depend on the quality of the crystal while predominance of dynamic (i.e., vibrational) disorder is intrinsic. See ref 45 for a detailed discussion. In the past it has been assumed quite often that the two contributions to measured atomic fluctuations can be discriminated on the basis of T dependent diffraction data. In many crystals it had been found empirically that the atomic fluctuations caused by vibrational motions are a linear function of T within a rather large temperature inter~al.4~But in Cm we have the situation that the spatial fluctuations of the carbon nuclei do not exhibit any strong T dependence between T = 0 K and room temperature. So by experimental measurements as a function of T, it is impossible to discriminate between static and dynamic processes causing the net measurable atomic fluctuations. In the literature it is well known that the determination of the crystal structure of hydrocarbon n systems by X-ray diffraction may be a difficult task as a result of strong disorder effects. In this context we want to mention the problems encountered in X-ray work on tetra-t-butylcy~lobutadiene~~ and 1,3,5,7-tetra-t-butyl-s-inda~ e n e . 4We ~ believe that zero-point vibrations as discussed in the present work may be of great importance in these compounds. They lead to strong dynamic disorder and may be a cause for the complications in the determination of high-quality structural data. On the basis of the present theoretical PI results, it is understandable that atomic motions in Cm may have a decisive influence on many observables. The largest effects are expected in metallic configurations of fulleride solids, where the electronic DOS at the Fermi energy EF or the shape of the Fermi surface is sensitively influenced by zero-point vibrations. In a recent ~ t u d y ~we , ~have * performed detailed band structure calculations on solid c60 and potassium-doped fullerides as a function of characteristic intramolecular degrees-of-freedom. For a theoretical tool, we have employed a crystal orbital (CO) formalism based on an INDO (intermediate neglect of differential overlap) H a m i l t ~ n i a n .In ~ ~these CO calculations we have recognized strong modifications of the above quantities as a response to deviations from a chosen “equilibrium” geometry along the normal coordinates of the zero-point vibrations. The corresponding vibrations are non-negligible for certain electronic properties such as the Fermi energy or the Fermi momentum. In ref 9 we have suggested that direct carrier attraction may contribute to the superconducting pairing in fullerides. The electronic prerequisite causing this direct attraction is a strongly localized hole pocket around the r point (center of the Brillouin zone). The zero-point motions allow an enhancement of the electronic factors which support direct Cooper pair formation. The totally symmetric radial breathing mode, for example, has a pronounced influence on the magnitude of an attractive nearest-neighbor “bond charge hopping” integral t h = (ab/aa) with a and b representing Wannier-type orbitals on adjacent fulleride units. We have used the “chemical” (11/22)convention in the two-electron element; see ref 9. fh is attractive in M3Cm fullerides (M = alkali atom) due to the u symmetry of the conduction band. The sign and magnitude of t h determine whether or not direct carrier attraction at EF is accessible in fullerides. Both quantities of the fh integral are influenced by zero-point motions. We have mentioned this dependence to emphasize that the evaluation of microscopic electronic param-

12408 J. Phys. Chem., Vol. 99, No. 33, 1995 eters in superconducting fullerides, such as the electronic DOS, FS properties, or the el-ph coupling constants, on the basis of a single geometry may contain remarkable uncertainties.

IV. Conclusions We have performed Feynman path-integral simulations to analyze finite-temperature properties of Cm. Comparison of classical and quantum simulations has shown the importance of quantum effects on the dynamics of the atoms over a large T interval. Below 200 K only zero-point motions are operative. The light mass of the C atoms in fullerenes leads to amplitudes of the zero-point vibrations that exceed 7 pm. This width is of particular importance since it is of the same order of magnitude as the difference in the length between the two types of CC bonds in c60. The PI simulations have demonstrated that the distribution functions confined to the two types of bonds show sizable overlap. In some detail we have discussed the PI results of the radial distribution function and the integrated number of CC bonds per atom as a function of the radial vector r. Such an overlap in the rdf for long and short CC bonds is even conserved in cyclobutadiene with a difference of roughly 24 pm between both CC bonds.50 In the C60 system quantum fluctuations cause overlapping peaks in the angular distribution function. It is our opinion that the theoretical results collected in Figures 3-8 are of some utility to develop a more realistic picture of the molecular and electronic structure of (hydro)carbon n systems, a picture which goes beyond the simple static Kekuli-like representations widely used in chemistry. As we already have mentioned, the precise knowledge of the microscopic origin of the spatial atomic fluctuations in solids should be of some importance in the determination of crystal structures. In the present PI simulations we had to employ an empirical model potential with its typical limitations in accuracy. This restriction has been necessary as a result of the prohibitive computational expenditure to derive quantum paths for any C atom of the c 6 0 system, a procedure yielding several thousands of different nuclear configurations. Our approach is different in philosophy from the Car-Paninello density functional work of ref 18. In contrast to our empirically fitted V(R) potential the potential energy used in the cited work had its origin in density functional calculations. The trajectories generated in ref 18, however, are classical while the present setup is an “exact” many-body description. The optimum setup to evaluate the finite-temperature properties of nuclear degrees-of-freedom in a quantum description would be the combination of a CarPaninello method and a path-integral formalism to derive the dynamics of the system,

Acknowledgment. This work has been supported by CICYT (Spain) under Contract No. PB93-1254 and by a DFG and CSIC travel grant to visit the TH Darmstadt (R.R.). M.C.B. thanks the Fonds der Chemischen Industrie for financial support. We are grateful to R. Schlogl for useful comments and to C.P. Herrero for stimulating discussions. We thank S . Philipp for critically reading the manuscript. References and Notes (1) Kroto, H. W.; Allaf, A. W.; Balm, S. P. Chem. Rev. 1991,91, 1213. (2) Buckminsterfullerenes; Billups, W. E., Ciutolini, M. A., Eds.; VCH: Weinheim, Germany, 1993. (3) Solid State Physics; Ehrenreich, H., Spaepen, F., Eds.; Academic: Boston, 1994; Vol. 48. (4) Jerome, D.; Schultz, H. J. Adv. Phys. 1982, 31, 399. (5) Bulaevskii, L. N. Adv. Phys. 1988, 37, 433. (6) Erwin, S. C.; Pickett, W. E. Science 1991, 254, 842. (7) Huang, M.-Z.; Xu, Y.-N.; Freeman, A. J. Physica C 1992, 191, 399.

Bohm and Ram’rez (8) Ramirez, A. P.; Kortan, A. R.; Rosseinsky, M. J.; Duclos, S. J.; Mujsce, A. M.; Haddon, R. C.; Murphy, D. W.; Makhija, A. V.; Zahurak, S. M; Lyons, K. B. Phys. Rev. Lett. 1992, 68, 1058. (9) BBhm, M. C., Schulte, J. Submitted for publication. (10) Chakravarty, S.; Kivelson, S. Europhys. Lett. 1992, 16, 751. (11) Jansen, L.; Block, R.; Lombardi, E. Physica C 1991, 182, 17. Jansen, L.; Chandran, L.; Block, R. Chem. Phys. 1993, 176, 1. (12) Zhang, Q.-M.; Yi, J.-Y.; Bemholc, J. Phys. Rev. Lett. 1991, 66, 2633. (1 3) Wang, C. Z.; Chan, C. T.; Ho, K. M. Phys. Rev. B 1992,46,976 1. (14) Kim, E.; Lee, Y. H.; Lee, J. Y. Phys. Rev. B 1993.48, 18230. (15) Robertson, D. H.; Brenner, D. W.; White, C. T. J . Phys. Chem. 1992, 96, 6133. (16) Zhang, B. L.; Wang, C. Z.; Chan, C. T.; Ho, K. M. J . Phys. Chem. 1993, 97, 3134. (17) Kim, S. G.; Tomhek, D. Phys. Rev. Lett. 1994, 72, 2418. (18) Kohanoff, J.; Andreoni, W.; Parrinello, M. Phys. Rev. B 1992, 46, 437 1. (19) Car, R.; Paninello, M. Phys. Rev. Lett. 1985, 55, 2471. (20) Pickett, W. E. In Solid State Physics; Academic: Boston, 1994; Vol. 48, p 225. (21) Feynman, R. P. Rev. Mod. Phys. 1948, 20, 367. (22) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill: New York, 1965. Feynman, R. P. Statistical Mechanics; Benjamin: New York, 1972. (23) Ceperly, D. M.; Kalos, M. H. Monte Carlo Methods in Statistical Physics; Binder, K., Ed.; Springer: Berlin, 1981. (24) Kalos, M. H. Monte Carlo Methods in Quantum Problems; Reidel: Dordrecht, 1984. (25) Gillan, M. J. Computer Modelling of Fluids, Polymers and Solids; Catlow, C. R. A., Parker, S. C., Allen, M. P., Eds.; Kluwer: Dordrecht, 1990. Gillan, M. J. Philos. Magn. A 1988, 58, 257. (26) von der Linden, W. Phys. Rep. 1992, 220, 53. (27) Rainirez, R.; Herrero, C. P. Phys. Rev. B 1993, 48, 14659. (28) Ramirez, R.; Herrero, C. P. Phys. Rev. Lett. 1994, 73, 126. (29) Bohm, M. C.; Schulte, J.; Utrera, L. Mol. Phys. 1993, 79, 1239. Bohm, M. C.; Schulte, J.; Philipp, S. Chem. Phys. 1993, 176, 109. (30) Ramirez, R.; Bohm, M. C. J . Phys. C, in press. (31) Coulson, C. A.; Rushbrooke, G. S. Proc. Cambridge Philos. SOC. 1940, 36, 193. (32) Axe, D.; MOSS,S. C.; Neumann, A. In Solid State Physics; Academic: Boston, 1994; Vol. 48, p 149. (33) Tersoff, J. Phys. Rev. B. 1988, 37, 6991; Phys. Rev. Lett. 1988, 61, 2879. (34) Brenner, D. W. Phys. Rev. B 1990, 42, 9458. (35) Chelikowsky, J. R. Phys. Rev. B 1992, 45, 12062. (36) Hall, R. W. J . Chem. Phys. 1989, 91, 1926. Ceperly, D. M. Phys. Rev. Lett. 1992, 69, 331. (37) Schiitt, J.; Schulte, J.; Bohm, M. C.; Soos, Z. G. Mol. Phys. 1995, 84, 1127. (38) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J . Chem. Phys. 1953, 21, 20. (39) Cappeletti, L.; Copley, R. D.; Kamitakahara, W. A,; Li, F.; Lannin, J. S.; Ramage, D. Phys. Rev. Lett. 1991, 66, 3261. (40) Wu, 2.C.; Jelski, D. A.; George, T. F. Chem. Phys. Lett. 1987, 137, 291. Cyvin, S. J.; Brensdal, E.; Cyvin, B. N.; Brunvoll, J. Chem. Phys. Lett. 1988, 143, 377. (41) Bethune, D. S.; Meijer, G.; Tang, W. C.; Rosen, H. J.; Golden, W. G.; Seki, H.; Brown, C. A.; de Vries, M. S. Chem. Phys. Lett. 1991, 179, 181. (42) Negri, F.; Orlandi, G.; Zerbetto, F. Chem. Phys. Lett. 1992, 190, 174. (43) Li, F.; Ramage, D.; Lannin, J. S.; Conceicao, J. Phys. Rev. B 1991, 44, 13167. (44) Bemholc, J.; Yi, J.-Y.; Zhang, Q.-M.; Brabec, C. J.; Anderson, E. B.; Davidson, B. N.; Kajihara, S. A. Z. Phys. D.1993, 26, 74. (45) Dunitz, J. D.; Maverick, E. F.; Trueblood, N. Angew. Chem. 1988, 100, 910. (46) Imgartinger, H.; Riegler, N.; Malsch, K.-D.; Schneider, K.-A,; Maier, G. Angew. Chem. 1980, 92, 214. Imgartinger, N.; Nixdorf, M. Angew. Chem. 1983, 95, 415. (47) Dunitz, J. D.; Kriiger, H.; Imgartinger, H.; Maverick, E. F.; Wang, Y.; Nixdorf, M. Angew. Chem. 1988, 100, 415. (48) Bohm, M. C.; Schulte, J. In preparation. (49) Ramirez, R.; Bohm, M. C. Int. J . Quantum Chem. 1988, 34, 47. (50) Facelli, J. C. THEOCHEM 1991, 82, 119. JP9502666