Self Avoiding Random Walk in a 3d lattice using pivot algorithm in python

Question

I was working on a problem from past few days, It is related to creating a self avoiding random walks using pivot algorithm and then to implement another code which places spherical inclusions in a 3d lattice. I have written two codes one code generates coordinates of the spheres by using Self avoiding random walk algorithm and then another code uses the co-ordinates generated by the first program and then creates a 3d lattice in abaqus with spherical inclusions in it.

The result which will be generated after running the codes:

Now the zoomed in portion (marked in red boundary)

Problem: The spheres generated are merging with each other, How to avoid this, i ran out of ideas. I just need some direction or algorithm to work on.

The code which 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-r-r-r)),(2*r)+(random.random()*(10.0-r-r-r)),(2*r)+(random.random()*(10.0-r-r-r))]

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

The other code would not be required because there is no geometry computation in the second code. Any help or some direction of how to avoid the collision of spheres is very helpful, Thank you all

Answer 1

To accommodate non-intersecting spheres with turns of 45,135,225,315 (really, only 45 and 315 are issues), you just need to make your spheres a little bit smaller. Take 3 consecutive sphere-centers, with a 45-degree turn in the middle. In the plane containing the 3 points, that makes an isosceles triangle, with a 45-degree center angle:

Note that the bottom circles (spheres) overlap. To avoid this, you shrink the radius by multiplying by a factor of 0.76: