/* applet No. 978
 *
 * NTP fcc(100)/bcc(110)/hcp(0001) crystal plane
 * - table function - periodic tight-binding potential
 *
 * Created by Ikeuchi Mitsuru on July 16 2006.
 * Copyright (c) 2006 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.07.16  created
 *
 */

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

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

/*-------------------------------------- define gloval -----*/

  Choice ch_tempMode, ch_mat, ch_mat2, ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_temp, sc_scale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  double xMax = 6.4*1.0e-9;
  double yMax = 6.4*1.0e-9;
  double zMax = 6.4*1.0e-9; 

  int resetSW = 0;
  int started = 1;
  int tempMode = 1;
  double contTemp = 300.0;
  double contPress = 0.0;
  int dispMode = 0;
  int MAT1 = 3;
  int MAT2 = 25;

  double dispWidth = 250.0;
  double dispScale = (dispWidth/xMax);
  double viewScale = 1.0;

  int Nmt = 10000;
  int NN = Nmt + 1;
  double t = 0.0;
  double dt = 1.0*1.0e-15;

  int Nsx = 16;
  int Nsy = 16;
  int Nsz = 16;


  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;

  double AU = 1.661e-27;
  double TB[][] = {
  /*  atom    0 mass    1 A(eV) 2 x(eV)  3 p     4 q     5 r0(A) 6 G/xx 7 color 8 f/bcc */
  /* 0 Al*/ { 26.98*AU, 0.1602, 1.5074,  7.5681, 2.7456, 2.8634, 3.552, 0.26, 0.0 },
  /* 1 Ca*/ { 40.08*AU, 0.0492, 0.6842, 11.2115, 2.6841, 3.9471, 3.557, 0.20, 0.0 },
  /* 2 Ni*/ { 58.71*AU, 0.0565, 1.4005, 14.0867, 1.7937, 2.4918, 3.655, 0.39, 0.0 },
  /* 3 Cu*/ { 63.54*AU, 0.0783, 1.2355, 11.1832, 2.3197, 2.5560, 3.589, 0.42, 0.0 },
  /* 4 Sr*/ { 87.62*AU, 0.0257, 0.5557, 12.3406, 1.8105, 4.3027, 3.652, 0.17, 0.0 },
  /* 5 Rh*/ {102.91*AU, 0.1086, 1.9776, 14.1315, 2.5555, 2.6901, 3.567, 0.60, 0.0 },
  /* 6 Pd*/ {106.40*AU, 0.1223, 1.5193, 11.3225, 3.0697, 2.7511, 3.532, 0.63, 0.0 },
  /* 7 Ag*/ {107.87*AU, 0.0812, 1.1081, 11.5597, 2.8316, 2.8890, 3.546, 0.66, 0.0 },
  /* 8 Ir*/ {192.20*AU, 0.2141, 2.7082, 12.8986, 3.4541, 2.7145, 3.513, 0.70, 0.0 },
  /* 9 Pt*/ {195.09*AU, 0.2906, 2.6715, 10.1423, 3.7878, 2.7746, 3.501, 0.73, 0.0 },
  /*10 Au*/ {196.97*AU, 0.1935, 1.7581, 10.4342, 3.9472, 2.8838, 3.497, 0.76, 0.0 },
  /*11 Pb*/ {207.19*AU, 0.0851, 0.8699, 10.0667, 3.3563, 3.5003, 3.517, 0.79, 0.0 },
  /*12 Th*/ {232.04*AU, 0.1200, 2.0937,  9.8344, 1.7972, 3.5951, 3.560, 0.82, 0.0 },

  /*13 Li*/ {  6.94*AU, 0.0488, 0.5729,  6.3675, 1.3969, 3.0391, 3.734, 0.11, 1.0 },
  /*14 Na*/ { 22.99*AU, 0.0353, 0.4083,  7.8536, 1.7477, 3.7158, 3.579, 0.09, 1.0 },
  /*15 K */ { 39.10*AU, 0.0230, 0.3170,  9.3093, 1.6143, 4.6073, 3.633, 0.06, 1.0 },
  /*16 V */ { 50.94*AU, 0.2572, 2.3126,  6.8543, 2.1886, 2.6223, 3.435, 0.30, 1.0 },
  /*17 Cr*/ { 52.00*AU, 0.0407, 1.1012, 13.1852, 0.8993, 2.4981, 4.048, 0.36, 1.0 },
  /*18 Fe*/ { 55.85*AU, 0.1184, 1.5418, 10.7613, 2.0379, 2.4824, 3.479, 0.33, 1.0 },
  /*19 Rb*/ { 85.47*AU, 0.0292, 0.3233,  8.1532, 1.9235, 4.9363, 3.516, 0.03, 1.0 },
  /*20 Nb*/ { 92.91*AU, 0.4546, 3.6302,  5.2702, 2.0552, 2.6033, 3.474, 0.45, 1.0 },
  /*21 Mo*/ { 95.94*AU, 0.2043, 2.5097, 10.0154, 2.0511, 2.7253, 3.475, 0.48, 1.0 },
  /*22 Cs*/ {132.91*AU, 0.0270, 0.3036,  8.4120, 1.9433, 5.3174, 3.509, 0.01, 1.0 },
  /*23 Ba*/ {137.34*AU, 0.0400, 0.6167, 10.1835, 1.5070, 4.3466, 3.680, 0.14, 1.0 },
  /*24 Ta*/ {180.95*AU, 0.3281, 3.3008,  8.2764, 2.2371, 2.8601, 3.422, 0.52, 1.0 },
  /*25 W */ {183.85*AU, 0.2490, 3.2055, 10.3715, 1.9916, 2.7410, 3.493, 0.55, 1.0 },

  /*26 Ti*/ { 47.90*AU, 0.1519, 1.8112,  8.6200, 2.3900, 2.9510, 3.600, 0.24, 2.0 },
  /*27 Zr*/ { 91.22*AU, 0.1934, 2.2792,  8.2500, 2.2490, 3.2320, 3.600, 0.26, 2.0 },
  /*28 Co*/ { 58.93*AU, 0.0950, 1.4880, 11.6040, 2.2860, 2.5070, 3.500, 0.30, 2.0 },
  /*29 Mg*/ { 24.31*AU, 0.0290, 0.4992, 12.8200, 2.2570, 3.1760, 3.600, 0.22, 2.0 },

  /*30 Ne*/ { 20.18*AU, 0.0031, 0.0302, 10.1900, 2.3900, 3.1320, 3.500, 0.90, 0.0 },
  /*31 Ar*/ { 39.95*AU, 0.0042, 0.0661, 14.6500, 2.3100, 3.7170, 3.500, 0.92, 0.0 },
  /*32 Kr*/ { 83.80*AU, 0.0038, 0.0830, 14.3500, 2.0300, 4.0350, 3.500, 0.94, 0.0 },
  /*33 Xe*/ {131.30*AU, 0.0075, 0.1182, 14.7700, 2.1500, 4.3860, 3.500, 0.96, 0.0 }

  };

  int kindMax = 34;


  int kind[] = new int[NN];
  double xx[] = new double[NN]; /* i-th position */
  double yy[] = new double[NN];
  double zz[] = new double[NN];
  double vx[] = new double[NN]; /* i-th velocity */
  double vy[] = new double[NN];
  double vz[] = new double[NN];
  double ffx[] = new double[NN]; /* i-th force    */
  double ffy[] = new double[NN];
  double ffz[] = new double[NN];

  double virxx[] = new double[NN];
  double viryy[] = new double[NN];
  double virzz[] = new double[NN];
  double virxy[] = new double[NN];
  double virxz[] = new double[NN];
  double viryz[] = new double[NN];

  double eDens[] = new double[NN];

  int regNr[][] = new int[NN][250];
  int section[][][][] = new int[Nsx][Nsy][Nsz][100];

  double pairTBL[][][] = new double[kindMax][kindMax][1000];
  double faiTBL[][][] = new double[kindMax][kindMax][1000];
  double ggRef[][] = new double[kindMax][kindMax];

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

  int rdf[][] = new int[4][260];
  int vdf[][] = new int[4][220];
  double meanvdf[][] = new double[4][220];
  double ffvMaxwell[] = new double[220];

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

  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 control functions -----*/

  public void init() {
    resize(630,420);
    setBackground(Color.white);

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

    ch_mat = new Choice();
    ch_mat.addItem("Al");
    ch_mat.addItem("Ca");
    ch_mat.addItem("Ni");
    ch_mat.addItem("Cu");
    ch_mat.addItem("Sr");
    ch_mat.addItem("Rh");
    ch_mat.addItem("Pd");
    ch_mat.addItem("Ag");
    ch_mat.addItem("Ir");
    ch_mat.addItem("Pt");
    ch_mat.addItem("Au");
    ch_mat.addItem("Pb");
    ch_mat.addItem("Th");
    ch_mat.addItem("Li");
    ch_mat.addItem("Na");
    ch_mat.addItem("K");
    ch_mat.addItem("V");
    ch_mat.addItem("Cr");
    ch_mat.addItem("Fe");
    ch_mat.addItem("Rb");
    ch_mat.addItem("Nb");
    ch_mat.addItem("Mo");
    ch_mat.addItem("Cs");
    ch_mat.addItem("Ba");
    ch_mat.addItem("Ta");
    ch_mat.addItem("W");
    ch_mat.addItem("Ti");
    ch_mat.addItem("Zr");
    ch_mat.addItem("Co");
    ch_mat.addItem("Mg");
    ch_mat.addItem("Ne");
    ch_mat.addItem("Ar");
    ch_mat.addItem("Kr");
    ch_mat.addItem("Xe");
    ch_mat.addItemListener(this); 
    ch_mat.select("Cu");

    ch_mat2 = new Choice();
    ch_mat2.addItem("Al");
    ch_mat2.addItem("Ca");
    ch_mat2.addItem("Ni");
    ch_mat2.addItem("Cu");
    ch_mat2.addItem("Sr");
    ch_mat2.addItem("Rh");
    ch_mat2.addItem("Pd");
    ch_mat2.addItem("Ag");
    ch_mat2.addItem("Ir");
    ch_mat2.addItem("Pt");
    ch_mat2.addItem("Au");
    ch_mat2.addItem("Pb");
    ch_mat2.addItem("Th");
    ch_mat2.addItem("Li");
    ch_mat2.addItem("Na");
    ch_mat2.addItem("K");
    ch_mat2.addItem("V");
    ch_mat2.addItem("Cr");
    ch_mat2.addItem("Fe");
    ch_mat2.addItem("Rb");
    ch_mat2.addItem("Nb");
    ch_mat2.addItem("Mo");
    ch_mat2.addItem("Cs");
    ch_mat2.addItem("Ba");
    ch_mat2.addItem("Ta");
    ch_mat2.addItem("W");
    ch_mat2.addItem("Ti");
    ch_mat2.addItem("Zr");
    ch_mat2.addItem("Co");
    ch_mat2.addItem("Mg");
    ch_mat2.addItem("Ne");
    ch_mat2.addItem("Ar");
    ch_mat2.addItem("Kr");
    ch_mat2.addItem("Xe");
    ch_mat2.addItemListener(this); 
    ch_mat2.select("W");

    ch_tempMode = new Choice();
    ch_tempMode.addItem("adiabatic");
    ch_tempMode.addItem("t scaling"); 
    //ch_tempMode.addItem("wall cont.");
    ch_tempMode.addItemListener(this);
    ch_tempMode.select("t scaling"); 

    ch_view = new Choice();
    ch_view.addItem("ball");
    ch_view.addItem("line");
    ch_view.addItem("both");
    ch_view.addItem("temp");
    ch_view.addItem("real");
    ch_view.addItem("velocity");
    ch_view.addItem("v-space");
    ch_view.addItem("Veff");
    ch_view.addItem("e-density");
    ch_view.addItem("surface");
    ch_view.addItemListener(this);
    ch_view.select("ball");

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

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

    sc_temp = new Scrollbar(Scrollbar.HORIZONTAL,30,10,1,210);
    sc_temp.addAdjustmentListener(this);

    sc_scale = new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,310);
    sc_scale.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
	//pnl.add(new Label("  "));
    pnl.add(ch_mat);
	pnl.add(new Label("  "));
    //pnl.add(ch_mat2);
    pnl.add(ch_tempMode);
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    pnl.add(ch_view);
    add(pnl,"North");

    setInitialCondition();
  }

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

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


  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_mat){
      MAT1 = ch_mat.getSelectedIndex();
      resetSW = 1;
    } else if (ev.getSource() == ch_mat2){
      MAT2 = ch_mat2.getSelectedIndex();
      resetSW = 1;
    } else if (ev.getSource() == ch_tempMode){
      tempMode = ch_tempMode.getSelectedIndex();
    } else if (ev.getSource() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

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

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_temp) { 
      contTemp = 10.0*(double)(sc_temp.getValue());
      vAjustment();
    } else if(ev.getSource() == sc_scale) { 
      viewScale = 0.01*(double)(sc_scale.getValue());
    }
  }

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

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

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

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

/* ----------------------------- offPaint -------------------- */

  void offPaint() {
    int gbx,gby;
    double tmp;

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

    tmp = temperature();

    gbx = 60;
    gby = 140;
    if (dispMode==0) {
      ballPlot(gbx,gby,0.6,0);
      rdfPlot(400,100,1.0);

    } else if (dispMode==1) {
      ballPlot(gbx,gby,0.0,2);
      rdfPlot(400,100,1.0);

    } else if (dispMode==2) {
      ballPlot(gbx,gby,0.5,2);
      rdfPlot(400,100,1.0);

    } else if (dispMode==3) {
      ballPlot(gbx,gby,0.6,1);
      rdfPlot(400,100,1.0);

    } else if (dispMode==4) {
      ballPlot(gbx,gby,1.0,0);
      rdfPlot(400,100,1.0);

    } else if (dispMode==5) {
      velocityPlot(gbx,gby, 50.0);
      vdfPlot(400,100,1.0);

    } else if (dispMode==6) {
      vSpacePlot(gbx,gby);
      vdfPlot(400,100,1.0);

    } else if (dispMode==7) {
      ballPlot(gbx,gby,0.6,0);
      drawEffectivePotential(400, 100 ,MAT1,MAT1,1.0);

    } else if (dispMode==8) {
      ballPlot(gbx,gby,1.0,3);
      drawEffectivePotential(400, 100 ,MAT1,MAT1,1.0);

    } else if (dispMode==9) {
      ballPlot(gbx,gby,1.0,4);
      drawEffectivePotential(400, 100 ,MAT1,MAT1,1.0);

    } else if (dispMode==10) {
      ;

    }

	gOff.setColor(Color.getHSBColor((float)(TB[MAT1][7]),0.95f,0.5f));
    gOff.drawString("atom 1",630/6*0+10,40);
    gOff.drawString("r0="+(int)(TB[MAT1][5]*1000+0.5)/1000.0+" A",630/6*0+10,60);
	if ((int)(TB[MAT1][8]+0.5)==0) {
      gOff.drawString("fcc(100)",630/6*0+10,80);
	} else if ((int)(TB[MAT1][8]+0.5)==1) {
      gOff.drawString("bcc(110)",630/6*0+10,80);
	} else if ((int)(TB[MAT1][8]+0.5)==2) {
      gOff.drawString("hcp(0001)",630/6*0+10,80);
	}

/*
    gOff.setColor(Color.getHSBColor((float)(TB[MAT2][7]),0.9f,0.7f));
    gOff.drawString("atom 2",630/6*1+10,40);
    gOff.drawString("r0="+(int)(TB[MAT2][5]*1000+0.5)/1000.0+" A",630/6*1+10,60);
	if ((int)(TB[MAT2][8]+0.5)==0) {
      gOff.drawString("fcc",630/6*1+10,80);
	} else if ((int)(TB[MAT2][8]+0.5)==1) {
      gOff.drawString("bcc",630/6*1+10,80);
	} else if ((int)(TB[MAT2][8]+0.5)==2) {
      gOff.drawString("hcp",630/6*1+10,80);
	}
*/

    gOff.setColor(Color.black);
    gOff.drawString("temp control",630/6*2+10,40);
    gOff.drawString("T_cont="+(int)(contTemp+0.5)+" K",630/6*3+10,40);
    gOff.drawString("scale="+(int)(viewScale*100.0+0.499)+"%",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+Math.floor(t*1.0e15+0.5)+" fs",630/6*2+10,60);
    gOff.drawString("T="+(int)(tmp*10.0+0.5)/10.0+" K",630/6*3+10,60);

    gOff.setColor(Color.blue);
    gOff.drawString("Lx="+(int)(xMax*1.0e13+0.5)/10.0+" pm",100,80);
    gOff.drawString("Ly="+(int)(yMax*1.0e13+0.5)/10.0+" pm",200,80);
    gOff.drawString("Lz="+(int)(zMax*1.0e13+0.5)/10.0+" pm",300,80);

    gOff.drawString("Pxx="+(int)(xxPress()/1.0e6+0.5)+"MPa",430,60);
    gOff.drawString("Pyy="+(int)(yyPress()/1.0e6+0.5)+"MPa",430,75);
    gOff.drawString("Pzz="+(int)(zzPress()/1.0e6+0.5)+"MPa",430,90);
    gOff.drawString("Pxy="+(int)(xyPress()/1.0e6+0.5)+"MPa",530,60);
    gOff.drawString("Pxz="+(int)(xzPress()/1.0e6+0.5)+"MPa",530,75);
    gOff.drawString("Pyz="+(int)(yzPress()/1.0e6+0.5)+"MPa",530,90);

  }

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

  void ballPlot(int gbx, int gby, double size,int mode) {
    int i,j,k,jj,gx,gy,gz,g2x,g2y;
    double sc,sz,bl,rt,rg;

    setWaku();
    rotate(theta,fai);
    zSort();
    if (mode>=3) setEDensity();

    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
	  sz = size*(1.732/2.0)*sc*TB[kind[j]][5]*1.1547*1.0e-10;
	  bl = 1.2*TB[kind[j]][5]*1.0e-10;
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      gz = (int)(sz*(0.2*pz[j]/zMax+0.8));

      if (mode==0) {
        gOff.setColor(Color.getHSBColor((float)(TB[kind[j]][7]),0.95f,(float)(0.5*pz[j]/zMax+0.2)));
      } else if (mode==1) {
	    rt = temp(j)/2000.0; if (rt>0.79) rt = 0.79;
	    gOff.setColor(Color.getHSBColor(
          (float)(TB[kind[j]][7]),0.95f,(float)(rt+0.2)));
	  } else if (mode==2) {
	    gOff.setColor(Color.getHSBColor((float)(TB[kind[j]][7]),0.3f,(float)(0.5*pz[j]/zMax+0.2)));
        for(k=1;k<regNr[i][0];k++) {
          jj = regNr[j][k];
          if (jj>j) {
            if (distance(j,jj)<bl) {
              g2x = (int)(sc*px[jj])+gbx;
              g2y = (int)(sc*py[jj])+gby;
              gOff.drawLine(gx,gy, g2x, g2y);
            }
          }
        }
		gOff.setColor(Color.getHSBColor((float)(TB[kind[j]][7]),0.95f,(float)(0.5*pz[j]/zMax+0.2)));
      } else if (mode==3) {
	    rg = eDens[j]/(ggRef[kind[j]][kind[j]]*ggRef[kind[j]][kind[j]]);
	    gz = (int)(gz*rg);
	    gOff.setColor(Color.getHSBColor((float)(0.66-0.33*rg),0.95f,(float)(0.5*pz[j]/zMax+0.2)));
	  } else if (mode==4) {
	    rg = ggRef[kind[j]][kind[j]];
	    rt = (eDens[j]-rg*rg)/(rg*rg);
		if (rt>0.0) {
		  gz = (int)(gz*1.0*rt);
		} else {
		  gz = (int)(gz*1.0*(-rt));
		}
	    gOff.setColor(Color.getHSBColor((float)(0.33+0.3*rt),0.95f,(float)(0.5*pz[j]/zMax+0.2)));
	  }
      if (gz>0) gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }

    drawWaku(gbx,gby,1);
  }

  double distance(int i,int j) {
    double x,y,z;

    x = xx[i]-xx[j];
    y = yy[i]-yy[j];
    z = zz[i]-zz[j];

    return(Math.sqrt(x*x+y*y+z*z));
  }

  void velocityPlot(int gbx, int gby, double mag) {
    int i,j,k,jj,gx,gy,gz,g2x,g2y;
    double sc,cosTh,sinTh,cosFi,sinFi;
    double p2x,p2y,p2z,rotp2x,rotp2y,rotp2z;

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

    cosTh = Math.cos(theta); sinTh = Math.sin(theta);
    cosFi = Math.cos(fai); sinFi = Math.sin(fai);

    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;

      p2x = xx[j]+mag*vx[j]*20*dt;
      p2y = yy[j]+mag*vy[j]*20*dt;
      p2z = zz[j]+mag*vz[j]*20*dt;
      rotp2y = cosFi*(p2y-cy)+sinFi*(p2z-cz) + cy; 
      rotp2z = -sinFi*(p2y-cy)+cosFi*(p2z-cz) + cz;
      rotp2x = cosTh*(p2x-cx)+sinTh*(rotp2z-cz) + cx; 
      rotp2z = -sinTh*(p2x-cx)+cosTh*(rotp2z-cz) + cz;
      g2x = (int)(sc*rotp2x)+gbx;
      g2y = (int)(sc*rotp2y)+gby;
      gOff.setColor(Color.getHSBColor((float)(TB[kind[j]][7]),0.5f,(float)(0.3*pz[j]/zMax+0.5))); 
	  if (vx[j]>0.0) gOff.setColor(Color.getHSBColor((float)(TB[kind[j]][7]),0.95f,(float)(0.3*pz[j]/zMax+0.5)));
      gOff.drawLine(gx,gy, g2x,g2y);
    }

    drawWaku(gbx,gby,1);
  }


  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);
          if (tempMode==2) {
            tmp = temperature();
            if (tmp>contTemp) {
              gOff.setColor(Color.black); 
            } else {
              gOff.setColor(Color.red);
            }
          }
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
  }

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

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

    i=0;
    sgm = 3.4e-10;

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

/*---------  rotate  ----------*/

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

    cosTh = Math.cos(th); sinTh = Math.sin(th);
    cosFi = Math.cos(fi); sinFi = Math.sin(fi);
    for (i=0; i<Nmt; i++) {
      py[i] = cosFi*(yy[i]-cy)+sinFi*(zz[i]-cz) + cy; 
      pz[i] = -sinFi*(yy[i]-cy)+cosFi*(zz[i]-cz) + cz;
    } 
    for (i=0; i<Nmt; i++) {
      px[i] = cosTh*(xx[i]-cx)+sinTh*(pz[i]-cz) + cx; 
      pz[i] = -sinTh*(xx[i]-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;
    }  
  }

/*----------------------------- v-space plot --------------*/

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

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

    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      gz = (int)(4.0*(0.2*pz[j]/zMax+0.8));
	  if (vx[j]>0) {
	    gOff.setColor(Color.getHSBColor((float)(TB[kind[j]][7]),0.95f,(float)(0.5*pz[j]/zMax+0.2)));
	  } else {
	    gOff.setColor(Color.getHSBColor((float)(TB[kind[j]][7]),0.5f,(float)(0.5*pz[j]/zMax+0.2)));
	  }
      gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
	}

    drawWaku(gbx,gby,1);
  }

/*---------  v-rotate  ----------*/

  void vRotate(double th,double fi) {
    int i;
    double cosTh,sinTh,cosFi,sinFi,sc;

    sc = xMax/2000.0;
    cosTh = Math.cos(th); sinTh = Math.sin(th);
    cosFi = Math.cos(fi); sinFi = Math.sin(fi);
    for (i=0; i<Nmt; i++) {
      py[i] = cosFi*(sc*vy[i]+yMax/2-cy)+sinFi*(sc*vz[i]+xMax/2-cz) + cy; 
      pz[i] = -sinFi*(sc*vy[i]+yMax/2-cy)+cosFi*(sc*vz[i]+zMax/2-cz) + cz;
    } 
    for (i=0; i<Nmt; i++) {
      px[i] = cosTh*(sc*vx[i]+xMax/2-cx)+sinTh*(pz[i]-cz) + cx; 
      pz[i] = -sinTh*(sc*vx[i]+xMax/2-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  ----------*/

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

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

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

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


/* ----------------------- radial distribution function ----------------*/

  void rdfPlot(int gbx,int gby,double mag) {
    int i,rd,rd2;
    double ascale;

    rdfSet();

    ascale = mag*(10000000.0/Nmt);
    gOff.setColor(Color.lightGray);
    for (i=0; i<200; i+=25) {
      gOff.drawLine(gbx+i,gby, gbx+i, gby+250);
    }
    for (i=0; i<250; i+=25) {
      gOff.drawLine(gbx,gby+i, gbx+200, gby+i);
    }
    gOff.setColor(Color.gray);
    gOff.drawRect(gbx,gby,200,250);
    gOff.drawString("1",gbx-5,gby+250+20);
    gOff.drawString("3A",gbx+90,gby+250+20);
    gOff.drawString("5A",gbx+190,gby+250+20);
    gOff.drawString("radial distribution",gbx+30,gby+300);

    gOff.setColor(Color.lightGray);
    for (i=50; i<250; i++) {
	  rd = (int)(rdf[0][i]*(ascale/(i*i))+0.5);
      if (rdf[0][i]>0) {
        gOff.drawLine(gbx+i-50,gby+250-rd, gbx+i-50, gby+250);
      }
    }

    gOff.setColor(Color.getHSBColor((float)(TB[MAT1][7]),0.9f,0.6f));
    for(i=50;i<250;i++) {
      rd = (int)(rdf[1][i]*(ascale/(i*i))+0.5);
	  rd2 = (int)(rdf[1][i+1]*(ascale/(i*i))+0.5);
      if (rd>0) {
        gOff.drawLine(gbx+i-50,gby+250-rd, gbx+i+1-50, gby+250-rd2);
      }
	}
    gOff.setColor(Color.getHSBColor((float)(TB[MAT2][7]),0.9f,0.6f));
    for(i=50;i<250;i++) {
      rd = (int)(rdf[2][i]*(ascale/(i*i))+0.5);
	  rd2 = (int)(rdf[2][i+1]*(ascale/(i*i))+0.5);
      if (rd>0) {
        gOff.drawLine(gbx+i-50,gby+250-rd, gbx+i+1-50, gby+250-rd2);
      }
	}
  }

  void rdfSet() {
    int i,j,k,ir;
    double xi,xj,yi,yj,zi,zj;
    double r2,rij;

    for (i=0; i<260; i++) {
      rdf[0][i]=0; rdf[1][i]=0; rdf[2][i]=0;
    }

    for (i=0; i<Nmt; i++) {
      for(k=1;k<regNr[i][0];k++) {
        j = regNr[i][k];
        r2=(xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])+(zz[i]-zz[j])*(zz[i]-zz[j]);
        rij = Math.sqrt(r2);
        if (rij<5.0e-10) {
          ir = (int)(rij/2.0e-12+0.5);
          rdf[0][ir] +=1;
		  if (kind[i]==MAT1 && kind[j]==MAT1) rdf[1][ir] +=1;
          if (kind[i]==MAT2 && kind[j]==MAT2) rdf[2][ir] +=1;
        }
      }
    }
  }

/*------------------------ verocity distribution function  --------------*/

  void vdfPlot(int gx, int gy, double mag) {
    int i,iy,iy2,hoganWidth,hoganHight;
    double ascale;

    setVelocityDistribution();
    setMaxwell(contTemp);

    ascale = 4000.0/Nmt*mag;
    hoganWidth = 200; hoganHight = 200;

    gOff.setColor(Color.lightGray);
    for (i=0; i<200; i++) {
      iy = (int)(ascale*vdf[0][i]);
      if (iy>0) {
        gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200);
      }
    }

    gOff.setColor(Color.getHSBColor((float)TB[MAT1][7],0.8f,0.6f));
    for (i=0; i<200-1; i++) {
      iy = (int)(4000*meanvdf[1][i]);
	  iy2 = (int)(4000*meanvdf[1][i+1]);
      if (iy>0) {
        gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200-iy2);
      }
    }

/*
    gOff.setColor(Color.getHSBColor((float)TB[MAT2][7],0.8f,0.6f));
    for (i=0; i<200-1; i++) {
      iy = (int)(4000*meanvdf[2][i]);
	  iy2 = (int)(4000*meanvdf[2][i+1]);
      if (iy>0) {
        gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200-iy2);
      }
    }
*/
/*
    gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.6f));
    for (i=0; i<200-1; i++) {
      iy = (int)(4000*ffvMaxwell[i]);
	  iy2 = (int)(4000*ffvMaxwell[i+1]);
      if (iy>0) {
        gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200-iy2);
      }
    }
*/
    gOff.setColor(Color.gray);
    for (i=0; i<200; i+=50) {
      gOff.drawLine(gx+i,gy, gx+i, gy+200);
    }
    for (i=0; i<200; i+=100) {
      gOff.drawLine(gx,gy+200-i, gx+200, gy+200-i);
    }
    gOff.setColor(Color.gray);
    gOff.drawRect(gx,gy,200,200);

    gOff.setColor(Color.black);
    gOff.drawString("0",gx-5,gy+hoganHight+20);
    gOff.drawString("2000 m/s",gx+hoganWidth-30,gy+hoganHight+20);
    gOff.drawString("velocity distribution",gx+30,gy+hoganHight+40);

    gOff.setColor(Color.lightGray);
    gOff.drawString("count",gx+120,gy+20);
    gOff.setColor(Color.getHSBColor((float)TB[MAT1][7],0.8f,0.6f));
    gOff.drawString("mean distr.",gx+120,gy+40);
    //gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.6f));
    //gOff.drawString("Maxwell(Tc)",gx+120,gy+60);
  }

  void setVelocityDistribution() {
    int i,iv;
    double v2,vv;

    for (i=0; i<220; i++) {
      vdf[0][i]=0; vdf[1][i]=0; vdf[2][i]=0;
    }

    for (i=0; i<Nmt; i++) {
      v2 = vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
      vv = Math.sqrt(v2);
      if (vv<2000.0) {
        iv = (int)(vv*0.1);
        vdf[0][iv] +=1;
		if (kind[i]==MAT1) vdf[1][iv] +=1;
		if (kind[i]==MAT2) vdf[2][iv] +=1;
      }
    }

    for (i=0; i<220; i++) {
      meanvdf[0][i] = 0.99*meanvdf[0][i]+0.01*(vdf[0][i]/(Nmt*1.0));
	  meanvdf[1][i] = 0.99*meanvdf[1][i]+0.01*(vdf[1][i]/(Nmt*1.0));
	  meanvdf[2][i] = 0.99*meanvdf[2][i]+0.01*(vdf[2][i]/(Nmt*1.0));
    }
  }

  void setMaxwell(double Tk) {
    int iv;
    double m,kT,a,v;

    m = TB[MAT1][0];
    kT = 1.38e-23*Tk;
    a = m/kT*Math.sqrt(2.0*m/(3.14159*kT));

    for (iv=0; iv<200; iv++) {
      v = (double)(iv*10);
      ffvMaxwell[iv] = a*v*v*Math.exp(-0.5*m*v*v/kT)*10.0;
    }
  }

/*-------------------------- effective potential -----------------*/

  void drawEffectivePotential(int gbx, int gby, int knd1, int knd2, double mag) {
    int i,ip,ipMax,ipMin,v1,v2;
    double ra,fp,fpMax,fpMin;

    gOff.setColor(Color.lightGray);
    for (i=0; i<200; i+=25) {
      gOff.drawLine(gbx+i,gby, gbx+i, gby+250);
    }
    for (i=0; i<250; i+=25) {
      gOff.drawLine(gbx,gby+i, gbx+200, gby+i);
    }
    gOff.setColor(Color.gray);
    gOff.drawRect(gbx,gby,200,250);
    gOff.drawString("1A",gbx-5,gby+250+20);
    gOff.drawString("3A",gbx+90,gby+250+20);
    gOff.drawString("5A",gbx+190,gby+250+20);
    gOff.drawString("0",gbx-20,gby+200+5);
    gOff.drawString("4eV",gbx-26,gby+200-50*4+5);
    gOff.drawString("effective potential",gbx+30,gby+300);

    gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.9f));
    gOff.drawString("pair 2*U",gbx+80,gby+60);
    for(i=0;i<200;i++) {
      ra = 0.02*i+1.0;
      v1 = (int)(50.0*mag*2.0*pairInter(ra,knd1,knd2)); if (v1>220) continue;
	  v2 = (int)(50.0*mag*2.0*pairInter(ra+0.02,knd1,knd2));
      if (v1<220) {
        gOff.drawLine(gbx+i,gby+200-v1, gbx+i, gby+200-v2);
      }
	}

    gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.9f));
    gOff.drawString("Veff",gbx+80,gby+80);
    for(i=0;i<200;i++) {
      ra = 0.02*i+1.0;
      v1 = (int)(50.0*mag*effectiveInter(ra,knd1,knd2)); if (v1>220) continue;
	  v2 = (int)(50.0*mag*effectiveInter(ra+0.02,knd1,knd2));
      if (v1<220) {
        gOff.drawLine(gbx+i,gby+200-v1, gbx+i, gby+200-v2);
      }
	}
  }

/*-------------------------------------  statictics  ----------*/

  double temperature() {
    int i;
    double s;

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += 0.5*TB[kind[i]][0]*(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    }
    return ( s/((double)Nmt*1.5*1.38e-23) );
  }

  double temp(int ii) {
    double s;

    s = 0.5*TB[kind[ii]][0]*(vx[ii]*vx[ii] + vy[ii]*vy[ii] + vz[ii]*vz[ii]);
    return (s/(1.5*1.38e-23));
  }

  double randomTemperature() {
    int i;
    double mvx,mvy,mvz,s;

    mvx=0.0; mvy=0.0; mvz=0.0;
    for (i=0; i<Nmt; i++) {
      mvx += vx[i]; mvy += vy[i]; mvz += vz[i];
    }
    mvx = mvx/Nmt; mvy = mvy/Nmt; mvz = mvz/Nmt;

    s = 0.0;
    for (i=0; i<Nmt; i++) {
      s += 0.5*TB[kind[i]][0]*((vx[i]-mvx)*(vx[i]-mvx) + (vy[i]-mvy)*(vy[i]-mvy) + (vz[i]-mvz)*(vz[i]-mvz));
    }
    return ((double)Math.floor(s/((double)Nmt*1.5*1.38e-23)*10.0)/10.0 );
  }

/*-------------------------------------  v ajustment  ---------*/

  void vAjustment() {
    int i;
    double tmp, r;

    tmp = temperature();
    r = Math.sqrt(contTemp/tmp);
    for (i=0; i<Nmt; i++) {
        vx[i] = r*vx[i];
        vy[i] = r*vy[i];
        vz[i] = r*vz[i];
    }
  }

/*-------------------- pressure  ----------*/

  double xxPress() {
    int i;
    double vir,nkt,vol;

    vir = 0.0; nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      vir += virxx[i];
      nkt += TB[kind[i]][0]*vx[i]*vx[i];
    }

    vol = xMax*yMax*zMax;
    return( (nkt+vir)/vol );
  }

  double yyPress() {
    int i;
    double vir,nkt,vol;

    vir = 0.0; nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      vir += viryy[i];
      nkt += TB[kind[i]][0]*vy[i]*vy[i];
    }

    vol = xMax*yMax*zMax;
    return( (nkt+vir)/vol );
  }

  double zzPress() {
    int i;
    double vir,nkt,vol;

    vir = 0.0; nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      vir += virzz[i];
      nkt += TB[kind[i]][0]*vz[i]*vz[i];
    }

    vol = xMax*yMax*zMax;
    return( (nkt+vir)/vol );
  }

  double xyPress() {
    int i;
    double vir,nkt,vol;

    vir = 0.0; nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      vir += virxy[i];
      nkt += TB[kind[i]][0]*vx[i]*vy[i];
    }

    vol = xMax*yMax*zMax;
    return( (nkt+vir)/vol );
  }

  double xzPress() {
    int i;
    double vir,nkt,vol;

    vir = 0.0; nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      vir += virxz[i];
      nkt += TB[kind[i]][0]*vx[i]*vz[i];
    }

    vol = xMax*yMax*zMax;
    return( (nkt+vir)/vol );
  }

  double yzPress() {
    int i;
    double vir,nkt,vol;

    vir = 0.0; nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      vir += viryz[i];
      nkt += TB[kind[i]][0]*vy[i]*vz[i];
    }

    vol = xMax*yMax*zMax;
    return( (nkt+vir)/vol );
  }

/*-------------------------------  vol control   -------------*/

  void volControl() {
    double d,a,x,y,z;

    d = dt*3.0;
    a = 1.0e-14;

    x = (1.0+a*(xxPress()-contPress))*xMax;
    if (x>xMax+d) {
      x = xMax + d;
    } else if (x<xMax-d) {
      x = xMax - d;
    }
    xMax = x;

    y = (1.0+a*(yyPress()-contPress))*yMax;
    if (y>yMax+d) {
      y = yMax + d;
    } else if (y<yMax-d) {
      y = yMax - d;
    }
    yMax = y;

    z = (1.0+a*(zzPress()-contPress))*zMax;
    if (z>zMax+d) {
      z = zMax + d;
    } else if (z<zMax-d) {
      z = zMax - d;
    }
    zMax = z;

    cx = 0.5*xMax;
    cy = 0.5*yMax;
    cz = 0.5*zMax;

  }

/*-------------------------  set initial condition   ---------------*/

  void setInitialCondition() {
    int i; 

    t = 0.0;

    setTable();

    i = 0;
	if (TB[MAT1][8]==0) { // fcc
      i = setFcc100Plane(i, MAT1, 10, 20, 20, 0.2*xMax);
	} else if (TB[MAT1][8]==1) { // bcc
      i = setBcc110Plane(i, MAT1, 10, 16, 24, 0.2*xMax);
	} else if (TB[MAT1][8]==2) { // hcp
	  i = setHcp0001Plane(i, MAT1, 10, 16, 14, 0.2*xMax);
	}

    Nmt = i-1;

    for (i=0; i<Nmt; i++) {
      if (xx[i] < 0.0) { 
        xx[i] += xMax; 
      }
      if (xx[i] > xMax) {
        xx[i] -= xMax; 
      }
      if (yy[i] < 0.0) { 
        yy[i] += yMax; 
      }
      if (yy[i] > yMax) {
        yy[i] -= yMax;
      }
      if (zz[i] < 0.0) { 
        zz[i] += zMax; 
      }
      if (zz[i] > zMax) {
        zz[i] -= zMax; 
      }
    }

    removeMacroMotion();
  }


  int setFcc100Plane(int ii, int knd, int nnLayer, int ny, int nz, double xPos) {
    int ip,ix,iy,iz; 
    double a,b,c,rm;

    ip = ii;

    a = 0.5*1.414*TB[knd][5]*1.0e-10;
    rm = 10.0/Math.sqrt(TB[knd][0]/AU);

    yMax = 1.0*a*ny;
    zMax = 1.0*a*nz;
    cy = 0.5*yMax;
    cz = 0.5*zMax;


    for (ix=0; ix<nnLayer; ix++) {
      for (iy=0; iy<ny; iy++) {
        for (iz=0; iz<nz; iz++) {
          if ((ix+iy+iz)%2==0) {
		    kind[ip] = knd;
            xx[ip] = a*(double)(ix)+xPos;
            yy[ip] = a*(double)(iy)+0.5*a;
            zz[ip] = a*(double)(iz)+0.5*a;
            vx[ip] = rm*(100.0*Math.random()-50.0);
            vy[ip] = rm*(100.0*Math.random()-50.0);
            vz[ip] = rm*(100.0*Math.random()-50.0);
            ffx[ip] = 0.0;
            ffy[ip] = 0.0;
            ffz[ip] = 0.0;
            ip += 1;
          }
        }
      }
    }
    return ( ip );
  }

  int setBcc110Plane(int ii, int knd, int nnLayer, int ny, int nz, double xPos) {
    int ip,ix,iy,iz; 
    double a,b,c,rm;

    ip = ii;

    a = 1.1547*TB[knd][5]*1.0e-10;
    b = 0.5*Math.sqrt(2.0)*a;
    c = 0.5*a;
    rm = 10.0/Math.sqrt(TB[knd][0]/AU);

    yMax = 1.02*b*ny;
    zMax = 1.02*c*nz;
    cy = 0.5*yMax;
    cz = 0.5*zMax;


    for (ix=0; ix<nnLayer; ix++) {
      for (iy=0; iy<ny; iy++) {
        for (iz=0; iz<nz; iz++) {
          if ((ix+iy+iz)%2==0) {
		    kind[ip] = knd;
            xx[ip] = b*(double)(ix)+xPos;
            yy[ip] = b*(double)(iy)+0.5*b;
            zz[ip] = c*(double)(iz)+0.5*b;
            vx[ip] = rm*(100.0*Math.random()-50.0);
            vy[ip] = rm*(100.0*Math.random()-50.0);
            vz[ip] = rm*(100.0*Math.random()-50.0);
            ffx[ip] = 0.0;
            ffy[ip] = 0.0;
            ffz[ip] = 0.0;
            ip += 1;
          }
        }
      }
    }
    return ( ip );
  }

  int setHcp0001Plane(int ii, int knd, int nnLayer, int ny, int nz, double xPos) {
    int ip,ix,iy,iz;
    double a,b,c,dy,dz,rm;

    ip = ii;
    a = 1.0*TB[knd][5]*1.0e-10;
    b = 0.8660*a;
    c = 0.8165*a;
    rm = 10.0/Math.sqrt(TB[knd][0]/AU);

    yMax = 1.0*b*ny;
    zMax = 1.0*a*nz;
    cy = 0.5*yMax;
    cz = 0.5*zMax;

    for (ix=0; ix<nnLayer; ix++) {
      for (iy=0; iy<ny; iy++) {
        for (iz=0; iz<nz; iz++) {
          kind[ip] = knd;

          xx[ip] = c*ix+xPos;
          yy[ip] = b*iy+0.333333*b*(ix%2) + 0.25*a;
          zz[ip] = a*iz+0.5*a*(iy%2)+0.5*a*(ix%2) + 0.25*a;
          vx[ip] = rm*(100.0*Math.random()-50.0);
          vy[ip] = rm*(100.0*Math.random()-50.0);
          vz[ip] = rm*(100.0*Math.random()-50.0);
          ffx[ip] = 0.0;
          ffy[ip] = 0.0;
          ffz[ip] = 0.0;
          ip += 1;
        }
      }
    }

    return ( ip );
  }


  void setTable() {
    int knd,ir,i,j;
    double TBmass,TBA,TBx,TBp,TBq,TBr0,TBG,rsw,rcut;
	double ra;

    for (knd=0; knd<kindMax; knd++) {
      TBmass = TB[knd][0];
      TBA = TB[knd][1];
      TBx = TB[knd][2];
      TBp = TB[knd][3];
      TBq = TB[knd][4];
      TBr0 = TB[knd][5];
      TBG = TB[knd][6]*TBx;

      rsw = 1.414*TBr0;
      rcut = 1.732*TBr0;
      if ((int)(TB[knd][8]+0.5)==1) {
        rsw = 1.2*TBr0;
        rcut = 1.6*TBr0;
      } else if ((int)(TB[knd][8]+0.5)==2) {
        rsw = 1.9149*TBr0;
        rcut = 2.2362*TBr0;
      }

      ggRef[knd][knd] = TBG;

      for (ir=1; ir<1000; ir++) {
        ra = 0.01*ir;
        pairTBL[knd][knd][ir] = cutoffTB(rsw,rcut,ra)*TBA*Math.exp(-TBp*(ra/TBr0-1.0));
        faiTBL[knd][knd][ir] = cutoffTB(rsw,rcut,ra)*TBx*TBx*Math.exp(-2.0*TBq*(ra/TBr0-1.0));
      }
    }
    for (i=0; i<kindMax; i++) {
      for (j=0; j<kindMax; j++) {
	    if (j>i) {
          for (ir=1; ir<1000; ir++) {
            ra = 0.01*ir;
            pairTBL[i][j][ir] = Math.sqrt(pairTBL[i][i][ir]*pairTBL[j][j][ir]);
            faiTBL[i][j][ir] = Math.sqrt(faiTBL[i][i][ir]*faiTBL[j][j][ir]);
			pairTBL[j][i][ir] = pairTBL[i][j][ir];
			faiTBL[j][i][ir] = faiTBL[i][j][ir];
		  }
		  ggRef[i][j] = Math.sqrt( ggRef[i][i]*ggRef[j][j] );
          ggRef[j][i] = ggRef[i][j];
        }
	  }
	}
  }

  double cutoffTB(double rsw, double rcut, double ra) {
    double x;

    if (ra>=rcut) return (0.0);
    if (ra<=rsw) return (1.0);
	x = (ra-rsw)/(rcut-rsw);
	return (1.0+x*x*x*(-10.0+x*(15.0-6.0*x)));
  }

  void removeMacroMotion() {
    int i;
    double mvx,mvy,mvz;

    mvx = 0.0; mvy = 0.0; mvz = 0.0;
    for(i=0;i<Nmt;i++) {
      mvx += vx[i]; mvy += vy[i]; mvz += vz[i]; 
    }
    mvx = mvx/Nmt; mvy = mvy/Nmt; mvz = mvz/Nmt; 

    for(i=0;i<Nmt;i++) {
      vx[i] -= mvx; vy[i] -= mvy; vz[i] -= mvz; 
    }
  }

/*-------------------------------  molecules motion   ----------------*/

  void timeEvolution(){
    int i;

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

    if (started==1) {

      if (tempMode==1) {
        vAjustment();
      }

      regNearTB();
      for (i=0; i<10; i++) {
        timeStep();
		volControl();
      }

	}
  }

  void timeStep() {
    int i; 
    double dtv2, a2, rr, tmp;

    dtv2 = dt/2.0;
    t = t + dt;

    rr = 1.0;
    if (tempMode==2) {
      tmp = temperature();
      rr = Math.sqrt(contTemp/tmp);
      if (rr<0.5) { rr = 0.5; } else if (rr>1.2){ rr=1.2; };
    }

    for (i=0; i<Nmt; i++) {
      a2 = dtv2/TB[kind[i]][0];
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
      vz[i] += a2*ffz[i];
      xx[i] += vx[i]*dt;
      yy[i] += vy[i]*dt;
      zz[i] += vz[i]*dt;
    }
    forceCalcTB();
    for (i=0; i<Nmt; i++) {
      a2 = dtv2/TB[kind[i]][0];
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
      vz[i] += a2*ffz[i];
    }

    for (i=0; i<Nmt; i++) {
      if (xx[i] < 0.0) { 
        xx[i] += xMax; 
      }
      if (xx[i] > xMax) {
        xx[i] -= xMax; 
      }
      if (yy[i] < 0.0) { 
        yy[i] += yMax; 
      }
      if (yy[i] > yMax) {
        yy[i] -= yMax;
      }
      if (zz[i] < 0.0) { 
        zz[i] += zMax; 
      }
      if (zz[i] > zMax) {
        zz[i] -= zMax; 
      }
    }
  }

  void forceCalcTB() {
    int i,j,k;
    double rij,f,fxij,fyij,fzij;
    double xi,xj,yi,yj,zi,zj;

    for(i=0;i<Nmt;i++) {
      ffx[i] = 0.0;
      ffy[i] = 0.0;
      ffz[i] = 0.0;

      virxx[i] = 0.0;
      viryy[i] = 0.0;
      virzz[i] = 0.0;
      virxy[i] = 0.0;
      virxz[i] = 0.0;
      viryz[i] = 0.0;
    }

    for(i=0;i<Nmt;i++) {
      for(k=1;k<regNr[i][0];k++) {
        j = regNr[i][k];

		  xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[j]; zi = zz[i]; zj = zz[j];
          if (xi-xj>0.5*xMax) xj += xMax;
          if (xi-xj<-0.5*xMax) xj -= xMax;
          if (yi-yj>0.5*yMax) yj += yMax;
          if (yi-yj<-0.5*yMax) yj -= yMax;
          if (zi-zj>0.5*zMax) zj += zMax;
          if (zi-zj<-0.5*zMax) zj -= zMax;
          rij = Math.sqrt((xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj));
		  f = force(rij,kind[i],kind[j]);
          fxij = f*(xi-xj)/rij;
          fyij = f*(yi-yj)/rij;
          fzij = f*(zi-zj)/rij;

          ffx[i] +=  fxij;
          ffy[i] +=  fyij;
          ffz[i] +=  fzij;

          ffx[j] += -fxij;
          ffy[j] += -fyij;
          ffz[j] += -fzij;

          virxx[i] +=  fxij*(xi-xj);
          viryy[i] +=  fyij*(yi-yj);
          virzz[i] +=  fzij*(zi-zj);
          virxy[i] +=  fxij*(yi-yj);
          virxz[i] +=  fxij*(zi-zj);
          viryz[i] +=  fyij*(zi-zj);

      }
    }
  }

  double force(double r, int ki, int kj) { 
    double h,ra,rt;

    h = 0.01;
    ra = r*1.0e10;
    rt = -0.5*(effectiveInter(ra+h,ki,kj)-effectiveInter(ra-h,ki,kj))/h;

    return ( rt*1.6e-19/1.0e-10 );
  }

  double effectiveInter(double ra, int ki, int kj) {
    double gg,y;

    gg = ggRef[ki][kj];
    y = faiInter(ra,ki,kj)/gg;
	return ( 2.0*pairInter(ra,ki,kj) - y + y*y/(4.0*gg) );
  }

  double pairInter(double ra, int ki, int kj) {
    int i;
    double a,hh;

    hh = 0.01;
    i = (int)(ra/hh);
    a = ra - i*hh;

    return ( ((hh-a)*pairTBL[ki][kj][i] + a*pairTBL[ki][kj][i+1])/hh );
  }

  double faiInter(double ra, int ki, int kj) {
    int i;
    double a,hh;

    hh = 0.01;
    i = (int)(ra/hh);
    a = ra - i*hh;

    return ( ((hh-a)*faiTBL[ki][kj][i] + a*faiTBL[ki][kj][i+1])/hh );
  }

/*--------------------  registration of near atoms   ---*/

  void regNearTB() {
    int ip,jp,kp,iq;
    int i,j,k,ii,jj,kk,i0,i1,j0,j1,k0,k1;
    double rrg,rrg2,r2,rcutoff,rcutoff2;
    double xi,xj,yi,yj,zi,zj;

    preRegistration();

    rcutoff = TB[MAT1][5]*1.732; 
    if ((int)(TB[MAT1][8]+0.5)==1) rcutoff = TB[MAT1][5]*1.6;
    if ((int)(TB[MAT1][8]+0.5)==2) rcutoff = TB[MAT1][5]*2.2362;
/*
    rcutoff2 = TB[MAT2][5]*1.732; 
	if ((int)(TB[MAT2][8]+0.5)==1) rcutoff2 = TB[MAT2][5]*1.6;
    if ((int)(TB[MAT2][8]+0.5)==2) rcutoff2 = TB[MAT2][5]*2.2362;
    if (rcutoff<rcutoff2) rcutoff = rcutoff2;
*/

    rrg = rcutoff*1.0e-10;
    rrg += 3000.0*dt*10;
    rrg2 = rrg*rrg;

    for(ip=0;ip<Nmt;ip++) {
      kp = 1;

      i0 = (int)(Nsx*(xx[ip]-rrg)/xMax+(double)Nsx)-Nsx;
      i1 = (int)(Nsx*(xx[ip]+rrg)/xMax);
      j0 = (int)(Nsy*(yy[ip]-rrg)/yMax+(double)Nsy)-Nsy;
      j1 = (int)(Nsy*(yy[ip]+rrg)/yMax);
      k0 = (int)(Nsz*(zz[ip]-rrg)/zMax+(double)Nsz)-Nsz;
      k1 = (int)(Nsz*(zz[ip]+rrg)/zMax);

      for(i=i0;i<=i1;i++) {
        ii = (i+Nsx)%Nsx;
        for(j=j0;j<=j1;j++) {
          jj = (j+Nsy)%Nsy;
          for(k=k0;k<=k1;k++) {
            kk = (k+Nsz)%Nsz;
            for(iq=1;iq<=section[ii][jj][kk][0];iq++) {
              jp = section[ii][jj][kk][iq];
              if (jp>ip) {
                xi = xx[ip]; xj = xx[jp]; yi = yy[ip]; yj = yy[jp]; zi = zz[ip]; zj = zz[jp];
                if (xi-xj>0.5*xMax) xj += xMax;
                if (xi-xj<-0.5*xMax) xj -= xMax;
                if (yi-yj>0.5*yMax) yj += yMax;
                if (yi-yj<-0.5*yMax) yj -= yMax;
                if (zi-zj>0.5*zMax) zj += zMax;
                if (zi-zj<-0.5*zMax) zj -= zMax;
                r2=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
                if (r2<rrg2) {
                  regNr[ip][kp]=jp;
                  kp += 1;
                }
              }
            }
          }
        }
      }
      regNr[ip][0]=kp;
    }

  }

  void preRegistration() {
    int i,j,k,ip,iq;

    for(i=0;i<Nsx;i++) {
      for(j=0;j<Nsy;j++) {
        for(k=0;k<Nsz;k++) {
          section[i][j][k][0] = 0;
        }
      }
    }

    for(ip=0;ip<Nmt;ip++) {
      i = (int)(Nsx*xx[ip]/xMax); if (i>=Nsx) i = Nsx-1;
      j = (int)(Nsy*yy[ip]/yMax); if (j>=Nsy) j = Nsy-1;
      k = (int)(Nsz*zz[ip]/zMax); if (k>=Nsz) k = Nsz-1;
      iq = section[i][j][k][0]+1;
      section[i][j][k][0] = iq;
      section[i][j][k][iq] = ip;
    }

  }

/*--------------------  tight-binding Potential Calculation   ---*/
/*
  double forceTB(double r) { 
    double h,ra,rt;

    h = 0.001;
    ra = r*1.0e10;
    rt = -0.5*(effectivePotential(ra+h)-effectivePotential(ra-h))/h;

    return ( rt*1.6e-19/1.0e-10 );
  }

  double effectivePotential(double ra) {
    double y;

    y = faiTB(ra)/TBG;
	return ( 2.0*pairPotential(ra) - y + y*y/(4.0*TBG) );
  }

  double pairPotential(double ra) {
	return ( cutoff(ra)*TBA*Math.exp(-TBp*(ra/TBr0-1.0)) );
  }

  double faiTB(double ra) {
	return ( cutoff(ra)*TBx*TBx*Math.exp(-2.0*TBq*(ra/TBr0-1.0)) );
  }

  double cutoff(double ra) {
    double x;

    if (ra>=rcut) return (0.0);
    if (ra<=rsw) return (1.0);
	x = (ra-rsw)/(rcut-rsw);
	return (1.0+x*x*x*(-10.0+x*(15.0-6.0*x)));
  }
*/

  void setEDensity() {
    int i,j,k;
    double xi,xj,yi,yj,zi,zj,rij,fij;

    for(i=0;i<Nmt;i++) {
      eDens[i] = 0.0;
	}

    for(i=0;i<Nmt;i++) {
      for(k=1;k<regNr[i][0];k++) {
        j = regNr[i][k];
        xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[j]; zi = zz[i]; zj = zz[j];
        if (xi-xj>0.5*xMax) xj += xMax;
        if (xi-xj<-0.5*xMax) xj -= xMax;
        if (yi-yj>0.5*yMax) yj += yMax;
        if (yi-yj<-0.5*yMax) yj -= yMax;
        if (zi-zj>0.5*zMax) zj += zMax;
        if (zi-zj<-0.5*zMax) zj -= zMax;
        rij = Math.sqrt((xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj));
		fij = faiInter(rij*1.0e10, kind[i], kind[j]);
		eDens[i] += fij;
        eDens[j] += fij;
	  }
	}
  }

/*----- end of molecules motion -----*/

}  

/*----- end of applet -----*/

/* Tight-binding potential
 *
 * original paper
 * M.A. Karolewski ; Radiation Effects and Defects in Solids, Vol. 153, p.239-255 (2001)
 *
 * total energy : E_tot
 *   E_tot = Sum[ U(r_ij), {i,j!=i} ] - Sum[ EB(rho_i), {i} ]
 *
 * replusive pair-potential U(r_ij)
 *   U(r_ij) = A exp(-p(r_ij/r0-1)),
 *
 * cohesive band energy term EB(rho)
 *   EB(rho) = -sqrt(rho),
 * electronic density at host atom i
 *   rho_i = Sum[ fai(r_ij), {j} ]
 *
 * electronic density function
 *   fai(r_ij) = xi^2 exp(-2q(r_ij/r0-1)) 
 *
 * effective pair-potential V_eff(r)
 *   V_eff(r_ij) = 2U(r_ij) - fai(r_ij)/G + fai^2(r_ij)/(4 G^3)
 *    G = Sum[ fai(r_kl), {k!=l} ] 
 *   (electeon density at host atom in reference state)
 *
 * cutoff
 * S(x) = 1 - 6 x^5 + 15 x^4 - 10 x^3,
 *      x = (r - r_sw)/(r_cut - r_sw), (r_sw < r < r_cut)
 *
 * for BCC crystal
 *  r_sw : 3rd neighbour - sqrt(8/3) r0
 *  r_cut : 4th neighbour - sqrt(11/3) r0
 *
 * for FCC crystal
 *  r_sw : 2nd neighbour - sqrt(2) r0
 *  r_cut : 3rd neighbour - sqrt(3) r0
 *
 * for HCP crystal
 *  data adapted :
 *    F. Cleri and V. Rosato; Phys. Rev. B vol.48 pp.22-33 (1993)
 *  r_sw : 1.9149 r0
 *  r_cut : 2.2362 r0
 *
 */