An essential parameter for crystal growth is the kinetic coefficient given by the proportionality between supercooling and average growth velocity. Here, we show that this coefficient can be computed in a single equilibrium simulation using the interface pinning method where two-phase configurations are stabilized by adding a spring-like bias field coupling to an order-parameter that discriminates between the two phases. Crystal growth is a Smoluchowski process and the crystal growth rate can, therefore, be computed from the terminal exponential relaxation of the order parameter. The approach is investigated in detail for the Lennard-Jones model. We find that the kinetic coefficient scales as the inverse square-root of temperature along the high temperature part of the melting line. The practical usability of the method is demonstrated by computing the kinetic coefficient of the elements Na and Si from first principles. A generalized version of the method may be used for computing the rates of crystal nucleation or other rare events.
I. INTRODUCTION
Crystal growth is of paramount importance in many branches of condensed matter physics.1,2 An important parameter in phase field equations3 describing crystal growth is the kinetic coefficient defined as the proportionality constant between supercooling (or superheating) and the average interface growth (or melting) velocity .4,5 We define the kinetic coefficient M using the difference in chemical potential between the two phases μsl = μs − μl as a measure of supercooling (or superheating),
In the spirit of the fluctuation-dissipation theorem,6 we suggest to learn about the interface dynamics by investigating spontaneous fluctuations when a bias potential is added to the Hamiltonian. Specifically, we determine M indirectly in a computation where trajectories are pinned to two-phase configurations (Fig. 1) by adding a harmonic energy coupling to an order-parameter of crystallinity.7,8 In effect, non-equilibrium crystal growth is converted into a well-defined equilibrium problem. We devise a stochastic model of fluctuations that assume Smoluchowski dynamics for crystal growth and show that the growth rate can be inferred from the terminal exponential relaxation of the order-parameter. The chemical potential difference between the solid and the liquid μsl is known from the average force exerted by the bias field on the system, and thus the proportionality constant M can be computed in a single equilibrium simulation.
Sketch of an elongated periodic simulation box containing a two-phase configuration. The projected interface position, ignoring capillary waves, is xs = Ns/2XYρs.
Sketch of an elongated periodic simulation box containing a two-phase configuration. The projected interface position, ignoring capillary waves, is xs = Ns/2XYρs.
The paper is organized as follows: Sec. II is a brief introduction to the interface pinning method. In Sec. III, we motivate and give the solution to a stochastic model of thermal fluctuations with and without a pinning potential. In Secs. IV and V, we compute crystal growth rates for the Lennard-Jones (LJ) model and combine the method with ab initio density functional theory (DFT) to compute growth rates of the elements Na and Si. The paper is finalized with a discussion of the method (Sec. VI). The Appendix gives a detailed analysis of our stochastic model.
II. INTERFACE PINNING
In the following, we give a brief review of the recently proposed interface pinning method7,8 for studying the solid-liquid phase transition. Consider a system of N particles in an elongated periodic simulation cell as sketched in Fig. 1. For the configuration R = (r1, r2, …, rN), let the system have the Hamiltonian . Hydrostatic pressure is ensured by allowing the box length in the z-direction to change and constant temperature is ensured by connecting the momenta of the particles to a heat bath. We refer to this as the NpzT-ensemble. When the initial R is a two-phase configuration, the thermodynamically stable phase will grow on the expense of the other phase. In an interface pinning computation, this is avoided by adding an auxiliary spring-like energy term with a spring constant κ and an anchor point to the Hamiltonian so that R is biased towards two-phase configurations,
Here, Q(R) is a measure of crystallinity of the system. In practice, we use long-range order as measured by the magnitude of the collective density field |ρ(k)| (Ref. 8) (unless otherwise stated)
where k is a wavevector that corresponds to a Bragg peak. The averaged force α the system exerts on Q is proportional to the chemical potential difference between the two phases
Here Qsl = Qs − Ql and Qs and Ql are the mean values of the order parameter when the system is entirely crystalline or liquid, respectively. When the system is in equilibrium with respect to , the relative position of the interface stops moving up to thermal fluctuations and the force α is balanced by the average force of the spring-like bias-field. Thus, μsl can be computed from the average spring force as follows:
Since, we have used a harmonic pinning potential, the distribution of the order parameter will be Gaussian with variance kBT/κ. The time evolution of the order parameter Q(t) depends on the trajectory R(t) that itself is determined by . The main idea of the method suggested in this paper is that if we can understand the time fluctuations of Q(t) when the interface pinning field is applied, we can deduce how the system evolves in the absence of the interface pinning potential. In Sec. III, we devise a stochastic model that assumes that crystal growth is an over-damped Smoluchowski process. From this model, we deduce the value of the kinetic growth coefficient M.
III. STOCHASTIC MODEL
This section suggests a simple model for the dynamics of Q(t) that enables us to determine the value of M from an interface pinning simulation. Assuming Newtonian dynamics of R, the Q(t)-trajectory is deterministic and uniquely defined from the initial values of R and . We are, however, only interested in the statistical behavior of Q(t), e.g., the autocorrelation function 〈ΔQ(0)ΔQ(t)〉. To this aim, we devise a stochastic model with an effective Hamiltonian that is a function of coarse-grained coordinates. The equations of motion of these coordinates include Langevin noise forces representing degrees of freedom that have been neglected.9 The essential model parameter for determining M is “the friction constant of crystal growth” γ. We determine the value of this parameter by fitting the model to simulation results.
A. Effective Hamiltonian
As mentioned earlier, the interface pinning method requires the definition of an order parameter Q(R) that quantifies how much of the system is crystalline at the configuration R, which gives the positions of all atoms. Ideally, the order parameter Q should be proportional to the number Nc of crystalline particles such that any increase in Q corresponds to a growth of the crystalline region. In practice, however, the crystallinity of the system is measured using the Fourier transform of the density field as order parameter, Eq. (3). In this case, Q changes not only due a decrease or increase of the number of crystalline particles but also due to other spatial fluctuations such as lattice vibrations. Thus, we write the order parameter Q as a sum of two contributions,
where q is assumed to be proportional to the number of crystalline atoms and f takes into account fluctuations that do not change the number of crystalline atoms. The argument t emphasizes that all three variables Q, q, and f evolve in time along the trajectory R(t) of the system.
In an interface pinning simulation, the order parameter Q is computed explicitly, while the variable q, which contains the relevant information about the growth of the crystalline region and the motion of the interface, is not directly accessible. To extract information about the interface dynamics from the simulation, we, therefore, construct a simplified two-dimensional model that separates the effect of q and f on the dynamics of Q. In this model, q and f are the sole dynamical variables. We postulate that their dynamics is governed by the effective Hamiltonian
The first term on the right hand side, αq, arises from the imbalance in chemical potential between solid and liquid, which drives the growth (or decrease) of the crystalline region. This term drives the interface motion in the absence of the pinning potential (κ = 0). The effective force α is proportional to the difference in chemical potential, μsl. The proportionality constant depends on the choice of crystallinity order parameter, Eq. (4). The second term, , couples to Q = q + f and holds the interface fixed near a given position. The parameters κ and appearing in the pinning potential denote the force constant and the position of the bias, respectively. The third term, (κf/2) f2, reflects the harmonic potential energy experienced by phonon fluctuations. Finally, the last term is the kinetic energy associated to the time derivative of the phonon degree of freedom f, which carries the effective mass mf. A kinetic energy term for the coordinate q is not necessary, because the dynamics of q is assumed to be overdamped.
B. Langevin equations of motion
Based on the effective Hamiltonian of Eq. (7), we next devise stochastic equations of motion to describe the coupled time evolution of q and f. Since, the growth of the interface separating the liquid from the crystalline phase is slow on the molecular time scale, inertial effects in the motion of the variable q are negligible, suggesting to describe its dynamics with an overdamped equation. In contrast, lattice vibrations are fast, which requires an underdamped description for the dynamics of f. Thus, on a coarse-grained level, we expect the dynamics of q and f to be governed by the following pair of coupled Langevin equations,
Here, γ and γf are friction constants associated with q and f, respectively, and ηq(t) and ηf(t) are δ-correlated Gaussian random forces. Both the friction and random forces reflect the degrees of freedom neglected in the simplified model.10 The variance of the random forces and the friction constants are related by the fluctuation-dissipation theorem,6 〈ηq(0) ηq(t)〉 = 2γkBTδ(t) and 〈ηf(0) ηf(t)〉 = 2γfkBTδ(t), where kB is the Boltzmann constant, T is the temperature, and the angular brackets 〈⋯〉 indicate a thermal average. Solving the above Langevin equations yields q(t) and f(t) from which the time evolution of Q(t) = q(t) + f(t) can be computed.
C. Solution of the stochastic model
To apply the model to simulation data, we first consider the autocorrelation function 〈ΔQ(0)ΔQ(t)〉, where ΔQ(t) = Q(t) − 〈Q〉. As shown in details in the Appendix, this autocorrelation function is a sum of three complex exponentials with rather complicated arguments. To obtain a practical expression, we consider the long-time limit of the model and utilize the “separation of time scales”-assumption γ(1/κ + 1/κf) ≫ γf/κf,
where A = κf/(κ + κf) is the strength and is the characteristic time of the terminal relaxation. Thus, the friction constant of crystal growth can be determined from a fit to the terminal relaxation time as
The kinetic coefficient of crystal growth M in Eq. (1) is inversely proportional to the friction coefficient γ. The proportionality constant depends on the specifics of the simulations cell used. It is found by first taking the average time-derivative of the definition of xs schematically depicted in Fig. 1, using that and ,
Since μsl is known from Eq. (5), the averaged crystal growth rate can also be deduced.
It is easier to write up the solution of Eqs. (8) and (9) in the frequency domain. Let us consider the complex admittance, i.e., the frequency dependent mobility μ(ω) of Q. From the fluctuations of Q(t), it is possible to compute μ(ω) using the fluctuation-dissipation theorem6,10 (see the Appendix)
The complex admittance of the model is given by
A derivation of this relation is provided in the Appendix. As shown in Subsection III D, this result can also be obtained using an analogy to an electric circuit that is mathematically equivalent to the stochastic model.
D. Electric network representation
A system of coupled linear equations of motion can be represented as an idealized electrical circuit. The circuit shown in Fig. 2 is equivalent to our stochastic model and we can use the rules of electrical networks to derive the solution. Order parameter values of Q, , and f are represented by electric charges while voltage drops correspond to forces acting on these order parameters. The electric circuit consists of two capacitors with the admittances iω/κ and iω/κf, two resistors with the admittances 1/γ and 1/γf, one inductor with the admittance −i/ωmf, and one battery with a voltage α. Additionally, the resistors include thermal Johnson-Nyquist noise ηq(t) and ηf(t). The Langevin equations of motion, Eqs. (8) and (9), are retrieved by applying Kirchhoff’s voltage law to the left and the right loops, respectively. The effective Hamiltonian Eq. (7) corresponds to the electric energy of the circuit.
Electric network model corresponding to the effective Hamiltonian and the Langevin equations for the interface motion, Eqs. (7), (8), and (9).
The frequency-dependent admittance (Eq. (14)) over the central capacitor can be obtained as follows: The admittances of the unconnected elements from left to right are μq(ω) = 1/γ, μQ(ω) = iω/κ, and μf(ω) = 1/(iωmf + γf − iκf/ω). For the latter, we used Ohm’s law: impedances (i.e., inverse admittances) connected in series are summed. When currents are running in parallel, however, the admittances are summed. With respect to the current though the central capacitor, the loops on the left and on the right are connected in parallel with each other and in series with the central capacitor. The frequency dependent admittance, we are interested in is, thus,
By inserting the admittances μq(ω), μf(ω), and μQ(ω), we arrive at Eq. (14).
IV. THE LENNARD-JONES SYSTEM
To validate the use of the interface pinning method to compute crystal growth rates, we investigate a Lennard-Jones system11 of N = 5120 particles (8 × 8 × 20 face centered cubic unit cells; solid-liquid interface in the 100 plane; see Ref. 8 for details). The potential part of the Hamiltonian is with u(r) = 4ε((σ/r)12 − (σ/r)6) − 4ε((1/2.5)12 − (1/2.5)6) for r < 2.5σ and zero, otherwise. Trajectories are generated with the LAMMPS software package12 using the Verlet integrator with a time step of where m is the particle mass. The NpzT-ensemble is realized using the Nose-Hoover thermostat13,14 with a coupling time of and the Parrinello-Rahman barostat15 with a coupling time of .
First, we compute μ(ω) using Eq. (13). The integral is evaluated numerically from a discrete time series Qi representing Q(t). The rate is computed from central difference of Qi: . Since the autocorrelation function is zero at t = 0 and t → ∞, we can replace the Fourier-Laplace transform with a regular Fourier transform, and use the efficient Fast Fourier-Transform (FFT) algorithm to evaluate the integral numerically (we note that the FFT algorithm assumes a periodic dataset and we would get an erroneous result if the integrand did not have the property of vanishing values in the limits). For the analysis, we ensure that the discrete Qi trajectory has a high sampling frequency so that aliasing is avoided (), and are sufficiently long so that slow interface fluctuations are represented (). The dots in Fig. 3 show the real and imaginary part of the computed μ(ω) at T = 0.8ε/kB and the coexistence pressure pm = 2.185ε/σ3.8 The solid lines show μ(ω)/ω of our stochastic model with the four parameters γ, γf, κf, and mf determined by a least square fit to the real part. The agreement with both the real and imaginary part, which was not used for the fit of the parameters, is excellent. The fit gives . Here, and throughout the paper, numbers in parenthesis indicate statistical uncertainties on the last digit estimated by dividing runs into statistically independent blocks.16 Fig. 4 validates that the terminal relaxation time of the autocorrelation function 〈ΔQ(0)ΔQ(t)〉 agrees with this γ, see Eq. (10). Using Eq. (12) with N = 5120, X = Y = 12.82σ, , and Qsl = 56.16,8 we arrive at a kinetic coefficient of . We note that the model does not give a perfect fit in the high-frequency part. This is, however, not important for the estimate of γ, as discussed in Sec. VI. Fig. 5 shows that changing κ over three decades yields consistent results with the model parameters determined from the data shown in Fig. 3.
Real and imaginary part of μ(ω)/ω in an interface pinning simulation with κ = 0.5ε. The solid lines correspond to the model shown in Fig. 2 with , κf = 1.22(3)ε, , and mf = 0.019(4) m. The parameters are determined by a least square fit to the real part of μ(ω)/ω. The inset zooms in on the high frequencies peak of the real part.
Real and imaginary part of μ(ω)/ω in an interface pinning simulation with κ = 0.5ε. The solid lines correspond to the model shown in Fig. 2 with , κf = 1.22(3)ε, , and mf = 0.019(4) m. The parameters are determined by a least square fit to the real part of μ(ω)/ω. The inset zooms in on the high frequencies peak of the real part.
Autocorrelation 〈ΔQ(0)ΔQ(t)〉 from an interface pinning simulation with κ = 0.5ε at T = 0.8ε/kB and the coexistence pressure. The blue dashed line is the predicted terminal relaxation, Eq. (10), using the parameters determined by fitting to μ(ω)/ω (see Fig. 3). The blue dashed line is the full solution: −0.0248e−34.0t + 0.321e−2.68t + 0.704e−0.0226t.
Autocorrelation 〈ΔQ(0)ΔQ(t)〉 from an interface pinning simulation with κ = 0.5ε at T = 0.8ε/kB and the coexistence pressure. The blue dashed line is the predicted terminal relaxation, Eq. (10), using the parameters determined by fitting to μ(ω)/ω (see Fig. 3). The blue dashed line is the full solution: −0.0248e−34.0t + 0.321e−2.68t + 0.704e−0.0226t.
Real part of μ(ω) with κ = {0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100} computed using Eq. (13). Solid lines correspond to the model with parameters (γ, γf, κf, and mf) determined in Fig. 3.
For comparison, we compute M in a direct simulation of crystal melting (κ = 0). The average interface velocity can be computed as , where ts is the average time to crystallize one box length Zs = N/ρsXY (negative values indicate melting). In a simulation, this time can be computed from the average rate of change of the order-parameter . The solid (black) line in Fig. 6 is the average trajectory when melting a crystal at T = 0.8ε/kB and p = 1.8ε/σ3 (obtained from 50 statistically independent runs). From this result, we find that . At this state point μsl = 0.0431ε and, thus, we get a kinetic coefficient of . This is within the statistical uncertainty of the value determined by the interface pinning method.
Time evolution of Q along 50 statistically independent LJ melting runs (gray) at T = 0.8ε/kB and p = 1.8ε/kB where the initial configurations are taken from equilibrated interface pinning simulations (κ = 4ε; ). At this state-point, the liquid is the thermodynamically stable phase (μsl = 0.043ε) and the crystal melts. From the average rate of (red dashed), we compute the average growth velocity to corresponding to (using that Zs = 32.1σ). The inset shows along the T = 0.8ε/kB isotherm.
Time evolution of Q along 50 statistically independent LJ melting runs (gray) at T = 0.8ε/kB and p = 1.8ε/kB where the initial configurations are taken from equilibrated interface pinning simulations (κ = 4ε; ). At this state-point, the liquid is the thermodynamically stable phase (μsl = 0.043ε) and the crystal melts. From the average rate of (red dashed), we compute the average growth velocity to corresponding to (using that Zs = 32.1σ). The inset shows along the T = 0.8ε/kB isotherm.
V. FIRST PRINCIPLE COMPUTATIONS
We used the interface pinning method in combination with ab initio simulations to compute the crystal growth rate for real materials. In Ref. 7, we implemented the interface pinning method in the Vienna Ab initio Simulation Package (VASP)17 and conducted DFT calculations of the melting temperature of the period three elements sodium (Na) and silicon (Si). These simulations allow us to compute the kinetic coefficients of these elements in the z-direction of the chosen crystal directions.
For computing the trajectories, we employed a Verlet integrator with a timestep of 4 and 3 fs for Na and Si, respectively. The NpzT-ensemble was realized using a Langevin thermostat with a coupling time of 1 ps and a Parrinello-Rahman barostat15 with a coupling time of 0.33 ps. The Perdew, Burke, Enzerhof (PBE) exchange correlation functional was used for Na, and the local density approximation (LDA) was used for Si.18 For the analysis, we use a 120 ps trajectory for Na and 53 ps for Si. In Fig. 7, we compute the kinetic coefficients by fitting our model to the autocorrelation function 〈ΔQ(0)ΔQ(t)〉. Table I lists the computed kinetic coefficients.
The Q(t) autocorrelation in an ab initio interface pinning simulation of Na (upper panel) and Si (lower panel). The solid lines are least square fits to the stochastic model (see the Appendix). The dashed lines indicate the terminal relaxation estimated by Eq. (10).
The Q(t) autocorrelation in an ab initio interface pinning simulation of Na (upper panel) and Si (lower panel). The solid lines are least square fits to the stochastic model (see the Appendix). The dashed lines indicate the terminal relaxation estimated by Eq. (10).
First principle kinetic coefficients of Na and Si.
Element . | Direction . | Unit cells . | T (K) . | N . | . | Z (Å) . | Qsl . | κ (eV) . | M (Å/(ps eV)) . |
---|---|---|---|---|---|---|---|---|---|
Na | bcc(100) | 5 × 5 × 10 | 400 | 500 | 39.86(4) | 43.04(2) | 0.31(2) | 500 | 250(50) |
Si | cd(100) | 3 × 3 × 6 | 1400 | 432 | 20.71(2) | 32.96(1) | 0.44(1) | 1000 | 30(15) |
Element . | Direction . | Unit cells . | T (K) . | N . | . | Z (Å) . | Qsl . | κ (eV) . | M (Å/(ps eV)) . |
---|---|---|---|---|---|---|---|---|---|
Na | bcc(100) | 5 × 5 × 10 | 400 | 500 | 39.86(4) | 43.04(2) | 0.31(2) | 500 | 250(50) |
Si | cd(100) | 3 × 3 × 6 | 1400 | 432 | 20.71(2) | 32.96(1) | 0.44(1) | 1000 | 30(15) |
VI. DISCUSSION
A. Dependence on stochastic model
We have devised a stochastic model of Q(t) fluctuations in an interface pinning simulation and used it to compute the crystal growth velocity. The critical reader might rightfully ask the question: How dependent is the method on the validity of the stochastic model? To answer this, we emphasize that it is the long-time ω → 0 limit of the model that is of importance. If we use a more complicated model of the fast f(t) fluctuations, the long-time solution is the same, except that κf is replaced by an effective spring constant. To motivate the use of a more complicated model, we note that the high-ω part of the μ(ω) spectrum (inset on Fig. 8) indicates a superposition of two peaks. The reason for this is that fluctuations in the bulk of the liquid and the crystal are different. Thus, a better model of Q(t) fluctuations would be to split the fast contribution into contributions of the two phases, f(t) = fl(t) + fs(t), and write the effective equations of motion as
The frequency dependent mobility (admittance) of this model is
The long-time limit of this model, however, is the same as in the simpler model we used above, except that κf is replaced by an effective spring constant [1/κl+1/κs]−1.
The kinetic coefficient M along the LJ melting line computed with the interface pinning method.
The kinetic coefficient M along the LJ melting line computed with the interface pinning method.
B. Other techniques
How does the interface pinning method compare to other methods for computing M? To answer this we note that suggested techniques can be grouped into classes:4 (i) free solidification simulations (see Fig. 6),19–21 (ii) fluctuations analysis of a crystallinity order-parameter,22–24 (iii) capillary wave analysis,4,25–27 and (iv) forced velocity simulations.28,29 This paper’s method is of class (ii). The interface pinning method can be viewed as a generalization of the approach suggested by Briels and Tepper.22,23 In the Briels-Tepper method two-phase configurations are stabilized by simulating the constant NVT ensemble. This can be viewed as a special case of the interface pinning method, where the order-parameter is the volume, and the bias potential has an infinitely large spring constant. Fluctuations in the number of crystalline particles can be monitored by pressure fluctuations or an order-parameter of crystallinity (e.g., Eq. (3)). An advantage that the interface pinning method inherits from the Briels-Tepper method is that it involves well-defined equilibrium computations that can be done ad infinitum. The following challenges of the Briels-Tepper method are solved with the interface pinning method: (a) In the Briels-Tepper method, the two-phase configurations are stabilized by keeping the volume constant. In effect, computations can only be done near coexistence. In the interface pinning method, the stabilization of two-phase configurations is done by connecting the system with a harmonic field that couples to an order-parameter Q that distinguishes between the crystal and the liquid, and simulations can be performed far into the supercooled or superheated regimes. Also, the size of the fluctuations in the number of crystalline particles is determined by the compressibility of the solid and liquid. In effect, they will be large for large systems, leading to long correlation times. With the interface pinning method, the size of the Ns fluctuations can be controlled by the value of κ and by choosing a good order-parameter making κf large. (b) μsl is directly computed with the interface pinning method. This is an essential property that otherwise had to be computed in separate calculations.
The interface pinning method also bears some similarities with the capillary wave analysis method4,25–27 as both approaches relate the interface mobility to the time correlations of spontaneous interface fluctuations based on the fluctuation-dissipation theorem. More specifically, in the capillary wave analysis method, one tracks the position and shape of the crystal-liquid interface and determines its k-space fluctuation spectrum. By combining the decay time of the correlations with the interface stiffness obtained from static correlations, one can then compute the interface mobility. The interface pinning method and the capillary wave analysis method, however, also differ in important ways. While in the interface pinning method, one considers only fluctuations of the global parameter Q(R) that quantifies the fraction of the system in the crystalline state, the capillary wave analysis method requires an exact localisation of the crystal-liquid interface. To do that one needs to be able to differentiate locally between solid-like and liquid-like particles, a distinction that is not necessary in the interface pinning method. Also, the capillary wave analysis method can be applied only to coexisting phases, thus, requiring knowledge of the melting temperature, while the interface pinning method provides the interface mobility away from coexistence too.
C. Latent heat and volume
In this paper, we consider the hydrostatic NpzT-ensemble7,8 near coexistence. Growth or melting of a solid results in latent heat and latent volume that must be removed to stay in the NpzT-ensemble. To avoid that the growth rate is trivially determined by the characteristic time of the thermostat and barostat, we investigated crystal growth dynamics with a strong coupling of the thermo- and barostat, i.e., short coupling times. We note that in a weak coupling or experimental situation latent heat and volume disperse away from the interface. In effect, the temperature and pressure at the interface are different from that at the boundaries of the system,24,29 and in general the temperature and strain tensor fields should be taken explicitly into consideration. Addressing this issue goes beyond the scope of this paper.
D. Computing crystal nucleation rates
Auxiliary harmonic potentials as the ones employed here are routinely used to compute free energies along reaction coordinates for crystal nucleation or other rare events. Specifically, the umbrella sampling method30,31 uses a series of auxiliary potentials referred to as umbrellas. By reweighting32 probability distributions generated in the presence of the umbrellas, it is possible to compute free energies along a given reaction coordinate. As we have demonstrated here, one can also extract crystal growth rates without further computations. By extending the bias potentials to the regime where only a small crystallite is present in the system, one may be able to compute crystal nucleation rates from IP-simulations. We leave such studies for future investigations.
E. Scaling of the kinetic coefficient along the Lennard-Jones melting line
Fig. 8 shows the kinetic coefficient M along the coexistence line. In the region T < 1.2, M is nearly constant. At higher temperatures, M decreases indicating that the crystal growth (or melting) rate is lower. The inset shows that M in reduced units, using the thermal energy for the energy scale, is roughly invariant at T > 1.2. In other words, M scales as along the high temperature (and pressure) part of the melting line. This scale invariance is expected if a mapping to the hard-sphere system is valid26 and is predicted by the “isomorph theory”33 of simple liquids.34 This theory states that there is a class of systems,35–37 including the LJ system,35 that have “isomorph”-lines in the dense and/or high temperature part of the phase diagram (i.e., not near the critical point). Along these lines structure, dynamics and some thermodynamic properties are invariant in reduced units. One prediction is that the melting line is an isomorph,33 and indeed this was shown for the LJ system in Refs. 8 and 38. Thus, we expect M to be invariant in reduced units along the melting line. This explains why M scales as . The scaling law does not apply near the triple point (note that the lowest temperature data-point on Fig. 8 is at negative pressure). Consistent with this, the melting line itself deviates from an isomorph near the triple point.8 This deviation is probably due to long-ranged attractive interactions.8 We will leave further investigation of isomorph-scale invariance of the crystallization to future studies.
Acknowledgments
This work was financially supported by the Austrian Science Fund FWF within the SFB ViCoM (F41). U.R.P. was funded by the Villum Foundation’s Grant No. VKR-023455. The computational results presented have been achieved using the Vienna Scientific Cluster (VSC). The authors are grateful for valuable comments and suggestions from Georg Kresse, Jeppe C. Dyre, Lorenzo Costigliola, Tage Christensen, and Peter Harrowell.
APPENDIX: ANALYSIS OF STOCHASTIC MODEL
In this appendix, we provide a detailed analysis of the stochastic model suggested in the main part of the paper.
1. Static averages and variances
The static averages and fluctuations of q and f (and hence also of Q) can be computed as averages over the equilibrium distribution
where β = 1/kBT and the partition function
normalizes the distribution. The moments of this multivariate Gaussian distribution can be determined analytically yielding the averages
and variances
Here, Δq = q − 〈q〉, Δf = f − 〈f〉, and are the deviations of q, f, and from their respective averages. While is uncorrelated to the other variables, q and f are correlated with covariance
From the above expressions, the average and variance of Q = q + f follow
Thus, the static fluctuations of Q depend only on the temperature and the force constant κ of the pinning potential.
2. Time correlation functions
To quantify the average dynamics of the model, we now introduce the time correlation functions
which correlate the state of the system at time 0 to its state at a time t later. It follows from the microscopic reversibility of the equations of motion that ϕqf(t) = ϕfq(t), such that we do not need to consider ϕqf(t) and ϕfq(t) separately. Thus, the time autocorrelation functions of the observable order parameter 〈ΔQ(0)ΔQ(t)〉 are
One can show, that the following equations hold for the time correlation function 〈α(0) β(t)〉, where α and β are two arbitrary functions of q and f,
Upon time reversal, the time correlation function transform as
where εα and εβ take values +1 or −1, depending on whether α and β, respectively, are even or odd under reversal of the momenta.
By averaging over the equations of motion, the following differential equations for the time correlations can be derived:
where we have omitted the argument of the time correlation functions for simplicity. Introducing the auxiliary function , we write these differential equations in matrix notation as
In a more compact form, this system of homogeneous linear first order differential equations with constant coefficients is expressed as
where and A is the constant 4 × 4 matrix on the right hand side of Eq. (A22). The formal solution of this equation is given by
with initial conditions
As can be seen in Eq. (A22), the time evolution equations for ϕqf(t), ϕff(t), and ψff(t) are independent of ϕqq(t). Hence, one can obtain the time correlations functions ϕqf(t), ϕff(t), and ψff(t) by solving
As a consequence, the three time correlation functions ϕqf(t), ϕff(t), and ψff(t) can be written as superpositions of three exponentials,
where xi(t) is component i of the time correlation function vector x(t). (Note that we start our numbering at 0 such that the 0-th component of x(t) is the time correlation function ϕqq(t)). The time constants λi are the eigenvalues of the 3 × 3 matrix on the right hand side of Eq. (A26) and the specific values of the constants aij can be determined by diagonalising this matrix and applying the initial conditions ϕqf(0), ϕff(0), and ψff(0). The eigenvalues λi, which are also eigenvalues of the matrix A and can be complex for certain parameter combinations, are the roots of the cubic equation
Real roots correspond to exponential decay, while complex roots lead to oscillatory behaviour. In general, the time correlation functions ϕqq(t), ϕqf(t), ϕff(t), and ψff(t) are superpositions of exponential and oscillatory terms. Whether the roots of this equation are real or complex can be determined by computing the discriminant D of the above polynomial in λ. This polynomial discriminant, as well as the roots of the equation, can be expressed explicitly in terms of the constants κ, κf, γ, γf, and mf, but the expressions are omitted here because they are complicated and do not provide much insight. If D > 0, one root is real and two are complex conjugate, if D < 0 all roots are real and different, and if D = 0, all roots are real and at least two of them are equal.
Once the time correlation function ϕqf(t) is known, ϕqq(t) can be determined by solving the remaining differential equation
Since ϕqf(t) is already given, this equation is an inhomogeneous linear differential equation with constant coefficients, which can be solved, for instance, by variation of constants yielding
Here, λ0 = κ/γ and the coefficient a0 depends on the initial conditions. Like λ1, λ2, and λ3, also λ0 is an eigenvalue of the matrix A. Thus, in general, the time correlation function ϕqq(t) is a sum of four exponentials, two of which can lead to oscillatory behaviour. However, explicit solutions of the differential equations for specific parameter sets indicate that the coefficient a0 is close to zero, such that also ϕqq(t) is, effectively, a sum of three exponentials.
While it is possible to compute the roots of Eq. (A28) analytically, the resulting expressions are exceedingly complicated. For the root closest to zero, which dominates the behaviour of the correlation functions for long times, a simple but accurate approximation can be derived. Neglecting the cubic and quadratic term, the root with the smallest magnitude is given by
As a further simplification we will use a “separation of time-scales” approximation to eliminate the last term. Thus, the terminal exponential relaxation time is
justifying Eq. (10) in the main part of the paper. We note that the term τf = γf/κf that we neglect is the characteristic time of the uncoupled (κ = 0) fast f vibrations. For the typical system studied in this paper, this term is less than 1% of τ.
3. The fluctuation-dissipation theorem and the solution in the frequency domain
According to the fluctuation-dissipation theorem,6,10 the response of a system to a weak perturbation can be related to the fluctuation properties of the system in equilibrium. To apply this theorem to the pinned interface, imagine a system with the Hamiltonian perturbed by a time dependent external field K(t) that couples linearly to the variable A(z). Here, z denotes a set of variables describing the state of the system, for instance, the positions and momenta of all particles in our simulation or the variables Q, q, f, and of our stochastic model. The time-dependent Hamiltonian of the perturbed system is then given by
The reaction of the system to the external perturbation can be monitored through the change in the variable B(z), which, like A, is also a function of z, with respect to its equilibrium average,
where the angular brackets denote an equilibrium average and B(t) is a shorthand notation for B[z(t)]. Assuming that the external force K(t) has been acting on the system from the infinite past, i.e., from t = − ∞, the average deviation 〈ΔB(t)〉ne at time t can be written in its most general form as
where the response function RBA(t) denotes the response of the system to an impulsive force, i.e., K(t) ∝ δ(t) applied to the system at time t = 0. The angular brackets 〈⋯〉ne imply an average over many realisations of the process. According to the fluctuation dissipation theorem, the response function is related to the fluctuation properties of the system by
where ΔA(t) = A(t) − 〈A〉 denotes the deviation of A(t) = A[z(t)] from its equilibrium average. Hence, in the linear regime, the response of a system to an external perturbation is related to the correlation of spontaneous fluctuations at different times in the equilibrium system.
The significance of the fluctuation dissipation theorem becomes particularly clear in the frequency domain. In this case, the linear relation between perturbation and response is expressed as
where
and
are the Fourier transforms of the response and the force, respectively. In the frequency representation, the fluctuation-dissipation theorem then links the complex admittance χBA(ω) to the Fourier-Laplace transform of the time correlation function ,
In the following, we will mainly use this frequency-dependent formulation of the fluctuation-dissipation theorem.
We next apply the fluctuation-theorem to our stochastic model with the Hamiltonian of Eq. (7) and evolving according to the Langevin equations (8) and (9). For this model, we will compute the frequency dependent complex admittance that describes the reaction of the system to an external perturbation acting on the variable Q = q + f. The response of the system to the perturbation is monitored using , the velocity associated with the variable Q. So what we would like to determine is the complex admittance μ(ω) that relates the average velocity to the strength of the external perturbation
By comparing the computed admittance μ(ω) with results of computer simulations carried out for the atomistic system, we will then determine the values of the model parameters. In particular, we will determine the parameter γ which describes the mobility of the interface driven by the difference in chemical potential between the phases.
Since these equations are linear, we can compute the response of the system analytically for an arbitrarily strong external force K(t). To do that, we carry out a Fourier transformation on the Langevin equations above and average over many realizations of the stochastic process. For ω ≠ 0, we obtain
where and are the Fourier transforms of Δq(t) = q(t) − 〈q〉 and Δf(t) = f(t) − 〈f〉, respectively. Note that in the above equations we have omitted the argument ω in the averages to simplify the notation. Since, we are interested in the response of the system in terms of the generalized velocity , we rewrite these equations for the Fourier transforms of the time derivatives of Δq and Δf,
Here, we have exploited that in frequency space taking a time derivative simply amounts to multiplication with iω, such that and .
To simplify the notation in the following, we introduce the complex admittances μq(ω) and μf(ω) for and separately in the absence of the coupling term in the Hamiltonian of Eq. (7), i.e., without pinning potential. In this case, which corresponds to κ = 0, the equations of motion yield
such that
By equating the right hand sides of the above equations, one obtains
implying that
The response of the system in terms of is then given by
such that the complex admittance μ(ω) is given by
equivalent with Eq. (14) in the main part of the paper. The fluctuation-dissipation theorem links the complex admittance to the Fourier-Laplace transform of the equilibrium time correlation function of ,
Thus, the complex admittance μ(ω) can be obtained by Fourier-Laplace transformation of the time correlation function and vice versa with an inverse Fourier-Laplace transformation.