# Demo entry 6716382

**HOMEWORK2 OF MSE 561**

Submitted by **Yanpeng**
on Feb 19, 2018 at 01:09

Language: Python. Code size: 4.8 kB.

#Homework2 MSE561 #Yanpeng QI from matplotlib import pyplot as plt #import this module for ploting def position(i, j): '''this function is used to get atom's position, x and y''' m = i % 20 #get m as postion x n = int(i) // 20 #get n as postion y if j == 1: #use j, we can judge use x position or y position position = 3.822 * m else: position = 3.822 * n #output the requested position return position def distance(n, m): '''this function is used to get the distance between atom n and atom m''' d = ((position(n, 1) - position(m, 1))**2 + (position(n, 2) - position(m, 2) )**2)**0.5 return d def potential(r): '''this function is according to the homework1 and get the potential between atoms according to d''' r_0 = 7.0 r_cut = 7.5 sigma = 3.405 epsi = 0.0103023 A = -6.8102128*10**(-3) B = -5.5640876*10**(-3) if r < r_0 and r > 0: #the potentail function is different between 0-7 and 7-7.5 u = 4 * epsi * ((sigma / r)**12 - (sigma / r)**6) elif r >= r_0 and r < r_cut: u = A * (r - r_cut)**3 + B * (r - r_cut)**2 else: u = 0 return u def d_potential(r): '''this function get us the derivate of potential with r, we will use it to get the stress in each direction''' r_0 = 7.0 r_cut = 7.5 sigma = 3.405 epsi = 0.0103023 A = -6.8102128*10**(-3) B = -5.5640876*10**(-3) if r < r_0 and r > 0: #this derivate of potentail is according to the function potentail above and there is still a difference in range d_u = 4 * epsi *(sigma**12 * (-12) * (1/r)**(13) + (sigma**6 * 6 * (1 / r)**(7))) elif r >= r_0 and r < r_cut: d_u = 3 * A * (r - r_cut)**2 + 2 * B * (r - r_cut) else: d_u = 0 return d_u def stress(i, j, alpha, beta): '''this function is used to get the stress in each direction between each atoms''' area = ((19 * 3.822) ** 2) / 400 #the whole area r_alpha = position(i, alpha) - position(j, alpha) #dstance in the direction alpha r_beta = position(i, beta) - position(j, beta) #distance in the direction beta r = distance(i, j) #distance between atom i and atom j d_p = d_potential(r) if i != j: #only when i and j are different atoms, we can use the formla stress = (1 / area) * d_p * r_alpha * r_beta / r else: #when i and j are the same atom, stress is equal to 0 stress = 0 return stress def main(): #open the txt file that we want to import data f = open('/Users/apple/Desktop/561/2/assignmentOut1.txt', 'w') r_cut = 7.5 #set space for these lists NEIGHBORS = [[0] * 20] * 400 E = [0.0] * 400 F = [[0] * 3] * 400 NOFNEIGHBORS = [0] * 400 for i in range(0, 400): #initilize the i_st element to zero and use the next for to add value to them NOFNEIGHBORS[i]=0 E[i]=0.0 F[i][0]=0 F[i][1]=0 F[i][2]=0 for j in range(0, 400): #F[i][0]:sigma 11 of atom i F[i][0] += stress(i, j, 1, 1) #F[i][1]:sigma 12 of atom i F[i][1] += stress(i, j, 1, 2) #F[i][2]:sigma 22 of atom i F[i][2] += stress(i, j, 2, 2) # r is the distance between i and j to judge whether it is neighbour r = distance(i, j) if r < r_cut and r > 0: #if r in the range(0, 7), create a neighbor, number of neighbours add one NOFNEIGHBORS[i] += 1 #recrod this number j into NEIGHBOS[i] NEIGHBORS[i][NOFNEIGHBORS[i]] = j E[i] += potential(r) f.write(str(position(i , 1)) + ' ' + str(position(i , 2)) + ' ' + str(E[i]) + ' ' + str(F[i][0]) + ' ' + str(F[i][1]) + ' ' + str(F[i][2]) + '\n') f.close() #plot RDF, this first neighbour should be 3.822, we choose 3.800 as the initial point r = 3.800 #create two list to get r and g(r) for RDF plot G_r = [] r_list = [] #use this while loop to check r from 3.800 to 15 with step 0.01 while r <= 15: N_r = 0 for i in range(0, 400): for j in range(0, 400): if distance(i, j) < r + 0.01 and distance(i, j) >= r: #use N_r to represent the number of atoms between r and 𝑟 + ∆𝑟 N_r += 1 pi = 3.1415926 p = 27.730012830676 #use this formula to get the value of g(r) G_r.append(1 / p * (N_r / (2 * pi * r * 0.01 * 400))) r_list.append(r) r += 0.01 #plot g(r)-r plt.plot(r_list, G_r) plt.xlabel('r 10^-10^-10 m') plt.ylabel('g(r) count/area') plt.show() if __name__ == '__main__': main()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.