Commit d26ad8ff authored by rym's avatar rym
Browse files

Compute relative angles by polar decomposition

parent 8b49594e

tests/.ipynb_checkpoints/vdW_structures_test-checkpoint.ipynb

0 → 100644
+179 −0
Original line number Diff line number Diff line 
%% Cell type:code id:3adc4ce9-130e-4a15-bf22-df5bccd44c94 tags:



``` python

from pymatgen.core.structure import Structure

from ase.geometry.geometry import get_layers

from vdW_structures.vdW_structure import VdWStructure

import numpy as np

from typing import List, Tuple

```



%% Cell type:code id:5889879d-17e3-4a1e-aa34-2bdd65684786 tags:



``` python

s = Structure.from_file('original_structures/BiSe2_slab.vasp')

#s = Structure.from_file('original_structures/Bi2Se3.vasp')

vdW_s = VdWStructure(s, minimum_vdW_gap=2.4)

vdW_s = vdW_s.z_solve_shift_vdW_layers(layer_inds=[0], solve_shift=[4], adjust_z_cell=True)

vdW_s.vdW_spacings

```



%% Output



    [4.053074184176202, 2.935904452246059]



%% Cell type:code id:0c76c871-6c23-496f-8c5b-c71be9c0b585 tags:



``` python

def get_vdW_layers(structure: Structure, minimum_vdW_gap=2.4):

    """Identify van der Waals layers in the structure based on spacing along the (0, 0, 1) Miller index."""



    vdW_layers, z_images, vdW_spacings = [[]], [[]], []

    atoms = structure.to_ase_atoms()

    layers, layer_dist = get_layers(atoms, (0, 0, 1))



    #print(layers, layer_dist)



    # Convert layers into the [[atomic layers], spacings]

    atomic_layers, atomic_spacings = [[] for i in range(max(layers)+1)], [None for i in range(max(layers)+1)]

    for site_ind, layer_ind in enumerate(layers):

        atomic_layers[layer_ind].append(site_ind)

        atomic_spacings[layer_ind] = layer_dist[layer_ind]



    #print(atomic_layers, atomic_spacings)



    # Get z-coordinates sorted by layer distance

    site_values = [layer_dist[layers[i]] for i in range(len(atoms))]

    sorted_indices = sorted(range(len(atoms)), key=lambda i: site_values[i])

    sorted_zs = [site_values[i] for i in sorted_indices]



    # Compute layer spacings, accounting for periodic boundaries

    diffs = [structure.lattice.c - sorted_zs[-1] + sorted_zs[0]] + [

        sorted_zs[i + 1] - sorted_zs[i] for i in range(len(sorted_zs) - 1)

    ]



    # Group atoms into vdW layers

    layer_index = 0

    for i, diff in enumerate(diffs):

        if diff > minimum_vdW_gap:

            vdW_spacings.append(diff)

            vdW_layers.append([])

            z_images.append([])

            layer_index += 1

        vdW_layers[layer_index].append(sorted_indices[i])

        z_images[layer_index].append(0)



    # Handle periodic boundary conditions

    if len(vdW_layers) > 1 and diffs[0] < minimum_vdW_gap:

        z_images[0] = [-1] * len(vdW_layers[0]) + [0] * len(vdW_layers[-1])

        vdW_layers[0] += vdW_layers.pop()

        z_images.pop()



    return [layer for layer in vdW_layers if layer], [z_images[i] for i, layer in enumerate(vdW_layers) if layer], vdW_spacings, atomic_layers, atomic_spacings





def shift_to_vdW_gap(

    structure: Structure,

    layers: List[List[int]],

    images: List[List[int]],

    spacings: List[float],

    atomic_layers: List[List[float]],

    atomic_spacings: List[float]

):



    def find_sublist_index(lists: list[list[int]], target: int) -> int:

        for i, sublist in enumerate(lists):

            if target in sublist:

                return i  # Return the index as soon as we find the target

        return -1



    # Identify the largest vdW gap and its index

    max_gap_idx = spacings.index(max(spacings))



    c_vector = structure.lattice.matrix[2]



    # Identify the largest site index from the below gap and

    if not max_gap_idx == 0: # Ignore if the largest gap is already across the PBC

        top_site_ind = layers[max_gap_idx][0]

        bottom_site_ind = layers[max_gap_idx - 1][-1]

        top_distance = atomic_spacings[find_sublist_index(atomic_layers, top_site_ind)]

        bottom_distance = atomic_spacings[find_sublist_index(atomic_layers, bottom_site_ind)]

        shift_distance = -1 * np.average([top_distance, bottom_distance])

        shift_vector = shift_distance * (c_vector / np.linalg.norm(c_vector))

        print(shift_vector)



        for site in structure:

            site.coords += shift_vector



        return Structure(

            lattice=structure.lattice,

            coords=[site.coords for site in structure],

            species=[site.specie for site in structure],

            coords_are_cartesian=True,

            to_unit_cell=True

            ).get_sorted_structure()



    else:

        return structure

```



%% Cell type:code id:b5ecc597-8f6d-4a2a-bd68-37032feb4011 tags:



``` python

layers, images, spacings, atomic_layers, atomic_spacings = get_vdW_layers(vdW_s.structure)

print(layers, images, spacings, atomic_layers, atomic_spacings)

'''

shifted_s = shift_to_vdW_gap(

    structure= vdW_s.structure,

    layers= layers,

    images= images,

    spacings= spacings,

    atomic_layers=atomic_layers,

    atomic_spacings=atomic_spacings

)

s_layers, s_images, s_spacings, s_atomic_layers, s_atomic_spacings = get_vdW_layers(shifted_s)

print(s_layers, s_images, s_spacings)

'''

```



%% Output



    [[8, 1, 4, 5, 0, 9], [10, 3, 6, 7, 2, 11]] [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]] [4.053074184176202, 2.935904452246059] [[8], [1], [4], [5], [0], [9], [10], [3], [6], [7], [2], [11]] [2.290107083779596, 2.994003526732217, 3.2876766256722334, 4.437926432997699, 4.731599531937719, 5.435495974890334, 8.371400427136393, 8.960167885072961, 9.429028389346314, 10.371894248570655, 10.840754752844, 11.429522210780569]



    '\nshifted_s = shift_to_vdW_gap(\n    structure= vdW_s.structure, \n    layers= layers, \n    images= images, \n    spacings= spacings,\n    atomic_layers=atomic_layers,\n    atomic_spacings=atomic_spacings\n)\ns_layers, s_images, s_spacings, s_atomic_layers, s_atomic_spacings = get_vdW_layers(shifted_s)\nprint(s_layers, s_images, s_spacings)\n'



%% Cell type:code id:63997040-f4ba-47e5-9197-16295b1920d3 tags:



``` python

shifted_s = shift_to_vdW_gap(

    structure= vdW_s.structure,

    layers= layers,

    images= images,

    spacings= spacings,

    atomic_layers=atomic_layers,

    atomic_spacings=atomic_spacings

)

s_layers, s_images, s_spacings, s_atomic_layers, s_atomic_spacings = get_vdW_layers(shifted_s)

print(s_layers, s_images, s_spacings, s_atomic_layers, s_atomic_spacings)

```



%% Output



    [[8, 1, 4, 5, 0, 9], [10, 3, 6, 7, 2, 11]] [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]] [4.053074184176202, 2.935904452246059] [[8], [1], [4], [5], [0], [9], [10], [3], [6], [7], [2], [11]] [2.290107083779596, 2.994003526732217, 3.2876766256722334, 4.437926432997699, 4.731599531937719, 5.435495974890334, 8.371400427136393, 8.960167885072961, 9.429028389346314, 10.371894248570655, 10.840754752844, 11.429522210780569]



%% Cell type:code id:017e136b-3f35-44b3-aa60-80f0a9577ddd tags:



``` python

vdW_s.structure.lattice.matrix

```



%% Output



    array([[ 2.09568326, -3.62982989,  0.        ],

           [ 2.09568326,  3.62982989,  0.        ],

           [ 0.        ,  0.        , 30.99278275]])



%% Cell type:code id:8c292c02-1da5-460b-bdf6-944dd1993af0 tags:



``` python

```

tests/original_structures/.ipynb_checkpoints/untitled-checkpoint.txt

0 → 100644
+20 −0
Original line number Diff line number Diff line 
Bi4 Se8

1.0

   4.1924904999999999    0.0000000000000000    0.0000000000000003

  -2.0962452500000013    8.3598153695497253   -0.7937824245145666

  -0.0000000000000000    0.0596282645423928   12.0751726993480037

Bi Se

4 8

direct

   0.4650013997226389    0.7282428648191984    0.3450828665136759 Bi

   0.7692165796508181    0.3366732246755562    0.2006051460916581 Bi

   0.0327581331200342    0.8637563316139885    0.8530472134485209 Bi

   0.7102548149997345    0.2187496953733891    0.6966800815621834 Bi

   0.0522071049778474    0.9026542753296150    0.2250234917583226 Se

   0.1820108743956097    0.1622618141651394    0.3206645208470112 Se

   0.4333817259382872    0.6650035172504946    0.7356649184577538 Se

   0.3096312221814812    0.4175025097368825    0.8140623765529508 Se

   0.3807804403509084    0.5598009460757367    0.1420775290930601 Se

   0.8534375390225485    0.5051151434190173    0.4036104835122735 Se

   0.0868523958005768    0.9719448569750740    0.6477252152258732 Se

   0.6561605523191913    0.1105611700123025    0.9020020797848313 Se
 
 No newline at end of file

tests/original_structures/BiSe2_slab.vasp

0 → 100644
+20 −0
Original line number Diff line number Diff line 
Bi4 Se8

1.0

   4.1924904999999999    0.0000000000000000    0.0000000000000003

  -2.0962452500000013    8.3598153695497253   -0.7937824245145666

  -0.0000000000000000    0.0596282645423928   12.0751726993480037

Bi Se

4 8

direct

   0.4650013997226389    0.7282428648191984    0.3450828665136759 Bi

   0.7692165796508181    0.3366732246755562    0.2006051460916581 Bi

   0.0327581331200342    0.8637563316139885    0.8530472134485209 Bi

   0.7102548149997345    0.2187496953733891    0.6966800815621834 Bi

   0.0522071049778474    0.9026542753296150    0.2250234917583226 Se

   0.1820108743956097    0.1622618141651394    0.3206645208470112 Se

   0.4333817259382872    0.6650035172504946    0.7356649184577538 Se

   0.3096312221814812    0.4175025097368825    0.8140623765529508 Se

   0.3807804403509084    0.5598009460757367    0.1420775290930601 Se

   0.8534375390225485    0.5051151434190173    0.4036104835122735 Se

   0.0868523958005768    0.9719448569750740    0.6477252152258732 Se

   0.6561605523191913    0.1105611700123025    0.9020020797848313 Se
 
 No newline at end of file

tests/vdW_structures_test.ipynb

0 → 100644
+179 −0
Original line number Diff line number Diff line 
%% Cell type:code id:3adc4ce9-130e-4a15-bf22-df5bccd44c94 tags:



``` python

from pymatgen.core.structure import Structure

from ase.geometry.geometry import get_layers

from vdW_structures.vdW_structure import VdWStructure

import numpy as np

from typing import List, Tuple

```



%% Cell type:code id:5889879d-17e3-4a1e-aa34-2bdd65684786 tags:



``` python

s = Structure.from_file('original_structures/BiSe2_slab.vasp')

#s = Structure.from_file('original_structures/Bi2Se3.vasp')

vdW_s = VdWStructure(s, minimum_vdW_gap=2.4)

vdW_s = vdW_s.z_solve_shift_vdW_layers(layer_inds=[0], solve_shift=[4], adjust_z_cell=True)

vdW_s.vdW_spacings

```



%% Output



    [4.053074184176202, 2.935904452246059]



%% Cell type:code id:0c76c871-6c23-496f-8c5b-c71be9c0b585 tags:



``` python

def get_vdW_layers(structure: Structure, minimum_vdW_gap=2.4):

    """Identify van der Waals layers in the structure based on spacing along the (0, 0, 1) Miller index."""



    vdW_layers, z_images, vdW_spacings = [[]], [[]], []

    atoms = structure.to_ase_atoms()

    layers, layer_dist = get_layers(atoms, (0, 0, 1))



    #print(layers, layer_dist)



    # Convert layers into the [[atomic layers], spacings]

    atomic_layers, atomic_spacings = [[] for i in range(max(layers)+1)], [None for i in range(max(layers)+1)]

    for site_ind, layer_ind in enumerate(layers):

        atomic_layers[layer_ind].append(site_ind)

        atomic_spacings[layer_ind] = layer_dist[layer_ind]



    #print(atomic_layers, atomic_spacings)



    # Get z-coordinates sorted by layer distance

    site_values = [layer_dist[layers[i]] for i in range(len(atoms))]

    sorted_indices = sorted(range(len(atoms)), key=lambda i: site_values[i])

    sorted_zs = [site_values[i] for i in sorted_indices]



    # Compute layer spacings, accounting for periodic boundaries

    diffs = [structure.lattice.c - sorted_zs[-1] + sorted_zs[0]] + [

        sorted_zs[i + 1] - sorted_zs[i] for i in range(len(sorted_zs) - 1)

    ]



    # Group atoms into vdW layers

    layer_index = 0

    for i, diff in enumerate(diffs):

        if diff > minimum_vdW_gap:

            vdW_spacings.append(diff)

            vdW_layers.append([])

            z_images.append([])

            layer_index += 1

        vdW_layers[layer_index].append(sorted_indices[i])

        z_images[layer_index].append(0)



    # Handle periodic boundary conditions

    if len(vdW_layers) > 1 and diffs[0] < minimum_vdW_gap:

        z_images[0] = [-1] * len(vdW_layers[0]) + [0] * len(vdW_layers[-1])

        vdW_layers[0] += vdW_layers.pop()

        z_images.pop()



    return [layer for layer in vdW_layers if layer], [z_images[i] for i, layer in enumerate(vdW_layers) if layer], vdW_spacings, atomic_layers, atomic_spacings





def shift_to_vdW_gap(

    structure: Structure,

    layers: List[List[int]],

    images: List[List[int]],

    spacings: List[float],

    atomic_layers: List[List[float]],

    atomic_spacings: List[float]

):



    def find_sublist_index(lists: list[list[int]], target: int) -> int:

        for i, sublist in enumerate(lists):

            if target in sublist:

                return i  # Return the index as soon as we find the target

        return -1



    # Identify the largest vdW gap and its index

    max_gap_idx = spacings.index(max(spacings))



    c_vector = structure.lattice.matrix[2]



    # Identify the largest site index from the below gap and

    if not max_gap_idx == 0: # Ignore if the largest gap is already across the PBC

        top_site_ind = layers[max_gap_idx][0]

        bottom_site_ind = layers[max_gap_idx - 1][-1]

        top_distance = atomic_spacings[find_sublist_index(atomic_layers, top_site_ind)]

        bottom_distance = atomic_spacings[find_sublist_index(atomic_layers, bottom_site_ind)]

        shift_distance = -1 * np.average([top_distance, bottom_distance])

        shift_vector = shift_distance * (c_vector / np.linalg.norm(c_vector))

        print(shift_vector)



        for site in structure:

            site.coords += shift_vector



        return Structure(

            lattice=structure.lattice,

            coords=[site.coords for site in structure],

            species=[site.specie for site in structure],

            coords_are_cartesian=True,

            to_unit_cell=True

            ).get_sorted_structure()



    else:

        return structure

```



%% Cell type:code id:b5ecc597-8f6d-4a2a-bd68-37032feb4011 tags:



``` python

layers, images, spacings, atomic_layers, atomic_spacings = get_vdW_layers(vdW_s.structure)

print(layers, images, spacings, atomic_layers, atomic_spacings)

'''

shifted_s = shift_to_vdW_gap(

    structure= vdW_s.structure,

    layers= layers,

    images= images,

    spacings= spacings,

    atomic_layers=atomic_layers,

    atomic_spacings=atomic_spacings

)

s_layers, s_images, s_spacings, s_atomic_layers, s_atomic_spacings = get_vdW_layers(shifted_s)

print(s_layers, s_images, s_spacings)

'''

```



%% Output



    [[8, 1, 4, 5, 0, 9], [10, 3, 6, 7, 2, 11]] [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]] [4.053074184176202, 2.935904452246059] [[8], [1], [4], [5], [0], [9], [10], [3], [6], [7], [2], [11]] [2.290107083779596, 2.994003526732217, 3.2876766256722334, 4.437926432997699, 4.731599531937719, 5.435495974890334, 8.371400427136393, 8.960167885072961, 9.429028389346314, 10.371894248570655, 10.840754752844, 11.429522210780569]



    '\nshifted_s = shift_to_vdW_gap(\n    structure= vdW_s.structure, \n    layers= layers, \n    images= images, \n    spacings= spacings,\n    atomic_layers=atomic_layers,\n    atomic_spacings=atomic_spacings\n)\ns_layers, s_images, s_spacings, s_atomic_layers, s_atomic_spacings = get_vdW_layers(shifted_s)\nprint(s_layers, s_images, s_spacings)\n'



%% Cell type:code id:63997040-f4ba-47e5-9197-16295b1920d3 tags:



``` python

shifted_s = shift_to_vdW_gap(

    structure= vdW_s.structure,

    layers= layers,

    images= images,

    spacings= spacings,

    atomic_layers=atomic_layers,

    atomic_spacings=atomic_spacings

)

s_layers, s_images, s_spacings, s_atomic_layers, s_atomic_spacings = get_vdW_layers(shifted_s)

print(s_layers, s_images, s_spacings, s_atomic_layers, s_atomic_spacings)

```



%% Output



    [[8, 1, 4, 5, 0, 9], [10, 3, 6, 7, 2, 11]] [[0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0]] [4.053074184176202, 2.935904452246059] [[8], [1], [4], [5], [0], [9], [10], [3], [6], [7], [2], [11]] [2.290107083779596, 2.994003526732217, 3.2876766256722334, 4.437926432997699, 4.731599531937719, 5.435495974890334, 8.371400427136393, 8.960167885072961, 9.429028389346314, 10.371894248570655, 10.840754752844, 11.429522210780569]



%% Cell type:code id:017e136b-3f35-44b3-aa60-80f0a9577ddd tags:



``` python

vdW_s.structure.lattice.matrix

```



%% Output



    array([[ 2.09568326, -3.62982989,  0.        ],

           [ 2.09568326,  3.62982989,  0.        ],

           [ 0.        ,  0.        , 30.99278275]])



%% Cell type:code id:8c292c02-1da5-460b-bdf6-944dd1993af0 tags:



``` python

```

vdW_structures/vdW_heterostructure.py

+20 −2
Original line number Diff line number Diff line 
from vdW_structures.vdW_structure import VdWStructure

from vdW_structures.unique_structures import UniqueStructureGetter

from typing import Union, List

from scipy.linalg import sqrtm

import numpy as np

import math

import os
@@ -17,6 +18,21 @@ class VdWHeterostructureGenerator(): 
        self.sm = StructureMatcher(**kwargs)

        pass



    def polar_decomposition(self, film_sl_vectors, substrate_sl_vectors):

        s1, s2 = substrate_sl_vectors[0][:2], substrate_sl_vectors[1][:2]

        f1, f2 = film_sl_vectors[0][:2], film_sl_vectors[1][:2]



        S = np.column_stack([s1, s2])   # shape (2,2)

        F = np.column_stack([f1, f2])   # shape (2,2)



        M = F @ np.linalg.inv(S)

        U = sqrtm(M.T @ M)    # shape (2,2), symmetric by construction

        R = M @ np.linalg.inv(U)

        angle_rad = np.arctan2(R[1,0], R[0,0])

        angle_deg = np.degrees(angle_rad)



        return angle_deg



    def compute_signed_film_substrate_angle(self, film_sl_vectors, substrate_sl_vectors):

        """

        Compute the signed angle between the film and substrate motifs, distinguishing rotations.
@@ -113,10 +129,12 @@ class VdWHeterostructureGenerator(): 
            interface, minimum_vdW_gap) for interface in interfaces_list]

        heterostructure_properties = [interface.interface_properties for interface in interfaces_list]



        # Compute the signed angle between interfaces and add to heterostructure_properties dictionary

        # Compute the signed and decomposition angles between interfaces and add to heterostructure_properties dictionary

        for heterostructure_property in heterostructure_properties:

            heterostructure_property['signed_angle'] = self.compute_signed_film_substrate_angle(heterostructure_property['film_sl_vectors'], 

                                                                                                heterostructure_property['substrate_sl_vectors'])  

            heterostructure_property['decomposition_angle'] = self.polar_decomposition(heterostructure_property['film_sl_vectors'],

                                                                                                heterostructure_property['substrate_sl_vectors'])



        return self.sort_by_structure_length(unique_vdW_heterostructures, heterostructure_properties)