def align_template_to_reference(msmseed, ref_msmseed):
    import modeller
    import tempfile
    import shutil
    import copy
    import os
    temp_dir = tempfile.mkdtemp()
    try:
        os.chdir(temp_dir)
        alignment_file = open('aln_tmp.pir','w')
        aln = _PIR_alignment(ref_msmseed.template_sequence, ref_msmseed.template_id, msmseed.template_sequence, msmseed.template_id)
        alignment_file.writelines(aln)
        alignment_file.close()
        template_file = open(msmseed.template_id + '.pdb','w')
        template_pdb = msmseed.template_structure
        template_pdb.writeFile(template_pdb.topology, template_pdb.positions, template_file)
        template_file.close()
        ref_pdb = ref_msmseed.template_structure
        ref_file = open(ref_msmseed.template_id + '.pdb', 'w')
        ref_pdb.writeFile(ref_pdb.topology, ref_pdb.positions, ref_file)
        ref_file.close()
        modeller.log.none()
        env = modeller.environ()
        env.io.atom_files_directory = temp_dir
        aln = modeller.alignment(env, file='aln_tmp.pir', align_codes=(ref_msmseed.template_id, msmseed.template_id))
        mdl  = modeller.model(env, file=ref_msmseed.template_id + '.pdb')
        mdl2 = modeller.model(env, file=msmseed.template_id+'.pdb')
        atmsel = modeller.selection(mdl).only_atom_types('CA')
        r = atmsel.superpose(mdl2, aln)
        msmseed.rmsd_to_reference = copy.deepcopy(r.rms)
    except Exception as e:
        msmseed.error_message = e.message
    finally:
        shutil.rmtree(temp_dir)
    return msmseed

def get_sas(pdb, probe):
    import modeller

    # Read the PDB file
    env = modeller.environ()
    mdl = modeller.model(env)
    mdl.read(file=pdb)

    # Calculate atomic accessibilities (in Biso) with appropriate probe_radius
    myedat = modeller.energy_data()
    myedat.radii_factor = 1.6
    mdl.write_data(edat=myedat, output="PSA ATOMIC_SOL", psa_integration_step=0.05, probe_radius=probe)

    mdl.write(file=pdb.rsplit(".", 1)[0] + ".sas")

    # read SAS
    with open("%s.sas" % (pdb.rsplit(".", 1)[0],)) as data:
        D = data.readlines()

    Sas = {}
    for d in D:
        d = d.strip()
        if d[:4] == "ATOM":
            atom, res, resid, cid = d[12:16], d[17:20], int(d[22:26]), d[21]
            if cid == " ":
                cid = "A"
            Sas[(atom, res, resid, cid)] = float(d[60:66])
    return Sas

 def test_script9(self):
     """Test step 9 (multiple fitting)"""
     # Get inputs (outputs from step 8)
     for i in ('top', 'bottom'):
         shutil.copy('precalculate_results/stage8_split_density/' \
                     'groel-11.5A.%s.mrc' % i, 'output')
     # Make sure the script runs without errors
     p = subprocess.check_call(['scripts/' \
                                'script9_symmetric_multiple_fitting.py'])
     e = modeller.environ()
     ref = modeller.model(e,
            file='precalculate_results/stage9_symmetric_multiple_fitting/' \
                 'model.top.0.pdb')
     sel = modeller.selection(ref).only_atom_types('CA')
     # At least one model in each ring should be close to the reference
     for side in ('top', 'bottom'):
         rms = []
         for i in range(6):
             fname = 'output/model.%s.%d.pdb' % (side, i)
             m =  modeller.model(e, file=fname)
             a = modeller.alignment(e)
             a.append_model(ref, align_codes='ref')
             a.append_model(m, align_codes='model')
             rms.append(sel.superpose(m, a).rms)
             os.unlink(fname)
         self.assertTrue(min(rms) < 10.0)
     os.unlink('output/intermediate_asmb_sols.out')
     for side in ('top', 'bottom'):
         os.unlink('output/multifit.%s.output' % side)
         os.unlink('output/multifit.%s.output.symm.ref' % side)
         os.unlink('output/multifit.%s.param' % side)

    def test_glyc(self):
        """Test glycosylation benchmark"""
        os.chdir(os.path.join(TOPDIR, 'benchmark', 'input_glyc'))
        # Cleanup anything left over from a previous run
        shutil.rmtree('pred_dECALCrAS1000', ignore_errors=True)

        subprocess.check_call(['allosmod', 'setup'])
        # Setup should generate ligand and script:
        for f in ['lig.pdb', 'qsub.sh']:
            self.assertTrue(os.path.exists(f))
        # Run the protocol
        subprocess.check_call(['/bin/sh', '--login', './qsub.sh'])
        # Should generate more files:
        os.chdir('pred_dECALCrAS1000/2AAS.pdb_0')
        for f in ['align2.ali', 'allosmod.py', 'converted.rsr',
                  'model_glyc.log', 'model_glyc.py', 'pm.pdb.B99990001.pdb',
                  'pm.pdb.D00000001', 'pm.pdb.V99990001', 'run.log']:
            self.assertTrue(os.path.exists(f))

        # Generated model should have sugars added to second chain
        e = modeller.environ()
        e.io.hetatm = True
        m = modeller.model(e, file='pm.pdb.B99990001.pdb')
        self.assertEqual(len(m.chains), 2)
        self.assertEqual(len(m.chains[0].residues), 124)
        self.assertEqual(len(m.chains[1].residues), 16)
        self.assertEqual([r.name for r in m.chains[1].residues],
                         ['NAG', 'NAG', 'BMA', 'MAN', 'MAN', 'MAN', 'MAN',
                          'MAN', 'NAG', 'NAG', 'BMA', 'MAN', 'MAN', 'MAN',
                          'MAN', 'MAN'])

 def test_script6(self):
     """Test step 6 (model building and assessment)"""
     # Get inputs (outputs from steps 3 and 5)
     shutil.copy('precalculate_results/stage3_density_segmentation/' \
                 'groel_subunit_11.mrc', 'output')
     shutil.copy('precalculate_results/stage5_template_alignment/' \
                 'groel-1iokA.ali', 'output')
     # Make sure the script runs without errors
     p = subprocess.check_call(['scripts/' \
                                'script6_model_building_and_assessment.py'])
     # Check output models
     e = modeller.environ()
     for i in range(1, 11):
         base = 'output/P0A6F5.B9999%04d' % i
         pdb = base + '.pdb'
         trunc_pdb = base + '.truncated.pdb'
         trunc_fit_pdb = base + '.truncated.fitted.pdb'
         m = modeller.model(e, file=pdb)
         self.assertEqual(len(m.residues), 548)
         m = modeller.model(e, file=trunc_pdb)
         self.assertEqual(len(m.residues), 524)
         m = modeller.model(e, file=trunc_fit_pdb)
         self.assertEqual(len(m.residues), 524)
         os.unlink(pdb)
         os.unlink(trunc_pdb)
         os.unlink(trunc_fit_pdb)
     scores = 'output/model_building.scores.output'
     wc = len(open(scores).readlines())
     # Should be one line for each of 10 models, plus a header
     self.assertEqual(wc, 11)
     os.unlink(scores)

def main():
    import modeller
    ff1, ff2, aln_file, opts = parse_args()
    env = modeller.environ()
    # Read in HETATM records from template PDBs
    env.io.hetatm = True
    aln = salign0(env, ff1, ff2)
    aln.write(file=aln_file)

def make_model(msmseed):
    """
    Use MODELLER from the Sali lab to create a model between the target and template specified in the input

    Parameters
    ----------
    msmseed : MSMSeed
        object containing the alignment between target and template and template structure

    Returns
    -------
    msmseed : MSMSeed
        object containing the homology model built from the input alignment and template structure
    """

    import tempfile
    import os
    import modeller
    import modeller.automodel
    import shutil
    import simtk.openmm.app as app
    #if the target and template are the same, modeller dies.
    if msmseed.template_id == msmseed.target_id:
        msmseed.target_model = msmseed.template_structure
        return msmseed
    #first, we need to make a temp directory where we can put the files MODELLER needs
    temp_dir = tempfile.mkdtemp()
    try:
        os.chdir(temp_dir)
        alignment_file = open('aln_tmp.pir','w')
        alignment_file.writelines(msmseed.alignment)
        alignment_file.close()
        template_file = open(msmseed.template_id + '.pdb','w')
        template_pdb = msmseed.template_structure
        template_pdb.writeFile(template_pdb.topology, template_pdb.positions, template_file)
        template_file.close()
        modeller.log.none()
        env = modeller.environ()
        env.io.atom_files_directory = temp_dir
        a = modeller.automodel.allhmodel(env,
                                         # file with template codes and target sequence
                                         alnfile  = 'aln_tmp.pir',
                                         # PDB codes of the template
                                         knowns   = msmseed.template_id,
                                         # code of the target
                                         sequence = msmseed.target_id)
        a.make()
        tmp_model_pdbfilename = a.outputs[0]['name']
        target_model = modeller.model(env, file=tmp_model_pdbfilename)
        msmseed.sequence_similarity = target_model.seq_id
        msmseed.target_model = app.PDBFile(tmp_model_pdbfilename)
        msmseed.target_restraints = open('%s.rsr' % msmseed.target_id, 'r').readlines()
    except:
        msmseed.error_state = -1

    finally:
        shutil.rmtree(temp_dir)
    return msmseed

 def test_script5(self):
     """Test step 5 (template alignment)"""
     # Make sure the script runs without errors
     p = subprocess.check_call(['scripts/script5_template_alignment.py'])
     # Check output alignment
     e = modeller.environ()
     a = modeller.alignment(e, file='output/groel-1iokA.ali')
     self.assertEqual([x.code for x in a], ['1iok', 'P0A6F5'])
     os.unlink('output/groel-1iokA.ali')

def perform_sequence_alignment():
    e = modeller.environ()
    m1 = modeller.model(e, file='experimental.pdb')
    m2 = modeller.model(e, file='rosetta.pdb')
    aln = modeller.alignment(e)
    aln.append_model(m1, align_codes='experimental', atom_files='experimental.pdb')
    aln.append_model(m2, align_codes='rosetta')
    aln.align2d()
    aln.write(file='align.ali', alignment_format='PIR')

    def test_integrative_modeling(self):
        """Test the entire integrative modeling run"""
        import modeller
        # Compile the clustering program
        subprocess.check_call(['gfortran', 'cluster.f', 'u3best.f',
                               '-o', 'cluster.x'],
                              cwd='integrative_modeling/bin')

        # Run sampling
        subprocess.check_call(['./run_modeling.py'],
                              cwd='integrative_modeling')

        # Analysis
        subprocess.check_call(['bin/get_frames.sh'],
                              cwd='integrative_modeling')

        # Make sure that at least two of the three "known good" clusters
        # are reproduced
        clusters = glob.glob('integrative_modeling/clustering/clus.*.pdb')
        clusters = [x for x in clusters if '-' not in x]
        exp_clusters = glob.glob('model_refinement/cluster*/model.pdb')

        env = modeller.environ()
        n_cluster = 0
        rms = []
        cluster_match = [0] * len(clusters)
        exp_cluster_match = [0] * len(exp_clusters)
        # Get a matrix of RMSD between all clusters and the expected clusters
        for ncluster, cluster in enumerate(clusters):
            per_cluster = []
            for nexp_cluster, exp_cluster in enumerate(exp_clusters):
                mc = modeller.model(env, file=cluster)
                s = modeller.selection(mc)
                a = modeller.alignment(env)
                me = modeller.model(env, file=exp_cluster)
                a.append_model(mc, align_codes='clus')
                a.append_model(me, align_codes='exp_clus')
                # We only care about the global (non-cutoff) RMSD, so use a
                # large cutoff so that refine_local doesn't increase the number
                # of equivalent positions at the expense of worsening the RMSD
                r = s.superpose(me, a, rms_cutoff=999.)
                if r.rms < 15.0:
                    cluster_match[ncluster] += 1
                    exp_cluster_match[nexp_cluster] += 1
                per_cluster.append(r.rms)
            rms.append(per_cluster)
        # Count the number of clusters which are close to an expected cluster
        ncluster_match = len(cluster_match) - cluster_match.count(0)
        # Count the number of expected clusters which are close to a cluster
        nexp_cluster_match = len(exp_cluster_match) - exp_cluster_match.count(0)
        # Make sure that at least 2 of the 3 expected clusters is close to one
        # of the clusters we produced (but not all the *same* cluster)
        self.assertTrue(ncluster_match >= 2 and nexp_cluster_match >= 2,
                        "Could not find any match between the %d clusters "
                        "found in this test and 2 of the 3 'known good' "
                        "clusters (match defined as all-atom RMSD less than "
                        "15.0A). RMSD matrix: %s" % (len(clusters), str(rms)))

 def make_model(self):
     import modeller
     env = modeller.environ()
     env.edat.dynamic_sphere= False
     with open('test.pdb', 'w') as fh:
         fh.write("ATOM      2  CA  ALA     1      27.449  14.935   5.140  1.00 29.87           C\n")
     m = modeller.model(env, file='test.pdb')
     os.unlink('test.pdb')
     return m

 def get_environ(self):
     """Get a Modeller environ object"""
     if not hasattr(self, '_modeller_environ'):
         # Speed tests up a little by only creating this object once
         env = modeller.environ()
         env.libs.topology.read('${LIB}/top_heav.lib')
         env.libs.parameters.read('${LIB}/par.lib')
         Tests._modeller_environ = env
     return self._modeller_environ

def setupENV():
    #setup Modeller
    env = modeller.environ()
    MPATH=epmv.__path__[0]+'/extension/Modeller/'
    env.io.atom_files_directory = [MPATH]
    env.edat.dynamic_sphere = True
    env.libs.topology.read(file='$(LIB)/top_heav.lib')
    env.libs.parameters.read(file='$(LIB)/par.lib')
    return env

def spline(pdb_file, in_restraints, out_restraints):
    import modeller
    # Needed to keep our custom form alive for restraints.read()
    from allosmod.modeller.forms import TruncatedGaussian

    e = modeller.environ()
    m = modeller.model(e, file=pdb_file)
    m.restraints.read(file=in_restraints)
    convert_restraints(m.restraints)
    m.restraints.write(file=out_restraints)

def mk_strct_al_modeller(strct_data1, strct_data2):
    _stdout = sys.stdout
    sys.stdout = open(os.devnull, 'w')

    tmp_file = tempfile.NamedTemporaryFile(suffix=".fasta", delete=False)
    env = m.environ()

    aln = m.alignment(env)
    code1 = 'pdb' + strct_data1['id']
    code2 = 'pdb' + strct_data2['id']
    chain1 = strct_data1['chain_id']
    chain2 = strct_data2['chain_id']
    env.io.atom_files_directory = ['.', PDB_DIR]
    result = {}
    try:
        for (code, chain) in ((code1, chain1), (code2, chain2)):
            mdl = m.model(env, file=code, model_segment=('FIRST:'+chain,
                                                         'LAST:'+chain))
            aln.append_model(mdl, atom_files=code, align_codes=code+chain)

        for (weights, write_fit, whole) in (((1., 0., 0., 0., 1., 0.), False,
                                             True),
                                            ((1., 0.5, 1., 1., 1., 0.), False,
                                             True),
                                            ((1., 1., 1., 1., 1., 0.), True,
                                             False)):
            r = aln.salign(rms_cutoff=3.5, normalize_pp_scores=False,
                           rr_file='$(LIB)/as1.sim.mat', overhang=30,
                           gap_penalties_1d=(-450, -50),
                           gap_penalties_3d=(0, 3), gap_gap_score=0,
                           gap_residue_score=0,
                           alignment_type='tree', # If 'progresive', the tree is not
                                                  # computed and all structures will be
                                                  # aligned sequentially to the first
                           #ext_tree_file='1is3A_exmat.mtx', # Tree building can be avoided
                                                             # if the tree is input
                           feature_weights=weights, # For a multiple sequence alignment only
                                                    # the first feature needs to be non-zero
                           improve_alignment=True, fit=True, write_fit=False,
                           write_whole_pdb=whole, output='ALIGNMENT QUALITY')
        if r.qscorepct > 70:
            aln.write(file=tmp_file.name, alignment_format='FASTA')
            with open(tmp_file.name) as a:
                alignment = unwrap(a.read().splitlines())

            for i in range(len(alignment[1])):
                if alignment[1] != '-' and alignment[3] != '-':
                    pos1 = get_real_position_al(alignment[1], i)
                    pos2 = get_real_position_al(alignment[3], i)
                    result[pos1] = pos2
    except:
        print 'Modeller failed'
    sys.stdout.close()
    sys.stdout = _stdout
    return result

 def test_script1(self):
     """Test step 1 (build profile)"""
     # Make sure the script runs without errors
     p = subprocess.check_call(['scripts/script1_build_profile.py'])
     # Make sure the profile contains the sequences we expect
     e = modeller.environ()
     a = modeller.alignment(e, file='output/build_profile.ali')
     self.assertEqual(sorted(s.code for s in a),
                      sorted(self.templates) + ['P0A6F5'])
     os.unlink('output/build_profile.prf')
     os.unlink('output/build_profile.ali')

def main():
    rsr_file, pdb_file, sclbreak, opts = parse_args()
    e = modeller.environ()
    e.io.hetatm = True
    a = ChargedContactFinder(e, rsr_file, pdb_file)
    a.find()
    a.print_contacts(open('contpres.dat', 'w'))
    if opts.cdensity_cutoff is not None:
        a.print_cdensity(opts.cdensity_cutoff, sclbreak, open('break.dat', 'a'))
    else:
        a.print_all_buried(sclbreak, open('break.dat', 'a'))

def main():
    import modeller
    file1, file2, rcut = parse_args()
    e = modeller.environ()
    e.io.hetatm = True
    mdl1 = modeller.model(e, file=file1)
    mdl2 = modeller.model(e, file=file2)
    for ri, rj, dist in get_inter_contacts(e, mdl1, mdl2, rcut):
        print("  %6s  %2s  %6s  %2s  %3s  %3s%3d%11.3f  %1d  %1d"
              % (ri.num, ri.chain.name, rj.num, rj.chain.name,
                 ri.pdb_name, rj.pdb_name, get_contact_type(ri, rj),
                 dist, 6, 6))

def count_alignments(aln_file, target, templates):
    import modeller
    modeller.log.none()
    env = modeller.environ()
    aln = modeller.alignment(env, file=aln_file)
    target = aln[target]
    templates = [aln[t] for t in templates]
    num_align = 0
    for r in target.residues:
        for template in templates:
            if r.get_aligned_residue(template) is not None:
                num_align += 1
    return num_align, len(target.residues)

def main():
    pdb1, ligand, pdb2, rcut, opts = parse_args()
    e = modeller.environ()
    e.io.hetatm = True
    a = AllostericSiteFinder(e, pdb1, ligand, pdb2, rcut)
    try:
        if opts.output_pdb:
            a.find().write(file=opts.output_pdb)
        if opts.atom_list:
            a.write_atom_list(open(opts.atom_list, 'w'))
    except AllostericSiteError as err:
        print(str(err), file=sys.stderr)
        sys.exit(1)