#A* ------------------------------------------------------------------- #B* This file contains source code for the PyMOL computer program #C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific. #D* ------------------------------------------------------------------- #E* It is unlawful to modify or remove this copyright notice. #F* ------------------------------------------------------------------- #G* Please see the accompanying LICENSE file for further information. #H* ------------------------------------------------------------------- #I* Additional authors of this source file include: #-* #-* #-* #Z* ------------------------------------------------------------------- import chempy import copy from cpv import * import operator 00020 class Base: #------------------------------------------------------------------------------ def update_index(self): if chempy.feedback['verbose']: print " "+str(self.__class__)+": updating indexes..." c = 0 self.index = {} idx = self.index for a in self.atom: idx[id(a)] = c c = c + 1 #------------------------------------------------------------------------------ def get_residues(self): list = [] if self.nAtom: last = self.atom[0] c = 0 start = 0 for a in self.atom: if not a.in_same_residue(last): list.append((start,c)) start = c last = a c = c + 1 if (c-start>1): list.append((start,c)) return list #------------------------------------------------------------------------------ def get_coord_list(self): lst = [] for a in self.atom: lst.append(a.coord) return lst #------------------------------------------------------------------------------ def get_mass(self): sm = 0.0 for a in self.atom: sm = sm + a.get_mass() return sm #------------------------------------------------------------------------------ def get_nuclear_charges(self): '''Return the sum of nuclear charges of all atoms in a molecule.''' sm = 0 for a in self.atom: sm = sm + a.get_number() return sm #------------------------------------------------------------------------------ def list(self): for a in self.atom: print a.symbol, a.name, a.coord for a in self.bond: print a.index #------------------------------------------------------------------------------ def get_implicit_mass(self): # mass calculation for implicit models valence = [0]*len(self.atom) implicit = [0]*len(self.atom) for a in self.bond: ai0 = a.index[0] ai1 = a.index[1] valence[ai0] = valence[ai0] + a.order valence[ai1] = valence[ai1] + a.order c = 0 for a in self.atom: valence[c] = valence[c] - a.formal_charge implicit[c] = a.get_implicit_valence()[valence[c]] c = c + 1 h_count = reduce(operator.__add__,implicit) hydrogen = chempy.Atom() hydrogen.symbol='H' return self.get_mass()+hydrogen.get_mass()*h_count #------------------------------------------------------------------------------ def assign_names(self,preserve=1): dct = {} cnt = {} if preserve: for a in self.atom: if a.has('name'): dct[a['name']] = 1 else: for a in self.atom: if hasattr(a,'name'): del a.name for a in self.atom: if not a.has('name'): if not cnt.has_key(a.symbol): c = 1 else: c = cnt[a.symbol] while 1: nam = a.symbol+str(c) c = c + 1 if not dct.has_key(nam): break cnt[a.symbol]=c a.name=nam dct[nam]=1 #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ 00131 class Indexed(Base): attr_value = { 'nAtom' : compile('len(self.atom)','Indexed','eval'), 'nBond' : compile('len(self.bond)','Indexed','eval'), } def __getattr__(self,attr): if Indexed.attr_value.has_key(attr): return eval(Indexed.attr_value[attr]) else: raise AttributeError(attr) #------------------------------------------------------------------------------ def __init__(self): self.reset() #------------------------------------------------------------------------------ def reset(self): self.index = None self.molecule = chempy.Molecule() self.atom = [] self.bond = [] #------------------------------------------------------------------------------ def get_min_max(self): if len(self.atom): mn = copy.deepcopy(self.atom[0].coord) mx = copy.deepcopy(self.atom[0].coord) for a in self.atom: ac = a.coord if mn[0]>ac[0]: mn[0]=ac[0] if mn[1]>ac[1]: mn[1]=ac[1] if mn[2]>ac[2]: mn[2]=ac[2] if mx[0]<ac[0]: mx[0]=ac[0] if mx[1]<ac[1]: mx[1]=ac[1] if mx[2]<ac[2]: mx[2]=ac[2] return [mn,mx] else: return [[0.0,0.0,0.0],[0.0,0.0,0.0]] #------------------------------------------------------------------------------ def merge(self,other): # steals atom objects from 'other' and resets 'other' if chempy.feedback['actions']: print " "+str(self.__class__)+": merging models..." nAtom=self.nAtom self.atom.extend(other.atom) for b in other.bond: b.index[0]=b.index[0]+nAtom b.index[1]=b.index[1]+nAtom self.bond.append(b) other.reset() if self.index: self.update_index() #-------------------------------------------------------------------------------- def delete_atom(self,index): if chempy.feedback['atoms']: print " "+str(self.__class__)+": deleting atom %d." % index nAtom=self.nAtom # update index if it exists if self.index: idx = self.index for k in idx.keys(): if idx[k] > index: idx[k] = idx[k] - 1 del idx[id(self.atom[index])] # delete atom del self.atom[index] # delete bonds associated with this atom nBond = len(self.bond) templist = [] for i in range(nBond): if index in self.bond[i].index: templist.append(i) for i in range(len(templist)): j = templist[i] - i del self.bond[j] # re-index bond table for b in self.bond: if b.index[0] > index: b.index[0] = b.index[0] - 1 if b.index[1] > index: b.index[1] = b.index[1] - 1 #-------------------------------------------------------------------------------- def delete_list(self,list): # delete a list of indexed atoms if chempy.feedback['atoms']: print " "+str(self.__class__)+": deleting atoms %s." % str(list) nAtom=self.nAtom lrev = copy.deepcopy(list) lrev.sort() lrev.reverse() # generate cross-reference tables o2n = {} # old to new if len(lrev): nxt = lrev.pop() else: nxt = None shft = 0 for i in range(nAtom): if i == nxt: o2n[i]=-1 if len(lrev): nxt = lrev.pop() else: nxt = None shft = shft - 1 else: ixs = i + shft o2n[i] = ixs if shft: # delete atoms new_atom = [] for i in range(nAtom): if o2n[i]>=0: new_atom.append(self.atom[i]) self.atom = new_atom # delete bonds new_bond = [] for b in self.bond: b0 = b.index[0] b1 = b.index[1] if (o2n[b0]>=0) and (o2n[b1]>=0): b.index[0] = o2n[b0] b.index[1] = o2n[b1] new_bond.append(b) self.bond = new_bond # update index if it exists if self.index: self.index = {} idx = self.index i = 0 for a in self.atom: idx[id(a)] = i i = i + 1 #------------------------------------------------------------------------------ def insert_atom(self,index,atom): if chempy.feedback['atoms']: print " "+str(self.__class__)+': inserting atom "%s" before %d.' % ( atom.name,index) nAtom=self.nAtom self.atom.insert(index,atom) # re-index bond table for b in self.bond: if b.index[0] >= index: b.index[0] = b.index[0] + 1 if b.index[1] >= index: b.index[1] = b.index[1] + 1 # update index if it exists if self.index: idx = self.index for k in idx.keys(): if idx[k] >= index: idx[k] = idx[k] + 1 idx[id(atom)] = index #------------------------------------------------------------------------------ def index_atom(self,atom): c = 0 id_at = id(atom) for a in self.atom: if id(a)==id_at: return c c = c + 1 return -1 #------------------------------------------------------------------------------ def add_atom(self,atom): if chempy.feedback['atoms']: print " "+str(self.__class__)+': adding atom "%s".' % atom.name self.atom.append(atom) index = self.nAtom - 1 if self.index: self.index[id(atom)] = index return index #------------------------------------------------------------------------------ def add_bond(self,bond): if chempy.feedback['bonds']: print " "+str(self.__class__)+": adding bond (%d,%d)." % \ (bond.index[0],bond.index[1]) self.bond.append(bond) #------------------------------------------------------------------------------ def remove_bond(self,index): if chempy.feedback['bonds']: print " "+str(self.__class__)+": removing bond %d." % index nBond=len(self.Bond) del self.bond[index] #------------------------------------------------------------------------------ def convert_to_connected(self): if chempy.feedback['verbose']: print " "+str(self.__class__)+": converting to connected model..." model = Connected() model.molecule = self.molecule model.atom = self.atom model.bond = [] model.index = None for a in model.atom: model.bond.append([]) for b in self.bond: model.bond[b.index[0]].append(b) # note two refs to same object model.bond[b.index[1]].append(b) # note two refs to same object self.reset() return model #------------------------------------------------------------------------------ def from_molobj(self,molobj): self.reset() mol = self.molecule if len(molobj.title): mol.title = molobj.title if len(molobj.comments): mol.comments = molobj.comments mol.chiral = molobj.chiral mol.dim_code = molobj.dimcode for a in molobj.atom: at = chempy.Atom() at.symbol = a.symbol at.name = a.name if a.resn != chempy.Atom.defaults['resn']: at.resn = a.resn if a.resn_code != chempy.Atom.defaults['resn_code']: at.resn_code = a.resn_code at.resi = a.resi at.resi_number = a.resi_number at.b = a.b at.q = a.q at.alt = a.alt at.hetatm = a.hetatm if a.segi != chempy.Atom.defaults['segi']: at.segi = a.segi if a.chain != chempy.Atom.defaults['chain']: at.chain = a.chain at.color_code = a.color_code at.coord = a.coord at.formal_charge = a.formal_charge at.partial_charge = a.partial_charge if a.numeric_type != -99: at.numeric_type = a.numeric_type if a.text_type != 'UNKNOWN': at.text_type = a.text_type at.stereo = a.stereo if hasattr(a,'flags'): at.flags = a.flags if hasattr(a,'vdw'): at.vdw = a.vdw self.atom.append(at) for b in molobj.bond: bnd = chempy.Bond() bnd.index = [b.atom[0],b.atom[1]] bnd.order = b.order bnd.stereo = b.stereo self.bond.append(bnd) #------------------------------------------------------------------------------ def sort(self): if chempy.feedback['verbose']: print " "+__name__+": sorting..." if not self.index: self.update_index() old_index = self.index self.atom.sort() self.update_index() xref = {} new_index = self.index for a in new_index.keys(): xref[old_index[a]] = new_index[a] for b in self.bond: b.index[0] = xref[b.index[0]] b.index[1] = xref[b.index[1]] del old_index del xref #------------------------------------------------------------------------------ def get_internal_tuples(self): # generates raw atom sets needed to construct an internal coordinate # description of the molecule model = self # get a connected version too cmodel = copy.deepcopy(model).convert_to_connected() center = [0.0,0.0,0.0] nAtom = model.nAtom to_go = nAtom done = {} if to_go<3: z_set = [(0),(1,0)] else: # get center of molecule for a in model.atom: center = add(center,a.coord) center = scale(center,1.0/nAtom) # find most central multivalent atom min_a = -1 c = 0 for a in model.atom: if len(cmodel.bond[c])>1: # must have at least two neighbors d = distance(a.coord,center) if min_a < 0: min_d = d min_a = c elif d<min_d: min_d=d min_a=c c = c + 1 fst = min_a # make this our first atom z_set = [( fst, )] done[fst] = 1 to_go = to_go - 1 # for the second atom, try to choose different multivalent neighbor nxt = -1 for b in cmodel.bond[fst]: nbr = b.index[0] if nbr == fst: nbr = b.index[1] if len(cmodel.bond[nbr])>1: nxt = nbr break # safety, choose any neighbor if nxt<0: nbr = b.index[0] if nbr == fst: nbr = b.index[1] nxt = nbr z_set.append((nxt,fst)) done[nxt] = 1 to_go = to_go - 1 # for the third atom, choose a different multivalent neighbor trd = -1 for b in cmodel.bond[fst]: nbr = b.index[0] if nbr == fst: nbr = b.index[1] if len(cmodel.bond[nbr])>1: if not done.has_key(nbr): trd = nbr break # safety, choose any unchosen neighbor if trd<0: for b in cmodel.bond[fst]: nbr = b.index[0] if nbr == fst: nbr = b.index[1] if not done.has_key(nbr): trd = nbr break z_set.append((trd,fst,nxt)) done[trd] = 1 result = 1 to_go = to_go - 1 if to_go: # now find all torsions in the molecule tors = {} for b in model.bond: # use bond as center of torsion a1 = b.index[0] a2 = b.index[1] for c in cmodel.bond[a1]: a0 = c.index[0] if a0 not in (a1,a2): # outside atom for d in cmodel.bond[a2]: a3 = d.index[0] if a3 not in (a0,a1,a2): # outside atom if a0 < a3: to = (a0,a1,a2,a3) else: to = (a3,a2,a1,a0) tors[to] = 1 a3 = d.index[1] if a3 not in (a0,a1,a2): # outside atom if a0 < a3: to = (a0,a1,a2,a3) else: to = (a3,a2,a1,a0) tors[to] = 1 a0 = c.index[1] if a0 not in (a1,a2): # outside atom for d in cmodel.bond[a2]: a3 = d.index[0] if a3 not in (a0,a1,a2): # outside atom if a0 < a3: to = (a0,a1,a2,a3) else: to = (a3,a2,a1,a0) tors[to] = 1 a3 = d.index[1] if a3 not in (a0,a1,a2): # outside atom if a0 < a3: to = (a0,a1,a2,a3) else: to = (a3,a2,a1,a0) tors[to] = 1 if len(tors.keys()): # choose remaining atoms based on existing atoms using torsion while to_go: for tor in tors.keys(): a0 = tor[0] a1 = tor[1] a2 = tor[2] a3 = tor[3] dh0 = done.has_key(a0) dh1 = done.has_key(a1) dh2 = done.has_key(a2) dh3 = done.has_key(a3) if ( (not dh0) and dh1 and dh2 and dh3 ): z_set.append((a0,a1,a2,a3)) done[a0] = 1 to_go = to_go - 1 elif ( dh0 and dh1 and dh2 and (not dh3) ): z_set.append((a3,a2,a1,a0)) done[a3] = 1 to_go = to_go - 1 else: # for molecules with no torsions (dichloromethane, etc.) # we have to generate torsions which include one virtual # bond for b in cmodel.bond[fst]: nbr = b.index[0] if nbr in [fst,nxt,trd]: nbr = b.index[1] if not done.has_key(nbr): z_set.append((nbr,trd,fst,nxt)) to_go = to_go - 1 done[nbr] = 1 return z_set #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ 00580 class Connected(Base): attr_value = { 'nAtom' : compile('len(self.atom)','Connected','eval'), } def __getattr__(self,attr): if Connected.attr_value.has_key(attr): return eval(Connected.attr_value[attr]) else: raise AttributeError(attr) #------------------------------------------------------------------------------ def __init__(self): self.reset() #------------------------------------------------------------------------------ def reset(self): self.index = None self.molecule = chempy.Molecule() self.atom = [] self.bond = [] #------------------------------------------------------------------------------ def convert_to_indexed(self): if chempy.feedback['verbose']: print " "+str(self.__class__)+": converting to indexed model..." indexed = Indexed() indexed.atom = self.atom indexed.molecule = self.molecule c = 0 for a in self.bond: for b in a: if b.index[0] == c: indexed.bond.append(b) c = c + 1 self.reset() return indexed #------------------------------------------------------------------------------ def insert_atom(self,index,atom): if chempy.feedback['atoms']: print " "+str(self.__class__)+': inserting atom "%s" before %d.' % ( atom.name,index) nAtom=self.nAtom self.atom.insert(index,atom) # re-index bond table for a in self.bonds: for b in a: if b.index[0] >= index: b.index[0] = b.index[0] + 1 if b.index[1] >= index: b.index[1] = b.index[1] + 1 # update index if it exists if self.index: idx = self.index for k in idx.keys(): if idx[k] >= index: idx[k] = idx[k] + 1 idx[id(atom)] = index #------------------------------------------------------------------------------ def delete_atom(self,index): if chempy.feedback['atoms']: print " "+str(self.__class__)+": deleting atom %d." % index nAtom=self.nAtom # update index if it exists if self.index: idx = self.index for k in idx.keys(): if idx[k] > index: idx[k] = idx[k] - 1 del idx[id(self.atom[index])] # delete atom del self.atom[index] # delete bonds associated with this atom nBond = len(self.bond) for a in self.bond: i = 0 templist = [] for b in a: if index in b.index: templist.append(i) i = i + 1 for i in range(len(templist)): j = templist[i] - i del a[j] # re-index bond table for b in self.bond: if b.index[0] > index: b.index[0] = b.index[0] - 1 if b.index[1] > index: b.index[1] = b.index[1] - 1 #------------------------------------------------------------------------------ def add_atom(self,atom): if chempy.feedback['atoms']: print " "+str(self.__class__)+': adding atom "%s".' % atom.name self.atom.append(atom) self.bond.append([]) index = self.nAtom - 1 if self.index: self.index[id(atom)] = index return index #------------------------------------------------------------------------------ def sort(self): if chempy.feedback['verbose']: print " "+__name__+": sorting..." if not self.index: self.update_index() old_index = self.index self.atom.sort() self.update_index() xref = {} new_index = self.index for a in new_index.keys(): xref[old_index[a]] = new_index[a] new_bond = [None] * self.nAtom c = 0 tmp_list = [] for a in self.bond: for b in a: if c==b.index[0]: tmp_list.append(b) new_bond[xref[c]] = a c = c + 1 for b in tmp_list: b.index[0] = xref[b.index[0]] b.index[1] = xref[b.index[1]] del self.bond self.bond = new_bond del old_index del xref