/** applet No. 1228
 *
 * sintering nickel(fcc) - molecular dynamics 3D
 * - selected fixed Gupta potential (SFG)
 * - multi thread  SFGMD :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on October 14 2007.
 * Copyright (c) 2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2007.10.14  created
 *
 */

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

public class sinteringNiSFG 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(2); // Ni(fcc)=2
  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 = 2;
  int MAT2 = 2;

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

  // 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(new Label("  "));
    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.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 = 0;
  private double xMax = md.getxMax();
  private double yMax = md.getyMax();
  private double zMax = md.getzMax();

  private int xImageLoc = 60;
  private int yImageLoc = 130;

  private int xPlot2dLoc = 400;
  private int yPlot2dLoc = 140;

  private double plotSize = 250.0;
  private double dispScale = plotSize/xMax;
  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 = md.getNN();
  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 electronDensityAt[] = new double[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 String xtals[] = {"fcc","bcc","hcp"};

  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();
    xMax = md.getxMax();
    yMax = md.getyMax();
    zMax = md.getzMax();
    MAT1 = md.getMAT1(); MAT2 = md.getMAT2();

    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));
    gOff.drawString(fg.getMatStr(MAT1)+" ("+xtals[fg.getXtalKindOf(MAT1)]+")",630/6*0+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.drawString("P="+(int)(md.getPressure()/1.0e6)+"MPa",630/6*4+10,90);

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

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

    gOff.setColor(Color.black);
    gOff.drawString("regMax="+md.regMax()+"",400,130);
    gOff.drawString("secMax="+md.secMax()+"",500,130);
  }


  // ------------------------------------ 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.setEDensity(electronDensityAt);

    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.getTemperatureOf(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 = electronDensityAt[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 = (electronDensityAt[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.getColorOfKind(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.getColorOfKind(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.getColorOfKind(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.getColorOfKind(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.getColorOfKind(MAT1),0.8f,0.6f));
    gOff.drawString("mean distr.",gx+120,gy+40);
    //gOff.setColor(Color.getHSBColor(md.getColorOfKind(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.getMassOfKind(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 ===========
//
// selected fixed Gupta potential
//

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][];
  private double [][][] faiTBL = new double[kindMax][kindMax][];
  private double [][] ggRef = new double[kindMax][kindMax];


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

  FixedGuptaAtom() {
    int i,j;

    for (i=0; i<kindMax; i++) {
      for (j=0; j<kindMax; j++) {
        malocTBL(i, j);
      }
    }

    for (i=0; i<kindMax; i++) {
      setTable(i);
    }
    for (i=0; i<kindMax; i++) {
      for (j=0; j<kindMax; j++) {
        if (j>i) {
          setTable(i, j);
        }
      }
    }
  }

  FixedGuptaAtom(int knd) {
    malocTBL(knd, knd);
    setTable(knd);
  }

  FixedGuptaAtom(int ki,int kj) {
    malocTBL(ki, ki);
    malocTBL(kj, kj);
    malocTBL(ki, kj);
    malocTBL(kj, ki);
    setTable(ki);
    setTable(kj);
    setTable(ki,kj);
    setTable(kj,ki);
  }

  FixedGuptaAtom(int ki,int kj,int kk) {
    malocTBL(ki, ki);
    malocTBL(kj, kj);
    malocTBL(kk, kk);
    malocTBL(ki, kj);
    malocTBL(kj, ki);
    malocTBL(ki, kk);
    malocTBL(kk, ki);
    malocTBL(kj, kk);
    malocTBL(kk, kj);
    setTable(ki);
    setTable(kj);
    setTable(kk);
    setTable(ki,kj);
    setTable(kj,ki);
    setTable(ki,kk);
    setTable(kk,ki);
    setTable(kj,kk);
    setTable(kk,kj);
  }

  private void malocTBL(int ki, int kj) {
    pairTBL[ki][kj] = new double[1000];
    faiTBL[ki][kj] = new double[1000];
  }

  //

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

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

  private void setTable(int i, int j) {
    int ir;
    double ra;

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


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

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

  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 = 20000;
  private static double xMax = 6.4*1.0e-9;
  private static double yMax = 6.4*1.0e-9;
  private static double zMax = 6.4*1.0e-9;
  private static int MAT1 = 2; // Ni
  private static int MAT2 = 2; // Ni

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

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

  private double xVirial = 0.0;
  private double yVirial = 0.0;
  private double zVirial = 0.0;

  private int Nsx = 24;
  private int Nsy = 24;
  private int Nsz = 24;

  private static final int NN = 20000;
  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 int regNr[][] = new int[NN][100];
  private int section[][][][] = new int[Nsx][Nsy][Nsz][50];

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

    t = 0.0;
    count = 0;

    nm = 1.0e-9;
    i = 0;
    i = setNanoParticle(i,MAT1, 1.7*nm, 1.7*nm, 3.2*nm, 1.4*nm, 0.0*Math.PI/4.0, 1.0*Math.PI/2.0);
    i = setNanoParticle(i,MAT1, 1.7*nm, 4.7*nm, 3.2*nm, 1.4*nm, 1.0*Math.PI/4.0, 0.0);
    i = setNanoParticle(i,MAT1, 4.7*nm, 1.7*nm, 3.2*nm, 1.4*nm, 2.0*Math.PI/4.0, 0.0);
    i = setNanoParticle(i,MAT1, 4.7*nm, 4.7*nm, 3.2*nm, 1.4*nm, 3.0*Math.PI/4.0, 1.0*Math.PI/2.0);
    Nmt = i;

    removeMacroMotion();
    initVirial();
  }

  private void initVirial() {
    xVirial = 0.0; yVirial = 0.0; zVirial = 0.0;
  }

  private int setNanoCube(int ii, int knd, double x0, double y0, double z0, double cLength, double th, double fi) {
    int xtal;

    xtal = fg.getXtalKindOf(knd);
	if (xtal==0) { // fcc
	  return setFCCNanoCube(ii, knd, x0, y0, z0, cLength, th, fi);
	} else if  (xtal==1) { // bcc
	  return setBCCNanoCube(ii, knd, x0, y0, z0, cLength, th, fi);
	} else if  (xtal==2) { // hcp
	  return setHCPNanoCube(ii, knd, x0, y0, z0, cLength, th, fi);
	}	
    return -1;
  }

  private int setNanoParticle(int ii, int knd, double x0, double y0, double z0, double radius, double th, double fi) {
    int xtal;

    xtal = fg.getXtalKindOf(knd);
	if (xtal==0) { // fcc
	  return setFCCNanoParticle(ii, knd, x0, y0, z0, radius, th, fi);
	} else if  (xtal==1) { // bcc
	  return setBCCNanoParticle(ii, knd, x0, y0, z0, radius, th, fi);
	} else if  (xtal==2) { // hcp
	  return setHCPNanoParticle(ii, knd, x0, y0, z0, radius, th, fi);
	}	
    return -1;
  }

  // fcc

  private int setFCCNanoCube(int ii, int knd, double x0, double y0, double z0, double cLength, double th, double fi) {
    int ip,i,j,k,n;
    double a,len, x,y,z, xp,yp,zp, sth,cth,sfi,cfi;

    ip = ii;

    a = 1.0/Math.sqrt(2.0)*fg.getDiaOf(knd)*1.0e-10;
    sth = Math.sin(th); cth = Math.cos(th);
    sfi = Math.sin(fi); cfi = Math.cos(fi);

    len = 1.8*cLength;
    n = (int)(len/a) + 1;

    for (i=0; i<n; i++) {
      x = a*i - len/2.0;
      for (j=0; j<n; j++) {
        y = a*j - len/2.0;
        for (k=0; k<n; k++) {
          z = a*k - len/2.0;
          if ( (i+j+k)%2==0 ) {
            xp = x0 + cLength/2.0 + cth*cfi*x + sth*cfi*y + sfi*z;
            yp = y0 + cLength/2.0 - sth*x + cth*y;
            zp = z0 + cLength/2.0 - cth*sfi*x - sth*sfi*y + cfi*z;
            if (xp>=x0 && xp<=x0+cLength && yp>=y0 && yp<=y0+cLength && zp>=z0 && zp<=z0+cLength) {
              kind[ip] = knd;
              xx[ip] = xp;
              yy[ip] = yp;
              zz[ip] = zp;
              vx[ip] = 500.0*(Math.random()-0.5);
              vy[ip] = 500.0*(Math.random()-0.5);
              vz[ip] = 500.0*(Math.random()-0.5);
              ffy[ip] = 0.0;
              ffy[ip] = 0.0;
              ffz[ip] = 0.0;
              ip += 1;
            }
          }
        }
      }
    }
    return ip;
  }

  private int setFCCNanoParticle(int ii, int knd, double x0, double y0, double z0, double radius, double th, double fi) {
    int ip,i,j,k,n;
    double a,len, x,y,z, xp,yp,zp, sth,cth,sfi,cfi;

    ip = ii;

    a = 1.0/Math.sqrt(2.0)*fg.getDiaOf(knd)*1.0e-10;
    sth = Math.sin(th); cth = Math.cos(th);
    sfi = Math.sin(fi); cfi = Math.cos(fi);

    len = 2.0*radius;
    n = (int)(len/a) + 1;

    for (i=0; i<n; i++) {
      x = a*i - len/2.0;
      for (j=0; j<n; j++) {
        y = a*j - len/2.0;
        for (k=0; k<n; k++) {
          z = a*k - len/2.0;
          if ( (i+j+k)%2==0 ) {
            if (x*x+y*y+z*z<=radius*radius) {
              xp = cth*cfi*x + sth*cfi*y + sfi*z;
              yp = - sth*x + cth*y;
              zp = - cth*sfi*x - sth*sfi*y + cfi*z;
              kind[ip] = knd;
              xx[ip] = x0 + xp;
              yy[ip] = y0 + yp;
              zz[ip] = z0 + zp;
              vx[ip] = 500.0*(Math.random()-0.5);
              vy[ip] = 500.0*(Math.random()-0.5);
              vz[ip] = 500.0*(Math.random()-0.5);
              ffy[ip] = 0.0;
              ffy[ip] = 0.0;
              ffz[ip] = 0.0;
              ip += 1;
            }
          }
        }
      }
    }
    return ip;
  }

  // bcc

  private int setBCCNanoCube(int ii, int knd, double x0, double y0, double z0, double cLength, double th, double fi) {
    int ip,i,j,k,n;
    double a,len, x,y,z, xp,yp,zp, sth,cth,sfi,cfi;

    ip = ii;

    a = 1.0/Math.sqrt(3.0)*fg.getDiaOf(knd)*1.0e-10;
    sth = Math.sin(th); cth = Math.cos(th);
    sfi = Math.sin(fi); cfi = Math.cos(fi);

    len = 1.8*cLength;
    n = (int)(len/a) + 1;

    for (i=0; i<n; i++) {
      x = a*i - len/2.0;
      for (j=0; j<n; j++) {
        y = a*j - len/2.0;
        for (k=0; k<n; k++) {
          z = a*k - len/2.0;
          if ( (i%2==0 && j%2==0 && k%2==0) || (i%2==1 && j%2==1 && k%2==1) ) {

            xp = x0 + cLength/2.0 + cth*cfi*x + sth*cfi*y + sfi*z;
            yp = y0 + cLength/2.0 - sth*x + cth*y;
            zp = z0 + cLength/2.0 - cth*sfi*x - sth*sfi*y + cfi*z;
            if (xp>=x0 && xp<=x0+cLength && yp>=y0 && yp<=y0+cLength && zp>=z0 && zp<=z0+cLength) {
              kind[ip] = knd;
              xx[ip] = xp;
              yy[ip] = yp;
              zz[ip] = zp;
              vx[ip] = 500.0*(Math.random()-0.5);
              vy[ip] = 500.0*(Math.random()-0.5);
              vz[ip] = 500.0*(Math.random()-0.5);
              ffy[ip] = 0.0;
              ffy[ip] = 0.0;
              ffz[ip] = 0.0;
              ip += 1;
            }
          }
        }
      }
    }
    return ip;
  }

  private int setBCCNanoParticle(int ii, int knd, double x0, double y0, double z0, double radius, double th, double fi) {
    int ip,i,j,k,n;
    double a,len, x,y,z, xp,yp,zp, sth,cth,sfi,cfi;

    ip = ii;

    a = 1.0/Math.sqrt(3.0)*fg.getDiaOf(knd)*1.0e-10;
    sth = Math.sin(th); cth = Math.cos(th);
    sfi = Math.sin(fi); cfi = Math.cos(fi);

    len = 2.0*radius;
    n = (int)(len/a) + 1;

    for (i=0; i<n; i++) {
      x = a*i - len/2.0;
      for (j=0; j<n; j++) {
        y = a*j - len/2.0;
        for (k=0; k<n; k++) {
          z = a*k - len/2.0;
          if ( (i%2==0 && j%2==0 && k%2==0) || (i%2==1 && j%2==1 && k%2==1) ) {
            if (x*x+y*y+z*z<=radius*radius) {
              xp = cth*cfi*x + sth*cfi*y + sfi*z;
              yp = - sth*x + cth*y;
              zp = - cth*sfi*x - sth*sfi*y + cfi*z;
              kind[ip] = knd;
              xx[ip] = x0 + xp;
              yy[ip] = y0 + yp;
              zz[ip] = z0 + zp;
              vx[ip] = 500.0*(Math.random()-0.5);
              vy[ip] = 500.0*(Math.random()-0.5);
              vz[ip] = 500.0*(Math.random()-0.5);
              ffy[ip] = 0.0;
              ffy[ip] = 0.0;
              ffz[ip] = 0.0;
              ip += 1;
            }
          }
        }
      }
    }
    return ip;
  }

  // hcp

  private int setHCPNanoCube(int ii, int knd, double x0, double y0, double z0, double cLength, double th, double fi) {
    int ip,i,j,k,nx,ny,nz;
    double a,b,c,len, x,y,z, xp,yp,zp, sth,cth,sfi,cfi;

    ip = ii;

    a = 1.0*fg.getDiaOf(knd)*1.0e-10;
    b = a*Math.sqrt(3.0)/2.0;
    c = a*Math.sqrt(2.0/3.0);
    sth = Math.sin(th); cth = Math.cos(th);
    sfi = Math.sin(fi); cfi = Math.cos(fi);

    len = 1.8*cLength;
    nx = (int)(len/a) + 1;
    ny = (int)(len/b) + 1;
    nz = (int)(len/c) + 1;

    for (i=0; i<nx; i++) {
      for (j=0; j<ny; j++) {
        for (k=0; k<nz; k++) {
          x = a*i + 0.5*a*(j%2) - len/2.0;
          y = b*j + 0.57735*a*(k%2) - len/2.0;
          z = c*k - len/2.0;

          xp = x0 + cLength/2.0 + cth*cfi*x + sth*cfi*y + sfi*z;
          yp = y0 + cLength/2.0 - sth*x + cth*y;
          zp = z0 + cLength/2.0 - cth*sfi*x - sth*sfi*y + cfi*z;
          if (xp>=x0 && xp<=x0+cLength && yp>=y0 && yp<=y0+cLength && zp>=z0 && zp<=z0+cLength) {
            kind[ip] = knd;
            xx[ip] = xp;
            yy[ip] = yp;
            zz[ip] = zp;
            vx[ip] = 500.0*(Math.random()-0.5);
            vy[ip] = 500.0*(Math.random()-0.5);
            vz[ip] = 500.0*(Math.random()-0.5);
            ffy[ip] = 0.0;
            ffy[ip] = 0.0;
            ffz[ip] = 0.0;
            ip += 1;
          }
        }
      }
    }
    return ip;
  }

  private int setHCPNanoParticle(int ii, int knd, double x0, double y0, double z0, double radius, double th, double fi) {
    int ip,i,j,k,nx,ny,nz;
    double a,b,c,len, x,y,z, xp,yp,zp, sth,cth,sfi,cfi;

    ip = ii;

    a = 1.0*fg.getDiaOf(knd)*1.0e-10;
    b = a*Math.sqrt(3.0)/2.0;
    c = a*Math.sqrt(2.0/3.0);
    sth = Math.sin(th); cth = Math.cos(th);
    sfi = Math.sin(fi); cfi = Math.cos(fi);

    len = 2.0*radius;
    nx = (int)(len/a) + 1;
    ny = (int)(len/b) + 1;
    nz = (int)(len/c) + 1;

    for (i=0; i<nx; i++) {
      for (j=0; j<ny; j++) {
        for (k=0; k<nz; k++) {
          x = a*i + 0.5*a*(j%2) - len/2.0;
          y = b*j + 0.57735*a*(k%2) - len/2.0;
          z = c*k - len/2.0;
          if (x*x+y*y+z*z<=radius*radius) {
              xp = cth*cfi*x + sth*cfi*y + sfi*z;
              yp = - sth*x + cth*y;
              zp = - cth*sfi*x - sth*sfi*y + cfi*z;
              kind[ip] = knd;
              xx[ip] = x0 + xp;
              yy[ip] = y0 + yp;
              zz[ip] = z0 + zp;
              vx[ip] = 500.0*(Math.random()-0.5);
              vy[ip] = 500.0*(Math.random()-0.5);
              vz[ip] = 500.0*(Math.random()-0.5);
              ffy[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(){

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

  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, xvir,yvir,zvir;

    for(i=0;i<Nmt;i++) {
      ffx[i] = 0.0;
      ffy[i] = 0.0;
      ffz[i] = 0.0;
    }
    xvir = 0.0; yvir = 0.0; zvir = 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];

        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;

        xvir += fxij*xij;
        yvir += fyij*yij;
        zvir += fzij*zij;
      }
    }
    xVirial = 0.995*xVirial + 0.005*xvir;
    yVirial = 0.995*yVirial + 0.005*yvir;
    zVirial = 0.995*zVirial + 0.005*zvir;
  }

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

  // 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 += 1000.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) {
                xij = xx[ip]-xx[jp]; yij = yy[ip]-yy[jp]; zij = zz[ip]-zz[jp];
                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];
    }
  }

  public int regMax() {
    int i,m;

    m = 0;
    for(i=0;i<Nmt;i++) {
      if (regNr[i][0]>m) {
        m = regNr[i][0];
      }
    }
    return(m);
  }

  public int secMax() {
    int i,j,k,m;

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

  //

  private double pressure() {
    double vol;

    vol = xMax*yMax*zMax;
    return (Nmt*1.38e-23*temperature()+(xVirial+yVirial+zVirial)/3.0)/vol;
  }

  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 setEDensity(double eDens[]) {
    int i,j,k;
    double fij;

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

    for(i=0;i<Nmt;i++) {
      for(k=1;k<regNr[i][0];k++) {
        j = regNr[i][k];
        fij = fg.faiInter(distance(i,j)*1.0e10, kind[i], kind[j]);
        eDens[i] += fij;
        eDens[j] += fij;
      }
    }
  }

  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 getNN() { return NN; }
  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 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 getColorOfKind(int knd) { return (float)fg.getColorOf(knd); }
  public double getMassOfKind(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 getTemperatureOf(int i) { return temp(i); }

  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 getPressure() { return pressure(); }
}


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