Quasiparticle approach to diffusional atomic scale self-assembly of complex structures: from disorder to complex crystals and double-helix polymers

Abstract

A self-organisation is an universal phenomenon in nature and, in particular, is highly important in materials systems. Our goal was to develop a new theory that provides a computationally effective approach to this problem. In this paper a quasiparticle theory of a diffusional self-organisation of atoms in continuum space during the diffusional time scale has been introduced. This became possible due to two novelties, a concept of quasiparticles, fratons, used for a description of dynamic degrees of freedom and model Hamiltonian taking into account a directionality, length and strength of interatomic bonds. To illustrate a predictive power and achievable level of complexity of self-assembled structures, the challenging cases of self-assembling of the diamond, zinc-blende, helix and double-helix structures, from a random atomic distribution, have been successfully modelled. This approach opens a way to model a self-assembling of complex atomic and molecular structures in the atomic scale during diffusional time.

Introduction

In the last 25 years, significant progress has been made by using molecular dynamic (MD) and Monte Carlo (MC) modelling for study of evolution of multi-atomic systems.^1–4 These methods supplement each other. The MD, is a straightforward approach of to a numerical solution of equations of motions for 6N dynamic variables where N is the number of atoms. However, a sheer number of variable limits the size of studied systems to N~10⁵ and time of its evolution to t~10⁻⁸ s. Recent modifications of MD by coarse graining improve the addressable time but at expense of the lost spatial resolution^5,6. However, the MD is still not well-suited to study slow evolving systems with the typical diffusion time scale.

The MC alternative complements the MD since it can be applicable to the diffusional time scale. The MC approximates the mechanics of atomic motion by a stochastic dynamics of the Markov chain evolution.^3,4,7 A stochastic sampling in the MC dynamics requires a generation of a Markov chain that takes a significant time because it requires a search and update of databases, time scale separation and one process-at-a time execution.

Therefore, there are still significant difficulties of atomic scale prototyping of a slow diffusional self-organisation of atoms in complex structures. This is especially the case if the evolution is in the continuum space, the system consists of comparatively large number of atoms, and evolution time is long, ranging from a fraction of seconds to years. In this paper, we propose such an approach that may be supplemental to MD and MC that addresses their aforementioned limitations of MD and MC in computationally very effective way.

This development turned out to be possible because we (i) introduced a characterisation of a multi-atomic system in terms of quasiparticles named fratons, (ii) proposed a new simple form of phenomenological model potentials describing a directionality, length and strength of atom–atom bonding and (iii) used the kinetic equations of the atomic density field (ADF) theory describing the atomic scale diffusion. The latter was obtained by extending the ADF theory that was first developed for the Ising lattice gas model^8,9 and for the Ising lattice sites diffusion.¹⁰ Recently, the extension of the ADF theory to the relevant case of atomic diffusion in the continuum space was also done.¹¹

To check how this theory works, we chose examples of diffusional self-assembling of several atomic systems. Criteria for these choices are the following:

• a) To estimate a predictive power of the theory, we should choose the initial atomic configuration that has no any resemblance to the expected self-assembled atomic structure. This condition was satisfied by a choice of an initial state described by a completely random distribution of fratons in continuum space. The ‘condensation’ of fratons into the desirable final atomic structure should be provided solely by the input parameters of the model potentials, which characterise the directional atomic bonding.

• b) The model potentials should have sufficient flexibility to be fitted not only to provide a desirable crystallography of a system, but also to approximate its thermodynamic and mechanic properties.

• c) To make an illustration of potentiality of the approach sufficiently convincing, we have to demonstrate that it allows one to successfully simulate a self-assembly of high-complexity structures.

We did successfully model a self-assembly of a randomly distributed atoms into complex atomic configurations of increasing complexity. The structures are ranged from single-component and multi-component crystals to single- and double-helix polymers. As far as we know, a self-assembling of the most complex of them that would started from the completely random state, has been never modelled before.

Model

A central part of the proposed theory is an introduction of non-traditional dynamic variables. Unlike the conventional approach describing the configuration of a classic multi-atomic system by coordinates of atomic centres, the proposed theory describes the atomic configurations by occupation numbers, c(r), of quasiparticles that we call fratons where c(r) is a stochastic function describing two possible events: c(r) is equal to 1 if the point, r, reside inside any atomic sphere and is equal to zero otherwise. In other words, these two events indicate two possible states of each point of the continuous space, r. In this description, the dynamics of the system is described by a creation or annihilation of a fraton at each point of the continuous system: a creation of a fraton at a point, r, indicates that atomic movements resulted in a situation wherein the point r, which was previously outside of any atomic spheres, turned out inside of one of them; annihilation of a fraton describes an opposite process wherein the point, r, which initially is within of an atomic sphere, becomes outside of it.

We also introduced an analogue to the Pauli exclusion principle assuming that two fratons cannot occupy the same point. This exclusion automatically forbids interpenetration of the atoms and thus provides a dynamic ‘exchange’ repulsion preventing the atomic overlap. The introduction of fratons, in many respect, is conceptually similar to a transition to the secondary quantisation for multi-particle Fermi systems.¹²

A m-component system characterised in terms of fratons is described by m stochastic numbers c_α(r) at each site r where α=1,2,..m labels the fratons related to the corresponding atom of the kind α. The averaging over the time-dependent Gibbs ensemble gives the occupation probability, ρ_α(r, t)=<c_α(r, t)>T≤1, where the symbol <…> implies averaging over the ensemble at temperature, T, and time, t. In this definition, the function ρ_α(r, t) is an occupation probability, that point, r, is at any point inside of atomic sphere of any atom of the kind α at the time t.

The temporal evolution of the density function of fratons of the multi-component system is described by the atomic scale kinetic equation of the ADF theory¹⁰ extended to the continuum space:

$$\begin{array}{}\text{(1)}& \frac{\mathrm{d}{\rho }_{\alpha }(\mathbf{r},t)}{\mathrm{d}t}=\sum _{\mathbf{r}\prime }\sum _{\beta =1}^{\beta =m}\frac{L_{\alpha \beta }(\mathbf{r}-\mathbf{r}\prime )}{k_{\mathrm{B}}T}\frac{\delta G}{\delta {\rho }_{\beta }(\mathbf{r}\prime ,t)}\end{array}$$

where indices, α and β, label fratons describing different kinds of atoms (α=1, 2, ..., m), k_B is the Boltzmann constant, G is the non-equilibrium Gibbs free energy functional, L_αβ(r) is a kinetic coefficients matrix. Summation is carried out over all points, r′ of the computational grid approximating the continuum space. The kinetic parameters employed in our microscopic model are related to the phenomenological diffusion coefficients in the continuum model. It was shown previously⁹ that in the continuous model (where k→0) the kinetic coefficients can be expressed as: L(k)=−M_ijk_ik_j where M is a diffusional mobility. The condition ${\sum }_{\mathbf{r}}{L}_{\alpha \beta }(\mathbf{r})=0$ guarantees the conservation of the total number of fratons of each kind during evolution (and thus the conservation of the total volume of the corresponding atoms).

The simplest Gibbs free energy functional entering equation (1) is:

$$\begin{array}{}\text{(2)}& G=\frac{1}{2}\sum _{\mathbf{r},\mathbf{r}\prime }\sum _{\alpha =1}^{\alpha =m}\sum _{\beta =1}^{\beta =m}{w}_{\alpha \beta }(\mathbf{r}-\mathbf{r}\prime ){\rho }_{\alpha }(\mathbf{r}){\rho }_{\beta }(\mathbf{r}\prime )+k_{\mathrm{B}}T\sum _{\mathit{r}}\left[\sum _{\alpha =1}^{\alpha =m}{\rho }_{\alpha }(\mathbf{r})\ln {\rho }_{\alpha }(\mathbf{r})+\left(1-\sum _{\alpha =1}^{\alpha =m}{\rho }_{\alpha }(\mathbf{r})\right)\ln \left(1-\sum _{\alpha =1}^{\alpha =m}{\rho }_{\alpha }(\mathbf{r})\right)\right]-\sum _{\mathit{r}}\sum _{\alpha =1}^{\alpha =m}{\mu }_{\alpha }{\rho }_{\alpha }(\mathbf{r})\end{array}$$

where w_αβ(r−r′) is the model potential of interaction of a pair of fratons of the components α and β, respectively, separated by a distance, r−r′, μ_α is the chemical potential of fratons of the kind α. Summation over r and r′ in equation (2) is carried out over all N₀ sites of the computational grid lattice introduced to discretise the continuous space. The free energy (see equation (2)) corresponds to the mean field approximation⁹. It is asymptotically accurate at low and high temperatures, and its accuracy asymptotically increases if the interaction radius is much greater than the distance between interacting particles¹³. The latter condition is automatically satisfied in our case because interacting particles are fratons, and the computational grid increment, which is the minimum permitted distance between fratons, is much smaller than the atomic radius. The latter is also a requirement of accuracy of a description of a continuous atomic movement. The computational grid can be also interpreted as an Ising lattice and a spacing of the grid as a crystal lattice parameter of this ‘Ising lattice’.

Equation (2) uses the Connolly–Williams approximation¹⁴ that maps the fraton–fraton interaction into the chosen model potential, w_αβ(r−r′). The Fourier transform (FT) of such a potential is:

$$\begin{array}{}\text{(3)}& {\tilde{w}}_{\alpha \beta }(\mathbf{k})=\frac{1}{N_0}\sum _{\mathbf{r}}{w}_{\alpha \beta }(\mathbf{r})e^{-i\mathbf{k}\mathbf{r}}\end{array}$$

where the summation is carried out over all sites of the computational grid, and the wave vector, k, is defined at all quasi-continuum points, k, of the first Brillouin zone of the computational grid, that is, at all N₀ the points in the k-space permitted by the periodical boundary conditions.

Structures that are really complex are usually formed in systems with directional covalent bonds between atoms. Therefore a theory whose goal is to study such systems by using the phenomenological model Hamiltonian should formulate the atom–atom interaction potentials, k, as directional functions of atomic separation distance, r − r'. This function should describe the directionality of these bonds, their length and strength. The parameters of the potentials can be fitted to those calculated in quantum chemistry, Density Function formalism and/or by fitting to the observed crystallographic, mechanic and thermodynamic properties of the system. To make modelling sufficiently efficient, the approximation of the potentials should be as simple as possible without sacrificing a predictive power of the model.

In this paper, such model potentials are proposed. We tested their validity for the most challenging cases of self-assembling, some of which are not being modelled before. A central idea of the formulation of these potentials is based on the explicit use of what we call a bonding star, a group of vectors with the common beginning that describes directions and lengths of interatomic bonds.

A structural cluster designated as αβ, is determined by a star of vectors numbered by the index j_αβ. The vectors, j_αβ, have the same origin, parallel to the corresponding bonds, j_αβ, between the atoms of the kinds α and β, and have the lengths of these bonds. Therefore, the cluster, αβ, is characterised by a set of geometrical points determined by the ends of the vectors of the star and its centre. The geometry of the stars, αβ, and the strength of the bonds characterising the star describe interaction in the system and should be the only factor determining the structure of a self-assembled multi-atomic system.

These definitions allow us to introduce a cluster function, ${\Psi }_{\alpha \beta }^{\mathrm{clstr}}(\mathit{r})$ whose FT is:

$$\begin{array}{}\text{(4)}& {\tilde{\Psi }}_{\alpha \beta }^{\mathrm{clstr}}(\mathbf{k})=\sum _{j_{\alpha \beta }}\omega (j_{\alpha \beta },\mathbf{k})e^{-i\mathbf{k}{\mathbf{r}}_{j_{\alpha \beta }}}\end{array}$$

where summation is carried out over all vectors, ${\mathbf{r}}_{j_{\alpha \beta }}$, of the αβ star of the cluster, ω(j_αβ, k) is a function characterising a strength of the bond j_αβ. The coefficients, ω(j_αβ, k) determine the thermodynamic and mechanical properties of the simulated system and can be used as fitting parameters to reproduce them.

Using these definitions, we present the FT of the model potential as sum of what we call the short-range and long-range interactions:

$$\begin{array}{}\text{(5)}& {\tilde{w}}_{\alpha \beta }(\mathit{k})={\lambda }_1{\tilde{\theta }}_{\alpha }(\mathrm{k}){\delta }_{\alpha \beta }+{\lambda }_2(\mathrm{k}){\Psi }_{\alpha \beta }^{\mathrm{cltr}}(\mathbf{k}){\Psi }_{\alpha \beta }^{\mathrm{cltr}}{(\mathit{k})}^{*}\end{array}$$

The first term in equation (5) describes the spherically symmetrical short-range fraton–fraton pair interaction. The function θ_α(r) is schematically presented in Figure 1a, where r₁ is a length parameter determining atomic radius, Δr is the width of repulsion part, $\xi =\frac{\left|{\theta }_{\max }(\mathbf{r})\right|}{\left|{\theta }_{\min }(\mathbf{r})\right|}$ is the ratio between the modules of minimum and maximum values of the shape function and λ₁ is a constant determining the strength of the short-range atomic repulsion. In particular, the parameter λ₁ characterises the rigidity of the atomic spheres and its value can be determined by a fitting of the elastic properties of a given system. The FT of the function, ${\tilde{\theta }}_{\alpha }(\mathbf{k})$, schematically shown in Figure 1b is:

$$\begin{array}{}\text{(6)}& \begin{array}{l}{\tilde{\theta }}_{\alpha }(\mathbf{k})=\frac{4\pi }{k^3}(-\sin (\mathbf{k}{\mathbf{r}}_1)-\mathbf{k}{\mathbf{r}}_1\cos (\mathbf{k}{\mathbf{r}}_1))+\xi (\sin (\mathbf{k}({\mathbf{r}}_1+\Delta \mathbf{r}))-\mathbf{k}({\mathbf{r}}_1+\Delta \mathbf{r})\cos (\mathbf{k}({\mathbf{r}}_1+\Delta \mathbf{r})\\ -\sin (\mathbf{k}{\mathbf{r}}_1)+\mathbf{k}{\mathbf{r}}_1\cos (\mathbf{k}{\mathbf{r}}_1)))\end{array}\end{array}$$

Figure 1

(a) Schematic representation of short-range potential θ(r). (b) Example of FT of θ(r) with the following input parameters: ξ=4, $\Delta \hat{r}=0.25$.

The second term of equation (5) is the long-range part of the fraton–fraton interaction describing a directional bonding of atoms of the kind α and β.

In fact, the long-range interaction in equation (5) is presented as a bilinear expansion in cluster functions, ${\Psi }_{\alpha \beta }^{\mathrm{clstr}}(\mathbf{k})$, λ₂(k) is a fitting parameter determining a strength of the long-range interaction.

The indexes, α and β can be dropped for a single-component system. Then equation (5) is simplified to:

$$\begin{array}{}\text{(7)}& \tilde{w}(\mathbf{k})={\lambda }_1\tilde{\theta }(\mathbf{k})+{\lambda }_2(\mathbf{k}){\left|{\Psi }^{\mathrm{clstr}}(\mathbf{k})\right|}^2\end{array}$$

In this paper we consider a particular case of application of the fratonic theory allowing to obtain a single crystalline state. This is a case of specifically oriented clusters. The constraint lifts the angular isotropy of the system and thus allows us to prevent the formation of a ‘polycrystalline state’ that is an atomic aggregate of grains with the same atomic structure but different orientation. A self-assembling producing such a ‘polycrystal’ would make it difficult to identify the equilibrium atomic structure. However, in the general case, in which the angular isotropy is not lifted, the potential described by equation (7) describes a growth of polycrystal (see Supplementary Figure 1S in Supplementary Information). In this case the interaction energy of a pair of fratons is independent from its orientation. It can be achieved using a rotational averaging of a cluster function ${\Psi }_{\alpha \beta }^{\mathrm{clstr}}(\mathbf{k})$.

Results

To illustrate the versatility and effectiveness of the fraton theory, we tested its application to the modelling of the self-assembly of three groups of three-dimensional structures of increasing complexity. They are single-component crystals, two-component crystals and a polymer with a double-helix structure mimicking biological macromolecules. The modelling was carried out by numerical solution of the FT representation of the kinetic equation (1).

In our simulations, we used the reduced parameters, and, in particular, average density, defined as ${\hat{\overline{\rho }}}_{\alpha }={\rho }_{\alpha }^{\mathrm{at}}\frac{4\pi R_{\alpha }^3}{3}$ where ${\rho }_{\alpha }^{\mathrm{at}}=\frac{N_{\alpha }}{V}$ is the density of α atoms in the ground state, N_α is number of the atoms of sort α, V is the total volume of the system and R_α is the atomic radius of this atom. According to this definition, the reduced density, ${\hat{\overline{\rho }}}_{\alpha }$, is also a fraction of all computational grid sites occupied by fratons of the kind α. The input parameter ξ of the energy θ_α(r) is measured in units of k_BT_o, where T_o is the solidification temperature. The lengths are measured in units of r₁, which is very close to the atomic radius; the grid lattice increment, $\hat{l}$ (the spacing of the underlying Ising lattice), is defined as a fraction of the atomic radius. The temperature $\hat{T}$ is also measured in units of T_o. The reduced time, $\hat{t}$, is measured in units of typical atomic migration time, τ_o. The reduced kinetic coefficients, $\hat{L}(\mathbf{r})$, are measured in units of ${\tau }_0^{-1}$ and $\tilde{L}(\mathbf{k})=D{k}^2$, where D is a constant. The numerical solution of the reduced form of kinetic equation (1) was obtained by using the semi-implicit Fourier-spectral method¹⁵. The number of computational grid sites in the simulation box was chosen in compromise between a necessity to reasonably reproduce the atomic structure and to optimise the simulation time. A time step of Δt=0.01 was used. The periodic boundary condition was imposed for all simulations. The first example that we will consider in this article is the self-assembly of a single-component crystal with several atoms in a Bravais lattice unit cell. As an example, we tested whether the model Hamiltonian describing tetrahedral orientation of the interatomic bonds results in a self-assembling producing the diamond structure. An image of the structural cluster in this case is convenient to choose as a star of eight bonding vectors whose ends are located at the sites of a cubic unit cell of the diamond lattice: (000), $\left(0\frac{1}{2}\frac{1}{2}\right)$, $\left(\frac{1}{2}0\frac{1}{2}\right)$, $\left(0\frac{1}{2}\frac{1}{2}\right)$, $\left(\frac{1}{4}\frac{1}{4}\frac{1}{4}\right)$, $\left(\frac{3}{4}\frac{1}{4}\frac{3}{4}\right)$, $\left(\frac{1}{4}\frac{3}{4}\frac{3}{4}\right)$ and $\left(\frac{3}{4}\frac{3}{4}\frac{1}{4}\right)$. The function Ψ^cltr(k) was constructed by using its definition (equation (4)) and the coordinates of the chosen structural cluster points:

$$\begin{array}{}\text{(8)}& {\Psi }^{\mathrm{cltr}}(\mathbf{k})=\left(1+e^{-i\frac{a}{4}\left(k_x+k_y+k_z\right)}\right)\left(1+e^{-i\frac{a}{2}\left(k_x+k_y\right)}+e^{-i\frac{a}{2}\left(k_y+k_z\right)}+e^{-i\frac{a}{2}\left(k_x+k_z\right)}\right)\end{array}$$

where a is a lattice constant of diamond structure, $k_i=\frac{2\pi m_i}{aN}$ (where i=x,y or z, m_i=1…N, N is the number of simulation grid in a given direction). With this definition, the FT of the long-range interaction, ${\tilde{w}}_{\mathrm{LR}}(\mathbf{k})$, can be written as:

$$\begin{array}{}\text{(9)}& {\tilde{w}}_{\mathrm{LR}}(\mathbf{k})={\lambda }_2(\mathbf{k}){\tilde{\Omega }}_{\mathrm{D}}(\mathbf{k})\end{array}$$

where

$$\begin{array}{}\text{(10)}& {\tilde{\Omega }}_{\mathrm{D}}\left(\mathbf{k}\right)={\Psi }^{\mathrm{cltr}}(\mathbf{k}){\left({\Psi }^{\mathrm{cltr}}(\mathbf{k})\right)}^{*}=\left(2+2\cos \left(\frac{a}{4}\left(k_x+k_y+k_z\right)\right)\right)\times \left(4+4\left(\cos \left(\frac{k_{x}a}{2}\right)\cos \left(\frac{k_{y}a}{2}\right)+\cos \left(\frac{k_{y}a}{2}\right)\cos \left(\frac{k_{z}a}{2}\right)+\cos \left(\frac{k_{x}a}{2}\right)\cos \left(\frac{k_{z}a}{2}\right)\right)\right)\end{array}$$

To choose for the model potential the form of λ(k), we used the following consideration. The function |Ψ(k)|² by definition (see equation (10)) is periodical in k-space. Its longest period is along the [111] direction and equal to $k_1=\frac{8\pi }{a}$. We assume that the function, λ₂(k), is a step function. It is defined as:

$$\begin{array}{}\text{(11)}& {\lambda }_2\left(\mathbf{k}\right)=\{\begin{array}{c}{\lambda }_2\mathrm{if}0\le \left(k_x+k_y+k_z\right)\le \frac{8\pi }{a}+\delta \\ 0\mathrm{otherwise}\end{array}\end{array}$$

where δ is a positive constant.

Then the FT of the model potential can be written as:

$$\begin{array}{}\text{(12)}& \tilde{w}(\mathbf{k})={\lambda }_1\tilde{\theta }(\mathbf{k})+{\lambda }_2(\mathbf{k}){\tilde{\Omega }}_{\mathrm{D}}(\mathbf{k})\end{array}$$

where the first term describes the short-range interaction. The energy parameters λ₁ and λ₂ were normalised by the absolute value of a difference between the maximum and minima values of the function ${\tilde{\Omega }}_{\mathrm{D}}(\mathbf{k})$ and $\tilde{\theta }(\mathbf{k})$, respectively.

In this simulation, the initial configuration was an embryo consisting of the small variation of the fratons density at the sites of the structural cluster of diamond structure embedded in the gas of disordered fratons (Figure 2a). However, the initial state can be also chosen as local random ‘infinitesimals’ statics fluctuations with respect to homogeneous state. The spontaneous self-organisation of fratons into the diamond structure is shown in Figure 2. The intermediate step in the pattern formation dynamics at the reduced time t=60,000 is shown in Figure 2b. A very interesting aspect of this growth is the development of the transient bcc structure (Figure 2c) in the early stages of growth. Its lattice parameter approximately is half that of the diamond structure.

Figure 2

Example of a self-assembly of fratons into diamond structure at reduced times $\hat{t}$ of (a) $\hat{t}=0$, (b) $\hat{t}=60,000$, (c) $\hat{t}=100,000$, (d) $\hat{t}=280,000$ and (e) $\hat{t}=300,000$. The parameters in these simulations are ${\hat{\lambda }}_1=14.085$, ${\hat{\lambda }}_2=-7.042$, $\hat{a}=4.57$, ξ=2, $\hat{D}=1$, $\hat{\overline{\rho }}=0.07$, $\hat{l}=0.286$, $\Delta \hat{r}=0.17$ and $\hat{T}=0.732$. The initial configuration is the atomic cluster of the diamond structure placed in the centre of the simulation box. The size of the simulation box is 64×64×64.

A probability of occupation of each sites of the obtained transient bcc lattice is less than unity. The latter indicates that the bcc lattice sites, in fact, are randomly occupied by atoms and their vacancies and form the disordered distribution in the bcc lattice. This result is fully consistent with the thermodynamics dictating a mandatory presence of a disordered distribution of vacancies in the bcc (or any other) lattice at finite temperature, and thus dictating the greater number of lattice sites than the number of occupying them atoms. It should be also noted that the total number of atoms in our simulation is automatically conserved because the kinetic equation (or more specifically, a choice of matrix of kinetic coefficients) guarantees the atomic conservation.

The formation of the diamond lattice, in fact, is a fratonic theory description of the bcc→B32 (NaTl-type) ordering of atoms and their vacancies over the preferential sites of the previously ‘condensed’ bcc host lattice. A result of such an ordering is a doubling of the crystal lattice parameters of the underlying bcc lattice. It should be also mentioned that a possibility of the crystallisation of the diamond structure through the transient bcc structure is a result that could hardly be predicted in advance for a system with tetrahedral atomic bonding.

To better visualise the final structure, presented in Figure 2e, the links between the first neighbours are shown. The growth dynamic of the diamond structure is shown in Supplementary Movie 1 (see Supplementary Information). To better characterise the different steps of the growth dynamic, the diffraction patterns at different time steps have been calculated. The diffraction pattern, which is a distribution of intensity of scattered radiation in the three-dimensional reciprocal space of the wave vectors, k=(k_x, k_y, k_z), was determined as the squared modulus of the FT of the density function. The strongest diffraction peaks that characterise the diamond structure are: (220), (111), (311) and (400). The peaks {200}, which are forbidden for the diamond lattice, are present in the diffraction pattern of the bcc structure. Then, following these peaks it is possible to distinguish the different stage of the growth dynamic. The evolution of the diffraction pattern during the growth dynamic of the diamond structure is shown in Supplementary Movie 2 (see Supplementary Information).

Continuing a gradual increase in the complexity of the modelled structure, we also considered the formation of a two-component crystalline phase in a system with tetrahedral direction of bonds expecting the formation of the zinc-blende structure. A two-component systems atomic arrangement is formed by a ‘condensation’ of two kinds of fratons, belonging to type A and B. This condensation should produce atoms A and B, respectively. The spontaneous arrangement caused by an equilibration of a disordered distribution of one sort of fratons is described by the kinetic equation (1). For a two-component system, equation (1) is reduced to two equations for two fraton densities ρ_A(r) and ρ_B(r). Thus in the reciprocal space these equations can be written as:

$$\begin{array}{}\text{(13a)}& \frac{\partial {\tilde{\rho }}_A(\mathbf{k},t)}{\partial t}=L_{AA}(\mathbf{k})\left({\tilde{w}}_{AA}(\mathbf{k}){\tilde{\rho }}_A(\mathbf{k},t)+{\tilde{w}}_{AB}(\mathbf{k}){\tilde{\rho }}_B(\mathbf{k},t)+{\left(\ln \frac{{\rho }_A(\mathbf{r}\mathbf{\prime },t)}{1-{\rho }_A(\mathbf{r}\mathbf{\prime },t)-{\rho }_B(\mathbf{r}\mathbf{\prime },t)}\right)}_{\mathbf{k}}\right)+L_{AB}(\mathbf{k})\left({\tilde{w}}_{BB}(\mathbf{k}){\tilde{\rho }}_B(\mathbf{k},t)+{\tilde{w}}_{AB}(\mathbf{k}){\tilde{\rho }}_A(\mathbf{k},t)+{\left(\ln \frac{{\rho }_B(\mathbf{r}\mathbf{\prime },t)}{1-{\rho }_A(\mathbf{r}\mathbf{\prime },t)-{\rho }_B(\mathbf{r}\mathbf{\prime },t)}\right)}_{\mathbf{k}}\right)\end{array}$$

$$\begin{array}{}\text{(13b)}& \frac{\partial {\tilde{\rho }}_B(\mathbf{k},t)}{\partial t}=L_{BB}(\mathbf{k})\left({\tilde{w}}_{BB}(\mathbf{k}){\tilde{\rho }}_B(\mathbf{k},t)+{\tilde{w}}_{AB}(\mathbf{k}){\tilde{\rho }}_A(\mathbf{k},t)+{\left(\ln \frac{{\rho }_B(\mathbf{r}\mathbf{\prime },t)}{1-{\rho }_A(\mathbf{r}\mathbf{\prime },t)-{\rho }_B(\mathbf{r}\mathbf{\prime },t)}\right)}_{\mathbf{k}}\right)+L_{AB}(\mathbf{k})\left({\tilde{w}}_{AA}(\mathbf{k}){\tilde{\rho }}_A(\mathbf{k},t)+{\tilde{w}}_{AB}(\mathbf{k}){\tilde{\rho }}_B(\mathbf{k},t)+{\left(\ln \frac{{\rho }_A(\mathbf{r}\mathbf{\prime },t)}{1-{\rho }_A(\mathbf{r}\mathbf{\prime },t)-{\rho }_B(\mathbf{r}\mathbf{\prime },t)}\right)}_{\mathbf{k}}\right)\end{array}$$

where A and B designate two sorts of fratons.

The FTs of the interaction energies, ${\tilde{w}}_{\alpha \beta }(\mathbf{k})$, determined by equation (5), are:

$$\begin{array}{}\text{(14a)}& {\tilde{w}}_{AA}(\mathbf{k})={\lambda }_{1A}{\tilde{\theta }}_A(\mathbf{k})+{\lambda }_{2A}(\mathbf{k}){\tilde{\Omega }}_{ZB}^{AA}(\mathbf{k})\end{array}$$

$$\begin{array}{}\text{(14b)}& {\tilde{w}}_{BB}(\mathbf{k})={\lambda }_{1B}{\tilde{\theta }}_B(\mathbf{k})+{\lambda }_{2B}(\mathbf{k}){\tilde{\Omega }}_{ZB}^{BB}(\mathbf{k})\end{array}$$

$$\begin{array}{}\text{(14c)}& {\tilde{w}}_{AB}(\mathbf{k})={\lambda }_{2AB}(\mathbf{k}){\tilde{\Omega }}_{ZB}^{AB}(\mathbf{k})\end{array}$$

where ${\tilde{\Omega }}_{ZB}^{\alpha \beta }(\mathbf{k})={\Psi }_{\alpha }(\mathbf{k}){{\Psi }^{*}}_{\beta }(\mathbf{k})$

The zinc-blende structure has two atoms, A and B in a primitive unit cell of the fcc Bravais lattice with positions (000) and $\left(\frac{1}{4}\frac{1}{4}\frac{1}{4}\right)$, correspondingly. To describe the model Hamiltonian providing evolution to this structure, we needed two structural clusters, which is, the clusters of type A and B. The cluster A consists of four points: the point (000) and its nearest neighbours in the fcc lattice, $a\left(\frac{1}{2}\frac{1}{2}0\right)$, $a\left(0\frac{1}{2}\frac{1}{2}\right)$, $a\left(\frac{1}{2}0\frac{1}{2}\right)$. The cluster B also consists of four point. They are obtained from the four points of the cluster A by the shift, $a\left[\frac{1}{4}\frac{1}{4}\frac{1}{4}\right]$. With this definition, the Ψ-functions for the two structural clusters are:

$$\begin{array}{}\text{(15a)}& {\Psi }_A^{\mathrm{cltr}}(\mathbf{k})=\left(1+e^{-i\frac{a}{2}\left(k_x+k_y\right)}+e^{-i\frac{a}{2}\left(k_y+k_z\right)}+e^{-i\frac{a}{2}\left(k_x+k_z\right)}\right)\end{array}$$

$$\begin{array}{}\text{(15b)}& {\Psi }_B^{\mathrm{cltr}}(\mathbf{k})=e^{-i\frac{a}{4}\left(k_x+k_y+k_z\right)}\left(1+e^{-i\frac{a}{2}\left(k_x+k_y\right)}+e^{-i\frac{a}{2}\left(k_y+k_z\right)}+e^{-i\frac{a}{2}\left(k_x+k_z\right)}\right)\end{array}$$

Using this definition in equation (7) and assuming that the functions λ_2A(k)=λ_2B(k)=λ_2AB(k)=λ₂(k), where λ₂(k) is defined by equation (11), we have:

$$\begin{array}{}\text{(16a)}& {\tilde{\Omega }}_{ZB}^{AA}\left(\mathbf{k}\right)={\tilde{\Omega }}_{ZB}^{BB}\left(\mathbf{k}\right)=4+4\left(\cos \left(\frac{k_{x}a}{2}\right)\cos \left(\frac{k_{y}a}{2}\right)+\cos \left(\frac{k_{y}a}{2}\right)\cos \left(\frac{k_{z}a}{2}\right)+\cos \left(\frac{k_{x}a}{2}\right)\cos \left(\frac{k_{z}a}{2}\right)\right)\end{array}$$

$$\begin{array}{}\text{(16b)}& \begin{array}{l}{\tilde{\Omega }}_{ZB}^{AA}\left(\mathbf{k}\right)=2\cos \left(\frac{a}{4}\left(k_x+k_y+k_z\right)\right)\times \left(4+4\right)\left(\cos \left(\frac{k_{x}a}{2}\right)\cos \left(\frac{k_{y}a}{2}\right)+\cos \left(\frac{k_{y}a}{2}\right)\cos \left(\frac{k_{z}a}{2}\right)+\cos \left(\frac{k_{x}a}{2}\right)\cos \left(\frac{k_{z}a}{2}\right)\right)\end{array}\end{array}$$

The difference in size of different species of atoms has been taken into account in the short-range potential. In these simulations, the ratio of two atomic radii was chosen to be 0.875. The size of the simulation grid, $\hat{l}=0.25$, was measured in the unities of r_1A. Therefore, the value of r_1B was chosen equal to $3.5\hat{l}$. The growth of the zinc-blende structure is shown in Figure 3 and in Supplementary Movie 3 (see Supplementary Information).

Figure 3

Example of a self-assembly of fratons into zinc-blende structure at reduced times $\hat{t}$ of (a) $\hat{t}=0$, (b) $\hat{t}=160,000$, (c) $\hat{t}=190,000$, (d) $\hat{t}=280,000$ and (e) $\hat{t}=3,000,000$. The parameters in this simulation are ξ=2, ${\hat{D}}_{AA}={\hat{D}}_{BB}=1$, ${\hat{D}}_{AB}=-0.5$, ${\hat{\lambda }}_{1A}=3.77$, ${\hat{\lambda }}_{2A}=-1.88$, ${\hat{\lambda }}_{1B}=5.84$, ${\hat{\lambda }}_{2B}=-2.92$, ${\hat{\lambda }}_{2AB}=-2.26$, $\hat{l}=0.25$, r_1A=1.143 r_1B, $\Delta \hat{r}=0.17$, $\hat{a}=4.0$, ${\hat{\overline{\rho }}}_A=0.07$, ${\hat{\overline{\rho }}}_B=0.045$ and $\hat{T}=0.235$. The initial configuration is the atomic cluster of a diamond structure placed in the centre of the simulation box. The size of the simulation box is 64×64×64. Two sorts of atoms with different atomic sizes are indicated in red and green.

Molecules with a helix architecture are observed in organic materials¹⁶, helix-shaped graphite nanotubes,^17,18 liquid crystal,^19,20 proteins, and, of course, DNA and RNA polymeric molecules.^21,22 However, the most challenging test of the potency of the fraton theory would be its ability to describe a spontaneous self-assembly for the most interesting case relevant to biology, that is, the self-assembly of a double-helix polymer from a ‘soup’ of randomly distributed monomers. We chosen this system because, as far as we know, a self-assembly of randomly distributed monomers into double-thread helix polymers was a too complex phenomenon to prototype by the existing methods.

To model a spontaneous self-assembling of a helix polymer consisting of two complimentary threads, we considered two types of mutually complementary fratons whose ‘condensation’ should produce two kinds of mutually complementary monomers.

The model fraton–fraton potential producing a helix structure should be directional and have a built-in chirality. This is achieved by a choice of a bonding star of the cluster that has the chirality corresponding to the desired helix geometry. The geometrical parameters of the configuration of the structural cluster are the pitch length, P, the number of coils per pitch, n₀=6, the distance between coils in z-direction, h, and the radius of the coil, u (see Figure 4). Then the coordinates of points of the helix occupied by molecules are:

$$\begin{array}{}\text{(17)}& {\mathbf{r}}_s=\left(u\cos \left(\frac{2\pi }{n_0}s\right),usin\left(\frac{2\pi }{n_0}s\right),\frac{h}{n_0}s\right)\end{array}$$

where s runs from 0 to n₀−1, n₀ is the number of coils in the pitch.

Figure 4

Illustration of the geometrical parameters of the helical structure: P is a pitch, h is a distance between the nearest coils along the z-axis and helix radius u.

We again start from a construction of the model Hamiltonian that should lead to a condensation of a randomly distributed fratons into the helix structure. To do this, the structural cluster that consists of two pitches has been chosen. We had to choose a two-period cluster because fratons in this model have no orientational degrees of freedom and thus have no built-in chirality. The second-pitch segment of the cluster is needed to introduce a chirality into the long-range part of the model potential. The size of the structural cluster would be drastically reduced if the chirality were built-in in the short-range part of the fraction–fraton interaction. This can be done by a straightforward modification of the short-range interaction.

Using this definition of cluster for the formulation of the function Ψ(k) for the helical structure, presented in Figure 4 gives:

$$\begin{array}{}\text{(18)}& {\Psi }^{\mathrm{cltr}}\left(\mathbf{k}\right)=\left(1+e^{-n_{0}h{k}_z}\right)\left(\sum _{n=0}^{n_0-1}{e}^{-i\left(x_n{k}_x+y_n{k}_y+z_n{k}_z\right)}\right)\end{array}$$

Then function ${\tilde{\Omega }}_H\left(\mathbf{k}\right)$ is:

$$\begin{array}{}\text{(19)}& {\tilde{\Omega }}_H\left(\mathbf{k}\right)=\left(2+2\cos \left(h{n}_0{k}_z\right)\right)\left(6+2\sum _{n>m=0}^{n_0-1}\cos \left(\varphi \left(n,m\right)\right)\right)\end{array}$$

where

$$\begin{array}{}\text{(20)}& \varphi \left(n,m\right)=k_{x}u\left(\cos \left(\frac{2\pi n}{n_0}\right)-\cos \left(\frac{2\pi m}{n_0}\right)\right)+k_{y}u\left(\sin \left(\frac{2\pi n}{n_0}\right)-\sin \left(\frac{2\pi m}{n_0}\right)\right)+k_{z}h\left(n-m\right)\end{array}$$

We also assumed that the desired helix has a single monomer in each coil and approximate these monomers by a spherical shape. The latter is not a critical assumption for the theory. It is just a simplification that reduces the computational time. The short-range interaction was described by equation (6). The parameters of the $\tilde{\theta }(\mathbf{k})$ function were chosen the same as for the diamond structure. The size of simulation box was 210×32×32. The initial embryo lifting the spatial and rotational energy degeneration was a one pitch inhomogeneity introduced in the centre of the simulation box. As was discussed previously, to construct a model potential we should choose the function λ₂(k). For the same reasons as before, we assumed that the function, λ₂(k), is a step function defined as:

$$\begin{array}{}\text{(21)}& {\lambda }_2\left(\mathbf{k}\right)=\{\begin{array}{l}{\lambda }_2\mathrm{if}-\delta \le \varphi \left(n,m\right)\le 2\pi +\delta \mathrm{for}n-m=1,\hfill \\ 0\mathrm{otherwise}\hfill \end{array}\end{array}$$

The condition n−m=1 defines the longest period of the function ${\tilde{w}}_{\mathrm{LR}}(\mathbf{k})$ for helix.

In the first step, we considered a random distribution of fratons of one kind. The solution of equation (1) with the model potential given by equations (19) and (20) in this case describes an evolution that eventually produces the single-stranded helix shown in Figure 5.

Figure 5

Example of self-assembly of the fratons into the helix structure at reduced time of (a) $\hat{t}=0$, (b) $\hat{t}=200,000$, (c) $\hat{t}=250,000$, (d) $\hat{t}=300,000$ and (e) $\hat{t}=700,000$. The input parameters in this simulation are: ${\hat{\lambda }}_1=61.14$, ${\hat{\lambda }}_2=-69.87$, $\hat{h}=\hat{u}=1.56$, n₀=6, ξ=2, $\hat{D}=1$, $\hat{\overline{\rho }}=0.0096$, $\hat{l}=0.22$, $\Delta \hat{r}=0.17$ and $\hat{T}=0.568$. The size of the simulation box is 32×32×210. The initial configuration shown in a is n₀+1 coils in the helix structure.

The introduction of the second kind of fratons, which is complementary to the first kind, results in their self-assembly of a complementary strand and the formation of a double-helix molecule. To model double-helix growth, we have to introduce (as for the zinc-blend structure) complimentary fratons of two types, A and B that form each helix thread. We used the same type of the structural clusters for the fratons of the kind A and B. Each of them consists of two pitches of helix structure. However, the clusters are rotated with respect to each other about z-axis by φ=π. Then, using the definition of the functions ${\Psi }_{\alpha }\left(\mathbf{k}\right)$ by equation (4) for α=A, B, gives:

$$\begin{array}{}\text{(22a)}& {\Psi }_A^{\mathrm{cltr}}\left(\mathbf{k}\right)=\left(1+e^{-i\frac{n_{0}h}{2}{k}_z}\right)\left(\sum _{n=0}^{n_0-1}{e}^{-i\left(x_n{k}_x+y_n{k}_y+z_n{k}_z\right)}\right)\end{array}$$

$$\begin{array}{}\text{(22b)}& {\Psi }_B^{\mathrm{cltr}}\left(\mathbf{k}\right)=e^{-i\frac{n_{0}h}{2}{k}_z}\left(1+e^{-i\frac{n_{0}h}{2}{k}_z}\right)\left(\sum _{n=0}^{n_0-1}{e}^{-i\left(x_n{k}_x+y_n{k}_y+z_n{k}_z\right)}\right)\end{array}$$

We assume that, where ${\tilde{\Omega }}_{DH}^{AA}\left(\mathbf{k}\right)={\tilde{\Omega }}_{DH}^{BB}\left(\mathbf{k}\right)={\tilde{\Omega }}_H\left(\mathbf{k}\right)$, ${\tilde{\Omega }}_H\left(\mathbf{k}\right)$ is defined by equation (19). Then the function ${\tilde{\Omega }}_{DH}^{AB}\left(\mathbf{k}\right)$ is:

$$\begin{array}{}\text{(23)}& {\tilde{\Omega }}_{DH}^{AB}\left(\mathbf{k}\right)=2\cos \left(h\frac{n_0}{2}{k}_z\right)\left(2+\mathrm{2cos}\left(h{n}_0{k}_z\right)\right)\left(6+2\sum _{n>m=0}^{n_0-1}\cos \left(\varphi \left(n,m\right)\right)\right)\end{array}$$

The simulation box size and input parameters for two kinds of complimentary fratons were chosen the same as for a single-thread helix and the interaction between helix is defined by ${\hat{\lambda }}_{2AB}=-1.78$, $\hat{T}=0.033$.

In spite of all these oversimplifications, this model describes some generic features relevant to the spontaneous formation of single-stranded polymeric molecule and the growth of the complementary strand of the monomers eventually producing a double-stranded helix configuration (see Figure 6). In this case, the first single-stranded helix is a template for the aggregation on it of complementary monomers to form a double-stranded helix. For clarity, we show the clusters of fratons (monomers) of the second strand in red. The spontaneous growth of the helix and double-helix structure is also shown in Supplementary Movie 4 and 5 (see Supplementary Information).

Figure 6

Self-assembly of the fratons into a double-helix structure at reduced time of (a) $\hat{t}=0$, (b) $\hat{t}=150,000$, (c) $\hat{t}=200,000$, (d) $\hat{t}=300,000$ and (e) $\hat{t}=1,500,000$. The input parameters in this simulation are ${\hat{\lambda }}_{1A}={\hat{\lambda }}_{1B}=4.07$, ${\hat{\lambda }}_{2A}={\hat{\lambda }}_{2B}=-4.07$, ${\hat{\lambda }}_{2AB}=-1.78$, r_1A=r_1B, $\hat{h}=\hat{u}=1.56$, n₀=6, ξ=2, ${\hat{D}}_{AA}={\hat{D}}_{BB}=1$, ${\hat{D}}_{AB}=-0.5$, ${\hat{\overline{\rho }}}_A={\hat{\overline{\rho }}}_B=0.0096$, $\hat{l}=0.22$, $\Delta \hat{r}=0.17$, and $\hat{T}=0.033$. The size of the simulation box is 32×32×210. The initial configuration shown in a is one helix and one coil of the second helix.

Discussion

In this paper, we selected the most difficult cases wherein the initial system is atomically disordered so that its configuration ‘knows’ nothing about the final atomic pattern that should be spontaneously self-assembled. This self-assembling is driven only by the chosen model Hamiltonian, and, specifically, by mutual orientation, length and strength of interatomic bonds. On a top of that, we considered situations wherein the self-assembling is diffusional and usually takes a long time, which may range from a fraction of a second to years. A typical time of this evolution is dictated by the typical time of evolution of time-dependent ensemble rather than typical times of atomic dynamics like time of atomic vibrations. Difficulty in addressing such slow evolving systems probably was a reason why a spontaneous formation of some of them (crystals and polymers) from a liquid solution of atoms or monomers has not been modelled yet in the diffusional time scale.

The developed approach opens a way to answer numerous outstanding questions concerning the atomistic mechanisms of the formation of defects (dislocations, grain boundaries, etc.), nucleation in solid–solid transformations, the formation of polymers due to aggregation of monomers in their solution, folding and crystallisation of polymers, and their responses to external stimuli. This list can be significantly extended.

Especially interesting are the modelling results describing the spontaneous self-assembly of monomers into a single-stranded polymeric helix and the formation of a double-helix structure obtained by aggregation of complementary monomers on the single-stranded helix playing the role of a template. This result may be also considered as an attempt to formulate and execute the simplest prototyping of the spontaneous formation of homopolymeric DNA from a liquid solution of monomers playing the role of nucleotides.

Finally, the use of the new model Hamiltonian formulated in terms of the structural clusters and proposed fraton model provide already a ready tool to address a general problem of spontaneous pattern formation by self-assembling of any randomly distributed building elements in the time scale ranging from sub-seconds to years. This approach can be also straightforwardly extended for the prototyping of self-assembly of elementary building block monomers with more complex molecular structures. In the latter case, we have to generalise the concept of fratons of atoms by introducing fraton of molecules and modify accordingly the model short-range part of the model fraton–fraton Hamiltonian. Then this modification should provide a ‘condensation’ of the molecular fratons into molecules and subsequent self-assembly of these molecules. In principle, this approach can be even used for the description of three-dimensional pattern formation by any macroscopic objects and optimisation of their properties. The ‘fratons’ in this case being fragments of these objects are also macroscopic.

Materials and methods

The proposed fraton theory rests on two novel conceptual premises: (a) the introduction of interacting pseudoparticles that we call fratons that described two configurational states of each point of continuum space. One is an event in which the point is inside the atomic sphere of any atom and the other is an event in which the point is outside of atomic sphere; the fratons are considered as a non-ideal gas whose ‘condensation’ describes a diffusional self-assembling of atomic system, and (b) a concept of a structural cluster function describing the directions, length and strength of interatomic bonds. The latter allows us to formulate a new and simple model Hamiltonian that is proportional to a bilinear expansion in these cluster functions. This model Hamiltonian provides the formation of a predetermined atomic structure and has a sufficient flexibility to describe the desired mechanic and thermodynamic properties of this structure.

The numerical solution of kinetic equation (1) for the density function of atomic fractons was carried out by using the semi-implicit Fourier-spectral method in which the time variable is discretised using semi-implicit schemes and space variables are discretised using the Fourier-spectral method¹⁵.