From 933f4f5b3bcc6553e9750f5241fb5e9f67736e2a Mon Sep 17 00:00:00 2001 From: Stefan Doerr Date: Wed, 1 Feb 2017 16:00:18 +0100 Subject: [PATCH] updated sequenceID to start from 0. updated all it's usages. Fixed bug in Molecule where writing files with MDtraj didn't delete the temp files --- htmd/builder/amber.py | 4 ++-- htmd/model.py | 4 ++-- htmd/molecule/molecule.py | 17 ++++++++++++++--- htmd/molecule/util.py | 4 ++-- htmd/projections/metricsecondarystructure.py | 2 +- 5 files changed, 21 insertions(+), 10 deletions(-) diff --git a/htmd/builder/amber.py b/htmd/builder/amber.py index aa576724d..d351319c0 100644 --- a/htmd/builder/amber.py +++ b/htmd/builder/amber.py @@ -215,7 +215,7 @@ def build(mol, ff=None, topo=None, param=None, prefix='structure', outdir='./bui if len(disulfide) != 0: for d in disulfide: # Convert to stupid amber residue numbering - uqseqid = sequenceID((mol.resid, mol.insertion, mol.segid)) + mol.resid[0] - 1 + uqseqid = sequenceID((mol.resid, mol.insertion, mol.segid)) + mol.resid[0] 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 @@ -262,7 +262,7 @@ def build(mol, ff=None, topo=None, param=None, prefix='structure', outdir='./bui f.write('# Adding disulfide bonds\n') for d in disulfide: # Convert to stupid amber residue numbering - uqseqid = sequenceID((mol.resid, mol.insertion, mol.segid)) + mol.resid[0] - 1 + uqseqid = sequenceID((mol.resid, mol.insertion, mol.segid)) + mol.resid[0] 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))])) f.write('bond mol.{}.SG mol.{}.SG\n'.format(uqres1, uqres2)) diff --git a/htmd/model.py b/htmd/model.py index 011db2e1a..f67607a75 100644 --- a/htmd/model.py +++ b/htmd/model.py @@ -546,7 +546,7 @@ def viewStates(self, states=None, statetype='macro', protein=None, ligand=None, viewer.send('start_sscache') def _viewStatesNGL(self, states, statetype, protein, ligand, mols, numsamples, gui=False): - from htmd.builder.builder import sequenceID + from htmd.molecule.util import sequenceID if states is None: states = range(self.macronum) if isinstance(states, int): @@ -568,7 +568,7 @@ def _viewStatesNGL(self, states, statetype, protein, ligand, mols, numsamples, g if ligand: mol = ref.copy() mol.remove(ligand, _logger=False) - mol.coords = np.atleast_3d(mol.coords[:, :, 0]) + mol.dropFrames(keep=0) mols[i].filter(ligand, _logger=False) mols[i].set('chain', '{}'.format(s)) tmpcoo = mols[i].coords diff --git a/htmd/molecule/molecule.py b/htmd/molecule/molecule.py index e205aa82f..3c7cd6716 100644 --- a/htmd/molecule/molecule.py +++ b/htmd/molecule/molecule.py @@ -222,6 +222,11 @@ def frame(self, value): raise NameError("Frame index out of range. Molecule contains {} frame(s). Frames are 0-indexed.".format(self.numFrames)) self._frame = value + @property + def numResidues(self): + from htmd.molecule.util import sequenceID + return sequenceID((self.resid, self.insertion, self.chain)) + def insert(self, mol, index, collisions=False, coldist=1.3): """Insert the contents of one molecule into another at a specific index. @@ -973,21 +978,27 @@ def view(self, sel=None, style=None, color=None, guessBonds=True, viewer=None, h self.write(xtc) # Call the specified backend + retval = None if viewer is None: from htmd.config import _config viewer = _config['viewer'] if viewer.lower() == 'notebook': - return self._viewMDTraj(psf, xtc) + retval = self._viewMDTraj(psf, xtc) elif viewer.lower() == 'vmd': self._viewVMD(psf, xtc, viewerhandle, name, guessBonds) + #retval = viewerhandle elif viewer.lower() == 'ngl' or viewer.lower() == 'webgl': - return self._viewNGL(gui=gui) + retval = self._viewNGL(gui=gui) else: + os.remove(xtc) + os.remove(psf) raise ValueError('Unknown viewer.') # Remove temporary files os.remove(xtc) os.remove(psf) + if retval is not None: + return retval def _viewVMD(self, psf, xtc, vhandle, name, guessbonds): if name is None: @@ -1186,7 +1197,7 @@ def sequence(self, oneletter=True): '?EEI' """ - from htmd.builder.builder import sequenceID + from htmd.molecule.util import sequenceID prot = self.atomselect('protein') segs = np.unique(self.segid[prot]) increm = sequenceID((self.resid, self.insertion, self.chain)) diff --git a/htmd/molecule/util.py b/htmd/molecule/util.py index bd1d102c7..c322be80c 100644 --- a/htmd/molecule/util.py +++ b/htmd/molecule/util.py @@ -44,7 +44,7 @@ def molTMscore(mol, ref, selCAmol, selCAref): from htmd.molecule.molecule import _residueNameTable def calculateVariables(currmol): - res = sequenceID((currmol.resid, currmol.insertion, currmol.segid, currmol.chain)) - 1 + res = sequenceID((currmol.resid, currmol.insertion, currmol.segid, currmol.chain)) caidx = currmol.name == 'CA' res = np.unique(res) reslen = len(res) @@ -144,7 +144,7 @@ def sequenceID(field, prepend=None): else: seq = np.empty(fieldlen, dtype=object) - c = int(1) + c = int(0) if prepend is None: seq[0] = c else: diff --git a/htmd/projections/metricsecondarystructure.py b/htmd/projections/metricsecondarystructure.py index 41e8a7e38..856e92feb 100644 --- a/htmd/projections/metricsecondarystructure.py +++ b/htmd/projections/metricsecondarystructure.py @@ -64,7 +64,7 @@ def _calcarrays(self, mol): mol = mol.copy() mol.filter(self.sel, _logger=False) - residues = sequenceID((mol.resid, mol.chain, mol.insertion)) - 1 # Start at 0 + residues = sequenceID((mol.resid, mol.chain, mol.insertion)) backbone = mol.atomselect('backbone') ca_indices = np.where(mol.name == 'CA')[0].astype(np.int32)