""" python/ui/model_symmetry.py: CCP4MG Molecular Graphics Program Copyright (C) 2009 University of York, CCLRC Copyright (C) 2012 STFC This library is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License version 3, modified in accordance with the provisions of the license to address the requirements of UK law. You should have received a copy of the modified GNU Lesser General Public License along with this library. If not, copies may be downloaded from http://www.ccp4.ac.uk/ccp4license.php This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. """ from global_definitions import * import mmdb2 as mmdb import math class model_symmetry: initialise = { 'symmetry_mates' : [], 'MolDisp_neighb_selection' : '', 'MolDisp_neighb_radius' : 5.0, 'search_radius' : 5.0, 'search_ncells' : 2 } MolDisp_selection_menu = ['all atoms','trace'] MolDisp_selection_command = ['','catrace or monomers'] #----------------------------------------------------------------------- #----------------------------------------------------------------------- #----------------------------------------------------------------------- def __init__(self,MolData_name,**kw): #----------------------------------------------------------------------- #----------------------------------------------------------------------- #----------------------------------------------------------------------- self.MolData_name = MolData_name self.biomol_list = [] self.previous_colour = MGCOLOUR().next('',exclude=['black','white','red','blue']) self.firstSymModel = -1 self.symmetry_MolDisp_objects = [] self.symmetry_colours = [] self.symmetry_visible = [] self.symmetry_mates = [] self.pending_symmetry_mates = [] self.MolDisp_visible = 1 self.MolDisp_selection = model_symmetry.MolDisp_selection_menu[0] #----------------------------------------------------------------------- def version_update_mates(self,mates=[]): #----------------------------------------------------------------------- # Check that format is not pre2.0 without the 'x'/'n' as first element # indicating its xtallographic or non-xtallographic symop # And add a 'visible' item updated_mates = [] for item in mates: # Version 1 did not support 'n' (non-xtallographic) if not ['x','n'].count(item[0]): item.insert(0,'x') # Pre 1.115 had extra 1 on end of definition for visibility which # is now handled by separate parameter if item[0]=='x' and len(item)>5: item = item[0:6] updated_mates.append(item) return updated_mates #----------------------------------------------------------------------- def molHnd(self): #----------------------------------------------------------------------- return data(self.MolData_name).molHnd #----------------------------------------------------------------------- def check_symmetry(self,pdbfile=''): #----------------------------------------------------------------------- ''' Get symmetry info from PDB and create menu of symops_titles ''' self.symops_titles = [] # Text list of symops #self.symops_repr = [] # parsed list of symops nSymOps = self.molHnd().GetNumberOfSymOps() #print "nSymOps",nSymOps if nSymOps <= 0: return [1,'No symmetry operations defined for space group',102] ''' if pdbfile: self.symops_titles = readPDBRemark290(pdbfile) print 'model_symmetry.check_symmetry',pdbfile,self.symops_titles else: self.symops_titles = [] ''' if not self.symops_titles: for n in range(nSymOps): # This uses the symmetry set up from the syminfo.lib #self.symops_titles.append(self.molHnd().GetSymOpTitle(n,0,0,0)) # This is using the mmdb symmetry from ???? self.symops_titles.append(self.molHnd().GetSymOp(n)) ''' rv,rep = parseSymOp(self.symops_titles[n]) if not rv: self.symops_repr.append(rep) else: self.symops_repr.append( [ [1,0,0,0],[0,1,0,0],[0,0,1,0]]) ''' #print 'model_symmetry.check_symmetry',self.symops_titles, self.symops_repr return [0,''] #----------------------------------------------------------------------- def load_biomolecule(self,pdbfile=''): #----------------------------------------------------------------------- rv = readPDBRemark350(pdbfile) if not rv[0]:self.biomol_list = rv[1] #print 'load_biomolecule readPDBRemark350',pdbfile,self.biomol_list #----------------------------------------------------------------------- def nofCrystalSymops(self): #----------------------------------------------------------------------- return len(self.symops_titles) #--------------------------------------------------------------------- def getParams(self): #--------------------------------------------------------------------- return {'symmetry_MolDisp_objects' : self.symmetry_MolDisp_objects, 'symmetry_colours' : self.symmetry_colours, 'symmetry_visible' : self.symmetry_visible, 'MolDisp_visible' : self.MolDisp_visible, 'MolDisp_selection' : self.MolDisp_selection } #--------------------------------------------------------------------- def setParams(self,params): #--------------------------------------------------------------------- for key in ['symmetry_MolDisp_objects','symmetry_colours','symmetry_visible','MolDisp_visible','MolDisp_selection']: if params.has_key(key): setattr(self,key,params[key]) #--------------------------------------------------------------------- def getGuiParams(self): #--------------------------------------------------------------------- #print 'model_symmetry.getParams',self.symmetry_mates symop=[] tran_a=[] tran_b=[] tran_c=[] visible = [] for item in self.symmetry_mates: if item[0] == 'x': symop.append(self.symops_titles[item[1]]) tran_a.append(str(item[2])) tran_b.append(str(item[3])) tran_c.append(str(item[4])) pars = { 'symop':symop, 'tran_a':tran_a, 'tran_b':tran_b, 'tran_c':tran_c , 'model_num' : self.get_symmetry_models(), 'visible' : self.symmetry_visible, 'colour' : self.symmetry_colours } #print 'model_symmetry.getParams',pars return pars #--------------------------------------------------------------------- def getSymmetryMates(self): #--------------------------------------------------------------------- #return self.symmetry_mates mates = [] for item in self.symmetry_mates: if item[0] == 'x': mates.append(item) elif item[0] == 'n': mates.append(('n',item[1].getParams(),item[2])) #print 'getSymmetryMates',mates return mates #--------------------------------------------------------------------- def setSymmetryMates(self,sym_mates=[]): #--------------------------------------------------------------------- # Ensure its updated from pre-2.0 format #print 'model_symmetry.setSymmetryMates',sym_mates sym_mates = self.version_update_mates(sym_mates) updated_mates = [] biomol = None for item in sym_mates: if item[0] == 'n' and not isinstance(item[1],BiomoleculeDefinition): biomol = BiomoleculeDefinition(params=item[1]) updated_mates.append(['n',biomol,item[2]]) else: updated_mates.append(item) # Reset the symmetry colors replace = 1 if replace: self.symmetry_colours = [] self.symmetry_visible = [] self.previous_colour = MGCOLOUR().next('',exclude=['black','white','red','blue']) self.symmetry_colours.extend(self.assignColours(updated_mates)) self.symmetry_visible.extend(self.assignVisible(updated_mates)) self.buildSymmetryMates(updated_mates,replace=replace) if biomol: # FIXME - A preference here for aromaticize. print "FIXME (model_symmetry.py) - A preference here for aromaticize." output = self.molHnd().GetMolBonds() #print 'setSymmetryMates GetMolBonds',output #for nm in range(1,self.molHnd().GetNumberOfModels()+1): #print 'model',nm,self.molHnd().GetModel(nm).GetNumberOfAtoms(0) #--------------------------------------------------------------------- def getPendingSymmetryMate(self,symop): #--------------------------------------------------------------------- # symop could be a sym_id of form '1_555' or list of form ['x',rot,tranx,trany,tranz] #print 'model_symmetrry.getPendingSymmetryMate',symop import types,string # Convert a sym_id to mg form if isinstance(symop,types.StringType): if symop == '1_555': return 1 split = string.split(symop,'_') if not len(split)==2 or len(split[0])<1 or len(split[0])>2 or len(split[1])!=3: return -1 try: symop = ['x',int(split[0]),int(split[1][0])-5,int(split[1][1])-5,int(split[1][2])-5] #print 'model_symmetry.getPendingSymmetryMate symop',symop except: return -1 if self.symmetry_mates.count(symop): return self.symmetry_mates.index(symop) + self.firstSymModel elif self.pending_symmetry_mates.count(symop): return self.pending_symmetry_mates.index(symop) + len(self.symmetry_mates) + self.firstSymModel else: self.pending_symmetry_mates.append(symop) return len(self.pending_symmetry_mates) + len(self.symmetry_mates) + self.firstSymModel - 1 #--------------------------------------------------------------------- def makePendingSymmetryMate(self,replace=0): #--------------------------------------------------------------------- #print 'model_symmetry.makePendingSymmetryMate',self.pending_symmetry_mates if self.pending_symmetry_mates: rv = self.buildSymmetryMates(self.pending_symmetry_mates,replace=replace) else: rv = 0 self.pending_symmetry_mates = [] return rv #--------------------------------------------------------------------- def setGuiParams(self,params,append=0): #--------------------------------------------------------------------- #print 'model_symmetry.setGuiParams',params mates = [] visible = [] colours = [] for ii in range(len(params['symop'])): if self.symops_titles.count(params['symop'][ii]): mates.append( [ 'x',self.symops_titles.index(params['symop'][ii]), int(params['tran_a'][ii]),int(params['tran_b'][ii]),int(params['tran_c'][ii])]) visible.append(params['visible'][ii]) colours.append(params['colour'][ii]) self.symmetry_visible = visible self.symmetry_colours = colours #print 'setGuiParams visible' , self.symmetry_visible self.buildSymmetryMates(mates,1) #--------------------------------------------------------------------- def findSymmetryNeighbours(self,selection='all',radius=5.0,xtl_limit=2): #--------------------------------------------------------------------- # Sort out if only want neighbours of limited selection of atoms import string neighb_selHnd = -1 if not ['all',''].count(string.strip(selection)): rv=data(self.MolData_name).parse_selection(command='{'+selection+'} and /1/*/*.*/*;*') if rv[0]: return [2,'Error parsing selection: '+selection] else: neighb_selHnd = rv[1] #print "findSymmetryNeighbours neighb_selHnd",selection,neighb_selHnd # Loop over all possible neighbouring symm mates and test a,b,c,alpha,beta,gamme = self.molHnd().GetCellInfo() allSelHnd = self.molHnd().NewSelection() self.molHnd().SelectAtoms(allSelHnd, 0,"*",mmdb.ANY_RES,"*",mmdb.ANY_RES,"*","*","*","*","*",mmdb.SKEY_NEW ) centre = self.molHnd().CentreOfCoordinatesAsCartesian(allSelHnd) xshifts = int(math.ceil(radius/a)); yshifts = int(math.ceil(radius/b)); zshifts = int(math.ceil(radius/c)); xoff = abs(int((centre.get_x())/a))+1; yoff = abs(int((centre.get_y())/b))+1; zoff = abs(int((centre.get_z())/c))+1; xtl_range_x = range(-xshifts-xoff,xshifts+2+xoff) xtl_range_y = range(-yshifts-yoff,yshifts+2+yoff) xtl_range_z = range(-zshifts-zoff,zshifts+2+zoff) symmetry_mates = [] interrupt = 0 n_new = 0 n_error = [] #xtl_range = range(-xtl_limit,xtl_limit+1) for i in xtl_range_x : if interrupt: break for j in xtl_range_y: if interrupt: break #MGGUI().UpdateGUIWindow('symm_mates_'+self.name, \ # [['PARAM','progress', \ # ['INITIAL','Testing tranform i='+str(i)+' j='+str(j)+' ....']]]) for k in xtl_range_z: if interrupt: break for nsym in range(0,len(self.symops_titles)): if interrupt: break if [nsym,i,j,k] != [0,0,0,0]: ncontacts = self.molHnd().IfSymmetryNeighbours(neighb_selHnd,\ 1,nsym,i,j,k,radius) #print "IfSymmetryNeighbour",tran,ncontacts if ncontacts>0: n_new = n_new + 1 symmetry_mates.append(['x',nsym,i,j,k]) elif ncontacts < 0: n_error.append(ncontacts) #interrupt=MGGUI().GetInterruptStatus() if interrupt: break return [interrupt,symmetry_mates] #---------------------------------------------------------------------- def buildSymmetryMates(self,symmetry_mates=[],replace=0): #---------------------------------------------------------------------- import pygl_coord #print "build_symmetry_mates symmetry_mates",symmetry_mates # First call of this method needs to find out how many models exist in the PDB if self.firstSymModel<0: self.firstSymModel = self.molHnd().GetNumberOfModels()+1 data(self.MolData_name).resetSelHandles() t = self.molHnd().GetTransform() oldM = pygl_coord.matrix(4,4,t) #print t self.molHnd().UnSetTransform(True) nBuilt = 0 if replace: if len(symmetry_mates) > 0 and len(self.symmetry_mates) > 0: # If necessary reset the transformation for n in range(min(len(symmetry_mates),len(self.symmetry_mates))): #print 'build_symmetry_mates',n,symmetry_mates[n] == self.symmetry_mates[n] if symmetry_mates[n] != self.symmetry_mates[n]: undo_transform = self.symmetry_mates.pop(n) transform = symmetry_mates[n] self.resetTransform(n+self.firstSymModel,undo_transform,transform) self.symmetry_mates.insert(n,transform) nBuilt = nBuilt + 1 if len(symmetry_mates) < len(self.symmetry_mates): # Delete some models for n in range(len(symmetry_mates),len(self.symmetry_mates)): #print 'buildSymmetryMates delete',n rc = self.molHnd().DeleteModel("/"+str(n+self.firstSymModel)) #rc = self.molHnd.DeleteModel(n) #print "deleting",n,self.symmetry_models[n],rc for n in range(len(symmetry_mates),len(self.symmetry_mates)): self.symmetry_mates.pop() nBuilt = nBuilt + 1 if (not replace) or (len(symmetry_mates) > len(self.symmetry_mates)): # Add some models first = 0 if replace: first = len(self.symmetry_mates) for item in symmetry_mates[first:]: self.makeSymmetryMate(item) nBuilt = nBuilt + 1 # Build bonds etc. #print 'buildSymmetryMates nBuilt',nBuilt if nBuilt: self.molHnd().FinishStructEdit() self.molHnd().MakeAllSelIndexes() data(self.MolData_name).GetMolBonds() #print 'done buildSymmetrymates',nBuilt,self.symmetry_mates, self.symmetry_colours self.molHnd().SetTransform(oldM,False) if nBuilt: data(self.MolData_name).redraw_all_objects() return nBuilt #------------------------------------------------------------------ def assignVisible(self,symmetry_mates=[]): #------------------------------------------------------------------ # Work out the mapping from existing mates list to new mates list visibility = [] for item in symmetry_mates: if self.symmetry_mates.count(item) and self.symmetry_mates.index(item)self.nofCrystalSymops(): return '' v = [] for ii in range(3): jj = utils.safeInt(ss[1][ii],-1) if jj<0: return -1 v.append(jj-5) return ['x',symop-1,v[0],v[1],v[2]] #-------------------------------------------------------------------- def modelSymid(self,model): #-------------------------------------------------------------------- symid = '' import types,utils if isinstance(model,types.StringType): splitmod = model.split() if len(splitmod): model = utils.safeInt(splitmod[0],-1) else: model = -1 if model == 1: symid = '1_000' else: ii = model-self.firstSymModel if ii>=0 and ii=0 and ii=0 and op[1]=0 and ii=0 and ii=0 and visible!=self.MolDisp_visible: self.MolDisp_visible = visible if self.MolDisp_visible: self.create_MolDisp() self.colour_MolDisp() else: self.delete_MolDisp() redraw = 1 if redraw and self.MolDisp_visible: self.select_MolDisp() return redraw #-------------------------------------------------------------------- #-------------------------------------------------------------------- #-------------------------------------------------------------------- class BiomoleculeTransform: #-------------------------------------------------------------------- #-------------------------------------------------------------------- #-------------------------------------------------------------------- IDENTITY_ROT = [ 1.0,0.0,0.0, 0.0,1.0,0.0, 0.0,0.0,1.0 ] IDENTITY_TRAN = [ 0.0,0.0,0.0 ] IDENTITY_TOLERANCE = 0.001 def __init__(self,params={},**kw): params.update(kw) self.setParams(params) def getParams(self): return { 'id' : self.id, 'chains' : self.chains, 'rotation' : self.rotation, 'translation' : self.translation } def setParams(self,params={}): self.id = params.get('id','') self.chains = params.get('chains',[]) self.chains.sort() self.rotation = params.get('rotation',BiomoleculeTransform.IDENTITY_ROT) self.translation = params.get('translation',BiomoleculeTransform.IDENTITY_TRAN) def setTransform(self,id,rotation,translation): self.id = id self.rotation = rotation self.translation = translation def isIdentity(self): for ii in range(9): if abs(self.rotation[ii]-BiomoleculeTransform.IDENTITY_ROT[ii])>BiomoleculeTransform.IDENTITY_TOLERANCE: return 0 for ii in range(3): if abs(self.translation[ii])>BiomoleculeTransform.IDENTITY_TOLERANCE: return 0 return 1 def getTMatrix(self): import mmdb2 as mmdb import pygl_coord v = pygl_coord.doublea(16) v[0] = self.rotation[0] v[1] = self.rotation[1] v[2] = self.rotation[2] v[3] = 0.0 v[4] = self.rotation[3] v[5] = self.rotation[4] v[6] = self.rotation[5] v[7] = 0.0 v[8] = self.rotation[6] v[9] = self.rotation[7] v[10] = self.rotation[8] v[11] = 0.0 v[12] = self.translation[0] v[13] = self.translation[1] v[14] = self.translation[2] v[15]= 1.0 #print 'BiomoleculeTransform.getTMatrix' #print v[0],v[1],v[2],v[3] #print v[4],v[5],v[6],v[7] #print v[8],v[9],v[10],v[11] #print v[12],v[13],v[14],v[15] return v def setBIOMT(self,lines): #print 'addBIOMT',lines if len(lines)!=3: print 'ERROR in BiomoleculeTransforms.addBIOMT input lines',len(lines) else: import string,utils rotmat = [] translate = [] for line in lines: fields = string.split(line) for ii in range(1,4): rotmat.append(utils.safeFloat(fields[ii],0.0)) translate.append(utils.safeFloat(fields[4],0.0)) self.setTransform(fields[0],rotmat,translate) def printInfo(self): print 'CHAINS: ',self.chains print 'id',self.id,'rotation: ',self.rotation,'translation',self.translation #-------------------------------------------------------------------- #-------------------------------------------------------------------- #-------------------------------------------------------------------- class BiomoleculeDefinition: #-------------------------------------------------------------------- #-------------------------------------------------------------------- #-------------------------------------------------------------------- def __init__(self,params={},**kw): params.update(kw) self.setParams(params) def setParams(self,params): self.label = params.get('label','') self.comments = params.get('comments','') self.transform_list = [] for item in params.get('transform_list',[]): if isinstance(item,BiomoleculeTransform): self.transform_list.append(item) else: self.transform_list.append(BiomoleculeTransform(params=item)) def getParams(self): transform_list = [] for item in self.transform_list: transform_list.append(item.getParams()) return { 'label' : self.label, 'comments' : self.comments, 'transform_list' : transform_list } def setData(self,biomolecule='',transform_list =[],comments=[]): self.label = biomolecule self.transform_list = transform_list self.comments = comments def printInfo(self): print 'BIOMOLECULE',self.label print 'comments', self.comments for item in self.transform_list: item.printInfo() def getNofTransforms(self): return len(self.transform_list) def isIdentity(self,nt): if nt<0 or nt>len(self.transform_list): return -1 else: return self.transform_list[nt].isIdentity() def getRotation(self,nt): if nt<0 or nt>len(self.transform_list): return None else: return self.transform_list[nt].rotation def getTranslation(self,nt): if nt<0 or nt>len(self.transform_list): return None else: return self.transform_list[mt].translation def getChains(self,nt): if nt>0 and ntlen(self.transform_list): return None else: return self.transform_list[nt].getTMatrix() def noTransformNeeded(self,chainsInModel): # If there is one (identity) transform applied to all chains # then this is identity # input chainsInModel is list of chains in PDB import string if self.getNofTransforms()>1: return 0 chn_list = [] for item in chainsInModel: chn_list.append(string.strip(item,'/')) if len(self.getChains(0)) == 0 or self.getChains(0) == chn_list.sort() : return 1 else: return 0 def readPDBRemark350(filename): try: f = open(filename) text = f.read() f.close() except: return [1,'Error reading file'] import string,utils lines = string.split(text,'\n') nLines = len(lines) #print 'readPDBRemark350 nLines',nLines biomol_def_list = [] biomolecule = 0 comments = [] transform_list = [] ii = 0 foundrem350 = 0 while ii < nLines-1: ii = ii +1 line = lines[ii] #print ii,'line',line if line[0:10] == 'REMARK 350': foundrem350 = foundrem350+1 if string.find(line,'BIOMOLECULE:')>=0: if transform_list: obj = BiomoleculeDefinition() rv = obj.setData(biomolecule,transform_list,comments) if not rv: biomol_def_list.append(obj) biomolecule = utils.safeInt(string.strip(line[string.find(line,'BIOMOLECULE:')+12]),0) comments = [] elif string.find(line,'CHAINS:')>=0: chains = [] while (string.find(lines[ii],'CHAINS:')>=0): #print 'reading chains' chns = string.split(lines[ii][string.find(line,'CHAINS:')+7:],',') #print 'readPDBRemark350',chns for item in chns: chains.append(string.strip(item)) ii = ii + 1 ii = ii - 1 elif string.find(line,'BIOMT1')>=0: #print 'reading BIOMT' tt = [] for jj in range(3): tt.append(lines[ii+jj][string.find(line,'BIOMT')+6:]) transform_list.append(BiomoleculeTransform(chains=chains)) transform_list[-1].setBIOMT(tt) ii = ii + 2 elif biomolecule and not transform_list: comments.append(string.strip(line[10:])) else: #print 'not rem350',foundrem350 if foundrem350>0: break #print 'out of loop ii',ii if transform_list: obj = BiomoleculeDefinition() rv = obj.setData(biomolecule,transform_list,comments) if not rv: biomol_def_list.append(obj) #for item in biomol_def_list: # item.printInfo() return [0,biomol_def_list] def readPDBRemark290(filename): import string,re try: f = open(filename) text = f.read() f.close() except: return [1,'Error reading file'] lines = string.split(text,'\n') nLines = len(lines) #print 'readPDBRemark290',nLines hits = [] foundrem290 = 0 ii = 0 while ii < nLines-1: ii = ii + 1 line = lines[ii] print line[0:10] if line[0:10] == 'REMARK 290': foundrem290 = foundrem290 + 1 print ii,'line',line m = re.match(r"([ ]*)(\d+)([ ]*)(.*)",line[11:]) if m: gps = m.groups() hits.append([gps[1],string.strip(gps[3])]) else: if foundrem290: break symops = [] for item in hits: symops.append(item[1]) return symops def parseSymOp(symop=''): import re chr_xyz = symop.split(',') #print 'chr_xyz',chr_xyz if len(chr_xyz) != 3 : return [1,'Symop should have three comma-separated fields'] rep_xyz = [] for ch0 in chr_xyz: ch = ch0.lower().strip() b = re.split('([+,-])',ch) if b[0] == '': del b[0] elif not ['+','-'].count(b[0]): b.insert(0,'+') rep = [0,0,0,0] ii = 0 #print 'parseSymOp b',b while ii+1