Dinstant Matrix

calcDisMat.py

# GPL v 3.0 
# script name: setupDisMat.py
# written by Naoyuki Miyashita on 09/07/2014 
# Naoyuki Miyashita (c) 2014
#
#-----------------------------------
# Use lProam, pylab, numpy
#-----------------------------------
# ussage:
# python calcDisMat.py 4 input.pdb input.dcd "name CA" output
#
# Output:
# output.ave.d : averaged Distance matrix
# output.div.d : averaged Diviation matrix
# output.dis.d : trajectories of Distance matrix
# output.dif.d : trajectories of Diviation matrix
#-----------------------------------
# Please use prody selection rules in the atom selections.
#
from prody import *
from pylab import *
#import pyopencl as cl
import numpy as np
import lProam
from multiprocessing import Process

#

def argv_read():
    cpuNumber    = int(sys.argv[1])    # number of CPU
    filename     = sys.argv[2]         # input.pdb
    dcdfile      = sys.argv[3]         # input.dcd
    selec1 = sys.argv[4]               #  selections usually "name CA"
    outPrefix = sys.argv[5]               # prefix
    argv_array=(cpuNumber,filename,dcdfile,selec1,outPrefix)
    return(argv_array)


def distances(range1, range2, ncoord, coord):
    allx=np.array([])
    for i in range(range1,range2):
        for j in range(ncoord):
            c=coord[i]-coord[j]
            norm=np.linalg.norm(c)
            allx=np.append(allx,norm)
            k=i*ncoord+j
    return(allx)

def main_calc():
    #input
    argv_array=argv_read()
    (cpuNumber,filename,dcdfile,select1,outPrefix)=argv_array
    #
    (pdb,dcd)=lProam.readFrame(filename,dcdfile)
    (pdb1,dcd1)=lProam.getSelectedFrames(pdb,dcd,select1)
    ncp=([])
    norms=np.array([])
    ndcd1=len(dcd1)
    print "dcd1=",ndcd1

    avefile=str(outPrefix) + ".dis.d"
    p=open(avefile,'w')
    for i, frame in enumerate(dcd1):
        coord=frame.getCoords()
        ncoord=len(coord)
        sw=1
        norms=distances(0,ncoord, ncoord, coord)
        print >> p
        if i == 0:
            stock=np.copy(norms)
        #    stocksq=np.copy(norms*2)
        else:
            stock=stock+norms
        #    stocksq=stocksq+norms**2
        print i,norms
        print >> p, '# ', i
        for s in range(ncoord):
            for j in range(ncoord):
                w=s*ncoord+j
                print >> p, '%d %d %d %f' % (i,s,j,norms[w])
            print >> p, ' '
    p.close()
    stock=stock/ndcd1
    #stocksq=stocksq/ndcd1
    #
    # Output
    avefile=str(outPrefix) + ".ave.d"
    p=open(avefile,'w')
    for i in range(ncoord):
        for j in range(ncoord):
            k=i*ncoord+j
            print >> p, '%d %d %f' %(i,j,stock[k])
        print >> p, ' '
    p.close()
    # v
    #
    difffile=str(outPrefix) + ".dif.d"
    p=open(difffile,'w')
    for i, frame in enumerate(dcd1):
        coord=frame.getCoords()
        ncoord=len(coord)
        sw=0
        print >> p
        norms=distances(0,ncoord, ncoord, coord)
        stock2=stock-norms
        print i,stock2
        if i == 0:
            stock3=np.copy(stock2**2)
        else:
            stock3=stock3+stock2**2
        print >> p, '# ',i
        for s in range(ncoord):
            for j in range(ncoord):
                w=s*ncoord+j
                print >> p, '%d %d %d %f' % (i,s,j,stock2[w])
            print >> p, ' '
    p.close()
    stock3=stock3/ndcd1
    # Output
    avefile=str(outPrefix) + ".div.d"
    p=open(avefile,'w')
    #outp=np.copy(stocksq)
    #outp=stocksq-stock**2
    for i in range(ncoord):
        for j in range(ncoord):
            k=i*ncoord+j
            print >> p, '%d %d %f' %(i,j,stock3[k])
        print >> p, ' '
    p.close()

    return()

if __name__ == '__main__':
    main_calc()