/* thermal expansion NaCl - PSX1S */
/* coded by Ikeuchi Mitsuru       */
/* ver 0.0.1   2003.11.03         */
/* ver 0.0.2   2004.02.11 improved pannel layout */

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

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

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

  Choice ctc,cvw;
  Scrollbar scr,scp;
  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 = 30.0/sgm;

  double xMax0 = 5.0*sgm;       /* x-Box size in (m)      */
  double yMax0 = 5.0*sgm;       /* y-Box size in (m)      */
  double zMax0 = 5.0*sgm;       /* z-Box size in (m)      */
  double xMax = xMax0;
  double yMax = yMax0;
  double zMax = zMax0;

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

  double contPress = 1.0*1.0e6;
  int tempMode = 1;
  double contTemp = 300.0;
  double scale = 1.0;
  int dispMode = 0;


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

  int dgX, dgY, dgXb,dgYb;

  double cx = 0.5*xMax;
  double cy = 0.5*yMax;
  double cz = 0.5*zMax;
  double theta = -10.0*3.14/180.0;
  double fai = 5.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 eForceConst = 1.0/(4.0*3.1416*8.8542e-12);

  double SXdata[][] = {
  /*            0 M   ,     1 z ,    2 a ,    3 b ,   4 c ,  5 r  */
  /* 0 B  */ { 10.81*AU,  3.0*EE, 0.720*AA, 0.080*AA,  0.0, 0.23 },
  /* 1 O  */ { 16.00*AU, -2.0*EE, 1.626*AA, 0.085*AA, 20.0, 1.40 },
  /* 2 F  */ { 19.00*AU, -1.0*EE, 1.565*AA, 0.085*AA, 20.0, 1.33 },
  /* 3 Na */ { 22.99*AU,  1.0*EE, 1.260*AA, 0.080*AA, 20.0, 1.02 },
  /* 4 Mg */ { 24.31*AU,  2.0*EE, 1.161*AA, 0.080*AA, 10.0, 0.72 },
  /* 5 Al */ { 26.98*AU,  3.0*EE, 1.064*AA, 0.080*AA,  2.0, 0.53 },
  /* 6 Si */ { 28.09*AU,  4.0*EE, 1.012*AA, 0.080*AA,  0.0, 0.40 },
  /* 7 Cl */ { 35.45*AU, -1.0*EE, 1.950*AA, 0.090*AA, 30.0, 1.81 },
  /* 8 K  */ { 39.10*AU,  1.0*EE, 1.595*AA, 0.080*AA, 15.0, 1.38 },
  /* 9 Ca */ { 40.08*AU,  2.0*EE, 1.414*AA, 0.080*AA, 10.0, 1.00 },
  /*10 Ti */ { 47.88*AU,  4.0*EE, 1.235*AA, 0.080*AA,  0.0, 0.61 },
  /*11 Sr */ { 87.62*AU,  2.0*EE, 1.632*AA, 0.080*AA, 15.0, 1.16 },
  /*12 Ba */ { 137.3*AU,  2.0*EE, 1.820*AA, 0.080*AA, 20.0, 1.36 }
  };

  int NmatKind = 13;
  double SXa[][] = new double[NmatKind][NmatKind];
  double SXb[][] = new double[NmatKind][NmatKind];
  double SXc[][] = new double[NmatKind][NmatKind];

  int MAT0 = 3;  /* Na */
  int MAT1 = 7;  /* Cl  */
  int MAT2 = 3;  /* Na */
  int MAT3 = 9;  /* Ca */
  int MAT4 = 7;  /* Cl */

  int kindT[] = new int[NmatKind];
  double forceTable[][][] = new double[5][5][1200];
  
  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 */

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

  int rdf[][] = new int[4][120];

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

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

    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("temp");
    cvw.addItemListener(this);

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

    scp= new Scrollbar(Scrollbar.HORIZONTAL,1,10,1,310);
    scp.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("  NaCl crystal"));
    pnl.add(ctc);
    pnl.add(cvw);
    pnl.add(new Label("  "));
    pnl.add(scr);
    pnl.add(scp);
    add(pnl,"North");

    setSX();
    setInitialPosition();
  }

  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() == ctc) {
      tempMode = ctc.getSelectedIndex();
      repaint(); 
    } else if (ev.getSource() == cvw){
      dispMode = cvw.getSelectedIndex();
      repaint(); 
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == scr) { 
      contTemp = 10.0*(double)(scr.getValue());
      vAjustment();
      repaint();
    } else if(ev.getSource() == scp) { 
      contPress = 1.0e6*(double)(scp.getValue());
      repaint();
    }
  }

  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,tmp,clr,tm;

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

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

    gbx = 80;
    gby = 100;
    sc = dispScale*scale;
    sz = 1.0*sc*1.0e-10*scale;

    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)*SXdata[kind[j]][5] )+2;
        clr = colorOf(kind[j]);
        gOff.setColor(Color.getHSBColor((float)clr,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)<3.3e-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;
        for (k=1; k<reg[j][0]; k++) {
          jj = reg[j][k];
          if (distance(j,jj)<3.3e-10) {
            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*(0.4*pz[j]/zMax+0.6)*SXdata[kind[j]][5] )+2;
        clr = colorOf(kind[j]);
        gOff.setColor(Color.getHSBColor((float)clr,0.95f,(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)*SXdata[kind[j]][5] );
        clr = colorOf(kind[j]);
        gOff.setColor(Color.getHSBColor((float)clr,0.95f,(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)*SXdata[kind[j]][5] )+2;
        tm = SXdata[kind[j]][0]*(vx[j]*vx[j]+vy[j]*vy[j]+vz[j]*vz[j])/(3.0*1.38e-23)/1000.0;
        if (tm>1.0) {
          tm = 1.0;
        }
        clr = colorOf(kind[j]);
        gOff.setColor(Color.getHSBColor((float)clr,0.95f,(float)(0.8*tm+0.2)));
        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<100; i++) {
      if (rdf[0][i]>0) {
        gOff.drawLine(gbx+i,gby+200-rdf[0][i], gbx+i, gby+200);
      }
    }


    gOff.setColor(Color.red);
    gOff.drawString("Na+",20,40);
    gOff.setColor(Color.blue);
    gOff.drawString("Cl-",50,40);
    gOff.setColor(Color.black);
    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("cont P="+(int)(contPress/1.0e6+0.5)+" MPa",630/6*5+10,40);

    gb = 630/6*5;
    gOff.setColor(Color.blue);
    gOff.drawString("t="+Math.floor(t*1.0e15)+"(fs)",gb,110);
    gOff.drawString("size="+(int)(xMax*1.0e12+0.5)+" pm",gb,130);
    gOff.drawString("T="+temperature()+" K",gb,150);
    gOff.drawString("P="+(int)(pressure()/1.0e6)+"MPa",gb,170);

  }

  double colorOf(int knd) {
    double c;

    c = 0.0;
    if (knd==MAT0) {
      c = 0.01;
    } else if (knd==MAT1) {
      c = 0.66;
    } else if (knd==MAT2) {
      c = 0.33;
    } else if (knd==MAT3) {
      c = 0.13;
    }
    return ( c );
  }

  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 marume(double x) {

    return(Math.floor(x*100.0)/100.0);
  }

/*----------------------------------  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++) {
      srtz[0][i] = i;
      srtz[1][i] = (int)(255.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 rmsV() {
    int i;
    double s;

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i];
    }
    return ((double)Math.floor(Math.sqrt(s/(double)Nmt)*10.0)/10.0);
  }

  double temperature() {
    int i;
    double s;

    s = 0.0; 
    for (i=0; i<Nmt; i++) {
      s += 0.5*(SXdata[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 density() {

    return ( 64.0*1.0e-28/(xMax*yMax*zMax) );
  }

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

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

    for (i=0; i<120; i++) {
      rdf[0][i]=0; rdf[1][i]=0; rdf[2][i]=0; rdf[3][i]=0;
    }

    for (i=0; i<Nmt; i++) {
        for(k=1;k<reg[i][0];k++) {
        j = reg[i][k];
        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);
        rij = Math.sqrt(r2);
        if (rij<5.0e-10) {
          ir = (int)(Math.floor(rij*2.0e11+0.49999));
          rdf[0][ir] +=1;
          if (kind[i]==MAT0 && kind[j]==MAT0) {
            rdf[1][ir] +=1;
          } else if (kind[i]==MAT1 && kind[j]==MAT1) {
            rdf[2][ir] +=1;
          } else  {
            rdf[3][ir] +=1;
          }
        }
      }
    }

  }

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

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

    gAjustment();
    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 gAjustment() {
    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++) {
      vx[i] = vx[i]-mvx;
      vy[i] = vy[i]-mvy;
      vz[i] = vz[i]-mvz;
    }
  }


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

  void initVirial() {
    int i;

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

  void setVirial() {
    int i,j,k;
    double xi,xj,yi,yj,zi,zj;
    double r2,rij,f,fxij,fyij,fzij;
    double xVir,yVir,zVir;

    xVir = 0.0; yVir = 0.0; zVir = 0.0;

    for(i=0;i<Nmt;i++) {
      for(k=1;k<reg[i][0];k++) {
        j = reg[i][k];
        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 = force(rij,kind[i],kind[j]);
          fxij = f*(xi-xj)/rij;
          fyij = f*(yi-yj)/rij;
          fzij = f*(zi-zj)/rij;

          xVir += fxij*(xi-xj);
          yVir += fyij*(yi-yj);
          zVir += fzij*(zi-zj);
        }
      }
    }

    virial[virialp] = xVir+yVir+zVir;
    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 );
  }

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

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

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

/*-------------------------------  vol control   -------------*/

  void volControl() {
    double d,pr,r,x,y,z;

    d = dt*30.0;
    pr = pressure1();
    if (pr<0.0) {
      pr = 0.001;
    }
    r = Math.pow(pr/contPress,0.3333);

    x = r*xMax;
    y = r*yMax;
    z = r*zMax;

    if (x>xMax+d) {
      x = xMax + d;
    } else if (x<xMax-d) {
      x = xMax - d;
    }
    xMax = x;

    if (y>yMax+d) {
      y = yMax + d;
    } else if (y<yMax-d) {
      y = yMax - d;
    }
    yMax = y;

    if (z>zMax+d) {
      z = zMax + d;
    } else if (z<zMax-d) {
      z = zMax - d;
    }
    zMax = z;

    cx = 0.5*xMax;
    cy = 0.5*yMax;
    cz = 0.5*zMax;

  }



/*----------------------  set SX data and ForceTable   ---------*/

  void setSX(){
    int i,j; 

    for (i=0; i<NmatKind; i++) {
      for (j=0; j<NmatKind; j++) {
        SXa[i][j] = SXdata[i][2]+SXdata[j][2];
        SXb[i][j] = SXdata[i][3]+SXdata[j][3];
        SXc[i][j] = SXdata[i][4]+SXdata[j][4];
      }
    }

  }

  void setForceTable(){
    int i,j,ir,k0,k1; 
    double a,b,c,zi,zj,fc,fr,r,sc;
    int kkk[] = {MAT0,MAT1,MAT2,MAT3,MAT4};

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

        a = SXa[k0][k1];
        b = SXb[k0][k1];
        c = SXc[k0][k1];
        zi = SXdata[k0][1];
        zj = SXdata[k1][1];    

        for (ir=10; ir<1000; ir++) {
          r = (double)ir*1.0e-12;
          fc = eForceConst*zi*zj*(Math.exp(-r/6.5e-10)/r)*(1.0/r+1.0/6.5e-10);
          fr = 6.9742e-11*Math.exp((a-r)/b);
          forceTable[i][j][ir]=(fc + fr);
        }

        for (ir=0; ir<10; ir++) {
          forceTable[i][j][ir]=forceTable[i][j][10];
        }

        for (ir=600; ir<1000; ir++) {
          forceTable[i][j][ir] = cutoff(ir)*forceTable[i][j][ir];
        }
      }
    }

  }

  double cutoff(int ir) {
    double rt;

    rt = 0.0;
    if (ir<600) {
      rt = 1.0;
    } else if (ir>=600 && ir<1000) {
      rt = 0.5+0.5*Math.cos(3.1415926536*(((double)ir)/1000.0-0.6)/0.4);
    } else {
      rt = 0.0;
    }
    return( rt );
  }

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

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

    kindT[MAT0] = 0; kindT[MAT1] = 1; kindT[MAT2] = 2;
    kindT[MAT3] = 3; kindT[MAT4] = 4;
    setForceTable();

    a = 0.5*5.64e-10;
    ip = 0;
    for (ix=0; ix<6; ix++) {
      for (iy=0; iy<6; iy++) {
        for (iz=0; iz<6; iz++) {
          kind[ip] = MAT0;
          if ((ix+iy+iz)%2==1 ) {
            kind[ip] = MAT1;
          }
          xx[ip] = a*(double)ix+0.5*a;
          yy[ip] = a*(double)iy+0.5*a;
          zz[ip] = a*(double)iz+0.5*a;

          vx[ip] = 200.0*Math.random()-100.0;
          vy[ip] = 200.0*Math.random()-100.0;
          vz[ip] = 200.0*Math.random()-100.0;

          ffx[ip] = 0.0;
          ffy[ip] = 0.0;
          ffz[ip] = 0.0;
          fbx[ip] = 0.0;
          fby[ip] = 0.0;
          fbz[ip] = 0.0;

          ip += 1;
        }
      }
    }
    xMax = 6*a;
    yMax = 6*a;
    zMax = 6*a;

    initVirial();
  }

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

  void registration() {
    int i,j,k,kk;
    double xi,xj,yi,yj,zi,zj;
    double r2,rreg,rreg2;

    rreg = rc+10000.0*dt*10;
    rreg2 = rreg*rreg;
    for(i=0;i<Nmt;i++) {
      k = 1;kk=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) {
          reg[i][k]=j;
          k += 1;
        }
      }
      reg[i][0]=k;
    }
  }

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

  void calcposition(){
    int i;

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

    registration();
    for (i=0; i<10; i++) {
      timeEvolution();
      volControl();
    }
    setVirial();

  }

  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]/SXdata[kind[i]][0];
      yy[i] += vy[i]*dt+0.5*dt*dt*ffy[i]/SXdata[kind[i]][0];
      zz[i] += vz[i]*dt+0.5*dt*dt*ffz[i]/SXdata[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])/SXdata[kind[i]][0];
      vy[i] += 0.5*dt*(ffy[i]+fby[i])/SXdata[kind[i]][0];
      vz[i] += 0.5*dt*(ffz[i]+fbz[i])/SXdata[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,rij,f,fxij,fyij,fzij;

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

    for(i=0;i<Nmt;i++) {
      for(k=1;k<reg[i][0];k++) {
        j = reg[i][k];
        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 = force(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 force(double r, int ki, int kj) {
    int i;
    double a,hh;

    hh = 1.0e-12;
    i = (int)(r/hh);
    a = r - i*hh;

    return ( ((hh-a)*forceTable[kindT[ki]][kindT[kj]][i] + 
      a*forceTable[kindT[ki]][kindT[kj]][i+1])/hh );
  }

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

}  

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