Analysis of Atomic Configuration Data of Copper Crystals¶

In this tutorial, we will analyze atomic configuration data of copper crystals.

Preparation¶

To perform the analysis, you'll need ASE and py3Dmol. You can install them using pip as follows:

pip3 install ase py3Dmol

Loading the Data¶

Let's analyze a copper crystal with a face-centered cubic (FCC) lattice structure.
We'll use data downloaded from the Materials Project.

The FCC copper crystal is identified by the ID mp-30, and detailed information is available at:
https://materialsproject.org/materials/mp-30

For the terms of use, please refer to:
https://materialsproject.org/about/terms
The data is licensed under CC-BY, which allows relatively flexible and open use.

Import the necessary libraries.

In [2]:
import ase
import ase.io
import homcloud.interface as hc
import numpy as np
import py3Dmol
import io
import homcloud.py3dmolhelper as py3dmolhelper

We'll use ase.io.read to load the CIF file.
Make sure the file is placed in the same directory beforehand.

In [3]:
cu_atoms = ase.io.read("Cu.cif")

By inspecting the atomic positions as shown below, we can see that the structure consists of four atoms.

In [4]:
cu_atoms.get_positions()
Out[4]:
array([[0.        , 0.        , 0.        ],
       [0.        , 1.78871534, 1.78871534],
       [1.78871534, 0.        , 1.78871534],
       [1.78871534, 1.78871534, 0.        ]])

Since analyzing only four atoms with HomCloud doesn't work well, we'll use a dataset in which the structure is repeated four times along each of the X, Y, and Z directions.

In [5]:
cu_atoms_4 = cu_atoms.repeat(4)

Let's visualize this data using py3Dmol. Since py3Dmol requires input in xyz format, we'll first convert the data accordingly. Here, we'll use Python's StringIO class to handle the conversion.

In [6]:
# Convert the data to xyz format
s = io.StringIO()
ase.io.write(s, cu_atoms_4, format="xyz")
xyz_data = s.getvalue()

# Setup the viewer
view = py3Dmol.view()
# Add the atomic data to the viewer
view.addModel(xyz_data, "xyz")
# Set the rendering style for atomic data
view.setStyle({'stick': {"colorscheme": "default"}})
view.addStyle({"sphere": {"scale": 0.3}})
# Display
view.zoomTo()
view.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

We'll extract the atomic coordinates using the get_positions method.

In [7]:
atoms_positions = cu_atoms_4.get_positions()

In HomCloud analysis, perfect periodic structures like crystals can cause problems.
The computed persistence diagram may become unreliable in cases such as:

  1. Three points lying on the same straight line
  2. Four points lying on the same plane
  3. Five points lying on the same sphere

To address this, we add a small amount of noise to the atomic positions.
Since the typical interatomic distance is about 2.53 Å, adding noise on the order of $10^{-4}$ will not significantly affect the analysis results.

In [8]:
atoms_positions += np.random.uniform(-1e-4, 1e-4,size=atoms_positions.shape)

We will compute the PDs and display both the 1-dimensional and 2-dimensional PDs. Since the units of both the X and Y axes are in square angstroms, we will include that information in the display.

In [9]:
hc.PDList.from_alpha_filtration(atoms_positions, save_boundary_map=True, save_to="cu.pdgm")
Out[9]:
PDList(path=cu.pdgm)
In [10]:
pdlist = hc.PDList("cu.pdgm")
pd1 = pdlist.dth_diagram(1)
pd1.histogram((1, 4), 32).plot(colorbar={"type": "log"}, title="1D", unit_name="$\\AA^2$")
pd2 = pdlist.dth_diagram(2)
pd2.histogram((1, 4), 32).plot(colorbar={"type": "log"}, title="2D", unit_name="$\\AA^2$")
No description has been provided for this image
No description has been provided for this image

Analysis of the 1-Dimensional Persistence Diagram¶

In the 1-dimensional persistence diagram, birth-death pairs that are significantly distant from the diagonal are all concentrated around the point (1.6, 2.1).
Therefore, we will extract these distinctive pairs for further analysis.

In [11]:
pairs = pd1.pairs_in_rectangle(1.5, 1.8, 2.0, 2.2)
In [12]:
print(len(pairs))
print(pairs[0])

921
Pair(1.5996103698063915, 2.1327238876468386)

We found that there are 921 pairs in total. Let's randomly select 10 of them and compute their stable volumes.

In [13]:
import random
stable_volumes = [
  pair.stable_volume(1e-3)
  for pair in random.sample(pairs, 10)
]

The result will be visualized.

In [14]:
# Prepare the viewer
view = py3Dmol.view()
# Add atomic configuration data to the viewer
view.addModel(xyz_data, "xyz")
# Set the rendering style for atomic data
view.setStyle({'stick': {"colorscheme": "default"}})
view.addStyle({"sphere": {"scale": 0.3}})
# Add stable volume information to the viewer
for stable_volume in stable_volumes:
    py3dmolhelper.add_edges(view, stable_volume.boundary(), "green", 0.4)

# Display the result
view.zoomTo()
view.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

This suggests that the shapes corresponding to these pairs are related to equilateral triangles in the face-centered cubic (FCC) lattice. To verify this, let's calculate the distances between the vertices of the rings obtained from the stable volumes. We'll use scipy's distance_matrix function to compute the distance matrix.

In [15]:
triangle = stable_volumes[0].boundary_points()
print(triangle)

[[10.732248193410797, 10.732372233083277, 7.1548608145419985], [10.732310908754183, 8.943489357587952, 8.943540459221547], [12.520933523409571, 10.732234429977877, 8.943653709043415]]
In [16]:
import scipy as sp
sp.spatial.distance_matrix(triangle, triangle)
Out[16]:
array([[0.        , 2.52971872, 2.52965912],
       [2.52971872, 0.        , 2.52958095],
       [2.52965912, 2.52958095, 0.        ]])

It is found that the distances between the vertices are all equal, measuring approximately 2.529.... This value is consistent with the known interatomic distance of 2.53 Å.

Let's extract only the atoms that form the triangle. Using boundary_points, you can obtain the coordinates of the vertices, but I will also explain how to retrieve other related information.
By specifying by="vindexes" in boundary_points, you can find out which atoms correspond to each vertex. Note that the indexing starts from 0.

In [17]:
stable_volumes[0].boundary_points(by="vindexes")
Out[17]:
[248, 233, 250]

By passing this array of indices to the ASE Atoms object, you can extract only the corresponding atoms. Using get_positions and get_chemical_symbols, you can retrieve the positions and types of the atoms. In this case, all atoms are Cu, so it may not be very interesting, but this approach can be quite useful when dealing with more complex materials.

In [18]:
triangle_atoms = cu_atoms_4[stable_volumes[0].boundary_points(by="vindexes")]
In [19]:
triangle_atoms.get_positions()
Out[19]:
array([[10.73229201, 10.73229201,  7.15486134],
       [10.73229201,  8.94357667,  8.94357667],
       [12.52100735, 10.73229201,  8.94357667]])
In [20]:
triangle_atoms.get_chemical_symbols()
Out[20]:
['Cu', 'Cu', 'Cu']

get_all_distances allows you to calculate interatomic distances. It returns the same results as sp.spatial.distance_matrix, but if you're familiar with ASE, this method might be more convenient.

In [21]:
triangle_atoms.get_all_distances()
Out[21]:
array([[0.        , 2.52962549, 2.52962549],
       [2.52962549, 0.        , 2.52962549],
       [2.52962549, 2.52962549, 0.        ]])

Let's convert the triangle data into XYZ format and save it as a file. Here's how to do it.

In [22]:
ase.io.write("cu_triangle_atoms.xyz", triangle_atoms)

he data has been saved in a file named cu_triangle_atoms.xyz, so you can analyze or visualize it using your preferred software, such as VESTA. Now, let's read and display the contents of the XYZ file.

In [23]:
with open("cu_triangle_atoms.xyz") as f:
    cu_triangle_atoms_xyz_data = f.read()

print(cu_triangle_atoms_xyz_data)

3
Lattice="14.30972268 0.0 0.0 0.0 14.30972268 0.0 0.0 0.0 14.30972268" Properties=species:S:1:pos:R:3:spacegroup_kinds:I:1 spacegroup="P 1" unit_cell=conventional occupancy="_JSON {\"0\": {\"Cu\": 1}, \"1\": {\"Cu\": 1}, \"2\": {\"Cu\": 1}, \"3\": {\"Cu\": 1}}" pbc="T T T"
Cu      10.73229201      10.73229201       7.15486134        0
Cu      10.73229201       8.94357667       8.94357667        1
Cu      12.52100735      10.73229201       8.94357667        2

We will pass this data to py3Dmol and render only the triangle.

In [24]:
view = py3Dmol.view()
view.addModel(cu_triangle_atoms_xyz_data, "xyz")
view.setStyle({'stick': {"colorscheme": "default"}})
view.addStyle({"sphere": {"scale": 0.3}})
view.zoomTo()
view.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Analysis of the 2D Persistence Diagram¶

Next, we will analyze the 2D persistence diagram. In this diagram, pairs are concentrated in two regions: around (2.1, 2.4) and around (2.1, 3.2).

We will begin by extracting the pairs located near (2.1, 2.4).．

In [25]:
pairs = pd2.pairs_in_rectangle(2.1, 2.3, 2.3, 2.5)
print(len(pairs)) # Print the number of pairs
print(pairs[0]) # Print the birth and death times

343
Pair(2.1328119801783947, 2.399367736598198)

As in the case of the 1D persistence diagram, let's randomly select 10 pairs and compute and visualize their stable volumes.

In [26]:
# Stable volumes are computed for ten random pairs
stable_volumes = [pair.stable_volume(1e-3) for pair in random.sample(pairs, 10)]
# 3D Visualization

view = py3Dmol.view()
view.addModel(xyz_data, "xyz")
view.setStyle({'stick': {"colorscheme": "default", "opacity": 0.9}}) # 棒の部分を半透明にする
view.addStyle({"sphere": {"scale": 0.3}})
# stable volumeの情報を青色でビューアに追加
for stable_volume in stable_volumes:
    py3dmolhelper.add_surface(view, stable_volume.boundary(), "blue") 

# 以下表示
view.zoomTo()
view.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

This reveals that the structure forms a regular tetrahedron. As with the triangle case, we will convert the tetrahedral data into XYZ format and render only that part.
Here, we also incorporate the following two improvements:

  • Encapsulating the process in a function to make it reusable
  • Using io.StringIO for a more elegant implementation
In [27]:
# Define a function
def show_atoms(atoms):
    f = io.StringIO()
    ase.io.write(f, atoms, format="xyz")

    view = py3Dmol.view()
    view.addModel(f.getvalue(), "xyz")
    view.setStyle({'stick': {"colorscheme": "default"}})
    view.addStyle({"sphere": {"scale": 0.3}})
    view.zoomTo()
    view.show()
In [28]:
show_atoms(cu_atoms_4[stable_volumes[0].boundary_points(by="vindexes")])

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

Next, let's examine the pairs located near (2.1, 3.2).

In [29]:
pairs = pd2.pairs_in_rectangle(2.0, 2.2, 3.1, 3.3)
print(len(pairs)) # Print the number of pairs
print(pairs[0]) # Print the birth and death times

108
Pair(2.1329866829461572, 3.1991477227674245)
In [30]:
stable_volumes = [pair.stable_volume(1e-3) for pair in random.sample(pairs, 10)]
# 3D Visualization
view = py3Dmol.view()
view.addModel(xyz_data, "xyz")
view.setStyle({'stick': {"colorscheme": "default", "opacity": 0.9}}) # Make the bond representations semi-transparent.
view.addStyle({"sphere": {"scale": 0.3}})
# Add the stable volume information to the viewer in blue
for stable_volume in stable_volumes:
    py3dmolhelper.add_surface(view, stable_volume.boundary(), "blue") 
# Display
view.zoomTo()
view.show()

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

This indicates that the structure corresponding to these pairs resembles a regular octahedron.

When extracting just one, it looks like as follows. We use the function defined earlier.

In [31]:
show_atoms(cu_atoms_4[stable_volumes[0].boundary_points(by="vindexes")])

3Dmol.js failed to load for some reason. Please check your browser console for error messages.

In other words, the two types of birth-death pairs that appear in the 2D PD correspond to the voids formed by tetrahedral and octahedral structures.

In fact, it is known that a face-centered cubic (FCC) lattice contains two types of voids—tetrahedral and octahedral sites—in a 2:1 ratio. Therefore, it is reasonable to assume that the two types of pairs observed in the PD reflect this structural feature.

However, the actual ratio of the number of pairs obtained is 343 to 108, which is approximately 3:1, deviating from the theoretical value. Why does this discrepancy occur?

This deviation arises because the periodic structure of the lattice is not fully taken into account at the boundaries of the data.

There are two main ways to address this issue:

  • Compute the PD using periodic boundary conditions: This is the ideal solution, but it can only be applied when the unit cell is a rectangular prism. Additionally, using periodic boundaries may disable certain features, such as those provided by PHTrees.

  • Increase the number of repetitions of the lattice: This is a more general approach, but it increases computational cost as the number of repetitions grows.

Let's check the lattice information.

In [32]:
cu_atoms_4.cell.angles()
Out[32]:
array([90., 90., 90.])

Since all angles between the lattice vectors are 90 degrees, the unit cell is rectangular. Therefore, we should adopt the first approach.

In [33]:
a, b, c = cu_atoms_4.cell.lengths()
hc.PDList.from_alpha_filtration(atoms_positions, periodicity=((0, a), (0, b), (0, c)), save_to="cu_8_pbc.pdgm")
pd2 = hc.PDList("cu_8_pbc.pdgm").dth_diagram(2)
len(pd2.pairs_in_rectangle(2.1, 2.3, 2.3, 2.5)) / len(pd2.pairs_in_rectangle(2.0, 2.2, 3.1, 3.3))
Out[33]:
2.007843137254902

There is a slight discrepancy, but the ratio is still approximately 2:1.

This concludes the current tutorial. In the next tutorial, we will analyze amorphous SiO₂, which has a more complex structure.

Exercise¶

Analyze the birth-death pairs located near the diagonal in the PD of a copper crystal.