-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathCentreOfMass.py
128 lines (103 loc) · 2.83 KB
/
CentreOfMass.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
#from pylab import *
import numpy
f=open('output.com')
cont=f.readlines()
radius=17.2
def Expectation(x,p):
tot=0.0 ## p is the y
for i in range(len(x)):
tot+=x[i]*p[i]
return tot/sum(p)
interval=numpy.linspace(0,radius+2,200)
chr_array=[]
for i in range(47):
array=[]
for i in range(len(interval)):
array.append(0)
chr_array.append(array)
timestep=[]
for i in range(3,len(cont),47):
l=cont[i].split()
timestep.append(int(l[0]))
if int(l[0])==2000000:
startPos=i
if int(l[0])==10000000:
endPos=i
print "no of frame", len(timestep),startPos,endPos
noOfFrame=0
for i in range(startPos,endPos,47):
noOfFrame+=1
for j in range(i+1,i+47):
l=cont[j].split()
index=int(l[0])
a=[float(l[1]),float(l[2]),float(l[3])]
dist=numpy.linalg.norm(a)
for k in range(1,len(interval)):
if(dist>=interval[k-1])&(dist<interval[k]):
chr_array[index][k]+=1 ##########
#chromosome[int(l[0])].append(dist)
print "used number of frame ",noOfFrame
for j in range(1,47):
for i in range(len(interval)):
chr_array[j][i]=chr_array[j][i]/float(noOfFrame)
binsize=interval[1]-interval[0]
for j in range(1,47):
sum1=0.0
for i in range(len(interval)):
sum1=sum1+(binsize*chr_array[j][i])
for i in range(len(interval)):
chr_array[j][i]=radius*chr_array[j][i]/sum1
xvec=[]
for i in range(len(interval)):
xvec.append(interval[i]/radius)
ff=open('COM_errorbar.dat','w')
xvec_square=[]
for i in xvec:
xvec_square.append(i**2)
for j in range(1,47):
y=Expectation(xvec,chr_array[j])
variance=Expectation(xvec_square,chr_array[j]) - (y**2)
z=numpy.sqrt(variance)
if j>23:
x=j-23
else:
x=j
ff.write(str(x)+'\t'+str(y)+'\t'+str(z)+'\n')
ff.close()
ff=open('CentreOfMass_distribution.dat','w')
for j in range(len(interval)):
ff.write(str(j+1)+'\t'+str(xvec[j])+'\t')
for i in range(1,47):
ff.write(str(chr_array[i][j])+'\t')
ff.write('\n')
'''
f, axarr = subplots(4, 6)
xx=-1
for i in range(1,24):
yy=(i-1)%6
if yy==0:
xx+=1
axarr[xx,yy].plot(xvec,chr_array[i],'r.-',label='chr '+str(i))
axarr[xx,yy].plot(xvec,chr_array[i+23],'b.-',label='chr '+str(i+23))
axarr[xx,yy].legend(loc='upper right',prop={'size':9})
if xx==3:
axarr[xx,yy].set_xlabel('R')
if yy==0:
axarr[xx,yy].set_ylabel('S(R)')
axarr[xx,yy].tick_params(axis='both', which='major', labelsize=10)
axarr[xx,yy].tick_params(axis='both', which='minor', labelsize=10)
axarr[xx,yy].locator_params(nbins=3)
setp([a.get_xticklabels() for a in axarr[3, 5:6]], visible=False)
setp([a.get_yticklabels() for a in axarr[3:4, 5]], visible=False)
#f.savefig('RadialDis'+str(i)+'_'+str(i+23)+'.png')
f.tight_layout()
#show()
for i in range(1,24):
f=figure()
plot(xvec,chr_array[i],'r.-',label='chr '+str(i))
plot(xvec,chr_array[i+23],'b.-',label='chr '+str(i+23))
legend(loc='upper right')
xlabel('Radius')
ylabel('COM')
f.savefig('COM'+str(i)+'_'+str(i+23)+'.png')
'''