namespace chill

Overview

CHILL and CHILL+ structure classification. More…

namespace chill {
 
// enums
 
enum CrystalKind;
 
// structs
 
struct BondClassifier;
struct LinearClassifier;
struct QlmAtom;
struct SteinhardtQl;
struct TemplateHit;
struct VoronoiWeights;
struct YlmAtom;
 
// global functions
 
std::vector<TemplateHit> classifyTemplates(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, const std::vector<std::vector<int>>& nList, int kNeigh = 12);
std::vector<double> soapSpectrum(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, int iatom, const std::vector<std::vector<int>>& nList, int nMax, int lMax, double rcut);
std::vector<std::vector<double>> soapSpectrumAll(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, const std::vector<std::vector<int>>& nList, int nMax, int lMax, double rcut);
std::vector<double> voronoiFeature(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, int iatom, double candidateCutoff);
std::vector<std::vector<double>> voronoiFeatures(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, double candidateCutoff);
std::vector<VoronoiWeights> voronoiFacetWeights(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, double candidateCutoff);
SteinhardtQl steinhardtQlVoronoi(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, double candidateCutoff, int orderL);
 
} // namespace chill

Detailed Documentation

CHILL and CHILL+ structure classification.

This namespace contains functions that are used in the CHILL/CHILL+ classification scheme, as well as a yodaCloud struct to hold all the information.

Both the CHILL and CHILL+ methods are based on the local bond order parameter method, developed by ten Wolde et al. for the identification of crystal nuclei in Lennard-Jones systems. The local order parameter method is based on order parameters introduced by Steinhardt et al.

The local environment of the water molecules is classified using an algorithm based on spherical harmonics, which is independent of the specific crystal structure and does not require the definition of a reference frame. The local order around each water molecule \(i\) is described by a local orientational bond order parameter \(q_{lm}(i)\), with \(2l+1\) complex components:

\[q_{lm}(i) = @frac{1}{N_b(i)} \Sigma_{j=1}^{N_b(i)} Y_{lm}(r_{ij})\]

Here, \(N_b(i)=4\) is the number of nearest neighbours for water molecule \(i\), \(l\) is a free integer parameter, \(m\) is an integer such that \(m \in [-l,l]\). The functions \(Y_{lm}(r_{ij})\) are the spherical harmonics and \(r_{ij}\) is the distance vector between molecule \(i\) and each one of its four nearest neighbours \(j\).

The alignment of the orientation of the local structures is measured by the normalized dot product of the local orientational bond order parameter, \(q_l(i)\), between each pair of neighbor molecules \(i\) and \(j\), given by:

\[a(i,j) = \frac{\sum_{m=-l}^{l} q_{lm}(i)\, q_{lm}^*(j)} {|q_l(i)|\,|q_l(j)|}\]

where, \(q_{lm}^*\) is the complex conjugate of \(q_{lm}\).

Depending on the values of \(a(i,j)\), the bond between each pair of molecules \(i\) and \(j\) is classified as being either staggered or eclipsed.

Algorithm

Eclipsed Bonds

Staggered Bonds

CHILL

CHILL+

Since each molecule \(i\) in deeply supercooled water has four nearest neighbours, the type of the resultant four bonds is used to identify phases according to the CHILL and CHILL+ algorithms.

The CHILL algorithm can classify water molecules as belonging to the cubic, hexagonal, interfacial or liquid (amorphous) phase:

Phase | E | S | Neighbours | Description | |———- |– |– |——— |—————————————————— | | Cubic | 0 | 4 | 4 | All 4 neighbours are staggered. | | Hexagonal | 1 | 3 | 4 | 3 staggered and 1 eclipsed bond. | | Interfacial | any | 2 | 4 | Belongs to the ice phase, but does not satisfy strict criteria of cubic and hexagonal | | | 0 | 3 | 4 | | | Liquid | N/A | N/A | any | Classification by exclusion. |

Here, E and S refer to eclipsed and staggered bonds. The neighbours column contains the number of nearest neighbours.

The CHILL+ algorithm, which is a modified version of CHILL, additionally identifies clathrate and interfacial clathrate phases. The criteria are enumerated in the table below:

Phase

E

S

Neighbours

Description

Cubic

0

4

4

no change was made to

staggered bond criterion from CHILL | | Hexagonal | 1 | 3 | 4 | wider eclipsed range compared to CHILL; identifies 99% hexagonal ice up to 270 K | | Interfacial | any | 2 | 4 | must have at least one first neighbor water with more than two staggered bonds | | | 0 | 3 | 4 | must have at least one first neighbor water with more than one staggered bond | | Clathrate | 4 | 0 | 4 | bulk clathrate and part of the interface of clathrates can be identified with four eclipsed bonds | | Interfacial clathrate | 3 | any | 4 | partial clathrate cages and threads of clathrate-like order | | Liquid | N/A | N/A | any | classifies as liquid if none of the above criteria are fulfilled |

Although both the CHILL and CHILL+ classification schemes take into account the local order of the environment of each particle, both schemes output a per-particle classification.

Changelog

Global Functions

std::vector<TemplateHit> classifyTemplates(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, const std::vector<std::vector<int>>& nList, int kNeigh = 12)

Overlay the k nearest neighbours of each particle onto FCC, HCP, BCC and SC shells.

Uses IRA when linked; otherwise Horn on a distance-sorted correspondence (high-symmetry shells only).

std::vector<double> soapSpectrum(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, int iatom, const std::vector<std::vector<int>>& nList, int nMax, int lMax, double rcut)

Bartok SOAP power spectrum for one particle (nMax radial Gaussians, spherical harmonics through lMax).

Length is nMax*nMax*(lMax+1).

std::vector<std::vector<double>> soapSpectrumAll(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, const std::vector<std::vector<int>>& nList, int nMax, int lMax, double rcut)

Bartok SOAP power spectrum for every particle.

Length is nop; each row is nMax*nMax*(lMax+1).

std::vector<double> voronoiFeature(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, int iatom, double candidateCutoff)

[q4, q6, q8] from the Voronoi-weighted Steinhardt path (Mickel).

std::vector<std::vector<double>> voronoiFeatures(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, double candidateCutoff)

[q4, q6, q8] for every particle. One tessellation per order (l = 4, 6, 8).

std::vector<VoronoiWeights> voronoiFacetWeights(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, double candidateCutoff)

Facet neighbours and area weights for every particle, with an a posteriori exactness certificate.

candidateCutoff seeds the bisector search; the cutoff is enlarged automatically until the certificate holds or the growth cap is reached.

Certificate. The cell of particle i clipped against every candidate within cutoff c is the true Voronoi cell whenever every vertex of the clipped cell lies within c/2 of i: a particle k beyond the cutoff has |d_k| > c, its bisector plane sits at distance |d_k|/2 > c/2 from i, and a half-space whose boundary is farther than c/2 cannot cut a region contained in the ball of radius c/2. An open or under-clipped cell keeps vertices on the seeding square at distance ~c and fails the certificate, which is exactly the failure the enlargement retries.

SteinhardtQl steinhardtQlVoronoi(const molSys::PointCloud<molSys::Point<double>, double>& yCloud, double candidateCutoff, int orderL)

Steinhardt parameters of degree orderL (3, 4, 6 or 8) with the bond sum weighted by Voronoi facet areas.

qlBar averages the weighted q_lm over the particle and its facet neighbours, after Lechner and Dellago.