提交: 5e8ffa5a9021639563f64874043fa18c211f95a7 - myprojects-hg-reborn (hg)

Commit MetaInfo

• 修訂: 5e8ffa5a9021639563f64874043fa18c211f95a7 (tree)
• 時間: 2008-01-13 12:15:23
• 作者: iselllo
• Commiter: iselllo

Log Message

Now plot_statistics2.py is also able to calculate the radius of gyration
for a given particle configuration.

Change Summary

• modified: Python-codes/plot_statistics2.py

差異

diff -r 13f04bfb4c63 -r 5e8ffa5a9021 Python-codes/plot_statistics2.py
--- a/Python-codes/plot_statistics2.py	Fri Jan 11 17:31:11 2008 +0000
+++ b/Python-codes/plot_statistics2.py	Sun Jan 13 03:15:23 2008 +0000
@@ -3,6 +3,22 @@
 import scipy as s
 import numpy as n
 import pylab as p
+from rpy import r
+import distance_calc as d_calc
+
+
+
+
+n_part=500
+n_config=1000
+
+density=0.01
+box_size=(n_part/density)**(1./3.)
+print "the box size is, ", box_size
+
+
+threshold=1.06
+
 
 n_clusters=p.load("number_cluster.dat")
 time=p.load("time.dat")
@@ -15,18 +31,6 @@
 p.savefig('n_cluster_loglog.pdf')
 p.hold(False)
 
-coag=n_clusters[0]/n_clusters-1.
-
-
-p.loglog(time,coag,"ro")
-p.xlabel('time')
-p.ylabel('N(0)/N(t)-1')
-#p.legend(('mean x','mean y','mean z',))
-p.title('Log-log plot of cluster population')
-p.grid(True)
-p.savefig('coagulation_loglog.pdf')
-p.hold(False)
-
 
 
 p.loglog(time,n_clusters,linewidth=2.)
@@ -40,12 +44,186 @@
 p.hold(False)
 
 #now I choose to plot the histogram of some cluster distributions
+# and then I investigate its structure
+choose_config=140
 
-choose_config=120
 cluster_name="cluster_dist%05d"%choose_config
 my_cluster=p.load(cluster_name)
 p.hist(my_cluster)
 p.show()
 
 
+#Now I read all the configurations
+
+tot_config=p.load("total_configuration.dat")
+
+	
+#and put them in a suitable shape	 
+	  
+	   
+print "the length of tot_config is, ", len(tot_config)
+tot_config=s.reshape(tot_config,(n_config,3*n_part))
+
+
+#Some arrays I will be using in the following
+
+x_arr=s.zeros((n_config,n_part))
+y_arr=s.zeros((n_config,n_part))
+z_arr=s.zeros((n_config,n_part))
+
+
+for i in xrange(0,n_part):
+        x_arr[:,i]=tot_config[:,3*i]
+        y_arr[:,i]=tot_config[:,(3*i+1)]
+        z_arr[:,i]=tot_config[:,(3*i+2)]
+
+print("OK before r.source")
+
+r.source("cluster_functions.R")
+
+dist_mat=s.zeros((n_part,n_part))
+
+# test_arr=tot_config[choose_config,:]
+# test_arr=s.reshape(test_arr,(n_part,3))
+
+x_pos=x_arr[choose_config,:]
+y_pos=y_arr[choose_config,:]
+z_pos=z_arr[choose_config,:]
+dist_mat=d_calc.distance_calc(x_pos,y_pos,z_pos, box_size)
+
+
+
+
+clust_struc= (r.mycluster2(dist_mat,threshold)) #a cumbersome
+
+#n_clusters[i]=clust_struc[0]
+#print 'the number of clusters in my configuration is, ', clust_struc[0]
+print 'and their composition is, ', clust_struc
+p.save("cluster_structure.dat",clust_struc)
+print 'saved cluster_structure.dat'
+print 'the length of cluster_struc is, ', len(clust_struc) #the number
+#of particles + 1 entry for the number of clusters
+
+#Now I try buiding the structure of the individual clusters
+
+
+
+#print 'initially clust_struc[0] is',clust_struc[0]
+
+#clust_struc=s.transpose(clust_struc)
+
+#cluster_struc=n.integer(clust_struc)
+
+list_part=s.arange(n_part)
+part_in_clust=s.arange(n_part)
+print 'the shape of list_part is, ', s.shape(list_part)
+print 'the shape of clust_struc is, ', s.shape(clust_struc)
+
+print 'len(my_cluster) is, ', len(my_cluster)
+
+
+#MMMMMhhhhh....maybe the following part could become some fortran code
+
+# copy_pos=0
+# for i in xrange(0,len(my_cluster)):
+#    # part_in_clust[:]= s.where(clust_struc[:]==(i*1.),list_part[:],part_in_clust[:])
+# #    count_corr=s.cumsum(my_cluster)
+# #    find_pos=n.int(0+i*(count_corr[i]-count_corr[0]))
+# #    print 'find_pos is, ', find_pos
+# #   find_pos=0
+ 
+#    for j in xrange(0,n_part):
+#        if(clust_struc[j]==i*1.):
+#            part_in_clust[copy_pos]=j
+#            copy_pos=copy_pos+1
+# #       find_pos=find_pos+1
+           
+part_in_clust=s.argsort(clust_struc)
+   
+print 'part_in_clust is, ', part_in_clust
+#Finally, part_in_clust groups together all the particles which are in the same cluster    
+
+
+r_gyr_dist=s.zeros(len(my_cluster)) #this will contain the distribution of the calculated
+#radia of gyration
+
+#pseudocode: loop on each cluster first.
+#define 3 arrays of dimensions related to the number of particles in each cluster
+#then for each cluster work out the radius of gyration as in the case of the
+#plot_statistics.py code.
+
+my_sum=s.cumsum(my_cluster)
+f=s.arange(1) #simply 0 but as an array!
+my_lim=s.concatenate((f,my_sum))
+
+
+#Now a function to calculate the radius of gyration
+def calc_radius(x_arr,y_arr,z_arr,Len):
+    #here x_arr is one-dimensional corresponding to a single configuration
+    r_0j=s.zeros((len(x_arr)-1))
+    for j in xrange(1,len(x_arr)): #so, particle zero is now the reference particle
+        r_0j[j-1]=x_arr[0]-x_arr[j]
+        r_0j[j-1]=-(r_0j[j-1]-Len*n.round(r_0j[j-1]/Len))
+        #r_ij[count_int]=-r_ij[count_int] #I have better reverse the signs now.
+        #print 'i and j are, ', (i+1), (j+1)
+
+
+    #Now I re-define the x_arr in order to be able to take tha variance correctly
+    x_arr[0]=0.
+    x_arr[1:n_part]=r_0j
+
+    #var_x_arr[:]=s.var(r_0j, axis=1)
+    var_x_arr=s.var(x_arr)
+
+    for j in xrange(1,len(y_arr)): #so, particle zero is now the reference particle
+        r_0j[j-1]=y_arr[0]-y_arr[j]
+        r_0j[j-1]=-(r_0j[j-1]-Len*n.round(r_0j[j-1]/Len))
+        #r_ij[count_int]=-r_ij[count_int] #I have better reverse the signs now.
+        #print 'i and j are, ', (i+1), (j+1)
+
+
+    #Now I re-define the x_arr in order to be able to take tha variance correctly
+    y_arr[0]=0.
+    y_arr[1:n_part]=r_0j
+
+    #var_x_arr[:]=s.var(r_0j, axis=1)
+    var_y_arr=s.var(y_arr)
+
+
+
+    for j in xrange(1,len(z_arr)): #so, particle zero is now the reference particle
+        r_0j[j-1]=z_arr[0]-z_arr[j]
+        r_0j[j-1]=-(r_0j[j-1]-Len*n.round(r_0j[j-1]/Len))
+        #r_ij[count_int]=-r_ij[count_int] #I have better reverse the signs now.
+        #print 'i and j are, ', (i+1), (j+1)
+
+
+    #Now I re-define the x_arr in order to be able to take tha variance correctly
+    z_arr[0]=0.
+    z_arr[1:n_part]=r_0j
+
+    #var_x_arr[:]=s.var(r_0j, axis=1)
+    var_z_arr=s.var(z_arr)
+
+    radius=s.sqrt(var_x_arr+var_y_arr+var_z_arr)
+    return radius
+
+    
+                     
+ 
+
+
+       
+for i in xrange(0,len(my_cluster)):
+    x_pos2=x_arr[choose_config,part_in_clust[my_lim[i]:my_lim[i+1]]]
+    y_pos2=x_arr[choose_config,part_in_clust[my_lim[i]:my_lim[i+1]]]
+    z_pos2=x_arr[choose_config,part_in_clust[my_lim[i]:my_lim[i+1]]]
+    r_gyr_dist[i]=calc_radius(x_pos2,y_pos2,z_pos2,box_size)
+
+print 'the distribution of Rg is, ', r_gyr_dist
+p.hist(r_gyr_dist)
+p.show()
+
+        
+        
 print("So far so good")