Category:Interface pinning: Difference between revisions

Revision as of 08:40, 7 April 2022

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 [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}) - Q_{6,\text{pinned}}\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{ Q_{6,\text{pinned}} }[/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{ Q_{6,\text{pinned}} }[/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 - Q_{6,\text{pinned}}) }[/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]. We use a smooth fading function [math]\displaystyle{ w(r) }[/math] to weight each pair of atoms at distance [math]\displaystyle{ r }[/math] for the calculation of the [math]\displaystyle{ Q_6 }[/math] order parameter

[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 given in the INCAR file 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.

How to

The interface pinning method uses the [math]\displaystyle{ Np_zT }[/math] ensemble where the barostat only acts in 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: This tag defines the near-fading distance [math]\displaystyle{ n }[/math].
OFIELD_Q6_FAR: This tag defines the far-fading distance [math]\displaystyle{ f }[/math].
OFIELD_KAPPA: This tag defines the coupling strength [math]\displaystyle{ \kappa }[/math] of the bias potential.
OFIELD_A: This tag defines the desired value of the order parameter [math]\displaystyle{ a }[/math].

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

↑ U. R. Pedersen, F. Hummel, G. Kresse, G. Kahl, and C. Dellago, Phys. Rev. B 88, 094101 (2013).

Contents

This category currently contains no pages or media.

[pedersen:prb:13-1] U. R. Pedersen, F. Hummel, G. Kresse, G. Kahl, and C. Dellago, Phys. Rev. B 88, 094101 (2013).

[1]

@@ Line 1: / Line 1: @@
-'''Interface pinning''' is a method for finding melting points from a [[:Category:Molecular dynamics|molecular-dynamics]] simulation of a system where the liquid and the solid phase are in contact.
+Use '''interface pinning''' to determine the melting point from a [[:Category: Molecular dynamics|molecular-dynamics]] simulation of the interface of a liquid and a solid phase.
+<!-- == Theory == -->
+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.
-== Theory ==
+The Steinhardt-Nelson order parameter <math>Q_6</math> discriminates between the solid and the liquid phase.
-To prevent melting or freezing at constant pressure and constant temperature, a bias potential applies penalty energy for deviations from the desired two-phase system.
+With the bias potential
-The Steinhardt-Nelson order parameter <math>Q_6</math> is used for discriminating the solid from the liquid phase and the bias potential is given by
+:<math>U_\text{bias}(\mathbf{R}) = \frac\kappa2 \left(Q_6(\mathbf{R}) - Q_{6,\text{pinned}}\right)^2 </math>
-:<math>U_\textrm{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>Q_6({\mathbf{R}})</math> and the one for the desired interface <math>Q_{6,\text{pinned}}</math>.
+<math>\kappa</math> is an adjustable parameter determining the strength of the pinning.
-where <math>Q_6({\mathbf{R}})</math> is the <math>Q_6</math> order parameter for the current configuration <math>\mathbf{R}</math> and <math>a</math> is the desired value of the order parameter close to the order parameter of the initial two-phase configuration.
+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>\langle Q_6 \rangle</math> in equilibrium and the desired order parameter <math>Q_{6,\text{pinned}}</math>.
+This difference relates to the the chemical potentials of the solid <math>\mu_\text{solid}</math> and the liquid <math>\mu_\text{liquid}</math> phase
-With the bias potential enabled, the system can equilibrate while staying in the two-phase configuration. From the difference between the average order parameter <math>\langle Q_6 \rangle</math> in equilibrium and the desired order
+:<math>
-parameter <math>a</math> one can directly compute the difference of the chemical potential of the solid and the liquid phase:
+N(\mu_\text{solid} - \mu_\text{liquid}) =
+\kappa (Q_{6,\text{solid}} - Q_{6,\text{liquid}})(\langle Q_6 \rangle - Q_{6,\text{pinned}})
-:<math> N(\mu_\textrm{solid} - \mu_\textrm{liquid}) =\kappa (Q_{6 \textrm{solid}} - Q_{6 \textrm{liquid}}) (\langle Q_6 \rangle - a) </math>
+</math>
 where <math>N</math> 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.
+Computing the forces requires a differentiable <math>Q_6(\mathbf{R})</math>.
+<!-- PLEASE REPHRASE - I did not understand this part and how it relates to Q_6(R) -->
-<math>Q_6(\mathbf{R})</math> needs to be continuous for computing the forces on the atoms originating from the bias potential. We use a smooth fading function <math>w(r)</math> to weight each pair of atoms at distance <math>r</math> for the calculation of the <math>Q_6</math> order parameter
+We use a smooth fading function <math>w(r)</math> to weight each pair of atoms at distance <math>r</math> for the calculation of the <math>Q_6</math> order parameter
 :<math> w(r) = \left\{ \begin{array}{cl} 1  &\textrm{for} \,\, r\leq n \\
@@ Line 25: / Line 32: @@
   &\textrm{for} \,\,f\leq r \end{array}\right. </math>
-Here <math>n</math> and <math>f</math> are the near- and far-fading distances given in the {{TAG|INCAR}} file respectively. A good choice for the fading range can be made from the radial distribution function <math>g(r)</math> of the crystal phase. We recommend to use the distance where <math>g(r)</math> goes below 1 after the first peak as the near fading distance <math>n</math> and the distance where <math>g(r)</math> goes above 1 again before the second peak as the far fading distance <math>f</math>. <math>g(r)</math> should be low where the fading function has a high derivative to prevent spurious stress.
+<!-- is w(r) equivalent to (1 - t)^2(1 + 2t) with t = (r - n) / (f - n)? -->
+Here <math>n</math> and <math>f</math> are the near- and far-fading distances given in the {{TAG|INCAR}} file respectively.
+<!-- END REPHRASE -->
+The radial distribution function <math>g(r)</math> of the crystal phase yields a good choice for the fading range.
+To prevent spurious stress, <math>g(r)</math> should be small where the derivative of <math>w(r)</math> is large.
+Set the near fading distance <math>n</math> to the distance where <math>g(r)</math> goes below 1 after the first peak.
+Set the far fading distance <math>f</math> to the distance where <math>g(r)</math> goes above 1 again before the second peak.
 == How to ==
-The '''interface pinning''' method uses the <math>Np_zT</math> 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 '''interface pinning''' method uses the <math>Np_zT</math> ensemble where the barostat only acts in 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: