# 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 = [ * 20] * 400
E = [0.0] * 400
F = [ * 3] * 400
NOFNEIGHBORS =  * 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
F[i]=0
F[i]=0
for j in range(0, 400):
#F[i]:sigma 11 of atom i
F[i] += stress(i, j, 1, 1)
#F[i]:sigma 12 of atom i
F[i] += stress(i, j, 1, 2)
#F[i]:sigma 22 of atom i
F[i] += 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]) + '  ' + str(F[i]) + '  ' + str(F[i]) + '\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.