Interface pinning
Interface pinning[1] is used to determine the melting point from a molecular-dynamics simulation of the interface between a liquid and a solid phase. The typical behavior of such a simulation is to freeze or melt, while the interface is pinned with a bias potential. This potential applies an energy penalty for deviations from the desired two-phase system. It is preferred simulating above the melting point because the bias potential prevents melting better than freezing.
The Steinhardt-Nelson[2] order parameter [math]\displaystyle{ Q_6 }[/math] discriminates between the solid and the liquid phase. With the bias potential
- [math]\displaystyle{ U_\text{bias}(\mathbf{R}) = \frac\kappa2 \left(Q_6(\mathbf{R}) - A\right)^2 }[/math]
penalizes differences between the order parameter for the current configuration [math]\displaystyle{ Q_6({\mathbf{R}}) }[/math] and the one for the desired interface [math]\displaystyle{ A }[/math]. [math]\displaystyle{ \kappa }[/math] is an adjustable parameter determining the strength of the pinning.
Under the action of the bias potential, the system equilibrates to the desired two-phase configuration. An important observable is the difference between the average order parameter [math]\displaystyle{ \langle Q_6\rangle }[/math] in equilibrium and the desired order parameter [math]\displaystyle{ A }[/math]. This difference relates to the the chemical potentials of the solid [math]\displaystyle{ \mu_\text{solid} }[/math] and the liquid [math]\displaystyle{ \mu_\text{liquid} }[/math] phase
- [math]\displaystyle{ N(\mu_\text{solid} - \mu_\text{liquid}) = \kappa (Q_{6,\text{solid}} - Q_{6,\text{liquid}})(\langle Q_6 \rangle - A) }[/math]
where [math]\displaystyle{ N }[/math] is the number of atoms in the simulation.
Computing the forces requires a differentiable [math]\displaystyle{ Q_6(\mathbf{R}) }[/math]. In the VASP implementation a smooth fading function [math]\displaystyle{ w(r) }[/math] is used to weight each pair of atoms at distance [math]\displaystyle{ r }[/math] for the calculation of the [math]\displaystyle{ Q_6(\mathbf{R},w) }[/math] order parameter. This fading function is given as
- [math]\displaystyle{ w(r) = \left\{ \begin{array}{cl} 1 &\textrm{for} \,\, r\leq n \\ \frac{(f^2 - r^2)^2 (f^2 - 3n^2 + 2r^2)}{(f^2 - n^2)^3} &\textrm{for} \,\, n\lt r\lt f \\ 0 &\textrm{for} \,\,f\leq r \end{array}\right. }[/math]
Here [math]\displaystyle{ n }[/math] and [math]\displaystyle{ f }[/math] are the near- and far-fading distances, respectively.
The radial distribution function [math]\displaystyle{ g(r) }[/math] of the crystal phase yields a good choice for the fading range.
To prevent spurious stress, [math]\displaystyle{ g(r) }[/math] should be small where the derivative of [math]\displaystyle{ w(r) }[/math] is large.
Set the near fading distance [math]\displaystyle{ n }[/math] to the distance where [math]\displaystyle{ g(r) }[/math] goes below 1 after the first peak.
Set the far fading distance [math]\displaystyle{ f }[/math] to the distance where [math]\displaystyle{ g(r) }[/math] goes above 1 again before the second peak.
References