Next: Local density profile
Up: Results: surface melting
Previous: Equilibration and calculation
In order to study the phenomenon of premelting,
we have to find the thermodynamic melting point
of our model. At least three different methods for estimation of the melting point
have been given in the literature. The most elegant and thermodynamically valid
is based on the calculation of the Gibbs free energies of the bulk solid and liquid
phases as a function of temperature at given pressure. The melting point is then found,
by definition,
by making the Gibbs free energies of solid and liquid equal [19].
Hoverer, this method is computationally quite elaborate and contains a
fairly large source of inaccuracy, due to the fact that the Gibbs free energies
usually depend quite weakly on temperature, and is then an intersection
point of two flat curves. Therefore, it is quite difficult to find the position of the intersection point
precisely and a small error in the obtained position of the point leads to large
error in estimation of .
In the second method, proposed by J.F. Lutchko et al. [77] one
determines the temperature dependent velocity of the moving solidtoliquid interface
at and thus, by extrapolating this velocity to the zero
one can obtain the solidliquid coexistence point, which is interpreted as .
This method has been shown to give a melting point comparable to what is obtained
from the free energy method.
The third method, is computationally the most straightforward one. In this method,
an excess amount of kinetic energy is supplied to the atoms of the surface region
of the sample. This results in a temperature gradient and in a nucleation of
liquid phase at the surface. Heating is then stopped and the sample is allowed to evolve back
to equilibrium. The temperature at which a solidliquid coexistence takes place
is interpreted as the thermodynamic melting point [89]. All these methods
can estimate with comparable accuracy.
We chose the second method and determined the velocity of the solidmelt
interface at different temperatures (in the range of 2300  2500 K) in the samples with
various lowindex free surfaces. The value of varies with the geometry
of the simulation cell, and especially with the number of the surface layers.
To test such effects we performed simulations
for the samples with the different number of layers, , using 64 atoms in a layer
(Va(111), Va(001)) and 72 atoms (Va(011)).
A temperature above
(the upper limit of was found in the bulk simulations)
was reached using the following method. We started the simulation at a low temperature ( usually at )
far below , and increased the temperature at a high rate (by
) after each 1000 steps.
In this way the system could be superheated up to . A very simple method of
the velocity rescaling in each MD step is used, e.g. the instantaneous kinetic energy is calculated as
and velocity of each atom is rescaled by the same factor to satisfy
. This simple
method of the velocity rescaling was chosen instead of the NoseHoover mechanism to significantly reduce
the fluctuations of temperature. When the desired temperature is reached,
the velocity of the solidmelt interface is
determined by monitoring the decrease of the structure order parameters,
, as a function of time (See Fig. 5.8).
Figure:
Order parameters of the superheated system vs. time (MD steps) at temperature .
The sample contains atoms.

The order parameter was calculated only over the top 10 surface layers
to avoid the ordering effects introduced by the static substrate.
The temperature dependence of the velocity is governed by the thermodynamic freeenergy
driving force, given by the free energy difference,
, between the solid and the liquid phases [90,91].

(6.1) 
where is the mean free path in the liquid phase, is the thickness of the liquid film,
is the selfdiffusion coefficient of atoms in the liquid phase. One can write with a good approximation [90]
, where is entropy difference between the solid and liquid phases,
and obtain at temperatures close to the melting point

(6.2) 
Using this expression for as a function of temperature, , we
calculate the thermodynamical melting point by a linear fit (See Fig 5.9).
The temperature at which the propagation velocity vanishes
is considered to be the thermodynamical melting point [89].
Figure:
Velocity of propagation of the solidtoliquid interface, , as a function of temperature and
extrapolation to zero velocity. The sample, atoms, layers.

The melting temperatures obtained
for the samples with different number of surface layers and the various types of lowindex faces are
summarized in Table 5.2. The smaller the number of the layers, the larger
the apparent . This is due to the influence of the underlying static substrate, which stabilizes and
orders the samples. Therefore, to obtain more accurate and reliable
estimation for the melting point we use the sample with the largest
number of free surface layers.
Table:
The melting temperatures for the various surfaces of vanadium,
calculated for the samples with different number of the surface layers, .

The calculated melting temperature of the sample with the least packed face Va(111)
is
is the close one to the experimental value .
For other surfaces is insignificantly larger due to superheating effects on the
solidliquid interface just above the thermodynamical melting point .
This effects are very significant for the close packed Va(011) surface.
This dependence of the thermodynamic melting point on the geometry of the crystal face
is also observed in MD simulations when temperature is increased gradually
up to the melting point. It was found the Va(111)
and Va(001) surfaces melt approximately at the same temperature
, while
the sample Va(111), with the close packed free surface, melts at more elevated temperature
about
.
We want to point out again that in the bulk simulations, where we employed periodic boundary conditions
in all directions, melting is observed to happen at a temperature
.
This melting temperature is not the true thermodynamic melting point ,
but a temperature at which mechanical instability arises, leading ultimately
to a mechanical collapse of the crystal lattice. The temperature at which mechanical
instability occurs can be determined by using the Born criterion at zeroexternal stress.
This shows that the superheating region, , in MD simulations of metals may
be large enough (250 K in case of vanadium) and great care should be
exercised while determining from the simulations results of bulk and surface melting.
Next: Local density profile
Up: Results: surface melting
Previous: Equilibration and calculation
20030115