/* Rosato potential periodic   */
/* coded by Ikeuchi Mitsuru    */
/* ver 0.0.1   2003.10.15      */
/* ver 0.0.2   2004.02.11 improved pannel layout */

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

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

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

  Choice cm1,ctc,cvw;
  Scrollbar scr,scs;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;
  double sgm = 3.4e-10;       /* L-J sigma for Ar (m)   */
  double dispScale = 20.0/sgm;

  double xMax = 8.0*sgm;       /* x-Box size in (m)      */
  double yMax = 8.0*sgm;       /* y-Box size in (m)      */
  double zMax = 8.0*sgm;       /* z-Box size in (m)      */


  int Nmt = 8*8*8/2;           /* number of molecules    */
  int NN = Nmt + 1;
  double t = 0.0;              /* time (s)               */
  double dt = 0.5*1.0e-15;     /* time division (s)      */

  int tempMode = 1;
  double contTemp = 10.0;
  double scale = 1.0;
  double scaleDens = 1.0;
  int dispMode = 0;


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

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

  double AA = 1.0e-10;
  double AU = 1.661e-27;
  double Rosato[][] = {
  /*  atom    0 mass      1 A     2 xi   3 p     4 q    5 a    6 Ec  7 cf */
  /* 0 Ni */ {  58.71*AU, 0.0376, 1.070, 16.999, 1.189, 3.523, 4.435, 0.98 },
  /* 1 Cu */ {  63.54*AU, 0.0855, 1.224, 10.960, 2.278, 3.615, 3.544, 0.90 },
  /* 2 Rh */ { 102.91*AU, 0.0629, 1.660, 18.450, 1.867, 3.803, 5.752, 0.97 },
  /* 3 Pd */ { 106.40*AU, 0.1746, 1.718, 10.867, 3.742, 3.887, 3.936, 0.88 },
  /* 4 Ag */ { 107.87*AU, 0.1028, 1.178, 10.928, 3.139, 4.085, 2.960, 0.90 },
  /* 5 Ir */ { 192.20*AU, 0.1156, 2.289, 16.980, 2.691, 3.839, 6.928, 0.95 },
  /* 6 Pt */ { 195.09*AU, 0.2975, 2.695, 10.612, 4.004, 3.924, 5.853, 0.88 },
  /* 7 Au */ { 196.97*AU, 0.2061, 1.790, 10.229, 4.036, 4.079, 3.779, 0.87 },
  /* 8 Pb */ { 207.19*AU, 0.0980, 0.914,  9.576, 3.648, 4.951, 2.043, 0.85 },
  /* 9 Al */ {  26.98*AU, 0.1221, 1.316,  8.612, 2.516, 4.050, 3.339, 0.86 }
  };
  int MAT1 = 0; /* Ni */

  int kind[] = new int[NN];       /* i-th molec kind */

  double xx[] = new double[NN];   /* i-th x-position */
  double yy[] = new double[NN];   /* i-th y-position */
  double zz[] = new double[NN];   /* i-th z-position */
  double vx[] = new double[NN];   /* i-th x-velocity */
  double vy[] = new double[NN];   /* i-th y-velocity */
  double vz[] = new double[NN];   /* i-th z-velocity */
  double ffx[] = new double[NN];  /* i-th x-force    */
  double ffy[] = new double[NN];  /* i-th y-force    */
  double ffz[] = new double[NN];  /* i-th z-force    */
  double fbx[] = new double[NN];  /* i-th x-force before dt */
  double fby[] = new double[NN];  /* i-th y-force before dt */
  double fbz[] = new double[NN];  /* i-th z-force before dt */


  int regNr[][] = new int[NN][1000];   /* registration */

  int reg[][] = new int[NN][1000];   /* registration */
  double rijReg[][] = new double[NN][1000];   /* registration */


  double mass = 68.71*1.661e-27;  /* Ni */
  double latticeConstant = 3.52e-10;
  double atomDia = 2.49; /* A */
  double rcut = 4.7895; /* A */

  double dr = 0.001; /* A */
  double droh = 0.01;

  double pair[] = new double[8000];
  double density[] = new double[8000];
  double embedded[] = new double[8000];

  int rdf[] = new int[120];


  int srtz[][] = new int[2][NN];   /* z-sort work space */

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

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

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

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

  int NFx = 28;
  int NFy = 28;
  int NFz = 28;
  double dens[][][] = new double[NFx][NFy][NFz];
  double adens[][][] = new double[13][13][13];

/*----------------------------- applet control functions -----*/

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

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

    cm1 = new Choice();
    cm1.addItem("Ni");  cm1.addItem("Cu"); cm1.addItem("Rh"); cm1.addItem("Pd");
    cm1.addItem("Ag"); cm1.addItem("Ir"); cm1.addItem("Pt"); cm1.addItem("Au");
    cm1.addItem("Pb"); cm1.addItem("Al");
    cm1.addItemListener(this);

    ctc = new Choice();
    ctc.addItem("adiabatic"); ctc.addItem("T scaling");
    ctc.addItemListener(this);
    ctc.select("T scaling");

    cvw = new Choice();
    cvw.addItem("ball"); cvw.addItem("line");
    cvw.addItem("both"); cvw.addItem("real");
    cvw.addItem("T color"); cvw.addItem("e density");
    cvw.addItemListener(this); 

    scr= new Scrollbar(Scrollbar.HORIZONTAL,1,10,1,110);
    scr.addAdjustmentListener(this);

    scs= new Scrollbar(Scrollbar.HORIZONTAL,100,10,20,210);
    scs.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(cm1);
    pnl.add(ctc);
    pnl.add(cvw);
    pnl.add(new Label("  "));
    pnl.add(scr);
    pnl.add(scs);
    add(pnl,"North");

    setEAM();
    setInitialPosition();
    setadens();
    setWaku();
  }

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

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

  public void itemStateChanged(ItemEvent ev){
     if (ev.getSource() == cm1){
      MAT1 = cm1.getSelectedIndex();
      t = 0.0;
      setInitialPosition();
      offPaint();
      repaint(); 
    } else if (ev.getSource() == ctc){
      tempMode = ctc.getSelectedIndex();
      offPaint();
      repaint(); 
    } else if (ev.getSource() == cvw){
      dispMode = cvw.getSelectedIndex();
      offPaint();
      repaint(); 
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if(ev.getSource() == scr) { 
      contTemp = 10.0*(double)(scr.getValue());
      vAjustment();
    } else if(ev.getSource() == scs) { 
      scale = 0.01*(double)(scs.getValue());
    }
  }

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

  public void mouseMoved  (MouseEvent ev){ }

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

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

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

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


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

  void offPaint() {
    int i,j,k,jj,farPoint;
    int gx,gy,gz,g2x,g2y,gbx,gby,gb;
    double sc,sz,tm,bl,tmp,clr;
    double rx,ry,rz;
    double AA = 1.0e-10;

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

    rotate(theta,fai);
    zSort();
    tmp = temperature();
    setRadialDistribution();

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

    gbx = 70;
    gby = 90;
    sc = dispScale*scale;
    sz = 0.5*atomDia*1.0e-10*sc;

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

    if (dispMode==0) {

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

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

    } else if (dispMode==1) {

      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
        gx = (int)(sc*px[j])+gbx;
        gy = (int)(sc*py[j])+gby;
        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
          if (distance(j,jj)<atomDia*1.2*1.0e-10) {
            g2x = (int)(sc*px[jj])+gbx;
            g2y = (int)(sc*py[jj])+gby;
            gOff.setColor(Color.gray);
            gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
      }

    } else if (dispMode==2) {

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

        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
          if (distance(j,jj)<atomDia*1.2*1.0e-10) {
            g2x = (int)(sc*px[jj])+gbx;
            g2y = (int)(sc*py[jj])+gby;
            gOff.setColor(Color.gray);
            gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
        gOff.setColor(Color.getHSBColor(0.33f,0.9f,(float)(0.4*pz[j]/zMax+0.4)));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
      }

    } else if (dispMode==3) {

      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
        gx = (int)(sc*px[j])+gbx;
        gy = (int)(sc*py[j])+gby;
        gz = (int)(2.0*sz*(0.4*pz[j]/zMax+0.6));

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

    } else if (dispMode==4) {

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

        tm = 0.5*mass*(vx[j]*vx[j]+vy[j]*vy[j]+vz[j]*vz[j])/(1.5*1.38e-23);
        tm = 0.66-tm/1000.0;
        if (tm<0.01) {
          tm = 0.01;
        }
        gOff.setColor(Color.getHSBColor((float)(tm),0.9f,0.7f));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
      }

    } else if (dispMode==5) {

      for (i=0; i<NFx; i++) {
        for (j=0; j<NFy; j++) {
          for (k=0; k<NFz; k++) {
            if (dens[i][j][k]>0.1) {

              ry = cosFi*((j)*AA-cy)+sinFi*((k)*AA-cz) + cy; 
              rz = -sinFi*((j)*AA-cy)+cosFi*((k)*AA-cz) + cz;
              rx = cosTh*((i)*AA-cx)+sinTh*(rz-cz) + cx; 
              rz = -sinTh*((i)*AA-cx)+cosTh*(rz-cz) + cz;

              gx = (int)(sc*rx)+gbx;
              gy = (int)(sc*ry)+gby;
              gz = (int)(0.2*(0.4*(scale*rz/zMax)+0.6)*scaleDens*dens[i][j][k]);
              if (gz>5) {
                gz = 5;
              }
              gOff.setColor(Color.getHSBColor(
               (float)(0.66-0.66*0.01*dens[i][j][k]),0.9f,0.9f));
              gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
            }
          }
        }
      }

    }


    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+gbx;
      gy = (int)(sc*pwky[boxp[i][0]])+gby;
      g2x = (int)(sc*pwkx[boxp[i][1]])+gbx;
      g2y = (int)(sc*pwky[boxp[i][1]])+gby;
      if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
        gOff.setColor(Color.gray);
        if (tempMode==2) {
          if (tmp>contTemp) {
            gOff.setColor(Color.black); 
          } else {
            gOff.setColor(Color.red);
          }
        }
        gOff.drawLine(gx,gy, g2x, g2y);
      }
    }


    gbx = 400;
    gby = 60;
    gOff.setColor(Color.gray);
    for (i=0; i<100; i+=20) {
      gOff.drawLine(gbx+i,gby, gbx+i, gby+200);
    }
    gOff.drawLine(gbx,gby+100, gbx+100, gby+100);
    gOff.setColor(Color.black);
    gOff.drawRect(gbx,gby,100,200);
    gOff.drawString("0",gbx-5,gby+215);
    gOff.drawString("5 A",gbx+90,gby+215);
    gOff.setColor(Color.cyan);
    for (i=0; i<50; i++) {
      if (rdf[i]>0) {
        gOff.fillRect(gbx+i*2,gby+200-rdf[i]/4, 2, rdf[i]/4);
      }
    }

    gOff.setColor(Color.black);
    gOff.drawString("atom",630/6*0+10,40);
    gOff.drawString("temp control",630/6*1+10,40);
    gOff.drawString("view",630/6*2+10,40);
    gOff.drawString("cont T="+contTemp+"K",630/6*4+10,40);
    gOff.drawString("scale="+(int)(scale*100.0+0.499)+"%",630/6*5+10,40);
    gbx = 630/6*5;
    gOff.setColor(Color.blue);
    gOff.drawString("t="+Math.floor(t*1.0e15+0.5)+"(fs)",gbx,80);
    gOff.drawString("T="+tmp+" K",gbx,100);
    gOff.drawString("Tr="+randomTemperature()+" K",gbx,120);


  }

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


/*----------------------------------  findMin()  ----------*/

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

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

/*-------------------------------------  sort z  ----------*/

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

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      srtz[1][i] = (int)(32767.0*pz[j]/zMax);
    } 
    for (i=0; i<100; i++) {
      for (j=0; j<Nmt-1; j++) {
        if (srtz[1][j]>srtz[1][j+1]) {
          w = srtz[0][j]; srtz[0][j]=srtz[0][j+1]; srtz[0][j+1]=w;
          w = srtz[1][j]; srtz[1][j]=srtz[1][j+1]; srtz[1][j+1]=w;
        }
      }
    }    
  }

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

    for (i=0; i<Nmt; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zMax);
    } 
    for (i=0; i<Nmt; i++) {
      for (j=0; j<Nmt-1; j++) {
        if (srtz[1][j]>srtz[1][j+1]) {
          w = srtz[0][j]; srtz[0][j]=srtz[0][j+1]; srtz[0][j+1]=w;
          w = srtz[1][j]; srtz[1][j]=srtz[1][j+1]; srtz[1][j+1]=w;
        }
      }
    }    
  }

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

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

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

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

/*-------------------------------------  setWaku  ---------*/

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

    i=0;

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


/*-------------------------------------  setadens  ---------*/

  void setadens() {
    int ix,iy,iz;
    double rA;

    for (ix=0; ix<13; ix++) {
      for (iy=0; iy<13; iy++) {
        for (iz=0; iz<13; iz++) {
          rA = Math.sqrt((double)((ix-6)*(ix-6)+(iy-6)*(iy-6)+(iz-6)*(iz-6)));
          if (rA<rcut) {
            adens[ix][iy][iz] = density[(int)(rA*1000)];
          } else {
            adens[ix][iy][iz] = 0.0;
          }
        }
      }
    }
  }


  void setdens() {
    int i,ix,iy,iz,bx,by,bz,px,py,pz;
    double rA;

    for (ix=0; ix<NFx; ix++) {
      for (iy=0; iy<NFy; iy++) {
        for (iz=0; iz<NFz; iz++) {
          dens[ix][iy][iz] = 0.0;
        }
      }
    }

    for (i=0; i<Nmt; i++) {
      bx = (int)(xx[i]*1.0e10);
      by = (int)(yy[i]*1.0e10);
      bz = (int)(zz[i]*1.0e10);

      for (ix=0; ix<13; ix++) {
        for (iy=0; iy<13; iy++) {
          for (iz=0; iz<13; iz++) {
            px = (bx+ix-6+NFx)%NFx;
            py = (by+iy-6+NFy)%NFy;
            pz = (bz+iz-6+NFz)%NFz;
            dens[px][py][pz] += adens[ix][iy][iz];
          }
        }
      }
    }
  }

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

  double temperature() {
    int i;
    double s;

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

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

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

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

    for (i=0; i<51; i++) {
      rdf[i]=0;
    }

    for (i=0; i<Nmt; i++) {
        for(k=1;k<reg[i][0];k++) {
        j = reg[i][k];
        rij=rijReg[i][k];
        if (rij<5.0e-10) {
          ir = (int)(Math.floor(rij*1.0e11+0.49999));
          rdf[ir] +=1;
        }
      }
    }
  }


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

  void vAjustment() {
    int i;
    double tmp, r;
    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*mass*((vx[i]-mvx)*(vx[i]-mvx) + (vy[i]-mvy)*(vy[i]-mvy) + (vz[i]-mvz)*(vz[i]-mvz));
    }
    tmp = s/((double)Nmt*1.5*1.38e-23);

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

/*-----------------------------  set EAM functions   ---------
*
*  pair[i]     : pair potential table  r = i*dr
*  density[i]  : electron density table  r = i*dr
*  embedded[i] : emdedded energy table  roh = i*droh
*
*--------------------------------------------------------------*/


  void setEAM() {
    int i;
    double x1, x2, rA, roh;

    x1 = 1.9149*0.707*Rosato[MAT1][5];
    x2 = 2.2362*0.707*Rosato[MAT1][5];

    for(i=0; i<8000; i++) {
      rA = dr*(double)(i);
      pair[i] = EAMpair(rA)*cutoff(x1,x2,rA);
      density[i] = EAMroh(rA)*cutoff(x1,x2,rA);
    }
    for(i=0; i<8000; i++) {
      roh = droh*(double)(i);
      embedded[i] = EAMff(roh);
    }
  }


  double EAMpair(double rA) {
    int i;
    double aa,p,r0;

    r0 = 0.707*Rosato[MAT1][5];
    aa = Rosato[MAT1][1];
    p = Rosato[MAT1][3];

    return( aa*Math.exp(-p*(rA/r0-1)) );
  }

  double EAMroh(double rA) {
    int i;
    double q,r0;

    r0 = 0.707*Rosato[MAT1][5];
    q = Rosato[MAT1][4];

    return( Math.exp(-2.0*q*(rA/r0-1)) );
  }

  double EAMff(double rh) {
    double xi;

    xi = Rosato[MAT1][2];

    return( -xi*Math.sqrt(rh) );
  }

  double cutoff(double x1, double x2, double x) {
    double r;

    r = 0.0;
    if (x < x1) {
      r = 1.0;
    } else if (x>=x1 && x<x2) {
      r = 0.5+0.5*Math.cos(2*3.1416*((x-x1)/(x2-x1)));
    }
    return ( r );
  }


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

  void setInitialPosition(){
    int i,j,ix,iy,iz; 
    double a;

    mass = Rosato[MAT1][0];
    latticeConstant = Rosato[MAT1][5];
    atomDia = 0.707*Rosato[MAT1][5];
    rcut = 2.236*atomDia;
    setEAM();

    a = 0.5*latticeConstant*1.0e-10*Rosato[MAT1][7];
    i = 0;
    for (ix=0; ix<8; ix++) {
      for (iy=0; iy<8; iy++) {
        for (iz=0; iz<8; iz++) {
          if ((ix+iy+iz)%2==0) {
            xx[i] = a*(double)(ix)+1.5*sgm;
            yy[i] = a*(double)(iy)+1.5*sgm;
            zz[i] = a*(double)(iz)+1.5*sgm;
            vx[i] = 100.0*Math.random()-50.0+300.0;
            vy[i] = 100.0*Math.random()-50.0;
            vz[i] = 100.0*Math.random()-50.0;
            ffx[i] = 0.0;
            ffy[i] = 0.0;
            ffz[i] = 0.0;
            fbx[i] = 0.0;
            fby[i] = 0.0;
            fbz[i] = 0.0;
            i += 1;
          }
        }
      }
    }

    zSortInit();
  }



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

  void calcposition(){
    int i;
    
    if (tempMode==1 || t<1000.0*dt) {
      vAjustment();
    }

    regNear();
    for (i=0; i<10; i++) {
      timeEvolution();
    }

    if (dispMode==5) {
      setdens();
    }
  }

  void timeEvolution(){
    int i; 
    double rr,tmp;

    t = t + dt;

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

    for (i=0; i<Nmt; i++) {
      xx[i] += vx[i]*dt+0.5*dt*dt*ffx[i]/mass;
      yy[i] += vy[i]*dt+0.5*dt*dt*ffy[i]/mass;
      zz[i] += vz[i]*dt+0.5*dt*dt*ffz[i]/mass;
      fbx[i] = ffx[i];
      fby[i] = ffy[i];
      fbz[i] = ffz[i];
    }
    forceCalc();
    for (i=0; i<Nmt; i++) {
      vx[i] += 0.5*dt*(ffx[i]+fbx[i])/mass;
      vy[i] += 0.5*dt*(ffy[i]+fby[i])/mass;
      vz[i] += 0.5*dt*(ffz[i]+fbz[i])/mass;
    }

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

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

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

    regProximity();
    for(i=0;i<Nmt;i++) {
      fpi = dFdroh(i);
      for(k=1;k<reg[i][0];k++) {
        j = reg[i][k];
        if (j>i) {
          xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[j]; zi = zz[i]; zj = zz[j];
          if (xi-xj>0.5*xMax) { xj += xMax; };
          if (xi-xj<-0.5*xMax) { xj -= xMax; };  
          if (yi-yj>0.5*yMax) { yj += yMax; };
          if (yi-yj<-0.5*yMax) { yj -= yMax; }; 
          if (zi-zj>0.5*zMax) { zj += zMax; };
          if (zi-zj<-0.5*zMax) { zj -= zMax; }; 
          rij = rijReg[i][k];
          f = force(rij,fpi);
          fxij = f*(xi-xj)/rij;
          fyij = f*(yi-yj)/rij;
          fzij = f*(zi-zj)/rij;

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

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

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

  void regNear() {
    int i,j,k;
    double xi,xj,yi,yj,zi,zj;
    double rrg,rrg2,r2;

    rrg = rcut*1.0e-10;
    rrg += 10000.0*dt*10;
    rrg2 = rrg*rrg;

    for(i=0;i<Nmt;i++) {
      k = 1;
      for(j=0;j<i;j++) {
        xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[j]; zi = zz[i]; zj = zz[j];
        if (xi-xj>0.5*xMax) { xj += xMax; };
        if (xi-xj<-0.5*xMax) { xj -= xMax; };  
        if (yi-yj>0.5*yMax) { yj += yMax; };
        if (yi-yj<-0.5*yMax) { yj -= yMax; }; 
        if (zi-zj>0.5*zMax) { zj += zMax; };
        if (zi-zj<-0.5*zMax) { zj -= zMax; }; 
        r2 = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
        if (r2<rrg2) {
          regNr[i][k]=j;
          k += 1;
        }
      }
      for(j=i+1;j<Nmt;j++) {
        xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[j]; zi = zz[i]; zj = zz[j];
        if (xi-xj>0.5*xMax) { xj += xMax; };
        if (xi-xj<-0.5*xMax) { xj -= xMax; };  
        if (yi-yj>0.5*yMax) { yj += yMax; };
        if (yi-yj<-0.5*yMax) { yj -= yMax; }; 
        if (zi-zj>0.5*zMax) { zj += zMax; };
        if (zi-zj<-0.5*zMax) { zj -= zMax; }; 
        r2 = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
        if (r2<rrg2) {
          regNr[i][k]=j;
          k += 1;
        }
      }
      regNr[i][0]=k;
    }

  }

/*--------------------  embedded atom method calculation   ---*/

  double force(double r,double fpi) { 
    double ra,rt;

    ra = r*1.0e10;
    rt = -dpairdr(ra)-2.0*fpi*ddensdr(ra);

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

  void regProximity() {
    int i,j,k,jj;
    double xi,xj,yi,yj,zi,zj;
    double rrg,rrg2,r2;

    rrg = rcut*1.0e-10;
    rrg2 = rrg*rrg;

    for(i=0;i<Nmt;i++) {
      k = 1;
      for(jj=1;jj<regNr[i][0];jj++) {
        j = regNr[i][jj];
        xi = xx[i]; xj = xx[j]; yi = yy[i]; yj = yy[j]; zi = zz[i]; zj = zz[j];
        if (xi-xj>0.5*xMax) { xj += xMax; };
        if (xi-xj<-0.5*xMax) { xj -= xMax; };  
        if (yi-yj>0.5*yMax) { yj += yMax; };
        if (yi-yj<-0.5*yMax) { yj -= yMax; }; 
        if (zi-zj>0.5*zMax) { zj += zMax; };
        if (zi-zj<-0.5*zMax) { zj -= zMax; }; 
        r2 = (xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
        if (r2<rrg2) {
          reg[i][k]=j;
          rijReg[i][k] = Math.sqrt(r2);
          k += 1;
        }
      }
      reg[i][0]=k;
    }
  }

  double dFdroh(int i) {
    int iroh;
    double rh,rt;

    rh = roh(i);
    iroh = (int)(rh/droh);
    rt = 0.0;
    if (rh>0.0 && rh<80.0) { 
      rt = (embedded[iroh+1] - embedded[iroh])/droh;
    }
    return ( rt );
  }


  double roh(int i) {
    int j,k;
    double s;

    s = 0.0;
    for(k=1;k<reg[i][0];k++) {
      j = reg[i][k];
      s += density[ (int)(rijReg[i][k]*1.0e10/dr) ];
    }
    return(s);
  }

  double dpairdr(double r) {
    int ir;

    if (r<rcut) {
      ir = (int)(r/dr);
      return ( (pair[ir+1] - pair[ir])/dr );
    } else {
      return ( 0.0 );
    }
  }

  double ddensdr(double r) {
    int ir;

    if (r<rcut) {
      ir = (int)(r/dr);
      return ( (density[ir+1] - density[ir])/dr );
    } else {
      return ( 0.0 );
    }
  }

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

}  

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