diff --git a/README.md b/README.md index feefdc87..9d33f7e9 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,13 @@ conda env create --file environment.yml ## Installation -Currently the best way to use pyCOFBuilder is to manually import it using the `sys` module, as exemplified below: +You can install pyCOFBuilder using pip: + +```Shell +pip install pycofbuilder +``` + +Alternativelly, you can use pyCOFBuilder by manually import it using the `sys` module, as exemplified below: ```python # importing module diff --git a/build/lib/__init__.py b/build/lib/__init__.py new file mode 100644 index 00000000..578f44f8 --- /dev/null +++ b/build/lib/__init__.py @@ -0,0 +1,48 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +import os +import sys + +# This code is written for Python 3 +if sys.version_info[0] != 3: + raise Exception("Sorry but pyCOFBuilder requires Python 3.") + sys.exit(1) + +# Import BuildingBlocks class +from pycofbuilder.building_block import BuildingBlock + +# Import Framework class +from pycofbuilder.framework import Framework + +# Import Tools +import pycofbuilder.tools as Tools +import pycofbuilder.io_tools as IO_Tools + +# Import ChemJSON +import pycofbuilder.chemjson as ChemJSON + +# Import Exceptions +import pycofbuilder.exceptions as Exceptions + +# Import Logger +import pycofbuilder.logger as Logger + +__all__ = [ + 'BuildingBlock', + 'Framework', + 'Tools', + 'IO_Tools', + 'ChemJSON', + 'Exceptions', + 'Logger' + ] + +_ROOT = os.path.abspath(os.path.dirname(__file__)) + +__author__ = "Felipe Lopes de Oliveira" +__license__ = "MIT" +__version__ = '0.0.6' +__email__ = "felipe.lopes@nano.ufrj.br" +__status__ = "Development" diff --git a/build/lib/building_block.py b/build/lib/building_block.py new file mode 100644 index 00000000..0b9c3b74 --- /dev/null +++ b/build/lib/building_block.py @@ -0,0 +1,686 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +""" +The BuildingBlock class is used to create the building blocks for the Framework class. +""" + +import os +import copy +import numpy as np +from scipy.spatial.transform import Rotation as R + +from pycofbuilder.tools import (rotation_matrix_from_vectors, + closest_atom, + closest_atom_struc, + find_index, + unit_vector) +from pycofbuilder.io_tools import (read_xyz, save_xyz, read_gjf) +from pycofbuilder.cjson import ChemJSON + +from pycofbuilder.logger import create_logger + +from pycofbuilder.exceptions import MissingXError + + +class BuildingBlock(): + + def __init__(self, name: str = '', **kwargs): + + self.name: str = name + + self.out_path: str = kwargs.get('out_dir', os.path.join(os.getcwd(), 'out')) + self.save_bb: bool = kwargs.get('save_bb', True) + self.bb_out_path: str = kwargs.get('bb_out_path', os.path.join(self.out_path, 'building_blocks')) + + self.logger = create_logger(level=kwargs.get('log_level', 'info'), + format=kwargs.get('log_format', 'simple'), + save_to_file=kwargs.get('save_to_file', False), + log_filename=kwargs.get('log_filename', 'pycofbuilder.log')) + + _ROOTDIR = os.path.abspath(os.path.dirname(__file__)) + self.main_path = os.path.join(_ROOTDIR, 'data') + + self.connectivity = kwargs.get('connectivity', None) + self.size = kwargs.get('size', None) + self.mass = kwargs.get('mass', None) + self.composition = kwargs.get('composition', None) + + self.atom_types = kwargs.get('atom_types', None) + self.atom_pos = kwargs.get('atom_pos', None) + self.atom_labels = kwargs.get('atom_labels', None) + + self.smiles = kwargs.get('smiles', None) + self.charge = kwargs.get('charge', 0) + self.multiplicity = kwargs.get('multiplicity', 1) + self.chirality = kwargs.get('chirality', 0) + self.symmetry = kwargs.get('symmetry', None) + + self.core = kwargs.get('core', None) + self.conector = kwargs.get('conector', None) + self.funcGroups = kwargs.get('funcGroups', None) + + self.available_symmetry = ['L2', + 'T3', + 'S4', 'D4', # 'R4' + 'H6', # 'O6', 'P6' + # 'C8', 'A8', 'E8' + # 'B12', 'I12', 'U12', 'X12' + ] + + # Check if bb_out_path exists and try to create it if not + if self.save_bb != '': + os.makedirs(self.bb_out_path, exist_ok=True) + + # If a name is provided create building block from this name + if self.name != '': + if '.' not in self.name: + self.from_name(self.name) + + # If a name is provided and it is a file, create building block from this file + if '.' in self.name: + self.from_file(self.out_path, self.name) + + def __str__(self): + return self.structure_as_string() + + def __repr__(self): + return 'BuildingBlock({}, {}, {}, {})'.format(self.symmetry, + self.core, + self.conector, + self.funcGroups) + + def copy(self): + '''Return a deep copy of the BuildingBlock object''' + return copy.deepcopy(self) + + def from_file(self, path, file_name): + '''Read a building block from a file''' + extension = file_name.split('.')[-1] + + read_func_dict = {'xyz': read_xyz, + 'gjf': read_gjf} + + self.name = file_name.rstrip(f'.{extension}') + self.atom_types, self.atom_pos = read_func_dict[extension](path, file_name) + self.atom_labels = ['C']*len(self.atom_types) + + if any([i == 'X' for i in self.atom_types]): + self.connectivity = len([i for i in self.atom_types if 'X' in i]) + else: + raise MissingXError() + + self.centralize(by_X=True) + self.calculate_size() + pref_orientation = unit_vector( + self.get_X_points()[1][0]) + + self.align_to(pref_orientation) + + def from_name(self, name): + '''Automatically read or create a buiding block based on its name''' + + # Check the existence of the building block files + symm_check, core_check, conector_check, funcGroup_check = self.check_existence(name) + + error_msg = "Building Block name is invalid!\n" + error_msg += "Symm: {}, Core: {}, Connector: {}, Functional Group:{}".format(symm_check, + core_check, + conector_check, + funcGroup_check) + assert all([symm_check, core_check, conector_check, funcGroup_check]), error_msg + + self.name = name + + BB_name = self.name.split('_') + self.symmetry = BB_name[0] + self.core = BB_name[1] + self.conector = BB_name[2] + possible_funcGroups = BB_name[3:] + ['H'] * (9 - len(BB_name[3:])) + + self.create_BB_structure(self.symmetry, + self.core, + self.conector, + *possible_funcGroups) + if self.save_bb: + self.save() + + def n_atoms(self): + ''' Returns the number of atoms in the unitary cell''' + return len(self.atom_types) + + def print_structure(self): + """ + Print the structure in xyz format: + `atom_label pos_x pos_y pos_z` + """ + + print(self.structure_as_string()) + + def centralize(self, by_X=True): + ''' Centralize the molecule on its geometrical center''' + + transposed = np.transpose(self.atom_pos) + if by_X is True: + x_transposed = np.transpose(self.get_X_points()[1]) + if by_X is False: + x_transposed = transposed + cm_x = transposed[0] - np.average(x_transposed[0]) + cm_y = transposed[1] - np.average(x_transposed[1]) + cm_z = transposed[2] - np.average(x_transposed[2]) + + self.atom_pos = np.transpose([cm_x, cm_y, cm_z]) + return np.transpose([cm_x, cm_y, cm_z]) + + def get_X_points(self): + '''Get the X points in a molecule''' + + if 'X' in self.atom_types: + X_labels, X_pos = [], [] + for i in range(len(self.atom_types)): + if self.atom_types[i] == 'X': + X_labels += [self.atom_types[i]] + X_pos += [self.atom_pos[i]] + + return X_labels, np.array(X_pos) + + else: + print('No X ponts could be found!') + return self.atom_types, self.atom_pos + + def get_Q_points(self, atom_types, atom_pos): + '''Get the Q points in a molecule''' + + Q_labels, Q_pos = [], [] + + for i in range(len(atom_types)): + if atom_types[i] == 'Q': + Q_labels += [atom_types[i]] + Q_pos += [atom_pos[i]] + + return Q_labels, np.array(Q_pos) + + def get_R_points(self, atom_types, atom_pos): + """ + Get the R points in a molecule + """ + + # Create a dict with the R points + R_dict = {'R': [], + 'R1': [], + 'R2': [], + 'R3': [], + 'R4': [], + 'R5': [], + 'R6': [], + 'R7': [], + 'R8': [], + 'R9': []} + + # Iterate over the R keys + for key in R_dict.keys(): + # Iterate over the atoms + for i in range(len(atom_types)): + if atom_types[i] == key: + R_dict[key] += [atom_pos[i]] + + return R_dict + + def calculate_size(self): + '''Calculate the size of the building block''' + _, X_pos = self.get_X_points() + self.size = [np.linalg.norm(i) for i in X_pos] + + def align_to(self, vec: list = [0, 1, 0], n: int = 0): + ''' + Align the first n-th X point to a given vector + + Parameters + ---------- + vec : list + Vector to align the molecule + n : int + Index of the X point to be aligned + align_to_y : bool + If True, the second point is aligned to the y axis + ''' + _, X_pos = self.get_X_points() + R_matrix = rotation_matrix_from_vectors(X_pos[n], vec) + + self.atom_pos = np.dot(self.atom_pos, np.transpose(R_matrix)) + + def rotate_around(self, rotation_axis: list = [1, 0, 0], angle: float = 0.0, degree: bool = True): + ''' + Rotate the molecule around a given axis + + Parameters + ---------- + rotation_axis : list + Rotation axis + angle : float + Rotation angle in degrees + degree : bool + If True, the angle is given in degrees. + Default is True. + ''' + + rotation_axis = unit_vector(rotation_axis) + rotation = R.from_rotvec(angle * rotation_axis, degrees=degree) + + self.atom_pos = rotation.apply(self.atom_pos) + + def rotate_to_xy_plane(self): + '''Rotate the molecule to the xy plane''' + _, X_pos = self.get_X_points() + + if len(X_pos) == 3: + + normal = np.cross(X_pos[0], X_pos[-1]) + if normal[0] != 0 and normal[1] != 0: + R_matrix = rotation_matrix_from_vectors(normal, [0, 0, 1]) + self.atom_pos = np.dot(self.atom_pos, np.transpose(R_matrix)) + + if len(X_pos) == 2: + normal = np.cross(X_pos[0], self.atom_pos[1]) + if normal[0] != 0 and normal[1] != 0: + R_matrix = rotation_matrix_from_vectors(normal, [0, 0, 1]) + self.atom_pos = np.dot(self.atom_pos, np.transpose(R_matrix)) + + def shift(self, shift_vector: list): + ''' + Shift the molecule by a given vector + + Parameters + ---------- + shift_vector : list + Shift vector + ''' + + self.atom_pos = np.array(self.atom_pos) + np.array(shift_vector) + + def structure_as_string(self): + struct_string = '' + for i, _ in enumerate(self.atom_types): + struct_string += '{:<5s}{:>10.7f}{:>15.7f}{:>15.7f}\n'.format(self.atom_types[i], + self.atom_pos[i][0], + self.atom_pos[i][1], + self.atom_pos[i][2]) + + return struct_string + + def add_connection_group(self, conector_name): + '''Adds the functional group by which the COF will be formed from the building blocks''' + + connector = ChemJSON() + connector.from_cjson(os.path.join(self.main_path, 'conector'), conector_name) + + self.smiles = self.smiles.replace('[Q]', + f"{connector.properties['smiles'].replace('[Q]', '')}") + + conector_label = connector.atomic_types + conector_pos = connector.cartesian_positions + + # Get the position of the Q points in the structure + location_Q_struct = self.get_Q_points(self.atom_types, self.atom_pos) + + for i in range(len(location_Q_struct[0])): + + n_conector_label = conector_label.copy() + n_conector_pos = conector_pos.copy() + + # Get the position of the closest atom to Q in the structure + close_Q_struct = closest_atom('Q', + location_Q_struct[1][i], + self.atom_types, + self.atom_pos)[1] + + # Get the position of Q in the conection group + location_Q_connector = self.get_Q_points(n_conector_label, n_conector_pos) + + # Get the position of the closest atom to Q in the conection group + close_Q_connector = closest_atom('Q', + location_Q_connector[1][0], + n_conector_label, + n_conector_pos)[1] + + # Create the vector Q in the structure + v1 = close_Q_struct - location_Q_struct[1][i] + # Create the vector Q in the conector + v2 = np.array(close_Q_connector) - np.array(location_Q_connector[1][0]) + + # Find the rotation matrix that align v2 with v1 + Rot_m = rotation_matrix_from_vectors(v2, v1) + + # Delete the "Q" atom position of the conector group and the structure + n_conector_pos = np.delete( + n_conector_pos, + find_index(np.array([0., 0., 0.]), n_conector_pos), + axis=0) + + self.atom_labels = np.delete( + self.atom_labels, + find_index(location_Q_struct[1][i], self.atom_pos), + axis=0 + ) + + self.atom_pos = np.delete( + self.atom_pos, + find_index(location_Q_struct[1][i], self.atom_pos), + axis=0 + ) + + # Rotate and translade the conector group to Q position in the strucutre + rotated_translated_group = np.dot(n_conector_pos, -np.transpose(Rot_m)) + location_Q_struct[1][i] + + # Add the position of conector atoms to the main structure + self.atom_pos = np.append(self.atom_pos, rotated_translated_group, axis=0) + + # Remove the Q atoms from structure + self.atom_types.remove('Q') + n_conector_label.remove('Q') + + self.atom_types = self.atom_types + n_conector_label + + self.atom_labels = np.append(self.atom_labels, [['Q'] * len(n_conector_label)]) + + def add_connection_group_symm(self, conector_name): + '''Adds the functional group by which the COF will be formed from the building blocks''' + + connector = ChemJSON() + connector.from_cjson(os.path.join(self.main_path, 'conector'), conector_name) + + self.smiles = self.smiles.replace('[Q]', + f"{connector.properties['smiles'].replace('[Q]', '')}") + + conector_types = connector.atomic_types + conector_pos = connector.cartesian_positions + + # Get the position of the Q points in the structure + _, Q_vec = self.get_Q_points(self.atom_types, self.atom_pos) + + # Remove the Q atoms from structure + self.atom_types = self.atom_types[:-4] + self.atom_pos = self.atom_pos[:-4] + self.atom_labels = self.atom_labels[:-4] + + # Create the vector Q in the structure + QS_vector = Q_vec[0] + + # Create the vector Q in the conector + QC_vector = np.array(conector_pos[1]) + + # Find the rotation matrix that align the connector with the structure + Rot_m = rotation_matrix_from_vectors(QC_vector, QS_vector) + + # Rotate and translade the conector group to Q position in the strucutre + conector_pos = np.dot(conector_pos, np.transpose(Rot_m)) + Q_vec[0] + conector_pos = R.from_rotvec( + -90 * unit_vector(conector_pos[1] - conector_pos[0]), degrees=True).apply(conector_pos) + + # Add the position of conector atoms to the main structure + self.atom_types = list(self.atom_types) + conector_types[1:] + self.atom_pos = list(self.atom_pos) + list(conector_pos[1:]) + self.atom_labels = np.append(self.atom_labels, [['Q'] * len(conector_types[1:])]) + + # Apply the improper rotations to match S4 symmetry + Q_vec = [unit_vector(i) for i in Q_vec] + + # First S4 axis is location_Q_struct[1] + R1 = R.from_rotvec(120 * Q_vec[1], degrees=True).apply(conector_pos) + R1 = R.from_rotvec(60 * unit_vector(R1[1] - R1[0]), degrees=True).apply(R1) + + # Add the position of conector atoms to the main structure + self.atom_types = list(self.atom_types) + conector_types[1:] + self.atom_pos = list(self.atom_pos) + list(R1[1:]) + self.atom_labels = np.append(self.atom_labels, [['Q'] * len(conector_types[1:])]) + + # Second S4 axis is location_Q_struct[2] + R2 = R.from_rotvec(120 * Q_vec[2], degrees=True).apply(conector_pos) + R2 = R.from_rotvec(-120 * unit_vector(R2[1]-R2[0]), degrees=True).apply(R2) + + # Add the position of conector atoms to the main structure + self.atom_types = list(self.atom_types) + conector_types[1:] + self.atom_pos = list(self.atom_pos) + list(R2)[1:] + self.atom_labels = np.append(self.atom_labels, [['Q'] * len(conector_types[1:])]) + + # Third S4 axis is location_Q_struct[3] + R3 = R.from_rotvec(120 * Q_vec[3], degrees=True).apply(conector_pos) + R3 = R.from_rotvec(60 * unit_vector(R3[1] - R3[0]), degrees=True).apply(R3) + + # Add the position of conector atoms to the main structure + self.atom_types = list(self.atom_types) + conector_types[1:] + self.atom_pos = list(self.atom_pos) + list(R3[1:]) + self.atom_labels = np.append(self.atom_labels, [['Q'] * len(conector_types[1:])]) + + def add_R_group(self, R_name, R_type): + '''Adds group R in building blocks''' + + rgroup = ChemJSON() + rgroup.from_cjson(os.path.join(self.main_path, 'func_groups'), R_name) + + self.smiles = self.smiles.replace(f'[{R_type}]', + f"{rgroup.properties['smiles'].replace('[Q]', '')}") + + group_label = rgroup.atomic_types + group_pos = rgroup.cartesian_positions + + # Get the position of the R points in the structure + location_R_struct = self.get_R_points(self.atom_types, self.atom_pos)[R_type] + + # Get the position of the R points in the R group + for i, _ in enumerate(location_R_struct): + n_group_label = group_label.copy() + n_group_pos = group_pos.copy() + + # Get the position of the closest atom to R in the structure + close_R_struct = closest_atom_struc(R_type, + location_R_struct[i], + self.atom_types, + self.atom_pos)[1] + + # Get the position of R in the R group + pos_R_group = self.get_R_points(n_group_label, n_group_pos)['R'] + + # Get the position of the closest atom to R in the R group + close_R_group = closest_atom('R', pos_R_group[0], n_group_label, n_group_pos)[1] + + # Create the vector R in the structure + v1 = close_R_struct - location_R_struct[i] + + # Create the vector R in the R group + v2 = np.array(close_R_group) - np.array(pos_R_group[0]) + + # Find the rotation matrix that align v2 with v1 + Rot_m = rotation_matrix_from_vectors(v2, v1) + + # Delete the "R" atom position of the R group and the structure + n_group_pos = np.delete( + n_group_pos, + find_index(np.array([0.0, 0.0, 0.0]), n_group_pos), + axis=0 + ) + + # Rotate and translade the R group to R position in the strucutre + rotated_translated_group = np.dot( + n_group_pos, + -np.transpose(Rot_m) + ) + location_R_struct[i] + + # Remove the R atoms from the labels list + # Remove the R atoms from structure + self.atom_labels = np.delete( + self.atom_labels, + find_index(location_R_struct[i], self.atom_pos), + axis=0 + ) + + # Remove the R atoms from structure + self.atom_pos = np.delete( + self.atom_pos, + find_index(location_R_struct[i], self.atom_pos), + axis=0 + ) + + # Add the position of rotated atoms to the main structure + self.atom_pos = np.append(self.atom_pos, rotated_translated_group, axis=0) + + # Remove the R atoms from structure + self.atom_types.remove(R_type) + + # Remove the R atoms from R group + n_group_label.remove('R') + + self.atom_types = self.atom_types + n_group_label + self.atom_labels = np.append(self.atom_labels, ['R'] * len(n_group_label)) + + def create_BB_structure(self, + symmetry='L2', + core_name='BENZ', + conector='CHO', + R1='H', + R2='H', + R3='H', + R4='H', + R5='H', + R6='H', + R7='H', + R8='H', + R9='H'): + '''Create a building block''' + + self.name = f'{symmetry}_{core_name}_{conector}' + + core = ChemJSON() + core.from_cjson(os.path.join(self.main_path, 'core', symmetry), core_name) + + self.smiles = core.properties['smiles'] + + self.atom_types = core.atomic_types + self.atom_pos = core.cartesian_positions + self.composition = core.formula + self.atom_labels = ['C']*len(self.atom_types) + + pref_orientation = unit_vector( + self.get_Q_points(core.atomic_types, core.cartesian_positions)[1][0]) + + if symmetry == 'D4': + self.add_connection_group_symm(conector) + else: + self.add_connection_group(conector) + + R_list_names = [R1, R2, R3, R4, R5, R6, R7, R8, R9] + R_list_labels = ['R1', 'R2', 'R3', 'R4', 'R5', 'R6', 'R7', 'R8', 'R9'] + + funcGroup_string = [] + for i in range(len(R_list_names)): + if R_list_labels[i] in self.atom_types: + self.add_R_group(R_list_names[i], R_list_labels[i]) + self.name += f'_{R_list_names[i]}' + funcGroup_string.append(R_list_names[i]) + + self.funcGroups = funcGroup_string + + self.connectivity = len([i for i in self.atom_types if 'X' in i]) + self.centralize() + self.align_to(pref_orientation) + self.calculate_size() + + def replace_X(self, target_type): + for i in range(len(self.atom_types)): + if self.atom_types[i] == "X": + self.atom_types[i] = target_type + + def remove_X(self): + + atom_types, atom_pos, atom_labels = [], [], [] + for i in range(len(self.atom_types)): + if self.atom_types[i] != "X": + atom_types.append(self.atom_types[i]) + atom_pos.append(self.atom_pos[i]) + atom_labels.append(self.atom_labels[i]) + + self.atom_types = atom_types + self.atom_pos = np.array(atom_pos) + self.atom_labels = atom_labels + + def save(self, extension='xyz'): + + if extension == 'xyz': + save_xyz(path=self.bb_out_path, + file_name=self.name + '.xyz', + atom_types=self.atom_types, + atom_pos=self.atom_pos) + + def get_available_core(self): + '''Get the list of available cores''' + + available_cores = {} + + for symm in self.available_symmetry: + symm_path = os.path.join(self.main_path, 'core', symm) + available_cores[symm] = [i.rstrip('.cjson') for i in os.listdir(symm_path) if '.cjson' in i] + + return available_cores + + def get_available_R(self): + '''Get the list of available functional groups''' + R_PATH = os.path.join(self.main_path, 'func_groups') + R_list = [i.rstrip('.cjson') for i in os.listdir(R_PATH) if '.cjson' in i] + + return R_list + + def get_available_conector(self): + '''Get the list of available conectores''' + C_PATH = os.path.join(self.main_path, 'conector') + C_list = [i.rstrip('.cjson') for i in os.listdir(C_PATH) if '.cjson' in i] + + return C_list + + def check_existence(self, name): + + symm_check = False + core_check = False + conector_check = False + funcGroup_check = True + + name = name.split('_') + symm = name[0] + core = name[1] + conector = name[2] + funcGroups = name[3:] + + BB_dict = self.get_available_core() + + if symm in self.available_symmetry: + symm_check = True + else: + print('ERROR!: Building Block symmetry must be L2, T3, S4, or H6.') + symm_check = False + + if core in BB_dict[symm]: + core_check = True + else: + print(f'ERROR!: {core} not available!') + print(f'Available cores with {symm} symmetry are {BB_dict[symm]}') + + if conector in self.get_available_conector(): + conector_check = True + else: + print(f'ERROR! {conector} is not a available conector.') + print(f'Available list: {self.get_available_conector()}') + + possible_funcGroups_list = self.get_available_R() + for func in funcGroups: + if func not in possible_funcGroups_list: + print(f'ERROR! Functional group {func} is not a available.') + print(f'Available list: {possible_funcGroups_list}') + funcGroup_check = False + + return symm_check, core_check, conector_check, funcGroup_check + + def get_buildingblock_list(self, shape, connector_group): + + files_list = os.listdir(self.bb_out_path) + + return [i.rstrip('.xyz') for i in files_list if shape == i.split('_')[0] and connector_group in i.split('_')[2]] diff --git a/build/lib/cjson.py b/build/lib/cjson.py new file mode 100644 index 00000000..f980126d --- /dev/null +++ b/build/lib/cjson.py @@ -0,0 +1,403 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +""" +The CJSON package implements functions to read, create and manipulate Chemical JSON objects. +""" + +import os +import simplejson +import numpy as np + +from pycofbuilder.tools import elements_dict + +import gemmi +from ase.cell import Cell + + +class ChemJSON: + ''' + Class to read, create and manupulate ChemJSON files. + + Attributes + ---------- + + file_name : str + The name of the file. + name : str + The name of the structure. + cell_parameters : list + The cell parameters of the structure as a (1,6) list. + cell_matrix : list + The cell matrix of the structure as a (3,3) list. + cartesian_positions : list + The cartesian positions of the structure as a (n,3) list. + fractional_positions : list + The fractional positions of the structure as a (n,3) list. + atomic_numbers : list + The atomic numbers of the structure as a (n,1) list. + atomic_types : list + The atomic types of the structure as a (n,1) list. + atomic_labels : list + The atomic labels of the structure as a (n,1) list. + formula : str + The formula of the structure. + properties : dict + The properties of the structure. + partial_charges : dict + A dictionary contaning the partial charges of the atoms on the structure. + Example: {'DDEC': [0.1, 0.2, 0.15], 'EQeq': [0.05, 0.15, 0.19]} + ''' + def __init__(self): + self.file_name = '' + self.name = '' + + # Structure properties + self.cell_parameters = None + self.cell_matrix = None + self.cartesian_positions = None + self.fractional_positions = None + self.atomic_numbers = None + self.atomic_types = None + self.atomic_labels = None + self.formula = '' + self.partial_charges = None + + self.properties = None + self.results = [] + + # Create a custom representation of the class + def __repr__(self): + ''' + Returns a custom representation of the class. + ''' + + repr_string = "ChemJSON(name='{}', formula='{}', number of atoms={}".format(self.name, + self.formula, + len(self.atomic_types)) + + return repr_string + + # Create a custom print of the class + def __str__(self): + ''' + Returns a custom print of the class. + ''' + + string_string = "ChemJSON(name='{}', formula='{}', number of atoms={})\n".format(self.name, + self.formula, + len(self.atomic_types)) + + if self.cell_parameters is not None: + string_string += f"""Cell parameters: + a = {self.cell_parameters[0]:>12.7f} Å + b = {self.cell_parameters[1]:>12.7f} Å + c = {self.cell_parameters[2]:>12.7f} Å + α = {self.cell_parameters[3]:>12.7f} ° + β = {self.cell_parameters[4]:>12.7f} ° + γ = {self.cell_parameters[5]:>12.7f} ° + +Cell matrix: +A {self.cell_matrix[0][0]:>12.7f} {self.cell_matrix[0][1]:>12.7f} {self.cell_matrix[0][2]:>12.7f} +B {self.cell_matrix[1][0]:>12.7f} {self.cell_matrix[1][1]:>12.7f} {self.cell_matrix[1][2]:>12.7f} +C {self.cell_matrix[2][0]:>12.7f} {self.cell_matrix[2][1]:>12.7f} {self.cell_matrix[2][2]:>12.7f} +""" + if self.cartesian_positions is not None: + string_string += "Cartesian positions:\n" + for i, position in enumerate(self.cartesian_positions): + string_string += " {:3} {:>9.5f} {:>9.5f} {:>9.5f}\n".format(self.atomic_types[i], + position[0], + position[1], + position[2] + ) + + if self.fractional_positions is not None: + string_string += "Fractional positions:\n" + for i, position in enumerate(self.fractional_positions): + string_string += " {:3} {:>9.5f} {:>9.5f} {:>9.5f}\n".format(self.atomic_types[i], + position[0], + position[1], + position[2]) + + return string_string + + def set_properties(self, properties): + ''' + Sets the properties of the structure. + ''' + self.properties = properties + + def set_results(self, results): + ''' + Sets the results of the structure. + ''' + self.results = results + + def set_cell_parameters(self, cell_parameters): + ''' + Sets the cell parameters of the structure. + ''' + self.cell_parameters = cell_parameters + + aseCell = Cell.fromcellpar(cell_parameters) + self.cell_matrix = np.array(aseCell) + + def set_cell_matrix(self, cell_matrix): + ''' + Sets the cell matrix of the structure. The cell + parameters will be calculated and also updated. + ''' + self.cell_matrix = cell_matrix + + aseCell = Cell(cell_matrix) + self.cell_parameters = aseCell.cellpar() + + def set_cartesian_positions(self, cartesian_positions): + ''' + Sets the cartesian positions of the structure. The fractional + positions will be calculated and also updated. + ''' + self.cartesian_positions = np.array(cartesian_positions).astype(float) + + if self.cell_parameters is not None: + aseCell = Cell.fromcellpar(self.cell_parameters) + self.fractional_positions = aseCell.scaled_positions(cartesian_positions) + + def set_fractional_positions(self, fractional_positions): + ''' + Sets the fractional positions of the structure. The cartesian + positions will be calculated and also updated. + ''' + self.fractional_positions = np.array(fractional_positions).astype(float) + + if self.cell_parameters is not None: + aseCell = Cell.fromcellpar(self.cell_parameters) + self.cartesian_positions = aseCell.cartesian_positions(fractional_positions) + + def set_atomic_types(self, atomic_types): + ''' + Sets the atomic labels of the structure. + ''' + self.atomic_types = atomic_types + + self.atomic_labels = [f"{atom}{i+1}" for i, atom in enumerate(self.atomic_types)] + + symbol_dict = elements_dict('atomic_number') + self.atomic_numbers = [symbol_dict[i] for i in atomic_types] + + self.formula = ''.join([f'{atom}{self.atomic_types.count(atom)}' for atom in set(self.atomic_types)]) + + def set_atomic_numbers(self, atomic_numbers): + ''' + Sets the atomic numbers of the structure. The atomic types and formula + will be calculated and also updated. + ''' + self.atomic_numbers = atomic_numbers + + symbol_dict = elements_dict('atomic_number') + number_dict = {j: i for i, j in zip(symbol_dict.keys(), symbol_dict.values())} + + self.atomic_types = [number_dict[i] for i in atomic_numbers] + + self.atomic_labels = [f"{atom}{i+1}" for i, atom in enumerate(self.atomic_types)] + + self.formula = ''.join([f'{atom}{self.atomic_types.count(atom)}' for atom in set(self.atomic_types)]) + + def from_cjson(self, path, file_name): + ''' + Reads a ChemJSON file from a given path and file_name. + ''' + self.file_name = os.path.join(path, file_name.split('.')[0] + '.cjson') + + with open(self.file_name, 'r') as file: + cjson_data = simplejson.load(file) + + if "name" in cjson_data: + self.name = cjson_data['name'] + + if "unitCell" in cjson_data: + if 'a' in cjson_data['unitCell']: + self.set_cell_parameters( + [cjson_data['unitCell'][i] for i in ['a', 'b', 'c', 'alpha', 'beta', 'gamma']] + ) + elif 'cellVectors' in cjson_data['unitCell']: + self.set_cell_matrix( + np.array(cjson_data['unitCell']['cellVectors']).reshape(3, 3) + ) + + if "atoms" in cjson_data: + if 'coords' in cjson_data['atoms']: + if '3d' in cjson_data['atoms']['coords']: + self.set_cartesian_positions( + np.array(cjson_data['atoms']['coords']['3d']).reshape(-1, 3) + ) + elif '3dFractional' in cjson_data['atoms']['coords']: + self.set_fractional_positions( + np.array(cjson_data['atoms']['coords']['3dFractional']).reshape(-1, 3) + ) + + if "elements" in cjson_data['atoms']: + if 'type' in cjson_data['atoms']['elements']: + self.set_atomic_types(cjson_data['atoms']['elements']['type']) + + elif 'number' in cjson_data['atoms']['elements']: + self.set_atomic_numbers(cjson_data['atoms']['elements']['number']) + + if 'label' in cjson_data['atoms']['elements']: + self.atomic_labels = cjson_data['atoms']['elements']['label'] + else: + self.atomic_labels = [f"{atom}{i+1}" for i, atom in enumerate(self.atomic_types)] + + if 'properties' in cjson_data: + self.set_properties(cjson_data['properties']) + if 'results' in cjson_data: + self.set_results(cjson_data['results']) + if 'partialCharges' in cjson_data: + self.partial_charges = cjson_data['partialCharges'] + + def from_xyz(self, path, file_name): + ''' + Reads a XYZ file from a given path and file_name. + ''' + + self.file_name = os.path.join(path, file_name.split('.')[0] + '.xyz') + self.name = file_name.split('.')[0] + + with open(self.file_name, 'r') as file: + xyz_data = file.read().splitlines() + n_atoms = int(xyz_data[0]) + + atomic_types = [] + cartesian_positions = [] + + for line in xyz_data[2: n_atoms + 3]: + atomic_types.append(line.split()[0]) + cartesian_positions.append([float(i) for i in line.split()[1:]]) + + self.set_atomic_types(atomic_types) + + self.set_cartesian_positions(np.array(cartesian_positions)) + + def from_gjf(self, path, file_name): + ''' + Reads a Gaussian input file from a given path and file_name. + ''' + + self.file_name = os.path.join(path, file_name.split('.')[0] + '.gjf') + self.name = file_name.split('.')[0] + + with open(self.file_name, 'r') as file: + gjf_data = file.read().splitlines() + + # Remove empty lines + gjf_data = [line for line in gjf_data if line != ''] + + atomic_types = [] + cartesian_positions = [] + + for line in gjf_data: + if line.split()[0] in elements_dict('atomic_number').keys(): + atomic_types.append(line.split()[0]) + cartesian_positions.append([float(i) for i in line.split()[1:4]]) + + cell_matrix = [] + for line in gjf_data: + if line.split()[0] == 'Tv': + cell_matrix.append([float(i) for i in line.split()[1:4]]) + + if cell_matrix != []: + self.set_cell_matrix(np.array(cell_matrix)) + + self.set_atomic_types(atomic_types) + + self.set_cartesian_positions(np.array(cartesian_positions)) + + def from_cif(self, path, file_name): + ''' + Reads a CIF file from a given path and file_name. + ''' + + # Read the cif file and get the lattice parameters and atomic positions + cif_filename = os.path.join(path, file_name.split('.')[0] + '.cif') + + cif = gemmi.cif.read_file(cif_filename).sole_block() + + a = float(cif.find_value('_cell_length_a').split('(')[0]) + b = float(cif.find_value('_cell_length_b').split('(')[0]) + c = float(cif.find_value('_cell_length_c').split('(')[0]) + beta = float(cif.find_value('_cell_angle_beta').split('(')[0]) + gamma = float(cif.find_value('_cell_angle_gamma').split('(')[0]) + alpha = float(cif.find_value('_cell_angle_alpha').split('(')[0]) + + CellParameters = [a, b, c, alpha, beta, gamma] + + AtomicTypes = list(cif.find_values('_atom_site_type_symbol')) + PosX = np.array(cif.find_values('_atom_site_fract_x')).astype(float) + PosY = np.array(cif.find_values('_atom_site_fract_y')).astype(float) + PosZ = np.array(cif.find_values('_atom_site_fract_z')).astype(float) + try: + charges = np.array(cif.find_values('_atom_site_charge')).astype(float) + charge_type = 'DDEC' + except Exception: + charges = None + charge_type = None + + self.set_cell_parameters(CellParameters) + + self.set_atomic_types(AtomicTypes) + + self.set_fractional_positions(np.array([PosX, PosY, PosZ]).T) + + if charges is not None: + self.partial_charges = {charge_type: charges} + + def as_dict(self): + ''' + Returns the structure as a dictionary. + ''' + structure_dict = { + 'chemical json': 1, + 'name': self.name, + 'formula': self.formula, + } + if self.cell_parameters is not None: + structure_dict['unit cell'] = { + 'a': self.cell_parameters[0], + 'b': self.cell_parameters[1], + 'c': self.cell_parameters[2], + 'alpha': self.cell_parameters[3], + 'beta': self.cell_parameters[4], + 'gamma': self.cell_parameters[5], + 'cellVectors': self.cell_matrix.flatten().tolist() + } + + structure_dict['atoms'] = { + 'elements': { + 'type': self.atomic_types, + 'number': self.atomic_numbers, + }, + 'coords': { + '3d': self.cartesian_positions.flatten().tolist(), + } + } + + if self.cell_parameters is not None: + structure_dict['atoms']['coords']['3dFractional'] = self.fractional_positions.flatten().tolist() + + if self.partial_charges is not None: + structure_dict['partialCharges'] = self.partial_charges + + structure_dict['properties'] = self.properties + + structure_dict['results'] = self.results + + return structure_dict + + def write_cjson(self, path, file_name): + ''' + Writes a ChemJSON file to a given path and file_name. + ''' + self.file_name = os.path.join(path, file_name.split('.')[0] + '.cjson') + with open(self.file_name, 'w') as file: + simplejson.dump(self.as_dict(), file, indent=4) diff --git a/build/lib/data/topology.py b/build/lib/data/topology.py new file mode 100644 index 00000000..2351f00f --- /dev/null +++ b/build/lib/data/topology.py @@ -0,0 +1,253 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +""" +The dictionary containing the definitions nets for a Framework buiding +""" + +import numpy as np + +TOPOLOGY_DICT = { + 'HCB': { + 'a': 2*np.cos(np.radians(30)), + 'b': 2*np.cos(np.radians(30)), + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 120, + 'vertice_connectivity': 3, + 'edge_connectivity': 0, + 'vertices': [ + {'position': [0, 0, 0], 'angle': 0}, + {'position': [0, np.sqrt(3)/3, 0], 'angle': 180} + ], + 'edges': [] + }, + 'HCB_A': { + 'a': 2*np.cos(np.radians(30))*2, + 'b': 2*np.cos(np.radians(30))*2, + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 120, + 'vertice_connectivity': 3, + 'edge_connectivity': 3, + 'vertices': [ + {'position': [0, 0, 0], 'angle': 0}, + {'position': [0, np.sqrt(3)/3, 0], 'angle': 180} + ], + 'edges': [ + {'position': [0, np.sqrt(3)/6, 0], 'angle': 0}, + {'position': [-1/4, 5*np.sqrt(3)/12, 0], 'angle': 120}, + {'position': [1/4, 5*np.sqrt(3)/12, 0], 'angle': 240} + ] + }, + 'SQL': { + 'a': 2/np.sqrt(2), + 'b': 2/np.sqrt(2), + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 90, + 'vertice_connectivity': 4, + 'edge_connectivity': 0, + 'vertices': [ + {'position': [0, 0, 0], 'angle': 0}, + {'position': [1/2, 1/2, 0], 'angle': 0} + ], + 'edges': [] + }, + 'SQL_A': { + 'a': 4/np.sqrt(2), + 'b': 4/np.sqrt(2), + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 90, + 'vertice_connectivity': 4, + 'edge_connectivity': 2, + 'vertices': [ + {'position': [0, 0, 0], 'angle': 0}, + {'position': [1/2, 1/2, 0], 'angle': 0} + ], + 'edges': [ + {'position': [1/4, 1/4, 0], 'angle': 45}, + {'position': [3/4, 1/4, 0], 'angle': 135}, + {'position': [3/4, 3/4, 0], 'angle': 225}, + {'position': [1/4, 3/4, 0], 'angle': 315}, + ] + }, + 'KGD': { + 'a': 2*np.cos(np.radians(30)), + 'b': 2*np.cos(np.radians(30)), + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 120, + 'vertice_connectivity': 6, + 'edge_connectivity': 3, + 'vertices': [ + {'position': [0, 0, 0], 'angle': 0}, + ], + 'edges': [ + {'position': [0, np.sqrt(3)/3, 0], 'angle': -180}, + {'position': [0.5, np.sqrt(3)/6, 0], 'angle': 0} + ] + }, + 'HXL_A': { + 'a': 2, + 'b': 2, + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 120, + 'vertice_connectivity': 6, + 'edge_connectivity': 0, + 'vertices': [ + {'position': [0, 0, 0], 'angle': 30}, + ], + 'edges': [ + {'position': [1/4, np.sqrt(3)/4, 0], 'angle': 30}, + {'position': [0.5, 0, 0], 'angle': 90}, + {'position': [-1/4, np.sqrt(3)/4, 0], 'angle': -30} + ] + }, + 'KGM': { + 'a': 4, + 'b': 4, + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 120, + 'vertice_connectivity': 4, + 'edge_connectivity': 0, + 'vertices': [ + {'position': [1/4, np.sqrt(3)/4, 0], 'angle': 30}, + {'position': [1/2, 0, 0], 'angle': -90}, + {'position': [-1/4, np.sqrt(3)/4, 0], 'angle': -30} + ], + 'edges': [] + }, + 'KGM_A': { + 'a': 4, + 'b': 4, + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 120, + 'vertice_connectivity': 4, + 'edge_connectivity': 2, + 'vertices': [ + {'position': [1/4, np.sqrt(3)/4, 0], 'angle': 30}, + {'position': [1/2, 0, 0], 'angle': -90}, + {'position': [-1/4, np.sqrt(3)/4, 0], 'angle': -30} + ], + 'edges': [ + {'position': [3/8, np.sqrt(3)/8, 0], 'angle': -30}, + {'position': [1/8, 3*np.sqrt(3)/8, 0], 'angle': -30}, + {'position': [5/8, np.sqrt(3)/8, 0], 'angle': 30}, + {'position': [-1/8, np.sqrt(3)/8, 0], 'angle': 30}, + {'position': [4/8, np.sqrt(3)/4, 0], 'angle': 90}, + {'position': [0, np.sqrt(3)/4, 0], 'angle': 90}, + ] + }, + 'FXT': { + 'a': 1, + 'b': 1, + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 120, + 'vertice_connectivity': 4, + 'edge_connectivity': 2, + 'vertices': [ + {'position': [1/4, 3*np.sqrt(3)/12, 0], 'angle': -15}, + {'position': [0.5, 0, 0], 'angle': 45}, + {'position': [-1/4, 3*np.sqrt(3)/12, 0], 'angle': 15} + ], + 'edges': [] + }, + 'FXT_A': { + 'a': 2, + 'b': 2, + 'c': 3.6, + 'alpha': 90, + 'beta': 90, + 'gamma': 120, + 'vertice_connectivity': 4, + 'edge_connectivity': 2, + 'vertices': [ + {'position': [1/4, 3*np.sqrt(3)/12, 0], 'angle': -15}, + {'position': [0.5, 0, 0], 'angle': 45}, + {'position': [-1/4, 3*np.sqrt(3)/12, 0], 'angle': 15} + ], + 'edges': [ + {'position': [22/64, 7*np.sqrt(3)/64, 0], 'angle': -30}, + {'position': [85/128, 7*np.sqrt(3)/64, 0], 'angle': 30}, + {'position': [4/8, 35*np.sqrt(3)/128, 0], 'angle': 90}, + {'position': [0, 29*np.sqrt(3)/128, 0], 'angle': 90}, + {'position': [21/128, 25*np.sqrt(3)/64, 0], 'angle': -30}, + {'position': [-21/128, 25*np.sqrt(3)/64, 0], 'angle': 30}, + ] + }, + 'DIA': { + 'a': 1, + 'b': 1, + 'c': 1, + 'alpha': 60, + 'beta': 60, + 'gamma': 60, + 'lattice': [[0, 1, 1], [1, 0, 1], [1, 1, 0]], + 'vertice_connectivity': 4, + 'edge_connectivity': 4, + 'vertices': [ + {'position': [0, 0, 0], 'angle': 55, 'align_v': [1, 1, 1]}, + {'position': [1/4, 1/4, 1/4], 'angle': -55, 'align_v': [-1, -1, -1]}, + ], + 'edges': [] + }, + 'DIA_A': { + 'a': 1, + 'b': 1, + 'c': 1, + 'alpha': 60, + 'beta': 60, + 'gamma': 60, + 'lattice': [[0, 1, 1], [1, 0, 1], [1, 1, 0]], + 'vertice_connectivity': 4, + 'edge_connectivity': 2, + 'vertices': [ + {'position': [0, 0, 0], 'angle': -7.5, 'align_v': [1, 1, 1]}, + {'position': [1/4, 1/4, 1/4], 'angle': 7.5, 'align_v': [-1, -1, -1]}, + ], + 'edges': [ + {'position': [1/8, 1/8, 1/8], 'angle': -8.8, 'align_v': [1, 1, 1]}, + {'position': [1/8, 3/8, 3/8], 'angle': 16, 'align_v': [-1/4, 1/4, 1/4]}, + {'position': [3/8, 1/8, 3/8], 'angle': -78, 'align_v': [1/4, -1/4, 1/4]}, + {'position': [3/8, 3/8, 1/8], 'angle': 16, 'align_v': [1/4, 1/4, -1/4]}, + ] + }, + 'BOR': { + 'a': 1, + 'b': 1, + 'c': 1, + 'alpha': 90, + 'beta': 90, + 'gamma': 90, + 'lattice': [[1, 0, 0], [0, 1, 0], [0, 0, 1]], + 'vertice_connectivity': 4, + 'edge_connectivity': 3, + 'vertices': [ + {'position': [1/2, 0, 0], 'angle': -10, 'align_v': [-1, 1, 1]}, + {'position': [0, 1/2, 0], 'angle': 40, 'align_v': [1, -1, 1]}, + {'position': [0, 0, 1/2], 'angle': -45, 'align_v': [1, 1, -1]}, + ], + 'edges': [ + {'position': [1/6, 1/6, 1/6], 'angle': 30, 'align_v': [1, 1, 1]}, + {'position': [1/6, 5/6, 5/6], 'angle': 0, 'align_v': [1/2, 1, 1]}, + {'position': [5/6, 1/6, 5/6], 'angle': 0, 'align_v': [1, 1/2, 1]}, + {'position': [5/6, 5/6, 1/6], 'angle': 10, 'align_v': [1, 1, 1]}, + ] + } + } diff --git a/build/lib/exceptions.py b/build/lib/exceptions.py new file mode 100644 index 00000000..613aac5c --- /dev/null +++ b/build/lib/exceptions.py @@ -0,0 +1,48 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +""" +This file implements custom exceptions for the pycofbuilder package. +""" + + +class BondLenghError(Exception): + """Exception raised when the bond length is shorter than 0.8 angstrom.""" + def __init__(self, atom_1=None, atom_2=None, dist=0, threshold=0.8): + self.message = 'WARNING: Atoms {} and {} are closer than {} A, {}'.format(atom_1, + atom_2, + dist, + threshold) + + def __str__(self): + return str(self.message) + + +class BBConnectivityError(Exception): + """Exception raised when the building block connectivity is not valid.""" + def __init__(self, connectivity=None, found_connectivity=None): + self.message = 'ERROR: The building block connectivity should be {} buy is {}'.format(connectivity, + found_connectivity) + + def __str__(self): + return str(self.message) + + +class ConnectionGroupError(Exception): + """Exception raised when the connection group is not valid.""" + def __init__(self, conn_1=None, conn_2=None): + self.message = 'ERROR: The connection group {} not compatible with {}'.format(conn_1, + conn_2) + + def __str__(self): + return str(self.message) + + +class MissingXError(Exception): + """Exception raised when the building block has missing X atoms.""" + def __init__(self): + self.message = 'No X points found in the structure!' + + def __str__(self): + return str(self.message) diff --git a/build/lib/framework.py b/build/lib/framework.py new file mode 100644 index 00000000..20376e51 --- /dev/null +++ b/build/lib/framework.py @@ -0,0 +1,4494 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +""" +The Framework class implements definitions and methods for a Framework buiding +""" + +import os +import copy +import numpy as np + +# Import pymatgen +from pymatgen.core import Lattice, Structure +from pymatgen.symmetry.analyzer import SpacegroupAnalyzer + +from scipy.spatial.transform import Rotation as R + +# Import pycofbuilder exceptions +from pycofbuilder.exceptions import (BondLenghError, + BBConnectivityError) + +# Import pycofbuilder building_block +from pycofbuilder.building_block import BuildingBlock + +# Import pycofbuilder topology data +from pycofbuilder.data.topology import TOPOLOGY_DICT + +# Import pycofbuilder tools +from pycofbuilder.tools import (get_bond_atom, + cell_to_cellpar, + cellpar_to_cell, + rotation_matrix_from_vectors, + unit_vector, + angle, + get_framework_symm_text) + +# Import pycofbuilder io_tools +from pycofbuilder.io_tools import (save_json, + save_chemjson, + save_cif, + save_xyz, + save_turbomole, + save_vasp, + save_xsf, + save_pdb, + save_pqr, + save_qe, + save_gjf) + +from pycofbuilder.logger import create_logger + + +class Framework(): + """ + A class used to represent a Covalent Organic Framework as a reticular entity. + + ... + + Attributes + ---------- + name : str + Name of the material + out_dir : str + Path to save the results. + If not defined, a `out` folder will be created in the current directory. + verbosity : str + Control the printing options. Can be 'none', 'normal', or 'debug'. + Default: 'normal' + save_bb : bool + Control the saving of the building blocks. + Default: True + lib_path : str + Path for saving the building block files. + If not defined, a `building_blocks` folder will be created in the current directory. + topology : str = None + dimention : str = None + lattice : str = None + lattice_sgs : str = None + space_group : str = None + space_group_n : str = None + stacking : str = None + mass : str = None + composition : str = None + charge : int = 0 + multiplicity : int = 1 + chirality : bool = False + atom_types : list = [] + atom_pos : list = [] + lattice : list = [[], [], []] + symm_tol : float = 0.2 + angle_tol : float = 0.2 + dist_threshold : float = 0.8 + available_2D_topologies : list + List of available 2D topologies + available_3D_topologies : list + List of available 3D topologies + available_topologies : list + List of all available topologies + available_stacking : list + List of available stakings for all 2D topologies + lib_bb : str + String with the name of the folder containing the building block files + Default: bb_lib + """ + + def __init__(self, name: str = None, **kwargs): + + self.name: str = name + + self.out_path: str = kwargs.get('out_dir', os.path.join(os.getcwd(), 'out')) + self.save_bb: bool = kwargs.get('save_bb', True) + self.bb_out_path: str = kwargs.get('bb_out_path', os.path.join(self.out_path, 'building_blocks')) + + self.logger = create_logger(level=kwargs.get('log_level', 'info'), + format=kwargs.get('log_format', 'simple'), + save_to_file=kwargs.get('save_to_file', False), + log_filename=kwargs.get('log_filename', 'pycofbuilder.log')) + + self.symm_tol = kwargs.get('symm_tol', 0.1) + self.angle_tol = kwargs.get('angle_tol', 0.5) + self.dist_threshold = kwargs.get('dist_threshold', 0.8) + + self.bb1_name = None + self.bb2_name = None + self.topology = None + self.stacking = None + self.smiles = None + + self.atom_types = [] + self.atom_pos = [] + self.atom_labels = [] + self.cellMatrix = np.eye(3) + self.cellParameters = np.array([1, 1, 1, 90, 90, 90]).astype(float) + + self.lattice_sgs = None + self.space_group = None + self.space_group_n = None + + self.dimention = None + self.n_atoms = self.get_n_atoms() + self.mass = None + self.composition = None + self.charge = 0 + self.multiplicity = 1 + self.chirality = False + + self.available_2D_top = ['HCB', 'HCB_A', + 'SQL', 'SQL_A', + 'KGD', + 'HXL', 'HXL_A', + 'FXT', 'FXT_A'] + + # To add: ['dia', 'bor', 'srs', 'pts', 'ctn', 'rra', 'fcc', 'lon', 'stp', 'acs', 'tbo', 'bcu', 'fjh', 'ceq'] + self.available_3D_top = ['DIA', 'DIA_A', 'BOR'] # Temporary + self.available_topologies = self.available_2D_top + self.available_3D_top + + # Define available stackings for all 2D topologies + self.available_stacking = { + 'HCB': ['A', 'AA', 'AB1', 'AB2', 'AAl', 'AAt', 'ABC1', 'ABC2'], + 'HCB_A': ['A', 'AA', 'AB1', 'AB2', 'AAl', 'AAt', 'ABC1', 'ABC2'], + 'SQL': ['A', 'AA', 'AB1', 'AB2', 'AAl', 'AAt', 'ABC1', 'ABC2'], + 'SQL_A': ['A', 'AA', 'AB1', 'AB2', 'AAl', 'AAt', 'ABC1', 'ABC2'], + 'KGD': ['A', 'AA', 'AB1', 'AB2', 'AAl', 'AAt', 'ABC1', 'ABC2'], + 'HXL_A': ['A', 'AA', 'AB1', 'AB2', 'AAl', 'AAt', 'ABC1', 'ABC2'], + 'FXT': ['A', 'AA', 'AB1', 'AB2', 'AAl', 'AAt', 'ABC1', 'ABC2'], + 'FXT_A': ['A', 'AA', 'AB1', 'AB2', 'AAl', 'AAt', 'ABC1', 'ABC2'], + 'DIA': [str(i + 1) for i in range(15)], + 'DIA_A': [str(i + 1) for i in range(15)], + 'BOR': [str(i + 1) for i in range(15)] + } + + if self.name is not None: + self.from_name(self.name) + + def __str__(self) -> str: + return self.as_string() + + def __repr__(self) -> str: + return f'Framework({self.bb1_name}, {self.bb2_name}, {self.topology}, {self.stacking})' + + def as_string(self) -> str: + """ + Returns a string with the Framework information. + """ + + fram_str = f'Name: {self.name}\n' + + # Get the formula of the framework + + if self.composition is not None: + fram_str += f'Full Formula ({self.composition})\n' + fram_str += f'Reduced Formula: ({self.composition})\n' + else: + fram_str += 'Full Formula ()\n' + fram_str += 'Reduced Formula: \n' + + fram_str += 'abc : {:11.6f} {:11.6f} {:11.6f}\n'.format(*self.cellParameters[:3]) + fram_str += 'angles: {:11.6f} {:11.6f} {:11.6f}\n'.format(*self.cellParameters[3:]) + fram_str += 'A: {:11.6f} {:11.6f} {:11.6f}\n'.format(*self.cellMatrix[0]) + fram_str += 'B: {:11.6f} {:11.6f} {:11.6f}\n'.format(*self.cellMatrix[1]) + fram_str += 'C: {:11.6f} {:11.6f} {:11.6f}\n'.format(*self.cellMatrix[2]) + + fram_str += f'Cartesian Sites ({self.n_atoms})\n' + fram_str += ' # Type a b c label\n' + fram_str += '--- ---- -------- -------- -------- -------\n' + + for i in range(len(self.atom_types)): + fram_str += '{:3d} {:4s} {:8.5f} {:8.5f} {:8.5f} {:>7}\n'.format(i, + self.atom_types[i], + self.atom_pos[i][0], + self.atom_pos[i][1], + self.atom_pos[i][2], + self.atom_labels[i]) + + return fram_str + + def get_n_atoms(self) -> int: + ''' Returns the number of atoms in the unitary cell''' + return len(self.atom_types) + + def get_available_topologies(self, dimensionality: str = 'all', print_result: bool = True): + """ + Get the available topologies implemented in the class. + + Parameters + ---------- + + dimensionality : str, optional + The dimensionality of the topologies to be printed. Can be 'all', '2D' or '3D'. + Default: 'all' + print_result: bool, optional + If True, the available topologies are printed. + + Returns + ------- + dimensionality_list: list + A list with the available topologies. + """ + + dimensionality_error = f'Dimensionality must be one of the following: all, 2D, 3D, not {dimensionality}' + assert dimensionality in ['all', '2D', '3D'], dimensionality_error + + dimensionality_list = [] + + if dimensionality == 'all' or dimensionality == '2D': + if print_result: + print('Available 2D Topologies:') + for i in self.available_2D_top: + if print_result: + print(i.upper()) + dimensionality_list.append(i) + + if dimensionality == 'all' or dimensionality == '3D': + if print_result: + print('Available 3D Topologies:') + for i in self.available_3D_top: + if print_result: + print(i.upper()) + dimensionality_list.append(i) + + return dimensionality_list + + def check_name_concistency(self, FrameworkName) -> tuple[str, str, str, str]: + """ + Checks if the name is in the correct format and returns a + tuple with the building blocks names, the net and the stacking. + + In case the name is not in the correct format, an error is raised. + + Parameters + ---------- + FrameworkName : str, required + The name of the COF to be created + + Returns + ------- + tuple[str, str, str, str] + A tuple with the building blocks names, the net and the stacking. + """ + + string_error = 'FrameworkName must be in the format: BB1_BB2_Net_Stacking' + assert isinstance(FrameworkName, str), string_error + + name_error = 'FrameworkName must be in the format: BB1_BB2_Net_Stacking' + assert len(FrameworkName.split('-')) == 4, name_error + + bb1_name, bb2_name, Net, Stacking = FrameworkName.split('-') + + net_error = f'{Net} not in the available list: {self.available_topologies}' + assert Net in self.available_topologies, net_error + + stacking_error = f'{Stacking} not in the available list: {self.available_stacking[Net]}' + assert Stacking in self.available_stacking[Net], stacking_error + + return bb1_name, bb2_name, Net, Stacking + + def from_name(self, FrameworkName, **kwargs) -> None: + """Creates a COF from a given FrameworkName. + + Parameters + ---------- + FrameworkName : str, required + The name of the COF to be created + + Returns + ------- + COF : Framework + The COF object + """ + bb1_name, bb2_name, Net, Stacking = self.check_name_concistency(FrameworkName) + + bb1 = BuildingBlock(name=bb1_name, bb_out_path=self.bb_out_path, save_bb=self.save_bb) + bb2 = BuildingBlock(name=bb2_name, bb_out_path=self.bb_out_path, save_bb=self.save_bb) + + self.from_building_blocks(bb1, bb2, Net, Stacking, **kwargs) + + def from_building_blocks(self, + bb1: BuildingBlock, + bb2: BuildingBlock, + net: str, + stacking: str, + **kwargs): + """Creates a COF from the building blocks. + + Parameters + ---------- + BB1 : BuildingBlock, required + The first building block + BB2 : BuildingBlock, required + The second building block + Net : str, required + The network of the COF + Stacking : str, required + The stacking of the COF + + Returns + ------- + COF : Framework + The COF object + """ + self.name = f'{bb1.name}-{bb2.name}-{net}-{stacking}' + self.bb1_name = bb1.name + self.bb2_name = bb2.name + self.topology = net + self.stacking = stacking + + # Check if the BB1 has the smiles attribute + if hasattr(bb1, 'smiles') and hasattr(bb2, 'smiles'): + self.smiles = f'{bb1.smiles}.{bb2.smiles}' + else: + print('WARNING: The smiles attribute is not available for the building blocks') + + net_build_dict = { + 'HCB': self.create_hcb_structure, + 'HCB_A': self.create_hcb_a_structure, + 'SQL': self.create_sql_structure, + 'SQL_A': self.create_sql_a_structure, + 'KGD': self.create_kgd_structure, + # 'HXL': self.create_hxl_structure, + 'HXL_A': self.create_hxl_a_structure, + 'FXT': self.create_fxt_structure, + 'FXT_A': self.create_fxt_a_structure, + 'DIA': self.create_dia_structure, + 'DIA_A': self.create_dia_a_structure, + 'BOR': self.create_bor_structure + } + + result = net_build_dict[net](bb1, bb2, stacking, **kwargs) + + return result + + def save(self, fmt: str = 'cif', supercell: list = [1, 1, 1], save_dir=None, primitive=False) -> None: + ''' + Save the structure in a specif file format. + + Parameters + ---------- + fmt : str, optional + The file format to be saved + Can be `json`, `cif`, `xyz`, `turbomole`, `vasp`, `xsf`, `pdb`, `pqr`, `qe`. + Default: 'cif' + supercell : list, optional + The supercell to be used to save the structure. + Default: [1,1,1] + save_dir : str, optional + The path to save the structure. By default, the structure is saved in a + `out` folder created in the current directory. + primitive : bool, optional + If True, the primitive cell is saved. Otherwise, the conventional cell is saved. + Default: False + ''' + + save_dict = { + 'json': save_json, + 'cjson': save_chemjson, + 'cif': save_cif, + 'xyz': save_xyz, + 'turbomole': save_turbomole, + 'vasp': save_vasp, + 'xsf': save_xsf, + 'pdb': save_pdb, + 'pqr': save_pqr, + 'qe': save_qe, + 'gjf': save_gjf + } + + file_format_error = f'Format must be one of the following: {save_dict.keys()}' + assert fmt in save_dict.keys(), file_format_error + + if primitive: + structure = self.prim_structure + + else: + structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + final_structure = structure.make_supercell(supercell, in_place=False) + + structure_dict = final_structure.as_dict() + + cell = structure_dict['lattice']['matrix'] + + atom_types = [site['species'][0]['element'] for site in structure_dict['sites']] + atom_labels = [site['properties']['source'] for site in structure_dict['sites']] + atom_pos = [site['xyz'] for site in structure_dict['sites']] + + if save_dir is None: + save_path = self.out_path + else: + save_path = save_dir + + save_dict[fmt](path=save_path, + file_name=self.name, + cell=cell, + atom_types=atom_types, + atom_labels=atom_labels, + atom_pos=atom_pos) + +# --------------- Net creation methods -------------------------- # + + def create_hcb_structure(self, + BB_T3_A, + BB_T3_B, + stacking: str = 'AA', + slab: float = 10.0, + shift_vector: list = [1.0, 1.0, 0], + tilt_angle: float = 5.0): + """Creates a COF with HCB network. + + The HCB net is composed of two tripodal building blocks. + + Parameters + ---------- + BB_T3_1 : BuildingBlock, required + The BuildingBlock object of the tripodal Buiding Block A + BB_T3_2 : BuildingBlock, required + The BuildingBlock object of the tripodal Buiding Block B + stacking : str, optional + The stacking pattern of the COF layers (default is 'AA') + slab : float, optional + Default parameter for the interlayer slab (default is 10.0) + shift_vector: list, optional + Shift vector for the AAl and AAt stakings (defatult is [1.0,1.0,0]) + tilt_angle: float, optional + Tilt angle for the AAt staking in degrees (default is 5.0) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + + connectivity_error = 'Building block {} must present connectivity {} not {}' + + if BB_T3_A.connectivity != 3: + self.logger.error(connectivity_error.format('A', 3, BB_T3_A.connectivity)) + raise BBConnectivityError(3, BB_T3_A.connectivity) + if BB_T3_B.connectivity != 3: + self.logger.error(connectivity_error.format('B', 3, BB_T3_B.connectivity)) + raise BBConnectivityError(3, BB_T3_B.connectivity) + + self.name = f'{BB_T3_A.name}-{BB_T3_B.name}-HCB-{stacking}' + self.topology = 'HCB' + self.staking = stacking + self.dimension = 2 + + self.charge = BB_T3_A.charge + BB_T3_B.charge + self.chirality = BB_T3_A.chirality or BB_T3_B.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_T3_A.conector, BB_T3_B.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_T3_A.conector, + BB_T3_B.conector)) + + # Replace "X" the building block + BB_T3_A.replace_X(bond_atom) + + # Remove the "X" atoms from the the building block + BB_T3_A.remove_X() + BB_T3_B.remove_X() + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = BB_T3_A.size[0] + BB_T3_B.size[0] + + # Calculate the delta size to add to the c parameter + delta_a = abs(max(np.transpose(BB_T3_A.atom_pos)[2])) + abs(min(np.transpose(BB_T3_A.atom_pos)[2])) + delta_b = abs(max(np.transpose(BB_T3_B.atom_pos)[2])) + abs(min(np.transpose(BB_T3_B.atom_pos)[2])) + + delta_max = max([delta_a, delta_b]) + + # Calculate the cell parameters + a = topology_info['a'] * size + b = topology_info['b'] * size + c = topology_info['c'] + delta_max + alpha = topology_info['alpha'] + beta = topology_info['beta'] + gamma = topology_info['gamma'] + + if self.stacking == 'A': + c = slab + + # Create the lattice + self.cellMatrix = Lattice.from_parameters(a, b, c, alpha, beta, gamma) + self.cellParameters = np.array([a, b, c, alpha, beta, gamma]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Add the A1 building blocks to the structure + vertice_data = topology_info['vertices'][0] + self.atom_types += BB_T3_A.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', + vertice_data['angle'], + degrees=True).as_matrix() + + rotated_pos = np.dot(BB_T3_A.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C1' if i == 'C' else i for i in BB_T3_A.atom_labels] + + # Add the A2 building block to the structure + vertice_data = topology_info['vertices'][1] + self.atom_types += BB_T3_B.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', + vertice_data['angle'], + degrees=True).as_matrix() + + rotated_pos = np.dot(BB_T3_B.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C2' if i == 'C' else i for i in BB_T3_B.atom_labels] + + # Creates a pymatgen structure + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + # Translates the structure to the center of the cell + StartingFramework.translate_sites( + range(len(StartingFramework.as_dict()['sites'])), + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + self.cellParameters = cell_to_cellpar(self.cellMatrix) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + + if stacking == 'A' or stacking == 'AA': + stacked_structure = StartingFramework + + if stacking == 'AB1': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [2/3, 1/3, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AB2': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/2, 0, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/2, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC1': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (2/3, 1/3, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (4/3, 2/3, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC2': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 0, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 0, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AAl': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + sv = np.array(shift_vector) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos + sv)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + # Create AA tilted stacking. + if stacking == 'AAt': + cell = StartingFramework.as_dict()['lattice'] + + # Shift the cell by the tilt angle + a_cell = cell['a'] + b_cell = cell['b'] + c_cell = cell['c'] * 2 + alpha = cell['alpha'] - tilt_angle + beta = cell['beta'] - tilt_angle + gamma = cell['gamma'] + + self.cellMatrix = cellpar_to_cell([a_cell, b_cell, c_cell, alpha, beta, gamma]) + self.cellParameters = np.array([a_cell, b_cell, c_cell, alpha, beta, gamma]).astype(float) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = stacked_structure.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + self.cellParameters = cell_to_cellpar(self.cellMatrix) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = stacked_structure.formula + + dist_matrix = StartingFramework.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(stacked_structure, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_refined_structure(keep_site_properties=True) + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_hcb_a_structure(self, + BB_T3: str, + BB_L2: str, + stacking: str = 'AA', + slab: float = 10.0, + shift_vector: list = [1.0, 1.0, 0], + tilt_angle: float = 5.0): + """Creates a COF with HCB-A network. + + The HCB-A net is composed of one tripodal and one linear building blocks. + + Parameters + ---------- + BB_T3 : BuildingBlock, required + The BuildingBlock object of the tripodal Buiding Block + BB_L2 : BuildingBlock, required + The BuildingBlock object of the linear Buiding Block + stacking : str, optional + The stacking pattern of the COF layers (default is 'AA') + c_parameter_base : float, optional + The base value for interlayer distance in angstroms (default is 3.6) + print_result : bool, optional + Parameter for the control for printing the result (default is True) + slab : float, optional + Default parameter for the interlayer slab (default is 10.0) + shift_vector: list, optional + Shift vector for the AAl and AAt stakings (defatult is [1.0,1.0,0]) + tilt_angle: float, optional + Tilt angle for the AAt staking in degrees (default is 5.0) + + Returns + ------- + list + A list of strings containing: + 1. the structure name + 2. lattice type + 3. hall symbol of the cristaline structure + 4. space group + 5. number of the space group, + 6. number of operation symmetry + """ + + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_T3.connectivity != 3: + self.logger.error(connectivity_error.format('A', 3, BB_T3.connectivity)) + raise BBConnectivityError(3, BB_T3.connectivity) + if BB_L2.connectivity != 2: + self.logger.error(connectivity_error.format('B', 3, BB_L2.connectivity)) + raise BBConnectivityError(2, BB_L2.connectivity) + + self.name = f'{BB_T3.name}-{BB_L2.name}-HCB_A-{stacking}' + self.topology = 'HCB_A' + self.staking = stacking + self.dimension = 2 + + self.charge = BB_L2.charge + BB_T3.charge + self.chirality = BB_L2.chirality or BB_T3.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_T3.conector, BB_L2.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_T3.conector, + BB_L2.conector)) + + # Replace "X" the building block + BB_L2.replace_X(bond_atom) + + # Remove the "X" atoms from the the building block + BB_T3.remove_X() + BB_L2.remove_X() + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = BB_T3.size[0] + BB_L2.size[0] + + # Calculate the delta size to add to the c parameter + delta_a = abs(max(np.transpose(BB_T3.atom_pos)[2])) + abs(min(np.transpose(BB_T3.atom_pos)[2])) + delta_b = abs(max(np.transpose(BB_L2.atom_pos)[2])) + abs(min(np.transpose(BB_L2.atom_pos)[2])) + + delta_max = max([delta_a, delta_b]) + + # Calculate the cell parameters + a = topology_info['a'] * size + b = topology_info['b'] * size + c = topology_info['c'] + delta_max + alpha = topology_info['alpha'] + beta = topology_info['beta'] + gamma = topology_info['gamma'] + + if self.stacking == 'A': + c = slab + + # Create the lattice + self.cellMatrix = Lattice.from_parameters(a, b, c, alpha, beta, gamma) + self.cellParameters = np.array([a, b, c, alpha, beta, gamma]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Add the building blocks to the structure + for vertice_data in topology_info['vertices']: + self.atom_types += BB_T3.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', vertice_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_T3.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C1' if i == 'C' else i for i in BB_T3.atom_labels] + + # Add the building blocks to the structure + for edge_data in topology_info['edges']: + self.atom_types += BB_L2.atom_types + + R_Matrix = R.from_euler('z', edge_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_L2.atom_pos, R_Matrix) + np.array(edge_data['position'])*a + + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C2' if i == 'C' else i for i in BB_L2.atom_labels] + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + # Translates the structure to the center of the cell + StartingFramework.translate_sites( + range(len(StartingFramework.as_dict()['sites'])), + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + + if stacking == 'A' or stacking == 'AA': + stacked_structure = StartingFramework + + if stacking == 'AB1': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [2/3, 1/3, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AB2': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/2, 0, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/2, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC1': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (2/3, 1/3, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (4/3, 2/3, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC2': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 0, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 0, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AAl': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + sv = np.array(shift_vector) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos + sv)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + # Create AA tilted stacking. + if stacking == 'AAt': + cell = StartingFramework.as_dict()['lattice'] + + # Shift the cell by the tilt angle + a_cell = cell['a'] + b_cell = cell['b'] + c_cell = cell['c'] * 2 + alpha = cell['alpha'] - tilt_angle + beta = cell['beta'] - tilt_angle + gamma = cell['gamma'] + + self.cellMatrix = cellpar_to_cell([a_cell, b_cell, c_cell, alpha, beta, gamma]) + self.cellParameters = np.array([a_cell, b_cell, c_cell, alpha, beta, gamma]).astype(float) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = stacked_structure.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = stacked_structure.formula + + dist_matrix = stacked_structure.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(stacked_structure, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_refined_structure(keep_site_properties=True) + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_sql_structure(self, + BB_S4_A: str, + BB_S4_B: str, + stacking: str = 'AA', + slab: float = 10.0, + shift_vector: list = [1.0, 1.0, 0], + tilt_angle: float = 5.0): + """Creates a COF with SQL network. + + The SQL net is composed of two tetrapodal building blocks. + + Parameters + ---------- + BB_S4_A : BuildingBlock, required + The BuildingBlock object of the tetrapodal Buiding Block A + BB_S4_B : BuildingBlock, required + The BuildingBlock object of the tetrapodal Buiding Block B + stacking : str, optional + The stacking pattern of the COF layers (default is 'AA') + print_result : bool, optional + Parameter for the control for printing the result (default is True) + slab : float, optional + Default parameter for the interlayer slab (default is 10.0) + shift_vector: list, optional + Shift vector for the AAl and AAt stakings (defatult is [1.0,1.0,0]) + tilt_angle: float, optional + Tilt angle for the AAt staking in degrees (default is 5.0) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_S4_A.connectivity != 4: + self.logger.error(connectivity_error.format('A', 4, BB_S4_A.connectivity)) + raise BBConnectivityError(4, BB_S4_A.connectivity) + if BB_S4_B.connectivity != 4: + self.logger.error(connectivity_error.format('B', 4, BB_S4_B.connectivity)) + raise BBConnectivityError(4, BB_S4_B.connectivity) + + self.name = f'{BB_S4_A.name}-{BB_S4_B.name}-SQL-{stacking}' + self.topology = 'SQL' + self.staking = stacking + self.dimension = 2 + + self.charge = BB_S4_A.charge + BB_S4_B.charge + self.chirality = BB_S4_A.chirality or BB_S4_B.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_S4_A.conector, BB_S4_B.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_S4_A.conector, + BB_S4_B.conector)) + + # Replace "X" the building block + BB_S4_A.replace_X(bond_atom) + + # Remove the "X" atoms from the the building block + BB_S4_A.remove_X() + BB_S4_B.remove_X() + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = BB_S4_A.size[0] + BB_S4_B.size[0] + + # Calculate the delta size to add to the c parameter + delta_a = abs(max(np.transpose(BB_S4_A.atom_pos)[2])) + abs(min(np.transpose(BB_S4_B.atom_pos)[2])) + delta_b = abs(max(np.transpose(BB_S4_A.atom_pos)[2])) + abs(min(np.transpose(BB_S4_B.atom_pos)[2])) + + delta_max = max([delta_a, delta_b]) + + # Calculate the cell parameters + a = topology_info['a'] * size + b = topology_info['b'] * size + c = topology_info['c'] + delta_max + alpha = topology_info['alpha'] + beta = topology_info['beta'] + gamma = topology_info['gamma'] + + if self.stacking == 'A': + c = slab + + # Create the lattice + self.cellMatrix = Lattice.from_parameters(a, b, c, alpha, beta, gamma) + self.cellParameters = np.array([a, b, c, alpha, beta, gamma]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Add the first building block to the structure + vertice_data = topology_info['vertices'][0] + self.atom_types += BB_S4_A.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', vertice_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_S4_A.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C1' if i == 'C' else i for i in BB_S4_A.atom_labels] + + # Add the second building block to the structure + vertice_data = topology_info['vertices'][1] + self.atom_types += BB_S4_B.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', vertice_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_S4_B.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C2' if i == 'C' else i for i in BB_S4_B.atom_labels] + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + # Translates the structure to the center of the cell + StartingFramework.translate_sites( + range(len(StartingFramework.as_dict()['sites'])), + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + + if stacking == 'A' or stacking == 'AA': + stacked_structure = StartingFramework + + if stacking == 'AB1': + + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/4, 1/4, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/4, 1/4, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AB2': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/2, 0, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/2, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC1': + + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (1/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 1/3, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 2/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 2/3, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC2': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (1/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 0, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 2/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 0, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + # Create AAl stacking. Tetragonal cell with two sheets + # per cell shifited by the shift_vector in angstroms. + if stacking == 'AAl': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + sv = np.array(shift_vector) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos + sv)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + # Create AA tilted stacking. + # Tilted tetragonal cell with two sheets per cell tilted by tilt_angle. + if stacking == 'AAt': + cell = StartingFramework.as_dict()['lattice'] + + # Shift the cell by the tilt angle + a_cell = cell['a'] + b_cell = cell['b'] + c_cell = cell['c'] * 2 + alpha = cell['alpha'] - tilt_angle + beta = cell['beta'] - tilt_angle + gamma = cell['gamma'] + + self.cellMatrix = cellpar_to_cell([a_cell, b_cell, c_cell, alpha, beta, gamma]) + self.cellParameters = np.array([a_cell, b_cell, c_cell, alpha, beta, gamma]).astype(float) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = stacked_structure.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = stacked_structure.formula + + dist_matrix = stacked_structure.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(stacked_structure, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_refined_structure(keep_site_properties=True) + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_sql_a_structure(self, + BB_S4: str, + BB_L2: str, + stacking: str = 'AA', + c_parameter_base: float = 3.6, + slab: float = 10.0, + shift_vector: list = [1.0, 1.0, 0], + tilt_angle: float = 5.0): + """Creates a COF with SQL-A network. + + The SQL-A net is composed of one tetrapodal and one linear building blocks. + + Parameters + ---------- + BB_S4 : BuildingBlock, required + The BuildingBlock object of the tetrapodal Buiding Block + BB_L2 : BuildingBlock, required + The BuildingBlock object of the bipodal Buiding Block + stacking : str, optional + The stacking pattern of the COF layers (default is 'AA') + slab : float, optional + Default parameter for the interlayer slab (default is 10.0) + shift_vector: list, optional + Shift vector for the AAl and AAt stakings (defatult is [1.0,1.0,0]) + tilt_angle: float, optional + Tilt angle for the AAt staking in degrees (default is 5.0) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_S4.connectivity != 4: + self.logger.error(connectivity_error.format('A', 4, BB_S4.connectivity)) + raise BBConnectivityError(4, BB_S4.connectivity) + if BB_L2.connectivity != 2: + self.logger.error(connectivity_error.format('B', 3, BB_L2.connectivity)) + raise BBConnectivityError(2, BB_L2.connectivity) + + self.name = f'{BB_S4.name}-{BB_L2.name}-SQL_A-{stacking}' + self.topology = 'SQL_A' + self.staking = stacking + self.dimension = 2 + + self.charge = BB_S4.charge + BB_L2.charge + self.chirality = BB_S4.chirality or BB_L2.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_S4.conector, BB_L2.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_S4.conector, + BB_L2.conector)) + + # Replace "X" the building block + BB_L2.replace_X(bond_atom) + + # Remove the "X" atoms from the the building block + BB_S4.remove_X() + BB_L2.remove_X() + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = BB_S4.size[0] + BB_L2.size[0] + + # Calculate the delta size to add to the c parameter + delta_a = abs(max(np.transpose(BB_S4.atom_pos)[2])) + abs(min(np.transpose(BB_S4.atom_pos)[2])) + delta_b = abs(max(np.transpose(BB_L2.atom_pos)[2])) + abs(min(np.transpose(BB_L2.atom_pos)[2])) + + delta_max = max([delta_a, delta_b]) + + # Calculate the cell parameters + a = topology_info['a'] * size + b = topology_info['b'] * size + c = topology_info['c'] + delta_max + alpha = topology_info['alpha'] + beta = topology_info['beta'] + gamma = topology_info['gamma'] + + if self.stacking == 'A': + c = slab + + # Create the lattice + self.cellMatrix = Lattice.from_parameters(a, b, c, alpha, beta, gamma) + self.cellParameters = np.array([a, b, c, alpha, beta, gamma]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Add the building blocks to the structure + for vertice_data in topology_info['vertices']: + self.atom_types += BB_S4.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', vertice_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_S4.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C1' if i == 'C' else i for i in BB_S4.atom_labels] + + # Add the building blocks to the structure + for edge_data in topology_info['edges']: + self.atom_types += BB_L2.atom_types + edge_pos = np.array(edge_data['position'])*a + + R_Matrix = R.from_euler('z', edge_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_L2.atom_pos, R_Matrix) + edge_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C2' if i == 'C' else i for i in BB_L2.atom_labels] + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + # Translates the structure to the center of the cell + StartingFramework.translate_sites( + range(len(StartingFramework.as_dict()['sites'])), + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + + if stacking == 'A' or stacking == 'AA': + stacked_structure = StartingFramework + + if stacking == 'AB1': + + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/4, 1/4, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/4, 1/4, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AB2': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/2, 0, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/2, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC1': + + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (1/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 1/3, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 2/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 2/3, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC2': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (1/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 0, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 2/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 0, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + # Create AAl stacking. Tetragonal cell with two sheets + # per cell shifited by the shift_vector in angstroms. + if stacking == 'AAl': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + sv = np.array(shift_vector) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos + sv)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + # Create AA tilted stacking. + # Tilted tetragonal cell with two sheets per cell tilted by tilt_angle. + if stacking == 'AAt': + cell = StartingFramework.as_dict()['lattice'] + + # Shift the cell by the tilt angle + a_cell = cell['a'] + b_cell = cell['b'] + c_cell = cell['c'] * 2 + alpha = cell['alpha'] - tilt_angle + beta = cell['beta'] - tilt_angle + gamma = cell['gamma'] + + self.cellMatrix = cellpar_to_cell([a_cell, b_cell, c_cell, alpha, beta, gamma]) + self.cellParameters = np.array([a_cell, b_cell, c_cell, alpha, beta, gamma]).astype(float) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = stacked_structure.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = stacked_structure.formula + + dist_matrix = stacked_structure.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(stacked_structure, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_refined_structure(keep_site_properties=True) + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_kgd_structure(self, + BB_H6: str, + BB_T3: str, + stacking: str = 'AA', + print_result: bool = True, + slab: float = 10.0, + shift_vector: list = [1.0, 1.0, 0], + tilt_angle: float = 5.0): + """Creates a COF with KGD network. + + The KGD net is composed of one hexapodal and one tripodal building blocks. + + Parameters + ---------- + BB_H6 : BuildingBlock, required + The BuildingBlock object of the hexapodal Buiding Block + BB_T3 : BuildingBlock, required + The BuildingBlock object of the tripodal Buiding Block + stacking : str, optional + The stacking pattern of the COF layers (default is 'AA') + c_parameter_base : float, optional + The base value for interlayer distance in angstroms (default is 3.6) + print_result : bool, optional + Parameter for the control for printing the result (default is True) + slab : float, optional + Default parameter for the interlayer slab (default is 10.0) + shift_vector: list, optional + Shift vector for the AAl and AAt stakings (defatult is [1.0,1.0,0]) + tilt_angle: float, optional + Tilt angle for the AAt staking in degrees (default is 5.0) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_H6.connectivity != 6: + self.logger.error(connectivity_error.format('A', 6, BB_H6.connectivity)) + raise BBConnectivityError(6, BB_H6.connectivity) + if BB_T3.connectivity != 3: + self.logger.error(connectivity_error.format('B', 3, BB_T3.connectivity)) + raise BBConnectivityError(3, BB_T3.connectivity) + + self.name = f'{BB_H6.name}-{BB_T3.name}-KGD-{stacking}' + self.topology = 'KGD' + self.staking = stacking + self.dimension = 2 + + self.charge = BB_H6.charge + BB_T3.charge + self.chirality = BB_H6.chirality or BB_T3.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_H6.conector, BB_T3.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_H6.conector, + BB_T3.conector)) + + # Replace "X" the building block + BB_H6.replace_X(bond_atom) + + # Remove the "X" atoms from the the building block + BB_H6.remove_X() + BB_T3.remove_X() + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = BB_H6.size[0] + BB_T3.size[0] + + # Calculate the delta size to add to the c parameter + delta_a = abs(max(np.transpose(BB_H6.atom_pos)[2])) + abs(min(np.transpose(BB_H6.atom_pos)[2])) + delta_b = abs(max(np.transpose(BB_T3.atom_pos)[2])) + abs(min(np.transpose(BB_T3.atom_pos)[2])) + + delta_max = max([delta_a, delta_b]) + + # Calculate the cell parameters + a = topology_info['a'] * size + b = topology_info['b'] * size + c = topology_info['c'] + delta_max + alpha = topology_info['alpha'] + beta = topology_info['beta'] + gamma = topology_info['gamma'] + + if self.stacking == 'A': + c = slab + + # Create the lattice + self.cellMatrix = Lattice.from_parameters(a, b, c, alpha, beta, gamma) + self.cellParameters = np.array([a, b, c, alpha, beta, gamma]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Add the building blocks to the structure + for vertice_data in topology_info['vertices']: + self.atom_types += BB_H6.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', + vertice_data['angle'], + degrees=True).as_matrix() + + rotated_pos = np.dot(BB_H6.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C1' if i == 'C' else i for i in BB_H6.atom_labels] + + # Add the building blocks to the structure + for edge_data in topology_info['edges']: + self.atom_types += BB_T3.atom_types + edge_pos = np.array(edge_data['position'])*a + + R_Matrix = R.from_euler('z', + edge_data['angle'], + degrees=True).as_matrix() + + rotated_pos = np.dot(BB_T3.atom_pos, R_Matrix) + edge_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C2' if i == 'C' else i for i in BB_T3.atom_labels] + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + # Translates the structure to the center of the cell + StartingFramework.translate_sites( + range(len(StartingFramework.as_dict()['sites'])), + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + + if stacking == 'A' or stacking == 'AA': + stacked_structure = StartingFramework + + if stacking == 'AB1': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [2/3, 1/3, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AB2': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/2, 0, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/2, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC1': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (2/3, 1/3, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (4/3, 2/3, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC2': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 0, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 0, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AAl': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + sv = np.array(shift_vector) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos + sv)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + # Create AA tilted stacking. + if stacking == 'AAt': + cell = StartingFramework.as_dict()['lattice'] + + # Shift the cell by the tilt angle + a_cell = cell['a'] + b_cell = cell['b'] + c_cell = cell['c'] * 2 + alpha = cell['alpha'] - tilt_angle + beta = cell['beta'] - tilt_angle + gamma = cell['gamma'] + + self.cellMatrix = cellpar_to_cell([a_cell, b_cell, c_cell, alpha, beta, gamma]) + self.cellParameters = np.array([a_cell, b_cell, c_cell, alpha, beta, gamma]).astype(float) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = stacked_structure.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = stacked_structure.formula + + dist_matrix = stacked_structure.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(stacked_structure, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_refined_structure(keep_site_properties=True) + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_hxl_a_structure(self, + BB_H6: str, + BB_L2: str, + stacking: str = 'AA', + print_result: bool = True, + slab: float = 10.0, + shift_vector: list = [1.0, 1.0, 0], + tilt_angle: float = 5.0): + """Creates a COF with HXL-A network. + + The HXL-A net is composed of one hexapodal and one linear building blocks. + + Parameters + ---------- + BB_H6 : BuildingBlock, required + The BuildingBlock object of the tetrapodal Buiding Block A + BB_L2 : BuildingBlock, required + The BuildingBlock object of the tetrapodal Buiding Block B + stacking : str, optional + The stacking pattern of the COF layers (default is 'AA') + print_result : bool, optional + Parameter for the control for printing the result (default is True) + slab : float, optional + Default parameter for the interlayer slab (default is 10.0) + shift_vector: list, optional + Shift vector for the AAl and AAt stakings (defatult is [1.0,1.0,0]) + tilt_angle: float, optional + Tilt angle for the AAt staking in degrees (default is 5.0) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_H6.connectivity != 6: + self.logger.error(connectivity_error.format('A', 6, BB_H6.connectivity)) + raise BBConnectivityError(6, BB_H6.connectivity) + if BB_L2.connectivity != 2: + self.logger.error(connectivity_error.format('B', 3, BB_L2.connectivity)) + raise BBConnectivityError(2, BB_L2.connectivity) + + self.name = f'{BB_H6.name}-{BB_L2.name}-HXL_A-{stacking}' + self.topology = 'HXL_A' + self.staking = stacking + self.dimension = 2 + + self.charge = BB_H6.charge + BB_L2.charge + self.chirality = BB_H6.chirality or BB_L2.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_H6.conector, BB_L2.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_H6.conector, + BB_L2.conector)) + + # Replace "X" the building block + BB_L2.replace_X(bond_atom) + + # Remove the "X" atoms from the the building block + BB_H6.remove_X() + BB_L2.remove_X() + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = BB_H6.size[0] + BB_L2.size[0] + + # Calculate the delta size to add to the c parameter + delta_a = abs(max(np.transpose(BB_H6.atom_pos)[2])) + abs(min(np.transpose(BB_H6.atom_pos)[2])) + delta_b = abs(max(np.transpose(BB_L2.atom_pos)[2])) + abs(min(np.transpose(BB_L2.atom_pos)[2])) + + delta_max = max([delta_a, delta_b]) + + # Calculate the cell parameters + a = topology_info['a'] * size + b = topology_info['b'] * size + c = topology_info['c'] + delta_max + alpha = topology_info['alpha'] + beta = topology_info['beta'] + gamma = topology_info['gamma'] + + if self.stacking == 'A': + c = slab + + # Create the lattice + self.cellMatrix = Lattice.from_parameters(a, b, c, alpha, beta, gamma) + self.cellParameters = np.array([a, b, c, alpha, beta, gamma]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Add the building blocks to the structure + for vertice_data in topology_info['vertices']: + self.atom_types += BB_H6.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', + vertice_data['angle'], + degrees=True).as_matrix() + + rotated_pos = np.dot(BB_H6.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C1' if i == 'C' else i for i in BB_H6.atom_labels] + + # Add the building blocks to the structure + for edge_data in topology_info['edges']: + self.atom_types += BB_L2.atom_types + edge_pos = np.array(edge_data['position'])*a + + R_Matrix = R.from_euler('z', + edge_data['angle'], + degrees=True).as_matrix() + + rotated_pos = np.dot(BB_L2.atom_pos, R_Matrix) + edge_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C2' if i == 'C' else i for i in BB_L2.atom_labels] + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + # Translates the structure to the center of the cell + StartingFramework.translate_sites( + range(len(StartingFramework.as_dict()['sites'])), + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + + if stacking == 'A' or stacking == 'AA': + stacked_structure = StartingFramework + + if stacking == 'AB1': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [2/3, 1/3, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AB2': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/2, 0, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/2, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC1': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (2/3, 1/3, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (4/3, 2/3, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC2': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 0, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 0, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AAl': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + sv = np.array(shift_vector) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos + sv)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + # Create AA tilted stacking. + if stacking == 'AAt': + cell = StartingFramework.as_dict()['lattice'] + + # Shift the cell by the tilt angle + a_cell = cell['a'] + b_cell = cell['b'] + c_cell = cell['c'] * 2 + alpha = cell['alpha'] - tilt_angle + beta = cell['beta'] - tilt_angle + gamma = cell['gamma'] + + self.cellMatrix = cellpar_to_cell([a_cell, b_cell, c_cell, alpha, beta, gamma]) + self.cellParameters = np.array([a_cell, b_cell, c_cell, alpha, beta, gamma]).astype(float) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = stacked_structure.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = stacked_structure.formula + + dist_matrix = stacked_structure.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(stacked_structure, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_refined_structure(keep_site_properties=True) + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_fxt_structure(self, + BB_S4_A: str, + BB_S4_B: str, + stacking: str = 'AA', + print_result: bool = True, + slab: float = 10.0, + shift_vector: list = [1.0, 1.0, 0], + tilt_angle: float = 5.0): + """Creates a COF with FXT network. + + The FXT net is composed of two tetrapodal building blocks. + + Parameters + ---------- + BB_S4_A : BuildingBlock, required + The BuildingBlock object of the tetrapodal Buiding Block A + BB_S4_B : BuildingBlock, required + The BuildingBlock object of the tetrapodal Buiding Block B + stacking : str, optional + The stacking pattern of the COF layers (default is 'AA') + print_result : bool, optional + Parameter for the control for printing the result (default is True) + slab : float, optional + Default parameter for the interlayer slab (default is 10.0) + shift_vector: list, optional + Shift vector for the AAl and AAt stakings (defatult is [1.0,1.0,0]) + tilt_angle: float, optional + Tilt angle for the AAt staking in degrees (default is 5.0) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_S4_A.connectivity != 4: + self.logger.error(connectivity_error.format('A', 4, BB_S4_A.connectivity)) + raise BBConnectivityError(4, BB_S4_A.connectivity) + if BB_S4_B.connectivity != 4: + self.logger.error(connectivity_error.format('B', 4, BB_S4_B.connectivity)) + raise BBConnectivityError(4, BB_S4_B.connectivity) + + self.name = f'{BB_S4_A.name}-{BB_S4_B.name}-FXT-{stacking}' + self.topology = 'FXT' + self.staking = stacking + self.dimension = 2 + + self.charge = BB_S4_A.charge + BB_S4_B.charge + self.chirality = BB_S4_A.chirality or BB_S4_B.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_S4_A.conector, BB_S4_B.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_S4_A.conector, + BB_S4_B.conector)) + + # Replace "X" the building block + BB_S4_A.replace_X(bond_atom) + + # Remove the "X" atoms from the the building block + BB_S4_A.remove_X() + BB_S4_B.remove_X() + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = 2 * (BB_S4_A.size[0] + BB_S4_B.size[0]) + + # Calculate the delta size to add to the c parameter + delta_a = abs(max(np.transpose(BB_S4_A.atom_pos)[2])) + abs(min(np.transpose(BB_S4_B.atom_pos)[2])) + delta_b = abs(max(np.transpose(BB_S4_A.atom_pos)[2])) + abs(min(np.transpose(BB_S4_B.atom_pos)[2])) + + delta_max = max([delta_a, delta_b]) + + # Calculate the cell parameters + a = topology_info['a'] * size + b = topology_info['b'] * size + c = topology_info['c'] + delta_max + alpha = topology_info['alpha'] + beta = topology_info['beta'] + gamma = topology_info['gamma'] + + if self.stacking == 'A': + c = slab + + # Create the lattice + self.cellMatrix = Lattice.from_parameters(a, b, c, alpha, beta, gamma) + self.cellParameters = np.array([a, b, c, alpha, beta, gamma]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Add the first building block to the structure + vertice_data = topology_info['vertices'][0] + self.atom_types += BB_S4_A.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', vertice_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_S4_A.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C1' if i == 'C' else i for i in BB_S4_A.atom_labels] + + # Add the second building block to the structure + for vertice_data in topology_info['vertices'][1:]: + self.atom_types += BB_S4_B.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', vertice_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_S4_B.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C2' if i == 'C' else i for i in BB_S4_B.atom_labels] + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + # Translates the structure to the center of the cell + StartingFramework.translate_sites( + range(len(StartingFramework.as_dict()['sites'])), + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + + if stacking == 'A' or stacking == 'AA': + stacked_structure = StartingFramework + + if stacking == 'AB1': + + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/4, 1/4, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/4, 1/4, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AB2': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/2, 0, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/2, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC1': + + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (1/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 1/3, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 2/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 2/3, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC2': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (1/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 0, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 2/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 0, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + # Create AAl stacking. Tetragonal cell with two sheets + # per cell shifited by the shift_vector in angstroms. + if stacking == 'AAl': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + sv = np.array(shift_vector) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos + sv)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + # Create AA tilted stacking. + # Tilted tetragonal cell with two sheets per cell tilted by tilt_angle. + if stacking == 'AAt': + cell = StartingFramework.as_dict()['lattice'] + + # Shift the cell by the tilt angle + a_cell = cell['a'] + b_cell = cell['b'] + c_cell = cell['c'] * 2 + alpha = cell['alpha'] - tilt_angle + beta = cell['beta'] - tilt_angle + gamma = cell['gamma'] + + self.cellMatrix = cellpar_to_cell([a_cell, b_cell, c_cell, alpha, beta, gamma]) + self.cellParameters = np.array([a_cell, b_cell, c_cell, alpha, beta, gamma]).astype(float) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = stacked_structure.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = stacked_structure.formula + + dist_matrix = stacked_structure.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(stacked_structure, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_refined_structure(keep_site_properties=True) + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_fxt_a_structure(self, + BB_S4: str, + BB_L2: str, + stacking: str = 'AA', + c_parameter_base: float = 3.6, + print_result: bool = True, + slab: float = 10.0, + shift_vector: list = [1.0, 1.0, 0], + tilt_angle: float = 5.0): + """Creates a COF with FXT-A network. + + The FXT-A net is composed of one tetrapodal and one linear building blocks. + + Parameters + ---------- + BB_S4 : BuildingBlock, required + The BuildingBlock object of the tetrapodal Buiding Block + BB_L2 : BuildingBlock, required + The BuildingBlock object of the bipodal Buiding Block + stacking : str, optional + The stacking pattern of the COF layers (default is 'AA') + print_result : bool, optional + Parameter for the control for printing the result (default is True) + slab : float, optional + Default parameter for the interlayer slab (default is 10.0) + shift_vector: list, optional + Shift vector for the AAl and AAt stakings (defatult is [1.0,1.0,0]) + tilt_angle: float, optional + Tilt angle for the AAt staking in degrees (default is 5.0) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_S4.connectivity != 4: + self.logger.error(connectivity_error.format('A', 4, BB_S4.connectivity)) + raise BBConnectivityError(4, BB_S4.connectivity) + if BB_L2.connectivity != 2: + self.logger.error(connectivity_error.format('B', 3, BB_L2.connectivity)) + raise BBConnectivityError(2, BB_L2.connectivity) + + self.name = f'{BB_S4.name}-{BB_L2.name}-FXT_A-{stacking}' + self.topology = 'FXT_A' + self.staking = stacking + self.dimension = 2 + + self.charge = BB_S4.charge + BB_L2.charge + self.chirality = BB_S4.chirality or BB_L2.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_S4.conector, BB_L2.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_S4.conector, + BB_L2.conector)) + + # Replace "X" the building block + BB_L2.replace_X(bond_atom) + + # Remove the "X" atoms from the the building block + BB_S4.remove_X() + BB_L2.remove_X() + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = 2 * (BB_S4.size[0] + BB_L2.size[0]) + + # Calculate the delta size to add to the c parameter + delta_a = abs(max(np.transpose(BB_S4.atom_pos)[2])) + abs(min(np.transpose(BB_S4.atom_pos)[2])) + delta_b = abs(max(np.transpose(BB_L2.atom_pos)[2])) + abs(min(np.transpose(BB_L2.atom_pos)[2])) + + delta_max = max([delta_a, delta_b]) + + # Calculate the cell parameters + a = topology_info['a'] * size + b = topology_info['b'] * size + c = topology_info['c'] + delta_max + alpha = topology_info['alpha'] + beta = topology_info['beta'] + gamma = topology_info['gamma'] + + if self.stacking == 'A': + c = slab + + # Create the lattice + self.cellMatrix = Lattice.from_parameters(a, b, c, alpha, beta, gamma) + self.cellParameters = np.array([a, b, c, alpha, beta, gamma]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Add the building blocks to the structure + for vertice_data in topology_info['vertices']: + self.atom_types += BB_S4.atom_types + vertice_pos = np.array(vertice_data['position'])*a + + R_Matrix = R.from_euler('z', vertice_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_S4.atom_pos, R_Matrix) + vertice_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C1' if i == 'C' else i for i in BB_S4.atom_labels] + + # Add the building blocks to the structure + for edge_data in topology_info['edges']: + self.atom_types += BB_L2.atom_types + edge_pos = np.array(edge_data['position'])*a + + R_Matrix = R.from_euler('z', edge_data['angle'], degrees=True).as_matrix() + + rotated_pos = np.dot(BB_L2.atom_pos, R_Matrix) + edge_pos + self.atom_pos += rotated_pos.tolist() + + self.atom_labels += ['C2' if i == 'C' else i for i in BB_L2.atom_labels] + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + # Translates the structure to the center of the cell + StartingFramework.translate_sites( + range(len(StartingFramework.as_dict()['sites'])), + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + + if stacking == 'A' or stacking == 'AA': + stacked_structure = StartingFramework + + if stacking == 'AB1': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [2/3, 1/3, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AB2': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [1/2, 0, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [1/2, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC1': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (2/3, 1/3, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (4/3, 2/3, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'ABC2': + self.cellMatrix *= (1, 1, 3) + self.cellParameters *= (1, 1, 3, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + _, B_list, C_list = np.split(np.arange(len(self.atom_types)), 3) + + # Translate the second sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + B_list, + (1/3, 0, 1/3), + frac_coords=True, + to_unit_cell=True + ) + + # Translate the third sheet by the vector (2/3, 1/3, 0) to generate the B positions + stacked_structure.translate_sites( + C_list, + (2/3, 0, 2/3), + frac_coords=True, + to_unit_cell=True + ) + + if stacking == 'AAl': + self.cellMatrix *= (1, 1, 2) + self.cellParameters *= (1, 1, 2, 1, 1, 1) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + sv = np.array(shift_vector) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos + sv)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + # Create AA tilted stacking. + if stacking == 'AAt': + cell = StartingFramework.as_dict()['lattice'] + + # Shift the cell by the tilt angle + a_cell = cell['a'] + b_cell = cell['b'] + c_cell = cell['c'] * 2 + alpha = cell['alpha'] - tilt_angle + beta = cell['beta'] - tilt_angle + gamma = cell['gamma'] + + self.cellMatrix = cellpar_to_cell([a_cell, b_cell, c_cell, alpha, beta, gamma]) + self.cellParameters = np.array([a_cell, b_cell, c_cell, alpha, beta, gamma]).astype(float) + + self.atom_types = np.concatenate((self.atom_types, self.atom_types)) + self.atom_pos = np.concatenate((self.atom_pos, self.atom_pos)) + self.atom_labels = np.concatenate((self.atom_labels, self.atom_labels)) + + stacked_structure = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ) + + # Get the index of the atoms in the second sheet + B_list = np.split(np.arange(len(self.atom_types)), 2)[1] + + # Translate the second sheet by the vector [2/3, 1/3, 0.5] to generate the B positions + stacked_structure.translate_sites( + B_list, + [0, 0, 0.5], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = stacked_structure.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = stacked_structure.formula + + dist_matrix = stacked_structure.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(stacked_structure, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_refined_structure(keep_site_properties=True) + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_dia_structure(self, + BB_D41: str, + BB_D42: str, + interp_dg: str = '1', + d_param_base: float = 7.2, + print_result: bool = True, + **kwargs): + """Creates a COF with DIA network. + + The DIA net is composed of two tetrapodal tetrahedical building blocks. + + Parameters + ---------- + BB_D41 : BuildingBlock, required + The BuildingBlock object of the tetrapodal tetrahedical Buiding Block + BB_D42 : BuildingBlock, required + The BuildingBlock object of the tetrapodal tetrahedical Buiding Block + interp_dg : str, optional + The degree of interpenetration of the framework (default is '1') + d_param_base : float, optional + The base value for interlayer distance in angstroms (default is 7.2) + print_result : bool, optional + Parameter for the control for printing the result (default is True) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_D41.connectivity != 4: + self.logger.error(connectivity_error.format('A', 4, BB_D41.connectivity)) + raise BBConnectivityError(4, BB_D41.connectivity) + if BB_D42.connectivity != 4: + self.logger.error(connectivity_error.format('B', 4, BB_D42.connectivity)) + raise BBConnectivityError(4, BB_D42.connectivity) + + self.name = f'{BB_D41.name}-{BB_D42.name}-DIA-{interp_dg}' + self.topology = 'DIA' + self.staking = interp_dg + self.dimension = 3 + + self.charge = BB_D41.charge + BB_D42.charge + self.chirality = BB_D41.chirality or BB_D42.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_D41.conector, BB_D42.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_D41.conector, + BB_D42.conector)) + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = np.average(BB_D41.size) + np.average(BB_D42.size) + + # Calculate the primitive cell vector assuming tetrahedical building blocks + a_prim = np.sqrt(2)*size*np.sqrt((1 - np.cos(1.9106316646041868))) + a_conv = np.sqrt(2)*a_prim + + # Create the primitive lattice + self.cellMatrix = Lattice(a_conv/2*np.array(topology_info['lattice'])) + self.cellParameters = np.array([a_prim, a_prim, a_prim, 60, 60, 60]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Align and rotate the building block 1 to their respective positions + BB_D41.align_to(topology_info['vertices'][0]['align_v']) + + # Determine the angle that alings the X[1] to one of the vertices of the tetrahedron + vertice_pos = unit_vector(np.array([1, 0, 1])) + Q_vertice_pos = BB_D41.get_X_points()[1][1] + + rotated_list = [ + R.from_rotvec( + angle * unit_vector(topology_info['vertices'][0]['align_v']), degrees=False + ).apply(Q_vertice_pos) + for angle in np.linspace(0, 2*np.pi, 360) + ] + + # Calculate the angle between the vertice_pos and the elements of rotated_list + angle_list = [angle(vertice_pos, i) for i in rotated_list] + + rot_angle = np.linspace(0, 360, 360)[np.argmax(angle_list)] + + BB_D41.rotate_around(rotation_axis=np.array(topology_info['vertices'][0]['align_v']), + angle=rot_angle, + degree=True) + + BB_D41.shift(np.array(topology_info['vertices'][0]['position'])*a_conv) + BB_D41.remove_X() + + # Add the building block 1 to the structure + self.atom_types += BB_D41.atom_types + self.atom_pos += BB_D41.atom_pos.tolist() + self.atom_labels += ['C1' if i == 'C' else i for i in BB_D41.atom_labels] + + # Align and rotate the building block 1 to their respective positions + BB_D42.align_to(topology_info['vertices'][0]['align_v']) + + # Determine the angle that alings the X[1] to one of the vertices of the tetrahedron + vertice_pos = unit_vector(np.array([1, 0, 1])) + Q_vertice_pos = BB_D42.get_X_points()[1][1] + + rotated_list = [ + R.from_rotvec( + angle * unit_vector(topology_info['vertices'][0]['align_v']), degrees=False + ).apply(Q_vertice_pos) + for angle in np.linspace(0, 2*np.pi, 360) + ] + + # Calculate the angle between the vertice_pos and the elements of rotated_list + angle_list = [angle(vertice_pos, i) for i in rotated_list] + + rot_angle = np.linspace(0, 360, 360)[np.argmax(angle_list)] + + BB_D42.rotate_around(rotation_axis=np.array(topology_info['vertices'][0]['align_v']), + angle=rot_angle, + degree=True) + + BB_D42.atom_pos = -BB_D42.atom_pos + + BB_D42.shift(np.array(topology_info['vertices'][1]['position'])*a_conv) + + BB_D42.replace_X(bond_atom) + BB_D42.remove_X() + + # Add the building block 2 to the structure + self.atom_types += BB_D42.atom_types + self.atom_pos += BB_D42.atom_pos.tolist() + self.atom_labels += ['C2' if i == 'C' else i for i in BB_D42.atom_labels] + + atom_types, atom_labels, atom_pos = [], [], [] + for n_int in range(int(self.stacking)): + int_direction = np.array([0, 1, 0]) * d_param_base * n_int + + atom_types += self.atom_types + atom_pos += (np.array(self.atom_pos) + int_direction).tolist() + atom_labels += self.atom_labels + + self.atom_types = atom_types + self.atom_pos = atom_pos + self.atom_labels = atom_labels + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + StartingFramework.to(os.path.join(os.getcwd(), 'TESTE_DIA.cif'), fmt='cif') + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = StartingFramework.formula + + dist_matrix = StartingFramework.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(StartingFramework, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_primitive_standard_structure(keep_site_properties=True) + + dict_structure = symm.get_refined_structure(keep_site_properties=True).as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + self.cellParameters = np.array([dict_structure['lattice']['a'], + dict_structure['lattice']['b'], + dict_structure['lattice']['c'], + dict_structure['lattice']['alpha'], + dict_structure['lattice']['beta'], + dict_structure['lattice']['gamma']]).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = self.prim_structure.formula + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_dia_a_structure(self, + BB_D4: str, + BB_L2: str, + interp_dg: str = '1', + d_param_base: float = 7.2, + print_result: bool = True, + **kwargs): + """Creates a COF with DIA-A network. + + The DIA net is composed of two tetrapodal tetrahedical building blocks. + + Parameters + ---------- + BB_D4 : BuildingBlock, required + The BuildingBlock object of the tetrapodal tetrahedical Buiding Block + BB_L2 : BuildingBlock, required + The BuildingBlock object of the dipodal linear Buiding Block + interp_dg : str, optional + The degree of interpenetration of the framework (default is '1') + d_param_base : float, optional + The base value for interlayer distance in angstroms (default is 7.2) + print_result : bool, optional + Parameter for the control for printing the result (default is True) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_D4.connectivity != 4: + self.logger.error(connectivity_error.format('A', 4, BB_D4.connectivity)) + raise BBConnectivityError(4, BB_D4.connectivity) + if BB_L2.connectivity != 2: + self.logger.error(connectivity_error.format('B', 2, BB_L2.connectivity)) + raise BBConnectivityError(2, BB_L2.connectivity) + + self.name = f'{BB_D4.name}-{BB_L2.name}-DIA_A-{interp_dg}' + self.topology = 'DIA_A' + self.staking = interp_dg + self.dimension = 3 + + self.charge = BB_D4.charge + BB_L2.charge + self.chirality = BB_D4.chirality or BB_L2.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_D4.conector, BB_L2.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_D4.conector, + BB_L2.conector)) + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + size = 2 * (np.average(BB_D4.size) + np.average(BB_L2.size)) + + # Calculate the primitive cell vector assuming tetrahedical building blocks + a_prim = np.sqrt(2)*size*np.sqrt((1 - np.cos(1.9106316646041868))) + a_conv = np.sqrt(2)*a_prim + + # Create the primitive lattice + self.cellMatrix = Lattice(a_conv/2*np.array(topology_info['lattice'])) + self.cellParameters = np.array([a_prim, a_prim, a_prim, 60, 60, 60]).astype(float) + + # Create the structure + self.atom_types = [] + self.atom_labels = [] + self.atom_pos = [] + + # Align and rotate the building block 1 to their respective positions + BB_D4.align_to(topology_info['vertices'][0]['align_v']) + + # Determine the angle that alings the X[1] to one of the vertices of the tetrahedron + vertice_pos = unit_vector(np.array([1, 0, 1])) + Q_vertice_pos = BB_D4.get_X_points()[1][1] + + rotated_list = [ + R.from_rotvec( + angle * unit_vector(topology_info['vertices'][0]['align_v']), degrees=False + ).apply(Q_vertice_pos) + for angle in np.linspace(0, 2*np.pi, 360) + ] + + # Calculate the angle between the vertice_pos and the elements of rotated_list + angle_list = [angle(vertice_pos, i) for i in rotated_list] + + rot_angle = np.linspace(0, 360, 360)[np.argmax(angle_list)] + + BB_D4.rotate_around(rotation_axis=np.array(topology_info['vertices'][0]['align_v']), + angle=rot_angle, + degree=True) + + BB_D4.shift(np.array(topology_info['vertices'][0]['position'])*a_conv) + BB_D4.remove_X() + + # Add the building block 1 to the structure + self.atom_types += BB_D4.atom_types + self.atom_pos += BB_D4.atom_pos.tolist() + self.atom_labels += ['C1' if i == 'C' else i for i in BB_D4.atom_labels] + + # Add the building block 1 to the structure + self.atom_types += BB_D4.atom_types + self.atom_pos += list(-np.array(BB_D4.atom_pos) + np.array(topology_info['vertices'][1]['position'])*a_conv) + self.atom_labels += ['C1' if i == 'C' else i for i in BB_D4.atom_labels] + + # Add the building blocks to the structure + for edge_data in topology_info['edges']: + # Copy the building block 2 object + BB = copy.deepcopy(BB_L2) + + # Align, rotate and shift the building block 2 to their respective positions + BB.align_to(edge_data['align_v']) + BB.rotate_around(rotation_axis=edge_data['align_v'], + angle=edge_data['angle']) + BB.shift(np.array(edge_data['position']) * a_conv) + + # Replace "X" the building block with the correct atom dicated by the connection group + BB.replace_X(bond_atom) + BB.remove_X() + + # Update the structure + self.atom_types += BB.atom_types + self.atom_pos += BB.atom_pos.tolist() + self.atom_labels += ['C2' if i == 'C' else i for i in BB.atom_labels] + + atom_types, atom_labels, atom_pos = [], [], [] + for n_int in range(int(self.stacking)): + int_direction = np.array([0, 1, 0]) * d_param_base * n_int + + atom_types += self.atom_types + atom_pos += (np.array(self.atom_pos) + int_direction).tolist() + atom_labels += self.atom_labels + + self.atom_types = atom_types + self.atom_pos = atom_pos + self.atom_labels = atom_labels + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + StartingFramework.translate_sites( + np.ones(len(self.atom_types)).astype(int).tolist(), + [0, 0, 0], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + self.cellParameters = np.array([dict_structure['lattice']['a'], + dict_structure['lattice']['b'], + dict_structure['lattice']['c'], + dict_structure['lattice']['alpha'], + dict_structure['lattice']['beta'], + dict_structure['lattice']['gamma']]).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = StartingFramework.formula + + StartingFramework.to(os.path.join(os.getcwd(), 'TESTE_DIA-A.cif'), fmt='cif') + + dist_matrix = StartingFramework.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(StartingFramework, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_primitive_standard_structure(keep_site_properties=True) + + dict_structure = symm.get_refined_structure(keep_site_properties=True).as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + self.cellParameters = np.array([dict_structure['lattice']['a'], + dict_structure['lattice']['b'], + dict_structure['lattice']['c'], + dict_structure['lattice']['alpha'], + dict_structure['lattice']['beta'], + dict_structure['lattice']['gamma']]).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = self.prim_structure.formula + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] + + def create_bor_structure(self, + BB_D4: str, + BB_T3: str, + interp_dg: str = '1', + d_param_base: float = 7.2, + print_result: bool = True, + **kwargs): + """Creates a COF with BOR network. + + The DIA net is composed of one tetrapodal tetrahedical building block and + one tripodal triangular building block. + + Parameters + ---------- + BB_D4 : BuildingBlock, required + The BuildingBlock object of the tetrapodal tetrahedical Buiding Block + BB_T3 : BuildingBlock, required + The BuildingBlock object of the tripodal triangular Buiding Block + interp_dg : str, optional + The degree of interpenetration of the framework (default is '1') + d_param_base : float, optional + The base value for interlayer distance in angstroms (default is 7.2) + print_result : bool, optional + Parameter for the control for printing the result (default is True) + + Returns + ------- + list + A list of strings containing: + 1. the structure name, + 2. lattice type, + 3. hall symbol of the cristaline structure, + 4. space group, + 5. number of the space group, + 6. number of operation symmetry + """ + connectivity_error = 'Building block {} must present connectivity {} not {}' + if BB_D4.connectivity != 4: + self.logger.error(connectivity_error.format('A', 4, BB_D4.connectivity)) + raise BBConnectivityError(4, BB_D4.connectivity) + if BB_T3.connectivity != 3: + self.logger.error(connectivity_error.format('B', 3, BB_T3.connectivity)) + raise BBConnectivityError(3, BB_T3.connectivity) + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + self.name = f'{BB_D4.name}-{BB_T3.name}-BOR-{interp_dg}' + self.topology = 'BOR' + self.staking = interp_dg + self.dimension = 3 + + self.charge = BB_D4.charge + BB_T3.charge + self.chirality = BB_D4.chirality or BB_T3.chirality + + self.logger.debug(f'Starting the creation of {self.name}') + + # Detect the bond atom from the connection groups type + bond_atom = get_bond_atom(BB_D4.conector, BB_T3.conector) + + self.logger.debug('{} detected as bond atom for groups {} and {}'.format(bond_atom, + BB_D4.conector, + BB_T3.conector)) + + # Get the topology information + topology_info = TOPOLOGY_DICT[self.topology] + + # Measure the base size of the building blocks + d_size = (np.array(BB_D4.size).mean() + np.array(BB_T3.size).mean()) + + # Calculate the primitive cell vector assuming tetrahedical building blocks + a_conv = np.sqrt(6) * d_size + + # Create the primitive lattice + self.cellMatrix = Lattice(a_conv * np.array(topology_info['lattice'])) + self.cellParameters = np.array([a_conv, a_conv, a_conv, 90, 90, 90]).astype(float) + + # Create the structure + atom_types = [] + atom_labels = [] + atom_pos = [] + + for D_site in topology_info['vertices']: + D4 = BB_D4.copy() + D4.align_to( + np.array(D_site['align_v']) + ) + + D4.rotate_around( + rotation_axis=D_site['align_v'], + angle=D_site['angle']) + + D4.shift(np.array(D_site['position'])*a_conv) + + atom_types += D4.atom_types + atom_pos += D4.atom_pos.tolist() + atom_labels += D4.atom_labels.tolist() + + # Translate all atoms to inside the cell + for i, pos in enumerate(atom_pos): + for j, coord in enumerate(pos): + if coord < 0: + atom_pos[i][j] += a_conv + + X_pos = [atom_pos[i] for i in np.where(np.array(atom_types) == 'X')[0]] + + T_site = topology_info['edges'][0] + + _, X = BB_T3.get_X_points() + BB_T3.rotate_around([0, 0, 1], T_site['angle'], True) + + R_matrix = rotation_matrix_from_vectors([0, 0, 1], + T_site['align_v']) + + BB_T3.atom_pos = np.dot(BB_T3.atom_pos, R_matrix.T) + + BB_T3.replace_X('O') + + # Get the 3 atoms that are closer to T_site['position'])*a_conv + X_pos_temp = sorted(X_pos, key=lambda x: np.linalg.norm(x - np.array(T_site['position'])*a_conv)) + + X_center = np.array(X_pos_temp[:3]).mean(axis=0) + + BB_T3.shift(X_center) + + atom_types += BB_T3.atom_types + atom_pos += BB_T3.atom_pos.tolist() + atom_labels += BB_T3.atom_labels.tolist() + + T4 = BB_T3.copy() + T4.rotate_around([0, 0, 1], 180, True) + + atom_types += T4.atom_types + atom_pos += T4.atom_pos.tolist() + atom_labels += T4.atom_labels.tolist() + + T2 = BB_T3.copy() + T2.rotate_around([0, 0, 1], 90, True) + T2.rotate_around([1, 0, 0], -90, True) + + atom_types += T2.atom_types + atom_pos += T2.atom_pos.tolist() + atom_labels += T2.atom_labels.tolist() + + T3 = BB_T3.copy() + T3.rotate_around([0, 0, 1], -90, True) + + T3.atom_pos *= np.array([1, 1, -1]) + + atom_types += T3.atom_types + atom_pos += T3.atom_pos.tolist() + atom_labels += T3.atom_labels.tolist() + + # Translate all atoms to inside the cell + for i, pos in enumerate(atom_pos): + for j, coord in enumerate(pos): + if coord < 0: + atom_pos[i][j] += a_conv + + # Remove the X atoms from the list + X_index = np.where(np.array(atom_types) == 'X')[0] + + self.atom_types = [atom_types[i] for i in range(len(atom_types)) if i not in X_index] + self.atom_pos = [atom_pos[i] for i in range(len(atom_pos)) if i not in X_index] + self.atom_labels = [atom_labels[i] for i in range(len(atom_labels)) if i not in X_index] + + atom_types, atom_labels, atom_pos = [], [], [] + for n_int in range(int(self.stacking)): + int_direction = np.array([0, 1, 0]) * d_param_base * n_int + + atom_types += self.atom_types + atom_pos += (np.array(self.atom_pos) + int_direction).tolist() + atom_labels += self.atom_labels + + self.atom_types = atom_types + self.atom_pos = atom_pos + self.atom_labels = atom_labels + + StartingFramework = Structure( + self.cellMatrix, + self.atom_types, + self.atom_pos, + coords_are_cartesian=True, + site_properties={'source': self.atom_labels} + ).get_sorted_structure() + + StartingFramework.translate_sites( + np.ones(len(self.atom_types)).astype(int).tolist(), + [0, 0, 0], + frac_coords=True, + to_unit_cell=True + ) + + dict_structure = StartingFramework.as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + self.cellParameters = np.array([dict_structure['lattice']['a'], + dict_structure['lattice']['b'], + dict_structure['lattice']['c'], + dict_structure['lattice']['alpha'], + dict_structure['lattice']['beta'], + dict_structure['lattice']['gamma']]).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = StartingFramework.formula + + StartingFramework.to('TESTE_BOR.cif', fmt='cif') + + dist_matrix = StartingFramework.distance_matrix + + # Check if there are any atoms closer than 0.8 A + for i in range(len(dist_matrix)): + for j in range(i+1, len(dist_matrix)): + if dist_matrix[i][j] < self.dist_threshold: + raise BondLenghError(i, j, dist_matrix[i][j], self.dist_threshold) + + # Get the simmetry information of the generated structure + symm = SpacegroupAnalyzer(StartingFramework, + symprec=self.symm_tol, + angle_tolerance=self.angle_tol) + + try: + self.prim_structure = symm.get_primitive_standard_structure(keep_site_properties=True) + + dict_structure = symm.get_refined_structure(keep_site_properties=True).as_dict() + + self.cellMatrix = np.array(dict_structure['lattice']['matrix']).astype(float) + self.cellParameters = np.array([dict_structure['lattice']['a'], + dict_structure['lattice']['b'], + dict_structure['lattice']['c'], + dict_structure['lattice']['alpha'], + dict_structure['lattice']['beta'], + dict_structure['lattice']['gamma']]).astype(float) + + self.atom_types = [i['label'] for i in dict_structure['sites']] + self.atom_pos = [i['xyz'] for i in dict_structure['sites']] + self.atom_labels = [i['properties']['source'] for i in dict_structure['sites']] + self.n_atoms = len(dict_structure['sites']) + self.composition = self.prim_structure.formula + + self.logger.debug(self.prim_structure) + + self.lattice_type = symm.get_lattice_type() + self.space_group = symm.get_space_group_symbol() + self.space_group_n = symm.get_space_group_number() + + symm_op = symm.get_point_group_operations() + self.hall = symm.get_hall() + + except Exception as e: + self.logger.exception(e) + + self.lattice_type = 'Triclinic' + self.space_group = 'P1' + self.space_group_n = '1' + + symm_op = [1] + self.hall = 'P 1' + + symm_text = get_framework_symm_text(self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)) + + self.logger.info(symm_text) + + return [self.name, + str(self.lattice_type), + str(self.hall[0:2]), + str(self.space_group), + str(self.space_group_n), + len(symm_op)] diff --git a/build/lib/io_tools.py b/build/lib/io_tools.py new file mode 100644 index 00000000..e0d79d2c --- /dev/null +++ b/build/lib/io_tools.py @@ -0,0 +1,1287 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +""" +This module contains tools for input and output file manipulation used by pyCOFBuilder. +""" + +import os +from datetime import date +import numpy as np + +from pymatgen.io.cif import CifParser + +import simplejson +from pycofbuilder.tools import (elements_dict, + cell_to_cellpar, + cellpar_to_cell, + get_fractional_to_cartesian_matrix, + get_cartesian_to_fractional_matrix, + get_kgrid, + formula_from_atom_list, + smiles_to_xsmiles, + cell_to_ibrav) + + +def save_csv(path, file_name, data, delimiter=',', head=False): + """ + Saves a file in format `.csv`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the `csv` file. Does not neet to contain the `.csv` extention. + data : list + Data to be saved. + delimiter: str + Delimiter of the columns. `,` is the default. + head : str + Names of the columns. + """ + + # Remove the extention if exists + file_name = file_name.split('.')[0] + file_name = os.path.join(path, file_name + '.csv') + + file_temp = open(file_name, 'w') + if head is not False: + file_temp.write(head) + + for i in range(len(data)): + file_temp.write(delimiter.join([str(j) for j in data[i]]) + '\n') + + file_temp.close() + + +def read_xyz(path, file_name): + """ + Reads a file in format `.xyz` from the `path` given and returns + a list containg the N atom labels and a Nx3 array contaning + the atoms coordinates. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the `xyz` file. Does not neet to contain the `.xyz` extention. + + Returns + ------- + atom_labels : list + List of strings containing containg the N atom labels. + atom_pos : numpy array + Nx3 array contaning the atoms coordinates + """ + + # Remove the extention if exists + file_name = file_name.split('.')[0] + + if os.path.exists(os.path.join(path, file_name + '.xyz')): + temp_file = open(os.path.join(path, file_name + '.xyz'), 'r').readlines() + + atoms = [i.split() for i in temp_file[2:]] + + atom_labels = [i[0] for i in atoms if len(i) > 1] + atom_pos = np.array([[float(i[1]), float(i[2]), float(i[3])] for i in atoms if len(i) > 1]) + + return atom_labels, atom_pos + else: + print(f'File {file_name} not found!') + return None + + +def read_pdb(path, file_name): + """ + Reads a file in format `.pdb` from the `path` given and returns + a list containg the N atom labels and a Nx3 array contaning + the atoms coordinates. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the `pdb` file. Does not neet to contain the `.pdb` extention. + + Returns + ------- + atom_labels : list + List of strings containing containg the N atom labels. + atom_pos : numpy array + Nx3 array contaning the atoms coordinates + """ + + # Remove the extention if exists + file_name = file_name.split('.')[0] + + if not os.path.exists(os.path.join(path, file_name + '.pdb')): + raise FileNotFoundError(f'File {file_name} not found!') + + temp_file = open(os.path.join(path, file_name + '.pdb'), 'r').read().splitlines() + + cellParameters = np.array([i.split()[1:] for i in temp_file if 'CRYST1' in i][0]).astype(float) + + AtomTypes = [i.split()[2] for i in temp_file if 'ATOM' in i] + CartPos = np.array([i.split()[4:7] for i in temp_file if 'ATOM' in i]).astype(float) + + return cellParameters, AtomTypes, CartPos + + +def read_gjf(path, file_name): + """ + Reads a file in format `.gjf` from the `path` given and returns + a list containg the N atom labels and a Nx3 array contaning + the atoms coordinates. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the `gjf` file. Does not neet to contain the `.gjf` extention. + + Returns + ------- + atom_labels : list + List of strings containing containg the N atom labels. + atom_pos : numpy array + Nx3 array contaning the atoms coordinates + """ + # Remove the extention if exists + file_name = file_name.split('.')[0] + + if os.path.exists(os.path.join(path, file_name + '.gjf')): + + temp_file = open(os.path.join(path, file_name + '.gjf'), 'r').readlines() + temp_file = [i.split() for i in temp_file if i != '\n'] + + atoms = [i for i in temp_file if i[0] in elements_dict()] + + atom_labels = [i[0] for i in atoms] + atom_pos = np.array([[float(i[1]), float(i[2]), float(i[3])] for i in atoms]) + + return atom_labels, atom_pos + else: + print(f'File {file_name} not found!') + return None + + +def read_cif(path, file_name): + """ + Reads a file in format `.cif` from the `path` given and returns + a list containg the N atom labels and a Nx3 array contaning + the atoms coordinates. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the `cif` file. Does not neet to contain the `.cif` extention. + + Returns + ------- + cell : numpy array + 3x3 array contaning the cell vectors. + atom_labels : list + List of strings containing containg the N atom labels. + atom_pos : numpy array + Nx3 array contaning the atoms coordinates + charges : list + List of strings containing containg the N atom partial charges. + """ + + # Remove the extention if exists + file_name = file_name.split('.')[0] + + if os.path.exists(os.path.join(path, file_name + '.cif')): + + temp_file = open(os.path.join(path, file_name + '.cif'), 'r').readlines() + cell = [] + atom_label = [] + atom_pos = [] + charges = [] + has_charges = False + for i in temp_file: + if 'cell_length_a' in i: + cell += [float(i.split()[-1])] + if 'cell_length_b' in i: + cell += [float(i.split()[-1])] + if 'cell_length_c' in i: + cell += [float(i.split()[-1])] + if 'cell_angle_alpha' in i: + cell += [float(i.split()[-1])] + if '_cell_angle_beta' in i: + cell += [float(i.split()[-1])] + if '_cell_angle_gamma' in i: + cell += [float(i.split()[-1])] + if '_atom_site_charge' in i: + has_charges = True + + for i in temp_file: + line = i.split() + if len(line) > 1 and line[0] in elements_dict().keys(): + atom_label += [line[0]] + atom_pos += [[float(j) for j in line[2:5]]] + if has_charges: + charges += [float(line[-1])] + cell = cellpar_to_cell(cell) + + return cell, atom_label, atom_pos, charges + else: + print(f'File {file_name} not found!') + return None + + +def save_xsf(path: str = None, + file_name: str = None, + cell: np.ndarray = np.eye(3), + atom_types: list = None, + atom_pos: list = None, + atom_labels: list = None, + frac_coords: bool = False): + """ + Save a file in format `.xsf` on the `path`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.xsf` extention. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + atom_label : list + List of strings containing containg the N atom label. + atom_pos : list + Nx3 array contaning the atoms coordinates. + """ + + file_name = file_name.split('.')[0] + + if len(cell) == 6: + cell = cellpar_to_cell(cell) + + if frac_coords: + # Convert to fractional coordinates + frac_matrix = get_fractional_to_cartesian_matrix(*cell_to_cellpar(cell)) + atom_pos = [np.dot(frac_matrix, [i[0], i[1], i[2]]) for i in atom_pos] + + xsf_file = open(os.path.join(path, file_name + '.xsf'), 'w') + xsf_file.write(' CRYSTAL\n') + xsf_file.write(' PRIMVEC\n') + + for i in range(len(cell)): + xsf_file.write(f' {cell[i][0]:>15.9f} {cell[i][1]:>15.9f} {cell[i][2]:>15.9f}\n') + + xsf_file.write(' PRIMCOORD\n') + xsf_file.write(f' {len(atom_pos)} 1\n') + + for i in range(len(atom_pos)): + xsf_file.write('{:3s} {:>15.9f} {:>15.9f} {:>15.9f}\n'.format(atom_types[i], + atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2])) + + xsf_file.close() + + +def save_pqr(path: str = None, + file_name: str = None, + cell: np.ndarray = np.eye(3), + atom_types: list = None, + atom_pos: list = None, + atom_labels: list = None, + partial_charges: list = None, + frac_coords: bool = False): + """ + Save a file in format `.pqr` on the `path`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.pqr` extention. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + atom_label : list + List of strings containing containg the N atom partial charges. + atom_pos : list + Nx3 array contaning the atoms coordinates. + partial_charges: list + List of strings containing containg the N atom partial charges. + """ + + file_name = file_name.split('.')[0] + + if len(cell) == 3: + cell = cell_to_cellpar(cell) + + if frac_coords: + # Convert to fractional coordinates + frac_matrix = get_fractional_to_cartesian_matrix(*cell) + atom_pos = [np.dot(frac_matrix, [i[0], i[1], i[2]]) for i in atom_pos] + + pqr_file = open(os.path.join(path, file_name + '.pqr'), 'w') + pqr_file.write(f'TITLE {file_name} \n') + pqr_file.write('REMARK 4\n') + pqr_file.write('CRYST1{:>9.3f}{:>9.3f}{:>9.3f}{:>7.2f}{:>7.2f}{:>7.2f} P1\n'.format(cell[0], + cell[1], + cell[2], + cell[3], + cell[4], + cell[5])) + + if partial_charges is None: + atom_line = 'ATOM {:>4} {:>2} MOL A 0 {:>8.3f}{:>8.3f}{:>8.3f} {:>15}\n' + for i in range(len(atom_pos)): + pqr_file.write(atom_line.format(i + 1, + atom_types[i], + atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2], + atom_types[i])) + else: + atom_line = 'ATOM {:>4} {:>2} MOL A 0 {:>8.3f}{:>8.3f}{:>8.3f}{:>8.5f} {:>15}\n' + for i in range(len(atom_pos)): + pqr_file.write(atom_line.format(i + 1, + atom_types[i], + atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2], + partial_charges[i], + atom_types[i])) + + pqr_file.close() + + +def save_pdb(path: str = None, + file_name: str = None, + cell: np.ndarray = np.eye(3), + atom_types: list = None, + atom_pos: list = None, + atom_labels: list = None, + frac_coords: bool = False): + """ + Save a file in format `.pdb` on the `path`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.pdb` extention. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + atom_label : list + List of strings containing containg the N atom partial charges. + atom_pos : list + Nx3 array contaning the atoms coordinates in cartesian form. + """ + + file_name = file_name.split('.')[0] + + if len(cell) == 3: + cell = cell_to_cellpar(cell) + + if frac_coords: + # Convert to fractional coordinates + frac_matrix = get_fractional_to_cartesian_matrix(*cell) + atom_pos = [np.dot(frac_matrix, [i[0], i[1], i[2]]) for i in atom_pos] + + pdb_file = open(os.path.join(path, file_name + '.pdb'), 'w') + pdb_file.write(f'TITLE {file_name} \n') + pdb_file.write('REMARK pyCOFBuilder\n') + pdb_file.write('CRYST1{:>9.3f}{:>9.3f}{:>9.3f}{:>7.2f}{:>7.2f}{:>7.2f} P1\n'.format(cell[0], + cell[1], + cell[2], + cell[3], + cell[4], + cell[5])) + + atom_line = 'ATOM {:>4} {:>2} MOL {:>13.3f}{:>8.3f}{:>8.3f} 1.00 0.00 {:>11}\n' + + for i in range(len(atom_pos)): + pdb_file.write(atom_line.format(i+1, + atom_types[i], + atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2], + atom_types[i])) + + pdb_file.close() + + +def save_gjf(path: str = None, + file_name: str = None, + cell: list = None, + atom_types: list = None, + atom_pos: list = None, + atom_labels: list = None, + frac_coords: bool = False, + text: str = 'opt pm6'): + + """ + Save a file in format `.gjf` on the `path`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.gjf` extention. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + atom_label : list + List of strings containing containg the N atom partial charges. + atom_pos : list + Nx3 array contaning the atoms coordinates. + text : str + Parameters for Gaussian calculations. + """ + + if len(cell) == 6: + cell = cellpar_to_cell(cell) + + if frac_coords: + # Convert to fractional coordinates + frac_matrix = get_fractional_to_cartesian_matrix(*cell_to_cellpar(cell)) + atom_pos = [np.dot(frac_matrix, [i[0], i[1], i[2]]) for i in atom_pos] + + file_name = file_name.split('.')[0] + + temp_file = open(os.path.join(path, file_name + '.gjf'), 'w') + temp_file.write(f'%chk={file_name}.chk \n') + temp_file.write(f'# {text}\n') + temp_file.write('\n') + temp_file.write(f'{file_name}\n') + temp_file.write('\n') + temp_file.write('0 1 \n') + + for i in range(len(atom_types)): + temp_file.write('{:<5s}{:>15.7f}{:>15.7f}{:>15.7f}\n'.format(atom_types[i], + atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2])) + if cell is not None: + for i in range(len(cell)): + temp_file.write('Tv {:>15.7f}{:>15.7f}{:>15.7f}\n'.format(*cell[i])) + + temp_file.write('\n\n') + temp_file.close() + + +def save_xyz(path: str = None, + file_name: str = None, + cell: np.ndarray = np.eye(3), + atom_types: list = None, + atom_pos: list = None, + atom_labels: list = None, + frac_coords: bool = False): + """ + Save a file in format `.xyz` on the `path`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.xyz` extention. + atom_types : list + List of strings containing containg the N atom types + atom_pos : list + Nx3 array contaning the atoms coordinates. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + """ + + if len(cell) == 3: + cell = cell_to_cellpar(cell) + + if frac_coords: + # Convert to fractional coordinates + frac_matrix = get_fractional_to_cartesian_matrix(cell) + atom_pos = [np.dot(frac_matrix, [i[0], i[1], i[2]]) for i in atom_pos] + + file_name = file_name.split('.')[0] + + temp_file = open(os.path.join(path, file_name + '.xyz'), 'w') + temp_file.write(f'{len(atom_types)}\n') + + if cell is None: + temp_file.write(f'{file_name}\n') + else: + temp_file.write(f'{cell[0]} {cell[1]} {cell[2]} {cell[3]} {cell[4]} {cell[5]}\n') + + for i in range(len(atom_types)): + temp_file.write('{:<5s}{:>15.7f}{:>15.7f}{:>15.7f}\n'.format(atom_types[i], + atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2])) + + temp_file.close() + + +def save_turbomole(path: str = None, + file_name: str = None, + cell: np.ndarray = np.eye(3), + atom_types: list = None, + atom_pos: list = None, + atom_labels: list = None, + frac_coords: bool = False): + """Save the structure in Turbomole .coord format + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.coord` extention. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + atom_label : list + List of strings containing containg the N atom partial charges. + atom_pos : list + Nx3 array contaning the atoms coordinates. + """ + + if cell.shape == (3, 3): + cell = cell_to_cellpar(cell) + + if frac_coords: + # Convert to fractional coordinates + frac_matrix = get_fractional_to_cartesian_matrix(*cell) + atom_pos = [np.dot(frac_matrix, [i[0], i[1], i[2]]) for i in atom_pos] + + with open(os.path.join(path, file_name + '.coord'), 'w') as temp_file: + temp_file.write('$coord angs\n') + + for i in range(len(atom_types)): + temp_file.write('{:>15.7f}{:>15.7f}{:>15.7f} {:<5s}\n'.format(atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2], + atom_types[i])) + + temp_file.write('$periodic 3\n') + temp_file.write('$cell\n') + temp_file.write('{} {} {} {} {} {}\n'.format(*cell)) + temp_file.write('$opt\n') + temp_file.write(' engine=inertial\n') + temp_file.write('$end\n') + + +def save_vasp(path: str = None, + file_name: str = None, + cell: np.ndarray = np.eye(3), + atom_types: list = None, + atom_pos: list = None, + atom_labels: list = None, + frac_coords: bool = False): + """Save the structure in VASP .vasp format + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.vasp` extention. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + atom_label : list + List of strings containing containg the N atom partial charges. + atom_pos : list + Nx3 array contaning the atoms coordinates. + coords_are_cartesian : bool + If True, the coordinates are in cartesian coordinates. + """ + + if cell.shape == 6: + cell = cellpar_to_cell(cell) + + unique_atoms = [] + for i in atom_types: + if i not in unique_atoms: + unique_atoms.append(i) + + composition_dict = {i: atom_types.count(i) for i in unique_atoms} + + with open(os.path.join(path, file_name + '.vasp'), 'w') as temp_file: + temp_file.write(f'{file_name}\n') + temp_file.write('1.0\n') + + for i in range(3): + temp_file.write('{:>15.7f}{:>15.7f}{:>15.7f}\n'.format(cell[i][0], + cell[i][1], + cell[i][2])) + + temp_file.write(' '.join(composition_dict.keys()) + '\n') + temp_file.write(' '.join([str(i) for i in composition_dict.values()]) + '\n') + + if frac_coords: + temp_file.write('Direct\n') + else: + temp_file.write('Cartesian\n') + + for i in range(len(atom_types)): + temp_file.write('{:>15.7f}{:>15.7f}{:>15.7f} {:<5s}\n'.format(atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2], + atom_types[i])) + + +def save_qe(path: str = None, + file_name: str = None, + cell: np.ndarray = np.eye(3), + atom_types: list = None, + atom_pos: list = None, + atom_labels: list = None, + frac_coords: bool = False, + calc_type: str = 'scf'): + ''' + Save the structure in Quantum Espresso .pwscf format. + + The `input_dict` can be used to specify the input parameters for the + QuantumESPRESSO calculation. + + This dictionary must contain the keys: `control`, `system`, `electrons`, and `ions`. + Each of these keys must contain a dictionary with the corresponding input parameters. + This dictionary can contain the kpoints item, with the kpoints grid as a list of 3 integers. + Additionally, it can contain the kspacing item, with the kpoints spacing as a float. In this + case the kpoints grid will be calculated automatically. By default, the kspacing is set to 0.3. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.pwscf` extention. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + atom_label : list + List of strings containing containg the N atom labels. + atom_pos : list + Nx3 array contaning the atoms coordinates. + frac_coords : bool + If True, the coordinates are in fractional coordinates. + calc_type : str + Type of pw.x calculation. Can be 'scf', 'nscf', 'bands', and 'vc-relax'. + ''' + + if len(cell) == 6: + cell_matrix = cellpar_to_cell(cell) + else: + cell_matrix = cell + + ibrav_dict = cell_to_ibrav(cell_matrix) + + input_dict = {} + + input_dict['control'] = { + 'prefix': f"'{file_name}'", + 'calculation': f"{calc_type}", + 'restart_mode': "'from_scratch'", + 'wf_collect': '.true.', + 'pseudo_dir': "'$PSEUDO_DIR'", + 'outdir': "'$SCRATCH_DIR'", + 'verbosity': "'high'", + 'tstress': '.true.', + 'tprnfor': '.true.', + 'etot_conv_thr': '1.0d-5', + 'forc_conv_thr': '1.0d-6', + 'nstep': 1000} + + input_dict['system'] = { + 'nat': len(atom_types), + 'ntyp': len(set(atom_types)), + 'ecutwfc': 40, + 'ecutrho': 360, + 'vdw_corr': "'grimme-d3'", + 'occupations': "'smearing'", + **ibrav_dict} + + input_dict['electrons'] = { + 'conv_thr': 1.0e-9, + 'electron_maxstep': 100, + 'mixing_beta': 0.3} + + if calc_type == 'vc-relax': + + input_dict['ions'] = { + 'ion_dynamics': "'bfgs'"} + + input_dict['cell'] = { + 'cell_dynamics': "'bfgs'", + 'cell_dofree': "'all'"} + + # If the kpoints grid is not specified, calculate it automatically + if 'k_points' not in input_dict.keys(): + if 'kspacing' not in input_dict.keys(): + input_dict['kspacing'] = 0.3 + input_dict['kpoints'] = get_kgrid(cell_matrix, input_dict['kspacing']) + + with open(os.path.join(path, file_name + '.pwscf'), 'w') as f: + f.write('&CONTROL\n') + for key in input_dict['control']: + f.write(f" {key} = {input_dict['control'][key]}\n") + f.write('/\n\n') + + f.write('&SYSTEM\n') + for key in input_dict['system']: + f.write(f" {key} = {input_dict['system'][key]}\n") + f.write('/\n\n') + + f.write('&ELECTRONS\n') + for key in input_dict['electrons']: + f.write(f" {key} = {input_dict['electrons'][key]}\n") + f.write('/\n\n') + + if calc_type == 'vc-relax': + f.write('&IONS\n') + for key in input_dict['ions']: + f.write(f" {key} = {input_dict['ions'][key]}\n") + f.write('/\n\n') + + f.write('&CELL\n') + for key in input_dict['cell']: + f.write(f" {key} = {input_dict['cell'][key]}\n") + f.write('/\n\n') + + f.write('ATOMIC_SPECIES\n') + for atom in set(atom_types): + f.write(f" {atom} {elements_dict()[atom]:>9.5f} {atom}.PSEUDO.UPF\n") + f.write('\n') + + # f.write('CELL_PARAMETERS (angstrom) \n') + # for v in cell_matrix: + # f.write(f'{v[0]:>15.9f} {v[1]:>15.9f} {v[2]:>15.9f}\n') + # f.write('\n') + + if frac_coords: + coords_type = 'crystal' + else: + coords_type = 'angstrom' + + f.write(f'ATOMIC_POSITIONS ({coords_type})\n') + + for i, atom in enumerate(atom_pos): + f.write('{:<5s}{:>15.9f}{:>15.9f}{:>15.9f} ! {:5}\n'.format(atom_types[i], + atom[0], + atom[1], + atom[2], + atom_labels[i])) + + f.write('\n') + f.write('K_POINTS automatic\n') + f.write(' {} {} {} 1 1 1\n'.format(*input_dict['kpoints'])) + + +def convert_cif_2_qe(out_path, file_name): + """ + Convert a cif file to a Quantum Espresso input file + + Parameters + ---------- + out_path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.cif` extention. + """ + + cell, atom_labels, atom_pos, _ = read_cif(out_path, file_name, has_charges=False) + + print(cell, atom_labels, atom_pos) + + save_qe(out_path, + file_name, + cell, + atom_labels, + atom_pos, + coords_are_cartesian=True, + supercell=False, + angs=False, + ecut=40, + erho=360, + k_dist=0.3) + + +def save_json(path, file_name, cell, atom_types, atom_pos, atom_labels, frac_coords=False): + """ + Save a file in format `.json` on the `path`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.cif` extention. + atom_label : list + List of strings containing containg the N atom partial charges. + atom_pos : list + Nx3 array contaning the atoms coordinates. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + """ + + file_name = file_name.split('.')[0] + + cof_json = create_COF_json(file_name) + + if len(cell) == 3: + cell_par = cell_to_cellpar(np.array(cell)).tolist() + cell_matrix = np.array(cell).astype(float).tolist() + + if len(cell) == 6: + cell_par = np.array(cell).astype(float).tolist() + cell_matrix = cellpar_to_cell(cell_par).tolist() + + cof_json['system']['geo_opt'] = False + + cof_json['geometry']['cell_matrix'] = cell_matrix + cof_json['geometry']['cell_parameters'] = cell_par + cof_json['geometry']['atom_labels'] = list(atom_types) + cof_json['geometry']['atom_pos'] = list(atom_pos) + + write_json(path, file_name, cof_json) + + +def save_chemjson(path, + file_name, + cell, + atom_types, + atom_pos, + atom_labels, + frac_coords=False): + """ + Save a file in format `.json` on the `path`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.cif` extention. + atom_label : list + List of strings containing containg the N atom partial charges. + atom_pos : list + Nx3 array contaning the atoms coordinates. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + """ + + file_name = file_name.split('.')[0] + if len(cell) == 6: + CellParameters = cell + CellMatrix = None + if len(cell) == 3: + CellParameters = None + CellMatrix = cell + + chemJSON = create_structure_CJSON(StructureName=file_name.split('.')[0], + CellParameters=CellParameters, + CellMatrix=CellMatrix, + AtomTypes=atom_types, + AtomPositions=atom_pos, + AtomLabels=atom_labels, + CartesianPositions=not frac_coords) + + write_json(path, file_name, chemJSON) + + +def save_cif(path, + file_name, + cell, + atom_types, + atom_pos, + atom_labels=None, + partial_charges=False, + frac_coords=False): + """ + Save a file in format `.cif` on the `path`. + + Parameters + ---------- + path : str + Path to the file. + file_name : str + Name of the file. Does not neet to contain the `.cif` extention. + atom_label : list + List of strings containing containg the N atom partial charges. + atom_pos : list + Nx3 array contaning the atoms coordinates. + cell : numpy array + Can be a 3x3 array contaning the cell vectors or a list with the 6 cell parameters. + """ + + file_name = file_name.split('.')[0] + + if len(cell) == 3: + a, b, c, alpha, beta, gamma = cell_to_cellpar(cell) + if len(cell) == 6: + a, b, c, alpha, beta, gamma = cell + + if atom_labels is None: + atom_labels = [''] * len(atom_types) + + cif_text = f"""\ +data_{file_name} + +_audit_creation_date {date.today().strftime("%Y-%d-%m")} +_audit_creation_method pyCOFBuilder +_audit_author_name 'Felipe Lopes de Oliveira' + +_chemical_name_common '{file_name}' +_cell_length_a {a:>10.6f} +_cell_length_b {b:>10.6f} +_cell_length_c {c:>10.6f} +_cell_angle_alpha {alpha:>6.2f} +_cell_angle_beta {beta:>6.2f} +_cell_angle_gamma {gamma:>6.2f} +_space_group_name_H-M_alt 'P 1' +_space_group_IT_number 1 + +loop_ +_symmetry_equiv_pos_as_xyz + 'x, y, z' + +loop_ + _atom_site_label + _atom_site_type_symbol + _atom_site_fract_x + _atom_site_fract_y + _atom_site_fract_z +""" + + if partial_charges is not False: + cif_text += ' _atom_site_charge\n' + + if frac_coords is False: + # Convert to fractional coordinates + frac_matrix = get_cartesian_to_fractional_matrix(a, b, c, alpha, beta, gamma) + atom_pos = [np.dot(frac_matrix, [i[0], i[1], i[2]]) for i in atom_pos] + + for i in range(len(atom_pos)): + u, v, w = atom_pos[i][0], atom_pos[i][1], atom_pos[i][2] + if partial_charges is not False: + cif_text += '{:<7} {} {:>15.9f} {:>15.9f} {:>15.9f} {:>10.5f}\n'.format( + f"{atom_types[i]}{str(i + 1)}_{atom_labels[i]}", + atom_types[i], + u, + v, + w, + partial_charges[i]) + else: + cif_text += '{:<7} {} {:>15.9f} {:>15.9f} {:>15.9f}\n'.format( + f"{atom_types[i]}{str(i + 1)}_{atom_labels[i]}", + atom_types[i], + u, + v, + w) + + # Write cif_text to file + cif_file = open(os.path.join(path, file_name + '.cif'), 'w') + cif_file.write(cif_text) + cif_file.close() + + +def convert_json_2_cif(origin_path, file_name, destiny_path, charge_type='None'): + """ + Convert a file in format `.json` to `.cif`. + + Parameters + ---------- + origin_path : str + Path to the '.json' file. + file_name : str + Name of the file. Does not neet to contain the `.json` extention. + destiny_path : str + path where the `.cif` file will be saved. + """ + + framework_JSON = read_json(origin_path, file_name) + + cell = framework_JSON['geometry']['cell_matrix'] + atom_labels = framework_JSON['geometry']['atom_labels'] + atom_pos = framework_JSON['geometry']['atom_pos'] + + if charge_type + '_charges' in list(framework_JSON['system'].keys()): + partial_charges = framework_JSON['geometry'][charge_type + '_charges'] + else: + partial_charges = False + + save_cif(destiny_path, + file_name, + cell, + atom_labels, + atom_pos, + partial_charges, + frac_coords=False) + + +def convert_gjf_2_xyz(path, file_name): + + file_name = file_name.split('.')[0] + + atom_labels, atom_pos = read_gjf(path, file_name + '.gjf') + + save_xyz(path, file_name + '.xyz', atom_labels, atom_pos) + + +def convert_xyz_2_gjf(path, file_name): + + file_name = file_name.split('.')[0] + + atom_labels, atom_pos = read_xyz(path, file_name + '.xyz') + + save_xyz(path, file_name + '.gjf', atom_labels, atom_pos) + + +def convert_cif_2_xyz(path, file_name, supercell=[1, 1, 1]): + + file_name = file_name.split('.')[0] + + structure = CifParser(os.path.join(path, file_name + '.cif')).get_structures(primitive=True)[0] + + structure.make_supercell([[supercell[0], 0, 0], [0, supercell[1], 0], [0, 0, supercell[2]]]) + + dict_sctructure = structure.as_dict() + + a, b, c = dict_sctructure['lattice']['a'] + b = dict_sctructure['lattice']['b'] + c = dict_sctructure['lattice']['c'] + + alpha = round(dict_sctructure['lattice']['alpha']) + beta = round(dict_sctructure['lattice']['beta']) + gamma = round(dict_sctructure['lattice']['gamma']) + + atom_labels = [i['label'] for i in dict_sctructure['sites']] + + atom_pos = [i['xyz'] for i in dict_sctructure['sites']] + + temp_file = open(os.path.join(path, file_name + '.xyz'), 'w') + temp_file.write(f'{len(atom_labels)} \n') + + temp_file.write(f'{a} {b} {c} {alpha} {beta} {gamma}\n') + + for i in range(len(atom_labels)): + temp_file.write('{:<5s}{:>15.7f}{:>15.7f}{:>15.7f}\n'.format(atom_labels[i], + atom_pos[i][0], + atom_pos[i][1], + atom_pos[i][2])) + + temp_file.close() + + +def write_json(path, name, COF_json): + + name = name.split('.')[0] + + if os.path.exists(path) is not True: + os.mkdir(path) + + save_path = os.path.join(path, name + '.cjson') + + with open(save_path, 'w', encoding='utf-8') as f: + simplejson.dump(COF_json, + f, + ensure_ascii=False, + separators=(',', ':'), + indent=2, + ignore_nan=True) + + +def read_json(path, name): + + cof_path = os.path.join(path, name + '.json') + + with open(cof_path, 'r') as r: + json_object = simplejson.loads(r.read()) + + return json_object + + +def create_COF_json(name) -> dict: + """ + Create a empty dictionary with the COF information. + """ + + system_info = 'Informations about the system.' + geometry_info = 'Informations about the geometry.' + optimization_info = 'Information about the optimization process.' + adsorption_info = 'Information about the adsorption simulation experiments on RASPA2' + textural_info = 'Information about the textural properties' + spectrum_info = 'Information about spectra simulation.' + experimental_info = 'Experimental data DRX, FTIR, ssNMR, UV-VIS...' + + COF_json = {'system': {'description': system_info, + 'name': name, + 'geo_opt': True, + 'execution_times_seconds': {}}, + 'geometry': {'description': geometry_info}, + 'optimization': {'description': optimization_info}, + 'adsorption': {'description': adsorption_info}, + 'textural': {'description': textural_info}, + 'spectrum': {'description': spectrum_info}, + 'experimental': {'description': experimental_info} + } + + return COF_json + + +def create_empty_CJSON() -> dict: + """ + Create a dictionary with the structure information to be saved using the + chemical JSON format. + """ + + chemJSON = { + "chemicalJson": 1, + "name": "", + "formula": "", + "unitCell": { + "a": 0.0, + "b": 0.0, + "c": 0.0, + "alpha": 0.0, + "beta": 0.0, + "gamma": 0.0, + "cellVectors": [] + }, + "atoms": { + "elements": { + "number": [], + "type": [] + }, + "coords": { + "3d": [], + "3dFractional": [] + }, + "formalCharges": [], + "labels": [] + }, + "bonds": { + "connections": { + "index": [] + }, + "order": [] + }, + "PartialCharges": {}, + "properties": { + "molecularMass": 0, + "totalCharge": 0, + "spinMultiplicity": 1, + "totalEnergy": 0, + "bandGap": 0, + }, + "spectra": {}, + "vibrations": {}, + "metadata": {}, + } + + return chemJSON + + +def create_structure_CJSON(StructureName: str, + CellParameters: list = None, + CellMatrix: list = None, + AtomTypes: list = None, + AtomPositions: list = None, + CartesianPositions: bool = False, + AtomLabels: list = [], + Bonds: list = [], + BondOrders: list = [], + PartialCharges: dict = {}, + ) -> dict: + """ + Creates a dictionary with the structure information to be saved using the + chemical JSON format. + + Parameters + ---------- + StructureName : str + Name of the structure. + CellParameters : list + List with the cell parameters. + CellMatrix : list + List with the cell matrix. Optional + AtomLabels : list + List with the atom labels. + AtomPositions : list + List with the atom positions. + + """ + + chemJSON = create_empty_CJSON() + + chemJSON['name'] = StructureName + chemJSON['formula'] = formula_from_atom_list(AtomTypes) + + if CellParameters is not None: + CellMatrix = cellpar_to_cell(CellParameters) + else: + CellParameters = cell_to_cellpar(CellMatrix) + CellMatrix = np.array(CellMatrix) + + chemJSON['unitCell']['a'] = CellParameters[0] + chemJSON['unitCell']['b'] = CellParameters[1] + chemJSON['unitCell']['c'] = CellParameters[2] + chemJSON['unitCell']['alpha'] = CellParameters[3] + chemJSON['unitCell']['beta'] = CellParameters[4] + chemJSON['unitCell']['gamma'] = CellParameters[5] + chemJSON['unitCell']['cellVectors'] = CellMatrix.flatten().tolist() + + AtomNumbers = [elements_dict(property="atomic_number")[i] for i in AtomTypes] + chemJSON['atoms']['elements']['number'] = AtomNumbers + chemJSON['atoms']['elements']['type'] = AtomTypes + + chemJSON['atoms']['elements']['labels'] = AtomLabels + + if CartesianPositions: + chemJSON['atoms']['coords']['3d'] = np.array(AtomPositions).flatten().tolist() + V_frac = get_cartesian_to_fractional_matrix(*CellParameters) + FracPosition = np.array([np.dot(V_frac, atom) for atom in AtomPositions]).flatten().tolist() + chemJSON['atoms']['coords']['3dFractional'] = FracPosition + + else: + chemJSON['atoms']['coords']['3dFractional'] = np.array(AtomPositions).flatten().tolist() + V_cart = get_fractional_to_cartesian_matrix(*CellParameters) + CartPosition = np.array([np.dot(V_cart, atom) for atom in AtomPositions]).flatten().tolist() + chemJSON['atoms']['coords']['3d'] = CartPosition + + chemJSON['atoms']['PartialCharges'] = PartialCharges + + chemJSON['bonds']['connections']['index'] = Bonds + chemJSON['bonds']['order'] = BondOrders + + return chemJSON + + +def generate_mol_dict(path, file_name, name, code, smiles): + + xsmiles, xsmiles_label, composition = smiles_to_xsmiles(smiles) + + if file_name.endswith('gjf'): + atom_types, atom_pos = read_gjf(path, file_name) + elif file_name.endswith('xyz'): + atom_types, atom_pos = read_xyz(path, file_name) + + mol_dict = { + "name": name, + "smiles": smiles, + "code": code, + "xsmiles": xsmiles, + "xsmiles_label": xsmiles_label, + "formula": composition, + "atoms": { + "elements": {"elementType": atom_types}, + "coords": {"3d": atom_pos.tolist()} + } + } + + print(mol_dict) + + write_json(path, file_name.split('.')[0], mol_dict) diff --git a/build/lib/logger.py b/build/lib/logger.py new file mode 100644 index 00000000..7a67aa74 --- /dev/null +++ b/build/lib/logger.py @@ -0,0 +1,87 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +""" +The logger creation function. +""" + +import logging.config +import logging.handlers + + +def create_logger(level: str = "DEBUG", + format: str = "detailed", + save_to_file: bool = False, + log_filename: str = 'pycofbuilder.log'): + """ + Build a logger with the given level and format. + + Parameters + ---------- + level : str + The logging level. Default is "DEBUG". + Can be one of "DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL". + format : str + The logging format. Default is "detailed". + Can be one of "simple", "detailed". + save_to_file : bool + Whether to save the logs to a file. Default is False. + log_filename : str + The file to save the logs to. Default is "pycofbuilder.log". + + Returns + ------- + logger : logging.Logger + The logger object. + """ + + # Check if the parameters are valid + allowed_levels = ["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"] + assert level.upper() in allowed_levels, "Invalid level, must be one of {}".format(allowed_levels) + + allowed_formats = ["simple", "detailed"] + assert format.lower() in allowed_formats, "Invalid format, must be one of {}".format(allowed_formats) + + # Set up logging + config_log = { + "version": 1, + "disable_existing_loggers": False, + "formatters": { + "simple": { + "format": "%(message)s" + }, + "detailed": { + "format": "[%(levelname)s|%(module)s|L%(lineno)d] %(asctime)s: %(message)s", + "datefmt": "%Y-%m-%dT%H:%M:%S%z" + } + }, + "handlers": { + "stderr": { + "class": "logging.StreamHandler", + "level": level.upper(), + "formatter": format.lower(), + "stream": "ext://sys.stderr" + }, + }, + "loggers": { + "root": {"level": level.upper(), "handlers": ["stderr"]} + } + } + + if save_to_file: + config_log["handlers"]["file"] = { + "class": "logging.FileHandler", + "level": level.upper(), + "formatter": format.lower(), + "filename": log_filename, + "mode": "a", + "encoding": "utf-8" + } + config_log["loggers"]["root"]["handlers"].append("file") + + logger = logging.getLogger("pycofbuilder") + + logging.config.dictConfig(config_log) + + return logger diff --git a/build/lib/tools.py b/build/lib/tools.py new file mode 100644 index 00000000..14f74f17 --- /dev/null +++ b/build/lib/tools.py @@ -0,0 +1,1006 @@ +# -*- coding: utf-8 -*- +# Created by Felipe Lopes de Oliveira +# Distributed under the terms of the MIT License. + +""" +This module contains the tools used by pyCOFBuilder. +""" + +import os +import simplejson +import numpy as np +from scipy.spatial import distance + + +def elements_dict(property='atomic_mass'): + '''Returns a dictionary containing the elements symbol and its selected property. + + Parameters + ---------- + prop : string + The desired property can be: + - "full_name" + - "atomic_number" + - "atomic_mass" + - "polarizability" + - "pauling_electronegativity" + - "thermo_electronegativity" + - "mulliken_electronegativity" + - "sanderson_electronegativity" + - "allen_electronegativity" + - "ghosh_electronegativity" + - "martynov_batsanov_electronegativity" + - "atomic_radius" + - "covalent_radius" + - "vdw_radius" + + Returns + ------- + prop_dic : dictionary + A dictionary containing the elements symbol and its respective property. + ''' + + file_name = os.path.join(os.path.dirname(__file__), 'data', 'periodic_table.json') + + with open(file_name, 'r') as f: + periodic_table = simplejson.load(f) + + prop_list = periodic_table['H'].keys() + + # Check if the property is valid + if property not in prop_list: + raise ValueError('Invalid property. Valid properties are: ' + ', '.join(prop_list)) + + prop_dic = {} + for element in periodic_table: + prop_dic[element] = periodic_table[element][property] + + return prop_dic + + +def unit_vector(vector): + """Return a unit vector in the same direction as x.""" + y = np.array(vector, dtype='float') + norm = y / np.linalg.norm(y) + return norm + + +def angle(v1, v2, unit='degree'): + """ + Calculates the angle between two vectors v1 and v2. + + Parameters + ---------- + v1 : array + (N,1) matrix with N dimensions + v2 : array + (N,1) matrix with N dimensions + unit : str + Unit of the output. Could be 'degree', 'radians' or 'cos'. + + Returns + ------- + angle : float + Angle in the selected unit. + """ + unit_vector1 = unit_vector(v1) + unit_vector2 = unit_vector(v2) + + dot_product = np.dot(unit_vector1, unit_vector2) + + if unit == 'degree': + angle = np.arccos(dot_product) * 180. / np.pi + if unit == 'radians': + angle = np.arccos(dot_product) + if unit == 'cos': + angle = dot_product + return angle + + +def rotation_matrix_from_vectors(vec1, vec2): + ''' + Find the rotation matrix that aligns vec1 to vec2 + + Parameters + ---------- + vec1 : array + (3,3) array + vec2 : array + (3,3) array + Returns + ------- + rotation_matrix : array + A transform matrix (3x3) which when applied to vec1, aligns it with vec2. + ''' + a, b = (vec1 / np.linalg.norm(vec1)).reshape(3), (vec2 / np.linalg.norm(vec2)).reshape(3) + v = np.cross(a, b) + c = np.dot(a, b) + s = np.linalg.norm(v) + if s != 0: + kmat = np.array([[0, -v[2], v[1]], + [v[2], 0, -v[0]], + [-v[1], v[0], 0]]) + + rotation_matrix = np.eye(3) + kmat + kmat.dot(kmat) * ((1 - c) / (s ** 2)) + return rotation_matrix + else: + return np.identity(3) + + +def rmsd(V, W): + """ + Calculate Root-mean-square deviation from two sets of vectors V and W. + Parameters + ---------- + V : array + (N,D) matrix, where N is points and D is dimension. + W : array + (N,D) matrix, where N is points and D is dimension. + Returns + ------- + rmsd : float + Root-mean-square deviation between the two vectors + """ + diff = np.array(V) - np.array(W) + N = len(V) + return np.sqrt((diff * diff).sum() / N) + + +def cell_to_cellpar(cell, radians=False): + """Returns the cell parameters [a, b, c, alpha, beta, gamma] + given a 3x3 cell matrix. + + Angles are in degrees unless radian=True is used. + + Parameters + ---------- + cell : array + (3,3) matrix of cell vectors v1, v2, and v3 + radians : bool + Return the cell angles in radians + + Returns + ------- + cellpar : array + (6,1) vector with the cell parameters + """ + lengths = [np.linalg.norm(v) for v in cell] + angles = [] + for i in range(3): + j = i - 1 + k = i - 2 + ll = lengths[j] * lengths[k] + if ll > 1e-16: + x = np.dot(cell[j], cell[k]) / ll + angle = 180.0 / np.pi * np.arccos(x) + else: + angle = 90.0 + angles.append(angle) + + # Corvet to radians if radians is True + if radians: + angles = [angle * np.pi / 180 for angle in angles] + + return np.array(lengths + angles) + + +def cellpar_to_cell(cellpar, ab_normal=(0, 0, 1), a_direction=None): + """Return a 3x3 cell matrix from cell parameters (a,b,c,alpha,beta, and gamma). + + Angles must be in degrees. + + The returned cell is orientated such that a and b are normal to `ab_normal` and a is + parallel to the projection of `a_direction` in the a-b plane. + + Default `a_direction` is (1,0,0), unless this is parallel to + `ab_normal`, in which case default `a_direction` is (0,0,1). + + The returned cell has the vectors va, vb and vc along the rows. The + cell will be oriented such that va and vb are normal to `ab_normal` + and va will be along the projection of `a_direction` onto the a-b + plane. + + Parameters + ---------- + cellpar : array + (6,1) vector with the cell parameters + ab_normal : array + Normal vector between a and b cell vectors. Default: (0, 0, 1) + a_direction : array + Specific direction for the a vector. Default: None + Returns + ------- + cell : array + (3,3) matrix of cell vectors v1, v2, and v3 + """ + if a_direction is None: + if np.linalg.norm(np.cross(ab_normal, (1, 0, 0))) < 1e-5: + a_direction = (0, 0, 1) + else: + a_direction = (1, 0, 0) + + # Define rotated X,Y,Z-system, with Z along ab_normal and X along the + # projection of a_direction onto the normal plane of Z. + ad = np.array(a_direction) + Z = unit_vector(ab_normal) + X = unit_vector(ad - np.dot(ad, Z) * Z) + Y = np.cross(Z, X) + + # Express va, vb and vc in the X,Y,Z-system + alpha, beta, gamma = 90., 90., 90. + if isinstance(cellpar, (int, float)): + a = b = c = cellpar + elif len(cellpar) == 1: + a = b = c = cellpar[0] + elif len(cellpar) == 3: + a, b, c = cellpar + else: + a, b, c, alpha, beta, gamma = cellpar + + # Handle orthorhombic cells separately to avoid rounding errors + eps = 2 * np.spacing(90.0, dtype=np.float64) # around 1.4e-14 + # alpha + if abs(abs(alpha) - 90) < eps: + cos_alpha = 0.0 + else: + cos_alpha = np.cos(alpha * np.pi / 180.0) + # beta + if abs(abs(beta) - 90) < eps: + cos_beta = 0.0 + else: + cos_beta = np.cos(beta * np.pi / 180.0) + # gamma + if abs(gamma - 90) < eps: + cos_gamma = 0.0 + sin_gamma = 1.0 + elif abs(gamma + 90) < eps: + cos_gamma = 0.0 + sin_gamma = -1.0 + else: + cos_gamma = np.cos(gamma * np.pi / 180.0) + sin_gamma = np.sin(gamma * np.pi / 180.0) + + # Build the cell vectors + va = a * np.array([1, 0, 0]) + vb = b * np.array([cos_gamma, sin_gamma, 0]) + cx = cos_beta + cy = (cos_alpha - cos_beta * cos_gamma) / sin_gamma + cz_sqr = 1. - cx * cx - cy * cy + assert cz_sqr >= 0 + cz = np.sqrt(cz_sqr) + vc = c * np.array([cx, cy, cz]) + + # Convert to the Cartesian x,y,z-system + abc = np.vstack((va, vb, vc)) + T = np.vstack((X, Y, Z)) + cell = np.dot(abc, T) + + return cell + + +def get_fractional_to_cartesian_matrix(cell_a: float, + cell_b: float, + cell_c: float, + alpha: float, + beta: float, + gamma: float, + angle_in_degrees: bool = True) -> np.array(float): + """ + Return the transformation matrix that converts fractional coordinates to + cartesian coordinates. + + Parameters + ---------- + a, b, c : float + The lengths of the edges. + alpha, gamma, beta : float + The angles between the sides. + angle_in_degrees : bool + True if alpha, beta and gamma are expressed in degrees. + Returns + ------- + T_matrix : ndarray + The 3x3 rotation matrix. ``V_cart = np.dot(T_matrix, V_frac)``. + """ + if angle_in_degrees: + alpha = np.deg2rad(alpha) + beta = np.deg2rad(beta) + gamma = np.deg2rad(gamma) + + cosa = np.cos(alpha) + cosb = np.cos(beta) + cosg = np.cos(gamma) + sing = np.sin(gamma) + + volume = np.sqrt(1.0 - cosa**2.0 - cosb**2.0 - cosg**2.0 + 2.0 * cosa * cosb * cosg) + + T_matrix = np.zeros((3, 3)) + + T_matrix[0, 0] = cell_a + T_matrix[0, 1] = cell_b * cosg + T_matrix[0, 2] = cell_c * cosb + T_matrix[1, 1] = cell_b * sing + T_matrix[1, 2] = cell_c * (cosa - cosb * cosg) / sing + T_matrix[2, 2] = cell_c * volume / sing + + return T_matrix + + +def get_cartesian_to_fractional_matrix(a: float, + b: float, + c: float, + alpha: float, + beta: float, + gamma: float, + angle_in_degrees: bool = True) -> np.array(float): + """ + Return the transformation matrix that converts cartesian coordinates to + fractional coordinates. + + Parameters + ---------- + a, b, c : float + The lengths of the edges. + alpha, gamma, beta : float + The angles between the sides. + angle_in_degrees : bool + True if alpha, beta and gamma are expressed in degrees. + + Returns + ------- + T_matrix : np.array + The 3x3 rotation matrix. ``R_frac = np.dot(T_matrix, R_cart)``. + """ + if angle_in_degrees: + alpha = np.deg2rad(alpha) + beta = np.deg2rad(beta) + gamma = np.deg2rad(gamma) + + cosa = np.cos(alpha) + cosb = np.cos(beta) + cosg = np.cos(gamma) + sing = np.sin(gamma) + + volume = np.sqrt(1.0 - cosa**2.0 - cosb**2.0 - cosg**2.0 + 2.0 * cosa * cosb * cosg) + + T_matrix = np.zeros((3, 3)) + + T_matrix[0, 0] = 1.0 / a + T_matrix[0, 1] = -cosg / (a * sing) + T_matrix[0, 2] = (cosa * cosg - cosb) / (a * volume * sing) + T_matrix[1, 1] = 1.0 / (b * sing) + T_matrix[1, 2] = (cosb * cosg - cosa) / (b * volume * sing) + T_matrix[2, 2] = sing / (c * volume) + + return T_matrix + + +def get_reciprocal_vectors(cell) -> tuple: + """ + Get the reciprocal vectors of a cell given in cell parameters of cell vectors. + + Parameters + ---------- + cell : array + (3,1) array for cell vectors or (6,1) array for cell parameters + + Returns + ------- + b1 : array + (3,1) array containing b_1 vector in the reciprocal space + b2 : array + (3,1) array containing b_2 vector in the reciprocal space + b3 : array + (3,1) array containing b_3 vector in the reciprocal space + """ + + if len(cell) == 3: + v1, v2, v3 = cell + if len(cell) == 6: + v1, v2, v3 = cellpar_to_cell(cell) + + vol = np.dot(v1, np.cross(v2, v3)) + + b1 = 2 * np.pi * np.cross(v2, v3) / vol + b2 = 2 * np.pi * np.cross(v3, v1) / vol + b3 = 2 * np.pi * np.cross(v1, v2) / vol + + return b1, b2, b3 + + +def get_kgrid(cell, distance=0.3) -> tuple: + """ + Get the k-points grid in the reciprocal space with a given distance for a + cell given in cell parameters of cell vectors. + + Parameters + ---------- + cell : array + (3,1) array for cell vectors or (6,1) array for cell parameters + distance : float + distance between the points in the reciprocal space + Returns + ------- + kx : int + Number of points in the x direction on reciprocal space + ky : int + Number of points in the y direction on reciprocal space + kz : int + Number of points in the z direction on reciprocal space + """ + + b1, b2, b3 = get_reciprocal_vectors(cell) + + b = np.array([np.linalg.norm(b1), np.linalg.norm(b2), np.linalg.norm(b3)]) + + kx = np.ceil(b[0]/distance).astype(int) + ky = np.ceil(b[1]/distance).astype(int) + kz = np.ceil(b[2]/distance).astype(int) + + return kx, ky, kz + + +def create_CellBox(A, B, C, alpha, beta, gamma): + """Creates the CellBox using the same expression as RASPA.""" + + tempd = (np.cos(alpha) - np.cos(gamma) * np.cos(beta)) / np.sin(gamma) + + ax = A + ay = 0 + az = 0 + bx = B * np.cos(gamma) + by = B * np.sin(gamma) + bz = 0 + cx = C * np.cos(beta) + cy = C * tempd + cz = C * np.sqrt(1 - np.cos(beta) ** 2 - tempd ** 2) + + CellBox = np.array([[ax, ay, az], + [bx, by, bz], + [cx, cy, cz]]) + + return CellBox + + +def calculate_UnitCells(cell, cutoff): + ''' + Calculate the number of unit cell repetitions so that all supercell lengths are larger than + twice the interaction potential cut-off radius. + + RASPA considers the perpendicular directions the directions perpendicular to the `ab`, `bc`, + and `ca` planes. Thus, the directions depend on who the crystallographic vectors `a`, `b`, + and `c` are and the length in the perpendicular directions would be the projections + of the crystallographic vectors on the vectors `a x b`, `b x c`, and `c x a`. + (here `x` means cross product) + + Parameters + ---------- + cell_matrix : array + (3,3) cell vectors or (6,1) + Returns + ------- + superCell + (3,1) list containg the number of repiting units in `x`, `y`, `z` directions. + ''' + + # Make sure that the cell is in the format of cell matrix + if len(cell) == 6: + cell_box = cellpar_to_cell(cell) + if len(cell) == 3: + cell_box = cell + + # Pre-calculate the cross products + axb = np.cross(cell_box[0], cell_box[1]) + bxc = np.cross(cell_box[1], cell_box[2]) + cxa = np.cross(cell_box[2], cell_box[0]) + + # Calculates the cell volume + V = np.dot(np.cross(cell_box[0], cell_box[1]), cell_box[2]) + + # Calculate perpendicular widths + cx = V / np.linalg.norm(bxc) + cy = V / np.linalg.norm(cxa) + cz = V / np.linalg.norm(axb) + + # Calculate UnitCells array + supercell = np.ceil(2.0 * cutoff / np.array([cx, cy, cz])).astype(int) + + return supercell + + +def cellpar_to_lammpsbox(a: float, + b: float, + c: float, + alpha: float, + beta: float, + gamma: float, + angle_in_degrees: bool = True): + """ + Return the box parameters lx, ly, lz, xy, xz, yz for LAMMPS data input. + Parameters + ---------- + a, b, c : float + The lengths of the edges. + alpha, gamma, beta : float + The angles between the sides. + angle_in_degrees : bool + True if alpha, beta and gamma are expressed in degrees. + Returns + ------- + r : array_like + The 1x6 array with the box parameters 'lx', 'ly', 'lz', 'xy', 'xz', 'yz'. + """ + if angle_in_degrees: + alpha = alpha*(np.pi/180) + beta = beta*(np.pi/180) + gamma = gamma*(np.pi/180) + + lx = a + xy = b * np.cos(gamma) + xz = c * np.cos(beta) + ly = np.sqrt(b ** 2 - xy ** 2) + yz = (b * c * np.cos(alpha) - xy * xz) / ly + lz = np.sqrt(c ** 2 - xz ** 2 - yz ** 2) + + return np.array([lx, ly, lz, xy, xz, yz]) + + +def find_index(element, e_list): + """ + Finds the index of a given element in a list + + Parameters + ---------- + element : string + String containing the label of the element in e_list + e_list : list + List with the atom labels + Returns + ---------- + i : int + The index of element in the e_list + """ + + index = None + for i in range(len(e_list)): + if np.array_equal(e_list[i], element): + index = i + break + return index + + +def change_X_atoms(atom_labels, atom_pos, bond_atom) -> tuple: + ''' + Changes the X atom for the desired bond_atom or remove it if bond_atom == 'R'. + + Parameters + ---------- + atom_labels : list + List containing the atom labels + atom_pos : list + List containing the atom position + Returns + ---------- + labels : list + List containing the processed atom labels + pos : list + List containing the processed atom position + ''' + label, pos = [], [] + + for i in range(len(atom_labels)): + if atom_labels[i] == 'X' and bond_atom != 'R': + label += [bond_atom] + pos += [atom_pos[i]] + if atom_labels[i] != 'X': + label += [atom_labels[i]] + pos += [atom_pos[i]] + + return label, pos + + +def closest_atom(label_1: str, pos_1: list, labels: list, pos: list): + ''' + Find the closest atom to a given atom + + Parameters + ---------- + label_1 : string + String containing the label of the atom + pos_1 : list + Array containing the position of the atom + labels : list + List containing the all the atom labels on the structure + pos : list + List containing the all the atom positions on the structure + + Returns + ---------- + closest_label : string + String containing the label of the closest atom + closest_position : array + Array containing the position of the closest atom + euclidian_distance : float + Euclidian distance between the two atoms + ''' + + list_labels = [] + list_pos = [] + + for i in range(len(labels)): + if labels[i] != label_1: + list_labels += [labels[i]] + list_pos += [pos[i]] + + if len(list_pos) == 0: + return None, np.array([0, 0, 0]), None + + closest_index = distance.cdist([pos_1], list_pos).argmin() + + closest_label = list_labels[closest_index] + closest_position = list_pos[closest_index] + euclidian_distance = np.linalg.norm(pos_1 - list_pos[closest_index]) + + return closest_label, closest_position, euclidian_distance + + +def closest_atom_struc(label_1, pos_1, labels, pos): + '''Finds the closest atom on the structure to a given atom''' + + list_labels = [] + list_pos = [] + for i in range(len(labels)): + if labels[i] != label_1: + if 'C' in labels[i]: + list_labels += [labels[i]] + list_pos += [pos[i]] + + closest_index = distance.cdist([pos_1], list_pos).argmin() + + closet_label = list_labels[closest_index] + closet_position = list_pos[closest_index] + euclidian_distance = np.linalg.norm(pos_1 - list_pos[closest_index]) + + return closet_label, closet_position, euclidian_distance + + +def get_bond_atom(connector_1: str, connector_2: str) -> str: + ''' + Get the atom that will be used to bond two building blocks. + ''' + + bond_dict = { + 'NH2': 'N', + 'NHOH': 'N', + 'COCHCHOH': 'N', + 'CONHNH2': 'N', + 'CHNNH2': 'N', + 'COOH': 'N', + 'BOH2': 'B', + 'OH2': 'B', + 'Cl': 'X', + 'Br': 'X', + 'CH2CN': 'C', + 'CH3': 'C' + } + + bond_atom = None + for group in list(bond_dict.keys()): + if group in [connector_1, connector_2]: + bond_atom = bond_dict[group] + + return bond_atom + + +def get_framework_symm_text(name, lattice, hall, space_group, space_number, symm_op): + '''Get the text for the framework symmop''' + text = '{:<60s} {:^12s} {:<4s} {:^4s} #{:^5s} {:^2} sym. op.'.format(name, + lattice, + hall.lstrip('-'), + space_group, + space_number, + symm_op) + return text + + +def print_framework_name(name, lattice, hall, space_group, space_number, symm_op): + '''Print the results of the created structure''' + print('{:<60s} {:^12s} {:<4s} {:^4s} #{:^5s} {:^2} sym. op.'.format(name, + lattice, + hall.lstrip('-'), + space_group, + space_number, + symm_op)) + + +def print_command(text, verbose, match): + if verbose in match: + print(text) + + +def formula_from_atom_list(AtomLabels: list) -> str: + """ + Create a string with the formula of the structure from the list of atoms. + + Parameters + ---------- + AtomLabels : list + List of strings containing the atom labels. + + Returns + ------- + formula : str + String with the formula of the structure. + """ + + formula = '' + for i in set(AtomLabels): + formula += i + str(AtomLabels.count(i)) + + return formula + + +def smiles_to_xsmiles(smiles_string: str) -> str: + ''' + Converts a SMILES string to an extended SMILES string with labels + + Parameters + ---------- + smiles_string : str + SMILES string to be converted + + Returns + ------- + xsmiles : str + Extended SMILES string with special labels + xsmiles_label : str + xsmiles labels for images with the special labels + composition : str + String containing the composition + ''' + SPECIAL_ATOMS = ['Q', 'R', 'X'] + REGULAR_ATOMS = ['C', 'N', 'H', 'O', 'S', 'B'] + + xsmiles = '' + labels = [] + atom_list = [] + + for i, letter in enumerate(smiles_string): + + if letter in SPECIAL_ATOMS: + xsmiles += '*' + labels.append(letter) + if letter == 'R': + atom_list.append(smiles_string[i:i+2]) + else: + atom_list.append(letter) + + elif letter.isnumeric(): + if smiles_string[i-1] == 'R': + labels[-1] = labels[-1] + letter + else: + xsmiles += letter + + elif letter in REGULAR_ATOMS: + xsmiles += letter + labels += [''] + atom_list.append(letter) + + else: + xsmiles += letter + + # Generate the xsmiles label + xsmiles_label = '|$' + ';'.join(labels) + '$|' + + # Generate the composition + composition = formula_from_atom_list(atom_list) + + return xsmiles, xsmiles_label, composition + + +def ibrav_to_cell(ibrav, celldm1, celldm2, celldm3, celldm4, celldm5, celldm6): + """ + Convert a value of ibrav to a cell. + + Parameters + ---------- + ibrav : int + celldmx: float + + Returns + ------- + cell : matrix + The cell as a 3x3 numpy array + """ + + alat = celldm1 * 0.5291772105638411 + + if ibrav == 1: + cell = np.identity(3) * alat + elif ibrav == 2: + cell = np.array([[-1.0, 0.0, 1.0], + [0.0, 1.0, 1.0], + [-1.0, 1.0, 0.0]]) * (alat / 2) + elif ibrav == 3: + cell = np.array([[1.0, 1.0, 1.0], + [-1.0, 1.0, 1.0], + [-1.0, -1.0, 1.0]]) * (alat / 2) + elif ibrav == -3: + cell = np.array([[-1.0, 1.0, 1.0], + [1.0, -1.0, 1.0], + [1.0, 1.0, -1.0]]) * (alat / 2) + elif ibrav == 4: + cell = np.array([[1.0, 0.0, 0.0], + [-0.5, 0.5 * 3**0.5, 0.0], + [0.0, 0.0, celldm3]]) * alat + elif ibrav == 5: + tx = ((1.0 - celldm4) / 2.0)**0.5 + ty = ((1.0 - celldm4) / 6.0)**0.5 + tz = ((1 + 2 * celldm4) / 3.0)**0.5 + cell = np.array([[tx, -ty, tz], + [0, 2 * ty, tz], + [-tx, -ty, tz]]) * alat + elif ibrav == -5: + ty = ((1.0 - celldm4) / 6.0)**0.5 + tz = ((1 + 2 * celldm4) / 3.0)**0.5 + a_prime = alat / 3**0.5 + u = tz - 2 * 2**0.5 * ty + v = tz + 2**0.5 * ty + cell = np.array([[u, v, v], + [v, u, v], + [v, v, u]]) * a_prime + elif ibrav == 6: + cell = np.array([[1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, celldm3]]) * alat + elif ibrav == 7: + cell = np.array([[1.0, -1.0, celldm3], + [1.0, 1.0, celldm3], + [-1.0, -1.0, celldm3]]) * (alat / 2) + elif ibrav == 8: + cell = np.array([[1.0, 0.0, 0.0], + [0.0, celldm2, 0.0], + [0.0, 0.0, celldm3]]) * alat + elif ibrav == 9: + cell = np.array([[1.0 / 2.0, celldm2 / 2.0, 0.0], + [-1.0 / 2.0, celldm2 / 2.0, 0.0], + [0.0, 0.0, celldm3]]) * alat + elif ibrav == -9: + cell = np.array([[1.0 / 2.0, -celldm2 / 2.0, 0.0], + [1.0 / 2.0, celldm2 / 2.0, 0.0], + [0.0, 0.0, celldm3]]) * alat + elif ibrav == 10: + cell = np.array([[1.0 / 2.0, 0.0, celldm3 / 2.0], + [1.0 / 2.0, celldm2 / 2.0, 0.0], + [0.0, celldm2 / 2.0, celldm3 / 2.0]]) * alat + elif ibrav == 11: + cell = np.array([[1.0 / 2.0, celldm2 / 2.0, celldm3 / 2.0], + [-1.0 / 2.0, celldm2 / 2.0, celldm3 / 2.0], + [-1.0 / 2.0, -celldm2 / 2.0, celldm3 / 2.0]]) * alat + elif ibrav == 12: + sinab = (1.0 - celldm4**2)**0.5 + cell = np.array([[1.0, 0.0, 0.0], + [celldm2 * celldm4, celldm2 * sinab, 0.0], + [0.0, 0.0, celldm3]]) * alat + elif ibrav == -12: + sinac = (1.0 - celldm5**2)**0.5 + cell = np.array([[1.0, 0.0, 0.0], + [0.0, celldm2, 0.0], + [celldm3 * celldm5, 0.0, celldm3 * sinac]]) * alat + elif ibrav == 13: + sinab = (1.0 - celldm4**2)**0.5 + cell = np.array([[1.0 / 2.0, 0.0, -celldm3 / 2.0], + [celldm2 * celldm4, celldm2 * sinab, 0.0], + [1.0 / 2.0, 0.0, celldm3 / 2.0]]) * alat + elif ibrav == 14: + sinab = (1.0 - celldm4**2)**0.5 + v3 = [celldm3 * celldm5, + celldm3 * (celldm6 - celldm5 * celldm4) / sinab, + celldm3 * ((1 + 2 * celldm6 * celldm5 * celldm4 + - celldm6**2 - celldm5**2 - celldm4**2)**0.5) / sinab] + cell = np.array([[1.0, 0.0, 0.0], + [celldm2 * celldm4, celldm2 * sinab, 0.0], + v3]) * alat + else: + raise NotImplementedError('ibrav = {0} is not implemented'.format(ibrav)) + + return cell + + +def equal_value(val1, val2, threshold=1e-3) -> bool: + ''' + Determine if two values are equal based on a given threshold. + ''' + return abs(val1 - val2) <= threshold + + +def classify_unit_cell(cell, thr=1e-3) -> str: + ''' + Determine the bravais lattice based on the cell lattice. + The cell lattice can be the cell parameters as (6,1) array or + the cell vectors as (3x3) array. + + Bravais lattice can be cubic, tetragonal, orthorhombic, hexagonal, + monoclinic, or triclinic. + + Parameters + ---------- + cell : array + Array with the cell vectors or parameters + threshold: float + Numeric threshold for the analysis. Default: 1e-3 + + Returns + ------- + cell_type : string + Bravais lattice. + ''' + + if len(cell) == 3: + a, b, c, alpha, beta, gamma = cell_to_cellpar(cell) + + cell_type = None + + if equal_value(alpha, 90, thr) and equal_value(beta, 90, thr) and equal_value(gamma, 90, thr): + if equal_value(a, b, thr) and equal_value(b, c, thr): + cell_type = "cubic" + if equal_value(a, b, thr) and not equal_value(a, c, thr): + cell_type = "tetragonal" + else: + cell_type = "orthorhombic" + elif equal_value(alpha, 90, thr) and equal_value(beta, 90, thr) and equal_value(gamma, 120, thr): + if equal_value(a, b, thr): + cell_type = "hexagonal" + elif equal_value(alpha, 90, thr) or equal_value(beta, 90, thr) or equal_value(gamma, 90, thr): + if not equal_value(a, b, thr) and not equal_value(b, c, thr) and not equal_value(a, c, thr): + cell_type = "monoclinic" + else: + cell_type = "triclinic" + + return cell_type + + +def cell_to_ibrav(cell): + ''' + Return the ibrav number for a given cell. + ''' + + if len(cell) == 3: + a, b, c, alpha, beta, gamma = cell_to_cellpar(cell) + else: + a, b, c, alpha, beta, gamma = cell + + cell_type = classify_unit_cell(cell) + + if cell_type == 'cubic': + celldm = {'ibrav': 1, + 'celldm(1)': a / 0.5291772105638411} + elif cell_type == 'hexagonal': + celldm = {'ibrav': 4, + 'celldm(1)': a / 0.5291772105638411, + 'celldm(3)': c / a} + elif cell_type == 'tetragonal': + celldm = {'ibrav': 6, + 'celldm(1)': a / 0.5291772105638411, + 'celldm(3)': c / a} + elif cell_type == 'orthorhombic': + celldm = {'ibrav': 8, + 'celldm(1)': a / 0.5291772105638411, + 'celldm(2)': b / a, + 'celldm(3)': c / a} + elif cell_type == 'monoclinic': + celldm = {'ibrav': 12, + 'celldm(1)': a / 0.5291772105638411, + 'celldm(2)': b / a, + 'celldm(3)': c / a, + 'celldm(4)': np.cos(np.deg2rad(beta))} + else: + celldm = {'ibrav': 14, + 'celldm(1)': a / 0.5291772105638411, + 'celldm(2)': b / a, + 'celldm(3)': c / a, + 'celldm(4)': np.cos(np.deg2rad(alpha)), + 'celldm(5)': np.cos(np.deg2rad(beta)), + 'celldm(6)': np.cos(np.deg2rad(gamma))} + + return celldm diff --git a/dist/pycofbuilder-0.0.7-py3-none-any.whl b/dist/pycofbuilder-0.0.7-py3-none-any.whl new file mode 100644 index 00000000..911a9ee3 Binary files /dev/null and b/dist/pycofbuilder-0.0.7-py3-none-any.whl differ diff --git a/dist/pycofbuilder-0.0.7.tar.gz b/dist/pycofbuilder-0.0.7.tar.gz new file mode 100644 index 00000000..0521ee81 Binary files /dev/null and b/dist/pycofbuilder-0.0.7.tar.gz differ diff --git a/pyproject.toml b/pyproject.toml index f73d4449..d34616b8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,18 +1,17 @@ [project] name = "pycofbuilder" -version = '0.0.6' +version = '0.0.7' authors = [ { name="Felipe Lopes", email="felipe.lopes@nano.ufrj.br" }, ] description = 'A package for Covalent Organic Frameworks sturcture creation based on the reticular approach.' readme = "README.md" -requires-python = ">=3.8" +requires-python = ">=3.10" classifiers = [ 'Programming Language :: Python', 'Programming Language :: Python :: 3', - 'Programming Language :: Python :: 3.9', + 'Programming Language :: Python :: 3.10', "Intended Audience :: Science/Research", - "License :: MIT License", "Operating System :: OS Independent", "Topic :: Scientific/Engineering :: Information Analysis", "Topic :: Scientific/Engineering :: Physics", @@ -22,4 +21,5 @@ classifiers = [ [project.urls] Homepage = "https://github.com/lipelopesoliveira/pyCOFBuilder" -Issues = "https://github.com/lipelopesoliveira/pyCOFBuilder/issues" \ No newline at end of file +Issues = "https://github.com/lipelopesoliveira/pyCOFBuilder/issues" +Documentation = "https://lipelopesoliveira.github.io/pyCOFBuilder/docs/_build/html/index.html" \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index f24a20cd..73a1b57b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ +python >= 3.10 os numpy >= 1.6.3 scipy >= 1.2 diff --git a/setup.py b/setup.py index 8c93e2fa..dad21e7c 100644 --- a/setup.py +++ b/setup.py @@ -6,7 +6,7 @@ license = f.read() -VERSION = '0.0.6' +VERSION = '0.0.7' DESCRIPTION = 'A package for Covalent Organic Frameworks sturcture assembly.' setup( diff --git a/src/pycofbuilder/pycofbuilder.egg-info/PKG-INFO b/src/pycofbuilder/pycofbuilder.egg-info/PKG-INFO new file mode 100644 index 00000000..787c1d87 --- /dev/null +++ b/src/pycofbuilder/pycofbuilder.egg-info/PKG-INFO @@ -0,0 +1,205 @@ +Metadata-Version: 2.1 +Name: pycofbuilder +Version: 0.0.7 +Summary: A package for Covalent Organic Frameworks sturcture creation based on the reticular approach. +Home-page: https://github.com/lipelopesoliveira/pyCOFBuilder +Author: Felipe Lopes Oliveira +Author-email: Felipe Lopes +License: MIT License + + Copyright (c) 2023, Felipe Lopes de Oliveira + + Permission is hereby granted, free of charge, to any person obtaining a copy + of this software and associated documentation files (the "Software"), to deal + in the Software without restriction, including without limitation the rights + to use, copy, modify, merge, publish, distribute, sublicense, and/or sell + copies of the Software, and to permit persons to whom the Software is + furnished to do so, subject to the following conditions: + + The above copyright notice and this permission notice shall be included in all + copies or substantial portions of the Software. + + THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR + IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, + FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE + AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER + LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, + OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE + SOFTWARE. +Project-URL: Homepage, https://github.com/lipelopesoliveira/pyCOFBuilder +Project-URL: Issues, https://github.com/lipelopesoliveira/pyCOFBuilder/issues +Project-URL: Documentation, https://lipelopesoliveira.github.io/pyCOFBuilder/docs/_build/html/index.html +Classifier: Programming Language :: Python +Classifier: Programming Language :: Python :: 3 +Classifier: Programming Language :: Python :: 3.10 +Classifier: Intended Audience :: Science/Research +Classifier: Operating System :: OS Independent +Classifier: Topic :: Scientific/Engineering :: Information Analysis +Classifier: Topic :: Scientific/Engineering :: Physics +Classifier: Topic :: Scientific/Engineering :: Chemistry +Classifier: Topic :: Software Development :: Libraries :: Python Modules +Requires-Python: >=3.10 +Description-Content-Type: text/markdown +License-File: LICENSE +Requires-Dist: os +Requires-Dist: math +Requires-Dist: simplejason +Requires-Dist: numpy>=1.6.3 +Requires-Dist: scipy>=1.2 +Requires-Dist: pymatgen>=2022.0.8 +Requires-Dist: jupyter +Requires-Dist: jupyterlab +Requires-Dist: pandas +Requires-Dist: tqdm +Requires-Dist: gemmi +Requires-Dist: ase + +# pyCOFBuilder + +![puCOFBuilder](docs/img/header.png) + +**pyCOFBuilder** is a simple and powerful python package to automatically assembly COF structures with specifics building blocks, topologies, and functionalizations following the reticular approach to build and represent COF structures. The project was developed to address the need for generation of COFs structures in a high-throughput style, based on a nomenclature tha allows direct sctructural feature interpretation from a simple name. The package uses [pymatgen](https://pymatgen.org/) to create the structures. + +This package is still under development and, but it is already possible to create a large number of COFs structures. + +Learn more at the [Documentation](https://lipelopesoliveira.github.io/pyCOFBuilder/index.html) + +## Requirements + +0. Python >= 3.10 +1. pymatgen >= 2022.0.0 +2. numpy >= 1.2 +3. scipy >= 1.6.3 +4. simplejson +5. ase +6. gemmi + +The Python dependencies are most easily satisfied using a conda +([anaconda](https://www.anaconda.com/distribution)/[miniconda](https://docs.conda.io/en/latest/miniconda.html)) +installation by running + +```Shell +conda env create --file environment.yml +``` + +## Installation + +Currently the best way to use pyCOFBuilder is to manually import it using the `sys` module, as exemplified below: + +```python +# importing module +import sys + +# appending a path +sys.path.append('{PATH_TO_PYCOFBUILDER}/pyCOFBuilder/src') + +import pycofbuilder as pcb +``` + +Just remember to change the `{PATH_TO_PYCOFBUILDER}` to the directory where you download the pyCOFBuilder package. + +## Basic Usage + +To create a specific COF, such as `T3_BENZ_NH2_OH-L2_BENZ_CHO_H-HCB_A-AA`: + +```python +# importing module +import sys + +# appending a path +sys.path.append('{PATH_TO_PYCOFBUILDER}/pyCOFBuilder/src') + +import pycofbuilder as pcb + +cof = pcb.Framework('T3_BENZ_CHO_OH-L2_BENZ_NH2_H-HCB_A-AA') +cof.save(fmt='cif', supercell = [1, 1, 2], save_dir = '.') +``` + +You should see an output such as: + +```python +T3_BENZ_NH2_OH-L2_BENZ_CHO_H_H-HCB_A-AA hexagonal P P6/m # 175 12 sym. op. +``` + +A `.cif` file (the default save format is CIF, but it can be easily changed by setting other value on the `fmt` option) will be created in the `out` folder. The code will print out some information about the structure created. + +Currently, it is possible to select the following formats: + +- `cif` +- `xsf` +- `pdb` +- `cjson` +- `vasp` +- `turbomole` +- `pqr` +- `qe` +- `gjf` +- `xyz` + +Besides, the variable `structure` now is a `Framework` object. This object has some attributes that can be accessed: + +```python +>>> cof.name +'T3_BENZ_NH2_OH-L2_BENZ_CHO_H-HCB_A-AA' +>>> cof.smiles +'(N)C1=C(O)C((N))=C(O)C((N))=C1O.(C([H])=O)C1=C([H])C([H])=C((C([H])=O))C([H])=C1[H]' +>>> cof.lattice +array([[ 22.49540055, 0. , 0. ], + [-11.24770028, 19.48158835, 0. ], + [ 0. , 0. , 3.6 ]]) +>>> cof.n_atoms +72 +>>> cof.space_group +'P6/m' +``` + +## COFs and Building Blocks nomenclature + +In order to ensure greater reproducibility as well as quickly and easily access to relevant information from the COFs, I've developed a simple nomenclature to name the structure. Generally speaking, a COF can be described as + +### `Building_Block_1`-`Building_Block_2`-`Net`-`Stacking` + +where: + +- `Building_Block_1`: The building block with the greater connectivity. +- `Building_Block_2`: The building block with the smaller connectivity. +- `Net`: The net describing the reticular structure. +- `Stacking`: The stacking (for 2D structures) or interpenetrating degree (for 3D structures) + +To name the building blocks I also developed a set of rules. The building block can be described as + +### `Symmetry`\_`Core`\_`Connector`\_`RadicalGroupR1`\_`RadicalGroupR2`\_`RadicalGroupR3`\_`...` + +where: + +- `Symmetry`: The general symmetry of the building block. Also represents the connectivity of the building block. For 2D building blocks can be `L2`, `T3` or `S4`, and for 3D building blocks can be `D4`. +- `Core`: The 4 letters code referring to the building block core. +- `Connector`: The type of functional group that will be used to assembly the COF structure. Ex.: `NH2`, `CHO`, `CONHNH2`, etc. +- `RadicalGroupRN`: The Nth radical group in the structure. The number of Radical groups will change according to the availability of the core. + +Note that every "card" for the building block name is separated by an underline (\_) and every "card" for the COF name is separated by a dash (-). This makes it easy to split the COF name into useful information. + +## Current available Building Blocks + +![Ditopic](docs/img/L2_1.png) +![Ditopic](docs/img/L2_2.png) +![Tritopic](docs/img/T3.png) +![Tetratopic](docs/img/S4.png) +![Hexatopic](docs/img/H6.png) + +## Current available Connector Groups + +![Connection groups](docs/img/Q_GROUPS.png) + +## Current available R Groups + +![Functional Groups](docs/img/R_GROUPS.png) + +## Citation + +If you find **pyCOFBuilder** useful in your research please consider citing the following paper: + +> F. L. Oliveira and P. M. Esteves, +> _pyCOFBuilder: A python package for automated creation of Covalent Organic Framework models based on the reticular approach_ +> +> _arxiv.org/abs/2310.14822_ [DOI](https://doi.org/10.48550/arXiv.2310.14822) diff --git a/src/pycofbuilder/pycofbuilder.egg-info/SOURCES.txt b/src/pycofbuilder/pycofbuilder.egg-info/SOURCES.txt new file mode 100644 index 00000000..e05419be --- /dev/null +++ b/src/pycofbuilder/pycofbuilder.egg-info/SOURCES.txt @@ -0,0 +1,21 @@ +LICENSE +MANIFEST.in +README.md +pyproject.toml +setup.py +src/pycofbuilder/__init__.py +src/pycofbuilder/building_block.py +src/pycofbuilder/cjson.py +src/pycofbuilder/exceptions.py +src/pycofbuilder/framework.py +src/pycofbuilder/io_tools.py +src/pycofbuilder/logger.py +src/pycofbuilder/tools.py +src/pycofbuilder/data/topology.py +src/pycofbuilder/pycofbuilder.egg-info/PKG-INFO +src/pycofbuilder/pycofbuilder.egg-info/SOURCES.txt +src/pycofbuilder/pycofbuilder.egg-info/dependency_links.txt +src/pycofbuilder/pycofbuilder.egg-info/requires.txt +src/pycofbuilder/pycofbuilder.egg-info/top_level.txt +tests/test_advanced.py +tests/test_basic.py \ No newline at end of file diff --git a/src/pycofbuilder/pycofbuilder.egg-info/dependency_links.txt b/src/pycofbuilder/pycofbuilder.egg-info/dependency_links.txt new file mode 100644 index 00000000..8b137891 --- /dev/null +++ b/src/pycofbuilder/pycofbuilder.egg-info/dependency_links.txt @@ -0,0 +1 @@ + diff --git a/src/pycofbuilder/pycofbuilder.egg-info/requires.txt b/src/pycofbuilder/pycofbuilder.egg-info/requires.txt new file mode 100644 index 00000000..a010dcd8 --- /dev/null +++ b/src/pycofbuilder/pycofbuilder.egg-info/requires.txt @@ -0,0 +1,12 @@ +os +math +simplejason +numpy>=1.6.3 +scipy>=1.2 +pymatgen>=2022.0.8 +jupyter +jupyterlab +pandas +tqdm +gemmi +ase diff --git a/src/pycofbuilder/pycofbuilder.egg-info/top_level.txt b/src/pycofbuilder/pycofbuilder.egg-info/top_level.txt new file mode 100644 index 00000000..01f7e5e9 --- /dev/null +++ b/src/pycofbuilder/pycofbuilder.egg-info/top_level.txt @@ -0,0 +1,9 @@ +__init__ +building_block +cjson +data +exceptions +framework +io_tools +logger +tools