def discrete_molecules(system, rebuild=None, tol=<NUM_LIT>):
<EOL>if rebuild is not None:<EOL><INDENT>mode = <NUM_LIT:3><EOL><DEDENT>else:<EOL><INDENT>if '<STR_LIT>' in system.keys():<EOL><INDENT>if system['<STR_LIT>'].shape == (<NUM_LIT:6>,):<EOL><INDENT>mode = <NUM_LIT:2><EOL><DEDENT>else:<EOL><INDENT>mode = <NUM_LIT:1><EOL><DEDENT><DEDENT>elif '<STR_LIT>' in system.keys():<EOL><INDENT>if system['<STR_LIT>'].shape == (<NUM_LIT:3>, <NUM_LIT:3>):<EOL><INDENT>mode = <NUM_LIT:2><EOL><DEDENT>else:<EOL><INDENT>mode = <NUM_LIT:1><EOL><DEDENT><DEDENT>else:<EOL><INDENT>mode = <NUM_LIT:1><EOL><DEDENT><DEDENT>try:<EOL><INDENT>elements = system['<STR_LIT>']<EOL>coordinates = system['<STR_LIT>']<EOL><DEDENT>except KeyError:<EOL><INDENT>raise _FunctionError(<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>"<STR_LIT>"<EOL>)<EOL><DEDENT>coordinates = system['<STR_LIT>']<EOL>args = (elements, coordinates)<EOL>adj = <NUM_LIT:0><EOL>if '<STR_LIT>' in system.keys():<EOL><INDENT>atom_ids = system['<STR_LIT>']<EOL>args = (elements, atom_ids, coordinates)<EOL>adj = <NUM_LIT:1><EOL><DEDENT>atom_list = compose_atom_list(*args)<EOL>atom_coor = decompose_atom_list(atom_list)[<NUM_LIT:1> + adj]<EOL>if mode == <NUM_LIT:2> or mode == <NUM_LIT:3>:<EOL><INDENT>origin = np.array([<NUM_LIT>, <NUM_LIT:0.>, <NUM_LIT:0.>])<EOL>if '<STR_LIT>' not in system.keys():<EOL><INDENT>matrix = unit_cell_to_lattice_array(system['<STR_LIT>'])<EOL><DEDENT>else:<EOL><INDENT>matrix = system['<STR_LIT>']<EOL><DEDENT>pseudo_origin_frac = np.array([<NUM_LIT>, <NUM_LIT>, <NUM_LIT>])<EOL>pseudo_origin = cartisian_from_fractional(pseudo_origin_frac, matrix)<EOL>if rebuild is not None:<EOL><INDENT>selements = rebuild['<STR_LIT>']<EOL>sids = rebuild['<STR_LIT>']<EOL>scoordinates = rebuild['<STR_LIT>']<EOL>satom_list = compose_atom_list(selements, sids, scoordinates)<EOL>satom_coor = decompose_atom_list(satom_list)[<NUM_LIT:1> + adj]<EOL><DEDENT>system_com = center_of_mass(elements, coordinates)<EOL>if np.allclose(system_com, origin, atol=<NUM_LIT>):<EOL><INDENT>boundary = np.array([-<NUM_LIT:0.5>, <NUM_LIT:0.5>])<EOL><DEDENT>else:<EOL><INDENT>boundary = np.array([<NUM_LIT:0.>, <NUM_LIT:1.>])<EOL><DEDENT><DEDENT>else:<EOL><INDENT>pseudo_origin = center_of_mass(<EOL>elements, coordinates) + np.array([<NUM_LIT>, <NUM_LIT:0.>, <NUM_LIT:0.>])<EOL><DEDENT>molecules = []<EOL>exceptions = ['<STR_LIT:H>', '<STR_LIT>', '<STR_LIT>', '<STR_LIT:F>', '<STR_LIT>', '<STR_LIT>', '<STR_LIT>', '<STR_LIT>', '<STR_LIT>', '<STR_LIT>']<EOL>set_of_elements = set(system['<STR_LIT>'])<EOL>max_r_cov = max([<EOL>atomic_covalent_radius[i.upper()] for i in set_of_elements])<EOL>max_dist = <NUM_LIT:2> * max_r_cov + tol<EOL>while atom_list:<EOL><INDENT>inside_atoms_heavy = [<EOL>i for i in atom_list if i[<NUM_LIT:0>].upper() not in exceptions<EOL>]<EOL>if inside_atoms_heavy:<EOL><INDENT>inside_atoms_coord_heavy = decompose_atom_list(inside_atoms_heavy)[<EOL><NUM_LIT:1> + adj]<EOL>dist_matrix = euclidean_distances(inside_atoms_coord_heavy,<EOL>pseudo_origin.reshape(<NUM_LIT:1>, -<NUM_LIT:1>))<EOL>atom_index_x, _ = np.unravel_index(dist_matrix.argmin(),<EOL>dist_matrix.shape)<EOL>potential_starting_point = inside_atoms_heavy[atom_index_x]<EOL>pot_arr = np.array(potential_starting_point[<NUM_LIT:1> + adj:])<EOL>dist_matrix = euclidean_distances(<EOL>atom_coor, pot_arr.reshape(<NUM_LIT:1>, -<NUM_LIT:1>)<EOL>)<EOL>idx = (dist_matrix > <NUM_LIT:0.1>) * (dist_matrix < max_dist)<EOL>if len(idx) < <NUM_LIT:1>:<EOL><INDENT>pass<EOL><DEDENT>else:<EOL><INDENT>working_list = [potential_starting_point]<EOL><DEDENT><DEDENT>else:<EOL><INDENT>break<EOL><DEDENT>final_molecule = []<EOL>while working_list:<EOL><INDENT>working_list_temp = []<EOL>try:<EOL><INDENT>atom_coor = decompose_atom_list(atom_list)[<NUM_LIT:1> + adj]<EOL><DEDENT>except _FunctionError:<EOL><INDENT>atom_coor = None<EOL><DEDENT>for i in working_list:<EOL><INDENT>if i[<NUM_LIT:0>].upper() not in exceptions:<EOL><INDENT>i_arr = np.array(i[<NUM_LIT:1> + adj:])<EOL>if atom_coor is not None:<EOL><INDENT>dist_matrix = euclidean_distances(<EOL>atom_coor, i_arr.reshape(<NUM_LIT:1>, -<NUM_LIT:1>)<EOL>)<EOL>idx = (dist_matrix > <NUM_LIT:0.1>) * (dist_matrix < max_dist)<EOL>neighbours_indexes = np.where(idx)[<NUM_LIT:0>]<EOL>for j in neighbours_indexes:<EOL><INDENT>j_arr = np.array(atom_coor[j])<EOL>r_i_j = distance(i_arr, j_arr)<EOL>r_cov_i_j = atomic_covalent_radius[<EOL>i[<NUM_LIT:0>].upper()] + atomic_covalent_radius[<EOL>atom_list[j][<NUM_LIT:0>].upper()]<EOL>if r_cov_i_j - tol < r_i_j < r_cov_i_j + tol:<EOL><INDENT>working_list_temp.append(atom_list[j])<EOL><DEDENT><DEDENT><DEDENT>if rebuild is not None:<EOL><INDENT>sdist_matrix = euclidean_distances(<EOL>satom_coor, i_arr.reshape(<NUM_LIT:1>, -<NUM_LIT:1>))<EOL>sidx = (sdist_matrix > <NUM_LIT:0.1>) * (sdist_matrix < max_dist)<EOL>sneighbours_indexes = np.where(sidx)[<NUM_LIT:0>]<EOL>for j in sneighbours_indexes:<EOL><INDENT>if satom_list[j] in atom_list:<EOL><INDENT>pass<EOL><DEDENT>else:<EOL><INDENT>j_arr = np.array(satom_coor[j])<EOL>r_i_j = distance(i_arr, j_arr)<EOL>r_cov_i_j = atomic_covalent_radius[<EOL>i[<NUM_LIT:0>].upper()<EOL>] + atomic_covalent_radius[<EOL>satom_list[j][<NUM_LIT:0>].upper()]<EOL>if r_cov_i_j - tol < r_i_j < r_cov_i_j + tol:<EOL><INDENT>working_list_temp.append(satom_list[j])<EOL><DEDENT><DEDENT><DEDENT><DEDENT>final_molecule.append(i)<EOL><DEDENT>else:<EOL><INDENT>final_molecule.append(i)<EOL><DEDENT><DEDENT>for i in working_list:<EOL><INDENT>try:<EOL><INDENT>atom_list.remove(i)<EOL><DEDENT>except ValueError:<EOL><INDENT>pass<EOL><DEDENT><DEDENT>working_list = []<EOL>working_list_temp = unique(working_list_temp)<EOL>for i in working_list_temp:<EOL><INDENT>if i not in final_molecule:<EOL><INDENT>working_list.append(i)<EOL><DEDENT><DEDENT><DEDENT>final_molecule_dict = {}<EOL>final_molecule_dict['<STR_LIT>'] = np.array(<EOL>[x[<NUM_LIT:0>] for x in final_molecule], dtype='<STR_LIT:str>')<EOL>final_molecule_dict['<STR_LIT>'] = np.array(<EOL>[[*xyz[<NUM_LIT:1> + adj:]] for xyz in final_molecule])<EOL>if adj == <NUM_LIT:1>:<EOL><INDENT>final_molecule_dict['<STR_LIT>'] = np.array(<EOL>[x[<NUM_LIT:1>] for x in final_molecule], dtype='<STR_LIT:str>')<EOL><DEDENT>bool_ = True<EOL>if rebuild is not None:<EOL><INDENT>com = center_of_mass(final_molecule_dict['<STR_LIT>'],<EOL>final_molecule_dict['<STR_LIT>'])<EOL>com_frac = fractional_from_cartesian(com, matrix)[<NUM_LIT:0>]<EOL>com_frac_round = np.around(com_frac, decimals=<NUM_LIT:8>)<EOL>bool_ = np.all(np.logical_and(com_frac_round >= boundary[<NUM_LIT:0>],<EOL>com_frac_round < boundary[<NUM_LIT:1>]),<EOL>axis=<NUM_LIT:0>)<EOL><DEDENT>if bool(bool_) is True:<EOL><INDENT>molecules.append(final_molecule_dict)<EOL><DEDENT><DEDENT>return molecules<EOL>
Decompose molecular system into individual discreet molecules.

Note
----
New formula for bonds: (26/07/17)
The two atoms, x and y, are considered bonded if the distance between
them, calculated with distance matrix, is within the ranges:
.. :math:

    Rcov(x) + Rcov(y) - t < R(x,y) <  Rcov(x) + Rcov(y) + t

where Rcov is the covalent radius and the tolarenace (t) is set to
0.4 Angstrom.