import sys
import numpy as np
import math
import copy
import os
import shutil
import gmsh
from geom.classes import molecule, parameters
from geom.functions import general, output, tools
from ase.cluster.cubic import FaceCenteredCubic, BodyCenteredCubic
from ase.build import graphene_nanoribbon
from ase.build import graphene as graphene_general_ase
from ase.io import write
# -------------------------------------------------------------------------------------
[docs]
def select_case(inp):
"""
Selects the appropriate geometry creation function based on user input.
Args:
inp (input_class): An instance of the input class containing user-defined parameters.
Returns:
None: Executes the corresponding geometry generation function.
"""
if (inp.gen_graphene): graphene(inp)
if (inp.gen_sphere): sphere(inp)
if (inp.gen_sphere_core_shell): sphere_core_shell(inp)
if (inp.gen_3d_mesh_sphere): sphere_3d_mesh(inp)
if (inp.gen_rod): rod(inp)
if (inp.gen_rod_core_shell): rod_core_shell(inp)
if (inp.gen_3d_mesh_rod): rod_3d_mesh(inp)
if (inp.gen_tip): tip(inp)
if (inp.gen_pyramid): pyramid(inp)
if (inp.gen_pentpyramid): pentpyramid(inp)
if (inp.gen_pentbipyramid): pentbipyramid(inp)
if (inp.gen_cone): cone(inp)
if (inp.gen_microscope): microscope(inp)
if (inp.gen_icosahedra): icosahedra(inp)
if (inp.gen_cto): cto(inp)
if (inp.gen_idh): idh(inp)
if (inp.gen_bipyramid): bipyramid(inp)
if (inp.gen_pencil): pencil(inp)
# Creation of dimer and bowtie structures
if (inp.create_dimer): tools.create_dimer(inp)
if (inp.create_bowtie): tools.create_bowtie(inp)
# Eliminate tmp folder containing bulk structure
if inp.create_ase_bulk: shutil.rmtree(inp.tmp_folder)
# -------------------------------------------------------------------------------------
[docs]
def graphene(inp):
"""
Generates a graphene structure (ribbon, disk, ring, or triangle).
Args:
inp (input_class): An instance containing the input parameters.
Returns:
None: Saves the generated geometry to a file.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Initialize bulk "molecule" and read geometry
mol = molecule.molecule()
mol.read_geom(inp.geom_file,False)
# Pick only atoms within the defined graphene geometry
if inp.graphene_structure=='rib': mol.filter_xyz_graphene_to_ribbon(inp)
if inp.graphene_structure=='disk': mol.filter_xyz_graphene_to_disk(inp)
if inp.graphene_structure=='ring': mol.filter_xyz_graphene_to_ring(inp)
if inp.graphene_structure=='triangle': mol.filter_xyz_graphene_to_triangle(inp)
# Remove dangling bonds
mol.remove_dangling_bonds_graphene(inp)
# Translate geometrical center to 000
mol.trans_geom_center_to_000()
# Save filtered geometry
if inp.graphene_structure=='rib': inp.xyz_output = f'graphene_ribbon_{inp.X_length}_{inp.Y_length}'
if inp.graphene_structure=='disk': inp.xyz_output = f'graphene_disk_{inp.radius}'
if inp.graphene_structure=='ring': inp.xyz_output = f'graphene_ring_out_{inp.radius_out}_in_{inp.radius_in}'
if inp.graphene_structure=='triangle': inp.xyz_output = f'graphene_triangle_{inp.graphene_edge_type}_{inp.side_length}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def sphere(inp):
"""
Generates a spherical nanoparticle geometry.
Args:
inp (input_class): An instance containing the input parameters.
Returns:
None: Saves the generated sphere geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Initialize bulk "molecule" and read geometry
mol = molecule.molecule()
mol.read_geom(inp.geom_file,False)
# Pick only atoms within the defined sphere
mol.filter_xyz_in_sphere(inp)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'sphere_{inp.atomtype}_r_{inp.radius}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def sphere_core_shell(inp):
"""
Generates a core-shell sphere structure.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the core-shell structure to a file.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Extract merge cutoff
param = parameters.parameters()
inp.merge_cutoff = param.min_dist.get(inp.atomtype)
# Initialize bulk "molecule" and read geometry
mol_out = molecule.molecule()
mol_out = mol_out.read_geom(inp.geom_file,False)
# Pick only atoms within the defined sphere
inp.radius = inp.radius_out
inp.atomtype = inp.atomtype_out
mol_out.filter_xyz_in_sphere(inp)
# Copy mol_out as initial guess, change atomtype, and create smaller sphere
mol_in = copy.deepcopy(mol_out)
mol_in.change_atomtype(inp.atomtype_in)
inp.radius = inp.radius_in
inp.atomtype = inp.atomtype_in
mol_in.filter_xyz_in_sphere(inp)
# Create shell by subtracting core geometry
mol_shell = tools.subtract_geoms(inp,mol_in,mol_out)
# Alloy core and shell
if inp.alloy:
inp.alloy_string = f"_alloy_{inp.alloy_perc}_perc"
inp.atomtype = inp.atomtype_out
inp.atomtype_alloy = inp.atomtype_in
mol_shell.create_alloy(inp)
inp.atomtype = inp.atomtype_in
inp.atomtype_alloy = inp.atomtype_out
mol_in.create_alloy(inp)
# Merge to create core-shell
mol_core_shell = tools.merge_geoms(inp,mol_in,mol_shell)
# Save filtered geometry
inp.xyz_output = f'sphere_core_{inp.atomtype_in}_r_{inp.radius_in}_shell_{inp.atomtype_out}_r_{inp.radius_out}{inp.alloy_string}'
output.print_geom(mol_core_shell, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def sphere_3d_mesh(inp):
"""
Generates a 3D mesh representation of a sphere using Gmsh in MeshFormat 2.2.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated 3D mesh structure.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Initialize Gmsh
gmsh.initialize([str(inp.radius), str(inp.mesh_size), str(inp.mesh_output)])
gmsh.model.add("sphere")
# Create a sphere centered at (0,0,0) with the given radius
sphere = gmsh.model.occ.addSphere(0, 0, 0, inp.radius)
# Synchronize the model with the OpenCASCADE representation
gmsh.model.occ.synchronize()
# Define mesh element size for the sphere
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), inp.mesh_size)
# Generate 3D mesh (only surface)
gmsh.model.mesh.generate(2)
# Save in MeshFormat 2.2
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.write(inp.mesh_output)
# Finalize Gmsh
gmsh.finalize()
# -------------------------------------------------------------------------------------
[docs]
def rod(inp):
"""
Generates a cylindrical rod geometry.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated rod geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Extract merge cutoff
param = parameters.parameters()
inp.merge_cutoff = param.min_dist.get(inp.atomtype)
# Initialize bulk "molecule" and read geometry
mol_sphere_1 = molecule.molecule()
mol_sphere_2 = molecule.molecule()
mol_cylinder = molecule.molecule()
mol_sphere_1.read_geom(inp.geom_file,False)
mol_sphere_2.read_geom(inp.geom_file,False)
mol_cylinder.read_geom(inp.geom_file,False)
# Create individual sphere at the extremes of the rods
tools.determine_sphere_center(inp,'+')
mol_sphere_1.filter_xyz_in_sphere(inp)
tools.determine_sphere_center(inp,'-')
mol_sphere_2.filter_xyz_in_sphere(inp)
# Create cylinder
mol_cylinder.filter_xyz_in_cylinder(inp)
# -------------------------------------
# Merge cylinder and spheres geometries
# -------------------------------------
mol_tmp = tools.merge_geoms(inp,mol_sphere_1,mol_sphere_2)
mol_rod = tools.merge_geoms(inp,mol_cylinder,mol_tmp)
# Alloy
if inp.alloy: mol_rod.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'rod_{inp.atomtype}_{inp.main_axis.upper()}_l_{inp.rod_length}_w_{inp.rod_width}{inp.alloy_string}'
output.print_geom(mol_rod, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def rod_core_shell(inp):
"""
Generates a core-shell rod structure.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated core-shell rod geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Extract merge cutoff
param = parameters.parameters()
inp.merge_cutoff = param.min_dist.get(inp.atomtype)
# Initialize bulk "molecule" and read geometry
mol_sphere_1 = molecule.molecule()
mol_sphere_2 = molecule.molecule()
mol_cylinder = molecule.molecule()
mol_sphere_1.read_geom(inp.geom_file,False)
mol_sphere_2.read_geom(inp.geom_file,False)
mol_cylinder.read_geom(inp.geom_file,False)
# ----------------------- #
# --------- Out --------- #
# Populate input class
inp.rod_length = inp.rod_length_out
inp.rod_width = inp.rod_width_out
inp.atomtype = inp.atomtype_out
# Create individual sphere at the extremes of the rods
tools.determine_sphere_center(inp,'+')
mol_sphere_1.filter_xyz_in_sphere(inp)
tools.determine_sphere_center(inp,'-')
mol_sphere_2.filter_xyz_in_sphere(inp)
# Create cylinder
mol_cylinder.filter_xyz_in_cylinder(inp)
# Merge cylinder and spheres geometries
mol_tmp = tools.merge_geoms(inp,mol_sphere_1,mol_sphere_2)
mol_out = tools.merge_geoms(inp,mol_cylinder,mol_tmp)
# ---------------------- #
# --------- In --------- #
# Copy mol_out as initial guess, change atomtype, and populate input class
mol_sphere_1 = copy.deepcopy(mol_out)
mol_sphere_1.change_atomtype(inp.atomtype_in)
mol_sphere_2 = copy.deepcopy(mol_out)
mol_sphere_2.change_atomtype(inp.atomtype_in)
inp.rod_length = inp.rod_length_in
inp.rod_width = inp.rod_width_in
inp.atomtype = inp.atomtype_in
# Create individual sphere at the extremes of the rods
tools.determine_sphere_center(inp,'+')
mol_sphere_1.filter_xyz_in_sphere(inp)
tools.determine_sphere_center(inp,'-')
mol_sphere_2.filter_xyz_in_sphere(inp)
# Create cylinder
mol_cylinder.filter_xyz_in_cylinder(inp)
# Merge cylinder and spheres geometries
mol_tmp = tools.merge_geoms(inp,mol_sphere_1,mol_sphere_2)
mol_in = tools.merge_geoms(inp,mol_cylinder,mol_tmp)
# Create shell by subtracting core geometry
mol_shell = tools.subtract_geoms(inp,mol_in,mol_out)
# Alloy core and shell
if inp.alloy:
inp.alloy_string = f"_alloy_{inp.alloy_perc}_perc"
inp.atomtype = inp.atomtype_out
inp.atomtype_alloy = inp.atomtype_in
mol_shell.create_alloy(inp)
inp.atomtype = inp.atomtype_in
inp.atomtype_alloy = inp.atomtype_out
mol_in.create_alloy(inp)
# Merge to create core-shell
mol_core_shell = tools.merge_geoms(inp,mol_in,mol_shell)
# Save filtered geometry
inp.xyz_output = f'rod_{inp.main_axis.upper()}_core_{inp.atomtype_in}_l_{inp.rod_length_in}_r_{inp.rod_width_in}_shell_{inp.atomtype_out}_l_{inp.rod_length_out}_r_{inp.rod_width_out}_shell_{inp.atomtype_out}{inp.alloy_string}'
output.print_geom(mol_core_shell, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def rod_3d_mesh(inp):
"""
Generates a 3D mesh representation of a rod using Gmsh in MeshFormat 2.2.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated 3D mesh structure.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Initialize Gmsh
gmsh.initialize([])
gmsh.model.add("rod")
# Define rod parameters
radius = inp.rod_width / 2.0 # Radius of spheres and cylinder
length = inp.rod_length - inp.rod_width # Total rod length
# Mapping for axes
axis_map = {'x': (length, 0, 0), 'y': (0, length, 0), 'z': (0, 0, length)}
axis = axis_map[inp.main_axis] # Choose correct alignment
# Define sphere center positions at rod ends
sphere1_pos = tuple(-coord / 2 for coord in axis)
sphere2_pos = tuple(coord / 2 for coord in axis)
# Create spheres at the ends
sphere1 = gmsh.model.occ.addSphere(*sphere1_pos, radius)
sphere2 = gmsh.model.occ.addSphere(*sphere2_pos, radius)
# Create the connecting cylinder **starting from the first sphere center**
cylinder_start = sphere1_pos
cylinder_direction = tuple(coord for coord in axis) # Direction of cylinder
cylinder = gmsh.model.occ.addCylinder(*cylinder_start, *cylinder_direction, radius)
# Synchronize before merging
gmsh.model.occ.synchronize()
# **Fuse all parts into a single surface representation**
rod, _ = gmsh.model.occ.fuse([(3, sphere1), (3, sphere2)], [(3, cylinder)])
# Synchronize after fusion
gmsh.model.occ.synchronize()
# **Remove any internal surfaces or duplicate points**
gmsh.model.occ.removeAllDuplicates()
gmsh.model.occ.synchronize()
# **Apply mesh size factor**
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), inp.mesh_size)
# **Generate only a 2D surface mesh**
gmsh.model.mesh.generate(2)
# Save mesh in MeshFormat 2.2
gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
gmsh.write(inp.mesh_output)
# Finalize Gmsh
gmsh.finalize()
# -------------------------------------------------------------------------------------
[docs]
def tip(inp):
"""
Generates a nanostructured round tip geometry.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated tip geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Initialize bulk "molecule" and read geometry
mol = molecule.molecule()
mol.read_geom(inp.geom_file,False)
# Pick only atoms within the defined paraboloid
mol.filter_xyz_in_elliptic_paraboloid(inp)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'tip_{inp.atomtype}_elliptic_paraboloid_a-{inp.elliptic_parabola_a}_b-{inp.elliptic_parabola_b}_zmin-{inp.z_min}_zmax-{inp.z_max}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def pyramid(inp):
"""
Generates a pyramid-shaped geometry with a square base.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated pyramid geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Initialize bulk "molecule" and read geometry
mol = molecule.molecule()
mol.read_geom(inp.geom_file, False)
# Define vertices with respect to the center (C[0,0,0]) of the XY-plane base
#
# XY base (center 5 in the +z direction)
#
# 4 ______ 1 ^ y
# | | |
# | C | o ---> x
# |______| z
# 3 2
centers = {
"center_1": [inp.side_length / 2.0, inp.side_length / 2.0, 0.0],
"center_2": [inp.side_length / 2.0, -inp.side_length / 2.0, 0.0],
"center_3": [-inp.side_length / 2.0, -inp.side_length / 2.0, 0.0],
"center_4": [-inp.side_length / 2.0, inp.side_length / 2.0, 0.0],
"center_5": [0.0, 0.0, inp.z_max]
}
# Define planes to cut the pyramid
planes = {
"n_125": tools.calculate_normal_and_rhs(centers["center_1"], centers["center_2"], centers["center_5"]),
"n_235": tools.calculate_normal_and_rhs(centers["center_2"], centers["center_3"], centers["center_5"]),
"n_345": tools.calculate_normal_and_rhs(centers["center_3"], centers["center_4"], centers["center_5"]),
"n_415": tools.calculate_normal_and_rhs(centers["center_4"], centers["center_1"], centers["center_5"])
}
# Pick only atoms within the defined pyramid
mol.filter_xyz_in_pyramid(inp,centers, planes)
# Remove dangling atom on extreme
mol.remove_dangling_atoms_metals(inp)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'pyramid_{inp.atomtype}_length-{inp.side_length}_zmin-{inp.z_min}_zmax-{inp.z_max}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def pentpyramid(inp):
"""
Generates a pyramid-shaped geometry with a pentagonal base.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated pyramid geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
param = parameters.parameters()
inp.merge_cutoff = param.min_dist.get(inp.atomtype)
mol = molecule.molecule()
mol.read_geom(inp.geom_file, False)
nearest_distance = param.min_dist.get(inp.atomtype)
z_max_half = inp.z_max / 2.0
cap_radius = 5.0
cap_start = max(0.0, z_max_half - cap_radius)
effective_base_width = inp.base_width + 2.5 * nearest_distance
pent_angles = [math.radians(90.0 + i * 72.0) for i in range(5)]
xs = mol.xyz[0]
ys = mol.xyz[1]
zs = mol.xyz[2]
keep_mask = np.zeros(mol.nAtoms, dtype=bool)
for idx in range(mol.nAtoms):
x, y, z = xs[idx], ys[idx], zs[idx]
if z < 0.0 or z > z_max_half + 1e-8:
continue
layer_fraction = min(z / cap_start, 1.0) if cap_start > 0.0 else 1.0
layer_width = effective_base_width * (1.0 - layer_fraction) + 2.0 * cap_radius * layer_fraction
if z >= cap_start:
dz = z - cap_start
cap_width = 2.0 * math.sqrt(max(0.0, cap_radius ** 2 - dz ** 2))
layer_width = min(layer_width, cap_width)
if layer_width <= nearest_distance:
if np.sqrt(x ** 2 + y ** 2) < nearest_distance:
keep_mask[idx] = True
continue
R_layer = layer_width / (2.0 * math.sin(2.0 * math.pi / 5.0))
apothem = R_layer * math.cos(math.pi / 5.0)
corner_rounding = min(3.0 * nearest_distance, 0.18 * R_layer)
rounded_radius = max(apothem + 0.35 * nearest_distance, R_layer - corner_rounding)
vertices = np.array([[R_layer * math.cos(a), R_layer * math.sin(a)] for a in pent_angles])
point = np.array([x, y])
inside = True
for vi in range(5):
edge = vertices[(vi + 1) % 5] - vertices[vi]
rel = point - vertices[vi]
if edge[0] * rel[1] - edge[1] * rel[0] < -1e-8:
inside = False
break
if inside and np.linalg.norm(point) <= rounded_radius:
keep_mask[idx] = True
keep_idx = np.where(keep_mask)[0]
mol.xyz = mol.xyz[:, keep_idx]
mol.nAtoms = len(keep_idx)
mol.atoms = [inp.atomtype] * mol.nAtoms
mol.xyz_center = np.mean(mol.xyz, axis=1)
mol.xyz_max = np.max(mol.xyz, axis=1)
mol.xyz_min = np.min(mol.xyz, axis=1)
if mol.nAtoms == 0:
output.error("Pentagonal pyramid generation produced no atoms. Increase z_max or base_width.")
# Align base to z=0
z_min = mol.xyz_min[-1]
if z_min > 0.0:
mol.translate_geom(z_min, [0.0, 0.0, -1.0])
# Remove dangling atoms on top-half core
mol.remove_dangling_atoms_metals(inp)
# Rotate 180 degrees to create bottom pyramid and merge
mol_rot = molecule.molecule()
mol_rot = tools.rotate(mol, 180.0, '+x', mol_rot)
mol_rot = tools.rotate(mol_rot, 180.0, '+z', mol_rot)
mol = tools.merge_geoms(inp, mol, mol_rot)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'pentpyramid_{inp.atomtype}_width-{inp.base_width}_zmax-{inp.z_max}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def bipyramid(inp):
"""
Generates a bipyramid-shaped geometry with a pentagonal base.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated bipyramidal geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Extract lattice constant and merge cutoff
param = parameters.parameters()
lattice_constant = param.lattice_constant.get(inp.atomtype)
inp.merge_cutoff = param.min_dist.get(inp.atomtype)
# Initialize bulk "molecule" and read geometry
mol = molecule.molecule()
mol.read_geom(inp.geom_file, False)
# Define vertices with respect to the center (C[0,0,0]) of the XY-plane base
#
# XY base (center 6 in the +z direction)
#
#
# 1
# ,'.
# ,' `.
# ,' `.
# 5 ,' `. 2 ^ y
# \ C / |
# \ / o ---> x
# \ / z
# \_______/
# 4 3
#
# ------------------------------
# Define pentagonal-base pyramid
# -------------------------------
# Recover true circumradius
R = inp.bipyramid_width
# Calculate pentagon vertices (angles offset by 90° to have a flat side at bottom)
angles = [math.radians(90 + i * 72) for i in range(5)] # 72° between vertices
centers = {}
# Pentagon base vertices
for i in range(5):
centers[f"center_{i+1}"] = [
R * math.cos(angles[i]), # x coordinate
R * math.sin(angles[i]), # y coordinate
0.0 # z coordinate (base)
]
# Apex point at the top (half of total length from base to apex)
centers["center_6"] = [0.0, 0.0, inp.bipyramid_length ] # z coordinate (top apex). Added extra length to ensure lateral sides are filled by atoms
# Define planes to cut the pyramid
planes = {
"n_126": tools.calculate_normal_and_rhs(centers["center_1"], centers["center_2"], centers["center_6"]),
"n_236": tools.calculate_normal_and_rhs(centers["center_2"], centers["center_3"], centers["center_6"]),
"n_346": tools.calculate_normal_and_rhs(centers["center_3"], centers["center_4"], centers["center_6"]),
"n_456": tools.calculate_normal_and_rhs(centers["center_4"], centers["center_5"], centers["center_6"]),
"n_516": tools.calculate_normal_and_rhs(centers["center_5"], centers["center_1"], centers["center_6"]),
}
# Include atoms within the defined pentagonal pyramid
mol.filter_xyz_in_pentagonal_pyramid(inp,centers, planes)
# If the pyramid base is not on the XY plane, translate it
z_min = mol.xyz_min[-1]
if (z_min > 0.0 ):
mol.translate_geom(z_min,[0.0,0.0, -1.0])
# Rotate 180 degrees to create bottom pyramid and merge
mol_rot = molecule.molecule()
mol_rot = tools.rotate(mol,180.0,'+x',mol_rot)
mol_rot = tools.rotate(mol_rot,180.0,'+z',mol_rot)
mol = tools.merge_geoms(inp,mol,mol_rot)
# Remove dangling atom on extreme
mol.remove_dangling_atoms_metals(inp)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.bipyramid_width = inp.bipyramid_width * 2.0 # Make it match with initial definition
inp.bipyramid_length = inp.bipyramid_length * 2.0 # Make it match with initial definition
inp.xyz_output = f'bipyramid_{inp.atomtype}_width-{inp.bipyramid_width}_length-{inp.bipyramid_length}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def pentbipyramid(inp):
"""
Generates a pentagonal bipyramid core-shell geometry.
Core uses the pentpyramid atom-by-atom stacking approach (atomtype_in).
Shell is an outer rod, either fullshell or halfshell (atomtype_out).
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
# Extract merge cutoff for outer rod
param = parameters.parameters()
inp.merge_cutoff = param.min_dist.get(inp.atomtype_out)
# -----------------------------------------------
# Build outer rod (atomtype_out) — same as pencil
# -----------------------------------------------
mol_sphere_1 = molecule.molecule()
mol_sphere_2 = molecule.molecule()
mol_cylinder = molecule.molecule()
mol_sphere_1.read_geom(inp.geom_file, False)
mol_sphere_2.read_geom(inp.geom_file, False)
mol_cylinder.read_geom(inp.geom_file, False)
inp.atomtype = inp.atomtype_out
tools.determine_sphere_center(inp, '+')
mol_sphere_1.filter_xyz_in_sphere(inp)
tools.determine_sphere_center(inp, '-')
mol_sphere_2.filter_xyz_in_sphere(inp)
mol_cylinder.filter_xyz_in_cylinder(inp)
mol_tmp = tools.merge_geoms(inp, mol_sphere_1, mol_sphere_2)
mol_out = tools.merge_geoms(inp, mol_cylinder, mol_tmp)
# -----------------------------------------------
# Build inner core (atomtype_in) — bulk lattice positions, rounded pentpyramid shape
# -----------------------------------------------
mol_in = molecule.molecule()
mol_in.read_geom(inp.geom_file, False)
mol_in.change_atomtype(inp.atomtype_in)
inp.atomtype = inp.atomtype_in
nearest_distance = param.min_dist.get(inp.atomtype_in) # match pentpyramid shape parameters
z_max_half = inp.z_max / 2.0
cap_radius = 5.0
cap_start = max(0.0, z_max_half - cap_radius)
effective_base_width = inp.base_width + 2.5 * nearest_distance
pent_angles = [math.radians(90.0 + i * 72.0) for i in range(5)]
xs = mol_in.xyz[0]
ys = mol_in.xyz[1]
zs = mol_in.xyz[2]
keep_mask = np.zeros(mol_in.nAtoms, dtype=bool)
for idx in range(mol_in.nAtoms):
x, y, z = xs[idx], ys[idx], zs[idx]
if z < 0.0 or z > z_max_half + 1e-8:
continue
layer_fraction = min(z / cap_start, 1.0) if cap_start > 0.0 else 1.0
layer_width = effective_base_width * (1.0 - layer_fraction) + 2.0 * cap_radius * layer_fraction
if z >= cap_start:
dz = z - cap_start
cap_width = 2.0 * math.sqrt(max(0.0, cap_radius ** 2 - dz ** 2))
layer_width = min(layer_width, cap_width)
if layer_width <= nearest_distance:
if np.sqrt(x ** 2 + y ** 2) < nearest_distance:
keep_mask[idx] = True
continue
R_layer = layer_width / (2.0 * math.sin(2.0 * math.pi / 5.0))
apothem = R_layer * math.cos(math.pi / 5.0)
corner_rounding = min(3.0 * nearest_distance, 0.18 * R_layer)
rounded_radius = max(apothem + 0.35 * nearest_distance, R_layer - corner_rounding)
vertices = np.array([[R_layer * math.cos(a), R_layer * math.sin(a)] for a in pent_angles])
point = np.array([x, y])
inside = True
for vi in range(5):
edge = vertices[(vi + 1) % 5] - vertices[vi]
rel = point - vertices[vi]
if edge[0] * rel[1] - edge[1] * rel[0] < -1e-8:
inside = False
break
if inside and np.linalg.norm(point) <= rounded_radius:
keep_mask[idx] = True
keep_idx = np.where(keep_mask)[0]
mol_in.xyz = mol_in.xyz[:, keep_idx]
mol_in.nAtoms = len(keep_idx)
mol_in.atoms = [inp.atomtype_in] * mol_in.nAtoms
mol_in.xyz_center = np.mean(mol_in.xyz, axis=1)
mol_in.xyz_max = np.max(mol_in.xyz, axis=1)
mol_in.xyz_min = np.min(mol_in.xyz, axis=1)
# Align base to z=0 and shift outer rod by same amount
z_min = mol_in.xyz_min[-1]
if z_min > 0.0:
mol_in.translate_geom(z_min, [0.0, 0.0, -1.0])
mol_out.translate_geom(z_min, [0.0, 0.0, -1.0])
# If halfshell: slice outer rod to keep only top half (z >= -0.1)
if inp.pentbipyramid_type == "halfshell":
inp.atomtype = inp.atomtype_out
mol_out.slice_xyz_by_z_threshold(inp, z_threshold=-0.1)
# Rotate core 180° to create bottom pyramid and merge into full bipyramid core
mol_rot = molecule.molecule()
mol_rot = tools.rotate(mol_in, 180.0, '+x', mol_rot)
mol_rot = tools.rotate(mol_rot, 180.0, '+z', mol_rot)
mol_in = tools.merge_geoms(inp, mol_in, mol_rot)
# Remove dangling atoms on merged bipyramid core
inp.atomtype = inp.atomtype_in
mol_in.remove_dangling_atoms_metals(inp)
# Create shell by subtracting core from outer rod
mol_shell = tools.subtract_geoms(inp, mol_in, mol_out)
# Alloy
if inp.alloy:
inp.alloy_string = f"_alloy_{inp.alloy_perc}_perc"
inp.atomtype = inp.atomtype_out
inp.atomtype_alloy = inp.atomtype_in
mol_shell.create_alloy(inp)
inp.atomtype = inp.atomtype_in
inp.atomtype_alloy = inp.atomtype_out
mol_in.create_alloy(inp)
# Merge core + shell
mol_core_shell = tools.merge_geoms(inp, mol_in, mol_shell)
# Save geometry
inp.xyz_output = f'pentbipyramid_{inp.pentbipyramid_type}_core_{inp.atomtype_in}_shell_{inp.atomtype_out}_in_width-{inp.base_width}_in_zmax-{inp.z_max}_out_width-{inp.rod_width}_out_length-{inp.rod_length}{inp.alloy_string}'
output.print_geom(mol_core_shell, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def cone(inp):
"""
Generates a conical geometry.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated cone geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Initialize bulk "molecule" and read geometry
mol = molecule.molecule()
mol.read_geom(inp.geom_file, False)
# Pick only atoms within the defined cone
mol.filter_xyz_in_cone(inp)
# Remove dangling atom on extreme
mol.remove_dangling_atoms_metals(inp)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'cone_{inp.atomtype}_radius-{inp.radius}_zmin-{inp.z_min}_zmax-{inp.z_max}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def microscope(inp):
"""
Generates a simulated microscope tip structure with a pyramidal and paraboloidal combination.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated microscope tip geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Extract lattice constant and merge cutoff
param = parameters.parameters()
lattice_constant = param.lattice_constant.get(inp.atomtype)
inp.merge_cutoff = param.min_dist.get(inp.atomtype)
# ---------------------------------------------
# Create paraboloid from bulk metallic geometry
# ---------------------------------------------
inp.z_max = inp.z_max_paraboloid
mol_paraboloid = molecule.molecule()
mol_paraboloid.read_geom(inp.geom_file,False)
mol_paraboloid.filter_xyz_in_elliptic_paraboloid(inp)
# ------------------------------------------
# Create pyramid from bulk metallic geometry
# ------------------------------------------
inp.z_max = inp.z_max_pyramid
mol_pyramid = molecule.molecule()
mol_pyramid.read_geom(inp.geom_file, False)
centers = {
"center_1": [inp.side_length / 2.0, inp.side_length / 2.0, 0.0],
"center_2": [inp.side_length / 2.0, -inp.side_length / 2.0, 0.0],
"center_3": [-inp.side_length / 2.0, -inp.side_length / 2.0, 0.0],
"center_4": [-inp.side_length / 2.0, inp.side_length / 2.0, 0.0],
"center_5": [0.0, 0.0, inp.z_max]
}
# Define planes to cut the pyramid
# Eqs. --> n_x * x + n_y * y + n_z * z + rhs = 0
planes = {
"n_125": tools.calculate_normal_and_rhs(centers["center_1"], centers["center_2"], centers["center_5"]),
"n_235": tools.calculate_normal_and_rhs(centers["center_2"], centers["center_3"], centers["center_5"]),
"n_345": tools.calculate_normal_and_rhs(centers["center_3"], centers["center_4"], centers["center_5"]),
"n_415": tools.calculate_normal_and_rhs(centers["center_4"], centers["center_1"], centers["center_5"])
}
# Pick only atoms within the defined pyramid
mol_pyramid.filter_xyz_in_pyramid(inp,centers, planes)
# Remove dangling atom on extreme
mol_pyramid.remove_dangling_atoms_metals(inp)
# Rotate pyramid 180 degrees along z
mol_pyramid_rot = molecule.molecule()
mol_pyramid_rot = tools.rotate(mol_pyramid,180.0,'+x',mol_pyramid_rot)
# Move 1/2*reticular_distance (lc) up in the z axis
shift = math.ceil((inp.z_max_pyramid/2.0)/lattice_constant) * lattice_constant
mol_pyramid_rot.translate_geom(shift,[0.0,0.0, 1.0])
# -----------------------------------------
# Merge paraboloid and pyramidal geometries
# -----------------------------------------
mol_microscope = tools.merge_geoms(inp,mol_paraboloid,mol_pyramid_rot)
# Alloy
if inp.alloy: mol_microscope.create_alloy(inp)
inp.xyz_output = f'microscope_{inp.atomtype}_parabola_{inp.z_max_paraboloid}_{inp.elliptic_parabola_a}_{inp.elliptic_parabola_b}_pyramid_{inp.z_max_pyramid}_{inp.side_length}{inp.alloy_string}'
output.print_geom(mol_microscope, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def icosahedra(inp):
"""
Generates an icosahedral nanoparticle geometry.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated icosahedral geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Check FCC lattice on selected atom type (requirement)
general.check_FCC(inp.atomtype,'icosahedron')
# Initialize bulk "molecule" and create icosahedra with ASE
mol = molecule.molecule()
mol.create_icosahedra(inp)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'icosahedron_{inp.atomtype}_r_{inp.radius}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def cto(inp):
"""
Generates a cuboctahedral nanoparticle geometry.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated cuboctahedral geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Check FCC lattice on selected atom type (requirement)
general.check_FCC(inp.atomtype,'cuboctahedron')
# Initialize bulk "molecule" and create cuboctahedron with ASE
mol = molecule.molecule()
mol.create_cuboctahedra(inp)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'cuboctahedron_{inp.atomtype}_r_{inp.radius}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def idh(inp):
"""
Generates a decahedral nanoparticle geometry.
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated decahedral geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Check FCC lattice on selected atom type (requirement)
general.check_FCC(inp.atomtype,'decahedron')
# Initialize bulk "molecule" and create decahedron with ASE
mol = molecule.molecule()
mol.create_decahedra(inp)
# Alloy
if inp.alloy: mol.create_alloy(inp)
# Save filtered geometry
inp.xyz_output = f'decahedron_{inp.atomtype}_r_{inp.radius}{inp.alloy_string}'
output.print_geom(mol, inp.xyz_output)
# -------------------------------------------------------------------------------------
[docs]
def pencil(inp):
"""
Generates a nanopencil nanoparticle geometry.
(Reference of experimental geometries: https://doi.org/10.1002/smll.202302302
Args:
inp (input_class): An instance containing input parameters.
Returns:
None: Saves the generated decahedral geometry.
"""
# Check input
inp.check_input_case()
general.create_results_geom()
#out_log = output.logfile_init()
# Extract lattice constant and merge cutoff
param = parameters.parameters()
lattice_constant = param.lattice_constant.get(inp.atomtype)
inp.merge_cutoff = param.min_dist.get(inp.atomtype)
# Initialize bulk "molecule" and read geometry
mol_sphere_1 = molecule.molecule()
mol_sphere_2 = molecule.molecule()
mol_cylinder = molecule.molecule()
mol_sphere_1.read_geom(inp.geom_file,False)
mol_sphere_2.read_geom(inp.geom_file,False)
mol_cylinder.read_geom(inp.geom_file,False)
# ----------------------- #
# --------- Out --------- #
# Populate input class
inp.atomtype = inp.atomtype_out
# Increase rod width and length to have a coating that
# is at least 3x(lattice constant) in the thinner part
minimum_coating = 3.0 * lattice_constant
if (abs(inp.rod_length - inp.bipyramid_length*2.0) < minimum_coating):
inp.rod_length = inp.bipyramid_length*2.0 + minimum_coating
if (abs(inp.rod_width - inp.bipyramid_width*2.0) < minimum_coating):
inp.rod_width = inp.bipyramid_width*2.0 + minimum_coating
# Create individual sphere at the extremes of the rods
tools.determine_sphere_center(inp,'+')
mol_sphere_1.filter_xyz_in_sphere(inp)
tools.determine_sphere_center(inp,'-')
mol_sphere_2.filter_xyz_in_sphere(inp)
# Create cylinder
mol_cylinder.filter_xyz_in_cylinder(inp)
# Merge cylinder and spheres geometries
mol_tmp = tools.merge_geoms(inp,mol_sphere_1,mol_sphere_2)
mol_out = tools.merge_geoms(inp,mol_cylinder,mol_tmp)
# ---------------------- #
# --------- In --------- #
# Copy mol_out as initial guess, change atomtype, and populate input class
mol_in = molecule.molecule()
mol_in.read_geom(inp.geom_file,False)
mol_in.change_atomtype(inp.atomtype_in)
inp.atomtype = inp.atomtype_in
# Define vertices with respect to the center (C[0,0,0]) of the XY-plane base
#
# XY base (center 6 in the +z direction)
#
#
# 1
# ,'.
# ,' `.
# ,' `.
# 5 ,' `. 2 ^ y
# \ C / |
# \ / o ---> x
# \ / z
# \_______/
# 4 3
#
# ------------------------------
# Define pentagonal-base pyramid
# -------------------------------
# Recover true circumradius
R = inp.bipyramid_width
# Calculate pentagon vertices (angles offset by 90° to have a flat side at bottom)
angles = [math.radians(90 + i * 72) for i in range(5)] # 72° between vertices
centers = {}
# Pentagon base vertices
for i in range(5):
centers[f"center_{i+1}"] = [
R * math.cos(angles[i]), # x coordinate
R * math.sin(angles[i]), # y coordinate
0.0 # z coordinate (base)
]
# Apex point at the top (half of total length from base to apex)
centers["center_6"] = [0.0, 0.0, inp.bipyramid_length ] # z coordinate (top apex). Added extra length to ensure lateral sides are filled by atoms
# Define planes to cut the pyramid
planes = {
"n_126": tools.calculate_normal_and_rhs(centers["center_1"], centers["center_2"], centers["center_6"]),
"n_236": tools.calculate_normal_and_rhs(centers["center_2"], centers["center_3"], centers["center_6"]),
"n_346": tools.calculate_normal_and_rhs(centers["center_3"], centers["center_4"], centers["center_6"]),
"n_456": tools.calculate_normal_and_rhs(centers["center_4"], centers["center_5"], centers["center_6"]),
"n_516": tools.calculate_normal_and_rhs(centers["center_5"], centers["center_1"], centers["center_6"]),
}
# Include atoms within the defined pentagonal pyramid
mol_in.filter_xyz_in_pentagonal_pyramid(inp,centers, planes)
# If the pyramid base is not on the XY plane, translate it
# To make it match with the outer rod, translate it as well
z_min = mol_in.xyz_min[-1]
if (z_min > 0.0 ):
mol_in.translate_geom(z_min,[0.0,0.0, -1.0])
mol_out.translate_geom(z_min,[0.0,0.0, -1.0])
# If halfshelf slice from z = 0 to z = -inf
if inp.pencil_type == "halfshell":
inp.atomtype = inp.atomtype_out
mol_out.slice_xyz_by_z_threshold(inp, z_threshold=-0.1)
# Rotate 180 degrees to create bottom pyramid and merge
mol_rot = molecule.molecule()
mol_rot = tools.rotate(mol_in,180.0,'+x',mol_rot)
mol_rot = tools.rotate(mol_rot,180.0,'+z',mol_rot)
mol_in = tools.merge_geoms(inp,mol_in,mol_rot)
# Remove dangling atom on extreme
mol_in.remove_dangling_atoms_metals(inp)
# Create shell by subtracting core geometry
mol_shell = tools.subtract_geoms(inp,mol_in,mol_out)
# Alloy core and shell
if inp.alloy:
inp.alloy_string = f"_alloy_{inp.alloy_perc}_perc"
inp.atomtype = inp.atomtype_out
inp.atomtype_alloy = inp.atomtype_in
mol_shell.create_alloy(inp)
inp.atomtype = inp.atomtype_in
inp.atomtype_alloy = inp.atomtype_out
mol_in.create_alloy(inp)
# Merge to create core-shell
mol_core_shell = tools.merge_geoms(inp,mol_in,mol_shell)
# Save filtered geometry
inp.bipyramid_width = inp.bipyramid_width * 2.0 # Make it match with initial definition
inp.bipyramid_length = inp.bipyramid_length * 2.0 # Make it match with initial definition
inp.xyz_output = f'pencil_{inp.pencil_type}_core_{inp.atomtype_in}_shell_{inp.atomtype_out}_in_width-{inp.bipyramid_width}_in_length-{inp.bipyramid_length}_out_width_{inp.rod_width}_out_length-{inp.rod_length}{inp.alloy_string}'
output.print_geom(mol_core_shell, inp.xyz_output)
# -------------------------------------------------------------------------------------
# -------------------------------------------------------------------------------------
[docs]
def get_layers(inp, lattice_constant):
"""
Dynamically calculates the required number of ASE layers for bulk structures
based on the structural shape.
Args:
inp (input_class): An instance containing input parameters.
lattice_constant (float): Lattice constant of the atomic structure.
Returns:
list[int]: Number of layers required along the x, y, and z axes.
"""
# Scaling factor to ensure enough layers are considered, independently of FCC of BCC structure
structure_scaling = 2.0
if inp.gen_sphere or inp.gen_sphere_core_shell:
# Find the maximum distance from the origin + radius
max_shift = max(abs(inp.sphere_center[0]), abs(inp.sphere_center[1]), abs(inp.sphere_center[2]))
R = inp.radius + max_shift
return [int(structure_scaling * R / lattice_constant) + 2] * 3 # Same layers for x, y, z
elif inp.gen_rod or inp.gen_rod_core_shell or inp.gen_pencil:
axis = inp.main_axis.lower()
R = inp.rod_width / 2.0
L = 2.0 * structure_scaling * inp.rod_length
layers = [max(3, int(structure_scaling * 2 * R / lattice_constant) + 2)] * 3
axis_index = {"x": 0, "y": 1, "z": 2}[axis]
layers[axis_index] = max(3, int(L / lattice_constant) + 2) # Adjust for rod length
return layers
elif inp.gen_bipyramid:
axis = "z"
R = inp.bipyramid_width / 2.0
L = 2.0 * structure_scaling * inp.bipyramid_length
layers = [max(3, int(structure_scaling * 2 * R / lattice_constant) + 2)] * 3
axis_index = {"x": 0, "y": 1, "z": 2}[axis]
layers[axis_index] = max(3, int(L / lattice_constant) + 2) # Adjust for rod length
return layers
elif inp.gen_pentpyramid:
R = inp.base_width / 2.0
L = structure_scaling * inp.z_max
layers = [max(3, int(structure_scaling * 2 * R / lattice_constant) + 2)] * 3
layers[2] = max(3, int(L / lattice_constant) + 2)
return layers
elif inp.gen_pentbipyramid:
R = inp.rod_width / 2.0
L = 2.0 * structure_scaling * inp.rod_length
layers = [max(3, int(structure_scaling * 2 * R / lattice_constant) + 2)] * 3
layers[2] = max(3, int(L / lattice_constant) + 2) # z axis
return layers
elif inp.gen_tip:
H = structure_scaling * inp.z_max * 1.8 # Tip height
return [int(H / lattice_constant) + 2] * 3 # Tip grows mostly in z-axis
elif inp.gen_cone:
H = inp.z_max
R = structure_scaling * inp.radius
return [int(2 * R / lattice_constant) + 2, int(2 * R / lattice_constant) + 2, int(H / lattice_constant) + 2]
elif inp.gen_pyramid:
if inp.z_max > inp.side_length:
H = 1.5 * structure_scaling * inp.z_max
L = inp.side_length
elif inp.side_length > inp.z_max:
H = inp.z_max
L = 1.5 * structure_scaling * inp.side_length
else:
# If both are equal, apply scaling to L
H = inp.z_max
L = 1.5 * structure_scaling * inp.side_length
return [int(L / lattice_constant) + 2, int(L / lattice_constant) + 2, int(H / lattice_constant) + 2]
elif inp.gen_microscope:
H_paraboloid = 1.5 * structure_scaling * inp.z_max_paraboloid
H_pyramid = 1.5 * structure_scaling * inp.z_max_pyramid
L = 1.5 * structure_scaling * inp.side_length
return [int(L / lattice_constant) + 2, int(L / lattice_constant) + 2, int((H_paraboloid + H_pyramid) / lattice_constant) + 2]
# -------------------------------------------------------------------------------------
[docs]
def create_ase_bulk_graphene(inp, base_dir):
"""
Creates a temporary bulk graphene XYZ file using ASE.
Args:
inp (input_class): An instance containing input parameters.
base_dir (str): Absolute path to the working directory.
Returns:
None: Saves the generated bulk graphene geometry.
"""
# Extract lattice constant from parameters dictionary
param = parameters.parameters()
lattice_constant = param.lattice_constant.get(inp.atomtype)
# Create tmp folder
inp.tmp_folder = os.path.join(base_dir,'tmp')
if os.path.exists(inp.tmp_folder): shutil.rmtree(inp.tmp_folder)
os.mkdir(inp.tmp_folder)
# Create initial graphene structure with ASE
if inp.graphene_structure == "rib" or inp.graphene_structure == "triangle":
# Increase size of bulk structure by a factor to ensure correct geometry creation
geom_scale = 1.5
# If triangle, create side_length X side_length ribbon to cut
if inp.graphene_structure == "triangle":
inp.X_length = inp.side_length
inp.Y_length = inp.side_length
geom_scale = 2.0
# Oversize by 1.5x
scaled_width = geom_scale * inp.X_length
scaled_length = geom_scale * inp.Y_length
# Convert dimensions to graphene_nanoribbon parameters
n = int(scaled_width / 2.13) # Number of dimer rows for width
m = int(scaled_length / 4.26) # Number of unit cells for length
# Create flat armchair nanoribbon in XY plane
# Armchair structure will then be managed
graphene_xyz = graphene_nanoribbon(n=n, m=m, type='armchair')
graphene_xyz.rotate(90, 'x', rotate_cell=True)
elif inp.graphene_structure == "disk" or inp.graphene_structure == "ring":
# Select radius to define structure. Both inp.radius and inp.radius_out are initialized to zero
radius_selected = max(inp.radius, inp.radius_out)
# Create a large graphene sheet in XY plane
radius_int = math.ceil(radius_selected) # Round to upper integer for function compatibility
graphene_xyz = graphene_general_ase(a=lattice_constant, size=(2*radius_int, 2*radius_int, 1)) # Large enough to extract disk
# Center the graphene at (0,0,0) and write on tmp/tmp_bulk.xyz
graphene_xyz.translate(-graphene_xyz.get_center_of_mass())
inp.geom_file = os.path.join(inp.tmp_folder,'tmp_bulk.xyz')
write(inp.geom_file, graphene_xyz)
# -------------------------------------------------------------------------------------