Updated example.ipynb and added new functionality to VdWHeterostructureGenerator (5bc552e8) · Commits · Morelock, Ryan / vdW_structures

tests/.ipynb_checkpoints/example-checkpoint.ipynb

+27 −38

Original line number	Diff line number	Diff line
		%% Cell type:code id:5647de5b-5ac4-41a1-b215-b319b8a75844 tags:

		``` python
		# package imports
		from pymatgen.analysis.interfaces.zsl import ZSLGenerator
		from pymatgen.core.structure import Structure
		from ase.visualize import view
		import numpy as np
		from datetime import datetime

		# class imports
		from vdW_structures.vdW_structure import VdWStructure
		from vdW_structures.vdW_heterostructure import VdWHeterostructureGenerator
		from vdW_structures.layer_shifter import LayerShifter
		from vdW_structures.unique_structures import UniqueStructureGetter
		```

		%% Cell type:markdown id:98fe58e1-bd90-43ee-8c96-cdfac46bf82d tags:

		## Input Bulk Structures

		%% Cell type:code id:4e69fc4f-4bee-4553-9018-1de3d3b190d8 tags:

		``` python
		bise2 = Structure.from_file('original_structures/BiSe2.vasp')
		bi2se3 = Structure.from_file('original_structures/Bi2Se3.vasp')
		bi3se4 = Structure.from_file('original_structures/Bi3Se4.vasp')
		bi4se3 = Structure.from_file('original_structures/Bi4Se3.vasp')
		test = bi3se4
		```

		%% Cell type:markdown id:039df08c-8f63-4d5a-8857-6bd2866a0e59 tags:

		## Generate Starting vdW Subunits for BiSe2, Bi2Se3, Bi3Se4

		%% Cell type:code id:fb3ed4e8-ecf3-4d15-a168-e9a979b8b20b tags:

		``` python
		vdW_bi2se3 = VdWStructure(bi2se3, minimum_vdW_gap=2.4)
		vdW_bi3se4 = VdWStructure(bi3se4, minimum_vdW_gap=2.4)
		vdW_bise2 = VdWStructure(bise2, minimum_vdW_gap=1.9)

		# Note: the vdW_bise2 structure has a smaller minimum vdW gap than vdW_bi2se3 and vdW_bi3se4, which could cause issues
		# To make this VdWStructure commensurate, create a vdW_gap in bise2 that is > 2.4

		#This can be done two ways (equivalent for this system):

		# Add vacuum to vdW_bise2. This works because bise2 is a single vdW unit.
		vdW_bise2_first = VdWStructure(vdW_bise2.set_vacuum(vacuum=2.41).structure, minimum_vdW_gap=2.4)
		print(vdW_bise2_first.vdW_layers, vdW_bise2_first.vdW_spacings)

		# Or, add padding to vdW layer 0 in vdW_bise2. This is how to change layer distances when multiple layers exist in a structure.
		vdW_bise2_second = VdWStructure(vdW_bise2.z_solve_shift_vdW_layers(layer_inds=0,
		solve_shift=2.41,
		adjust_z_cell=True).structure,
		minimum_vdW_gap=2.4)
		print(vdW_bise2_second.vdW_layers, vdW_bise2_second.vdW_spacings)
		```

		%% Output

		[[1, 0, 2]] [2.41]
		[[1, 0, 2]] [2.41]

		%% Cell type:code id:a2e9d9d3-16bb-4741-a8fa-222e7fc03ab8 tags:

		``` python
		# Look at the arrangement of vdW layers (by structure indices) and interlayer distances between the layers
		print(vdW_bi2se3.vdW_layers, vdW_bi2se3.vdW_spacings)
		print(vdW_bi3se4.vdW_layers, vdW_bi3se4.vdW_spacings)
		print(vdW_bise2_second.vdW_layers, vdW_bise2_second.vdW_spacings)
		```

		%% Output

		[[10, 3, 9, 4, 14], [13, 5, 12, 0, 8], [7, 1, 6, 2, 11]] [2.9359044522460795, 2.9359044522460938, 2.9359044522460884]
		[[9, 1, 12, 3, 15, 5, 14], [13, 4, 16, 6, 19, 8, 18], [17, 7, 20, 0, 11, 2, 10]] [3.0050793718446416, 3.0050793718446407, 3.0050793718446407]
		[[1, 0, 2]] [2.41]

		%% Cell type:code id:e850ce93-92fb-4e41-b253-31feeaaee138 tags:

		``` python
		# View the vdW_bise2_second structure
		view(vdW_bise2_second.structure.to_ase_atoms(), viewer='ngl')
		```

		%% Output


		HBox(children=(NGLWidget(), VBox(children=(Dropdown(description='Show', options=('All', 'Bi', 'Se'), value='Al…

		%% Cell type:code id:f2231b75-695b-4efd-be44-d91bfcf0e2e4 tags:

		``` python
		# Choose which layers you'd like to include in the heterostructure construction
		# For the current example, it's sufficient to use a single vdW layer from each structure

		subunit_vdW_bi2se3_0 = vdW_bi2se3.extend_vdW_layers([0])
		subunit_vdW_bi3se4_1 = vdW_bi3se4.extend_vdW_layers([1])
		subunit_vdW_bise2_0 = vdW_bise2_second.extend_vdW_layers([0])
		```

		%% Cell type:code id:237268a4-e426-43b5-8a5e-01e64b9bb055 tags:

		``` python
		# Look at the Bi2Se3 subunit, with the appropriate vdW gap across the PBC
		view(subunit_vdW_bi2se3_0.structure.to_ase_atoms(), viewer='ngl')
		```

		%% Output

		HBox(children=(NGLWidget(), VBox(children=(Dropdown(description='Show', options=('All', 'Bi', 'Se'), value='Al…

		%% Cell type:markdown id:54ad7b33-ecd1-4ec3-b5e7-25335651a5e8 tags:

		## Generate vdW Heterostructures

		%% Cell type:code id:3e53c9c7-8847-4bce-b865-1681a1c6fdde tags:

		``` python
		# This requires you to specify a ZSLGenerator structure, with the tolerances passed dictating the size(s) of vdW heterostructures generated

		# This ZSL object yields four unique heterostructures for Bi2Se3 on Bi3Se4

		bi2se3_bi3se4_zsl = ZSLGenerator(max_area_ratio_tol=0.5,
		max_area=150,
		max_area=100,
		max_length_tol=0.05,
		max_angle_tol=0.05)

		vsg = VdWHeterostructureGenerator()
		start = datetime.now()
		vdW_heterostructures_bi2se3_bi3se4 = vsg.generate_unique_vdW_heterostructures(film=subunit_vdW_bi2se3_0,
		vdW_heterostructures_bi2se3_bi3se4, bi2se3_bi3se4_heterostructure_properties = vsg.generate_vdW_heterostructures(film=subunit_vdW_bi2se3_0,
		substrate=subunit_vdW_bi3se4_1,
		zsl_generator=bi2se3_bi3se4_zsl,
		n_chunks=4,
		ltol=0.2, stol=0.3, angle_tol=5)
		n_chunks=4, get_unique=True,
		ltol=0.001, stol=0.001, angle_tol=0.01)
		end = datetime.now()
		print(end - start)
		# Print the lengths of the unique heterostructures
		s_vdW_heterostructures_bi2se3_bi3se4 = sorted(vdW_heterostructures_bi2se3_bi3se4, key=lambda obj: len(obj.structure))


		print([len(heterostructure.structure) for heterostructure in s_vdW_heterostructures_bi2se3_bi3se4])
		print([len(heterostructure.structure) for heterostructure in vdW_heterostructures_bi2se3_bi3se4])
		print([heterostructure_properties['signed_angle'] for heterostructure_properties in bi2se3_bi3se4_heterostructure_properties])
		```

		%% Output

		/Users/rym/miniconda3/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['wyckoffs']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
		warnings.warn(
		/Users/rym/miniconda3/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['equivalent_atoms']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
		warnings.warn(

		Unique structures: 2
		Unique structures: 2
		Unique structures: 4
		Unique structures: 4
		Unique structures: 4
		0:00:12.435890
		[12, 24, 84, 84]
		Unique structures: 9
		Unique structures: 10
		Unique structures: 11
		Unique structures: 11
		Unique structures: 11
		0:00:15.889802
		[12, 24, 24, 24, 24, 24, 24, 36, 48, 60, 72]
		[-2.8029520988333133e-15, -60.00000000000001, -60.00000000000001, -120.00000000000001, 60.00000000000001, 120.00000000000001, 2.8029520988333133e-15, 120.00000000000001, 3.7372694651110844e-15, 60.00000000000001, 0.0]

		%% Cell type:code id:43caf0bd-7abe-42e7-8fba-5f1146767bf0 tags:

		``` python
		# We'll use the first heterostructure returned, which has the fewest # of atoms
		view(vdW_heterostructures_bi2se3_bi3se4[0].structure.to_ase_atoms(), viewer='ngl')
		view(vdW_heterostructures_bi2se3_bi3se4[1].structure.to_ase_atoms(), viewer='ngl')
		```

		%% Output

		HBox(children=(NGLWidget(), VBox(children=(Dropdown(description='Show', options=('All', 'Bi', 'Se'), value='Al…

		%% Cell type:code id:ce0f96a9-4fe7-4760-bc73-c7bfedc730c0 tags:

		``` python
		# Add a BiSe2 film to this substrate
		bise2_bi2se3_bi3se4_zsl = ZSLGenerator(max_area_ratio_tol=0.5,
		max_area=150,
		max_length_tol=0.05,
		max_angle_tol=0.05)

		vdW_heterostructures_bise2_bi2se3_bi3se4 = vsg.generate_unique_vdW_heterostructures(film=subunit_vdW_bise2_0,
		vdW_heterostructures_bise2_bi2se3_bi3se4, bise2_bi2se3_bi3se4_heterostructure_properties = vsg.generate_vdW_heterostructures(
		film=subunit_vdW_bise2_0,
		substrate=vdW_heterostructures_bi2se3_bi3se4[0],
		zsl_generator=bise2_bi2se3_bi3se4_zsl)

		# Print the lengths of the unique heterostructures
		print([len(heterostructure.structure) for heterostructure in vdW_heterostructures_bise2_bi2se3_bi3se4])
		```

		%% Output

		/Users/rym/miniconda3/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['wyckoffs']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
		warnings.warn(
		/Users/rym/miniconda3/lib/python3.9/site-packages/spglib/spglib.py:115: DeprecationWarning: dict interface (SpglibDataset['equivalent_atoms']) is deprecated.Use attribute interface ({self.__class__.__name__}.{key}) instead
		warnings.warn(

		Unique structures: 6
		Unique structures: 6
		Unique structures: 7
		Unique structures: 7
		Unique structures: 7
		[15, 30, 30, 45, 45, 60, 105]
		Unique structures: 2
		Unique structures: 2
		Unique structures: 3
		Unique structures: 3
		Unique structures: 3
		[15, 45, 105]

		%% Cell type:code id:8957dbcd-3039-4927-b582-c29b373eec41 tags:

		``` python
		# And view the first heterostructure generated
		view(vdW_heterostructures_bise2_bi2se3_bi3se4[0].structure.to_ase_atoms(), viewer='ngl')
		```

		%% Output

		HBox(children=(NGLWidget(), VBox(children=(Dropdown(description='Show', options=('All', 'Bi', 'Se'), value='Al…

		%% Cell type:code id:d44dfd9c-1d4a-41a9-9ae4-76d4fab4b50b tags:

		``` python
		# Generate a cubic supercell for easier viewing

		from pymatgen.transformations.advanced_transformations import CubicSupercellTransformation
		cst = CubicSupercellTransformation(force_90_degrees=True)

		# Choose the heterostruture at index 3
		index = 1
		index = 0
		sc_heterostructure = VdWStructure(cst.apply_transformation(vdW_heterostructures_bise2_bi2se3_bi3se4[index].structure),
		vdW_heterostructures_bise2_bi2se3_bi3se4[index].minimum_vdW_gap)
		```

		%% Cell type:code id:63157904-a7fa-4597-baf4-f39e1b374c79 tags:

		``` python
		view(sc_heterostructure.structure.to_ase_atoms(), viewer='ngl')
		```

		%% Output

		HBox(children=(NGLWidget(), VBox(children=(Dropdown(description='Show', options=('All', 'Bi', 'Se'), value='Al…

		%% Cell type:code id:d73e497e-abd2-4f7b-aad0-4dec4159e287 tags:

		``` python
		len(sc_heterostructure.structure)
		```

		%% Output

		240

		%% Cell type:markdown id:c9f2ac4e-4d08-411b-8f74-d753efe3e621 tags:

		## Modifying the Heterostructure

		%% Cell type:markdown id:a4a15e07-dcc8-45dc-8e27-f4e5d78097b7 tags:

		### Shift a layer within the fixed box

		%% Cell type:code id:a14c9758-94a5-4159-8fb4-df8636c499df tags:

		``` python
		# Can be performed for x, y , or z directions, but we'll shift relative to z because it's easiest to visualize

		z_shifted_layer_1 = sc_heterostructure.z_shift_vdW_layers(layer_inds=1,
		shift=-1.0,
		coords_are_cartesian=True)
		view(z_shifted_layer_1.structure.to_ase_atoms(), viewer='ngl')
		```

		%% Output

		HBox(children=(NGLWidget(), VBox(children=(Dropdown(description='Show', options=('All', 'Bi', 'Se'), value='Al…

		%% Cell type:markdown id:e2d16126-4ae8-46c0-9984-5952312e7527 tags:

		### Increase interlayer distance, preserving other vdW spacings

		%% Cell type:code id:59ea4315-6bfa-4b52-b4da-7e59efa05252 tags:

		``` python
		# Increase a layer and increase box to preserve vdW spacing
		increased_gap_layer_2 = sc_heterostructure.z_solve_shift_vdW_layers(layer_inds=[2, 0],
		solve_shift=[2.8, 3.4],
		adjust_z_cell=True)
		view(increased_gap_layer_2.structure.to_ase_atoms(), viewer='ngl')
		```

		%% Output

		HBox(children=(NGLWidget(), VBox(children=(Dropdown(description='Show', options=('All', 'Bi', 'Se'), value='Al…

		%% Cell type:code id:f6ae9624-8e65-486c-a097-78dc06ff9516 tags:

		``` python
		print(increased_gap_layer_2.vdW_spacings)
		```

		%% Output

		[3.400000000000006, 2.9359044522460724, 2.7999999999999936]
		[3.4000000000000026, 2.9359044522460813, 2.7999999999999865]

		%% Cell type:markdown id:beb807a4-5527-44c4-b621-8ceee62f6dc3 tags:

		### Change the layering, so that it is Bi3Se4 - Bi2Se3 - BiSe2 - Bi2Se3

		%% Cell type:code id:76d06136-4e29-4322-8ba0-d52fa2da6f33 tags:

		``` python
		modified_layering = sc_heterostructure.extend_vdW_layers([0, 1, 2, 1])
		view(modified_layering.structure.to_ase_atoms(), viewer='ngl')
		```

		%% Output

		HBox(children=(NGLWidget(), VBox(children=(Dropdown(description='Show', options=('All', 'Bi', 'Se'), value='Al…

		%% Cell type:code id:295e9f76-e8cc-4c12-952a-0d74369aac8b tags:

		``` python
		print(modified_layering.vdW_spacings)
		```

		%% Output

		[3.005079371844629, 2.935904452246076, 2.4099999999999824, 2.93590445224606]
		[3.005079371844629, 2.935904452246074, 2.4099999999999824, 2.93590445224606]

		%% Cell type:code id:e0a0322a-fdc5-4e98-b299-9a9012ffeaed tags:
		%% Cell type:code id:b8d8d764-a293-41b7-a8dc-57caa5093910 tags:

		``` python
		```

tests/example.ipynb

+54 −187

File changed.

Preview size limit exceeded, changes collapsed.

vdW_structures/.vdW_heterostructure.py.swp

0 → 100644

+24 KiB

File added.

No diff preview for this file type.

View file

vdW_structures/vdW_heterostructure.py

+82 −3

Original line number	Diff line number	Diff line
		@@ -3,6 +3,8 @@ from vdW_structures.unique_structures import UniqueStructureGetter
		from typing import Union, List
		import numpy as np
		import math
		import os
		import json
		from pymatgen.core.structure import Structure
		from pymatgen.analysis.structure_matcher import StructureMatcher
		from pymatgen.core.surface import SlabGenerator
		@@ -10,12 +12,44 @@ from pymatgen.analysis.interfaces.zsl import ZSLGenerator
		from pymatgen.analysis.interfaces.coherent_interfaces import CoherentInterfaceBuilder



		class VdWHeterostructureGenerator():
		def __init__(self, **kwargs):
		self.sm = StructureMatcher(**kwargs)
		pass

		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.

		Args:
		film_sl_vectors (list of list): Superlattice vectors of the film (2D list of 3D vectors).
		substrate_sl_vectors (list of list): Superlattice vectors of the substrate (2D list of 3D vectors).

		Returns:
		float: Signed angle in degrees between the film and substrate repeating motifs.
		"""
		# Convert to numpy arrays and take the first in-plane vectors
		film_v1 = np.array(film_sl_vectors[0][:2]) # Only consider x, y components
		sub_v1 = np.array(substrate_sl_vectors[0][:2])

		# Compute angle using atan2 for proper signed angle determination
		dot_product = np.dot(film_v1, sub_v1)
		cross_product = np.cross(np.append(film_v1, 0), np.append(sub_v1, 0)) # Append z=0 for cross product
		angle_radians = np.arctan2(cross_product[2], dot_product) # Use z-component of cross product

		# Convert to degrees
		angle_degrees = np.degrees(angle_radians)

		return angle_degrees # Returns values from -180 to 180 degrees

		def sort_by_structure_length(self, vdW_heterostructures: List[VdWStructure], heterostructure_properties: List):
		'''Sorts structures by the number of atoms present'''
		zipped_heterostructures = zip(vdW_heterostructures, heterostructure_properties)
		sorted_zipped_heterostructures = sorted(zipped_heterostructures, key=lambda x: len(x[0].structure))
		s_vdW_heterostructures, s_heterostructure_properties = zip(*sorted_zipped_heterostructures)

		return s_vdW_heterostructures, s_heterostructure_properties

		def get_unique_structures(self, structures: List[Structure], **kwargs):
		'''Function that returns only the unique structures (based on PMG StructureMatcher) in a list of pymatgen Structure objects'''
		sm = StructureMatcher(**kwargs)
		@@ -32,7 +66,7 @@ class VdWHeterostructureGenerator():
		unique_structures.append(structure)
		return unique_structures

		def generate_unique_vdW_heterostructures(self,
		def generate_vdW_heterostructures(self,
		film: VdWStructure,
		substrate: VdWStructure,
		zsl_generator: ZSLGenerator,
		@@ -77,5 +111,50 @@ class VdWHeterostructureGenerator():
		minimum_vdW_gap = np.min(film.vdW_spacings + substrate.vdW_spacings + [interface_spacing]) - threshold_tolerance
		unique_vdW_heterostructures = [VdWStructure(
		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
		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'])

		return self.sort_by_structure_length(unique_vdW_heterostructures, heterostructure_properties)

		def write_heterostructures_data(self, directory, vdW_heterostructures, heterostructure_properties):

		def convert_ndarrays_to_lists(d):
		"""
		Recursively convert all numpy ndarray objects in a dictionary to lists.

		Args:
		d (dict): The input dictionary.

		Returns:
		dict: The dictionary with all np.ndarray objects converted to lists.
		"""
		if isinstance(d, dict):
		return {key: convert_ndarrays_to_lists(value) for key, value in d.items()}
		elif isinstance(d, list):
		return [convert_ndarrays_to_lists(item) for item in d]
		elif isinstance(d, np.ndarray):
		return d.tolist() # Convert numpy array to list
		else:
		return d # Return other types unchanged

		len_dct = {}
		for i, vdW_heterostructure in enumerate(vdW_heterostructures):
		length = len(vdW_heterostructure.structure)
		if length in len_dct.keys():
		len_dct[length] += 1
		else:
		len_dct[length] = 0

		write_path = os.path.join(directory, str(length), str(len_dct[length]))
		os.makedirs(write_path, exist_ok=True)
		vdW_heterostructure.structure.to(filename=os.path.join(write_path, 'POSCAR'))

		with open(os.path.join(write_path, 'interface.json'), 'w') as f:
		json.dump(convert_ndarrays_to_lists(heterostructure_properties[i]),
		f, indent=4)

		return unique_vdW_heterostructures
		return

vdW_structures/vdW_structure.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -223,7 +223,7 @@ class VdWStructure:
		]))

		for new_layer_number_ind, new_layer_number in enumerate(new_layer_numbers):
		print(self.vdW_layers[new_layer_number])
		#print(self.vdW_layers[new_layer_number])
		for site_number in self.vdW_layers[new_layer_number]:
		site_shift = np.multiply(new_layer_number_ind, self.structure.lattice.matrix[2])
		coords = self.structure[site_number].coords + site_shift