'''Adapted from helix bundle rotation code written by Nathalie for TMM@''' import MMTK from MMTK import * from MMTK.PDB import PDBConfiguration from MMTK.Proteins import Protein from MMTK.NormalModes import NormalModes from numpy import sqrt, add, array, argmax, argsort, argmin,linalg, dot from Scientific import * from Scientific.IO.ArrayIO import writeArray from Scientific.Visualization import VMD import sys #from rpy import * #import PlotTableR, util def massWeightedNormalVector(v): return v/sqrt(v.massWeightedDotProduct(v)) filename=sys.argv[1].replace(".mode","") modes = MMTK.load(sys.argv[1]) #should have the input modes file for the TIM barrel structure in question universe = modes.universe protein = universe.protein protein.normalizeConfiguration('Ir') protein.writeToFile(sys.argv[2]+"_"+sys.argv[3]+"_from_modes.pdb")#sys.argv[2] should be filename, sys.argv[3] should be motion_type ##sys.argv[4] should be mode type ##the structure annotation part ALPHA_annot = open("ALPHA_from_pdb.txt","r") BETA_annot = open("BETA_from_pdb.txt","r") def pickout(file1): name_list={} name='' first=0 for line in file1: line=line.rstrip() row=line.split("\t") if len(row[0]) ==4: name=row[0] first = int(row[1]) name_list[(name,first)]=[] old_line=[] if len(row[0])<4 and "#" not in line: start=int(row[0])-first end=int(row[1])-first name_list[(name,first)]=old_line old_line.append((start,end)) # print name_list return name_list def make_collection(sse_annot,filename): sse_dict = pickout(sse_annot) for key in sse_dict.keys(): if key[0] in filename: sse_list=[] chain1 = protein[0] for x in sse_dict[key]: #print filename #for r in (chain1[x[0]:x[1]+1]): sse_list.append(chain1[x[0]:x[1]+1]) return Collection(sse_list) ## helix is changed to sse_collection def calculatingthings(sse_collection,suffix, motion_type): cm, mi = sse_collection.centerAndMomentOfInertia() ev, axes = mi.diagonalization() # print argmin(ev), argmax(ev) sse_axis = axes[argmax(ev)].asVector().normal() vector = ParticleVector(universe) #print vector #print sse_axis # graphics = VMD # scene = graphics.Scene(scale=1./Units.Ang) # scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*vector[atom],0.01)) # # scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, # cm+2.*Units.nm*sse_axis, # 0.05)) # scene.writeToFile("test_scene.vrml") # scene.view() # vector=massWeightedNormalVector(vector) p = [] if motion_type == "indiv_translation": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) #print "7", len(sse_collection[7]) for counter in xrange(len(sse_collection)): element=sse_collection[counter] element_cm, element_mi = element.centerAndMomentOfInertia() ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmin(ele_ev)].asVector().normal() indiv_trans_vector = ParticleVector(universe) #print indiv_trans_vector for atom in element.atomList(): indiv_trans_vector[atom] = ele_axis scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*indiv_trans_vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) indiv_trans_vector=massWeightedNormalVector(indiv_trans_vector) for i in range(6, len(modes)): ## changed from 0 to 6 to ignore trivial mode = massWeightedNormalVector(modes[i]) # mode = modes.rawMode(i) projection = mode.massWeightedDotProduct(indiv_trans_vector) # projection = mode.dotProduct(indiv_trans_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_translation_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_translation_%s_%s.esc.dat'%(suffix,counter+1)) p = [] scene.writeToFile("indiv_translation_scene_%s.vmd"% suffix) if motion_type=="indiv_slide": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for counter in xrange(len(sse_collection)): element=sse_collection[counter] element_cm, element_mi = element.centerAndMomentOfInertia() # print element_mi ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmax(ele_ev)].asVector().normal() indiv_slide_vector = ParticleVector(universe) indiv_rot_vector = sse_axis.cross(cm - element_cm) # print cm - element_cm #indiv_rot_vector = massWeightedNormalVector(indiv_rot_vector) indiv_slide_v = sse_axis.cross(indiv_rot_vector) if (element_cm - cm) * indiv_slide_v <0: indiv_slide_v = indiv_slide_v * -1 #indiv_slide_v = massWeightedNormalVector(indiv_slide_v) for atom in element.atomList(): indiv_slide_vector[atom] = indiv_slide_v scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*indiv_slide_vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) indiv_slide_vector =massWeightedNormalVector(indiv_slide_vector) for i in range(6,len(modes)): ## changed from 0 to 6 to exclude first mode = massWeightedNormalVector(modes[i]) # mode = modes.rawMode(i) # projection = mode.dotProduct(indiv_slide_vector) projection = mode.massWeightedDotProduct(indiv_slide_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_slide_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_slide_%s_%s.esc.dat'% (suffix, counter+1)) p=[] scene.writeToFile("indiv_slide_scene_%s.vmd"% suffix) if motion_type=="indiv_tilt": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for counter in xrange(len(sse_collection)): #print counter element=sse_collection[counter] element_cm, element_mi = element.centerAndMomentOfInertia() ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmax(ele_ev)].asVector().normal() indiv_tilt_vector = ParticleVector(universe) indiv_rot_vector = sse_axis.cross(cm - element_cm) #print cm - element_cm indiv_tilt_v = sse_axis.cross(indiv_rot_vector) '''if (element_cm - cm) * indiv_tilt_v <0: indiv_slide_v = indiv_tilt_v * -1''' for atom in element.atomList(): tilt_disp = atom.position()[1] - element_cm[1] #tilt_disp = linalg.norm(atom.position()-element_cm) #if (atom.position()-element_cm) * indiv_tilt_v < 0: #tilt_disp = tilt_disp * -1 #print tilt_disp indiv_tilt_vector[atom] = indiv_tilt_v * tilt_disp scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*indiv_tilt_vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) indiv_tilt_vector =massWeightedNormalVector(indiv_tilt_vector) for i in range(6,len(modes)): ## changed from 0 to 6 to exclude first mode = massWeightedNormalVector(modes[i]) # mode = modes.rawMode(i) # projection = mode.dotProduct(indiv_tilt_vector) projection = mode.massWeightedDotProduct(indiv_tilt_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_tilt_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_tilt_%s_%s.esc.dat'% (suffix, counter+1)) p=[] scene.writeToFile("indiv_tilt_scene_%s.vmd"% suffix) if motion_type=="indiv_tilt_opposite": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for counter in xrange(len(sse_collection)): #print counter element=sse_collection[counter] element_cm, element_mi = element.centerAndMomentOfInertia() ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmax(ele_ev)].asVector().normal() indiv_tilt_opp_vector = ParticleVector(universe) indiv_rot_vector = sse_axis.cross(cm - element_cm) #print cm - element_cm indiv_tilt_opp_v = sse_axis.cross(indiv_rot_vector) '''if (element_cm - cm) * indiv_tilt_v <0: indiv_slide_v = indiv_tilt_v * -1''' for atom in element.atomList(): tilt_opp_disp = -1 * (atom.position()[1] - element_cm[1]) #tilt_disp = linalg.norm(atom.position()-element_cm) #if (atom.position()-element_cm) * indiv_tilt_v < 0: #tilt_disp = tilt_disp * -1 #print tilt_disp indiv_tilt_opp_vector[atom] = indiv_tilt_opp_v * tilt_opp_disp scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*indiv_tilt_opp_vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) indiv_tilt_opp_vector =massWeightedNormalVector(indiv_tilt_opp_vector) for i in range(6,len(modes)): ## changed from 0 to 6 to exclude first mode = massWeightedNormalVector(modes[i]) # mode = modes.rawMode(i) # projection = mode.dotProduct(indiv_tilt_vector) projection = mode.massWeightedDotProduct(indiv_tilt_opp_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_tilt_opp_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_tilt_opp_%s_%s.esc.dat'% (suffix, counter+1)) p=[] scene.writeToFile("indiv_tilt_opp_scene_%s.vmd"% suffix) if motion_type=="indiv_bend": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for counter in xrange(len(sse_collection)): #print counter element=sse_collection[counter] element_cm, element_mi = element.centerAndMomentOfInertia() ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmax(ele_ev)].asVector().normal() indiv_bend_vector = ParticleVector(universe) indiv_rot_vector = sse_axis.cross(cm - element_cm) #print cm - element_cm indiv_bend_v = sse_axis.cross(indiv_rot_vector) '''if (element_cm - cm) * indiv_bend_v <0: indiv_slide_v = indiv_bend_v * -1''' for atom in element.atomList(): bend_disp = (atom.position()[1] - element_cm[1]) if bend_disp < 0: bend_disp = bend_disp * 0.000001 #print bend_disp indiv_bend_vector[atom] = indiv_bend_v * bend_disp scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*indiv_bend_vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) indiv_bend_vector =massWeightedNormalVector(indiv_bend_vector) for i in range(6,len(modes)): ## changed from 0 to 6 to exclude first mode = massWeightedNormalVector(modes[i]) # mode = modes.rawMode(i) # projection = mode.dotProduct(indiv_bend_vector) projection = mode.massWeightedDotProduct(indiv_bend_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_bend_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_bend_%s_%s.esc.dat'% (suffix, counter+1)) p=[] scene.writeToFile("indiv_bend_scene_%s.vmd"% suffix) if motion_type=="indiv_bend_alt": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for counter in xrange(len(sse_collection)): #print counter element=sse_collection[counter] element_cm, element_mi = element.centerAndMomentOfInertia() ele_ev, ele_axes = element_mi.diagonalization() ele_axis = ele_axes[argmax(ele_ev)].asVector().normal() indiv_bend_alt_vector = ParticleVector(universe) indiv_rot_vector = sse_axis.cross(cm - element_cm) #print cm - element_cm indiv_bend_alt_v = sse_axis.cross(indiv_rot_vector) '''if (element_cm - cm) * indiv_bend_alt_v <0: indiv_slide_v = indiv_bend_alt_v * -1''' step = float(1)/len(element) weight=1 for atom in element.atomList(): bend_alt_disp = -1*(atom.position()[1] - element_cm[1]) if bend_alt_disp < 0: bend_alt_disp = bend_alt_disp * 0.000001 #print bend_alt_disp indiv_bend_alt_vector[atom] = indiv_bend_alt_v * bend_alt_disp * weight weight-=step scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*indiv_bend_alt_vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) indiv_bend_alt_vector =massWeightedNormalVector(indiv_bend_alt_vector) for i in range(6,len(modes)): ## changed from 0 to 6 to exclude first mode = massWeightedNormalVector(modes[i]) # mode = modes.rawMode(i) # projection = mode.dotProduct(indiv_bend_alt_vector) projection = mode.massWeightedDotProduct(indiv_bend_alt_vector) p.append((i+1,projection**2)) p=array(p) writeArray(p, 'indiv_bend_alt_%s_%s.pics.dat'% (suffix, counter+1)) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc, 'indiv_bend_alt_%s_%s.esc.dat'% (suffix, counter+1)) p=[] scene.writeToFile("indiv_bend_alt_scene_%s.vmd"% suffix) if motion_type=="rotation": # graphics = VMD # scene = graphics.Scene(scale=1./Units.Ang) for atom in sse_collection.atomList(): #print atom vector[atom] = sse_axis.cross(atom.position()-cm) # scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*vector[atom],0.01)) # scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, # cm+2.*Units.nm*sse_axis, # 0.05)) #for atom2 in modes.universe.atomList(): ## vector[atom2]=(cm-atom2.position()) #print vector vector=massWeightedNormalVector(vector) #for i in vector: # print i for i in range(6, len(modes)): ##changed from 0 to 6 to exclude trivial modes mode = massWeightedNormalVector(modes[i]) projection = mode.massWeightedDotProduct(vector) p.append((i+1,projection**2)) p = array(p) writeArray(p, 'rotation_%s.pics.dat' %suffix) pacc = add.accumulate(p) pacc[:,0] = p[:,0] writeArray(pacc,'rotation_%s.esc.dat' %suffix) # scene.writeToFile("rotation_nontriv_scene_%s.vmd"%suffix) if motion_type == "translation": graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) for atom in sse_collection.atomList(): vector[atom] = sse_axis scene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*vector[atom],0.01)) scene.addObject(graphics.Cylinder(cm-2.*Units.nm*sse_axis, cm+2.*Units.nm*sse_axis, 0.05)) # for atom2 in modes.universe.atomList(): # if atom2 not in sse_collection.atomList(): #vector[atom2]=(cm-atom2.position()) #print vector vector=massWeightedNormalVector(vector) #print vector p = [] pacc = 0 for i in range(6,len(modes)): ##changed from 0 mode=massWeightedNormalVector(modes[i]) # mode = modes.rawMode(i) projection=mode.massWeightedDotProduct(vector) # projection = mode.dotProduct(vector) p.append((i+1, projection**2)) p=array(p) pacc=add.accumulate(p) writeArray(p, 'translation_%s.pics.dat' %suffix) pacc[:,0] = p[:,0] writeArray(pacc, 'translation_%s.esc.dat' %suffix) # scene.writeToFile("translation_nontriv_raw_scene_%s.vmd"%suffix) return #print alpha_collection #print alpha_collection.centerAndMomentOfInertia() if "indiv" in sys.argv[3]: alpha_collection = make_collection(ALPHA_annot, filename) calculatingthings(alpha_collection, sys.argv[2]+"_ALPHA_"+sys.argv[4], sys.argv[3]) beta_collection = make_collection(BETA_annot, filename) calculatingthings(beta_collection, sys.argv[2]+"_BETA_"+sys.argv[4], sys.argv[3]) else: alpha_collection = make_collection(ALPHA_annot, filename) calculatingthings(alpha_collection, sys.argv[2]+"_ALPHA_"+sys.argv[4], sys.argv[3]) beta_collection = make_collection(BETA_annot, filename) calculatingthings(beta_collection, sys.argv[2]+"_BETA_"+sys.argv[4], sys.argv[3]) calculatingthings(protein, sys.argv[2]+"_WHOLE_"+sys.argv[4],sys.argv[3]) #print make_collection(ALPHA_annot,filename)