Skip to content

Uploaded WIP terminations.py #9

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
kevinlar2 wants to merge 10 commits into
base: main
Choose a base branch
from
Draft

Uploaded WIP terminations.py #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
ed3cf9d
Uploaded WIP terminations.py
kevinlar2 Jun 3, 2024
cdd375d
Updated WIP terminations.py
kevinlar2 Jul 3, 2024
7f40bc1
Update and rename find_terminations_3_1_1.py to OgreInterface/surface…
kevinlar2 Jul 3, 2024
35ce7a0
Delete OgreInterface/surfaces/find_terminations_3_0_0.py
kevinlar2 Jul 5, 2024
c52d199
Refactored and protected mutables
kevinlar2 Jul 5, 2024
d3d73e9
Integrate terminations into base_surface_generator.py
kevinlar2 Aug 21, 2024
a3a360a
Added terminations.py
kevinlar2 Aug 21, 2024
c921a7e
Update terminations.py
kevinlar2 Aug 29, 2024
9784185
Update terminations.py
kevinlar2 Aug 29, 2024
536bd9c
Update terminations.py
kevinlar2 Aug 30, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
326 changes: 17 additions & 309 deletions OgreInterface/generate/base_surface_generator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
from OgreInterface.surfaces.oriented_bulk import OrientedBulk
from OgreInterface.surfaces.surface import Surface
from OgreInterface.surfaces.molecular_surface import MolecularSurface
from OgreInterface.surfaces.terminations import Terminator

SelfBaseSurfaceGenerator = TypeVar(
"SelfBaseSurfaceGenerator", bound="BaseSurfaceGenerator"
Expand Down Expand Up @@ -105,6 +106,7 @@ def __init__(
refine_structure: bool = True,
make_planar: bool = True,
generate_all: bool = True,
smoothest_only: bool = False,
lazy: bool = False,
suppress_warnings: bool = False,
layer_grouping_tolarence: Optional[float] = None,
Expand All @@ -126,6 +128,7 @@ def __init__(

self.vacuum = vacuum
self.generate_all = generate_all
self.smoothest_only = smoothest_only
self.lazy = lazy

self.obs = OrientedBulk(
Expand All @@ -143,6 +146,15 @@ def __init__(
(minimum_thickness // self.obs.layer_thickness) + 1
)

self.terminator = Terminator(
bulk=self.bulk_structure,
plane=self.miller_index,
molecular=True if self._surface_type == MolecularSurface else False,
generate_all=self.generate_all,
num_layers=self.layers,
vacuum=self.vacuum
)

if not self.lazy:
self._slabs = self._generate_slabs()
else:
Expand All @@ -159,6 +171,7 @@ def from_file(
refine_structure: bool = True,
make_planar: bool = True,
generate_all: bool = True,
smoothest_only: bool = False,
lazy: bool = False,
suppress_warnings: bool = False,
layer_grouping_tolarence: Optional[float] = None,
Expand Down Expand Up @@ -193,6 +206,7 @@ def from_file(
refine_structure=refine_structure,
make_planar=make_planar,
generate_all=generate_all,
smoothest_only=smoothest_only,
lazy=lazy,
suppress_warnings=suppress_warnings,
layer_grouping_tolarence=layer_grouping_tolarence,
Expand Down Expand Up @@ -240,158 +254,6 @@ def _get_point_group_operations(self):

return unique_operations

def _get_slab(
self,
slab_base: OrientedBulk,
shift: float = 0.0,
) -> Tuple[Structure, Structure, float, Tuple[int, ...]]:
"""
This method takes in shift value for the c lattice direction and
generates a slab based on the given shift. You should rarely use this
method. Instead, it is used by other generation algorithms to obtain
all slabs.

Args:
slab_base: Oriented bulk structure used to generate the slab
shift: A shift value in fractional c-coordinates that determines
how much the slab_base should be shifted to select a given
termination.
tol: Optional tolarance for grouping the atomic layers together

Returns:
Returns a tuple of the shifted slab base, orthogonalized slab,
non-orthogonalized slab, actual value of the vacuum in angstroms,
and the tuple of layer indices and bulk equivalents that is used
to filter out duplicate surfaces.
"""
# Shift the slab base to the termination defined by the shift input
slab_base.translate_sites(
vector=[0, 0, -shift],
frac_coords=True,
)

# Round and mod the structure
slab_base.round(tol=6)

# Get the fractional c-coords
c_coords = slab_base._oriented_bulk_structure.frac_coords[:, -1]

# Calculate the shifts again on the shifted structure to get the upper
# and lower bounds of where an atomic layer should be.
shifts = self._calculate_possible_shifts(
structure=slab_base._oriented_bulk_structure,
)
shifts += [1.0]

# Group the upper and lower bounds into a list of tuples
atomic_layer_bounds = [
(shifts[i], shifts[i + 1]) for i in range(len(shifts) - 1)
]

# Define an array of -1's that will get filled in later with atomic
# layer indices
atomic_layers = -np.ones(len(c_coords))
for i, (bottom_bound, top_bound) in enumerate(atomic_layer_bounds):
# Find atoms that have c-position between the top and bottom bounds
layer_mask = (c_coords > bottom_bound) & (c_coords < top_bound)

# Set the atomic layer index to i
atomic_layers[layer_mask] = i

# Add the atomic layer site property to the slab base
slab_base.add_site_property(
"atomic_layer_index",
np.round(atomic_layers).astype(int).tolist(),
)

# Get the bulk equivalent to create the key associated with the given
# surface so that we can extract unique surface terminations later on
bulk_equiv = np.array(slab_base.site_properties["bulk_equivalent"])

# The surface key is sorted by atomic layer and the bulk equivalent
# i.e. [(0, 0), (0, 1), (0, 2), (1, 0), (1, 2), ...]
surf_key = sorted(
[(idx, eq) for idx, eq in zip(atomic_layers, bulk_equiv)],
key=lambda x: (x[0], x[1]),
)

# Concatenate the key and turn it into one long tuple of ints
surf_key = tuple(np.concatenate(surf_key).astype(int))

# Get the top c-coord
top_c = c_coords.max()

# Get the bottom c-coord
# bot_c = c_coords.min()

# Get the inds of the top atoms so we can shift of so the top atom
# has a and b positions of zero.
max_c_inds = np.where(np.isclose(top_c, c_coords))[0]

dists = []
for i in max_c_inds:
# Get distance from the origin
dist, image = slab_base[i].distance_and_image_from_frac_coords(
fcoords=[0.0, 0.0, 0.0]
)
dists.append(dist)

# Get the atom index of the top atom that is closest to a=0, b=0
horiz_shift_ind = max_c_inds[np.argmin(dists)]

# Find the planar shift required to plane the top atom at a=0, b=0
horiz_shift = -slab_base[horiz_shift_ind].frac_coords
horiz_shift[-1] = 0

# Shift the slab base (this is mostly just for aesthetics)
slab_base.translate_sites(
vector=horiz_shift,
frac_coords=True,
)

# Round and mod the structure
slab_base.round(tol=6)

# Calculate number of empty unit cells are needed for the vacuum
# Make sure the number is even so the surface can be nicely centered
# in the vacuum region.
vacuum_scale = self.vacuum // self.obs.layer_thickness

if vacuum_scale % 2:
vacuum_scale += 1

if vacuum_scale == 0:
vacuum_scale = 2

# Get the actuall vacuum in angstroms
vacuum = self.obs.layer_thickness * vacuum_scale

# Create the non-orthogonalized surface
non_orthogonal_slab = utils.get_layer_supercell(
structure=slab_base._oriented_bulk_structure,
layers=self.layers,
vacuum_scale=vacuum_scale,
)
utils.sort_slab(non_orthogonal_slab)
# non_orthogonal_slab.sort()

# Center the surfaces within the vacuum region by shifting along c
center_shift = 0.5 * (vacuum_scale / (vacuum_scale + self.layers))

non_orthogonal_slab.translate_sites(
indices=range(len(non_orthogonal_slab)),
vector=[0, 0, center_shift],
frac_coords=True,
to_unit_cell=True,
)

return (
slab_base,
non_orthogonal_slab,
vacuum,
surf_key,
)

def _generate_slabs(self) -> List[Union[Surface, MolecularSurface]]:
"""
This function is used to generate slab structures with all unique
Expand All @@ -400,162 +262,8 @@ def _generate_slabs(self) -> List[Union[Surface, MolecularSurface]]:
Returns:
A list of Surface classes
"""
# Determine if all possible terminations are generated
slab_base = self._get_slab_base()
possible_shifts = self._calculate_possible_shifts(
structure=slab_base._oriented_bulk_structure
)
shifted_slab_bases = []
non_orthogonal_slabs = []
surface_keys = []

if not self.generate_all:
(
shifted_slab_base,
non_orthogonal_slab,
actual_vacuum,
surf_key,
) = self._get_slab(
slab_base=deepcopy(slab_base),
shift=possible_shifts[0],
)
non_orthogonal_slab.sort_index = 0
shifted_slab_bases.append(shifted_slab_base)
non_orthogonal_slabs.append(non_orthogonal_slab)
surface_keys.append((surf_key, 0))
else:
for i, possible_shift in enumerate(possible_shifts):
(
shifted_slab_base,
non_orthogonal_slab,
actual_vacuum,
surf_key,
) = self._get_slab(
slab_base=deepcopy(slab_base),
shift=possible_shift,
)
non_orthogonal_slab.sort_index = i
shifted_slab_bases.append(shifted_slab_base)
non_orthogonal_slabs.append(non_orthogonal_slab)
surface_keys.append((surf_key, i))

surfaces = []

sorted_surface_keys = sorted(surface_keys, key=lambda x: x[0])

groups = groupby(sorted_surface_keys, key=lambda x: x[0])

unique_inds = []
for group_key, group in groups:
_, inds = list(zip(*group))
unique_inds.append(min(inds))

unique_inds.sort()

# Loop through slabs to ensure that they are all properly oriented and reduced
# Return Surface objects
for i in unique_inds:
# Create the Surface object
surface = self._surface_type(
slab=non_orthogonal_slabs[i], # KEEP
oriented_bulk=shifted_slab_bases[i], # KEEP
miller_index=self.miller_index, # KEEP
layers=self.layers, # KEEP
vacuum=actual_vacuum, # KEEP
termination_index=i, # KEEP
)
surfaces.append(surface)

return surfaces

def _calculate_possible_shifts(
self,
structure: Structure,
):
"""
This function calculates the possible shifts that need to be applied to
the oriented bulk structure to generate different surface terminations

Args:
structure: Oriented bulk structure
tol: Grouping tolarence in angstroms.
If None, it will automatically be calculated based on the input
structure.

Returns:
A list of fractional shift values along the c-vector
"""
frac_coords = structure.frac_coords[:, -1]

# Projection of c lattice vector in
# direction of surface normal.
h = self.obs.layer_thickness

if self._layer_grouping_tolarence is None:
cart_coords = structure.cart_coords
projected_coords = np.dot(cart_coords, self.obs.surface_normal)
extended_projected_coords = np.round(
np.concatenate(
[
projected_coords - h,
projected_coords,
projected_coords + h,
]
),
5,
)
unique_cart_coords = np.sort(np.unique(extended_projected_coords))
diffs = np.diff(unique_cart_coords)
max_diff = diffs.max()
self._layer_grouping_tolarence = 0.15 * max_diff

n = len(frac_coords)

if n == 1:
# Clustering does not work when there is only one data point.
shift = frac_coords[0] + 0.5
return [shift - math.floor(shift)]

# We cluster the sites according to the c coordinates. But we need to
# take into account PBC. Let's compute a fractional c-coordinate
# distance matrix that accounts for PBC.
dist_matrix = np.zeros((n, n))

for i, j in combinations(list(range(n)), 2):
if i != j:
cdist = frac_coords[i] - frac_coords[j]
cdist = abs(cdist - round(cdist)) * h
dist_matrix[i, j] = cdist
dist_matrix[j, i] = cdist

condensed_m = squareform(dist_matrix)
z = linkage(condensed_m)
clusters = fcluster(
z,
self._layer_grouping_tolarence,
criterion="distance",
)
if self.smoothest_only:
return [self.terminator.smoothest_surface]

# Generate dict of cluster to c val - doesn't matter what the c is.
c_loc = {c: frac_coords[i] for i, c in enumerate(clusters)}

# Put all c into the unit cell.
possible_c = [c - math.floor(c) for c in sorted(c_loc.values())]

# Calculate the shifts
nshifts = len(possible_c)
shifts = []
for i in range(nshifts):
if i == nshifts - 1:
# There is an additional shift between the first and last c
# coordinate. But this needs special handling because of PBC.
shift = (possible_c[0] + 1 + possible_c[i]) * 0.5
if shift > 1:
shift -= 1
else:
shift = (possible_c[i] + possible_c[i + 1]) * 0.5
shifts.append(shift - math.floor(shift))

shifts = sorted(shifts)

return shifts
return self._terminator.surfaces
Loading