Investigation of the properties of
a perfect crystal of vanadium

The aim of this project is to investigate the rôle of point defects in the bulk melting transition. With this aim in mind we examine the properties of a perfect crystal of vanadium before introducing defects into the bulk of the solid. We have chosen a set of physical properties of the solid, which are relevant in studying of melting transition.

First of all, a structure order parameter $\eta$ was considered which can help in distinguishing between the crystalline and liquid phases, as well as in the detection of the melting transition. The structure order parameter is defined as:

$$\eta=\left < \frac{1}{N}\left \vert \sum_{i=1}^N exp(i\vec k \vec r_i) \right \vert \right >$$ (5.1)

where: $\vec k=[0,0,\frac{2 \pi}{a}]$ is a vector of the reciprocal lattice, $\vec r_i$ is a vector pointing on the atom $i$, $N$ is the number of the atoms in the system, $<>$ is ensemble average. The structure order parameter is unity at zero temperature, when all atoms are at rest at their lattice sites (See ), but at non-zero temperatures $T \le T_m$ its value is less than unity, $0 <\eta < 1$ due to the atom vibrations and defect formation (See Fig. 4.3). This parameter fluctuates around zero ($\eta \simeq 0$) for a molten crystal $T \ge T_m$, when the bcc structure is completely lost (See Fig. ).

Figure: Evolution with time of the structural order parameter at melting transition, initial conditions are those of a solid at $T=2500K$.

Figure: A perfect bcc lattice of vanadium at $T=0K$, $2000$ atoms.

Figure: A bcc lattice of vanadium at T=2100K, 2000 atoms.

The shape and the volume of the sample at various temperatures are other properties of interest to us. Those properties could be investigated by means of the Parinello - Rahman (PR) method [68]. Using this method we inspected the diagonal as well as the off-diagonal elements of the $H_{\alpha \beta}$ matrix which describes the shape and the volume of the sample. The elements of $H_{\alpha \beta}$ were recalculated in each MD step, and their mean values were obtained by averaging over more than $4,000,000$ steps in a steady-state. The variation of the diagonal elements, of a sample containing 2000 atoms is shown in Figure , and the variation of the off-diagonal elements are represented in Figure .

Figure: Variation of the diagonal elements of $H_{\alpha \beta}$ in time at $T= 2450$ K.

Figure: Variation of the off - diagonal elements of $H_{\alpha \beta}$ in time at $T= 2450$ K.

We found that the shape of the computational cell, at zero stress and pressure, is a cubic one. The values of the off-diagonal elements fluctuate around zero, and the difference between the diagonal elements of the $H_{\alpha \beta}$ matrix is negligible. Knowledge of the elements of the matrix $H_{\alpha \beta}$ allows us to calculate the volume of the sample $V=det(H_{\alpha \beta})$ as a function of temperature (See Fig. 4.6). Relying on these data the thermal expansion coefficient is estimated

$$\alpha_t=\frac{1}{V}\frac{\partial V}{\partial T}$$ (5.2)

Figure: Specific volume of a sample with 2000 atoms as function of temperature in comparison with extrapolated experimental data.

We found that the thermal expansion coefficient at low temperatures is
$\alpha_{calc}=(18 \pm 6)\times10^{-6}~[1/K]$ and its value is the same order as the experimental value measured at room temperature $\alpha_{exp}=8.6 \times 10^{-6}~~[1/K]$. The thermal expansion of vanadium was studied experimentally only at the room temperatures, far from the melting point. Thus, in order to compare the results of our simulations at high temperatures with the experiment we would have to extrapolate the available experimental data to the high temperature region. The simplest approximate extrapolation formula for the thermal expansion was applied:

$$v(T)=v_0\exp(\alpha_{\exp}(T-T_0))$$ (5.3)

where $v_0$ is the molar volume at room temperature $T_0$ (See Fig. ). The calculated values are closer to the experimental data at low temperatures. This result could be expected, since the lattice constant $a$ measured at that temperatures was used as an input parameter for the FS potential. At higher temperatures deviation from the experiment is larger, possibly due to the large anharmonicity introduced by the FS potential.

The above mentioned quantities: the structure order parameter, the shape and volume as well as the total energy and pressure of the system were calculated in the framework of the Parinello Rahman method ($NtT$ ensemble). The bulk melting temperature (strictly speaking its upper limit) was also estimated by slowly heating the perfect sample up to the temperature at which the crystal lattice becomes unstable and collapses. At the same time the order parameter jumps to zero and the shape of the sample changes from cubic to tetragonal. In addition, we calculated the shear elastic moduli $C_{\alpha \beta}$ using method proposed by Parinello and Ray ($NVT$ ensemble) (See Appendix C). The stress tensor fluctuations as well as the mean value of the Born term were calculated to find the shear elastic coefficients. The Born term converges very fast in comparison with the stress tensor fluctuations which reach their equilibrium value at a given temperature after a very long time (about $1,500,000$ steps in average). The variation of the elastic coefficients as a function of MD time steps is shown in Fig.

Figure: The shear elastic moduli vs. time at temperature $2200$ K.