分子動力学による粒子のシミュレーション

粒子の動きを分子動力学で計算しmatplotlibで可視化してみる。
ただし各粒子間には相互作用は働かないシンプルな系として考える。
速度の計算はベルレ法(参考)。

各粒子はお互いが見えていないので、初速のまま真っ直ぐに飛び続ける。境界を越えたら反対側へ。

#-*- coding:utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
import random

a  = 10
Dt = 0.1

#Define xy index
PX = 0
PY = 1
VX = 2
VY = 3
FX = 4
FY = 5

def integ(xy):
  #verlet integration
  particle_num = len(xy)
  pnum = 0
  while pnum < particle_num:
    #update velocity
    xy[pnum][VX] = xy[pnum][VX] + 0.5*Dt*xy[pnum][FX]
    xy[pnum][VY] = xy[pnum][VY] + 0.5*Dt*xy[pnum][FY]
    #update position
    xy[pnum][PX] = xy[pnum][PX] + Dt*xy[pnum][VX]
    xy[pnum][PY] = xy[pnum][PY] + Dt*xy[pnum][VY]
    #boundary condition
    if xy[pnum][PX] > a:
       xy[pnum][PX] -= a
    if xy[pnum][PY] > a:
       xy[pnum][PY] -= a
    if xy[pnum][PX] < 0:
       xy[pnum][PX] += a
    if xy[pnum][PY] < 0:
       xy[pnum][PY] += a
    pnum += 1


def update(i, fig_title, xy, total_step):
    if i != 0:
        plt.cla()
    ax.set_xlim(0,10)
    ax.set_ylim(0,10)
    integ(xy)
    nnum = 0
    while nnum < len(xy):
       plt.plot(xy[nnum][PX], xy[nnum][PY], "o")
       nnum = nnum + 1
    plt.title(fig_title + ' step = ' + str(i) + '/' + str(total_step))


def init_set(xy):
  pnum = 0
  xv_sum = 0
  yv_sum = 0
  particle_num = len(xy)
  while pnum < particle_num:
    xy[pnum][PX] = 0  # x-position
    xy[pnum][PY] = 0  # y-position
    xy[pnum][VX] = random.uniform(-1,1) # x-velocity
    xy[pnum][VY] = random.uniform(-1,1) # y-velocity
    xv_sum += xy[pnum][VX]
    yv_sum += xy[pnum][VY]
    xy[pnum][FX] = 0   # x-force
    xy[pnum][FY] = 0   # y-force
    pnum += 1
  #Total momentum to zero
  xv_sum = xv_sum / particle_num
  yv_sum = yv_sum / particle_num
  pnum = 0
  while pnum < particle_num:
    xy[pnum][VX] = xy[pnum][VX] - xv_sum
    xy[pnum][VY] = xy[pnum][VY] - yv_sum
    pnum += 1


if __name__ == "__main__":
    fig = plt.figure(figsize = (10, 6))
    fig.set_size_inches(10,10)
    ax = fig.add_subplot(111)
    perticle_num = 100
    xy = [[0 for i in range(6)] for j in range(perticle_num)]
    init_set(xy)
    total_step =400
    ani = anm.FuncAnimation(fig, \
                        update, fargs = ('Molecular Dynamics ', xy, total_step), \
                        interval = 10, frames =  total_step)
    plt.show()
    #ani.save("Sample_MD.gif", writer = 'imagemagick', dpi=50)


  

実行の結果は下記。