Python Molecule.read方法代码示例

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
def _processSim(sim, projectionlist, uqmol, skip):
    pieces = sim.trajectory
    try:
        if uqmol is not None:
            mol = uqmol.copy()
        else:
            mol = Molecule(sim.molfile)
        logger.debug(pieces[0])

        mol.read(pieces, skip=skip)
        #Gianni testing
        #_highfreqFilter(mol,10)
 
        data = []
        for p in projectionlist:
            pj = _project(p, mol)
            if pj.ndim == 1:
                pj = np.atleast_2d(pj).T
            data.append(pj)
        data = np.hstack(data)
        if data.dtype == np.float64:
            data = data.astype(np.float32)
    except Exception as e:
        logger.warning('Error in simulation with id: ' + str(sim.simid) + '. "' + e.__str__() + '"')
        return None, None, None, True

    return data, _calcRef(pieces, mol.fileloc), mol.fstep, False

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
def _filtSim(i, sims, outFolder, filterSel):
    name = _simName(sims[i].trajectory[0])
    directory = path.join(outFolder, name)
    if not path.exists(directory):
        makedirs(directory)

    logger.debug('Processing trajectory ' + name)

    fmolfile = path.join(outFolder, 'filtered.pdb')
    (traj, outtraj) = _renameSims(sims[i].trajectory, name, outFolder)
    if not traj:
        ftrajectory = _listXTCs(path.join(outFolder, name))
        return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile)

    try:
        mol = Molecule(sims[i].molfile)
    except:
        logger.warning('Error! Skipping simulation ' + name)
        return

    sel = mol.atomselect(filterSel)

    for j in range(0, len(traj)):
        try:
            mol.read(traj[j])
        except IOError as e:
            logger.warning(e.strerror + ', skipping trajectory')
            break

        mol.write(outtraj[j], sel)

    ftrajectory = _listXTCs(path.join(outFolder, name))
    #bar.progress()
    return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
def _writeInputsFunction(i, f, epoch, inputpath, coorname):
    regex = re.compile('(e\d+s\d+)_')
    frameNum = f.frame
    piece = f.piece
    if f.sim.parent is None:
        currSim = f.sim
    else:
        currSim = f.sim.parent

    traj = currSim.trajectory[piece]
    if currSim.input is None:
        raise NameError('Could not find input folder in simulation lists. Cannot create new simulations.')

    wuName = _simName(traj)
    res = regex.search(wuName)
    if res:  # If we are running on top of adaptive, use the first name part for the next sim name
        wuName = res.group(1)

    # create new job directory
    newName = 'e' + str(epoch) + 's' + str(i + 1) + '_' + wuName + 'p' + str(piece) + 'f' + str(frameNum)
    newDir = path.join(inputpath, newName, '')

    # copy previous input directory including input files
    copytree(currSim.input, newDir, symlinks=False, ignore=ignore_patterns('*.coor', '*.rst', '*.out', *_IGNORE_EXTENSIONS))

    # overwrite input file with new one. frameNum + 1 as catdcd does 1 based indexing

    mol = Molecule(currSim.molfile)  # Always read the mol file, otherwise it does not work if we need to save a PDB as coorname
    mol.read(traj)
    mol.dropFrames(keep=frameNum)  # Making sure only specific frame to write is kept
    mol.write(path.join(newDir, coorname))

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
def _loadMols(self, i, rel, molfile, wrapsel, alignsel, refmol):
    frames = self.data.rel2sim(rel)
    mol = Molecule(molfile)
    trajs = np.empty(0, dtype=str)
    frs = np.empty(0, dtype=int)
    for f in frames:
        trajs = np.append(trajs, f.sim.trajectory[f.piece])
        frs = np.append(frs, f.frame)
    mol.read(trajs, frames=frs)
    if len(wrapsel) > 0:
        mol.wrap(wrapsel)
    if refmol is not None:
        mol.align(alignsel, refmol=refmol)
    return mol

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
    def test_tmscore(self):
        from htmd.molecule.molecule import Molecule
        expectedTMscore = np.array([0.21418524, 0.2367377, 0.23433833, 0.21362964, 0.20935164,
                                    0.20279461, 0.27012895, 0.22675238, 0.21230793, 0.2372011])
        expectedRMSD = np.array([3.70322128, 3.43637027, 3.188193, 3.84455877, 3.53053882,
                                 3.46781854, 2.93777629, 2.97978692, 2.70792428, 2.63051318])

        mol = Molecule(os.path.join(home(dataDir='tmscore'), 'filtered.pdb'))
        mol.read(os.path.join(home(dataDir='tmscore'), 'traj.xtc'))
        ref = Molecule(os.path.join(home(dataDir='tmscore'), 'ntl9_2hbb.pdb'))
        tmscore, rmsd = molTMscore(mol, ref, mol.atomselect('protein'), ref.atomselect('protein'))

        self.assertTrue(np.allclose(tmscore, expectedTMscore))
        self.assertTrue(np.allclose(rmsd, expectedRMSD))

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
    def _writeInputs(self, simsframes, epoch=None):
        if epoch is None:
            epoch = self._getEpoch() + 1

        test = glob(path.join(self.inputpath, 'e' + str(epoch) + '*'))
        if len(test) != 0:
            raise NameError('Input dirs of epoch ' + str(epoch) + ' already exists.')

        if path.exists(path.join(self.inputpath, 'e' + str(epoch) + '_writeinputs.log')):
            raise NameError('Epoch logfile already exists. Cant overwrite it.')

        fid = open(path.join(self.inputpath, 'e' + str(epoch) + '_writeinputs.log'), 'w')

        regex = re.compile('(e\d+s\d+)_')
        for i, f in enumerate(simsframes):
            frameNum = f.frame
            piece = f.piece
            #print(frameNum)
            if f.sim.parent is None:
                currSim = f.sim
            else:
                currSim = f.sim.parent

            traj = currSim.trajectory[piece]
            if currSim.input is None:
                raise NameError('Could not find input folder in simulation lists. Cannot create new simulations.')

            wuName = _simName(traj)
            res = regex.search(wuName)
            if res:  # If we are running on top of adaptive, use the first name part for the next sim name
                wuName = res.group(1)

            # create new job directory
            newName = 'e' + str(epoch) + 's' + str(i+1) + '_' + wuName + 'p' + str(piece) + 'f' + str(frameNum)
            newDir = path.join(self.inputpath, newName, '')
            # copy previous input directory including input files
            copytree(currSim.input, newDir, symlinks=False, ignore=ignore_patterns('*.dcd', '*.xtc', '*.coor'))
            # overwrite input file with new one. frameNum + 1 as catdcd does 1 based indexing
            mol = Molecule()
            mol.read(traj)
            mol.frame = frameNum
            mol.write(path.join(newDir, 'input.coor'))

            # write nextInput file
            fid.write('# {0} \n{1} {2}\n'.format(newName, traj, frameNum))

        fid.close()

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
    def test_project(self):
        from htmd.molecule.molecule import Molecule
        from htmd.home import home
        from os import path
        mol = Molecule(path.join(home(), 'data', 'metricdistance', 'filtered.pdb'))
        mol.read(path.join(home(), 'data', 'metricdistance', 'traj.xtc'))
        ref = mol.copy()
        ref.coords = np.atleast_3d(ref.coords[:, :, 0])
        metr = MetricCoordinate('protein and name CA')
        data = metr.project(mol)

        lastcoors = np.array([-24.77000237, -27.76000023, -30.44000244, -33.65000153,
                              -33.40999985, -36.32000351, -36.02000427, -36.38000107,
                              -39.61000061, -41.01000214, -43.80000305, -45.56000137,
                              -45.36000061, -47.13000488, -49.54000473, -50.6000061 ,
                              -50.11000061, -52.15999985, -55.1400032 , -55.73000336], dtype=np.float32)
        assert np.all(np.abs(data[-1, -20:] - lastcoors) < 0.001), 'Coordinates calculation is broken'

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
    def test_project_align(self):
        from htmd.molecule.molecule import Molecule
        from htmd.home import home
        from os import path
        mol = Molecule(path.join(home(), 'data', 'metricdistance', 'filtered.pdb'))
        mol.read(path.join(home(), 'data', 'metricdistance', 'traj.xtc'))
        ref = mol.copy()
        ref.coords = np.atleast_3d(ref.coords[:, :, 0])
        metr = MetricCoordinate('protein and name CA', ref)
        data = metr.project(mol)

        lastcoors = np.array([6.79283285, 5.55226946, 4.49387407, 2.94484425,
                              5.36937141, 3.18590879, 5.75874281, 5.48864174,
                              1.69625032, 1.58790839, 0.57877392, -2.66498065,
                              -3.70919156, -3.33702421, -5.38465405, -8.43286991,
                              -8.15859032, -7.85062265, -10.92551327, -13.70733166], dtype=np.float32)
        assert np.all(np.abs(data[-1, -20:] - lastcoors) < 0.001), 'Coordinates calculation is broken'

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
    def _analyse(self, mol, pdb, rtf, prm, traj, ftraj):
        t = Molecule(pdb)
        t.read(traj)
        t.filter('not water')
        t.write(ftraj)
        m = FFMolecule(filename=mol, rtf=rtf, prm=prm)
        m.read(ftraj)
        torsions = m.getRotatableDihedrals()
        # For each torsion
        for i in range(len(torsions)):
            # Create title
            title = '{}-{}-{}-{}'.format(m.name[torsions[i][0]], m.name[torsions[i][1]], m.name[torsions[i][2]],
                                         m.name[torsions[i][3]])
            # Measure
            (r, theta) = self._measure_torsion(torsions[i], m.coords)

            self._plot_scatter(r, theta, title)
            self._plot_hist(theta, title)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
    def _fb_potential2restraints(self, inputdir):
        from htmd.molecule.molecule import Molecule
        restraints = list()

        fb_box = np.array(self.fb_box)
        # convert fb_box to width
        width = list(np.concatenate(np.diff(np.array([fb_box[::2], fb_box[1::2]]), axis=0)))

        # If fb_box is not symmetrical
        if not np.all(fb_box[::2] == -fb_box[1::2]):
            # convert fb_box and fb_reference to fbcentre and width
            mol = Molecule(os.path.join(inputdir, self.acemd.structure))
            mol.read(os.path.join(inputdir, self.acemd.coordinates))
            fb_refcentre = mol.get('coords', sel=self.fb_reference).mean(axis=0).squeeze()

            fbcentre = list(np.around(np.mean(np.array([fb_box[::2], fb_box[1::2]]), axis=0) + fb_refcentre, 3))
            restraints.append(GroupRestraint(self.fb_selection, width, [(self.fb_k, 0)], fbcentre=fbcentre))
        else:
            restraints.append(GroupRestraint(self.fb_selection, width, [(self.fb_k, 0)], fbcentresel=self.fb_reference))

        return restraints

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
def _filtSim(i, sims, outFolder, filterSel):
    name = _simName(sims[i].trajectory[0])
    directory = path.join(outFolder, name)
    if not path.exists(directory):
        makedirs(directory)

    logger.debug('Processing trajectory ' + name)

    fmolfile = path.join(outFolder, 'filtered.pdb')
    (traj, outtraj) = _renameSims(sims[i].trajectory, name, outFolder)
    if not traj:
        ftrajectory = _autoDetectTrajectories(path.join(outFolder, name))
        numframes = _getNumFrames(sims[i], ftrajectory)
        return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile, numframes=numframes)

    try:
        from htmd.molecule.molecule import Molecule
        mol = Molecule(sims[i].molfile)
    except:
        logger.warning('Error! Skipping simulation ' + name)
        return

    sel = mol.atomselect(filterSel)

    for j in range(0, len(traj)):
        try:
            mol.read(traj[j])
        except IOError as e:
            logger.warning('{}, skipping trajectory'.format(e))
            break

        mol.write(outtraj[j], sel)

    ftrajectory = _autoDetectTrajectories(path.join(outFolder, name))
    numframes = _getNumFrames(sims[i], ftrajectory)
    return Sim(simid=sims[i].simid, parent=sims[i], input=None, trajectory=ftrajectory, molfile=fmolfile, numframes=numframes)

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
    plt.imshow(contactmap, interpolation="nearest", aspect="equal")
    plt.colorbar()
    # plt.axis('off')
    # plt.tick_params(axis='x', which='both', bottom='off', top='off', labelbottom='off')
    # plt.tick_params(axis='y', which='both', left='off', right='off', labelleft='off')
    plt.show()


if __name__ == "__main__":
    from htmd.molecule.molecule import Molecule
    from htmd.home import home
    from os import path

    mol = Molecule(path.join(home(), "data", "metricdistance", "filtered.pdb"))
    mol.read(path.join(home(), "data", "metricdistance", "traj.xtc"))
    metr = MetricDistance("protein and name CA", "resname MOL and noh", metric="contacts")
    data = metr.project(mol)
    contframes = np.array(
        [
            46,
            47,
            48,
            49,
            50,
            51,
            52,
            53,
            54,
            55,
            75,

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]

#.........这里部分代码省略.........
    #_checkProteinGaps(mol)
    if patches is None:
        patches = []
    if isinstance(patches, str):
        patches = [patches]
    # Find protonated residues and add patches for them
    patches += _protonationPatches(mol)

    f = open(path.join(outdir, 'build.vmd'), 'w')
    f.write('# psfgen file generated by charmm.build\n')
    f.write('package require psfgen;\n')
    f.write('psfcontext reset;\n\n')

    # Copying and printing out the topologies
    for i in range(len(topo)):
        if topo[i][0] != '.' and path.isfile(path.join(charmmdir, topo[i])):
            topo[i] = path.join(charmmdir, topo[i])
        localname = '{}.'.format(i) + path.basename(topo[i])
        shutil.copy(topo[i], path.join(outdir, localname))
        f.write('topology ' + localname + '\n')
    f.write('\n')

    _printAliases(f)

    # Printing out segments
    logger.info('Writing out segments.')
    segments = _getSegments(mol)
    for seg in segments:
        pdbname = 'segment' + seg + '.pdb'
        mol.write(path.join(outdir, pdbname), sel='segid '+seg)

        segatoms = mol.atomselect('segid {}'.format(seg))
        segwater = mol.atomselect('segid {} and water'.format(seg))

        f.write('segment ' + seg + ' {\n')
        if np.all(segatoms == segwater):  # If segment only contains waters, set: auto none
            f.write('\tauto none\n')
        f.write('\tpdb ' + pdbname + '\n')
        if caps is not None and seg in caps:
            for c in caps[seg]:
                f.write('\t' + c + '\n')
        f.write('}\n')
        f.write('coordpdb ' + pdbname + ' ' + seg + '\n\n')

    # Printing out patches for the disulfide bridges
    if disulfide is None:
        disulfide = detectDisulfideBonds(mol)

    if len(disulfide) != 0:
        for d in disulfide:
            f.write('patch DISU {}:{} {}:{}\n'.format(d.segid1, d.resid1, d.segid2, d.resid2))
        f.write('\n')

    # Printing out extra patches
    if len(patches) != 0:
        for p in patches:
            f.write(p + '\n')
        f.write('\n')

    f.write('guesscoord\n')
    f.write('writepsf ' + prefix + '.psf\n')
    f.write('writepdb ' + prefix + '.pdb\n')
    #f.write('quit\n')
    f.close()

    if param is not None:
        _charmmCombine(param, path.join(outdir, 'parameters'))

    molbuilt = None
    if execute:
        logpath = os.path.abspath('{}/log.txt'.format(outdir))
        logger.info('Starting the build.')
        currdir = os.getcwd()
        os.chdir(outdir)
        f = open(logpath, 'w')
        #call([vmd, '-dispdev', 'text', '-e', './build.vmd'], stdout=f)
        call([psfgen, './build.vmd'], stdout=f)
        f.close()
        _logParser(logpath)
        os.chdir(currdir)
        logger.info('Finished building.')

        if path.isfile(path.join(outdir, 'structure.pdb')) and path.isfile(path.join(outdir, 'structure.psf')):
            molbuilt = Molecule(path.join(outdir, 'structure.pdb'))
            molbuilt.read(path.join(outdir, 'structure.psf'))
        else:
            raise NameError('No structure pdb/psf file was generated. Check {} for errors in building.'.format(logpath))

        if ionize:
            shutil.move(path.join(outdir, 'structure.pdb'), path.join(outdir, 'structure.noions.pdb'))
            shutil.move(path.join(outdir, 'structure.psf'), path.join(outdir, 'structure.noions.psf'))
            totalcharge = np.sum(molbuilt.charge)
            nwater = np.sum(molbuilt.atomselect('water and noh'))
            anion, cation, anionatom, cationatom, nanion, ncation = ionizef(totalcharge, nwater, saltconc=saltconc, ff='charmm', anion=saltanion, cation=saltcation)
            newmol = ionizePlace(mol, anion, cation, anionatom, cationatom, nanion, ncation)
            # Redo the whole build but now with ions included
            return build(newmol, topo=topo, param=param, prefix=prefix, outdir=outdir, ionize=False, caps=caps,
                         execute=execute, saltconc=saltconc, disulfide=disulfide, patches=patches, psfgen=psfgen)
    _checkFailedAtoms(molbuilt)
    return molbuilt

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]

#.........这里部分代码省略.........
        param = []
    if caps is None:
        caps = _defaultCaps(mol)

    _missingSegID(mol)
    _checkMixedSegment(mol)

    logger.info('Converting CHARMM membranes to AMBER.')
    mol = _charmmLipid2Amber(mol)

    #_checkProteinGaps(mol)
    _applyCaps(mol, caps)

    f = open(path.join(outdir, 'tleap.in'), 'w')
    f.write('# tleap file generated by amber.build\n')

    # Printing out the forcefields
    for force in ff:
        f.write('source ' + force + '\n')
    f.write('\n')

    # Loading TIP3P water parameters
    f.write('# Loading ions and TIP3P water parameters\n')
    f.write('loadamberparams frcmod.ionsjc_tip3p\n\n')

    # Printing out topologies
    logger.info('Writing prepi files.')
    f.write('# Loading prepi topologies\n')
    for t in topo:
        shutil.copy(t, outdir)
        f.write('loadamberprep ' + path.basename(t) + '\n')
    f.write('\n')

    # Printing and loading the PDB file. AMBER can work with a single PDB file if the segments are separate by TER
    logger.info('Writing PDB file for input to tleap.')
    pdbname = path.join(outdir, 'input.pdb')
    mol.write(pdbname)
    if not os.path.isfile(pdbname):
        raise NameError('Could not write a PDB file out of the given Molecule.')
    f.write('# Loading the system\n')
    f.write('mol = loadpdb input.pdb\n\n')

    # Printing out patches for the disulfide bridges
    '''if disulfide is None and not ionize:
        logger.info('Detecting disulfide bonds.')
        disulfide = detectDisulfideBonds(mol)

    if not ionize and len(disulfide) != 0:  # Only make disu bonds after ionizing!
        f.write('# Adding disulfide bonds\n')
        for d in disulfide:
            # Convert to stupid amber residue numbering
            uqseqid = sequenceID(mol.resid)
            uqres1 = int(np.unique(uqseqid[mol.atomselect('segid {} and resid {}'.format(d.segid1, d.resid1))]))
            uqres2 = int(np.unique(uqseqid[mol.atomselect('segid {} and resid {}'.format(d.segid2, d.resid2))]))
            # Rename the CYS to CYX if there is a disulfide bond
            mol.set('resname', 'CYX', sel='segid {} and resid {}'.format(d.segid1, d.resid1))
            mol.set('resname', 'CYX', sel='segid {} and resid {}'.format(d.segid2, d.resid2))
            f.write('bond mol.{}.SG mol.{}.SG\n'.format(uqres1, uqres2))
        f.write('\n')'''

    f.write('# Writing out the results\n')
    f.write('savepdb mol ' + prefix + '.pdb\n')
    f.write('saveamberparm mol ' + prefix + '.prmtop ' + prefix + '.crd\n')
    f.write('quit')
    f.close()

    molbuilt = None
    if execute:
        logpath = os.path.abspath('{}/log.txt'.format(outdir))
        logger.info('Starting the build.')
        currdir = os.getcwd()
        os.chdir(outdir)
        f = open(logpath, 'w')
        try:
            call([tleap, '-f', './tleap.in'], stdout=f)
        except:
            raise NameError('tleap failed at execution')
        f.close()
        os.chdir(currdir)
        logger.info('Finished building.')

        if path.getsize(path.join(outdir, 'structure.pdb')) != 0 and path.getsize(path.join(outdir, 'structure.prmtop')) != 0:
            molbuilt = Molecule(path.join(outdir, 'structure.pdb'))
            molbuilt.read(path.join(outdir, 'structure.prmtop'))
            molbuilt.bonds = []  # Causes problems in ionization mol.remove and mol._removeBonds
        else:
            raise NameError('No structure pdb/prmtop file was generated. Check {} for errors in building.'.format(logpath))

        if ionize:
            shutil.move(path.join(outdir, 'structure.pdb'), path.join(outdir, 'structure.noions.pdb'))
            shutil.move(path.join(outdir, 'structure.crd'), path.join(outdir, 'structure.noions.crd'))
            shutil.move(path.join(outdir, 'structure.prmtop'), path.join(outdir, 'structure.noions.prmtop'))
            totalcharge = np.sum(molbuilt.charge)
            nwater = np.sum(molbuilt.atomselect('water and noh'))
            anion, cation, anionatom, cationatom, nanion, ncation = ionizef(totalcharge, nwater, saltconc=saltconc, ff='amber', anion=saltanion, cation=saltcation)
            newmol = ionizePlace(molbuilt, anion, cation, anionatom, cationatom, nanion, ncation)
            # Redo the whole build but now with ions included
            return build(newmol, ff=ff, topo=topo, param=param, prefix=prefix, outdir=outdir, caps={}, ionize=False,
                         execute=execute, saltconc=saltconc, disulfide=disulfide, tleap=tleap)
    return molbuilt

# 需要导入模块: from htmd.molecule.molecule import Molecule [as 别名]
# 或者: from htmd.molecule.molecule.Molecule import read [as 别名]
    return angles, dihedrals


# A test method
if __name__ == "__main__":
    from htmd.molecule.molecule import Molecule
    from htmd.home import home
    import numpy as np
    from os import path
    import doctest
    #doctest.testmod(extraglobs={"mol" : Molecule("3PTB")})

    expectedTMscore = np.array([ 0.21418524,  0.2367377 ,  0.23433833,  0.21362964,  0.20935164,
        0.20279461,  0.27012895,  0.22675238,  0.21230793,  0.2372011 ])
    expectedRMSD = np.array([ 3.70322128,  3.43637027,  3.188193  ,  3.84455877,  3.53053882,
        3.46781854,  2.93777629,  2.97978692,  2.70792428,  2.63051318])

    mol = Molecule(path.join(home(), 'data', 'tmscore', 'filtered.pdb'))
    mol.read(path.join(home(), 'data', 'tmscore', 'traj.xtc'))
    ref = Molecule(path.join(home(), 'data', 'tmscore', 'ntl9_2hbb.pdb'))
    tmscore, rmsd = molTMscore(mol, ref, mol.atomselect('protein'), ref.atomselect('protein'))

    assert np.allclose(tmscore, expectedTMscore)
    assert np.allclose(rmsd, expectedRMSD)

    #
    # rhodopsin = Molecule('1F88')
    # d3r = Molecule('3PBL')
    # alnmol = sequenceStructureAlignment(rhodopsin, d3r)