Demo entry 6706016

Homework1 of MSE561

   

Submitted by Yanpeng on Jan 30, 2018 at 19:53
Language: Python 3. Code size: 2.5 kB.

#Yanpeng QI
#qiyp5
#Homework1 of 561

from matplotlib import pyplot as plt

r_cut_off = 7.50
r_0 = 7.00

''' E = E1(r) + E2(r)'''

'when distance is less than r_0, we get E1()'
def E1(r):
    sigma = 3.405
    e = 0.010323
    E_1 = 4 * e * ((sigma / r) ** 12-(sigma / r) ** 6)
    return E_1

'when distance is less than r_0, we get E2()'
def E2(r):
    A = -0.0068102128
    B = -0.0055640876
    global r_cut_off
    global r_0
    E_2 = A * (r - r_cut_off) ** 3 + B * (r - r_cut_off) ** 2
    return E_2

'initialize E in every condition'
E_r = []

E_r_1 = []
E_r_2 = []
E_r_3 = []
E_r_4 = []
E_r_5 = []
E_r_6 = []

E = []
R = []

r = 7.5

'''r change from 7.5 to 0 and for each stage, we get different value'''
while r >= r_0 and r <= r_cut_off:
    E_r_1.append([6 * E2(r), r])
    E.append(6 * E2(r))
    R.append(r)
    r -= 0.01
#print(min(E_1))    for test

while r >= r_cut_off / (3 ** 0.5) and r < r_0:
    E_r_2.append([6 * E1(r), r])
    E.append(6 * E1(r))
    R.append(r)
    r -= 0.01
#print(min(E_2))    for test

while r >= r_0 / (3 ** 0.5) and r < r_cut_off / (3 ** 0.5):
    E_r_3.append([6 * E1(r)+ 6 * E2((3 ** 0.5) * r), r])
    E.append(6 * E1(r)+ 6 * E2((3 ** 0.5) * r))
    R.append(r)
    r -= 0.01
#print(min(E_3))    for test

while r >= r_cut_off / 2 and r < r_0 / (3 ** 0.5):
    E_r_4.append([6 * E1(r) + 6 * E1(3 ** 0.5 * r), r])
    E.append(6 * E1(r) + 6 * E1(3 ** 0.5 * r))
    R.append(r)
    r -= 0.01
#print(min(E_4))    for test

while r >= r_0 / 2 and r < r_cut_off / 2:
    E_r_5.append([6 * E1(r) + 6 * E1(3 ** 0.5 * r) + 6 * E2(2 * r), r])
    E.append(6 * E1(r) + 6 * E1(3 ** 0.5 * r) + 6 * E2(2 * r))
    R.append(r)
    r -= 0.01
#print(min(E_5))    for test

''' Use 3 as the minmal value of r. 
I tried to use 0 as the minmal value, but I found it is not suitable to plot. 
When r=0, the potential energy will be so large that we can not find the trend.  '''
while r >= 3 and r < r_0 / 2:
    E_r_6.append([6 * E1(r) + 6 * E1(3 ** 0.5 * r) + 6 * E1(2 * r), r])
    E.append(6 * E1(r) + 6 * E1(3 ** 0.5 * r) + 6 * E1(2 * r))
    R.append(r)
    r -= 0.01
#print(min(E_6))    for test

#sum all data for different stages up and use print(E) for test
E_r = E_r_1 + E_r_2 + E_r_3 + E_r_4 + E_r_5 + E_r_6

#print(E_r) #for test
E_min = 0
r_min = 0
for y, x in E_r:
    if y < E_min:
        E_min = y
        r_min = x

print(E_min, r_min)
print(E)
print(R)
plt.plot(R, E)
plt.xlabel('Lattice Parameter')
plt.ylabel('Potential Energy')
plt.show()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).