/* applet No. 963
 *
 * fcc crystals - table function - tight-binding potential
 *
 * Created by Ikeuchi Mitsuru on June 25 2006.
 * Copyright (c) 2006 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.06.25  created
 *
 */

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

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

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

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

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  double xMax = 8.0*1.0e-9;
  double yMax = 6.0*1.0e-9;
  double zMax = 6.0*1.0e-9; 

  int resetSW = 0;
  int started = 1;
  int tempMode = 1;
  double contTemp = 300.0;
  int dispMode = 0;
  int MAT1 = 12;
  int MAT2 = 5;

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

  int Ncube = 10;
  // for fcc int Nmt = Ncube*Ncube*Ncube/2;
  int Nmt = Ncube*Ncube*Ncube*2;
  int NN = Nmt + 1;
  double t = 0.0;
  double dt = 1.0*1.0e-15;

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

  double TBmass = 26.98*1.661e-27;
  double TBA = 0.1602;
  double TBx = 1.5074;
  double TBp = 7.5681;
  double TBq = 2.7456;
  double TBr0 = 2.8634;
  double TBG = 3.552*TBx;
  // for fcc double rsw = 1.414*TBr0;
  // for fcc double rcut = 1.732*TBr0;
  double rsw = 1.2*TBr0;
  double rcut = 1.6*TBr0;

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

  double TB[][] = {
  /*  atom    0 mass    1 A(eV) 2 x(eV)  3 p     4 q     5 r0(A) 6 G/xx */
  /* 0 Li*/ {  6.94*AU, 0.0488, 0.5729,  6.3675, 1.3969, 3.0391, 3.734, 0.11 },
  /* 1 Na*/ { 22.99*AU, 0.0353, 0.4083,  7.8536, 1.7477, 3.7158, 3.579, 0.09 },
  /* 2 K */ { 39.10*AU, 0.0230, 0.3170,  9.3093, 1.6143, 4.6073, 3.633, 0.06 },
  /* 3 V */ { 50.94*AU, 0.2572, 2.3126,  6.8543, 2.1886, 2.6223, 3.435, 0.30 },
  /* 4 Cr*/ { 52.00*AU, 0.0407, 1.1012, 13.1852, 0.8993, 2.4981, 4.048, 0.36 },
  /* 5 Fe*/ { 55.85*AU, 0.1184, 1.5418, 10.7613, 2.0379, 2.4824, 3.479, 0.33 },
  /* 6 Rb*/ { 85.47*AU, 0.0292, 0.3233,  8.1532, 1.9235, 4.9363, 3.516, 0.03 },
  /* 7 Nb*/ { 92.91*AU, 0.4546, 3.6302,  5.2702, 2.0552, 2.6033, 3.474, 0.45 },
  /* 8 Mo*/ { 95.94*AU, 0.2043, 2.5097, 10.0154, 2.0511, 2.7253, 3.475, 0.48 },
  /* 9 Cs*/ {132.91*AU, 0.0270, 0.3036,  8.4120, 1.9433, 5.3174, 3.509, 0.01 },
  /*10 Ba*/ {137.34*AU, 0.0400, 0.6167, 10.1835, 1.5070, 4.3466, 3.680, 0.14 },
  /*11 Ta*/ {180.95*AU, 0.3281, 3.3008,  8.2764, 2.2371, 2.8601, 3.422, 0.52 },
  /*12 W */ {183.85*AU, 0.2490, 3.2055, 10.3715, 1.9916, 2.7410, 3.493, 0.55 }
  };

  int kindMax = 13;

  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 virx[] = new double[NN];
  double viry[] = new double[NN];
  double virz[] = 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 virial[] = new double[256];
  int virialp = 0;

  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.addItemListener(this); 
    ch_mat.select("Ir");

    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.addItemListener(this); 
    ch_mat2.select("Ni");
*/

    ch_mat = new Choice();
    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.addItemListener(this); 
    ch_mat.select("W");

    ch_mat2 = new Choice();
    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.addItemListener(this); 
    ch_mat2.select("Fe");

    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 1-1");
    ch_view.addItem("Veff 2-2");
    ch_view.addItem("Veff 1-2");
    ch_view.addItemListener(this);
    ch_view.select("ball");

    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(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 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 = 30;
    gby = 100;
    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);

    } else if (dispMode==8) {
      ballPlot(gbx,gby,0.6,0);
      drawEffectivePotential(400, 100 ,MAT2,MAT2);

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

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

    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.drawString("<P>="+(int)(pressure()/1.0e6)+"MPa",630/6*4+10,60);

    gOff.setColor(Color.black);
    gOff.drawString("box="+(int)(xMax*1.0e10+0.5)/10.0+"x"+(int)(yMax*1.0e10+0.5)/10.0+"x"+(int)(zMax*1.0e10+0.5)/10.0+"nm",430,80);
  }

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

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

    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)));
      }
      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.getHSBColor(0.6f,0.8f,0.8f));
    gOff.drawString("count",gx+120,gy+20);
    gOff.setColor(Color.getHSBColor(0.15f,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) {
    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("5eV",gbx-30,gby+200-125+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)(25.0*2.0*pairInter(ra,knd1,knd2)); if (v1>220) continue;
	  v2 = (int)(25.0*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)(25.0*effectiveInter(ra,knd1,knd2)); if (v1>220) continue;
	  v2 = (int)(25.0*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];
    }
  }

/*-------------------- virial  ----------*/

  void initVirial() {
    int i;

    for (i=0; i<256; i++) {
      virial[i] = 0.0;
    }
    virialp = 0;
  }

  void setVirial() {
    int i;
    double vir;

    vir = 0.0;
    for (i=0; i<Nmt; i++) {
      vir += virx[i]+viry[i]+virz[i];
    }

    virial[virialp] = vir;
    virialp = (virialp+1)%256;
  }

  double pressure() {
    int i;
    double s,vir,tm,vol,p;

    s = 0.0;
    for (i=0; i<256; i++) {
      s += virial[i];
    }
    vir = s/256.0;

    tm = temperature();
    vol = xMax*yMax*zMax;
    p = (Nmt*1.38e-23*tm+vir/3.0)/vol;
    return( p );
  }


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

  void setInitialCondition() {
    int i; 
    double a;

    t = 0.0;

    setTable();

/* for fcc 
    a = 0.5*1.414*TB[MAT1][5]*1.0e-10;

    i = 0;
	i = setFccCrystal(i, MAT1, Ncube/2, 0.48*xMax-a*Ncube/2);
    i = setFccCrystal(i, MAT2, Ncube/2, 0.52*xMax);

    Nmt = i-1;
*/

    a = 1.1547*TB[MAT1][5]*1.0e-10;

    i = 0;
	i = setBccCrystal(i, MAT1, 6, 0.48*xMax-a*6);
    i = setBccCrystal(i, MAT2, 6, 0.52*xMax);

    Nmt = i-1;


    removeMacroMotion();
    initVirial();
  }

  int setFccCrystal(int ip, int knd, int nnCube, double xPos) {
    int ix,iy,iz;
    double a,dy,dz,rm;

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

    for (ix=0; ix<nnCube; ix++) {
      for (iy=0; iy<nnCube; iy++) {
        for (iz=0; iz<nnCube; iz++) {
          if ((ix+iy+iz)%2==0) {
		    kind[ip] = knd;
            xx[ip] = a*(double)(ix)+xPos;
            yy[ip] = a*(double)(iy)+dy;
            zz[ip] = a*(double)(iz)+dz;
            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 setBccCrystal(int ip, int knd, int nnCube, double xPos) {
    int i,j,ix,iy,iz,n2,n3;
    double a,dy,dz,rm;

    a = 1.1547*TB[knd][5]*1.0e-10;
    dy = (yMax-a*nnCube)/2.0;
    dz = (zMax-a*nnCube)/2.0;
    rm = 10.0/Math.sqrt(TB[knd][0]/AU);
    n2 = nnCube*nnCube; n3 = nnCube*nnCube*nnCube;

    for (i=0; i<n3; i++) {
      ix = i/n2;
      iy = (i%n2)/nnCube;
      iz = (i%n2)%nnCube;

	  kind[ip] = knd;
      xx[ip] = a*(double)(ix)+xPos;
      yy[ip] = a*(double)(iy)+dy;
      zz[ip] = a*(double)(iz)+dz;
      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;

	  kind[ip] = knd;
      xx[ip] = a*(double)(ix)+0.5*a+xPos;
      yy[ip] = a*(double)(iy)+0.5*a+dy;
      zz[ip] = a*(double)(iz)+0.5*a+dz;
      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 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;
      // for fcc rsw = 1.414*TBr0;
      // for rcut = 1.732*TBr0;
      rsw = 1.2*TBr0;
      rcut = 1.6*TBr0;

      ggRef[knd][knd] = TBG;

      for (ir=1; ir<1000; ir++) {
        ra = 0.01*ir;
        pairTBL[knd][knd][ir] = pairPotential(ra);
        faiTBL[knd][knd][ir] = faiTB(ra);
      }
    }
    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];
        }
	  }
	}
  }

  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 (tempMode==1) {
      vAjustment();
    }

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

  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] = 0.0; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i];  vz[i] = rr*vz[i]; 
      }
      if (xx[i] > xMax) {
        xx[i] = xMax; vx[i] = -rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = rr*vz[i];  
      }
      if (yy[i] < 0.0) { 
        yy[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i]; 
      }
      if (yy[i] > yMax) {
        yy[i] = yMax; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; vz[i] = rr*vz[i];
      }
      if (zz[i] < 0.0) { 
        zz[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i]; 
      }
      if (zz[i] > zMax) {
        zz[i] = zMax; vx[i] = rr*vx[i]; vy[i] = rr*vy[i]; vz[i] = -rr*vz[i];
      }
    }
  }

  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;

      virx[i] = 0.0;
      viry[i] = 0.0;
      virz[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];

          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;

          virx[i] +=  fxij*(xi-xj);
          viry[i] +=  fyij*(yi-yj);
          virz[i] +=  fzij*(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;
    double xi,xj,yi,yj,zi,zj;

    preRegistration();

    // for fcc rcutoff = TB[MAT1][5]*1.732;
    // for fcc if (TB[MAT1][5]<TB[MAT2][5]) rcutoff = TB[MAT2][5]*1.732;
    rcutoff = TB[MAT1][5]*1.6;
    if (TB[MAT1][5]<TB[MAT2][5]) rcutoff = TB[MAT2][5]*1.6;
    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); if (i0<0) i0 = 0;
      i1 = (int)(Nsx*(xx[ip]+rrg)/xMax); if (i1>=Nsx) i1 = Nsx-1;
      j0 = (int)(Nsy*(yy[ip]-rrg)/yMax); if (j0<0) j0 = 0;
      j1 = (int)(Nsy*(yy[ip]+rrg)/yMax); if (j1>=Nsy) j1 = Nsy-1;
      k0 = (int)(Nsz*(zz[ip]-rrg)/zMax); if (k0<0) k0 = 0;
      k1 = (int)(Nsz*(zz[ip]+rrg)/zMax); if (k1>=Nsz) k1 = Nsz-1;

      for(i=i0;i<=i1;i++) {
        for(j=j0;j<=j1;j++) {
          for(k=k0;k<=k1;k++) {
            for(iq=1;iq<=section[i][j][k][0];iq++) {
              jp = section[i][j][k][iq];
              if (jp>ip) {
                xi = xx[ip]; xj = xx[jp]; yi = yy[ip]; yj = yy[jp]; zi = zz[ip]; zj = zz[jp];
                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)));
  }


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