Scripts:Overlap Script

From Harlem

#  Script to superimpose two molecules 
#  using a Monte Carlo trajectory.
#  Molecule set can consist of several molecules.
#  Harlem file should exist in the same directory. 
#  Make pdb file from harlem file.
#  Load molecule in HLM format.
#  Load script and execute it. 
#
#  Carnegie Mellon University
#  Department of Chemistry
#  Kurnikova's lab

import os
nm1 = []

# name of the trajectory file, format (snaphoot,x,y,z,q1,q2,q3,q4)
trj_file = "all_mc_pts_May13_100K.dat"
#name of file with rmsd values?
out_file1=open("rmsd_1py6_SurfEne.out","w")
# what HLM file to use for fitting?
nm1.append("Crystal_Jose.hlm") # HLM file 

# load HLM file on pmset1 object
ifile = open(trj_file, "r")
pmset1 = HaMolSet()
pmset1.FetchFile(FormatHarlem,nm1[0])

# change the identification number of residues
n=0
aitr_m1 = AtomIteratorMolSet(pmset1)
aptr = aitr_m1.GetFirstAtom()
while( aptr != None):
    at_name = aptr.GetName()
    if(at_name =="N"):	
	    n= n+1
    pres = aptr.GetHostRes()
    pres.serno = n    
    aptr = aitr_m1.GetNextAtom()

#save the structure in PDB format
pmset1.SavePDBFile("temp.pdb")

#creates two HaMolset objects from the PDB file
pmset = HaMolSet()
pmset2 =HaMolSet()
pmset.FetchFile(FormatPDB,"temp.pdb")
pmset2.FetchFile(FormatPDB,"temp.pdb")

#gets the number of molecules in pmset1
nmol = pmset1.GetNMol()

#initializes the Empirical Potential method 
empmod = pmset1.GetEmpiricalMod(1)
empmod.Initialize()

#creates auxiliary arrays
atl1 = AtomGroup()
atm_ca_arr_set1 = []
atm_ca_arr_set2 = []
new_at = HaAtom()

#creates the first atom group with the template structure atl1
#the template atom group in this case will contains CA atoms
#the identity of atoms can be modified ie: ("CA" and "N")
aitr_m1 = AtomIteratorMolSet(pmset)
aptr = aitr_m1.GetFirstAtom()
while( aptr != None):
    at_name = aptr.GetName()
    if(at_name =="CA") :
        atl1.InsertAtom(aptr)
    aptr = aitr_m1.GetNextAtom()

#creates and array of CA atoms from pmset2 object(atm_ca_arr_set2)
aitr_m1 = AtomIteratorMolSet(pmset2)
aptr = aitr_m1.GetFirstAtom()
while( aptr != None):
	at_name = aptr.GetName()
 	if(at_name =="CA") :
		atm_ca_arr_set2.append( aptr)
	aptr = aitr_m1.GetNextAtom()

#creates an array of CA atoms from pmset1 object(atm_ca_arr_set1)
aitr_m1 = AtomIteratorMolSet(pmset1)
aptr = aitr_m1.GetFirstAtom()
while( aptr != None):
	at_name = aptr.GetName()
	if(at_name =="CA") :
		atm_ca_arr_set1.append( aptr)
	aptr = aitr_m1.GetNextAtom()

#creates the translation and quaternion vectors
#creates qang,qx,qy,qz components
trans = Vec3D()
quat = Quaternion()
qang = doublep()
qx = doublep()
qy = doublep()
qz = doublep()

#start to read the trajectory file
snap = 0 #count the number of snapshoots
str = ifile.readline() #read each line of the trajectory file 

while (1):
         snap = snap +1
         coord_arr=[]
         imol=0
         # Set positions of rigid bodies in pmset1 objects
         while ( imol < nmol):   
            nums = str.split() #read the coordinate values for each rigid body
            qang.assign(float(nums[4]))
            qx.assign(float(nums[5]))
            qy.assign(float(nums[6]))
            qz.assign(float(nums[7]))
            trans.SetX( float(nums[1]) )
            trans.SetY( float(nums[2]))
            trans.SetZ( float(nums[3]))
            #get the molecule imol from pmset1 object
            pmol1 = pmset1.GetMoleculeNum(imol)
            quat.SetQuaternion(qang,qx,qy,qz)
            #translate and rotate the molecule imol
            pmol1.SetQuaternionTrans(quat, trans)
            str = ifile.readline()
            imol = imol+1 #while loop
            
         if(snap > 0 and snap%1 == 0):

          #creates atl2 Atomgroup that will contains the CA atoms of the current snapshoot
          atl2 = AtomGroup()

          #gets the x,y,z coordinates for each CA atom in pmset1 (already translated and rotated) 
          for i in range (len(atm_ca_arr_set1) ):
                new_at = atm_ca_arr_set1[i]
                coord_arr.append(new_at.GetX() )
                coord_arr.append(new_at.GetY() )
                coord_arr.append(new_at.GetZ() )

          #update the x,y,z coordinates for each CA atom in pmset2 (pmset2 now translated and rotated)
          #add all CA atoms from pmset2 to atl2 group
          ind =0
          for i in range (len(atm_ca_arr_set2) ):
                new_at = atm_ca_arr_set2[i]
                new_at.SetX(coord_arr[ind])
                ind = ind+1
                new_at.SetY(coord_arr[ind])
                ind = ind+1
                new_at.SetZ(coord_arr[ind])
                ind = ind+1
                atl2.InsertAtom(new_at)

          #calculates the RMSD and energy value
          eps = pmset.OverlapMol(atl1, atl2)
          enep = empmod.PenaltyPairwise()

          #saves a temporal snapshoot
          pmset1.SavePDBFile("temp.pdb")

          #print the snaphoot number, RMSD and energy
          print >> out_file1, snap, eps, enep
         if str == '':
                break

out_file1.close()
ifile.close()
print "###END###"