#GLUIでリアルタイム3D描画＋GUI化

原子構造が今どうなってんのか？をリアルタイムで見られるような分子動力学法(MD)のコードを書きたかった。ついでにGUI機能もついていると使いやすいよね。ってことで色々探したけれど、この２つの要求を同時に満たすsolutionはGLUI利用しか見つからなかった。

###0. コンパイル環境構築(Linux)

#####freeglutのインストール
CentOSなどyumによるパッケージ管理の環境では

$ sudo yum install freeglut
$ sudo yum install freeglut-devel

ubuntuなどのapt-getによるパッケージ管理環境では

$ sudo apt-get install freeglut3-dev

#####gluiのインストール
http://glui.sourceforge.net/ などからソース(最新版はglui-2.36.tgz)を入手。
srcディレクトリに行き，makefileを下記のように修正
　(1) FreeGLUTと書いてある場所の下の2行のコメントアウトを外し、
　　/usr/X11R6/lib -lfreeglut を /usr/lib -lglut に、
　　/usr/X11R6/include を /usr/include に直す。
　(2) GLUT の下の2行をコメントアウト。
make し、所定の場所に libglui.a と glui.h をコピー．（例えば、/ust/lib/libglui.a, /usr/include/GL/glui.h)

#####コンパイル用のスクリプト

gccgl

#! /bin/sh
file=$1
fileh=${file%%.*}
g++ $fileh.cpp -o $fileh -lGL -lglut -lglui -lGLU -lm -L/usr/X11R6/lib -I/usr/X11R6/include -lX11 -lXmu -lpthread

上記のようなスクリプトを書いておけば、例えばtest.cppに対して

$ ./gccgl test.cpp

とすればtestという名のバイナリが作成される。

###1. 最小機能のMD（3原子系）

Morseポテンシャルで相互作用する原子３個のMDを書いてみた。速度Verlet法で時間発展をしている。Morseポテンシャルのパラメータ、原子の質量などは適当である。(一応Alを想定)

コードを見る

threeatomsGLUI_1.cpp

#include <string.h>
#include <iostream>
#include <GL/glui.h>
#include<fstream>
#include<math.h>
#define NATOM 3
double rx[NATOM], ry[NATOM], rz[NATOM]; //Position of atoms
double rx_org[NATOM], ry_org[NATOM], rz_org[NATOM]; //Initial position of atoms
double vx[NATOM], vy[NATOM], vz[NATOM]; //Velocity
double fx[NATOM], fy[NATOM], fz[NATOM]; //Force
double epot[NATOM], epotsum, ekinsum; //Potential and kinetic energy
double v(double rmeter);  //Potential function
double vp(double rmeter); //Force function
void md();
void idle();

int istep = 1; //MD step counter
int mdmotion = 0; //to switch MD On/Off

float rotate[16] = { //Rotation matrix
  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}; //Object position
float scl = 1.0; //Scaling factor for magnification of object

int   main_window; //"Atom viewer" window

/** These are the live variables passed into GLUI ***/
int   wireframe = 0; //Wireframe or Solid drawing
int   segments = 10; //Smoothness of sphere
int   radius = 10; //Radius of each sphere

GLUI *glui;
GLUI_Checkbox   *checkbox;
GLUI_Spinner    *spinner;
GLUI_Spinner    *spinner2;

/* 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(); //To apply renewed settings to glui widgets
    break;
  case 'r': // Reset rotation matrix and initialize atom position
    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++) {
      rx[i]=rx_org[i];ry[i]=ry_org[i];rz[i]=rz_org[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': // Termination
    exit(0);
    break;
  };
  glutPostRedisplay(); //Redraw objects
}

/***************************************** 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 {
      glutIdleFunc(0);
      mdmotion = 0;
    }      
    break;
  }
}

/***************************************** myGlutMotion() **********/
/* Actions for mouse drag must be written here */
void myGlutMotion(int x, int y )
{
  //glutPostRedisplay(); 
}

/**************************************** 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();
}

// Draw atom as sphere
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();
  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");

  /*   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 );

  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); // To arrange wigets horizontally
  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 = 2.0e-15;
  double wm = 1.6726e-27 * 26.98;
  // Velet[1]
  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;
      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;
}

Atom viewerのウィンドウ内でマウスを左クリックするとMDが開始して原子がうにょうにょ動き出す。マウス右クリックで停止。

####コード中身の解説

#####(1)冒頭部分

#include <string.h>
#include <iostream>
#include <GL/glui.h>
#include<fstream>
#include<math.h>
#define NATOM 3
double rx[NATOM], ry[NATOM], rz[NATOM]; //Position of atoms
double rx_org[NATOM], ry_org[NATOM], rz_org[NATOM]; //Initial position of atoms
double vx[NATOM], vy[NATOM], vz[NATOM]; //Velocity
double fx[NATOM], fy[NATOM], fz[NATOM]; //Force
double epot[NATOM], epotsum, ekinsum; //Potential and kinetic energy
double v(double rmeter);  //Potential function
double vp(double rmeter); //Force function
void md();
void idle();

お決まりのinclude文が並んだあと、defineで原子数（NATOM）を指定。
原子配置とその初期値を格納する配列、速度・力ベクトルの配列、エネルギー保存用配列と変数を定義。
Morse型のポテンシャルエネルギー(v)とその微分値(vp)を計算する関数を定義。
MDを1ステップ進める関数をmd()、アイドル状態（何も操作が行われない際）に呼ばれる関数をidle()と定義しておく。

#####(2) 準備の続き

int istep = 1; //MD step counter
int mdmotion = 0; //to switch MD On/Off

float rotate[16] = { //Rotation matrix
  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}; //Object position
float scl = 1.0; //Scaling factor for magnification of object

MDのステップカウントをistep、MDのOn/Off切り替え状態を保存する変数をmdmotionとしておく。
オブジェクト（あるいは視点と見ても良い）の回転角を保存するマトリックス（回転処理に必要）、オブジェクトの位置座標を示すベクトル（並進処理に必要）、拡大率を示す変数（拡大・縮小処理に必要）を定義。
（これらはgluiを使う際の定型と考えよう）

#####(3)ウィンドウ、ウィジェット等の準備

int   main_window; //"Atom viewer" window

/** These are the live variables passed into GLUI ***/
int   wireframe = 0; //Wireframe or Solid drawing
int   segments = 10; //Smoothness of sphere
int   radius = 10; //Radius of each sphere

GLUI *glui;
GLUI_Checkbox   *checkbox;
GLUI_Spinner    *spinner;
GLUI_Spinner    *spinner2;

メインウィンドウをintとして定義しておく。
GLUIでの原子（球）の3D描画に使う変数を定義。wireframeは、球をワイヤーフレームだけにするか或いは表面を塗りつぶすかを切り替え。segmentsは球をどのくらいの数のポリゴンで書くかを指定（数が多いほうが滑らかな球面に見える）。radiusは球の半径。

後に生成するウィンドウやウィジェットのポインタをここで与えておく。
必要ならmain関数以外からも呼べるように、グローバルとして定義している。

#####(4)コールバックの処理

何かイベント（キー入力、ボタン押す等）が発生したときに呼ばれる関数を作っておく。
各ウィジェットの状態を呼び出す例として書いたのみで、ここではコールバック処理を積極的に利用して何かをしているわけでは無い（つまり、無くても良い）。

/* 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() );
}

それぞれのパーツの状態を取得して表示している。

#####(5) キー入力の処理
いくつかのキーボード入力を受け付けるようにしておく。
wは球の描画モード（wireframe/solid）の切り替え、rは初期状態に戻す（視点も）、qはプログラム修了。
GLUI_Master.sync_live_all(); は、値の変更をウィジェットに反映させるために必要。
また最後の方にglutPostRedisplay(); があるが、これが無いとオブジェクトが再描画されない。

/**************************************** 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(); //To apply renewed settings to glui widgets
    break;
  case 'r': // Reset rotation matrix and initialize atom position
    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++) {
      rx[i]=rx_org[i];ry[i]=ry_org[i];rz[i]=rz_org[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': // Termination
    exit(0);
    break;
  };
  glutPostRedisplay(); //Redraw objects
}

#####(6) マウス操作の処理
マウスイベント処理の例。左クリックでMDスタート、右クリックでストップ。
MDスタートした場合はアイドル状態でidle関数が常に呼ばれねばならない（そうしないとMDが動き続けてくれない）ので、その宣言をしている。
中央ボタンの処理は・・これで合ってるんだっけ？（笑）

/***************************************** 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 {
      glutIdleFunc(0);
      mdmotion = 0;
    }      
    break;
  }
}

#####(7) マウスドラッグの処理
マウスドラッグに対して何かしたい場合には、ここに処理を書く。

/***************************************** myGlutMotion() **********/
/* Actions for mouse drag must be written here */
void myGlutMotion(int x, int y )
{
  //glutPostRedisplay(); 
}

#####(8) Reshapeの処理
ウィンドウのサイズが変えられた場合の処理。
ウィンドウサイズ変化に対してオブジェクトを拡大・縮小したい場合などはここを書き換える必要がある（と思われる）。

/**************************************** 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();
}

#####(9) 球描画の関数定義
指定された場所に球を描画する関数を用意しておく。
座標を平行移動し、球を描いて、また座標を戻す、という（奇妙な？）操作となっているが、こうするしかないようである。

// Draw atom as sphere
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);
}

#####(10) アイドル関数の定義
上記（7)で、マウス操作によってアイドル時にどの関数が呼ばれるかを定義した。つまり左クリックするとアイドル時にはidleを常に呼ぶように、と指定したのであるが、ここでその処理の内容を書いておく。
以下のように、idle()では画面の再描画を繰り返すことになっている。再描画の際にmd()が呼ばれるように書いておくことで（後述）、MDの時間発展が繰り返される仕様になっている。

void idle( void )
{
  glutPostRedisplay();
  glui->sync_live();
}

#####(11) 画面再描画の関数定義
画面再描画の際に呼ばれる関数を用意。
glClearColorで、viewportで指定した範囲を指定の色で塗りつぶす（引数はr,g,b,a(透明度)）
glClearで、色バッファ(GL_COLOR_BUFFER_BIT)とデプスバッファ(GL_DEPTH_BUFFER_BIT)をglClearColorで指定した色で塗りつぶす。このあたり、詳細はGLUTの文献を参照のこと。

glMatrixMode(GL_MODELVIEW)として、モデルビュー変換行列を設定。
glPushMatrix〜glTranslatef〜glMultMatrixf〜オブジェクト描画（putSphere）〜glPopMatrix〜glSwapBuffersの流れは座標を回転させてオブジェクトを描画する定形のプロセスとなっており、これについても詳しくはGLUTの文献を参照されたい。

すべての原子を球として描画する前にMDのステップを１つ進め（md()を呼ぶ）、ステップのカウンタを一つ増やしている。

/***************************************** myGlutDisplay() *****************/
void myGlutDisplay( void )
{
  glClearColor( .9f, .9f, .9f, 1.0f );
  glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );

  glMatrixMode( GL_MODELVIEW );

  glPushMatrix();
  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(); 
}

#####(12-1) main関数（冒頭部）
原子の初期状態を与え、結果出力のためのファイルを開いている。

/**************************************** 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");

#####(12-2) main関数（中間部：GLUT関連）
OpenGL描画のための初期設定など。こちらも詳細はGLUT関連の文献をあたっていただきたい。
初期化、ディスプレイモードの設定に続き、ウィンドウポジションとサイズを決定。
その後で、イベントに対してどの関数を呼ぶのかの設定をしている。
その次にライトの設定があり、z-バッファリングをenable（陰面消去をする設定）している。

  /*   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);

#####(12-3) main関数（後半部：GLUI関連）
ここがGLUI関連の記述で、ウィジェットを設定していく部分になっている。

  /*         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 );
  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); // To arrange wigets horizontally
  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にはこのように、オブジェクト回転・並進・拡大縮小のウィジェットが用意されているので便利。
new GLUI_Column(panel, false); をするとウィジェットを横並びにできる（そうでないと下に下にと配置されていく）。

  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 );

MDステップカウンタの表示のためGLUI_EdtTextを用いている。表示させたいだけで値を入力する必要は無いのだが、EditTextとして用意するしかないようである。
終了ボタンも配置。こう書いておくと、押されたときにexitする。

#####(12-4) main関数（終結部）

  glui->set_main_gfx_window( main_window );

  GLUI_Master.set_glutIdleFunc( idle );

  // Loop

  glutMainLoop();
  fout.close();
  return EXIT_SUCCESS;
}

各GLUIウィンドウにメインウィンドウのIDを知らせた後、アイドル時の関数としてidle()を設定。これらはgluiのおまじないと考えてよいかも。
最後にglutMainLoopをしておくことで、プログラムは無限ループを実行し、イベントを待ち続けることになる。

#####(13) MD
速度Verlet法で時間発展をしている。

void md()
{
  double rr, drx, dry, drz;
  double dt = 2.0e-15;
  double wm = 1.6726e-27 * 26.98;
  // Velet[1]
  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;
      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]);
    }
}

#####(14) ポテンシャル関数の処理

Morseポテンシャルで２体間のエネルギーと力を計算する関数。

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;
}