mdapy package

Submodules

mdapy.ackland_jones_analysis module

class mdapy.ackland_jones_analysis.AcklandJonesAnalysis(pos, box, boundary=[1, 1, 1], verlet_list=None, distance_list=None)

Bases: object

This class applies Ackland Jones Analysis (AJA) method to identify the lattice structure.

The AJA method can recgonize the following structure:

1. Other

2. FCC

3. HCP

4. BCC

5. ICO

Note

If you use this module in publication, you should also cite the original paper. Ackland G J, Jones A P. Applications of local crystal structure measures in experiment and simulation[J]. Physical Review B, 2006, 73(5): 054104..

Hint

The module uses the legacy algorithm in LAMMPS.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• verlet_list (np.ndarray, optional) – (\(N_p\), >=14), first 14 neighbors is sorted, if not given, use kdtree to obtain it. Defaults to None.

• distance_list (np.ndarray, optional) – (\(N_p\), >=14), first 14 neighbors is sorted, if not given, use kdtree to obtain it. Defaults to None.

Outputs:

• aja (np.ndarray) - (\(N_p\)) AJA value per atoms.

• structure (list) - the corresponding structure to each aja value.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> AJA = mp.AcklandJonesAnalysis(FCC.pos, FCC.box, [1, 1, 1]) # Initialize AJA class

>>> AJA.compute() # Calculate the aja per atoms

>>> AJA.aja # Check the aja value

compute()

Do the real AJA calculation.

mdapy.atomic_strain module

class mdapy.atomic_strain.AtomicStrain(ref_pos, ref_box, cur_pos, cur_box, verlet_list, neighbor_number, boundary=[1, 1, 1], affi_map='off')

Bases: object

This class is used to calculate the atomic shear strain. More details can be found here. https://www.ovito.org/docs/current/reference/pipelines/modifiers/atomic_strain.html

Parameters:

• ref_pos (np.ndarray) – position of the reference configuration.

• ref_box (np.ndarray) – box of the reference configuration.

• cur_pos (np.array) – position of the current configuration.

• cur_box (np.array) – box of the current configuration.

• verlet_list (np.ndarray) – neighbor atom index for the reference configuration.

• distance_list (np.ndarray) – neighbor atom number for the reference configuration..

• boundary (list, optional) – boundary condition. Defaults to [1, 1, 1].

• affi_map (str, optional) – selected in [‘off’, ‘ref’]. If use to ‘ref’, the current position will affine to the reference frame. Defaults to ‘off’.

Outputs:

• shear_strain (np.ndarray) - shear strain value per atoms.

compute()

Do the real compute.

mdapy.bond_analysis module

class mdapy.bond_analysis.BondAnalysis(pos, box, boundary, rc, nbins, verlet_list=None, distance_list=None, neighbor_number=None)

Bases: object

This class calculates the distribution of bond length and angle based on a given cutoff distance.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(3, 2\)) or (\(4, 3\)) system box.

• boundary (list) – boundary conditions, 1 is periodic and 0 is free boundary.

• rc (float) – cutoff distance.

• nbins (int) – number of bins.

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number.

Outputs:

• bond_length_distribution (np.ndarray) - (nbins) bond length distribution.

• bond_angle_distribution (np.ndarray) - (nbins) bond angle distribution.

• r_length (np.ndarray) - (nbins) bond length (x axis).

• r_angle (np.ndarray) - (nbins) bond angle (x axis).

compute()

Do the real computation.

plot_bond_angle_distribution()

plot_bond_length_distribution()

mdapy.box module

mdapy.box.init_box(box)

This function is used to obtain box array.

• case 1: a float/int number, such as 10, it will generate a rectangle box with box length of 10 A.

• case 2: a str, such as “10”, the result is same with case 1.

• case 3: a list/tuple/np.ndarray.

If it includs three elements, such as [10, 20, 30], it will build a rectangle box with box length of x(10), y(20) and z(30). If the input is 2 2-D np.array or nested list. We accept the shape of (3, 2), (3, 3) and (4, 3). The shape of (3, 2) indicates the first colume is the origin and the second column is the maximum points of box. The shape of (3, 3) indicates the three box vector. The shape of (4, 3) indicates the three box vector and the fourth row is the origin position. The defaults origin is [0, 0, 0]. The final box is a np.ndarray with shape of (4, 3):

• ax ay az (x axis)

• bx by bz (y axis)

• cx cy cz (z axis)

• ox oy oz (origin)

Parameters:

box (float | str | list | np.ndarray) – The input box.

Returns:

box (4, 3), inverse box (3, 3), rec. The rec indicates if the box is reactangle.

Return type:

tuple[np.ndarray, np.ndarray, bool]

mdapy.cell_opt module

mdapy.centro_symmetry_parameter module

class mdapy.centro_symmetry_parameter.CentroSymmetryParameter(N, pos, box, boundary=[1, 1, 1], verlet_list=None)

Bases: object

This class is used to compute the CentroSymmetry Parameter (CSP), which is heluful to recgonize the structure in lattice, such as FCC and BCC. The CSP is given by:

\[p_{\mathrm{CSP}} = \sum_{i=1}^{N/2}{|\mathbf{r}_i + \mathbf{r}_{i+N/2}|^2},\]

where \(r_i\) and \(r_{i+N/2}\) are two neighbor vectors from the central atom to a pair of opposite neighbor atoms. For ideal centrosymmetric crystal, the contributions of all neighbor pairs will be zero. Atomic sites within a defective crystal region, in contrast, typically have a positive CSP value.

This parameter \(N\) indicates the number of nearest neighbors that should be taken into account when computing the centrosymmetry value for an atom. Generally, it should be a positive, even integer. Note that larger number decreases the calculation speed. For FCC is 12 and BCC is 8.

Note

If you use this module in publication, you should also cite the original paper. Kelchner C L, Plimpton S J, Hamilton J C. Dislocation nucleation and defect structure during surface indentation[J]. Physical review B, 1998, 58(17): 11085..

Hint

The CSP is calculated by the same algorithm as LAMMPS. First calculate all \(N (N - 1) / 2\) pairs of neighbor atoms, and the summation of the \(N/2\) lowest weights is CSP values.

Parameters:

• N (int) – Neighbor number.

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• verlet_list (np.ndarray, optional) – (\(N_p\), >=N), first N neighbors is sorted, if not given, use kdtree to obtain it. Defaults to None.

Outputs:

• csp (np.ndarray) - (\(N_p\)) CSP value per atoms.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> CSP = mp.CentroSymmetryParameter(12, FCC.pos, FCC.box, [1, 1, 1]) # Initialize CSP class

>>> CSP.compute() # Calculate the csp per atoms

>>> CSP.csp # Check the csp value

compute()

Do the real CSP calculation.

mdapy.cluser_analysis module

class mdapy.cluser_analysis.ClusterAnalysis(rc, verlet_list=None, distance_list=None, neighbor_number=None, pos=None, box=None, boundary=None, type_list=None)

Bases: object

This class is used to divide atoms connected within a given cutoff distance into a cluster. It is helpful to recognize the reaction products or fragments under shock loading.

Parameters:

• rc (float | dict) – cutoff distance. One can also assign multi cutoff for different elemental pair, such as {‘1-1’:1.5, ‘1-6’:1.7}. The unassigned elemental pair will default use the maximum cutoff distance.

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1. Defaults to None.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom. Defaults to None.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor number per atoms. Defaults to None.

• pos (np.ndarray, optional) – (\(N_p, 3\)) particles positions. Defaults to None.

• box (np.ndarray, optional) – (\(3, 2\)) or (\(4, 3\)) system box. Defaults to None.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Such as [1, 1, 1]. Defaults to None.

• type_list (np.ndarray, optional) – (\(N_p\)) atom type. It is needed only if rc is a Dict.

Outputs:

• particleClusters (np.ndarray) - (\(N_p\)) cluster ID per atoms.

• cluster_number (int) - cluster number.

Note

This class use the same method as in Ovito.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> FCC.pos = FCC.pos[((FCC.pos[:, 0]>1*3.615) & (FCC.pos[:, 0]<4*3.615))
                    | ((FCC.pos[:,0]>6*3.615) & (FCC.pos[:,0]< 9*3.615))] # Do a slice operation.

>>> neigh = mp.Neighbor(FCC.pos, FCC.box,
                        4., max_neigh=50) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> Clus = mp.ClusterAnalysis(4., neigh.verlet_list, neigh.distance_list, neigh.neighbor_number) # Initilize Cluster class.

>>> Clus.compute() # Do cluster calculation.

>>> Clus.cluster_number # Check cluster number, should be 2 here.

>>> Clus.particleClusters # Check atom in which cluster.

>>> Clus.get_size_of_cluster(1) # Obtain the atom number in cluster 1.

compute()

Do the real cluster analysis.

get_size_of_cluster(cluster_id)

This method can obtain the number of atoms in cluster \(cluster\_id\).

Parameters:

cluster_id (int) – cluster ID, it should be larger than 0 and lower than total \(cluster\_number\).

Returns:

the number of atoms in given cluster.

Return type:

int

mdapy.cluser_analysis.filter_by_type(verlet_list: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None), distance_list: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None), neighbor_number: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None), type_list: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None), type1: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None), type2: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None), r: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None))

mdapy.common_neighbor_analysis module

class mdapy.common_neighbor_analysis.CommonNeighborAnalysis(pos, box, boundary=[1, 1, 1], rc=None, verlet_list=None, neighbor_number=None)

Bases: object

This class use Common Neighbor Analysis (CNA) method to recgonize the lattice structure, based on which atoms can be divided into FCC, BCC, HCP and Other structure.

Note

If one use this module in publication, one should also cite the original paper. Stukowski, A. (2012). Structure identification methods for atomistic simulations of crystalline materials. Modelling and Simulation in Materials Science and Engineering, 20(4), 045021..

Hint

We use the same algorithm as in OVITO.

The original CNA method is sensitive to the given cutoff distance. The suggesting cutoff can be obtained from the following formulas:

\[r_{c}^{\mathrm{fcc}} = \frac{1}{2} \left(\frac{\sqrt{2}}{2} + 1\right) a \approx 0.8536 a,\]

\[r_{c}^{\mathrm{bcc}} = \frac{1}{2}(\sqrt{2} + 1) a \approx 1.207 a,\]

\[r_{c}^{\mathrm{hcp}} = \frac{1}{2}\left(1+\sqrt{\frac{4+2x^{2}}{3}}\right) a,\]

where \(a\) is the lattice constant and \(x=(c/a)/1.633\) and 1.633 is the ideal ratio of \(c/a\) in HCP structure.

Prof. Alexander Stukowski has improved this method using adaptive cutoff distances based on the atomic neighbor environment, which is the default method in mdapy from version 0.11.1.

The CNA method can recgonize the following structure:

1. Other

2. FCC

3. HCP

4. BCC

5. ICO

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• rc (float, optional) – cutoff distance, if not given, will using adaptive cutoff.

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1. Defaults to None.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number. Defaults to None.

Outputs:

• pattern (np.ndarray) - (\(N_p\)) CNA results.

• structure (list) - the corresponding structure to each pattern.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> neigh = mp.Neighbor(FCC.pos, FCC.box,
                        3.615*0.8536, max_neigh=20) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> CNA = mp.CommonNeighborAnalysis(FCC.pos, FCC.box, [1, 1, 1], 3.615*0.8536,
            neigh.verlet_list, neigh.neighbor_number) # Initialize CNA class

>>> CNA.compute() # Calculate the CNA per atoms

>>> CNA.pattern # Check results, should be 1 here.

>>> CNA.structure[CNA.pattern[0]] # Structure of atom 0, should be fcc here.

compute()

Do the real CNA calculation.

mdapy.common_neighbor_parameter module

class mdapy.common_neighbor_parameter.CommonNeighborParameter(pos, box, boundary, rc, verlet_list=None, distance_list=None, neighbor_number=None)

Bases: object

This class use Common Neighbor Parameter (CNP) method to recgonize the lattice structure.

Note

If one use this module in publication, one should also cite the original paper. Tsuzuki H, Branicio P S, Rino J P. Structural characterization of deformed crystals by analysis of common atomic neighborhood[J]. Computer physics communications, 2007, 177(6): 518-523..

CNP method is sensitive to the given cutoff distance. The suggesting cutoff can be obtained from the following formulas:

\[r_{c}^{\mathrm{fcc}} = \frac{1}{2} \left(\frac{\sqrt{2}}{2} + 1\right) a \approx 0.8536 a,\]

\[r_{c}^{\mathrm{bcc}} = \frac{1}{2}(\sqrt{2} + 1) a \approx 1.207 a,\]

\[r_{c}^{\mathrm{hcp}} = \frac{1}{2}\left(1+\sqrt{\frac{4+2x^{2}}{3}}\right) a,\]

where \(a\) is the lattice constant and \(x=(c/a)/1.633\) and 1.633 is the ideal ratio of \(c/a\) in HCP structure.

Some typical CNP values:

• FCC : 0.0

• BCC : 0.0

• HCP : 4.4

• FCC (111) surface : 13.0

• FCC (100) surface : 26.5

• FCC dislocation core : 11.

• Isolated atom : 1000. (manually assigned by mdapy)

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(3, 2\)) or (\(4, 3\)) system box.

• boundary (list) – boundary conditions, 1 is periodic and 0 is free boundary.

• rc (float) – cutoff distance.

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number.

Outputs:

• cnp (np.ndarray) - (\(N_p\)) CNP results.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> neigh = mp.Neighbor(FCC.pos, FCC.box,
                        3.615, max_neigh=20) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> CNP = mp.CommonNeighborParameter(FCC.pos, FCC.box, [1, 1, 1], 3.615*0.8536, neigh.verlet_list, neigh.distance_list,
            neigh.neighbor_number) # Initialize CNP class

>>> CNP.compute() # Calculate the CNP per atoms

>>> CNP.cnp[0] # Check results, should be 0.0 here.

compute()

Do the real CNP calculation.

mdapy.create_polycrystalline module

class mdapy.create_polycrystalline.Cell(face_vertices, vertices, volume, cavity_radius, face_areas, pos)

Bases: object

cavity_radius()

face_areas()

face_vertices()

vertices()

volume()

class mdapy.create_polycrystalline.Container(pos, box, boundary, num_t)

Bases: list

class mdapy.create_polycrystalline.CreatePolycrystalline(box, seednumber, filename=None, metal_latttice_constant=None, metal_lattice_type=None, randomseed=None, metal_overlap_dis=2.0, add_graphene=False, gra_lattice_constant=1.42, metal_gra_overlap_dis=3.1, gra_overlap_dis=1.2, seed=None, if_rotation=True, theta_list=None, face_threshold=5.0, output_name=None, num_t=None)

Bases: object

This class is used to create polycrystalline structure with random grain orientation based on the Voronoi diagram. It provides an interesting feature to replace the metallic grain boundary into graphene, which can generate a three-dimensional graphene structure. The metallic matrix can be FCC, BCC and HCP structure.

Parameters:

• box (np.ndarray) – (\(3, 2\)), system box.

• seednumber (int) – number of initial seed to generate the Voronoi polygon.

• metal_latttice_constant (float) – lattice constant.

• metal_lattice_type (str) – metallic matrix structure type, supporting FCC, BCC and HCP.

• randomseed (int, optional) – randomseed to generate the random number, this is helpful to reproduce the same structure. If not given, a random seed will be generated.

• metal_overlap_dis (float, optional) – minimum distance between metallic particles. If not given, this parameter will be determined by the given metal_lattice_type and metal_latttice_constant.

• add_graphene (bool, optional) – whether use graphene as grain boundary. Defaults to False.

• gra_lattice_constant (float, optional) – C-C bond length in graphene. Defaults to 1.42.

• metal_gra_overlap_dis (float, optional) – minimum distance between metallic particles and graphene. Defaults to 3.1.

• gra_overlap_dis (float, optional) – minimum distance between atoms in graphene. Defaults to 1.2.

• seed (np.ndarray, optional) – (seednumber, 3) initial position of seed to generate Voronoi polygon. If not given, it will be generated by the given randomseed.

• if_rotation (bool, optional) – whether rotate the grain orientation. Defaults to True.

• theta_list (np.ndarray, optional) – (seednumber, 3) rotation degree of each grain along x, y, z axis. If not given, it will be generated by the given randomseed.

• face_threshold (float, optional) – minimum voronoi polygon face area to add graphene. Defaults to 5.0.

• output_name (str, optional) – filename of DUMP file. Defaults to None.

• num_t (int, optional) – threads number to generate Voronoi diagram. If not given, use all avilable threads.

Outputs:

• generate an polycrystalline DUMP file with grain ID.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> box = np.array([[0.0, 200.0], [0.0, 200.0], [0.0, 200.0]]) # Generate a box.

>>> polycry = mp.CreatePolycrystalline(
              box, 20, 3.615, "FCC", add_graphene=True) # Initilize a Poly class.

>>> polycry.compute() # Generate a graphene/metal structure with 20 grains.

compute(save_dump=True)

Do the real polycrystalline structure building.

class mdapy.create_polycrystalline.DeleteOverlap(pos, box, rc)

Bases: object

compute()

mdapy.dft2nepxyz module

class mdapy.dft2nepxyz.DFT2NEPXYZ(filename_list, fmt='CP2K-SCF', interval=10, energy_shift=None, save_virial=True, force_max=None, stress_max=None, mode='w')

Bases: object

This class is used to generate NEP training XYZ file from many DFT data. Now we only support SCF calculation in CP2K. In the future the AIMD and VASP may also be implemented.

Parameters:

• filename_list (list) – all DFT output file you want to save, such as [‘a/output.log’, ‘b/output.log’]

• fmt (str, optional) – DFT calculation code. Defaults to “CP2K-SCF”.

• interval (int, optional) – if provided, we will save it to test.xyz per interval. Defaults to 10.

• energy_shift (dict, optional) – unit is eV. Ff provided, the energy will substract the base energy, such as {‘Fe’:89.0, ‘O’:50.0}. Defaults to None.

• save_virial (bool, optional) – if set False, the virial information will not be saved. Defaults to True.

• force_max (float, optional) – if system’s abusolute maximum force (in eV/A) is larger than this value, it will not be saved to train.xyz.

• stress_max (float, optional) – if system’s abusolute maximum stress (in GPa) is larger than this value, it will not be saved to train.xyz.

• mode (str, optional) – if mode is ‘w’, it will generate new train.xyz/test.xyz. If mode is ‘a’, it will append in current train.xyz/test.xyz. Defaults to ‘w’.

Outputs:

• Generate train.xyz and test.xyz in current folder (if interval provides).

class mdapy.dft2nepxyz.LabeledSystem(filename, fmt='CP2K-SCF')

Bases: object

This class is used to read first principle calculation data, obataining the energy, force, box and virial information. Those information can be saved as initial database for deep learning training, aiming to develop high accurancy potential function.

The units are listed as below:

• energy : eV (per-cell)

• force : eV/Å (per-atom)

• virial : eV (per-cell)

• stress : GPa (per-cell)

• pos : Å (per-atom)

• box : Å

Now we only support SCF calculation in CP2K. In the future the AIMD and VASP may also be implemented.

Parameters:

• filename (str) – filename of CP2K SCF output file.

• fmt (str, optional) – DFT calculation code. Defaults to “CP2K-SCF”.

Outputs:

• content (str) - the whole content of input file.

• structure (dict) - a dict contains energy, pos, box, force, type_list and virial (if computed).

Examples

>>> import mdapy as mp

>>> mp.init()

>>> LS = mp.LabeledSystem('output.log')

>>> LS.data['energy'] # check energy per cell.

get_data()

mdapy.eam_average module

class mdapy.eam_average.EAMAverage(filename, concentration, output_name=None)

Bases: EAM

This class is used to generate the average EAM (A-atom) potential, which is useful in alloy investigation. The A-atom potential has the similar formula with the original EAM potential:

\[E_{i}=\sum_{i} F^{A}\left(\bar{\rho}_{i}\right)+\frac{1}{2} \sum_{i, j \neq i} \phi_{i j}^{A A},\]

\[F^{A}\left(\bar{\rho}_{i}\right)=\sum_{\alpha} c_{\alpha} F^{\alpha}\left(\bar{\rho}_{i}\right),\]

\[\phi_{i j}^{A A}=\sum_{\alpha, \beta } c_{\alpha} c_{\beta} \phi_{i j}^{\alpha \beta},\]

\[\quad \bar{\rho}_{i}=\sum_{j \neq i} \sum_{\alpha} c_{\alpha} \rho_{i j}^{\alpha},\]

where \(A\) denotes an average-atom.

Note

If you use this module in publication, you should also cite the original paper. Average-atom interatomic potential for random alloys

Parameters:

• filename (str) – filename of eam.alloy file.

• concentration (list) – atomic ratio list, such as [0.5, 0.5] and the summation should be equal to 1.

• output_name (str, optional) – filename of generated average EAM potential.

Outputs:

• generate an averaged eam.alloy potential file with A element.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> potential = mp.EAMAverage("./example/CoNiFeAlCu.eam.alloy",
                [0.2, 0.2, 0.2, 0.2, 0.2]) # Generate the EAMAverage class.

>>> potential.plot() # plot the results.

mdapy.eam_generate module

class mdapy.eam_generate.EAMGenerate(elements_list, output_name=None)

Bases: object

This class is used to create EAM potential including one or more elements. This is the python version translated from the fortran version by Zhou et. al.

Note

If you use this module in publication, you should also cite the original paper. Misfit-energy-increasing dislocations in vapor-deposited CoFe/NiFe multilayers

Parameters:

• elements_list (list) – elements list, such as [‘Co’, ‘Ni’], wchich should choose from [“Cu”,”Ag”,”Au”,”Ni”,”Pd”,”Pt”,”Al”,”Pb”,”Fe”,”Mo”,”Ta”,”W”,”Mg”,”Co”,”Ti”,”Zr”].

• output_name (str, optional) – filename of generated EAM file.

Outputs:

• generate an eam.alloy potential file.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> potential = mp.EAMGenerate(["Co", "Ni", "Fe"]) # Generate the EAMGenerate class.

>>> potential.plot() # plot the results.

mdapy.entropy module

class mdapy.entropy.AtomicEntropy(vol=None, verlet_list=None, distance_list=None, pos=None, box=None, boundary=None, rc=5.0, sigma=0.2, use_local_density=False, compute_average=False, average_rc=None)

Bases: object

This class is used to calculate the entropy fingerprint, which is useful to distinguish between ordered and disordered environments, including liquid and solid-like environments, or glassy and crystalline-like environments. The potential application could identificate grain boundaries or a solid cluster emerging from the melt. One of the advantages of this parameter is that no a priori information of structure is required.

This parameter for atom \(i\) is computed using the following formula:

\[s_S^i=-2\pi\rho k_B \int\limits_0^{r_m} \left [ g(r) \ln g(r) - g(r) + 1 \right ] r^2 dr,\]

where \(r\) is a distance, \(g(r)\) is the radial distribution function, and \(\rho\) is the density of the system. The \(g(r)\) computed for each atom \(i\) can be noisy and therefore it can be smoothed using:

\[g_m^i(r) = \frac{1}{4 \pi \rho r^2} \sum\limits_{j} \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(r-r_{ij})^2/(2\sigma^2)},\]

where the sum over \(j\) goes through the neighbors of atom \(i\) and \(\sigma\) is a parameter to control the smoothing. The average of the parameter over the neighbors of atom \(i\) is calculated according to:

\[\left< s_S^i \right> = \frac{\sum_j s_S^j + s_S^i}{N + 1},\]

where the sum over \(j\) goes over the neighbors of atom \(i\) and \(N\) is the number of neighbors. The average version always provides a sharper distinction between order and disorder environments.

Note

If you use this module in publication, you should also cite the original paper. Entropy based fingerprint for local crystalline order

Note

This class uses the same algorithm with LAMMPS.

Tip

Suggestions for FCC, the \(rc = 1.4a\) and \(average\_rc = 0.9a\) and for BCC, the \(rc = 1.8a\) and \(average\_rc = 1.2a\), where the \(a\) is the lattice constant.

Parameters:

• vol (float, optional) – system volume. Defaults to None.

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1. Defaults to None.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom. Defaults to None.

• pos (np.ndarray, optional) – (\(N_p, 3\)) particles positions. Defaults to None.

• box (np.ndarray, optional) – (\(3, 2\)) or (\(4, 3\)) system box. Defaults to None.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Such as [1, 1, 1]. Defaults to None.

• rc (float, optional) – cutoff distance. Defaults to 5.0.

• sigma (float, optional) – smoothing parameter. Defaults to 0.2.

• use_local_density (bool, optional) – whether use local atomic volume. Defaults to False.

• compute_average (bool, optional) – whether compute the average version. Defaults to False.

• average_rc (float, optional) – cutoff distance for averaging operation, if not given, it is equal to rc. This parameter should be lower than rc.

Outputs:

• entropy (np.ndarray) - (\(N_p\)), atomic entropy.

• entropy_average (np.ndarray) - (\(N_p\)), averaged atomic entropy if compute_average is True.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> import numpy as np

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> neigh = mp.Neighbor(FCC.pos, FCC.box,
                        3.615*1.4, max_neigh=80) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> vol = np.product(FCC.box[:, 1] - FCC.box[:, 0]) # Calculate system volume.

>>> Entropy = mp.AtomicEntropy(
        vol,
        neigh.verlet_list,
        neigh.distance_list,
        None,
        None,
        None,
        neigh.rc,
        sigma=0.2,
        use_local_density=False,
        compute_average=True,
        average_rc=3.615*0.9,
    ) # Initilize the entropy class.

>>> Entropy.compute() # Calculate the atomic entropy.

>>> Entropy.entropy # Check atomic entropy.

>>> Entropy.entropy_average # Check averaged atomic entropy.

compute()

Do the real entropy calculation.

mdapy.identify_SFs_TBs module

class mdapy.identify_SFs_TBs.IdentifySFTBinFCC(structure_types, ptm_indices)

Bases: object

This class is used to identify the stacking faults (SFs) and coherent twin boundaries (TBs) in FCC structure based on the Polyhedral Template Matching (PTM). It can identify the following structure:

1. 0 = Non-hcp atoms (e.g. perfect fcc or disordered)

2. 1 = Indeterminate hcp-like (isolated hcp-like atoms, not forming a planar defect)

3. 2 = Intrinsic stacking fault (two adjacent hcp-like layers)

4. 3 = Coherent twin boundary (one hcp-like layer)

5. 4 = Multi-layer stacking fault (three or more adjacent hcp-like layers)

6. 5 = Extrinsic stacking fault

Note

This class is translated from that implementation in Ovito but optimized to be run parallely. And so-called multi-layer stacking faults maybe a combination of intrinsic stacking faults and/or twin boundary which are located on adjacent {111} plane. It can not be distiguished by the current method.

Parameters:

• structure_types (np.ndarray) – (\(N_p\)) structure type from ptm method.

• ptm_indices (np.ndarray) – (\(N_p, 18\)) ptm ordered verlet_list for all atoms.

Outputs:

• fault_types (np.ndarray) - (\(N_p\)) planar faults types.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> system = mp.System(r"./example/ISF.dump") # Read dump

>>> ptm = mp.PolyhedralTemplateMatching(
            system.pos, system.box, system.boundary, "default", 0.1, None, True
            ) # Initilize ptm method.

>>> ptm.compute() # Compute ptm per atoms.

>>> structure_type = np.array(ptm.output[:, 0], int) # Obtain ptm structure

>>> SFTB = mp.IdentifySFTBinFCC(structure_type, ptm.ptm_indices) # Initialize SFTB class

>>> SFTB.compute() # Compute the planar faults type

>>> SFTB.fault_types # Check results

compute()

Do the real computation.

mdapy.identify_diamond_structure module

class mdapy.identify_diamond_structure.IdentifyDiamondStructure(pos, box, boundary=[1, 1, 1], verlet_list=None)

Bases: object

This class is used to identify the Diamond structure. The results and algorithm should be the same in Ovito. More details can be found in https://www.ovito.org/manual/reference/pipelines/modifiers/identify_diamond.html .

Parameters:

• pos (np.ndarray) – atom position.

• (np.ndarray (box) – system box.

• boundary (list, optional) – boundary. Defaults to [1, 1, 1].

• verlet_list (np.ndarray, optional) – first 12 neighbors is sorted, if not given, use kdtree to obtain it. Defaults to None.

Outputs:

• pattern (np.ndarray) - structure pattern per atoms.

The identified structures include:

• 0 “other”,

• 1 “cubic_diamond”,

• 2 “cubic_diamond_1st_neighbor”,

• 3 “cubic_diamond_2st_neighbor”,

• 4 “hexagonal_diamond”,

• 5 “hexagonal_diamond_1st_neighbor”,

• 6 “hexagonal_diamond_2st_neighbor”

compute()

mdapy.lattice_maker module

class mdapy.lattice_maker.LatticeMaker(lattice_constant=None, lattice_type=None, x=1, y=1, z=1, crystalline_orientation=None, basis_vector=None, basis_atoms=None, type_list=None)

Bases: object

This class is used to create some standard lattice structure.

Parameters:

• lattice_constant (float) – lattice constant \(a\).

• lattice_type (str) – lattice type, seleted in [‘FCC’, ‘BCC’, ‘HCP’, ‘GRA’]. Here the HCP is ideal structure and the \(c/a=1.633\).

• x (int) – repeat times along \(x\) axis.

• y (int) – repeat times along \(y\) axis.

• z (int) – repeat times along \(z\) axis.

• crystalline_orientation (np.ndarray, optional) – (\(3, 3\)). Crystalline orientation, only support for ‘FCC’ and ‘BCC’ lattice. If not given, the orientation is x[1, 0, 0], y[0, 1, 0], z[0, 0, 1].

• basis_vector (np.ndarray) – (\(4, 3\)) or (\(3, 2\)) repeat vector.

• basis_atoms (np.ndarray) – (\(N_p, 3\)) basis atom positions.

• type_list (np.ndarray) – (\(N_p\)) type for basis atoms.

Outputs:

• box (np.ndarray) - (\(4, 3\)), system box.

• pos (np.ndarray) - (\(N_p, 3\)), particle position.

• type_list (np.ndarray) - (\(N_p\)), replicated type list.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure.

>>> FCC.compute() # Get atom positions.

>>> FCC.write_data() # Save to DATA file.

>>> GRA = mp.LatticeMaker(1.42, 'GRA', 10, 20, 1) # Create a graphene structure.

>>> GRA.write_data(output_name='graphene.data') # Save to graphene.data file.

>>> BCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10,
          crystalline_orientation=np.array([[1, 1, 1], [1, -1, 0], [1, 1, -2]])) # Create a BCC lattice with special orientations.

>>> BCC.compute() # Get atom positions.

>>> BCC.write_dump() # Save to dump file.

property N

The particle number.

compute()

Do the real lattice calculation.

property vol

The box volume.

write_POSCAR(output_name=None, type_name=None, reduced_pos=False)

This function writes position into a POSCAR file.

Parameters:

• output_name (str, optional) – filename of generated POSCAR file.

• type_name (list, optional) – species name. Such as [‘Al’, ‘Fe’].

• reduced_pos (bool, optional) – whether save directed coordination. Defaults to False.

write_cif(output_name=None, type_name=None)

This function writes position into a cif file.

Parameters:

• output_name (str, optional) – filename of generated cif file.

• type_name (list, optional) – species name. Such as [‘Al’, ‘Fe’].

write_data(output_name=None, data_format='atomic', type_list=None, num_type=None, type_name=None)

This function writes position into a DATA file.

Parameters:

• output_name (str, optional) – filename of generated DATA file.

• data_format (str, optional) – data format, selected in [‘atomic’, ‘charge’].

• type_list (np.ndarray, optional) – one can mannually assign the type_list.

• num_type (int, optional) – explictly assign a number of atom type. Defaults to None.

• type_name (list, optional) – explictly assign elemantal name list, such as [‘Al’, ‘C’]. Defaults to None.

write_dump(output_name=None, compress=False)

This function writes position into a DUMP file.

Parameters:

• output_name (str, optional) – filename of generated DUMP file.

• compress (bool, optional) – whether compress the DUMP file.

write_xyz(output_name=None, type_name=None, classical=False)

This function writes position into a xyz file.

Parameters:

• output_name (str, optional) – filename of generated xyz file. Defaults to None.

• type_name (list, optional) – assign the species name. Such as [‘Al’, ‘Cu’]. Defaults to None.

• classical (bool, optional) – whether save with classical format. Defaults to False.

mdapy.lindemann_parameter module

class mdapy.lindemann_parameter.LindemannParameter(pos_list, only_global=False)

Bases: object

This class is used to calculate the Lindemann index, which is useful to distinguish the melt process and determine the melting points of nano-particles. The Lindemann index is defined as the root-mean-square bond-length fluctuation with following mathematical expression:

\[\left\langle\sigma_{i}\right\rangle=\frac{1}{N_{p}(N_{p}-1)} \sum_{j \neq i} \frac{\sqrt{\left\langle r_{i j}^{2}\right\rangle_t-\left\langle r_{i j}\right\rangle_t^{2}}}{\left\langle r_{i j}\right\rangle_t},\]

where \(N_p\) is the particle number, \(r_{ij}\) is the distance between atom \(i\) and \(j\) and brackets \(\left\langle \right\rangle_t\) represents an time average.

Note

This class is partly referred to a work on calculating the Lindemann index.

Note

This calculation is high memory requirement. One can estimate the memory by: \(2 * 8 * N_p^2 / 1024^3\) GB.

Tip

If only global lindemann index is needed, the class can be calculated in parallel. The local Lindemann index only run serially due to the dependencies between different frames. Here we use the Welford method to update the varience and mean of \(r_{ij}\).

Parameters:

• pos_list (np.ndarray) – (\(N_f, N_p, 3\)), \(N_f\) frames particle position, which need to be unwrapped for periodic boundary.

• only_global (bool, optional) – whether only calculate the global index. Defaults to False.

Outputs:

• lindemann_atom (np.ndarray) - (\(N_f, N_p\)), local Lindemann index per atoms.

• lindemann_frame (np.ndarray) - (\(N_f\)), Lindemann index per frames.

• lindemann_trj (float) - global Lindemann index for an entire trajectory.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> import numpy as np

>>> pos_list = np.cumsum(
    np.random.choice([-1.0, 0.0, 1.0], size=(200, 1000, 3)), axis=0
    ) # Generate a random walk trajectory with 200 frames and 1000 particles.

>>> LDMG = mp.LindemannParameter(pos_list, only_global=True) # Generate a Lindemann class.

>>> LDMG.compute() # Only calculate the global Lindemann index, it will be much faster.

>>> LDML = mp.LindemannParameter(pos_list) # Generate a Lindemann class.

>>> LDML.compute() # Calculate the global and local Lindemann index.

>>> np.isclose(LDML.lindemann_trj, LDMG.lindemann_trj) # Should return True.

>>> LDML.lindemann_frame # Check Lindemann index per frame.

>>> LDML.plot() # Plot the evolution of Lindemann index per frame.

compute()

Do the real Lindemann index calculation.

plot()

Plot the evolution of Lindemann index per frame.

Raises:

Exception – One should compute lidemann_frame first!

Returns:

(fig, ax) matplotlib figure and axis class.

Return type:

tuple

mdapy.load_save_data module

class mdapy.load_save_data.BuildSystem

Bases: object

static fromarray(pos, box, boundary, vel, type_list)

static fromdata(data, box, boundary)

classmethod fromfile(filename, fmt)

static getformat(filename, fmt=None)

static read_POSCAR(filename)

static read_cif(filename)

static read_data(filename)

static read_dump(filename, fmt)

static read_xyz(filename)

class mdapy.load_save_data.SaveFile

Bases: object

static write_POSCAR(output_name, box, data, type_name=None, reduced_pos=False, selective_dynamics=False, save_velocity=False)

static write_cif(output_name, box, data, type_name=None)

static write_cp2k(output_name, box, boundary, data=None, pos=None, type_name=None)

static write_data(output_name, box, data=None, pos=None, type_list=None, num_type=None, type_name=None, data_format='atomic')

static write_dump(output_name, box, boundary, data=None, pos=None, type_list=None, timestep=0, compress=False)

static write_xyz(output_name, box, data, boundary, classical=False, **kargs)

mdapy.main module

mdapy.main.ackland_jones_analysis()

mdapy.main.atomic_entropy()

mdapy.main.atomic_temperature()

mdapy.main.box2axis(box)

mdapy.main.box2lines(box)

mdapy.main.build_lattice()

mdapy.main.build_polycrystal()

mdapy.main.callback()

mdapy.main.centro_symmetry_parameter()

mdapy.main.cluster_analysis()

mdapy.main.common_neighbor_analysis()

mdapy.main.common_neighbor_parameter()

mdapy.main.identify_SFs_TBs()

mdapy.main.identify_species()

mdapy.main.init_global_parameters()

mdapy.main.load_file_gui()

mdapy.main.loadfile(filename)

mdapy.main.main()

mdapy.main.polyhedral_template_matching()

mdapy.main.radiul_distribution_function()

mdapy.main.replicate()

mdapy.main.save_file_POSCAR()

mdapy.main.save_file_cif()

mdapy.main.save_file_data()

mdapy.main.save_file_dump()

mdapy.main.save_file_gui()

mdapy.main.save_file_xyz()

mdapy.main.show_system_info()

mdapy.main.spatial_binning()

mdapy.main.steinhardt_bond_orientation()

mdapy.main.void_analysis()

mdapy.main.voronoi_analysis()

mdapy.main.warren_cowley_parameter()

mdapy.main.whether_init_file(filename)

mdapy.mean_squared_displacement module

class mdapy.mean_squared_displacement.MeanSquaredDisplacement(pos_list, mode='windows')

Bases: object

This class is used to calculate the mean squared displacement MSD of system, which can be used to reflect the particle diffusion trend and describe the melting process. Generally speaking, MSD is an average displacement over all windows of length \(m\) over the course of the simulation (so-called ‘windows’ mode here) and defined by:

\[MSD(m) = \frac{1}{N_{p}} \sum_{i=1}^{N_{p}} \frac{1}{N_{f}-m} \sum_{k=0}^{N_{f}-m-1} (\vec{r}_i(k+m) - \vec{r}_i(k))^2,\]

where \(r_i(t)\) is the position of particle \(i\) in frame \(t\). It is computationally extensive while using a fast Fourier transform can remarkably reduce the computation cost as described in nMoldyn - Interfacing spectroscopic experiments, molecular dynamics simulations and models for time correlation functions and discussion in StackOverflow.

Note

One can install pyfftw to accelerate the calculation, otherwise mdapy will use scipy.fft to do the Fourier transform.

Sometimes one only need the following atomic displacement (so-called ‘direct’ mode here):

\[MSD(t) = \dfrac{1}{N_p} \sum_{i=1}^{N_p} (r_i(t) - r_i(0))^2.\]

Parameters:

• pos_list (np.ndarray) – (\(N_f, N_p, 3\)), \(N_f\) frames particle position, which need to be unwrapped for periodic boundary.

• mode (str, optional) – calculation mode, selected from [‘windows’, ‘direct’]. Defaults to “windows”.

Outputs:

• msd (np.ndarray) - (\(N_f\)), mean squared displacement per frames.

• particle_msd (np.ndarray) - (\(N_f, N_p\)), mean squared displacement per atom per frames.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> import numpy as np

>>> pos_list = np.cumsum(
    np.random.choice([-1.0, 0.0, 1.0], size=(200, 1000, 3)), axis=0
    ) # Generate a random walk trajectory with 200 frames and 1000 particles.

>>> MSD = mp.MeanSquaredDisplacement(pos_list, mode="windows") # Initilize MSD class.

>>> MSD.compute() # Calculate the MSD in 'windows' mode.

>>> MSD.msd # Check msd.

>>> MSD.particle_msd.shape # Check msd per particle, should be (200, 1000) here.

>>> MSD.plot() # Plot the evolution of msd per frame.

compute()

Do the real MSD calculation.

plot()

Plot the evolution of MSD per frame.

Returns:

(fig, ax) matplotlib figure and axis object.

Return type:

tuple

mdapy.minimizer module

class mdapy.minimizer.Minimizer(pos, box, boundary, potential, elements_list, type_list, mini_type='FIRE', fmax=1e-05, max_itre=200, volume_change=False, hydrostatic_strain=False)

Bases: object

This function use the fast inertial relaxation engine (FIRE) method to minimize the system, including optimizing position and box. More details can be found in paper:

Guénolé, Julien, et al. “Assessment and optimization of the fast inertial relaxation engine (fire) for energy minimization in atomistic simulations and its implementation in lammps.” Computational Materials Science 175 (2020): 109584.

Parameters:

• pos (np.ndarray) – atom position.

• box (np.ndarray) – system box.

• boundary (list) – boundary, such as [1, 1, 1].

• potential (BasePotential) – mdapy potential.

• elements_list (list) – elemental name, such as [‘Al’, ‘C’].

• type_list (np.ndarray) – atom type list.

• mini_type (str, optional) – only support ‘FIRE’ now. Defaults to “FIRE”.

• fmax (float, optional) – maximum force per atom to consider as converged. Defaults to 1e-5.

• max_itre (int, optional) – maximum iteration times. Defaults to 200.

• volume_change (bool, optional) – whether change the box to optimize the pressure. Defaults to False.

• hydrostatic_strain (bool, optional) – sonstrain the cell by only allowing hydrostatic deformation. Defaults to False.

compute()

mdapy.nearest_neighbor module

class mdapy.nearest_neighbor.NearestNeighbor(pos, box, boundary=[1, 1, 1])

Bases: object

This class is used to query the nearest neighbor with fixed number. For rectangle box, this class is a wrapper of kdtree of scipy and helful to obtain the certain nearest atom neighbors considering the periodic/free boundary. If you want to access the atom neighbor within a spherical distance, the Neighbor class is suggested. For triclinic box, this class use cell-list to find the nearest neighbors.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)), particles positions.

• box (np.ndarray) – (\(4, 3\) or \(3, 2\)), system box.

• boundary (list) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> kdt = mp.kdtree(FCC.pos, FCC.box, [1, 1, 1]) # Build a kdtree.

>>> dis, index = kdt.query_nearest_neighbors(12) # Query the 12 nearest neighbors per atom.

query_nearest_neighbors(K, workers=-1)

Query the \(n\) nearest atom neighbors.

Parameters:

• K (int) – number of neighbors to query.

• worker (int) – maximum cores used. Defaults to -1, indicating use all aviliable cores. Only works for scipy backend.

Returns:

(distance, index), distance of atom \(i\) to its neighbor atom \(j\), and the index of atom \(j\).

Return type:

tuple

class mdapy.nearest_neighbor.kdtree(pos, box, boundary)

Bases: object

This class is a wrapper of kdtree of scipy and helful to obtain the certain nearest atom neighbors considering the periodic/free boundary. If you want to access the atom neighbor within a spherical distance, the Neighbor class is suggested.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)), particles positions.

• box (np.ndarray) – (\(3, 2\)), system box.

• boundary (list) – boundary conditions, 1 is periodic and 0 is free boundary. Such as [1, 1, 1].

Outputs:

• shift_pos (np.ndarray) - (\(N_p, 3\)), shifted position, making the lower cornor is zero.

• kdt (scipy.spatial.KDTree) - a KDTree class.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> kdt = mp.kdtree(FCC.pos, FCC.box, [1, 1, 1]) # Build a kdtree.

>>> dis, index = kdt.query_nearest_neighbors(12) # Query the 12 nearest neighbors per atom.

query_nearest_neighbors(K, workers=-1)

Query the \(n\) nearest atom neighbors.

Parameters:

• K (int) – number of neighbors to query.

• worker (int) – maximum cores used. Defaults to -1, indicating use all aviliable cores. Only works for scipy backend.

Returns:

(distance, index), distance of atom \(i\) to its neighbor atom \(j\), and the index of atom \(j\).

Return type:

tuple

mdapy.neighbor module

class mdapy.neighbor.Neighbor(pos, box, rc, boundary=[1, 1, 1], max_neigh=None)

Bases: object

This module is used to cerate neighbor of atoms within a given cutoff distance. Using linked-cell list method makes fast neighbor finding possible.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• rc (float) – cutoff distance.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• max_neigh (int, optional) – maximum neighbor number. If not given, will estimate atomatically. Default to None.

Outputs:

• verlet_list (np.ndarray) - (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1.

• distance_list (np.ndarray) - (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom.

• neighbor_number (np.ndarray) - (\(N_p\)) neighbor atoms number.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> import numpy as np

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> neigh = mp.Neighbor(FCC.pos, FCC.box,
                        5.) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> neigh.verlet_list # Check the neighbor index.

>>> neigh.distance_list # Check the neighbor distance.

>>> neigh.neighbor_number # Check the neighbor atom number.

compute()

Do the real neighbor calculation.

sort_verlet_by_distance(N: int)

This function sorts the first N-th verlet_list by distance_list.

Parameters:

N (int) – number of sorted values

mdapy.orthogonal_box module

class mdapy.orthogonal_box.OrthogonalBox(pos, box, type_list=None)

Bases: object

This class try to change the box to rectangular.

Parameters:

• pos (np.ndarray) – atom position.

• box (np.ndarray) – system box.

• type_list (np.ndarray, optional) – type list. Defaults to None.

Outputs:

• rec_pos (np.ndarray) - new position.

• rec_box (np.ndarray) - new box.

• rec_type_list (np.ndarray) - new type_list.

compute(N=10)

Do the real computation.

Parameters:

N (int, optional) – search limit. If you can’t found rectangular box, increase N. Defaults to 10.

mdapy.pair_distribution module

class mdapy.pair_distribution.PairDistribution(rc, nbin, box, boundary, verlet_list=None, distance_list=None, neighbor_number=None, pos=None, type_list=None)

Bases: object

This class is used to calculate the radiul distribution function (RDF),which reflects the probability of finding an atom at distance r. The seperate pair-wise combinations of particle types can also be computed:

\[g(r) = c_{\alpha}^2 g_{\alpha \alpha}(r) + 2 c_{\alpha} c_{\beta} g_{\alpha \beta}(r) + c_{\beta}^2 g_{\beta \beta}(r),\]

where \(c_{\alpha}\) and \(c_{\beta}\) denote the concentration of two atom types in system and \(g_{\alpha \beta}(r)=g_{\beta \alpha}(r)\).

Parameters:

• rc (float) – cutoff distance.

• nbin (int) – number of bins.

• box (np.ndarray) – (\(3, 2\)) or (\(4, 3\)) system box. Defaults to None.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Such as [1, 1, 1]. Defaults to [1, 1, 1].

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number.

• pos (np.ndarray, optional) – (\(N_p, 3\)) particles positions. Defaults to None.

• type_list (np.ndarray, optional) – (\(N_p\)) atom type list. If not given, all atoms types are set as 1.

Outputs:

• r (np.ndarray) - (nbin), distance.

• g_total (np.ndarray) - (nbin), global RDF.

• Ntype (int) - number of species.

• g (np.ndarray): (\(N_{type}, N_{type}, nbin\)), partial RDF.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure.

>>> FCC.compute() # Get atom positions.

>>> neigh = mp.Neighbor(FCC.pos, FCC.box,
                        5., max_neigh=80) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> gr = mp.PairDistribution(neigh.rc, 200, FCC.box, [1, 1, 1],
                          neigh.verlet_list, neigh.distance_list, neigh.neighbor_number,) # Initilize the RDF class.

>>> gr.compute() # Calculate the RDF.

>>> gr.g_total # Check global RDF.

>>> gr.g # Check partial RDF, here is the same due to only one species.

>>> gr.r # Check the distance.

>>> gr.plot() # Plot global RDF.

>>> gr.plot_partial() # Plot partial RDF.

compute()

Do the real RDF calculation.

plot()

Plot the global RDF.

Returns:

(fig, ax) matplotlib figure and axis class.

Return type:

tuple

plot_partial(elements_list=None)

Plot the partial RDF.

Parameters:

elements_list (list, optional) – species element names list, such as [‘Al’, ‘Ni’].

Returns:

(fig, ax) matplotlib figure and axis class.

Return type:

tuple

mdapy.perturb_model module

class mdapy.perturb_model.PerturbModel(filename=None, lattice_type=None, lattice_constant=None, x=1, y=1, z=1, crystalline_orientation=None, scale_list=None, pert_num=None, pert_box=None, pert_atom=None, type_list=None, type_name=None, save_path='res', save_type='cp2k', fmt=None)

Bases: object

This class is used to generate atomic model with small geometry perturb, which is helpful to preparing database for deep learning training.

Parameters:

• filename (str, optional) – filename one wants to perturb. Defaults to None.

• lattice_type (str, optional) – lattice type, selected in [‘FCC’, ‘BCC’, ‘HCP’, ‘GRA’]. This parameter will be ignored if filename is provides. Defaults to None.

• lattice_constant (float, optional) – lattice constant. Defaults to None.

• x (int, optional) – replicate along x direction. Defaults to 1.

• y (int, optional) – replicate along y direction. Defaults to 1.

• z (int, optional) – replicate along z direction. Defaults to 1.

• crystalline_orientation (np.ndarray, optional) – (\(3, 3\)). Crystalline orientation, only support for ‘FCC’ and ‘BCC’ lattice. If not given, the orientation is x[1, 0, 0], y[0, 1, 0], z[0, 0, 1]. Defaults to None.

• scale_list (list, optional) – one can scale system isotropicly, such as [0.9, 1.0, 1.1]. Defaults to None.

• pert_num (int, optional) – number of models per-scale, such as 20. Defaults to None.

• pert_box (float, optional) – perturb on box, such as 0.03. Defaults to None.

• pert_atom (float, optional) – perturb on atom, such as 0.03. Defaults to None.

• type_list (list[int], optional) – type list only for generated lattice, such as [1, 1, 1, 2] for FCC lattice. Defaults to None.

• type_name (list[str], optional) – element name, such as [‘Al’, ‘Fe’]. Defaults to None.

• save_path (str, optional) – save path. Defaults to “res”.

• save_type (str, optional) – selected in [‘cp2k’, ‘vasp’]. For ‘vasp’, POSCAR will be saved. For ‘cp2k’, cif and xyz will be saved. Defaults to “cp2k”.

• fmt (str, optional) – this explitly gives the file format for filename, selected in [“data”, “lmp”, “POSCAR”, “xyz”, “cif”]. Defaults to None.

Outputs:

• Generate a series of perturb model in save_path.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> pert = mp.PerturbModel(
        lattice_type="FCC",
        lattice_constant=4.05,
        x=3,
        y=3,
        z=3,
        scale_list=[0.8, 1.0, 1.2],
        pert_num=20,
        pert_atom=0.03,
        pert_box=0.03,
        type_name=["Al"],
        save_path="Al_FCC",
        save_type="cp2k"
    )

>>> pert.compute() # check results in './Al_FCC'

compute(init_type='init_bulk', direction='z', distance=10.0)

Generate perturb models.

Parameters:

• init_type (str, optional) – selected in [“init_bulk”, “init_surf”]. Defaults to “init_bulk”.

• direction (str, optional) – the direction to generate vacuum layer. This parameter will be ignored for init_bulk. Defaults to “z”.

• distance (float, optional) – this length of vacuum layer. This parameter will be ignored for init_bulk. Defaults to 10.0.

mdapy.phonon module

mdapy.pigz module

Functions and classes to speed up compression of files by utilizing multiple cores on a system.

class mdapy.pigz.PigzFile(compression_target, output_filename=None, compresslevel=9, blocksize=128, workers=2)

Bases: object

Class to implement Pigz functionality in Python

calculate_chunk_check(chunk: bytes)

Calculate the check value for the chunk.

clean_up()

Close the output file. Clean up the processing pool.

process_compression_target()

Setup output file. Start read and write threads. Join to write thread.

write_file_trailer()

Write the trailer for the compressed data.

mdapy.pigz.compress_file(source_file, output_file=None, inplace=False, compresslevel=9, blocksize=128, workers=2)

This function provides a interface to compress a file to .gz format parallelly.

Parameters:

• source_file (str) – filename you want to compress.

• output_file (str, optional) – output compressed filename. If not give, mdapy will automatically append a .gz postfix.

• inplace (bool, optional) – whether inplace the original file. Defaults to False.

• compresslevel (int, optional) – 1 is fastest but worst, 9 is slowest but best. Defaults to 9.

• blocksize (int, optional) – blocksize, generally do not need change. Defaults to 128 KB.

• workers (int, optional) – number of threads to be used to do compression. Defaults to all your CPU cores.

mdapy.plotset module

mdapy.plotset.cm2inch(value)

Centimeters to feet.

Parameters:

value (float) – centimeters.

Returns:

feet.

Return type:

float

mdapy.plotset.pltset(color_cycler=None, **kargs)

This is a drawing template function, which refers to the science style from SciencePlots and color also comes from Paul Tot’s website.

Parameters:

color_cycler (str | list[str] | tuple[str], optional) – One can use other color cycler belongs to [‘bright’, ‘high-contrast’, ‘high-vis’, ‘light’, ‘muted’, ‘retro’, ‘std-colors’, ‘vibrant’], or directly provide a list or tuple contains colors, such as [‘#FFAABB’, ‘#99DDFF’, ‘#44BB99’]. Defaults to None.

One also can modify the default parameters by provide a dict. For example, if you want to use ‘Times New Roman’ font and hide the minor ticks: you can try:

pltset(**{‘font.serif’:’Times New Roman’, ‘xtick.minor.visible’:False, ‘ytick.minor.visible’:False})

If you want to change the axes linewidth and major tick width, and make bolded font, you can try:

pltset(**{“xtick.major.width”:1., “ytick.major.width”:1., “axes.linewidth”:1., “font.weight”:’bold’, “axes.labelweight”:’bold’})

More parameters can be found in plt.rcParams.keys().

mdapy.plotset.pltset_old(color_cycler=None, **kargs)

This function used to generate the same style with mdapy<0.9.9. Note that the ‘Times New Roman’ font is required.

Parameters:

color_cycler (str | list[str] | tuple[str], optional) – One can use other color cycler belongs to [‘bright’, ‘high-contrast’, ‘high-vis’, ‘light’, ‘muted’, ‘retro’, ‘std-colors’, ‘vibrant’], or directly provide a list or tuple contains colors, such as [‘#FFAABB’, ‘#99DDFF’, ‘#44BB99’]. Defaults to None.

One also can modify the default parameters by provide a dict. For example, if you want to let the font normal and set the font size, you can try:

pltset_old(**{“font.weight”:’normal’, “axes.labelweight”:’normal’, “font.size” : 10})

More parameters can be found in plt.rcParams.keys().

mdapy.plotset.set_figure(figsize=(10, 6), figdpi=300, top=0.935, bottom=0.18, left=0.115, right=0.965, hspace=0.2, wspace=0.2, nrow=1, ncol=1, use_pltset=False, use_pltset_old=False)

This function can generate a Figure and a Axes object easily.

Parameters:

• figsize (tuple, optional) – figsize with units of Centimeters. Defaults to (10, 6).

• figdpi (int, optional) – figure dpi to show. Defaults to 300.

• top (float, optional) – axes range in figure at top, should be [0., 1.]. Defaults to 0.935.

• bottom (float, optional) – axes range in figure at bottom, should be [0., 1.]. Defaults to 0.18.

• left (float, optional) – axes range in figure at left, should be [0., 1.]. Defaults to 0.115.

• right (float, optional) – axes range in figure at right, should be [0., 1.]. Defaults to 0.965.

• hspace (float, optional) – space between two subplots along height direction, should be [0., 1.]. Defaults to 0.2.

• wspace (float, optional) – space between two subplots along width direction, should be [0., 1.]. Defaults to 0.2.

• nrow (int, optional) – the rows number. Defaults to 1.

• ncol (int, optional) – the columns number. Defaults to 1.

• use_pltset (bool, optional) – whether use the pltset. Defaults to False.

• use_pltset_old (bool, optional) – whether use the old pltset. This will be used only if use_pltset is False. Defaults to False.

Returns:

Figure and Axes object.

Return type:

tuple

mdapy.polyhedral_template_matching module

class mdapy.polyhedral_template_matching.PolyhedralTemplateMatching(pos, box, boundary=[1, 1, 1], structure='fcc-hcp-bcc', rmsd_threshold=0.1, verlet_list=None, return_verlet=False)

Bases: object

This class identifies the local structural environment of particles using the Polyhedral Template Matching (PTM) method, which shows greater reliability than e.g. Common Neighbor Analysis (CNA). It can identify the following structure:

1. other = 0

2. fcc = 1

3. hcp = 2

4. bcc = 3

5. ico (icosahedral) = 4

6. sc (simple cubic) = 5

7. dcub (diamond cubic) = 6

8. dhex (diamond hexagonal) = 7

9. graphene = 8

Hint

If you use this class in publication, you should cite the original papar:

Larsen P M, Schmidt S, Schiøtz J. Robust structural identification via polyhedral template matching[J]. Modelling and Simulation in Materials Science and Engineering, 2016, 24(5): 055007.

Note

The present version is translated from that in LAMMPS, which is fully parallel via openmp from mdapy verison 0.8.3.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) or (\(3, 2\)) system box, must be rectangle.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• structure (str, optional) – the structure one want to identify, one can choose from [“fcc”,”hcp”,”bcc”,”ico”,”sc”,”dcub”,”dhex”,”graphene”,”all”,”default”], such as ‘fcc-hcp-bcc’. ‘default’ represents ‘fcc-hcp-bcc-ico’. Defaults to ‘fcc-hcp-bcc’.

• rmsd_threshold (float, optional) – rmsd threshold. Defaults to 0.1.

• verlet_list (np.ndarray, optional) – (\(N_p, >=18\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1. Defaults to None.

• return_verlet (bool, optional) – whether return ptm_indicis for pre-processing, if you do not need, set it to False. Defaults to False.

Outputs:

• output (np.ndarray) - (\(N_p, 7\)) the columns represent [‘type’, ‘rmsd’, ‘interatomic distance’, ‘qw’, ‘qx’, ‘qy’, ‘qz’].

• ptm_indices (np.ndarray) - (\(N_p, 18\)) ptm neighbors.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> import numpy as np

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> ptm = mp.PolyhedralTemplateMatching(FCC.pos, FCC.box) # Initilize ptm class.

>>> ptm.compute() # Compute ptm per atoms.

>>> ptm.output[:, 0] # Should be 1, that is fcc.

compute()

Do the real ptm computation.

mdapy.potential module

class mdapy.potential.BasePotential

Bases: ABC

abstract compute(pos, box, elements_list, type_list, boundary=[1, 1, 1])

Interface function.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• elements_list (list[str]) – elements to be calculated, such as [‘Al’, ‘Ni’].

• type_list (np.ndarray) – (\(N_p\)) atom type list.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

Returns:

energy, force, virial.

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray]

class mdapy.potential.EAM(filename)

Bases: BasePotential

This class is used to read/write a embedded-atom method (EAM) potentials. The energy of atom \(i\) is given by:

\[E_i = F_\alpha \left(\sum_{j \neq i}\rho_\beta (r_{ij})\right) + \frac{1}{2} \sum_{j \neq i} \phi_{\alpha\beta} (r_{ij}),\]

where \(F\) is the embedding energy, \(\rho\) is the electron density, \(\phi\) is the pair interaction, and \(\alpha\) and \(\beta\) are the elements species of atom \(i\) and \(j\). Both summation go over all neighbor atom \(j\) of atom \(i\) withing the cutoff distance.

Note

Only support eam.alloy definded in LAMMPS format now.

Parameters:

filename (str) – filename of eam.alloy file.

Outputs:

• Nelements (int) - number of elements.

• elements_list (list) - elements list.

• nrho (int) - number of \(\rho\) array.

• nr (int) - number of \(r\) array.

• drho (float) - spacing of electron density \(\rho\) space.

• dr (float) - spacing of real distance \(r\) space. Unit is Angstroms.

• rc (float) - cutoff distance. Unit is Angstroms.

• r (np.ndarray) - (nr), distance space.

• rho (np.ndarray) - (nrho), electron density space.

• aindex (np.ndarray) - (Nelements), element serial number.

• amass (np.ndarray) - (Nelements), element mass.

• lattice_constant (np.ndarray) - (Nelements), lattice constant. Unit is Angstroms.

• lattice_type (list) - (Nelements), lattice type, such as [fcc, bcc].

• embedded_data (np.ndarray) - (Nelements, nrho), embedded energy \(F\).

• elec_density_data (np.ndarray) - (Nelements, nr), electron density \(\rho\).

• rphi_data (np.ndarray) = (Nelements, Nelements, nr), \(r*\phi\).

• d_embedded_data (np.ndarray) - (Nelements, nrho), derivatives of embedded energy \(dF/d\rho\).

• d_elec_density_data (np.ndarray) - (Nelements, nr), derivatives of electron density \(d\rho/dr\).

• phi_data (np.ndarray) = (Nelements, Nelements, nr), pair interaction \(\phi\).

• d_phi_data (np.ndarray) = (Nelements, Nelements, nr), derivatives of pair interaction \(d\phi/dr\).

Examples

>>> import mdapy as mp

>>> mp.init()

>>> potential = mp.EAM("./example/CoNiFeAlCu.eam.alloy") # Read a eam.alloy file.

>>> potential.embedded_data # Check embedded energy.

>>> potential.phi_data # Check pair interaction.

>>> potential.plot() # Plot information of potential.

>>> FCC = LatticeMaker(4.05, "FCC", x, y, z)

>>> FCC.compute()

>>> Compute energy, force and virial.

>>> energy, force, virial = potential.compute(FCC.pos, FCC.box, ["Al"], np.ones(FCC.pos.shape[0], dtype=np.int32))

compute(pos, box, elements_list, type_list, boundary=[1, 1, 1], verlet_list=None, distance_list=None, neighbor_number=None)

This function is used to calculate the energy, force and virial.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• elements_list (list[str]) – elements to be calculated, such as [‘Al’, ‘Ni’].

• type_list (np.ndarray) – (\(N_p\)) atom type list.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1. Defaults to None.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom. Defaults to None.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number. Defaults to None.

Returns:

energy, force, virial.

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray]

Units & Shape:

energy : eV (\(N_p\)) force : eV/A (\(N_p, 3\)). The order is fx, fy and fz. virial : eV*A^3 (\(N_p, 9\)). The order is xx, yy, zz, xy, xz, yz, yx, zx, zy.

plot()

Plot \(F\), \(\rho\), \(\phi\).

plot_embded_rho()

Plot the \(F\) with \(\rho\).

Returns:

(fig, ax) matplotlib figure and axis class.

Return type:

tuple

plot_phi_r()

Plot the \(\phi\) with \(r\).

Returns:

(fig, ax) matplotlib figure and axis class.

Return type:

tuple

plot_rho_r()

Plot the \(\rho\) with \(r\).

Returns:

(fig, ax) matplotlib figure and axis class.

Return type:

tuple

write_eam_alloy(output_name=None)

Save to a eam.alloy file.

Parameters:

output_name (str, optional) – filename of generated eam file.

class mdapy.potential.EAMCalculator(potential, pos, boundary, box, elements_list, init_type_list, verlet_list=None, distance_list=None, neighbor_number=None)

Bases: object

This class is used to calculate the atomic energy and force based on the given embedded atom method EAM potential. Multi-elements alloy is also supported.

Parameters:

• potential (mp.EAM) – A EAM class.

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• boundary (list) – boundary conditions, 1 is periodic and 0 is free boundary. Such as [1, 1, 1].

• box (np.ndarray) – (\(3, 2\)) or (\(4, 3\)) system box.

• elements_list (list) – elements need to be calculated. Such as [‘Al’, ‘Fe’].

• init_type_list (np.ndarray) – (\(N_p\)) per atom type.

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1. Defaults to None.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom. Defaults to None.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number. Defaults to None.

Outputs:

• energy (np.ndarray) - (\(N_p\)) atomic energy (eV).

• force (np.ndarray) - (\(N_p, 3\)) atomic force (eV/A).

• virial (np.ndarray) - (\(N_p, 9\)) atomic force (eV*A^3).

Examples

>>> import mdapy as mp

>>> mp.init()

>>> potential = mp.EAM("./example/CoNiFeAlCu.eam.alloy") # Read a eam.alloy file.

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> neigh = mp.Neighbor(FCC.pos, FCC.box,
                        potential.rc, max_neigh=100) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> Cal = EAMCalculator(
        potential,
        FCC.pos,
        [1, 1, 1],
        FCC.box,
        ["Al"],
        np.ones(FCC.pos.shape[0], dtype=np.int32),
        neigh.verlet_list,
        neigh.distance_list,
        neigh.neighbor_number,
    ) # Initialize Calculator class.

>>> Cal.compute() # Calculate the atomic energy and force.

>>> Cal.energy # Check the energy.

>>> Cal.force # Check the force.

>>> Cal.virial # Check the virial.

compute()

Do the real energy, force and virial calculation.

class mdapy.potential.LammpsPotential(pair_parameter, units='metal', atomic_style='atomic', extra_args=None, conversion_factor=None)

Bases: BasePotential

This class provide a interface to use potential supported in lammps to calculate the energy, force and virial.

Parameters:

• pair_parameter (str) – including pair_style and pair_coeff.

• units (str, optional) – lammps units, such as metal, real etc. Defaults to “metal”.

• atomic_style (str, optional) – atomic_style, such as atomic, charge etc. Defaults to “atomic”.

• extra_args (str, optional) – any lammps commond. Defaults to None.

• conversion_factor (dict, optional) – units conversion. It must be {‘energy’:float, ‘force’:float, ‘virial’:float}. The float can be any number, while the key is fixed. Defaults to None.

compute(pos, box, elements_list, type_list, boundary=[1, 1, 1], centroid_stress=False)

This function is used to calculate the energy, force and virial.

# If one use NEP potential, set centroid_stress to True to obtain right 9 indices virial.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• elements_list (list[str]) – elements to be calculated, such as [‘Al’, ‘Ni’].

• type_list (np.ndarray) – (\(N_p\)) atom type list.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• centroid_stress (bool, optional) – if Ture, use compute stress/atomm. If False, use compute centroid/stress/atom. Defaults to False.

Returns:

energy, force, virial.

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray]

Shape:

energy : (\(N_p\)) force : (\(N_p, 3\)). The order is fx, fy and fz. virial : (\(N_p, 9\)). If centroid_stress is True, the order is xx, yy, zz, xy, xz, yz, yx, zx, zy. Otherwise the order is xx, yy, zz, xy, xz, yz.

to_lammps_box(box)

class mdapy.potential.NEP(filename)

Bases: BasePotential

This class is a python interface for NEP_CPU version, which can be used to evaluate the energy, force and virial of a given system.

Parameters:

filename (str) – filename of a NEP potential file, such as nep.txt.

Outputs:

• rc (float) - cutoff distance. Unit is Angstroms.

• info (dict) - information for NEP potential.

Example

>>> import mdapy as mp

>>> mp.init()

>>> FCC = LatticeMaker(4.05, "FCC", x, y, z)

>>> FCC.compute()

>>> nep = NEP("nep.txt")

>>> energy, force, virial = nep.compute(
        FCC.pos, FCC.box, ["Al"], np.ones(FCC.pos.shape[0], dtype=np.int32)
        )

>>> des = nep.get_descriptors(
            FCC.pos, FCC.box, ["Al"], np.ones(FCC.pos.shape[0], dtype=np.int32)
            ) # obtain the descriptor.

compute(pos, box, elements_list, type_list, boundary=[1, 1, 1])

This function is used to calculate the energy, force and virial.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• elements_list (list[str]) – elements to be calculated, such as [‘Al’, ‘Ni’].

• type_list (np.ndarray) – (\(N_p\)) atom type list.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

Returns:

energy, force, virial.

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray]

Units & Shape:

energy : eV (\(N_p\)) force : eV/A (\(N_p, 3\)). The order is fx, fy and fz. virial : eV*A^3 (\(N_p, 9\)). The order is xx, yy, zz, xy, xz, yz, yx, zx, zy.

fps_sample(n_sample, des_total=None, filename_list=None, elements_list=None, start_idx=None, fmt=None)

This function is used to sample the configurations using farthest point sampling method, based on the NEP descriptors. It is helpful to select the structures during active learning process.

Parameters:

• n_sample (int) – number of structures one wants to select.

• des_total (np.ndarray) – two dimensional ndarray, it actually can be any descriptors. If this parameter is given, the filename_list, elements_list and fmt will be ignored. Defaults to None.

• filename_list (list) – filename list, such as [‘0.xyz’, ‘1.xyz’]. Defaults to None.

• elements_list (list) – elements list, such as [‘Al’, ‘C’]. Defaults to None.

• start_idx (int, optional) – for deterministic results, fix the first sampled point index. Defaults to None.

• fmt (str, optional) – selected in [‘data’, ‘lmp’, ‘dump’, ‘dump.gz’, ‘poscar’, ‘xyz’, ‘cif’], One can explicitly assign the file format or mdapy will handle it with the postsuffix of filename. Defaults to None.

Returns:

selected index.

Return type:

np.ndarray

get_descriptors(pos, box, elements_list, type_list, boundary=[1, 1, 1])

This function is used to calculate the descriptor.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(4, 3\)) system box.

• elements_list (list[str]) – elements to be calculated, such as [‘Al’, ‘Ni’].

• type_list (np.ndarray) – (\(N_p\)) atom type list.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

Returns:

descriptor.

Return type:

np.ndarray

mdapy.replicate module

class mdapy.replicate.Replicate(pos, box, x=1, y=1, z=1, type_list=None)

Bases: object

This class used to replicate a position with np.ndarray format.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)), initial position to be replicated.

• box (np.ndarray) – (\(4, 3\)), initial system box.

• x (int, optional) – replication number (positive integer) along x axis. Defaults to 1.

• y (int, optional) – replication number (positive integer) along y axis. Defaults to 1.

• z (int, optional) – replication number (positive integer) along z axis. Defaults to 1.

• type_list (np.ndarray, optional) – (\(N_p\)), type for each atom. Defaults to None.

Outputs:

• box (np.ndarray) - (\(4, 3\)), replicated box.

• pos (np.ndarray) - (\(N, 3\)), replicated particle position.

• N (int) - replicated atom number.

• type_list (np.ndarray) - (\(N\)), replicated type list.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure.

>>> FCC.compute() # Get atom positions.

>>> repli = mp.Replicate(FCC.pos, FCC.box, 5, 5, 5) # Initilize a replication class.

>>> repli.compute() # Get replicated positions.

>>> repli.pos # Check replicated positions.

>>> repli.write_data() # Save to DATA file.

property N

The particle number.

compute()

Do the real replicate calculation.

write_POSCAR(output_name=None, type_name=None, reduced_pos=False)

This function writes position into a POSCAR file.

Parameters:

• output_name (str, optional) – filename of generated POSCAR file.

• type_name (list, optional) – species name. Such as [‘Al’, ‘Fe’].

• reduced_pos (bool, optional) – whether save directed coordination. Defaults to False.

write_cif(output_name=None, type_name=None)

This function writes position into a cif file.

Parameters:

• output_name (str, optional) – filename of generated cif file.

• type_name (list, optional) – species name. Such as [‘Al’, ‘Fe’].

write_data(output_name=None, data_format='atomic', num_type=None, type_name=None)

This function writes position into a DATA file.

Parameters:

• output_name (str, optional) – filename of generated DATA file.

• data_format (str, optional) – data format, selected in [‘atomic’, ‘charge’]

• num_type (int, optional) – explictly assign a number of atom type. Defaults to None.

• type_name (list, optional) – explictly assign elemantal name list, such as [‘Al’, ‘C’]. Defaults to None.

write_dump(output_name=None, compress=False)

This function writes position into a DUMP file.

Parameters:

• output_name (str, optional) – filename of generated DUMP file.

• compress (bool, optional) – whether compress the DUMP file.

write_xyz(output_name=None, type_name=None, classical=False)

This function writes position into a xyz file.

Parameters:

• output_name (str, optional) – filename of generated xyz file. Defaults to None.

• type_name (list, optional) – assign the species name. Such as [‘Al’, ‘Cu’]. Defaults to None.

• classical (bool, optional) – whether save with classical format. Defaults to False.

mdapy.spatial_binning module

class mdapy.spatial_binning.SpatialBinning(pos, direction, vbin, wbin=5.0, operation='mean')

Bases: object

This class is used to divide particles into different bins and operating on each bin. One-dimensional to Three-dimensional binning are supported.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• direction (str) – binning direction, selected in [‘x’, ‘y’, ‘z’, ‘xy’, ‘xz’, ‘yz’, ‘xyz’]

• vbin (np.ndarray) – (\(N_p, x\)), values to be operated, \(x\) means arbitrary columns.

• wbin (float, optional) – width of each bin. Defaults to 5.0.

• operation (str, optional) – operation on each bin, selected from [‘mean’, ‘sum’, ‘min’, ‘max’]. Defaults to “mean”.

Outputs:

• res (np.ndarray) - binning results, res[:, …, 0] is the number of atoms each bin.

• coor (dict) - coordination along binning direction, such as coor[‘x’].

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(4.05, 'FCC', 100, 50, 50) # Create a FCC structure.

>>> FCC.compute() # Get atom positions.

>>> binning = SpatialBinning(
        FCC.pos,
        "xz",
        FCC.pos[:, 0],
        operation="mean",
    ) # Initilize the binning class along 'xz' plane.

>>> binning.compute() # Do the binnning calculation.

>>> binning.res # Check the binning results.

>>> binning.coor['x'] # Check the x coordination.

>>> binning.plot(bar_label='x coordination')

compute()

Do the real binning calculation.

plot(value_label=None)

Plot the binning results multidimensional binning. For multi values, only the first value will be plotted.

Parameters:

value_label (str, optional) – y-lable for One-dimensional binning, or colorbar label for Two/Three-dimensional binning. Defaults to None.

Returns:

(fig, ax) matplotlib figure and axis class.

Return type:

tuple

mdapy.steinhardt_bond_orientation module

class mdapy.steinhardt_bond_orientation.SteinhardtBondOrientation(pos, box, boundary=[1, 1, 1], verlet_list=None, distance_list=None, neighbor_number=None, rc=0.0, qlist=array([4, 6, 8, 10, 12]), nnn=12, wlflag=False, wlhatflag=False, max_neigh=60, use_weight=False, weight=None, voronoi=False, average=False)

Bases: object

This class is used to calculate a set of bond-orientational order parameters \(Q_{\ell}\) to characterize the local orientational order in atomic structures. We first compute the local order parameters as averages of the spherical harmonics \(Y_{\ell m}\) for each neighbor:

\[\bar{Y}_{\ell m} = \frac{1}{nnn}\sum_{j = 1}^{nnn} Y_{\ell m}\bigl( \theta( {\bf r}_{ij} ), \phi( {\bf r}_{ij} ) \bigr),\]

where the summation goes over the \(nnn\) nearest neighbor and the \(\theta\) and the \(\phi\) are the azimuthal and polar angles. Then we can obtain a rotationally invariant non-negative amplitude by summing over all the components of degree \(l\):

\[Q_{\ell} = \sqrt{\frac{4 \pi}{2 \ell + 1} \sum_{m = -\ell }^{m = \ell } \bar{Y}_{\ell m} \bar{Y}^*_{\ell m}}.\]

For a FCC lattice with \(nnn=12\), \(Q_4 = \sqrt{\frac{7}{192}} \approx 0.19094\). More numerical values for commonly encountered high-symmetry structures are listed in Table 1 of J. Chem. Phys. 138, 044501 (2013), and all data can be reproduced by this class.

If \(wlflag\) is True, this class will compute the third-order invariants \(W_{\ell}\) for the same degrees as for the \(Q_{\ell}\) parameters:

\[\begin{split}W_{\ell} = \sum \limits_{m_1 + m_2 + m_3 = 0} \begin{pmatrix}\ell & \ell & \ell \\m_1 & m_2 & m_3\end{pmatrix}\bar{Y}_{\ell m_1} \bar{Y}_{\ell m_2} \bar{Y}_{\ell m_3}.\end{split}\]

For FCC lattice with \(nnn=12\), \(W_4 = -\sqrt{\frac{14}{143}} \left(\frac{49}{4096}\right) \pi^{-3/2} \approx -0.0006722136\).

If \(wlhatflag\) is true, the normalized third-order invariants \(\hat{W}_{\ell}\) will be computed:

\[\begin{split}\hat{W}_{\ell} = \frac{\sum \limits_{m_1 + m_2 + m_3 = 0} \begin{pmatrix}\ell & \ell & \ell \\m_1 & m_2 & m_3\end{pmatrix}\bar{Y}_{\ell m_1} \bar{Y}_{\ell m_2} \bar{Y}_{\ell m_3}}{\left(\sum \limits_{m=-l}^{l} |\bar{Y}_{\ell m}|^2 \right)^{3/2}}.\end{split}\]

For FCC lattice with \(nnn=12\), \(\hat{W}_4 = -\frac{7}{3} \sqrt{\frac{2}{429}} \approx -0.159317\). More numerical values of \(\hat{W}_{\ell}\) can be found in Table 1 of Phys. Rev. B 28, 784, and all data can be reproduced by this class.

Hint

If you use this class in your publication, you should cite the original paper:

Steinhardt P J, Nelson D R, Ronchetti M. Bond-orientational order in liquids and glasses[J]. Physical Review B, 1983, 28(2): 784.

Note

This class is translated from that in LAMMPS.

We also further implement the bond order to identify the solid or liquid state for lattice structure. For FCC structure, one can compute the normalized cross product:

\[s_\ell(i,j) = \frac{4\pi}{2\ell + 1} \frac{\sum_{m=-\ell}^{\ell} \bar{Y}_{\ell m}(i) \bar{Y}_{\ell m}^*(j)}{Q_\ell(i) Q_\ell(j)}.\]

According to J. Chem. Phys. 133, 244115 (2010), when \(s_6(i, j)\) is larger than a threshold value (typically 0.7), the bond is regarded as a solid bond. Id the number of solid bond is larger than a threshold (6-8), the atom is considered as solid phase.

Hint

If you use identifySolidLiquid function in this class in your publication, you should cite the original paper:

Filion L, Hermes M, Ni R, et al. Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling, and forward flux sampling: A comparison of simulation techniques[J]. The Journal of chemical physics, 2010, 133(24): 244115.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(3, 2\)) system box, must be rectangle.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1. Defaults to None.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom. Defaults to None.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number. Defaults to None.

• rc (float, optional) – cutoff distance to find neighbors. Defaults to 0.0.

• qlist (list|int, optional) – the list of order parameters to be computed, which should be a non-negative integer. Defaults to np.array([4, 6, 8, 10, 12], int).

• nnn (int, optional) – the number of nearest neighbors used to calculate \(Q_{\ell}\). If \(nnn > 0\), the \(rc\) has no effects, otherwise the summation will go over all neighbors within \(rc\). Defaults to 12.

• wlflag (bool, optional) – whether calculate the third-order invariants \(W_{\ell}\). Defaults to False.

• wlhatflag (bool, optional) – whether calculate the normalized third-order invariants \(\hat{W}_{\ell}\). If \(wlflag\) is False, this parameter has no effect. Defaults to False.

• max_neigh (int, optional) – a given maximum neighbor number per atoms. Defaults to 60.

• use_weight (bool, optional) – whether use weight. Defaults to False.

• weight (np.ndarray, optional) – a given weight to calculate the qlm. Such as a voronoi face area array. Defaults to None.

• voronoi (bool, optional) – indicate that whether the neighbor is Voronoi neighbor. Defaults to False.

• average (bool, optional) – whether average the qlm within the neighbor atoms. Defaults to False.

Outputs:

• qnarray (np.ndarray) - (math:N_p, len(qlist)*(1+wlflag+wlhatflag)) consider the \(qlist=[4, 6]\) and \(wlflag\) and \(wlhatflag\) is True, the columns of \(qnarray\) are [\(Q_4, Q_6, W_4, W_6, \hat{W}_4, \hat{W}_6\)].

• solidliquid (np.ndarray) - (math:N_p), 1 indicates solid state and 0 indicates liquid state.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure

>>> FCC.compute() # Get atom positions

>>> BO = SteinhardtBondOrientation(
            FCC.pos,
            FCC.box,
            [1, 1, 1],
            None,
            None,
            None,
            0.0,
            [4, 6, 8, 10, 12],
            12,
            wlflag=False,
            wlhatflag=False,
        ) # Initialize BondOrder class

>>> BO.compute() # Do the BondOrder computation.

>>> BO.qnarray[0] # Check qnarray, it should be [0.19094067 0.57452428 0.40391458 0.01285704 0.60008306].

>>> BO.identifySolidLiquid() # Identify solid/liquid state.

>>> BO.solidliquid[0] # Should be 1, that is solid.

compute()

Do the Steinhardt Bondorder calculation.

identifySolidLiquid(threshold=0.7, n_bond=7)

Identify the solid/liquid phase. Make sure you computed the 6 in qlist.

Parameters:

• threshold (float, optional) – threshold value to determine the solid bond. Defaults to 0.7.

• n_bond (int, optional) – threshold to determine the solid atoms. Defaults to 7.

mdapy.system module

class mdapy.system.MultiSystem(filename_list, unwrap=True, sorted_id=True)

Bases: list

This class is a collection of mdapy.System class and is helful to handle the atomic trajectory.

Parameters:

• filename_list (list) – ordered filename list, such as [‘melt.0.dump’, ‘melt.100.dump’].

• unwrap (bool, optional) – make atom positions do not wrap into box due to periotic boundary. Using minimum image criterion, see https://en.wikipedia.org/wiki/Periodic_boundary_conditions#Practical_implementation:_continuity_and_the_minimum_image_convention. Defaults to True.

• sorted_id (bool, optional) – sort data by atoms id. Defaults to True.

Outputs:

• pos_list (np.ndarray): (\(N_f, N_p, 3\)), \(N_f\) frames particle position.

• Nframes (int): number of frames.

cal_lindemann_parameter(only_global=False)

Calculate the Lindemann index, which is useful to distinguish the melt process and determine the melting points of nano-particles. The Lindemann index is defined as the root-mean-square bond-length fluctuation with following mathematical expression:

\[\left\langle\sigma_{i}\right\rangle=\frac{1}{N_{p}(N_{p}-1)} \sum_{j \neq i} \frac{\sqrt{\left\langle r_{i j}^{2}\right\rangle_t-\left\langle r_{i j}\right\rangle_t^{2}}}{\left\langle r_{i j}\right\rangle_t},\]

where \(N_p\) is the particle number, \(r_{ij}\) is the distance between atom \(i\) and \(j\) and brackets \(\left\langle \right\rangle_t\) represents an time average.

Note

This class is partly referred to a work on calculating the Lindemann index.

Note

This calculation is high memory requirement. One can estimate the memory by: \(2 * 8 * N_p^2 / 1024^3\) GB.

Tip

If only global lindemann index is needed, the class can be calculated in parallel. The local Lindemann index only run serially due to the dependencies between different frames. Here we use the Welford method to update the varience and mean of \(r_{ij}\).

Parameters:

only_global (bool, optional) – whether only calculate the global index. Defaults to False.

Outputs:

• The result adds a Lindemann (mdapy.LindemannParameter) class.

Tip

One can check the results by:

>>> MS = mp.MultiSystem(filename_list) # Generate a MultiSystem class.

>>> MS.cal_lindemann_parameter() # Calculate Lindemann index.

>>> MS.Lindemann.plot() # Plot the Lindemann index.

>>> MS[0].data['lindemann'] # Check Lindemann index per particles.

cal_mean_squared_displacement(mode='windows')

Calculate the mean squared displacement MSD of system, which can be used to reflect the particle diffusion trend and describe the melting process. Generally speaking, MSD is an average displacement over all windows of length \(m\) over the course of the simulation (so-called ‘windows’ mode here) and defined by:

\[MSD(m) = \frac{1}{N_{p}} \sum_{i=1}^{N_{p}} \frac{1}{N_{f}-m} \sum_{k=0}^{N_{f}-m-1} (\vec{r}_i(k+m) - \vec{r}_i(k))^2,\]

where \(r_i(t)\) is the position of particle \(i\) in frame \(t\). It is computationally extensive while using a fast Fourier transform can remarkably reduce the computation cost as described in nMoldyn - Interfacing spectroscopic experiments, molecular dynamics simulations and models for time correlation functions and discussion in StackOverflow.

Note

One can install pyfftw to accelerate the calculation, otherwise mdapy will use scipy.fft to do the Fourier transform.

Sometimes one only need the following atomic displacement (so-called ‘direct’ mode here):

\[MSD(t) = \dfrac{1}{N_p} \sum_{i=1}^{N_p} (r_i(t) - r_i(0))^2.\]

Parameters:

mode (str, optional) – ‘windows’ or ‘direct’. Defaults to “windows”.

Outputs:

• The result adds a MSD (mdapy.MeanSquaredDisplacement) class

Tip

One can check the results by:

>>> MS = mp.MultiSystem(filename_list) # Generate a MultiSystem class.

>>> MS.cal_mean_squared_displacement() # Calculate MSD.

>>> MS.MSD.plot() # Plot the MSD results.

>>> MS[0].data['msd'] # Check MSD per particles.

write_dumps(output_col=None, compress=False)

Write all data to a series of DUMP files.

Parameters:

• output_col (list, optional) – columns to be saved, such as [‘id’, ‘type’, ‘x’, ‘y’, ‘z’].

• compress (bool, optional) – whether compress the DUMP file.

class mdapy.system.System(filename=None, fmt=None, data=None, box=None, pos=None, boundary=[1, 1, 1], vel=None, type_list=None, type_name=None, sorted_id=False)

Bases: object

This class can generate a System class for rapidly accessing almost all the analysis method in mdapy.

Note

• mdapy now supports both rectangle and triclinic box from version 0.9.0.

• mdapy only supports the simplest DATA format, atomic and charge, which means like bond information will cause an error.

• We recommend you use DUMP as input file format or directly give particle positions and box.

Parameters:

• filename (str, optional) – DATA/DUMP filename. Defaults to None.

• fmt (str, optional) – selected in [‘data’, ‘lmp’, ‘dump’, ‘dump.gz’, ‘poscar’, ‘xyz’, ‘cif’], One can explicitly assign the file format or mdapy will handle it with the postsuffix of filename. Defaults to None.

• data (polars.Dataframe, optional) – all particles information. Defaults to None.

• box (np.ndarray, optional) – (\(4, 3\) or \(3, 2\)) system box. Defaults to None.

• pos (np.ndarray, optional) – (\(N_p, 3\)) particles positions. Defaults to None.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• vel (np.ndarray, optional) – (\(N_p, 3\)) particles velocities. Defaults to None.

• type_list (np.ndarray, optional) – (\(N_p\)) type per particles. Defaults to 1.

• type_name (list, optional) – one can assign the type name per type, such as [‘Al’, ‘C’], indicate type1 is Al and type2 is C. Defaults to None.

• sorted_id (bool, optional) – whether sort system data by the particle id. Defaults to False.

Note

• Mdapy supports load/save POSCAR format from version 0.9.6. We will convert the box vector to be compatiable with that defined in lammps.

Note

• Mdapy supports load/save xyz file with classical and extended format from version 0.9.6.

Examples

There are two ways to create a System class. The first is directly reading from a DUMP/DATA file generated from LAMMPS.

>>> import mdapy as mp

>>> mp.init('cpu')

>>> system = mp.System('example.dump')

One can also create a System by giving pos, box manually.

>>> import numpy as np

>>> box = np.array([[0, 100], [0, 100], [0, 100.]])

>>> pos = np.random.random((100, 3))*100

>>> system = mp.System(box=box, pos=pos)

Then one can access almost all the analysis method in mdapy with uniform API.

>>> system.cal_atomic_entropy() # calculate the atomic entropy

One can check the calculation results:

>>> system.data

And easily save it into disk with DUMP/DATA format.

>>> system.write_dump()

property N

particle number.

Returns:

particle number.

Return type:

int

atom_distance(i, j)

Calculate the distance of atom \(i\) and atom \(j\) considering the periodic boundary.

Parameters:

• i (int) – atom \(i\).

• j (int) – atom \(j\).

Returns:

distance between given two atoms.

Return type:

float

atoms_colored_by(values, vmin=None, vmax=None, cmap='rainbow')

Rendered atoms by the given values.

Parameters:

• values (str) – column names in data. Atoms will be colored base on it. A numpy.ndarray or polars.Series is also accept.

• vmin (float, optional) – color range min, if not give, use the values.min(). Defaults to None.

• vmax (float, optional) – color range max, if not give, use the values.max(). Defaults to None.

• cmap (str, optional) – colormap name from matplotlib, check https://matplotlib.org/stable/users/explain/colors/colormaps.html. Defaults to “rainbow”.

property boundary

boundary information.

Returns:

boundary information.

Return type:

list

property box

box information.

Returns:

box information.

Return type:

np.ndarray

build_neighbor(rc=5.0, max_neigh=None)

Build neighbor withing a spherical distance based on the mdapy.Neighbor class.

Parameters:

• rc (float, optional) – cutoff distance. Defaults to 5.0.

• max_neigh (int, optional) – maximum neighbor number (highly recommened to assign a number). If not given, will estimate atomatically. Default to None.

Outputs:

• verlet_list (np.ndarray) - (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1.

• distance_list (np.ndarray) - (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom.

• neighbor_number (np.ndarray) - (\(N_p\)) neighbor atoms number.

• rc (float) - cutoff distance.

build_neighbor_voronoi(ncore=-1)

Build voronoi neighbor. Only support rectangular and periodic boundary. After building, the system classs will have voro_neighbor_number, voro_verlet_list, voro_distance_list, voro_face_areas. This neighbor information is mainly used to calculate the steinhardt bond order.

Parameters:

ncore (int, optional) – parallel CPU cores, -1 means use all cores. Defaults to -1.

cal_ackland_jones_analysis()

Using Ackland Jones Analysis (AJA) method to identify the lattice structure.

The AJA method can recgonize the following structure:

1. Other

2. FCC

3. HCP

4. BCC

5. ICO

Note

If you use this module in publication, you should also cite the original paper. Ackland G J, Jones A P. Applications of local crystal structure measures in experiment and simulation[J]. Physical Review B, 2006, 73(5): 054104..

Hint

The module uses the legacy algorithm in LAMMPS.

Outputs:

• The result is added in self.data[‘aja’].

cal_atomic_entropy(rc=5.0, sigma=0.2, use_local_density=False, compute_average=False, average_rc=None, max_neigh=80)

Calculate the entropy fingerprint, which is useful to distinguish between ordered and disordered environments, including liquid and solid-like environments, or glassy and crystalline-like environments. The potential application could identificate grain boundaries or a solid cluster emerging from the melt. One of the advantages of this parameter is that no a priori information of structure is required.

This parameter for atom \(i\) is computed using the following formula:

\[s_S^i=-2\pi\rho k_B \int\limits_0^{r_m} \left [ g(r) \ln g(r) - g(r) + 1 \right ] r^2 dr,\]

where \(r\) is a distance, \(g(r)\) is the radial distribution function, and \(\rho\) is the density of the system. The \(g(r)\) computed for each atom \(i\) can be noisy and therefore it can be smoothed using:

\[g_m^i(r) = \frac{1}{4 \pi \rho r^2} \sum\limits_{j} \frac{1}{\sqrt{2 \pi \sigma^2}} e^{-(r-r_{ij})^2/(2\sigma^2)},\]

where the sum over \(j\) goes through the neighbors of atom \(i\) and \(\sigma\) is a parameter to control the smoothing. The average of the parameter over the neighbors of atom \(i\) is calculated according to:

\[\left< s_S^i \right> = \frac{\sum_j s_S^j + s_S^i}{N + 1},\]

where the sum over \(j\) goes over the neighbors of atom \(i\) and \(N\) is the number of neighbors. The average version always provides a sharper distinction between order and disorder environments.

Note

If you use this module in publication, you should also cite the original paper.

Entropy based fingerprint for local crystalline order

Note

This class uses the same algorithm with LAMMPS.

Tip

Suggestions for FCC, the \(rc = 1.4a\) and \(average\_rc = 0.9a\) and

for BCC, the \(rc = 1.8a\) and \(average\_rc = 1.2a\), where the \(a\) is the lattice constant.

Parameters:

• rc (float, optional) – cutoff distance. Defaults to 5.0.

• sigma (float, optional) – smoothing parameter. Defaults to 0.2.

• use_local_density (bool, optional) – whether use local atomic volume. Defaults to False.

• compute_average (bool, optional) – whether compute the average version. Defaults to False.

• average_rc (_type_, optional) – cutoff distance for averaging operation, if not given, it is equal to rc. This parameter should be lower than rc.

• max_neigh (int, optional) – maximum number of atom neighbor number. Defaults to 80.

Outputs:

• The entropy is added in self.data[‘atomic_entropy’].

• The averaged entropy is added in self.data[‘ave_atomic_entropy’].

cal_atomic_strain(ref, rc=5.0, affi_map='off')

This class is used to calculate the atomic shear strain. More details can be found here. https://www.ovito.org/docs/current/reference/pipelines/modifiers/atomic_strain.html

Parameters:

• ref (mdapy.System) – a reference system object.

• rc (float, optional) – cutoff distance to determine the neighbor environments. Defaults to 5. A.

• affi_map (str, optional) – selected in [‘off’, ‘ref’]. If use to ‘ref’, the current position will affine to the reference frame. Defaults to ‘off’.

Outputs:

• The result is added in self.data[‘shear_strain’].

cal_atomic_temperature(amass=None, elemental_list=None, rc=5.0, units='metal', max_neigh=None)

Calculate an average thermal temperature per atom, wchich is useful at shock simulations. The temperature of atom \(i\) is given by:

\[T_i=\sum^{N^i_{neigh}}_0 m^i_j(v_j^i -v_{COM}^i)^2/(3N^i_{neigh}k_B),\]

where \(N^i_{neigh}\) is neighbor atoms number of atom \(i\), \(m^i_j\) and \(v^i_j\) are the atomic mass and velocity of neighbor atom \(j\) of atom \(i\), \(k_B\) is the Boltzmann constant, \(v^i_{COM}\) is the center of mass COM velocity of neighbor of atom \(i\) and is given by:

\[v^i_{COM}=\frac{\sum _0^{N^i_{neigh}}m^i_jv_j^i}{\sum_0^{N^i_{neigh}} m_j^i}.\]

Here the neighbor of atom \(i\) includes itself.

Parameters:

• amass (np.ndarray, optional) – (\(N_{type}\)) atomic mass per species, such as np.array([1., 12.]).

• elemental_list (list, optional) – (\(N_{type}\)) elemental name, such as [‘H’, ‘C’].

• rc (float, optional) – cutoff distance. Defaults to 5.0.

• units (str, optional) – units selected from [‘metal’, ‘charge’]. Defaults to “metal”.

• max_neigh (int, optional) – maximum neighbor number. If not given, will estimate atomatically. Default to None.

Outputs:

• The result is added in self.data[‘atomic_temp’].

cal_bond_analysis(rc, nbins=100, max_neigh=None)

This function calculates the distribution of bond length and angle based on a given cutoff distance.

Parameters:

• rc (float) – cutoff distance.

• nbins (int, optional) – number of bins. Defaults to 100.

• max_neigh (int, optional) – maximum number of atom neighbor number. Defaults to None.

Outputs:

• The result adds a BA (mdapy.BondAnalysis) class.

Tip

One can check the results by:

>>> system.BA.plot_bond_length_distribution() # Plot the bond length distribution.

>>> system.BA.plot_bond_angle_distribution() # Plot the bond angle distribution.

>>> system.BA.bond_length_distribution # Check the data.

>>> system.BA.bond_angle_distribution # Check the data.

cal_centro_symmetry_parameter(N=12)

Compute the CentroSymmetry Parameter (CSP), which is heluful to recgonize the structure in lattice, such as FCC and BCC. The CSP is given by:

\[p_{\mathrm{CSP}} = \sum_{i=1}^{N/2}{|\mathbf{r}_i + \mathbf{r}_{i+N/2}|^2},\]

where \(r_i\) and \(r_{i+N/2}\) are two neighbor vectors from the central atom to a pair of opposite neighbor atoms. For ideal centrosymmetric crystal, the contributions of all neighbor pairs will be zero. Atomic sites within a defective crystal region, in contrast, typically have a positive CSP value.

This parameter \(N\) indicates the number of nearest neighbors that should be taken into account when computing the centrosymmetry value for an atom. Generally, it should be a positive, even integer. Note that larger number decreases the calculation speed. For FCC is 12 and BCC is 8.

Note

If you use this module in publication, you should also cite the original paper.

Kelchner C L, Plimpton S J, Hamilton J C. Dislocation nucleation and defect structure during surface indentation[J]. Physical review B, 1998, 58(17): 11085..

Hint

The CSP is calculated by the same algorithm as LAMMPS.

First calculate all \(N (N - 1) / 2\) pairs of neighbor atoms, and the summation of the \(N/2\) lowest weights is CSP values.

Parameters:

N (int, optional) – neighbor atom number considered, should be a positive and even number. Defaults to 12.

Outputs:

• The result is added in self.data[‘csp’].

cal_cluster_analysis(rc=5.0, max_neigh=80)

Divide atoms connected within a given cutoff distance into a cluster. It is helpful to recognize the reaction products or fragments under shock loading.

Note

This class use the same method as in Ovito.

Parameters:

• rc (float | dict) – cutoff distance. One can also assign multi cutoff for different elemental pair, such as {‘1-1’:1.5, ‘1-6’:1.7}. The unassigned elemental pair will default use the maximum cutoff distance. Defaults to 5.0.

• max_neigh (int, optional) – maximum number of atom neighbor number. Defaults to 80.

Returns:

cluster number.

Return type:

int

Outputs:

• The cluster id per atom is added in self.data[‘cluster_id’]

cal_common_neighbor_analysis(rc=None, max_neigh=30)

Use Common Neighbor Analysis (CNA) method to recgonize the lattice structure, based on which atoms can be divided into FCC, BCC, HCP and Other structure.

Note

If one use this module in publication, one should also cite the original paper. Stukowski, A. (2012). Structure identification methods for atomistic simulations of crystalline materials. Modelling and Simulation in Materials Science and Engineering, 20(4), 045021..

Hint

We use the same algorithm as in OVITO.

CNA method is sensitive to the given cutoff distance. The suggesting cutoff can be obtained from the following formulas:

\[r_{c}^{\mathrm{fcc}} = \frac{1}{2} \left(\frac{\sqrt{2}}{2} + 1\right) a \approx 0.8536 a,\]

\[r_{c}^{\mathrm{bcc}} = \frac{1}{2}(\sqrt{2} + 1) a \approx 1.207 a,\]

\[r_{c}^{\mathrm{hcp}} = \frac{1}{2}\left(1+\sqrt{\frac{4+2x^{2}}{3}}\right) a,\]

where \(a\) is the lattice constant and \(x=(c/a)/1.633\) and 1.633 is the ideal ratio of \(c/a\) in HCP structure.

Prof. Alexander Stukowski has improved this method using adaptive cutoff distances based on the atomic neighbor environment, which is the default method in mdapy from version 0.11.1.

The CNA method can recgonize the following structure:

1. Other

2. FCC

3. HCP

4. BCC

5. ICO

Parameters:

• rc (float, optional) – cutoff distance, if not given, will using adaptive cutoff. Defaults to None.

• max_neigh (int, optional) – maximum number of atom neighbor number. Defaults to 30.

Outputs:

• The CNA pattern per atom is added in self.data[‘cna’].

cal_common_neighbor_parameter(rc=3.0, max_neigh=30)

Use Common Neighbor Parameter (CNP) method to recgonize the lattice structure.

Note

If one use this module in publication, one should also cite the original paper. Tsuzuki H, Branicio P S, Rino J P. Structural characterization of deformed crystals by analysis of common atomic neighborhood[J]. Computer physics communications, 2007, 177(6): 518-523..

CNP method is sensitive to the given cutoff distance. The suggesting cutoff can be obtained from the following formulas:

\[r_{c}^{\mathrm{fcc}} = \frac{1}{2} \left(\frac{\sqrt{2}}{2} + 1\right) a \approx 0.8536 a,\]

\[r_{c}^{\mathrm{bcc}} = \frac{1}{2}(\sqrt{2} + 1) a \approx 1.207 a,\]

\[r_{c}^{\mathrm{hcp}} = \frac{1}{2}\left(1+\sqrt{\frac{4+2x^{2}}{3}}\right) a,\]

where \(a\) is the lattice constant and \(x=(c/a)/1.633\) and 1.633 is the ideal ratio of \(c/a\) in HCP structure.

Some typical CNP values:

• FCC : 0.0

• BCC : 0.0

• HCP : 4.4

• FCC (111) surface : 13.0

• FCC (100) surface : 26.5

• FCC dislocation core : 11.

• Isolated atom : 1000. (manually assigned by mdapy)

Parameters:

• rc (float, optional) – cutoff distance. Defaults to 3.0.

• max_neigh (int, optional) – maximum number of atom neighbor number. Defaults to 30.

Outputs:

• The CNP pattern per atom is added in self.data[‘cnp’].

cal_energy_force_virial(potential, elements_list, centroid_stress=False)

Calculate the atomic energy and force based on the given potential.

Parameters:

• potential (BasePotential) – base potential class defined in mdapy, which must including a compute method to calculate the energy, force, virial.

• elements_list (list) – elements to be calculated, such as [‘Al’, ‘Ni’] indicates setting type 1 as ‘Al’ and type 2 as ‘Ni’.

• centroid_stress (bool, optional) – Only for LammpsPotential. If Ture, use compute stress/atomm. If False, use compute centroid/stress/atom. Defaults to False.

Returns:

energy, force, virial.

Return type:

tuple[np.ndarray, np.ndarray, np.ndarray]

Units & Shape (Only for our originally supported EAM and NEP potential.):

energy : eV (\(N_p\)) force : eV/A (\(N_p, 3\)). The order is fx, fy and fz. virial : eV*A^3 (\(N_p, 9\)). The order is xx, yy, zz, xy, xz, yz, yx, zx, zy.

cal_identify_SFs_TBs(rmsd_threshold=0.1)

This function is used to identify the stacking faults (SFs) and coherent twin boundaries (TBs) in FCC structure based on the Polyhedral Template Matching (PTM). It can identify the following structure:

1. 0 = Non-hcp atoms (e.g. perfect fcc or disordered)

2. 1 = Indeterminate hcp-like (isolated hcp-like atoms, not forming a planar defect)

3. 2 = Intrinsic stacking fault (two adjacent hcp-like layers)

4. 3 = Coherent twin boundary (one hcp-like layer)

5. 4 = Multi-layer stacking fault (three or more adjacent hcp-like layers)

6. 5 = Extrinsic stacking fault

Note

This class is translated from that implementation in Ovito but optimized to be run parallely. And so-called multi-layer stacking faults maybe a combination of intrinsic stacking faults and/or twin boundary which are located on adjacent {111} plane. It can not be distiguished by the current method.

Parameters:

rmsd_threshold (float, optional) – rmsd_threshold for ptm method. Defaults to 0.1.

Outputs:

• The result is added in self.data[‘ptm’].

• The result is added in self.data[‘fault_types’].

cal_identify_diamond_structure()

This class is used to identify the Diamond structure. The results and algorithm should be the same in Ovito. More details can be found in https://www.ovito.org/manual/reference/pipelines/modifiers/identify_diamond.html .

Outputs:

• The pattern per atom is added in self.data[‘ids’].

The identified structures include:

• 0 “other”,

• 1 “cubic_diamond”,

• 2 “cubic_diamond_1st_neighbor”,

• 3 “cubic_diamond_2st_neighbor”,

• 4 “hexagonal_diamond”,

• 5 “hexagonal_diamond_1st_neighbor”,

• 6 “hexagonal_diamond_2st_neighbor”

cal_local_warren_cowley_parameter(pair_list, rc=3.0, use_voronoi=False, max_neigh=50, default_value=-10000.0)

Calculate the local sro parameter per atom:

\[WCP_{mn} = 1 - Z_{mn}/(X_n Z_{m}),\]

where \(Z_{mn}\) is the number of \(n\)-type atoms around \(m\)-type atoms, \(Z_m\) is the total number of atoms around \(m\)-type atoms, and \(X_n\) is the atomic concentration of \(n\)-type atoms in the alloy.

Parameters:

• pair_list (list[str]) – the atom pair you want to compute, such as [‘1-2’, ‘2-3’].

• rc (float, optional) – cutoff distance. Defaults to 3.0.

• use_voronoi (bool, optional) – If set to True, we will use Voronoi neighbor to calculate the results. Defaults to False.

• max_neigh (int, optional) – maximum number of atom neighbor number. Defaults to 50.

• default_value (float, optional) – defaults values for non-reference atoms, one can assign any number here. Defaults to -10000.

cal_pair_distribution(rc=5.0, nbin=200, max_neigh=80)

Calculate the radiul distribution function (RDF),which reflects the probability of finding an atom at distance r. The seperate pair-wise combinations of particle types can also be computed:

\[g(r) = c_{\alpha}^2 g_{\alpha \alpha}(r) + 2 c_{\alpha} c_{\beta} g_{\alpha \beta}(r) + c_{\beta}^2 g_{\beta \beta}(r),\]

where \(c_{\alpha}\) and \(c_{\beta}\) denote the concentration of two atom types in system and \(g_{\alpha \beta}(r)=g_{\beta \alpha}(r)\).

Parameters:

• rc (float, optional) – cutoff distance. Defaults to 5.0.

• nbin (int, optional) – number of bins. Defaults to 200.

• max_neigh (int, optional) – maximum number of atom neighbor number. Defaults to 80.

Outputs:

• The result adds a PairDistribution (mdapy.PairDistribution) class

Tip

One can check the results by:

>>> system.PairDistribution.plot() # Plot gr.

>>> system.PairDistribution.g # Check partial RDF.

>>> system.PairDistribution.g_total # Check global RDF.

cal_phono_dispersion(path, labels, potential, elements_list, symprec=1e-05, replicate=None, displacement=0.01, cutoff_radius=None)

This function can be used to calculate the phono dispersion based on Phonopy (https://phonopy.github.io/phonopy/). We support NEP and eam/alloy potential now.

Parameters:

• path (str | np.ndarray| nested list) – band path, such as ‘0. 0. 0. 0.5 0.5 0.5’, indicating two points.

• labels (str | list) – high symmetry points label, such as [“$Gamma$”, “K”, “M”, “$Gamma$”].

• potential (BasePotential) – base potential class defined in mdapy, which must including a compute method to calculate the energy, force, virial.

• elements_list (list[str]) – element list, such as [‘Al’]

• pair_style (str, optional) – pair style, selected in [‘nep’, ‘eam/alloy’]. Defaults to “eam/alloy”.

• symprec (float) – this is used to set geometric tolerance to find symmetry of crystal structure. Defaults to 1e-5.

• replicate (list, optional) – replication to pos, such as [3, 3, 3]. If not given, we will replicate it exceeding 15 A per directions. Defaults to None.

• displacement (float, optional) – displacement distance. Defaults to 0.01.

• cutoff_radius (float, optional) – Set zero to force constants outside of cutoff distance. If not given, the force constant will consider the whole supercell. This parameter will not reduce the computation cost. Defaults to None.

Outputs:

Phon : a Phonon object, which can be used to plot phonon dispersion.

cal_polyhedral_template_matching(structure='fcc-hcp-bcc', rmsd_threshold=0.1, return_rmsd=False, return_atomic_distance=False, return_wxyz=False)

This function identifies the local structural environment of particles using the Polyhedral Template Matching (PTM) method, which shows greater reliability than e.g. Common Neighbor Analysis (CNA). It can identify the following structure:

1. other = 0

2. fcc = 1

3. hcp = 2

4. bcc = 3

5. ico (icosahedral) = 4

6. sc (simple cubic) = 5

7. dcub (diamond cubic) = 6

8. dhex (diamond hexagonal) = 7

9. graphene = 8

Hint

If you use this class in publication, you should cite the original papar:

Larsen P M, Schmidt S, Schiøtz J. Robust structural identification via polyhedral template matching[J]. Modelling and Simulation in Materials Science and Engineering, 2016, 24(5): 055007.

Note

The present version is translated from that in LAMMPS and only can run serially, we will try to make it parallel.

Parameters:

• structure (str, optional) – the structure one want to identify, one can choose from [“fcc”,”hcp”,”bcc”,”ico”,”sc”,”dcub”,”dhex”,”graphene”,”all”,”default”], such as ‘fcc-hcp-bcc’. ‘default’ represents ‘fcc-hcp-bcc-ico’. Defaults to ‘fcc-hcp-bcc’.

• rmsd_threshold (float, optional) – rmsd threshold. Defaults to 0.1.

• return_rmsd (bool, optional) – whether return rmsd. Defaults to False.

• return_atomic_distance (bool, optional) – whether return interatomic distance. Defaults to False.

• return_wxyz (bool, optional) – whether return local structure orientation. Defaults to False.

Outputs:

• The result is added in self.data[‘ptm’].

cal_species_number(element_list, search_species=None, check_most=10)

This function can recgnized the species based on the atom connectivity. For atom i and atom j, if rij <= (vdwr_i + vdwr_j) * 0.6, we think two atoms are connected. Similar method can be found in OpenBabel.

Parameters:

• element_list (list) – elemental name for your system. Such as [‘C’, ‘H’, ‘O’].

• search_species (list, optional) – the molecular formula you want to find. Such as [‘H2O’, ‘CO2’, ‘Cl2’, ‘N2’].

• check_most (int, optional) – if search_species is not given, we will give the N-most species.

Returns:

search species and the corresponding number.

Return type:

dict

cal_steinhardt_bond_orientation(nnn=12, qlist=[6], rc=0.0, wlflag=False, wlhatflag=False, use_voronoi=False, use_weight=False, average=False, solidliquid=False, max_neigh=60, threshold=0.7, n_bond=7)

This function is used to calculate a set of bond-orientational order parameters \(Q_{\ell}\) to characterize the local orientational order in atomic structures. We first compute the local order parameters as averages of the spherical harmonics \(Y_{\ell m}\) for each neighbor:

\[\bar{Y}_{\ell m} = \frac{1}{nnn}\sum_{j = 1}^{nnn} Y_{\ell m}\bigl( \theta( {\bf r}_{ij} ), \phi( {\bf r}_{ij} ) \bigr),\]

where the summation goes over the \(nnn\) nearest neighbor and the \(\theta\) and the \(\phi\) are the azimuthal and polar angles. Then we can obtain a rotationally invariant non-negative amplitude by summing over all the components of degree \(l\):

\[Q_{\ell} = \sqrt{\frac{4 \pi}{2 \ell + 1} \sum_{m = -\ell }^{m = \ell } \bar{Y}_{\ell m} \bar{Y}^*_{\ell m}}.\]

For a FCC lattice with \(nnn=12\), \(Q_4 = \sqrt{\frac{7}{192}} \approx 0.19094\). More numerical values for commonly encountered high-symmetry structures are listed in Table 1 of J. Chem. Phys. 138, 044501 (2013), and all data can be reproduced by this class.

If \(wlflag\) is True, this class will compute the third-order invariants \(W_{\ell}\) for the same degrees as for the \(Q_{\ell}\) parameters:

\[\begin{split}W_{\ell} = \sum \limits_{m_1 + m_2 + m_3 = 0} \begin{pmatrix}\ell & \ell & \ell \\m_1 & m_2 & m_3\end{pmatrix}\bar{Y}_{\ell m_1} \bar{Y}_{\ell m_2} \bar{Y}_{\ell m_3}.\end{split}\]

For FCC lattice with \(nnn=12\), \(W_4 = -\sqrt{\frac{14}{143}} \left(\frac{49}{4096}\right) \pi^{-3/2} \approx -0.0006722136\).

If \(wlhatflag\) is true, the normalized third-order invariants \(\hat{W}_{\ell}\) will be computed:

\[\begin{split}\hat{W}_{\ell} = \frac{\sum \limits_{m_1 + m_2 + m_3 = 0} \begin{pmatrix}\ell & \ell & \ell \\m_1 & m_2 & m_3\end{pmatrix}\bar{Y}_{\ell m_1} \bar{Y}_{\ell m_2} \bar{Y}_{\ell m_3}}{\left(\sum \limits_{m=-l}^{l} |\bar{Y}_{\ell m}|^2 \right)^{3/2}}.\end{split}\]

For FCC lattice with \(nnn=12\), \(\hat{W}_4 = -\frac{7}{3} \sqrt{\frac{2}{429}} \approx -0.159317\). More numerical values of \(\hat{W}_{\ell}\) can be found in Table 1 of Phys. Rev. B 28, 784, and all data can be reproduced by this class.

Hint

If you use this class in your publication, you should cite the original paper:

Steinhardt P J, Nelson D R, Ronchetti M. Bond-orientational order in liquids and glasses[J]. Physical Review B, 1983, 28(2): 784.

Note

This class is translated from that in LAMMPS.

We also further implement the bond order to identify the solid or liquid state for lattice structure. For FCC structure, one can compute the normalized cross product:

\[s_\ell(i,j) = \frac{4\pi}{2\ell + 1} \frac{\sum_{m=-\ell}^{\ell} \bar{Y}_{\ell m}(i) \bar{Y}_{\ell m}^*(j)}{Q_\ell(i) Q_\ell(j)}.\]

According to J. Chem. Phys. 133, 244115 (2010), when \(s_6(i, j)\) is larger than a threshold value (typically 0.7), the bond is regarded as a solid bond. Id the number of solid bond is larger than a threshold (6-8), the atom is considered as solid phase.

Hint

If you set solidliquid is True in your publication, you should cite the original paper:

Filion L, Hermes M, Ni R, et al. Crystal nucleation of hard spheres using molecular dynamics, umbrella sampling, and forward flux sampling: A comparison of simulation techniques[J]. The Journal of chemical physics, 2010, 133(24): 244115.

Parameters:

• nnn (int, optional) – the number of nearest neighbors used to calculate \(Q_{\ell}\). If \(nnn > 0\), the \(rc\) has no effects, otherwise the summation will go over all neighbors within \(rc\). Defaults to 12.

• qlist (list, optional) – the list of order parameters to be computed, which should be a non-negative integer. Defaults to [6].

• rc (float, optional) – cutoff distance to find neighbors. Defaults to 0.0.

• wlflag (bool, optional) – whether calculate the third-order invariants \(W_{\ell}\). Defaults to False.

• wlhatflag (bool, optional) – whether calculate the normalized third-order invariants \(\hat{W}_{\ell}\). If \(wlflag\) is False, this parameter has no effect. Defaults to False.

• use_voronoi (bool, optional) – whether use Voronoi neighbor. Defaults to False.

• use_weight (bool, optional) – whether use Voronoi face area to weight the qlm. Defaults to False.

• average (bool, optional) – whether avcerage the qlm with the neighbor atoms. Defaults to False.

• solidliquid (bool, optional) – whether identify the solid/liquid phase. Defaults to False.

• max_neigh (int, optional) – a given maximum neighbor number per atoms. Defaults to 60.

• threshold (float, optional) – threshold value to determine the solid bond. Defaults to 0.7.

• n_bond (int, optional) – threshold to determine the solid atoms. Defaults to 7.

Outputs:

• The result is added in self.data[[‘ql’, ‘wl’, ‘whl’, ‘solidliquid’]].

cal_void_distribution(cell_length, out_void=False, out_name='void.dump')

This class is used to detect the void distribution in solid structure. First we divid particles into three-dimensional grid and check the its neighbors, if all neighbor grid is empty we treat this grid is void, otherwise it is a point defect. Then useing clustering all connected ‘void’ grid into an entire void. Excepting counting the number and volume of voids, this class can also output the spatial coordination of void for analyzing the void distribution.

Note

The results are sensitive to the selection of \(cell\_length\), which should be illustrated if this class is used in your publication.

Parameters:

• cell_length (float) – length of cell, larger than lattice constant is okay.

• out_void (bool, optional) – whether outputs void coordination. Defaults to False.

• out_name (str, optional) – output filename. Defaults to void.dump.

Returns:

(void_number, void_volume), number and volume of voids.

Return type:

tuple

cal_voronoi_volume(num_t=None)

This class is used to calculate the Voronoi polygon, which can be applied to estimate the atomic volume. The calculation is conducted by the voro++ package and this class only provides a wrapper. From mdapy v0.8.6, we use extended parallel voro++ to improve the performance, the implementation can be found in An extension to VORO++ for multithreaded computation of Voronoi cells.

Parameters:

num_t (int, optional) – threads number to generate Voronoi diagram. If not given, use all avilable threads.

Outputs:

• The atomic Voronoi volume is added in self.data[‘voronoi_volume’].

• The atomic Voronoi neighbor is added in self.data[“voronoi_number”].

• The atomic Voronoi cavity radius is added in self.data[“cavity_radius”].

cal_warren_cowley_parameter(rc=3.0, max_neigh=50)

This class is used to calculate the Warren Cowley parameter (WCP), which is useful to analyze the short-range order (SRO) in the 1st-nearest neighbor shell in alloy system and is given by:

\[WCP_{mn} = 1 - Z_{mn}/(X_n Z_{m}),\]

where \(Z_{mn}\) is the number of \(n\)-type atoms around \(m\)-type atoms, \(Z_m\) is the total number of atoms around \(m\)-type atoms, and \(X_n\) is the atomic concentration of \(n\)-type atoms in the alloy.

Note

If you use this class in publication, you should also cite the original paper: X‐Ray Measurement of Order in Single Crystals of Cu3Au.

Parameters:

• rc (float, optional) – cutoff distance. Defaults to 3.0.

• max_neigh (int, optional) – maximum number of atom neighbor number. Defaults to 50.

Outputs:

• The result adds a WarrenCowleyParameter (mdapy.WarrenCowleyParameter) class

Tip

One can check the results by:

>>> system.WarrenCowleyParameter.plot() # Plot WCP.

>>> system.WarrenCowleyParameter.WCP # Check WCP.

cell_opt(pair_parameter, elements_list, units='metal', atomic_style='atomic', extra_args=None, conversion_factor=None)

This function can be used to optimize box and position using lammps.

Args:

pair_parameter (str): including pair_style and pair_coeff, such as “pair_style eam/alloy

pair_coeff * * example/Al_DFT.eam.alloy Al”.

elements_list (list[str]): elements to be calculated, such as [‘Al’, ‘Ni’]. units (str, optional): lammps units, such as metal, real etc. Defaults to “metal”. atomic_style (str, optional): atomic_style, such as atomic, charge etc. Defaults to “atomic”. extra_args (str, optional): any lammps commond. Defaults to None. conversion_factor (float, optional): units conversion. Make sure converse the length units to A. Defaults to None.

Returns:

System: an optimized system.

change_filename(filename)

change the filename.

Parameters:

filename (str) – the new filename.

property data: DataFrame

check particles information.

Returns:

particles information.

Return type:

polars.Dataframe

display()

Visualize the System.

property distance_list

distance neighbor information. Each row indicates the neighbor atom’s distance.

Returns:

distance information.

Return type:

np.ndarray

property filename

obtain filename.

Returns:

filename.

Return type:

str

property fmt

obtain file format.

Returns:

file format.

Return type:

str

property global_info

obtaine global info, such as energy, virial and stress.

Returns:

global information.

Return type:

dict

minimize(elements_list, potential, fmax=1e-05, max_itre=200, volume_change=False, hydrostatic_strain=False)

This function use the fast inertial relaxation engine (FIRE) method to minimize the system, including optimizing position and box.

Parameters:

• elements_list (list) – element name, such as [‘Al’, ‘C’]

• potential (BasePotential) – a BasePotential

• fmax (float, optional) – maximum force per atom to consider as converged. Defaults to 1e-5.

• max_itre (int, optional) – maximum iteration times. Defaults to 200.

• volume_change (bool, optional) – whether change the box to optimize the pressure. Defaults to False.

• hydrostatic_strain (bool, optional) – sonstrain the cell by only allowing hydrostatic deformation. Defaults to False.

Returns:

optimized system.

Return type:

System

property neighbor_number

neighbor number information. Each row indicates the neighbor atom’s number.

Returns:

neighbor number.

Return type:

np.ndarray

orthogonal_box(N=10)

This function try to change the box to rectangular.

Parameters:

N (int, optional) – search limit. If you can’t found rectangular box, increase N. Defaults to 10.

Returns:

a new system with reactangular box. The atoms number may be changed.

Return type:

System

property pos

particle position information. Do not change it directly. If you want to modify the positions, modify the data and call self.update_pos()

Returns:

position information.

Return type:

np.ndarray

property rc

current cutoff distance.

Returns:

cutoff distance.

Return type:

float

replicate(x=1, y=1, z=1)

Replicate the system.

Parameters:

• x (int, optional) – replication number (positive integer) along x axis. Defaults to 1.

• y (int, optional) – replication number (positive integer) along y axis. Defaults to 1.

• z (int, optional) – replication number (positive integer) along z axis. Defaults to 1.

property rho

system number density.

Returns:

system number density.

Return type:

float

select(data: DataFrame)

Generate a subsystem.

Parameters:

data (polars.DataFrame) – a new dataframe. Such as system.data.filter((pl.col(‘x’)>50) & (pl.col(‘y’)<100))

Returns:

a new subsystem.

Return type:

System

spatial_binning(direction, vbin, wbin=5.0, operation='mean')

This class is used to divide particles into different bins and operating on each bin. One-dimensional to Three-dimensional binning are supported.

Parameters:

• direction (str) – binning direction, selected in [‘x’, ‘y’, ‘z’, ‘xy’, ‘xz’, ‘yz’, ‘xyz’].

• vbin (str/list) – values to be operated, such as ‘x’ or [‘vx’, ‘pe’].

• wbin (float, optional) – width of each bin. Defaults to 5.0.

• operation (str, optional) – operation on each bin, selected from [‘mean’, ‘sum’, ‘min’, ‘max’]. Defaults to “mean”.

Outputs:

• The result adds a Binning (mdapy.SpatialBinning) class

Tip

One can check the results by:

>>> system.Binning.plot() # Plot the binning results.

>>> system.Binning.res # Check the binning results.

>>> system.Binning.coor # Check the binning coordination.

update_data(data: DataFrame, update_pos=False, update_vel=False)

Provide a interface to directly update the particle information. If you are not sure, do not use it.

Parameters:

• data (pl.DataFrame) – a new dataframe

• update_pos (bool, optional) – if new data change the position, set it to True. Defaults to False.

• update_vel (bool, optional) – if new data change the velocity, set it to True. Defaults to False.

update_pos()

Call it only if you modify the positions information by modify the data.

update_vel()

Call it only if you modify the velocities information by modify the data.

property vel

particle velocity information. Do not change it directly. If you want to modify the velocities, modify the data and call self.update_vel()

Returns:

velocity information.

Return type:

np.ndarray

property verlet_list

verlet neighbor information. Each row indicates the neighbor atom’s indice.

Returns:

verlet information.

Return type:

np.ndarray

property vol

system volume.

Returns:

system volume.

Return type:

float

wrap_pos()

Wrap atom position into box considering the periodic boundary.

write_POSCAR(output_name=None, type_name=None, reduced_pos=False, selective_dynamics=False, save_velocity=False)

This function writes particles information into a POSCAR file.

Parameters:

• output_name (str, optional) – output filename. Defaults to None.

• type_name (list, optional) – species name. Such as [‘Al’, ‘Fe’].

• selective_dynamics (bool, optional) – whether do selective dynamics. Defaults to False.

• reduced_pos (bool, optional) – whether save directed coordination. Defaults to False.

• save_velocity (bool, optional) – whether save velocities information. Defaults to False.

write_cif(output_name=None, type_name=None)

This function writes particles information into a cif file.

Parameters:

• output_name (str, optional) – output filename. Defaults to None.

• type_name (list, optional) – species name. Such as [‘Al’, ‘Fe’].

write_cp2k(output_name=None, type_name=None)

This function writes particles information for cp2k calculation.

Parameters:

• output_name (str, optional) – output filename. Defaults to None.

• type_name (list, optional) – species name. Such as [‘Al’, ‘Fe’].

write_data(output_name=None, data_format='atomic', num_type=None, type_name=None)

This function writes particles information into a DATA file.

Parameters:

• output_name (str, optional) – output filename. Defaults to None.

• data_format (str, optional) – selected in [‘atomic’, ‘charge’]. Defaults to “atomic”.

• num_type (int, optional) – explictly assign a number of atom type. Defaults to None.

• type_name (list, optional) – explictly assign elemantal name list, such as [‘Al’, ‘C’]. Defaults to None.

write_dump(output_name=None, output_col=None, compress=False)

This function writes position into a DUMP file.

Parameters:

• output_name (str, optional) – filename of generated DUMP file.

• output_col (list, optional) – which columns should be saved.

• compress (bool, optional) – whether compress the DUMP file.

write_xyz(output_name=None, type_name=None, classical=False)

This function writes position into a XYZ file. Classical model only saves [type x y z] information and do not contain the box information. Extended model can includes any particles information, just like dump format.

Parameters:

• output_name (str, optional) – filename of generated XYZ file. Defaults to None.

• type_name (list, optional) – assign the species name. Such as [‘Al’, ‘Cu’]. Defaults to None.

• classical (bool, optional) – whether save with classical format. Defaults to False.

mdapy.temperature module

class mdapy.temperature.AtomicTemperature(amass, vel, atype_list, rc, verlet_list=None, distance_list=None, pos=None, box=None, boundary=None, units='metal')

Bases: object

This class is used to calculated an average thermal temperature per atom, wchich is useful at shock simulations. The temperature of atom \(i\) is given by:

\[T_i=\sum^{N^i_{neigh}}_0 m^i_j(v_j^i -v_{COM}^i)^2/(3N^i_{neigh}k_B),\]

where \(N^i_{neigh}\) is neighbor atoms number of atom \(i\), \(m^i_j\) and \(v^i_j\) are the atomic mass and velocity of neighbor atom \(j\) of atom \(i\), \(k_B\) is the Boltzmann constant, \(v^i_{COM}\) is the center of mass COM velocity of neighbor of atom \(i\) and is given by:

\[v^i_{COM}=\frac{\sum _0^{N^i_{neigh}}m^i_jv_j^i}{\sum_0^{N^i_{neigh}} m_j^i}.\]

Here the neighbor of atom \(i\) includes itself.

Parameters:

• amass (np.ndarray) – (\(N_{type}\)) atomic mass.

• vel (np.ndarray) – (\(N_p, 3\)), atomic velocity.

• atype_list (np.ndarray) – (\(N_p\)), atomic type.

• rc (float) – cutoff distance to average.

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1.

• distance_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) distance_list[i, j] means distance between i and j atom.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number.

• pos (np.ndarray, optional) – (\(N_p, 3\)) particles positions. Defaults to None.

• box (np.ndarray, optional) – (\(3, 2\)) or (\(4, 3\)) system box. Defaults to None.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Such as [1, 1, 1]. Defaults to None.

• units (str, optional) – units defined in LAMMPS, supports metal and charge. Defaults to “metal”.

Outputs:

• T (np.ndarray) - (\(N_p\)), atomic temperature.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(3.615, 'FCC', 10, 10, 10) # Create a FCC structure.

>>> FCC.compute() # Get atom positions.

>>> neigh = mp.Neighbor(FCC.pos, FCC.box,
                        5., max_neigh=50) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> def init_vel(N, T, Mass=1.0):
        # Generate random velocity at T K.
        Boltzmann_Constant = 8.617385e-5
        np.random.seed(10086)
        x1 = np.random.random(N * 3)
        x2 = np.random.random(N * 3)
        vel = (
            np.sqrt(T * Boltzmann_Constant / Mass)
            * np.sqrt(-2 * np.log(x1))
            * np.cos(2 * np.pi * x2)
        ).reshape(N, 3)
        vel -= vel.mean(axis=0)
        return vel * 100  # A/ps

>>> vel = init_vel(FCC.N, 300.0, 1.0) # Generate random velocity at 300 K.

>>> Temp = AtomicTemperature(
        np.array([1.0]),
        vel,
        np.ones(FCC.N, dtype=int),
        5.0,
        neigh.verlet_list,
        neigh.distance_list,
    ) # Initilize the temperature class.

>>> Temp.compute() # Do the temperature calculation.

>>> Temp.T.mean() # Average temperature should be close to 300 K.

compute()

Do the real temperature calculation.

mdapy.tool_function module

mdapy.tool_function.atomic_numbers = {'Ac': 89, 'Ag': 47, 'Al': 13, 'Am': 95, 'Ar': 18, 'As': 33, 'At': 85, 'Au': 79, 'B': 5, 'Ba': 56, 'Be': 4, 'Bh': 107, 'Bi': 83, 'Bk': 97, 'Br': 35, 'C': 6, 'Ca': 20, 'Cd': 48, 'Ce': 58, 'Cf': 98, 'Cl': 17, 'Cm': 96, 'Cn': 112, 'Co': 27, 'Cr': 24, 'Cs': 55, 'Cu': 29, 'Db': 105, 'Ds': 110, 'Dy': 66, 'Er': 68, 'Es': 99, 'Eu': 63, 'F': 9, 'Fe': 26, 'Fl': 114, 'Fm': 100, 'Fr': 87, 'Ga': 31, 'Gd': 64, 'Ge': 32, 'H': 1, 'He': 2, 'Hf': 72, 'Hg': 80, 'Ho': 67, 'Hs': 108, 'I': 53, 'In': 49, 'Ir': 77, 'K': 19, 'Kr': 36, 'La': 57, 'Li': 3, 'Lr': 103, 'Lu': 71, 'Lv': 116, 'Mc': 115, 'Md': 101, 'Mg': 12, 'Mn': 25, 'Mo': 42, 'Mt': 109, 'N': 7, 'Na': 11, 'Nb': 41, 'Nd': 60, 'Ne': 10, 'Nh': 113, 'Ni': 28, 'No': 102, 'Np': 93, 'O': 8, 'Og': 118, 'Os': 76, 'P': 15, 'Pa': 91, 'Pb': 82, 'Pd': 46, 'Pm': 61, 'Po': 84, 'Pr': 59, 'Pt': 78, 'Pu': 94, 'Ra': 88, 'Rb': 37, 'Re': 75, 'Rf': 104, 'Rg': 111, 'Rh': 45, 'Rn': 86, 'Ru': 44, 'S': 16, 'Sb': 51, 'Sc': 21, 'Se': 34, 'Sg': 106, 'Si': 14, 'Sm': 62, 'Sn': 50, 'Sr': 38, 'Ta': 73, 'Tb': 65, 'Tc': 43, 'Te': 52, 'Th': 90, 'Ti': 22, 'Tl': 81, 'Tm': 69, 'Ts': 117, 'U': 92, 'V': 23, 'W': 74, 'X': 0, 'Xe': 54, 'Y': 39, 'Yb': 70, 'Zn': 30, 'Zr': 40}

Van der Waals radii in [A] taken from: A cartography of the van der Waals territories S. Alvarez, Dalton Trans., 2013, 42, 8617-8636 DOI: 10.1039/C3DT50599E

mdapy.tool_function.run_cmd(command)

mdapy.tool_function.split_dump(input_file, outputdir='res', output_file_prefix=None, fix_N=True)

This function can be used to split a dump trajectory into seperate dump files, which can be read by mp.System.

Parameters:

• input_file (str) – xyz filename.

• outputdir (str, optional) – output folder. Defaults to ‘res’.

• output_file_prefix (str, optional) – output file prefix. If not given, mdapy will generate one based on the input_file. Defaults to None.

• fix_N (bool, optional) – if the whole trajectory has the same atoms, set this to True, it will improve the performance. Defaults to True.

mdapy.tool_function.split_xyz(input_file, outputdir='res', output_file_prefix=None, fix_N=True)

This function can be used to split a xyz trajectory into seperate xyz files, which can be read by mp.System.

Parameters:

• input_file (str) – xyz filename.

• outputdir (str, optional) – output folder. Defaults to ‘res’.

• output_file_prefix (str, optional) – output file prefix. If not given, mdapy will generate one based on the input_file. Defaults to None.

• fix_N (bool, optional) – if the whole trajectory has the same atoms, set this to True, it will improve the performance. Defaults to True.

mdapy.tool_function.timer(function)

Decorators function for timing.

mdapy.visualize module

class mdapy.visualize.Visualize(data, box)

Bases: object

atom_colored_by(values, vmin=None, vmax=None, cmap='rainbow')

atom_colored_by_atom_type()

atom_colored_by_atom_type_name()

atom_colored_by_structure_type(method, show_label=False)

box2lines(box)

close()

delete_color_bar()

display()

hide_object_by_group_name(name, remove=False)

init_color()

init_plot(vertices, indices)

init_radius()

show_object_by_group_name(name)

mdapy.visualize.value2color(colors_rgb: NdarrayType(dtype=<taichi.lang.matrix.VectorType object at 0x7f785cbbc070>, ndim=None, layout=Layout.AOS, needs_grad=None), value: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None), vmin: float, vmax: float, colors: NdarrayType(dtype=None, ndim=None, layout=Layout.AOS, needs_grad=None))

mdapy.void_distribution module

class mdapy.void_distribution.VoidDistribution(pos, box, cell_length, boundary=[1, 1, 1], out_void=False, out_name='void.dump')

Bases: object

This class is used to detect the void distribution in solid structure. First we divid particles into three-dimensional grid and check the its neighbors, if all neighbor grid is empty we treat this grid is void, otherwise it is a point defect. Then useing clustering all connected ‘void’ grid into an entire void. Excepting counting the number and volume of voids, this class can also output the spatial coordination of void for analyzing the void distribution.

Note

The results are sensitive to the selection of \(cell\_length\), which should be illustrated if this class is used in your publication.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(3, 2\)) system box.

• cell_length (float) – length of cell, larger than lattice constant is okay.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Defaults to [1, 1, 1].

• out_void (bool, optional) – whether outputs void coordination. Defaults to False.

• out_name (str, optional) – filename of generated void DUMP file. Defaults to “void.dump”.

Outputs:

• void_number (int) - total number of voids.

• void_volume (float) - total volume of voids

Examples

>>> import mdapy as mp

>>> mp.init()

>>> import numpy as np

>>> FCC = mp.LatticeMaker(4.05, 'FCC', 50, 50, 50) # Create a FCC structure.

>>> FCC.compute() # Get atom positions.

Generate four voids.

>>> pos = FCC.pos.copy()

>>> pos = pos[np.sum(np.square(pos - np.array([50, 50, 50])), axis=1) > 100]

>>> pos = pos[np.sum(np.square(pos - np.array([100, 100, 100])), axis=1) > 100]

>>> pos = pos[np.sum(np.square(pos - np.array([150, 150, 150])), axis=1) > 400]

>>> pos = pos[np.sum(np.square(pos - np.array([50, 150, 50])), axis=1) > 400]

>>> void = mp.VoidDistribution(pos, FCC.box, 5., out_void=True) # Initilize void class.

>>> void.compute() # Calculated the voids and generate a file named void.dump.

>>> void.void_number # Check the void number, should be 4 here.

>>> void.void_volume # Check the void volume.

compute()

Do the real void calculation.

mdapy.voronoi_analysis module

class mdapy.voronoi_analysis.VoronoiAnalysis(pos, box, boundary=[1, 1, 1], num_t=None)

Bases: object

This class is used to calculate the Voronoi polygon, which can be applied to estimate the atomic volume. The calculation is conducted by the voro++ package and this class only provides a wrapper. From mdapy v0.8.6, we use extended parallel voro++ to improve the performance, the implementation can be found in An extension to VORO++ for multithreaded computation of Voronoi cells.

Parameters:

• pos (np.ndarray) – (\(N_p, 3\)) particles positions.

• box (np.ndarray) – (\(3, 2\)) system box.

• boundary (list) – boundary conditions, 1 is periodic and 0 is free boundary, default to [1, 1, 1].

• num_t (int, optional) – threads number to generate Voronoi diagram. If not given, use all avilable threads.

Outputs:

• vol (np.ndarray) - (\(N_p\)), atom Voronoi volume.

• neighbor_number (np.ndarray) - (\(N_p\)), atom Voronoi neighbor number.

• cavity_radius (np.ndarray) - the distance from the particle to the farthest vertex of its Voronoi cell.

Examples

>>> import mdapy as mp

>>> mp.init()

>>> FCC = mp.LatticeMaker(4.05, 'FCC', 10, 10, 10) # Create a FCC structure.

>>> FCC.compute() # Get atom positions.

>>> avol = mp.VoronoiAnalysis(FCC.pos, FCC.box, [1, 1, 1]) # Initilize the Voronoi class.

>>> avol.compute() # Calculate the Voronoi volume.

>>> avol.vol # Check atomic Voronoi volume.

>>> avol.neighbor_number # Check neighbor number.

>>> avol.cavity_radius # Check the cavity radius.

compute()

Do the real Voronoi volume calculation.

mdapy.warren_cowley_parameter module

class mdapy.warren_cowley_parameter.WarrenCowleyParameter(type_list, verlet_list=None, neighbor_number=None, rc=None, pos=None, box=None, boundary=None)

Bases: object

This class is used to calculate the Warren Cowley parameter (WCP), which is useful to analyze the short-range order (SRO) in the 1st-nearest neighbor shell in alloy system and is given by:

\[WCP_{mn} = 1 - Z_{mn}/(X_n Z_{m}),\]

where \(Z_{mn}\) is the number of \(n\)-type atoms around \(m\)-type atoms, \(Z_m\) is the total number of atoms around \(m\)-type atoms, and \(X_n\) is the atomic concentration of \(n\)-type atoms in the alloy.

Note

If you use this class in publication, you should also cite the original paper: X‐Ray Measurement of Order in Single Crystals of Cu3Au.

Parameters:

• type_list (np.ndarray) – (\(N_p\)) atom type.

• verlet_list (np.ndarray, optional) – (\(N_p, max\_neigh\)) verlet_list[i, j] means j atom is a neighbor of i atom if j > -1.

• neighbor_number (np.ndarray, optional) – (\(N_p\)) neighbor atoms number.

• pos (np.ndarray, optional) – (\(N_p, 3\)) particles positions. Defaults to None.

• box (np.ndarray, optional) – (\(3, 2\)) or (\(4, 3\)) system box. Defaults to None.

• boundary (list, optional) – boundary conditions, 1 is periodic and 0 is free boundary. Such as [1, 1, 1]. Defaults to None.

Outputs:

• WCP (np.ndarray) - (\(N_{type}, N_{type}\)) WCP for all pair elements

Examples

>>> import mdapy as mp

>>> mp.init()

>>> system = mp.System("./example/CoCuFeNiPd-4M.dump") # Read a alloy DUMP file.

>>> neigh = mp.Neighbor(system.pos, system.box, 3.0, max_neigh=30) # Initialize Neighbor class.

>>> neigh.compute() # Calculate particle neighbor information.

>>> wcp = mp.WarrenCowleyParameter(
    system.data["type"].values, neigh.verlet_list, neigh.neighbor_number
    ) # Initilize the WCP class.

>>> wcp.compute() # Calculate the WCP.

>>> wcp.WCP # Check the WCP matrix.

This should get the same results with that in paper: Simultaneously enhancing the ultimate strength and ductility of high-entropy alloys via short-range ordering.

>>> wcp.plot(["Co", "Cu", "Fe", "Ni", "Pd"]) # Plot the results.

compute()

Do the real WCP calculation.

plot(elements_list=None, vmin=-2, vmax=1, cmap='GnBu')

Plot the WCP matrix.

Parameters:

• elements_list (list, optional) – elements list, such as [‘Al’, ‘Ni’, ‘Co’].

• vmin (int, optional) – vmin for colorbar. Defaults to -2.

• vmax (int, optional) – vmax for colorbar. Defaults to 1.

• cmap (str, optional) – cmap name. Defaults to “GnBu”.

Returns:

(fig, ax) matplotlib figure and axis class.

Return type:

tuple

Module contents

mdapy.init(arch='cpu', cpu_max_num_threads=-1, offline_cache=False, debug=False, device_memory_fraction=0.9, kernel_profiler=False)

Initilize the mdapy calculation. One should call this function after import mdapy. This is a simple wrapper function of taichi.init().

Parameters:

• arch (str, optional) – run on CPU or GPU. Defaults to “cpu”, choose in ‘cpu’ and ‘gpu’.

• cpu_max_num_threads (int, optional) – maximum CPU core to use in calculation. Defaults to -1, indicating using all available CPU cores.

• offline_cache (bool, optional) – whether save compile cache. Defaults to False.

• debug (bool, optional) – whether use debug mode. Defaults to False.

• device_memory_fraction (float, optional) – preallocate available GPU memory fraction. Defaults to 90%.

• kernel_profiler (bool, optional) – whether enable profiler. Defaults to False.

Raises:

ValueError – Unrecognized arch, please choose in [‘cpu’, ‘gpu’].