files/fluct_plot.py

   1 #! /usr/bin/env python
   2 import sys
   3 import math
   4
   5 import matplotlib
   6 #matplotlib.use('GTK')
   7 matplotlib.use('Agg')
   8 import matplotlib.pyplot as plt
   9 import os.path
  10
  11 iatom=0
  12 nmol=0
  13 sx={}
  14 sy={}
  15 sz={}
  16 s2x={}
  17 s2y={}
  18 s2z={}
  19 with open(sys.argv[1]) as f:
  20    for line in f:
  21       if line[0:4]=='ATOM' and line[13:15]=='CA':
  22 #         print line
  23          x=float(line[30:38])
  24          y=float(line[38:46])
  25          z=float(line[46:54])
  26 #         print x,y,z
  27          iatom=iatom+1
  28          sx[iatom]=sx.get(iatom,0)+x
  29          sy[iatom]=sy.get(iatom,0)+y
  30          sz[iatom]=sz.get(iatom,0)+z
  31          s2x[iatom]=s2x.get(iatom,0)+x*x
  32          s2y[iatom]=s2y.get(iatom,0)+y*y
  33          s2z[iatom]=s2z.get(iatom,0)+z*z
  34
  35
  36       if line[0:6]=='ENDMDL':
  37          natom=iatom
  38          iatom=0
  39          nmol=nmol+1
  40
  41 x=[]
  42 y=[]
  43 for i in range(1,natom+1):
  44    fluct=math.sqrt((s2x.get(i)-sx.get(i)*sx.get(i)/nmol\
  45                    +s2y.get(i)-sy.get(i)*sy.get(i)/nmol\
  46                    +s2z.get(i)-sz.get(i)*sz.get(i)/nmol  )/nmol)
  47    print fluct
  48    x.append(i)
  49    y.append(fluct)
  50
  51 b=[]
  52 if os.path.exists('plik.pdb'):
  53  with open('plik.pdb') as f:
  54    for line in f:
  55       if line[0:4]=='ATOM' and line[13:15]=='CA':
  56 #         print line
  57          b.append(math.sqrt(float(line[60:66])*3/8/math.pi/math.pi))
  58
  59
  60
  61 plt.xlabel('residue')
  62 plt.ylabel('fluctuations')
  63 plt.xlim(0,natom+1)
  64 plt.plot(x,y,'-')
  65 if len(b)!=0:
  66  plt.plot(x,b,'-',c='red')
  67  plt.legend(['fluctuations','sqrt(3*bfactor/8*pi^2)'])
  68 plt.savefig('fluct_plot.png')
  69