Category:Interface pinning: Difference between revisions

Theory

Interface Pinning is a method for finding melting points from an MD simulation of a system where the liquid and the solid phase are in contact. To prevent melting or freezing at constant pressure and constant temperature, a bias potential applies a penalty energy for deviations from the desired two-phase system.

The Steinhardt-Nelson order parameter ${\displaystyle Q_{6}}$ is used for discriminating the solid from the liquid phase and the bias potential is given by

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

where ${\displaystyle Q_{6}({\mathbf {R} })}$ is the ${\displaystyle Q_{6}}$ order parameter for the current configuration ${\displaystyle \mathbf {R} }$ and ${\displaystyle a}$ is the desired value of the order parameter close to the order parameter of the initial two-phase configuration.

With the bias potential enabled, the system can equilibrate while staying in the two-phase configuration. From the difference of the average order parameter ${\displaystyle \langle Q_{6}\rangle }$ in equilibrium and the desired order parameter ${\displaystyle a}$ one can directly compute the difference of the chemical potential of the solid and the liquid phase:

${\displaystyle N(\mu _{\textrm {solid}}-\mu _{\textrm {liquid}})=\kappa (Q_{6{\textrm {solid}}}-Q_{6{\textrm {liquid}}})(\langle Q_{6}\rangle -a)}$

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

It is preferable to simulate in the super-heated regime, as it is easier for the bias potential to prevent a system from melting than to prevent a system from freezing.

${\displaystyle Q_{6}(\mathbf {R} )}$ needs to be continuous for computing the forces on the atoms originating from the bias potential. 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. A good choice for the fading range can be made from the radial distribution function ${\displaystyle g(r)}$ of the crystal phase. We recommend to use the distance where ${\displaystyle g(r)}$ goes below 1 after the first peak as the near fading distance ${\displaystyle n}$ and the distance where ${\displaystyle g(r)}$ goes above 1 again before the second peak as the far fading distance ${\displaystyle f}$. ${\displaystyle g(r)}$ should be low where the fading function has a high derivative to prevent spurious stress.

How to

The interface pinning method uses the ${\displaystyle Np_{z}T}$ ensemble where the barostat only acts on the direction of the lattice that is perpendicular to the solid liquid interface. This uses a Langevin thermostat and a Parrinello-Rahman barostat with lattice constraints in the remaining two dimensions.

The following variables need to be set for the interface pinning method:

• OFIELD_Q6_NEAR:
• OFIELD_Q6_FAR:
• OFIELD_KAPPA: This tag defines the coupling strength ${\displaystyle \kappa }$ of the bias potential.
• OFIELD_A: This tag 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

Contents