mwe_surface_sites.py

#!/usr/bin/env python

import pprint
try:
    import io as StringIO
except:
    import StringIO


import ase.io
import catkit.surface

bulk_cif = (
    'data_image0\n'
    '_cell_length_a       3.91\n'
    '_cell_length_b       3.91\n'
    '_cell_length_c       3.91\n'
    '_cell_angle_alpha    90\n'
    '_cell_angle_beta     90\n'
    '_cell_angle_gamma    90\n'
    '\n'
    '_symmetry_space_group_name_H-M    "P 1"\n'
    '_symmetry_int_tables_number       1\n'
    '\n'
    'loop_\n'
    '  _symmetry_equiv_pos_as_xyz\n'
    "  'x, y, z'\n"
    '\n'
    'loop_\n'
    '  _atom_site_label\n'
    '  _atom_site_occupancy\n'
    '  _atom_site_fract_x\n'
    '  _atom_site_fract_y\n'
    '  _atom_site_fract_z\n'
    '  _atom_site_thermal_displace_type\n'
    '  _atom_site_B_iso_or_equiv\n'
    '  _atom_site_type_symbol\n'
    '  Pd1      1.0000 0.00000  0.00000  0.00000  Biso   1.000  Pd\n'
    '  Cu1      1.0000 0.00000  0.50000  0.50000  Biso   1.000  Cu\n'
    '  Cu2      1.0000 0.50000  0.00000  0.50000  Biso   1.000  Cu\n'
    '  Cu3      1.0000 0.50000  0.50000  0.00000  Biso   1.000  Cu\n')

slab_cif = (
    ('data_image0\n'
     '_cell_length_a       8.74303\n'
     '_cell_length_b       5.52958\n'
     '_cell_length_c       22.3944\n'
     '_cell_angle_alpha    90\n'
     '_cell_angle_beta     90\n'
     '_cell_angle_gamma    50.7685\n'
     '\n'
     '_symmetry_space_group_name_H-M    "P 1"\n'
     '_symmetry_int_tables_number       1\n'
     '\n'
     'loop_\n'
     '  _symmetry_equiv_pos_as_xyz\n'
     "  'x, y, z'\n"
     '\n'
     'loop_\n'
     '  _atom_site_label\n'
     '  _atom_site_occupancy\n'
     '  _atom_site_fract_x\n'
     '  _atom_site_fract_y\n'
     '  _atom_site_fract_z\n'
     '  _atom_site_thermal_displace_type\n'
     '  _atom_site_B_iso_or_equiv\n'
     '  _atom_site_type_symbol\n'
     '  Pd1      1.0000 0.66667  0.83333  0.48218  Biso   1.000  Pd\n'
     '  Cu1      1.0000 0.66667  0.33333  0.48218  Biso   1.000  Cu\n'
     '  Cu2      1.0000 0.33333  0.91667  0.44654  Biso   1.000  Cu\n'
     '  Cu3      1.0000 0.33333  0.41667  0.44654  Biso   1.000  Cu\n'
     '  Pd2      1.0000 0.33333  0.66667  0.55346  Biso   1.000  Pd\n'
     '  Cu4      1.0000 0.33333  0.16667  0.55346  Biso   1.000  Cu\n'
     '  Cu5      1.0000 0.00000  0.75000  0.51782  Biso   1.000  Cu\n'
     '  Cu6      1.0000 0.00000  0.25000  0.51782  Biso   1.000  Cu\n')
)


slab_cifs = [slab_cif] * 2  # just for demo purpose


if __name__ == '__main__':
    mem_file = StringIO.StringIO()
    mem_file.write(bulk_cif)
    mem_file.seek(0)
    bulk_atoms = ase.io.read(mem_file, format='cif')

    gen = catkit.surface.SlabGenerator(
        bulk=bulk_atoms,
        miller_index=[1, 1, 1],
        layers=8,
    )

    in_mem_files = []
    images = []
    for cif_image in slab_cifs:
        mem_file = StringIO.StringIO()
        mem_file.write(cif_image)
        mem_file.seek(0)
        atoms = ase.io.read(mem_file, format='cif')
        images.append(atoms)

    sites_list = []
    for atoms in images:
        # As a WORKAROUND one can create a new slab
        # and comment out the following lines
        # generator but this seems counter-intuitive
        #gen = catkit.surface.SlabGenerator(
            #bulk=bulk_atoms,
            #miller_index=[1, 1, 1
                         #],
            #layers=4,)

        sites = gen.get_adsorption_sites(atoms)
        sites_list.append(sites)

    pprint.pprint(sites_list)

The surface sites you are loading from the CIF file are of a (2, 1, 1) surface with 4 layers. That's too thin to properly identify the surface sites using current methods in CatKit. I'll work on making this more robust.

There are some safeguards in place now to prevent identification of surface sites when slabs are too thin.

I also added functionality to produce slabs of a minimum desired width. Might help make usage more robust.