/* effective potential change due to emvironment */
/* coded by Ikeuchi Mitsuru    */
/* ver 0.0.1   2004.03.13      */
/* ver 0.0.2   2004.03.18  improved display method */

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

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

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

  Choice cvw;
  Scrollbar scr,scs,scd;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;
  int atomSizeX = 4; 
  int atomSizeY = 4; 
  int atomSizeZ = 4; 
  int dispR = 40;         /* display molecular radious */
  int paintSizeX = atomSizeX*dispR; 
  int paintSizeY = atomSizeY*dispR;
  int wakuX = paintSizeX + dispR;
  int wakuY = paintSizeY + dispR;


  int Nmt = 8;           /* number of molecules    */
  int NN = Nmt + 1;
  double t = 0.0;              /* time (s)               */
  double dt = 2.0*1.0e-15;     /* time division (s)      */
  double sgm = 3.40e-10;       /* L-J sigma for Ar (m)   */
  int gerNrCount = 0;

  double contTemp = 300.0;
  double scale = 1.0;
  double scaleDens = 1.0;
  int dispMode = 4;
  int iptr = 1;

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

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

  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 AU = 1.661e-27;
  double FS[][] = {
  /*  atom      0 mass    1 d    2 A(eV)  3 b  4 c     5 c0       6 c1       7 c2    8 a */
  /* 0 W */ {183.85 *AU,4.400224,1.896373,0.0,3.25,47.1346499,-33.7665655, 6.2541999, 3.16, 1.01},
  /* 1 Mo*/ { 95.94 *AU,4.114825,1.887117,0.0,3.25,43.4475218,-31.9332978, 6.0804249, 3.15, 1.01},
  /* 2 Ta*/ {180.948*AU,4.076980,2.591061,0.0,4.20, 1.2157373,  0.0271471,-0.1217350, 3.30, 1.05},
  /* 3 Fe*/ { 55.847*AU,3.569745,1.828905,1.8,3.40, 1.2371147, -0.3592185,-0.0385607, 2.87, 0.97},
  /* 4 Cr*/ { 51.996*AU,3.915720,1.453418,1.8,2.90,29.1429813,-23.3975027, 4.7578297, 2.88, 0.98},
  /* 5 Nb*/ { 92.906*AU,3.915354,3.013789,0.0,4.20,-1.5640104,  2.0055779,-0.4663764, 3.30, 1.08},
  /* 6 V */ { 50.942*AU,3.692767,2.010637,0.0,3.80,-0.8816318,  1.4907756,-0.3976370, 3.03, 1.06}
  };
  int MAT1 = 0; /* W */

  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 reg[][] = new int[NN][250];   /* registration */
  double rijReg[][] = new double[NN][250];   /* registration */

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

  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 = 44;
  int NFy = 44;
  int NFz = 44;
  double dens[][][] = new double[NFx][NFy][NFz];
  double adens[][][] = new double[17][17][17];

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

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

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

    cvw = new Choice();
    cvw.addItem("e density"); cvw.addItem("real"); 
    cvw.addItem("ball"); cvw.addItem("line"); cvw.addItem("ball+line");
    cvw.addItem("dens+ball");
    cvw.select("ball+line");
    cvw.addItemListener(this);

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

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

    addMouseListener(this);  
    addMouseMotionListener(this);

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

    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() == cvw){
      dispMode = cvw.getSelectedIndex();
    }
  }

  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());
    } else if(ev.getSource() == scd) { 
      scaleDens = 0.01*(double)(scd.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;
    int bx,by,bz;
    double tm;
    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();

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

    gbx = 80;
    gby = 100;

    farPoint = findMin();
    for (i=0; i<12; i++) {
      gx = (int)(paintSizeX*scale*pwkx[boxp[i][0]]/xMax)+gbx;
      gy = (int)(paintSizeY*scale*pwky[boxp[i][0]]/yMax)+gby;
      g2x = (int)(paintSizeX*scale*pwkx[boxp[i][1]]/xMax)+gbx;
      g2y = (int)(paintSizeY*scale*pwky[boxp[i][1]]/yMax)+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=4; i<NFx-4; i++) {
        for (j=4; j<NFy-4; j++) {
          for (k=4; k<NFz-4; k++) {
            if (dens[i][j][k]>1.0) {

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

              gx = (int)(paintSizeX*scale*rx/xMax)+gbx;
              gy = (int)(paintSizeY*scale*ry/yMax)+gby;
              gz =  (int)(0.1*(5.0*(scale*rz/zMax)+5.0)*scaleDens*dens[i][j][k]/6.0);

              if (gz>5) {
                gz = 5;
              }

              gOff.setColor(Color.getHSBColor(
               (float)(0.66-1.5*dens[i][j][k]/100.0),0.9f,0.9f));
              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)(paintSizeX*scale*px[j]/xMax)+gbx;
        gy = (int)(paintSizeY*scale*py[j]/yMax)+gby;
        gz = (int)(14.0*scale*(0.6+0.3*(pz[j]/zMax))*FS[MAT1][8]);

        gOff.setColor(Color.getHSBColor((float)(1.0/8.0*j),0.95f,0.95f));
        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)(paintSizeX*scale*px[j]/xMax)+gbx;
        gy = (int)(paintSizeY*scale*py[j]/yMax)+gby;
        gz = (int)(20.0*(scale*pz[j]/zMax))+6;

        gOff.setColor(Color.getHSBColor((float)(1.0/8.0*j),0.95f,0.95f));
        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)(paintSizeX*scale*px[j]/xMax)+gbx;
        gy = (int)(paintSizeY*scale*py[j]/yMax)+gby;
        gz = (int)(20.0*(scale*pz[j]/zMax))+6;

        for(k=1;k<reg[i][0];k++) {
          jj = reg[j][k];
          if (jj>j) {
            if (distance(j,jj)<1.3*FS[MAT1][8]*1.0e-10) {
              g2x = (int)(paintSizeX*scale*px[jj]/xMax)+gbx;
              g2y = (int)(paintSizeY*scale*py[jj]/yMax)+gby;
              gOff.setColor(Color.gray);
              gOff.drawLine(gx,gy, g2x, g2y);
            }
          }
        }
      }

    } else if (dispMode==4) {

      for (i=0; i<Nmt; i++) {
        j = srtz[0][i];
        gx = (int)(paintSizeX*scale*px[j]/xMax)+gbx;
        gy = (int)(paintSizeY*scale*py[j]/yMax)+gby;
        gz = (int)(20.0*(scale*pz[j]/zMax))+6;

        for(k=1;k<reg[i][0];k++) {
          jj = reg[j][k];
          if (jj>j) {
            if (distance(j,jj)<1.3*FS[MAT1][8]*1.0e-10) {
              g2x = (int)(paintSizeX*scale*px[jj]/xMax)+gbx;
              g2y = (int)(paintSizeY*scale*py[jj]/yMax)+gby;
              gOff.setColor(Color.gray);
              gOff.drawLine(gx,gy, g2x, g2y);
            }
          }
        }

        gOff.setColor(Color.getHSBColor((float)(1.0/8.0*j),0.95f,0.95f));
        gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);

      }

    } else if (dispMode==5) {

      for (i=0; i<Nmt; i++) {
        bx = (int)(xx[i]/5.0e-11);
        by = (int)(yy[i]/5.0e-11);
        bz = (int)(zz[i]/5.0e-11);
        dens[bx+8][by+8][bz+8] += 200.0;
      }

      for (i=4; i<NFx-4; i++) {
        for (j=4; j<NFy-4; j++) {
          for (k=4; k<NFz-4; k++) {
            if (dens[i][j][k]>1.0) {

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

              gx = (int)(paintSizeX*scale*rx/xMax)+gbx;
              gy = (int)(paintSizeY*scale*ry/yMax)+gby;
              gz =  (int)(0.1*(5.0*(scale*rz/zMax)+5.0)*scaleDens*dens[i][j][k]/6.0);

              gOff.setColor(Color.getHSBColor(
                (float)(0.66-1.5*dens[i][j][k]/100.0),0.9f,0.9f));
              if (gz>5 && gz<20.0) {
                gz = 5;
              } else if (gz>=20.0) {
                gz = (int)(12.0*(scale*rz/zMax))+6;
                gOff.setColor(Color.green);
              }
              gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
            }
          }
        }
      }

    }


    for (i=0; i<12; i++) {
      gx = (int)(paintSizeX*scale*pwkx[boxp[i][0]]/xMax)+gbx;
      gy = (int)(paintSizeY*scale*pwky[boxp[i][0]]/yMax)+gby;
      g2x = (int)(paintSizeX*scale*pwkx[boxp[i][1]]/xMax)+gbx;
      g2y = (int)(paintSizeY*scale*pwky[boxp[i][1]]/yMax)+gby;
      if (boxp[i][0]!=farPoint && boxp[i][1]!=farPoint) { 
        gOff.setColor(Color.black); 
        gOff.drawLine(gx,gy, g2x, g2y);
      }
    }

    drawEffectivePotential(iptr);
	
    gOff.setColor(Color.black);
    gOff.drawString("view",630/6*0+10,40);
    gOff.drawString("cont T="+contTemp+"K",630/6*3+10,40);
    gOff.drawString("scale="+(int)(scale*100.0+0.5)+"%",630/6*4+10,40);
    gOff.drawString("dens="+(int)(scaleDens*100.0+0.5)+"%",630/6*5+10,40);

    gOff.setColor(Color.blue);
    gOff.drawString("t="+Math.floor(t*1.0e15+0.5)+"(fs)",630/6*2+10,40);
    gOff.drawString("T="+temperature()+" K",630/6*3+10,60);

  }

  double distance(int i,int j) {
    double x,y,z;

    x = xx[i]-xx[j];
    y = yy[i]-yy[j];
    z = zz[i]-zz[j];

    return(Math.sqrt(x*x+y*y+z*z));
  }

  void drawEffectivePotential(int ip) {
    int i,ipp,gbx,gby;
    int v1,v2;
    double ra,fp;
    
    gbx = 400;
    gby = 80;
    gOff.setColor(Color.gray);
    for (i=0; i<200; i+=25) {
      gOff.drawLine(gbx+i,gby, gbx+i, gby+200);
    }
    for (i=0; i<200; i+=25) {
      gOff.drawLine(gbx,gby+i, gbx+200, gby+i);
    }
    gOff.drawRect(gbx,gby,200,200);
    gOff.drawString("1A",gbx-5,gby+220);
    gOff.drawString("3A",gbx+90,gby+220);
    gOff.drawString("5A",gbx+190,gby+220);
    gOff.drawString("0",gbx-20,gby+200-25-20);
    gOff.drawString("5eV",gbx-30,gby+200-20-150);
    
    for (ipp=0; ipp<8; ipp++) {
      gOff.setColor(Color.getHSBColor((float)(1.0/8.0*ipp),0.95f,0.95f));
      fp = dFdroh(ipp);
      for(i=0;i<200;i++) {
        ra = 0.02*i+1.0;
        v1 = (int)(25.0*effectivePotential(ra,fp));
        v2 = (int)(25.0*effectivePotential(ra+0.02,fp));
        if (v1<180) {
          gOff.drawLine(gbx+i,gby+200-50-v1, gbx+i, gby+200-50-v2);
        }
      }
    }
  }
  
  double effectivePotential(double ra, double fp) {
    return( potV(ra)-2.0*FSA*fp*fai(ra) );
  }


/*----------------------------------  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<17; ix++) {
      for (iy=0; iy<17; iy++) {
        for (iz=0; iz<17; iz++) {
          rA = 0.5*Math.sqrt((double)((ix-8)*(ix-8)+(iy-8)*(iy-8)+(iz-8)*(iz-8)));
          if (rA<FSd) {
            adens[ix][iy][iz] = (rA-FSd)*(rA-FSd);
          } else {
            adens[ix][iy][iz] = 0.0;
          }
        }
      }
    }
  }


  void setdens() {
    int i,ix,iy,iz,bx,by,bz;
    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]/5.0e-11);
      by = (int)(yy[i]/5.0e-11);
      bz = (int)(zz[i]/5.0e-11);

      for (ix=0; ix<17; ix++) {
        for (iy=0; iy<17; iy++) {
          for (iz=0; iz<17; iz++) {
            dens[bx+ix][by+iy][bz+iz] += adens[ix][iy][iz];
          }
        }
      }
    }
  }


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

  double temperature() {
    int i;
    double s;

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += 0.5*FSmass*(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 ii) {
    double s;

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

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

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

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

    FSmass = FS[MAT1][0];
    FSd = FS[MAT1][1];
    FSA = FS[MAT1][2];
    FSb = FS[MAT1][3];
    FSc = FS[MAT1][4];
    FSc0 = FS[MAT1][5];
    FSc1 = FS[MAT1][6];
    FSc2 = FS[MAT1][7];

    a = FS[MAT1][8]*1.0e-10;

    d = 1.0*sgm;

    for (i=0; i<8; i++) {
      ix = i/4;
      iy = (i%4)/2;
      iz = (i%4)%2;
      xx[i] = 2.0*sgm*(double)(ix)+d;
      yy[i] = 2.0*sgm*(double)(iy)+d;
      zz[i] = 2.0*sgm*(double)(iz)+d;
      vx[i] = 400.0*Math.random()-200.0;
      vy[i] = 400.0*Math.random()-200.0;
      vz[i] = 400.0*Math.random()-200.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;

    vAjustment();

    regNear();
    for (i=0; i<10; i++) {
      timeEvolution();
    }
    if (dispMode==0 || dispMode==5) {
      setdens();
    }
  }

  void timeEvolution(){
    int i; 
    double rr;

    t = t + dt;

    rr = 1.0;

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

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

  void forceCalc() {
    int i,j,k;
    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) {
          rij = rijReg[i][k];
          f = force(rij,fpi);
          fxij = f*(xx[i]-xx[j])/rij;
          fyij = f*(yy[i]-yy[j])/rij;
          fzij = f*(zz[i]-zz[j])/rij;
        } else {
          fxij = 0.0;
          fyij = 0.0;
          fzij = 0.0;
        }
        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 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++) {
        r2=(xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])+(zz[i]-zz[j])*(zz[i]-zz[j]);
        if (r2<rrg2) {
          regNr[i][k]=j;
          k += 1;
        }
      }
      for(j=i+1;j<Nmt;j++) {
        r2=(xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])+(zz[i]-zz[j])*(zz[i]-zz[j]);
        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 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];
        r2=(xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])+(zz[i]-zz[j])*(zz[i]-zz[j]);
        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;
    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)*(1.0+FSb*(r-FSd)/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)+3.0*FSb*(r-FSd)*(r-FSd)/FSd );
    } else {
      return ( 0.0 );
    }
  }

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

}  

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