Commit: 7a7995b1bc4a93c6f5a7c2d4692edcd722a2bd0f - myprojects-hg-reborn (hg)

Commit MetaInfo

• Revision: 7a7995b1bc4a93c6f5a7c2d4692edcd722a2bd0f (tree)
• Zeit: 2008-04-16 00:27:54
• Autor: iselllo
• Commiter: iselllo

Log Message

I added the code diffusion_many_monomers.py which is useful to study the
diffusion results of the espresso simulations involving monomers only.

Ändern Zusammenfassung

• add: Python-codes/diffusion_many_monomers.py

Diff

diff -r b4909edfb539 -r 7a7995b1bc4a Python-codes/diffusion_many_monomers.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Python-codes/diffusion_many_monomers.py	Tue Apr 15 15:27:54 2008 +0000
@@ -0,0 +1,1327 @@
+#! /usr/bin/env python
+
+import scipy as s
+import numpy as n
+import pylab as p
+#from rpy import r
+import distance_calc as d_calc
+import igraph as ig
+import hydrodynamic_radius_calc as h
+import calc_rad as c_rad
+
+from scipy import stats #I need this module for the linear fit
+
+
+box_size=10000.  #here the box size does not depend on the total particle number
+#but on the total number of particles in a box
+print "the box size is, ", box_size
+
+#NB: now the box size is fixed a priori and is no longer a derivate object
+
+n_part=3000
+#density=0.03
+
+size_cluster=1 #number of monomers in each cluster.
+
+types=n_part/size_cluster #different kinds of particles i.e. total number of clusters when coagulation
+# is done
+
+print "types is, ", types
+
+
+min_clus_stat=250  # this defines the minimum number of clusters which need to have formed before I start taking some statistics
+
+
+# x_cm=s.zeros(types)   
+# y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
+# z_cm=s.zeros(types)
+
+# v_x_cm=s.zeros(types)
+# v_y_cm=s.zeros(types)  #velocities of the centres of mass of each cluster
+# v_z_cm=s.zeros(types)
+
+
+#Again, I prefer to give the number of particles and the density as input parameters
+# and work out the box size accordingly
+
+
+ini_config=0
+fin_config=250 #for large data post-processing I need to declare an initial and final
+
+               #configuration I want to read and post-process
+
+
+by=1 #this tells how many configurations there are in the file I am reading
+
+figure=0 #whether I sould print many figures or not
+
+
+n_config=(fin_config-ini_config)/by # total number of configurations, which are numbered from 0 to n_config-1
+my_selection=s.arange(ini_config,fin_config,by)
+
+
+threshold=1.06 #threshold to consider to particles as directly connected
+
+
+
+
+time=p.load("time.dat")
+
+time=time[my_selection] #I adapt the time to the initial and final configuration
+
+p.save("time_red.dat",time)
+
+# x_mean=s.zeros(len(time))     # mean position (average on all the clusters) of the cluster CM
+# y_mean=s.zeros(len(time))
+# z_mean=s.zeros(len(time))
+
+# #y_sq_mean=s.zeros(len(time))
+# #z_sq_mean=s.zeros(len(time))
+# #z_sq_mean=s.zeros(len(time))
+
+
+# delta_x_sq=s.zeros(len(time))
+# delta_y_sq=s.zeros(len(time))
+# delta_z_sq=s.zeros(len(time))
+
+
+
+# # v_x_sq_mean=s.zeros(len(time))
+# # v_y_sq_mean=s.zeros(len(time))
+# # v_z_sq_mean=s.zeros(len(time))
+
+# delta_v_x_sq=s.zeros(len(time))
+# delta_v_y_sq=s.zeros(len(time))
+# delta_v_z_sq=s.zeros(len(time))
+
+
+
+#dist_mat=s.zeros((size_cluster,size_cluster)) # I need to find out when the particles of each kind
+#coagulate into a single cluster...I am not interested in looking at at all the particles in the system now
+
+
+
+def R_cm_distr(test_arr,types,n_config,n_part,size_cluster):  #I can use this function both for the statistics on position
+    #and velocity
+
+    count_left=s.arange(0,n_part,size_cluster)
+    count_right=count_left+size_cluster         
+
+
+    
+#     x_cm=s.zeros(types)   
+#     y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
+#     z_cm=s.zeros(types)
+
+    
+#     #print "OK here"
+#     x_pos_temp=test_arr[:,0]
+#     y_pos_temp=test_arr[:,1]
+#     z_pos_temp=test_arr[:,2]
+
+    R_cm_arr=s.zeros((types,3)) #this array will contain the distribution of hydrodynamic radia
+    #calculated for a certain snapshot of the system
+
+
+    for i in xrange(types):
+        r_pos=test_arr[count_left[i]:count_right[i],:]
+#         x_pos=x_pos_temp[count_left[i]:count_right[i]]
+#         y_pos=y_pos_temp[count_left[i]:count_right[i]]
+#         z_pos=z_pos_temp[count_left[i]:count_right[i]]
+
+     
+
+#        print "s.shape(r_pos) is, ", s.shape(r_pos)
+        
+        R_cm_arr[i,:]=r_pos.mean(axis=0)
+
+
+    return R_cm_arr
+
+
+
+
+
+
+def R_h_distr(test_arr,types,n_config,n_part,size_cluster):  #I can use this function both for the statistics on position
+    #and velocity
+
+    count_left=s.arange(0,n_part,size_cluster)
+    count_right=count_left+size_cluster         
+
+
+    
+    x_cm=s.zeros(types)   
+    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
+    z_cm=s.zeros(types)
+
+    
+    #print "OK here"
+    x_pos_temp=test_arr[:,0]
+    y_pos_temp=test_arr[:,1]
+    z_pos_temp=test_arr[:,2]
+
+    R_hydro_arr=s.zeros(types) #this array will contain the distribution of hydrodynamic radia
+    #calculated for a certain snapshot of the system
+
+
+    for i in xrange(types):
+        x_pos=x_pos_temp[count_left[i]:count_right[i]]
+        y_pos=y_pos_temp[count_left[i]:count_right[i]]
+        z_pos=z_pos_temp[count_left[i]:count_right[i]]
+
+        #print "x_pos is, ", x_pos
+        
+        R_hydro_arr[i]=h.hydro_radius_calc(x_pos,y_pos,z_pos,box_size)
+   #     print "R_hydro_arr[i] is, ", R_hydro_arr[i]
+        
+    R_h=R_hydro_arr.mean()
+
+
+    return R_h
+
+
+
+
+
+def R_gyr_distr(test_arr,types,n_config,n_part,size_cluster):  #I can use this function both for the statistics on position
+    #and velocity
+
+    count_left=s.arange(0,n_part,size_cluster)
+    count_right=count_left+size_cluster         
+
+
+    
+    x_cm=s.zeros(types)   
+    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
+    z_cm=s.zeros(types)
+
+    
+    #print "OK here"
+    x_pos_temp=test_arr[:,0]
+    y_pos_temp=test_arr[:,1]
+    z_pos_temp=test_arr[:,2]
+
+    R_gyr_arr=s.zeros(types) #this array will contain the distribution of hydrodynamic radia
+    #calculated for a certain snapshot of the system
+
+
+    for i in xrange(types):
+        x_pos=x_pos_temp[count_left[i]:count_right[i]]
+        y_pos=y_pos_temp[count_left[i]:count_right[i]]
+        z_pos=z_pos_temp[count_left[i]:count_right[i]]
+
+        #print "x_pos is, ", x_pos
+        
+        R_gyr_arr[i]=c_rad.rad_gyr(x_pos,y_pos,z_pos,box_size)
+   #     print "R_hydro_arr[i] is, ", R_hydro_arr[i]
+        
+    R_gyr=R_gyr_arr.mean()
+
+
+    return R_gyr
+
+
+
+
+def calc_time_stat(test_arr,types,n_config,n_part,size_cluster,box_size):  #I can use this function both for the statistics on position
+    #and velocity
+
+    count_left=s.arange(0,n_part,size_cluster)
+    count_right=count_left+size_cluster         
+
+
+    
+    x_cm=s.zeros(types)   
+    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
+    z_cm=s.zeros(types)
+
+
+    #fold everything inside the box
+
+    test_arr=s.remainder(test_arr,box_size)
+    
+    #print "OK here"
+    x_pos=test_arr[:,0]
+    y_pos=test_arr[:,1]
+    z_pos=test_arr[:,2]
+
+
+    results=s.zeros(7)
+
+
+    for j in xrange(types):
+        x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
+        y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
+        z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])
+
+
+
+
+
+    x_mean=x_cm.mean()
+    y_mean=y_cm.mean()
+    z_mean=z_cm.mean()
+
+    delta_x_sq=x_cm.var()
+    delta_y_sq=y_cm.var()
+    delta_z_sq=z_cm.var()
+
+    delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq
+
+
+    results[0]=x_mean
+    results[1]=y_mean
+    results[2]=z_mean
+
+    results[3]=delta_x_sq
+    results[4]=delta_y_sq
+    results[5]=delta_z_sq
+    results[6]=delta_r_sq
+
+
+    return results
+
+def calc_time_stat_vel(test_arr,types,n_config,n_part,size_cluster,box_size):  #I can use this function both for the statistics on position
+    #and velocity
+
+    count_left=s.arange(0,n_part,size_cluster)
+    count_right=count_left+size_cluster         
+
+
+    
+    x_cm=s.zeros(types)   
+    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
+    z_cm=s.zeros(types)
+
+#here I have the only difference wrt the position fluctuations: I do not have to fold the velocity!!!
+ #    #fold everything inside the box
+
+#     test_arr=s.remainder(test_arr,box_size)
+    
+    #print "OK here"
+    x_pos=test_arr[:,0]
+    y_pos=test_arr[:,1]
+    z_pos=test_arr[:,2]
+
+
+    results=s.zeros(7)
+
+
+    for j in xrange(types):
+        x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
+        y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
+        z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])
+
+
+
+
+
+    x_mean=x_cm.mean()
+    y_mean=y_cm.mean()
+    z_mean=z_cm.mean()
+
+    delta_x_sq=x_cm.var()
+    delta_y_sq=y_cm.var()
+    delta_z_sq=z_cm.var()
+
+    delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq
+
+
+    results[0]=x_mean
+    results[1]=y_mean
+    results[2]=z_mean
+
+    results[3]=delta_x_sq
+    results[4]=delta_y_sq
+    results[5]=delta_z_sq
+    results[6]=delta_r_sq
+
+
+    return results
+
+
+
+
+
+
+
+
+def calc_time_stat_clusters_together(test_arr,types,n_config,n_part,size_cluster):  #I can use this function both for the statistics on position
+    #and velocity
+
+    count_left=s.arange(0,n_part,size_cluster)
+    count_right=count_left+size_cluster         
+
+
+    
+    x_cm=s.zeros(types)   
+    y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
+    z_cm=s.zeros(types)
+
+    
+    #print "OK here"
+    x_pos=test_arr[:,0]
+    y_pos=test_arr[:,1]
+    z_pos=test_arr[:,2]
+
+
+    results=s.zeros(7)
+
+
+
+
+
+
+
+    x_mean=x_pos.mean()
+    y_mean=y_pos.mean()
+    z_mean=z_pos.mean()
+
+    delta_x_sq=x_pos.var()
+    delta_y_sq=y_pos.var()
+    delta_z_sq=z_pos.var()
+
+    delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq
+
+
+    results[0]=x_mean
+    results[1]=y_mean
+    results[2]=z_mean
+
+    results[3]=delta_x_sq
+    results[4]=delta_y_sq
+    results[5]=delta_z_sq
+    results[6]=delta_r_sq
+
+
+    return results
+
+
+
+
+
+
+
+
+
+
+
+
+def calc_time_stat_selection(test_arr,types,n_config,n_part,size_cluster,clust_choice):
+#I can use this function both for the statistics on position
+    #and velocity
+
+
+    
+    x_cm=s.zeros(len(s.ravel(clust_choice)))   
+    y_cm=s.zeros(len(s.ravel(clust_choice)))    #positions of the centres of mass of each cluster
+    z_cm=s.zeros(len(s.ravel(clust_choice)))
+
+
+   # print "len(clust_choice) is, ", len(clust_choice)
+    count_left=s.arange(0,(len(s.ravel(clust_choice))*size_cluster),size_cluster)
+    count_right=count_left+size_cluster         
+
+
+
+    
+    #print "OK here"
+    x_pos_temp=test_arr[:,0]
+    y_pos_temp=test_arr[:,1]
+    z_pos_temp=test_arr[:,2]
+
+    #print "len(x_pos_temp) is, ", len(x_pos_temp)
+
+    x_pos=s.zeros(len(s.ravel(clust_choice))*size_cluster)
+    y_pos=s.zeros(len(s.ravel(clust_choice))*size_cluster)
+    z_pos=s.zeros(len(s.ravel(clust_choice))*size_cluster)
+
+    #print "len(x_pos) is,", len(x_pos)
+
+    for l in xrange(len(s.ravel(clust_choice))):
+     #   print "count_left and count_right are, ", count_left[l],count_right[l]
+
+        #print "clust_choice[l] is, ", clust_choice[l]
+        x_pos[count_left[l]:count_right[l]]=\
+        x_pos_temp[(s.ravel(clust_choice)[l]*size_cluster):(s.ravel(clust_choice)[l]*size_cluster+size_cluster)]
+        y_pos[count_left[l]:count_right[l]]=\
+                 y_pos_temp[(s.ravel(clust_choice)[l]*size_cluster):(s.ravel(clust_choice)[l]*size_cluster+size_cluster)]
+
+        z_pos[count_left[l]:count_right[l]]=\
+                 z_pos_temp[(s.ravel(clust_choice)[l]*size_cluster):(s.ravel(clust_choice)[l]*size_cluster+size_cluster)]
+
+
+
+
+
+    
+    #Then all the rest can stay the same
+
+    results=s.zeros(7)
+
+
+
+    for j in xrange(len(s.ravel(clust_choice))): #I iterate on the existing fully-coagulated clusters
+        #print "j is, ", j
+        x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
+        y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
+        z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])
+
+
+
+
+
+    x_mean=x_cm.mean()
+    y_mean=y_cm.mean()
+    z_mean=z_cm.mean()
+
+    delta_x_sq=x_cm.var()
+    delta_y_sq=y_cm.var()
+    delta_z_sq=z_cm.var()
+
+    delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq
+
+
+    results[0]=x_mean
+    results[1]=y_mean
+    results[2]=z_mean
+
+    results[3]=delta_x_sq
+    results[4]=delta_y_sq
+    results[5]=delta_z_sq
+    results[6]=delta_r_sq
+
+
+    return results
+
+
+
+
+
+
+
+def coagulation_clusters(test_arr,types,n_config,box_size,size_cluster,threshold):
+
+    count_left=s.arange(0,n_part,size_cluster)
+    count_right=count_left+size_cluster         
+    
+    coag_arr=s.zeros(types)
+
+    
+#     x_cm=s.zeros(types)   
+#     y_cm=s.zeros(types)    #positions of the centres of mass of each cluster
+#     z_cm=s.zeros(types)
+
+    
+    #print "OK here"
+    x_pos=test_arr[:,0]
+    y_pos=test_arr[:,1]
+    z_pos=test_arr[:,2]
+
+
+    results=s.zeros(7)
+
+
+    for j in xrange(types):
+        x_pos_clu=x_pos[count_left[j]:count_right[j]]
+        y_pos_clu=y_pos[count_left[j]:count_right[j]]
+        z_pos_clu=z_pos[count_left[j]:count_right[j]]
+
+        dist_mat=d_calc.distance_calc(x_pos_clu,y_pos_clu,z_pos_clu, box_size)
+        cluster_obj=ig.Graph.Adjacency((dist_mat <= threshold).tolist(),\
+					       ig.ADJ_UNDIRECTED)
+
+        cluster_obj.simplify()
+        clustering=cluster_obj.clusters()
+        
+        #n_clusters[i]=p.load("nc_temp.dat")
+#        print "clustering is, ", clustering
+
+        coag_arr[j]=len(clustering)
+
+   # print "coag_arr is, ", coag_arr
+    return coag_arr
+
+
+######################################################################
+######################################################################
+######################################################################
+######################################################################
+
+
+
+
+
+
+
+# #print "Ok here"
+
+# for i in xrange(0,n_config):
+
+#     read_pos="read_pos_%1d"%my_selection[i]
+#     read_vel="read_vel_%1d"%my_selection[i]
+
+#     #print "OK here 2"
+#     #cluster_name="cluster_dist%05d"%my_selection[i]
+#     test_arr=p.load(read_pos)
+    
+#     #print "OK here"
+#     x_pos=test_arr[:,0]
+#     y_pos=test_arr[:,1]
+#     z_pos=test_arr[:,2]
+
+#     for j in xrange(types):
+#         x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
+#         y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
+#         z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])
+
+#         v_x_cm[j]=s.mean(x_pos[count_left[j]:count_right[j]])
+#         v_y_cm[j]=s.mean(y_pos[count_left[j]:count_right[j]])
+#         v_z_cm[j]=s.mean(z_pos[count_left[j]:count_right[j]])
+
+
+
+#     x_mean[i]=x_cm.mean()
+#     y_mean[i]=y_cm.mean()
+#     z_mean[i]=z_cm.mean()
+
+#     delta_x_sq[i]=x_cm.var()
+#     delta_y_sq[i]=y_cm.var()
+#     delta_z_sq[i]=z_cm.var()
+
+
+
+# delta_r_sq=delta_x_sq+delta_y_sq+delta_z_sq
+
+
+# fig = p.figure()
+# axes = fig.gca()
+
+
+
+# axes.plot(time,x_mean,"ro",\
+#           label="<x>",linewidth=2.)
+# axes.plot(time,y_mean,"k+",\
+#           label="<y>")
+# axes.plot(time,z_mean,"mx",\
+#           label="<z> ",linewidth=2.)
+# #p.ylim(1.4,2.6)
+# #axes.legend()
+# p.xlabel('Time')
+# p.ylabel('Mean position')
+# p.savefig("mean_x_y_z.pdf")
+ 
+# p.clf()    
+
+        
+   
+
+# fig = p.figure()
+# axes = fig.gca()
+
+
+
+# axes.plot(time,delta_r_sq,"ro",\
+#           label="<x>",linewidth=2.)
+# # axes.plot(time,y_mean,"k+",\
+# #           label="<y>")
+# # axes.plot(time,z_mean,"mx",\
+# #           label="<z> ",linewidth=2.)
+# #p.ylim(1.4,2.6)
+# #axes.legend()
+# p.xlabel('Time')
+# p.ylabel('delta position squared')
+# p.savefig("delta_r_squared.pdf")
+ 
+# p.clf()    
+
+
+    
+#     #print "the length of x_pos is, ", len(x_pos)
+
+        
+#########################################################################################
+#########################################################################################
+#########################################################################################
+
+#initialization of some arrays I will have to take in & out the loop
+
+
+#pos_results=s.zeros(4) #just to initialize it
+#pos_results_selection=s.zeros(4) #just to initialize it
+#pos_statistics_selection=s.zeros(2)  # I just need to initialize it!
+time_selection=s.zeros(2)
+
+
+
+
+flag_ini=0 # I will need this flag to check is some other array has been initialized
+
+R_h_arr=s.zeros(n_config)
+
+R_gyr_arr=s.zeros(n_config)
+
+
+for i in xrange(0,n_config):
+
+    read_pos="read_pos_%1d"%my_selection[i]
+    read_vel="read_vel_%1d"%my_selection[i]
+
+    #print "OK here 2"
+    #cluster_name="cluster_dist%05d"%my_selection[i]
+    pos_arr=p.load(read_pos)
+    vel_arr=p.load(read_vel)
+        
+    pos_statistics=calc_time_stat(pos_arr,types,n_config,n_part,size_cluster, box_size)        
+
+#    pos_statistics_clusters_together=calc_time_stat_clusters_together(pos_arr,types,n_config,n_part,size_cluster)        
+
+
+    vel_statistics=calc_time_stat_vel(vel_arr,types,n_config,n_part,size_cluster,box_size)        
+    clust_stat=coagulation_clusters(pos_arr,types,n_config,box_size,size_cluster,threshold)
+    #print "clust_stat is, ", clust_stat
+    clust_choice=s.where(clust_stat[:]==1.) #i.e. I select the clusters which have already finished
+    # coagulating
+
+#     R_h_arr[i]=R_h_distr(pos_arr,types,n_config,n_part,size_cluster)
+#     R_gyr_arr[i]=R_gyr_distr(pos_arr,types,n_config,n_part,size_cluster)
+#     R_cm=R_cm_distr(pos_arr,types,n_config,n_part,size_cluster)
+
+#     if (i==0):
+#         R_cm_arr=R_cm
+#     elif (i>0):
+#         R_cm_arr=s.hstack((R_cm_arr,R_cm))
+
+    print "R_gyr_arr[i] is, ", R_gyr_arr[i]
+    print "R_h_arr[i] is, ", R_h_arr[i]
+    
+    #print "len(clust_choice) is, ", len(s.ravel(clust_choice))
+    #print "clust_choice is, ", clust_choice
+
+    if (len(s.ravel(clust_choice))>=min_clus_stat):
+        pos_statistics_selection=\
+             calc_time_stat_selection(pos_arr,types,n_config,n_part,size_cluster,clust_choice)
+     #   print "pos_statistics_selection is, ", pos_statistics_selection
+
+    
+
+
+
+    if (i==0):
+        
+        pos_results=pos_statistics
+        #R_h=R_h_temp
+
+#        pos_results_clusters_together=pos_statistics_clusters_together
+
+
+#         print "pos_results is, ", pos_results
+#         print "the len of pos_results is, ", len(pos_results)
+        vel_results=vel_statistics
+        clust_results=clust_stat
+        #The following if condition is academic since I'll never have the right number of
+        #clusters at the beginning
+        if (len(s.ravel(clust_choice))>=min_clus_stat):
+            pos_results_selection=pos_statistics_selection
+            time_selection=[time[i]]
+            
+        
+            
+    elif(i!=0):
+#         print "now pos_results is, ", pos_results
+#         print "now the shape of pos_results is, ", s.shape(pos_results)
+#         print "now pos_statistics is, ", pos_statistics
+
+        pos_results=s.concatenate((pos_results,pos_statistics))
+        #R_h=s.concatenate((R_h,R_h_temp))
+#        pos_results_clusters_together=s.concatenate((pos_results_clusters_together,pos_statistics_clusters_together))
+
+
+
+#         print "OK here"
+        vel_results=s.concatenate((vel_results,vel_statistics))
+        clust_results=s.concatenate((clust_results,clust_stat))
+
+
+        if ((len(s.ravel(clust_choice))>=min_clus_stat) and (flag_ini==0)):
+            pos_results_selection=pos_statistics_selection
+            time_selection=[time[i]]
+            flag_ini=1 #to be sure I am doing this only the fist time, then I can concatenate
+
+
+        elif ((len(s.ravel(clust_choice))>=min_clus_stat) and (flag_ini!=0)):
+            pos_results_selection=s.concatenate((pos_results_selection,pos_statistics_selection))
+            time_selection=s.concatenate((time_selection,[time[i]] ))
+
+
+
+pos_results.shape=(-1,7)
+vel_results.shape=(-1,7)  #where -1 means "do whatever is needed to reshape it right"
+
+
+
+#pos_results_clusters_together.shape=(-1,7)
+
+
+
+
+
+clust_results.shape=(-1,types)
+
+
+#pos_results_selection.shape=(-1,7)
+
+
+p.save("position_cluster_centres_of_mass.dat", pos_results)
+p.save("velocity_cluster_centres_of_mass.dat", vel_results)
+p.save("evolution_coagulation_of_each_particle_kind.dat", clust_results)   
+
+
+
+fig = p.figure()
+axes = fig.gca()
+
+
+axes.plot(time,pos_results[:,6],"ro",\
+          label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
+# axes.plot(time,y_mean,"k+",\
+#           label="<y>")
+# axes.plot(time,z_mean,"mx",\
+#           label="<z> ",linewidth=2.)
+#p.ylim(1.4,2.6)
+axes.legend()
+p.xlabel('Time')
+p.ylabel('delta position squared')
+p.savefig("delta_r_squared.pdf")
+ 
+p.clf()    
+
+
+
+
+
+
+#print "pos_results_clusters_together is, ", pos_results_clusters_together
+
+
+# fig = p.figure()
+# axes = fig.gca()
+
+
+# axes.plot(time,pos_results_clusters_together[:,6],"ro",\
+#           label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
+# # axes.plot(time,y_mean,"k+",\
+# #           label="<y>")
+# # axes.plot(time,z_mean,"mx",\
+# #           label="<z> ",linewidth=2.)
+# #p.ylim(1.4,2.6)
+# axes.legend()
+# p.xlabel('Time')
+# p.ylabel('delta position squared')
+# p.savefig("delta_r_squared_clusters_together.pdf")
+ 
+# p.clf()    
+
+
+
+
+
+
+# fig = p.figure()
+# axes = fig.gca()
+
+
+# axes.plot(time_selection,pos_results_selection[:,6],"ro",\
+#           label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
+# # axes.plot(time,y_mean,"k+",\
+# #           label="<y>")
+# # axes.plot(time,z_mean,"mx",\
+# #           label="<z> ",linewidth=2.)
+# #p.ylim(1.4,2.6)
+# axes.legend()
+# p.xlabel('Time')
+# p.ylabel('delta position squared (selection)')
+# p.savefig("delta_r_squared_selection.pdf")
+ 
+# p.clf()    
+
+
+
+
+
+
+fig = p.figure()
+axes = fig.gca()
+
+axes.plot(time,pos_results[:,0],"ro",\
+          label="<x>",linewidth=2.)
+axes.plot(time,pos_results[:,1],"k+",\
+          label="<y>")
+axes.plot(time,pos_results[:,2],"mx",\
+          label="<z> ",linewidth=2.)
+#p.ylim(1.4,2.6)
+#axes.legend()
+p.xlabel('Time')
+p.ylabel('Mean position')
+p.savefig("mean_x_y_z.pdf")
+ 
+p.clf()    
+
+
+
+
+
+
+   
+
+fig = p.figure()
+axes = fig.gca()
+
+
+axes.plot(time,vel_results[:,6],"ro",\
+          label="<v_x^2+v_y^2+v_z^2>",linewidth=2.)
+# axes.plot(time,y_mean,"k+",\
+#           label="<y>")
+# axes.plot(time,z_mean,"mx",\
+#           label="<z> ",linewidth=2.)
+#p.ylim(1.4,2.6)
+axes.legend()
+p.xlabel('Time')
+p.ylabel('delta velocity squared')
+p.savefig("delta_v_squared.pdf")
+ 
+p.clf()    
+
+
+
+
+
+fig = p.figure()
+axes = fig.gca()
+
+axes.plot(time,vel_results[:,0],"ro",\
+          label="<v_x>",linewidth=2.)
+axes.plot(time,vel_results[:,1],"k+",\
+          label="<v_y>")
+axes.plot(time,vel_results[:,2],"mx",\
+          label="<v_z> ",linewidth=2.)
+#p.ylim(1.4,2.6)
+#axes.legend()
+p.xlabel('Time')
+p.ylabel('Mean velocity')
+p.savefig("mean_vel_x_vel_y_vel_z.pdf")
+ 
+p.clf()    
+
+
+
+fig = p.figure()
+axes = fig.gca()
+
+
+n_clus=s.sum(clust_results,axis=1)
+
+print "the number of clusters is, ", n_clus
+
+
+axes.plot(time,n_clus,"ro",\
+          label="evolution total number of clusters")
+types_arr=s.zeros(len(time))
+types_arr[:]=types
+
+axes.plot(time,types_arr,"b",\
+          label="coagulation of all the individual clusters",linewidth=2.)
+
+
+# axes.plot(time,y_mean,"k+",\
+#           label="<y>")
+# axes.plot(time,z_mean,"mx",\
+#           label="<z> ",linewidth=2.)
+#p.ylim(1.4,2.6)
+axes.legend()
+p.xlabel('Time')
+p.ylabel('Number of clusters')
+p.savefig("Coagulation_individual_clusters.pdf")
+ 
+p.clf()    
+
+#Now I think about a linear fit of delta_r_sq vs time
+
+#         lin_fit=stats.linregress(s.log10(n_time_small),s.log10(r_time_small))
+#         my_fra_low[count_config]=1./lin_fit[0]
+#         r_stat_low[count_config]=lin_fit[2]
+
+#I select the results only for the case when all the clusters are formed
+#print "n_clus is, ", n_clus
+#print "the shape of n_clus is, ", s.shape(n_clus)
+cluster_formed=s.where(n_clus[:]==types)
+#print "cluster_formed is, ", cluster_formed
+#print "the shape of cluster_formed is, ", s.shape(cluster_formed)
+
+time_formed=time[cluster_formed]
+#print "time_formed is, ", time_formed
+my_delta_r_sq=pos_results[:,6]
+delta_r_sq_formed=my_delta_r_sq[cluster_formed]
+
+
+lin_fit=stats.linregress(time_formed,delta_r_sq_formed)
+
+delta_r_sq_fitted=lin_fit[0]*time_formed+lin_fit[1]
+
+
+
+Diff_coeff=lin_fit[0]
+
+print "the linear coefficient [after cluster formation] is, ", lin_fit[0]
+
+
+
+
+# fig = p.figure()
+# axes = fig.gca()
+
+
+# axes.plot(time,pos_results[:,6],"ro",\
+#           label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
+# # axes.plot(time,y_mean,"k+",\
+# #           label="<y>")
+# # axes.plot(time,z_mean,"mx",\
+# #           label="<z> ",linewidth=2.)
+# #p.ylim(1.4,2.6)
+
+# axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)
+
+# axes.legend()
+# p.xlabel('Time')
+# p.ylabel('delta position squared and fit')
+# p.savefig("delta_r_squared_and_fit.pdf")
+ 
+# p.clf()
+
+
+
+
+# my_delta_r_sq=pos_results_clusters_together[:,6]
+# delta_r_sq_formed=my_delta_r_sq[cluster_formed]
+
+
+# lin_fit=stats.linregress(time_formed,delta_r_sq_formed)
+
+# delta_r_sq_fitted=lin_fit[0]*time_formed+lin_fit[1]
+
+
+
+# fig = p.figure()
+# axes = fig.gca()
+
+
+# axes.plot(time,pos_results_clusters_together[:,6],"ro",\
+#           label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
+# # axes.plot(time,y_mean,"k+",\
+# #           label="<y>")
+# # axes.plot(time,z_mean,"mx",\
+# #           label="<z> ",linewidth=2.)
+# #p.ylim(1.4,2.6)
+
+# axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)
+
+# axes.legend()
+# p.xlabel('Time')
+# p.ylabel('delta position squared and fit')
+# p.savefig("delta_r_squared_and_fit_clusters_together.pdf")
+ 
+# p.clf()    
+
+# print "the linear coefficient is [clusters_together], ", lin_fit[0]
+
+
+#Now the dynamics before the formation of the clusters
+
+
+cluster_unformed=s.where(n_clus[:]!=types)
+#print "cluster_formed is, ", cluster_formed
+#print "the shape of cluster_formed is, ", s.shape(cluster_formed)
+
+time_unformed=time[cluster_unformed]
+#print "time_formed is, ", time_formed
+#my_delta_r_sq=pos_results[:,6]
+delta_r_sq_unformed=my_delta_r_sq[cluster_unformed]
+
+
+lin_fit=stats.linregress(time_unformed,delta_r_sq_unformed)
+
+delta_r_sq_fitted_unformed=lin_fit[0]*time_unformed+lin_fit[1]
+
+
+print "the linear coefficient before forming the clusters is, ", lin_fit[0]
+
+
+
+coag_choice=s.where((n_clus/n_clus[0])>0.2)
+
+
+
+time_coag=time[coag_choice]
+#print "time_formed is, ", time_formed
+#my_delta_r_sq=pos_results[:,6]
+delta_r_sq_coag=my_delta_r_sq[coag_choice]
+
+
+print "delta_r_sq_coag is, ", delta_r_sq_coag
+
+
+lin_fit=stats.linregress(time_coag,delta_r_sq_coag)
+
+delta_r_sq_coag_fitted=lin_fit[0]*time_coag+lin_fit[1]
+
+
+
+
+print "the linear coefficient at the early stages of coagulation is, ", lin_fit[0]
+
+
+
+
+
+
+
+
+fig = p.figure()
+axes = fig.gca()
+
+
+axes.plot(time,pos_results[:,6],"ro",\
+          label="<delta_x^2+delta_y^2+delta_z^2>",linewidth=2.)
+# axes.plot(time,y_mean,"k+",\
+#           label="<y>")
+# axes.plot(time,z_mean,"mx",\
+#           label="<z> ",linewidth=2.)
+#p.ylim(1.4,2.6)
+
+axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)
+
+axes.plot(time_unformed, delta_r_sq_fitted_unformed,"k",linewidth=2.)
+
+axes.plot(time_coag, delta_r_sq_coag_fitted,"m",linewidth=2.)
+
+
+
+axes.legend()
+p.xlabel('Time')
+p.ylabel('delta position squared and fit')
+p.savefig("delta_r_squared_and_fit.pdf")
+ 
+p.clf()    
+
+
+
+
+fig = p.figure()
+axes = fig.gca()
+
+
+axes.plot(time,R_h_arr,"ro",\
+          label="hydrodynamic radius",linewidth=2.)
+# axes.plot(time,y_mean,"k+",\
+#           label="<y>")
+# axes.plot(time,z_mean,"mx",\
+#           label="<z> ",linewidth=2.)
+#p.ylim(1.4,2.6)
+
+# axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)
+
+# axes.plot(time_unformed, delta_r_sq_fitted_unformed,"k",linewidth=2.)
+
+# axes.plot(time_coag, delta_r_sq_coag_fitted,"m",linewidth=2.)
+
+
+
+axes.legend()
+p.xlabel('Time')
+p.ylabel('R_h')
+p.savefig("hydrodynamic_radius.pdf")
+ 
+p.clf()    
+
+
+
+
+
+fig = p.figure()
+axes = fig.gca()
+
+
+axes.plot(time,R_gyr_arr,"ro",\
+          label="gyration radius",linewidth=2.)
+# axes.plot(time,y_mean,"k+",\
+#           label="<y>")
+# axes.plot(time,z_mean,"mx",\
+#           label="<z> ",linewidth=2.)
+#p.ylim(1.4,2.6)
+
+# axes.plot(time_formed, delta_r_sq_fitted,"b",linewidth=2.)
+
+# axes.plot(time_unformed, delta_r_sq_fitted_unformed,"k",linewidth=2.)
+
+# axes.plot(time_coag, delta_r_sq_coag_fitted,"m",linewidth=2.)
+
+
+
+axes.legend()
+p.xlabel('Time')
+p.ylabel('R_gyr')
+p.savefig("gyration_radius.pdf")
+ 
+p.clf()    
+
+
+
+
+#print "R_h_arr is, ", R_h_arr
+
+#Now I evaluate the ratio between the hydrodynamic radius and the radius of gyration as time goes by and I
+# evaluate the 
+
+
+# R_h_over_R_gyr=s.mean(R_h_arr[cluster_formed]/R_gyr_arr[cluster_formed])
+
+
+# R_h_mean=R_h_arr[cluster_formed].mean()
+
+
+# R_gyr_mean=R_gyr_arr[cluster_formed].mean()
+
+
+# print "the ratio hydrodynamic radius over radius of gyration is, " ,R_h_over_R_gyr
+
+
+
+
+
+
+#Now I take some time averages for a single cluster.
+
+if (types==1): #i.e. I have a single cluster
+    delta_r_sq_time_aver=s.var(pos_results[:,0:3], axis=0)
+    delta_r_sq_time_aver=delta_r_sq_time_aver.sum()
+    v_mean_time=s.mean(vel_results[:,0:3],axis=0)
+
+    print "delta_r_sq_time_aver is, ", delta_r_sq_time_aver
+    print "v_mean_time is, ", v_mean_time
+    print "vel results are, ", vel_results
+
+    print "the shape of vel_results is, ", s.shape(vel_results)
+
+
+    print "<delta x squared>, ", pos_results[:,0].var()
+    print "<delta y squared>, ", pos_results[:,1].var()
+    print "<delta z squared>, ", pos_results[:,2].var()
+
+    print "<x>, ", pos_results[:,0].mean()
+    print "<y>, ", pos_results[:,1].mean()
+    print "<z>, ", pos_results[:,2].mean()
+
+
+
+
+    print "<v_x>, ", vel_results[:,0].mean()
+    print "<v_y>, ", vel_results[:,1].mean()
+    print "<v_z>, ", vel_results[:,2].mean()
+    
+
+    delta_v_sq_time_aver=s.var(vel_results[:,0:3], axis=0)
+    delta_v_sq_time_aver=delta_v_sq_time_aver.sum()
+
+    print "<delta v_x squared>, ", vel_results[:,0].var()
+    print "<delta v_y squared>, ", vel_results[:,1].var()
+    print "<delta v_z squared>, ", vel_results[:,2].var()
+    
+
+    
+
+    print "delta_v_sq_time_aver is, ", delta_v_sq_time_aver
+
+
+#Now I calculate the diffusion coefficient from the collected cluster positions
+
+print "s.shape(R_cm_arr) is, ", s.shape(R_cm_arr)
+
+#displacement=(R_cm_arr[:,0]-R_cm_arr)**2.
+
+
+my_len=s.shape(R_cm_arr)[1]/3
+
+p.save("centre_of_mass_positions.dat", R_cm_arr)
+
+
+displacement=s.zeros((s.shape(R_cm_arr)[0],my_len))
+
+print "R_cm_arr is, ", R_cm_arr[:,0:10]
+
+for i in xrange(my_len):
+    displacement[:,i]=s.sum(((R_cm_arr[:,0:3]-R_cm_arr[:,(i*3):((i+1)*3)])**2.),axis=1)
+
+print "displacement is, ", displacement[:, 0:3]
+
+p.save("displacement.dat", displacement)
+
+# x_pos=x.zeros(my_len)
+# y_pos=x.zeros(my_len)
+# z_pos=x.zeros(my_len)
+
+# for i in xrange(my_len):
+#     x_pos=
+    
+
+displacement=displacement/6./time.max()
+
+D_distr=displacement.mean(axis=1)
+
+print "the trajectory-averaged distribution of diffusion coeff is, ", D_distr
+
+D_mean=D_distr.mean()
+
+print "the trajectory-averaged diffusion coeff is, ", D_mean
+
+print "and its std is, ", D_distr.std()
+
+
+
+
+
+def Diff(kT,mu,r_eff):
+    D=kT/(6.*s.pi*mu*r_eff)
+    return D
+
+def R_mob(kT,mu, D):
+    r_eff=kT/(6.*s.pi*mu*D)
+
+    return r_eff
+
+
+
+
+kT=0.5
+mu=1./(3.*s.pi)
+
+# Diff_R_h=Diff(kT,mu,R_h_mean)
+
+
+# Diff_R_gyr=Diff(kT,mu,R_gyr_mean)
+
+# print "D calculated with R_gyr is, ", Diff_R_gyr
+
+# print "D calculated with R_h is, ", Diff_R_h
+
+# R_mobility=R_mob(kT,mu,D_mean)
+
+# print "The effective mobility radius is ", R_mobility
+
+# print "the ratio R_mobility/R_gyration is, ", R_mobility/R_gyr_mean
+
+# print "the ratio R_mobility/R_hydro_mean is, ", R_mobility/R_h_mean
+
+
+
+print "So far so good"