/** applet No. 1041
 *
 * adatom on fcc(111) / bcc(110) / hcp(0001) crystal plane 
 * - table function - Fixed Gupta potential
 * - multi thread  FGMD :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on October 15 2006.
 * Copyright (c) 2006 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.10.15  created
 * ver 0.0.2   2006.12.02  improved code
 * ver 0.0.3   2007.05.20  added setBoundary()
 * ver 0.0.4   2007.05.27  improved code
 * ver 0.0.5   2007.07.14  improved code
 * ver 0.0.6   2007.09.15  changed virial calculation method
 *
 */

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

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

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

  Thread th = null; // for run()-paint() thread
  FixedGuptaMD md = null;
  FixedGuptaAtom fg = new FixedGuptaAtom();
  DrawGraph3d dg3d = new DrawGraph3d();

  // for event
  Choice ch_tempMode, ch_mat, ch_mat2, ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_temp, sc_scale;

  // for off-paint buffer
  Dimension dim;
  Image imgOff;
  Graphics gOff;
  
  int MAT1 = md.getMAT1();
  int MAT2 = md.getMAT2();

  // mouse 
  int dgX, dgY, dgXb, dgYb;

  // display
  int sleepTime = 100;
  int dispMode = 0;
  double viewScale = 1.0;
  int thCount = 0;


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

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

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

    ch_mat = new Choice();
    for (int i=0; i<fg.getkindMax(); i++) {
      ch_mat.add(fg.getMatStr(i));
    }
    ch_mat.addItemListener(this); 
    ch_mat.select(fg.getMatStr(MAT1));

    ch_mat2 = new Choice();
    for (int i=0; i<fg.getkindMax(); i++) {
      ch_mat2.add(fg.getMatStr(i));
    }
    ch_mat2.addItemListener(this); 
    ch_mat2.select(fg.getMatStr(MAT2));

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

    ch_view = new Choice();
    ch_view.add("ball");
    ch_view.add("line");
    ch_view.add("both");
    ch_view.add("temp");
    ch_view.add("real");
    ch_view.add("velocity");
    ch_view.add("v-space");
    ch_view.add("Veff");
    ch_view.add("e-density");
    ch_view.add("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); 

    bt_step = new Button("step");
    bt_step.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(2,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(bt_step);
    pnl.add(ch_tempMode);
    pnl.add(new Label("  "));
    pnl.add(ch_view);

    pnl.add(ch_mat);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(sc_temp);
    pnl.add(new Label("  "));
    pnl.add(sc_scale);
    add(pnl,"North");
  }

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

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


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

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_mat){
      MAT1 = ch_mat.getSelectedIndex();
      md.setMAT1(MAT1);
      md.reset();
      thCount = 0;
    } else if (ev.getSource() == ch_mat2){
      MAT2 = ch_mat2.getSelectedIndex();
      md.setMAT2(MAT2);
      md.reset();
      thCount = 0;
    } else if (ev.getSource() == ch_tempMode){
      md.setTempMode( ch_tempMode.getSelectedIndex() );
    } else if (ev.getSource() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      md.reset();
      thCount = 0;
    } else if (ev.getSource() == bt_startStop){
      if (md.getStartSW()==0) {
        md.setStartSW(1);
      } else {
        md.setStartSW(0);
      }
    } else if (ev.getSource() == bt_step){
      if (md.getStartSW()==0) {
        md.setStepSW();
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_temp) { 
      md.setContTemp( 10.0*(double)(sc_temp.getValue()) );
    } else if(ev.getSource() == sc_scale) { 
      viewScale = 0.01*(double)(sc_scale.getValue());
      dg3d.setScale(viewScale);
    }
  }

  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;
    }
    dg3d.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
  }



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

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

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

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

  public void offPaint() {
    if (md != null) {
      dg3d.plotGraph3d(gOff, md, fg, dispMode, thCount);
    }
  }

}  



// =============================== draw graph 3D class ===========

class DrawGraph3d {

  private Graphics gOff;
  FixedGuptaMD md;
  FixedGuptaAtom fg;

  private int nMax;
  private double xMax = md.getxMax();
  private double yMax = md.getyMax();
  private double zMax = md.getzMax();

  private int xImageLoc = 60;
  private int yImageLoc = 160;

  private int xPlot2dLoc = 400;
  private int yPlot2dLoc = 150;

  private double plotSize = 200.0;
  private double dispScale = plotSize/yMax;
  private double viewScale = 1.0;

  private double viewTheta = -20.0*Math.PI/180.0;
  private double viewFai = 10.0*Math.PI/180.0;
  private double dtheta = 0.0*Math.PI/180.0;

  private double cosTh = Math.cos(viewTheta); 
  private double sinTh = Math.sin(viewTheta);
  private double cosFi = Math.cos(viewFai);
  private double sinFi = Math.sin(viewFai);

  private double ccx = xMax/2.0;
  private double ccy = yMax/2.0;
  private double ccz = zMax/2.0;


  private int MAT1 = -1;
  private int MAT2 = -1;
  private int rdf[][] = new int[4][260];
  private int vdf[][] = new int[4][220];
  private double meanvdf[][] = new double[4][220];
  private double ffvMaxwell[] = new double[220];

  private int NN = 10000;
  private double px[] = new double[NN]; 
  private double py[] = new double[NN]; 
  private double pz[] = new double[NN]; 
  private int srtz[][] = new int[2][NN];
  private float col[] = new float[NN];

  private double wkx[] = new double[8]; 
  private double wky[] = new double[8]; 
  private double wkz[] = new double[8]; 
  private double pwkx[] = new double[8]; 
  private double pwky[] = new double[8]; 
  private double pwkz[] = new double[8]; 

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

  private int drawnCount = 0;


  // ------------------------------- class constructor -----------

  DrawGraph3d() {
  }


  // ----------------------------------- class methods -----------

  public void plotGraph3d(Graphics gg, FixedGuptaMD fgmd, FixedGuptaAtom fga, int dMode, int thCount) {
    int mdCount;

    gOff = gg;
    fg = fga;
    md = fgmd;
    if (md==null) return;
    mdCount = md.getCount();
    //if (mdCount==drawnCount) return;
    drawnCount = mdCount;

    nMax = md.getNmt();
    MAT1 = md.getMAT1(); MAT2 = md.getMAT2();
    xMax = md.getxMax(); yMax = md.getyMax(); zMax = md.getzMax();
    ccx = xMax/2.0;
    ccy = yMax/2.0;
    ccz = zMax/2.0;

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,630,460);

    if (dMode==0) { // ball
      ball3dPlot(0.6, 0);
      rdfPlot(1.0);

    } else if (dMode==1) { // line
      ball3dPlot(0.0, 2);
      rdfPlot(1.0);

    } else if (dMode==2) { // line + ball
      ball3dPlot(0.5, 2);
      rdfPlot(1.0);

    } else if (dMode==3) { // temperature
      ball3dPlot(0.6, 1);
      rdfPlot(1.0);

    } else if (dMode==4) { // real
      ball3dPlot(1.0, 0);
      rdfPlot(1.0);

    } else if (dMode==5) { // velocity
      ball3dPlot(0.0, 3);
      vdfPlot(1.0);

    } else if (dMode==6) { // v-space
      ball3dPlot(1.0, 6);
      vdfPlot(1.0);

    } else if (dMode==7) { // effective potential
      ball3dPlot(0.6, 0);
      effectivePotentialPlot(MAT1,MAT1,1.0);

    } else if (dMode==8) { // electron density
      ball3dPlot(0.6, 4);
      effectivePotentialPlot(MAT1,MAT1,1.0);

    } else if (dMode==9) { // surface
      ball3dPlot(1.0, 5);
      effectivePotentialPlot(MAT1,MAT1,1.0);

    } else if (dMode==10) { //
      ;

    }

    gOff.setColor(Color.getHSBColor((float)(fg.getColorOf(MAT1)),0.95f,0.5f));
    if (fg.getXtalKindOf(MAT1)==0) {
      gOff.drawString("fcc(111)",630/6*0+10,70);
    } else if (fg.getXtalKindOf(MAT1)==1) {
      gOff.drawString("bcc(110)",630/6*0+10,70);
    } else if (fg.getXtalKindOf(MAT1)==2) {
      gOff.drawString("hcp(0001)",630/6*0+10,70);
    }
    gOff.drawString("r0="+(int)(fg.getDiaOf(MAT1)*1000+0.5)/1000.0+" A",630/6*0+10,90);

    gOff.setColor(Color.getHSBColor((float)(fg.getColorOf(MAT2)),0.95f,0.5f));
    gOff.drawString("adatom="+fg.getMatStr(MAT2)+"",630/6*1+10,70);

    gOff.setColor(Color.black);
    gOff.drawString("T_cont="+(int)(md.getContTemp()+0.5)+" K",630/6*3+10,70);
    gOff.drawString("scale="+(int)(viewScale*100.0+0.499)+"%",630/6*5+10,70);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(md.getTime()*1.0e16+0.5)/10.0+" fs",630/6*1+10,90);
    gOff.drawString("N="+md.getNmt()+" ",630/6*2+10,90);
    gOff.drawString("T="+(int)(md.getTemperature()*10.0+0.5)/10.0+" K",630/6*3+10,90);

    gOff.setColor(Color.black);
    gOff.drawString("box="+(int)(md.getxMax()*1.0e12+0.5)/1000.0+"x" 
	  +(int)(md.getyMax()*1.0e12+0.5)/1000.0+"x" 
	  +(int)(md.getzMax()*1.0e12+0.5)/1000.0+"nm",430,90);

    gOff.setColor(Color.black);
    gOff.drawString("thread MD = "+md.getCount()+" ",20,440);
    gOff.drawString("thread disp = "+thCount+" ",220,440);

    gOff.setColor(Color.blue);
    gOff.drawString("Pxx="+(int)(md.getxxPressure()/1.0e6+0.5)+"MPa",430,110);
    gOff.drawString("Pyy="+(int)(md.getyyPressure()/1.0e6+0.5)+"MPa",430,125);
    gOff.drawString("Pzz="+(int)(md.getzzPressure()/1.0e6+0.5)+"MPa",430,140);
    gOff.drawString("Pxy="+(int)(md.getxyPressure()/1.0e6+0.5)+"MPa",530,110);
    gOff.drawString("Pxz="+(int)(md.getxzPressure()/1.0e6+0.5)+"MPa",530,125);
    gOff.drawString("Pyz="+(int)(md.getyzPressure()/1.0e6+0.5)+"MPa",530,140);
  }


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

  // ball 3d plot
  //
  // ball3dPlot(gOff, md, double size, int mode)
  // size: ball size 0.0 - 1.0(real atom dia)
  // mode: 0-ball, 1-temperature, 2-line, 3-velocity, 
  //       4-electron density, 5-surface, 6-v space
  //

  private void ball3dPlot(double size, int mode) {
    int i,j,k,jj,gx,gy,gz,g2x,g2y;
    double sc,sz,rt,rg,deltat;
    float col;

    if (mode==4 || mode==5) { 
      md.setStatMode(1);
    } else {
      md.setStatMode(0);
    }

    if (mode!=6) { // set position
      for (i=0; i<nMax; i++) {
        px[i] = md.getxx(i);
        py[i] = md.getyy(i);
        pz[i] = md.getzz(i);
      }
    } else { // set velocity
      deltat = 0.5*xMax/1000.0;
      for (i=0; i<nMax; i++) {
        px[i] = 0.5*xMax + md.getvx(i)*deltat;
        py[i] = 0.5*yMax + md.getvy(i)*deltat;
        pz[i] = 0.5*zMax + md.getvz(i)*deltat;
      }
    }

    rotate();
    zSort();

    drawBox(0);

    sc = dispScale*viewScale;
    sz = size*1.0e-10*sc;
    for (i=0; i<nMax; i++) {
      j = srtz[0][i];

      gx = (int)(sc*px[j]) + xImageLoc;
      gy = (int)(sc*py[j]) + yImageLoc;
      gz = (int)(md.getDia(j)*sz*(0.4*pz[j]/zMax+0.6));
      col = (float)(md.getColor(j));

      if (mode==0) { // 0-ball
        gOff.setColor(Color.getHSBColor( col, 0.9f,(float)(0.4*pz[j]/zMax+0.3)));
      } else if (mode==1) {  // 1-temperature
        rt = md.getTemperature(j)/2000.0; if (rt>0.79) rt = 0.79;
        gOff.setColor(Color.getHSBColor( col, 0.9f,(float)(rt+0.2)));
      } else if (mode==2) { // 2-line
        gOff.setColor(Color.getHSBColor( col, 0.3f,(float)(0.4*pz[j]/zMax+0.3)));
        for(k=1;k<md.getnnRegNr(j);k++) {
          jj = md.getRegNr(j,k);
          if (jj>j) {
            if (md.distance(j,jj)<0.6*(md.getDia(j)+md.getDia(jj))*1.0e-10) {
              g2x = (int)(sc*px[jj]) + xImageLoc;
              g2y = (int)(sc*py[jj]) + yImageLoc;
              gOff.drawLine(gx,gy, g2x, g2y);
            }
          }
        }
        gOff.setColor(Color.getHSBColor( col, 0.9f,(float)(0.4*pz[j]/zMax+0.3)));

      } else if (mode==3) { // 3-velocity
        drawVLineOf(j, sc, 20.0);
        gOff.setColor(Color.getHSBColor( col, 0.9f,(float)(0.4*pz[j]/zMax+0.3)));

      } else if (mode==4) { // 4-electron density at atom
        rg = md.getElectronDensityAt(j)/(md.getggRef(j)*md.getggRef(j));
        gz = (int)(gz*rg);
        gOff.setColor(Color.getHSBColor((float)(0.66-0.33*rg),0.9f,(float)(0.4*pz[j]/zMax+0.3)));

      } else if (mode==5) { // 5-surface
        rg = md.getggRef(j);
        rt = (md.getElectronDensityAt(j)-rg*rg)/(rg*rg);
        if (rt>0.0) {
          gz = (int)(gz*1.0*rt);
        } else {
          gz = (int)(gz*1.0*(-rt));
        }
	if (gz>20) gz = 20;
        gOff.setColor(Color.getHSBColor(
          (float)(0.33+0.3*rt),0.95f,(float)(0.5*pz[j]/zMax+0.2)));

      } else if (mode==6) { // 6-v space
        sz = 8.0*size*viewScale;
        gz = (int)(sz*(0.4*pz[j]/zMax+0.3));
        gOff.setColor(Color.getHSBColor(0.05f,0.8f,(float)(0.3*pz[j]/zMax+0.5))); 
        if (md.getvx(j)>0.0) gOff.setColor(Color.getHSBColor(0.66f,0.8f,(float)(0.3*pz[j]/zMax+0.5)));
      }

      if (gz>0) gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }

    drawBox(1);
  }

  private void drawVLineOf(int j, double sc, double mag) {
    int gx,gy,g2x,g2y;
    double p2x,p2y,p2z,rotp2x,rotp2y,rotp2z;

    gx = (int)(sc*px[j]) + xImageLoc;
    gy = (int)(sc*py[j]) + yImageLoc;
    p2x = md.getxx(j)+mag*md.getvx(j)*mag*1.0e-15;
    p2y = md.getyy(j)+mag*md.getvy(j)*mag*1.0e-15;
    p2z = md.getzz(j)+mag*md.getvz(j)*mag*1.0e-15;
    rotp2y = cosFi*(p2y-ccy)+sinFi*(p2z-ccz) + ccy; 
    rotp2z = -sinFi*(p2y-ccy)+cosFi*(p2z-ccz) + ccz;
    rotp2x = cosTh*(p2x-ccx)+sinTh*(rotp2z-ccz) + ccx; 
    rotp2z = -sinTh*(p2x-ccx)+cosTh*(rotp2z-ccz) + ccz;
    g2x = (int)(sc*rotp2x) + xImageLoc;
    g2y = (int)(sc*rotp2y) + yImageLoc;
    gOff.setColor(Color.getHSBColor(0.05f,0.8f,(float)(0.3*pz[j]/zMax+0.5))); 
    if (md.getvx(j)>0.0) gOff.setColor(Color.getHSBColor(0.66f,0.8f,(float)(0.3*pz[j]/zMax+0.5)));
    gOff.drawLine(gx,gy, g2x,g2y);
  }

  // 3d plot utility 

  private void rotate() {
    int i;
    double x,y,z;

    for (i=0; i<nMax; i++) {
      x = px[i]; y = py[i]; z = pz[i];
      py[i] = cosFi*(y-ccy)+sinFi*(z-ccz) + ccy; 
      pz[i] = -sinFi*(y-ccy)+cosFi*(z-ccz) + ccz;
      px[i] = cosTh*(x-ccx)+sinTh*(pz[i]-ccz) + ccx; 
      pz[i] = -sinTh*(x-ccx)+cosTh*(pz[i]-ccz) + ccz;
    }  
  }

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

    nn = nMax;
    for (i=0; i<nn; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zMax);
    }
    qSort(0,nn-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 );
  }

 // 

  private int drawBox(int imode) {
    int i,gx,gy,g2x,g2y,farPoint;
    double sc,tmp;

    if (imode==0) {
      setBox();
      rotateBox();
    }

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

    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+xImageLoc;
      gy = (int)(sc*pwky[boxp[i][0]])+yImageLoc;
      g2x = (int)(sc*pwkx[boxp[i][1]])+xImageLoc;
      g2y = (int)(sc*pwky[boxp[i][1]])+yImageLoc;
      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.darkGray);
          if (md.getTempMode()==2) {
            tmp = md.getTemperature();
            if (tmp>md.getContTemp()) {
              gOff.setColor(Color.black); 
            } else {
              gOff.setColor(Color.red);
            }
          }
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
    return farPoint;
  }

  private int findFarPoint() {
    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);
  }

  private void setBox() {
    int i,ix,iy,iz;
    double ballSize;

    i=0;
    ballSize = 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*ballSize+(xMax+ballSize)*(double)(ix);
          wky[i] = -0.5*ballSize+(yMax+ballSize)*(double)(iy);
          wkz[i] = -0.5*ballSize+(zMax+ballSize)*(double)(iz);
          i += 1;
        }
      }
    }
  }

  private void rotateBox() {
    int i;

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

  // radial distribution function

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

    md.rdfSet(rdf);
    gbx = xPlot2dLoc; gby = yPlot2dLoc;

    ascale = mag*(10000000.0/nMax);
    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(md.getColorOf(MAT1),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(md.getColorOf(MAT2),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);
      }
    }
  }

  // verocity distribution function

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

    setVelocityDistribution();
    setMaxwell( md.getContTemp() );

    gx = xPlot2dLoc; gy = yPlot2dLoc;

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

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

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

    // MAT1
    gOff.setColor(Color.getHSBColor(md.getColorOf(MAT1),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);
      }
    }

    // MAT2
    gOff.setColor(Color.getHSBColor(md.getColorOf(MAT2),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);
      }
    }

    // grid
    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(md.getColorOf(MAT1),0.8f,0.6f));
    gOff.drawString("mean distr.",gx+120,gy+40);
    //gOff.setColor(Color.getHSBColor(md.getColorOf(MAT2),0.8f,0.6f));
    //gOff.drawString("mean distr.",gx+120,gy+60);
    gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.6f));
    gOff.drawString("Maxwell(Tc)",gx+120,gy+80);
  }

  private void setVelocityDistribution() {
    int nn,i,iv,knd;
    double vx,vy,vz,v2,vv;

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

    nn = md.getNmt();
    for (i=0; i<nn; i++) {
      vx = md.getvx(i); vy = md.getvy(i); vz = md.getvz(i);
      knd = md.getKind(i);
      v2 = vx*vx+vy*vy+vz*vz;
      vv = Math.sqrt(v2);
      if (vv<2000.0) {
        iv = (int)(vv*0.1);
        vdf[0][iv] +=1;
        if (knd==MAT1) vdf[1][iv] +=1;
        if (knd==MAT2) vdf[2][iv] +=1;
      }
    }

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

  }

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

    m = md.getMassOf(MAT1);
    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

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

    gbx = xPlot2dLoc; gby = yPlot2dLoc;

    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*md.getPairInter(ra,knd1,knd2)); if (v1>220) continue;
      v2 = (int)(50.0*mag*2.0*md.getPairInter(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*md.getEffectiveInter(ra,knd1,knd2)); if (v1>220) continue;
      v2 = (int)(50.0*mag*md.getEffectiveInter(ra+0.02,knd1,knd2));
      if (v1<220) {
        gOff.drawLine(gbx+i,gby+200-v1, gbx+i, gby+200-v2);
      }
    }
  }


  // ------------------------------------- I/O methods -----------

  public void setViewDelta(double xDisplacement, double yDisplacement) {
    viewTheta += 0.5*3.14159/180.0*xDisplacement;
    viewFai += 0.5*3.14159/180.0*yDisplacement;
    cosTh = Math.cos(viewTheta); 
    sinTh = Math.sin(viewTheta);
    cosFi = Math.cos(viewFai);
    sinFi = Math.sin(viewFai);
  }

  public void setLoc(int gx, int gy) { xImageLoc = gx; yImageLoc = gy; }
  public void setScale(double sc) { viewScale = sc; }

}



// ================================== Gupta atom class ===========

class FixedGuptaAtom {

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

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

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

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

  };
  private static final int kindMax = 34;

  private static final String matStr[] = {
    "Al","Ca","Ni","Cu","Sr","Rh","Pd","Ag","Ir","Pt","Au","Pb","Th",
    "Li","Na","K","V","Cr","Fe","Rb","Nb","Mo","Cs","Ba","Ta","W",
    "Ti","Zr","Co","Mg","Ne","Ar","Kr","Xe"
  };

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


  // ------------------------------- class constructor -----------

  FixedGuptaAtom() {
    setTable();
  }

  //

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

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


  // ----------------------------------- class methods -----------

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

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

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

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

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


  // ------------------------------------- I/O methods -----------

  public double getAU() { return AU; } 
  public int getkindMax() { return kindMax; }
  public String getMatStr(int knd) { return matStr[knd]; }

  public double getggRefOf(int knd) { return ggRef[knd][knd]; }
  public double getMassOf(int knd) { return TB[knd][0]; }
  public double getDiaOf(int knd) { return TB[knd][5]; }
  public double getColorOf(int knd) { return TB[knd][7]; }
  public int getXtalKindOf(int knd) { return (int)(TB[knd][8]+0.5); }
  
}



// ======== molecular dynamics - fixed Gupta potential ===========

class FixedGuptaMD extends Thread {
  private FixedGuptaAtom fg;

  private static int Nmt = 10000;
  private static double xMax = 5.0*1.0e-9;
  private static double yMax = 5.0*1.0e-9;
  private static double zMax = 5.0*1.0e-9;
  private static int MAT1 = 25; // W
  private static int MAT2 = 23; // Ba

  private double t = 0.0;
  private double dt = 1.0*1.0e-15;

  private int statMode = 1;
  private int tempMode = 1;
  private double contTemp = 300.0;
  private double contPress = 0.0;

  private double xxVirial = 0.0;
  private double yyVirial = 0.0;
  private double zzVirial = 0.0;
  private double xyVirial = 0.0;
  private double xzVirial = 0.0;
  private double yzVirial = 0.0;

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

  private static final int NN = 10000;
  private int kind[] = new int[NN];
  private double xx[] = new double[NN];
  private double yy[] = new double[NN];
  private double zz[] = new double[NN];
  private double vx[] = new double[NN];
  private double vy[] = new double[NN];
  private double vz[] = new double[NN];
  private double ffx[] = new double[NN];
  private double ffy[] = new double[NN];
  private double ffz[] = new double[NN];

  private double electronDensityAt[][] = new double[NN][2];

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

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 0;
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int count = 0;


  // ------------------------------- class constructor -----------

  FixedGuptaMD(FixedGuptaAtom fga) {
    fg = fga;
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      count += 1;
      timeEvolution();
      if (sleepTime>0) {
        try {
          Thread.sleep(sleepTime);
        } catch (InterruptedException e) { }
      }
    }
  }


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

  private void setInitialCondition() {
    int i; 
    double a;

    count = 0;
    t = 0.0;

    i = 0;

    // base crystal plane
    if (fg.getXtalKindOf(MAT1)==0) {
      i = setFcc111YBase(i, MAT1);

    } else if (fg.getXtalKindOf(MAT1)==1) {
      i = setBcc110YBase(i, MAT1);

    } else if (fg.getXtalKindOf(MAT1)==2) {
      i = setHcp0001YBase(i, MAT1);

    }

    // adatom
    a = 0.6*(fg.getDiaOf(MAT1)+fg.getDiaOf(MAT2))*1.0e-10;
    i = setAdatom(i, MAT2, 0.5*yMax-a);

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

    removeMacroMotion();
    initVirial();
  }

  private void initVirial() {
    xxVirial = 0.0; yyVirial = 0.0; zzVirial = 0.0;
    xyVirial = 0.0; xzVirial = 0.0; yzVirial = 0.0;
  }

  // fcc(111) plane

  int setFcc111YBase(int ii, int knd) {
    int nx,ny,nz;
    double a,b,c,size;

    a = 1.01*0.5*fg.getDiaOf(knd)*1.0e-10;
    b = (Math.sqrt(3.0)/3.0)*a;
    c = (Math.sqrt(6.0)/3.0)*2.0*a;
    size = 5.5e-9;

    nx = 2*(int)(size/a/4.0);
    ny = (int)(0.48*yMax/c);
    nz = 2*(int)(size/b/6.0);

    return( setFcc111YPlane(ii, knd, nx, ny, nz, 0.5*yMax) );
  }

  int setFcc111YPlane(int ii, int knd, int nx, int ny, int nz, double yPos) {
    int ip,ix,iy,iz,nnx,nnz; 
    double a,b,c,rm;

    ip = ii;

    a = 1.01*0.5*fg.getDiaOf(knd)*1.0e-10;
    b = (Math.sqrt(3.0)/3.0)*a;
    c = (Math.sqrt(6.0)/3.0)*2.0*a;
    rm = 10.0/Math.sqrt(fg.getMassOf(knd)/fg.getAU());

    nnx = nx*2;
    nnz = nz*3;

    xMax = a*nnx;
    zMax = b*nnz;

    for (ix=0; ix<nnx; ix++) {
      for (iy=0; iy<ny; iy++) {
        for (iz=0; iz<nnz; iz++) {
          if (fcc111Q(ix,iy,iz)==1) {
            kind[ip] = knd;
            xx[ip] = a*(double)(ix)+0.5*a;
            yy[ip] = c*(double)(iy)+yPos;
            zz[ip] = b*(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 fcc111Q(int ix, int iy, int iz) {
    int iix, iiy, iiz;

    iix = ix%2; iiy = iy%3; iiz = iz%6;

    if (iiy==0 && iix==0 && iiz==0) return(1);
    if (iiy==0 && iix==1 && iiz==3) return(1);
    if (iiy==1 && iix==0 && iiz==4) return(1);
    if (iiy==1 && iix==1 && iiz==1) return(1);
    if (iiy==2 && iix==0 && iiz==2) return(1);
    if (iiy==2 && iix==1 && iiz==5) return(1);
	
    return(0);
  }

  // bcc(110) plane

  int setBcc110YBase(int ii, int knd) {
    int nx,ny,nz;
    double a,b,c,size;

    a = 1.1547*fg.getDiaOf(knd)*1.0e-10;
    b = 0.5*Math.sqrt(2.0)*a;
    c = 0.5*a;
    size = 5.5e-9;

    nx = 2*(int)(size/b/2.0);
    ny = (int)(0.48*yMax/b);
    nz = 2*(int)(size/c/2.0);
	
    return( setBcc110YPlane(ii, knd, nx, ny, nz, 0.5*yMax) );
  }

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

    ip = ii;

    a = 1.1547*fg.getDiaOf(knd)*1.0e-10;
    b = 0.5*Math.sqrt(2.0)*a;
    c = 0.5*a;
    rm = 10.0/Math.sqrt(fg.getMassOf(knd)/fg.getAU());


    xMax = 1.03*b*nnLayer;
    zMax = 1.02*c*nz;

    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] = 1.03*b*(double)(ix)+0.5*b;
            yy[ip] = 1.02*b*(double)(iy)+yPos;
            zz[ip] = 1.02*c*(double)(iz)+0.5*c;
            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 );
  }

  // hcp(0001) plane

  private int setHcp0001YBase(int ii, int knd) {
    int nx,ny,nz;
    double a,b,c,size;

    a = 1.0*fg.getDiaOf(knd)*1.0e-10;
    b = 0.8660*a;
    c = 0.8165*a;
    size = 5.5e-9;

    nx = 2*(int)(size/b/2.0);
    ny = (int)(0.48*yMax/c);
    nz = 2*(int)(size/a/2.0);

    return( setHcp0001YPlane(ii, knd, nx, ny, nz, 0.5*yMax) );
  }

  private int setHcp0001YPlane(int ii, int knd, int nx, int ny, int nz, double yPos) {
    int ip,ix,iy,iz;
    double a,b,c,dy,dz,rm;

    ip = ii;
    a = 1.0*fg.getDiaOf(knd)*1.0e-10;
    b = 0.8660*a;
    c = 0.8165*a;
    rm = 10.0/Math.sqrt(fg.getMassOf(knd)/fg.getAU());

    xMax = 1.00*b*nx;
    zMax = 1.00*a*nz;

    for (ix=0; ix<nx; ix++) {
      for (iy=0; iy<ny; iy++) {
        for (iz=0; iz<nz; iz++) {
          kind[ip] = knd;
          xx[ip] = b*ix+0.333333*b*(iy%2) + 0.25*b;
          yy[ip] = c*iy+yPos;
          zz[ip] = a*iz+0.5*a*(ix%2)+0.5*a*(iy%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 );
  }

  //

  private int setAdatom(int ii, int knd, double yPos) {
    int ip,ix,iy,iz,nx,ny,nz; 
    double a,b,c,rm;

    ip = ii;

    a = 1.2*fg.getDiaOf(knd)*1.0e-10;
    rm = 10.0/Math.sqrt(fg.getMassOf(knd)/fg.getAU());

    nx = (int)(xMax/a);
    ny = 1;
    nz = (int)(zMax/a);

    for (ix=0; ix<nx; ix++) {
      for (iy=0; iy<ny; iy++) {
        for (iz=0; iz<nz; iz++) {
          kind[ip] = knd;
          xx[ip] = a*(double)(ix)+0.5*(xMax-nx*a);
          yy[ip] = a*(double)(iy)+yPos;
          zz[ip] = a*(double)(iz)+0.5*(zMax-nz*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 );
  }

  //

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


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

  private void timeEvolution(){
    int i;

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

    if (startSW==1) {
      timeStep();

    } else {
      if (stepSW==1) {
        timeStep();
        stepSW = 0;
      }
      try {
        Thread.sleep(stopSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void timeStep() {
    int i;

    if (tempMode==1) {
      vAjustment();
    }
    regNearTB();
    for (i=0; i<10; i++) {
      timeStep1();
      if (statMode==1) setElectronDensity();
      volControl();
    }
  }

  private void timeStep1() {
    int i; 
    double dtv2, a2;

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

    for (i=0; i<Nmt; i++) {
      a2 = dtv2/fg.getMassOf(kind[i]);
      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/fg.getMassOf(kind[i]);
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
      vz[i] += a2*ffz[i];
    }

    setBoundary();
  }

  private void forceCalcTB() {
    int i,j,k;
    double rij,f,fxij,fyij,fzij;
    double xij,yij,zij, xxvir,yyvir,zzvir,xyvir,xzvir,yzvir;

    for(i=0;i<Nmt;i++) {
      ffx[i] = 0.0;
      ffy[i] = 0.0;
      ffz[i] = 0.0;
    }
    xxvir = 0.0; yyvir = 0.0; zzvir = 0.0;
    xyvir = 0.0; xzvir = 0.0; yzvir = 0.0;

    for(i=0;i<Nmt;i++) {
      for(k=1;k<regNr[i][0];k++) {
        j = regNr[i][k];
        xij = xx[i]-xx[j]; yij = yy[i]-yy[j]; zij = zz[i]-zz[j];
        if (xij>0.5*xMax) xij -= xMax;
        if (xij<-0.5*xMax) xij += xMax;
        //if (yij>0.5*yMax) yij -= yMax;
        //if (yij<-0.5*yMax) yij += yMax;
        if (zij>0.5*zMax) zij -= zMax;
        if (zij<-0.5*zMax) zij += zMax;
        rij = Math.sqrt(xij*xij+yij*yij+zij*zij);
        f = fg.force(rij,kind[i],kind[j]);
        fxij = f*xij/rij;
        fyij = f*yij/rij;
        fzij = f*zij/rij;

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

        xxvir +=  fxij*xij;
        yyvir +=  fyij*yij;
        zzvir +=  fzij*zij;
        xyvir +=  fxij*yij;
        xzvir +=  fxij*zij;
        yzvir +=  fyij*zij;
      }
    }
    xxVirial = 0.995*xxVirial + 0.005*xxvir;
    yyVirial = 0.995*yyVirial + 0.005*yyvir;
    zzVirial = 0.995*zzVirial + 0.005*zzvir;
    xyVirial = 0.995*xyVirial + 0.005*xyvir;
    xzVirial = 0.995*xzVirial + 0.005*xzvir;
    yzVirial = 0.995*yzVirial + 0.005*yzvir;
  }

  private void setBoundary() {
    int i; 
    double rr, tmp;

    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++) {
      if (xx[i] < 0.0) { 
        xx[i] += xMax; 
      }
      if (xx[i] > xMax) {
        xx[i] -= xMax; 
      }
      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] += zMax; 
      }
      if (zz[i] > zMax) {
        zz[i] -= zMax; 
      }
    }
  }

  // registration of near atoms

  private 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 xij,yij,zij;

    preRegistration();

    rcutoff = fg.getDiaOf(MAT1)*1.732; 
    if (fg.getXtalKindOf(MAT1)==1) rcutoff = fg.getDiaOf(MAT1)*1.6;
    if (fg.getXtalKindOf(MAT1)==2) rcutoff = fg.getDiaOf(MAT1)*2.2362;

    rcutoff2 = fg.getDiaOf(MAT2)*1.732; 
    if (fg.getXtalKindOf(MAT2)==1) rcutoff2 = fg.getDiaOf(MAT2)*1.6;
    if (fg.getXtalKindOf(MAT2)==2) rcutoff2 = fg.getDiaOf(MAT2)*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); 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+(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;
          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) {
                xij = xx[ip]-xx[jp]; yij = yy[ip]-yy[jp]; zij = zz[ip]-zz[jp];
                if (xij>0.5*xMax) xij -= xMax;
                if (xij<-0.5*xMax) xij += xMax;
                //if (yij>0.5*yMax) yij -= yMax;
                //if (yij<-0.5*yMax) yij += yMax;
                if (zij>0.5*zMax) zij -= zMax;
                if (zij<-0.5*zMax) zij += zMax;
                r2 = xij*xij+yij*yij+zij*zij;
                if (r2<rrg2) {
                  regNr[ip][kp]=jp;
                  kp += 1;
                }
              }
            }
          }
        }
      }
      regNr[ip][0]=kp;
    }
  }

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

  }


  // ----------------------------------------- utility -----------

  private 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

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

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      nkt += fg.getMassOf(kind[i])*vx[i]*vx[i];
    }

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

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

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      nkt += fg.getMassOf(kind[i])*vy[i]*vy[i];
    }

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

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

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      nkt += fg.getMassOf(kind[i])*vz[i]*vz[i];
    }

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

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

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      nkt += fg.getMassOf(kind[i])*vx[i]*vy[i];
    }

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

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

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {
      nkt += fg.getMassOf(kind[i])*vx[i]*vz[i];
    }

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

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

    nkt= 0.0;
    for (i=0; i<Nmt; i++) {

      nkt += fg.getMassOf(kind[i])*vy[i]*vz[i];
    }

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

  // vol control

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

  //

  private void setElectronDensity() {
    int i,j,k;
    double xij,yij,zij,rij,fij;

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

    for(i=0;i<Nmt;i++) {
      for(k=1;k<regNr[i][0];k++) {
        j = regNr[i][k];
        xij = xx[i]-xx[j]; yij = yy[i]-yy[j]; zij = zz[i]-zz[j];
        if (xij>0.5*xMax) xij -= xMax;
        if (xij<-0.5*xMax) xij += xMax;
        //if (yij>0.5*yMax) yij -= yMax;
        //if (yij<-0.5*yMax) yij += yMax;
        if (zij>0.5*zMax) zij -= zMax;
        if (zij<-0.5*zMax) zij += zMax;
        rij = Math.sqrt(xij*xij+yij*yij+zij*zij);
        fij = fg.faiInter(rij*1.0e10, kind[i], kind[j]);
        electronDensityAt[i][1] += fij;
        electronDensityAt[j][1] += fij;
      }
    }

    for(i=0;i<Nmt;i++) {
      electronDensityAt[i][0] = electronDensityAt[i][1];
    }
  }

  //

  private double temperature() {
    int i;
    double s;

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

  private double temp(int ii) {
    double s;

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

  private 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*fg.getMassOf(kind[i])*((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 );
  }

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

  public void rdfSet(int rdf[][]) {
    int nn,i,j,k,ir;
    double 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];
        rij = distance(i,j);
        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;
        }
      }
    }
  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public int getCount() { return count; }

  public static int getNmt() { return Nmt; }
  public static double getxMax() { return xMax; }
  public static double getyMax() { return yMax; }
  public static double getzMax() { return zMax; }
  public void setMAT1(int m) { MAT1 = m; }
  public static int getMAT1() { return MAT1; }
  public void setMAT2(int m) { MAT2 = m; }
  public static int getMAT2() { return MAT2; }

  public void setStatMode(int m) { statMode = m; }
  public void setTempMode(int tmode) { tempMode = tmode; }
  public int getTempMode() { return tempMode; }
  public void setContTemp(double tmp) {
    contTemp = tmp;
    vAjustment();
  }
  public double getContTemp() { return contTemp; }

  // get data
  public double getTime() { return t; }
  public double getdt() { return dt; }

  public float getColorOf(int knd) { return (float)fg.getColorOf(knd); }
  public double getMassOf(int knd) { return fg.getMassOf(knd); }
  public double getEffectiveInter(double ra, int ki, int kj) {
    return fg.effectiveInter(ra,ki,kj); }
  public double getPairInter(double ra, int ki, int kj) {
    return fg.pairInter(ra,ki,kj); }

  public double getDia(int i) { return fg.getDiaOf(kind[i]); }
  public double getColor(int i) { return fg.getColorOf(kind[i]); }
  public double getggRef(int i) { return fg.getggRefOf(kind[i]); }

  public int getKind(int i) { return kind[i]; }
  public double getxx(int i) { return xx[i]; }
  public double getyy(int i) { return yy[i]; }
  public double getzz(int i) { return zz[i]; }
  public double getvx(int i) { return vx[i]; }
  public double getvy(int i) { return vy[i]; }
  public double getvz(int i) { return vz[i]; }
  public double getTemperature(int i) { return temp(i); }
  public double getElectronDensityAt(int i) { return electronDensityAt[i][0]; }

  public int getnnRegNr(int i) { return regNr[i][0]; }
  public int getRegNr(int i, int k) { return regNr[i][k]; }

  public double getTemperature() { return temperature(); }
  public double getxxPressure() { return xxPress(); }
  public double getyyPressure() { return yyPress(); }
  public double getzzPressure() { return zzPress(); }
  public double getxyPressure() { return xyPress(); }
  public double getxzPressure() { return xzPress(); }
  public double getyzPressure() { return yzPress(); }

}


/* fixed G (electron density = reference crystal one) Gupta potential
 *
 * Tight-binding second-moment approximation
 *
 * 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
 *
 */