How to use the pybel.readfile function in pybel

def read_mol(mol, fmt="mol2"):
    """Read a molecular object either from a file or string"""
    if not os.path.isfile(mol):
        name = '/tmp/structure.{}'.format(fmt)
        with open(name, 'w') as f:
            f.write(mol)
        mol = name

    return pybel.readfile(fmt, mol).next()

name (str): name to assign to molecule
        format (str): File format: pdb, sdf, mol2, bbll, etc.

    Returns:
        moldesign.Molecule: parsed result
    """
    # TODO: check for openbabel molecule name?
    if format is None:
        format = filename.split('.')[-1]

    if force_remote:
        with open(filename, 'r') as infile:
            mol = read_string(infile.read(), format, name=name)
        return mol
    else:
        pbmol = pb.readfile(format=format, filename=filename).next()
        if name is None: name = filename
        mol = pybel_to_mol(pbmol, name=os.path.basename(name))
        mol.filename = filename
        return mol

def standardizeMol2(path, reference, output):
    traj = next(pybel.readfile("mol2", path))
    ref = next(pybel.readfile("mol2", reference))
    for iatom in traj.atoms:
        ob = iatom.OBAtom
        idx = ob.GetIndex()

        jatom = ref.OBMol.GetAtomById(idx)
        ob.SetType(jatom.GetType())
        ob.SetAtomicNum(jatom.GetAtomicNum())

    traj.write("xyz", output, True)

# select only those that have a hydrogen atom attached 
                    
                    min_num_vs = numpy.min(num_vs)
                    self.min_acidic_pKa.append(min_num_vs) # lowest value
                    # get index corresponding to the lowest acidic pKa in 
                    index_acidic_atom = numpy.argmin(num_vs) # get the index of the minimum value     
                    alternate_index_acidic_atom = numpy.argmax(num_vs)
                    msg = "Lowest acidic pKa = " + str(min_num_vs)
                    logger.info(msg)
                    
                    
                    atoms_acidpka = atom_nums[0:nacidpka]
                    i1 = numpy.int(atoms_acidpka[index_acidic_atom])
                    
                    mol = pybel.readfile("mol",mfile).next()
          
                    iatom = mol.atoms[i1-1]
                    nbatoms = openbabel.OBAtomAtomIter(iatom.OBAtom)
                    
                    found = 0
                    for nb in nbatoms:
                        nb_atomnum = nb.GetAtomicNum()
                        if(nb_atomnum==1):
                            found = 1
                    if (found == 0):
                        msg = "Atom number with lowest acidic pKa has no hydrogen atom attached"
                        logger.info(msg)
                        msg = "Checking the other site"
                        logger.info(msg)
                        
                        i1 = numpy.int(atoms_acidpka[alternate_index_acidic_atom])

ppdb = 'tmp/p.pdb'
    __, intype = os.path.splitext(inprot)
    if intype[1:].lower() != 'pdb':
        prot = pybel.readfile(intype[1:], inprot).__next__()
        output = pybel.Outputfile("pdb", ppdb, overwrite=True)
        output.write(prot)
        output.close()
    else:
        # change possible HETATM to ATOM in pdb
        os.system("""sed 's/HETATM/ATOM\ \ /g' """ + inprot + " > " + ppdb)

    # convert ligand file to be PDB by openbabel
    lpdb = 'tmp/l.pdb'
    __, intype = os.path.splitext(inlig)
    if intype[1:].lower() != 'pdb':
        lig = pybel.readfile(intype[1:], inlig).__next__()
        output = pybel.Outputfile("pdb", lpdb, overwrite=True)
        output.write(lig)
        output.close()
    else:
        # change possible HETATM to ATOM in pdb
        os.system("""sed 's/HETATM/ATOM\ \ /g' """ + inlig + " > " + lpdb)

    os.chdir('tmp')
    #copy atom typefiel into directory
    #os.system("cp $DXGB/atmtypenumbers .")
    # Process p.pdb/l.pdb to be p_sa.pdb/l_sa.pdb after pharma assignment
    # get full atom idx list and pharma
    ppdb2 = 'p_sa.pdb'
    lpdb2 = 'l_sa.pdb'

    pidx, ppharm = pharma('p.pdb').assign(write=True, outfn = ppdb2)

def _load_pybel(self, input, input_type, **kwargs):
        """
        The internal function to load a molecule using rdkit engine.

        """
        # available variables
        pybel_mol = None
        creator = None

        if input_type == 'xyz':
            if os.path.isfile(input):
                creator = ('XYZ', input)
                gen = pybel.readfile("xyz", input)
                mols = list(gen)
                if len(mols) == 1:
                    pybel_mol = mols[0]
                else:
                    warnings.warn('More than one Molecule object is created and is returned as a generator.')
                    return self._multiple_molecules(mols, creator)
            else:
                msg = "The input '%s' is not a valid XYZ input file."%input
                raise ValueError(msg)

        if pybel_mol is None:
            msg = 'The input is not a legit %s string. Refere to the error message that was already displayed by Pybel library!' % \
              creator[0]
            raise ValueError(msg)

        # if the molecule is already being created warn the user

def run_pdbqt_python(pdb_path):
    import pybel
    mol = next(pybel.readfile("pdb", pdb_path))

    mol.addh()

    pdbqt = mol.write("pdbqt")
    autodock_features = {}
    for atom_index in range(mol.OBMol.NumAtoms()):
        a = mol.OBMol.GetAtom(atom_index + 1)

        if a.IsCarbon() and a.IsAromatic():
            element = 'A'
        elif a.IsOxygen():
            element = 'OA'
        elif a.IsNitrogen() and a.IsHbondAcceptor():
            element = 'NA'
        elif a.IsSulfur() and a.IsHbondAcceptor():
            element ='SA'

def rmsd(mol1, mol2):
    a = next(pybel.readfile("xyz", mol1))
    b = next(pybel.readfile("xyz", mol2))

    align = openbabel.OBAlign(False, True)

    align.SetRefMol(a.OBMol)
    align.SetTargetMol(b.OBMol)
    align.Align()
    return align.GetRMSD()

import numpy.linalg as la
    from ase.calculators import mopac
    from ase import Atoms
    import pybel

    os.chdir('../test')

    MOPAC = os.path.join(os.getcwd(), 'MOPAC')
    os.environ['MOPAC_LICENSE'] = MOPAC
    os.environ['LD_LIBRARY_PATH'] = MOPAC

    mopac_calc = mopac.MOPAC()
    mopac_calc.command = 'MOPAC/MOPAC2016.exe PREFIX.mop 2> /dev/null'
    mopac_calc.set(method='pm3')

    mol = next(pybel.readfile('xyz', 'ts2.xyz'))
    atoms = Atoms(numbers=[a.atomicnum for a in mol.atoms],
                  positions=[a.coords for a in mol.atoms])
    atoms.set_positions(atoms.positions - atoms.get_center_of_mass())
    atoms.set_calculator(mopac_calc)

    irccalc = IRC(atoms, stride=0.15, mw=True, forward=True, trajectory='ts1.traj')
    for _ in irccalc.run():
        pass

def read_mol(path, fmt='mol2'):
    return Mol(next(pybel.readfile(fmt, path)))