Source code for ase.build.general_surface

import numpy as np
from numpy.linalg import norm, solve

from ase.utils import gcd, basestring
from ase.build import bulk
from copy import deepcopy


[docs]def surface(lattice, indices, layers, vacuum=None, tol=1e-10, termination=0): """Create surface from a given lattice and Miller indices. lattice: Atoms object or str Bulk lattice structure of alloy or pure metal. Note that the unit-cell must be the conventional cell - not the primitive cell. One can also give the chemical symbol as a string, in which case the correct bulk lattice will be generated automatically. indices: sequence of three int Surface normal in Miller indices (h,k,l). layers: int Number of equivalent layers of the slab. vacuum: float Amount of vacuum added on both sides of the slab. termination: int The termination "number" for your crystal. The same value will not produce the same termination for different symetrically identical bulk structures, but changing this value allows your to explore all the possible terminations for the bulk structure you provide it. note: this code is not well tested """ indices = np.asarray(indices) if indices.shape != (3,) or not indices.any() or indices.dtype != int: raise ValueError('%s is an invalid surface type' % indices) if isinstance(lattice, basestring): lattice = bulk(lattice, cubic=True) h, k, l = indices h0, k0, l0 = (indices == 0) if termination != 0: #changing termination import warnings warnings.warn('Work on changing terminations is currently in ' 'progress. Code may not behave as expected.') lattice1 = deepcopy(lattice) cell = lattice1.get_cell() pt = [0,0,0] millers = list(indices) for index,item in enumerate(millers): if item == 0: millers[index] = 10**9 #make zeros large numbers elif pt == [0,0,0]: #for numerical stability pt = list(cell[index]/float(item)/np.linalg.norm(cell[index])) h1,k1,l1 = millers N = np.array(cell[0]/h1+cell[1]/k1+cell[2]/l1) n = N/np.linalg.norm(N)# making a unit vector normal to cut plane d = [np.round(np.dot(n,(a-pt)),4) for a in lattice.get_scaled_positions()] d = set(d) d = sorted(list(d)) d = [0]+d #distances of atoms from cut plane displacement = (h*cell[0]+k*cell[1]+l*cell[2])*d[termination] lattice1.positions += displacement lattice = lattice1 if h0 and k0 or h0 and l0 or k0 and l0: # if two indices are zero if not h0: c1, c2, c3 = [(0, 1, 0), (0, 0, 1), (1, 0, 0)] if not k0: c1, c2, c3 = [(0, 0, 1), (1, 0, 0), (0, 1, 0)] if not l0: c1, c2, c3 = [(1, 0, 0), (0, 1, 0), (0, 0, 1)] else: p, q = ext_gcd(k, l) a1, a2, a3 = lattice.cell # constants describing the dot product of basis c1 and c2: # dot(c1,c2) = k1+i*k2, i in Z k1 = np.dot(p * (k * a1 - h * a2) + q * (l * a1 - h * a3), l * a2 - k * a3) k2 = np.dot(l * (k * a1 - h * a2) - k * (l * a1 - h * a3), l * a2 - k * a3) if abs(k2) > tol: i = -int(round(k1 / k2)) # i corresponding to the optimal basis p, q = p + i * l, q - i * k a, b = ext_gcd(p * k + q * l, h) c1 = (p * k + q * l, -p * h, -q * h) c2 = np.array((0, l, -k)) // abs(gcd(l, k)) c3 = (b, a * p, a * q) surf = build(lattice, np.array([c1, c2, c3]), layers, tol) if vacuum is not None: surf.center(vacuum=vacuum, axis=2) return surf
def build(lattice, basis, layers, tol): surf = lattice.copy() scaled = solve(basis.T, surf.get_scaled_positions().T).T scaled -= np.floor(scaled + tol) surf.set_scaled_positions(scaled) surf.set_cell(np.dot(basis, surf.cell), scale_atoms=True) surf *= (1, 1, layers) a1, a2, a3 = surf.cell surf.set_cell([a1, a2, np.cross(a1, a2) * np.dot(a3, np.cross(a1, a2)) / norm(np.cross(a1, a2))**2]) # Change unit cell to have the x-axis parallel with a surface vector # and z perpendicular to the surface: a1, a2, a3 = surf.cell surf.set_cell([(norm(a1), 0, 0), (np.dot(a1, a2) / norm(a1), np.sqrt(norm(a2)**2 - (np.dot(a1, a2) / norm(a1))**2), 0), (0, 0, norm(a3))], scale_atoms=True) surf.pbc = (True, True, False) # Move atoms into the unit cell: scaled = surf.get_scaled_positions() scaled[:, :2] %= 1 surf.set_scaled_positions(scaled) surf.cell[2] = 0.0 return surf def ext_gcd(a, b): if b == 0: return 1, 0 elif a % b == 0: return 0, 1 else: x, y = ext_gcd(b, a % b) return y, x - y * (a // b)