Self Avoiding Random Walk in 3d grid using pivot algorithm in python

Question

I have been working on a problem in the last few days. This is due to the creation of self-avoidance of random walks using a summary algorithm, and then to implement another code that puts spherical inclusions in a 3d lattice. I wrote two codes, one code generates the coordinates of the spheres using the Self algorithm, which avoids random walks, and then another code uses the coordinates generated by the first program, and then creates a 3d lattice in abaqus with spherical inclusions in it.

The result that will be generated after the codes are run:

Now the enlarged part (marked with a red border)

Problem: Generated spheres merge with each other. To avoid this, I ran out of ideas. I just need some direction or algorithm of work .

The code I wrote is as follows: code 1: generates the coordinates of the spheres

# Functions of the code: 1) This code generates the coordinates for the spherical fillers using self avoiding random walk which was implemented by pivot algorithm # Algorithm of the code: 1)Prepare an initial configuration of a N steps walk on lattice(equivalent N monomer chain) # 2)Randomly pick a site along the chain as pivot site # 3)Randomly pick a side(right to the pivot site or left to it), the chain on this side is used for the next step. # 4)Randomly apply a rotate operation on the part of the chain we choose at the above step. # 5)After the rotation, check the overlap between the rotated part of the chain and the rest part of the chain. # 6)Accept the new configuration if there is no overlap and restart from 2th step. # 7)Reject the configuration and repeat from 2th step if there are overlaps. ################################################################################################################################################################ # Modules to be imported are below import numpy as np import timeit from scipy.spatial.distance import cdist import math import random import sys import ast # define a dot product function used for the rotate operation def v_dot(a):return lambda b: np.dot(a,b) def distance(x, y): #The Euclidean Distance to Spread the spheres initiallly if len(x) != len(y): raise ValueError, "vectors must be same length" sum = 0 for i in range(len(x)): sum += (x[i]-y[i])**2 return math.sqrt(sum) x=1/math.sqrt(2) class lattice_SAW: # This is the class which creates the self avoiding random walk coordinates def __init__(self,N,l0): self.N = N #No of spheres self.l0 = l0 #distance between spheres # initial configuration. Usually we just use a straight chain as inital configuration self.init_state = np.dstack((np.arange(N),np.zeros(N),np.zeros(N)))[0] #initially set all the coordinates to zeros self.state = self.init_state.copy() # define a rotation matrix # 9 possible rotations: 3 axes * 3 possible rotate angles(45,90,135,180,225,270,315) self.rotate_matrix = np.array([[[1,0,0],[0,0,-1],[0,1,0]],[[1,0,0],[0,-1,0],[0,0,-1]] ,[[1,0,0],[0,0,1],[0,-1,0]],[[0,0,1],[0,1,0],[-1,0,0]] ,[[-1,0,0],[0,1,0],[0,0,-1]],[[0,0,-1],[0,1,0],[-1,0,0]] ,[[0,-1,0],[1,0,0],[0,0,1]],[[-1,0,0],[0,-1,0],[0,0,1]] ,[[0,1,0],[-1,0,0],[0,0,1]],[[x,-x,0],[x,x,0],[0,0,1]] ,[[1,0,0],[0,x,-x],[0,x,x]] ,[[x,0,x],[0,1,0],[-x,0,x]],[[-x,-x,0],[x,-x,0],[0,0,1]] ,[[1,0,0],[0,-x,-x],[0,x,-x]],[[-x,0,x],[0,1,0],[-x,0,-x]] ,[[-x,x,0],[-x,-x,0],[0,0,1]],[[1,0,0],[0,-x,x],[0,-x,-x]] ,[[-x,0,-x],[0,1,0],[x,0,-x]],[[x,x,0],[-x,x,0],[0,0,1]] ,[[1,0,0],[0,x,x],[0,-x,x]],[[x,0,-x],[0,1,0],[x,0,x]]]) # define pivot algorithm process where t is the number of successful steps def walk(self,t): #this definitions helps to start walking in 3d acpt = 0 # while loop until the number of successful step up to t while acpt <= t: pick_pivot = np.random.randint(1,self.N-1) # pick a pivot site pick_side = np.random.choice([-1,1]) # pick a side if pick_side == 1: old_chain = self.state[0:pick_pivot+1] temp_chain = self.state[pick_pivot+1:] else: old_chain = self.state[pick_pivot:] # variable to store the coordinates of the old chain temp_chain = self.state[0:pick_pivot]# for the temp chain # pick a symmetry operator symtry_oprtr = self.rotate_matrix[np.random.randint(len(self.rotate_matrix))] # new chain after symmetry operator new_chain = np.apply_along_axis(v_dot(symtry_oprtr),1,temp_chain - self.state[pick_pivot]) + self.state[pick_pivot] # use cdist function of scipy package to calculate the pair-pair distance between old_chain and new_chain overlap = cdist(new_chain,old_chain) #compare the old chain and the new chain to check if the new chain overlaps on any previous postion overlap = overlap.flatten() # just to combine the coordinates in a list which will help to check the overlap # determinte whether the new state is accepted or rejected if len(np.nonzero(overlap)[0]) != len(overlap): continue else: if pick_side == 1: self.state = np.concatenate((old_chain,new_chain),axis=0) elif pick_side == -1: self.state = np.concatenate((new_chain,old_chain),axis=0) acpt += 1 # place the center of mass of the chain on the origin, so then we can translate these coordinates to the initial spread spheres self.state = self.l0*(self.state - np.int_(np.mean(self.state,axis=0))) #now write the coordinates of the spheres in a file sys.stdout = open("C:\Users\haris\Desktop\CAME_MINE\HIWI\Codes\Temp\Sphere_Positions.txt", "w") n = 30 #text_file.write("Number of Spheres: " + str(n) + "\n") #Radius of one sphere: is calculated based on the volume fraction here it 2 percent r = (3*0.40)/(4*math.pi*n) #print r r = r**(1.0/3.0) #print "Sphere Radius is " + str(r) sphereList = [] # list to maintain track of the Spheres using their positions sphereInstancesList = [] # to maintain track of the instances which will be used to cut and or translate the spheres sphere_pos=[] #Create n instances of the sphere #After creating n instances of the sphere we trie to form cluster around each instance of the sphere for i in range(1, n+1): InstanceName = 'Sphere_' + str(i) # creating a name for the instance of the sphere #print InstanceName #text_file.write(InstanceName) #Maximum tries to distribute sphere maxTries = 10000000 while len(sphereList) < i: maxTries -= 1 if maxTries < 1: print "Maximum Distribution tries exceded. Error! Restart the Script!" break; # Make sure Spheres dont cut cube sides: this will place the spheres well inside the cube so that they does'nt touch the sides of the cube # This is done to avoid the periodic boundary condition: later in the next versions it will be used vecPosition = [(2*r)+(random.random()*(10.0-rrr)),(2*r)+(random.random()*(10.0-rrr)),(2*r)+(random.random()*(10.0-rrr))] sphere_pos.append(vecPosition) for pos in sphereList: if distance(pos, vecPosition) < 2*r: # checking whether the spheres collide or not break else: sphereList.append(vecPosition) print vecPosition #text_file.write(str(vecPosition) + "\n") cluster_Size=[10,12,14,16,18,20,22,24,26,28,30] #list to give the random number of spheres which forms a cluster for i in range(1,n+1): Number_of_Spheres_Clustered=random.choice(cluster_Size) #selecting randomly from the list cluster_Size radius_sphr=2*r #Distance between centers of the spheres pivot_steps=1000 # for walking the max number of steps chain = lattice_SAW(Number_of_Spheres_Clustered,radius_sphr) #initializing the object chain.walk(pivot_steps) # calling the walk function to walk in the 3d space co_ordinates=chain.state # copying the coordinates into a new variable for processing and converting the coordinates into lists for c in range(len(co_ordinates)): temp_cds=co_ordinates[c] temp_cds=list(temp_cds) for k in range(len(temp_cds)): temp_cds[k]=temp_cds[k]+sphere_pos[i-1][k] #text_file.write(str(temp_cds) + "\n") print temp_cds #sys.stdout redirected into the file which stores the coordinates as lists sys.stdout.flush() f2=open("C:\Users\haris\Desktop\CAME_MINE\HIWI\Codes\Temp\Sphere_Positions.txt", "r") remove_check=[] for line in f2: temp_check=ast.literal_eval(line) if (temp_check[0]>10 or temp_check[0]<-r or temp_check[1]>10 or temp_check[1]<-r or temp_check[2]>10 or temp_check[2]<-r): remove_check.append(str(temp_check)) f2.close() flag=0 f2=open("C:\Users\haris\Desktop\CAME_MINE\HIWI\Codes\Temp\Sphere_Positions.txt", "r") f3=open("C:\Users\haris\Desktop\CAME_MINE\HIWI\Codes\Temp\Sphere_Positions_corrected.txt", "w") for line in f2: line=line.strip() if any(line in s for s in remove_check): flag=flag+1 else: f3.write(line+'\n') f3.close() f2.close() 

Another code is not required because the second code does not have geometric calculations. Any help or some direction on how to avoid a collision of spheres is very useful , thanks to everyone

Answer 1

To place disjoint spheres with 45,135,225,315 orbits (indeed, only 45 and 315 are problems), you just need to make your spheres a little smaller. Take 3 consecutive spherical centers with a 45-degree rotation in the middle. In a plane containing 3 points, forms an isosceles triangle with a center of 45 degrees:

Note that the lower circles (spheres) overlap. To avoid this, you reduce the radius by multiplying the coefficient by 0.76: