Analytical carbon-oxygen reactive potential
A. Kutana and K. P. Giapis
a
Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena,
California 91125, USA
Received 22 January 2008; accepted 9 May 2008; published online 19 June 2008
We present a reactive empirical potential with environment-dependent bond strengths for the
carbon-oxygen
CO
system. The distinct feature of the potential is the use of three adjustable
parameters characterizing the bond: the strength, length, and force constant, rather than a single
bond order parameter, as often employed in these types of potentials. The values of the parameters
are calculated by fitting results obtained from density functional theory. The potential is tested in a
simulation of oxidative unzipping of graphene sheets and carbon nanotubes. Previous higher-level
theoretical predictions of graphene unzipping by adsorbed oxygen atoms are confirmed. Moreover,
nanotubes with externally placed oxygen atoms are found to unzip much faster than flat graphene
sheets. ©
2008 American Institute of Physics
.
DOI:
10.1063/1.2940329
I. INTRODUCTION
Empirical interatomic potentials offer substantial advan-
tages in calculating the properties of large statistical entities
as compared to more rigorous
ab initio
or tight binding
methods, which typically produce more accurate results al-
beit for much smaller systems. In particular, reactive empiri-
cal bond order
REBO
potentials have been very successful
in simulating covalently bonded materials such as Si
Ref.
1
and hydrocarbons.
2
,
3
They have also facilitated the descrip-
tion of chemical reactions and barriers in large systems
4
and
established molecular dynamics simulations as an important
tool in designing new materials, estimating their properties,
and understanding materials processing. This success moti-
vated the development of REBO potentials for chemically
complex systems of different atoms.
5
–
7
The REBO approach is based on a fixed relationship
between bond energy and bond order, which is a function of
the coordination numbers of atom pairs forming a bond. This
empirical relationship is only approximately satisfied in real
molecules and may yield imprecise results especially in sys-
tems encompassing two or more chemically distinct atoms.
Each new atom considered introduces a new set of param-
eters, which must be optimized by fitting the equilibrium
distance and energy of bonds formed between the new atom
and the previous set of atoms. This fitting strategy is not
specific to the REBO potential. It can be implemented to
other potentials with their own sets of parameters, provided
they are complex enough to capture an increasing number of
bonding configurations.
We present here a new approach to generating empirical
reactive potentials for covalent systems. We focus on the
carbon-oxygen system as a model. Carbon-oxygen chemical
reactions are ubiquitous, from combustion to metabolic reac-
tions in cells.
8
The importance of carbon-oxygen chemistry
has led to many electronic structure and potential energy
surface calculations for oxygen-hydrocarbon interactions
using
ab initio
methods.
9
–
12
Yet, fast analytical interatomic
potentials capable of capturing the chemical interactions
between these elements have been scarce. In fact, there exists
currently only one such potential
6
to our best knowledge. In
order to fill this gap, a carbon-oxygen potential based on an
extensive fitting database of
ab initio
calculations is pre-
sented here. The motivation for this work has also been an
attempt to avoid any predefined fixed relationships between
bond energy, length, and force constants present in
traditional bond order potentials.
In Sec. II, we describe the formalism used to represent
the binding energy as a function of the environment of a
chemical bond. In Sec. III, the parameters of the formulas
given in Sec. II are fitted to results of
ab initio
calculations.
In Sec. IV, the C
u
O potential developed is applied to the
problem of unzipping of graphene sheets and carbon
nanotubes by adsorbed oxygen. Conclusions are given in
Sec. V.
II. REACTIVE POTENTIAL
One of the most successful approaches to modeling re-
active chemistry in atomistic simulations has been the bond
order method.
1
It introduces environment-dependent correc-
tions to the strengths of chemical bonds based on the local
coordination of atoms. Unlike multibody potential expan-
sions, where each successive term depends on one extra
atomic coordinate, the bond order formalism retains the form
of a pair potential; instead, all the corrections are included
into the coefficients of the potential function. Previous reac-
tive potentials for silicon
1
and hydrocarbons
C+H system
2
used a bond order parameter
b
ij
b
ij
1
to account for the
environment-dependent bond order between atoms
i
and
j
as
follows:
a
Author to whom correspondence should be addressed. Electronic mail:
giapis@cheme.caltech.edu.
THE JOURNAL OF CHEMICAL PHYSICS
128
, 234706
2008
0021-9606/2008/128
23
/234706/8/$23.00
© 2008 American Institute of Physics
128
, 234706-1
Downloaded 20 Jun 2008 to 131.215.225.9. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp
U
ij
=
f
c
r
ij
V
R
r
ij
−
b
ij
V
A
r
ij
.
1
Here,
r
ij
is the distance between atoms
i
and
j
,
f
c
r
ij
is the
cutoff function used for smooth potential truncation, and
V
R
r
ij
and
V
A
r
ij
are the repulsive and attractive
interactions usually given by Morse-type expressions
V
R
r
=
A
exp
−
1
r
,
2
V
A
r
=
B
exp
−
2
r
,
where
A
,
B
,
1
, and
2
are all positive
1
=2
2
=
in the
original Morse potential
. The total bonding energy is
obtained by summing
U
ij
over all bonds as follows:
E
bond
=
j
i
U
ij
.
3
In a more recent REBO potential for hydrocarbons,
3
a
repulsive term of the form
V
R
r
=
1+
Q
/
r
A
exp
−
1
r
4
was introduced, where
Q
is a fitting parameter. The modified
pre-exponential factor 1+
Q
/
r
accounts for Coulombic repul-
sion, thus rectifying the unphysical behavior of the Morse
potential at small
r
.
Morse-type functions are well suited for capturing the
universal energetics of covalent and metallic bonds.
13
,
14
The
bond order parameter
b
ij
depends on the atomic coordination
in a molecule or solid and is central to the potential’s ability
to reproduce correct binding energies. The form of Eqs.
1
and
2
with
1
=2
2
=
requires that bond energies and
bond distances follow the Pauling rule,
15
stating that the
logarithm of the bond order
bond strength
U
0
is a linear
function of the equilibrium bond distance
r
e
as follows:
log
U
0
= log
A
−2
r
e
.
5
Figure
1
shows the energies of CO
n
n
=1,4
molecules calcu-
lated using density functional theory
DFT
. The Pauling rule
is followed only approximately: the CO
3
and CO
4
molecules
are more stable than what would follow from a linear fit for
the CO and CO
2
molecules. Usually, a best fit to the Pauling
rule is performed in order to minimize the errors over differ-
ent coordinations. The fitting could also be improved by in-
troducing corrections to the angular energy of the molecule,
at the cost of a more complex expression for it. As an alter-
native to these approaches, we remove the requirement that
the approximate Pauling rule be satisfied, and allow all three
parameters
A
,
B
, and
in Eqs.
1
and
2
to vary indepen-
dently. Using the parameters
U
0
=
B
2
/
4
A
and
r
e
=
1
/
ln
2
A
/
B
, which are the poten-
tial well depth and equilibrium interatomic distance, respec-
tively, instead of
A
and
B
in Eq.
2
, the expression for the
bond energy becomes
U
ij
=
f
c
r
ij
U
0
exp
−2
r
ij
−
r
e
− 2 exp
−
r
ij
−
r
e
,
6
where
U
0
,
r
e
, and
are functions of the numbers of neigh-
bors on both sides of the bond between atoms
i
and
j
as
follows:
U
0
=
U
0
N
C
i
,
N
O
i
,
N
C
j
,
N
O
j
,
r
e
=
r
e
N
C
i
,
N
O
i
,
N
C
j
,
N
O
j
,
7
=
N
C
i
,
N
O
i
,
N
C
j
,
N
O
j
.
Here,
N
C
i
is the number of carbon atoms connected to atom
i
of the
ij
bond,
N
O
i
is the number of oxygen atoms con-
nected to this atom, while
N
C
j
and
N
O
j
are similar quantities
for atom
j
. Each parameter is a function of 2
K
variables,
where
K
is the total number of the atomic species in the
system. The coefficients
U
0
,
r
e
, and
are calculated for
stable molecules where bonds are assigned integers
N
C,O
i
and
N
C,O
j
. Use of fractional
N
values permits smooth switching
through different coordinations and allows chemical reac-
tions to occur.
The particular form of the smooth interpolating
functions used here for fitting the parameters in Eq.
7
as
functions of the number of neighbors is a multivariate
polynomial of
m
variables of the order
n
as follows:
P
n
x
1
,
x
2
, ...,
x
m
=
i
=0
n
x
1
i
j
=0
n
−
i
x
2
j
k
=0
n
−
i
−
j
x
3
k
̄
l
=0
n
−
i
−
j
−
k
−...
a
ijk
...
l
x
m
l
,
8
where
x
i
are the variables corresponding to
N
C
i
,...,
N
O
j
in
Eq.
7
, and
a
ijk
̄
l
are the coefficients of the polynomial. The
FIG. 1. The Pauling plot for the potential well depth and equilibrium bond
length of the CO bond represented by a Morse potential with variable
coefficients as a function of the number of the O atoms connected to a C
atom. The molecular energies are obtained using
GAMESS
DFT, 6-311G
*
,
and B3LYP density functional
.
234706-2 A. Kutana and K. P. Giapis
J. Chem. Phys.
128
, 234706
2008
Downloaded 20 Jun 2008 to 131.215.225.9. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp
number of polynomial variables
m
is equal to twice the num-
ber of the present reactive species, which is 2 for the current
CO system, hence,
m
=4. The degree
n
of the polynomial is
n
=3. The 35 unknown coefficients
a
ijkl
in
P
3
are found by
solving a linear system of equations as follows:
P
3
q
=
i
=0
3
x
1
i
q
j
=0
n
−
i
x
2
j
q
k
=0
n
−
i
−
j
x
3
k
q
l
=0
n
−
i
−
j
−
k
a
ijkl
x
4
l
q
,
9
where
q
=1,2,...,35
enumerates points in the four-
dimensional fitting space. The numbers of neighbors
x
1
,...,
x
4
on each side of the bond are calculated according
to
x
1
N
C
i
=
k
C near
i
f
c
r
ik
1−
f
c
r
jk
,
x
2
N
O
i
=
k
O near
i
f
c
r
ik
1−
f
c
r
jk
,
10
x
3
N
C
j
=
k
C near
j
f
c
r
jk
1−
f
c
r
ik
,
x
4
N
O
j
=
k
O near
j
f
c
r
jk
1−
f
c
r
ik
,
where
f
c
r
ik
is the cutoff function
f
c
r
=
1,
r
D
min
1 + cos
r
−
D
min
/
D
max
−
D
min
/
2,
D
min
r
D
max
0,
r
D
max
,
11
for which the appropriate parameters
D
max
and
D
min
are
taken for an atomic pair
ik
. The numerical values of
D
max
and
D
min
for different atomic pairs are listed in Table
I
. Note
that according to Eq.
10
, an atom, which is connected to
both atoms constituting the bond, is not considered to be a
neighbor of either of these atoms.
In addition to the bonding energy, the full form of the
total potential energy also includes the angular, dihedral, and
van der Walls
vdW
components as follows:
E
total
=
E
bond
+
E
angle
+
E
tor
+
E
vdW
=
j
i
U
ij
+
j
,
i
,
k
E
jik
angle
+
k
,
i
,
j
,
l
E
kijl
tor
+
E
vdW
.
12
The three-atom angular energy
E
jik
angle
is presented by
E
jik
angle
jik
=
f
c
r
ij
f
c
r
ik
1
2
k
angle
jik
−
0
2
,
13
where
jik
is the angle formed by the adjacent bonds
ij
and
ik
, and
0
is the equilibrium angle. The cutoff function
factors ensure smooth truncation of the angular term upon
dissociation of either of the bonds in the angle. The dihedral,
or torsional energy
E
kijl
tor
is described as
E
kijl
tor
kijl
=
f
c
r
ij
f
c
r
ik
f
c
r
jl
1
2
k
tor
1 − cos
2
kijl
,
14
where
kijl
is the torsional angle between the two planes
formed by triplets of atoms
kij
and
ijl
centered on the
ij
bond. Again, the cutoff functions are added for eliminating
the dihedral term when one or more of the bonds constituting
the angle dissolve. In the current version of the potential, the
spring constants and equilibrium angle in Eqs.
13
and
14
are insensitive to the local environment. Since the equilib-
rium angle is fixed at 180° in most cases, the angular com-
ponents in nonlinear molecules at equilibrium may not be
zero and will contribute to the total energy of these mol-
ecules. As a result, the bonding and angular energies are
coupled via the total energy. The coupling mandates that
bonding and angular terms be fitted simultaneously to yield
the correct total energy of the molecule.
At very short interatomic distances, the Morse-type
potential
1
approaches a finite value, while in reality the
interaction follows the Coulombic repulsion with 1
/
r
scaling. For a more realistic description of the interaction at
short
r
, a screened Coulomb potential function has been
introduced into the total potential function in that region. A
two-body Moliere potential
16
U
M
r
=
A
/
r
0.35 exp
− 0.3
r
/
a
s
+ 0.55 exp
− 1.2
r
/
a
s
+ 0.1 exp
− 6.0
r
/
a
s
15
represents the interaction between a pair of atoms as given
by the Thomas–Fermi model. Here,
A
/
r
is the Coulomb
TABLE I. Values of the parameters
D
min
and
D
max
used in the cutoff
function.
D
min
Å
D
max
Å
C
u
C
1.7
2.0
C
u
O
1.7
2.0
O
u
O
1.5
1.9
TABLE II. Parameters of the angular function for the CO potential.
Angle
0
deg
k
ang
eV
O
u
C
u
O
180.0
4.56
O
u
C
u
C
180.0
1.85
C
u
O
u
C
180.0
1.89
C
u
C
u
C
180.0
3.04
O
u
O
u
O
118.2
14.0
O
u
O
u
C
180.0
1.37
234706-3 Analytical carbon-oxygen reactive potential
J. Chem. Phys.
128
, 234706
2008
Downloaded 20 Jun 2008 to 131.215.225.9. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp
repulsion between a pair of bare nuclei, while
a
s
is the
Firsov
17
screening length
a
s
= 0.885 34
a
Bohr
Z
1
1
/
2
+
Z
2
1
/
2
−2
/
3
,
16
which is expressed as a function of the atomic numbers
Z
1
and
Z
2
of the two atoms, and the Bohr distance
a
Bohr
=0.529 Å. The total potential energy can then be written as
U
tot
r
=
f
c
r
U
M
r
+
1−
f
c
r
U
r
,
17
where
U
r
stands for the reactive many-body part of the
potential described by Eq.
6
. The switching range of the
cutoff function
f
c
r
is the same for all atomic pairs, namely,
0.3 Å
r
0.9 Å.
III. THE FITTING PROCEDURE
The calculations of energies, force constants, equilib-
rium bond distances, and angles in isolated molecules were
performed with the
GAMESS
program.
18
All calculations were
done within DFT using B3LYP density functional and the
6-311G
*
basis set. The force constants were calculated from
seminumerical Hessians in
GAMESS
. The properties of
graphene and diamond were calculated with the
ABINIT
program
19
within DFT, using the Perdew–Burke–Ernzerhof
exchange-correlation generalized gradient approximation
functional.
20
Different bond types were labeled by listing the symbols
of the elements forming the bond followed by the numbers
of neighbors of each type on each side of the bond in paren-
theses. For instance, CO
2110
designates a carbon-oxygen
bond with two carbons and one oxygen connected to a
carbon atom and one carbon and no oxygens connected to an
oxygen atom. When denoting bonds between identical at-
oms, the left and right parts in parentheses can be exchanged,
e.g., both CC
1000
and CC
0010
refer to the same bond.
The angular spring constants in Eq.
13
were obtained
from the calculated Hessians and equilibrium bond distances
of linear triatomic molecules. The values of the coefficients
in Eq.
13
for various combinations of atoms are given in
Table
II
. The dihedral coefficient in Eq.
14
for C atoms
TABLE III. Calculated parameters of the molecules used in the fitting database. Molecules are in a singlet state,
unless indicated otherwise. Bond types in column 3 are labeled by the symbols of the elements forming the
bond, i.e., CO, followed by numbers in parentheses signifying how many neighbors of each atomic species each
side of the bond has. For instance, CO
2110
designates a carbon-oxygen bond with two carbons and one
oxygen connected to a carbon atom and one carbon and no oxygen connected to an oxygen atom.
Molecule
Atomization
E
eV
Bond Bond order
r
e
Å
k
eV
/
Å
2
O
u
C
carbon monoxide
−10.978 7
CO
0000
2.408 1.1265 123.4
O
u
C
u
O
carbon dioxide
−16.645 8
CO
0100
2.097 1.1609 102.0
O
u
CO
2
−16.083 7
CO
0200
1.521 1.2487
60.4
O
u
CO
3
triplet
T
d
−12.088 1
CO
0300
1.085 1.3913
12.7
C
2
O
triplet
D
h
−10.837 3
CO
0010
1.259 1.2586
34.0
C
3
O
D
3
h
−5.727 22
CO
0020
0.802 1.4683
1.3
C
u
O
u
O
*
singlet
C
v
−9.332 51
CO
0001
1.992 1.1540
96.5
OO
1000
0.779 1.3592
17.8
C
u
O
3
−7.330 92
OO
1100
0.725 1.4118
18.0
CO
0002
1.423 1.2280
59.9
O
u
C
u
C
u
O
triplet
D
h
−20.258 4
CC
0101
1.611 1.2806
72.5
CO
1000
2.075 1.1871
86.4
C
2
O
2
C
2
v
−16.003 3
CC
0002
1.011 1.4208
32.4
CO
1100
1.645 1.2452
35.5
O
u
C
u
C
u
C
C
h
−21.502 7
CC
0110
1.448 1.2952
66.3
CO
1000
2.247 1.1497 108.3
O
u
C
u
C
triplet
C
v
−13.656 2
CO
1000
2.162 1.1627
96.3
CC
0100
1.595 1.3595
40.7
C
4
O
triplet
C
3
v
all C atoms in one plane
−12.977 5
CC
0021
0.948 1.4385
23.9
CO
3000
0.757 1.4184
33.9
C
u
C
u
O
u
C
singlet
a
−16.846 2
CO
1010
0.903 1.3029
41.4
C
u
C
dicarbon
−5.089 28
CC
0000
3.470 1.2514
77.2
C
u
C
u
C
linear tricarbon
−13.711 8
CC
1000
1.715 1.2915
67.5
O
u
C
u
C
u
O
triplet
−20.258 4
CC
0101
1.611 1.281
72.5
C
4
triplet
D
3
h
−14.161 9
CC
2000
1.152 1.4103
28.0
C
5
D
4
h
b
triplet
−18.064 6
CC
3000
0.854 1.4109
26.4
C
4
linear
−17.993 3
CC
1010
1.845 1.2961
63.4
O
u
O
triplet
−5.268 6
OO
0000
1.748 1.2065
79.1
O
3
ozone, singlet
C
2
v
−5.811 0
OO
0100
1.248 1.2584
45.8
graphene
c
−7.917 3
CC
2020
1.303 1.4189
45.0
diamond
c
−7.773 7
CC
3030
1.000 1.5417
24.0
a
Triplet calculation diverges.
b
Singlet state diverges; pentet state is higher in energy.
c
Energy per atom.
234706-4 A. Kutana and K. P. Giapis
J. Chem. Phys.
128
, 234706
2008
Downloaded 20 Jun 2008 to 131.215.225.9. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp
k
CCCC
tor
=1.583 eV was deduced by comparing the energy of a
graphene sheet with that of a hypothetical H-6 carbon
structure,
21
in which atoms have the same threefold local
coordination as in graphene, but in which one-third of dihe-
dral angles are 60° instead of 0°. All other dihedral interac-
tions were not considered. The bond energies
U
ij
are found
by subtracting the angular and dihedral components from the
total energy of the molecule. Since bonds of the same type
can occur in more than one molecule, there is an ambiguity
in the choice of molecules for the fitting. Here, molecules
with the least number of atoms representing a given bond
type were employed.
In symmetric bonding configurations, when all bonds are
of the same type
e.g., a CO
0200
bond in a
D
3
h
CO
3
, the
assignment of the energy to each of the bonds is obvious. In
other molecules, such as a linear C
u
O
u
O, where only a
sum of the energies of the two bonds is known, the bond
energies were defined to be proportional to the bond orders
as calculated by
GAMESS
.
CO bond
. The results of calculations for the CO bond are
shown in Table
III
. Since the CO compounds are usually
gases at normal conditions, their binding energy must
quickly decrease with coordination. The calculated molecu-
lar energies listed in the table demonstrate that this is indeed
the case. In diatomic carbon monoxide, carbon and oxygen
form a strong covalent bond with the energy of −11 eV.
Connecting a second oxygen atom to the carbon atom results
in a carbon dioxide molecule, predicted to be 5.7 eV lower in
energy than CO. Being more stable than carbon monoxide,
CO
2
nevertheless has less energy per each of its two
CO
0100
bonds. The molecular stability starts decreasing
when more oxygen atoms are added on the C side of the CO
bond. The calculated energy of
D
3
h
carbon trioxide CO
3
Ref.
22
is slightly above CO
2
, and its bond energy is also
smaller. The energy of
T
d
CO
4
is above that of CO
3
by
4.0 eV, although the formal bond strength is higher than in
CO
3
due to high contributions from the angular terms in this
molecule. The molecular stability decreases even more rap-
idly when adding carbon atoms on the O side of the CO
bond. The linear C
2
O molecule in its
X
3
−
ground state
23
is
less stable than the carbon monoxide, while
D
3
h
C
3
O has the
energy of only −5.7 eV.
OO bond
. The calculations for the singlet ozone mol-
ecule predict an equilibrium O
u
O
u
O angle of 118.2°, and
an atomization energy of −5.81 eV, which is 0.54 eV lower
than the energy of the triplet ground state of O
2
. The param-
eters of the OO
0101
bond were determined from the cal-
culation of the
D
4
h
square
O
4
molecule, predicted to be less
stable than O
3
−5.58 eV
. Larger O atom clusters were not
considered due to their instability. The energies of the OO
bonds in configurations involving neighboring C atoms are
listed in Table
III
.
CC bond
. Due to the greater electron delocalization in
various spatially extended structures formed by carbon, local
potentials may be expected to be less suitable for description
of these structures than various CO molecules. While delo-
calization effects must be most pronounced in periodic solids
exemplified by graphene or diamond, they are also found to
be evident in small structures such as uniform rings of
carbon atoms.
Carbon
C
n
rings and linear carbon chains
. Uniform car-
bon rings C
n
have
n
equivalent CC
1010
bonds contributing
equally to the total energy of the molecule. In order to
maximize the scope of the fit, C
n
rings of various sizes with
n
between 5 and 50 were calculated. From Eq.
13
with
0
=
, it follows that
E
total
C
n
ring
=
nU
CC
1010
+
2
k
CCC
2
n
.
18
According to Eq.
18
, the plot of
E
total
C
n
ring
/
n
vs 1
/
n
2
should be linear, yielding
k
CCC
and
U
CC
1010
from the slope
and intercept, correspondingly. Figure
2
shows this plot for
n
=5–33. The validity of the approximation that the total
energy can be represented as a sum of the bond energy and
angle energy can now be examined. Three distinct lines are
generated, corresponding to
n
=4
k
+2,
n
=4
k
, and
n
=2
k
+1.
Within each of these three groups, it is permissible to de-
scribe the total energy with Eq.
18
. The most stable mol-
ecules are aromatic rings with
n
=4
k
+2, while the least
stable are the rings with the even number of atoms not sat-
isfying the aromaticity condition. The odd-membered rings
have intermediate stability. The stability of ring molecules
with
n
=4
k
+2 originates from the formation of closed shell
configurations, as described by Huckel’s rule. The fitted val-
ues are given in Table
IV
, giving different types of mol-
ecules’ different angular stiffness, while the bond energies
are the same in all molecules. These calculations demon-
strate the limitations of a local potential in approximating the
energies of these molecules. Since the local potential cannot
estimate the total number of atoms in the ring, a global quan-
TABLE IV. Coefficients for the linear fittings of the energies of C
n
rings as
described by Eq.
18
.
U
CC
1010
eV
k
CCC
eV
4
n
−6.308 219 187
3.872 626 914
4
n
+2
−6.327 616 194
2.116 749 568
Odd
−6.300 275 152
3.087 086 693
FIG. 2. Calculated energy of
n
-atom carbon rings C
n
normalized by the
number of atoms
n
as a function of 1
/
n
2
.
Circles
n
=4
k
+2;
squares
n
=4
k
;
triangles
odd
n
. The lines are linear fits to the calculated points.
234706-5 Analytical carbon-oxygen reactive potential
J. Chem. Phys.
128
, 234706
2008
Downloaded 20 Jun 2008 to 131.215.225.9. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp
tity, it cannot account for the differences between the lines in
Fig.
2
. For the purposes of fitting the CC
1010
bond in the
current potential, the values of
U
CC
1010
and
k
CCC
were
averaged over the corresponding values for
n
=4
k
+2,
n
=4
k
,
and
n
=2
k
+1 molecules.
Unlike homogeneous C
n
rings, linear carbon molecules
have
n
−3
CC
1010
bonds, plus two terminating CC
1000
bonds. The energy of a linear carbon chain with
n
atoms is
E
total
C
n
linear
=
n
−3
U
CC
1010
+2
U
CC
1000
.
19
According to Eq.
19
,
E
total
C
n
linear
yields a line
when plotted as a function of
n
. Figure
3
shows the calcu-
lated
E
total
C
n
linear
as a function of
n
for
n
=3–50, as well
as a linear fit, yielding
U
CC
1010
=−6.327 eV and
U
CC
1000
=−0.6403 eV. As expected, the obtained value of
U
CC
1010
is
very close to the value
U
CC
1010
=−6.321 eV obtained from
the fitting of C
n
rings.
Graphene structures
. The energy and equilibrium dis-
tance of the CC
2020
bond in a single graphene sheet ob-
tained with
ABINIT
simulations are given in Table
III
. The
force constant of 45.0 eV
/
Å
2
based on the experimentally
measured elastic stiffness coefficients C
11
and C
12
of
graphite
24
was used for this bond.
The energies of defects in graphene were calculated with
GAMESS
using a C
56
H
18
molecule shown in Fig.
4
. After
optimizing the geometry of the molecule, atoms 1–6 were
removed and the energy of the resulting structures was cal-
culated. Removal of atoms 1–6 is equivalent to the elimina-
tion of 12 CC
2020
bonds, and conversion of another 12
CC
2020
bonds into CC
1020
. The obtained energy of the
CC
1020
bond is given in Table
III
. The parameters of the
CC
3030
bond were obtained from the calculation of the
equilibrium structure of diamond. Other carbon-carbon
bonds are also given in Table
III
.
IV. MD SIMULATION OF GRAPHENE UNZIPPING
BY OXYGEN
As a test of our potential, a classical molecular dynamics
MD
simulation was performed to describe the ordered oxi-
dation and subsequent unzipping and cracking of graphene
sheets by chemisorbed oxygen. Using
ab initio
DFT calcula-
tions, Li
et al.
9
first described the reaction mechanism that
causes oxygen-driven unzipping of graphitic materials. They
found that cooperative unzipping ensues as a result of the
formation of linear arrays of epoxy groups on opposite sites
of hexagons in graphene. In our simulation, we adopted the
same setup: oxygen atoms were placed above C
u
C bonds
in graphene along the opposite sides of hexagons, as shown
in Fig.
5
a
. Periodic boundary conditions were assumed.
The system was attached to a 1200 K thermostat and then
allowed to evolve according to Newtonian mechanics. Fig-
ures
5
b
–
5
d
depict snapshots of the unzipping process of
graphene by adsorbed oxygen atoms. The first C
u
C bond is
cleaved by an O atom at
t
=349 ps, creating a defect site.
Then, unzipping is initiated at this site and proceeds serially
in both directions along the linear array of adjacent oxygen
atoms. Figure
5
c
captures the breakage of the second
carbon-carbon bond next to the first epoxy group, which oc-
curs at
t
=459 ps. Eight C
u
C bonds unravel after 1061 ps,
as shown in Fig.
5
d
. All C
u
C bonds become disconnected
at
t
=1163 ps. All bond cleavages subsequent to the first one
occur in a chain fashion next to the already cleaved bonds,
taking on average 90 ps per each C
u
C bond. This unzip-
ping process is irreversible: No restoration of the CC bonds
has been seen in the simulation up to
t
=1800 ps. After cleav-
ingaC
u
C bond, each oxygen atom remains connected to
two carbons, oscillating between two buckled positions on
either side of the graphene plane. Our results agree well with
the observations and predictions made by Li
et al.
9
for a
FIG. 3. Calculated energy of n-atom linear carbon molecules C
n
as a func-
tion of
n
. The line gives a linear fit to the calculated points.
FIG. 4. A
GAMESS
-optimized structure of the C
56
H
18
molecule. Numbers
1–6 mark the atoms that were removed in order to evaluate the energy of the
CC
1020
bond. The illustrations of molecules have been prepared using the
package
VMD
Ref.
26
.
FIG. 5.
Color online
Snapshots from the simulation of graphene unzipping
by oxygen atoms at
T
=1200 K. The snapshots are taken at
a
t
=0 ps,
b
t
=349 ps,
c
t
=459 ps, and
d
t
=1061 ps.
234706-6 A. Kutana and K. P. Giapis
J. Chem. Phys.
128
, 234706
2008
Downloaded 20 Jun 2008 to 131.215.225.9. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/jcp/copyright.jsp