Skip to content

Commit

Permalink
updated sequenceID to start from 0. updated all it's usages. Fixed bu…
Browse files Browse the repository at this point in the history
…g in Molecule where writing files with MDtraj didn't delete the temp files
  • Loading branch information
Stefan Doerr committed Feb 1, 2017
1 parent 21040be commit 933f4f5
Show file tree
Hide file tree
Showing 5 changed files with 21 additions and 10 deletions.
4 changes: 2 additions & 2 deletions htmd/builder/amber.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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))
Expand Down
4 changes: 2 additions & 2 deletions htmd/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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
Expand Down
17 changes: 14 additions & 3 deletions htmd/molecule/molecule.py
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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))
Expand Down
4 changes: 2 additions & 2 deletions htmd/molecule/util.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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:
Expand Down
2 changes: 1 addition & 1 deletion htmd/projections/metricsecondarystructure.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit 933f4f5

Please sign in to comment.