Distinguishing noisy crystal symmetries in coarse-grained computer simulations: New procedures for noise reduction and lattice reconstruction
Evgeniia Filimonova Viktor Ivanov Timur Shakirov


twoside=false, left=4cm, right=2.5cm, top=1.5cm, bottom=2.5cm
  gobble
  plain
  
  
  A Unified Framework for Total Variation Regularized Optimization in Fluid Dynamics and Related Physical Systems
    mailto:vvarsha@purdue.eduvvarsha@purdue.edu
    April 29, 2024
===============================================================================================================

fancy
  

  gobble
  fancy
  
  arabic
  
    
  
  

  

	
	evgeniia.filimonova@physik.uni-halle.de
	Institute of Physics, Martin Luther University, 06099 Halle (Saale), Germany

	
	Institute of Physics, Martin Luther University, 06099 Halle (Saale), Germany
	
	
	Institute of Physics, Martin Luther University, 06099 Halle (Saale), Germany

Distinguishing noisy crystal symmetries in coarse-grained computer simulations]Distinguishing noisy crystal symmetries in coarse-grained computer simulations: New procedures for noise reduction and lattice reconstruction






We suggest new modification (we call it a noise reduction procedure) for Steinhardt parameters which are often used for detecting crystalline structures in computer simulation of solids and soft matter systems. We have also developed a new methodology how to reconstruct "ideal" lattice structure in the whole simulation box that would be most close to a real noisy crystalline symmetry, when it is defined locally and then averaged over the whole box. For this second procedure, which we call lattice reconstruction procedure, we have developed an algorithm for finding the lattice vectors from the values of Steinhardt parameters obtained after the noise reduction procedure. We apply noise to the classical crystalline structures (sc, bcc, fcc, hcp), and use both procedures to detect the crystalline structures in these classical but noisy systems. We demonstrate advantages of our procedures in comparison with existing methods and discuss their applicability limits. 


[
    T. Shakirov
    April 29, 2024
==================



§ INTRODUCTION



The model of hard spheres (HS) is widely used in theory and simulations of colloidal systems and it describes rather well experimental results on crystallization in colloidal systems <cit.>. The model of tangent hard spheres (THS) is very popular for polymer systems. It describes flexible polymers if tangent hard spheres are freely jointed and semiflexible polymers if some bending potential on the angles between rigid bonds adjacent along the chain is applied. Phase transitions in such systems, in particular crystallization and liquid crystalline transition, are entropy driven. The Monte Carlo method is very convenient for computer simulation of such models and has been widely used for more than three decades <cit.>.

A very important (if not the most important) question in the studies of crystallization phenomena in soft matter systems is how to distinguish structural motifs, in particular ordered symmetries/morphologies, e.g., crystalline and liquid crystalline ones, on local and global scales in conformations which are usually quite noisy and often inhomogeneous, i.e., the local particle density and values of order parameters can fluctuate strongly, especially in the vicinity of phase transitions. 

Usually, the radial distribution function (RDF) and static structure factor are calculated first. For example, for freely jointed hard-sphere chains one can see more pronounced maxima in the solid phase in comparison to the liquid one <cit.>, but one cannot reliably distinguish different crystalline symmetries by data on RDF.

Another parameters and methods of analysis for crystals are Voronoi tesselation <cit.>, common neighbor analysis <cit.>, nematic order parameter, chain segments local alignment parameter <cit.>.

Bond orientational order parameters or Steinhardt parameters <cit.> have been suggested long ago to distinguish crystals with different symmetries. However, these parameters work well only for ideal crystals or for crystals with very weak fluctuations in the particle positions, but they need modifications for typical cases in real soft matter systems, and such modifications have been suggested <cit.>. 

The characteristic crystallographic element (CCE) norm as a powerful descriptor of local structure in atomistic and particulate systems was introduced in Ref. LasoCCE2009. The CCE-norm is sensitive both to radial and orientational deviations from perfect local order. Recently, the CCE-norm has been revised and extended for reliable identification of local structure in 2d and 3d atomic systems <cit.>.
The local and global order in dense packings of semi-flexible chains consisting of tangent hard spheres has been studied in Ref. LasoPolym2023-1. It was analyzed how the packing density and chain stiffness influence the self-organization of chains at the local and global levels. The local order corresponds to crystallinity, and it was  quantified by the CCE-norm descriptor <cit.>, while the global order was computed through the scalar orientational order parameter <cit.>.

Recently, we have started investigation of surface phenomena in crystallization of semiflexible tangent hard spheres' chains by means of flat histogram Monte Carlo simulations <cit.>. We use the same model for which crystallization in the bulk has been studied before <cit.>, but consider melts at different substrates modelled by either purely repulsive walls or by attracting walls with various depth of well potential. In our simulations we use all above mentioned parameters, except CCE-norm.  
However, we have found that all these parameters, including modifications of Steinhardt parameters, are still not entirely accurate for the systems which we get in our simulations. Therefore, we have suggested another new modification for Steinhardt parameters and called it a noise reduction procedure. We have also developed a novel methodology how to reconstruct "ideal" lattice structure in the whole simulation box that would be most close to a real local noisy crystalline symmetry when averaged over the whole box. For this second procedure, which we call lattice reconstruction procedure, we have developed an algorithm for finding the lattice vectors from the values of Steinhardt parameters obtained after the noise reduction procedure. We apply noise to the classical crystalline structures (sc, bcc, fcc, hcp), and afterwards use both procedures to determine the crystalline structures in classical but noisy systems. 

Our paper is organized as follows. In Section <ref> we describe our model which we have used in our simulations of semiflexible chains of tangent hard spheres and present some examples of conformations which motivated us to developed the new modifications for structural analysis. We present here also our estimates of the noise amplitude in those systems, which we will use to apply noise to classical crystalline symmetries, It should be emphasized here, that our results on phase transitions in semiflexible THS systems will be published in a separate paper, and here we concentrate only on the methodological developments, which are (in our opinion) quite important for the field of computer simulations of soft matter systems, in particular for polymers, colloids and liquid crystals. In Section <ref> we present our results on structure analysis in noisy sc, bcc, fcc, hcp lattices by previously suggested methods and demonstrate the problems with correct detecting of the symmetries, Then, in Section <ref> we present our new noise reduction procedure and discuss the results of structure analysis by this method. A Section on our new lattice reconstruction procedure will be added soon. Section <ref> contains our conclusions. 
  

§ MODEL



Flexible tangent hard-sphere chain model used in this work is well established in works devoted to the study of crystallization in short chain melts <cit.>. 
The hard-sphere-type interactions of non-bonded beads:


    U_nb (r)=
    ∞,     r ≤σ ,
    
    0,     r > σ ,


where r is the distance between the centers of the two beads. 
This interaction does not give a numerical contribution to the energy, but imposes constraints on the available configuration space. 
As a result, possible angles between the bonds of neighboring spheres θ in the chain cannot exceed 120^o.
The lengths of chains in all studied systems are N=10 beads. 



Stiffness potential of the chains: 


    U_θ (θ)=
    
    - ϵ,    θ≤θ_s ,
    
    0,    θ_s < θ < 120^o  , 
    ∞,     120^o ≤θ ,


The ground-state stiffness energy for systems of
N_c chains of length N is defined as E^s_min = -N_c(N - 2) ϵ, while the maximum energy equal to E^s_max = 0. We consider several boxes with constant volume fraction ϕ≈0.496. The fixed value of cos(θ_s)=0.9 used in this study corresponds to  θ_s≈ 26^o, which at a volume fraction ϕ≈0.496 makes it possible to observe an ordered structure. 
During the simulation, the size of the box (L_x, L_y, and L_z) does not change.
Energies and temperatures are measured in units of well depth ϵ=1, while all lengths are measured in units of solid sphere diameter σ=1.



 §.§ System in the absence of walls

We simulate the system that we will use as a reference system of the polymer in the unconstrained bulk. Since in small size systems, the finite size effects of the system have a great impact, we use 2 simulating boxes: with sizes  L_x = L_y = 20 and L_z = 19 (N_c = 720) and  L_x = L_y = 20 and L_z = 39 (N_c = 1440). This ratio of the number of particles and the size of the box provides the desired volume fraction:


    ϕ=NN_c/L_x L_y L_zπσ^3/6≈0.496


Since there are no constraints and there are periodic boundary conditions in all directions, the total energy E is composed only of the stiffness energy E^s of the chains:

    E = E^s = ∑ U_θ.




 §.§ System with two purely repulsive walls

In the case of uniformly repulsive walls, a potential U_rep is applied between the bead and the wall :

    U_rep (z)=
    ∞,     |z|  ≥  (L_z - σ) / 2  ,
    
    0,     |z|   <   (L_z - σ) / 2,


where z represents a coordinate of a bead along z axis. 
In fact, this potential hinders chain units from passing through the wall and does not contribute numerically to the total energy value. That is, the contribution to the energy is determined by the same summand (Eq. <ref>) as in the case of the unconstrained system. 
As the centres of the beads cannot come closer to the surface than their own radius of σ/2, we increase the size of the box along the z axis to the value of σ in order to ensure that we preserve the bulk volume fraction of the polymer, as in the system without walls (ϕ≈ 0.496).
Therefore, in the presence of walls, the dimensions are L_x = L_y = L_z = 20 (N_c = 720) and L_x = L_y = 20 and L_z = 40 (N_c = 1440). 
Thus, two parallel walls are positioned on the planes at z=± 10 and z=± 20 respectively. Along the x and y axes, periodic boundary conditions are implemented.



 §.§ System with purely repulsive and attractive walls

We also examine two scenarios involving a small system ( L_x = L_y = L_z = 20, N_c = 720 ) where one wall exhibits repulsion and is conserved at (z = -10), while the second wall has an attracting potential at (z = 10): 


    U_at_1 (z)=
    
    
    - ϵ,    L_z / 2 - σ ≤  z   <   (L_z - σ) / 2 ,
    
    0,    - (L_z - σ) / 2   <   z   <   L_z / 2 - σ,
    ∞,    |z|  ≥  (L_z - σ) / 2



    U_at_2 (z)=
    
    
    - 4 ϵ,    L_z / 2 - σ ≤  z   <   (L_z - σ) / 2  ,
    
    0,    - (L_z - σ) / 2   <   z   <   L_z/ 2 - σ,
    ∞,    |z|  ≥  (L_z - σ) / 2



  

§ STRUCTURE ANALYSIS



In general, a crystal means a solid body having a three-dimensional long-range translational order. 
The arrangement of atoms in a crystal is characterized by its unit cell containing one or more atoms in a certain spatial arrangement. 
The crystal structure of the substance is reproduced by periodic translation of the unit cell. 
In this study, we do not pretend to determine the exact class of the crystal, since this would represent a separate in-depth research in the field of crystallography. 
In addition, in polymer systems, we can only talk about local translational ordering and the formation of a set of crystallites. 
However, we claim to be able to detect the coexistence of various crystal symmetries.
Despite the fact that a significant number of methods for analyzing crystal structures have been developed, this task is still the cornerstone of computer simulation research. In this chapter, we will look at the most common methods of analysis applied to frequently considered systems such as simple cubic (sc), body-centered cubic (bcc), face-centered cubic (fcc), hexagonal close-packed (hcp), as well as to our system and propose another method of analysis that we consider the most successful in this study.
  

 §.§ Noise estimation



Before starting to analyze the crystal structure, it is necessary to pay attention to the fact that one cannot expect to obtain an ideal crystal within the framework of the model used.
In the studied systems, noise will be observed due to the breadth of the potential.
Since the condition for adding the value -1 is an angle ranging from 0 to 26 degrees, then the lowest energy of E=-5760 in a small system will be given by both, as well as fully elongated chains with all angles between the chains equal to 0, and elongated chains with all angles close to 26 from below. 
Thus, the width of the potential is the cause of noise in the resulting crystal structures. 
This paragraph will describe the procedure for evaluating noise in our system.

Let's choose a system with a sufficiently low energy. For example, consider a system with energy E=-5727 (Fig.<ref>), for other energies and systems, estimates give a similar result. 



While in the XY and YZ planes one can notice a variety of patterned structures of particles, in the XZ plane one can observe uniting into groups. 
In the above case, the organization of projections into groups in this plane is explained by the elongation of the chains along the y axis (Fig.<ref> a). Thus, the center of the groups has the meaning of the atom of the crystal lattice, and the size of the group characterizes the noise in the crystal.
At this stage, clustering was carried out using the k-means method <cit.> (with an input number of 360 groups) and a method similar to k-means, however, not the number of groups was set as an input parameter, but the initial roughened size of the groups (R=0.8). Both methods give the same result for groups positions.
As can be seen from the inserts in Fig. <ref>, the distribution of deviations of points from the centers of groups resembles the Gaussian distribution. 
One can use quantile-quantile plot (QQ plot) to check normality visually. 
QQ plot draw a correlation between the sample and the normal distribution. 
In a QQ plot, each observation is displayed as a single point. 
If the data correspond to the norm, then the dots should form a straight line. 
Thus, in our case one can talk about the normality of the distribution. 
In addition to the ability to visually assess whether the distribution belongs to the normal distribution, this method allows estimating the median value of the distribution and the standard deviation. 
The value on the y axis at x=0 corresponds to the median. 
We see that it is equal to zero; this correctly reflects our construction: we translated the center of each group to the origin before counting the deviation of the position of the particle of each group from its center. 
The tangent of the slope of the approximating lines corresponds to the standard deviation of the distribution. 
From the resulting plots we can extract √(σ^2)≈ 0.112. 
This value characterizes the noise in our system. 



Another way to estimate the parameter √(σ^2), assuming that the distribution of deviations of particles from the centers of groups D(x) is normal, is to construct a lnD(x) :



    D(x) ∼exp{- (x-x_0)^2/2σ^2}


    ln D(x) = const_1 - (x-x_0)^2/2σ^2

Since we assume that the centers of all groups have already been translated to the origin, then x_0 = 0.

    ∂ln D(x)/∂ x = - x/σ^2

As can be calculated from the coefficients obtained in the linear model √(σ^2) is in both cases: 

    √(σ^2)≈ 0.112,   (σ^2 ≈ 0.0125 ).


Taking into consideration that this is an estimate for one group along one direction, for a pair of noisy groups in three-dimensional space we get 6σ^2 ≈ 0.075.
We will use this estimation in the following parts of the work. 



  

 §.§ Test structures preparation



In order to make sure that methods used are reliable, we will perform tests on well-known structures: sc, bcc, fcc, hcp.
For these structures primitive vectors are known, by which it is possible to reproduce the crystal lattice. 
In order that the results obtained for the test lattices could be compared with the structures studied in this work, two conditions were taken into account: 
(a) the minimum distance between the lattice atoms of the perfect test structures should be equal to the length of the rigid bond in the polymer model under study, i.e. d=1; 
(b) the atomic packing factor (APF) of the test structures should coincide with the volume fraction of the polymer in the system under study (APF = ϕ = 0.496), since this parameter is the same for all the systems under study and does not change during simulating process.

Since the dependence of various parameters on noise √(σ^2) is investigated and special attention is paid to the value estimated in the section <ref> (√(σ^2)=0.11), here we will also describe the procedure for introducing noise √(σ^2) into the test structures taking into account the excluded volume.

Example

As an example, let us consider the construction of bcc structure on the basis of which comparisons of the investigated parameters will be made. 
The construction of other test structures is carried out in full analogy.

Primitive lattice vectors for bcc: 

    a = -a/2e_1 + a/2e_2 + a/2e_3, 
    b = a/2e_1 - a/2e_2 + a/2e_3, 
    c = a/2e_1 + a/2e_2 - a/2e_3,


where a is the lattice parameter, or the side of the cube, related to the distance between the nearest lattice atoms d by the relation

    a=2/√(3)d.

Usually, to calculate the  APF one assume that d=2r, where r is the radius of the particle.
However, to adjust the APF, we fix d=1, according to assumption (a) and change the particle radius r. APF_bcc is expressed in terms of the volume of a single particle V_1, the number of particles in a unit cell N_bcc and the volume of a unit cell V_0:


    APF_bcc=V_1 ·  N_bcc/V_0 = 4/3π r^3 · 2/a^3.


After using the condition Eq. <ref>, and the requirements (a) and (b) that APF_bcc = ϕ = 0.496 from Eq. <ref>, one obtains:


    r = (ϕ/π√(3))^1/3≈ 0.45.


This value is particularly important in the procedure for the noise introducing to the hard sphere system.

Procedure for introducing noise

We assume that the noise along any direction in our system has a Gaussian distribution. 

Therefore, we choose orthogonal vectors e_1, e_2, e_3, and we will drive noise along each of these directions. 

Step 1

At the first step in the constructed ideal lattice, each particle acquires a normally distributed Gaussian displacement 


    Δ x_1 = 𝒩 (0,σ^2)

Since at the initial stage we did not control the overlap of particles, now we need to eliminate it.  

Step 2

We are looking for a pair of points i_0, j_0 that form the minimum pair distance d_min in the resulting noisy system. 

Step 3

Randomly select one particle i ^0 from i_0 and j_0. With a probability of 0.1%, the particle i ^0 is randomly selected from all the particles of the system. For the particle i ^0, the displacement introduced earlier is replaced by:


    Δ x_2 ⟶αΔ x_2 + 𝒩 (0,σ_1^2)


where α∈ [0;1]. 
According to the properties of the normal distribution, the distribution Δ x_2 will be distributed normally with a width of √(σ^2) if one sets:


    σ_1^2 = (1-α^2)σ^2.


An attempt to change the offset (Eq. <ref>) is accepted if the paired distances to the particle i ^0 taking into account this step are not less than d_min, otherwise i ^0 retains its previous offset Δ x_1. 

Step 4

If a new displacement Δ x_2 is accepted, a new d_min and particles i_0, j_0 are searched.

Step 5

Repeat Step 2 - Step 5 until d_min < 2 r, where r was evaluated in Eq. <ref>.

Remark

After completing the procedure, the resulting structures were checked for the normality of the distribution of the position of the particles relative to the noiseless positions, similar to how it was done in the section <ref>.
The test showed that this method works well in our task for small values of (√(σ^2) < 0.07). 
When √(σ^2)> 0.07, the dependence of the observed noise on the √(σ^2) value we introduce during the procedure is no longer linear. 
For example, for hcp, the observed noise becomes comparable to that observed in our system (Eq. <ref>) at induced √(σ^2) = 0.16. 
Later in the text, speaking about noise and using the notation σ, we will use the actually observed values obtained during this check.




  

 §.§ Radial distribution function



Let's start the study of structures by calculating the radial distribution function (RDF). 
Let G(r) be the probability density of the presence of a particle at a distance r from a given particle. 
Then the probability of the presence of a particle at a distance r from this particle is defined as G(r)dr. 
Then the radial distribution function g_2(r) is calculated as:


    g_2(r) = G(r)dr/4π r^2 dr V,

where V is the volume under study. 




As it was established in the section <ref>, the formation of a perfect crystal cannot be expected in the structure under the study.
Noise can contribute to the displacement and blurring of the radial distribution function.
To study the behavior of the g_2(r) depending on noise, sc, bcc, fcc, hcp structures were selected as test structures.
Gaussian noise √(σ^2) is introduced into the initially perfect crystal lattice in each of the three directions, as it was described in the section <ref>. 
The red curves (Fig.<ref>) correspond to noise similar to the system under study, since noise √(σ^2)≈0.11 is introduced in this case. It can be seen from the data obtained that the number of peaks on the RDF decreases and the determination of the symmetry of the structure becomes impossible. 
The position of the first maxima shifts to the left on all structures with increasing noise. For small √(σ^2) values, the maxima of all structures are located at the point r=1.0, which corresponds to the method of constructing lattices described in section <ref>. 
This value corresponds to the minimum distance between the lattice nodes. With a sufficiently large √(σ^2), the peak corresponds to the minimum possible distance that the lattice atoms can approach, that is, the diameter of solid spheres. 

On the graphs of g_2(r) (Fig.<ref>) for the system under study (at energies E =-2040 and E=-5727, which corresponds to a melt and an ordered structure respectively), peaks at r=1 and r=2 can be observed. These maxima correspond to neighboring particles along the chain in our model. 

The calculation of the RDF is a necessary step for calculating the Steinhardt parameters. 
In general, Steinhardt and co-authors <cit.> recommend using particles that fall into the sphere of the cutoff radius R_c=1.2r_0 as the nearest neighbors, where r_0 is the position of the first peak of the g_2(r). 
Such a choice should ensure that all particles in the first coordination sphere are taken into account. 
However, since the noise estimation in our system is 6σ^2=0.27, it was decided to use the R_c=1.3r_0. Since r_0=1 coincides with the diameter of the hard spheres in further text we will immediately write R_c=1.3. Our experience has shown that this value is optimal for calculating parameters. 
If we consider the distributions for the model under study (Fig.<ref>), one can see that the selected value of the R_c is to the left of both minima. 
As this value is close to the minima of both our structure and the test ones, this choice is the most optimal, since it is most likely to capture all particles from the first coordination sphere (R_c=1.3 magenta in the Fig. <ref>).


  

 §.§ Local bond order parameters




  §.§.§ P. J. Steinhardt, R. Nelson and M. Ronchetti parameters



One of the most discussed problems of simulation of the crystallization process is how to attribute a particle to a liquid or crystal and also determine the type of crystal. 
Recently, a widely used method of distinguishing the type of particles is local bond order parameters, also known as Steinhardt parameters <cit.>. 
To calculate these parameters, it is not necessary to have a reference structure for comparison with the studied one, as is required in, for example, common neighbor analysis <cit.>. 
Another convenience of using this method is the absence of reference to the coordinate system, since the calculation of these parameters is based on spherical harmonics. 

Local bond order parameter of particle i is defined as:


    q_l(i) = √(4π/2l+1∑_m=-l^l q_lm(i)^2 ),


where l is the order of parameter, q_lm - complex vector which is defined as


    q_lm(i) = 1/N_b(i)∑_j=1^N_b(i) Y_lm(r_ij).


Here, N_b(i) is the number of nearest neighbors of the particle i. There is no strict rule for determining the particles that are considered to be the nearest neighbors. 
Most studies use the concept of the cutoff radius R_c (however, different studies use different values), and some use the Voronoi cell.
We prefer to use the concept of the cutoff radius with R_c=1.3, as discussed in section <ref>. 
The functions Y_lm(r_ij) are the spherical harmonics, r_ij is the vector connecting the particles i and j; m is an integer that runs from m=-l to m=l.
Along with q_l, w_l parameters are also often used, which can be computed according to formulas:


    w_l(i) = ∑_m_1+m_2+m_3=0[   l   l   l; m_1 m_2 m_3 ]
        q_lm_1(i)q_lm_2(i)q_lm_3(i)
         /( ∑_m=-l^l q_lm(i)^2 )^3/2 ,


where the summation is carried out by the integers m_1, m_2, m_3 from -l to l, which satisfy the condition m_1+m_2+m_3=0. 
The expression in parentheses is the Wigner 3-j symbol. 
Since no significant conclusions can be drawn from the calculation of the parameters w_i, the results for these parameters are presented in supplementary material (subsection <ref>).
To analyze the structure, q_4, q_6, q_8 are most often used. 
In this section, we will present the results for the parameters of test structures with noise, and also apply this method to analyze the structures obtained during our simulation.



The study of Steinhardt parameters for test structures with the value √(σ^2)=0.11 showed that the structures hcp, fcc, bcc become absolutely indistinguishable (Fig. <ref>). 
This fact does not allow us to apply parameters for our system. 

The study of the dependence of the average parameters on √(σ^2) showed that the mean values with increasing noise deviate greatly from the values for ideal crystal lattices (horizontal lines in Fig. <ref>). 
Starting from √(σ^2)≈ 0.08, the average number of nearest neighbors also becomes indistinguishable for bcc and hcp, fcc and for none of these values corresponds to the true number of neighbors.





Thus, for this model, the application of classical Steinhardt bond order parameters turns out to be impossible. 
In all non-ideal systems, where there is a deviation of the positions of particles (for example, due to thermal fluctuations) from perfect crystal lattices, difficulties arise with determining the crystal structure in simulation. Nevertheless, the potential of the Steinhardt approach is great. 
Spherical invariants are reasonable functions for describing the symmetry of particle clusters.
A number of works have been devoted to modifications of the local bond order parameters to improve the accuracy of determining the type of structures in the presence of thermal fluctuations. 
In this chapter, we will look at several well-known modifications of this method.



  §.§.§ W. Lechner and C. Dellago parameters



In 2008, W. Lechner and C. Dellago <cit.> proposed a procedure for averaging Steinhardt parameters. 
The authors tested their results on two different systems of soft spheres.  
In this study , they propose to average local bond parameters as follows:


    q̅_l(i)=√(4π/2l+1∑^l_m=-l|q̅_lm(i)|^2 ),


where 


    q̅_lm(i) = 1/Ñ_b(i)∑^Ñ_b(i)_k=0 q_lm(k).


The equation involves a sum over all neighbors of particle i, encompassing the particle itself, ranging from k=0 to Ñ_b(i). 
To determine the local orientational order vectors for particle i, one takes the average of q_lm across both the particle i itself and its nearby environment Ñ_b(i). 
While q_l(i) reveals the structural intricacies of the first shell around particle i, its averaged counterpart q̅_l(i) accounts for the influence of the second shell. 
The effective considering of the second particle shell is critical in this context.

In the supplementary material (subsection <ref>) you can also find the results for the updated parameters w̅_̅l̅ calculated using q̅_lm:


    w̅_l(i) = ∑_m_1+m_2+m_3=0[   l   l   l; m_1 m_2 m_3 ]q̅_lm_1(i)q̅_lm_2(i)q̅_lm_3(i)
         /( ∑_m=-l^l q̅_lm(i)^2 )^3/2 ,


Authors <cit.> conducted their tests on 2 types of soft spheres: with the Lennard-Jones potential and the Gaussian core model. 
We have calculated these parameters for our test structures consisting of solid spheres.




From the results obtained (Fig. <ref>) it can be seen that the average values still vary greatly with the growth of √(σ^2), however, averaging definitely reduces the spread in values. 
On the Fig. <ref> (d) it can be seen that the average values of the parameter q_4 are able to distinguish bcc, fcc, hcp, although it is still unreliable. 
The average value for the bcc structure has stabilized and slightly deviates from the undisturbed state, which was not observed with the original method of calculating the Steinhardt parameters Fig. <ref> (a). 
The parameters q_6, q_8 still do not allow us to distinguish structures with noise values comparable to those observed in the system under study (√(σ)^2≈ 0.11).



  §.§.§ H. Eslami, P. Sedaghat and F. Müller-Plathe parameters



A more advanced method was proposed by H. Eslami et al. <cit.>. 
In this paper, using the example of the Lennard-Jones system (12-6), a comparison was made with the parameters of previous authors <cit.>. 
These order parameters examine the ratio of the orientational orders of the second-shell to the first-shell neighbors of a central particle, so the parameters change from 0 (disordered structure) to 1 (in crystal). 
The authors note that despite the fact that the parameters in crystal structures are close to 1, thermal fluctuations reduce the parameters. 
Moreover, it is stronger in bcc than, for example, in fcc, due to the large amount of free space and lower density. 
At the same time, no studies have been conducted on the dependence of parameter changes on temperature fluctuations.

These parameters can be calculated using the following formulas:


    q̃_l(i)=q_l(i)·q_l(j)=1/N_b(i)∑_j∈ N_b(i)∑^l_m=-lq̂_lm(i)q̂_lm^*(j),


where


    q̂_lm(i)=q_lm(i)/( ∑^l_m=-l | q_lm(i)|^2 )^1/2.



After the averaging q̃_l(i) over the first coordination shell neighbors of particle i one get the following:


    q̅̃̅_l(i)=1/Ñ_b(i)∑_j∈Ñ_b(i)q̃_l(j).



Despite the fact that these parameters are not intended to distinguish ideal (noiseless) crystal structures, since they are equal to 1, it is interesting to investigate the dependence of these parameters on noise for systems of solid spheres.
There is a possibility that using a combination of parameters q_4, q_6, q_8 will allow distinguishing different types of structures by the degree of deviation of the parameters from 1. 




From the obtained data (Fig. <ref>), it can be seen that the parameter q̅̃̅_4 still seems to be the most sensitive for separating structures. 
Depending on the type of structure, noise has a different effect on the deviation of parameters from undisturbed values.
Such a wide range of values, unfortunately, does not allow using combinations of parameters q_4, q_6, q_8 to distinguish between different types of structures.

Thus, since the parameters discussed in the section <ref> do not allow us to reliably distinguish test structures with a noise value comparable to that estimated in the system under study, we cannot rely on them when analyzing the system under study. 
In the next section, another way of using Steinhardt parameters (Eq. <ref>) will be proposed, the efficiency of the method on test structures will be checked, and the use of the method on the example of the system under study will also be demonstrated.
  

 §.§ Noise reduction procedure



As we have established in the previous sections, the methods used for analyzing structures do not allow us to distinguish satisfactorily different symmetries in the presence of noise. 
Therefore, at this stage of the work, an attempt was made to reduce the influence of noise by averaging the position of coordinates in space.



We propose the following procedure (Fig. <ref>). 
Let's choose the particle i for which one will compute Steinhardt order parameters. 
We draw a sphere of radius R_s=2 around the particle (Fig. <ref> a). 
Let it have N_i neighbors of particle i. Next, select particle j (green triangle on Fig. <ref> b) from the list of neighbors {N_i}. 
We also draw a sphere of radius R_s=2 around it, into which N_j particles fall (blue particles on Fig. <ref>b). 
Then we translate particle j together with all its neighbors into particle i (Fig. <ref> c). 
We perform the same operation with all the particles from the list {N_i} (Fig. <ref> d). 

The next step after moving the points inside the sphere is averaging their positions. 
To do this, we set the positions of the neighbors of the initial particle as the starting position of the groups r'_gr (black dots inside the black sphere on the Fig.<ref> a)). 
Then we iteratively refine the present position of the centers of the groups r_gr.  To do this, add the position of each of the above points (blue on the Fig. <ref> d)) to each of the existing groups as an gaussian term. To get the updated coordinates of the groups, it is necessary to weigh the amounts received.


    r_gr = 1/Z_gr∑_k r_k  exp-(r_k-r'_gr)^2/2· 2σ^2


    Z_gr = ∑_k  exp-(r_k-r'_gr)^2/2· 2σ^2,


where r_k is the position of the particle relative to the particle under consideration i, r'_gr are group positions at previous iteration step, r_gr are updated group positions, 2σ^2 = 0.0242, as defined in the section <ref>.
We repeat this procedure a number of times, achieving convergence of the position of the groups. 
Thus, we obtained the averaged positions of the neighbors of the initial point i, which we use to calculate the Steinhardt parameters (Eq. <ref>). 

It should be noted that for averaging coordinates, the radius of the sphere (R_s=2) is set obviously larger than the cutoff radius for calculating the Steinhardt parameters, which is still equal to r_c=1.3. 
During iterations, the coordinates of the centers of the groups r_gr, which will be used as the coordinates of neighbors at the end, may shift. 
Thus, if the original point i had an insufficient or excessive number of points inside the sphere R_c=1.3, in the process of iterations, this amount changes and becomes less sensitive to noise around the original particle. Thus, the choice of R_s=2 provides a more reliable averaging.



As we can see from the results obtained (Fig. <ref>), the procedure of averaging coordinates in space makes it possible to make the mean values of all parameters more stable with increasing noise. 
At the same time, the spread of order parameters naturally decreases. 
On the planes q_4-q_6, q_4-q_8, q_6-q_8 (Fig. <ref> a - c) it is also clearly seen that when introducing displacements in the lattice coordinates of the order √(σ^2) = 0. 11 the structures are clearly distinguishable. 
The same pictures show the division of sc into 3 groups. 
It is established that the reason for this is the different number of neighbors. 
Thus, one can conclude that the combination of three parameters q_4, q_6, q_8 allows to distinguish different symmetries after the averaging procedure, even if there is a sufficiently large noise comparable to that observed in the system under study.

Another modification of this method was also tested. Instead of combining the particles of the black sphere during averaging procedure by translations (Fig.<ref> a), the coordinates of the centers of mass for each of the spheres (black and all blue ones) were superimposed. 
Thus, one can say that the displacement of the central red particle is also taken into account. 
However, our experience has shown that this modification does not provide significant improvements. 
Initially, the described method can be considered optimal, and it will be used in further analysis. 
The results for the superposition of the centers of mass can be found in supplementary materials (<ref>).

Since the proposed approach has proven itself well in relation to test structures, we have every reason to apply it to the system under study. 
To begin with, let's take a system without walls with low energy E=-5727 (Fig. <ref>). 
The obtained parameters after the averaging procedure are shown in Fig. <ref>.



From the data obtained, it can be seen that the parameters q_4, q_6, q_8 are clearly divided into several groups. 
Two dense clusters are allocated on the plane q_6-q_8. 
On the other two planes, smaller clumps are also noticeable. 
However, since their size is small and they are located close to larger ones, we will consider them part of larger clusters. 
It was also found that most of the particles in two biggest the most dense clusters have 12 neighbors, while scattered points have a predominantly different number of neighbors than 12. 
Since the predominance of the correct number of neighbors can be observed on the test structures after the averaging procedure, then at this stage we can assert that there is a coexistence of 2 types of symmetries with 12 neighbors in this structure, while a different number of neighbors is formed due to the noise of the structure. 
Here we make the assumption that these 2 crystal structures originate at higher energies, where they coexist with the disordered phase, which we also will call Melt. 
We present an algorithm for isolating co-existing structures at different energies.

Procedure for structure separation

Step 1

At the first stage, we calculated the parameters for the averaged coordinates of a low-energy system (E=-5727) and divided the particles into 2 clusters using the k-means method <cit.> built into the R environment. 
q_4, q_6, q_8 were used as input data for clustering to increase accuracy. 
At the same time, our experience has shown that additional consideration of neighbors as input data does not significantly affect the result.
The number of clusters (2) and the number of launches (n_start = 500, 1000) are used as external parameters for the method. 
The main disadvantage of this method is a strong dependence on the starting points chosen randomly. 
To get around this problem, one should do several launches, and it is also better to visually observe the results of the cluster separation on the graphs. 
Despite the aforementioned drawback, within the framework of our task, the results are well reproducible and do not change when restarted. 
One can also use other clustering methods (for example, various variations of the k-means method <cit.>, the c-means method, EM-clustering). 
Each of these methods has its own advantages and disadvantages, however, in our work we will not describe and compare clustering algorithms in detail.
We prefer to use k-means method because of its simplicity, but testing other methods has yielded similar results.
After the initial separation of the parameters of the low-energy structure into 2 clusters, it is necessary to clarify the positions of the clusters centers in the space q_4, q_6, q_8. 
To take into account the asymmetry of the clusters, we calculate the covariance matrix:


    Ĉ_k = [  cov ( q_4^k,  q_4^k )  cov ( q_4^k,  q_6^k )  cov ( q_4^k,  q_8^k );  cov ( q_6^k,  q_4^k )  cov ( q_6^k,  q_6^k ) cov  ( q_6^k,  q_8^k );  cov ( q_4^k,  q_8^k )  cov ( q_6^k,  q_8^k )  cov ( q_8^k,  q_8^k ) ],


here k = {1,2} is the number of cluster; q_j^k - the parameter of the particle related to the cluster k; j = {4,6,8} is the order of the parameter. 
Initially, the mean value of the parameter in the initial cluster is taken as the center of the cluster, then this value will be refined during the iterative process. 
In the iterative process the weight with which each particle i of the structure enters each cluster k is calculated:


    Ω_i^k = √(det  Ĉ_k^-1)exp -1/2 (q_i,4 -q̅_4^k,  q_i,6 - q̅_6^k,  q_i,8 - q̅_8^k )   Ĉ_k^-1[ q_i,4 -q̅_4^k; q_i,6 -q̅_6^k; q_i,8 -q̅_8^k ],


where Ĉ_k^-1 is the inverse matrix of Ĉ_k; q̅_j^k - the center of cluster k. 
Then the particles are reassigned to clusters: each particle i belongs to the cluster k for which the weight turned out to be greater. 
After that, the positions of cluster centers are updated as:


    q̅_j^k = ∑_i q_i,j^k Ω_i^k/∑_i Ω_i^k,


where j={4,6,8}; i - the number of particle form the cluster k. 
The iterative process ends at the moment when q̅_j^k stops to change. 

Note: 
The Ĉ_k and Ĉ_k^-1 do not change during the iterative process. 
In the case of matrix recalculation, we will get nested clusters with close centers.

Thus, at the Step 1, the centers of 2 crystal structures (q̅_j^k, k={1,2}, j={4,6,8} ) and covariance matrices (Ĉ_k, k={1,2}) for them were obtained. 
These values are fixed and do not change in further steps.

Step 2.

After 2 types of particles with different symmetries have been determined for this structure, we propose that these types are formed at higher energies and coexist with the polymer melt. 
For the 2 found crystal states we fix the found values of cluster centers, as well as covariance matrices. 
At this step, our task is to find the values (q̅_4^3, q̅_6^3, q̅_8^3, Ĉ_3) corresponding to the polymer melt. 
To do this, we choose a high-energy structure, such that only a melt is observed. 
For example, let's take the energy E = -2040. 
For simplicity, assuming that there are no crystalline particles in the sample, we find q̅_4^3, q̅_6^3, q̅_8^3 in the sample as the mean parameters for the entire system:


    q̅_j^3 = 1/N · N_c∑_i q_i,j ,

where N· N_c is the number of particles in the system (N · N_c = 7200 in small and N· N_c = 14400 in large systems, respectively).

And based on these data, we build a covariance matrix (Eq. <ref>).

Thus, at the Step 2, the center of disordered structure (q̅_j^3, j={4,6,8} ) and covariance matrix (Ĉ_3) for it was obtained. 
These values are fixed and do not change in further step. 
One can find results for q̅_j^k (k = {1,2,3}) in the Tab. <ref>.




Step 3.

Since we know the coordinates in the {q_4, q_6, q_8} space for the three phases and their covariance matrices, it is now possible to distribute the particles in the system under study of any energy into 3 types. 
First of all, we can return to the low-energy structure and redistribute using weights (Eq. <ref>) the particles over the 3 sets of {q̅_4, q̅_6, q̅_8 } found. 
Since the system we have chosen does not have the lowest of all possible energies (E = -5760), the presence of several isotropic phase particles is possible. 





Remark

To make sure that the sets obtained in the Tab. <ref> for these 3 clusters do not depend significantly on the choice of specific structures, additional tests were conducted. 
The procedure of searching q̅_j^k and Ĉ_k (j={4,6,8}, k={1,2,3}) was repeated for several sets of low and high energy conformations. 
The set of conformations obtained during SAMC can be considered independent, since they were taken with a large time difference (on the order of several months of computer calculations). 
The results (Step 1, Step 2) of several tests for the system described here can be found in the supplementary material <ref>. The values q̅_j^k and Ĉ_k do not depend significantly on the choice of conformations. 
The assignment of a part to different clusters (Step 3) also does not change significantly.



In the Fig. <ref>, the distributions of the parameters q_4, q_6, q_8 for the system under study are shown in blue. 
Black, red and green represent the distributions of parameters weighted by the fraction of particles of the corresponding cluster. 
The dotted line of the cyan color represents the sum of these individual contributions. 
It is obvious that the distribution of parameters for all particles is actually a set of distributions of particles belonging to individual clusters, since the blue curve exactly coincides with the cyane line. 
The shoulder in the black curve at q_8 distribution for Crystal 1 (Fig. <ref> c) corresponds to the number of neighbors other than 12 (detailed distributions can be seen in the additional materials <ref>).








The method was also tested on a system with repulsive walls ( Tab. <ref>). 
In these two studied systems, structures with identical sets of q̅_j^k are formed. 
If we compare the spatial arrangement of the structures of Tab. <ref> and Tab. <ref> (see Supplementary materials below), the dominant crystalline phase in both systems is the Crystal 2 (red). 
In a system with repulsive walls, the alternation of crystal phases is observed parallel to the walls, while in a system without walls, alternation occurs in the planes located at an angle to the planes of the simulation box.
The conformations are dominated by layered crystal phases, with the Crystal 2 dominating in the vicinity of the walls.

The received data Tab. <ref> it is said that the structures found do not exactly define any of the test structures (sc, bcc, fcc, hcp). 
However, taking into account the estimated spread of the values obtained, it can be seen that the parameters of the Crystal 2 (red structure) are closest to fcc, while the Crystal 1 (black structure) is close in parameters to hcp. 
The only significant difference is observed for the parameter q_4 of the Crystal 1 and hcp. 
The deviation of this parameter from the reference value indicates the deformation of the structure. 
Indeed, when considering the projections of the centers of the particles of the system on different planes, irregular, deformed hexagons were observed. 
The reason for the deformation is the connectivity in the chain and the dimensions of the box. So, for example, if in low energy chains are stretched along the y axis one after the other, then in the perpendicular plane XZ we see sections of only 360 chains. 
By no such transformations, it is impossible to stretch particles stacked in a regular polygon on the plane to fit it into a rectangle of size L_x · L_z = 20 · 19. 
Note that although in reality the y direction is not highlighted, the final, small size of the system has a strong influence on the laying in our case. 
Nevertheless, the organization of short chains into layers was tested on a similar system in the work of T. Shakirov<cit.>.


Thus, it can be argued that the approaches developed in this work are valid within the framework of this model.



  

§ CONCLUSIONS



We have suggested two new procedures of analysis of noisy crystalline structures which appear in computer simulations of soft matter systems. 

The first of these two procedures – a noise reduction – is a special averaging over neighbors of a chosen particle before calculation of bond local order parameters (Steinhardt parameters). This procedure leads to an essential reduction of noise in particles' positions and allows to determine the local crystalline symmetry more reliably. 

Another procedure – lattice reconstruction – allows to reconstruct the “ideal” lattice structure in the whole simulation box that would be most close to a real noisy crystalline symmetry when it is first detected locally and then averaged over the whole box.

We plan to apply both these procedures to analysis of crystallization transitions in the melts and in thin films of short semiflexible chains of tangent hard spheres. 


  
  
  
  

  
  =10000
  
  
  
    unsrtnat 






twoside=false, left=2cm, right=1.5cm, top=2cm, bottom=2cm



    
