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 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 find(self):
        """Return a Modeller selection corresponding to the allosteric site.
           @raise AllostericSiteError on error."""
        if self.__allosteric_site is None:
            # align PDB2 to PDB1 and superimpose antigen
            try:
                salign0(self.env, self.pdb1, self.pdb2)
            except modeller.ModellerError as err:
                raise AllostericSiteError("Could not align %s with %s: %s. "
                               "This is usually due to a poor alignment."
                               % (self.pdb2, self.pdb1, str(err)))
            pmfit = get_fit_filename(self.pdb2)

            # determine residues in PDB2 that contact LIG1
            self.__pmfit = modeller.model(self.env, file=pmfit)
            lig1 = modeller.model(self.env, file=self.ligand)
            self.__allosteric_site = \
                 modeller.selection([ri for ri, rj, dist \
                            in get_inter_contacts(self.env, self.__pmfit, lig1,
                                                  self.rcut)])
            os.unlink(pmfit)
            os.unlink(get_fit_filename(self.pdb1))
            if len(self.__allosteric_site) == 0:
                raise AllostericSiteError("No allosteric site found")
        return self.__allosteric_site

 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 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 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 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_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_feature_sidechain_biso(self):
     """Check average sidechain Biso feature"""
     env = self.get_environ()
     mlib = self.get_mdt_library()
     self.assertRaises(ValueError, mdt.features.SidechainBiso, mlib, bins=mdt.uniform_bins(5, 0, 10), protein=3)
     sidechain_biso = mdt.features.SidechainBiso(mlib, bins=mdt.uniform_bins(5, 0, 10))
     mdl = modeller.model(env)
     mdl.build_sequence("A")
     aln = modeller.alignment(env)
     aln.append_model(mdl, align_codes="test")
     s = aln[0]
     # Mainchain atom Biso should be ignored:
     for mainchain in ("N:1", "C:1", "O:1", "OXT:1", "CA:1"):
         s.atoms[mainchain].biso = 1000
     for (biso, bin) in (
         (22, 2),
         (32, 3),  # Map regular values to bins
         (0, -1),  # Zero Biso should be "undefined"
         (1, 3),
     ):  # Biso < 2 is multiplied by 4pi^2
         s.atoms["CB:1"].biso = biso
         m = mdt.Table(mlib, features=sidechain_biso)
         m.add_alignment(aln)
         self.assertEqual(m.shape, (6,))
         self.assertEqual(m.sum(), 1)
         self.assertEqual(m[bin], 1)

 def test_disulfide(self):
     """Test handling of disulfide bonds"""
     mlib = self.get_all_libraries()
     bsep = mdt.features.AtomBondSeparation(mlib,
                                     bins=mdt.uniform_bins(20, 0, 1.0))
     bsep_ss = mdt.features.AtomBondSeparation(mlib,
                                     bins=mdt.uniform_bins(20, 0, 1.0),
                                     disulfide=True)
     env = self.get_environ()
     mdl = modeller.model(env)
     mdl.build_sequence('CC')
     # When SG-SG distance is small enough, an extra bond
     # (separation feature = 1) should be detected, but only with
     # disulfide=True
     for (dist, num) in [(2.6, 11.0), (2.4, 12.0)]:
         sg1 = mdl.residues[0].atoms['SG']
         sg2 = mdl.residues[1].atoms['SG']
         sg1.x = sg1.y = sg1.z = 0.
         sg2.x = sg2.y = 0.
         sg2.z = dist
         a = modeller.alignment(env)
         a.append_model(mdl, atom_files='test', align_codes='test')
         m = mdt.Table(mlib, features=bsep)
         m.add_alignment(a, residue_span_range=(-999,0,0,999))
         self.assertEqual(m[1], 11.0)
         m2 = mdt.Table(mlib, features=bsep_ss)
         m2.add_alignment(a, residue_span_range=(-999,0,0,999))
         self.assertEqual(m2[1], num)

 def setup_atoms(self, env):
     self.m = modeller.model(env, file=self.pdb_file)
     self.atoms = [Atom(a) for a in self.m.atoms]
     self.contacts = get_contacts(self.contacts_pdbs, self.rcut)
     if self.break_file:
         self.breaks = get_breaks(open(self.break_file))
     else:
         self.breaks = {}
     self.beta_structure = get_beta(self.pdb_file)
     NUCLEIC_ACIDS = dict.fromkeys(['ADE', 'A', 'DA', 'THY', 'T', 'DT',
                                    'URA', 'U', 'DU', 'GUA', 'G', 'DG',
                                    'CYT', 'C', 'DC'])
     BACKBONE_ATOMS = dict.fromkeys(['CA', 'CB', 'O', 'N', 'C', 'OT', 'NA',
                                     'NB', 'NC', 'ND', 'C1A', 'C2A', 'C3A',
                                     'C4A', 'C1B', 'C2B', 'C3B', 'C4B',
                                     'C1C', 'C2C', 'C3C', 'C4C', 'C1D',
                                     'C2D', 'C3D', 'C4D'])
     for a in self.atoms:
         r = a.a.residue
         if r.pdb_name in NUCLEIC_ACIDS:
             a.isNUC = True
             a.torestr = get_nuc_restrained(a.a.name, r.pdb_name)
             for rj in range(1, len(self.m.residues) + 1):
                 self.contacts[(r.index,rj)] = True
         if a.a.name in BACKBONE_ATOMS or r.pdb_name in NUCLEIC_ACIDS:
             a.isSC = False
             a.isCA = a.a.name == 'CA'
             a.isCB = a.a.name == 'CB'
         else:
             a.isSC = a.a.name != 'H'
     for a, asrs in zip(self.atoms,
                        parse_atomlist_asrs(open(self.atomlist_asrs))):
         a.isAS = asrs

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 copy_to_modeller(env, particles):
    fh = open("temp_particles.pdb", "w")
    for i in range(len(particles)):
        fh.write("ATOM  %5d  N   ALA     0       0.000   0.000   0.000  "
                 "1.00  0.00           C  \n" % (i))
    fh.close()
    mdl = modeller.model(env, file='temp_particles.pdb')
    os.unlink('temp_particles.pdb')
    return mdl

 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 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 build_mdt_from_sequence(self, mlib, features, seq, **keys):
        """Build a simple test MDT for a given sequence"""
        env = self.get_environ()
        mdl = modeller.model(env)
        mdl.build_sequence(seq)

        m = mdt.Table(mlib, features=features)
        a = modeller.alignment(env)
        a.append_model(mdl, atom_files='test', align_codes='test')
        m.add_alignment(a, **keys)
        return m

    def get_test_mdt(self, mlib, features):
        env = self.get_environ()
        mdl = modeller.model(env)
        mdl.build_sequence('C')

        m = mdt.Table(mlib, features=features)
        a = modeller.alignment(env)
        a.append_model(mdl, atom_files='test', align_codes='test')
        m.add_alignment(a)
        m = m.reshape(features, [0] * len(features), [-1] * len(features))
        return m

def salign0(env, ff1, ff2):
    import modeller
    aln = modeller.alignment(env)
    code = ff1
    mdl = modeller.model(env, file=code, model_segment=('FIRST:@', 'END:'))
    fit_atoms = determine_fit_atoms(mdl)
    aln.append_model(mdl, atom_files=code, align_codes=code)
    code = ff2
    mdl = modeller.model(env, file=code, model_segment=('FIRST:@', 'END:'))
    aln.append_model(mdl, atom_files=code, align_codes=code)

    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)):
        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, fit_atoms=fit_atoms,
                   alignment_type='tree', feature_weights=weights,
                   improve_alignment=True, fit=True, write_fit=write_fit,
                   write_whole_pdb=whole, output='ALIGNMENT QUALITY')
    return aln