The interatomic potential

The heart of the MD simulations is the inter-atomic potential. In classical simulations the atoms are most often represented by point-like centers, which interact through many-body interactions potential. In that way the highly complex description of electron dynamics is abandoned and effective picture is adopted. In this picture the main features like the hard core of particles and internal degrees of freedom are modeled by a set of parameters and analytical functions, which depend on the mutual positions of the atoms in the configuration. These parameters and functions give a complete information about the system energy, as well as about the forces acting on each particle.

The best choice of a potential for simulations of metals is a many-body potential. It is well known fact that pair-wise potentials, like Lenard - Jones (LJ) potentials :

$$\Phi_{LJ}(r)=4\epsilon(\frac{\sigma^{12}}{r^{12}}-\frac{\sigma^6}{r^6})$$ (4.4)

do not give adequate description of all the properties of metals. For example, the LJ potential imposes the Cauchy relation $C_{12}=C_{44}$ ( $C_{12},~C_{44}$ are elastic constants). This relation is proved to be wrong for most of the metals. Pair-wise potentials fail to estimate the structure relaxation and reconstruction around point defects (vacancies and self-interstitials) in metals. The vacancy formation energy obtained by means of pair - wise potentials is overestimated, and is found to be about equal to the bulk cohesive energy.

The solution of the problem is an introduction of a many-body potential, which include a pair-wise interactions only as part of the full potential. This first part of the many-body potential accounts for the core - core interactions (or ion - ion interactions), while the second part incorporates the complex nature of metallic cohesion by an additional term:

$$V = V_1+V_2 = \frac {1}{2} \sum_{i:~i \neq j}^N \Phi(r_{ij}) +\sum_{i=1}^N U(n _i)$$ (4.5)

where the $\Phi(r)$ is the two-body part, and the many-body part $U(n _ i)$ depends on electronic charge density $n _ i$ around the atom $i$:

$$n_i=\sum_{i \neq j}^{N_g} \rho (r_{ij})$$ (4.6)

$N_g$ is the number of nearest neighbors of the atom $i$.

Various many-body potentials for metals, include the Finnis-Sinclair potentials (FS) [44], the effective medium potentials [45], tight-binding type potentials [47,48], as well as potentials based on the embedded atom methods (EAM) [43] and the glue model [46]. Many-body potentials for fcc metals are well developed and thoroughly tested in numerous simulations. The situation with the many-body potentials for bcc metals is not so good. That can be explained by the more complicated nature of the bcc metals (especially transition bcc metals), where d-orbitals have to be taken into account. Different versions of the many-body potentials for bcc metals are often adopted for some specific problems, for example, to describe defect configurations in transition metals at low temperatures. Therefore, it is not clear if these potentials could be applied to simulate bcc metals at high temperatures up to the melting point. One way to find the answer is to perform simulations, using these potentials, and compare the results with available experimental data and theoretical predictions.

In our project, the Finnis - Sinclair (FS) potential for a bcc metal vanadium [42] was chosen, for the following reason:
First, the FS potential gives good description of the properties of transition bcc metals, and it is widely used in Molecular Dynamics and Monte - Carlo simulations. Second, the FS potential is an analytical one, which is more convenient for to use in MD simulations, than the numerical EAM type potentials. It is especially important when we calculate high-order derivatives of the potential, for example, in calculation of the shear elastic moduli.

The FS potential includes the pair-wise and many-body parts:

$$U =U_1+U_2=\frac {1}{2} \sum_{i,j:~i \neq j}^N V_1(r_{ij}) +\sum_{i=1}^N V_2(n _i)$$ (4.7)

The pair-wise part $V_1(r_{ij})$ is represented by a quartic form (see below), and the many-body part $V_2(n _i)$ of the potential is given by:

$$V_2(n _i)=f(n_i)=-A\sqrt{n_i}$$ (4.8)

where $n _ i$ is the local electronic charge density around the atom $i$:

$$n_i=\sum_{i:~i \neq j}^{Ng} \Phi(r_{ij})$$ (4.9)

the function $\Phi(r_{ij})$ is given by:

$$\Phi(r_{ij})=\left\{ \begin{array}{rcl} (r_{ij}-d)^2,~~r_{ij} \le d\\ 0,~~r_{ij}> d\\ \end{array}\right.$$ (4.10)

where $r_{ji}$ is distance between the atoms $i$ and $j$,
and $A,d$ are the fitting parameters.

The function $f(n_i)$ is obviously taken to mimic the results of the tight-binding theory. Indeed, $n _ i$ can be considered as the local electronic charge density at a site $i$, which is constructed by a rigid superposition of atomic charge densities $\phi(r)$. This is actually a viewpoint of Daw and Baskin [49]. To a first approximation, the energy of an atom at a site $i$ is assumed to be the same as if the atom is placed into an uniform electron gas of that density. The remaining part of the potential is based on empirical fitting. The pair-wise potential represents repulsive core - core interactions. The pair potential can be written in the form:

$$U_1 =\frac {1}{2} \sum_{i:~i \neq j}^N V_1(r_{ij})$$ (4.11)

where the pair potential $V_1$ is a quartic polynomial:

$$V_1(r_{ij})=\left\{ \begin{array}{rcl} (r_{ij}-c)^2\times(c_... ...{ij}^2),~~r_{ij} \le c\\ 0,~~r_{ij}> c\\ \end{array}\right.$$ (4.12)

It is assumed that $c$ lies between the second and the third nearest neighbors interatomic distance. The three fitting parameters $c_0, c_1, c_2$ can be found using experimental data. Values of a lattice spacing $a_0$, cohesive energy $E_g$, bulk modulus $B$, and the three independent elastic constants $C_{11}$,$C_{12}$ and $C_{44}$ have to be taken as an input to find these fitting parameters. In the case of vanadium, which has a bcc structure, the values of the above mentioned physical quantities measured at room temperature, are listed in Table .

Table: Experimental data used for a semi-empirical FS potential of vanadium, from ref. [50,51].

Exp. quantities	Numerical values
$a_0$	$3.0399~ \dot A$
$E_g$	$5.31~ eV$
$C_{11}$	$2.279~ eV/\dot A^3$
$C_{12}$	$1.187~ eV/\dot A^3$
$C_{44}$	$0.426~ eV/\dot A^3$
$B$	$1.551~ eV/\dot A^3$

The elastic coefficients $C_{11}$, $C_{12}$ and $C_{44}$ are taken from the compilation of Bujard [50] and the cohesive energy $E_g$ of vanadium is taken from Kittel [51]. The resulting values for the two-body part of the FS potential are summarized in Table .

Table: Parameters of the $FS$ potential of vanadium, taken from ref. [42].

Fitting parameters	Numerical values
$d$	$3.692767~ \dot A$
$A$	2.010637  eV
$c$	$3.8~ \dot A$
$c_0$	-0.8816318
$c_1$	1.4907756
$c_2$	-0.397637

The most important feature of the FS potential is its many-body part, which is explicitly non-central. Therefore this potential is capable to reproduce correctly the elastic constants. In addition to this property, it reproduces exactly the lattice parameters, bulk modulus, cohesive energy as well as the phonon dispersion, surface and vacancy formation energy in good agreement with experiment [42,44]. However, the chief shortcoming of the potential for us is that it was fitted to the low temperature experimental data, i.e. far away from the melting point.

Despite the fact that the FS potentials were explicitly constructed in order to cope with highly defective systems, such as crystals with point defects, dislocations, cracks, etc., it was found by Rebonato et al.[52], that some modification must to be introduced to the original potential for bcc metals. It was discovered in series of MD and MC simulations that the original FS potentials [42], appear to give unphysical results for properties involving small interatomic separations, namely the pressure versus volume relation [52, 53]. Marshese et al. [53], have shown that the FS potentials predict thermal expansion that is too low in comparison with experiment and in some cases even negative. This flaw implies that the FS potentials may be inaccurate at large distortions of the lattice of the perfect crystal. In order to overcome these problems Rebonato et al.[53], suggested some modifications of the original FS potentials. The repulsive part of theFS potential was changed in the following way:

$$V_1(r) \longrightarrow V_1(r)+h(r)$$ (4.13)

where:

$$h(r)=\left\{ \begin{array}{rcl} K(FND-r)^n,~~r \le FND\\ 0,~~r> FND\\ \end{array}\right.$$ (4.14)

with FND as the first-neighbor distance at zero pressure.

$$FND=a_0\sqrt3/2$$ (4.15)

The many-body part was left unaltered.

These fitting parameters for vanadium are in Table

Table: Parameters of the modified FS potential of vanadium, from ref. [53].

K $eV/ \dot A^3$	n	FND $\dot A$
3.3	3.0	2.63271

This modification improves substantially the original FS potentials to get adequate description of the pressure versus volume behavior, as well as overcome high-compression instability of these potentials. The energies of the stable configurations of vacancies and self - interstitials, are in good agreement with experimental data [53].

Once a potential is selected the energy and the forces can be calculated. Usually the potential is cut-off at some distance to speed up the calculations. The cut-off determines the accuracy of the simulation results, the larger the cut-off length, the more accurate the results. Besides that, the length of simulation box has to be larger than the cut-off length. An atom $i$ interacts only with the atoms, whose distance from the atom $i$ is less than the cut-off length of the potential. The neighbors of the $i$ atom are accounted by means of a neighbor list, which contains information about the environment of the $i$th atom. The neighbor list is usually extended beyond the cutoff length to avoid the need to update it in every simulation time step.