myprojects-hg-reborn

Fork

#! /usr/bin/env python

import scipy as s
import numpy as n
import pylab as p

import visual as v
#from visual import *


#v.sphere()


n_part=500
n_config=20

density=0.01

box_size=(n_part/density)**(1./3.)
print "the box size is, ", box_size




already_ordered=0



if (already_ordered==0):


    tot_config=p.load("total_configuration.dat")
    tot_config=s.reshape(tot_config,(n_config,3*n_part))

    #r_gyr=s.zeros(n_config)

    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)]



    # #Now the arrays which will be containing the variance
    # var_x_arr=s.zeros(n_config)
    # var_y_arr=s.zeros(n_config)
    # var_z_arr=s.zeros(n_config)



    # 	len_int=(n_part*n_part-n_part)/2
    # 	print 'the number of interactions to be counted is, ', len_int

    r_0j=s.zeros((n_config,(n_part-1)))

    count_int=0

            #The following commented loop is what I would need to evaluate all the interaction
            #distances among the particles and it goes like n_part**2. I cannot run it as it is
            #since the array x_arr are 2D, but it gives an idea of the mechanism.
            #If I want to reconstruct the position of the particles wrt particle zero, then I
            #need only a loop on n-1 particles for each coord (and each stored configuration)


    #  	for i in xrange(0,(n_part-1)):
    #  		for j in xrange((i+1),n_part):
    # 			r_ij[count_int]=x_arr[i]-x_arr[j]
    # 			r_ij[count_int]=r_ij[count_int]-Len*n.round(r_ij[count_int]/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)

    # 			count_int=count_int+1

    # 	print 'r_ij is, ', r_ij




    Len=box_size


    #First I am going to calculate the correct variance of x

    for i in xrange(0,n_config):
        for j in xrange(1,n_part): #so, particle zero is now the reference particle
            r_0j[i,j-1]=x_arr[i,0]-x_arr[i,j]
            r_0j[i,j-1]=-(r_0j[i,j-1]-Len*n.round(r_0j[i,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, axis=1)


    # p.save("test_config_after_manipulation.dat",r_0j[600,:])

    #Now I am going to do the same along y and z


    for i in xrange(0,n_config):
        for j in xrange(1,n_part): #so, particle zero is now the reference particle
            r_0j[i,j-1]=y_arr[i,0]-y_arr[i,j]
            r_0j[i,j-1]=-(r_0j[i,j-1]-Len*n.round(r_0j[i,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)

    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, axis=1)






    for i in xrange(0,n_config):
        for j in xrange(1,n_part): #so, particle zero is now the reference particle
            r_0j[i,j-1]=z_arr[i,0]-z_arr[i,j]
            r_0j[i,j-1]=-(r_0j[i,j-1]-Len*n.round(r_0j[i,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)


            #print "r_0j is, ", r_0j


    z_arr[:,0]=0.
    z_arr[:,1:n_part]=r_0j


    #Now I save the new ordered configuration re-using the total_conf array

    for i in xrange(0,n_part):
        tot_config[:,(3*i)]=x_arr[:,i]
        tot_config[:,(3*i+1)]=y_arr[:,i]
        tot_config[:,(3*i+2)]=z_arr[:,i]
 
 




    p.save("tot_config_ordered.dat",tot_config)

    

            #var_x_arr[:]=s.var(r_0j, axis=1)
    # var_z_arr[:]=s.var(z_arr, axis=1)

    #	print 'z_arr is,', z_arr


    # r_gyr=s.sqrt(var_x_arr+var_y_arr+var_z_arr)



    # print 'the radius of gyration is, ', r_gyr




    # sphere(pos=(-5,0,0),radius=0.5,color=color.red)
    # #box(pos=(6,0,0),size=(0.2,4,4),color=color.green)

    # pos_list=s.arange(3)*1.3

    # sphere(pos=pos_list,radius=0.5,color=color.blue)


    #Now I select a certain configuration and plot it

    my_conf=12

    x_list=x_arr[my_conf,:]
    y_list=y_arr[my_conf,:]
    z_list=z_arr[my_conf,:]


    pos_list=s.zeros(3) # I initialize the array to contain the single particle's coordinates 

    my_rad=1.06/2.

    for i in xrange(0,n_part):
        pos_list[0]=x_list[i]
        pos_list[1]=y_list[i]
        pos_list[2]=z_list[i]
        particle=v.sphere(pos=pos_list,radius=my_rad,color=v.color.blue)

        

elif(already_ordered==1):
    
    tot_config=p.load("tot_config_ordered.dat")


    my_conf=0





#     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)]




    pos_list=s.zeros(3) # I initialize the array to contain the single particle's coordinates 

    my_rad=1.06/2.


    x_list=s.zeros(n_part)
    y_list=s.zeros(n_part)
    z_list=s.zeros(n_part)


    for i in xrange(0,n_part):


        x_list[i]=tot_config[my_conf,3*i]
        y_list[i]=tot_config[my_conf,(3*i+1)]
        z_list[i]=tot_config[my_conf,(3*i+2)]
        
        
   
#     v.scene = v.display(title="System Snapshot", width=box_size, height=box_size, x=0, y=0,
#                 range=box_size, center=(0.,0.,0.))

    particles=[v.sphere(pos=loc,radius=my_rad,color=v.color.blue)\
               for loc in zip(x_list,y_list,z_list)]    

    print 'particles[10].pos is, ', particles[10].pos






    for my_conf in xrange(1,n_config):
        
        v.rate(0.5)
        
        for i in xrange(0,n_part):
                #v.rate(5000)
                x_list[i]=tot_config[my_conf,3*i]
                y_list[i]=tot_config[my_conf,(3*i+1)]
                z_list[i]=tot_config[my_conf,(3*i+2)]
                particles[i].pos=(x_list[i],y_list[i],z_list[i])    

 

        



print 'So far so good'