double rx[natom], ry[natom], rz[natom]; 原子の座標速度力ポテンシャルエ double vx[natom], vy[natom], vz[natom]; ネルギーを受ける配列を準備 double fx[natom], fy[natom], fz[natom

Size: px

Start display at page:

Download "double rx[natom], ry[natom], rz[natom]; 原子の座標速度力ポテンシャルエ double vx[natom], vy[natom], vz[natom]; ネルギーを受ける配列を準備 double fx[natom], fy[natom], fz[natom"

ともなりおおばま
5 years ago
Views:

GLUI による MD の GUI 化前提条件 :GLUI のプログラミング環境が整っていること 3 原子の MD コード ( 下図 ) viewer ウィンドウ内のマウス左クリックで MD 開始右クリックで MD 停止 control パネルは solid/wireframe を切り替えるチェックボタン球の滑らかさと半径を決める窓 ( スピナー ) オブジェクトを回転移動

1 GLUI による MD の GUI 化前提条件 :GLUI のプログラミング環境が整っていること 3 原子の MD コード ( 下図 ) viewer ウィンドウ内のマウス左クリックで MD 開始右クリックで MD 停止 control パネルは solid/wireframe を切り替えるチェックボタン球の滑らかさと半径を決める窓 ( スピナー ) オブジェクトを回転移動拡大縮小させるコントローラ MD ステップの表示窓 ( テキストボックス ) ソフトを修了させるボタンからなる以下にソースコードを示す #include <string.h> #include <iostream> #include <GL/glui.h> #include<fstream> #include<math.h> #define NATOM 3 色々な関数を呼ぶための include 文原子の総数を NATOM と define

2 double rx[natom], ry[natom], rz[natom]; 原子の座標速度力ポテンシャルエ double vx[natom], vy[natom], vz[natom]; ネルギーを受ける配列を準備 double fx[natom], fy[natom], fz[natom]; 系全体のポテンシャル運動エネルギー double epot[natom], epotsum, ekinsum; の変数も用意 double v(double rmeter); double vp(double rmeter); ポテンシャル力を求める関数を宣言その他 void md(); MD を進める関数アイドル関数 ( 何もしてい void idle(); ないときに呼ばれる関数 ) も型宣言する int istep = 1; int mdmotion = 0; MD ステップ数と MD on/off 状態を入れる変数 float rotate[16] = 1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1 ; float obj_pos[] = 0.0, 0.0, 0.0; オブジェクトの回転状態を表す行列などこの辺りは深く解説しないのでそういうものと思って受け入れていただきたい float scl = 1.0; オブジェクトの拡大度合いを示す変数 int main_window; メインウィンドウ ( Atom viewer ウィンドウ : オブジェクト ( 原子 ) を表示するウィンドウ ) の宣言をしておく /** These are the live variables passed into GLUI ***/ int wireframe = 0; オブジェクトの描画方法などを決める変数 Wireframe int segments = 10; が 1 ならワイヤーフレームで 0 ならソリッドで描画 int radius = 10; segment は球の滑らかさ radius は球の半径 GLUI *glui; GLUI_Checkbox GLUI_Spinner GLUI_Spinner *checkbox; *spinner; *spinner2; glui ウィンドウ ( Control ウィンドウ :glui のボタンなどが配置される ) をポインタ宣言しておくその他 glui のパーツ ( ウィジェットと呼ばれる ) も main 関数以外の関数からも参照したい場合にはここでポインタ宣言しておくと良い

3 /* GLUI control callback コールバック関数つまり glui のウィジェット */ ( チェックボックスやボタンなど ) が操作された場合に呼ば void control_cb( int control ) れる関数ここではそれぞれのパーツの状態を取得して表示する方法を例示している printf( "callback: ID = %d, values are below n", control ); printf( " checkbox: %d n", checkbox->get_int_val() ); printf( " spinner: %d n", spinner->get_int_val() ); printf( " spinner2: %d n", spinner2->get_int_val() ); /**************************************** myglutkeyboard() **********/ void myglutkeyboard(unsigned char Key, int x, int y) viewer ウィンドウ上でキーボード入力がされた場合に switch(key) 実行される動作内容 // A few keys here to test the sync_live capability. case 'w': // Toggle wireframe mode wireframe =!wireframe; GLUI_Master.sync_live_all(); この命令は値の変更を glui ウィジェットに反映させるために必要 case 'r': for (int i=0; i<=15; i++) rotate[i]=0; rotate[0]=1;rotate[5]=1;rotate[10]=1;rotate[15]=1; for (int i=0; i<natom; i++) vx[i]=0.0;vy[i]=0.0;vz[i]=0.0; fx[i]=0.0;fy[i]=0.0;fz[i]=0.0; case 'q': exit(0); ; glutpostredisplay(); オブジェクトを描画しなおす命令これがないと r キーを押しても視点がリセットされない ( 再描画されない )

4 /***************************************** myglutmouse() **********/ void myglutmouse(int button, int button_state, int x, int y ) switch (button) マウスボタン押下のイベントに対して行 case GLUT_LEFT_BUTTON: う命令を記述 glutidlefunc(idle); mdmotion = 1; case GLUT_RIGHT_BUTTON: glutidlefunc(0); mdmotion = 0; case GLUT_MIDDLE_BUTTON: if (button_state == GLUT_DOWN) glutidlefunc(0); glutpostredisplay(); mdmotion = 1; else // glutpostredisplay(); glutidlefunc(0); mdmotion = 0; /***************************************** myglutmotion() **********/ void myglutmotion(int x, int y ) マウスの移動イベントに対して行う命令を記述マ //glutpostredisplay(); ウスのドラッグなどに対してアクションを起こしたければここに書く

5 /**************************************** myglutreshape() *************/ void myglutreshape( int x, int y ) マウスの移動イベントに対して行う命令を記述マウスのドラッグなどに対してアクションを起こした glviewport( 0, 0, x, y); ければここに書く glloadidentity(); glortho(-x/200.0, x/200.0, -y/200.0, y/200.0, , 100.0); glutpostredisplay(); void putsphere( float x, float y, float z, float rad, int seg, int solid) gltranslated(x*scl, y*scl, z*scl); 指定された座標に球を書く関数を用意してお if (solid==0) く座標を平行移動し球を書いてまた座標 glutwiresphere( rad*scl, seg, seg ); を戻すという奇妙な操作となっているがと else にかくこうするしかないようである glutsolidsphere( rad*scl, seg, seg ); gltranslated(-x*scl, -y*scl, -z*scl); void idle( void ) glutpostredisplay(); glui->sync_live(); アイドル時に呼ばれる関数常に再描画し続け glui ウィジェットも更新するという命令 /***************************************** myglutdisplay() *****************/ void myglutdisplay( void ) viewer 窓にオブジェクト ( ここでは原子 ) を描画する命令を関数として用意 glclearcolor(.9f,.9f,.9f, 1.0f ); glclear( GL_COLOR_BUFFER_BIT GL_DEPTH_BUFFER_BIT ); glmatrixmode( GL_MODELVIEW ); glpushmatrix(); OpenGL で 3 次元描画を行う際の色設定や座標変換およびそれに関連するモード設定詳細は OpenGL の教科書を参照されたいがこのまま流用できると思われる

6 gltranslatef( obj_pos[0], obj_pos[1], -obj_pos[2] ); glmultmatrixf( rotate ); float rad = (float)radius/100.0; int solid = 1; if (wireframe == 1) solid = 0; if (mdmotion == 1) istep++; md(); for (int i=0; i<natom; i++) float xx=rx[i]*1e9; ここで原子を描画 ( 球を置く ) している float yy=ry[i]*1e9; putsphere(xx, yy, 0.0, rad, segments, solid); glpopmatrix(); glutswapbuffers(); /**************************************** main() ********************/ int main(int argc, char* argv[]) プログラムの本体 (main 関数 ) rx[0] = 0.0e-10; ry[0] = 0.0e-10; rz[0] = 0.0e-10; rx[1] = 5.0e-10; ry[1] = 1.0e-10; rz[1] = 0.0e-10; rx[2] = 5.0e-10; ry[2] = 5.0e-10; rz[2] = 0.0e-10; vx[0] = 0.0e3; vy[0] = 0.0e3 ; vz[0] = 0.0e3; vx[1] = 0.0e3; vy[1] = 0.0e3 ; vz[1] = 0.0e3; vx[2] = 0.0e3; vy[2] = 0.0e3 ; vz[2] = 0.0e3; for (int i=0; i<natom; i++) rx_org[i]=rx[i];ry_org[i]=ry[i];rz_org[i]=rz[i]; std::ofstream fouttraj("trajectory.d"); std::ofstream fout("energy.d");

7 glut を使用するためのおまじな /* Initialize GLUT and create window */ い glutinit(&argc, argv); glutinitdisplaymode( GLUT_RGB GLUT_DOUBLE GLUT_DEPTH ); glutinitwindowposition( 50, 50 ); glutinitwindowsize( 500, 500 ); main_window = glutcreatewindow( "Atom viewer" ); glutdisplayfunc( myglutdisplay ); ウィンドウの定義とイベント glutreshapefunc( myglutreshape ); に対してどの関数が呼ばれるか glutkeyboardfunc( myglutkeyboard ); の定義 glutmotionfunc( myglutmotion ); glutmousefunc( myglutmouse ); 光の当て方の設定 /* Set up OpenGL lights */ GLfloat light0_ambient[] = 0.1f, 0.1f, 0.3f, 1.0f; GLfloat light0_diffuse[] =.6f,.6f, 1.0f, 1.0f; GLfloat light0_position[] = 200.0f, 200.0f, -4.0f, 1.0f; glenable(gl_lighting); glenable(gl_light0); gllightfv(gl_light0, GL_AMBIENT, light0_ambient); gllightfv(gl_light0, GL_DIFFUSE, light0_diffuse); gllightfv(gl_light0, GL_POSITION, light0_position); /* Enable z-buferring */ glenable(gl_depth_test); 陰面消去をするとの設定 /* Here's the GLUI code */ これ以降にウィジェットを書いていく glui = GLUI_Master.create_glui( "Control", 0, 600, 50 ); // name, flags, x, and y new GLUI_StaticText( glui, "Crude MD with GLUI" ); new GLUI_Separator( glui ); checkbox = new GLUI_Checkbox( glui, "Wireframe", &wireframe, 1, control_cb ); spinner = new GLUI_Spinner( glui, "Segments:", &segments, 2, control_cb ); spinner->set_int_limits( 3, 30 );

8 spinner2 = new GLUI_Spinner( glui, "Radius:", &radius, 5, control_cb ); spinner2->set_int_limits( 1, 20 ); GLUI_Panel *panel = new GLUI_Panel(glui,"", false); GLUI_Rotation *view_rot = new GLUI_Rotation(panel, "Rotation", rotate); new GLUI_Column(panel, false); GLUI_Translation *trans_xy = new GLUI_Translation(panel, "Objects XY", GLUI_TRANSLATION_XY, obj_pos ); trans_xy->set_speed(0.005); new GLUI_Column(panel, false); GLUI_Translation *trans_z = new GLUI_Translation(panel, "Objects Z", GLUI_TRANSLATION_Z, &scl ); trans_z->set_speed(0.005); GLUI_Panel *panel2 = new GLUI_Panel(glui,"", false); GLUI_EditText *text = new GLUI_EditText(panel2,"Step", &istep); new GLUI_Button( glui, "Quit", 0,(GLUI_Update_CB)exit ); 各 GLUI ウィンドウにメインウィンドウ glui->set_main_gfx_window( main_window ); の ID を知らせたのちアイドル時の関 GLUI_Master.set_glutIdleFunc( idle ); 数として idle() を設定これらは glui のおまじないと考えてよいかと // Loop glutmainloop(); これによりプログラムは無限ループを実行しイベ fout.close(); ントを待ち続けることになる return EXIT_SUCCESS; void md() double rr, drx, dry, drz; double dt = 1.0e-15; double wm = e-27; // Velet[1]

9 for (int i=0; i<natom; i++) rx[i] = rx[i] + dt * vx[i] + (dt*dt/2.0) * fx[i] / wm; ry[i] = ry[i] + dt * vy[i] + (dt*dt/2.0) * fy[i] / wm; rz[i] = rz[i] + dt * vz[i] + (dt*dt/2.0) * fz[i] / wm; vx[i] = vx[i] + dt/2.0 * fx[i] / wm; vy[i] = vy[i] + dt/2.0 * fy[i] / wm; vz[i] = vz[i] + dt/2.0 * fz[i] / wm; // Force and energy calculation for (int i=0; i<natom; i++) fx[i] = 0.0; fy[i] = 0.0; fz[i] = 0.0; epot[i] = 0.0; for (int i=0; i<natom; i++) for (int j=0; j<natom; j++) if (j!=i) drx = rx[i] - rx[j]; dry = ry[i] - ry[j]; drz = rz[i] - rz[j]; rr = sqrt(drx*drx+dry*dry+drz*drz); fx[i] = fx[i]-vp(rr)/rr*drx; fy[i] = fy[i]-vp(rr)/rr*dry; fz[i] = fz[i]-vp(rr)/rr*drz; epot[i]=epot[i]+v(rr)/2.0; // Velet[2] for (int i=0; i<natom; i++) vx[i] = vx[i] + dt/2.0 * fx[i] / wm; vy[i] = vy[i] + dt/2.0 * fy[i] / wm;

10 vz[i] = vz[i] + dt/2.0 * fz[i] / wm; // Potential and kinetic energies epotsum=0.0; ekinsum=0.0; for (int i=0; i<natom; i++) epotsum=epotsum+epot[i]; ekinsum=ekinsum+wm/2.0*(vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i]); double v(double rmeter) double vval, rang, ep, al, ro, ev; rang = rmeter * 1.0e10; ep = ; al = ; ro = 3.253; ev = e-19; vval=ep*(exp(-2.0*al*(rang-ro))-2.0*exp(-al*(rang-ro))) *ev; return vval; double vp(double rmeter) double vpval, rang, ep, al, ro, ev; rang = rmeter * 1.0e10; ep = ; al = ; ro = 3.253; ev = e-19; vpval=-2.0*al*ep*(exp(-2.0*al*(rang-ro))-exp(-al*(rang-ro))) *ev*1.0e10; return vpval;

#include <stdio.h> 2 #include <stdlib.h> 3 #include <GL/glut.h> 4 Program 1 (OpenGL GameSample001) 5 // 6 static bool KeyUpON = false; // 7 sta

#include <stdio.h> 2 #include <stdlib.h> 3 #include <GL/glut.h> 4 Program 1 (OpenGL GameSample001) 5 // 6 static bool KeyUpON = false; // 7 sta 1 1. 1 #include 2 #include 3 #include 4 Program 1 (OpenGL GameSample001) 5 // 6 static bool KeyUpON = false; // 7 static bool KeyDownON = false; // 8 static bool KeyLeftON

double rx[natom], ry[natom], rz[natom]; 原子の座標 速度 力 ポテンシャルエ double vx[natom], vy[natom], vz[natom]; ネルギーを受ける配列を準備 double fx[natom], fy[natom], fz[natom

double rx[natom], ry[natom], rz[natom]; 原子の座標速度力ポテンシャルエ double vx[natom], vy[natom], vz[natom]; ネルギーを受ける配列を準備 double fx[natom], fy[natom], fz[natom