Category:Interface pinning: Difference between revisions

Use interface pinning to determine the melting point from a molecular-dynamics simulation of the interface of a liquid and a solid phase. Because the typical behavior of such a simulation is to freeze or melt, the interface is pinned with a bias potential. This potential applies an energy penalty for deviations from the desired two-phase system. Prefer simulating above the melting point because the bias potential prevents melting better than freezing.

The Steinhardt-Nelson order parameter ${\displaystyle Q_{6}}$ discriminates between the solid and the liquid phase. With the bias potential

${\displaystyle U_{\text{bias}}(\mathbf {R} )={\frac {\kappa }{2}}\left(Q_{6}(\mathbf {R} )-A\right)^{2}}$

penalizes differences between the order parameter for the current configuration ${\displaystyle Q_{6}({\mathbf {R} })}$ and the one for the desired interface ${\displaystyle A}$. ${\displaystyle \kappa }$ 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 ${\displaystyle \langle Q_{6}\rangle }$ in equilibrium and the desired order parameter ${\displaystyle A}$. This difference relates to the the chemical potentials of the solid ${\displaystyle \mu _{\text{solid}}}$ and the liquid ${\displaystyle \mu _{\text{liquid}}}$ phase

${\displaystyle N(\mu _{\text{solid}}-\mu _{\text{liquid}})=\kappa (Q_{6,{\text{solid}}}-Q_{6,{\text{liquid}}})(\langle Q_{6}\rangle -A)}$

where ${\displaystyle N}$ is the number of atoms in the simulation.

Computing the forces requires a differentiable ${\displaystyle Q_{6}(\mathbf {R} )}$. We use a smooth fading function ${\displaystyle w(r)}$ to weight each pair of atoms at distance ${\displaystyle r}$ for the calculation of the ${\displaystyle Q_{6}}$ order parameter

${\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<r<f\\0&{\textrm {for}}\,\,f\leq r\end{array}}\right.}$

Here ${\displaystyle n}$ and ${\displaystyle f}$ are the near- and far-fading distances given in the INCAR file respectively. The radial distribution function ${\displaystyle g(r)}$ of the crystal phase yields a good choice for the fading range. To prevent spurious stress, ${\displaystyle g(r)}$ should be small where the derivative of ${\displaystyle w(r)}$ is large. Set the near fading distance ${\displaystyle n}$ to the distance where ${\displaystyle g(r)}$ goes below 1 after the first peak. Set the far fading distance ${\displaystyle f}$ to the distance where ${\displaystyle g(r)}$ goes above 1 again before the second peak.

How to

Interface pinning uses the ${\displaystyle Np_{z}T}$ ensemble where the barostat only acts along the ${\displaystyle z}$ direction. This uses a Langevin thermostat and a Parrinello-Rahman barostat with lattice constraints in the remaining two dimensions. The solid-liquid interface must be in the ${\displaystyle x-y}$ plane perpendicular to the action of the barostat.

Set the following tags for the interface pinning method:

OFIELD_Q6_NEAR

Defines the near-fading distance ${\displaystyle n}$.

OFIELD_Q6_FAR

Defines the far-fading distance ${\displaystyle f}$.

OFIELD_KAPPA

Defines the coupling strength ${\displaystyle \kappa }$ of the bias potential.

OFIELD_A

Defines the desired value of the order parameter ${\displaystyle A}$.

The following is a sample INCAR file for interface pinning of sodium^[1]:

TEBEG = 400                   # temperature in K
POTIM = 4                     # timestep in fs
IBRION = 0                    # do MD
ISIF = 3                      # use Parrinello-Rahman barostat for the lattice
MDALGO = 3                    # use Langevin thermostat
LANGEVIN_GAMMA = 1.0          # friction coef. for atomic DoFs for each species
LANGEVIN_GAMMA_L = 3.0        # friction coef. for the lattice DoFs
PMASS = 100                   # mass for lattice DoFs
LATTICE_CONSTRAINTS = F F T   # fix x&y, release z lattice dynamics
OFIELD_Q6_NEAR = 3.22         # fading distances for computing a continuous Q6
OFIELD_Q6_FAR = 4.384         # in Angstrom
OFIELD_KAPPA = 500            # strength of bias potential in eV/(unit of Q)^2
OFIELD_A = 0.15               # desired value of the Q6 order parameter

References

Contents