Skip to content

Oxygen Displacement

This tutorial is used to calculate the displacement of each oxygen atom in HfO₂ systems relative to the geometric center of its neighboring regular tetrahedron, and to determine whether the structure has global polarization.

Code

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
from ase import Atoms
from ase.io import read
import numpy as np


def signed_volume(a, b, c, d):
    """Calculate the signed volume of a tetrahedron defined by points a, b, c, d"""
    return np.dot(np.cross(b - a, c - a), d - a) / 6.0


def is_point_in_tetrahedron(P, A, B, C, D, tol=0.02):
    """Determine if point P is inside the tetrahedron defined by points A, B, C, D"""
    # Calculate signed volumes of sub-tetrahedra
    V0 = signed_volume(A, B, C, D)
    V1 = signed_volume(P, B, C, D)
    V2 = signed_volume(A, P, C, D)
    V3 = signed_volume(A, B, P, D)
    V4 = signed_volume(A, B, C, P)

    # Consistent signs → point is inside (allowing ±tol error)
    signs = np.array([V0, V1, V2, V3, V4])
    positive = signs > -tol
    negative = signs <  tol
    return np.all(positive) or np.all(negative)


# ---------- Main program ----------

# Read POSCAR
atoms = read("POSCAR")

# Get indices of Hf / O separately
Hf_indices = [i for i, at in enumerate(atoms) if at.symbol == "Hf"]
O_indices  = [i for i, at in enumerate(atoms) if at.symbol == "O"]

# Generate "nth atom of same element" numbering for each atom
element_specific_indices = {}
counter = {}
for i, at in enumerate(atoms):
    counter[at.symbol] = counter.get(at.symbol, 0) + 1
    element_specific_indices[i] = counter[at.symbol]

sum_dx = sum_dy = sum_dz = 0.0
valid_O_count = 0
O_displacements = {}

with open("B4ID.dat", "w") as f_b4id, open("O-disp.dat", "w") as f_disp:

    for i in O_indices:
        O_pos = atoms.get_positions()[i]

        # O-Hf distances (considering periodicity)
        distances = [(j, atoms.get_distance(i, j, mic=True)) for j in Hf_indices]
        distances.sort(key=lambda x: x[1])

        first_three = [idx for idx, _ in distances[:3]]
        found = False

        # Try 4th Hf in sequence
        for fourth in (idx for idx, _ in distances[3:]):
            cand = first_three + [fourth]

            # Get minimum image coordinates of 4 Hf atoms
            Hf_pos = [
                O_pos + atoms.get_distance(i, idx, mic=True, vector=True)
                for idx in cand
            ]

            # All four edges must be < 5 Å
            edges = [
                atoms.get_distance(cand[m], cand[n], mic=True)
                for m in range(4) for n in range(m + 1, 4)
            ]
            if np.all(np.array(edges) < 5.0):
                # Write to B4ID.dat
                f_b4id.write(" ".join(map(str, cand)) + "\n")
                found = True

                # Tetrahedron centroid
                center = np.mean(np.array(Hf_pos), axis=0)

                # Displacement Δ⃗
                dx, dy, dz = O_pos - center
                f_disp.write(f"{dx:.6f} {dy:.6f} {dz:.6f}\n")

                O_displacements[i] = (dx, dy, dz)
                sum_dx += dx; sum_dy += dy; sum_dz += dz
                valid_O_count += 1
                break

        if not found:
            # Not coordinated
            num = element_specific_indices[i]
            f_b4id.write(f"O atom {num} did not find suitable 4th Hf atom\n")
            O_displacements[i] = (0.0, 0.0, 0.0)


# Average displacement
if valid_O_count:
    avg_dx = sum_dx / valid_O_count
    avg_dy = sum_dy / valid_O_count
    avg_dz = sum_dz / valid_O_count
    with open("avgdisp.dat", "w") as f:
        f.write(f"{avg_dx:.6f} {avg_dy:.6f} {avg_dz:.6f}\n")
    print(f"Average displacement: ({avg_dx:.6f}, {avg_dy:.6f}, {avg_dz:.6f}) Å")
else:
    print("No valid O atom tetrahedra found.")


# Generate POSCAR.xsf (with displacement vectors)
with open("POSCAR.xsf", "w") as f:
    f.write("CRYSTAL\nPRIMVEC\n")
    for v in atoms.get_cell():
        f.write(f" {v[0]:.16f} {v[1]:.16f} {v[2]:.16f}\n")

    f.write("PRIMCOORD\n")
    f.write(f"{len(atoms)}\n")

    for i, at in enumerate(atoms):
        x, y, z = atoms.get_positions()[i]
        if at.symbol == "Hf":
            f.write(f"Hf {x:.8f} {y:.8f} {z:.8f} 0.00000000 0.00000000 0.00000000\n")
        else:
            dx, dy, dz = O_displacements[i]
            f.write(f"O  {x:.8f} {y:.8f} {z:.8f} {dx:.6f} {dy:.6f} {dz:.6f}\n")

Tetrahedron-Displacement Algorithm Explanation

1. Determining "whether a point is inside a tetrahedron" — Signed Volume Criterion

Given a tetrahedron formed by four vertices A B C D, its signed volume is: 

1
V₀ = [(B−A) × (C−A)] · (D−A) / 6

For any test point P, construct 4 "smaller" tetrahedra: 

1
2
V₁ = V(P,B,C,D)   V₂ = V(A,P,C,D)
V₃ = V(A,B,P,D)   V₄ = V(A,B,C,P)
  • If V₁ … V₄ have consistent signs with V₀, P is inside;
  • If some Vᵢ≈0, then P is on a face/edge;
  • If signs are not all the same, P is outside.

The code uses a ±tol tolerance, considering five volumes to have the same sign if all are > −tol or < tol.

2. Finding 4 Coordinating Hf Atoms for O

  1. For each O atom, calculate minimum image distances to all Hf atoms;
  2. Sort distances in ascending order, always use the nearest 3;
  3. Try the 4th nearest Hf in sequence, accept only when all 6 edges of the 4 Hf atoms are < 5 Å;
  4. Centroid: 
1
r₀ = (r₁ + r₂ + r₃ + r₄) / 4

Displacement (polarization vector): 

1
Δ⃗ = r(O) − r₀

Write to O-disp.dat. Record as "not coordinated" if no suitable fourth vertex is found.

3. Interactive 3D Visualization

Save the following complete HTML as view_tetra.html to view in a browser:

  • Left dropdown menu can select any O atom
  • Scene elements:
  • Gray spheres: all Hf
  • Pink spheres: all O
  • Red sphere: selected O
  • Gold spheres: 4 Hf atoms forming the tetrahedron
  • Yellow lines: 6 edges of tetrahedron
  • Green sphere: centroid
  • Blue arrow: polarization vector Δ⃗
  • Supports drag rotation and zoom

🔗 Click to view interactive demo

The above content provides both algorithm explanation and visualization examples for further debugging or demonstration.