commit 8ab593eecce6e63e883bcb6b746b543d4ce77225
Author: Carles <carles.marti-aliod@ens-lyon.fr>
Date:   Tue Apr 7 21:01:26 2020 +0200

    Implemented the read coordinates module working to read xyz and mol files and obtain rdkit and ase objects to work with

diff --git a/modules/read_coords.py b/modules/read_coords.py
new file mode 100644
index 0000000..12e0308
--- /dev/null
+++ b/modules/read_coords.py
@@ -0,0 +1,23 @@
+def read_coords(requirement, coord_file):
+    from ase.io.formats import filetype
+
+    req_vals = ['rdkit', 'ase']
+    file_type_vals = ['xyz', 'mol']
+    assert requirement in req_vals, ValueError('Not an adequate requirement.\n'
+                                    f'adequate requirements: {req_vals}')
+    assert filetype(coord_file) in file_type_vals, TypeError('The file formnat '
+                                                             'is not supported')
+
+    if requirement == 'rdkit':
+        if filetype(coord_file) == 'xyz':
+            from xyz2mol import xyz2mol, read_xyz_file
+            atoms, charge, xyz_coordinates = read_xyz_file(coord_file)
+            rd_mol_obj = xyz2mol(atoms, xyz_coordinates, charge=charge)
+            return rd_mol_obj
+        elif filetype(coord_file) == 'mol':
+            from rdkit.Chem import MolFromMolFile
+            return MolFromMolFile(coord_file)
+
+    if requirement == 'ase':
+        import ase.io
+        return ase.io.read(coord_file)
diff --git a/modules/xyz2mol.py b/modules/xyz2mol.py
new file mode 100755
index 0000000..a08eb8a
--- /dev/null
+++ b/modules/xyz2mol.py
@@ -0,0 +1,794 @@
+"""
+Module for generating rdkit molobj/smiles/molecular graph from free atoms
+
+Implementation by Jan H. Jensen, based on the paper
+
+    Yeonjoon Kim and Woo Youn Kim
+    "Universal Structure Conversion Method for Organic Molecules: From Atomic Connectivity
+    to Three-Dimensional Geometry"
+    Bull. Korean Chem. Soc. 2015, Vol. 36, 1769-1777
+    DOI: 10.1002/bkcs.10334
+
+"""
+
+import copy
+import itertools
+
+from rdkit.Chem import rdmolops
+try:
+    from rdkit.Chem import rdEHTTools #requires RDKit 2019.9.1 or later
+except ImportError:
+    rdEHTTools = None
+
+from collections import defaultdict
+
+import numpy as np
+import networkx as nx
+
+from rdkit import Chem
+from rdkit.Chem import AllChem, rdmolops
+
+global __ATOM_LIST__
+__ATOM_LIST__ = \
+    ['h',  'he',
+     'li', 'be', 'b',  'c',  'n',  'o',  'f',  'ne',
+     'na', 'mg', 'al', 'si', 'p',  's',  'cl', 'ar',
+     'k',  'ca', 'sc', 'ti', 'v ', 'cr', 'mn', 'fe', 'co', 'ni', 'cu',
+     'zn', 'ga', 'ge', 'as', 'se', 'br', 'kr',
+     'rb', 'sr', 'y',  'zr', 'nb', 'mo', 'tc', 'ru', 'rh', 'pd', 'ag',
+     'cd', 'in', 'sn', 'sb', 'te', 'i',  'xe',
+     'cs', 'ba', 'la', 'ce', 'pr', 'nd', 'pm', 'sm', 'eu', 'gd', 'tb', 'dy',
+     'ho', 'er', 'tm', 'yb', 'lu', 'hf', 'ta', 'w',  're', 'os', 'ir', 'pt',
+     'au', 'hg', 'tl', 'pb', 'bi', 'po', 'at', 'rn',
+     'fr', 'ra', 'ac', 'th', 'pa', 'u',  'np', 'pu']
+
+
+global atomic_valence
+global atomic_valence_electrons
+
+atomic_valence = defaultdict(list)
+atomic_valence[1] = [1]
+atomic_valence[6] = [4]
+atomic_valence[7] = [3,4]
+atomic_valence[8] = [2,1]
+atomic_valence[9] = [1]
+atomic_valence[14] = [4]
+atomic_valence[15] = [5,3] #[5,4,3]
+atomic_valence[16] = [6,3,2] #[6,4,2]
+atomic_valence[17] = [1]
+atomic_valence[32] = [4]
+atomic_valence[35] = [1]
+atomic_valence[53] = [1]
+
+atomic_valence_electrons = {}
+atomic_valence_electrons[1] = 1
+atomic_valence_electrons[6] = 4
+atomic_valence_electrons[7] = 5
+atomic_valence_electrons[8] = 6
+atomic_valence_electrons[9] = 7
+atomic_valence_electrons[14] = 4
+atomic_valence_electrons[15] = 5
+atomic_valence_electrons[16] = 6
+atomic_valence_electrons[17] = 7
+atomic_valence_electrons[32] = 4
+atomic_valence_electrons[35] = 7
+atomic_valence_electrons[53] = 7
+
+
+def str_atom(atom):
+    """
+    convert integer atom to string atom
+    """
+    global __ATOM_LIST__
+    atom = __ATOM_LIST__[atom - 1]
+    return atom
+
+
+def int_atom(atom):
+    """
+    convert str atom to integer atom
+    """
+    global __ATOM_LIST__
+    print(atom)
+    atom = atom.lower()
+    return __ATOM_LIST__.index(atom) + 1
+
+
+def get_UA(maxValence_list, valence_list):
+    """
+    """
+    UA = []
+    DU = []
+    for i, (maxValence, valence) in enumerate(zip(maxValence_list, valence_list)):
+        if not maxValence - valence > 0:
+            continue
+        UA.append(i)
+        DU.append(maxValence - valence)
+    return UA, DU
+
+
+def get_BO(AC, UA, DU, valences, UA_pairs, use_graph=True):
+    """
+    """
+    BO = AC.copy()
+    DU_save = []
+
+    while DU_save != DU:
+        for i, j in UA_pairs:
+            BO[i, j] += 1
+            BO[j, i] += 1
+
+        BO_valence = list(BO.sum(axis=1))
+        DU_save = copy.copy(DU)
+        UA, DU = get_UA(valences, BO_valence)
+        UA_pairs = get_UA_pairs(UA, AC, use_graph=use_graph)[0]
+
+    return BO
+
+
+def valences_not_too_large(BO, valences):
+    """
+    """
+    number_of_bonds_list = BO.sum(axis=1)
+    for valence, number_of_bonds in zip(valences, number_of_bonds_list):
+        if number_of_bonds > valence:
+            return False
+
+    return True
+
+
+def BO_is_OK(BO, AC, charge, DU, atomic_valence_electrons, atoms, valances,
+    allow_charged_fragments=True):
+    """
+    Sanity of bond-orders
+
+    args:
+        BO -
+        AC -
+        charge -
+        DU - 
+
+
+    optional
+        allow_charges_fragments - 
+
+
+    returns:
+        boolean - true of molecule is OK, false if not
+    """
+
+    if not valences_not_too_large(BO, valances):
+        return False
+
+    # total charge
+    Q = 0
+
+    # charge fragment list
+    q_list = []
+
+    if allow_charged_fragments:
+
+        BO_valences = list(BO.sum(axis=1))
+        for i, atom in enumerate(atoms):
+            q = get_atomic_charge(atom, atomic_valence_electrons[atom], BO_valences[i])
+            Q += q
+            if atom == 6:
+                number_of_single_bonds_to_C = list(BO[i, :]).count(1)
+                if number_of_single_bonds_to_C == 2 and BO_valences[i] == 2:
+                    Q += 1
+                    q = 2
+                if number_of_single_bonds_to_C == 3 and Q + 1 < charge:
+                    Q += 2
+                    q = 1
+
+            if q != 0:
+                q_list.append(q)
+
+    check_sum = (BO - AC).sum() == sum(DU)
+    check_charge = charge == Q
+    # check_len = len(q_list) <= abs(charge)
+
+    if check_sum and check_charge:
+        return True
+
+    return False
+
+
+def get_atomic_charge(atom, atomic_valence_electrons, BO_valence):
+    """
+    """
+
+    if atom == 1:
+        charge = 1 - BO_valence
+    elif atom == 5:
+        charge = 3 - BO_valence
+    elif atom == 15 and BO_valence == 5:
+        charge = 0
+    elif atom == 16 and BO_valence == 6:
+        charge = 0
+    else:
+        charge = atomic_valence_electrons - 8 + BO_valence
+
+    return charge
+
+
+def clean_charges(mol):
+    """
+    This hack should not be needed anymore, but is kept just in case
+
+    """
+
+    Chem.SanitizeMol(mol)
+    #rxn_smarts = ['[N+:1]=[*:2]-[C-:3]>>[N+0:1]-[*:2]=[C-0:3]',
+    #              '[N+:1]=[*:2]-[O-:3]>>[N+0:1]-[*:2]=[O-0:3]',
+    #              '[N+:1]=[*:2]-[*:3]=[*:4]-[O-:5]>>[N+0:1]-[*:2]=[*:3]-[*:4]=[O-0:5]',
+    #              '[#8:1]=[#6:2]([!-:6])[*:3]=[*:4][#6-:5]>>[*-:1][*:2]([*:6])=[*:3][*:4]=[*+0:5]',
+    #              '[O:1]=[c:2][c-:3]>>[*-:1][*:2][*+0:3]',
+    #              '[O:1]=[C:2][C-:3]>>[*-:1][*:2]=[*+0:3]']
+
+    rxn_smarts = ['[#6,#7:1]1=[#6,#7:2][#6,#7:3]=[#6,#7:4][CX3-,NX3-:5][#6,#7:6]1=[#6,#7:7]>>'
+                  '[#6,#7:1]1=[#6,#7:2][#6,#7:3]=[#6,#7:4][-0,-0:5]=[#6,#7:6]1[#6-,#7-:7]',
+                  '[#6,#7:1]1=[#6,#7:2][#6,#7:3](=[#6,#7:4])[#6,#7:5]=[#6,#7:6][CX3-,NX3-:7]1>>'
+                  '[#6,#7:1]1=[#6,#7:2][#6,#7:3]([#6-,#7-:4])=[#6,#7:5][#6,#7:6]=[-0,-0:7]1']
+
+    fragments = Chem.GetMolFrags(mol,asMols=True,sanitizeFrags=False)
+
+    for i, fragment in enumerate(fragments):
+        for smarts in rxn_smarts:
+            patt = Chem.MolFromSmarts(smarts.split(">>")[0])
+            while fragment.HasSubstructMatch(patt):
+                rxn = AllChem.ReactionFromSmarts(smarts)
+                ps = rxn.RunReactants((fragment,))
+                fragment = ps[0][0]
+                Chem.SanitizeMol(fragment)
+        if i == 0:
+            mol = fragment
+        else:
+            mol = Chem.CombineMols(mol, fragment)
+
+    return mol
+
+
+def BO2mol(mol, BO_matrix, atoms, atomic_valence_electrons,
+           mol_charge, allow_charged_fragments=True):
+    """
+    based on code written by Paolo Toscani
+
+    From bond order, atoms, valence structure and total charge, generate an
+    rdkit molecule.
+
+    args:
+        mol - rdkit molecule
+        BO_matrix - bond order matrix of molecule
+        atoms - list of integer atomic symbols
+        atomic_valence_electrons -
+        mol_charge - total charge of molecule
+
+    optional:
+        allow_charged_fragments - bool - allow charged fragments
+
+    returns
+        mol - updated rdkit molecule with bond connectivity
+
+    """
+
+    l = len(BO_matrix)
+    l2 = len(atoms)
+    BO_valences = list(BO_matrix.sum(axis=1))
+
+    if (l != l2):
+        raise RuntimeError('sizes of adjMat ({0:d}) and Atoms {1:d} differ'.format(l, l2))
+
+    rwMol = Chem.RWMol(mol)
+
+    bondTypeDict = {
+        1: Chem.BondType.SINGLE,
+        2: Chem.BondType.DOUBLE,
+        3: Chem.BondType.TRIPLE
+    }
+
+    for i in range(l):
+        for j in range(i + 1, l):
+            bo = int(round(BO_matrix[i, j]))
+            if (bo == 0):
+                continue
+            bt = bondTypeDict.get(bo, Chem.BondType.SINGLE)
+            rwMol.AddBond(i, j, bt)
+
+    mol = rwMol.GetMol()
+
+    if allow_charged_fragments:
+        mol = set_atomic_charges(
+            mol,
+            atoms,
+            atomic_valence_electrons,
+            BO_valences,
+            BO_matrix,
+            mol_charge)
+    else:
+        mol = set_atomic_radicals(mol, atoms, atomic_valence_electrons, BO_valences)
+
+    return mol
+
+
+def set_atomic_charges(mol, atoms, atomic_valence_electrons,
+                       BO_valences, BO_matrix, mol_charge):
+    """
+    """
+    q = 0
+    for i, atom in enumerate(atoms):
+        a = mol.GetAtomWithIdx(i)
+        charge = get_atomic_charge(atom, atomic_valence_electrons[atom], BO_valences[i])
+        q += charge
+        if atom == 6:
+            number_of_single_bonds_to_C = list(BO_matrix[i, :]).count(1)
+            if number_of_single_bonds_to_C == 2 and BO_valences[i] == 2:
+                q += 1
+                charge = 0
+            if number_of_single_bonds_to_C == 3 and q + 1 < mol_charge:
+                q += 2
+                charge = 1
+
+        if (abs(charge) > 0):
+            a.SetFormalCharge(int(charge))
+
+    mol = clean_charges(mol)
+
+    return mol
+
+
+def set_atomic_radicals(mol, atoms, atomic_valence_electrons, BO_valences):
+    """
+
+    The number of radical electrons = absolute atomic charge
+
+    """
+    for i, atom in enumerate(atoms):
+        a = mol.GetAtomWithIdx(i)
+        charge = get_atomic_charge(
+            atom,
+            atomic_valence_electrons[atom],
+            BO_valences[i])
+
+        if (abs(charge) > 0):
+            a.SetNumRadicalElectrons(abs(int(charge)))
+
+    return mol
+
+
+def get_bonds(UA, AC):
+    """
+
+    """
+    bonds = []
+
+    for k, i in enumerate(UA):
+        for j in UA[k + 1:]:
+            if AC[i, j] == 1:
+                bonds.append(tuple(sorted([i, j])))
+
+    return bonds
+
+
+def get_UA_pairs(UA, AC, use_graph=True):
+    """
+
+    """
+
+    bonds = get_bonds(UA, AC)
+
+    if len(bonds) == 0:
+        return [()]
+
+    if use_graph:
+        G = nx.Graph()
+        G.add_edges_from(bonds)
+        UA_pairs = [list(nx.max_weight_matching(G))]
+        return UA_pairs
+
+    max_atoms_in_combo = 0
+    UA_pairs = [()]
+    for combo in list(itertools.combinations(bonds, int(len(UA) / 2))):
+        flat_list = [item for sublist in combo for item in sublist]
+        atoms_in_combo = len(set(flat_list))
+        if atoms_in_combo > max_atoms_in_combo:
+            max_atoms_in_combo = atoms_in_combo
+            UA_pairs = [combo]
+
+        elif atoms_in_combo == max_atoms_in_combo:
+            UA_pairs.append(combo)
+
+    return UA_pairs
+
+
+def AC2BO(AC, atoms, charge, allow_charged_fragments=True, use_graph=True):
+    """
+
+    implemenation of algorithm shown in Figure 2
+
+    UA: unsaturated atoms
+
+    DU: degree of unsaturation (u matrix in Figure)
+
+    best_BO: Bcurr in Figure
+
+    """
+
+    global atomic_valence
+    global atomic_valence_electrons
+
+    # make a list of valences, e.g. for CO: [[4],[2,1]]
+    valences_list_of_lists = []
+    AC_valence = list(AC.sum(axis=1))
+    for atomicNum in atoms:
+        valences_list_of_lists.append(atomic_valence[atomicNum])
+
+    # convert [[4],[2,1]] to [[4,2],[4,1]]
+    valences_list = itertools.product(*valences_list_of_lists)
+
+    best_BO = AC.copy()
+
+    for valences in valences_list:
+
+        UA, DU_from_AC = get_UA(valences, AC_valence)
+
+        check_len = (len(UA) == 0)
+        if check_len:
+            check_bo = BO_is_OK(AC, AC, charge, DU_from_AC,
+                atomic_valence_electrons, atoms, valences,
+                allow_charged_fragments=allow_charged_fragments)
+        else:
+            check_bo = None
+
+        if check_len and check_bo:
+            return AC, atomic_valence_electrons
+
+        UA_pairs_list = get_UA_pairs(UA, AC, use_graph=use_graph)
+        for UA_pairs in UA_pairs_list:
+            BO = get_BO(AC, UA, DU_from_AC, valences, UA_pairs, use_graph=use_graph)
+            status = BO_is_OK(BO, AC, charge, DU_from_AC,
+                        atomic_valence_electrons, atoms, valences,
+                        allow_charged_fragments=allow_charged_fragments)
+
+            if status:
+                return BO, atomic_valence_electrons
+
+            elif BO.sum() >= best_BO.sum() and valences_not_too_large(BO, valences):
+                best_BO = BO.copy()
+
+    return best_BO, atomic_valence_electrons
+
+
+def AC2mol(mol, AC, atoms, charge, allow_charged_fragments=True, use_graph=True):
+    """
+    """
+
+    # convert AC matrix to bond order (BO) matrix
+    BO, atomic_valence_electrons = AC2BO(
+        AC,
+        atoms,
+        charge,
+        allow_charged_fragments=allow_charged_fragments,
+        use_graph=use_graph)
+
+    # add BO connectivity and charge info to mol object
+    mol = BO2mol(
+        mol,
+        BO,
+        atoms,
+        atomic_valence_electrons,
+        charge,
+        allow_charged_fragments=allow_charged_fragments)
+
+    return mol
+
+
+def get_proto_mol(atoms):
+    """
+    """
+    mol = Chem.MolFromSmarts("[#" + str(atoms[0]) + "]")
+    rwMol = Chem.RWMol(mol)
+    for i in range(1, len(atoms)):
+        a = Chem.Atom(atoms[i])
+        rwMol.AddAtom(a)
+
+    mol = rwMol.GetMol()
+
+    return mol
+
+
+def read_xyz_file(filename, look_for_charge=True):
+    """
+    """
+
+    atomic_symbols = []
+    xyz_coordinates = []
+    charge = 0
+    title = ""
+
+    with open(filename, "r") as file:
+        for line_number, line in enumerate(file):
+            if line_number == 0:
+                num_atoms = int(line)
+            elif line_number == 1:
+                title = line
+                if "charge=" in line: #TODO: More flexible charge declaration
+                    charge = int(line.split("=")[1])
+            else:
+                atomic_symbol, x, y, z = line.split()
+                atomic_symbols.append(atomic_symbol)
+                xyz_coordinates.append([float(x), float(y), float(z)])
+
+    atoms = [int_atom(atom) for atom in atomic_symbols]
+
+    return atoms, charge, xyz_coordinates
+
+
+def xyz2AC(atoms, xyz, charge, use_huckel=False):
+    """
+
+    atoms and coordinates to atom connectivity (AC)
+
+    args:
+        atoms - int atom types
+        xyz - coordinates
+        charge - molecule charge
+
+    optional:
+        use_huckel - Use Huckel method for atom connecitivty
+
+    returns
+        ac - atom connectivity matrix
+        mol - rdkit molecule
+
+    """
+
+    if use_huckel:
+        return xyz2AC_huckel(atoms, xyz, charge)
+    else:
+        return xyz2AC_vdW(atoms, xyz)
+
+
+def xyz2AC_vdW(atoms, xyz):
+
+    # Get mol template
+    mol = get_proto_mol(atoms)
+
+    # Set coordinates
+    conf = Chem.Conformer(mol.GetNumAtoms())
+    for i in range(mol.GetNumAtoms()):
+        conf.SetAtomPosition(i, (xyz[i][0], xyz[i][1], xyz[i][2]))
+    mol.AddConformer(conf)
+
+    AC = get_AC(mol)
+
+    return AC, mol
+
+
+def get_AC(mol, covalent_factor=1.3):
+    """
+
+    Generate adjacent matrix from atoms and coordinates.
+
+    AC is a (num_atoms, num_atoms) matrix with 1 being covalent bond and 0 is not
+
+
+    covalent_factor - 1.3 is an arbitrary factor
+
+    args:
+        mol - rdkit molobj with 3D conformer
+
+    optional
+        covalent_factor - increase covalent bond length threshold with facto
+
+    returns:
+        AC - adjacent matrix
+
+    """
+
+    # Calculate distance matrix
+    dMat = Chem.Get3DDistanceMatrix(mol)
+
+    pt = Chem.GetPeriodicTable()
+    num_atoms = mol.GetNumAtoms()
+    AC = np.zeros((num_atoms, num_atoms), dtype=int)
+
+    for i in range(num_atoms):
+        a_i = mol.GetAtomWithIdx(i)
+        Rcov_i = pt.GetRcovalent(a_i.GetAtomicNum()) * covalent_factor
+        for j in range(i + 1, num_atoms):
+            a_j = mol.GetAtomWithIdx(j)
+            Rcov_j = pt.GetRcovalent(a_j.GetAtomicNum()) * covalent_factor
+            if dMat[i, j] <= Rcov_i + Rcov_j:
+                AC[i, j] = 1
+                AC[j, i] = 1
+
+    return AC
+
+
+def xyz2AC_huckel(atomicNumList,xyz,charge):
+    """
+
+    args
+        atomicNumList - atom type list
+        xyz - coordinates
+        charge - molecule charge
+
+    returns
+        ac - atom connectivity
+        mol - rdkit molecule
+
+    """
+    mol = get_proto_mol(atomicNumList)
+
+    conf = Chem.Conformer(mol.GetNumAtoms())
+    for i in range(mol.GetNumAtoms()):
+        conf.SetAtomPosition(i,(xyz[i][0],xyz[i][1],xyz[i][2]))
+    mol.AddConformer(conf)
+
+    num_atoms = len(atomicNumList)
+    AC = np.zeros((num_atoms,num_atoms)).astype(int)
+
+    mol_huckel = Chem.Mol(mol)
+    mol_huckel.GetAtomWithIdx(0).SetFormalCharge(charge) #mol charge arbitrarily added to 1st atom    
+
+    passed,result = rdEHTTools.RunMol(mol_huckel)
+    opop = result.GetReducedOverlapPopulationMatrix()
+    tri = np.zeros((num_atoms, num_atoms))
+    tri[np.tril(np.ones((num_atoms, num_atoms), dtype=bool))] = opop #lower triangular to square matrix
+    for i in range(num_atoms):
+        for j in range(i+1,num_atoms):
+            pair_pop = abs(tri[j,i])   
+            if pair_pop >= 0.15: #arbitry cutoff for bond. May need adjustment
+                AC[i,j] = 1
+                AC[j,i] = 1
+
+    return AC, mol
+
+
+def chiral_stereo_check(mol):
+    """
+    Find and embed chiral information into the model based on the coordinates
+
+    args:
+        mol - rdkit molecule, with embeded conformer
+
+    """
+    Chem.SanitizeMol(mol)
+    Chem.DetectBondStereochemistry(mol, -1)
+    Chem.AssignStereochemistry(mol, flagPossibleStereoCenters=True, force=True)
+    Chem.AssignAtomChiralTagsFromStructure(mol, -1)
+
+    return
+
+
+def xyz2mol(atoms, coordinates,
+    charge=0,
+    allow_charged_fragments=True,
+    use_graph=True,
+    use_huckel=False,
+    embed_chiral=True):
+    """
+    Generate a rdkit molobj from atoms, coordinates and a total_charge.
+
+    args:
+        atoms - list of atom types (int)
+        coordinates - 3xN Cartesian coordinates
+        charge - total charge of the system (default: 0)
+
+    optional:
+        allow_charged_fragments - alternatively radicals are made
+        use_graph - use graph (networkx)
+        use_huckel - Use Huckel method for atom connectivity prediction
+        embed_chiral - embed chiral information to the molecule
+
+    returns:
+        mol - rdkit molobj
+
+    """
+
+    # Get atom connectivity (AC) matrix, list of atomic numbers, molecular charge,
+    # and mol object with no connectivity information
+    print("atoms variable: ", atoms)
+    AC, mol = xyz2AC(atoms, coordinates, charge, use_huckel=use_huckel)
+
+    # Convert AC to bond order matrix and add connectivity and charge info to
+    # mol object
+    new_mol = AC2mol(mol, AC, atoms, charge,
+        allow_charged_fragments=allow_charged_fragments,
+        use_graph=use_graph)
+
+    # Check for stereocenters and chiral centers
+    if embed_chiral:
+        chiral_stereo_check(new_mol)
+
+    return new_mol
+
+
+def main():
+
+
+    return
+
+
+if __name__ == "__main__":
+
+    import argparse
+
+    parser = argparse.ArgumentParser() # usage='%(prog)s [options] molecule.xyz')
+    parser.add_argument('structure', metavar='structure', type=str)
+    parser.add_argument('-s', '--sdf',
+        action="store_true",
+        help="Dump sdf file")
+    parser.add_argument('--ignore-chiral',
+        action="store_true",
+        help="Ignore chiral centers")
+    parser.add_argument('--no-charged-fragments',
+        action="store_true",
+        help="Allow radicals to be made")
+    parser.add_argument('--no-graph',
+        action="store_true",
+        help="Run xyz2mol without networkx dependencies")
+
+    # huckel uses extended Huckel bond orders to locate bonds (requires RDKit 2019.9.1 or later)
+    # otherwise van der Waals radii are used
+    parser.add_argument('--use-huckel',
+        action="store_true",
+        help="Use Huckel method for atom connectivity")
+    parser.add_argument('-o', '--output-format',
+        action="store",
+        type=str,
+        help="Output format [smiles,sdf] (default=sdf)")
+    parser.add_argument('-c', '--charge',
+        action="store",
+        metavar="int",
+        type=int,
+        help="Total charge of the system")
+
+    args = parser.parse_args()
+
+    # read xyz file
+    filename = args.structure
+
+    # allow for charged fragments, alternatively radicals are made
+    charged_fragments = not args.no_charged_fragments
+
+    # quick is faster for large systems but requires networkx
+    # if you don't want to install networkx set quick=False and
+    # uncomment 'import networkx as nx' at the top of the file
+    quick = not args.no_graph
+
+    # chiral comment
+    embed_chiral = not args.ignore_chiral
+
+    # read atoms and coordinates. Try to find the charge
+    atoms, charge, xyz_coordinates = read_xyz_file(filename)
+
+    # huckel uses extended Huckel bond orders to locate bonds (requires RDKit 2019.9.1 or later)
+    # otherwise van der Waals radii are used
+    use_huckel = args.use_huckel
+
+    # if explicit charge from args, set it
+    if args.charge is not None:
+        charge = int(args.charge)
+
+    # Get the molobj
+    mol = xyz2mol(atoms, xyz_coordinates,
+        charge=charge,
+        use_graph=quick,
+        allow_charged_fragments=charged_fragments,
+        embed_chiral=embed_chiral,
+        use_huckel=use_huckel)
+
+    # Print output
+    if args.output_format == "sdf":
+        txt = Chem.MolToMolBlock(mol)
+        print(txt)
+
+    else:
+        # Canonical hack
+        isomeric_smiles = not args.ignore_chiral
+        smiles = Chem.MolToSmiles(mol, isomericSmiles=isomeric_smiles)
+        m = Chem.MolFromSmiles(smiles)
+        smiles = Chem.MolToSmiles(m, isomericSmiles=isomeric_smiles)
+        print(smiles)
