QM and MD Studies on Formation Mechanisms of Fullerenes

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

^{a}*Materials and Process Simulation Center,
California Institute of Technology,
Pasadena, California, 91125*

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
C_{60} 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 (C_{60}, C_{70}, 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 (C_{2}, C_{3})?
Is the C_{60} fullerene formed by adding C_{1}, C_{2}, or
C_{3} to some
smaller fullerene or is C_{60} formed by isomerization of some type of
precursor molecule C_{60}?
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 C_{60} 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.

- {For n=4m, the minimum energy structure has a
polyacytelene geometry of alternating single and triple bonds.
The bond length difference is from 0.5 A to $0.9 A.
We find that inclusion of correlation reduces the dimerization amplitude,
simular to the case in polyacetylene [Konig 1990].
Comparing to the DFT geometry, Hartree-Fock (HF)gives too large of a bond
alternation, along with too large of angle alternations.
Our HF calculation gives a bond difference of $0.16 A, in
agreement with that of Feyereisen {\it et al}.[Feyereisen 1992]
As for angle alternation, for C
_{20}, HF gives 160^{o}-164^{o}[Raghavachari 1993] while DFT gives 161.5^{o}-162.5^{o} - For n=4m+2, the minimum energy structure has the polyallene geometry with equal bond lengths. This is due to the resonance between the two structures, involving the pi-bond perpendicular to the plane and pi-bond parallel to the plane.

C_{4m+2} is more stable than C_{4m}. But as n --> infinity,
the difference in E_{coh} decreases to zero,
leading to E_{infinity}^{sp1}=6.56eV.

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

Indeed we found strain energy increase linearly with 1/n^{2} with slopes of
63.3 eV/n^{2} for 4m and 40.1eV/n^{2} for the 4m+2, respectively.
Both converg to E_{infinity}^{sp1}=6.56eV$

The force field took the following form:

Here q_{r1}(l)=R_{i}(l)-R_{i0}(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
q_{theta}(l)=theta(l)-theta_{0}(l)
We use the periodic boundary condition so that, with theta_{0} = 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.
E_{0} is a reference energy corresponding to zero strain energy structure.
(infinite linear chain)
Comparing E_{MSX} with E_{DFT} for several structures,
we can derive the force field parameters.

In a simular fasion we can derive the force field for the sp^{2}
bonded carbons.
The optimum structure for bulk carbon is graphite, which has each carbon
bonded to three others (sp^{2} 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 C_{n} 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 sp^{2} carbon of
E_{coh}^{sp2} = 7.71 eV.
This can be compared to the experimental cohesive energy of a single graphitic
sheet of E^{sheet}_{coh}=7.74eV.
This is derived from the experimental cohesive
energy [CRC Handbook] of graphite of E^{graphite}=7.8eV plus total
Van der Waals attraction of E^{vdw}=0.056eV between sheets calculated
using the graphite force field. [Guo 1992]

Now that we have the energy and force field of both sp^{1} and
sp^{2}
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.

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

- The reaction from two C
_{30}ring (**I**) to C_{60}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. - The energy of
**I**is determined directly from Figure 1. - The energy difference between
**I**and**III**is a strain energy which can be calculated using the MSX FF. - 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 C_{6}H_{2}molecules to form to the 4-membered ring, C_{12}H_{2}. (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**. - Thus the energy of C
_{60}bycyclic ring (**II**) can be calculated by**II--III--I**.

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:

We have chosen E_{0}=60 epsilon_{1},
as zero point. Here, n_{2} is the sp^{2} bonded carbons,
n_{2} (epsilon_{1} - epsilon_{2}) gives the energy gained by
converting sp^{1} bonded carbon into spn^{2} bonded carbon, with
epsilon_{1}=-6.56eV$ and epsilon_{2}=-7.71 eV.
d_{1} is the energy of a dangling bond relative to the sigma-bonded state,
n_{R} number of of such dangling bond(radicals);
d_{2} is the energy of an atom participating bended planar
pi-bond relative to the sigma-bonded state and
n_{sigma pi} is the number of such atoms.
We use the Benson-like scheme to evaluate d_{1} and d_{2} [Guo 1992]
and found d_{1}=2.32eV and d_{2}=1.64eV.
E^{str} (n_{2}) 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.

At the beginning atomic carbons combine themselves to form dimers and trimer,
C_{2}, C_{3}. These would then grow into linear chain of carbons
C_{n}, 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 C_{60} formation, as suggested by
Jarrolds experiments, [Bowers 1993a, 1993b, Jarrold 1994]
is to combine two C_{30} rings to form a bycyclic C_{60} ring, which
in turn isomerized into a C_{60} 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 C_{60} 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 E_{o}=60 epsilon_{1},
where epsilon_{1}=-6.56eV.

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

- (i)
**1'**={30+4+30} --->**2**={30+4+6+30}. This is a Bergman diyne cyclization which forms a 6-membered 1,4 benzyne-like ring from two triple bonds. This leads to two isolated radical sites (sp^{2}-like orbtials in the plane, that cannot form a bond), and we find that this increases the energy by about 0.7eV. - (ii)
**2**={30+4+6+30} --->**3**={30+8+30}. This process kills two dangling bond by breaking one sigma bond and forming two pi bonds. This process is down hill by about 1.3eV. - (iii)
**3**={30+8+30} --->**4**={30+8+6+22}. This involves breaking an in-plane $\pi$ bond and forming a $\sigma$ bond. In the process there is bending of one triple bond to from a 1,2-benzyne-like ring which includes a new radical site. This process is uphill by 1.66eV. Jarrold postulated**4'**which is 2.1eV above the bycyclic rings from our calculation. - (iv)
**4**={30+8+6+22} --->**5**={6+6+55}. This involves twisting open the original 4-membered ring. Then it is followed by relaxing the 50 carbon chain to reduce the strain energy. This {6+6+55} contains two dangling bonds. This process is downhill by about 0.67eV. - (v)Spiral growth around the {6+6}. As a first step
**5**={6+6+55} --->**6**={6+6+53+5}. This uses one of the sp^{2}orbitals of the 1,2-benzyne-like ring to attack a triple bond and form a new 5-membered ring. This process is down hill by 0.13eV. - (vi) Continue the spiral growth to form C
_{60}fullerene. The energies calculated using the extended MSX FF on these systems are shown in Figure 5 where we see that they are monotonically downhill. The overall gain of energy from {6+6+53+5} to C_{60}is about 30eV, so that no barriers are expected to impede these steps. Figure 6 illustrates some of the intermediates between**5**and the fullerene.Figure 5. Energies calculated by MSXX FF. The driving force for the growth is the gain in forming sp

^{2}sigma bond. The opposing forces are the energy lost by the radicals created along the way and the increasing strain energies.The Jarrold mechanism represents an innovative major step forward in understanding the formation of C

_{60}fullerene. Our energetic analysis shows that some of the reactions pathways have large energy barriers, however they never exceed the energy available to the unimolecular reaction.The similar approach could be used to study the formation of other fullerenes, like C

_{40}, C_{50}, C_{70}, etc..Figure 6. Intermediates between the 5 and the fullerene. ### Summary

Why C

_{60}is so stable and how C_{60}fullerenes are formed, these are the two most interesting problems in basic fullerene research. We have studied the formation mechanism of C_{60}fullerenes using first principles calculation and molecular dynamics simulations. We have derived a force field (MSX FF) that is suitable to describe both the sp^{1}hybrid and sp^{2}hybrid carbons. Combining DFT and MD with MSX FF we found the relative thermal stability of various neutral isomers at each cluster size n and predicted the relative abundancy of these neutral species for thermal equilibrium. We identified a complete path to form a C_{60}fullerene from atomic carbons and calculated its energetics. Our approach is fully applicable to other possible reaction paths and other fullerenes.## References

- Joanna M Hunter, James L Fye, Eric j. Roskamp, and Martin F.Jarrold
*J. Phys. Chem.***98**, 1810-1818 (1994); - Joanna M Hunter, James L Fye, and Martin F.Jarrold
*J. Chem. Phys.***99**, 1785-1795 (1993) - Gert Von Helden, Nigel G. Gotts and Michael T. Bowers,
*Nature***363**, 60 (1993); Gert von Helden, Ming-Teh Hsu, Nigel Gotts, and Michael T. Bowers,*J. Phys. Chem***97**8182-8192, (1993) - Benny G. Johnson, Peter M. BW. Gill, and John A. Pople,
*J. Chem. Phys.***98**5612 (1993). - Murco N. Ringnalda, Jean-Marc Langlois, Burnham H. Greeley, Robert B Murphy, Thomas V. Russo, Christian Cortis, Richard P. Muller, Bryan Marten, Robert E. Donnelly, Jr., Daniel T. Mainz, Julie R. Wright, W. Thomas Pollard, Yixiang Cao, Youngdo Won, Gregory H. Miller, William A. Goddard III, and Richard A. Friesner,
*PS-GVB v2.2*, Schrodinger, Inc., (1995) - G. Konig, G. Stollhoff
*Phys. Rev. Lett.***65**, 1239 (1990) - Martin Feyereisen, Maciej Gutowski, and Jack Simons, and Jan Alm,
*J. Chem. Phys.***96**, 2926 (1992) - K. Raghavachari
*el al., Chem. Phys. Lett.***214**, 357 (1993). -
*CRC Handbook* - Yuejing Guo,
*Ph.D. Thesis*, California Institute of Technology. (1992) - Jurg Hutter, Hans Peter Luthi, and Francois Diederich
*J. Am. Chem. Soc.***116**750-756 (1994) - S. C. O'Brien, J. R. Heath, H. W. Kroto, R. F. Curl, and R. E. Smalley,
*Chem. Phys. Lett.***132**99-102 (1986)

- Joanna M Hunter, James L Fye, Eric j. Roskamp, and Martin F.Jarrold