VASP - suggestion on optimization of interfaced slabs

Question

I am currently trying to create a CeO2 (100) / CoO (200) interface for a bulk oxide. Unlike what I usually see, slabs stacked on the z-axis ("vertically"), I would like to stack the two "horizontally", in this case along the x-axis, so that both desired split surfaces are exposed to possible adsorption. I show here an image (z-view) of the system generated via the ASE. The slab was generated from .cif files of CoO and CeO2, cleaved to the desired indices and merged.

Although I have been working with VASP for a few years now, this is the first time I have approached systems of this type. In fact, I soon realised that it is not at all easy, as CoO has a smaller lattice constant than CeO2, so you can clearly see that the structure undergoes strain. I thought of tackling the problem by

• 1. create a zero vacuum cell (emulating bulk) with the slabs to optimise the positions and thus the cell parameters (perhaps a few cycles of calculations)
• 2. at the end, add the vacuum in the z-direction to obtain the real slabs.
• 3. Increase precision of the calculations (spin polarization, Hubbard, etc...)

My main concern is the validity of step 1. and how could i overcome it in case it is wrong. Because i am not sure surfaces are meant to be "bulk-like" stacked from the VASP PBC point of view.

At this point, is this the right way to proceed or is it a total dead end? Could anyone suggest a different workflow based on their previous experience, perhaps one that overcomes the problems of ion optimisation (always critical...) ? Unfortunately I have found relatively little information about this in the literature.

Thanks in advance.

Answer 1

Thanks to @ShaunHa's reply, which cleverly allows slabs to be stacked in the z-direction, I created a customised (simplified) add_adsorbate function that can stack slabs in the 'x', 'y' or 'z' direction. Note that all bond lengths and distances for each slab are those for the CIF file used, so no deformation is inserted.

Here are the code-snippet and the generated structures.

from ase import Atoms, Atom
import numpy as np

def my_add_adsorbate(slab, adsorbate, height, ax='z'):
    axes = {
        'x' : 0,
        'y' : 1,
        'z' : 2
    }
    indexes = {
        'x' : (1,3,1), # (y,z, step)
        'y' : (0,3,2), # (x,z, step)
        'z' : (0,2,1)  # (x,z, step)
    }

    i,j,s = indexes[ax]
    idx = axes[ax]

    pos = np.array([0.0, 0.0])  
    spos = np.array([0.0, 0.0])  

    cell = slab.get_cell()[i:j:s, i:j:s]

    pos += np.dot(spos, cell)
    ads = adsorbate

    # Get desired ax
    a = slab.positions[:, idx]

    # Max value position...
    a = a.argmax()

    # ...and value+height
    v = slab.positions[a, idx] + height

    # Move adsorbate into position
    vector = [0,0,0]
    vector[idx] = v
    print(vector)
    ads.translate(vector)

    # Attach the adsorbate
    slab.extend(ads)


from ase.visualize import view

# Build CeO2
CeO2 = read('CeO2_mp-1018664_primitive.cif')
# Build CoO 
CoO = read('CoO_fm3m.cif')
# Match the two unit cells in the z direction
ax = 'y'
axes = {'x' : 0,
        'y' : 1,
        'z' : 2}

ax_ = axes[ax]

CeO2.cell[ax_] = CoO.cell[ax_]
indexes = [1,0,0]

CeO2 = surface(
    CeO2, 
    indices=indexes, 
    layers=3, #3 
    periodic=True
    )

CoO = surface(
    CoO, 
    indices=indexes, 
    layers=3, #3 
    periodic=True
    )


fracs = CeO2.get_scaled_positions()

CeO2.cell = CoO.cell

CeO2.positions = fracs @ CoO.cell

# Stack the two structures
my_add_adsorbate(CeO2, adsorbate=CoO, height=3., ax=ax)
CeO2.center(vacuum=0, axis=ax_)

view(CeO2)