How to read a pdb file and perform FFT based docking using python?
Question: Suppose I have two PDB files (one of them is as follows)

    
    
    ATOM      1  N   MET A   1      66.104  56.583 -35.505  
    ATOM      2  CA  MET A   1      66.953  57.259 -36.531  
    ATOM      3  C   MET A   1      67.370  56.262 -37.627  
    ATOM      4  O   MET A   1      67.105  55.079 -37.531  
    ATOM      5  CB  MET A   1      68.227  57.852 -35.867  
    ATOM      6  CG  MET A   1      67.848  58.995 -34.899  
    ATOM      7  SD  MET A   1      66.880  58.593 -33.421  
    ....      .  ..  ... .   .      ......  ......  ......
    ....      .  ..  ... .   .      ......  ......  ......
    

This file can be read in python using following script.

    
    
    import sys
    x=[];y=[];z=[]
    res=[]
    Nr=0
    for fn in sys.argv[1:]:
       f=open(fn,'r')
       while 1:
          line=f.readline()
          if not line: break
          if line[0:6]=='ATOM  ' :
             rx=float(line[30:38]);ry=float(line[38:46]);rz=float(line[46:54])
             if line[21]=='A' :
                x.append(rx); y.append(ry); z.append(rz)
                Nr=Nr+1
                res.append(line[17:20])
       for i in range(1,Nr-1):
          print fn, i, res[i], x[i], y[i], z[i]
    f.close
    

Now I would like to generate the grid of `N*N*N` dimension and rotate and
translate the molecule on the grid. The rotation and translation can be done
by using FFT (Fast Fourier Transform).

I tried to write something like follows

    
    
    import numpy as np
    import fftw as fft
    
    class Grid3D(object):
    
       def __init__(self, grid_dimension):
          x = y = z = grid_dimension
          self.grid = np.zeros([x, y, z], dtype=float)
    

All this is actually to perform docking of two molecules using 3d grid and
FFT. I wanted to know how to proceed further or any better way?


Answer: An answer to you first question, **"How to read a pdb file"**

If you want to end up with a numpy array, you can use `numpy.genfromtxt` which
is very nice, and much easier to implement and use than your looped reading.
It is also much more robust to spacing of the files, etc.

    
    
    import numpy as np
    data = np.genfromtxt('filename.txt',
            names = 'ATOM,index,res,MET,A,count,x,y,z',
            dtype=['S4',int,'S2','S3','S1',int,float,float,float])
    

Now `data` is a numpy "structured array", which can easily be accessed as
follows:

    
    
    In [13]: data
    Out[13]: 
    array([('ATOM', 1, 'N', 'MET', 'A', 1, 66.104, 56.583, -35.505),
           ('ATOM', 2, 'CA', 'MET', 'A', 1, 66.953, 57.259, -36.531),
           ('ATOM', 3, 'C', 'MET', 'A', 1, 67.37, 56.262, -37.627),
           ('ATOM', 4, 'O', 'MET', 'A', 1, 67.105, 55.079, -37.531),
           ('ATOM', 5, 'CB', 'MET', 'A', 1, 68.227, 57.852, -35.867),
           ('ATOM', 6, 'CG', 'MET', 'A', 1, 67.848, 58.995, -34.899),
           ('ATOM', 7, 'SD', 'MET', 'A', 1, 66.88, 58.593, -33.421)], 
          dtype=[('ATOM', 'S4'), ('index', '<i8'), ('el', 'S2'), ('MET', 'S3'), ('A', 'S1'), ('count', '<i8'), ('x', '<f8'), ('y', '<f8'), ('z', '<f8')])
    
    In [14]: data['x']
    Out[14]: array([ 66.104,  66.953,  67.37 ,  67.105,  68.227,  67.848,  66.88 ])
    
    In [15]: data['y']
    Out[15]: array([ 56.583,  57.259,  56.262,  55.079,  57.852,  58.995,  58.593])
    
    In [16]: data['index']
    Out[16]: array([1, 2, 3, 4, 5, 6, 7])
    
    In [17]: data[3]
    Out[17]: ('ATOM', 4, 'O', 'MET', 'A', 1, 67.105, 55.079, -37.531)
    

