#はじめに

[Pythonによる科学・技術計算]初歩的な分子動力学シミュレーション, 2次元系, 熱・統計物理学

で作成した分子動力学(MD)シミュレーションの応用の一つとして，ここでは原子団に原子(弾丸)を撃ち込んでみます。
計算サイズは小さいですが，破壊現象などを原子レベルでシミュレートする際の基礎になると思います。

#内容

2次元原子系のMDシミュレーションを行い，平面に並んだ粒子に高速粒子(弾丸)を打ち込むアニメーションを作成します。
シミュレーションではレナード-ジョーンズポテンシャルによる2体相互作用を考慮します。

打ち込む弾丸の速度として以下の3種類を考えます。速度によって衝突過程がどのように変わるのか調べます。
(1) ~140 (高速)
(2) ~14 (低速)

#どんな人向け？

下記の方たちにとって少なからず有益だと思います。

● MDシミュレーションの応用の一つを見てみたい方
● 原子衝突によって木っ端みじんになる原子団に興味がある方

#計算コード

時間間隔: delta_t = 0.002
原子数: 49原子

ターゲットとなる原子段の初期速度をゼロとします。
左下斜めからそれらにめがけて弾丸を撃ち込みます。

なお，本シミュレーションでは力のカットオフや周期境界条件を使いません。

弾丸の速度と初期位置を変更する場合は
x[0]とy[0]，そしてvx[0]とvy[0]の値を変更します。

なお，このコードは，基本的には[1]にあるコードで初期条件を変えたものになります。
各関数については[1]をご参照ください。

%%time
"""
レナード-ジョーンズ系のMDシミュレーション
弾丸衝突: 2D
Sept. 2017
"""
%matplotlib nbagg

import numpy as np
import random
import matplotlib.pyplot as plt
import matplotlib.animation as animation



fig = plt.figure(figsize = (8, 6))
ims = []



lat_x = 1 # 格子定数: x
lat_y = 1
Lx = 7 # 格子点: x
Ly = 7 
N_atom = Lx*Ly

delta_t = 0.002

T_initial = 0 

# 粒子生成

Max_step =500

x = np.zeros([N_atom],float)
y = np.zeros([N_atom],float)

vx = np.zeros([N_atom],float) #原子団の初速を0
vy = np.zeros([N_atom],float)

fx = np.zeros([N_atom,2],float)#
fy = np.zeros([N_atom,2],float)   




def initialposvel():
    i = -1
    for ix in range(Lx):
        for iy in range(Ly):
            i += 1
            x[i] = ix*lat_x
            y[i] = iy*lat_y
    
    # 弾丸作成
    x[0] = -3*lat_x
    y[0] = -3*lat_y
    vx[0] = 100
    vy[0] = 100
    

    
def Forces(t,w,PE,PEorW):

    PotEnergy = 0

    
    fx[:,t] = 0
    fy[:,t] = 0
    for i in range(N_atom-1):
        for j in range(i+1, N_atom):
            dx = x[i] - x[j]
            dy = y[i] - y[j]
            r2 = dx**2+dy**2
                
         
            if (r2 == 0):
                r2 = 0.0001
            invr2 = 1/r2
                    
            wij = 48*(invr2**3-0.5)*(invr2**3)
            fijx = wij*invr2*dx
            fijy = wij*invr2*dy
                    
            fx[i,t] += fijx
            fy[i,t] += fijy
            fx[j,t] -= fijx
            fy[j,t] -= fijy
                    
            PotEnergy += 4*(invr2**3)*((invr2**3)-1)
            w += wij
    if (PEorW == 1) :
        return PotEnergy
    else :
        return w
         
# メイン: 時間発展
def time_evolution():
    
 
    
    t1 = 0
    PE = 0.0
    
 #   hover2 = delta_t/2.0
    
    # 初期化

    initialposvel()

    
    
    KE = (np.sum(vx**2)+np.sum(vy**2))/2
    w=0
    PE = Forces(t1, w, PE, 1)
    time = 1
        
    n_step=0
    

    
    
    while n_step < Max_step:
        if n_step % 2 ==0:
            im1=plt.scatter(x[1:], y[1:], s=200, c="pink", alpha=0.9, linewidths="1.5",edgecolors="black")
            im2=plt.scatter(x[0], y[0], s=200, c="blue", alpha=0.9, linewidths="1.5",edgecolors="black")
    
        ims.append([im1,im2])
       
        

        n_step += 1
        if n_step % 50==0: 
            print ("step=",n_step)
            
        for i in range(N_atom):  
            PE = Forces(t1, w, PE, 1)    
            # 速度ベルレー法による更新
            x[i] +=delta_t*(vx[i]+delta_t*fx[i,t1]/2)
            y[i] +=delta_t*(vy[i]+delta_t*fy[i,t1]/2)


                
        PE = 0.
        t2=1
        PE = Forces(t2, w, PE, 1)


        
        
        # 速度ベルレー法による更新
        vx[:] += delta_t*(fx[:,t1]+fx[:,t2])/2
        vy[:] += delta_t*(fy[:,t1]+fy[:,t2])/2


time_evolution()

plt.axis([-4*lat_x, (Lx+2)*lat_x, -4*lat_y,(Ly+2)*lat_y])
ani = animation.ArtistAnimation(fig, ims, interval =1)
plt.show() 
ani.save("output.gif", writer="imagemagick")

#結果(1):高速: 速度 140

この速度だと原子団が切り裂かれ原子団はバラバラになります。再起不能です。

#結果(2):低速: 速度 14

低速の場合はどうでしょうか...。

低速衝突の場合は，部分的な破壊が生じています。
わずかな擾乱でも原子団からすれば原子間結合がブチブチと切られてしまい，不安定な構造となってしまいます。

#参考文献

[1][Pythonによる科学・技術計算]初歩的な分子動力学シミュレーション, 2次元系, 熱・統計物理学