from MMTK import * from MMTK.PDB import PDBConfiguration from MMTK.Proteins import Protein from MMTK.NormalModes import NormalModes from Numeric import sqrt, add, array, argmax, argsort, argmin from Scientific import * from Scientific.IO.ArrayIO import writeArray from Scientific.Visualization import VMD import PlotTableR, util import os def massWeightedNormalVector(v): return v/sqrt(v.massWeightedDotProduct(v)) # # # # Vectors for kink-motion was added august 2007 by Edvin Fuglebakk. # # Code that does similar work on the kinkvectors and the previously implemented vectors might be scattered around the file # so be careful not to alter the meaning (or values) of: ## move_vector ## datafile ## path ## modes ## bundle ## names ## slide_pre_norm #... possibly more ... # from MotionMaker import * # # # Use this to see where the time is spent # ## def Rotation_wrapper(*args): ## import profile ## def rp(callableTest): ## prof = profile.Profile() ## prof.runcall(callableTest) ## import pstats ## stats = pstats.Stats(prof) ## stats.strip_dirs().sort_stats('cumulative').print_stats() ## def a(): ## return Rotation_wrapper0(*args) ## return rp(a) ## In this file we calculate the Overlap between the displacement vectors and the normal modes. def Rotation_wrapper(items, helices, bundle, residues, chains, names, kinkpoints): ## print "items" ## print items ## print "helices" ## print helices ## print "bundle" ## print bundle ## print "residues" ## print residues ## print "names" ## print names ## print "kinkpoints" ## print kinkpoints #print chains datafile = 'scratch/' + items['jobid'] graphics = VMD kink = graphics.Scene(scale=1./Units.Ang) tilt = graphics.Scene(scale=1./Units.Ang) rotation = graphics.Scene(scale=1./Units.Ang) translation = graphics.Scene(scale=1./Units.Ang) slide = graphics.Scene(scale=1./Units.Ang) bundle_rot = graphics.Scene(scale=1./Units.Ang) bundle_trans = graphics.Scene(scale=1./Units.Ang) ## Load the mode file modefile = items['modefile_name'] modes = load(modefile) universe = modes.universe # Then we have our protein protein = universe[0] natoms = len(protein.residues()) num_helices = len(helices) all_data = [] graphics = VMD scene = graphics.Scene(scale=1./Units.Ang) # There are three types of movement we wish to investigate. # We store all projections and cumulative values into these arrays. cum_rot_data = [] cum_trans_data = [] cum_move_data = [] cum_tilt_data = [] # cum_kink_data = [] # #Added aug 2007 #tables for plotting kink-vectors, similar to cum_rot_data, above # bananaPlot=[] midFlexPlot=[] upPlot=[] lowPlot=[] #Scenes for makeing vmd-files, similar to scenes above # bananaScene = graphics.Scene(scale=1./Units.Ang) midFlexScene = graphics.Scene(scale=1./Units.Ang) upScene = graphics.Scene(scale=1./Units.Ang) downScene = graphics.Scene(scale=1./Units.Ang) #/tables #/Scenes # # #Table for legendes for plot of kinkvectors # #This makes sure that every helix has got a name in kinkvectors. #For other vector unnamed helices ( name[helix]=='' ) are handled in PlotR kinknames=names[0:len(names)] #/legends rot_data = [] trans_data = [] move_data = [] tilt_data = [] # kink_data = [] helix_data = [] ## In this first part, we calculate the center of mass of the TM helical bundle ## We also calculate the bundle axis helix_bundle = Collection() all_axes=[] for i in range (0, len(bundle)): chain_nr = util.get_chain_nr(helices[bundle[i]][0], chains) chain = protein[chain_nr] start_res = helices[bundle[i]][1] start_res = util.findIndex(start_res, residues, chain_nr) h_len = len(helices[bundle[i]])-1 end_res = helices[bundle[i]][h_len] end_res = util.findIndex(end_res, residues, chain_nr) helix = chain[start_res:end_res] helix_bundle.addObject(helix) cm, mi = helix.centerAndMomentOfInertia() ev, axes = mi.diagonalization() helix_axis = axes[argmin(ev)].asVector().normal() all_axes.append(helix_axis) b_cm, b_mi = helix_bundle.centerAndMomentOfInertia() ev, axes = b_mi.diagonalization() s_min = 0 best_index = 0 for i in range(0, len(axes)): for j in range (0, len(all_axes)): s = all_axes[j]*(axes[i].asVector().normal()) s = sqrt(s*s) if s > s_min: s_min = s best_index = i #bundle_axis = axes[argmin(ev)].asVector().normal() bundle_axis = axes[best_index].asVector().normal() ## This loop goes once for each helix in the TM region. ## Calculating the overlap for each one. for i in range (0, len(bundle)): print 'calculating helix ' + str(bundle[i]) # The three displacement vectors are first set according to the initial configuration rot_vector = ParticleVector(universe) trans_vector = ParticleVector(universe) move_vector = ParticleVector(universe) tilt_vector = ParticleVector(universe) # kink_vector = ParticleVector(universe) chain_nr = util.get_chain_nr(helices[bundle[i]][0], chains) chain = protein[chain_nr] start_res = helices[bundle[i]][1] start_res = util.findIndex(start_res, residues, chain_nr) h_len = len(helices[bundle[i]])-1 end_res = helices[bundle[i]][h_len] end_res = util.findIndex(end_res, residues, chain_nr) helix = chain[start_res:end_res] helix_bundle.addObject(helix) # Calculating the axis of the helix cm, mi = helix.centerAndMomentOfInertia() ev, axes = mi.diagonalization() helix_axis = axes[argmin(ev)].asVector().normal() #slide.addObject(graphics.Cylinder(cm-2.*Units.nm*helix_axis,cm+2.*Units.nm*helix_axis,0.05)) # The slide vector is a vector moving away from the center of mass of the bundle (b_cm) # we normalize it first... #slide_v = cm - b_cm #slide_norm = sqrt( slide_v.x()*slide_v.x() + slide_v.y()*slide_v.y() + slide_v.z()*slide_v.z() ) #slide_v = slide_v / slide_norm rot = bundle_axis.cross(b_cm-cm) rot_norm = sqrt( rot.x()*rot.x() + rot.y()*rot.y() + rot.z()*rot.z() ) rot = rot / rot_norm slide_v = helix_axis.cross( rot ) if (cm-b_cm) * slide_v < 0: slide_v = slide_v * -1 slide_norm = sqrt( slide_v.x()*slide_v.x() + slide_v.y()*slide_v.y() + slide_v.z()*slide_v.z() ) slide_v = slide_v / slide_norm tilt_v = rot.cross(helix_axis) tilt_norm = sqrt( tilt_v.x()*tilt_v.x() + tilt_v.y()*tilt_v.y() + tilt_v.z()*tilt_v.z() ) tilt_v = tilt_v / tilt_norm # kink_v = tilt_v # Loop going for each atom in the helix # This is where we set the displacement vector of each atom in the helix length = float(int(len(helix))) if int(length)%2 == 0.0: counter = float(int(-1*(length/2))) else: counter = float(int(-1*((length-1)/2))) for atom in helix.atomList(): #rot_vector[atom] = helix_axis.cross(atom.position()-cm) t = helix_axis.cross(atom.position()-cm) rot_norm = sqrt( t.x()*t.x() + t.y()*t.y() + t.z()*t.z() ) rot_vector[atom] = t / rot_norm trans_vector[atom] = helix_axis move_vector[atom] = slide_v #t = rot_vector[atom] #print sqrt( t.x()*t.x() + t.y()*t.y() + t.z()*t.z() ) if counter == 0.0: if int(length)%2 == 0.0: counter = float(1.0)+float(counter) #kink_vector[atom] = kink_v / ( (length/2) / sqrt(counter*counter) ) tilt_vector[atom] = ( tilt_v / ((length/2) / counter )) else: #kink_vector[atom] = Vector(0.00001,0.00001,0.00001) tilt_vector[atom] = Vector(0.00001,0.00001,0.00001) else: #kink_vector[atom] = kink_v / ( (length/2) / sqrt(counter*counter) ) tilt_vector[atom] = ( tilt_v / ((length/2) / counter )) #print sqrt( tilt_vector[atom].x()*tilt_vector[atom].x() + tilt_vector[atom].y()*tilt_vector[atom].y() + tilt_vector[atom].z()*tilt_vector[atom].z() ) counter+=1.0 #kink.addObject(graphics.Arrow(atom.position(),atom.position()+1.*kink_vector[atom],0.01)) tilt.addObject(graphics.Arrow(atom.position(),atom.position()+1.*tilt_vector[atom],0.01)) rotation.addObject(graphics.Arrow(atom.position(),atom.position()+1.*rot_vector[atom],0.01)) slide.addObject(graphics.Arrow(atom.position(), atom.position()+1.*move_vector[atom], 0.01)) translation.addObject(graphics.Arrow(atom.position(), atom.position()+1.*trans_vector[atom], 0.01)) # #Aug07: #Keep a reference to the unormalized move_vector as we need it to make the kink-vectors # # slide_pre_norm=move_vector # Normalizing the displacement vectors rot_vector=massWeightedNormalVector(rot_vector) trans_vector=massWeightedNormalVector(trans_vector) move_vector=massWeightedNormalVector(move_vector) tilt_vector=massWeightedNormalVector(tilt_vector) # kink_vector=massWeightedNormalVector(kink_vector) p_rot = [] p_trans = [] p_move = [] p_tilt = [] # p_kink = [] rot_pacc = 0 trans_pacc = 0 move_pacc = 0 tilt_pacc = 0 # kink_pacc = 0 # Let us project these displacement vectors on to the full set of normal modes for j in range(6, len(modes)): # since we want to start counting from 1, which is the first non-trivial mode, we define an output mode-number output_mode_nr = j-5 mode = massWeightedNormalVector(modes[j]) projection1 = mode.massWeightedDotProduct(rot_vector) projection2 = mode.massWeightedDotProduct(trans_vector) projection3 = mode.massWeightedDotProduct(move_vector) projection4 = mode.massWeightedDotProduct(tilt_vector) # projection5 = mode.massWeightedDotProduct(kink_vector) p_rot.append((output_mode_nr, (projection1**2))) p_trans.append((output_mode_nr, (projection2**2))) p_move.append((output_mode_nr, (projection3**2))) p_tilt.append((output_mode_nr, (projection4**2))) # p_kink.append((output_mode_nr, (projection5**2))) p_rot = array(p_rot) p_trans = array(p_trans) p_move = array(p_move) p_tilt = array(p_tilt) # p_kink = array(p_kink) # The cumulative values is simply summing the projections rot_pacc = add.accumulate(p_rot) trans_pacc = add.accumulate(p_trans) move_pacc = add.accumulate(p_move) tilt_pacc = add.accumulate(p_tilt) # kink_pacc = add.accumulate(p_kink) rot_pacc[:,0] = p_rot[:,0] trans_pacc[:,0] = p_trans[:,0] move_pacc[:,0] = p_move[:,0] tilt_pacc[:,0] = p_tilt[:,0] # kink_pacc[:,0] = p_kink[:,0] # Raw data is stored in the files writeArray(rot_pacc, datafile+'.ROT.CUM.' + str(bundle[i]) + '.dat') writeArray(trans_pacc, datafile+'.TRANS.CUM.' + str(bundle[i]) + '.dat') writeArray(move_pacc, datafile+'.MOVE.CUM.' + str(bundle[i]) + '.dat') writeArray(tilt_pacc, datafile+'.TILT.CUM.' + str(bundle[i]) + '.dat') # writeArray(kink_pacc, datafile+'.KINK.CUM.' + str(bundle[i]) + '.dat') writeArray(p_rot, datafile+'.ROT.' + str(bundle[i]) + '.dat') writeArray(p_trans, datafile+'.TRANS.' + str(bundle[i]) + '.dat') writeArray(p_move, datafile+'.MOVE.' + str(bundle[i]) + '.dat') writeArray(p_tilt, datafile+'.TILT.' + str(bundle[i]) + '.dat') # writeArray(p_kink, datafile+'.KINK.' + str(bundle[i]) + '.dat') d = cm, str(start_res)+':'+str(end_res), rot_pacc, trans_pacc, bundle[i] helix_data.append(d) # We need to store the overlap data because we are not ready for plotting just yet. cum_rot_data.append(rot_pacc) cum_trans_data.append(trans_pacc) cum_move_data.append(move_pacc) cum_tilt_data.append(tilt_pacc) # cum_kink_data.append(kink_pacc) rot_data.append(p_rot) trans_data.append(p_trans) move_data.append(p_move) tilt_data.append(p_tilt) # kink_data.append(p_kink) # # # Aug 2007 # Added kinkvectors # # # # #slice off trivial modes nontrivials=modes[6:len(modes)] #partition helix #find kinkpoint #calculate the index in helix from the given pdb identifiers (residuenumbers) pdbIndex=kinkpoints[bundle[i]] helixIndex=None for j in range(0, len(helices[bundle[i]])): if helices[bundle[i]][j] == pdbIndex: helixIndex=j-1 #the tables in helices has chain-letter prepended to them if helixIndex==None: print "pdbIndex: " + str(pdbIndex) print helices[bundle[i]] raise Exception("kinkpoints pdbindex not in this helix") lower=helix[0:helixIndex] kinkpoint=helix[helixIndex:helixIndex+1] upper=helix[helixIndex+1:len(helix)] #make the kinkvectors #the slide vector is called move_vector in the above code banana=makeBananaKink(lower, kinkpoint, upper, slide_pre_norm) midFlex=makeMidFlexKink(lower, kinkpoint, upper, slide_pre_norm) upKink=makeUpKink(upper, helix, slide_pre_norm) lowKink=makeLowKink(lower, helix, slide_pre_norm) #find the best rotation bestBanana, bananaAngle, bananaCumOverlap = bestRotation(banana, 0.6, nontrivials, helix) bestMidFlex, midFlexAngle, midFlexCumOverlap = bestRotation(midFlex, 0.6, nontrivials, helix) bestUpKink, upAngle, upCumOverlap = bestRotation(upKink, 0.6, nontrivials, helix) bestLowKink, lowAngle, lowCumOverlap = bestRotation(lowKink, 0.6, nontrivials, helix) # #make all files and prepare data for the plots # #the cumulative overlaps writeArray(array(bananaCumOverlap), datafile+'.BANANA_KINK.CUM.'+str(bundle[i])+'.dat') writeArray(array(midFlexCumOverlap), datafile+'.MIDFLEX_KINK.CUM.'+str(bundle[i])+'.dat') writeArray(array(upCumOverlap), datafile+'.UP_KINK.CUM.'+str(bundle[i])+'.dat') writeArray(array(lowCumOverlap), datafile+'.DOWN_KINK.CUM.'+str(bundle[i])+'.dat') #the non-cumulative overlaps ... recalculates for now #--- recalculation can probably be avoided --- writeArray(array( overlapProfile(bestBanana, nontrivials) ), datafile+'.BANANA_KINK.'+str(bundle[i])+'.dat') writeArray(array( overlapProfile(bestMidFlex, nontrivials) ), datafile+'.MIDFLEX_KINK.'+str(bundle[i])+'.dat') writeArray(array( overlapProfile(bestUpKink, nontrivials) ), datafile+'.UP_KINK.'+str(bundle[i])+'.dat') writeArray(array( overlapProfile(bestLowKink, nontrivials) ), datafile+'.DOWN_KINK.'+str(bundle[i])+'.dat') #add data for this helix to the plots bananaPlot.append(bananaCumOverlap) midFlexPlot.append(midFlexCumOverlap) upPlot.append(upCumOverlap) lowPlot.append(lowCumOverlap) #add legend to kinkvectors # #This should be consistent with how PlotR deals with '' names # if kinknames[bundle[i]] == '': kinknames[bundle[i]] = str(bundle[i]) kinknames[bundle[i]]+= " k:" + str(pdbIndex) + " r:" + str(bananaAngle) #add vectors to scenes for vmd-files for atom in helix.atomList(): bananaScene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*bestBanana[atom],0.01)) midFlexScene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*bestMidFlex[atom],0.01)) upScene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*bestUpKink[atom],0.01)) downScene.addObject(graphics.Arrow(atom.position(),atom.position()+1.*bestLowKink[atom],0.01)) # # # #/added kinkvectors # # # print 'done for all helices! yayey!!!' print 'start plotting' #Having calculated the overlap for all helices, we can go ahead with the plotting #suffix = 'png' #suffix = 'pdf' rot_file = datafile+'.ROT.' trans_file = datafile+'.TRANS.' move_file = datafile+'.MOVE.' tilt_file = datafile+'.TILT.' # kink_filee = datafile+'.KINK.' cum_rot_file = datafile+'.CUM.ROT.' cum_trans_file = datafile+'.CUM.TRANS.' cum_move_file = datafile+'.CUM.MOVE.' cum_tilt_file = datafile+'.CUM.TILT.' # cum_kink_file = datafile+'.CUM.KINK.' # R handles the plotting for us. The R code is saved in PlotTableR PlotTableR.Plot( cum_rot_data, 'Cumulative Rotation of helices', 'mode number', 'cumulative squared overlap', cum_rot_file, modes, bundle, 'pdf', names, ) PlotTableR.Plot( cum_trans_data, 'Cumulative Translation of helices', 'mode number', 'cumulative squared overlap', cum_trans_file, modes, bundle, 'pdf', names,) PlotTableR.Plot( cum_move_data, 'Cumulative Slide of helices', 'mode number', 'cumulative squared overlap', cum_move_file, modes, bundle, 'pdf', names,) PlotTableR.Plot( cum_tilt_data, 'Cumulative Tilt of helices', 'mode number', 'cumulative squared overlap', cum_tilt_file, modes, bundle, 'pdf', names ) # PlotTableR.Plot( cum_kink_data, 'Cumulative Kink of helices', 'mode number', 'cumulative squared overlap', cum_kink_file, modes, bundle, 'pdf' ) #PlotTableR.Plot( cum_rot_data, 'Cumulative Rotation of helices', 'mode number', 'cumulative squared overlap', cum_rot_file, modes, bundle, 'png' ) #PlotTableR.Plot( cum_trans_data, 'Cumulative Translation of helices', 'mode number', 'cumulative squared overlap', cum_trans_file, modes, bundle, 'png' ) #PlotTableR.Plot( cum_move_data, 'Cumulative Slide of helices', 'mode number', 'cumulative squared overlap', cum_move_file, modes, bundle, 'png' ) #PlotTableR.Plot( cum_tilt_data, 'Cumulative Tilt of helices', 'mode number', 'cumulative squared overlap', cum_tilt_file, modes, bundle, 'png' ) # PlotTableR.Plot( cum_kink_data, 'Cumulative Kink of helices', 'mode number', 'cumulative squared overlap', cum_kink_file, modes, bundle, 'png' ) print 'done plotting' #### HERE COMES THE CALCULATION OF THE MOBILITY OF THE BUNDLE ############ print 'starting bundle computing' whole_vector_rot = ParticleVector(universe) # whole_vector_trans = ParticleVector(universe) helix_bundle = Collection() #all_axes = [] for i in range (0, len(bundle)): rot_vector = ParticleVector(universe) # trans_vector = ParticleVector(universe) #chain_nr = helices[bundle[i]][0] chain_nr = util.get_chain_nr(helices[bundle[i]][0], chains) chain = protein[chain_nr] #start_res = helices[bundle[i]][1]-1 start_res = helices[bundle[i]][1] start_res = util.findIndex(start_res, residues, chain_nr) h_len = len(helices[bundle[i]])-1 #end_res = helices[bundle[i]][h_len] end_res = helices[bundle[i]][h_len] end_res = util.findIndex(end_res, residues, chain_nr) helix = chain[start_res:end_res] helix_bundle.addObject(helix) tilt.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) rotation.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) translation.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) slide.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) bundle_rot.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) # bundle_trans.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) # #Added aug2007 # bananaScene.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) midFlexScene.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) upScene.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) downScene.addObject(graphics.Cylinder(b_cm-2.*Units.nm*bundle_axis,b_cm+2.*Units.nm*bundle_axis,0.05)) #/ # for atom in helix_bundle.atomList(): # # Aug07: # Endret parameter bm til b_cm # # rot = bundle_axis.cross(atom.position()-b_cm) rot_norm = sqrt( rot.x()*rot.x() + rot.y()*rot.y() + rot.z()*rot.z() ) rot = rot / rot_norm whole_vector_rot[atom] = rot # whole_vector_trans[atom] = bundle_axis bundle_rot.addObject(graphics.Arrow(atom.position(),atom.position()+1.*whole_vector_rot[atom],0.01)) # bundle_trans.addObject(graphics.Arrow(atom.position(),atom.position()+1.*whole_vector_trans[atom],0.01)) p_rot = [] # p_trans = [] whole_vector_rot = massWeightedNormalVector(whole_vector_rot) # whole_vector_trans= massWeightedNormalVector(whole_vector_trans) for j in range(6, len(modes)): # since we want to start counting from 1, which is the first non-trivial mode, we define an output mode-number output_mode_nr = j-5 mode = massWeightedNormalVector(modes[j]) projection1 = mode.massWeightedDotProduct(whole_vector_rot) # projection2 = mode.massWeightedDotProduct(whole_vector_trans) p_rot.append((output_mode_nr, (projection1**2))) # p_trans.append((j+1, (projection2**2))) p_rot = array(p_rot) rot_pacc2 = add.accumulate(p_rot) # trans_pacc2 = add.accumulate(p_trans) rot_pacc2[:,0] = p_rot[:,0] # trans_pacc2[:,0] = p_trans[:,0] writeArray(p_rot, datafile+'.ROT.BUNDLE.dat') # writeArray(p_trans, datafile+'.TRANS.BUNDLE.dat') rot_file = datafile+'.ROT.BUNDLE.' # trans_file = datafile+'.TRANS.BUNDLE.' cum_rot_file = datafile+'.CUM.ROT.BUNDLE.' # cum_trans_file = datafile+'.CUM.TRANS.BUNDLE.' writeArray(rot_pacc2, datafile+'.CUM.ROT.BUNDLE.dat') # writeArray(trans_pacc2, datafile+'.CUM.TRANS.BUNDLE.dat') PlotTableR.Plot2( p_rot, 'Rotation of the bundle', 'mode number', 'squared overlap', rot_file, modes, 'pdf' ) # PlotTableR.Plot2( p_rot, 'Rotation of the bundle', 'mode number', 'squared overlap', rot_file, modes, 'png' ) # PlotTableR.Plot2( p_trans, 'Translation of the bundle', 'mode number', 'squared overlap', trans_file, modes, 'pdf' ) # PlotTableR.Plot2( p_trans, 'Translation of the bundle', 'mode number', 'squared overlap', trans_file, modes, 'png' ) PlotTableR.Plot3( rot_pacc2, 'Cumulative rotation of the bundle', 'mode number', 'cumulative squared overlap', cum_rot_file, modes, 'pdf' ) #PlotTableR.Plot3( rot_pacc2, 'Cumulative rotation of the bundle', 'mode number', 'cumulative squared overlap', cum_rot_file, modes, 'png' ) # PlotTableR.Plot3( trans_pacc2, 'Cumulative translation of the bundle', 'mode number', 'cumulative squared overlap', cum_trans_file, modes, 'pdf' ) # PlotTableR.Plot3( trans_pacc2, 'Cumulative translation of the bundle', 'mode number', 'cumulative squared overlap', cum_trans_file, modes, 'png' ) print 'done bundle computing' path = './scratch/'+items['jobid']+'/' os.system('mkdir ' + path) os.system('mkdir ' + path + 'HelixRotation') os.system('mkdir ' + path + 'HelixTranslation') os.system('mkdir ' + path + 'HelixSlide') os.system('mkdir ' + path + 'HelixTilt') ## os.system('mkdir ' + path + 'HelixKink') os.system('mkdir ' + path + 'CumulativeHelixRotation') os.system('mkdir ' + path + 'CumulativeHelixTranslation') os.system('mkdir ' + path + 'CumulativeHelixSlide') os.system('mkdir ' + path + 'CumulativeHelixTilt') ## os.system('mkdir ' + path + 'CumulativeHelixKink') os.system('mkdir ' + path + 'BundleRotation') # os.system('mkdir ' + path + 'BundleTranslation') os.system('mkdir ' + path + 'CumulativeBundleRotation') # os.system('mkdir ' + path + 'CumulativeBundleTranslation') os.system('cd ./scratch && cp '+ items['jobid']+'.ROT.?.dat ' + items['jobid']+'/' + 'HelixRotation/') os.system('cd ./scratch && cp '+ items['jobid']+'.TRANS.?.dat ' +items['jobid']+'/' + 'HelixTranslation/') os.system('cd ./scratch && cp '+ items['jobid']+'.MOVE.?.dat ' + items['jobid']+'/' + 'HelixSlide/') os.system('cd ./scratch && cp '+ items['jobid']+'.TILT.?.dat ' + items['jobid']+'/' + 'HelixTilt/') ## os.system('cd ./scratch && cp '+ items['jobid']+'.KINK.?.dat ' + items['jobid']+'/' + 'HelixKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.ROT.??.dat ' + items['jobid']+'/' + 'HelixRotation/') os.system('cd ./scratch && cp '+ items['jobid']+'.TRANS.??.dat ' +items['jobid']+'/' + 'HelixTranslation/') os.system('cd ./scratch && cp '+ items['jobid']+'.MOVE.??.dat ' + items['jobid']+'/' + 'HelixSlide/') os.system('cd ./scratch && cp '+ items['jobid']+'.TILT.??.dat ' + items['jobid']+'/' + 'HelixTilt/') ## os.system('cd ./scratch && cp '+ items['jobid']+'.KINK.??.dat ' + items['jobid']+'/' + 'HelixKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.ROT.CUM.*.dat ' +items['jobid']+'/' + 'CumulativeHelixRotation/') os.system('cd ./scratch && cp '+ items['jobid']+'.TRANS.CUM.*.dat ' + items['jobid']+'/' + 'CumulativeHelixTranslation/') os.system('cd ./scratch && cp '+ items['jobid']+'.MOVE.CUM.*.dat ' + items['jobid']+'/' + 'CumulativeHelixSlide/') os.system('cd ./scratch && cp '+ items['jobid']+'.TILT.CUM.*.dat ' + items['jobid']+'/' + 'CumulativeHelixTilt/') ## os.system('cd ./scratch && cp '+ items['jobid']+'.KINK.CUM.*.dat ' + items['jobid']+'/' + 'CumulativeHelixKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.ROT.BUNDLE.dat ' + items['jobid']+'/' + 'BundleRotation/') # os.system('cd ./scratch && cp '+ items['jobid']+'.TRANS.BUNDLE.dat ' + items['jobid']+'/' + 'BundleTranslation/') os.system('cd ./scratch && cp '+ items['jobid']+'.CUM.ROT.BUNDLE.dat ' +items['jobid']+'/' + 'CumulativeBundleRotation/') # os.system('cd ./scratch && cp '+ items['jobid']+'.CUM.TRANS.BUNDLE.dat ' + items['jobid']+'/' + 'CumulativeBundleTranslation/') # # Aug 2007: This was moved to the end to include kink-vector data as well # # os.system('cd ./scratch && zip -r '+items['jobid']+'.RAW_DATA.zip ' + items['jobid']+'/') # os.system('cd ./scratch && rm -r '+ items['jobid']) # # #kink.addObject(graphics.Molecules(items['pdbfilename'])) #tilt.addObject(graphics.Molecules(items['pdbfilename'])) #rotation.addObject(graphics.Molecules(items['pdbfilename'])) #translation.addObject(graphics.Molecules(items['pdbfilename'])) #slide.addObject(graphics.Molecules(items['pdbfilename'])) #bundle_rot.addObject(graphics.Molecules(items['pdbfilename'])) #bundle_trans.addObject(graphics.Molecules(items['pdbfilename'])) ## kink.writeToFile('./scratch/'+items['jobid']+'.KINK.vmd') tilt.writeToFile('./scratch/'+items['jobid']+'.TILT.vmd') rotation.writeToFile('./scratch/'+items['jobid']+'.ROTATION.vmd') translation.writeToFile('./scratch/'+items['jobid']+'.TRANSLATION.vmd') slide.writeToFile('./scratch/'+items['jobid']+'.SLIDE.vmd') bundle_rot.writeToFile('./scratch/'+items['jobid']+'.BUNDLE_ROT.vmd') # bundle_trans.writeToFile('./scratch/'+items['jobid']+'.BUNDLE_TRANS.vmd') #kink.view() #slide.view() ## ## ## Aug 2007 Added plotting and storage of kinkvectors ## print 'plot kink-vectors' #plot kinkvectors PlotTableR.Plot( bananaPlot , 'Cumulative kink of helices (kink 1)', 'mode number', 'cumulative squared overlap', datafile+'.CUM.BANANA_KINK.', modes, bundle, 'pdf', kinknames, 3.0 ) PlotTableR.Plot( midFlexPlot , 'Cumulative kink of helices (kink 2)', 'mode number', 'cumulative squared overlap', datafile+'.CUM.MIDFLEX_KINK.', modes, bundle, 'pdf', kinknames, 3.0 ) PlotTableR.Plot( upPlot , 'Cumulative kink of helices (kink 3)', 'mode number', 'cumulative squared overlap', datafile+'.CUM.UP_KINK.', modes, bundle, 'pdf', kinknames, 3.0 ) PlotTableR.Plot( lowPlot , 'Cumulative kink of helices (kink 4)', 'mode number', 'cumulative squared overlap', datafile+'.CUM.DOWN_KINK.', modes, bundle, 'pdf', kinknames, 3.0 ) print 'really done plotting' #add stuff to the scratch/ directory #basicly copied from above code, I don't know the reason for every line os.system('mkdir ' + path + 'HelixBananaKink') os.system('mkdir ' + path + 'HelixMidFlexKink') os.system('mkdir ' + path + 'HelixUpKink') os.system('mkdir ' + path + 'HelixDownKink') os.system('mkdir ' + path + 'CumulativeHelixBananaKink') os.system('mkdir ' + path + 'CumulativeHelixMidFlexKink') os.system('mkdir ' + path + 'CumulativeHelixUpKink') os.system('mkdir ' + path + 'CumulativeHelixDownKink') os.system('cd ./scratch && cp '+ items['jobid']+'.BANANA_KINK.?.dat ' + items['jobid']+'/' + 'HelixBananaKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.MIDFLEX_KINK.?.dat ' + items['jobid']+'/' + 'HelixMidFlexKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.UP_KINK.?.dat ' + items['jobid']+'/' + 'HelixUpKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.DOWN_KINK.?.dat ' + items['jobid']+'/' + 'HelixDownKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.BANANA_KINK.??.dat ' + items['jobid']+'/' + 'HelixBananaKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.MIDFLEX_KINK.??.dat ' + items['jobid']+'/' + 'HelixMidFlexKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.UP_KINK.??.dat ' + items['jobid']+'/' + 'HelixUpKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.DOWN_KINK.??.dat ' + items['jobid']+'/' + 'HelixDownKink/') os.system('cd ./scratch && cp '+ items['jobid']+'.BANANA_KINK.CUM.*.dat ' +items['jobid']+'/' + 'CumulativeHelixBananaKink') os.system('cd ./scratch && cp '+ items['jobid']+'.MIDFLEX_KINK.CUM.*.dat ' + items['jobid']+'/' + 'CumulativeHelixMidFlexKink') os.system('cd ./scratch && cp '+ items['jobid']+'.UP_KINK.CUM.*.dat ' + items['jobid']+'/' + 'CumulativeHelixUpKink') os.system('cd ./scratch && cp '+ items['jobid']+'.DOWN_KINK.CUM.*.dat ' + items['jobid']+'/' + 'CumulativeHelixDownKink') #this was moved down here from above code to include kink-vector data in zip-file os.system('cd ./scratch && zip -r '+items['jobid']+'.RAW_DATA.zip ' + items['jobid']+'/') os.system('cd ./scratch && rm -r '+ items['jobid']) #write vmd-files bananaScene.writeToFile('./scratch/'+items['jobid']+'.BANANA_KINK.vmd') midFlexScene.writeToFile('./scratch/'+items['jobid']+'.MIDFLEX_KINK.vmd') upScene.writeToFile('./scratch/'+items['jobid']+'.UP_KINK.vmd') downScene.writeToFile('./scratch/'+items['jobid']+'.DOWN_KINK.vmd') # # #/plotting and storage of kinkvectors # # return 1 def sort_by_value(d): """ Returns the keys of dictionary d sorted by their values """ items=d.items() backitems=[ [v[1],v[0]] for v in items] backitems.sort() backitems.reverse() return [ int(backitems[i][1]) for i in range(0,len(backitems)) ] # #Aug 2007: #some aux-functions # # def bestRotation(vector, goal, modes, helix, degrees=[0,30,60,90,120,150]): #0 could be taken away and the unrotated vector prependend to the list of rotations """ Determines which of the rotations of vector that first reaches goal in a cumulative squared projections plot over the normalmodes mode returns (vector, rotation in degreees, cumulative projections-plot """ #parameters neede for rotations axis, __, __ = helixParams(helix) radians=[] for d in degrees: radians.append(math.radians(d)) rotations=makeRotations(vector, axis, radians) #project all the rotations and sort them by modes needed to reach goal projections=[] for i in range(0, len(rotations)): projection = cumulativeOverlapProfile(rotations[i], modes) goalIndex=lowestAbove(projection, goal) projections.append( (goalIndex, rotations[i], degrees[i], projection) ) projections.sort() return projections[0][1:4] def lowestAbove(table, delim): """ returns the first value higher than delim in the sortet array table """ for i in range (0, len(table)): if table[i] > delim: return i if __name__ == '__main__': initiate()