Bulk melting and point defects

After investigation of a perfect crystal of vanadium, point defects, i.e. the simplest structural imperfection in solids, are introduced either by removal of atoms (vacancies) from lattice sites or by insertion additional atoms of the same kind (self - interstitials) between the lattice sites (See Fig. ). Initially these point defects are distributed homogeneously inside the bulk of the solid. In our simulations we introduced point defects of one type only to avoid their mutual annihilation or recombination.

Figure: An initial configuration of 5 interstitials in a sample with $2000$ atoms.

The off - lattice Monte Carlo method, namely simulated tempering, was implemented to find the most stable configuration of atoms in the vicinity of a point defect inside the bulk at low temperatures. Our simulations were carried out for a sample containing $128$ atoms plus a self-interstitials ($128+1$), and the temperature set was chosen to be $\{ 75~K,80~K,85~K,90~K\}$. The most energetically favored configuration was found to be the $<011>$ dumb-bell split - interstitial with a formation energy of $E_f=4.18 \pm 0.02~eV$ (See Fig. ).

Figure: A $<011>$ dumb - bell split - interstitial in a bcc metal vanadium.

The defect formation energy was calculated in following way:

$$E_f=E(N-1,1)-E(N,0)$$ (5.4)

where $E(N,N_{def})$ is the potential energy of a sample which contains $N$ lattice atoms and $N_{def}$ point defects. The potential energy of a perfect crystal with $N$ atoms is given by:

$$E(N,0)=NE_{coh}$$ (5.5)

where $E_{coh}$ is the cohesion energy per atom calculated for a pure sample. The calculated value of the defect formation energy $E_f=4.18 eV$ is close enough to the results obtained by Ackland et. al [73] using the DEVIL program which based on the conjugate-gradient method. (See Table 4.1).

Table: Defect formation energy of various split-interstitial defects, from ref. [73].

Type of split - interstitial	Formation energy, eV
$<001>$	4.963
$<011>$	4.163
$<111>$	4.608
crowdion	4.6

Other possible configurations (octahedral, tetrahedral, crowdion) possess larger defect formation energy, and therefore they are less energetically favored and less stable. In our simulations we implemented various initial configurations:either we started very close to the most stable configuration, i.e. $<011>$ dumb-bell split - interstitial, or inserted an additional atom in a random fashion between the lattice sites (See Fig. ) In each case, the configuration with lowest formation energy was found (See Fig. ).

Figure: An initial configuration: an interstitial (white color) and its neighbors.

Figure: Equilibrium configuration of the $< 110 >$ split interstitial.

When the most stable configuration of point defects inside the bulk of vanadium was found, we began to study how point defects influence the various properties of vanadium which is interested to us. To simulate bulk properties of vanadium we prepared various samples with different concentrations of point defects. The MD simulations were performed in the NtT ensemble by using the PR method. We found that introduction of point defects leads to the structural disordering (see Fig. ).

Figure: Structure order parameter as a function of concentration of self-interstitials at several temperatures.

Increase in the concentration of self - interstitials results in noticeable decrease of the structure order parameter $\eta$, while the same effect of vacancies is relatively weak (See Fig. ).

Figure: Structure order parameter as a function of point defect concentration. A comparison between interstials and vacancies. The lines to guide the eye.

Self - interstitials expand the volume of the sample as it is shown in Fig. 4.14 where the lattice parameter $a \simeq (V/N)^{1/3}$, while vacancies lead to decrease of the volume (See Fig. 4.15).

Figure: The lattice parameter $a$ as a function of the concentration self - interstitials.

Figure: The lattice parameter $a$ as a function of the point defect concentration. A comparison between self - interstitials and vacancies.

It is interesting to compare the dependence of the specific volume on the concentration of self - interstitials for vanadium (bcc lattice) and copper (fcc lattice). In both cases the specific volume is increased, but for fcc lattice the effect is more noticeable. This effect can be attributed to the more compact structure of the fcc lattice, where even a small concentration of self - interstitials lead to a large distortion of the fcc lattice, and therefore is more noticeable in comparison with the bcc lattice.

The next stage of our bulk simulations is the most important one - investigation of the rôle of point defects in bulk melting transition. According to the Born criterion [9, 16] bulk melting transition takes place when the specific volume of the crystal reaches a critical value. As shown by Kanigel et. al [17,67] and it does not matter in which way the critical value is reached. The critical volume at which crystal melts could be attained either by heating of the sample or by doping it with point defects at a constant temperature which leads to the expansion of the sample and in the end to melting. The temperature at which melting occurs can lower than the bulk melting point of a perfect sample, i.e. point defects lower the bulk melting temperature! In this sense the mechanical melting process is universal, e.g. it determined only by the sample expansion up to the critical volume.

In our simulations we verified that theoretical prediction. We prepared samples with a specified concentrations of point defects, and heated them gradually up to the melting point. In this way the value of the critical volume and the bulk melting temperature $T_b$ was obtained. By repeating this procedure for various defect concentrations we found the dependence of the bulk melting temperature on the concentration of point defects.

The initial temperature is far below the melting point of a perfect sample $T \simeq 0.7T_{b}$, (the bulk melting temperature $T_b=2500\pm5~K$ is calculated in the simulations of perfect sample). After each $100,000$ MD steps we increased the temperature by $100$ K until we reached $2100$ K. After that the temperature was increased by in a smaller step of $50$ K followed by $200,000$ MD steps. In the end we reached the temperature of $2400$ K, and from this point and onward we increased the temperature incrementally by $10$ K. In this region each sample configuration (positions and velocities of all atoms) was saved before the elevation of the temperature, in order to use the stored configurations again if needed. The number of MD steps between two successive temperature changes was increased to $2,000,000$ MD steps. At some temperature we observed an abrupt decrease of the structure order parameter, and a drastic increase of the total energy and the volume of the sample (See Fig. and Fig. ).

Figure: Increase of the total energy at the melting point.

Figure: Jump of the sample volume at the melting transition.

Figure: Variation of the diagonal elements of the $H_{\alpha \beta}$ matrix at the bulk melting transition.

At that temperature one sees a sharp bifurcation in the lattice dimension where the system elongates in two directions and contracts in the third (See Fig. ). This is a clear sign of symmetry change, from cubic to tetragonal. The same effects were observed at melting transition of fcc metals [63]. Bulk melting transition occurs during a very short time scale corresponding approximately to the several vibration periods of atoms.

It is not improbable that we do not encounter the true bulk melting temperature, but find only its upper limit. It is not known in advance how long the simulation has to be carried out before the expected phenomenon will be observed. There is a possibility that we missed the melting point during the heating of the sample and the melting transition would occur at a lower temperature, provided we could run our simulations for a longer time. Therefore, after the upper limit is detected, we returned to one of the previous configurations. The temperature of the recovered configuration is lower than the bulk melting temperature, but close enough (actually we took the closest one). The simulation were repeated for a quite long time up $15,000,000$ MD steps, in the hope to observe a possible melting transition. If the transition was observed, we repeated the procedure again.

The results of the various simulations performed at the different temperatures and the defect concentrations can be summarized in a phase diagram (See Fig. ).

Figure: The influence of point defects on the melting temperature of vanadium obtaining using periodic boundary conditions.

We see that increase in concentration of the self - interstitials leads to decrease of the bulk melting temperature, while the vacancies almost do not affect the bulk melting temperature, at least if their concentration is small. The same effect of decrease of the bulk melting temperature induced by point defects was obtained by A. Kanigel [67] for the fcc metal copper (See Fig. 4.20).

Figure: The influence of point defects on the melting temperature of Cu, a sample of 1372 atoms, from ref. [67].

If we compare the phase diagram of vanadium with the phase diagram of copper, then it can be immediately seen that the melting point of vanadium is lowered by $50$K at a concentration of self-interstitials about $0.006 \%$, while the melting point of copper is lowered by $80$K at the same defect concentration. Interstitials in copper induce larger distortion of its close packed fcc lattice, than interstititials in vanadium. Therefore, the specific volume of copper is increases more than that of vanadium, and a smaller temperature increase is needed to expand copper up to the critical volume, at which melting occurs.

In conclusion, some additional remarks are necessary:
First, the concentration of point defects in the simulations is increased up to a high enough value $n_{def} \simeq 10^{-2} ~~[defects]/[atoms]$ to see their possible effect on the bcc lattice of vanadium. This is unrealistically large value in comparison with the typical laboratory values $n_{def} \simeq 10^{-6}~~[defects]/[atoms]$. At these high defect concentrations the effect of self - interstitials and vacancies could not be considered as a simple lattice expansion in analogy with the thermal expansion. One has take interactions between these defects into account, which alter in some way the physics of the phenomenon Second, we can not neglect the fact that MD simulation with the Noose - Hoover thermostat is plagued by temperature fluctuations due to the small sample size. That means it is hard to approach $T_b$ to an accuracy of better than about $\sim 1\%$. In summary, it could be noted that the calculated phase diagram is qualitative, because of the finite sample size and the finite time of our computer simulations.