Local crystalline order ¶

Asap contains three methods for analysing the local crystalline order of a simulation: Coordination Numbers, Polyhedral Template Matching and Common Neighbor Analysis.

Unless otherwise noted, these methods work in both serial and parallel simulations. All objects on this page are imported like this:

from asap3.analysis import XXX

Coordination Numbers ¶

The coordination number of an atom is the number of nearest neighbor atoms. In a realistic system, is not necessarily well defined if two atoms are nearest neighbors, so the coordination number is defined as the number of neighbors within a cutoff distance. Radial Distribution Functions can be used to choose a good value for this cutoff (the position of the minimum between the first and second peak in the RDF is a good choice).

Asap provides a function for calculating the coordination numbers.

CoordinationNumbers(atoms, rCut=None):: Returns the coordination numbers of the atoms. The cutoff can be given, otherwise it tries to guess a value between the first and second coordination shell, based on the atomic numbers.

Note that for pure elements, not providing a cutoff distance (rCut) typically works well since the function will get a reasonable value based on the crystalline lattice constant of the material.

Polyhedral Template Matching (PTM)¶

Available from Asap version 3.8.17.

PTM is a new alternative to the popular Common Neigbor Analysis, providing much of the same advantages, but with a greater robustness against thermal vibrations, and does not depend critically on a cutoff.

The PTM classifies the local crystalline order, and identifies local simple cubic (SC), face-centered cubic (FCC), body-centered cubic (FCC), hexagonal closed-packed (HCP) and icosahedral (ICO) order. In addition, some ordered alloys based on the FCC and BCC structures are also detected, namely L1_0, L1_2 and B2 structures.

The PTM also return the local rotation of the lattice (if a crystalline order is found) and optionally the elastic strain tensor.

The method is described by P.M. Larsen, S. Schmidt, and J. Schiøtz, Modelling Simul. Mater. Sci. Eng. 24, 055007 (2016). It is also available as a preprint here.

The PTM version in Asap is dated, a more modern version is available in the project mdapy (the Asap authors are not associated with that project).

PTM(atoms, cutoff=10.0, rmsd_max=None, target_structures=None, calculate_strains=False, quick=False, return_nblist=False):

Run Complex Hull Analysis on an Atoms object.

Parameters:

atoms:

The atoms object

cutoff=10.0:

A cutoff used for the neighborlist. Must be large enough that all nearest neighbors are returned (second-nearest for BCC). Using a too large value may impact performance, but not results.

rmsd_max=None:

If set, matches with a RMSD above this threshold are reclassified as unknown structure (type 0).

This is done as post-processing, and is completely equivalent to:

data['structure'][data['rmsd'] > rmsd_max] = 0

target_structures=None:

A tuple of structures to be investigated. It defaults to (‘sc’, ‘fcc’, ‘hcp’, ‘ico’, ‘bcc’). It MUST be a tuple, not a list or other sequence (bug?).

calculate_strains=False:

Set to True to calculate strains.

quick=False:

Set to True to use Euclidian ordering instead of topological ordering of the neighbor atoms. It speeds up the analysis by approximately a factor 2, but is less robust at high temperatures.

return_nblist=False:

Return the NeighborList object used internally. This may be used for preprocessing. Do not enable unless you need it, as this object can use a lot of memory.

The PTM function returns a dictionary with NumPy arrays as values. Each of these are a NumPy array of the same length as the number of atoms.

Key names and descriptions of returned data. In the description, i is the first index of the array and identifies the atom.

‘structure’:: The local crystal structure around atom i, if any. 0 = none; 1 = FCC; 2 = HCP; 3 = BCC; 4 = Icosahedral; 5 = SC.
‘alloytype’:: The alloy structure identified. 0 = unidentified; 1 = pure element; 2 = L1_0; 3 = L1_2 majority atom; 4 = L1_2 minority atom; 5 = B2. (0 is returned if structure is not 1 or 3; or if no known alloy structure is recognized)
‘rmsd’:: The RMSD error in the fitting to the template, or INF if no structure was identified.
‘scale’:: The average distance to the nearest neighbors for structures 1-4; or the average distance to nearest and next-nearest neighbors for structure 5 (BCC); or INF if no structure was identified.
‘orientation’:: The orientation of the crystal lattice, expressed as a unit quaternion. If no structure was found, the illegal value (0, 0, 0, 0) is returned.
‘strain’: (only present if calculate_strains=True). The strain tensor as a symmetric 3x3 matrix. The trace of the matrix is 1.0, since a hydrostatic component of the strain cannot be determined without a-priori knowledge of the reference lattice parameter. If such knowledge is available, the hydrostatic component of the strain can be calculated from scales[i].
‘info’:: A tuple of two integers, the number of atoms analyzed by PTM and the number of atoms skipped due to insufficient neighbors. Mostly for debugging purposes, and may change without warning.
‘nblist’:: A Neighbor lists object representing a full neighbor list with the cutoff specified when calling PTM. Only returned on request.

Example ¶

Assuming that atoms is an atoms object, you can visualize the different local crystal structures like this:

from asap3.analysis import PTM
from ase.visualize import view

ptmdata = PTM(atoms, rmsd_max=0.2)
atoms.set_tags(ptmdata['structure']
view(atoms)

Note that the rmsd_max argument is not necessary, but it prevents that some surface atoms are misclassified as highly distorted simple cubic structures. A value of 0.2 is high enough to only exclude spurious classifications.

PTM Observer ¶

A wrapper object which does PTM analysis can be attached to the dynamics, in order to automatically do PTM analysis during a simulation. For an example, see Running Polyhedral Template Matching under the Examples tab on these pages.

PTMobserver(atoms, cutoff=10.0, rmsd_max=None, target_structures=None, calculate_strains=False, quantity='structure', analyze_first=True):

Per default, sets the tags of the atoms according to their local crystal structure. 0 = none; 1 = FCC; 2 = HCP; 3 = BCC; 4 = Icosahedral; 5 = SC.

An optional argument \(quantity\) can be set to ‘alloytype’ or None, then the tags are set to the alloytype calculated by the PTM function (see that function’s documentation); or not set if None is given.

This class is intended to use as an observer, so its analyze() function is called automatically by e.g. the dynamics. It can itself act as a subject, so a Plotter or a Trajectory can be called just after the calculations.

It has a method get_data which returns the full dictionary of PTM data from the last calculation.

Parameters:

\(atoms\):

The atoms object

\(cutoff\):

A cutoff used for the neighborlist. Must be large enough that all nearest neighbors are returned (second-nearest for BCC). Using a too large value may impact performance, but not results.

\(rmsd_max=None\):

If set, matches with a RMSD above this threshold are reclassified as unknown structure (type 0).

This is done as post-processing, and is completely equivalent to:

data['structure'][data['rmsd'] > rmsd_max] = 0

\(target_structures=None\):

A tuple of structures to be investigated. It defaults to (‘sc’, ‘fcc’, ‘hcp’, ‘ico’, ‘bcc’). It MUST be a tuple, not a list or other sequence (bug?).

\(calculate_strains=False\):

Set to True to calculate strains.

\(quantity='structure'\):

The quantity used to set the tags.

\(analyze_first=True\):

Should an analysis be made as soon as this object is created? Leave as true if you plan to attach a Trajectory as an observer to this object, and if the Trajectory will save the initial state.

\(quick=False\):

Set to True to use Euclidian ordering instead of topological ordering of the neighbor atoms. It speeds up the analysis by approximately a factor 2, but is less robust at high temperatures.

PTMdislocations ¶

PTMdislocations is an analysis module build on top of PTM. It identifies stacking faults and dislocation lines in systems with an FCC crystal structure. It is used exactly like PTM Observer above, but identifies the following structures:

0 = none; 1 = FCC; 2 = HCP; 3 = BCC; 4 = Icosahedral; 5 = SC; 6 = Twin boundary; 7 = Stacking Fault; 8 = Dislocation line.

It works by looking at the neighbors to atoms identified as HCP.

If an HCP atom has 6 HCP neighbors centered around the atom then it is assumed that it is part of a single (111) layer of HCP atoms, i.e. a twin boundary.
If an HCP atom has 9 HCP neighbors then it is assumed that it is part of a double (111) layer of HCP atoms, i.e. an intrinsic stacking fault.
If an HCP atom has between 5 and 8 neighbors, and they are not evenly distributed around it, then it is likely to the edge of a stacking fault, i.e. the core of a partial dislocation.

Common Neighbor Analysis (CNA)¶

Common Neighbor Analysis (CNA) is a way to get information about the local crystalline order of a system. It can for example be used to find stacking faults and dislocations in a simulation of a metal. CNA works by examining all bonds (nearest-neighbor pairs of atoms), each such bond is classified according to how the common neighbors of the two atoms are located, and finally all bonds related to an atom is used to classify the crystal structure of the atom. Asap provides two impementations of CNA:

Restricted CNA which classifies atoms in one of three groups, FCC, HCP and Other according to the local crystal structure but skips doing a full CNA analysis
Full CNA, which does a full CNA analysis, and can return the data in several ways.

RestrictedCNA(atoms, rCut=None, analyze_first=True):

Creates an object for doing restricted CNA analysis. This object can be used as-is or to create an observer object. The object will store a reference to the atoms, and use it to calculate the CNA when the analyze method is called.

It requires a value of rCut between the first and second shell. It tries to guess a value for rCut from the atomic numbers if none is specified.

If the optional parameter analyze_first is true (the default), a CNA analysis is run already when the object is created, instead of waiting for the first call of analyze().

The result of the restricted CNA analysis is returned by setting the tags of the atoms. FCC atoms are given tag 0, HCP atoms 1, and everything else 2.

FullCNA(atoms, rCut=None, normal=False, per_z=False, total=False, total_z=False, raw=False):

Creates an object for doing full CNA analysis. This object can be used as-is or to create an observer object. The object will store a reference to the atoms, and use it to calculate the CNA when one of the get methods is called.

update(atoms=None, cutoff=None):: Update the stored reference to the atoms and/or the cutoff distance.
get_normal_cna():: CNA is returned as a dictionary for each atom, with CNA 3-tuples as keys and number of bonds as values.
get_total_cna():: A single dictionary with CNA 3-tuples as keys, and total number of bonds as values. In a parallel simulation, data for the entire system is returned on all nodes. This function is not supported in parallel simulations.
get_total_perz_cna():: A single dictionary with pairs of atomic numbers as keys and dictionaries for values. These dictionaries contain the same data as for the “total” method, but only refers to atoms with the two atomic numbers specified. The pairs of atomic numbers (Z1, Z2) are always chosen so Z1 <= Z2. This function is not supported in parallel simulations.
get_normal_and_total_cna():: Returns a tuple with two elements, the first is the same as get_normal_cna() would return, the second is what get_total_cna() would return. However, unlike calling the two methods separately, the actual CNA calculation is only done once. This function is not supported in parallel simulations.
get_raw_cna():: Raw CNA data as a NumPy array of shape (X, 3) and type int32, where X is the number of bond. Each row consists of the indices of the two atoms involved in the bond, and the CNA encoded as 65536*i1 + 256*i2 + i3 for the CNA tuple (i1,i2,i3). Mostly for internal use.
get_cna_tuple(cnaint):: Translate CNA encoded as an integer (as described under get_raw_cna() to CNA as a 3-tuple. If called multiple times with the same cnaint, it returns references to the same tuple. Mostly for internal use..

Implementation details: The CNA calculation resulting in the raw CNA data is implemented in C++. The access methods may be implemented in Python, but some of them are also implemented in C++ for performance. The object will reuse (and possibly update) the neighbor list of any associated calculator, if the neighbor list is usable, or create its own if not.

Local crystalline order¶

Coordination Numbers¶