def do_for_each(name, num_threads):
    _format = name[name.rfind(".")+1:]
    basename = name[:name.rfind(".")]
    out_name = basename + ".num"
    pdb_id = basename.split("_")[-1]
    pdb_name = "/home/xmluo/jlpeng/cMet/pdb_protein/%s_protein.pdb"%pdb_id
    global ligand,protein,outf,ap_iter
    ligand = pybel.readfile(_format, name).next()
    protein = pybel.readfile("pdb", pdb_name).next()
    outf = open(out_name, "w")
    #make_pairs()
    #print "after creating pair_queue"
    ap_iter = AtomPairs()
    threads = []
    for i in xrange(num_threads):
        threads.append(AtomPairWalker())
    print "len(threads)=%d"%len(threads)
    for t in threads:
        t.start()
    for t in threads:
        t.join()
        if t.isAlive():
            print "thread %s is alive"%t.getName()
        else:
            print "thread %s is dead"%t.getName()
    outf.close()

def mp_helper(file_one, file_two, iteration_depth = 1):
    """
        Helper function for the multiprocessing library.
        Two fragment files gets passed and we merge all against all molecules in that two files.
    """
    results = list()
    fragments = list()
    for mol_one in pybel.readfile( 'smi', file_one ):
        for i,mol_two in enumerate(pybel.readfile( 'smi', file_two )):
            #print 'merge:', str(mol_one).strip(), str(mol_two).strip()
            result, fragment = merge(mol_two,mol_one, options, iteration_depth)
            if result:
                #print '\tr', result
                results.extend( result )
            if fragment:
                #print '\tf', fragment
                fragments.extend( fragment )

    fragment_return, molecule_return = None, None
    if fragments:
        fragment_file = tempfile.NamedTemporaryFile(dir=temp_dir, delete=False)
        fragment_file.writelines( fragments )
        fragment_file.close()
        fragment_return = fragment_file.name
    if results:
        result_file = tempfile.NamedTemporaryFile(dir=temp_dir, delete=False)
        result_file.writelines( results )
        result_file.close()
        molecule_return = result_file.name
    return molecule_return, fragment_return

def dock_ligands_to_receptors(docking_dir, worker_pool=False, exhaustiveness=None, chosen_receptor=None, restrict_box=True):
  subdirs = glob.glob(os.path.join(docking_dir, '*/'))
  for subdir in subdirs:
    subdir = subdir.rstrip('/')
    receptor_name = os.path.basename(subdir)
    if chosen_receptor is not None and chosen_receptor != receptor_name:
        continue
    print("receptor name = %s" % receptor_name)
    receptor_filename = os.path.join(subdir, "%s.pdbqt" % receptor_name)
    if not os.path.exists(receptor_filename):
      continue

    print("Examining %s" % receptor_filename)

    receptor_pybel = next(pybel.readfile("pdb", 
        os.path.join(subdir, "%s.pdb" % receptor_name)))
    protein_centroid, protein_range = get_molecule_data(receptor_pybel)

    box_dims = protein_range + 5.0

    ligands = sorted(glob.glob(os.path.join(subdir, '*_prepared.pdbqt')))
    print("Num ligands = %d" % len(ligands))


    dock_ligand_to_receptor_partial = partial(dock_ligand_to_receptor, receptor_filename=receptor_filename,
                                              protein_centroid=protein_centroid, box_dims=box_dims,
                                              subdir=subdir, exhaustiveness=exhaustiveness)

    if restrict_box:
      active_ligand = ""
      for ligand in ligands:
        if "CHEM" in ligand:
          active_ligand = ligand
          break

      print("Docking to %s first to ascertain centroid and box dimensions" % active_ligand)

      out_pdb_qt = dock_ligand_to_receptor_partial(active_ligand)
      ligand_pybel = next(pybel.readfile("pdbqt", 
                                    out_pdb_qt))
      ligand_centroid, _ = get_molecule_data(ligand_pybel)
      print("Protein centroid = %s" %(str(protein_centroid)))
      print("Ligand centroid = %s" %(str(ligand_centroid)))
      box_dims = np.array([20., 20., 20.])
      dock_ligand_to_receptor_partial = partial(dock_ligand_to_receptor, receptor_filename=receptor_filename,
                                          protein_centroid=ligand_centroid, box_dims=box_dims,
                                          subdir=subdir, exhaustiveness=exhaustiveness)

      print("Finished docking to %s, docking to remainder of ligands now." % active_ligand)

    if worker_pool is False:
      for i, ligand_file in enumerate(ligands):
        a = time.time()
        dock_ligand_to_receptor_partial(ligand)
        print("took %f seconds to dock single ligand." %(time.time() - a))
    else:
      print("parallelizing docking over worker pool")

      worker_pool.map_sync(dock_ligand_to_receptor_partial, ligands)

    def helper_run(self, result):
        result_pdbqt = result.output().path
        result_lig = pybel.readfile("pdbqt", result_pdbqt).next()
        native_lig = pybel.readfile("pdbqt",
                                    LigPdbqt(self.tname).output().path).next()
        result_rmsd = rmsd_between(native_lig, result_lig)

        return result_rmsd

def create_inchi(job, chemid, xyz_file=''):
    if xyz_file == '':
        xyz_file = os.path.expanduser(job) + 'xyz/' + chemid + '.xyz'
    obmol = list(pybel.readfile('xyz', xyz_file))[0]
    try:
        obmol = list(pybel.readfile('xyz', xyz_file))[0]
    except NameError:
        logging.error('Pybel is not installed or loaded correctly.')
        sys.exit()

    return obmol.write("inchi", opt={'T': 'nostereo'}).split()[0]

def prepare_dsstox_dataset(root, name, dest=None, overwrite=False):
    """ This method bootstraps the analysis of DSSTox data.
       - Rename the compounds
       - Merge train/test
       - Generate 3D conformations
       - Save "master" and "saliviewer" tables
       - Redirects stdout/stderr to a "prepare.log" file
    """
    if not dest: dest = root

    dataset_root = op.join(dest, name)
    dest_sdf = op.join(dataset_root, name + '.sdf')
    if op.exists(dest_sdf) and not overwrite:
        print '%s is already there and not overwriting requested' % dest_sdf
        return

    print 'Reading %s' % name
    train_mols = list(pybel.readfile('sdf', op.join(root, name + '_training.sdf')))
    test_mols = list(pybel.readfile('sdf', op.join(root, name + '_prediction.sdf')))

    print '\tCreating dataset root: %s' % dataset_root
    if not op.exists(dataset_root):
        os.makedirs(dataset_root)

    print '\tRenaming the compounds to keep track of the provenance'
    rename_mols_by_index(train_mols, name + '-train-')
    rename_mols_by_index(test_mols, name + '-test-')

    print '\tGenerating conformations'
    for mol in train_mols + test_mols:
        #Some molecules from mutagenicity produce segfault on make3D
        #See bug report at https://sourceforge.net/tracker/?func=detail&aid=3374324&group_id=40728&atid=428740
        #Train 3988: OC(=O)[[email protected]]1(C)CCC[[email protected]]2(C1CC[[email protected]]13C2CC[[email protected]](C3)([[email protected]]2(C1)OC2)O)C
        #Train 4205: CC(CCC[[email protected]]([[email protected]]1CC[[email protected]@H]2[[email protected]]1(C)CC[[email protected]]1[[email protected]]2CC2([[email protected]@H]3[[email protected]]1(C)CC[[email protected]@H](C3)Br)S(=O)(=O)CCS2(=O)=O)C)C
        #This kind of fatal errors are worrying, is there any robust way of controlling them in python/java? Will need to create one
        if not any(name in mol.title for name in ('train-3988', 'train-4205')):
            try:
                print 'Conformation for %s' % mol.title
                mol.make3D()
            except Exception:
                print 'Error computing a 3D conformation for %s' % mol.title

    print '\tSaving compounds'
    save_mols(train_mols + test_mols, dest_sdf)

    master_table = op.join(dataset_root, name + '-master.csv')
    print '\tCreating \"master\" table: %s' % master_table
    create_master_table(dest_sdf, master_table)

    sali_table = op.join(dataset_root, name + '-saliviewer.csv')
    print '\tCreating \"saliviewer\" table: %s' % sali_table
    create_saliviewer_input(master_table, sali_table)

 def __init__(self, filename, ext):
     self.data = {'name':'','index':'', 'metal':'', 'topology':'', 'parent':'',
             'atomic_info':'', 'bond_table':'', 'connectivity':'',
             'connect_flag':'', 'connect_sym':''}
     name = os.path.split(filename)[-1]
     self.name = clean(name, ext) 
     self.update(name=self.name)
     # may be a source of error.. untested
     if version_info.major >= 3:
         self.mol = next(pybel.readfile(ext, filename))
     else:
         self.mol = pybel.readfile(ext, filename).next()
     self._reset_formal_charges()

    def pred_foo(tname):
        try:
            geaux_sdf = os.path.join(GEAUX_OUTPUT, tname,
                                     tname + '_pred_0.sdf')
            native_sdf = VinaPath(tname).lig_sdf

            geaux_lig = pybel.readfile("sdf", geaux_sdf).next()
            native_lig = pybel.readfile("sdf", native_sdf).next()
            result_rmsd = rmsd_between(native_lig, geaux_lig)
            return result_rmsd
        except Exception as e:
            print(e)
            return None

    def caculateRMSD(self):
        vina_task = self.requires()
        predicted_pdbqt = vina_task.output().path
        predicted_mol = pybel.readfile('pdbqt', predicted_pdbqt).next()
        crystal_pdbqt = vina_task.lig_pdbqt
        crystal_mol = pybel.readfile('pdbqt', crystal_pdbqt).next()

        def rmsd(m1, m2):
            c1 = [a.coords for a in m1 if not a.OBAtom.IsHydrogen()]
            c2 = [a.coords for a in m2 if not a.OBAtom.IsHydrogen()]
            return dockedpose.rmsd(c1, c2)

        return rmsd(predicted_mol, crystal_mol)

 def helper_geauxdock(tname, version="0.7"):
     try:
         native_sdf = VinaPath(tname).lig_sdf
         geaux_sdf = os.path.join(
             ModelPath(tname,
                       version=version).work_dir,
             "{}_{}_0.sdf".format(tname, version))
         geaux_lig = pybel.readfile("sdf", geaux_sdf).next()
         native_lig = pybel.readfile("sdf", native_sdf).next()
         result_rmsd = rmsd_between(native_lig, geaux_lig)
         return result_rmsd
     except Exception as e:
         print(e)
         return None

 def run(self):
     data = {}
     for sdf_id in self.getSdfs():
         path = Path(sdf_id)
         lig = pybel.readfile('sdf', path.astex_sdf()).next()
         lig.removeh()
         prt = pybel.readfile('pdb', path.astex_pdb()).next()
         prt.removeh()
         lig_sz = len(lig.atoms)
         prt_sz = len(prt.atoms)
         data[sdf_id] = {"lig_sz": lig_sz,
                         "prt_sz": prt_sz}
     dset = pd.DataFrame(data)
     dset.to_csv(self.output().path)

def main():
    file_list=GetFileList()
    dc_fn='dist_splited_train_0525.dat'
    dc=dist.UnformatedInput(dc_fn)
    pts=dist.PTS(dc)
    for i in range(len(file_list)):
        lig_fn=file_list[i][0]
        pro_fn=file_list[i][1]
        pdbid=file_list[i][2]
        if os.path.isfile(lig_fn):
            lig=pybel.readfile('sdf',lig_fn).next()
            pro=pybel.readfile('pdb',pro_fn).next()
            score,scorev=GetFinger(lig,pro,pts,pdbid)
            print score,scorev

def main( args ):

    for infile in args.infile:
        file_extension = args.format or os.path.splitext( infile )[-1].lstrip('.')

        if not args.format and file_extension not in ['smi', 'sdf', 'inchi', 'mol']:
            sys.exit('Could not guess the format from the file extension please specify with the --format option.')

        molecules = pybel.readfile(file_extension, infile)
        for mol in molecules:
            args.outfile.write( 't # id %s\n' % mol.title.strip() )
            for atom in openbabel.OBMolAtomIter( mol.OBMol):
                label = atom.GetAtomicNum()
                vertex_index = atom.GetIdx()
                args.outfile.write('v %s %s\n' % (vertex_index, label))

            for bond in openbabel.OBMolBondIter( mol.OBMol):
                src_index = bond.GetBeginAtomIdx()
                dest_index = bond.GetEndAtomIdx()
                assert(src_index > 0)
                assert(dest_index > 0)
                if bond.IsAromatic():
                    label = 'a'
                elif bond.IsSingle():
                    label = 's'
                elif bond.IsDouble():
                    label = 'd'
                elif bond.IsTriple():
                    label = 't'
                atom1 = bond.GetBeginAtom()
                atom2 = bond.GetEndAtom()
                args.outfile.write('e %s %s %s\n' % (src_index, dest_index, label))

def run():
	
	inputfile=pybel.readfile(sys.argv[1].split(".")[-1],sys.argv[1])
	value=()
	for mol in inputfile:
		descvalues=mol.calcdesc()
		value= value+(descvalues.get('TPSA'),)
		value= value+(descvalues.get('HBD'),)
		value= value+(descvalues.get('logP'),)
		value= value+(descvalues.get('MW'),)
		value= value+(descvalues.get('tbonds'),)
		value= value+(descvalues.get('nF'),)
		value= value+(descvalues.get('bonds'),)
		value= value+(descvalues.get('atoms'),)
		value= value+(descvalues.get('HBA1'),)
		value= value+(descvalues.get('HBA2'),)
		value= value+(descvalues.get('sbonds'),)
		value= value+(descvalues.get('dbonds'),)
		value= value+(descvalues.get('MR'),)
		value= value+(descvalues.get('abonds'),)
	
		smarts = pybel.Smarts("[+]")
		num=smarts.findall(mol)				
		value= value+(len(num),)			
	
		smarts = pybel.Smarts("[-]")
		num=smarts.findall(mol)				
		value= value+(len(num),)


	model=joblib.load('volume_model/volume.pkl')
	for result in model.predict(value):
		return round(result,2)

def main():
    fa = sys.stdin.read()
    fa = fa.rstrip()
    of = open(args['outfile'], 'w')
    of.write('cid,')
    labels = ''
    for i in properties:
        labels = labels + i + ','
    labels = re.match(r"^(.*),", labels).group(1)
    of.write(labels + '\n')
    inputTemp = tempfile.NamedTemporaryFile(suffix='.sdf', delete=False)
    inputTemp.write(fa)
    inputTempName = inputTemp.name
    inputTemp.close()
    for mol in pybel.readfile('sdf', inputTempName):
        info = ''
        myid = mol.title
        mol.addh()
        mol.make3D()
        desc = mol.calcdesc()
        for thisdesc in properties:
            info = info + str(desc[thisdesc]) + ','
        info = re.match(r"^(.*),", info).group(1)
        of.write(myid.strip() + ',' + info + '\n')
    of.close()
    os.unlink(inputTempName)

def read_file(filename, name=None, format=None):
    """ Read a molecule from a file

    Note:
        Currently only reads the first conformation in a file

    Args:
        filename (str): path to file
        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 main():
    fa = sys.stdin.read()
    fa = fa.rstrip()
    of = open(args['outfile'], 'w')
    of.write('cid\t')
    labels = 'InChI\tInChIkey'
    of.write(labels + '\n')
    inputTemp = tempfile.NamedTemporaryFile(suffix='.sdf', delete=False)
    inputTemp.write(fa)
    inputTempName = inputTemp.name
    inputTemp.close()
    for mol in pybel.readfile('sdf', inputTempName):
        info = ''
        myid = mol.title
        mol.addh()
        mol.make3D()
        conv = ob.OBConversion()
        conv.SetInAndOutFormats("sdf", "inchi")
        inchi = conv.WriteString(mol.OBMol)
        conv.SetInAndOutFormats("sdf", "inchikey")
        inchikey = conv.WriteString(mol.OBMol)
        info = re.match(r"^InChI=(.*)\n", inchi).group(1)
        of.write(myid.strip() + '\t' + info + '\t' + inchikey + '\n')
    of.close()
    os.unlink(inputTempName)

def step(context, action=None, projkey=None):
    # something here
    from cbh_core_model.models import Project, CBHCompoundBatch
    from rdkit import Chem
    from rdkit.Chem import AllChem, inchi

    path = "/dev/cbh_compound_batches/"
    resp = context.api_client.get(
        path,
        format='json',
        data=context.post_data,
    )

    reg_cmpds = context.ser.deserialize(resp.content)["objects"]
    reg_inchis = []
    # get a list of inchis from the reponse
    for cmpd in reg_cmpds:
        reg_inchis.append(cmpd['standardInchi'].strip())

    fn = os.path.join(os.path.dirname(__file__), 'files/inchi-list.txt')
    inchis = [mol.write("inchi").split("\t")[0].strip()
              for mol in readfile('inchi', fn)]

    # do an array subtraction of the hardcoded inchis from the registered inchis
    # print(set(inchis))
    print(len(inchis))
    # print(set(reg_inchis))
    print(len(reg_inchis))
    diff = list(set(inchis) - set(reg_inchis))
    print(len(diff))
    # print(diff)
    assert len(diff) == 0

def readfile(format, filename, opt=None, lazy=False):
    if lazy and format == 'mol2':
        return _filereader_mol2(filename, opt=opt)
    elif lazy and format == 'sdf':
        return _filereader_sdf(filename, opt=opt)
    else:
        return pybel.readfile(format, filename, opt=opt)

def main(input_ext, inputfile, output_ext, outputfilename,
         nconfs, rmsd_cutoff, energy_cutoff):
    
    ff = pybel._forcefields['mmff94']
    outputfile = pybel.Outputfile(output_ext, outputfilename, overwrite=True)
    for i, mol in enumerate(pybel.readfile(input_ext, inputfile)):
        t = time.time()
        
        print "**Molecule %d\n..title = %s" % (i, mol.title)
        print "..number of rotatable bonds = %d" % mol.OBMol.NumRotors()
        mol.addh()
        ff.Setup(mol.OBMol)
        ff.DiverseConfGen(rmsd_cutoff, nconfs, energy_cutoff)

        ff.GetConformers(mol.OBMol)
        confdata = pybel.ob.toConformerData(mol.OBMol.GetData(pybel.ob.ConformerData))
        energies = confdata.GetEnergies()

        N = mol.OBMol.NumConformers()
        assert N == len(energies)
        print "..generated %d conformers"
        
        u = time.time()
        data = []
        for i in range(N):
            mol.OBMol.SetConformer(i)
            outputfile.write(mol)

        print "..(overall time = %.1fs  writing results = %.1fs)" % (time.time() - t,
                                                                  time.time() -u)
        print "\n"
        
    outputfile.close()