pyscal Reference

pyscal.core module

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.

class pyscal.core.System(filename=None, format='lammps-dump', compressed=False, customkeys=None)

Bases: object

Python class for holding the properties of an atomic configuration

__init__(filename=None, format='lammps-dump', compressed=False, customkeys=None)

add_atoms(atoms)

Cleanly add a given list of atoms

Parameters:

atoms (dict) –

Return type:

None

apply_mask(mask_type='primary', ids=None, indices=None, condition=None, selection=False)

Notes

Masks can be used to exclude atoms from neighbor calculations. An atom for which mask is set to True is excluded from the calculation. There are two types of masks, primary or secondary. For example, neighbors are being calculated for a central atom i. The neighbor atom is denoted as j. If primary mask of i is True, no neighbor calculation is carried out for i. If it is False, i is considered. Now if secondary mask of j is True, it will not included in the list of neighbors of i even if it is within the cutoff distance. The primary mask of j has no effect in this situation.

An example situation can be to calculate the local concentration around Ni atoms in a NiAl structure. In this case, the primary mask of all Al atoms can be set to True so that only Ni atoms are considered. Now, in a second case, the task is to count the number of Al atoms around each Ni atom. For this case, the primary mask of all Al atoms can be set to True, and the secondary mask of all Ni atoms can be set to True.

The masks for ghost atoms are copied from the corresponding mask for real atoms.

apply_selection(ids=None, indices=None, condition=None)

property atoms

average_over_neighbors(key, include_self=True)

Perform a simple average over neighbor atoms

Parameters:

• key (string) – atom property

• include_self (bool, optional) – If True, include the host atom in the calculation

property box

Wrap for inbuilt box

property box_dimensions

calculate_angularcriteria()

Calculate the angular criteria for each atom

Parameters:

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,

\[A = \sum_{i=1}^6 (\cos(\theta_i) + \frac{1}{3})^2\]

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 can be accessed by system.angular

References

[1]

Uttormark, MJ, Thompson, MO, Clancy, P, Phys. Rev. B 47, 1993

calculate_chiparams(angles=False)

Calculate the chi param vector for each atom

Parameters:

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.

References

[1]

Ackland, Jones, Phys. Rev. B 73, 2006

calculate_disorder(averaged=False, q=6)

Calculate the disorder criteria for each atom

Parameters:

• averaged (bool, optional) – If True, calculate the averaged disorder. Default False.

• q (int, optional) – The Steinhardt parameter value over which the bonds have to be calculated. Default 6.

Return type:

None

Notes

Calculate the disorder criteria as introduced in [1]. The disorder criteria value for each atom is defined by, .. math:

D_j = \frac{1}{N_b^j} \sum_{i=1}^{N_b} [ S_{jj} + S_{kk} -2S_{jk}]

where .. math:: S_{ij} = sum_{m=-6}^6 q_{6m}(i) q_{6m}^*(i)

Any q value other than six can also be used. This can be specified using the q argument.

The keyword averaged is True, the disorder value is averaged over the atom and its neighbors. For ordered systems, the value of disorder would be zero which would increase and reach one for disordered systems.

This function creates two new attributes for this class: disorder and avg_disorder.

References

[1]

Kawasaki, T, Onuki, A, J. Chem. Phys. 135, 2011

calculate_q(q, averaged=False, continuous_algorithm=False)

Find the Steinhardt parameter q_l for all atoms.

Parameters:

• q (int or list of ints) – A list of all Steinhardt parameters to be found.

• averaged (bool, optional) – If True, return the averaged q values, default False

• continuous_algorithm (bool, optional) – See Notes for description.

Returns:

q – calculated q values

Return type:

list of floats

Notes

Enables calculation of the Steinhardt parameters [1] q. The type of q values depend on the method used to calculate neighbors. See the description find_neighbors() for more details.

The option continuous_algorithm specifies which algorithm to use for calculations. If False, an algorithm [3] is used. The C++ algorithm is faster is a large, consecutive number of q values (> 200) are to be calculated.

This function creates three new attributes for this class: qx, qx_real and qx_imag, where stands for the q number.

References

[1]

Steinhardt, PJ, Nelson, DR, Ronchetti, M. Phys Rev B 28, 1983

[2]

Lechner, W, Dellago, C, J Chem Phys, 2013

calculate_rdf(rmin=0, rmax=5.0, bins=100)

Calculate the radial distribution function.

Parameters:

• rmin (float, optional) – minimum value of the distance histogram. Default 0.0.

• rmax (float, optional) – maximum value of the distance histogram. Default 5.0

• bins (int) – number of bins in the histogram

Returns:

• rdf (array of ints) – Radial distribution function

• r (array of floats) – radius in distance units

cluster_atoms(condition, largest=True, cutoff=0)

Cluster atoms based on a property

Parameters:

• 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 should be a boolean array the same length as number of atoms in the system.

property composition

property concentration

delete(ids=None, indices=None, condition=None, selection=False)

property direct_coordinates

embed_in_cubic_box(input_box=None, return_box=False)

Embedded the triclinic box in a cubic box

find_largest_cluster()

Find largest cluster among given clusters

Parameters:

None –

Returns:

lc – Size of the largest cluster.

Return type:

int

find_neighbors(method='cutoff', cutoff=0, threshold=2, voroexp=1, padding=1.2, nlimit=6, cells=None, nmax=12, assign_neighbor=True)

Find neighbors of all atoms in the System.

Parameters:

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.

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.

cellsbool, optional

If True, always use cell lists. Default None.

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 thefind_neighbors adaptive methods, pyscal sets an initial guess for a cutoff distance. This is calculated as,

\[r_{initial} = threshold * (simulation~box~volume/ number~of~particles)^{(1/3)}\]

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,

\[r_c(i) = padding * ((1/nlimit) * \sum_{j=1}^{nlimit}(r_{ij}))\]

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 is 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.

If cells is None, cell lists are used if number of atoms are higher than 2500. If True, cell lists are always used.

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.

References

[1]

Stukowski, A, Model Simul Mater SC 20, 2012

[2]

van Meel, JA, Filion, L, Valeriani, C, Frenkel, D, J Chem Phys 234107, 2012

[3]

Haeberle, J, Sperl, M, Born, P, arxiv 2019

find_solids(bonds=0.5, threshold=0.5, avgthreshold=0.6, cluster=True, q=6, cutoff=0, right=True)

Distinguish solid and liquid atoms in the system. :param bonds: 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.

Parameters:

• 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], .. math:: s_{ij} = sum_{m=-6}^6 q_{6m}(i) q_{6m}^*(i) geq threshold 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.

References

[1]

Auer, S, Frenkel, D. Adv Polym Sci 173, 2005

classmethod from_structure(structure, lattice_constant=1.0, repetitions=None, ca_ratio=1.633, noise=0, element=None, chemical_symbol=None)

get_concentration()

Return a dict containing the concentration of the system

Parameters:

None –

Returns:

condict – dict of concentration values

Return type:

dict

get_distance(pos1, pos2, vector=False)

Get the distance between two atoms.

Parameters:

• pos1 (list) – first atom position

• pos2 (list) – second atom position

• vector (bool, optional) – If True, return the vector between two atoms

Returns:

distance – distance between the first and second atom.

Return type:

double

Notes

Periodic boundary conditions are assumed by default.

iter_atoms()

property natoms

read_inputfile(filename, format='lammps-dump', compressed=False, customkeys=None)

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 can be accessed directly as self.atoms[‘customkey’]

remap_atoms_into_box()

Go through atoms in the list and remap them into the bix

remove_mask(mask_type='primary', ids=None, indices=None, condition=None, selection=False)

Remove applied masks

Parameters:

mask_type (string, optional) – type of mask to be applied, either primary, secondary or all

Return type:

None

remove_selection(ids=None, indices=None, condition=None)

repeat(repetitions, atoms=None, ghost=False, scale_box=True, assign=False, return_atoms=False)

reset_neighbors()

Reset the neighbors of all atoms in the system.

Parameters:

None –

Return type:

None

Notes

It is used automatically when neighbors are recalculated.

to_ase(species=None)

Convert system to an ASE Atoms object

Parameters:

species (list of string) – The chemical species

Return type:

None

to_file(outfile, format='lammps-dump', customkeys=None, customvals=None, compressed=False, timestep=0, species=None)

Save the system instance to a trajectory file.

Parameters:

• outfile (string) – name of the output file

• 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 .

property volume

Volume of box

pyscal.crystal_structures module

pyscal module for creating crystal structures.

class pyscal.crystal_structures.ElementCreator(element_dict)

Bases: object

Create an elementary structure

__init__(element_dict)

class pyscal.crystal_structures.LatticeCreator(element_dict)

Bases: ElementCreator

Create a lattice

class pyscal.crystal_structures.Structure

Bases: object

A class for structure creation

element

create elementary structures

lattice

create structures by specifying lattice

__init__()

structure_dict(structure)

pyscal.crystal_structures.general_lattice(species, positions, scaling_factors=[1.0, 1.0, 1.0], lattice_constant=1.0, repetitions=None, noise=0, element=None)

Create a general lattice structure.

species: list

list of per-atom species

positions:

list of relative positions positions of reach atom (between 0-1)

scaling_fractors:

factors with which the unit cell should be scaled, for example hcp could have [1,1.73, 1.63]. Default [1,1,1]

lattice_constantfloat, optional

lattice constant of the crystal structure, default 1

repetitionslist 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].

noisefloat, optional

If provided add normally distributed noise with standard deviation noise to the atomic positions.

elementstring, optional

The chemical element

pyscal.crystal_structures.structure_creator(structure, lattice_constant=1.0, repetitions=None, ca_ratio=1.633, noise=0, element=None)

Create a crystal structure and return it as a System object.

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.

• element (string, optional) – The chemical element

Returns:

System – system will be populated with given atoms and simulation box

Return type:

pyscal System

Examples

>>> sys = structure_creator('bcc', lattice_constant=3.48, repetitions=[2,2,2])

pyscal.trajectory module

class pyscal.trajectory.Timeslice(trajectory, blocklist)

Bases: object

Timeslice containing info about a single time slice Timeslices can also be added to each

__init__(trajectory, blocklist)

Initialize instance with data

to_ase(species=None)

Get block as Ase objects

Parameters:

blockno (int) – number of the block to be read, starts from 0

Returns:

sys

Return type:

ASE object

to_dict()

Get the required block as data

to_file(outfile, mode='w')

Get block as outputfile

Parameters:

• outfile (string) – name of output file

• mode (string) – write mode to be used, optional default “w” write also can be “a” to append.

Return type:

None

to_system(customkeys=None)

Get block as pyscal system

Parameters:

blockno (int) – number of the block to be read, starts from 0

Returns:

sys – pyscal System

Return type:

System

class pyscal.trajectory.Trajectory(filename)

Bases: object

A Trajectory class for LAMMPS

__init__(filename)

Initiaze the class

Parameters:

• filename (string) – name of the inputfile

• customkeys (list of string) – keys other than position, id that needs to be read in from the input file

get_block(blockno)

Get a block from the file as raw data

Parameters:

blockno (int) – number of the block to be read, starts from 0

Returns:

data – list of strings containing data

Return type:

list

load(blockno)

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.

unload(blockno)

Unload the data that is loaded to memory using load method

Parameters:

blockno (int) – number of the block to be read, starts from 0

Return type:

None

pyscal.traj_process module

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.

pyscal.traj_process.read_file(filename, format='lammps-dump', compressed=False, customkeys=None)

Read input file

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

pyscal.traj_process.split_trajectory(infile, format='lammps-dump', compressed=False)

Read in a trajectory file and convert it to individual time slices.

Parameters:

• filename (string) – name of input file

• format (format of the input file) – only lammps-dump is supported now.

• compressed (bool, optional) – force to read a gz zipped file. If the filename ends with .gz, use of this keyword is not necessary.

Returns:

snaps – a list of filenames which contain individual frames from the main trajectory.

Return type:

list of strings

Notes

This is a wrapper function around split_traj_lammps_dump function.

pyscal.traj_process.write_file(sys, outfile, format='lammps-dump', compressed=False, customkeys=None, customvals=None, timestep=0, species=None)

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.

Return type:

None

pyscal.misc module

pyscal.misc.compare_atomic_env(infile, atomtype=2, precision=2, format='poscar', print_results=True, return_system=False)

Compare the atomic environment of given types of atoms in the inputfile. The comparison is made in terms of Voronoi volume and Voronoi fingerprint.

Parameters:

• infile (string) – name of the inputfile

• atomtype (int, optional) – type of the atom default 2

• precision (float, optional) – precision for comparing Voronoi volumes default 3

• format (string, optional) – format of the input file default poscar

• print_results (bool, optional) – if True, print the results. If False, return the data instead. default True

• return_system (bool, optional) – if True, return the system object. default False

Returns:

• vvx (list of floats) – unique Voronoi volumes. Returned only if print results is False

• vrx (list of strings) – unique Voronoi polyhedra. Returned only if print results is False

• vvc (list of ints) – number of unique quantities specified above. Returned only if print results is False

pyscal.misc.find_tetrahedral_voids(infile, format='poscar', print_results=True, return_system=False, direct_coordinates=True, precision=0.1)

Check for tetrahedral voids in the system

Parameters:

• infile (string) – name of the input file

• format (string) – format of the input file, optional default poscar

• print_results (bool, optional) – if True, print the results. If False, return the data instead. default True

• return_system (bool, optional) – if True, return the system object. default False

• direct_coordinates (bool, optional) – if True, results are provided in direct coordinates default False

• precision (int, optional) – the number of digits to check for distances. default 1

Returns:

• types (list of atom types)

• volumes (list of atom volumes)

• pos (list of atom positions)

• sys (system object, returns only if return_sys is True)