Inconsistency in oriented_unit_cell when generating (001)/(010)/(100) surface

[hezhengda]

Bug description

When I generated the slab of TiNb(010) surface, I have found that the thickness of the slab (6.95Å) is smaller than min_slab_size, although primitive=False has been applied.

Reproduce the bug

from pymatgen.core.surface import generate_all_slabs
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen.core import Structure

bulk_structure = Structure.from_file('TiNb_conv.cif')
sga = SpacegroupAnalyzer(bulk_structure)
conv_structure = sga.get_conventional_standard_structure()

slabs = generate_all_slabs(structure=conv_structure, 
                           max_index=2, 
                           min_slab_size=8,
                           min_vacuum_size=15, 
                           center_slab=True,
                           primitive=False)

The cif file for TiNb_conv.cif is:

# generated using pymatgen
data_TiNb
_symmetry_space_group_name_H-M   'P 1'
_cell_length_a   3.28200886
_cell_length_b   4.63609090
_cell_length_c   4.64594597
_cell_angle_alpha   90.00000000
_cell_angle_beta   90.00000000
_cell_angle_gamma   90.00000000
_symmetry_Int_Tables_number   1
_chemical_formula_structural   TiNb
_chemical_formula_sum   'Ti2 Nb2'
_cell_volume   70.69128019
_cell_formula_units_Z   2
loop_
 _symmetry_equiv_pos_site_id
 _symmetry_equiv_pos_as_xyz
  1  'x, y, z'
loop_
 _atom_site_type_symbol
 _atom_site_label
 _atom_site_symmetry_multiplicity
 _atom_site_fract_x
 _atom_site_fract_y
 _atom_site_fract_z
 _atom_site_occupancy
  Ti  Ti0  1  0.50000000  0.00000000  0.50000000  1.0
  Ti  Ti1  1  0.00000000  0.50000000  0.50000000  1.0
  Nb  Nb2  1  0.00000000  0.00000000  0.00000000  1.0
  Nb  Nb3  1  0.50000000  0.50000000  0.00000000  1.0

---

Try to pin point the bug

I have done some digging to understand the source of the bug.

The lattice constant of oriented_unit_cell for (010) surface is:

The lattice constants of oriented_unit_cell is:
-3.28,0.00,-0.00
0.00,0.00,-4.65
0.00,-4.64,-0.00

But the fractional coordinates of the oriented_unit_cell is (after dividing nlayers):

Coordination of the oriented_unit_cell is:
0.50, 0.50, 0.04
0.00, 0.50, 0.12
0.00, 0.00, 0.04
0.50, 0.00, 0.12

lattice for generated (010) slab
3.28, 0.00, 0.00
-0.00, 4.65, 0.00
0.00, 0.00, 27.82

Since the length of c-direction is 27.82Å. The length of c lattice in oriented_unit_cell calculated from fractional coordinate is: (0.12-0.04)*27.82 = 2.23Å, which is not consistent with the 4.64Å shown in the lattice matrix of oriented_unit_cell.

I think it could be related to the fact that in (010) surface, the oriented_unit_cell has FCC structure, it should have 6 atoms per unit cell, but 4 atoms are found. A schematic is shown in below, where yellow highlight is the oriented_unit_cell from lattice matrix and blue highlight is the oriented_unit_cell from fractional coordinates.

[hongyi-zhao]

I confirmed this problem as follows:

from numpy.linalg import norm
from numpy import cross
from pymatgen.core.structure import Structure
from pymatgen.core.surface import generate_all_slabs
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
import numpy as np

def calculate_slab_thickness(slab):
    # Get the first two lattice vectors of the slab, which are in the slab's plane
    a = slab.lattice.matrix[0]
    b = slab.lattice.matrix[1]

    # Calculate the normal vector as the unit vector of the cross product of vectors a and b
    normal_vector = cross(a, b)
    normal_vector = normal_vector / norm(normal_vector)

    # Calculate the thickness of the slab
    distances = [np.dot(site.coords, normal_vector) for site in slab.sites]
    thickness = max(distances) - min(distances)
    return thickness

# Load the structure from a CIF file
bulk_structure = Structure.from_file('TiNb_conv.cif')
sga = SpacegroupAnalyzer(bulk_structure)
conv_structure = sga.get_conventional_standard_structure()

# Generate slabs with specified parameters
slabs = generate_all_slabs(structure=conv_structure,
                           max_index=2,
                           min_slab_size=8,
                           min_vacuum_size=15,
                           center_slab=False, lll_reduce=True,
                           primitive=False)

# Filter slabs based on their calculated thickness
filtered_slabs = [slab for slab in slabs if calculate_slab_thickness(slab) < 8.0]

# Calculate and print the thickness of filtered slabs
for i, slab in enumerate(filtered_slabs):
    thickness = calculate_slab_thickness(slab)
    print(f"Filtered Slab {i+1} thickness: {thickness:.2f} Å")

The result is as follows:

Filtered Slab 1 thickness: 6.96 Å
Filtered Slab 2 thickness: 6.70 Å
Filtered Slab 3 thickness: 6.95 Å
Filtered Slab 4 thickness: 7.26 Å

The TiNb_conv.cif used above is also attached: TiNb_conv.zip

See here for the related discussion.