next up previous contents
Next: The tight binding model Up: Interatomic interactions Previous: Interatomic interactions

The Tersoff potential

The molecular dynamics technique is a computer calculation that solves numerically a multi-body problem of mechanics. Here, it is applied to the description of velocities and positions of carbon atoms disrupted in the diamond structure, as a function of time. From the initial position and velocity of the atoms, and the forces of interaction between them, the microscopic behavior of the system can be computed by solving differential equations of motion. The molecular dynamics technique applied here is a deterministic process (in the limits of the computer capabilities) since no stochastic factors are involved [91].

The basic assumption made in the present calculation is that the dynamics can be treated classically and that the atoms are spherical and chemically inert. Many mathematical models were proposed to simulate the interatomic potential energy, and from it, the interaction forces. A common feature of these models is the resemblance to the Taylor expansion of the energy as a function of the atomic positions:

\begin{displaymath}E = \sum_{i} V_1(r_i) + \sum_{i<j} V_2(r_i,r_j) + \sum_{i<j<k}
V_3(r_i,r_j,r_k) + \ldots
\end{displaymath} (7.1)

V1, V2, V3 are the one body, two body and three body potentials and depend on the coordinates of each atom. The first term of the expansion is related to an external force and thus disappears if we consider just the interatomic forces. The pair potential (V2) alone may be appropriate for rather closed packed structures (liquefied noble gases Ar, Kr, Xe are typical for such systems) but unsuitable for describing strongly covalent systems with more a open structure. Indeed, no reasonable pair potential will stabilize the diamond structure. For this purpose, the next term in the expansion (the three body potential V3) should be taken into account. Added to the two body term, it can give a sufficiently accurate description of the real physical system of interest [18,92].

The specific form of the V2 term (in fact, only V2 and V3 are significant because any external interactions are not usually taken into account and adding more terms in the expansion will make the computation impracticable) varies from the 1/rn interaction (a ``Lennard-Jones'' type), to the $e^{-\alpha r}$ interaction (a Morse type) [92,93], or a combination of them [18]. In these cases, a cutoff function is added to limit the range of the potential and permit a reduction in computational time.

J. Tersoff [93] abandoned the use of N-body potential form and proposed a new approach by effectively coupling two body and higher multi atom correlations into the model. The central idea is that in real systems, the strength of each bond depends on the local environment, i.e. an atom with many neighbors forms weaker bonds than an atom with few neighbors. Then, J. Tersoff developed a pair potential the strength of which depends on the environment. It was calibrated firstly for silicon [93] and later for carbon [47]. As for the Biswas and Hamann potential [92] , the Morse form is adopted, related to the exponential decay dependence of the electronic density. It is written in the following form :

  
$\displaystyle E = \sum_{i} E_i = \frac{1}{2} \sum_{i\neq j} V_{ij}$     (7.2)
       
Vij = fC(rij)[fR(rij)+bijfA(rij)]     (7.3)

where the potential energy is decomposed into a site energy Ei and a bonding energy Vij, rij is the distance between the atoms i and j, fA and fR are the attractive and repulsive pair potential respectively, and fC is a smooth cutoff function.
$\displaystyle f_R(r) = Ae^{(-{\lambda_1 r})}$     (7.4)
$\displaystyle f_A(r) = -Be^{(-{\lambda_2 r})}$      


\begin{displaymath}f_C(r) = \cases {1, & $r < R-D$ \cr \frac{1}{2} - \frac{1}{2}...
...rac{\pi}{2}(r-R)/D], & $R-D<r<R+D$ \cr 0, & $r>R+D$\space \cr}
\end{displaymath} (7.5)

It has to be noticed that the parameters R and D are not systematically optimized but are chosen so as to include the first-neighbor shell only for several selected high-symmetry bulk structure of silicon, namely for Si2, graphite, diamond, simple cubic, and face-centered cubic structures. The fC function, thus, decreases from 1 to 0 in the range R-D<r<R+D. The main feature of this potential is the presence of the bijterm. As explained before, the basic idea is that the strength of each bond depends upon the local environment and is lowered when the number of neighbors is relatively high. This dependence is expressed by bij, which can accentuate or diminish the attractive force relative to the repulsive force, according to the environment, such that

$\displaystyle b_{ij} = \frac{1}{(1+\beta ^n \zeta_{ij} ^n)^{1/2n}}$     (7.6)
       
$\displaystyle \zeta_{ij} = \sum_{k\neq i,j} f_C(r_{ij})g(\theta_{ijk})
e^{[\lambda_3 ^3(r_{ij} - r_{ik})^3]}$      
$\displaystyle g(\theta) = 1+\frac{c^2}{d^2}-\frac{c^2}{[d^2+(h-\cos\theta)^2]}$      

The term $\zeta_{ij}$ defines the effective coordination number of atom i, i.e. the number of nearest neighbors, taken into account the relative distance of two neighbors rij - rik and the bond-angle $\theta$. The function $g(\theta)$ has a minimum for $h =
\cos (\theta)$, the parameter d determines how sharp the dependence on angle is, and c expresses the strength of the angular effect.

This potential and the parameters were chosen to fit theoretical and experimental data obtained for realistic and hypothetical silicon configurations, namely the cohesive energy of several high-symmetry bulk structures mentioned above, the lattice constant and bulk modulus of the silicon lattice in the diamond configuration. It has later been calibrated for carbon atoms by J. Tersoff [47] in the same manner, constraining the vacancy formation energy in diamond to have at least 4 eV, close to the value found by Bernholc et al[55]. The following parameters (table 7.1) are reported from paper [14] that shows a model of interatomic potentials for multicomponent systems, taking $\lambda_3 = 0$.


 
Table 7.1: Parameters for the Tersoff potential, for atoms of carbon and silicon.
  C Si
A (eV) 1393.6 1830.8
B (eV) 346.7 471.18
$\lambda (\AA^{-1})$ 3.4879 2.4799
$\mu (\AA^{-1})$ 2.2119 1.7322
$\beta$ $1.5724\times 10^{-7}$ $1.1\times 10^{-6}$
n 0.72751 0.78734
c $3.8049\times 10^4$ $1.0039\times 10^5$
d 4.384 16.217
h -0.57058 -0.59825
$R (\AA)$ 1.8 2.7
$S (\AA)$ 2.1 3.0
 

The potential for carbon atoms only [47] was tested by calculating the cohesive energy and the structure of diverse carbon geometries, the elastic constants and phonon frequencies, defect energies and migration barriers in diamond and graphite (the energy required for interstitial and vacancy migration). The results obtained with the parameters given in table 7.1 are in good agreement with experiment (for elastic constants and phonon dispersion) and with ab initio calculations [55] (for defect energies). The first nearest neighbor distance for graphite calculated with this potential is 1.46 Å, rather close to the theoretical value of 1.42 Å[1], keeping in mind that the diamond lattice constant was the only length used for calibration. The cohesive energy calculated is -7.37 eV/atom for diamond and -7.39 eV/atom for graphite, a difference of 0.02 eV/atom, the same as found by others [1]. This ability to describe the small energy difference between diamond and graphite is one of the attractive features of this potential, and shows that it is able to distinguish among different carbon environments, fourfold sp3 bond as well as threefold sp2bond [48].


next up previous contents
Next: The tight binding model Up: Interatomic interactions Previous: Interatomic interactions
David Saada
2000-06-22