/* W(110) - Morse atom sputtering */
/* coded by Ikeuchi Mitsuru */
/* ver 0.0.1   2002.11.23   */
/* ver 0.0.2   2002.11.25  bug fixed in forceCalc() */
/* ver 0.0.3   2003.08.08 rewrited based on (316) template */
/* ver 0.0.4   2003.09.05  registoration bug fixed */
/* ver 0.0.5   2003.11.21  improved display code */
/* ver 0.0.6   2004.02.08  improved pannel layout */

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

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

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

  Button bt_reset,bt_start;
  Choice cm1,cvw;
  Scrollbar scr,sce;
  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 = 25.0/sgm;

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

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

  int tempMode = 0;
  double contTemp = 100.0;
  double scale = 1.0;
  double ESet = 100.0;
  int dispMode = 4;
  int started = 0;

  double rc = 1.0e-9;
  double rc2 = rc*rc;

  double rreg = 1.5e-9;
  double rreg2 = rreg*rreg;

  double FSmass = 183.85*1.661e-27;
  double FSd = 4.400224;
  double FSA = 1.896373;
  double FSc = 3.25;
  double FSc0 = 47.1346499;
  double FSc1 = -33.7665655;
  double FSc2 = 6.2541999;


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


  double AU = 1.661e-27;
  double EE = 1.602e-19;
  double AA = 1.0e-10;
  double Morse[][] = {
  /*         mass         charge   D(eV)     A(A^-1)    r0(A)   */
  /* 0 W */ {183.85 *AU,  0.0*EE, 0.9906*EE, 1.4116/AA, 3.032*AA},
  /* 1 Mo*/ { 95.94 *AU,  0.0*EE, 0.8032*EE, 1.5079/AA, 2.976*AA},
  /* 2 Cr*/ { 51.996*AU,  0.0*EE, 0.4414*EE, 1.5721/AA, 2.754*AA},
  /* 3 Fe*/ { 55.847*AU,  0.0*EE, 0.4174*EE, 1.3885/AA, 2.845*AA},
  /* 4 Ni*/ { 58.71 *AU,  0.0*EE, 0.4205*EE, 1.4199/AA, 2.780*AA},
  /* 5 Al*/ { 26.98 *AU,  0.0*EE, 0.2703*EE, 1.1646/AA, 3.253*AA},
  /* 6 Pb*/ {207.19 *AU,  0.0*EE, 0.2348*EE, 1.1836/AA, 3.733*AA},
  /* 7 Cu*/ { 63.54 *AU,  0.0*EE, 0.3429*EE, 1.3588/AA, 2.866*AA},
  /* 8 Ag*/ {107.87 *AU,  0.0*EE, 0.3323*EE, 1.3690/AA, 3.115*AA},
  /* 9 Ca*/ { 40.08 *AU,  0.0*EE, 0.1623*EE, 0.8054/AA, 4.569*AA},
  /*10 Sr*/ { 87.62 *AU,  0.0*EE, 0.1513*EE, 0.7878/AA, 4.988*AA},
  /*11 Ba*/ {137.34 *AU,  0.0*EE, 0.1416*EE, 0.6570/AA, 5.373*AA},
  /*12 Na*/ { 22.99 *AU,  0.0*EE, 0.0633*EE, 0.5900/AA, 5.336*AA},
  /*13 K */ { 39.102*AU,  0.0*EE, 0.0542*EE, 0.4977/AA, 6.369*AA},
  /*14 Rb*/ { 85.47 *AU,  0.0*EE, 0.0464*EE, 0.4298/AA, 7.207*AA},
  /*15 Cs*/ {132.905*AU,  0.0*EE, 0.0449*EE, 0.4157/AA, 7.557*AA},
  /*16 Ne*/ { 20.183*AU,  0.0*EE, 0.0031*EE, 1.6500/AA, 3.076*AA},
  /*17 Ar*/ { 39.948*AU,  0.0*EE, 0.0104*EE, 1.3400/AA, 3.816*AA},
  /*18 Kr*/ { 83.80 *AU,  0.0*EE, 0.0141*EE, 1.2500/AA, 4.097*AA},
  /*19 Xe*/ {131.30 *AU,  0.0*EE, 0.0200*EE, 1.2400/AA, 4.467*AA},
  /*20 Hg*/ {200.59 *AU,  0.0*EE, 0.0734*EE, 1.4900/AA, 3.255*AA}
  };
  int NmatKind = 21;
  double MorseD[][] = new double[NmatKind][NmatKind];
  double MorseA[][] = new double[NmatKind][NmatKind];
  double Morser0[][] = new double[NmatKind][NmatKind];

  int MAT0 = 0; /* W fixed */
  int MAT1 = 11;/* Ba */
  int MAT2 = 17;/* Ar */


  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 regMorse[][] = new int[NN][860];

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

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

  int kindT[] = new int[NmatKind];
  double forceTable[][][] = new double[3][3][10200];

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

/*----------------------------- 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("W");  cm1.addItem("Mo"); cm1.addItem("Cr");
    cm1.addItem("Fe"); cm1.addItem("Ni"); cm1.addItem("Al");
    cm1.addItem("Pb"); cm1.addItem("Cu"); cm1.addItem("Ag");
    cm1.addItem("Ca"); cm1.addItem("Sr"); cm1.addItem("Ba");
    cm1.addItem("Na"); cm1.addItem("K");  cm1.addItem("Rb");
    cm1.addItem("Cs");
    cm1.addItem("Ne"); cm1.addItem("Ar");
    cm1.addItem("Kr"); cm1.addItem("Xe"); cm1.addItem("Hg");
    cm1.addItemListener(this);
    cm1.select("Ar");

    cvw = new Choice();
    cvw.addItem("ball"); cvw.addItem("line");
    cvw.addItem("both"); cvw.addItem("real");
    cvw.addItem("temp");
    cvw.addItemListener(this);
    cvw.select("temp");

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

    bt_start = new Button("start");
    bt_start.addActionListener(this);

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

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

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,7,5,0));
    pnl.add(cm1);
    pnl.add(cvw);
    pnl.add(bt_reset);
    pnl.add(bt_start);
    pnl.add(new Label("  "));
    pnl.add(scr);
    pnl.add(sce);
    add(pnl,"North");

    setMorse();
    setInitialPosition();
    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 actionPerformed(ActionEvent ev){ 
    if (ev.getSource() == bt_start) {
      if (started==0) {
        started=1;
        vx[Nmt-1] -= Math.sqrt(2.0*ESet*1.6e-19/FSmass);
        repaint();
      }
    } else if (ev.getSource() == bt_reset) {
      started=0;
      t=0.0;
      setInitialPosition();
      repaint();
    }
  }

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == cm1){
      MAT2 = cm1.getSelectedIndex();
      kind[Nmt-1] = MAT2;
      if (started==1) {
        started=0;
        t = 0.0;
        setInitialPosition();
      }
      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() == sce) { 
      ESet = 1.0*(double)(sce.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,gbz,gb;
    double sc,sz,tmp,tm,clr;

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

    rotate(theta,fai);
    zSort();
    tmp = temperature();
    gbx = 70;
    gby = 90;
    sc = dispScale*scale;
    sz = 4.5*scale/AA;

    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*Morse[kind[j]][4]*(0.4*(pz[j]/zMax)+0.6));
        clr = ballColor(kind[j]);
        gOff.setColor(Color.getHSBColor(
          (float)clr,0.9f,(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)<0.55*(Morse[kind[j]][4]+Morse[kind[jj]][4]) ) {
              g2x = (int)(sc*px[jj])+gbx;
              g2y = (int)(sc*py[jj])+gby;
              gOff.setColor(Color.gray);
              gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
        if (kind[j]!=MAT0) {
          gz = (int)(sz*Morse[kind[j]][4]*(0.4*(pz[j]/zMax)+0.6));
          clr = ballColor(kind[j]);
          gOff.setColor(Color.getHSBColor(
            (float)clr,0.9f,(float)(0.4*pz[j]/zMax+0.4)));
          gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
        }
      }

    } 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;
        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
          if ( distance(j,jj)<0.55*(Morse[kind[j]][4]+Morse[kind[jj]][4]) ) {
              g2x = (int)(sc*px[jj])+gbx;
              g2y = (int)(sc*py[jj])+gby;
              gOff.setColor(Color.gray);
              gOff.drawLine(gx,gy, g2x, g2y);
          }
        }
        gz = (int)(sz*Morse[kind[j]][4]*(0.4*(pz[j]/zMax)+0.6));
        clr = ballColor(kind[j]);
        gOff.setColor(Color.getHSBColor(
          (float)clr,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)(1.6*sz*Morse[kind[j]][4]*(0.4*(pz[j]/zMax)+0.6));
        clr = ballColor(kind[j]);
        gOff.setColor(Color.getHSBColor(
          (float)clr,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*Morse[kind[j]][4]*(0.4*(pz[j]/zMax)+0.6));
        clr = ballColor(kind[j]);
        tm = 0.5*Morse[kind[j]][0]*(vx[j]*vx[j]+vy[j]*vy[j]+vz[j]*vz[j])/(1.5*1.38e-23)/2000.0 + 0.3;
        if(tm>0.99) {
          tm = 0.99;
        }
        gOff.setColor(Color.getHSBColor( (float)clr, 0.99f, (float)tm ));
        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);
      }
    }

    gOff.setColor(Color.black);
    gOff.drawString("atom",630/7*0+10,40);
    gOff.drawString("view",630/7*1+10,40);
    gOff.drawString("T="+(int)(contTemp+0.5)+"K",630/7*5+10,40);
    gOff.drawString("E="+(int)(ESet+0.5)+"eV",630/7*6+10,40);
    gb = 520;
    gOff.setColor(Color.blue);
    gOff.drawString("t="+Math.floor(t*1.0e15+0.5)+"(fs)",gb,110);
    gOff.drawString("T="+temperature()+" K",gb,130);
    gOff.setColor(Color.green);
    gOff.drawString("crystal=W(110)",gb,150);

  }

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

  double ballColor(int knd) {
    double c;

    c = 0.0;
    if (knd==MAT0) {
      c = 0.33;
    } else if (knd==MAT1) {
      c = 0.15;
    } else if (knd==MAT2) {
      c = 0.85;
    }
    return ( c );
  }


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


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

  double temperature() {
    int i;
    double s;

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += 0.5*Morse[kind[i]][0]*(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 temp(int k) {
    int i,n;
    double s;

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

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

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

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

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

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

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

  void setMorse(){
    int i,j; 

    for (i=0; i<NmatKind; i++) {
      for (j=0; j<NmatKind; j++) {
        MorseD[i][j] = Math.sqrt(Morse[i][2]*Morse[j][2]);
        MorseA[i][j] = 0.5*(Morse[i][3]+Morse[j][3]);
        Morser0[i][j] = 0.5*(Morse[i][4]+Morse[j][4]);
      }
    }
  }

  void setForceTable(){
    int i,j,ir,k0,k1; 
    double d,a,r0,y,r,sc;
    int kkk[] = {MAT0,MAT1,MAT2};

    for (i=0; i<3; i++) {
      k0=kkk[i];
      for (j=0; j<3; j++) {
        k1=kkk[j];

        d = Math.sqrt(Morse[k0][2]*Morse[k1][2]);
        a = 0.5*(Morse[k0][3]+Morse[k1][3]);
        r0 = 0.5*(Morse[k0][4]+Morse[k1][4]);

        for (ir=100; ir<10000; ir++) {
          r = (double)ir*1.0e-13;
          y = Math.exp(-a*(r-r0));
          forceTable[i][j][ir]=2.0*d*a*y*(y-1.0);
        }
        for (ir=0; ir<100; ir++) {
          forceTable[i][j][ir]=forceTable[i][j][100];
        }

        for (ir=6000; ir<10000; ir++) {
          sc = 1.0-((double)(ir-6000))/4000.0;
          forceTable[i][j][ir] = sc*forceTable[i][j][ir];
        }
      }
    }
  }

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

    kindT[MAT0] = 0; kindT[MAT1] = 1; kindT[MAT2] = 2;
    setForceTable();

    a = 3.16e-10;
    b = 0.5*Math.sqrt(2.0)*a;
    c = 0.5*a;

    i = 0;
    for (ix=0; ix<8; ix++) {
      for (iy=0; iy<10; iy++) {
        for (iz=0; iz<12; iz++) {
          j = ix+iy+iz;
          if (j%2==0) {
            kind[i] = MAT0;
            xx[i] = b*(double)(ix)+0.3*sgm;
            yy[i] = b*(double)(iy)+0.5*b;
            zz[i] = c*(double)(iz)+0.5*c;
            vx[i] = 100.0*Math.random()-50.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;
          }
        }
      }
    }

    kind[i] = MAT2;
    xx[i] = 0.95*xMax;
    yy[i] = 0.5*yMax+1.0*sgm*(Math.random()-0.5);
    zz[i] = 0.5*zMax+1.0*sgm*(Math.random()-0.5);
    vx[i] = 0.0;
    vy[i] = 0.0;
    vz[i] = 0.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;

    zSortInit();

  }

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

  void calcposition(){
    int i;
/*
    if (t<100*dt) {
      vAjustment();
    }
*/
    registrationMorse();
    regNear();
    for (i=0; i<10; i++) {
      timeEvolution();
    }
  }

  void timeEvolution(){
    int i; 

    t = t + dt;

    for (i=0; i<Nmt; i++) {
      xx[i] += vx[i]*dt+0.5*dt*dt*ffx[i]/Morse[kind[i]][0];
      yy[i] += vy[i]*dt+0.5*dt*dt*ffy[i]/Morse[kind[i]][0];
      zz[i] += vz[i]*dt+0.5*dt*dt*ffz[i]/Morse[kind[i]][0];
      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])/Morse[kind[i]][0];
      vy[i] += 0.5*dt*(ffy[i]+fby[i])/Morse[kind[i]][0];
      vz[i] += 0.5*dt*(ffz[i]+fbz[i])/Morse[kind[i]][0];
    }

    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 r2,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++) {
      if(kind[i]==MAT0) {
        fpi = dFdroh(i);
        for(k=1;k<reg[i][0];k++) {
          j = reg[i][k];
          if (j>i && kind[j]==MAT0) {
            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;
          }
        }
      }
    }

    for(i=0;i<Nmt;i++) {
      for(k=1;k<regMorse[i][0];k++) {
        j = regMorse[i][k];
        if (kind[i]!=MAT0 || kind[j]!=MAT0) {
          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<rc2) {
            rij = Math.sqrt(r2);
            f = forceMorse(rij,kind[i],kind[j]);
            fxij = f*(xi-xj)/rij;
            fyij = f*(yi-yj)/rij;
            fzij = f*(zi-zj)/rij;

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

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

  }

  double forceMorse(double r, int knd0, int knd1) {
    return ( forceTable[kindT[knd0]][kindT[knd1]][(int)(r*1.0e13)] );
  }

/*-------------------------------  registration   -------------*/

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

    for(i=0;i<Nmt;i++) {
      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<rreg2) {
          regMorse[i][k]=j;
          k += 1;
        }
      }
      regMorse[i][0]=k;
    }
  }

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

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

    rrg = FSc*1.0e-10;
    if (rrg<FSd*1.0e-10) {
      rrg = FSd*1.0e-10;
    }
    rrg += 3000.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;
    }

  }

/*--------------------  Finnis-Sinclair Potential Calculation   ---*/

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

    ra = r*1.0e10;
    rt = -dpotVdr(ra)+2.0*FSA*fpi*dfaidr(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 = FSc*1.0e-10;
    if (rrg<FSd*1.0e-10) {
      rrg = FSd*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) {
    double rh,rt;

    rh = roh(i);
    rt = 0.0;
    if (rh>0.0) { 
      rt = 0.5*Math.sqrt(1.0/rh);
    }
    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 += fai( rijReg[i][k]*1.0e10 );
    }
    return(s);
  }


  double potV(double r) {
    if (r<FSc) {
      return ( (r-FSc)*(r-FSc)*(FSc0+FSc1*r+FSc2*r*r) );
    } else {
      return ( 0.0 );
    }
  }

  double fai(double r) {
    if (r<FSd) {
      return ( (r-FSd)*(r-FSd) );
    } else {
      return ( 0.0 );
    }
  }

  double dpotVdr(double r) {
    if (r<FSc) {
      return ( 2.0*(r-FSc)*(FSc0+FSc1*r+FSc2*r*r)+(r-FSc)*(r-FSc)*(FSc1+2.0*FSc2*r) );
    } else {
      return ( 0.0 );
    }
  }

  double dfaidr(double r) {
    if (r<FSd) {
      return ( 2.0*(r-FSd) );
    } else {
      return ( 0.0 );
    }
  }

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

}  

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