OSはubunts 18
%pip3 install PyOpenGL PyOpenGL_accelerate
分子動力学を使用したイオントラップシミュレーションとOpenGLによる可視化(multiprocessingによる処理並列化版)
分子動力学を使用したイオントラップシミュレーションとOpenGLによる可視化
３D格子状に配置した粒子をopenGLで可視化、マウスで視点移動
ブラウン運動のシミュレーション可視化
以前ここで試した単純な粒子の中に大きいサイズの粒子を放り込む。周囲から無数の粒子がぶつかりランダムに動くブラウン運動のような動きをみることができる。
■結果
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 |
import random import math from OpenGL.GL import * from OpenGL.GLU import * from OpenGL.GLUT import * from PIL import Image from PIL import ImageOps a = 2 Dt = 0.01 CUT = 5 #Define xy index PX = 0 PY = 1 VX = 2 VY = 3 FX = 4 FY = 5 pot = [0,0] perticle_num = 1000 xy = [[0 for i in range(6)] for j in range(perticle_num)] BP = [0 for i in range(6)] Radi = 0.1 step = 1 def capture(): global step pad_step = '{0:04d}'.format(step) print pad_step width = glutGet(GLUT_WINDOW_WIDTH) height = glutGet(GLUT_WINDOW_HEIGHT) glReadBuffer(GL_FRONT) glPixelStorei(GL_UNPACK_ALIGNMENT, 1) data = glReadPixels(0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE) image = Image.fromstring("RGBA", (width, height), data) image = ImageOps.flip(image) image.save("test"+pad_step+".png") def find_pair(): global Radi particle_num = len(xy) pnum = 0 w_ratio = 0.1 while pnum < particle_num: dx = xy[pnum][PX] - BP[PX] dy = xy[pnum][PY] - BP[PY] if dx > 0.5*a: dx = dx + a if dy > 0.5*a: dy = dy + a rr = math.sqrt(dx*dx + dy*dy) if rr < Radi: v2 = math.sqrt(xy[pnum][VX]**2 + xy[pnum][VY]**2) cos = (xy[pnum][VX]*dx + xy[pnum][VY]*dy)/(v2*rr) sin = math.sqrt(1 - cos**2) tmpx= xy[pnum][VX]*(cos**2 - sin**2) - xy[pnum][VY]*(2*cos*sin) tmpy= xy[pnum][VX]*(2*cos*sin) + xy[pnum][VY]*(cos**2 - sin**2) xy[pnum][VX] = abs(tmpx) * dx / (abs(dx)) xy[pnum][VY] = abs(tmpy) * dy / (abs(dy)) BP[VX] = BP[VX] - w_ratio * abs(tmpx) * dx / (abs(dx)) BP[VY] = BP[VY] - w_ratio * abs(tmpy) * dy / (abs(dy)) pnum += 1 def integ(): find_pair() particle_num = len(xy) BP[PX] = BP[PX] + Dt*BP[VX] BP[PY] = BP[PY] + Dt*BP[VY] #boundary condition if BP[PX] > a: BP[PX] -= a if BP[PY] > a: BP[PY] -= a if BP[PX] < 0: BP[PX] += a if BP[PY] < 0: BP[PY] += a pnum = 0 while pnum < particle_num: #update position xy[pnum][PX] = xy[pnum][PX] + Dt*xy[pnum][VX] xy[pnum][PY] = xy[pnum][PY] + Dt*xy[pnum][VY] 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 calc_ena(): ena = 0.0 particle_num = len(xy) for i in range(0,particle): vx2 = xy[i][VX]*xy[i][VX] vy2 = xy[i][VY]*xy[i][VY] ena = ena + vx2 + vy2 ena = 0.5*ena return ena def update(): global Radi PI = 3.14159263 global step step = step + 1 # If you want to get image... #capture() integ() #print step nnum = 0 while nnum < len(xy): glPointSize(3.0) glColor3f(0.3, 0.3, 1.0) glBegin(GL_POINTS) glVertex2d(xy[nnum][PX], xy[nnum][PY]) glEnd() nnum = nnum + 1 N_t = 30 glColor3f(1.0, 0.0, 0.0) glBegin(GL_POLYGON) for j in range(0,N_t): glVertex2d(Radi*math.cos(2.0*PI*j/N_t)+BP[PX], Radi*math.sin(2.0*PI*j/N_t)+BP[PY]) glEnd() def init_set(): pnum = 0 xv_sum = 0 yv_sum = 0 particle_num = len(xy) while pnum < particle_num: xy[pnum][PX] = 0 # 0 x-position xy[pnum][PY] = 0 # 1 y-position xy[pnum][VX] = random.uniform(-1,1) # 2 x-velocity xy[pnum][VY] = random.uniform(-1,1) # 3 y-velocity xv_sum += xy[pnum][VX] yv_sum += xy[pnum][VY] pnum += 1 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 BP[PX] = a/2 BP[PY] = a/2 BP[VX] = random.uniform(-0.5,0.5) BP[VY] = random.uniform(-0.5,0.5) def draw(): glClearColor(0.0, 0.0, 0.0, 0.0) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) glLoadIdentity() glColor3f(1.0,1.0,0.0) gluLookAt(a/2, a/2, 4.0 ,a/2 ,a/2 ,0 , 0, 1.0, 0.0) glBegin(GL_LINE_LOOP) glVertex2d(0,0) glVertex2d(0,a) glVertex2d(a,a) glVertex2d(a,0) glEnd() update() glutSwapBuffers() def resize(w, h): glViewport(0, 0, w, h) glLoadIdentity() glOrtho(0, a, 0, a, -a, a) 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) if __name__ == "__main__": init_set() glutInitWindowPosition(50, 50); glutInitWindowSize(340, 340); glutInit(sys.argv) glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE ) glutCreateWindow("pyOpenGL TEST") glutDisplayFunc(draw) glutReshapeFunc(reshape) init() glutIdleFunc(idle) glutMainLoop() |
brownian motion simulation visualization by python and openGL
pyOpenGLを使って分子動力学を可視化
以前ここで試したmatplotlibを使った簡単な分子動力学計算の可視化処理ですが、描画が遅いのでOpenGLを使って書き直してみた。python用にはPyOpenGLのパッケージが用意されている。結果だがとても速い！
■結果をgifアニメに落としたものが↓。
gifアニメだと実際の描画より遅いので感動が伝わりにくい。。。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 |
import random from OpenGL.GL import * from OpenGL.GLU import * from OpenGL.GLUT import * from PIL import Image from PIL import ImageOps a = 10 Dt = 0.1 #Define xy index PX = 0 PY = 1 VX = 2 VY = 3 FX = 4 FY = 5 perticle_num = 1000 xy = [[0 for i in range(6)] for j in range(perticle_num)] step = 1 def capture(): global step step = step + 1 pad_step = '{0:04d}'.format(step) print pad_step width = glutGet(GLUT_WINDOW_WIDTH) height = glutGet(GLUT_WINDOW_HEIGHT) glReadBuffer(GL_FRONT) glPixelStorei(GL_UNPACK_ALIGNMENT, 1) data = glReadPixels(0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE) image = Image.fromstring("RGBA", (width, height), data) image = ImageOps.flip(image) image.save("test"+pad_step+".png") def integ(): #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(): # If you want to get image... # capture() integ() nnum = 0 while nnum < len(xy): glPointSize(3.0) glColor3f(0.3, 0.3, 1.0) glBegin(GL_POINTS) glVertex2d(xy[nnum][PX], xy[nnum][PY]) glEnd() nnum = nnum + 1 def init_set(): pnum = 0 xv_sum = 0 yv_sum = 0 particle_num = len(xy) while pnum < particle_num: xy[pnum][PX] = 0 # 0 x-position xy[pnum][PY] = 0 # 1 y-position xy[pnum][VX] = random.uniform(-1,1) # 2 x-velocity xy[pnum][VY] = random.uniform(-1,1) # 3 y-velocity xv_sum += xy[pnum][VX] yv_sum += xy[pnum][VY] xy[pnum][FX] = 0 # 0 x-force xy[pnum][FY] = 0 # 1 y-force pnum += 1 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 #exit() def draw(): glClearColor(0.0, 0.0, 0.0, 0.0) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) glLoadIdentity() glColor3f(1.0,1.0,0.0) gluLookAt(a/2, a/2, 20.0 ,a/2 ,a/2 ,0 , 0, 1.0, 0.0) glBegin(GL_LINE_LOOP) glVertex2d(0,0) glVertex2d(0,a) glVertex2d(a,a) glVertex2d(a,0) glEnd() 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) if __name__ == "__main__": init_set() glutInitWindowPosition(50, 50); glutInitWindowSize(340, 340); glutInit(sys.argv) glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE ) glutCreateWindow("pyOpenGL TEST") glutDisplayFunc(draw) glutReshapeFunc(reshape) init() glutIdleFunc(idle) glutMainLoop() |
pythonでopenGLを試す
先日matplotlibを使った簡単な分子動力学の可視化を試したが非常に描画が遅い。ググってみるとどうもmatplotlibはプロッティングの遅さに定評があるらしい。
そんなわけで違う方法での可視化を試す。
試してみるのはopenGL。python用にはpyOpenGLというパッケージが用意されている。当方の環境centos7では下記コマンドで導入が出来た。
%easy_install PyOpenGL
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
from OpenGL.GL import * from OpenGL.GLUT import * def draw(): glClearColor(1.0, 0.0, 0.0, 0.0) glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT) glFlush() glutSwapBuffers() if __name__ == "__main__": glutInit(sys.argv) glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH) glutInitWindowSize(320, 240) glutCreateWindow("openGL test") glutDisplayFunc(draw) glutMainLoop() |
■実行結果
