/** applet No. 1005
 *
 * gas particles in 3D glass network - 64 x 64 x 64
 * - fast kinetic Monte-Carlo simulation
 *
 * Created by Ikeuchi Mitsuru on August 27 2006.
 * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.08.27  created
 * ver 0.0.2   2007.05.27  improved code
 *
 */

import java.awt.*;
import java.awt.event.*;
import java.applet.*;

public class gasInGlass64KMC extends Applet
  implements MouseListener, MouseMotionListener,
    ItemListener, ActionListener, AdjustmentListener, Runnable {

  // ------------------------------------ preset field -----------

  Thread th = null; // for run()-paint() thread

  // for event
  Choice ch_view;
  Button bt_reset, bt_restart, bt_startStop;
  Scrollbar sc_nn, sc_qq1, sc_qq2, sc_qq3, sc_kT, sc_zview;

  // for off-paint buffers
  Dimension dim;
  Image imgOff,imgOffInit;
  Graphics gOff,gOffInit;

  // KMC3D
  double t = 0.0;
  double ddh = 1.0e-9;

  int resetSW = 0;
  int restartSW = 0;
  int started = 1;
  int changeFlag = 1;

  double kT = 0.025;

  int NNa = 1000;
  int NNlayer = 8;

  int NNx = 64;
  int NNy = 64;
  int NNz = 64;

  double xMax = ddh*NNx;
  double yMax = ddh*NNy;
  double zMax = ddh*NNz; 

  double totalRate = 1.0e18;

  int NNN = 60000;
  int lattice[][][] = new int[NNx][NNy][NNz];
  int doorwayKind[][][][] = new int[NNx][NNy][NNz][6];
  int atom[][] = new int[NNN][3];
  double trans[] = new double[NNN*6];
  double transitionAtom[] = new double[NNN];
  double activationEnergy[] = { 1000.0, 0.40, 0.30, 0.20 };
  int doorKindInfo[][] = new int[2][4];

  // 3D control ( mouse )
  int dgX, dgY, dgXb, dgYb;
  double cx = 0.5*xMax;
  double cy = 0.5*yMax;
  double cz = 0.5*zMax;
  double theta = -20.0*3.14/180.0;
  double fai = 10.0*3.14/180.0;
  double dtheta = 0.0*3.14/180.0;
  double pai = 3.1415926536;

  // display 3D
  int dispMode = 0;
  int sleepTime = 50;
  int dispWidth = 630;
  int dispHeight = 420;
  int zview = 32;

  double dispScale = (240.0/xMax);
  double viewScale = 1.0;

  int srtz[][] = new int[2][NNN];

  double px[] = new double[NNN]; 
  double py[] = new double[NNN]; 
  double pz[] = new double[NNN]; 

  double wkx[] = new double[8]; 
  double wky[] = new double[8]; 
  double wkz[] = new double[8]; 

  double pwkx[] = new double[8]; 
  double pwky[] = new double[8]; 
  double pwkz[] = new double[8]; 

  int boxp[][] = {
   {0,1},{0,2},{0,4},{1,3},{1,5},{2,3},{2,6},{4,5},{4,6},{3,7},{5,7},{6,7}
  };


  // ------------------------------ applet main thread -----------

  public void init() {

    resize(dispWidth,dispHeight);
    setBackground(Color.white);

    dim = getSize();
    imgOffInit = createImage(dim.width,dim.height);
    gOffInit = imgOffInit.getGraphics();

    imgOff = createImage(dim.width,dim.height);
    gOff = imgOff.getGraphics();

    ch_view = new Choice();
    ch_view.add("gas 2d at z");
    ch_view.add("way 2d at z");
    ch_view.add("gas 3d");
    ch_view.addItemListener(this);
    ch_view.select("gas 2d at z");

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_restart = new Button("restart");
    bt_restart.addActionListener(this); 

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 

    sc_nn= new Scrollbar(Scrollbar.HORIZONTAL,300,10,1,610);
    sc_nn.addAdjustmentListener(this);

    sc_qq1= new Scrollbar(Scrollbar.HORIZONTAL,40,10,30,60);
    sc_qq1.addAdjustmentListener(this);

    sc_qq2= new Scrollbar(Scrollbar.HORIZONTAL,30,10,20,50);
    sc_qq2.addAdjustmentListener(this);

    sc_qq3= new Scrollbar(Scrollbar.HORIZONTAL,20,10,10,40);
    sc_qq3.addAdjustmentListener(this);

    sc_kT= new Scrollbar(Scrollbar.HORIZONTAL,25,10,20,110);
    sc_kT.addAdjustmentListener(this);

    sc_zview= new Scrollbar(Scrollbar.HORIZONTAL,32,10,0,73);
    sc_zview.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(2,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_restart);
    pnl.add(bt_startStop);
    pnl.add(new Label("  "));
    pnl.add(new Label("    temp"));
    pnl.add(ch_view);

    pnl.add(sc_qq1);
    pnl.add(sc_qq2);
    pnl.add(sc_qq3);
    pnl.add(new Label("  "));
    pnl.add(sc_kT);
    pnl.add(sc_zview);
    add(pnl,"North");

    setInitialCondition();
  }

  public void start() {
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

  public void stop() {
    if (th != null) th = null;
  }


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
      if (dispMode==0 || dispMode==1) setoffImageInit(dispMode);
    }
  }

  public void actionPerformed(ActionEvent ev){
    if (ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_restart){ 
      restartSW = 1;
    } else if (ev.getSource() == bt_startStop){
      if (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_nn) { 
      NNa = sc_nn.getValue();
      resetSW = 1;
    } else if (ev.getSource() == sc_qq1) { 
      activationEnergy[1] = 0.01*(double)(sc_qq1.getValue());
      changeFlag = 1;
    } else if (ev.getSource() == sc_qq2) { 
      activationEnergy[2] = 0.01*(double)(sc_qq2.getValue());
      changeFlag = 1;
    } else if (ev.getSource() == sc_qq3) { 
      activationEnergy[3] = 0.01*(double)(sc_qq3.getValue());
      changeFlag = 1;
    } else if (ev.getSource() == sc_kT) { 
      kT = 0.001*(double)(sc_kT.getValue());
      changeFlag = 1;
    } else if (ev.getSource() == sc_zview) { 
      zview = sc_zview.getValue();
      if (dispMode==0 || dispMode==1) setoffImageInit(dispMode);
    }
  }

  public void mousePressed(MouseEvent ev){ }
  public void mouseReleased(MouseEvent ev){
    dgXb = 0; dgYb = 0;
    dgX = 0; dgY = 0;
  }
  public void mouseClicked(MouseEvent ev){ }
  public void mouseEntered(MouseEvent ev){ } 
  public void mouseExited (MouseEvent ev){ }

  public void mouseMoved  (MouseEvent ev){ }

  public void mouseDragged(MouseEvent ev){
    dgXb = dgX; dgYb = dgY;
    dgX=ev.getX();
    dgY=ev.getY();
    if (dgXb==0 && dgYb==0) {
      dgXb = dgX; dgYb = dgY;
    }
    theta += 0.5*3.14/180.0*(dgX-dgXb);
    fai += 0.5*3.14/180.0*(dgY-dgYb);
  }



  // ========================= run() - paint() loop ==============

  public void run() {
    while (th != null) {
      timeEvolution();
      offPaint();
      repaint();
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }


  // --------------------------------------- off-paint -----------

  private void setoffImageInit(int mode) {
    int gx,gy,sz;

    gx = 30; gy = 90; sz = 5;
    gOffInit.setColor(Color.white);
    gOffInit.fillRect(0,0,dim.width,dim.height);

    if (mode==0) {
      dispInitDoorway(gx, gy, sz, zview);
    } else if (mode==1) {
      dispInitDoorway(gx, gy, sz, zview);
      dispInitEasyway(gx, gy, sz, zview);
    }
  }

  private void offPaint() {
    int gx,gy,sz,gbx,k;
    double col;

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,dim.width,dim.height);

    gx = 30; gy = 90; sz = 5;

    if (dispMode==0) {
     gOff.drawImage(imgOffInit,0,0,this);
      dispLattice(gx, gy, sz, zview);

    } else if (dispMode==1) {
      gOff.drawImage(imgOffInit,0,0,this);
      dispLattice(gx, gy, sz, zview);

    } else if (dispMode==2) {
      gasAtomPlot(60,120);

    }

    col = colorOf(1);
    gOff.setColor(Color.getHSBColor((float)col,0.8f,0.8f));
    gOff.drawString("Q1="+(int)(activationEnergy[1]*100.0+0.5)/100.0+" eV",dispWidth/6*0+10,70);
    col = colorOf(2);
    gOff.setColor(Color.getHSBColor((float)col,0.8f,0.8f));
    gOff.drawString("Q2="+(int)(activationEnergy[2]*100.0+0.5)/100.0+" eV",dispWidth/6*1+10,70);
    col = colorOf(3);
    gOff.setColor(Color.getHSBColor((float)col,0.8f,0.8f));
    gOff.drawString("Q3="+(int)(activationEnergy[3]*100.0+0.5)/100.0+" eV",dispWidth/6*2+10,70);

    gOff.setColor(Color.black);
    gOff.drawString("kT="+(int)(kT*1000+0.5)/1000.0+" eV",dispWidth/6*4+10,70);
    gOff.drawString("  ="+(int)(kT*11600+0.5)+" K",dispWidth/6*4+10,90);

    if (dispMode==0 || dispMode==1) {
      gOff.drawString("view  z = "+zview+" ",dispWidth/6*5+10,70);
    } else {
      gOff.drawString("  view",dispWidth/6*5+10,70);
    }

    gbx = 400;
    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*1e12+0.5)/1000.0+" ns",gbx,120);
    //gOff.drawString("No. of gas="+numberOfgasatom()+"",gbx,140);
    gOff.drawString("No. of gas="+NNa+"",gbx,140);

    gOff.setColor(Color.black);
    gOff.drawString("hopping counts / No. of door",gbx,200);
    for (k=1;k<4;k++) {
      col = colorOf(k);
      gOff.setColor(Color.getHSBColor((float)col,0.8f,0.8f));
      gOff.drawString("Q"+k+"="+doorKindInfo[1][k]+"",gbx,200+20*k);
      gOff.drawString(" / "+doorKindInfo[0][k]+"",gbx+100,200+20*k);
    }
  }


  // ------------------------------------ plot methods -----------

  private void dispLattice(int gx, int gy, int span,int k0) {
    int i,j,ir;
    double col;

    gOff.setColor(Color.getHSBColor(0.6f,0.9f,0.9f));
    ir = span-2;
    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        if (lattice[i][j][k0]>=0) {
          gOff.fillOval(gx+span*i+1,gy+span*j+1, ir, ir);
        }
      }
    }
  }

  private double colorOf(int knd) {
    double col;

    if (knd==0) {
      col = 0.01;
    } else if (knd==1) {
      col = 0.05;
    } else if (knd==2) {
      col = 0.15;
    } else if (knd==3) {
      col = 0.33;
    } else {
      col = 0.85;
    }
    return (col);
  }

  // disp init

  private void dispInitDoorway(int gx, int gy, int span,int k0) {
    int i,j,knd;
    double col;

    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        knd = doorwayKind[i][j][k0][2];
        col = colorOf(knd);
        gOffInit.setColor(Color.getHSBColor((float)col,0.8f,0.8f));

        gOffInit.drawLine(gx+span*i,gy+span*j+span, gx+span*i+span-1, gy+span*j+span);

        knd = doorwayKind[i][j][k0][0];
        col = colorOf(knd);
        gOffInit.setColor(Color.getHSBColor((float)col,0.8f,0.8f));
        gOffInit.drawLine(gx+span*(i+1),gy+span*j, gx+span*(i+1), gy+span*j+span-1);
      }
    }

    i=0;
      for(j=0;j<NNy;j++) {
        knd = doorwayKind[i][j][k0][1];
        col = colorOf(knd);
        gOffInit.setColor(Color.getHSBColor((float)col,0.8f,0.8f));
        gOffInit.drawLine(gx,gy+span*j, gx, gy+span*j+span-1);
      }

    for(i=0;i<NNx;i++) {
      j=0;
        knd = doorwayKind[i][j][k0][3];
        col = colorOf(knd);
        gOffInit.setColor(Color.getHSBColor((float)col,0.8f,0.8f));
        gOffInit.drawLine(gx+span*i,gy, gx+span*i+span-1, gy);
    }
  }

  private void dispInitEasyway(int gx, int gy, int span, int k0) {
    int i,j,ir,k,knd;
    double col;

    ir = span-1;
    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        knd = -1;
        for(k=0;k<6;k++) {
          if (doorwayKind[i][j][k0][k]>knd) knd = doorwayKind[i][j][k0][k];
        }
        col = colorOf(knd);
        gOffInit.setColor(Color.getHSBColor((float)col,0.2f,0.8f));
        gOffInit.fillRect(gx+span*i+1,gy+span*j+1, ir, ir);
      }
    }
  }

  // 3d plot

  private void gasAtomPlot(int gbx, int gby) {
    int i,j,gx,gy,gz;
    double sc,sz;

    setWaku();
    rotate(theta,fai);
    zSort();

    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;

    for (i=0; i<NNa; i++) {
      j = srtz[0][i];
	  sz = 4;
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      gz = (int)(sz*(0.2*pz[j]/zMax+0.8));

      gOff.setColor(Color.getHSBColor(0.60f,0.60f,(float)(0.3*pz[j]/zMax+0.5)));
      gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }

    drawWaku(gbx,gby,1);
  }

  // draw waku

  private void drawWaku(int gbx, int gby, int imode) {
    int i,gx,gy,g2x,g2y,farPoint;
    double sc,tmp;

    sc = dispScale*viewScale;
    farPoint = findMin();

    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+gbx;
      gy = (int)(sc*pwky[boxp[i][0]])+gby;
      g2x = (int)(sc*pwkx[boxp[i][1]])+gbx;
      g2y = (int)(sc*pwky[boxp[i][1]])+gby;
      if (imode==0) {
        if (boxp[i][0]==farPoint || boxp[i][1]==farPoint) { 
          gOff.setColor(Color.lightGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      } else if (imode==1) {
        if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
          gOff.setColor(Color.gray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
  }

  private int findMin() {
    int i,im;
    double m;

    im = 0; m = pwkz[im];
    for (i=0; i<8; i++) {
      if (pwkz[i]<m) {
        m = pwkz[i];
        im = i;
      }
    }
    return(im);
  }

  // set waku

  private void setWaku() {
    int i,ix,iy,iz;
    double sgm;

    i=0;
    sgm = ddh;

    for (ix=0; ix<2; ix++) {
      for (iy=0; iy<2; iy++) {
        for (iz=0; iz<2; iz++) {
          wkx[i] = -0.5*sgm+(xMax)*(double)(ix);
          wky[i] = -0.5*sgm+(yMax)*(double)(iy);
          wkz[i] = -0.5*sgm+(zMax)*(double)(iz);
          i += 1;
        }
      }
    }
  }

  // rotate

  private void rotate(double th,double fi) {
    int i;
    double x,y,z,cosTh,sinTh,cosFi,sinFi;

    cosTh = Math.cos(th); sinTh = Math.sin(th);
    cosFi = Math.cos(fi); sinFi = Math.sin(fi);

    for (i=0; i<NNa; i++) {
      y = ddh*atom[i][1]; z = ddh*atom[i][2];
      py[i] = cosFi*(y-cy)+sinFi*(z-cz) + cy; 
      pz[i] = -sinFi*(y-cy)+cosFi*(z-cz) + cz;
    } 
    for (i=0; i<NNa; i++) {
      x = ddh*atom[i][0];
      px[i] = cosTh*(x-cx)+sinTh*(pz[i]-cz) + cx; 
      pz[i] = -sinTh*(x-cx)+cosTh*(pz[i]-cz) + cz;
    } 

    for (i=0; i<8; i++) {
      pwky[i] = cosFi*(wky[i]-cy)+sinFi*(wkz[i]-cz) + cy; 
      pwkz[i] = -sinFi*(wky[i]-cy)+cosFi*(wkz[i]-cz) + cz;
    } 
    for (i=0; i<8; i++) {
      pwkx[i] = cosTh*(wkx[i]-cx)+sinTh*(pwkz[i]-cz) + cx; 
      pwkz[i] = -sinTh*(wkx[i]-cx)+cosTh*(pwkz[i]-cz) + cz;
    }  
  }

  // sort z

  private void zSort() {
    int i,j,w;

    for (i=0; i<NNa; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zMax);
    }
    qSort(0,NNa-1);    
  }

  private void qSort(int i,int j) {
    int k,p;

    if(i==j) return;
    p = i+1;
    while (p<=j && srtz[1][i]==srtz[1][p]) p++;
    if (p<=j) {
      if (srtz[1][i]>=srtz[1][p]) p = i;
      k = qSortPartition(i,j,srtz[1][p]);
      qSort(i,k-1);
      qSort(k,j);
    }
  }

  private int qSortPartition(int i,int j, int x) {
    int l,r,w;

    l = i; r = j;
    while (l<=r) {
      while (l<=j && srtz[1][l]<x) l++;
      while (r>=i && srtz[1][r]>=x) r--;
      if (l>r) break;
      w = srtz[0][l]; srtz[0][l] = srtz[0][r]; srtz[0][r] = w;
      w = srtz[1][l]; srtz[1][l] = srtz[1][r]; srtz[1][r] = w;
      l++; r--;
    }
    return ( l );
  }



  // ================ kinetic Monte-Carlo Method 3D ==============

  //  set initial condition

  private void setInitialCondition() {
    int ia;

    t = 0.0;

    setDoorwayKind();
    initLattice();
    NNa = setGasAtom();

    setAllTransition();
    changeFlag = 0;

    setDoorKindInfo();
    setoffImageInit(dispMode);
  }

  private void initLattice() {
    int i,j,k;

    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        for(k=0;k<NNz;k++) {
          lattice[i][j][k] = -1;
        }
      }
    }
  }

  private int setGasAtom() {
    int i,j,k,ia,nx;

    nx = NNlayer;
    ia = 0;
    for(i=0;i<nx;i++) {
      for(j=0;j<NNy;j++) {
        for(k=0;k<NNz;k++) {
          atom[ia][0] = i; atom[ia][1] = j; atom[ia][2] = k;
          lattice[i][j][k] = ia;
          ia += 1;
        }
      }
    }
    return ( ia );
  }

  private void setDoorwayKind() {
    int i,j,k,ik;

    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        for(k=0;k<NNz;k++) {
          doorwayKind[i][j][k][0] = 0;
          doorwayKind[i][j][k][1] = 0;
          doorwayKind[i][j][k][2] = 0;
          doorwayKind[i][j][k][3] = 0;
          doorwayKind[i][j][k][4] = 0;
          doorwayKind[i][j][k][5] = 0;
        }
      }
    }

    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        for(k=0;k<NNz;k++) {
          ik = 1 + (int)(1.4999*(Math.random()+Math.random()));
          doorwayKind[i][j][k][0] = ik;
          doorwayKind[(i+1)%NNx][j][k][1] = doorwayKind[i][j][k][0];

          ik = 1 + (int)(1.4999*(Math.random()+Math.random()));
          doorwayKind[i][j][k][2] = ik;
          doorwayKind[i][(j+1)%NNy][k][3] = doorwayKind[i][j][k][2];

          ik = 1 + (int)(1.4999*(Math.random()+Math.random()));
          doorwayKind[i][j][k][4] = ik;
          doorwayKind[i][j][(k+1)%NNz][5] = doorwayKind[i][j][k][4];
        }
      }
    }

    i=0;
      for(j=0;j<NNy;j++) {
        for(k=0;k<NNz;k++) {
          doorwayKind[i][j][k][1] = 0;
        }
      }

    i=NNx-1;
      for(j=0;j<NNy;j++) {
        for(k=0;k<NNz;k++) {
          doorwayKind[i][j][k][0] = 0;
        }
      }
  }


  // ---------------------------------- time evolution -----------

  private void timeEvolution() {

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
      started = 1;
    }

    if (restartSW==1) {
      t = 0.0;
      initLattice();
      NNa = setGasAtom();

      setAllTransition();
      setDoorKindInfo();
      changeFlag = 0;
      restartSW = 0;
      started = 1;
    }

    if (started==1) {
      timeStep();
    }
  }

  private void timeStep() {
    int i,nn;
    double deltat;

    nn = NNa/10+1; if (nn>1000) nn = 1000;

    for(i=0;i<nn;i++) {
      deltat = kmcStep();
      t = t + deltat;
    }
    correctTotalRate();
  }

  private void correctTotalRate() {
    int i;
    double s;

    s = 0.0;
    for(i=0;i<NNa;i++) {
      s += transitionAtom[i];
    }
    totalRate = s;
  }

  private double kmcStep() {
    int i,ip;
    double r;

    if (changeFlag==1) {
      totalRate = setAllTransition();
      changeFlag = 0;
    }

    ip = selectTransion(totalRate);
    setMoveDoorKindInfo(ip);
    move(ip);

    r = (1.0-1.0e-12)*Math.random()+1.0e-12;
    return(-Math.log(r)/totalRate);
  }


  // set transition , and calc. sum of transition

  private double setAllTransition() {
    int ia;
    double rsum;

    rsum = 0.0;
    for(ia=0;ia<NNa;ia++) {
      rsum += setTransition(ia);
    }
    return (rsum);
  }

  private double setTransition(int ia) {
    int i,j,k,id,nb;
    double s;

    i = atom[ia][0]; j = atom[ia][1]; k = atom[ia][2];

    if (lattice[(i+1)%NNx][j][k]>=0) {
      trans[ia*6+0] = 0.0;
    } else {
      trans[ia*6+0] = transition(i,j,k,0);
    }

    if (lattice[(i-1+NNx)%NNx][j][k]>=0) {
      trans[ia*6+1] = 0.0;
    } else {
      trans[ia*6+1] = transition(i,j,k,1);
    }

    if (lattice[i][(j+1)%NNy][k]>=0) {
      trans[ia*6+2] = 0.0;
    } else {
      trans[ia*6+2] = transition(i,j,k,2);
    }

    if (lattice[i][(j-1+NNy)%NNy][k]>=0) {
      trans[ia*6+3] = 0.0;
    } else {
      trans[ia*6+3] = transition(i,j,k,3);
    }

    if (lattice[i][j][(k+1)%NNz]>=0) {
      trans[ia*6+4] = 0.0;
    } else {
      trans[ia*6+4] = transition(i,j,k,4);
    }

    if (lattice[i][j][(k-1+NNz)%NNz]>=0) {
      trans[ia*6+5] = 0.0;
    } else {
      trans[ia*6+5] = transition(i,j,k,5);
    }

    s = 0.0;
    for (id=0; id<6; id++) {
      s += trans[ia*6+id];
    }
    transitionAtom[ia] = s;
    return ( s );
  }

  private double transition(int i, int j, int k, int dir) {
    return( transitionRate(activationEnergy[doorwayKind[i][j][k][dir]]) );
  }

  private double transitionRate(double energy) {
    double h,nue0;

    h = 6.626e-34;
    nue0 = kT*1.602e-19/h;

    return( nue0*Math.exp(-energy/kT) );
  }

  // select transion

  private int selectTransion(double rsum) {
    int ia,iap;
    double rnd;

    rnd = rsum*Math.random();

    for (ia=0; ia<NNa; ia++) {
      rnd -= transitionAtom[ia];
      if (rnd<0.0) {
        rnd += transitionAtom[ia];
        for (iap=0; iap<6; iap++) {
          rnd -= trans[ia*6+iap];
          if (rnd<0.0) return (ia*6+iap);
        }
        return(ia*6+5);
      }
    }
    return((NNa-1)*6+5);
  }

  // move (update adatom[][],lattice[][])

  private void move(int ip) {
    int ia,id,i,j,k;
    double rs0,rs1;

    if (trans[ip]>0.0) {
      ia = ip/6;
      id = ip%6;
      i = atom[ia][0]; j = atom[ia][1]; k = atom[ia][2];

      lattice[i][j][k] = -1;
      if (id==0) i = (i+1)%NNx;
      if (id==1) i = (i-1+NNx)%NNx;
      if (id==2) j = (j+1)%NNy;
      if (id==3) j = (j-1+NNy)%NNy;
      if (id==4) k = (k+1)%NNz;
      if (id==5) k = (k-1+NNz)%NNz;
      atom[ia][0] = i; atom[ia][1] = j; atom[ia][2] = k;
      lattice[i][j][k] = ia;
	
      rs0 = rateAround(i,j,k);
      rs1 = setRateAround(i,j,k);
      totalRate = totalRate - rs0 + rs1;
    }
  }

  private double setRateAround(int i0, int j0, int k0) {
    int i,j,k,ii,jj,kk;
    double srate;

    srate = 0.0;
    for (i=i0-2; i<=i0+2; i++) {
      ii = (i+NNx)%NNx;
      for (j=j0-2; j<=j0+2; j++) {
        jj =  (j+NNy)%NNy;
        for (k=k0-2; k<=k0+2; k++) {
          kk =  (k+NNz)%NNz;
          if (lattice[ii][jj][kk]>=0) {
            srate += setTransition(lattice[ii][jj][kk]);
          }
        }
      }
    }
    return ( srate );
  }

  private double rateAround(int i0, int j0, int k0) {
    int i,j,k,ii,jj,kk;
    double srate;

    srate = 0.0;
    for (i=i0-2; i<=i0+2; i++) {
      ii = (i+NNx)%NNx;
      for (j=j0-2; j<=j0+2; j++) {
        jj =  (j+NNy)%NNy;
        for (k=k0-2; k<=k0+2; k++) {
          kk =  (k+NNz)%NNz;
          if (lattice[ii][jj][kk]>=0) {
            srate += transitionAtom[lattice[ii][jj][kk]];
          }
        }
      }
    }
    return ( srate );
  }


  // --------------------------------------- utilities -----------

  private int numberOfgasatom() {
    int i,j,k,n;

    n = 0;
    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        for(k=0;k<NNz;k++) {
          if (lattice[i][j][k]>=0) n += 1;
        }
      }
    }
    return(n);
  }

  private void setDoorKindInfo() {
    int i,j,k,knd,id;

    for(knd=0;knd<4;knd++) {
      doorKindInfo[0][knd] = 0;
      doorKindInfo[1][knd] = 0;
    }

    for(i=0;i<NNx;i++) {
      for(j=0;j<NNy;j++) {
        for(k=0;k<NNz;k++) {
          for(id=0;id<6;id++) {
            knd = doorwayKind[i][j][k][id];
            doorKindInfo[0][knd] += 1;
          }
        }
      }
    }
  }

  private void setMoveDoorKindInfo(int ip) {
    int ia,id,i,j,k,knd;

    ia = ip/6;
    id = ip%6;
    i = atom[ia][0]; j = atom[ia][1]; k = atom[ia][2];

    knd = doorwayKind[i][j][k][id];
    doorKindInfo[1][knd] += 1;
  }

}


/* ----- kinetic Monte-Carlo algorithm -----
 * ( or the Bortz-Kalos-Liebowitz (BKL) algorithm )
 *
 * i-th transition rate in the system : ri
 * total transition rate : RN = sum(ri, {i=1,2,...,N})
 *
 * 1) t = 0
 *
 * 2) select i : randomly (proportional to ri/RN)
 *  calculate R(i) = sum(rj, {j=1,2,...,i})
 *  get a uniform random number u in [0,1]
 *  find i  R(i-1) < u RN <= R(i)
 *
 * 3) execute i-th transition
 *
 * 4) update the time
 *  get a uniform random numner u in (0,1]
 *  dt = -log(u)/R
 *  t = t + dt
 *
 * 5) goto 2)
 *
 * ----------
 *
 * transion rate k
 *  k = nue0*exp(-dE/kT)
 *  nue0 = kT/h , h: Plank constant
 *
 */