Positional stability of some diamondoid and graphitic
nanomechanical structures: A molecular dynamics study

G. I. Leach, R. E. Tuzun, D. W. Noid, B. G. Sumpter

Abstract:

Molecular dynamics simulations indicate positional stability to be an
important issue in a wide variety of molecular nanotechnology
applications. It can determine the difference between the success and
failure of mechanical nanodevice designs. Diamondoid materials are
proposed for many such designs in part because of the
stiffness and strength of diamond compared to other materials. These
properties should allow the problems of positional stability to be
minimized while simultaneously minimizing atom count (bulk, molecular
weight, etc.). Because of their synthetic availability and desirable
mechanical properties, graphitic materials, in particular carbon
nanotubes and buckyballs, are also components in many proposed
nanomachine designs.

We present results of a molecular dynamics study of the positional
stability of several diamondoid and carbon nanotube structures,
including comparing similar functional designs made from both these
materials. In particular, we investigate hollow tubular and solid
rectangular ``struts'' -- which, for example, could form part of the
linear actuators required in a Stewart platform molecular positioning
device. The effects of constraints on positional stability and
nanomachine performance are examined. Quantum mechanical
results for some of these structures are briefly compared with
molecular dynamics results.

Finally, we discuss the use of the web-based virtual reality modeling
language -- VRML 2 -- to publish the results of molecular dynamics
simulations on the world wide web as three-dimensional movies.

For the purpose of feasibility studies of proposed
nanomachine designs, which have
not yet been built, it is necessary to perform some sort of
computational analysis [1,2].
One common, if not the most common,
technique is molecular dynamics (MD) because it is
quick and versatile and because it can be used to study
processes at the atomic scale which cannot be directly
observed experimentally. Within the last several years,
several groups have performed MD simulations on bearings
[3], motors [4], gears [5,6],
and other systems in order to study their mechanical properties
and to explore nanomachine design considerations.

One of the most important issues in nanomachine design is
positional stability. Nanomachine components must be fairly
stiff in order for bearings to rotate smoothly, for gears to
mesh well, for sub-angstrom resolution on positional control
devices, and so on. This issue can be investigated with MD.
Results so far indicate that the geometry of a nanostructure
plays as crucial role in its positional stability as its
potential energy surface. For example,
using MD Sumpter and Noid [7] recently found the
dynamical stability of carbon nanotubes to depend in large part on
the aspect ratio.

More recently, the correspondence of classical and
quantum mechanics has begun to emerge as an area of interest to
nanotechnology. Of particular interest to us is the well-known
zero point energy problem. Classically, the flow of energy is
unrestricted. Quantum mechanically, however, some energy must
always remain locked within each vibrational mode. This
could mean that classical simulations may sometimes contain
extraneous vibrational motion from the flow of what should be
zero point energy. Results of a new method for
calculating vibrational ground states of systems with many atoms
[8,9] indicate that this is indeed a serious issue,
especially for systems with sparsely connected bond networks and
with no constraints or external forces. This is actually good
news for nanotechnology because it implies that many nanomachine
designs may perform better than MD results indicate and/or
require fewer atoms for the same positional stability.

Adding constraints or collapsing degrees of freedom can help
reduce extraneous vibrational motion. This can be done in several
ways using methods which are already used for computational savings or
other reasons. Internal cordinate methods [10] have been used
by the Goddard group [11] for several biological and
nanotechnology applications. In recent nano-fluid dynamics
simulations [12,13] the ends of a nanotube were frozen in
order to not allow an axial fluid flow to drag the tube along. In the
gear simulations performed at NASA Ames [6] the ends of the
carbon nanotube shaft were constrained to not elongate but were still
allowed to move within a plane transverse to the tube symmetry axis.
Finally, in two recent papers rigid-body (quaternion) dynamics methods
were compared to fully atomistic MD for the operation or assembly of
molecular bearings [14,15].

In this paper we investigate the positional stability of
diamondoid blocks and of carbon nanotubes by MD for several geometries,
constraints, and initial conditions and using several potential energy
surfaces. We include both small aspect ratios,
which are so far the most commonly modeled, and large aspect
ratiso, which have not been modeled as extensively and which would
appear in some nano-hydraulic components or support struts for
Stewart platforms or other positioning devices. As an illustration
of the considerations for positional stability, we begin to
address the question
of how thick a support strut should be for a Stewart platform,
given a target compliance on the order of a bond length (atomic
precision).

The details of the molecular dynamics methds used in this study
have been well-described elsewhere and are not repeated here.
Essentially, MD consists of integrating the classical equations of
motion over small time steps (in our case, 1 fs). Initial momenta
are set to correspond to desired temperatures, rotational velocities, or
other desired initial conditions. Overall translational and
rotational motion are removed before the beginning of simulation
so that only random thermal motion
remains. External forces maybe applied at any
time. The simulations presented here were performed at constant
energy, using symplectic integration [17] and recent improvements
to the bonded interaction portion of the calculation [18,19,20].

For the sake of comparison, several potential energy surfaces were
used in this study. For terminated diamond blocks, we used a modified
MM2 potential [21] with and without non-bonded interactions.
The unmodified potential, which is most suitable for molecular mechanics
calculations, has a cubic bond stretch term, which is physically
unreasonable at moderate to large bond distances. For our MD simulations
we used a Morse functional form, which is well-behaved at
all bond distances, with the same equilibrium position and curvature
as the cubic potential. For unterminated diamond blocks, we used
a special potential with stretch and bend interactions. For carbon
nanotubes, we used the graphite potential energy surface of Guo
et. al. [22] with and without torsion interactions.

Studying positional stability using MD requires starting at the
equilibrium geometry or
at least near a local minimum in the potential energy surface. Two
methods may be used to obtain initial geometry: molecular mechanics
[23], or, as in this work, annealing
(also called dynamical steepest descent). If a structure is away
from equilibrium and at rest, it will move toward equilibrium.
Annealing consists of performing an MD simulation in which kinetic
energy is
gradually removed until the structure is nearly motionless and almost
completely relaxed. We demonstrate later the dramatic effects of
incomplete annealing on subsequent MD simulation.

Unterminated diamondoid rectangular blocks with approximately
square cross section were generated for three aspect
(length to width) ratios: 1, 10, and 100 (hereafter referred to
as short, medium length, and long) using the CrystalSketchpad [24]
package. In all of these blocks, the surfaces were either (110) or
(100) planes. To test the effects of bond network connectivity
two thicknesses were chosen: one unit cell (about 0.4 nm) and 1 nm.

Figure 1:
Unterminated diamondoid blocks, side view.

(a) AR=1, CS=uc.

(b) AR=10, CS=uc.

(c) AR=100, CS=uc.

(d) AR=1, CS=nm.

(e) AR=10, CS=nm.

(f) AR=100, CS=nm.

Figure 2:
Unterminated diamondoid blocks, end view.

(a) CS=uc

(b) CS=nm

These blocks were then used to generate terminated blocks (no dangling
bonds). The (110) surface was terminated with hydrogens. Rather than
terminating with hydrogens on the (100) surface using the postulated
surface rearrangement [28], the dangling bonds were bridged
with oxygen atoms. Similar strategies have been used in several
proposed nanomachine designs. When inserting the bridging oxygen
atoms, due to the geometrical constraints from the diamond surface it
was possible to insert the oxygen at a point where either the C-O bond
lengths or the O-C-O bond angle was correct, but not both. We chose
to begin with correct bond lengths. The final adjustment of the
geometry was performed by annealing, as described above.

Figure 3:
O-bridged, H-terminated diamondoid blocks.

(a) Side view.

(b) End view.

We use the following abbreviations to describe simulation conditions:

Figure 3:
O-bridged, H-terminated diamondoid blocks.

Symbol

Quantity

Possible Values

AR

aspect ratio

1, 10, or 100

CS

cross section

uc (unit cell) or nm

S

surface termination

C (unterminated) or O,H (oxygen,hydrogen)

T

temperature

150K or 300K

NB

non-bondeds

on or off

NC

constraints

0, 1, or 2 constrained ends

All of the above blocks were run for 20 ps. In some cases,
simulations were run with one or both ends of the block frozen. The
two key indicators of position stability of interest, end-to-end
distance and maximum transverse displacement, were obtained from the
saved trajectory files.

Similar runs were performed for carbon nanotubes with 20 atoms
per ring (about 8 Åin diameter). This diameter was chosen because
it has been used in many of our previous simulations. Tubes with aspect
ratios of 1 (5 rings), 10 (39 rings), and 100 (375 rings) were generated
and annealed with and without torsion interactions.

Figure 4 shows profiles of maximum transverse displacement
for several diamondoid block simulations. Figure 4(a)
shows the effects of temperature and aspect ratio for blocks with an
unterminated nm cross section. Initial conditions for the 150K and
300K simulations are exactly the same except for a proportionality
factor for the momenta. Gross features of the dynamics, and some fine
features, are therefore similar for both temperatures. The maximum
transverse displacement is about 6 Åin this set of simulations, for
the longest block and highest temperatures. Constraining the ends
changes this very little; when both ends are constrained, the maximum
transverse displacement drops to about 4 Å.

Figure 4:
Maximum transverse displacements for diamondoid blocks.

Similar trends are observed for end-to-end distance (Figure
5). Results are shown for medium and long
unterminated diamondoid blocks at temperatures of 150K and 300K,
for unit cell and 1 nm cross sections. Here, refers to
the change in end-to-end distance from the equilibrium value.

Unlike the nm cross section, the unit cell diamond blocks build up
appreciable torsional motion. The blocks with aspect ratios of 1 and
10 remain more or less straight. However, the longest block partially
coils, and so the maximum transverse displacement continually
increases throughout the simulation. Movies of this simulation show
that several flexion modes become excited. Unsurprisingly,
constraining both ends (not just a single end) reduces the maximum
transverse displacement considerably, down to about 5 Å
(Figure 4(b)).

For the 1 nm cross section blocks, surface termination appears to have
very little effect on positional stability (Figure 4(c)).

Turning non-bonded interactions on or off also appears to have
little effect. This is important because the non-bonded calculations
are the most expensive portion of the simulation and turning off
non-bonded interactions can save considerable time in initial
feasibility studies.

As observed previously, the long unit cell cross section block
partially coils. This shows up in the end-to-end distance as well as
the transverse displacement. In addition, as long as the block
doesn't coil, the end-to-end distance is more obviously periodic than
the transverse displacement. It is interesting to note that the
dominant vibrational frequency is inversely proportional to the
length, as observed previously for carbon nanotubes [7].

Trends in the classical dynamics of carbon nanotubes are similar
to those for the diamondoid struts.
Figure 7 shows profiles of maximum transverse displacement for
nanotubes with diameters of about 8 Å. Only aspect ratios of 10 and 100
are shown for the sake of clarity.
In the first 25 ps, the maximum transverse displacement reaches
a maximum of about 2.5 Åfor the 300 K simulation. After this, it reaches
4.3 Å(a 5% strain).
The aperiodicity of the transverse displacement profile
appears to indicate that multiple modes are excited.

The end to end distance profiles are consistent with this interpretation.
For the longest tube, the deviation from equilibrium end to end
distance has a period of about 8.2 ps. However, the minimum and maximum
deviation change between periods, meaning that some energy is transferring
between modes.

Figure 7:
profiles for carbon nanotubes.

In these simulations, torsion interactions were included. Simulations
in which torsion interactions were excluded showed similar dynamics as the
above results.

In Figure 8
we illustrate the effects of incomplete annealing on subsequent MD
results. When the terminated, 1 nm thick, medium length block was
only partially annealed, a spurious flexion (S) mode developed. This
occurred because in the partially annealed structure all of the O-C-O
bond angles were too large. At the start of the MD simulation, these
bond angles begin to decrease, forcing the tube into an elongation mode
and then a flexion mode.

One reliable indicator of incomplete
annealing is too high a temperature in the simulation. Initially, for
an average target temperature of 150K, the initial momenta are set to
give a temperature of 300K. This is done because it is normally expected
for energy to be evenly distributed among potential and kinetic
energy on the average. After an initial transient, the temperature
decreases, but settles about an average of about 180K
(Figure 9).
When the structure is completely annealed, the average temperature settles
around an average of 150K, and no flexion mode develops.

The quantum mechanical ground state wavefunction of a molecule
contains probability distributions for bond lengths and angles and
for overall features such as end-to-end distances. From
these distributions the dynamical stability at the low temperature
limit can be estimated.

We would expect the positional stability of
diamondoid blocks to be greater than that for
more sparsely connected bond networks such as polyethylene chains.
The ground state of model polyethylene chains has been explored
elsewhere [8]. The quantum mechanical
half-width in end-to-end distance is on the order
of 0.1 Å, which is far less than the classical results
for diamondoid blocks -- most especially for those with a unit cell
cross section. Although the end-to-end distance
oscillations from the diamondoid simulations are larger by more than an
order of magnitude,
they would probably have little effect on the overall performance of
nanomachines using these blocks -- at least for those with
a 1 nm cross section. Transverse displacement, however, may
still be problematically large in some MD simulations.

The classical and quantum simulation results for carbon nanotubes
are compared in detail in [9]. Essentially, the disagreement
between classical and quantum end-to-end distance distributions appears
to not seriously affect overall simulation results as long as the tubes
are well enough constrained.

The virtual reality modeling language (VRML) is a mechanism for
publishing three-dimensional (3D) models on the world wide web. The
most recent version of VRML (VRML 2.0 [25]) has added
mechanisms for describing motion and behaviour of 3D objects,
including specialised nodes and a general purpose scripting facility
(supporting Javascript and/or Java). VRML 2.0 browsing capability
for the latest versions of both Netscape Navigator and
Microsoft Internet Explorer is either built into the browser or
available as plug-ins.

So far the MPEG digital movie format has been the most common
mechanism for publishing MD trajectories on the web. The use of VRML
2.0 to publish MD trajectories offers the same advantages over MPEG
movies that
providing 3D object data has over 2D image data in general: it gives
back to the viewer control of viewing, lighting, culling and other
parameters and allows them to dynamically position and orient the
objects for improved comprehension and focus.

Our simulation results have been published as VRML 2 trajectories on
the web [26]. A comparison of the use of VRML 2 and MPEG for
this purpose may be found in [27].

The positional stability of diamondoid blocks and
carbon nanotubes in a wide variety
of sizes and initial conditions has been investigated using MD.
Two major criteria for positional stability are end-to-end
distance and maximum transverse displacement. Both of these showed
reasonable trends for temperature, cross-sectional thickness, and
length. Details such as surface termination and functional form of
the potential energy surface played a minor role compared to the
arrangement of the bond network.

So far, MD has been used for initial feasibility studies for several
types of nanomachine components. Later, when the designs are further
refined, it will be necessary to study finer features of the dynamics.
At this level of detail, the simulation data may reflect more about
the limitations of MD than about the dynamics of the problem. In the
nanotechnology simulation results reported so far, constraints and
external forces have tended to limit the magnitude of structural changes
that would harm simulation results for overall trends in nanomachine
performance. However, in future simulations new
correction strategies may be needed.

For the Stewart platform, MD simulations indicate
that a support strut should be
at least 1 nm square in cross section. However, previous QMC results appear
to indicate that this thickness could serve to overengineer the design
requirements, and a smaller thickness may suffice. To refine this estimate
it would be necessary to perform additional simulations, under different
compressive and transverse loads and other operating
conditions.