分子動力学を使用したイオントラップシミュレーションとOpenGLによる可視化

イオン粒子を座標原点向きに力を加えトラップしつつ速度をゆっくりと落としていくイオントラップの分子動力学シミュレーション。
粒子どうしはクーロン力により反発しつつ、中心向きの力を受け摂動しながら中央に集まっていく。分子動力学1ステップ毎に粒子の速度を減速させる。最終的には球体状の形状で安定構造となる。

import random
import math
import itertools
import time
from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *

#from PIL import Image
#from PIL import ImageOps

#delta time
Dt = 0.01

step = 0
#Define xyz index
PX = 0
PY = 1
PZ = 2
VX = 3
VY = 4
VZ = 5
FX = 6
FY = 7
FZ = 8

FF = [FX, FY, FZ]
PP = [PX, PY, PZ]

#Coulomb's Force parameter
CF = 1000.0

#mouse move
mx = 0.0
my = 0.0

#lattice size
lsize = 1.0

#number of particles in a line
line_num = 8

#total particle num
PN = line_num * line_num * line_num

th = math.pi * 0.5
ph = math.pi * 0.5

#ready to 9 parameters for particle 
#(PX, PY, PZ, VX, VY, VZ, FX, FY, FZ)
xyz = [[0 for i in range(9)] for j in range(PN)] 

def find_pair():
  global PN
  for element in itertools.combinations(range(PN), 2):
    i  = element[0]
    j  = element[1]
    dx = xyz[i][PX] - xyz[j][PX]
    dy = xyz[i][PY] - xyz[j][PY]
    dz = xyz[i][PZ] - xyz[j][PZ]
    r  = math.sqrt(dx*dx + dy*dy + dz*dz)

    xyz[i][FX] = xyz[i][FX] + dx/(r*r*r)
    xyz[i][FY] = xyz[i][FY] + dy/(r*r*r)
    xyz[i][FZ] = xyz[i][FZ] + dz/(r*r*r)
    xyz[j][FX] = xyz[j][FX] - dx/(r*r*r)
    xyz[j][FY] = xyz[j][FY] - dy/(r*r*r)
    xyz[j][FZ] = xyz[j][FZ] - dz/(r*r*r)


def update():
  global step
  step += 1

  find_pair()

  pnum = 0
  while pnum < PN:
    # set force that collects particles
    xyz[pnum][FX] = xyz[pnum][FX] + CF*(lsize/2.0 - xyz[pnum][PX])
    xyz[pnum][FY] = xyz[pnum][FY] + CF*(lsize/2.0 - xyz[pnum][PY])
    xyz[pnum][FZ] = xyz[pnum][FZ] + CF*(lsize/2.0 - xyz[pnum][PZ])

    # if too large force is making , reset.
    CUT = 10000
    for f in FF:
      if xyz[pnum][f] > abs(CUT):
        xyz[pnum][f] = 0    

    xyz[pnum][VX] = xyz[pnum][VX] + 0.5*Dt*xyz[pnum][FX]
    xyz[pnum][VY] = xyz[pnum][VY] + 0.5*Dt*xyz[pnum][FY]
    xyz[pnum][VZ] = xyz[pnum][VZ] + 0.5*Dt*xyz[pnum][FZ]
    xyz[pnum][FX] = 0.0
    xyz[pnum][FY] = 0.0
    xyz[pnum][FZ] = 0.0
    pnum += 1
 

  print "- step =", step," -"
  pnum = 0
  while pnum < PN:
    #update position
    xyz[pnum][PX] = xyz[pnum][PX] + Dt*xyz[pnum][VX]
    xyz[pnum][PY] = xyz[pnum][PY] + Dt*xyz[pnum][VY]
    xyz[pnum][PZ] = xyz[pnum][PZ] + Dt*xyz[pnum][VZ]
    pnum += 1

def init_lattice():

  #PN = line_num^3 
  delt = lsize/line_num;
  xv_sum = 0
  yv_sum = 0
  zv_sum = 0

  pnum = 0
  k = 0
  while k < line_num:
    j = 0
    while j < line_num:
      i = 0
      while i < line_num:
        xyz[pnum][PX] = delt * i
        xyz[pnum][PY] = delt * j
        xyz[pnum][PZ] = delt * k
        pnum += 1
        i += 1
      j += 1
    k += 1

  pnum = 0
  while pnum < PN:
    xyz[pnum][VX] = random.uniform(-1,1) # 1 x-velocity
    xyz[pnum][VY] = random.uniform(-1,1) # 2 y-velocity
    xyz[pnum][VZ] = random.uniform(-1,1) # 3 z-velocity
    xv_sum += xyz[pnum][VX]
    yv_sum += xyz[pnum][VY]
    zv_sum += xyz[pnum][VZ]
    pnum += 1

  xv_sum = xv_sum / PN
  yv_sum = yv_sum / PN
  zv_sum = zv_sum / PN

  pnum = 0
  while pnum < PN:
    xyz[pnum][VX] = xyz[pnum][VX] - xv_sum
    xyz[pnum][VY] = xyz[pnum][VY] - yv_sum
    xyz[pnum][VZ] = xyz[pnum][VZ] - zv_sum
    pnum += 1


def draw():
    cx = cy = cz = lsize * 0.5
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glClearColor(0.0, 0.0, 0.0, 1.0)

    glLoadIdentity()
    glColor3f(1.0,1.0,0.0)

    gluLookAt(5*math.sin(th)*math.cos(ph)+cx, 5*math.cos(th)+cy, 5*math.sin(th)*math.sin(ph)+cz,\
            cx,  cy,  cz, -math.cos(th)*math.cos(ph),math.sin(th),-math.cos(th)*math.sin(ph) )

    glBegin(GL_LINE_LOOP)
    glVertex3d(0,0,0)
    glVertex3d(0,lsize,0)
    glVertex3d(lsize,lsize,0)
    glVertex3d(lsize,0,0)
    glEnd()
    glBegin(GL_LINE_LOOP)
    glVertex3d(0,0,lsize)
    glVertex3d(0,lsize,lsize)
    glVertex3d(lsize,lsize,lsize)
    glVertex3d(lsize,0,lsize)
    glEnd()
    glBegin(GL_LINES)
    glVertex3d(0,0,0)
    glVertex3d(0,0,lsize)
    glEnd()
    glBegin(GL_LINES)
    glVertex3d(0,lsize,0)
    glVertex3d(0,lsize,lsize)
    glEnd()
    glBegin(GL_LINES)
    glVertex3d(lsize,0,0)
    glVertex3d(lsize,0,lsize)
    glEnd()
    glBegin(GL_LINES)
    glVertex3d(lsize,lsize,0)
    glVertex3d(lsize,lsize,lsize)
    glEnd()

    pnum = 0
    while pnum < PN:
       glPointSize(3.0)
       glColor3f(0.3, 0.3, 1.0)

       glBegin(GL_POINTS)
       glVertex3d(xyz[pnum][PX], xyz[pnum][PY], xyz[pnum][PZ])
       glEnd()

       #Reduce particle velocity
       xyz[pnum][VX] *= 0.99 
       xyz[pnum][VY] *= 0.99 
       xyz[pnum][VZ] *= 0.99 
       pnum = pnum + 1

    update()
    glutSwapBuffers()


def init():
    glClearColor(0.7, 0.7, 0.7, 0.7)


def idle():
    glutPostRedisplay()


def reshape(w, h):
    glViewport(0, 0,w,h)
    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()
    gluPerspective(30.0, w/h, 1.0, 100.0)
    glMatrixMode (GL_MODELVIEW)


def motion(x, y):
  global ph, th, mx, my

  dltx = mx -x
  dlty = my -y

  #Invalid large motion
  if 10 < abs(dltx):
   dltx = 0 
  if 10 < abs(dlty):
   dlty = 0 
  
  ph = ph - 0.01*dltx
  th = th + 0.01*dlty
  glutPostRedisplay()
  glutSwapBuffers()
  mx = x
  my = y


if __name__ == "__main__":
  init_lattice()

  glutInitWindowPosition(50, 50);
  glutInitWindowSize(500, 500);
  glutInit(sys.argv)
  glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE )
  glutCreateWindow("pyOpenGL TEST")
  glutDisplayFunc(draw)
  glutReshapeFunc(reshape)
  glutMotionFunc(motion);
  init()
  glutIdleFunc(idle)
  glutMainLoop()

# END ----

Slowly down the speed while trapping the ion particles.
Particles are repulsed by the Coulomb force, and they are gathered in the center while being perturbed by the center-directed force. Molecular Dynamics Reduce the velocity of the particle every step.