Measuring crystallinity with Spherical Harmonics
In the study of molecular dynamics simulations, having methods to characterize local structural order of the constituents of the system is required to investigate phase transitions. In order to study processes like nucleation, dynamical arrest or crystallization that happen when molecular systems are cooled, or supercooled, bond-order parameters are used as a measure of local order or “crystallinity” in the system.
Say we want a reliable and quantitative way to differentiate a system that has undergone crystallization from one that has no structural order anywhere (like a liquid or a glass). In a crystalline solid the way in which neighbors arrange around a given particle should be the same for all the particles in the structure.
In analogy, one can think of a process like the calculation of the radial distribution function, where instead of measuring the distances to the n-nearest neighbors and averaging for every particle in the sample, one measures the directions in which the nearest neighbors of each particle arrange around it. When a given set of directions for the neighbors of each particle is found to be repeated throughout the sample, we say the system is crystalline.
SHPERICAL HARMONICS AS BOND-ORDER PARAMETERS
It is clear that as a first requirement, the parameters to be used should be rotationally invariant, as the orientation of the sample with respect to a chosen frame of reference shouldn’t bear on the measure of order between particles in the sample.
When there’s high order in the sample the parameter should take a high value, and when it is averaged over a sample with no order, it should be zero. Some combinations of spherical harmonics make good candidates for these purposes. For example, combinations of spherical harmonics with l=6 can be related to icosahedral symmetry, which can be related to nucleation processes; at least for L-J-like systems, when small clusters form, the particles tend to arrange in lower-energy non-crystallographic structures (like the icosahedron) rather than in crystallographic structures like fcc or hcp. The crystallographic structures don’t commonly appear until the clusters get hundreds of atoms in size.
To characterize the order around a particle i in the sample, one averages the spherical harmonics taken over the directions of its near neighbors j. In the definition of near neighbor it is customary to take particles that are at a distance of <1.2 times the diameter of the particles. The term “bonded” is commonly used (although not in the chemical sense). We define this average for the particle i as
with
Using such definitions, a scalar product that measures the correlation between the structures that surround two near particles i and j can be defined as
If this correlation takes a large value (higher than an adequate threshold), it means that the particles share a part of the same surroundings and one therefore says that the two particles i and j are connected. If a certain particle is connected to more than a certain number of particles (around 6 or 8), it is said to be solid-like.
In the figure above it can be seen how correlation between the surroundings of two particles decreases as these separate (into a more liquid-like structure).
Although this parameter can very easily distinguish between solid-like and liquid-like structures, it doesn’t discriminate between different crystal structures. In general, to get a rotationally-invariant combination of spherical harmonics we can use the addition theorem
where the value of the sum depends only on the angle between the vectors as
and if the product is taken for the same vector
we then get a normalized and rotationally-invariant factor
By using this combination with the averages over the bonds of a given particle i as previously expressed and taking the square root, we get the Steinhardt bond-order parameters
These parameters (specially l=4 and l=6) do hold the information of the local structure and take different specific values for fcc, bcc, hcp, sc and icosahedral structures, which allows to completely characterize the local order.
One can go further and define averaged Steinhardt order parameters as
where
where k runs over neighbors of the particle i and the average is taken over its neighbors plus the particle i itself. So the averaged versions of the Steinhardt order parameters take into account not only the first shell of neighbors, but also the second.
These definitions can be used to increase the accuracy when thermal fluctuations may smear out the order parameter distributions and make difficult to distinguish between crystalline structures; however it reduces the spatial resolution.
CODE & TEST Shown below are two different MD simulations. The one on the left is pure L-J argon, which crystalizes at low temperatures (crystalline planes start to appear on the faces of the simulation box). The one on the right is a system called Kob-Andersen mixture, which consists of two different types of argon-like particles used to avoid crystallization and study systems in dynamical arrest.
It can be seen that the first system crystallizes while the one on the right doesn’t get any apparent order.
The Steinhardt order parameter obtained for this Kob Anderson mixture is on the order of ~10^-2, which is way below from what is obtained for any crystalline system (~0.3-0.6). The reason the parameter isn’t completely zero is because cubic periodic boundaries in the simulation induce certain symmetry in the sample. This MD simulations were made using in GROMACS (to learn to program your own MD simulations see Computational simulations of dynamical systems).
Below is the fortran code used for the Steinhardt order parameter averaged over a complete system.
bond=1.2*sigma !neighbor definition with diameter sigma l=6 yp=0.0 do m= -l, l nb=0 !total number of bonds ymp=(0.0,0.0) !sum of parameters for each m (complex) do j= 1, nump !number of particles in the sample open(1,file=’input positions.dat’) !file with the positions of the particles in the sample do p= 1, nump read(1,*)x0,y0,z0 if(p .eq. j) then close(1) go to 10 end if end do 10 continue open(1,file=’input positions.dat’) do i= 1, nump read(1,*)x,y,z r=sqrt((x-x0)**2 + (y-y0)**2 + (z-z0)**2) if(r .le. bond) .and. r .gt. 0) then nb=nb + 1 theta=acos((z-z0)/r) if(sin(theta) .eq. 0) then phi=1.0 !if theta=0, phi is arbitrary else phi=atan((y-y0)/(x-x0)) end if Y6m=spherical_h(l,m,theta,phi) !call your spherical harmonics function ymp=ymp+Y6m !sum over the bonds end if end do mymp=abs(ymp) close(1,*) end do yp=yp + (mymp/np)**2 !sum over m end do Ql=sqrt(4*pi*yp/(2*l+1)) write(*,*)’Q_l=’,Ql,’, l=’,l END
References
P.J. Steinhardt, D. R. Nelson, and M. Ronchetti (1983) “Bond-orientatonal order in liquids and glasses”. https://doi.org/10.1103/PhysRevB.28.784. J. S. van Duijneveldt et al (1992). “Computer simulation study of free energy barriers in crystal nucleation”. http://dx.doi.org/10.1063/1.462802. W. Lechner and C. Dellago (2008). “Accurate determination of crystal structures based on averaged local bond order parameters”. http://dx.doi.org/10.1063/1.2977970.
W. Kob and H. C. Andersen (1995). "Testing mode-coupling theory for a supercooled binary Lennard-Jones: The van Hove correlation function" https://doi.org/10.1103/PhysRevE.51.4626.
Plotting the Spherical Harmonics in Python: http://balbuceosastropy.blogspot.mx/2015/06/spherical-harmonics-in-python.html Spherical Harmonics and platonic solids: https://nukephysik101.wordpress.com/2017/08/14/spherical-harmonics-and-platonic-solids/
Comentarios