Roman
Ryltsev
*ab and
Nikolay
Chtchelkatchev
cdef
aInstitute of Metallurgy, UB RAS, 101 Amundsen str., 620016, Ekaterinburg, Russia. E-mail: rrylcev@mail.ru
bUral Federal University, 19 Mira str., 620002, Ekaterinburg, Russia
cL.D. Landau Institute for Theoretical Physics, RAS, Ac. Semenov 1-A, 142432, Chernogolovka, Russia
dAll-Russia Research Institute of Automatics, 22 Suschevskaya, Moscow 127055, Russia
eMoscow Institute of Physics and Technology, Institutskiy per. 9, Dolgoprudny, 141700, Moscow Region, Russia
fInstitute for High Pressure Physics, Russian Academy of Sciences, 108840, Moscow (Troitsk), Russia
First published on 14th June 2017
Using molecular dynamics simulations, we study computational self-assembly of one-component three-dimensional dodecagonal (12-fold) quasicrystals in systems with two-length-scale potentials. Existing criteria for three-dimensional quasicrystal formation are quite complicated and rather inconvenient for particle simulations. So to localize numerically the quasicrystal phase, one should usually simulate over a wide range of system parameters. We show how to universally localize the parameter values at which dodecagonal quasicrystal order may appear for a given particle system. For that purpose, we use a criterion recently proposed for predicting decagonal quasicrystal formation in one-component two-length-scale systems. The criterion is based on two dimensionless effective parameters describing the fluid structure which are extracted from the radial distribution function. The proposed method allows reduction of the time spent for searching the parameters favoring a certain solid structure for a given system. We show that the method works well for dodecagonal quasicrystals; this result is verified on four systems with different potentials: the Dzugutov potential, the oscillating potential which mimics metal interactions, the repulsive shoulder potential describing effective interactions for the core/shell model of colloids and the embedded-atom model potential for aluminum. Our results suggest that the mechanism of dodecagonal quasicrystal formation is universal for both metallic and soft-matter systems and it is based on competition between interparticle scales.
The stability of three-dimensional (3D) one-component QCs has been theoretically predicted using density functional theory10,11 and then confirmed by molecular dynamics simulations.12–16 The results obtained in these papers suggest a general idea explaining QC formation as due to the existence of two or more interparticle length-scales. This idea is supported by the fact that effective interactions for metallic17–21 and soft matter systems22–26 are often described by multi-length-scale potentials.
A general problem of computer simulation of 3D QCs is the lack of simple geometrical criteria of QC formation. So to localize numerically the QC phase, one should usually simulate over a wide range of system parameters.13,15 A similar problem exists for complex crystal phases in systems with multi-scale interactions.13 Thus, a universal procedure allowing us to predict somehow the formation of complex solid structures (including QCs) is extremely urgent.
Recently, we have proposed a method to predict the self-assembly of decagonal QCs in one-component two-length-scale systems.14 The method suggests that the formation of QCs from the fluid phase is mostly determined by the values of two dimensionless structural parameters of the fluid. The parameters reflect the existence of two effective interparticle distances (bond lengths) originated from the two-length-scale nature of the interaction potential. These are the ratio between effective bond lengths, λ, and the fraction of short-bonded particles ϕ. It has been shown that the criterion proposed is robust under change of potential and may be applicable to any system with two-length-scale interactions.
Here we show that the criterion works well for the case of dodecagonal (12-fold) quasicrystals (DDQCs). In order to show that, we use four different two-length-scale potentials: the Dzugutov potential12 and the oscillating pair potential (OPP)13 which mimic oscillating metal interactions, the repulsive shoulder system (RSS) potential27 corresponding to the core/shell model of colloids and the embedded-atom model (EAM) potential for aluminum proposed in ref. 28. The values of effective parameters favoring dodecagonal order are determined from the system with Dzugutov potential for which the temperature–density domain of DDQC formation is known.12 On adjusting the states of RSS and OPP fluids to obtain the same values of effective parameters, we observe self-assembly of the same DDQC phases at cooling. The values of the parameters for the EAM model28 of liquid aluminum near the liquid–DDQC transition reported in ref. 29 are also the same. This result suggests the common nature of both metallic and soft matter DDQCs arising from competition between length scales.
The proposed method allows reduction of the time spent for searching the parameters favoring a certain solid structure for two-length scale systems. Given the value of effective parameters favoring the formation of some structure, we can predict if this structure self-assembles from the fluid at cooling (see the general scheme in Fig. 1 for an explanation).
Udz = U1(r) + U2(r), | (1) |
(2) |
(3) |
The second potential we use is the repulsive shoulder system (RSS) potential:27
Urss(r) = ε(d/r)n + εnf[2k0(r − σ)], | (4) |
The hard-core analog of the RSS potential was developed by Adler and Yong to explain melting curve extrema.33 Later it was used by Stishov to discuss the possibility of liquid–liquid phase transitions.34 The smooth form presented by eqn (4) was later used to describe a broad range of phenomena, particularly glassy dynamics and formation of decagonal QCs (see ref. 14 and 27 and references therein). Here we take n = 14, k0 = 10, and σ = 1.75 to produce the same values of effective parameters as those for the Dzugutov potential (see Section 3).
The third potential used is the modified oscillating pair potential (OPPm):
Uopp(r) = 1/r15 + aexp(−(r/b)m)cos(kr − φ) | (5) |
Finally, the fourth potential used is the EAM potential proposed in ref. 28 for aluminum. Within the frameworks of EAM, the potential energy of the system Epot is represented as the sum of the pair interaction contribution Upair and the embedding energy F(ρ) depending on the local electron density ρ. So an EAM potential is effectively a many-body one. To compare visually the aluminum EAM potential with two-length-scale pair potentials described above, we use the effective pair format for EAM Ueff = Upair(r) + F(ρ(r)) taking into account the distance dependence of electron density.28 Note that simulations were performed with the original many-body formulation of aluminum EAM.
For molecular dynamics simulations, we use the LAMMPS package.35,36 The system of N = 20000 particles was simulated under periodic boundary conditions using the Nose–Hoover NVT ensemble. This amount of particles is enough to obtain satisfactory diffraction patterns to study (quasi)crystal symmetry (see Fig. 4). Larger systems require too much calculation time necessary for QC equilibration. The molecular dynamics time step was δt = 0.003–0.01 depending on the system temperature.37,38
To study solid phases, we cooled the system starting from a fluid in a stepwise manner and completely equilibrated at each step. The time dependencies of temperature, pressure and configurational energy were analyzed to control equilibration.37
To study the structure of both fluid and solid phases we use radial distribution functions g(r), bond order parameters ql,39–41 diffraction analysis and visual analysis of the snapshots. A detailed description of these methods as well as the procedure for preparing and relaxing the solid phases is presented in ref. 14.
The effective parameters are well defined at 1.2 < λ < 1.6. In this case g(r) subpeaks corresponding to short- and long-bonded particles are perfectly separated at arbitrary ϕ values. For example, in Fig. 3a, we show g(r) of RSS fluids for the effective parameters λ = 1.37 and ϕ = 0.474 corresponding to decagonal QCs.14 But the situation is more complicated in the case of DDQCs considered here. Indeed, in Fig. 3b, we show g(r) for the system with Dzugutov potential with the parameters corresponding to fluids slightly above the fluid–DDQC transition.12 As seen from the figure, the value of λ is about 1.7, which means the first coordination shell of long-bonded particles overlaps with the second coordination shell of short-bonded particles (see splitting of the second g(r) peak in Fig. 3b). To determine effective parameters in this case, we use the method of peak separation widely used in spectroscopy.42,43 The method is based on using high order (2 and 4th) derivatives to separate overlapped peaks. In Fig. 3b we show the second derivative of g(r) for the Dzugutov potential. As seen from the figure, the maximum of d2g(r)/dr2 allows estimation of the distance rs corresponding to intersection of the subpeaks. So the effective numbers of short- and long-bonded particles can be estimated as and .
Fig. 3b can cause a feeling that the separation of the second g(r) peak into two subpeaks is an artificial, non-physical and mathematically fragile procedure. We describe below a method that makes one sure that the second g(r) peak indeed has two subparts. We consider g(r) as the reduction of some complex-valued function of complex variable r on real axes. Using the analytical continuation procedure, we can reconstruct this complex-valued function using g(r) as the “source”. We perform analytical continuation of g(r) into the complex (Rer, Imr) plain numerically using Pade-approximants (see ref. 44–46 for details of the procedure). In Fig. 3c, we show the three dimensional graph with the density-plot projection of analytical continuation of g(r) presented in Fig. 3b where g(r) at real r is shown as a green “tube”. We see that, in the complex-r plain, the second g(r) peak transforms into two peaks. Detailed investigation shows that these peaks correspond to poles where the complex-valued function diverges as 1/(r − Ri), where Ri, i = 1, 2, are the complex coordinates of the poles. If we return to real-number “physical” coordinates, then 1/(r − Ri)-contribution transforms into a Lorentzian peak: 1/[(r − ReRi)2 + (ImRi)2], so ReRi and ImRi become the center and the width of the Lorentzian, respectively.46 The two poles as shown in Fig. 3c transform into a superposition of two Lorentzian peaks at real values of r. Thus the second peak of g(r) in Fig. 3b is indeed the superposition of two distinct subpeaks.
Note that, for two-length systems with finite pair potentials having well defined Fourier transform, the length scale ratio can be extracted from the positions of potential minima in the Fourier space.47
The very fact that the pair correlation function reduces in the main approximation to the superposition of Lorentzian functions rather than, for example, the Gaussian functions is quite an interesting observation that can be generalized to many other systems.48,49 This issue and technical details will be described in a separate paper.
Using the methods described above, the effective parameters of the system with Dzugutov potential slightly above the fluid–DDQC transition have been estimated to be λ = 1.74 and ϕ = 0.42. These values will be further used as reference values to obtain DDQCs in other two-length-scale systems under consideration.
The systems with the parameters, chosen as described above, were cooled down from the fluid phase till the fluid–solid transition occurs. The resulting solid state in all cases consisted of a few highly ordered quasicrystalline grains with pronounced dodecagonal symmetry. Typical snapshots of such DDQC grains in the plain orthogonal to the 12-fold axis are presented in Fig. 4. We see that all the systems demonstrate the same dodecagonal structure. The diffraction patterns of each structure are also shown to demonstrate the identical 12-fold symmetry of the samples.
Hereafter we use the term DDQC having in mind that the system may also fall into a crystalline approximant with local QC symmetry. Moreover, any QC-like configuration constrained by periodic boundary conditions is in fact a periodic approximant in the sense of global order. It should also be noted that the DDQC phase observed may not be the thermodynamically stable one for the systems under consideration. For example, it is known that, for the Dzugutov potential system, the DDQC phase is thermodynamically metastable with respect to the σ-phase periodic approximant.32 The study of thermodynamic stability of DDQC phases observed as well as the investigation of subtle structural features like difference between true QC and approximant phases is out of the framework of this paper. Anyway, the observed DDQC structures are physically stable over the time scale available for simulation and have QC structure on the mesoscale of the simulation box; it is enough for the purposes of this work.
Thus the values of effective parameters favorable to DDQC formation are estimated to be λ ∼ 1.7 and ϕ ∼ 0.4. These values are only estimates; in fact the DDQC phase can form in certain intervals of effective parameters around estimates of this type.14 The exact determination of these intervals for each system under consideration is a matter of separate work. It should be also noted that the obtained values of effective parameters can be only used to predict the formation of certain type of DDQCs presented in Fig. 4. Other types of one-component 3D DDQCs recently observed in computer simulations15,16 have different structure and so different values of effective parameters. The same holds true for the recently reported three dimensional decagonal QCs whose structure differs from that of decagonal QCs obtained in ref. 14. The applicability of the effective parameters method for these new types of one-component QCs is a matter of separate work.
Note that the DDQC structure obtained for aluminum is of much worse quality than that for other systems investigated. Indeed, in Fig. 4d we see a lot of structure defects disturbing the QC structure. This is because we did not tune either the parameters of the EAM potential or the thermodynamic state for aluminum. The system with original parameters proposed in ref. 28 and 29 generates a fluid whose ϕ value is slightly less than the optimal one obtained from the Dzugutov potential (see also Fig. 6).
As reported in ref. 12, the DDQC structure of the system with Dzugutov potential is locally icosahedral. We have checked that the structures of other systems studied are the same as those obtained by Dzugutov. For example, in Fig. 5a we show the typical fragment of DDQC tiling for RSS; the particles which are the centers of icosahedra are colored red and one such icosahedron is marked by interparticle bonds. Such icosahedron is the screen plain projection of the spatial tube structure made of edge-shared icosahedra12 (see Fig. 5b). So the red particles in Fig. 5a represent the axes of these tubes. As seen from Fig. 4 and 5a, there are two joining mechanisms of dodecagonal rings: triple and quadruple junctions. In Fig. 5c we show the local structure of the triple one made of three face-shared icosahedra.
Fig. 6 The potentials of mean force Upmf(r) = −kTln(g(r)) for the systems under investigation. The radial distribution functions g(r) used to calculate Upmf are the same as those presented in Fig. 2b. |
Note that earlier we reported the decagonal QC formation for RSS as well as for other two-length scale systems.14 The building block of such quasicrystals is an icosahedral tube similar to that shown in Fig. 5b but made of face-shared icosahedra.
This journal is © The Royal Society of Chemistry 2017 |