QM and MD Studies on Formation Mechanisms of Fullerenes

Xinlei Hua, Tahir Cagin, Jianwei Che and William A. Goddard III

aMaterials and Process Simulation Center, California Institute of Technology,
Pasadena, California, 91125

E-mail: tahir@wag.caltech.edu

This is a draft paper for the
First ElBA Foresight European Conference on Molecular Nanotechnology.


One of the most puzzling aspects of fullerenes is how such complicated symmetric molecules are formed from a gas of atomic carbons, namely, the atomistic or chemical mechanisms. Are the atoms added one by one or as molecules (C2, C3)? Is there a critical nucleus beyond which formation proceeds at gas kinetic rates? What determines the balance between forming buckyballs, buckytubes, graphite and soot? The answer to these questions is extremely important in manipulating the systems to achieve particular products. A difficulty in current experiments is that the products can only be detected on time scales of microseconds long after many of the important formation steps have been completed. Consequently, it is necessary to use simulations, quantum mechanics and molecular dynamics, to determine these initial states. Experiments serve to provide the boundary conditions that severely limit the Using quantum mechanical methods (DFT) we derived a force field (MSXX FF) to describe one-dimensional (rings) and two-dimensional (fullerene) carbon molecules. Combining DFT with the MSXX FF, we calculated the energetics for the ring fusion spiral zipper (RFSZ) mechanism for formation of C60 fullerenes. Our results shows that the RFSZ mechanism is consistent with the quantum mechanics (with a slight modification for some of the intermediates).


One of the most puzzling aspects of fullerenes (C60, C70, etc.) is how such complicated symmetric molecules are formed from a gas of atomic carbons [atomcarbon], namely, the atomistic or chemical mechanism. Are the atoms added one by one or as molecules (C2, C3)? Is the C60 fullerene formed by adding C1, C2, or C3 to some smaller fullerene or is C60 formed by isomerization of some type of precursor molecule C60? Is there a critical nucleus beyond which formation proceeds at gas kinetic reates? What determines the balance bwtween forming buckyballs, buckytubes, graphite, and soot? The answer to these questions might lead to means of manipulating the systems to achieve particular products.

A difficulty in current experiments [Jarrold 1993, Jarrold 1994, Bowers 1993a, 1993b] is that the products can only be detected on times scales of $\mu$s, long after many of the important formation steps have been completed. Consequently, it is necessary to use first principles quantum mechanical theory to determine these initial states; however, the experiments serve to provide boundary condition that severely limint the possibilities, making the use of first principles theory practicable.


We selected density functinal theory (DFT) as the best compromise between accuracy and speed for studying these systems. We use the Becke gradient corrected exchange and the gradient corrected correlation functional of Lee, Yang, and Parr. [Johnson 1993] The calculations were carried out using PS-GVB with the 6-31G* basis set [Rignalda 1995].

Because the carbon rings play a central role, we studied how the structures and energetics of such rings changed with size and extracted a force field (denoted as the MSX FF) that would reproduce the DFT energetics and structures. This MSX FF would be used later in conjunction with the DFT calculations on various multiring systems to estimate the energetics of the full 60 atom systems without the necessity of DFT on the complete system.

The calculations on ring systems up to C60 are shown in Figure 1. The energies quoted here are cohesive energy par carbon atom. In calculating these energies we used as our reference the triplet C atom, calculated by LSDA.

Figure 1. Cohesive energy per carbon atom.
We found that

C4m+2 is more stable than C4m. But as n --> infinity, the difference in Ecoh decreases to zero, leading to Einfinitysp1=6.56eV.

Both the polyacetylene and polyallene structures involve sigma-bonds that are sp1 hybrids, which prefer linear geometries. Thus we expect a strain energy proportional to

Indeed we found strain energy increase linearly with 1/n2 with slopes of 63.3 eV/n2 for 4m and 40.1eV/n2 for the 4m+2, respectively. Both converg to Einfinitysp1=6.56eV$

The force field took the following form:

Here qr1(l)=Ri(l)-Ri0(l)$ is the bond strain term, where for n=4m, i=1 is the triple bond and i=2 is the single bonds; for n=4m+2 their bonds are equivalent. The angle strain term is qtheta(l)=theta(l)-theta0(l) We use the periodic boundary condition so that, with theta0 = pi, n/2+1=1, where n is the total number of atoms in the system and n/2 is the number of unit cells. E0 is a reference energy corresponding to zero strain energy structure. (infinite linear chain) Comparing EMSX with EDFT for several structures, we can derive the force field parameters.

In a simular fasion we can derive the force field for the sp2 bonded carbons. The optimum structure for bulk carbon is graphite, which has each carbon bonded to three others (sp2 bonding) to form hexagonal sheets stacked on each other. The fullerenes structures can be considered as finite two dimensional anologs, in which each carbon is distorted (strained) from its preferred planar configuration. Since the strain should be proportinal to the square of the planar distortion angle, delta psi, we expect that the strain energy should scale as 1/n

We have performed the DFT(Becke/LYP) calculations on Cn fullerenes with n=20, 32 and 60. Figure 1 shows the cohesive energies per carbon atom.

Extrapolating the calculated cohesive energy to n ---> infinity leads to a cohesive energy per sp2 carbon of Ecohsp2 = 7.71 eV. This can be compared to the experimental cohesive energy of a single graphitic sheet of Esheetcoh=7.74eV. This is derived from the experimental cohesive energy [CRC Handbook] of graphite of Egraphite=7.8eV plus total Van der Waals attraction of Evdw=0.056eV between sheets calculated using the graphite force field. [Guo 1992]

Now that we have the energy and force field of both sp1 and sp2 hybridized carbon we can get the energetics of any carbon clusters. Adding the entropic controbution within harmonic approximation using FF, we get the free energy of various species at different temperature, which dictate the thermal equilibrium distribution of these species. Our population analysis is shown in Figure 2.

Figure 2. Population analysis of species.

For studying formation reaction sequence we adopt two level of models, a fine one and a coarse one, as explained below.

(A)Fine model:

The energies of a structure were computed via the following procedures that combine the DFT with MD as illustrated in Figure 3.

Figure 3. Paths.
  1. The reaction from two C30 ring (I) to C60 bycyclic ring molecules (II) is achieved via an intermediate III. For each of I, II and III, the system is partitioned into two parts: part A involving bond lengths changes, and part B involving continium deformation.
  2. The energy of I is determined directly from Figure 1.
  3. The energy difference between I and III is a strain energy which can be calculated using the MSX FF.
  4. The energy difference between III and II is calculated in two parts. (a)part A: DFT calculations are carried out on the reaction of two C6H2 molecules to form to the 4-membered ring, C12H2. (b)part B: we calculate the corresponding change in going from III to II using MSX FF. Then we combine A and B to get the energy difference between III and II.
  5. Thus the energy of C60 bycyclic ring (II) can be calculated by II--III--I.

(B)Coarse model:

We extend the MSX FF to include terms capable of describing the different bonding schemes. The key components are the additive energy terms for the dangling bond and the energy cost for bending a triple bond to form a 1,2-benzyne. Our FF are defined as follows:

Etot n2 = Ebond + Eradical + Estrain = n2 (epsilon1 - epsilon2) + d1 nR + d2 nsigma pi + Estr (n2)

We have chosen E0=60 epsilon1, as zero point. Here, n2 is the sp2 bonded carbons, n2 (epsilon1 - epsilon2) gives the energy gained by converting sp1 bonded carbon into spn2 bonded carbon, with epsilon1=-6.56eV$ and epsilon2=-7.71 eV. d1 is the energy of a dangling bond relative to the sigma-bonded state, nR number of of such dangling bond(radicals); d2 is the energy of an atom participating bended planar pi-bond relative to the sigma-bonded state and nsigma pi is the number of such atoms. We use the Benson-like scheme to evaluate d1 and d2 [Guo 1992] and found d1=2.32eV and d2=1.64eV. Estr (n2) is the strain energy and it is evaluated at the minimum energy structure.

We would use the fine model for the initial steps in the $C_{60}$ formation. As the reaction take off and begin to release more and more energy, we switch to the coase one.

The spiral model of fullerene formation

At the beginning atomic carbons combine themselves to form dimers and trimer, C2, C3. These would then grow into linear chain of carbons Cn, etc., for n<10. [Hutter 1994] When n>10 the carbon cluster prefer ring structure [Hutter 1994] because beyond n>10 the energy gain in killing the dangling bonds at the two ends over compensate for the strain energy incured by folding up the chain. At around n>30 the ring structures give way to fullerene structures,[Bowers 1993a, 1993b] because replacing more pi-bond by sigma-bond over compensate for the strain of folding the 2-D net.

One process of C60 formation, as suggested by Jarrolds experiments, [Bowers 1993a, 1993b, Jarrold 1994] is to combine two C30 rings to form a bycyclic C60 ring, which in turn isomerized into a C60 fullerene. This unimolecular reaction will be the focus of our study.

As a mnemonic for referring to the various structures, we will simply denote the ring sizes of a structure. Thus the simple C60 ring is denoted as 60, while the double ring system, 1, is 30+4+30. This notation does not uniquely describe a structure, but it is for the species we will consider. We will take the reference energy to be Eo=60 epsilon1, where epsilon1=-6.56eV.

Following Jarrold, the first few steps in the reaction are as follows: (see Figure 4)

Figure 4. First few steps of reactins.