Main module of pyscal. This module contains definitions of the two major
classes in pyscal - the System and Atom.
Atom is a pure pybind11 class whereas System is a hybrid class with additional
python definitions. For the ease of use, Atom class should be imported from the core
module. The original pybind11 definitions of Atom and System can be found in catom
and csystem respectively.
A System consists of two
major components - the simulation box and the atoms. All the associated variables
are then calculated using this class.
Note
atoms can be accessed or set as atoms. However, due to
technical reasons individual atoms should be accessed using the
get_atom() method. An atom can be assigned
to the atom using the set_atom() method.
Calculate the angular criteria for each atom
:param None:
Return type:
None
Notes
Calculates the angular criteria for each atom as defined in [1]_. Angular criteria is
useful for identification of diamond cubic structures. Angular criteria is defined by,
where cos(theta) is the angle size suspended by each pair of neighbors of the central
atom. A will have a value close to 0 for structures if the angles are close to 109 degrees.
The calculated A parameter for each atom is stored in angular.
nmax (int, optional) – number of neighbors to be considered for centrosymmetry
parameters. Has to be a positive, even integer. Default 12
Return type:
None
Notes
Calculate the centrosymmetry parameter for each atom which can be accessed by
centrosymmetry attribute. It calculates the degree of inversion
symmetry of an atomic environment. Centrosymmetry recalculates the neighbor using
the number method as specified in ¬pyscal.core.System.find_neighbors() method. This
is the ensure that the required number of neighbors are found for calculation of the parameter.
The Greedy Edge Selection (GES) [1] as specified in [2] in used in this method.
GES algorithm is implemented in LAMMPS and Ovito. Please see [2] for
a detailed description of the algorithms.
angles (bool, optional) – If True, return the list of cosines of all neighbor pairs
Returns:
angles – list of all cosine values, returned only if angles is True.
Return type:
array of floats
Notes
This method tries to distinguish between crystal structures by finding the cosines of angles
formed by an atom with its neighbors. These cosines are then historgrammed with bins
[-1.0, -0.945, -0.915, -0.755, -0.705, -0.195, 0.195, 0.245, 0.795, 1.0] to find a vector for
each atom that is indicative of its local coordination. Compared to chi parameters from chi_0 to
chi_7 in the associated publication, the vector here is from chi_0 to chi_8. This is due to an additional
chi parameter which measures the number of neighbors between cosines -0.705 to -0.195.
Parameter nlimit specifies the number of nearest neighbors to be included in the analysis to find the cutoff.
If parameter angles is true, an array of all cosine values is returned. The publication further provides
combinations of chi parameters for structural identification which is not implemented here. The calculated
chi params can be accessed using chiparams.
lattice_constant (float, optional) – lattice constant to calculate CNA. If not specified,
adaptive CNA will be used
Returns:
cna – dict containing the cna signature of the system
Return type:
dict
Notes
Performs the common neighbor analysis [1][2] or the adaptive common neighbor
analysis [2] and assigns a structure to each atom.
If lattice_constant is specified, a convential common neighbor analysis is
used. If lattice_constant is not specified, adaptive common neighbor analysis is used.
The assigned structures can be accessed by structure.
The values assigned for stucture are 0 Unknown, 1 fcc, 2 hcp, 3 bcc, 4 icosahedral.
where .. math:: S_{ij} = sum_{m=-6}^6 q_{6m}(i) q_{6m}^*(i)
The keyword averaged is True, the disorder value is averaged over the atom and its neighbors. The disorder value
can be accessed using disorder and the averaged version can be accessed using
avg_disorder. For ordered systems, the value of disorder would be zero which would increase
and reach one for disordered systems.
If averaged is True, the energy is averaged over the neighbors of an
atom. If neighbors were calculated before calling this method, those neighbors
are used for averaging. Otherwise neighbors are calculated on the fly
with an adaptive cutoff method.
Calculate pairwise multicomponent short range order
Parameters:
reference_type (int, optional) – type of the atom to be used a reference. default 1
compare_type (int, optional) – type of the atom to be used to compare. default 2
average (bool, optional) – if True, average over all atoms of the reference type in the system.
default True.
delta (bool, optional) – if True, SRO calculation contain the Kronecker delta in the definition.
if False, delta is always 0, and the function return the Cowley-SRO value.
default True.
Returns:
vec – The short range order averaged over the whole system for atom of
the reference type. Only returned if average is True. First value is SRO
of the first neighbor shell and the second value corresponds to the second
nearest neighbor shell.
Return type:
list of float
Notes
Calculates the pairwise multicomponent short range order for a higher-dimensional systems alloy using the approach by
Fontaine [1]. Pairwise multicomponent short range order is calculated as,
where i refers to reference type. n_j is the number of atoms of the non reference type among the c_j atoms
in the ith shell. m_A is the concentration of the non reference atom. delta_{ij}” = 1 if i = j. Please
note that the value is calculated for shells 1 and 2 by default. In order for
this to be possible, neighbors have to be found first using the find_neighbors()
method. The selected neighbor method should include the second shell as well. For this
purpose method=cutoff can be chosen with a cutoff long enough to include the second
shell. In order to estimate this cutoff, one can use the calculate_rdf()
method.
q_l (int or list of ints) – A list of all Steinhardt parameters to be found from 2-12.
averaged (bool, optional) – If True, return the averaged q values, default False
only_averaged (bool, optional) – If True, only calculate the averaged part. default False
condition (callable or atom property) – Either function which should take an Atom object, and give a True/False output
or an attribute of atom class which has value or 1 or 0.
clear_condition (bool, optional) – clear the condition variable for all atoms
Return type:
None
Notes
Enables calculation of the Steinhardt parameters [1] q from 2-12. The type of
q values depend on the method used to calculate neighbors. See the description
find_neighbors() for more details. If the keyword average is set to True,
the averaged versions of the bond order parameter [2] is returned. If only the averaged
versions need to be calculated, only_averaged keyword can be set to False.
The neighbors over which the q values are calculated can also be filtered. This is done
through the argument condition which is passed as a parameter.
condition can be of two types. The first type is a function which takes an
Atom object and should give a True/False value. condition can also be an
Atom attribute or a value from custom values stored in an atom. See
cluster_atoms() for more details. If the
condition is equal for both host atom and the neighbor, the neighbor is considered for
calculation of q parameters. This is slightly different from cluster_atoms()
where the condition has to be True for both atoms. condition is only cleared when neighbors are
recalculated. Additionally, the keyword clear_condition can also be used to clear the condition
and reset it to 0. By default, condition is applied to both unaveraged and averaged q parameter
calculation. If condition is needed for only averaged q parameters, this function can be called
twice, initially without condition and averaged=False, and then with a condition specified
and averaged=True. This way, the condition will only be applied to the averaged q calculation.
reference_type (int, optional) – type of the atom to be used a reference. default 1
average (bool, optional) – if True, average over all atoms of the reference type in the system.
default True.
Returns:
vec – The short range order averaged over the whole system for atom of
the reference type. Only returned if average is True. First value is SRO
of the first neighbor shell and the second value corresponds to the second
nearest neighbor shell.
Return type:
list of float
Notes
Calculates the short range order for an AB alloy using the approach by
Cowley [1]. Short range order is calculated as,
\[\alpha_i = 1 - \frac{n_i}{m_A c_i}\]
where n_i is the number of atoms of the non reference type among the c_i atoms
in the ith shell. m_A is the concentration of the non reference atom. Please
note that the value is calculated for shells 1 and 2 by default. In order for
this to be possible, neighbors have to be found first using the find_neighbors()
method. The selected neighbor method should include the second shell as well. For this
purpose method=cutoff can be chosen with a cutoff long enough to include the second
shell. In order to estimate this cutoff, one can use the calculate_rdf()
method.
edge_cutoff (float, optional) – cutoff for edge length. Default 0.05.
area_cutofffloat, optional
cutoff for face area. Default 0.01.
edge_lengthbool, optional
if True, a list of unrefined edge lengths are returned. Default false.
Returns:
vorovector – array of the form (n3, n4, n5, n6)
Return type:
array like, int
Notes
Returns a vector of the form (n3, n4, n5, n6), where n3 is the number
of faces with 3 vertices, n4 is the number of faces with 4
vertices and so on. This can be used to identify structures [1] [2].
The keywords edge_cutoff and area_cutoff can be used to tune the values to minimise
the effect of thermal distortions. Edges are only considered in the analysis if the
edge_length/sum(edge_lengths) is at least edge_cutoff. Similarly, faces are only
considered in the analysis if the face_area/sum(face_areas) is at least face_cutoff.
condition (callable or atom property) – Either function which should take an Atom object, and give a True/False output
or an attribute of atom class which has value or 1 or 0.
largest (bool, optional) – If True returns the size of the largest cluster. Default False.
cutoff (float, optional) – If specified, use this cutoff for calculation of clusters. By default uses the cutoff
used for neighbor calculation.
Returns:
lc – Size of the largest cluster. Returned only if largest is True.
Return type:
int
Notes
This function helps to cluster atoms based on a defined property. This property
is defined by the user through the argument condition which is passed as a parameter.
condition can be of two types. The first type is a function which takes an
Atom object and should give a True/False value. condition can also be an
Atom attribute or a value from custom values stored in an atom.
When clustering, the code loops over each atom and its neighbors. If the
condition is true for both host atom and the neighbor, they are assigned to
the same cluster. For example, a condition to cluster solid atoms would be,
defcondition(atom):#if both atom is solidreturn(atom1.solid)
The same can be done by passing “solid” as the condition argument instead of the above
function. Passing a function allows to evaluate complex conditions, but is slower than
passing an attribute.
This method finds in the underlying fcc/hcp lattice in diamond. It works
by the method described in this publication .
For each atom, 4 atoms closest to it are identified. The neighbors of the its neighbors
are further identified and the common neighbors shared with the host atom are selected.
These atom will fall in the underlying fcc lattice for cubic diamond or hcp lattice
for hexagonal lattice.
If neighbors are previously calculated, they are reset when this method is used.
method ({'cutoff', 'voronoi', 'number'}) – cutoff method finds neighbors of an atom within a specified or adaptive cutoff distance from the atom.
voronoi method finds atoms that share a Voronoi polyhedra face with the atom. Default, cutoff
number method finds a specified number of closest neighbors to the given atom. Number only populates
cutoff{ float, ‘sann’, ‘adaptive’}
the cutoff distance to be used for the cutoff based neighbor calculation method described above.
If the value is specified as 0 or adaptive, adaptive method is used.
If the value is specified as sann, sann algorithm is used.
thresholdfloat, optional
only used if cutoff=adaptive. A threshold which is used as safe limit for calculation of cutoff.
filter{‘None’, ‘type’, ‘type_r’}, optional
apply a filter to nearest neighbor calculation. If the filter keyword is set to
type, only atoms of the same type would be included in the neighbor calculations.
If type_r, only atoms of a different type will be included in the calculation. Default None.
voroexpint, optional
only used if method=voronoi. Power of the neighbor weight used to weight the contribution of each atom towards
Steinhardt parameter values. Default 1.
paddingdouble, optional
only used if cutoff=adaptive or cutoff=number. A safe padding value used after an adaptive cutoff is found. Default 1.2.
nlimitint, optional
only used if cutoff=adaptive. The number of particles to be considered for the calculation of adaptive cutoff.
Default 6.
nmaxint, optional
only used if cutoff=number. The number of closest neighbors to be found for each atom. Default 12
Return type:
None
Raises:
RuntimeWarning – raised when threshold value is too low. A low threshold value will lead to ‘sann’ algorithm not converging
when finding a neighbor. This function will try to automatically increase threshold and check again.
RuntimeError – raised when neighbor search was unsuccessful. This is due to a low threshold value.
Notes
This function calculates the neighbors of each particle. There are several ways to do this. A complete description of
the methods can be found here.
Method cutoff and specifying a cutoff radius uses the traditional approach being the one in which the neighbors of an atom
are the ones that lie in the cutoff distance around it.
In order to reduce time during the distance sorting during the adaptive methods, pyscal sets an initial guess for a cutoff distance.
This is calculated as,
threshold is a safe multiplier used for the guess value and can be set using the threshold keyword.
In Method cutoff, if cutoff='adaptive', an adaptive cutoff is found during runtime for each atom [1].
Setting the cutoff radius to 0 also uses this algorithm. The cutoff for an atom i is found using,
padding is a safe multiplier to the cutoff distance that can be set through the keyword padding. nlimit keyword sets the
limit for the top nlimit atoms to be taken into account to calculate the cutoff radius.
In Method cutoff, if cutoff='sann', sann algorithm is used [2]. There are no parameters to tune sann algorithm.
The second approach is using Voronoi polyhedra which also assigns a weight to each neighbor in the ratio of the face area between the two atoms.
Higher powers of this weight can also be used [3]. The keyword voroexp
can be used to set this weight.
If method os number, instead of using a cutoff value for finding neighbors, a specified number of closest atoms are
found. This number can be set through the argument nmax.
Warning
Adaptive and number cutoff uses a padding over the intial guessed “neighbor distance”. By default it is 2. In case
of a warning that threshold is inadequate, this parameter should be further increased. High/low value
of this parameter will correspond to the time taken for finding neighbors.
bonds (int or float, optional) – Minimum number of solid bonds for an atom to be identified as
a solid if the value is an integer. Minimum fraction of neighbors
of an atom that should be solid for an atom to be solid if the
value is float between 0-1. Default 0.5.
threshold (double, optional) – Solid bond cutoff value. Default 0.5.
avgthreshold (double, optional) – Value required for Averaged solid bond cutoff for an atom to be identified
as solid. Default 0.6.
cluster (bool, optional) – If True, cluster the solid atoms and return the number of atoms in the largest
cluster.
q (int, optional) – The Steinhardt parameter value over which the bonds have to be calculated.
Default 6.
cutoff (double, optional) – Separate value used for cluster classification. If not specified, cutoff used
for finding neighbors is used.
right (bool, optional) – If true, greater than comparison is to be used for finding solid particles.
default True.
Returns:
solid – Size of the largest solid cluster. Returned only if cluster=True.
Return type:
int
Notes
The neighbors should be calculated before running this function.
Check find_neighbors() method.
bonds define the number of solid bonds of an atom to be identified as solid.
Two particles are said to be ‘bonded’ if [1],
where threshold values is also an optional parameter.
If the value of bonds is a fraction between 0 and 1, at least that much of an atom’s neighbors
should be solid for the atom to be solid.
An additional parameter avgthreshold is an additional parameter to improve solid-liquid distinction.
In addition to having a the specified number of bonds,
\[\langle s_{ij} \rangle > avgthreshold\]
also needs to be satisfied. In case another q value has to be used for calculation of S_ij, it can be
set used the q attribute. In the above formulations, > comparison for threshold and avgthreshold
can be changed to < by setting the keyword right to False.
If cluster is True, a clustering is done for all solid particles. See find_clusters()
for more details.
find_neighbors (bool, optional) – If True, find 4 closest neighbors
Returns:
diamondstructure – dict of structure signature
Return type:
dict
Notes
Identify diamond structure using the algorithm mentioned in [1]. It is an
extended CNA method. The integers 5, 6, 7, 8, 9 and 10 are assigned to the
structure variable of the atom. 5 stands for cubic diamond, 6 stands for first
nearest neighbors of cubic diamond and 7 stands for second nearest neighbors
of cubic diamond. 8 signifies hexagonal diamond, the first nearest neighbors
are marked with 9 and second nearest neighbors with 10.
Read input file that contains the information of system configuration.
Parameters:
filename (string) – name of the input file.
format ({'lammps-dump', 'poscar', 'ase', 'mdtraj'}) – format of the input file, in case of ase the ASE Atoms object
compressed (bool, optional) – If True, force to read a gz compressed format, default False.
customkeys (list) – A list containing names of headers of extra data that needs to be read in from the
input file.
Return type:
None
Notes
format keyword specifies the format of the input file. Currently only
a lammps-dump and poscar files are supported. Additionaly, the widely
use Atomic Simulation environment (https://wiki.fysik.dtu.dk/ase/ase/ase.html).
mdtraj objects (http://mdtraj.org/1.9.3/) are also supported by using the keyword
‘mdtraj’ for format. Please note that triclinic boxes are not yet supported for
mdtraj format.
Atoms object can also be used directly. This function uses the
traj_process() module to process a file which is then assigned to system.
compressed keyword is not required if a file ends with .gz extension, it is
automatically treated as a compressed file.
Triclinic simulation boxes can also be read in.
If custom_keys are provided, this extra information is read in from input files if
available. This information is not passed to the C++ instance of atom, and is stored
as a dictionary. It can be accessed directly as atom.custom[‘customval’]
For example, an Atom at location i in the list of all atoms in
System can be queried by,
atom=System.get_atom(i), then any kind of modification, for example, the
position of the Atom can done by, atom.pos=[2.3,4.5,4.5]. After
modification, the Atom can be set back to its position in System by
set_atom().
Although the complete list of atoms can be accessed or set using atoms=sys.atoms,
get_atom and set_atom functions should be used for accessing individual atoms.
If an atom already exists at that index in the list, it will be overwritten and will
lead to loss of information.
factor (float, optional) – factor for multiplication of cutoff value.
default 1
Return type:
None
Notes
Assign cutoffs for each atom based on the nearest
neighbor distance. The cutoff assigned is the average nearest
neighbor distance multiplied by factor.
colorby (string, optional) – property over which the atoms are to be colored. It can be any
attributed of Atom, a custom attribute, or calculated q values which can be accessed
as qx or aqx where x stands for the q number.
filterby (string, optional) – property over which the atoms are to be filtered before plotting.
It can be any attribute of atom, or a custom value of atom. It should provide
a True or False value.
format (string, {'lammps-dump', 'lammps-data', 'poscar'}) – format of the output file, default lammps-dump
Currently only lammps-dump format is supported.
customkeys (list of strings, optional) – a list of extra atom wise values to be written in the output file.
customvals (list or list of lists, optional) – If customkey is specified, customvals take an array of the same length
as number of atoms, which contains the values to be written out.
compressed (bool, optional) – If true, the output is written as a compressed file.
timestep (int, optional) – timestep to be written to file. default 0
species (None, optional) – species of the atoms. Required if any format other than ‘lammps-dump’ is used. Required
for convertion to ase object.
Return type:
None
Notes
to_file method can handle a number of file formats. The most customizable format is the
lammps-dump which can take a custom options using customkeys and customvals. customkeys
will be the header written to the dump file. It can be any Atom attribute, any property
stored in custom variable of the Atom, or calculated q values which can be given by q4,
aq4 etc. External values can also be provided using customvals option. customvals array
should be of the same length as the number of atoms in the system.
For all other formats, ASE is used to write out the file, and hence the species keyword
needs to be specified. If initially, an ASE object was used to create the System, species
keyword will already be saved, and need not be specified. In other cases, species should
be a list of atomic species in the System. For example [“Cu”] or [“Cu”, “Al”], depending
on the number of species in the System. In the above case, atoms of type 1 will be mapped to
Cu and of type 2 will be mapped to Al. For a complete list of formats that ASE can handle,
see here .
pos (list of floats of length 3) – position of the Atom, default [0,0,0]
id (int) – id of the Atom, default 0
type (int) – type of the Atom, default 1
Notes
A pybind11 class for holding the properties of a single atom. Various properties of the atom
can be accessed through the attributes and member functions which are described below in detail. Atoms can
be created individually or directly by reading a file. Check the examples for more
details on how atoms are created. For creating atoms directly from an input file check
the documentation of System class.
Although an Atom object can be created independently, Atom should be thought of
inherently as members of the System class. All the properties that define an atom are
relative to the parent class. System has a list of all atoms. All the properties of an
atom, hence should be calculated through System.
Examples
>>> #method 1 - individually>>> atom=Atom()>>> #now set positions of the atoms>>> atom.pos=[23.0,45.2,34.2]>>> #now set id>>> atom.id=23>>> #now set type>>> atom.type=1>>> #Setting through constructor>>> atom=Atom([23.0,45.2,34.2],23,1)
Float.
The value of angular parameter A of an atom. The angular parameter measures the tetrahedral coordination of an atom.
Meaningful values are only returned if the property is calculated using calculate_angularcriteria().
float. Averaged version of the Voronoi volume which is calculated as an average over
itself and its neighbors. Only calculated when the find_neighbors()
using the method=’voronoi’ option is used.
Float.
The value of chiparameter of an atom. The return value is a vector of length 8.
Meaningful values are only returned if chi params are calculated using calculate_chiparams().
int.
condition that specifies if an atom is included in the clustering algorithm or not.
Only atoms with the value of condition=1 will be used for clustering in
cluster_atoms().
dict.
dictionary specfying custom values for an atom. The module only stores the id, type and
position of the atom. If any extra values need to be stored, they can be stored in custom
using atom.custom = {“velocity”:12}. read_inputfile() can also
read in extra atom information. By default, custom values are treated as string.
list of floats. For each face, this vector contains the lengths of edges
that make up the Voronoi polyhedra of the atom. Only calculated when the find_neighbors()
using the method=’voronoi’ option is used.
list of floats. List consisting of the perimeters of each Voronoi face of an atom.
Only calculated when the find_neighbors()
using the method=’voronoi’ option is used.
list of floats. A list of the number of vertices shared between an atom and its
neighbors. Only calculated when the find_neighbors()
using the method=’voronoi’ option is used.
List of floats.
Used to weight the contribution of each neighbor atom towards the value of
Steinhardt’s parameters. By default, each atom has a weight of 1 each. However,
if find_neighbors() is used with method=’voronoi’,
each neighbor gets a weight proportional to the area shared between the neighboring
atom and host atom.
list of floats. For each Voronoi face of the atom, this values includes a List
of vertices that constitute the face. Only calculated when the find_neighbors()
using the method=’voronoi’ option is used.
list of floats. A list of positions of each vertex of the Voronoi polyhedra of
the atom. Only calculated when the find_neighbors()
using the method=’voronoi’ option is used.
float. Voronoi volume of the atom. The Voronoi volume is only calculated if neighbors
are found using the find_neighbors() using the method=’voronoi’
option.
list of ints. A vector of the form (n3, n4, n5, n6) where n3 is the number of faces with 3 vertices,
n4 is the number of faces with 4 vertices and so on. This can be used to identify structures [1][2].
Vorovector is calculated if the calculate_vorovector() method is used.
Create a basic crystal structure and return it as a list of Atom objects
and box dimensions.
Parameters:
structure ({'sc', 'bcc', 'fcc', 'hcp', 'diamond', 'a15' or 'l12'}) – type of the crystal structure
lattice_constant (float, optional) – lattice constant of the crystal structure, default 1
repetitions (list of ints of len 3, optional) – of type [nx, ny, nz], repetions of the unit cell in x, y and z directions.
default [1, 1, 1].
ca_ratio (float, optional) – ratio of c/a for hcp structures, default 1.633
noise (float, optional) – If provided add normally distributed noise with standard deviation noise to the atomic positions.
Returns:
atoms (list of Atom objects) – list of all atoms as created by user input
box (list of list of floats) – list of the type [[xlow, xhigh], [ylow, yhigh], [zlow, zhigh]] where each of them are the lower
and upper limits of the simulation box in x, y and z directions respectively.
Load the data of a block into memory as a dictionary
of numpy arrays
Parameters:
blockno (int) – number of the block to be read, starts from 0
Return type:
None
Notes
When the data of a block is loaded, it is accessible
through Trajectory.data[x]. This data can then be
modified. When the block is written out, the modified
data is written instead of existing one. But, loaded
data is kept in memory until unloaded using unload
method.
pyscal module containing methods for processing of a trajectory. Methods for
reading of input files formats, writing of output files etc are provided in
this module.
Write the state of the system to a trajectory file.
Parameters:
sys (System object) – the system object to be written out
outfile (string) – name of the output file
format (string, optional) – format of the output file
compressed (bool, default false) – write a .gz format
customkey (string or list of strings, optional) – If specified, it adds this custom column to the dump file. Default None.
customvals (list or list of lists, optional) – If customkey is specified, customvals take an array of the same length
as number of atoms, which contains the values to be written out.
timestep (int, optional) – Specify the timestep value, default 0
species (None, optional) – species of the atoms. Required if any format other than ‘lammps-dump’ is used. Required
for convertion to ase object.