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.

Delete this entry (admin only).