/* applet No. 951
 *
 * bcc metal -  AFS potential
 *
 * Created by Ikeuchi Mitsuru on May 27 2006.
 * Copyright (c) 2006 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.05.27  created
 * ver 0.0.2   2006.06.17  bug fixed in velocity plot
 *
 */

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

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

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

  Choice ch_tempMode,ch_mat,ch_view;
  Scrollbar sc_temp,sc_scale;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  double xMax = 6.0*1.0e-9;
  double yMax = 6.0*1.0e-9;
  double zMax = 6.0*1.0e-9; 

  int resetSW = 0;
  int started = 1;
  int tempMode = 1;
  double contTemp = 300.0;
  int dispMode = 7;

  double dispWidth = 250.0;
  double dispScale = (dispWidth/xMax);
  double viewScale = 1.0;

  int Ncube = 16;
  int Nmt = Ncube*Ncube*Ncube*2; // number of atoms
  int NN = Nmt + 1;
  double t = 0.0;
  double dt = 1.0*1.0e-15;

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

  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 FSbb = 90.3;
  double FSalpha = 1.2;
  double FSb0 = 2.7411;


  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 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.00},
  /* 1 Mo*/ { 95.94 *AU,4.114825,1.887117,0.0,3.25,43.4475218,-31.9332978, 6.0804249, 3.15, 1.00},
  /* 2 Fe*/ { 55.847*AU,3.569745,1.828905,1.8,3.40, 1.2371147, -0.3592185,-0.0385607, 2.87, 1.00},
  /* 3 Cr*/ { 51.996*AU,3.915720,1.453418,1.8,2.90,29.1429813,-23.3975027, 4.7578297, 2.88, 1.00},
  /* 4 Ta*/ {180.948*AU,4.076980,2.591061,0.0,4.20, 1.2157373,  0.0271471,-0.1217350, 3.30, 1.00},
  /* 5 Nb*/ { 92.906*AU,3.915354,3.013789,0.0,4.20,-1.5640104,  2.0055779,-0.4663764, 3.30, 0.99},
  /* 6 V */ { 50.942*AU,3.692767,2.010637,0.0,3.80,-0.8816318,  1.4907756,-0.3976370, 3.03, 1.00}
  };
  int MAT1 = 0; /* W */

  double Ackland[][] = {
  /*  atom   0 B(eV/A3)  1 alpha    2 b0(A) */
  /* 0 W */ {   90.3,     1.2,      2.7411 },
  /* 1 Mo*/ { 1223.0,     3.9,      2.7255 },
  /* 2 Fe*/ {    0.0,     1.0,      0.0    },
  /* 3 Cr*/ {    0.0,     1.0,      0.0    },
  /* 4 Ta*/ {   91.0,     1.05,     2.8629 },
  /* 5 Nb*/ {   48.0,     0.8,      2.8585 },
  /* 6 V */ {   23.0,     0.5,      2.6320 }
  };


  int kind[] = new int[NN];
  double xx[] = new double[NN]; /* i-th position */
  double yy[] = new double[NN];
  double zz[] = new double[NN];
  double vx[] = new double[NN]; /* i-th velocity */
  double vy[] = new double[NN];
  double vz[] = new double[NN];
  double ffx[] = new double[NN]; /* i-th force    */
  double ffy[] = new double[NN];
  double ffz[] = new double[NN];

  double virx[] = new double[NN];
  double viry[] = new double[NN];
  double virz[] = new double[NN];

  int reg[][] = new int[NN][50];
  double rijReg[][] = new double[NN][50];

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

  int srtz[][] = new int[2][NN];

  int rdf[] = new int[260];
  int vdf[] = new int[220];
  double meanvdf[] = new double[220];
  double ffvMaxwell[] = new double[220];

  double virial[] = new double[256];
  int virialp = 0;

  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,420);
    setBackground(Color.white);

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

    ch_mat = new Choice();
    ch_mat.addItem("W");
    ch_mat.addItem("Mo");
    ch_mat.addItem("Fe");
    ch_mat.addItem("Cr");
    ch_mat.addItem("Ta");
    ch_mat.addItem("Nb");
    ch_mat.addItem("V");
    ch_mat.addItemListener(this); 

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

    ch_view = new Choice();
    ch_view.addItem("ball");
    ch_view.addItem("line");
    ch_view.addItem("both");
    ch_view.addItem("temp");
    ch_view.addItem("real");
    ch_view.addItem("velocity");
    ch_view.addItem("v-space");
    ch_view.addItem("ball+Veff");
    ch_view.addItem("rho@atom");
    ch_view.addItem("dF/drho");
    ch_view.addItemListener(this);
    ch_view.select("ball+Veff");

    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(1,6,5,0));
	/* pnl.add(new Label("  ")); */
    pnl.add(ch_mat);
    pnl.add(new Label("  "));
    pnl.add(ch_tempMode);
    pnl.add(sc_temp);
    pnl.add(sc_scale);
    pnl.add(ch_view);
    add(pnl,"North");

    setInitialCondition();
  }

  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() == ch_mat){
      MAT1 = ch_mat.getSelectedIndex();
      resetSW = 1;
    } else if (ev.getSource() == ch_tempMode){
      tempMode = ch_tempMode.getSelectedIndex();
    } else if (ev.getSource() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_temp) { 
      contTemp = 10.0*(double)(sc_temp.getValue());
      vAjustment();
    } else if(ev.getSource() == sc_scale) { 
      viewScale = 0.01*(double)(sc_scale.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 {
        timeEvolution();
        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 gbx,gby;
    double tmp;

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

    tmp = temperature();

    gbx = 60;
    gby = 100;
    if (dispMode==0) {
      ballPlot(gbx,gby,0.6,0);
      rdfPlot(400,100,1.0);

    } else if (dispMode==1) {
      ballPlot(gbx,gby,0.0,2);
      rdfPlot(400,100,1.0);

    } else if (dispMode==2) {
      ballPlot(gbx,gby,0.5,2);
      rdfPlot(400,100,1.0);

    } else if (dispMode==3) {
      ballPlot(gbx,gby,0.6,1);
      rdfPlot(400,100,1.0);

    } else if (dispMode==4) {
      ballPlot(gbx,gby,1.0,0);
      rdfPlot(400,100,1.0);

    } else if (dispMode==5) {
      velocityPlot(gbx,gby, 50.0);
      vdfPlot(400,100,1.0);

    } else if (dispMode==6) {
      vSpacePlot(gbx,gby);
      vdfPlot(400,100,1.0);

    } else if (dispMode==7) {
      ballPlot(gbx,gby,0.6,0);
      drawEffectivePotential(400, 100);

    } else if (dispMode==8) {
      rhoPlot(gbx,gby);
      drawEffectivePotential(400, 100);

    } else if (dispMode==9) {
      fpPlot(gbx,gby);
      drawEffectivePotential(400, 100);

    } else if (dispMode==10) {
      ;

    }


    gOff.setColor(Color.black);
    gOff.drawString("atom",630/6*0+10,40);
    gOff.drawString("temp control",630/6*2+10,40);
    gOff.drawString("T_cont="+(int)(contTemp+0.5)+" K",630/6*3+10,40);
    gOff.drawString("scale="+(int)(viewScale*100.0+0.499)+"%",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gOff.setColor(Color.blue);
    gOff.drawString("a="+FS[MAT1][8]+"A",630/6*0+10,60);
    gOff.drawString("t="+Math.floor(t*1.0e15+0.5)+" fs",630/6*1+10,60);
    gOff.drawString("T="+(int)(tmp*10.0+0.5)/10.0+" K",630/6*3+10,60);
    gOff.drawString("<P>="+(int)(pressure()/1.0e6)+"MPa",630/6*4+10,60);

    gOff.setColor(Color.black);
    gOff.drawString("box="+(int)(xMax*1.0e10+0.5)/10.0+"x"+(int)(yMax*1.0e10+0.5)/10.0+"x"+(int)(zMax*1.0e10+0.5)/10.0+"nm",430,80);
  }

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

  void ballPlot(int gbx, int gby, double size,int mode) {
    int i,j,k,jj,gx,gy,gz,g2x,g2y;
    double sc,sz,bl;

    setWaku();
    rotate(theta,fai);
    zSort();

    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;
    sz = size*(1.732/2.0)*sc*FS[MAT1][8]*1.0e-10;
    bl = 0.935*FS[MAT1][8]*1.0e-10;

    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.2*pz[j]/zMax+0.8));

      if (mode==0) {
        gOff.setColor(Color.getHSBColor(0.15f,0.95f,(float)(0.5*pz[j]/zMax+0.2)));
      } else if (mode==1) {
	    gOff.setColor(Color.getHSBColor(
          (float)(0.66-temp(j)/3000.0),0.95f,(float)(0.5*pz[j]/zMax+0.2)));
	  } else if (mode==2) {
	    gOff.setColor(Color.getHSBColor(0.15f,0.3f,(float)(0.5*pz[j]/zMax+0.2)));
        for(k=1;k<reg[i][0];k++) {
          jj = reg[j][k];
          if (jj>j) {
            if (distance(j,jj)<bl) {
              g2x = (int)(sc*px[jj])+gbx;
              g2y = (int)(sc*py[jj])+gby;
              gOff.drawLine(gx,gy, g2x, g2y);
            }
          }
        }
		gOff.setColor(Color.getHSBColor(0.15f,0.95f,(float)(0.5*pz[j]/zMax+0.2)));
      }
      if (gz>0) gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }

    drawWaku(gbx,gby,1);
  }

  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 velocityPlot(int gbx, int gby, double mag) {
    int i,j,k,jj,gx,gy,gz,g2x,g2y;
    double sc,cosTh,sinTh,cosFi,sinFi;
    double p2x,p2y,p2z,rotp2x,rotp2y,rotp2z;

    setWaku();
    rotate(theta,fai);
    zSort();

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

    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;

      p2x = xx[j]+mag*vx[j]*20*dt;
      p2y = yy[j]+mag*vy[j]*20*dt;
      p2z = zz[j]+mag*vz[j]*20*dt;
      rotp2y = cosFi*(p2y-cy)+sinFi*(p2z-cz) + cy; 
      rotp2z = -sinFi*(p2y-cy)+cosFi*(p2z-cz) + cz;
      rotp2x = cosTh*(p2x-cx)+sinTh*(rotp2z-cz) + cx; 
      rotp2z = -sinTh*(p2x-cx)+cosTh*(rotp2z-cz) + cz;
      g2x = (int)(sc*rotp2x)+gbx;
      g2y = (int)(sc*rotp2y)+gby;
      gOff.setColor(Color.getHSBColor(0.01f,0.95f,(float)(0.3*pz[j]/zMax+0.5))); 
	  if (vx[j]>0.0) gOff.setColor(Color.getHSBColor(0.60f,0.95f,(float)(0.3*pz[j]/zMax+0.5)));
      gOff.drawLine(gx,gy, g2x,g2y);
    }

    drawWaku(gbx,gby,1);
  }

  void drawWaku(int gbx, int gby, int imode) {
    int i,gx,gy,g2x,g2y,farPoint;
    double sc,tmp;

    sc = dispScale*viewScale;
    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 (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.gray);
          if (tempMode==2) {
            tmp = temperature();
            if (tmp>contTemp) {
              gOff.setColor(Color.black); 
            } else {
              gOff.setColor(Color.red);
            }
          }
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
  }

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

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

    i=0;
    sgm = 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*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;
        }
      }
    }
  }

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

/*----------------------------- dF/drho plot --------------*/

  void fpPlot(int gbx, int gby) {
    int i,j,gx,gy,gz;
    double sc,sz,fp0;

    setWaku();
    rotate(theta,fai);
    zSort();

    drawWaku(gbx,gby,0);

    fp0 = dFdroh(Ncube*Ncube*Ncube/2+Ncube*Ncube/2+Ncube/2);
    sc = dispScale*viewScale;
    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      sz = 0.5*dFdroh(j)/fp0; if (sz>0.999) sz = 0.999;
      gz = (int)(10.0*sz+0.5);
      gOff.setColor(Color.getHSBColor((float)(0.66-0.66*sz),0.7f,(float)(0.5*pz[j]/zMax+0.2)));
      if (gz>0) gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }

    drawWaku(gbx,gby,1);
  }

/*----------------------------- rho plot --------------*/

  void rhoPlot(int gbx, int gby) {
    int i,j,gx,gy,gz;
    double sc,sz,rh0;

    setWaku();
    rotate(theta,fai);
    zSort();

    drawWaku(gbx,gby,0);

    rh0 = roh(Ncube*Ncube*Ncube/2+Ncube*Ncube/2+Ncube/2);
    sc = dispScale*viewScale;
    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      sz = 0.5*roh(j)/rh0; if (sz>0.999) sz = 0.999;
      gz = (int)(10.0*sz+0.5);
      gOff.setColor(Color.getHSBColor((float)(0.66-0.66*sz),0.7f,(float)(0.5*pz[j]/zMax+0.2)));
      if (gz>0) gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
    }

    drawWaku(gbx,gby,1);
  }

/*----------------------------- v-space plot --------------*/

  void vSpacePlot(int gbx, int gby) {
    int i,j,gx,gy,gz;
    double sc;


    setWaku();
    vRotate(theta,fai);
    zSort();

    drawWaku(gbx,gby,0);

    sc = dispScale*viewScale;

    for (i=0; i<Nmt; i++) {
      j = srtz[0][i];
      gx = (int)(sc*px[j])+gbx;
      gy = (int)(sc*py[j])+gby;
      gz = (int)(4.0*(0.2*pz[j]/zMax+0.8));
	  if (vx[j]>0) {
	    gOff.setColor(Color.getHSBColor(0.66f,0.95f,(float)(0.5*pz[j]/zMax+0.2)));
	  } else {
	    gOff.setColor(Color.getHSBColor(0.01f,0.95f,(float)(0.5*pz[j]/zMax+0.2)));
	  }
      gOff.fillOval(gx-gz/2,gy-gz/2, gz, gz);
	}

    drawWaku(gbx,gby,1);
  }

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

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

    sc = xMax/2000.0;
    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*(sc*vy[i]+yMax/2-cy)+sinFi*(sc*vz[i]+xMax/2-cz) + cy; 
      pz[i] = -sinFi*(sc*vy[i]+yMax/2-cy)+cosFi*(sc*vz[i]+zMax/2-cz) + cz;
    } 
    for (i=0; i<Nmt; i++) {
      px[i] = cosTh*(sc*vx[i]+xMax/2-cx)+sinTh*(pz[i]-cz) + cx; 
      pz[i] = -sinTh*(sc*vx[i]+xMax/2-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;
    }  
  }

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

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

    for (i=0; i<Nmt; i++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(32767.0*pz[i]/zMax);
    }
    qSort(0,Nmt-1);    
  }

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

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


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

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

    rdfSet();

    ascale = mag*(5000000.0/Nmt);
    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("0",gbx-5,gby+250+20);
    gOff.drawString("2A",gbx+90,gby+250+20);
    gOff.drawString("4A",gbx+190,gby+250+20);
    gOff.drawString("radial distribution",gbx+30,gby+300);

    gOff.setColor(Color.cyan);
    for (i=0; i<200; i++) {
	  rd = (int)(rdf[i]*(ascale/(i*i))+0.5);
      if (rdf[i]>0) {
        gOff.drawLine(gbx+i,gby+250-rd, gbx+i, gby+250);
      }
    }
  }

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

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

    for (i=0; i<Nmt; i++) {
      for(k=1;k<regNr[i][0];k++) {
        j = regNr[i][k];
        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]);
        rij = Math.sqrt(r2);
        if (rij<5.0e-10) {
          ir = (int)(rij/2.0e-12+0.5);
          rdf[ir] +=1;
        }
      }
    }
  }

/*------------------------ verocity distribution function  --------------*/

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

    setVelocityDistribution();
    setMaxwell(contTemp);

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

    gOff.setColor(Color.getHSBColor(0.6f,0.8f,0.6f));
    for (i=0; i<200; i++) {
      iy = (int)(ascale*vdf[i]);
      if (iy>0) {
        gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200);
      }
    }

    gOff.setColor(Color.getHSBColor(0.15f,0.8f,0.6f));
    for (i=0; i<200-1; i++) {
      iy = (int)(4000*meanvdf[i]);
	  iy2 = (int)(4000*meanvdf[i+1]);
      if (iy>0) {
        gOff.drawLine(gx+i,gy+200-iy,gx+i,gy+200-iy2);
      }
    }

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

    gOff.setColor(Color.lightGray);
    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.getHSBColor(0.6f,0.8f,0.8f));
    gOff.drawString("count",gx+120,gy+20);
    gOff.setColor(Color.getHSBColor(0.15f,0.8f,0.6f));
    gOff.drawString("mean distr.",gx+120,gy+40);
    gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.6f));
    gOff.drawString("Maxwell(Tc)",gx+120,gy+60);
  }

  void setVelocityDistribution() {
    int i,iv;
    double v2,vv;

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

    for (i=0; i<Nmt; i++) {
      v2 = vx[i]*vx[i]+vy[i]*vy[i]+vz[i]*vz[i];
      vv = Math.sqrt(v2);
      if (vv<2000.0) {
        iv = (int)(vv*0.1);
        vdf[iv] +=1;
      }
    }

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

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

    m = FS[MAT1][0];
    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 -----------------*/

  void drawEffectivePotential(int gbx, int gby) {
    int i,ip,ipMax,ipMin,v1,v2;
    double ra,fp,fpMax,fpMin;

    ipMax=0; ipMin=0; fpMax=-1.0e99; fpMin=1.0e99;
    for (ip=0;ip<Nmt;ip++) {
      fp = dFdroh(ip);
      if (fp>fpMax) {
        fpMax = fp; ipMax = ip;
      }
      if (fp<fpMin) {
        fpMin = fp; ipMin = ip;
      }
    }

    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("5eV",gbx-30,gby+200-125+5);
    gOff.drawString("effective potential",gbx+30,gby+300);

    gOff.setColor(Color.getHSBColor(0.10f,0.8f,0.9f));
    gOff.drawString("max Veff",gbx+80,gby+80);
    for(i=0;i<200;i++) {
      ra = 0.02*i+1.0;
      v1 = (int)(25.0*effectivePotential(ra,fpMax)); if (v1>220) continue;
	  v2 = (int)(25.0*effectivePotential(ra+0.02,fpMax));
      if (v1<220) {
        gOff.drawLine(gbx+i,gby+200-v1, gbx+i, gby+200-v2);
      }
    }

    gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.9f));
    gOff.drawString("V(r)",gbx+80,gby+40);
    for(i=0;i<200;i++) {
      ra = 0.02*i+1.0;
      v1 = (int)(25.0*potV(ra));if (v1>220) continue;
	  v2 = (int)(25.0*potV(ra+0.02));
      if (v1<220) {
        gOff.drawLine(gbx+i,gby+200-v1, gbx+i, gby+200-v2);
      }
    }

    gOff.setColor(Color.getHSBColor(0.33f,0.8f,0.9f));
    gOff.drawString("-2AF'fai(r)",gbx+80,gby+60);
    for(i=0;i<200;i++) {
      ra = 0.02*i+1.0;
      v1 = (int)(25.0*embeddedPotential(ra,fpMin));if (v1<-50) continue;
	  v2 = (int)(25.0*embeddedPotential(ra+0.02,fpMin));
      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("effective potential",gbx+80,gby+20);
    for(i=0;i<200;i++) {
      ra = 0.02*i+1.0;
      v1 = (int)(25.0*effectivePotential(ra,fpMin)); if (v1>220) continue;
	  v2 = (int)(25.0*effectivePotential(ra+0.02,fpMin));
      if (v1<220) {
        gOff.drawLine(gbx+i,gby+200-v1, gbx+i, gby+200-v2);
      }
    }
  }

  double effectivePotential(double ra, double fp) {
    return( potV(ra)-2.0*FSA*fp*fai(ra) );
  }

  double embeddedPotential(double ra, double fp) {
    return( -2.0*FSA*fp*fai(ra) );
  }

/*-------------------------------------  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 ( s/((double)Nmt*1.5*1.38e-23) );
  }

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

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

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

/*-------------------- virial  ----------*/

  void initVirial() {
    int i;

    for (i=0; i<256; i++) {
      virial[i] = 0.0;
    }
    virialp = 0;
  }

  void setVirial() {
    int i;
    double vir;

    vir = 0.0;
    for (i=0; i<Nmt; i++) {
      vir += virx[i]+viry[i]+virz[i];
    }

    virial[virialp] = vir;
    virialp = (virialp+1)%256;
  }

  double pressure() {
    int i;
    double s,vir,tm,vol,p;

    s = 0.0;
    for (i=0; i<256; i++) {
      s += virial[i];
    }
    vir = s/256.0;

    tm = temperature();
    vol = xMax*yMax*zMax;
    p = (Nmt*1.38e-23*tm+vir/3.0)/vol;
    return( p );
  }


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

  void setInitialCondition() {
    int i,j,ix,iy,iz,n2,n3; 
    double a,d;

    t = 0.0;

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

    FSbb = Ackland[MAT1][0];
    FSalpha = Ackland[MAT1][1];
    FSb0 = Ackland[MAT1][2];

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

    d = 5.0e-10;


    n2 = Ncube*Ncube; n3 = Ncube*Ncube*Ncube;
    for (i=0; i<n3; i++) {
      ix = i/n2;
      iy = (i%n2)/Ncube;
      iz = (i%n2)%Ncube;
      xx[i] = a*(double)(ix)+d; //-3.0*d; if (xx[i]<0.0) xx[i] += xMax;
      yy[i] = a*(double)(iy)+d;
      zz[i] = a*(double)(iz)+d;
      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;
    }
    for (i=n3; i<Nmt; i++) {
      j = i-n3;
      ix = j/n2;
      iy = (j%n2)/Ncube;
      iz = (j%n2)%Ncube;
      xx[i] = a*(double)(ix)+0.5*a+d; //-3.0*d; if (xx[i]<0.0) xx[i] += xMax;
      yy[i] = a*(double)(iy)+0.5*a+d;
      zz[i] = a*(double)(iz)+0.5*a+d;
      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;
    }

    initVirial();
  }


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

  void timeEvolution(){
    int i;

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

    if (tempMode==1) {
      vAjustment();
    }

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

  void timeStep() {
    int i; 
    double dtv2, a2, rr, tmp;

    dtv2 = dt/2.0;
    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++) {
      a2 = dtv2/FSmass;
      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;
    }
    forceCalc();
    for (i=0; i<Nmt; i++) {
      a2 = dtv2/FSmass;
      vx[i] += a2*ffx[i];
      vy[i] += a2*ffy[i];
      vz[i] += a2*ffz[i];
    }

    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;
    double xi,xj,yi,yj,zi,zj;

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

      virx[i] = 0.0;
      viry[i] = 0.0;
      virz[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];

          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;

          virx[i] +=  fxij*(xi-xj);
          viry[i] +=  fyij*(yi-yj);
          virz[i] +=  fzij*(zi-zj);

        }
      }
    }
  }

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

  void regNear() {
    int ip,jp,kp,iq;
    int i,j,k,ii,jj,kk,i0,i1,j0,j1,k0,k1;
    double rrg,rrg2,r2;
    double xi,xj,yi,yj,zi,zj;

    preRegistration();

    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(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) {
                xi = xx[ip]; xj = xx[jp]; yi = yy[ip]; yj = yy[jp]; zi = zz[ip]; zj = zz[jp];
                r2=(xi-xj)*(xi-xj)+(yi-yj)*(yi-yj)+(zi-zj)*(zi-zj);
                if (r2<rrg2) {
                  regNr[ip][kp]=jp;
                  kp += 1;
                }
              }
            }
          }
        }
      }
      regNr[ip][0]=kp;
    }

  }

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

  }

/*--------------------  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;
    double xi,xj,yi,yj,zi,zj;

    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];
        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,ret;

    rh = roh(i);
    ret = 0.0;
    if (rh>0.0) { 
      ret = 0.5*Math.sqrt(1.0/rh);
    }
    return ( ret );
  }

  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) {
      if (r>=FSb0) {
        return ( (r-FSc)*(r-FSc)*(FSc0+FSc1*r+FSc2*r*r) );
      } else {
        return ( (r-FSc)*(r-FSc)*(FSc0+FSc1*r+FSc2*r*r) +FSbb*(FSb0-r)*(FSb0-r)*(FSb0-r)*Math.exp(-FSalpha*r) );
      }
    } else {
      return ( 0.0 );
    }
  }

  double dpotVdr(double r) {
    if (r<FSc) {
      if (r>=FSb0) {
        return ( 2.0*(r-FSc)*(FSc0+FSc1*r+FSc2*r*r)+(r-FSc)*(r-FSc)*(FSc1+2.0*FSc2*r) );
      } else {
        return ( 2.0*(r-FSc)*(FSc0+FSc1*r+FSc2*r*r)+(r-FSc)*(r-FSc)*(FSc1+2.0*FSc2*r) + FSbb*(FSb0-r)*(FSb0-r)*Math.exp(-FSalpha*r)*(-3.0-FSalpha*(FSb0-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 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 -----*/

/* Finnis-Sinclair potential
 *
 * original paper
 * M. W. Finnis and J. E. Sinclair ; Phil. Mag. A, Vol. 50, No. 1, p.45-55 (1984)
 * core correction 
 * G. J. Ackland and R. Thetford ; Phil. Mag. A, Vol. 56, No. 1, p.15-30 (1987)
 *
 * total energy : E_tot
 *   E_tot = (1/2) Sum[ V(r_ij), {i,j} ] - Sum[ F(rho_i), {i} ]
 *
 * pair-potential V(r_ij)
 *   V(r_ij) = (r-c)^2 (c0 + c1 r + c2 r^2), (r<=c else 0)
 * (Ackland core correction)
 *   V(r_ij) = (r-c)^2 (c0 + c1 r + c2 r^2) + B (b0-r)^3 exp(-alpha r), (r<b0<c)
 *
 * n-body (embedding) function F(rho)
 *   F(rho) = A sqrt(rho)
 * electronic density at host atom i
 *   rho_i = Sum[ fai(r_ij), {j} ]
 *
 * electronic density function
 *   fai(r) = (r-d)^2 (r<=d else 0) 
 * (Ackland correction)
 *   fai(r) = (r-d)^2 + B(r-d)^3/d (r<=d else 0) 
 *
 * effective pair-potential V_eff(r)
 *   (F"=0 apploximation)
 *   V_eff(r) = V(r) - 2 A F' fai(r)
 *    ( F' = d F(rho)/drho )
 *
 */