Tahir Cagin, Jinwei Che, Yanhua Zhou, Stephen J. Walch+ and

William A. Goddard III,

Materials and Process Simulation Center, 139-74

Division of Chemistry and Chemical Engineering

California Institute of Technology

Pasadena, CA 91125, USA

+) ELORET, 690 W. Fremont Ave, Suite 8, Sunnyvale, CA 94087-4202


Advances in theory and methods making it practical to consider fully first principles (de novo) predictions of surface and interface structures and the tribological properties of materials. However, despite the progress there remain enormous challenges in bridging the vast range of distances and time scales between de novo atomistic simulations and the quantitative continuum models for the macroscopic systems essential in industrial design. Recent advances relevant to such developments include: Quantum Chemistry including continuum solvation, de novo Force Fields to describe reactions and phase transitions, Molecular Dynamics to represent a steady state system, Nonequilibrium MD for studying rheology of lubricants. To provide some flavor for the opportunities we will illustrate some of the progress and challenges by summarizing some recent developments in methods and their applications. Four different applications will be covered. 1) Friction and wear mechanisms on clean diamond surfaces, 2) Rheology of confined lubricants, 3) Nano-tubes as nano-tribological probes, and 4) Reactions on surfaces of diamond.



1. Introduction

In order to develop new materials and compositions with designed new properties, it is essential that these properties be predicted before preparation, processing, and experimental characterization. Despite the tremendous advances made in the modeling of the structural, thermal, mechanical and transport properties of materials at the macroscopic level (finite element analysis of complicated structures) there remains tremendous uncertainty about how to predict many critical properties related to performance. The fundamental problem here is that these properties depend on the atomic level interactions and chemistry (eg. making and breaking of bonds) dealing with the electronic and atomic level description at the level of nanometers and picoseconds, while the materials designer needs answers from macroscopic modeling (finite element paradigm) of components having scales of cm and milliseconds or larger. To dramatically advance the ability to design useful high performance materials, it is essential that we insert the chemistry into the mesoscopic and macroscopic (finite element) modeling.

The difficulties in doing this are illustrated in Figure 1, where we see that vast length and time scales separate the Quantum Mechanics (QM) from the macroscopic world of engineering design. Tremendous advances have been made recently in first principles QM predictions of chemical reactions, but the state of the art can handle accurately reactions with only ~50 atoms. There is no practical approach to carrying out a QM calculation on the initiation and propagation of a crack through a stabilized zirconia ceramic. Despite this difficulty, the computations MUST be based on accurate first-principles quantum mechanics if we are to predict the properties of new materials.


Figure 1 Multiscale modeling hierarchy



Our strategy for accomplishing this objective is to develop an overlapping array of successively coarser modeling techniques where at each plateau (a range of length and time scales), the parameters of the coarse description are based on the first principles based parameters of the immediately finer description, as illustrated in Figure 1. Thus based on accurate QM calculations we find a Force Field (FF) including charges, force constants, polarization, van der Waals interactions etc that accurately reproduces the QM. With the FF, the dynamics is described with Newton's equations [Molecular Dynamics (MD)], instead of the Schrodinger Equation. The MD level allows one to predict the structures and properties for systems ~ 105 times larger than for QM, allowing direct simulations for the properties of many interesting systems. This leads to many results relevant and useful in materials design, however, many critical problems in materials design require time and length scales for too large for practical MD.

Thus we need to develop methods treating the mesoscale in between the atomic length and time scales of MD and the macroscopic length and time scales (microns to mm and m sec to sec) of Finite Element Analysis (FEA). This linking through the mesoscale in which we can describe microstructure is probably the greatest challenge to developing reliable first principles methods for practical materials design applications.

Only by establishing this connection from microscale to mesoscale will it be possible to build first principles methods for describing the properties of new materials and composites (the domain of materials science and engineering) in terms of fundamental principles of physics and chemistry. Thus, for fundamental prediction is to play a direct role in materials innovation and design, it is essential to bridge the micro-meso gap The problem here is that the methods of coarsening the description from atomistic to mesoscale or mesoscale to continuum is not so obvious as it was in going from electrons to atoms.

Given the concepts, it is necessary to carry out calculations for realistic time scales fast enough to be useful in design. This requires developing software tools useful by design engineers, by incorporating the methods and results of the QM to MD to mesoscale simulations. To accomplish the goals of developing methods for accurate calculations of materials and properties, we focus on (i) implementations that make full use of modern highly parallel computers and (ii) building in knowledge based heuristic methods of accessing this information automatically so that designers can focus on the macroscopic issues without concern for the details of theory and simulation. At this point, we expect a revolution in materials design applications where the use of first principles multiscale simulations allows an every increasing amount of the design to be done on the computer before experiment.

2. Progress in methods developments

Our strategy is to transcend from the most fundamental theory (quantum mechanics, QM) to practical engineering designs in a sequence of several levels as indicated in Figure 1.

2.1 Quantum Mechanics

2.1.1 Ab initio Quantum Chemistry Applications

It is important to use QM to describe systems in which bonds are being broken and formed. Only then can we be sure to obtain accurate barrier heights and bond energies. The modern methods of QM (GVB1, PS-GVB2, MR-CI3, and GDS-DFT4) can give accurate barriers for reactions useful in studying the properties of nanoscale materials. However, despite the progress in first principles electronic structure theory, the calculations are often far too slow for studying the dynamics of nanotechnology applications. Thus is necessary to have general approaches for averaging the electrons from QM to obtain a FF in terms of atomic positions. For finite molecules; A new methodology (PS-GVB) combining pseudo-spectral (PS) multi-grid and de-aliasing strategies with the sophisticated many-body wave functions [generalized valence bond (GVB)] were implemented and applied to large scale problems. PS-GVB led to considerably better scaling with size (N2 rather than the N3, N4, N5, N6. characteristic of alternative methods) and simpler parallelization. PS-GVB has been extended to treat all atoms of the periodic table (using core effective potentials), handle new sophisticated wave functions (GVB-RCI, MP2), and describe important properties (solvation energies, reaction rates, activation/reaction barriers). We will optimize these methods for parallel implementations, and extend the methodology to include GVB-RCI-MP2 and self-consistent GVB-RCI.

2.1.2 Density Functional Theory Applications

Most practical materials properties require a description of infinite systems using periodic boundary conditions (PBC). This is three-dimensional (3D) for bulk properties or two-dimensional (2D) for surface growth and interfaces. For this purpose, we have developed a new method, Gaussian dual space density functional theory (GDS-DFT), in which most parts scale linearly with N. In implementing this, we have developed a new separable pseudo-potential that can be applied to all atoms of the periodic table. We reformulated the theories for electronic structure calculations of periodic systems in a way suitable for large-scale calculations using Gaussian basis functions. An accurate grid is introduced for efficient calculation of matrix elements. A dual-space approach is used to calculate the Coulomb potential with computational cost that scales linearly with the size of basis set. Preconditioned generalized conjugate gradients approach is introduced for rapidly converging wave functions expressed in terms of Gaussian basis functions. This method is applied to a variety of systems with excellent results.

2.2 Force Fields (FF)

Using quantum mechanical results we develop FF descriptions to provide the energetics needed for the simulations of the nano-phase materials and their properties. The FF must even be accurate enough to obtain the proper energy differences for representing phase behavior of the materials and transferable so that one can apply it to phase transformations and interface phenomena. Standard FF generally use simple springs to represent bonds and angles in describing structures and vibrations of molecules. This is not generally sufficiently accurate to obtain a FF that accurately describes the properties of a specific class of molecules or polymers. For better FF, we fit to the QM using the Hessian-Biased FF (HBFF)5 approach, which combines normal mode information from HF theory with the frequency information from theory or experiment. This HBFF approach has been used to develop accurate FF for polymers (e.g., PE, PVDF, nylon, POM, SiH) 6-10, ceramics (e.g., Si3N4, C3N4), 11,12 semiconductors13 and metals.14

On the other hand, for fast qualitatively considerations of new systems, we find that generic FF suitable for general classes of systems are most useful. DREIDING FF15 (for the main group elements) and the Universal force field (UFF)16 (all elements, including inorganics and organometallics) are such FF. In studying the surface and interface properties of materials, especially in the study of friction and wear mechanisms in materials, one needs to have accurate environment dependent, and/or reactive force fields derived from first principles. In recent years we made critical advances made in developing ab QM based FF for describing

  1. metals where many-body interactions play critical role on their physical properties;17-18 In these force fields we have described cohesion in metals through many body interactions or through an embedding energy. The parameterization of the force fields is based on DFT level quantum mechanical calculations on pure metals.
  2. oxides, ceramics, and zeolites where competition between ionic and covalent bonding is often very important, especially in describing polymorphic phase transitions, reactions, surface and interface properties. 19-21 Here, we have used a Morse stretch two body pairwise interaction to account for short range interactions, whereas the electrostatic interactions are environment dependent. The charges of the ions are assumed variable and reevaluated at each configuration using charge equilibration.
  3. covalent-bonded system such as carbon, hydrocarbons, silicon, germanium and their behavior far from equilibrium where the description of bond breaking and forming must be a part of an accurate classical description. 21-22 Bond order dependent force fields are truly reactive force fields and extremely useful in studying the friction and wear properties of such systems.

2.3 Molecular Mechanics and Molecular Dynamics: MPSim

Using these force fields (FF) in large-scale Molecular Dynamics simulations allows practical calculations on up to several millions of atoms. Our objective is developing new strategies and algorithms in addition to taking advantage evolving hardware and software technologies to extend the time and distance scale to 100s of ns and close to microns. This will involve using fast multi-pole techniques, multiple time step approaches, NEIMO (Newton Euler Inverse Mass Operator) method, CFA (Constrained Force Algorithms) method, hyper molecular dynamics approaches where they are suitable.

The focus here is on extending the methods of MD to physical systems of molecules, polymers, liquids, and inorganic materials with up to 100 million atoms while accurately treating long-range interactions using the cell multi-pole method (CMM). For fast internal coordinate dynamics on a million atoms, we have developed the Newton-Euler Inverse Mass Operator (NEIMO) method. This methodology handles periodic systems and will be extended (Gibbs-Ensemble molecular dynamics. The applications for this method will enable us to investigate the long time dynamics of liquid polymer and solid interfaces; which has tremendous impact on broad range of technological applications; such as wetting, adhesion, phase separation, coatings.

In MD simulations, the FF is used to predict the equations of motion. This leads to trajectories [xi(t), vi(t):× × × =1, × × × ,3N] that can be analyzed (using statistical mechanics and thermodynamics23 principles) to obtain macroscopic properties. MD simulations of heterogeneous nano-phase materials may require millions of atoms to be considered explicitly (a 25 nm cube of polyethylene has 1 million atoms). The most time-consuming aspect of the MD simulations of large systems is the accurate evaluation of the long-range interactions (electrostatic and dispersion), which decrease slowly with distance. Without cutoffs, this cost scales as order (N2) for N particles (a system of 10 million atoms leads to 1014 terms to be evaluated each step). However, cutoffs can lead to excessive errors. For a periodic system, the Ewald procedure leads to accurate summations for these interactions, but the problem scales as N1.5, totally impractical for systems with million atoms.24,25 In order to simulate systems with millions of atoms, we developed methods and optimized parallel computer programs efficient for high capacity MD with the following advanced features:

  1. Cell Multipole Method 26 (CMM) which dramatically reduces the cost of long-range Coulomb and van der Waals interactions while retaining high accuracy. The cost scales linearly with size, allowing atomic-level simulations for million atom systems.27-29
  2. Reduced Cell Multipole Method30 (RCCM) which handles the special difficulties with long-range Coulomb interactions for crystals by combining a reduced unit cell plus CMM for interaction of the unit cell with its adjacent cells. The cost scales linearly with size while retaining high accuracy, allowing simulation of crystals having a million atoms per unit cell (the major use is for models of amorphous and semi-crystalline materials).
  3. Newton Euler Inverse Mass Operator method (NEIMO) 31-33 for internal coordinate dynamics (e.g., treats torsions only). This allows the solution of the dynamical equations for internal coordinates without inverting the mass tensor (moment of inertia tensor). The cost of NEIMO is linear in the number of degrees of freedom and small compared to other costs for million atom systems. More recently we also developed a new constrained force algorithm (CFA) for massively parallel MD simulation of polymers and dendrimer.34,35
  4. Advanced MD algorithms to simulate systems under constant temperature and constant pressure conditions.36, 37
  5. Nonequilibrium MD: We have implemented synthetic equations of motions to simulate various nonequilibrium conditions to predict transport properties such as viscosity, thermal conductivity of materials.38, 39 Using these methods, we have studied the effect of molecular topology of liquid alkanes on their measured viscosity indices40,41
  6. Steady State MD Methods are used to simulate non-equilibrium processes such as friction and wear in nanoscopically confined lubricants and diamond surfaces. Here, the external work is dissipated through the material and coupled to a thermal bath using the Langevin equation.42- 45 This method is specifically developed for studying tribological properties of materials.



3. Applications

3.1 Friction and wear properties of diamond

Microelectromechanical Systems (MEMS) is an emerging enabling technology that merges advances in information processing, storage, and display with advances in sensors and actuators to bring about a revolution in the way we perceive and control the environment 46-52. MEMS are natural progression in the capabilities of semiconductor devices. The ability of MEMS to gather and process information, decide on a course of action, and control the environment through actuators increases the affordability, functionality, and number of smart systems. A very significant role is attributed to MEMS for the technologies of 21st century, as applied to the following areas: Optical switches and aligners, parts handling in manufacturing, inertial navigation, pumps and valves, data storage, computer display, structural control of aircrafts.

The design of MEMS presents challenges at various length and time scales. The length scale ranges from Angstroms to millimeters where theories such as quantum mechanics, classical mechanics, solid (continuum) mechanics may be employed to improve practical performance.

Most of the MEMS devices explored and designed to date are based on Silicon, due to technological know-how accumulated on this particular semiconductor on manipulating, machining, manufacturing etc. on the microscale. Recent tribometric experiments carried out at Hughes by Gardos 53-56 strengthen the hypothesis in that the magnitude of adhesion and thus the adhesive friction between silicon and diamond surfaces are essentially defined by the number of dangling (high-friction), reconstructed (reduced-friction) or adsorbate-passivated (low-friction) surface bonds. Incipient linking of the sliding counterfaces by unsaturated (dangling) bonds on heating and their passivation by benign, chemi-sorbed gasses (e.g., hydrogen) were suggested as the primary causes of radically increased and reduced adhesion and friction, respectively. Hence to a performance degradation in Silicon based MEMS through a higher sliding friction and a resulting short wear life for devices with moving parts.

With this observation one must look into alternative materials, such as single crystal and polycrystalline diamond where sliding friction measurements lead to considerably lower values. As previously mentioned the wear rate of PCD is about four orders of magnitude lower than those of variety of Si crystallinites (poly-Si, Si(100) and Si(111)). In designing the next generation MEMS using PCD, there is a lot to be done. The hypothesis brought forth by Gardos to explain his experimental results on the basis of the chemical state of the surface bonds needs to be explored and rationalized theoretically.

Replacing conventional devices with their MEMS/NEMS analogs presents both opportunities and challenges. The evolution of small machines and sensors demonstrates that ever-increasing number of these devices with electronic circuits will yield a window to the world of motion, sound, heat and other physical forces. The current MEMS devices contain about 10 components MEMS devices. The anticipated applications such as pumps and valves require 103 - 105 components, optical aligners may have 102 - 104 components. In computer displays this range goes up to 106 or 107. The design, processing and fabrication methodology for such applications require research not only in the materials area but also on the operational conditions of such multi components systems, where simulation codes developed for robotics application will play a critical role33-34.


3.1.1 Problem of Friction and Wear

Friction is an unavoidable, parasitic force in flea-powered MEMS-MMA's. It is a crucial factor that determines not only efficiency but durability. Gardos' SEM-tribometric data in 53-58 indicated that high adhesive friction within the Hertzian contact zones of both Silicon and PCD causes shear-induced microcracking in the wake of the rolling/sliding contact. The ensuing wear and tolerance losses alone will fail these devices, especially when they are fabricated from high friction and high wear rate bearing materials such as Silicon. Although the hardness and chemical stability of PCD renders it a difficult material to shape into MEMS-MMA components and assembly of these components into complex micromachinery is even more difficult than those made of Silicon, the properties of PCD are far better for MEMS-MMA applications. Diamond is the hardest material known in nature, with a surprisingly high strength and fracture toughness in its polycrystalline form. Its thermal conductivity is high, readily dissipating both frictional and microelectronics-generated heat from functional areas. Due to the higher strength of the covalent carbon-hydrogen bonds, as compared to that of the silicon-hydrogen bonds, hydride-passivated PCD surfaces remain lubricative at higher operating temperatures.

3.1.2 Theoretical Background

Molecular dynamics have been widely used to study structures and properties of variety substance. It provides the information on atomic level directly. Therefore, people can test their theories, explain the experimental results on the bases of atomic motion. In molecular dynamics simulation, one of the most important steps is to apply a proper empirical potential to the system of interest. For a system with a given initial condition, its subsequent dynamics are completely determined by the potential surface that it evolves on. Choosing a proper classical potential is crucial to reproduce the real physical behavior of the system. In our study, we employ a modified bond order hydrocarbon potential 61-63. It allows chemical reactions during dynamics, therefore, it is a natural candidate to simulate both wearless friction and general friction.

3.1.3 Simulation Results

With the modified bond order potential, we study the friction process on reconstructed diamond 100 surfaces. Both bare and hydrogenated surfaces are investigated. In our simulation, two diamond crystals are put in contact. The bottom two layers of the lower block are held still, a constant velocity is assigned to the top two layers of the upper block. The external forces that are required to maintain the steady motion are recorded. As a result, the ratio between the external driving force and the normal force gives the friction coefficient. In order to prevent the system from being heated up, a stochastic thermal bath has been applied to maintain the average kinetic energy.

Figure 2 The top view of Hydrogenated 100 surfaces of two diamond crystal. Only the atoms at the interface are shown here. The gray and blue atoms are carbon atoms of moving block and fixed block, respectively. The red and green atoms are hydrogen atoms of upper and lower surfaces, respectively.



In Figure 2, we show the top view of two hydrogenated diamond 100 surfaces. Both surfaces are C 2 x 1 reconstructed. However, they are placed 90o respect to each other. The atoms on each surface create the surface potential with mountains and valleys. One can expect that different friction forces might arise from different movement on the microscopic level. The general experimental friction coefficients are the average value of many different microscopic results. In particular, the friction of poly crystalline diamond (PCD) is an average quantity of friction along many different crystal orientations. In our simulation, we are able to extract each individual friction coefficient for each direction on a specific surface. We call it differential friction coefficient. For the system shown in Fig. 2, we have calculated the differential friction for three different sliding directions. The three characteristic directions are x-direction, y-direction, and xy-direction. Obviously, sliding along xy-direction will have lowest resistance, because the path is lying along a potential valley. On the other hand, moving in the x-direction and y-direction has to climb up the potential barriers created by the presence of hydrogen atoms along the sliding direction.

Figure 3. Running average of normal forces. The force from x-direction movement is shown as solid line, the force from y-direction movement is shown as dash line, and the force from xy-direction movement is shown as dotted line.



In Fig. 3, we show the running averages of forces in the normal direction respect to the movement. The initial oscillations represent the process of reaching a steady state motion. Usually, the system arrives at the steady state within 10 ps. The constant velocity is maintained at 1 A/ps (100 m/s). Fig. 3 illustrates that the normal forces are very similar despite different moving directions. This result indicates that any difference in the differential friction coefficients can only come from different motion arrangements in the regime of our investigation. Our normal force is around 850 kcal/A/mol. Fig. 4 depicts the running average of the friction coefficients calculated for each sliding direction. As we expected, xy-direction has lowest friction coefficient, and x-direction and y-direction have high friction coefficients. If the two surfaces are perfectly aligned, x-direction and y-direction would have same friction coefficient. However, in our simulation, the model set up such that the top block is shifted in y-direction, to create a slightly different potential energy surface. In turn, this gives rise to smaller friction coefficient along x-direction.

Figure 4. Running average of differential friction coefficients. The solid line shows x-direction friction coefficient, the dash line shows y-direction friction coefficient, and the dotted line shows xy-direction friction coefficient.




One phenomenon that always accompanied with friction is energy dissipation. The friction coefficient can also be extracted from energy consumption curve. In Fig. 5, we plot the external work done to the system as a function of time. The slope of each curve gives the friction coefficient, and the detailed features of each curve reflects the potential surface contour along the each slide path. For instance, the three peaks in each period of xy-direction curve represent each hydrogen atom on the lower surface passing through three reconstructed C--C bonds of the upper surface.

Figure 5. External work done to the system as a function of time. The solid line is calculated from x-direction motion, the dash line is calculated from x-direction motion, the dash line is calculated from y-direction motion, and the dotted line is calculated from xy-direction motion.




In general, clean crystal surfaces are less stable because of the dangling bonds. This is also the case for diamond. The pi-chain reconstructed diamond (100) surface is known to be energetically favorable. However, when two such surfaces are brought together, they can easily react, which increases the friction force more than 2 orders of magnitude. In this case, frictional force comes from not only nonbond interaction, but chemical binding forces. This is the main cause of large frictional force. Fig. 6 shows the running average of the friction coefficient and the external work done on the system. Compared to the hydrogenated surface, the friction coefficient is much larger. Our simulation reveals the cost of formation and breaking of chemical bonds on the interface during structures of our simulation.

Figure 6 Running average of a differential friction coefficient and external work as a function of time.





3.1.4 Conclusions and Discussions

Friction and wear are complicated physical and chemical processes. Our investigation mainly focused on atomic scale origin of friction. For crystal surfaces, sliding in different directions gives different friction coefficients. The surface potential energy landscape (which is periodic for single crystal diamond) is the main source of this. Also, dangling bonds on surfaces can create large frictional force. In this case, it is necessary to break these bond in order to have relative movement. Our calculations clearly show how these aspects contribute to the atomic friction and wear. In all cases, surface passivation has to be done to have small friction and wear rate. If possible, a specific sliding direction is preferred. Another way of decreasing friction forces is to add proper lubricants to the interface. Gardos[11] has reported that water is a perfect lubricant for both diamond and silicon. It also needs to be noted that there are other factors giving rise to overall frictional force. For instance, surface asperities and contact area can dramatically modify total friction coefficients.

3.2 Structural and Dynamic Properties of Lubricants Under Shear Flow in a Confined Geometry

3.2.1 Introduction

Advances in experimental and computer simulation techniques have led to an increasing understanding of fluids in ultra-thin films at the atomic scale.65, 66 It has been found that when a liquid, confined by solid surfaces, as the separation becomes progressively narrower, its physical properties change at first gradually then more dramatically.67, 68 In this study we investigate the dynamics of nanoscale hexadecane film under the influence of a shear flow. The primary objective of the study is to enhance our understanding on how lubricants behave when they past engine surfaces and to assist in designing engine surfaces with reduced engine wear.69-71

The model of the present study builds upon our previous static studies of iron-oxide surfaces covered by a self-assembled monolayer (SAM) of adsorbed wear inhibitors (``brush''). 72 A lubricant film, namely liquid hexadecane, is sandwiched between two such protected iron-oxide surfaces. A planar flow is generated by moving the top surface and keeping the bottom surface stationary. Our model system consists of realistic lubricants, realistic solid surfaces, and realistic wear inhibitor molecules adsorbed on the surfaces.

The wear inhibitor molecules are dithiophosphates (DTP= S2P(OR)2, where R represents one of the three types of organic groups: R=isopropyl (iPr), isobutyl (iBu), and phenyl (Ph)).

3.2.2 Simulation Model

The details of simulations performed in this study are described in our recent paper, 73 and a typical structure of the system is shown in Figure 7. In brief, the fluid consisting of hexadecane molecules was confined between two iron-oxide surfaces, each of which is covered by 8 dithiophosphate (DTP) molecules (acting as wear inhibitors). The Fe2O3 surface is oriented in such a way that the (0001) surface is perpendicular to z-axis. A planar flow is generated in xz-plane by moving the top surface at a constant velocity along the x-axis while keeping the bottom surface fixed.

The initial velocities were chosen from a Maxwell-Boltzmann distribution at T = 500 K. We added a constant velocity, v, to all atoms of the top slab of Fe2O3 and to the top SAM layer. The lubricant and bottom layer started with the initial set of velocities from the Maxwell-Boltzmann distribution with zero center of mass velocity. In the simulations, the iron-oxide surfaces are sheared uniformly with respect to each other at a constant velocity and the DTP plus lubricant are allowed readjusting to the shearing motion until they attain a steady state shear velocity profile. We carried out ~200 ps of shear dynamics simulations. The time step is 1 fs. The trajectories were stored every 0.05 ps.

Figure 7. The structure and setup of the system for the shear dynamics simulations. Shown here is the snapshot after 200 ps for the case of iPr. The values for H, L, etc. are in Table 1. For the middle region, the 32 C16 molecules are partitioned into 8 different layers according to their zcm; each layer is plotted with a different color.




The force field (FF) used in these calculation was obtained by combining the QM calculations on Fe-DTP clusters with the Dreiding force field74 to describe the DTP and iron oxide interactions. The derivation and the parameters for this FF were described elsewhere.71-72

The results presented below are from four different simulation runs:

1) R = iPr, v = 1 A/ps, and H ~45 A;

2) R = iBu, v = 1 A/ps, and H ~45 A;

3) R = Ph, v = 1 A/ps, and H ~43 A;

4) R = iPr, v = 0.5 A}/ps, and H ~20 A.

A list of parameters of these simulations is given in Table 1. For each system the cell length in the z direction, L, was determined by minimizing the energy of the entire system (with fixed cell parameters in C16H34 layer is ~ 20 A.

Table 1. Values of the geometry parameters shown in Figure 7 for various SAMs.


La (A)

Hb (A)

Dc (A)

h (A)d

Runtime (ps)

(a) 44 A Film



















(b) 20 A Film







a. Total length of periodic cell in z direction

b. Thickness of the lubricant region.

c. Thickness of the ``brush'' wall consisting of Fe2O3 and DTP monolayers.

d. Thickness of wear inhibitor SAM.

3.2.3 Results and Discussions

3.2.3.a Density Oscillations

Figure 8 shows the distribution of the density of lubricants along z-axis from the four simulations. We observed distinct density oscillations, which is consistent with the previous work.75 The density oscillations are quite significant near the walls (within 15 to 20 A). In the central region of the fluid, the layering decreases, with the density approaching to that of the bulk lubricant. For the separations of 44 A we find 9-10 layers in the density. For the thinner system with just ~ 20 A separation (Figure 8d), we found 5 layers, but here the density oscillations are easily noticeable even in the central region of the film. In an earlier study73 we have looked at the distribution of carbon atoms and that of hydrogen atoms separately. There we found that the distribution of hydrogen atoms is quite uniform and that the density oscillations are mostly resulted from the oscillation in the number density of carbon atoms.

Figure 8. The density distribution of the lubricant as a function of z (perpendicular to the surface). We used 1000 bins to count the number of atoms falling in a zone between z and z+dz at each time-step. This density was averaged over the last 100 ps of the 200 ps runs: (a) iPr, H ~ 45 A; (b) iBu, H ~ 45 A; (c) Ph, H ~ 43 A; and (d) iPr, H ~ 20 A.




Here we further analyzed the contributions from each individual molecule to the oscillatory distribution of total carbon atoms in the lubricant film. The results are shown in Figure 9 [(a) for iPr, (b) for iBu, and (c) for Ph], where the carbon atom distribution for each of the 32 molecules (the dark solid curves) is plotted against the distribution of total carbon atoms (the light solid curve). The results clearly indicate the layering occurring in the lubricant film. In other words, the motions of the molecules in the z-direction are mostly confined within 1-2 layers during the entire 200~ps dynamics. Only occasionally a molecule can move across three layers. The molecules closest to the bottom and top interfaces are the most confined, due to the interactions with the ``brush'' surfaces. Their motions are restrained to one layer only.

Figure 9 The distribution of carbon atoms for each lubricant molecule (dark solid curve). The background (light solid) curve is the distribution of total carbon atoms of all molecules. These are results averaged over last 100~ps. (a) For iPr, (b) for iBu, and (c) for Ph.





For a comparison, we did a (static) simulation of 140~ps at zero shear velocity for the case of iPr with the exactly same setup.73 There we still observed the strong density oscillations. From that, we concluded that the density oscillations observed in the lubricant films are induced by geometric confinement, rather than shearing.

3.2.3.b Stick-Slip Motions

Figure 10 shows the displacement of the center-of-mass of each lubricant molecule, along the shear direction (x-axis) as a function of time. In the figures the diagonal black line and the horizontal black line represent the displacements of the top and bottom SAM molecules, which clearly indicate a velocity of 1 A/ps for the top SAM and a zero velocity for the bottom SAM. To distinguish the motion of each molecule, we partitioned the 32 lubricant molecules into 8 layers according to the similarity of their motion in the x direction. Each layer is indicated with a distinct color. The label for each molecule in the figures is consistent with the numbering of the molecule shown in Figure 9.

Figure 10. The displacement of each lubricant molecule in the shearing (x) direction as a function of time, for iPr, iBu, and Ph. For easy viewing, the 32 lubricant molecules are partitioned into 8 layers by the final positions in x; each layer is indicated with a distinct color. The label for each molecule is in agreement with number of molecule in Figure~3.






We can see stick-slip motion in the molecules near the bottom interface in all three cases. The stick-slip motion, often observed in macroscopic studies on solid surfaces, 68,76 is characterized by a sequence of move (slip), pause (stick), and move (slip) again. The case of iBu shows the most regular character, with a periodicity of ~ 40 ps whereas the stick-slip motions for the iPr and Ph cases are more irregular.

Normally, in observing the stick-slip motion, the experimental setup is to place a solid sliding block on top of a solid surface and drive the top block elastically. Here, in our simulation, the top ``brush'' surface was driven smoothly at a constant velocity, and we observed the stick-slip motion in the bottom layers of the lubricant. In this case the stick-slip motion of the interfacial boundary propagates into the lubricant. The result indicates that the stick-slip phenomena is rather universal, observable on the nano- or molecular-scales as well.

What is the origin of the lubricant stick-slip motion? In the case of stick-slip sliding on solid surfaces, there is a one-to-one correspondence between the stick-slip motion and the frictional force or shear stress (the ratio of the frictional force to the contact area). The friction is highest during the stick period (referred as static friction) and lowest during the slip period (referred as kinetic or dynamic friction). We have evaluated the frictional forces in our systems, but unfortunately, the fluctuations of the friction forces themselves are so large that it is impossible to differentiate them from the difference between the static friction and kinetic friction.9 However, the fluctuations of frictional forces in our system are too large to be explained in terms of pure thermal fluctuations. We think that the frictional forces are in some way responsible for the observed stick-slip motion in the lubricant.

Since the flow behavior of the lubricant changes dramatically with the nature of the wear inhibitor molecules at the interface, we concluded that the properties of the SAM molecules must play a role. Since the SAM molecules stay at the same relative position on the iron oxide surface, they probably play the role on modulating (elastically) the frictional forces on the lubricant. Thus, although the top wall is moving at a constant velocity, the wear inhibitor molecules lead to elastic response causing an oscillatory driving forces to the lubricant. That is, the lubricant might compress or push the SAM molecules, which then rebounds to press on the neighboring lubricant molecules.

3.2.3.C Structures at Various Stick States

To examine the structure of the lubricant molecules when they are stuck, we took the snapshots of molecule `2' for the iPr case, molecule `1' for the iBu case, and molecule `1' for the Ph case at five different times, each of which corresponds to a stick state. Shown in Figure 11a are the structures of molecule `2' in iPr at 45 ps (white), 80 ps (light brown), 105 ps (yellow), 145 ps (orange), and 190 ps (dark red); in Figure 5b are those of molecule `1' in iBu at 20 ps (white), 70 ps (light brown), 130 ps (yellow), 150 ps (orange), and 185 ps (dark red); and in Figure 511 c are those of molecule `1' in Ph at 20 ps (white), 70 ps (light brown), 105 ps (yellow), 160 ps (orange), and 190 ps (dark red).


Figure 11. Snap shots of lubricant molecules at various stick-states. Top -- top view (looking down along the z-axis); middle -- side view along the $y$-axis (perpendicular to the shearing direction); bottom -- side view of along the $a$-axis of hexagonal unit cell of iron oxide. (a) Molecule `2' in the iPr case, (b) molecule `1' in the iBu case, and (c) molecule `1' in the case Ph. Times at which the snap shots are taken are indicated at the top.





The top panels are the top view (looking down along the z-axis), the middle ones are the side view parallel to the y-axis (perpendicular to the shearing direction), and the bottom ones are the side view along the a-axis of the hexagonal unit cell of iron oxide.

From Figure 11a, we can see that this lubricant molecule did not move much during the 200 ps dynamics (~20 A). Indeed, the molecule is mostly stuck to the SAM interface for the last 80 ps. We think that this is because the iPr is rather rigid and interactions between iPr and the lubricant molecule are strong. On the other hand, we can see small "hills" and "valleys" on the surface of R=iBu case formed by the slightly longer hydrocarbon chains of iBu's. It is this surface roughness that causes the stick-slip motion of the lubricant molecule the most distinct in the iBu case. The slip states in this case is also more profound than the other cases. It appears that the stick takes place when the molecule is in the "valley" between the two hydrocarbon chains extended from the dithiophosphate branches. In the Ph case, we found the surface is quite uniformly covered by the phenyl groups, as shown in the bottom panel of Figure 11c. The stick-slip motion of the molecule in this case is the least manifested. The period for molecule changing from slip to stick is short. We may speculate that as the surface becomes smoother, the stick-slip motion will become less profound and may eventually disappear.

3.2.3.D Radius of Gyration, End-To-End Distance, and Torsional Conformation Distribution

We also calculated the radius of gyration (R0), end-to-end distance (L), and number of trans-conformers (Nt) for each lubricant molecule. The results are presented in Figure 12a for iPr, Figure 12b for iBu, and Figure 12c for Ph. We see strong correlations between these properties in all three cases. When the molecule is stretched, the radius of gyration increases, so do the end-to-end distance and the number of trans-conformers in the carbon chains.


Figure 12. The radius of gyration (top), end-to-end distance (middle), and number of {\it trans}-conformers (bottom) for each lubricant molecule as a function of time. (a) For iPr, (b) for iBu, and (c) for Ph The numbering of the molecules is the same as in Figure 3.



However, by comparing Figure 10 with Figure 12, we found no correlation between these properties and the stick-slip motion of the molecules. Apparently, when the molecule gets stuck, it is equally possible for it to be stretched or to be recoiled. So, changing from slip to stick state, the end-to-end distance of the molecule can either increase or decrease. In addition, there are situations where the molecule appears in Stick State in the shearing direction but it actually moves in one of the other two perpendicular directions. In these circumstances, the information on the motion in the shearing direction along (which defines the stick-slip motion) cannot determine whether the molecule has been stretched or shrunk. Thus, we did not observe a one-to-one relationship between the stick-slip motion and the variation of the radius of gyration or the end-to-end distance of the molecules.

Figure 13. The average radius of gyration (left), end-to-end distance (center), and number of {\it trans}-conformers (right) for each molecule as a function of its average position in z. (a) for iPr, (b) for iBu, and (c) for Ph.




Furthermore, we did not observe systematic trends in radius of gyration, end-to-end distance, and number of trans-conformers across the lubricant films. Intuitively, one may think that the top molecules may be stretched out more than the bottom molecules since they have higher shear velocities. In Figure 13, we plotted the radius of gyration, end-to-end distance, and number of trans-conformers for each molecule averaged over the last 100 ps of the simulations against the average z-position of the center-of-mass of the molecule. Here, besides the fluctuations, there are no systematic changes in these quantities across the lubricant films (i.e., along the z-axis).

3.2.4 Remarks

With the advances of computer simulation techniques we are able to simulate rather realistic and complex systems. Here we simulated nanoscale hexadecane lubricant films confined between two ``brush'' surfaces consisting of iron-oxide surface and adsorbed wear inhibitor molecules. We found that in such a confined geometry, the lubricant films exhibit strong density oscillation. For the separations of 44 A we found 9-10 layers in the density, and for a 20 A separation we found 5 layers. In addition, we found that the motion of each individual molecule in the direction perpendicular to the surfaces is very limited. They are confined within 1-2 layers.

Under a shear flow, the lubricant molecules near the bottom surface boundary show stick-slip motion in the shear direction. However, the stick-slip motion cannot be simply characterized by the stretching or compressing of molecule. The characteristics of stick-slip motion is very sensitive to the type of ``brush'' molecules on the surfaces that are in direct interactions with the lubricant molecules. The stick-slip motion should disappear on an ideally smoothened surface.

3.3 Carbon Nanotube as a Tribological Probe

Over the past decade since their discovery, carbon nanotubes have been one of the focal points of both scientific research and technological developments. This is mainly because of their natural beauty of structures and unique properties. Large quantities of carbon nanotubes are now available commercially, although it is still a challenge to synthesize and pack nanotube into single crystal form. Carbon nanotubesí electronic and mechanical properties are two major areas of current research interests. Single transistor has been built based on a single walled carbon nanotube, and electronic conductance of single walled and multi walled carbon nanotube have also been measured. In addition, it is also demonstrated that carbon nanotubes can be used as field emission materials. Besides remarkable electronic properties, very unique mechanical properties of carbon nanotubes have been predicted by various theoretical methods and measured by experiments. It is predicted that tensile strength of carbon nanotubes is about 100 times stronger than steel and the bending moduli are also very high. Recently, nanotubes have been used in various Atomic Force Microscope (AFM) experiments and shown great advantages of their fine structures, mechanical strength, and stability/reversibility.

Recent progress in tribology studies demonstrated that ancient view of friction and wears needed to be greatly modified at the microscopic level, where detailed surface structures at the nanometer scale are very important. At the microscopic level, the macroscopic rule of friction that is the result of statistical averaging does not hold. The wear itself by natural is the event happening at microscopic level (i.e. chemical reactions). In order to address the relationship between the fine surface structure and the nanoscale friction behavior. It is important to design new experiments that can simultaneously probe the surface structure and the atomic friction at the nanometer scale. Based on the unique properties and structures, carbon nanotubes can serve as ideal candidates to achieve our goal.

Figure 14 Configuration of functionalized nanotube tip and active site of diamond reconstructed surface (100).



In this section, we carry out some prototype simulations to illustrate the potential applications of nanotube tip as a nanotribology probe. One nice property of carbon nanotube is that the tube tip can be functionalized by various active species, so that it can be used to initiate surface reactions at a precise position. By measuring the force/energy on the nanotube, we can probe the specific reaction between the tube tip and the surface. In Fig. 14, we show a very simple example in which the nanotube tip is functionalized by a hydrogen atom and intentionally left with an unpaired electron. Similarly, the reconstructed diamond surface (100) has one hydrogen atom attached and one dangling bond. By pushing down the nanotube tip, we can expect some chemical reactions between the tip and the diamond surface active site since both tube tip and surface would like to be passivated. Three possible reaction pathways here are: 1) hydrogen extraction from the nanotube tip, 2) hydrogen extraction from the diamond surface, and 3) pairing two dangling bonds between tip and the surface. Our simulations show that the third pathway is the most probable reaction since it is the most energetically favored route on the deformed tip and surface. In Fig. 15, we show that the chemical bond (s - bond) is formed between the diamond surface and the tip.

Figure 15 Chemical bond formation between the tip and the surface.




In Fig. 16, we show the external energy consumption during the loading and unloading processes. As we can see from the figure, the unloading curve has a clear turnover that indicates the formation of the new chemical bond. By measuring the slope of the curve, the bonding strength can be obtained.

Figure 16 External energy consumed during the loading and unloading process.



In addition to monitoring and probing chemical reactions on atomic surfaces. Carbon nanotube can also be used to probe the frictional behavior at a very fine resolution without introduce chemical reactions. In Fig. 17, a specially constructed diamond rough surface is used to test the sensitivity of the nanotube "tribology machine". As we can see in the picture, two atomic layers are removed periodically to resemble the asperity at the atomic scale. When the nanotube is moved across the surface, it will feel the stronger or weaker interactions due to hills and valleys on the surface. Therefore, it can simultaneously measure the forces both in the normal direction and lateral directions. From the normal direction, it can be used to extract the information of the surface pattern, and its relationship with lateral forces, usually referred as frictional forces, can be drawn subsequently. In Fig. 18 and Fig. 19, we show some snap shots of the configuration during the process of dynamic movements.

Figure 17 Nanotube tip and rough diamond surface



Figure 18 Nanotube tip encounters the hill on diamond surface.




Figure 19 The upper panel shows the nanotube tip at the top of a surface hill. The Lower panel shows the nanotube tip at the top of a surface valley.





In Fig. 20 and Fig. 21, normal force and lateral frictional forces are depicted. From the normal force envelope, the surface roughness is clearly seen, and the shape of the surface roughness is also shown in the picture. The lateral forces show very clear stick-slip motion as we can see from the low pass filtered force envelopes. The use of Fourier low pass filter is to eliminate the high frequency oscillations that resemble the atomic vibrational motions.

Figure 20 Normal force acting on the nanotube tip. The black curve is the actual force, and the red curve is the Fourier filtered force



Figure 21 Lateral forces on the nanotube tip. The blue and black curves are actual forces, and the pink and yellow curves are Fourier filtered forces.



From the averaged lateral forces and normal forces, it is possible to deduce the macroscopic friction coefficient. Also, this quantity can be obtained by monitoring the energy dissipation or external work done to maintain the constant motion. In Fig. 22, we show the external energy consumed as a function of time during the process of nanotube movements. The energy curve not only illustrates the overall energy consumption due to friction between the tip and surface, but also shows the detailed information about surface roughness. From the linear fit of the energy curve, we can deduce the effective friction coefficient, given that the friction coefficient does not vary appreciably for different speed of tip movements. In our simulations, both friction coefficients calculated from the ratio between lateral force and normal force and from energy dissipation give consistent value, about 0.25. Of course, the accuracy of the value depends on the accuracy of the atomic interactions.

Figure 22 Energy dissipation due to friction between diamond surface and nanotube tip.



3.4 Reactions of Hydrogen on the surface of diamond

A problem in nanotechnology is the need to tie off surface dangling bonds with hydrogen atoms. One approach to this problem involves the use of tools to add H atoms one at a time. While this approach is useful for demonstration of the concept and could be used to construct small pieces of nanomachinery or possibly an assembler, a more efficient method would be useful in many cases. We have considered the possibility of reacting a bare surface with gas phase H2/H as a means to passify the surface. In this study, we have considered the diamond 100, 110, and 111 surfaces.

Cluster models for the three surfaces of diamond are shown in Fig. 23. The unrelaxed 100 surface has carbene like surface carbon atoms; however, for the relaxed surface these dimerize to give rows of surface dimers and there is a considerable amount of p bonding between the radical orbitals of the dimer. The 110 surface has zigzag rows of carbon atoms with a dangling bond on each carbon atom. These dangling bonds are hybridized away from each other and thus interact less strongly than for the 100 surface. Finally, the 111 surface has surface carbon atoms arranged in a triangular pattern and the surface dangling bonds are well separated from each other (second nearest neighbor distance) leading to almost no interaction between adjacent dangling bonds. These qualitative features may be quantified by computing the overlap of adjacent dangling bonds in a GVB(pp) calculation (Bobrowicz and Goddard, 1977). The overlaps are 0.462, 0.292, and 0.016 for the diamond 100, 110, and 111 surfaces, respectively, and 0.322 for the silicon 100 surface. (Note that these are with a 6-31G basis set and all the overlaps would be increased for a larger basis set, but these results clearly indicate the trends.) As discussed elsewhere (Walch and Merkle, 1998), the different overlaps between adjacent radical orbitals result in different reactivity for the various surfaces. For example, in the reaction with a carbene, the diamond 111 surface behaves like a radical, but the diamond 100 surface behaves like a p bond. Based on the GVB(pp) overlaps, the 110 surface is expected to be somewhere in between the 100 and 111 surfaces in reactivity.

Figure 23 Cluster models for the 100, 110, and 111 surfaces of diamond.



Fig. 24 shows some reactions of H2/H with the reconstructed diamond (100) surface. Reaction (1) is a 2+2 cycloaddition and is expected to have a barrier. The size of the barrier will decrease as the amount of p bonding in the surface dimer decreases. However, this process is expected to have a significant barrier for the diamond 100 and 110 surfaces. Reaction (2a) is addition of H2 to a surface carbene atom. This process may have a low barrier; however, it requires activating the surface to create the reactive site. Reaction (2b) is one reaction which may activate the surface by creating a surface carbene. Reaction (3) is a hydrogen abstraction process. This process requires that the dimer have some biradical character. Reaction (4) is a reaction of atomic H with the surface. This process will have only a small barrier but requires that the H2 be dissociated (as e.g. by reaction with a hot tungsten filament).

Figure 24 Some reactions between H2/H and diamond (100) considered in this work.




A review of diamond surface chemistry is given in (Wei and Yates, 1995) and a discussion of H/silicon interactions is given in (Murty and Atwater, 1995).

3.4.1 Calculational Details.

The density functional theory (DFT) calculations used the B3LYP functional (Becke, 1993) and the 6-31G basis set (Ditchfield, Hehre, and Pople, 1971). the calculations were carried out with Gaussian94 (M. J. Frisch et al., 1995). In the case of Si, the Si1s22s22p6 core was replaced by an effective core potential (Stephens/Basch/Krauss ECP split valence-CEP-31G) (Stevens, Basch, and Krauss, 1984 ). The smallest basis set used here (denoted by small) is equivalent to the 6-31G basis set that was used for diamond. Because polarization functions are expected to be more important for silicon than for carbon, some larger basis sets were used for the silicon clusters. The medium basis set adds a set of d functions on silicon, and the large basis set also adds a set of p functions on H.

The calculations for ethylene + H2 were ab initio calculations. The stationary points were located with complete-active-space-self-consistent-field (CASSCF)/derivative methods (as implemented in SIRIUS/ABACUS) using the cc-pVDZ basis set (Dunning, 1989) and the energetics were determined with internally contracted CI (ICCI) (Werner and Knowles, 1988) (Knowles and Werner,1988) using the cc-pVTZ basis set (Dunning, 1989).

3.4.2 Discussion.

3.4.2.A Diamond 100 Surface.

a. Reaction (1).

Attempts to study reaction (1) using DFT were not initially successful. It is clear that reaction (1) requires a complete-active-space-self-consistent field (CASSCF) zero-order description since it involves a crossing of two configurations which differ by a double excitation. In order to gain insight into the problem, the gas phase reaction of ethylene with H2 was studied first using an ab initio approach. The starting point for these calculations was a CASSCF calculation with 4 electrons in 4 obitals; this was followed by an internally-contracted-configuration-interaction (ICCI) calculation with the cc-pVTZ basis set. The active orbitals correspond to the CC p bond in ethylene and the HH bond. Fig. 3 shows the results of this calculation. The process of adding H2 to ethylene in a C2v constrained geometry is seen to have a barrier in excess of 0.69 aJ (100 kcal/mol). This is a second order saddle point and allowing the symmetry to break leads to the other saddle point, which is basically an H abstraction process. (Note that the saddle point is above ethyl radical + H at the CASSCF level but ICCI places the saddle point below ethyl radical + H and this process probably occurs with no barrier.)

In order to model the analagous reaction with a dimer on the diamond (100) surface, a cluster was used which was a distorted ethylene with the CC distance fixed at the distance obtained for the singlet state in the DFT calculations. In one case the distorted ethylene was planar and in the other case the distorted ethylene had a tetrahedral configuration around the two carbons. The saddle points for H2 addition in a C2v constrained geometry were obtained from a CASSCF grid and the energetics were obtained from ICCI calculations as for the ethylene + H2. This gave barrier heights of 105.7 kcal/mol for the planar cluster and 73.7 kcal/mol for the tetrahedral cluster. As expected the tetrahedral cluster weakens the p bonding and lowers the barrier to addition. However, it is clear that this process will be unfavorable for a dimer on the diamond (100) surface. Table I shows energetics for the same process obtained with the DFT method. Here the distance above the surface was fixed and the surface dimer and H2 were allowed to relax at each distance in order to generate a C2v symmetry constrained minimum energy path. From Table I it is seen that DFT predicts a barrier of about 50 kcal/mol for reaction (1).

b. Reaction (2).

Reaction (2) should have a low barrier since it is a carbene insertion process. As a first model a small cluster (See Fig. 4) was used which had one surface carbon and two second layer carbon atoms. (With the second layer carbon atoms fixed at bulk values and tied off with H atoms.) The barrier for reaction (2) was found to be 1.4 kcal/mol and the exothermicity for this reaction is 113.3 kcal/mol. (See Table II.)

Several other processes were also studied for this cluster. The inversion barrier for a H atom bonded to the carbene atom is less than 0.1 kcal/mol. The HCC2 moiety is only 10 deg. from planar, which is reasonable by analogy to CH3, and consistent with the small computed barrier. The barrier to abstraction of an H atom from a surface CH2 by H atom was found to be to be 9.6 kcal/mol and the heat of reaction is 4.4 kcal/mol (endothermic). The final question was the barrier for adding H atom to a surface CH. No barrier was found for this ( doublet + doublet reaction) and the bond strength is found to be 105.4 kcal/mol.


Fig. 5 and Table III show energetics for reaction (2b). (The cluster is shown in Fig. 1.) This process is assuming an entrance channel saddle point that leads to the structure labeled by min2 in Table III. This structure has two H atoms bonded to one C and the other carbon is a carbene site. An interesting feature of Fig. 5 is the low barrier for migration of a H onto a carbene center. This process is well known in gas phase chemistry (e.g. the isomerization of vinylidene to acetylene). (Walch and Taylor, 1995)

Another issue which was considered is the formation of a dihydride phase (i.e. two H atoms per surface C atom).

The structures which were considered here are shown in Fig. 7 and the energetics are shown in Table IVa. One question that hadn't been considered before is the structure of the two minima for three and four H atoms on the 100 surface dimer. One speculation was that these structures would twist out-of-plane to minimize repulsions. However, calculations that were started with twisted geometries reverted to planar geometries when optimized. From the structure of the 4 H atom minimum it is clear that there would be significant interactions with an adjacent dimer site with 4 H atoms. Thus, for the extended surface with 2H per surface carbon there would be large non-bonded repulsions with the CH bonds of adjacent monohydrogenated dimers. The main effects of this non-bonded repulsion were included by placing CH2 groups in the locations of the nearer carbon of the adjacent surface dimers. (Here one CH bond is oriented along the CC bond of the adjacent dimer and the other CH bond is in the location of the CH bond of the adjacent surface dimer. See Fig. 6.) This is similar to what was done earlier in our studies of the monofluorinated silicon surface (Walch and Halicioglu, to be published). Table IVa gives computed energetics for the case without adjacent CH2 groups. Here it is seen that the binding energy for a third H atom to the mono hydrogenated dimer is 50.3 kcal/mol or about half a CH bond strength. This is due to disruption of the CC p bond of the dimer. The barrier to H migration is only 13.7 kcal/mol, which is much less than the barrier of 63.8 kcal/mol for H atom migration of a single H atom on a surface dimer. Adding another H atom to the dimer plus three H atoms leads to a CH bond strength of 98.4 kcal/mol in the dimer plus 4 H atom cluster.


A number of other results are also summarized in Table IVa. These include the barrier for carbene insertion of H2 into the structure (2b) of Fig. 2. This barrier is 6.3 kcal/mol, which is 5 kcal/mol larger than for reaction (2a).

From Table IVb it is seen that the presence of adjacent CH2 groups reduces the binding energies of H atoms to a surface dimer. The reductions are 2.0 kcal/mol, 12.4 kcal/mol, and 18.1 kcal/mol for 2H, 3H, and 4H, respectively.

c. Reaction (3).

Table V and Fig. 8 show the reaction of H2 with a surface dimer along an abstraction pathway. From Fig. 8 it is seen that the abstraction reaction is nearly thermoneutral and has a small barrier. The calculations were carried out for the triplet spin state. Fig. 8 also shows a singlet pathway. This pathway is constructed assuming that the barrier is comparable on the singlet surface. Note that for the reverse reaction either singlet or triplet spin is possible, and the saddle point region on the singlet surface is assumed to be an open shell singlet. Here it is seen that on the singlet surface, the barrier including the endothermicity is about 15 kcal/mol.

Table VI and Fig. 9 show the process of adding a second H atom to the surface dimer via the abstraction pathway. This process occurs only on a doublet surface and has a small (about 6 kcal/mol ) barrier.

d. Reaction (4).

Table VII and Fig. 10 show energetics for adding an H atom to a dimer on diamond (100). DFT calculations show no barrier to addition of H to the singlet state of the dimer. Addition of a second H atom should also be a barrierless process. The exchange of an H atom between the two C atoms is also shown. This process has a barrier of 63.8 kcal/mol.

3.4.2.B. Diamond 110 Surface.

We now consider reactions on the diamond 110 surface. Fig. 11 shows the clusters that were used for this surface. This cluster has 4 surface and 6 second layer C atoms with the dangling bonds tied of with H atoms. Here the position of the 4 surface atoms was partially optimized (with the constraint that the three CC bond lengths remain the same). (Calculations in which the CC bond lengths were allowed to be different indicated a slight energetic preference for an alternating structure.) The center two surface atoms of this cluster will be considered as a dimer pair. The computed energetics are given in Table VIII. From Table VIII it is seen that the CH bond strength for adding one H atom to the dimer pair (cluster110b+h3) is 109.6 kcal/mol) and the exchange barrier is 72.0 kcal/mol as compared to 63.8 kcal/mol for the 100 surface. This larger barrier probably results from more hydridization of the radical orbitals away from each other as compared to the 100 surface. The overlap of adjacent singlet coupled radical orbitals is 0.292 in a GVB(pp) calculation and the singlet to triplet excitation energy is 7.3 kcal/mol for the diamond 110 surface, compared to 0.462 and 32 kcal/mol for the diamond 100 surface. The abstraction barriers were also computed for abstracting an H atom from the mono and dihydrogenated dimer (cluster110b+h2+h2.abs and cluster110b+h3+h2.abs, respectively). The abstraction barriers are 5.2 kcal/mol and 6.8 kcal/mol and the heats of reaction are 0.1 kcal/mol and 1.9 kcal/mol for the mono and dihydrogenated dimer, respectively.

Calculations were also carried out for a symmetry constrained addition of H2 across the dimer. These results are given in Table IX. Here it is seen that the barrier for the symmetry constrained approach is ~43 kcal/mol. For diamond 100 the barrier was 49 kcal/mol. This is consistent with the smaller overlap of the p bond for the 110 surface compared to the 100 surface.


3.4.2.C Diamond 111 Surface.

Calculations were also carried out for several processes on the diamond 111 surface. The cluster used here has two surface carbons and a second layer carbon. (See Fig. 12 and Table X) The two surface dangling bonds are separated by second nearest neighbor distances and interact only very weakly. This weak interaction is evident in the first and second CH bond energies which are 116.1 kcal/mol and 115.9 kcal/mol, respectively. Also these binding energies are significantly larger than on the 100 and 110 surfaces, which also is expected from the free radical character of the surface dangling bond on the 111 surface, compared to the other surfaces, where there is some p bonding. The barrier for abstraction is calculated to be 4.0 kcal/mol and the reaction is exothermic by 6.1 kcal/mol. The exothermic nature of the abstraction reaction is due to the greater CH bond strength for the 111 surface

3.4.3. Concludung remarks.

We have looked at the reaction of H2 with the diamond 100, 110, and 111 surfaces using a cluster model and DFT with the B3LYP functional. For each of the four surfaces we have studied the symmetry constrained and non-symmetry constrained addition of H2, the addition of H atom, and abstraction of an H atom from a surface CH by a gas phase H atom.

For the dimer reconstructed 100 surfaces we considered the reactions shown in Fig. 24. The symmetry constrained addition of H2 to the surface dimer (reaction (1)) is found to have barriers of 50 kcal/mol for diamond 100 surface. respectively.

The non-symmetry constrained addition of H2 to the surface dimer (reaction (2)) differs for diamond and silicon. In the diamond case, one possibility is activation of the surface dimer by moving both H atoms onto one dimer carbon atom (reaction (2b) ). (CHCH --> CH2C) This process is endothermic by 59.5 kcal/mol and requires 62.3 kcal/mol to surmount the barrier. The resulting CH2C structure is 40.4 kcal/mol below the surface dimer plus H2 and there is a small barrier ( about 6 kcal/mol ) to addition of H2 to the localized carbene center to give a dihydride structure. The other possibility involves a competition between the dimerization energy per surface dimer and the small barrier for addition of H2 to a localized carbene center. However the dimerization energy for diamond 100 is computed to be 69.3 kcal/mol per dimer and therefore the barrier to addition would be predicted to be greater than 70 kcal/mol. (This saddle point has not been located.)

For hydrogen abstraction by a gas phase H atom (reaction(3) ), there is a barrier of about 6 kcal/mol with respect to H atom and the surface CHC species. For the process H2 + CC --> H + HCC the reaction is approximately thermoneutral on the triplet surface, but is uphill by about 10 kcal/mol on the singlet surface in addition to the 6 kcal/mol barrier. For the silicon triplet surface, this reaction is uphill by 0.15 aJ (22 kcal/mol), but there is no barrier other than the endothermicity. For the reaction H2 + CCH --> H + HCCH, the reaction is exothermic by 0.4 kcal/mol and has an approximately (6 kcal/mol barrier.

For the addition of a H atom to a surface dimer, there is no barrier for the diamond surfaces.

Another process which was studied for the 100 surface is the barrier to exchange of an H atom ( HCC --> CCH). The barriers here are 63.8 kcal/mol for diamond surfaces.

Important differences are seen between diamond and silicon with respect to the abstraction process. For the diamond surface, the abstraction process has the lowest barrier (other than the H atom addition, which is barrierless on both surfaces) but the C2v constrained addition of H2 has a large barrier. On the silicon surface, this ordering is reversed and the symmetry constrained addition of H2 has a lower barrier than the abstraction reaction. These differences are related to the properties of the surface dimer, especially the lower amount of p bonding for the silicon surface as compared to the diamond surface.

Calculations were also carried out for the analagous reactions on diamond 110 and 111 surfaces. For the diamond 110 surface the exchange barrier is larger 0.50 aJ (72.0 kcal/mol) than the 0.44 aJ (63.8 kcal/mol) barrier obtained for the 100 surface. The barrier to symmetry constrained addition of H2 is 0.30 aJ (43 kcal/mol) as compared to 0.34 aJ (49 kcal/mol) for the 100 surface. Both of these results are related to the lower overlap of the two orbitals of the dimer pair p bond for the 110 surface as compared to the 100 surface. Abstraction barriers were also computed. These are 0.04 aJ (5.2 kcal/mol) and 0.04 aJ (5.8 kcal/mol) for the reactions H + CHC --> H2 + CC and H + CHCH --> H2 + CCH, respectively, and the reactions are approximately thermoneutral. For the 110 surface the non-symmetry constrained addition of H2 would be expected to be unfavorable, since breaking the dimer bond involves breaking at least 3 CC bonds.

For the 111 surface of diamond, the surface dangling bonds are only very weakly interacting, and the only reaction that was considered likely on this surface is abstraction of a H atom by a gas phase H atom. If the reaction is written in the direction H2 + CCH --> H + HCCH the barrier is 4.0 kcal/mol and the reaction is exothemic by 0.04 aJ (6.1 kcal/mol). Note that for the diamond 100 surface the analogous reaction is endothermic by 0.07 aJ (10 kcal/mol) in addition to the 0.04 aJ (6 kcal/mol) barrier (on the singlet surface), while for the 110 surface there is a 0.3-0.4 aJ (5-6 kcal/mol) barrier and the reaction is essentially thermoneutral, and for the 111 surface the barrier is reduced to 0.03 aJ (4 kcal/mol) and the reaction is exothermic. These differences relate to the increase in CH bond strength with 111 > 110 > 100, and these differences in turn relate to the differences in overlap of the dangling bonds of a dimer pair.



This research was supported by NASA Computational Nanotechnology Grant, NASA-JPL-Hughes-Raytheon TCA, and the Chevron Chemical Company (Oronite Technology group). The development of the methodologies used here was supported by grants from the DOE-ASCI-ASAP and NSF (CHE 95-12279 and ASC 92-17368). The facilities of the MSC used in this research are also supported by grants from ARO-DURIP, BP Amoco, ARO-MURI (Kiserow), Beckman Institute, Seiko-Epson, Exxon, Avery-Dennison Corp., 3M, Dow, and Chevron Research Technology Co.



  1. W. A. Goddard III, T. H. Dunning Jr., W. J. Hunt and P. J. Hay, Accts. Chem. Res. 6, 368 (1973)
  2. D. J. Tannor, B. Marten, R. Murphy, R. A. Friesner, D. Sitkoff, A. Nicholls, M. Ringnalda, W. A. Goddard III, and B. Honig, J. Am. Chem. Soc. 116, 11875 (1994)
  3. B. H. Greeley, T. V. Russo, D. T. Mainz, R. A. Friesner, J-M. Langlois, W. A. Goddard III, R. E. Donnelly, and M. N. Ringnalda, J. Chem. Phys. 101, 4028 (1994).
  4. X. J. Chen. J-M. Langlois, and W. A. Goddard, Phys. Rev. B 52, 2348 (1995)
  5. S. Dasgupta, T. Yamasaki, and W. A. Goddard III, J. Chem. Phys. 104, 2898 (1996
  6. N. Karasawa, S. Dasgupta, and W. A. Goddard III, J. Phys. Chem. 95, 2260 (1990
  7. S. Dasgupta, W. B. Hammond, and W. A. Goddard III, J. Am. Chem. Soc., 118, 12291-12301 (1996).
  8. N. Karasawa and W. A. Goddard III, Macromolecules 28, 6765 (1995)
  9. N. Karasawa and W. A. Goddard III, Macromolecules 25, 7268 (1992)
  10. S. Dasgupta, K.A. Brameld, C. F. Fan and W. A. Goddard, Spect. Acta. A 53, 1347-1363 (1997).
  11. S. Dasgupta, K. A. Smith, and W. A. Goddard III, J. Phys. Chem. 97, 10891 (1993)
  12. C. B. Musgrave, S. Dasgupta, and W. A. Goddard III, J. Phys. Chem. 99, 13321 (1995)
  13. C. B. Musgrave, S. J. Harris, and W. A. Goddard III Chem. Phys. Lett. 247, 359 (1995)
  14. M. McAdon and W. A. Goddard, III, J. Phys. Chem. 91,2607 (1987).
  15. S. L. Mayo, B. D. Olafson, and W. A. Goddard III, J. Phys. Chem. 94, 8897 (1990).
  16. A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard III, and W. M. Skiff, J. Am. Chem. Soc. 114, 10024 (1992)
  17. T. Cagin, Y. Qi, H. Li, Y. Kimura, H. Ikeda, W. L. Johnson, and W. A. Goddard, III. Bulk Metallic Glasses, eds. A. Inoue, W. L. Johnson, C. T. Liu, MRS Symp. Ser. 554, 43-48 (1999); Y. Kimura, T. Cagin, Y. Qi and W. A. Goddard, III, Phys. Rev. B, submitted; A. Strachan, T. Cagin, O. Gulseren, R. E. Cohen, W. A. Goddard, III, "qEAM FF for Tanatalum," unpublished.
  18. T. Cagin, E. Demiralp, and W. A. Goddard, III, MRS Symp. Proc. Vol 492, (1998) 287-292.; E. Demiralp, T. Cagin and W. A. Goddard, Phys. Rev. Lett. 82, 1708-1711 (1999);
  19. N. T. Huff, E. Demiralp, T. Cagin and W. A. Goddard, III, J. Noncrystal. Solids 253, 133-142 (1999);
  20. O. Kitao, E. Demiralp, T. Cagin, S. Dasgupta, M. Mikami, K. Tanabe, S. Ono, W. A. Goddard III, Comput. Mat. Sci. 14, 135-138 (1999).
  21. A. Strachan, T. Cagin, W. A. Goddard, III, Phys. Rev. B 60, 15084 (1999);
  22. J.W. Che, T. Cagin, W. A. Goddard, III, Theo. Chem. Acct, 102, 346-354 (1999);
  23. A. Duin and W. A. Goddard, III to be submitted,
  24. N. Karasawa and W. A. Goddard III, J. Phys. Chem. 93, 7320 (1989)
  25. Z. M. Chen, T. Cagin and W. A. Goddard, III, J. Comp. Chem. 18, 1365 (1997).
  26. H. Q. Ding, N. Karasawa and W. A. Goddard III, J. Chem. Phys. 97, 6 (1992)
  27. Lim, K.T., Brunett, S., Iotov, M., McClurg, R.B., Vaidehi, N., Dasgupta, S., Taylor, S. and Goddard III, W.A. (1997) J. Comp. Chem. 18, 501-521.
  28. M. Iotov, S. Kashihara, S. D. Dasgupta, W. A. Goddard, III, "Diffusion of Gases in Polymers", unpublished; M. Iotov, Ph. D. Thesis, California Institute of Technology, December 1997.
  29. P. Miklis, T. Cagin, and W. A. Goddard III, J. Am. Chem. Soc. 119, 7458 (1997).
  30. H. Ding, N. Karasawa, and W. A. Goddard III, Chem. Phys. Lett. 196, 6 (1992)
  31. A. M. Mathiowetz, A. Jain, N. Karasawa, and W. A. Goddard III Proteins 20, 227 (1994)
  32. N. Vaidehi, A. Jain, and W. A. Goddard III, J. Phys. Chem. 100, 10508 (1996).
  33. A. Fijany, T. Cagin, A.J-Botero, and W. A. Goddard, III, Adv. Eng. Software, 29, 441-450 (1998).
  34. A. Fijany, A. J-Botero, T. Cagin, in Parallel Computing: Fundamentals, Applications and New Directions, Eds. D'Hollander, G. Joubert, F. Peters, and U. Trottenberg, 505-515 (1998).
  35. G. Gao, M. Iotov, T. Cagin, and W. A. Goddard III, "MPSim97: Massively Parallel molecular dynamics Simulation program", in preparation
  36. T. Cagin, N. Karasawa, S. Dasgupta, and W. A. Goddard III, Mat. Res. Soc. Symp. Proc. 278-283 (1992)
  37. T. Cagin, W. A. Goddard III, and M. L. Ary, Comp. Polymer Sci. 1 (1991) 241.
  38. M. Parrinello and A. Rahman, Phys. Rev. Lett. 45, 1196 (1980); M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
  39. S. Nose, Mol. Phys. 52, 255 (1984); S. Nose, J. Chem. Phys. 81, (1984) 511.
  40. Yue Qi, T. Cagin. Y. Kimura, W. A. Goddard III, "Shear Viscosity of a Liquid Metal Alloy from NEMD: Au-Cu", unpublished.
  41. J. Che, T. Cagin and W. A. Goddard, III, "Thermal conductivity of diamond from NEMD", in preparation.
  42. T. Cagin, P. Miklis, and W. A. Goddard, III, "Effect of molecular topology on the shear viscosity: linear, star and hyper-branched alkanes," in preparation.
  43. T. Cagin, J. Che, M. N. Gardos, A. Fijany, and W. A. Goddard, III, Nanotech. 10, 278-284 (1999).
  44. T. Cagin, Y. Zhou, E. S. Yamaguchi, R. Frazier, A. Ho, Y. Tang, and W. A. Goddard III, in Dynamics in Small Confining Systems V, MRS Symp. Ser. 543, 79-84 (1999). Eds. J. M. Drake, G. S. Grest, J. Klafter, and R. Kopelman; Y. Zhou, T. Cagin, E.S. Yamaguchi, A. Ho, Y-C. Tang, and W. A. Goddard III, "Dynamical Shear Studies of Lubricant Films Confined to Nanoscale Separation Between Iron Oxide Surfaces Covered with Wear Inhibitors," J. Phys. Chem. submitted.
  45. Y. Qi, T. Cagin and W. A. Goddard, III, "Glass formation and Crystallization in Ni:Cu and Cu:Ag alloys'" submitted (1998)
  46. Engineering Microscopic Machines, by Ken Gabriel, Scientific American, Sept. 1995, Vol. 273, No. 3, 118-121.
  47. Silicon micromechanics: sensors and actuators on a chip, by Roger T. Howe, Richard S. Muller, Kaigham J. Gabriel, William S. N. Trimmer. IEEE Spectrum, July 1990.
  48. 3.Fabrication Technology for an Integrated Surface-Micromachined Sensor, by Theresa A. Core, W. K. Tsang, Steven J. Sherman. Solid State Technology, October 1993.
  49. Micro-machines hold promise for aerospace, by William B. Scott. Aviation Week and Space Technology, March 1993.
  50. Mighty Mites Hit it Big, by Peter Coy, John Carey, and Nei Gross. Business Week, April 1993.
  51. Micron Machinations, by Gary Stix. Scientific American, November 1992.
  52. Mirrors on a Chip, by Jack M. Younse. IEEE Spectrum, November 1993.
  53. M. N. Gardos (1996) Trib. Lett. Volume 2, pages 173-187, Surface chemistry-controlled tribological behavior of Si and diamond.
  54. M.N. Gardos, Tribol. Lett. 2 (1996) 355.
  55. M. N. Gardos (1998) Trib. Lett. Volume 4, pages 175-188, Re(de)construction-induced friction signatures of polished polycrystalline diamond films in vacuum and Hydrogen
  56. M.N. Gardos, in "Protective Coatings and Thin Films," Proc. NATO Adv. Res. Workshop, May 30 - June 5, 1996, NATO ARW Series, eds. Y. Pauleau and P.B. Barna, (Kluwer Academic Publishers, Dordrecht, The Netherlands, 1997), p.185.
  57. M.N. Gardos and B.L.Soriano, J. Mater. Res. 5 (1990) 2599.
  58. M.N. Gardos, in Synthetic Diamond: Emerging CVD Science and Technology, Electrochem. Soc. Monograph, eds. K.E. Spear, and J.P. Dismukes (Wiley, New York, 1994), ch. 12, p. 419.
  59. M.N. Gardos and K.V. Ravi, Dia. Films & Technol. 4, (1994) 139. (1993) 674.
  60. A. Fijany, T. Cagin, A. J. Botero, and W. A. Goddard, III. "Novel Algorithms for massively parallel, long term simulation of molecular dynamics systems," Advances in Engineering Software, 29, 441-450 (1998).
  61. D. W. Brenner, ``Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films,'' Phys. Rev. B 42, 9458 (1990).
  62. J. Che, T. Cagin, and W. A. Goddard, III, Extenion of Bond Order Dependent Potentials to include Long Range Interactions Theochim. Act. (in press).
  63. J. Che, T. Cagin, and W. A. Goddard, III, "Studies of Fullerenes and Carbon Nanotubes by an Extended Bond Order Potential" this conference.
  64. S. P. Walch, W. A. Goddard III, and T. Cagin "Computational Studies of the Interaction of H/H2 with Diamond and Silicon Surfaces", this conference
  65. B. Bhushan, J. N. Isrealachvili, U. Landman, Nature, 1995, 374, 607.
  66. Thompson, P. A.; Troian, S. M. Nature, 1997, 389, 360.
  67. 3 Granick, S. Science, 1991, 253, 1374.
  68. 4 Yoshizawa, H.; Israelachvili, J. N. J. Phys. Chem., 1993, 97, 11300.
  69. 5 Georges, J. M.; Tonick, A., Poletti; S., Yamaguchi, E. S.; Ryason, P. R. Trib. Trans. 1998, 41, 543.
  70. 6 Jiang, S.; Dasgupta, S.; Blanco, M.; Frazier, R.; Yamaguchi, E. S.; Tang, Y., Goddard, W. A., III J. Phys. Chem. 1996, 100, 15760.
  71. 7 Jiang, S.; Frazier, R.; Yamaguchi, E. S.; Blanco, M.; Dasgupta, S.;Zhou, Y.; Cagin, T.; Tang, Y.; and Goddard, W. A., III J. Phys. Chem. B, 1997, 101, 7702.
  72. 8 Zhou, Y.; Jiang, S.; Cagin, T.; Yamaguchi, E. S.; Frazier, R.; Ho, A.; Tang, Y., Goddard, W. A., III, J. Phys. Chem. A, 2000, 104, 2508.
  73. 9 Zhou, Y.; Cagin, T.; Yamaguchi, E. S.; Frazier, R.; Ho, A.; Tang, Y., Goddard, W. A., III ``Dynamical Shear Studies of Nanoscale Lubricant Films Confined Between Iron Oxide Surfaces Covered with DTP Wear Inhibitors". J. Phys. Chem., 2000, submitted.
  74. 10 Mayo, S. L.; Olafson, B. D.; Goddard, W. A., III J. Phys. Chem. 1990}, 94, 8897.
  75. 11 Rhykerd, C. L. Jr.; Schoen, M.; Diestler, D. J.; Cushman, J. H. Nature, 1987, 330, 461. Stevens, M. J.; Mondello, M.; Grest, G. S.; Cui, S. T.; Cochran, H. D.; Cummings, P. T. J. Chem. Phys., 1997, 106, 7303.
  76. Yoshizawa, H.; Chen, Y. L; Israelachvili, J. N. J. Phys. Chem., 1993, 97, 4128.
  77. A. D. Becke, J. Chem. Phys., 98, 5648 (1993).
  78. F. W. Bobrowicz and W. A. Goddard III, in Methods of Electronic Structure Theory, Modern Theoretical Chemistry, Plenum, New York, 1977, pp 79-126
  79. R. Ditchfield, W. J. Hehre, and J. A. Pople, J. Chem. Phys., 54, 724(1971).
  80. M. J. Frisch, G. W. Trucks, H. B. Schlegel, P. M. W. Gill, B. G. Johnson, M. A. Robb, J. R. Cheeseman, T. Keith, G. A. Petersson, J. A. Montgomery, K. Raghavachari, M. A. Al-Laham, V. G. Zakrzewski, J. V. Ortiz, J. B. Foresman, J. Cioslowski, B. B. Stefanov, A. Nanayakkara, M. Challacombe, C. Y. Peng, P. Y. Ayala, W. Chen, M. W. Wong, J. L. Andres, E. S. Replogle, R. Gomperts, R. L. Martin, D. J. Fox, J. S. Binkley, D. J. Defrees, J. Baker, J. P. Stewart, M. Head-Gordon, C. Gonzalez, and J. A. Pople, Gaussian 94, Revision D.1, Gaussian, Inc., Pittsburgh PA, 1995.
  81. T. H. Dunning, Jr., J. Chem. Phys., 90, 1007 (1989).
  82. P. J. Knowles and H.-J. Werner, Chem. Phys. Lett., 145, 514(1988).
  83. M. V. Ramana Murty and H. A. Atwater, Phys. Rev. B, 51, 4889(1995).
  84. SIRIUS is an MCSCF program written by H. J. Jensen and H. Agren and ABACUS is an MCSCF derivatives program written by T. Helgaker, H. J. Jensen, P. Jørenson, J. Olsen, and P. R. Taylor.
  85. W. Stevens, H. Basch, and J. Krauss, J. Chem. Phys., 81, 6026(1984).
  86. S. P. Walch and R. C. Merkle, Nanotechnology, 1998
  87. J. Wei and J. T. Yates, Jr., Critical reviews in Surface Chemistry, 5, 1 (1995).
  88. H. -J. Werner and P. J. Knowles, J. Chem. Phys. , 89, 5803(1988).
  89. C. T. Wu and E. A. Carter, Chem. Phys. Lett, 185, 172 (1991).