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

(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(); void idle(); int istep = 1; int mdmotion = 0; 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;

/** These are the live variables passed into GLUI ***/ int wireframe = 0; int segments = 10; int radius = 10; GLUI *glui; GLUI_Checkbox *checkbox; GLUI_Spinner *spinner; GLUI_Spinner *spinner2; 原子の座標、速度、力、ポテンシャルエネルギーを受ける配列を準備。系全体のポテンシャル、運動エネルギーの変数も用意ポテンシャル・力を求める関数を宣言。その他、 MD を進める関数、アイドル関数（何もしていないときに呼ばれる関数）も型宣言する MD ステップ数と MD on/off 状態を入れる変数オブジェクトの回転状態を表す行列など。この辺りは深く解説しないので、「そういうもの」と思って受け入れていただきたいオブジェクトの拡大度合いを示す変数メインウィンドウ（”Atom viewer”ウィンドウ：オブジェクト（原子）を表示するウィンドウ）の宣言をしておく。オブジェクトの描画方法などを決める変数。Wireframe が 1 ならワイヤーフレームで、0 ならソリッドで描画。 segment は球の滑らかさ、radius は球の半径。 glui ウィンドウ（”Control”ウィンドウ：glui のボタンなどが配置される）をポインタ宣言しておく。その他 glui のパーツ（ウィジェットと呼ばれる）も、 main 関数以外の関数からも参照したい場合には、ここでポインタ宣言しておくと良い。

(3)

/* GLUI control callback */

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)

{

switch(Key) {

// A few keys here to test the sync_live capability. case 'w':

// Toggle wireframe mode wireframe = !wireframe; GLUI_Master.sync_live_all(); break;

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; } break; case 'q': exit(0); break; }; glutPostRedisplay(); コールバック関数、つまり glui のウィジェット（チェックボックスやボタンなど）が操作された場合に呼ばれる関数。ここではそれぞれのパーツの状態を取得して表示する方法を例示している。 viewer ウィンドウ上でキーボード入力がされた場合に実行される動作内容。この命令は、値の変更を glui ウィジェットに反映させるために必要。オブジェクトを描画しなおす命令。これがないと r キーを押しても視点がリセットされない（再描画されない）

(4)

}

/***************************************** myGlutMouse() **********/

void myGlutMouse(int button, int button_state, int x, int y ) { switch (button) { case GLUT_LEFT_BUTTON: glutIdleFunc(idle); mdmotion = 1; break; case GLUT_RIGHT_BUTTON: glutIdleFunc(0); mdmotion = 0; break; case GLUT_MIDDLE_BUTTON: if (button_state == GLUT_DOWN) { glutIdleFunc(0); glutPostRedisplay(); mdmotion = 1; } else { // glutPostRedisplay(); glutIdleFunc(0); mdmotion = 0; } break; } } /***************************************** 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, 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(); }

/***************************************** myGlutDisplay() *****************/ void myGlutDisplay( void )

{

glClearColor( .9f, .9f, .9f, 1.0f );

glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );

glMatrixMode( GL_MODELVIEW ); glPushMatrix(); マウスの移動イベントに対して行う命令を記述。マウスのドラッグなどに対してアクションを起こしたければここに書く。指定された座標に球を書く関数を用意しておく。座標を平行移動し、球を書いて、また座標を戻す、という奇妙な操作となっているが、とにかくこうするしかないようである。アイドル時に呼ばれる関数。常に再描画し続け、glui ウィジェットも更新する、という命令。 viewer 窓にオブジェクト（ここでは原子）を描画する命令を関数として用意。 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[])

{

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"); ここで原子を描画（球を置く）しているプログラムの本体（main 関数）。

(7)

/* 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 ); glut を使用するためのおまじないウィンドウの定義と、イベントに対してどの関数が呼ばれるかの定義。光の当て方の設定。陰面消去をするとの設定。これ以降にウィジェットを書いていく。

(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->set_main_gfx_window( main_window ); GLUI_Master.set_glutIdleFunc( idle ); // Loop glutMainLoop(); fout.close(); return EXIT_SUCCESS; } void md() { double rr, drx, dry, drz; double dt = 1.0e-15; double wm = 1.6726e-27; // Velet[1] これにより、プログラムは無限ループを実行し、イベントを待ち続けることになる。各 GLUI ウィンドウにメインウィンドウの ID を知らせたのち、アイドル時の関数として idle()を設定。これらは glui のおまじないと考えてよいかと。

(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 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]

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 = 0.2703; al = 1.1646; ro = 3.253; ev = 1.6021892e-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 = 0.2703; al = 1.1646; ro = 3.253; ev = 1.6021892e-19; vpval=-2.0*al*ep*(exp(-2.0*al*(rang-ro))-exp(-al*(rang-ro))) *ev*1.0e10; return vpval; }

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

GLUI による MD の GUI 化

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