# Microscopic Rate Constants of Crystal Growth from Molecular Dynamic Simulations Combined with Metadynamics

Authors: Dániel Kozma; Gergely Tóth
Journal: Advances in Physical Chemistry
 (2012)

---

## Abstract

Atomistic simulation of crystal growth can be decomposed into two steps: the determination of the microscopic rate constants and a mesoscopic kinetic Monte Carlo simulation. We proposed a method to determine kinetic rate constants of crystal growth. We performed classical molecular dynamics on the equilibrium liquid/crystal interface of argon. Metadynamics was used to explore the free energy surface of crystal growth. A crystalline atom was selected at the interface, and it was displaced to the liquid phase by adding repulsive Gaussian potentials. The activation free energy of this process was calculated as the maximal potential energy density of the Gaussian potentials. We calculated the rate constants at different interfacial structures using the transition state theory. In order to mimic real crystallization, we applied a temperature difference in the calculations of the two opposite rate constants, and they were applied in kinetic Monte Carlo simulation. The novelty of our technique is that it can be used for slow crystallization processes, while the simple following of trajectories can be applied only for fast reactions. Our method is a possibility for determination of elementary rate constants of crystal growth that seems to be necessary for the long-time goal of computer-aided crystal design.


---

## Body

## 1. Introduction

The macroscopic shape of a crystal reflects the net effect of complex processes. The number of the important processes and the dependence of the results on many factors mean a large task for theoretical description and computational modelling. The most important factors are the properties of the solute and the solvent, the temperature, the pressure, the impurities, and the form of the possible seeds. Crystallization is usually divided into two steps, the formation of the seeds and the growth of the crystal [1]. There are numerous ways to model both steps. A possible grouping of the methods is according to the size scale of the methods terming as microscopic, mesoscopic, and macroscopic ones. In detailed modelling, we can use quantum mechanics where the electrons of the system are explicitly included. If the effect of the electrons is averaged in classical mechanical interaction potentials among the atoms or molecules, we can use atomic/molecular modelling with classical mechanics. We can use mesoscopic lattice models, where the atoms or molecules occupy lattice points defined in advance. We may even simplify the model up to limit, where the model includes only the information, that a lattice point is occupied or it is vacant. If the basics of the model are planes of the crystal, we approach to macroscopic description of the processes [2]. In these macroscopic cases, the increase or decrease of the interfacial planes is mostly described by differential equations or by stochastic methods, and thermodynamic and kinetic properties of the crystal planes are the parameters of the model. Recently, classical density functional method is used in the simulations as well [3, 4].Despite the recent development of computational power, the microscopic modelling is still not feasible for complex systems in the case of crystal growth from solution. Here, three points have to be fulfilled simultaneously: accurate description of the system, reasonable simulation time, and large system size. Most crystallization processes need longer simulation time than the accessible 10−9–10−6 s in atomistic approaches. Although there are few systems [5–7] that can be reasonablly studied within this time domain, the most practically important crystallizations take hours or days in the reality. Considering the size limitation, quantum mechanical calculations are restricted to systems, where only very simple interfaces can be found, and the solution phase can be highly biased by artificial periodicity.These model systems are usually too small to be treated as an adequate representation of the interfacial problem [1]. In the case of classical mechanical models, it is seldom possible to include the whole crystal with the surrounding solution into the simulation. Usually, only the modelling of clearly defined interfaces are performed one at a time and some interfacial structures, like kinks, edges, or even screw dislocations are included. While the time domain of these simulations can be extended by crude simplification of the classical mechanical interaction potential, it leads to a wrong representation of reality. In the case of the classical density functional method, up to now, the system consists of only simple model particles, and the application on realistic molecules has not been elaborated yet [3, 4].The time and size range restrictions can be eliminated, if we use mesoscopic methods. A feasible method is the kinetic Monte Carlo simulation [1, 8–11]. Here, elementary kinetic processes were performed at randomly selected crystal positions of a supercell. The supercell is constructed from the three-dimensional periodic images of the unit cell. The supercell contains hundred thousands to several millions of unit cells. Initially, a small portion of the lattice points in the middle of the supercell are occupied creating a seed (crystal phase), while the other lattice points are unoccupied (solution phase). We take into consideration two elementary reactions, a crystallization process where an unoccupied lattice point changes to be occupied, and a dissolution process, where an occupied lattice point becomes unoccupied. If these processes are performed at randomly selected positions of the supercell with a probability reflecting the rate constant of the process in the given lattice point environment, we can reasonablly simulate the growth of the crystal. The exact knowledge of the rate constants is crucial. Using the transition state theory, we can calculate these values from the free energy surface of the system. Free energy calculations are not straightforward, but they can be approximated by potential-energy-surface calculations extended with vibration data of the important phase-space points. In the case of quantum mechanical calculations, we are able to calculate it for a few phase space points, but we need a very good guess, which are the important ones. There are different methods to determine the transition pathways, but their adaptation to complex systems containing both crystal and solution is not trivial. The given kinetic parameter set is seldom complete, because the small systems are inadequate representations of the reality. In the case of classical mechanical models, where reasonable system size can be used, the vibration calculations are not trustworthy, because classical mechanical models are seldom suitable for accurate vibration calculations.Piana et al. [5–7] applied another way to determine the rate constants. They performed classical mechanical molecular dynamics on different interfaces of urea crystal and urea solution in water and in methanol. They count the number of the elementary processes (crystallization and dissolution) at different interfacial positions. They differentiated the positions by the occupation of the neighbouring lattice points. If we take into account only the first 6 hydrogen-bonded neighbours in a urea crystal, and we use some symmetry considerations, there are 36 different positions. It means 72 rate constants for the two elementary processes. Piana et al. counted the visually detectable crystallization and dissolution processes for each type of positions. They got rate constants with acceptable statistics, because the crystal growth of urea is fast. Thereafter, they performed kinetic Monte Carlo simulations with the rate constant sets for both solvents (water, methanol). They obtained micrometer size crystals with macroscopic shapes corresponding well to experimental data. Their investigation clearly showed in a numerical experiment, that microscopic interactions determine the macroscopic shape of the crystals. The reverse direction was investigated in our previous paper: how can we estimate the microscopic rate constant on the basis of the macroscopic shape of the crystal [12]. Using a combined technique of kinetic Monte Carlo simulation and genetic algorithm for parameter estimation, we obtained kinetic rate constant sets for urea in water and in methanol. The sets were qualitatively similar to the ones of Piana et al. Our investigation showed that macroscopic shape contains important and relevant information to clarify the microscopic crystal growth processes. There is a lack of complete rate constant sets on crystal growth except the studies of Piana et al. and our reverse work in the literature. The absence of data does not mean disinterest, because these sets are maybe mandatory to develop computer aided crystal design on mesoscopic scale. The lack of the data is due to the difficulties of obtaining the rate constants. Up to now, these microscopic rate constants cannot be determined experimentally. The method of counting the events in a molecular dynamic simulation can be performed for very fast crystal growth processes, but the most of the scientifically and technically important processes are slow. The parameterization on the macroscopic shape supplies uncontrollable data, because the physical background is missing from our multiparameter fitting algorithm. It turned out to be feasible for urea, but it should be tested on other systems as well. The determination of the potential energy surface augmented with vibration data is usable for a few phase space points and for the corresponding few rate constants, but the determination of all rate constants exceeds the actual computation capacity. Furthermore, the systems are very small in the quantum mechanical calculations.During the preparation of our paper, two interesting papers have been published on combined techniques to identify microscopic rate constants of crystallization processes [13, 14]. Liu et al. used classical and quantum mechanical simulations with metadynamics and transition path sampling to investigate the free energy surface for the initial steps of dissolution of NaCl in water [13]. They focused on the possible reaction pathways. Stack et al. [14] described combination of MD simulation and rare event methods to find the reaction pathway and rate constants of attachment and detachment of a barium ion onto a stepped barite surface. They determined the rate constants of elementary steps, like kink to bidentate, bidentate to inner shell, inner shell to outer shell, and outer shell to dissolution processes.In our paper, we present a similar way to determine the rate constants of crystal growth in microscopic simulations, but here we perform a kMC simulation also using our microscopic rate constants. We combined classical mechanical molecular dynamics with metadynamics [15] to explore the free energy surface of crystal/liquid interfaces. Our method keeps the advantage of classical mechanical simulations, that the simulated system is large enough to be adequate representation of the interface. It explores directly the free energy surface: there is no need of classical mechanical vibration calculations. The main issue is that our method can be used for slow processes as well, because it uses the transition state theory as opposed to the counting method of Piana et al. In this feasibility study, we applied the method on the liquid/crystal interface of liquid argon modelled by Lennard-Jones interaction. The crystallization of Lennard-Jones systems was studied several times with computer simulations [16–18], but their aim was different than the determination of rate constants. This one-component study is only our first step. Here, we were able to mimic temperature-induced crystallization instead of the more important concentration-induced cases. The rate constants obtained in the combined molecular dynamics and metadynamics method were used as inputs of mesoscopic kinetic Monte Carlo simulations. We set up this simplified system only for studying the usability of our method, but unfortunately it turned to be more complicated than the original concentration-induced two-component situation. Despite the unexpected difficulties, we could show that rate constants can be determined this way.


## 2. Details of Calculation/Methods

### 2.1. Molecular Dynamics

The basic system of the simulations was a three-dimensional periodic Lennard-Jones model of argon containing two stable crystal/liquid interfaces (Figure1). The crystal was hexagonal close-packed and the interfaces were the HCP(0001) ones. Each interfacial crystal atom had 6 neighbours in the plane of the interface and three neighbours in the lower layer of the crystal at equal distances. The space of the three missing crystal neighbours was filled by the liquid. The density of the phases, the temperature, the simulation box size, and the crystal parameters were set to obtain stable interfaces in equilibrium.Figure 1
Simulation cell. Atoms: grey—crystal phase, blue—liquid phase.We performed a systematic search close to the triple-point of Ar [19, 20]. The parameters were determined and optimized in several simulations, the reliability of the parameters was checked by comparison to data (Figure 2) in the literature [21]. The stability of the system was checked visually on the intact structure of the solid phase and by calculating the square of the individual displacement of the atoms in the different phases [21], because this value was characteristically larger for the atoms in the liquid phase than for those ones in the solid phase. The isotropic feature of the pressure was checked also, because the cell contained a crystalline phase. It means the three axes of the simulation cell were set to obtain the same pressure in the three directions. The system consisted of 1400 atoms (700-700 atoms in the crystal and in the liquid phases). The final simulation parameters were: T=82 K, ρcryst=0.021837 Å−3, ρliquid=0.021519 Å−3. There were 10-10 rows of atoms at the interfacial crystal in the x and y directions, and altogether 7 (4+3) layers of crystal in the zdirection. The Lennard-Jones parameters were ε=956.11 J/mol and σ=3.405 Å [19]. The equilibrium distance in the crystal was 3.822 Å between two neighbours. It corresponded to the minimum of the Lennard-Jones potential. The equation of motions was integrated with the velocity Verlet method using a time step of 10−15 s. Nose-Hoover thermostat was used to maintain the temperature of the simulations in the canonical ensemble [19]. We used an NVT ensemble instead of NPT, that would be a more sophisticated solution, but the “phase transition” of a single atom does not change the density values significantly, so using NVT is acceptable.Figure 2
State variables of our system and the coexistence curve of Ar liquid/solid phase [19]. The values are in reduced units, where T* is the reduced temperature, ρ* is the reduced density. The dotted line with × denotes the solid, the dotted line with black square denotes the fluid system. Δ and ○ represent the state variables of our fluid and solid phase.The equilibrated basic system was used in the construction of different interfacial structures. Usually, atoms were taken out from the crystalline part at the interface to create vacancies, kinks, and steps. The density was maintained by addition of particles onto the other solid/liquid interface.


### 2.2. Metadynamics

Metadynamics is a numerical technique to determine free energy surfaces in computer simulations [15]. We define an s collective variable, which depends on the phase space variables of the system. The dynamics of the system is biased by auxiliary potentials placed on the trajectory of the collective variable. These cumulatively introduced potentials repel the system of the previously occupied regions of the collective variable. The dynamics of the system are artificial, but the auxiliary potentials contain information on the free energy surface of the system. Usually, the auxiliary potentials are Gaussian functions centered on the average value of sin a previous time interval. The functions are cumulatively added to the system at each T period. It can be shown, that the s-dependent free energy of the system can be approximated for an infinitely long simulation as
(1)F(s)=-w∑t′=T,2T,3T,…exp(-(s-s(t′))22δ2).w is the height and δ is the standard deviation parameter of the Gaussian. Metadynamics is a widely used technique for the determination of different free energy barriers and possible reaction pathways. There are many technical articles as well about the algorithm, the three fundamental articles of this method can be found in [15, 22, 23].We used metadynamics to calculate the free energy barrier of displacing an atom from the crystal to the liquid phase (Figure3). We defined the s variable as the actual Cartesian coordinates of the selected interfacial atom in the crystal. During the simulation, the repelling Gaussian potentials were added to the average position of this particle in every T period. A simulation was finished, if the chosen particle left the crystal phase and entered into the liquid phase, or it was finished after 70000 time steps that time was long enough to displace the investigated particle. We approximated the activation energy of the dissolution to the maxima of F(s).Figure 3
Stimulated displacement of an atom from the crystal to the liquid phase. Atoms: grey—crystal phase, red—atom to be displaced, black background—liquid phase.The activation energy of the crystallization was approximated by this one as well, because our two phases were in equilibrium. We note here that the obtained free energies were maybe lower limits of the activation energies, because the kinetic energy of the given particles contributed to pass the activation barrier. For the sake of simplicity, we did not apply any correction term. The parameters of the metadynamics were determined in test calculations. An important issue was that the total calculation had to be short enough to avoid spontaneous changes at the interface, like dissolution of other interfacial atoms or crystallization of atoms from the liquid phase. At first, we optimized the parameters of the auxiliary Gaussian functions. In our final choice, the height of the Gaussian (w) was 0.4 kBT and the standard deviation (δ) was 0.2 Å, where kB denotes the Boltzmann constant. The relatively large height and standard deviation of the Gaussians were necessary, because there was a fast rearrangement of the interface. The height of the Gaussian was comparable to the average kinetic energy of one degree of freedom. We note here, that despite these large values, the auxiliary potential surface was relatively smooth, because the few hundred Gaussians were put along a three-dimensional s and the probability of putting two Gaussians very close to each other was small. Anyway, to reduce the uncertainty of the average free energy barrier, we performed at least 30 simulations for every case. The height of the Gaussians was kept fixed during the simulations. Perhaps, the well-tempered metadynamics method [24] would increase the accuracy of the free energy barriers, but the possible gain would be unimportant due to the large uncertainty caused by the interface reconstruction (detailed later). The time period (T) of adding the repelling potential and the time length of averaging were determined by checking the periodic motions of the phase at the interface. We found that the periodicity of an average in cage motion of the selected atom was about 600–800 time steps. We used different setups for the optimization of the number of time steps to wait after the addition of Gaussian and also for the position averaging of the selected particle to calculate s(t′). Finally, we concluded that T=300 time steps was a feasible choice to add a new Gaussian. We did not apply waiting period, so all of the positions of the selected particles were used during this T time steps to calculate the position of the next Gaussian.The activation energy for a given simulation was defined as the maximal potential density of the auxiliary Gaussians after finishing the simulation. The points were identified with global extreme search adapted for the task. At first, we calculated at eachs(t′) points, how many other Gaussian functions could be found within 3δ radius. Then we started a randomly driven simulated annealing [25] search from the s(t′) coordinate having the highest number of points. In this way, we usually found a global or local maxima around this s(t′). We repeated the annealing process several times, and we started simulating annealing also for all other Gaussian centers, where the number of the neighbouring Gaussians was more than 90% of the maximal one. The activation energy was calculated, where the highest Gaussian density was found.


### 2.3. The Simulated Lattice Environments

The lattice points were distinguished according the occupation of the neighbouring lattice points. In the case of anisotropy, taking into account the 12 first neighbours and demanding at least one occupied neighbour position, there were 211 = 2048 cases. If we supposed isotropic symmetry relations, this number reduced to 328. For example, the 12 cases with one missing neighbour were identical and the cases with only one existing neighbour were identical, as well. In the case of two missing neighbours (or 10 missing neighbours), there were 5 different subcases. We did not calculate all of the possible cases, but we simplly selected some examples.One set of different lattice environments were created by deleting one or more neighbouring crystal atoms at the interface of the selected atom. A schematic view of this can be seen in Figure4. The second set of lattice environments were created by forming different structures at the interface (Figure 5). Some details of the systems and the calculations can be found in Table 1. Few structures of the first and the second set had identical first-neighbour occupation, but their second-neighbour occupations were not the same. This redundancy was used to check the effect of the second neighbours.Table 1
Details of the interfacial positions and activation energies. Columns are the following: interfacial structure, number of neighbours of the pushed particle, number of calculation performed on a given interfacial type, calculated average activation free energies with 95% confidence interval, median, minima, and maxima of the calculated activation free energy values, and the activation free energy values determined by the probabilistic method.
No. of neighbours
No. of calculations
〈
F
‡
〉
F
‡
median
F
‡
min
F
‡
max
F
‡
Gaussian
(J/mol + 95% conf.)
(J/mol)
(J/mol)
(J/mol)
(J/mol)
Full
9
30
4653 ± 246
4749
2976
6109
4757
Interfacial positions with vacancies at the interface
1 Vacancy
8
59
4517 ± 232
4542
2649
6620
4600
2 Vac—orto
7
30
4177 ± 357
4030
2501
5886
4259
2 Vac—meta
7
60
4405 ± 226
4389
2344
6713
2 Vac—para
7
30
4409 ± 312
4324
2369
6412
3 Vac—side by side
6
60
4412 ± 205
4484
2594
5853
4090
4 Vac—side by side
5
60
4203 ± 266
4162
1897
6756
3957
5 Vacancies
4
60
3758 ± 284
4006
1387
5603
3652
6 Vacancies
3
60
3770 ± 246
3726
1852
6193
3295
Interfacial structures
Corner minus
8
90
4127 ± 188
4235
2007
5950
Line
6
90
3488 ± 176
3446
1574
5483
Front
6
90
3293 ± 186
3114
1425
5339
Corner plus
5
90
2014 ± 96
1918
711
3423
Top
3
90
2021 ± 120
1964
847
3512Figure 4
Interfacial positions with vacancies at the interface. Atoms: light grey: upper layer in the crystal.Figure 5
Interfacial structures. Atoms: light grey—upper layer in the crystal phase at the interface, dark grey—second layer in the crystal phase, red—atom be displaced.


### 2.4. A Simple Probabilistic Model for Evaluating the Calculated Free Energies

The probability distribution of the activation free energies at a type of interfacial structure can be supposed to be unimodal, as it can be seen in Figure6 for some positions. If reconstruction of the interface occurred during the simulation, the histogram split to bi- or multimodal distribution. The different peaks might be identified as other structures, especially for the cases of the planar interfaces with different number of vacancies. The observed activation free energy histogram of an initial structure can be treated as a linear combination of the activation free energy histograms of the ideal cases of all of the structures developed during the interfacial reconstruction. If we assumed Gaussian distributions for the ideal cases of the planar interfacial structures with vacancies, the observed histograms of the full and of 1 to 6 vacancies could be described as linear combinations of the same 7 ideal Gaussian distributions. Numerically, it meant the decomposition of the curves to linear combinations of 7 Gaussian functions, where the fitted parameters were the means and standard deviations of the Gaussian functions and the linear coefficients. Our choice to use Gaussians as bases functions of the decompositions was according to its success in spectroscopic curve decompositions and that the surface rearrangement of free full and top positions looked like normally distributed.Figure 6
Free energy distributions for the different interfacial structures.There were too many independent linear coefficients; therefore, we applied a simple time-independent, probabilistic approach to reduce the number of the free parameters. The model consisted of two free parameters: probability of the occupation of a vacancy and probability of the dissolution of a neighbour. We used the formula of binomial distribution (2) to calculate the other linear coefficients using the fitted probability parameters. The shape of the equation is the same in the direction of crystallization and dissolution as well. 								(2)ζprocess(n,k)=(nk)pk(1-p)(n-k).In the case of crystallizationp denotes the probability of occupation of a vacancy and n denotes the number of vacancies, k is the number of crystallized particles. In the case of dissolution, these values are the probability of dissolution, number of occupied positions, number of dissolved neighbours. So, using this formula we can build a coefficient matrix, the upper and lower triangle contains the coefficients of the reconstructed interfaces with more or less neighbours, while the diagonal is the proportion of the original surface structure.Using the free energy values determined in the numerous calculations of each surface type, we created a normalized energy histogram between 100 and 7100 J/mol with a 200 J/mol grid size. We fitted these histograms with the corresponding linear combination of Gaussian functions in the reduced parameter space. The parameters were the means and the widths of the 7 ideal Gaussian functions and the twop-s (for occupation and dissolution). The criteria of the fitting was to minimize the sum of the squared deviation of histograms and linear combination of the pure Gaussian functions. To ensure the convergence, we predefined the order of the activation energies and kept the probability values in the positive region. We got interpretable free energy values for the interfacial structures shown in Table 1.


## 2.1. Molecular Dynamics

The basic system of the simulations was a three-dimensional periodic Lennard-Jones model of argon containing two stable crystal/liquid interfaces (Figure1). The crystal was hexagonal close-packed and the interfaces were the HCP(0001) ones. Each interfacial crystal atom had 6 neighbours in the plane of the interface and three neighbours in the lower layer of the crystal at equal distances. The space of the three missing crystal neighbours was filled by the liquid. The density of the phases, the temperature, the simulation box size, and the crystal parameters were set to obtain stable interfaces in equilibrium.Figure 1
Simulation cell. Atoms: grey—crystal phase, blue—liquid phase.We performed a systematic search close to the triple-point of Ar [19, 20]. The parameters were determined and optimized in several simulations, the reliability of the parameters was checked by comparison to data (Figure 2) in the literature [21]. The stability of the system was checked visually on the intact structure of the solid phase and by calculating the square of the individual displacement of the atoms in the different phases [21], because this value was characteristically larger for the atoms in the liquid phase than for those ones in the solid phase. The isotropic feature of the pressure was checked also, because the cell contained a crystalline phase. It means the three axes of the simulation cell were set to obtain the same pressure in the three directions. The system consisted of 1400 atoms (700-700 atoms in the crystal and in the liquid phases). The final simulation parameters were: T=82 K, ρcryst=0.021837 Å−3, ρliquid=0.021519 Å−3. There were 10-10 rows of atoms at the interfacial crystal in the x and y directions, and altogether 7 (4+3) layers of crystal in the zdirection. The Lennard-Jones parameters were ε=956.11 J/mol and σ=3.405 Å [19]. The equilibrium distance in the crystal was 3.822 Å between two neighbours. It corresponded to the minimum of the Lennard-Jones potential. The equation of motions was integrated with the velocity Verlet method using a time step of 10−15 s. Nose-Hoover thermostat was used to maintain the temperature of the simulations in the canonical ensemble [19]. We used an NVT ensemble instead of NPT, that would be a more sophisticated solution, but the “phase transition” of a single atom does not change the density values significantly, so using NVT is acceptable.Figure 2
State variables of our system and the coexistence curve of Ar liquid/solid phase [19]. The values are in reduced units, where T* is the reduced temperature, ρ* is the reduced density. The dotted line with × denotes the solid, the dotted line with black square denotes the fluid system. Δ and ○ represent the state variables of our fluid and solid phase.The equilibrated basic system was used in the construction of different interfacial structures. Usually, atoms were taken out from the crystalline part at the interface to create vacancies, kinks, and steps. The density was maintained by addition of particles onto the other solid/liquid interface.


## 2.2. Metadynamics

Metadynamics is a numerical technique to determine free energy surfaces in computer simulations [15]. We define an s collective variable, which depends on the phase space variables of the system. The dynamics of the system is biased by auxiliary potentials placed on the trajectory of the collective variable. These cumulatively introduced potentials repel the system of the previously occupied regions of the collective variable. The dynamics of the system are artificial, but the auxiliary potentials contain information on the free energy surface of the system. Usually, the auxiliary potentials are Gaussian functions centered on the average value of sin a previous time interval. The functions are cumulatively added to the system at each T period. It can be shown, that the s-dependent free energy of the system can be approximated for an infinitely long simulation as
(1)F(s)=-w∑t′=T,2T,3T,…exp(-(s-s(t′))22δ2).w is the height and δ is the standard deviation parameter of the Gaussian. Metadynamics is a widely used technique for the determination of different free energy barriers and possible reaction pathways. There are many technical articles as well about the algorithm, the three fundamental articles of this method can be found in [15, 22, 23].We used metadynamics to calculate the free energy barrier of displacing an atom from the crystal to the liquid phase (Figure3). We defined the s variable as the actual Cartesian coordinates of the selected interfacial atom in the crystal. During the simulation, the repelling Gaussian potentials were added to the average position of this particle in every T period. A simulation was finished, if the chosen particle left the crystal phase and entered into the liquid phase, or it was finished after 70000 time steps that time was long enough to displace the investigated particle. We approximated the activation energy of the dissolution to the maxima of F(s).Figure 3
Stimulated displacement of an atom from the crystal to the liquid phase. Atoms: grey—crystal phase, red—atom to be displaced, black background—liquid phase.The activation energy of the crystallization was approximated by this one as well, because our two phases were in equilibrium. We note here that the obtained free energies were maybe lower limits of the activation energies, because the kinetic energy of the given particles contributed to pass the activation barrier. For the sake of simplicity, we did not apply any correction term. The parameters of the metadynamics were determined in test calculations. An important issue was that the total calculation had to be short enough to avoid spontaneous changes at the interface, like dissolution of other interfacial atoms or crystallization of atoms from the liquid phase. At first, we optimized the parameters of the auxiliary Gaussian functions. In our final choice, the height of the Gaussian (w) was 0.4 kBT and the standard deviation (δ) was 0.2 Å, where kB denotes the Boltzmann constant. The relatively large height and standard deviation of the Gaussians were necessary, because there was a fast rearrangement of the interface. The height of the Gaussian was comparable to the average kinetic energy of one degree of freedom. We note here, that despite these large values, the auxiliary potential surface was relatively smooth, because the few hundred Gaussians were put along a three-dimensional s and the probability of putting two Gaussians very close to each other was small. Anyway, to reduce the uncertainty of the average free energy barrier, we performed at least 30 simulations for every case. The height of the Gaussians was kept fixed during the simulations. Perhaps, the well-tempered metadynamics method [24] would increase the accuracy of the free energy barriers, but the possible gain would be unimportant due to the large uncertainty caused by the interface reconstruction (detailed later). The time period (T) of adding the repelling potential and the time length of averaging were determined by checking the periodic motions of the phase at the interface. We found that the periodicity of an average in cage motion of the selected atom was about 600–800 time steps. We used different setups for the optimization of the number of time steps to wait after the addition of Gaussian and also for the position averaging of the selected particle to calculate s(t′). Finally, we concluded that T=300 time steps was a feasible choice to add a new Gaussian. We did not apply waiting period, so all of the positions of the selected particles were used during this T time steps to calculate the position of the next Gaussian.The activation energy for a given simulation was defined as the maximal potential density of the auxiliary Gaussians after finishing the simulation. The points were identified with global extreme search adapted for the task. At first, we calculated at eachs(t′) points, how many other Gaussian functions could be found within 3δ radius. Then we started a randomly driven simulated annealing [25] search from the s(t′) coordinate having the highest number of points. In this way, we usually found a global or local maxima around this s(t′). We repeated the annealing process several times, and we started simulating annealing also for all other Gaussian centers, where the number of the neighbouring Gaussians was more than 90% of the maximal one. The activation energy was calculated, where the highest Gaussian density was found.


## 2.3. The Simulated Lattice Environments

The lattice points were distinguished according the occupation of the neighbouring lattice points. In the case of anisotropy, taking into account the 12 first neighbours and demanding at least one occupied neighbour position, there were 211 = 2048 cases. If we supposed isotropic symmetry relations, this number reduced to 328. For example, the 12 cases with one missing neighbour were identical and the cases with only one existing neighbour were identical, as well. In the case of two missing neighbours (or 10 missing neighbours), there were 5 different subcases. We did not calculate all of the possible cases, but we simplly selected some examples.One set of different lattice environments were created by deleting one or more neighbouring crystal atoms at the interface of the selected atom. A schematic view of this can be seen in Figure4. The second set of lattice environments were created by forming different structures at the interface (Figure 5). Some details of the systems and the calculations can be found in Table 1. Few structures of the first and the second set had identical first-neighbour occupation, but their second-neighbour occupations were not the same. This redundancy was used to check the effect of the second neighbours.Table 1
Details of the interfacial positions and activation energies. Columns are the following: interfacial structure, number of neighbours of the pushed particle, number of calculation performed on a given interfacial type, calculated average activation free energies with 95% confidence interval, median, minima, and maxima of the calculated activation free energy values, and the activation free energy values determined by the probabilistic method.
No. of neighbours
No. of calculations
〈
F
‡
〉
F
‡
median
F
‡
min
F
‡
max
F
‡
Gaussian
(J/mol + 95% conf.)
(J/mol)
(J/mol)
(J/mol)
(J/mol)
Full
9
30
4653 ± 246
4749
2976
6109
4757
Interfacial positions with vacancies at the interface
1 Vacancy
8
59
4517 ± 232
4542
2649
6620
4600
2 Vac—orto
7
30
4177 ± 357
4030
2501
5886
4259
2 Vac—meta
7
60
4405 ± 226
4389
2344
6713
2 Vac—para
7
30
4409 ± 312
4324
2369
6412
3 Vac—side by side
6
60
4412 ± 205
4484
2594
5853
4090
4 Vac—side by side
5
60
4203 ± 266
4162
1897
6756
3957
5 Vacancies
4
60
3758 ± 284
4006
1387
5603
3652
6 Vacancies
3
60
3770 ± 246
3726
1852
6193
3295
Interfacial structures
Corner minus
8
90
4127 ± 188
4235
2007
5950
Line
6
90
3488 ± 176
3446
1574
5483
Front
6
90
3293 ± 186
3114
1425
5339
Corner plus
5
90
2014 ± 96
1918
711
3423
Top
3
90
2021 ± 120
1964
847
3512Figure 4
Interfacial positions with vacancies at the interface. Atoms: light grey: upper layer in the crystal.Figure 5
Interfacial structures. Atoms: light grey—upper layer in the crystal phase at the interface, dark grey—second layer in the crystal phase, red—atom be displaced.


## 2.4. A Simple Probabilistic Model for Evaluating the Calculated Free Energies

The probability distribution of the activation free energies at a type of interfacial structure can be supposed to be unimodal, as it can be seen in Figure6 for some positions. If reconstruction of the interface occurred during the simulation, the histogram split to bi- or multimodal distribution. The different peaks might be identified as other structures, especially for the cases of the planar interfaces with different number of vacancies. The observed activation free energy histogram of an initial structure can be treated as a linear combination of the activation free energy histograms of the ideal cases of all of the structures developed during the interfacial reconstruction. If we assumed Gaussian distributions for the ideal cases of the planar interfacial structures with vacancies, the observed histograms of the full and of 1 to 6 vacancies could be described as linear combinations of the same 7 ideal Gaussian distributions. Numerically, it meant the decomposition of the curves to linear combinations of 7 Gaussian functions, where the fitted parameters were the means and standard deviations of the Gaussian functions and the linear coefficients. Our choice to use Gaussians as bases functions of the decompositions was according to its success in spectroscopic curve decompositions and that the surface rearrangement of free full and top positions looked like normally distributed.Figure 6
Free energy distributions for the different interfacial structures.There were too many independent linear coefficients; therefore, we applied a simple time-independent, probabilistic approach to reduce the number of the free parameters. The model consisted of two free parameters: probability of the occupation of a vacancy and probability of the dissolution of a neighbour. We used the formula of binomial distribution (2) to calculate the other linear coefficients using the fitted probability parameters. The shape of the equation is the same in the direction of crystallization and dissolution as well. 								(2)ζprocess(n,k)=(nk)pk(1-p)(n-k).In the case of crystallizationp denotes the probability of occupation of a vacancy and n denotes the number of vacancies, k is the number of crystallized particles. In the case of dissolution, these values are the probability of dissolution, number of occupied positions, number of dissolved neighbours. So, using this formula we can build a coefficient matrix, the upper and lower triangle contains the coefficients of the reconstructed interfaces with more or less neighbours, while the diagonal is the proportion of the original surface structure.Using the free energy values determined in the numerous calculations of each surface type, we created a normalized energy histogram between 100 and 7100 J/mol with a 200 J/mol grid size. We fitted these histograms with the corresponding linear combination of Gaussian functions in the reduced parameter space. The parameters were the means and the widths of the 7 ideal Gaussian functions and the twop-s (for occupation and dissolution). The criteria of the fitting was to minimize the sum of the squared deviation of histograms and linear combination of the pure Gaussian functions. To ensure the convergence, we predefined the order of the activation energies and kept the probability values in the positive region. We got interpretable free energy values for the interfacial structures shown in Table 1.


## 3. Results and Discussion

A typical trajectory ofs is shown for the full structure in Figure 7. More exactly, the coordinates of the Gaussian potentials are shown corresponding to the averaged coordinates of the selected particle in the subsequent T time periods. The particles left the initial position after the placement of 50–150 Gaussians. The particles entered to the liquid phase in many cases, but rather often they squeezed out a neighbouring particle to the liquid phase and annexed the neighbouring position. Sometimes the selected particle left to the direction of the bulk crystal and destroyed locally the crystal structure. In these cases, the interface decomposed fast. We express here that the trajectories of metadynamics are not related to real trajectories. The main goal of metadynamics is the determination of the F(s) free energy surface.Figure 7
Typical “trajectory” of the Gaussians. Trajectories: black—particle to be displaced, blue—neighbouring atoms in the same crystal layer, green—neighbouring atoms in the next crystal layer, red—second neighbours in the same crystal layer.The main aim of our investigation was to develop a method to determine kinetic constants for crystal growth processes for crystallization from solution. We chose a one-component system in our feasibility study. Here, instead of investigating concentration-induced crystallization, we studied temperature-induced crystallization. We considered it would be a simplification, but unfortunately it turned to be more complicated than the original concentration-induced two-component situation. In the case of two components, the empty lattice points of the crystal at the interface are mostly occupied by solvent molecules. In our one-component case, the unoccupied lattice points could be real voids or they were filled by the liquid. If they were voids, the atoms of the liquid tended to occupy the cavities. In the case of simple spherical model of Lennard-Jones argon, it was hard to distinguish an atom at a given position, that it was a part of the crystal or the liquid phase. The occupation of voids and rearrangement of the interface seemed to be a process comparable to the time scale of our simulation with metadynamics. We got difficulties in the interpretation of our data, and these difficulties could have been avoided, if we had chosen a two-component solvent + solute system.In the case of the full interface, the activation energies were reliable. In the case of the simulations started with one vacancy at the surface, a large part of our simulations finished with different neighbourhood occupation than the initial one. We were able to remove the selected particle from the central position, but in the meanwhile the initially vacant neighbouring position was often occupied by a particle. Usually, an atom of the liquid phase filled up the hole. In the case of the planar interface with two neighbouring vacancies, the final neighbouring occupation was identical to the initial one only in 10 percent of the simulations. In the case of the other group of different interfacial structures (5 line, corner plus, etc), the main characteristics remained mostly intact during the simulations. Crystallization at the vacant positions was observed mostly only in the 5 corner minus situation. In the other cases, especially in the 5 corner plus case, the dissolution of the crystal was more probable. The structure of the interface far from the vacancies or interfacial structures remained usually unaltered. It means the differentiation of the initial liquid and crystalline phases remained still possible.We tried different simulation constraints to avoid the reconstruction of the interface, but we were unable to find a method which did not influence the numerical results drastically. Therefore, we postponed the problem to the data evaluation process. At first, we selected visually the cases, where there was no surface reconstruction. We built a sphere with 3δ around the smax value of the maximal free energy density in each simulation, and we searched the simulation time of the last added Gauss potential within this sphere. The snapshot of the particles at this time was checked for the initial condition. The simulations satisfying the condition were used in the calculation of average activation energies. The others with surface reconstruction were not used. It was a rather time-consuming method and did not provide clear cut results. The distinction of the altered/unaltered structures was fuzzy, and we got quantitatively different results depending on the subjective selecting. Therefore, we elaborated a method (see Details of Calculation Methods) on the histograms of the activation free energies. The results of the calculations are summarized in Table 1. The calculations provided data on interfacial structures, where the number of the occupied neighbour positions was between 3 and 9. The simple average of the activation free energies does not correspond well to a simplified chemical intuition based on a linear relation to the number of the neighbours. Let us discuss at first the planar interfacial structures with vacancies. Here the simple averages were not meaningfully caused by the enhanced reconstruction of the initial surface, mostly by the occupation of the vacancies. If we used medians for the estimation of the expected values, we obtained similarly biased activation free energies as for the averages. The distributions are rather broad (Figure 6) as numerically can be seen in the minimum and maximum values of the activation free energies (Table 1). If we used the Gaussian curve deconvolution described previously, we obtained meaningful activation energies for the ideal (pure) surface situations. On the other hand, also these activation energies could not be approximated by simple linear relationship to respect with the number of the neighbours. The activation process was not a clear cut energetic process based on the interaction among the neighbours. The entropic terms unpredictable differed for each structure. If we compared the case of the interfacial structures to the planar interfacial structures with vacancies, there was a large difference in the activation energies.As it was mentioned, here the interfaces remained mostly unaltered causing more relevant activation free energies. More exactly, in some cases (e.g., 5 line, front, and corner plus) the opposite process, the melting of neighbours, took place. In the case of the planar structures with vacancies, the mean activation energies were larger than the ideal ones, but here the activation energies were maybe smaller than the ideal ones.There was a significant difference among the data belonging to the same number of occupied neighbours. This difference could be easily explained by the difference in the second neighbour shells. In the case of the interfacial structures, there were missing second neighbours as well in the plane. There were lower activation energies (meaning faster crystal growth) for the crystallization or melting of partially existing planes in the parallel direction. These numerical data corresponded well to the general view of crystal growth processes. We did not perform the Gaussian curve decomposition for the interfacial structures. The necessary incorporation of the second neighbours would increase the number of the Gaussian functions up to a numerically underdetermined task.The rate constants of crystal growth and melting were calculated using the transition state theory as follows:(3)k=RThexp(-F‡RT),
where h is the Planck’s constant. We applied the simplest RT/h transmission coefficient, although it is not correct in the case of our reactions due to the liquid and solid phases and the non-elementary feature of the reactions. In the case of the kMC simulations, the difference of the opposite processes and their variability on the crystalline positions determine the shape of a crystal. If the same transmission coefficient is supposed in all rates, the magnitude of the transmission coefficient does not change the shape of the crystals, it scales only the timescale. We did not put effort in a better estimation of transition coefficient, because it would extend the study extremely. Therefore, our rates can be termed as relatively accurate rates only. If T was the temperature of the molecular dynamic simulation of the two-phase equilibrium system, we supposed that the two constants (kcryst and kmelt) were equal. If we wanted to mimic temperature-induced crystallization, we defined different temperatures for the crystallization and the melting process (i.e., in the case of cooling the solid phase). In this approach, secondary nucleations can be neglected, because at this temperature, seeds cannot be formed in the liquid phase. In Table 2, we show rate constant of a system, where the temperature of the liquid was equal to the temperature of the equilibrium system, but the temperature of the crystal was 5 K less. This setup caused that the kcryst-s were equal to the equilibrium ones, but the rate constant of melting decreased. The net effect can be interpreted as a temperature-induced crystal growth. In the calculation of the rate constants, we used the activation free energies printed in boldface in Table 1. It means for the interfacial positions with vacancies, we used the mean values of the Gaussian values of the curve resolution process. In the case of the interfacial structures, we used the average activation energies.Table 2
Rate constants for temperature-induced crystallization ofΔT=5 K.
k
cryst
k
melt
T
=
82 K (mol−1 s−1)
T
=
77 K (mol−1 s−1)
Full
9.61 · 1032
5.73 · 1032
1 Vacancy
1.21 · 1033
7.33 · 1032
2 Vac—orto
1.99 · 1033
1.25 · 1033
2 Vac—meta
1.99 · 1033
1.25 · 1033
2 Vac—para
1.99 · 1033
1.25 · 1033
3 Vac—side by side
2.56 · 1033
1.63 · 1033
4 Vac—side by side
3.11 · 1033
2.00 · 1033
5 Vacancies
4.86 · 1033
3.22 · 1033
6 Vacancies
8.20 · 1033
5.63 · 1033
Corner minus
2.42 · 1033
1.53 · 1033
Line
6.18 · 1033
4.16 · 1033
Front
8.23 · 1033
5.65 · 1033
Corner plus
5.37 · 1034
4.17 · 1034
Top
5.32 · 1034
4.12 · 1034To present the usability of the rate constants, we performed kinetic Monte Carlo simulations. The detailed algorithm of the simulations can be found in [1, 12]. It was briefly described in the introduction as well. The supercell of the kMC simulation contained 1.5625·107 lattice points. The form of the supercell was close to a cube. The initial crystal occupied 0.5% of the supercell. The initially occupied lattice points were put in the middle of the supercell in a cubic form. The simulations were finished, if the crystal reached the border of the supercell or completely dissolved.In the case of argon crystal, the number of the lattice points with different first-neighbourhood occupancy is 212. If we use symmetry relations and isotropic approximation, we might reduce it to few hundreds of different environments. We did not calculate all of the different environments in this study. We calculated some typical ones and used them for the similar lattice positions. This similarity search was performed with comparing distance histograms calculated among the occupied neighbours for all type of lattice environments. These histograms were calculated for the lattice points for which the rate constants were determined in the simulations and for the other type of lattice points as well. The rate constants of a calculated case were assigned to that uncalculated environment which provided the highest overlap integral between the two histograms.A typical kMC simulation is shown in Figure8 with the set of Table 2. The temperature-induced simulation provided a featureless (close to spherical) crystal. It was not surprising because the simple spherical interactions of liquid argon and the HCP crystal structure do not cause necessarily characteristic crystals. Furthermore, it seems to be requisite to include second neighbourhood information in the simulation, for example, to mimic the preferred in plane growth of planes. Without the incorporation of the second neighbours in the classification of the positions, it was impossible to compare our data to other results on argon crystal growth [16–18].Figure 8
Argon crystal obtained in the kinetic Monte Carlo simulation.


## 4. Conclusion

The atomistic simulation of crystal growth can be decomposed into two steps: the determination of the microscopic rate constants and a mesoscopic kinetic Monte Carlo simulation. In our study, we proposed a method to determine kinetic rate constants of crystal growth. We performed classical mechanical molecular dynamic simulation on the liquid/crystal interface of argon at the equilibrium. The technique of metadynamics was used during the simulation to explore the free energy surface of crystal growth. A crystalline atom was selected at the interface, and it was displaced to the liquid phase by adding repulsive Gaussian potentials. The activation free energy of this process was calculated as the maximal potential-energy density of the Gaussian potentials in accordance with the basic idea of metadynamics. The calculations were performed for different interfacial structures and vacancies at the neighbouring positions of the selected particle.The activation free energies were transformed to kinetic rate constants using the transition state theory. We supposed equal rate constants for the crystallization and the melting processes at a crystal position, since our system was in equilibrium. In order to mimic real crystallization, we applied a temperature difference in the transition state theory calculation of the two opposite reactions. The temperature perturbed set of kinetic rate constants was applied in mesoscopic kinetic Monte Carlo simulation of argon crystal growth.The main advantage of this determination of rate constants is that it can be used for slow processes, while the simple following of trajectories in molecular dynamics and counting of elementary reactions taken place can be applied only for fast reactions. Our method provides a possibility for routine determination of elementary rate constants of crystal growth that seems to be necessary for the long-time goal of computer aided crystal design.We should admit that our system choice was awkward for a feasibility study. We chose a one-component spherical system with high number of first neighbours both in the crystal and liquid phase. Here we were able to model crystal growth and crystal melting elementary processes of one selected particle instead of crystallization and dissolution processes. The determination of the activation free energy was not a simple task for some interfacial structures. The initial structure was often destroyed during the simulations with metadynamics. Therefore, the initial characterization of the lattice environment changed. It caused uncertainty in the interpretation of the obtained activation free energies. We think in the more important crystal growth from solutions, there is no need to take care with this bias, because in the two-component solute-solvent case, the solvent occupies the most of the liquid phase at the interface providing stable metadynamics [13]. The natural continuation of this feasibility study is the application of the method on scientifically important systems.



---

*Source: 135172-2012-08-30.xml*