/* La on W(100) slab 3x3x3 - y,z-periodic, V=0,psi=0 at x boundary */ 
/* external charge density - real space DFT-LDA */
/* LDA:  J. P. Perdew and A. Zunger; Phys. Rev., B23, 5048 (1981) */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2005.01.16       */

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

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

/* -------------------------------------- set global ------ */

  Choice ch_atom,ch_noe,cvw, ch_md;
  Button bt_reset, bt_start, bt_stop;
  Scrollbar sc_dist,sc_mixing,sc_broadening,sc_orbit;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int dispMode = 10;
  int dispObj = 76;

  double ppDistance = 6.0;
  double mixing = 0.2;
  double broadening = 0.001;
  double SDdump = 0.2;

  int iterCount = 0;
  int started = 1;

  double NOE = 165.0;
  double Znuc = 165.0;
  int Norbit = 92;
  int setVhQ = 0;

  int dgX, dgY, dgXb,dgYb;
  double viewTheta = -15.0*3.14/180.0;
  double viewFai = -72.0*3.14/180.0;
  double dtheta = 0.0*3.14/180.0;
  double pai = 3.1415926536;

  double cosTh = Math.cos(viewTheta); 
  double sinTh = Math.sin(viewTheta);
  double cosFi = Math.cos(viewFai);
  double sinFi = Math.sin(viewFai);


  double t = 0.0;

  double dx = 1.0;
  double dy = 1.0;
  double dz = 1.0;
  double dt = 1.0*(dx*dx);

  int NNx = 32;
  int NNy = 18;
  int NNz = 18;
  int NN = NNx+1; /* max(NNx,NNy,NNz)+1 */

  double vv[][][] = new double[NNx][NNy][NNz];

  double vvext[][][] = new double[NNx][NNy][NNz];
  double vvh[][][] = new double[NNx][NNy][NNz];
  double vvx[][][] = new double[NNx][NNy][NNz];
  double vvc[][][] = new double[NNx][NNy][NNz];

  double rhoExt[][][] = new double[NNx][NNy][NNz];
  double rho[][][] = new double[NNx][NNy][NNz];

  int NNs = 92;
  double sdEnergy[] = new double[NNs];
  double sdState[][][][] = new double[NNs][NNx][NNy][NNz];
  double occ[] = new double[NNs];

  double nuc[][] = {
    { 3.0, 12.0+6.0,  9.0,  9.0  },

    { 6.0,  6.0,  0.0,  0.0  },
    { 6.0,  6.0,  0.0,  6.0  },
    { 6.0,  6.0,  0.0, 12.0  },
    { 6.0,  6.0,  6.0,  0.0  },
    { 6.0,  6.0,  6.0,  6.0  },
    { 6.0,  6.0,  6.0, 12.0  },
    { 6.0,  6.0, 12.0,  0.0  },
    { 6.0,  6.0, 12.0,  6.0  },
    { 6.0,  6.0, 12.0, 12.0  },
  
    { 6.0,  9.0,  3.0,  3.0  },
    { 6.0,  9.0,  3.0,  9.0  },
    { 6.0,  9.0,  3.0, 15.0  },
    { 6.0,  9.0,  9.0,  3.0  },
    { 6.0,  9.0,  9.0,  9.0  },
    { 6.0,  9.0,  9.0, 15.0  },
    { 6.0,  9.0, 15.0,  3.0  },
    { 6.0,  9.0, 15.0,  9.0  },
    { 6.0,  9.0, 15.0, 15.0  },

    { 6.0, 12.0,  0.0,  0.0  },
    { 6.0, 12.0,  0.0,  6.0  },
    { 6.0, 12.0,  0.0, 12.0  },
    { 6.0, 12.0,  6.0,  0.0  },
    { 6.0, 12.0,  6.0,  6.0  },
    { 6.0, 12.0,  6.0, 12.0  },
    { 6.0, 12.0, 12.0,  0.0  },
    { 6.0, 12.0, 12.0,  6.0  },
    { 6.0, 12.0, 12.0, 12.0  },

    { 0.0,  0.0 , 0.0 , 0.0  }
  };

  int pxBox[] = new int[8];
  int pyBox[] = new int[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 xpts[] = new int[NN];
  int ypts[] = new int[NN];

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

  public void init() {

    resize(630,320);
    setBackground(Color.white);

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

    ch_atom = new Choice();
    ch_atom.addItem("2");ch_atom.addItem("3");
    ch_atom.addItem("4");ch_atom.addItem("5");
    ch_atom.addItem("6");ch_atom.addItem("7");
    ch_atom.addItem("8");ch_atom.addItem("9");
    ch_atom.addItem("10");ch_atom.addItem("11");
    ch_atom.addItem("12");ch_atom.addItem("13");
    ch_atom.addItem("14");ch_atom.addItem("15");
    ch_atom.addItem("16");ch_atom.addItem("17");
    ch_atom.addItem("18");ch_atom.addItem("19");
    ch_atom.addItem("20");ch_atom.addItem("21");
    ch_atom.addItem("22");
    ch_atom.addItemListener(this);
    ch_atom.select("20");

    ch_noe = new Choice();
    ch_noe.addItem("ne=Z-1");ch_noe.addItem("ne=Z-0.5");
    ch_noe.addItem("ne=Z");ch_noe.addItem("ne=Z+0.5");ch_noe.addItem("ne=Z+1");
    ch_noe.addItemListener(this);
    ch_noe.select("ne=Z");

    cvw = new Choice();
    cvw.addItem("ph-0");cvw.addItem("ph-1");cvw.addItem("ph-2");
    cvw.addItem("ph-3");cvw.addItem("ph-4");cvw.addItem("ph-5");
    cvw.addItem("ph-6");cvw.addItem("ph-7");cvw.addItem("ph-8");
    cvw.addItem("ph-9");cvw.addItem("ph-10");cvw.addItem("ph-11");
    cvw.addItem("ph-12");cvw.addItem("ph-13");cvw.addItem("ph-14");
    cvw.addItem("ph-15");
    cvw.addItemListener(this);
    cvw.select("ph-0");

    ch_md = new Choice();
    ch_md.addItem("density");ch_md.addItem("phase");ch_md.addItem("ph-grid");
    ch_md.addItem("Vext/4"); ch_md.addItem("VH/4"); ch_md.addItem("Vx");
    ch_md.addItem("Vc"); ch_md.addItem("Veff"); ch_md.addItem("rho");
    ch_md.addItem("rho 3D");ch_md.addItem("rho+charge");
    ch_md.addItemListener(this);
    ch_md.select("rho+charge");


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

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

    bt_stop = new Button("stop");
    bt_stop.addActionListener(this);

    sc_dist= new Scrollbar(Scrollbar.HORIZONTAL,6,10,0,20);
    sc_dist.addAdjustmentListener(this);

    sc_mixing= new Scrollbar(Scrollbar.HORIZONTAL,20,10,1,100);
    sc_mixing.addAdjustmentListener(this);

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

    sc_orbit= new Scrollbar(Scrollbar.HORIZONTAL,76,10,0,101);
    sc_orbit.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(sc_dist);
/*
    pnl.add(new Label("  "));
*/
    pnl.add(ch_noe);
    pnl.add(sc_mixing);
    pnl.add(sc_broadening);
    pnl.add(sc_orbit);
    pnl.add(ch_md);
    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_atom){
      Znuc = 2.0 + ch_atom.getSelectedIndex();
	  ch_noe.select("ne=Z");
	  NOE = Znuc;
	  iterCount = 0;
	  setVhQ = 0;
	  setPotential();
    } else if (ev.getSource() == ch_noe){
      NOE = Znuc - 1.0 + 0.5*ch_noe.getSelectedIndex();
/*
	  iterCount = 0;
	  setVhQ = 0;
	  setPotential();
*/
    } else if (ev.getSource() == cvw){
      dispObj = cvw.getSelectedIndex();
    } else if (ev.getSource() == ch_md){
      dispMode = ch_md.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      started = -1;
    } else if(ev.getSource() == bt_start){ 
      started = 1;
    } else if(ev.getSource() == bt_stop){
      started = 0;
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_dist) {
	  /* iterCount = 0; */
      ppDistance = dx*sc_dist.getValue();
      nuc[0][1] = 12.0+ppDistance;
	  setRhoExt(ppDistance);
    } else if (ev.getSource() == sc_broadening) {
	  broadening = 0.001*sc_broadening.getValue();
    } else if (ev.getSource() == sc_mixing) { 
      mixing = 0.01*(double)(sc_mixing.getValue());
    } else if (ev.getSource() == sc_orbit) { 
      dispObj = sc_orbit.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;
    }
    viewTheta += 0.5*3.14/180.0*(dgX-dgXb);
    viewFai += 0.5*3.14/180.0*(dgY-dgYb);
    cosTh = Math.cos(viewTheta); 
    sinTh = Math.sin(viewTheta);
    cosFi = Math.cos(viewFai);
    sinFi = Math.sin(viewFai);
  }

  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 i,gb,iorbit;
    double ev,ek;
    double dens;

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

    dens = 1500.0;

    if (dispMode==0 || dispMode==1) {

      boxPlot();
      dens3DPlotSD(dens,dispMode,dispObj);
      boxPlot2();

    } else if (dispMode==2) {

      gridFnPlot(sdState[dispObj], NNz/2, 50.0);

    } else if (dispMode==3) {

      gridFnPlot(vvext, NNz/2, 0.25);

    } else if (dispMode==4) {

      gridFnPlot(vvh, NNz/2, 0.25);

    } else if (dispMode==5) {

      gridFnPlot(vvx, NNz/2, 1.0);

    } else if (dispMode==6) {

      gridFnPlot(vvc, NNz/2, 1.0);

    } else if (dispMode==7) {

      gridFnPlot(vv, NNz/2, 1.0);

    } else if (dispMode==8) {

      gridFnPlot(rho, NNz/2, 50.0);

    } else if (dispMode==9) {

      boxPlot();
      rho3DPlot(10.0);
      boxPlot2();

    } else if (dispMode==10) {

      boxPlot();
      charge3DPlot(10.0);
      boxPlot2();

    }

    gOff.setColor(Color.black);
    gOff.drawString("W|-La="+(int)(ppDistance*10.0+0.5)/10.0+" ",630/6*0+10,40);
    gOff.drawString("Z="+(int)(Znuc+0.5)+", ne="+(int)(NOE*10+0.5)/10.0+" ",630/6*1+0,40);
    gOff.drawString("mix="+(int)(mixing*100.0+0.5)/100.0+" ",630/6*2+30,40);
    gOff.drawString("kT0="+(int)(broadening*1000+0.5)/1000.0+" ",630/6*3+10,40);
    gOff.drawString("orbit="+dispObj+" ",630/6*4+10,40);
    gOff.drawString("view",630/6*5+20,40);

    gb = 400;
    gOff.setColor(Color.blue);
    gOff.drawString("kT="+(float)(levelWidth())+" au",gb,60);
    gOff.drawString("iteration="+iterCount+" ",gb+120,60);
    gOff.setColor(Color.black);
    gOff.drawString("box="+(int)(NNx*dx*10.0)/10.0+"x"+ +(int)(NNy*dy*10.0)/10.0+"x"+(int)(NNz*dz*10.0)/10.0+", dx="+(int)(dx*1000.0)/1000.0+" au",gb,80);
    gOff.setColor(Color.blue);
	gOff.drawString("total energy="+(float)totalEnergy()+" au",gb,100);
	
    gOff.setColor(Color.black);
    gOff.drawString("No.",gb,120);
    gOff.drawString("occupation",gb+40,120);
    gOff.drawString("energy",gb+130,120);
	iorbit = dispObj;
	if (iorbit>76) iorbit=76;
    for (i=iorbit; (i<Norbit && i<iorbit+16); i++) {
      gOff.setColor(Color.getHSBColor((float)(0.66-0.33*occ[i]), 0.9f,0.8f));
      gOff.drawString(""+i+"",gb,135+12*(i-iorbit));
      gOff.drawString(" "+(float)(occ[i])+"",gb+30,135+12*(i-iorbit));
      gOff.drawString(" "+(float)(sdEnergy[i])+"",gb+130,135+12*(i-iorbit));
    }
  }

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

  void rho3DPlot(double weight) {
    int i,j,k,ix,iy,ir;
    int gbx,gby;
    double cxx,cyy,czz;
    double z,ppx,ppy,ppz;
    double mag,dens,col;

    gbx = 30;
    gby = 80;
    mag = 10.0;
    cxx = mag*(NNx/2.0); 
    cyy = mag*(NNy/2.0);
    czz = mag*(NNz/2.0);

    for (i=0; i<2*NNx-1; i+=1) {
      for(j=0; j<2*NNy-1; j +=1) {
        for (k=0; k<2*NNz-1; k+=1) {

          dens = 1.25*(rho[i/2][j/2][k/2]+rho[i/2+i%2][j/2][k/2]+rho[i/2][j/2+j%2][k/2]+rho[i/2][j/2][k/2+k%2]+rho[i/2+i%2][j/2+j%2][k/2]+rho[i/2+i%2][j/2][k/2+k%2]+rho[i/2][j/2+j%2][k/2+k%2]+rho[i/2+i%2][j/2+j%2][k/2+k%2]);
		  
          ir = (int)(weight*dens+0.8);
          if (ir>0) {
            if (ir>6) {
              ir = 6;
            }
            z = 0.5*k*mag;
            ppy = cosFi*(0.5*j*mag-cyy)+sinFi*(z-czz) + cyy; 
            ppz = -sinFi*(0.5*j*mag-cyy)+cosFi*(z-czz) + czz;
            ppx = cosTh*(0.5*i*mag-cxx)+sinTh*(ppz-czz) + cxx; 
            ppz = -sinTh*(0.5*i*mag-cxx)+cosTh*(ppz-czz) + czz;

            ix = (int)(ppx)+gbx;
            iy = (int)(ppy)+gby;
			
            gOff.setColor(Color.getHSBColor((float)(0.66-0.66*0.2*dens), 0.9f,0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
        }
      }
    }
  }

  void charge3DPlot(double weight) {
    int i,j,k,ix,iy,ir,inuc;
    int gbx,gby;
    double cxx,cyy,czz;
    double z,ppx,ppy,ppz;
    double mag,dens,col;

    gbx = 30;
    gby = 80;
    mag = 10.0;
    cxx = mag*(NNx/2.0); 
    cyy = mag*(NNy/2.0);
    czz = mag*(NNz/2.0);

    for (i=0; i<2*NNx-1; i+=1) {
      for(j=0; j<2*NNy-1; j +=1) {
        for (k=0; k<2*NNz-1; k+=1) {

          dens = 1.25*(rho[i/2][j/2][k/2]+rho[i/2+i%2][j/2][k/2]+rho[i/2][j/2+j%2][k/2]+rho[i/2][j/2][k/2+k%2]+rho[i/2+i%2][j/2+j%2][k/2]+rho[i/2+i%2][j/2][k/2+k%2]+rho[i/2][j/2+j%2][k/2+k%2]+rho[i/2+i%2][j/2+j%2][k/2+k%2]);
          ir = (int)(weight*dens+0.8); if (ir>1) ir = 1;

          z = 0.5*k*mag;
          ppy = cosFi*(0.5*j*mag-cyy)+sinFi*(z-czz) + cyy; 
          ppz = -sinFi*(0.5*j*mag-cyy)+cosFi*(z-czz) + czz;
          ppx = cosTh*(0.5*i*mag-cxx)+sinTh*(ppz-czz) + cxx; 
          ppz = -sinTh*(0.5*i*mag-cxx)+cosTh*(ppz-czz) + czz;

          ix = (int)(ppx)+gbx;
          iy = (int)(ppy)+gby;

          if (ir>0) {
            gOff.setColor(Color.getHSBColor(0.6f, 0.9f,0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
		  if (i%2==0 && j%2==0 && k%2==0 && rhoExt[i/2][j/2][k/2]<0.0) {
		    ir = 3;
			gOff.setColor(Color.getHSBColor(0.01f, 0.9f,0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
        }
      }
    }
  }

  void dens3DPlotSD(double weight,int mode, int istate) {
    int i,j,k,ix,iy,ir;
    int gbx,gby;
    double cxx,cyy,czz;
    double z,ppx,ppy,ppz;
    double mag,dens,col,d0,d1,d2,d3,d4,d5,d6,d7;

    gbx = 30;
    gby = 80;
    mag = 10.0;
    cxx = mag*(NNx/2.0); 
    cyy = mag*(NNy/2.0);
    czz = mag*(NNz/2.0);

    for (i=0; i<2*NNx-1; i+=1) {
      for(j=0; j<2*NNy-1; j +=1) {
        for (k=0; k<2*NNz-1; k+=1) {
		
          d0 = sdState[istate][i/2][j/2][k/2];
		  d1 = sdState[istate][i/2+i%2][j/2][k/2];
		  d2 = sdState[istate][i/2][j/2+j%2][k/2];
		  d3 = sdState[istate][i/2][j/2][k/2+k%2];
		  d4 = sdState[istate][i/2+i%2][j/2+j%2][k/2];
		  d5 = sdState[istate][i/2+i%2][j/2][k/2+k%2];
		  d6 = sdState[istate][i/2][j/2+j%2][k/2+k%2];
		  d7 = sdState[istate][i/2+i%2][j/2+j%2][k/2+k%2];
		  
		  dens = 0.125*(d0*d0+d1*d1+d2*d2+d3*d3+d4*d4+d5*d5+d6*d6+d7*d7);
          ir = (int)(weight*dens+0.8);
          if (ir>0) {
            if (ir>6) {
              ir = 6;
            }
            z = 0.5*k*mag;
            ppy = cosFi*(0.5*j*mag-cyy)+sinFi*(z-czz) + cyy; 
            ppz = -sinFi*(0.5*j*mag-cyy)+cosFi*(z-czz) + czz;
            ppx = cosTh*(0.5*i*mag-cxx)+sinTh*(ppz-czz) + cxx; 
            ppz = -sinTh*(0.5*i*mag-cxx)+cosTh*(ppz-czz) + czz;

            ix = (int)(ppx)+gbx;
            iy = (int)(ppy)+gby;
            col = 0.0;
            if (mode==0) {
              col = 0.66 -0.05*dens;
              if (col<0.0) {
                col = 0.0;
              }
            } else if (mode==1) {
              if (d0+d1+d2+d3+d4+d5+d6+d7>=0) {
                col = 0.33;
              } else {
                col = 0.01;
              }
            }
			
            gOff.setColor(Color.getHSBColor((float)col, (float)(0.3+0.3*occ[istate]),0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
        }
      }
    }
  }

  void gridFnPlot(double fn[][][], int zPos, double weight) {
    int i,j,k,gbx,gby;
    int cx,cy,cz;
    double z,px,py,pz,sc;
    double ph2;

    cx = NNx/2;
    cy = NNy/2;
    cz = 0;
    sc = 10.0;

    gbx = 30+(int)(cx/sc);
    gby = 80+(int)(cy/sc);
    k = zPos;

    for(j=0; j<NNy; j += 1) {
      for (i=0; i<NNx; i++) { 
        z = weight*fn[i][j][k] + cz;
        py = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        pz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
        px = cosTh*(i-cx)+sinTh*(pz-cz) + cx; 
        pz = -sinTh*(i-cx)+cosTh*(pz-cz) + cz;
        xpts[i] = (int)(sc*px)+gbx;
        ypts[i] = (int)(sc*py)+gby;
      }
      for (i=0; i<NNx-1; i++) {
        gOff.setColor(Color.getHSBColor(0.66f,0.6f,0.9f));
        gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
      }
    }

    for(i=0; i<NNx; i += 1) {
      for (j=0; j<NNy; j++) {
        z = weight*fn[i][j][k] + cz;
        py = cosFi*(j-cy)+sinFi*(z-cz) + cy; 
        pz = -sinFi*(j-cy)+cosFi*(z-cz) + cz;
        px = cosTh*(i-cx)+sinTh*(pz-cz) + cx; 
        pz = -sinTh*(i-cx)+cosTh*(pz-cz) + cz;
        xpts[j] = (int)(sc*px)+gbx;
        ypts[j] = (int)(sc*py)+gby;
      }
      for (j=0; j<NNy-1; j++) { 
        gOff.setColor(Color.getHSBColor(0.66f,0.6f,0.9f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }


  void boxPlot() {
    int i,j,k,ip,i1,i2;
    int gbx,gby;
    double cxx,cyy,czz;
    double z,ppx,ppy,ppz;
    double mag,dens;

    gbx = 30;
    gby = 80;
    mag = 10.0;
    cxx = mag*(NNx/2.0); 
    cyy = mag*(NNy/2.0);
    czz = mag*(NNz/2.0);

    ip = 0;
    for (i=0; i<NNx; i+=NNx-1) {
      for(j=0; j<NNy; j+=NNy-1) {
        for (k=0; k<NNz; k+=NNz-1) {
          z = k*mag;
          ppy = cosFi*(j*mag-cyy)+sinFi*(z-czz) + cyy; 
          ppz = -sinFi*(j*mag-cyy)+cosFi*(z-czz) + czz;
          ppx = cosTh*(i*mag-cxx)+sinTh*(ppz-czz) + cxx; 
          ppz = -sinTh*(i*mag-cxx)+cosTh*(ppz-czz) + czz;

          pxBox[ip] = (int)(ppx)+gbx;
          pyBox[ip] = (int)(ppy)+gby;
          ip += 1;
        }
      }
    }
    gOff.setColor(Color.gray);
    for (i=0; i<12; i++) {
      i1= boxp[i][0]; i2= boxp[i][1];
      if (i1==0) {
        gOff.setColor(Color.gray);
        gOff.drawLine(pxBox[i1],pyBox[i1],pxBox[i2],pyBox[i2]);
      }
    }
  }

  void boxPlot2() {
    int i,i1,i2;

    for (i=0; i<12; i++) {
      i1= boxp[i][0]; i2= boxp[i][1];
      if (i1 != 0) {
        gOff.setColor(Color.black);
        gOff.drawLine(pxBox[i1],pyBox[i1],pxBox[i2],pyBox[i2]);
      }
    }
  }

/*-------------------------- statistics ----------*/

  double meanKinetic(double ph[][][]) {
    int i,j,k,ip,im,jp,jm,kp,km;
    double p,d2ph,h2;

    h2 = dx*dx;
    p = 0.0;
    for (i=0;i<NNx;i++) {
      ip = i+1; if (ip>=NNx) ip = 0;
      im = i-1; if (im<0) im = NNx-1;
      for (j=0; j<NNy; j++) {
        jp = j+1; if (jp>=NNy) jp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;
        for (k=0; k<NNz; k++) {
          kp = k+1; if (kp>=NNz) kp = 0;
          km = k-1; if (km<0) km = NNz-1;
          d2ph = (ph[ip][j][k]+ph[im][j][k]+ph[i][jp][k]+ph[i][jm][k]+ph[i][j][kp]+ph[i][j][km]-6.0*ph[i][j][k])/h2;
          p += ph[i][j][k]*d2ph;
        }
      }
    }
    return ( -0.5*p*dx*dy*dz );
  }

  double meanValue(double ph[][][],double a[][][]) {
    int i,j,k;
    double s;

    s = 0.0;
    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          s += a[i][j][k]*ph[i][j][k]*ph[i][j][k];
        }
      }
    }
    return ( s*dx*dy*dz );
  }

  double totalEnergy() {
    int i,j,k;
	double x,y,z;
    double s,sei,enn;

    sei = 0.0;
    for (i=0;i<Norbit;i++) {
      if (occ[i]<=0.0) break;
	  sei += occ[i]*sdEnergy[i];
    }

    enn = 0.0;
    for (i=0;i<10000;i++) {
      if (nuc[i][0]<=0.0) break;
      for (j=i+1;j<10000;j++) {
        if (nuc[j][0]<=0.0) break;
        x = nuc[i][1]-nuc[j][1];
        y = nuc[i][2]-nuc[j][2];
        z = nuc[i][3]-nuc[j][3];
        sei += nuc[i][0]*nuc[j][0]/Math.sqrt(x*x+y*y+z*z);
      }
    }

    s = 0.0;
    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          s += (-0.5*vvh[i][j][k]-0.25*vvx[i][j][k]+eeCorrelation(rho[i][j][k])-vvc[i][j][k])*rho[i][j][k];
        }
      }
    }
	s = s*dx*dy*dz;
	
    return ( enn+sei+s );
  }


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

  void setInitialCondition() {

    t = 0.0;
    iterCount = 0;

    setPotential();
    setOcc(Norbit,NOE);
    initSD(NNs);
  }

/* ----------------------------- set potential -------------- */

  void setPotential() {
    int i,j,k;
    double  x,y,z,r2,s;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
		  vvext[i][j][k] = 0.0;
		  vvh[i][j][k] = 0.0;
		  rho[i][j][k] = 0.0;
		  rhoExt[i][j][k] = 0.0;
        }
      }
    }
    setRhoExt(ppDistance);
  }

  void setRhoExt(double dd) {
    int i,j,k,inuc;
    double  x,y,z,s,x0,y0,z0;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          rhoExt[i][j][k] = 0.0;
        }
      }
    }

    inuc = 0;
    setLaIon(nuc[inuc][1]*dx, nuc[inuc][2]*dy, nuc[inuc][3]*dz);

    for (inuc=1;inuc<10000;inuc++) {
      if (nuc[inuc][0]<=0.0) break;
      setWIon(nuc[inuc][1]*dx, nuc[inuc][2]*dy, nuc[inuc][3]*dz);
    }
  }


  void setWIon(double x0, double y0, double z0) {
    setIon(74.0, 6.0, 2.2, 2.2, x0,y0,z0);
  }

  void setLaIon(double x0, double y0, double z0) {
    setIon(57.0, 3.0, 2.7, 2.3, x0,y0,z0);
  }

  void setIon(double zz, double zv, double rp, double rm, double x0, double y0, double z0) {
    int i,j,k,i0,j0,k0;
    double  x,y,z,r2,sp,sm;

    i0 = (int)(x0/dx+0.5); j0 = (int)(y0/dy+0.5); k0 = (int)(z0/dz+0.5);

    sp = 0.0; sm = 0.0;
    for (i=i0-5;i<=i0+5;i++) {
      x = dx*(i-i0);
      for (j=j0-5;j<=j0+5;j++) {
        y = dy*(j-j0);
        for (k=k0-5;k<=k0+5;k++) {
          z = dz*(k-k0);
          r2 = x*x+y*y+z*z;
          if (r2<rp*rp) {
            sp += 1.0;
          }
		  if (r2<rm*rm) {
			sm += 1.0;
		  }
        }
      }
    }
	
    for (i=i0-5;i<=i0+5;i++) {
      x = dx*(i-i0);
      for (j=j0-5;j<=j0+5;j++) {
        y = dy*(j-j0);
        for (k=k0-5;k<=k0+5;k++) {
          z = dz*(k-k0);
          r2 = x*x+y*y+z*z;
          if (r2<rp*rp) {
            rhoExt[(i+NNx)%NNx][(j+NNy)%NNy][(k+NNz)%NNz] += -zz/sp/(dx*dy*dz);
          }
		  if (r2<rm*rm) {
			rhoExt[(i+NNx)%NNx][(j+NNy)%NNy][(k+NNz)%NNz] += +(zz-zv)/sm/(dx*dy*dz);
		  }
        }
      }
    }
  }

/* ----------------------------- timeEvolution -------------- */

  void timeEvolution() {
    int i;

    if (started==1) {

      iterCount += 1;
      setVeff();
      SD(Norbit, vv, 2, SDdump);
      sortState(Norbit);
      setOcc(Norbit,NOE);

    } else if (started==-1) {

      t = 0.0;
      iterCount = 0;
      setInitialCondition();
      started = 0;

    }
  }

  void sortState(int maxState){
    int ist,i,j,k;
    double w;

    for (ist=maxState-2; ist>=0; ist--) {
      if (sdEnergy[ist]>sdEnergy[ist+1]+0.0001) {
        for (i=0;i<NNx;i++) {
          for (j=0;j<NNy;j++) {
            for (k=0;k<NNz;k++) {
			   w = sdState[ist][i][j][k];
               sdState[ist][i][j][k] = sdState[ist+1][i][j][k];
			   sdState[ist+1][i][j][k] = w;
            }
          }
        }
        w = sdEnergy[ist]; sdEnergy[ist]= sdEnergy[ist+1]; sdEnergy[ist+1]= w;
      }
    }
  }

/* ----------------------------- occupation --------------------- */

  void setOcc(int maxState, double nElectron) {
    int i;
    double eUpper,eLower,eFermi,ntrial;

    eUpper = sdEnergy[maxState-1]+1.0;
    eLower = sdEnergy[0]-1.0;
    for (i=0;i<maxState;i++) {
	  if(sdEnergy[i]>eUpper) eUpper = sdEnergy[i];
	  if(sdEnergy[i]<eLower) eLower = sdEnergy[i];
    }

    while (eUpper-eLower>1.0e-12) {
      eFermi = (eUpper+eLower)/2.0;
	  ntrial = trialOcc(maxState, eFermi);
	  if (ntrial<nElectron) {
	    eLower = eFermi;
	  } else {
        eUpper = eFermi;
	  }
	}
    eFermi = (eUpper+eLower)/2.0;

    for (i=0;i<maxState;i++) {
      occ[i] = 2.0*FermiDirac(sdEnergy[i], eFermi);
	  if (occ[i]<0.0001) occ[i] = 0.0;
      if (2.0-occ[i]<0.0001) occ[i] = 2.0;
	}

  }

  double trialOcc(int maxState, double eFermi) {
    int i;
    double s;

    s = 0.0;
    for (i=0;i<maxState;i++) {
      s += 2.0*FermiDirac(sdEnergy[i], eFermi);
    }
    return ( s );
  }

  double FermiDirac(double ee, double ef) {
    double et;
    et = levelWidth();
    return ( 1.0/(Math.exp((ee-ef)/et)+1.0) );
  }

  double levelWidth() {
    double kT;

    kT = broadening;
	if (iterCount<10) {
	  kT=0.1;
	} else {
      kT = 0.1-0.005*(iterCount-9);
	  if (kT<broadening) kT = broadening;
	}
	return ( kT );
  }

/* ----------------- effective potential <-- electron density ------- */

  void setVeff() {
    int i,j,k;

    setRho();
    poisson(5);
    setVxc();

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          vv[i][j][k] = vvext[i][j][k]+vvh[i][j][k]+vvx[i][j][k]+vvc[i][j][k];
        }
      }
    }
  }

  void setRho() {
    int i,j,k,ie;
    double mx;

    mx = mixing; if (iterCount<20) mx = 1.0;

	for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          rho[i][j][k] *= (1.0-mx);
          for (ie=0; ie<Norbit; ie++) {
			if (occ[ie]<=0.0) break;
            rho[i][j][k] += mx*occ[ie]*(sdState[ie][i][j][k]*sdState[ie][i][j][k]);
          }
        }
      }
    }
  }

  void poisson(int iterMax) {
    int i,j,k,ip,im,jp,jm,kp,km,iter;
    double h2,w;

    h2 = 4.0*3.1416*dx*dx;
    w = 0.2;

	for (iter=1; iter<iterMax; iter++) {
      for (i=1; i<NNx-1; i++) {
	  /*
        ip = i+1; if (ip>=NNx) ip = 0;
        im = i-1; if (im<0) im = NNx-1;*/
        for (j=0; j<NNy; j++) {
          jp = j+1; if (jp>=NNy) jp = 0;
          jm = j-1; if (jm<0) jm = NNy-1;
          for (k=0; k<NNz; k++) {
            kp = k+1; if (kp>=NNz) kp = 0;
            km = k-1; if (km<0) km = NNz-1;
            vvh[i][j][k] += w*(vvh[i+1][j][k]+vvh[i-1][j][k]+vvh[i][jp][k]+vvh[i][jm][k]+vvh[i][j][kp]+vvh[i][j][km]-6.0*vvh[i][j][k] +h2*rho[i][j][k]);

			vvext[i][j][k] += w*(vvext[i+1][j][k]+vvext[i-1][j][k]+vvext[i][jp][k]+vvext[i][jm][k]+vvext[i][j][kp]+vvext[i][j][km]-6.0*vvext[i][j][k] +h2*rhoExt[i][j][k]);

          }
        }
      }
    }
  }

  /* LDA:  J. P. Perdew and A. Zunger; Phys. Rev., B23, 5048 (1981) */

  void setVxc() {
    int i,j,k;
    double c1,rh,rh3,rs,sqrtrs,ec;

    c1 = -0.984745022;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          rh = rho[i][j][k];
          rh3 = Math.pow(rh,0.33333333);
          vvx[i][j][k] = c1*rh3;

          rs = 0.6204/(rh3+1.0e-20);
          if (rs>=1.0) {
		    sqrtrs = Math.sqrt(rs);
            ec = -0.1423/(1.0+1.0529*sqrtrs+0.3334*rs);
			vvc[i][j][k] = ec*(1.0+1.22838*sqrtrs+0.4445*rs)/(1.0+1.0529*sqrtrs+0.3334*rs);
          } else {
            vvc[i][j][k] = -0.05837-0.0084*rs +(0.0311+0.00133*rs)*Math.log(rs);
          }
        }
      }
    }
  }

  double eeCorrelation(double rh) {
    double r,ec;

    r = 0.6204/(Math.pow(rh,0.33333333)+1.0e-20);
    if (r>=1.0) {
      ec = -0.1423/(1.0+1.0529*Math.sqrt(r)+0.3334*r);
    } else {
      ec = -0.0480-0.0116*r+(0.0311+0.0020*r)*Math.log(r);
    }
    return ( ec );
  }


/* ----------------- Gram-Schmidt --------- */

  void GramSchmidt(int stateMax) {
    int i,j,k,ist,istate;
    double s;

    normalize(sdState[0]);

    for (istate=1; istate<stateMax; istate++) {
      for (ist=0; ist<istate; ist++) {
        s = innerProduct(sdState[ist],sdState[istate]);
        for (i=0; i<NNx; i++) {
          for (j=0; j<NNy; j++) {
            for (k=0; k<NNz; k++) {
              sdState[istate][i][j][k] -= s*sdState[ist][i][j][k];
            }
          }
        }
      }
      normalize(sdState[istate]);
    }
  }

  double innerProduct(double f[][][],double g[][][]) {
    int i,j,k;
    double s,dv;

    dv = dx*dy*dz;
    s = 0.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          s += f[i][j][k]*g[i][j][k]*dv;
        }
      }
    }
    return ( s );
  }

  void normalize(double ph[][][]) {
    int i,j,k;
    double s,a,dv;

    dv = dx*dy*dz;
    s = 0.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          s += ph[i][j][k]*ph[i][j][k]*dv;
        }
      }
    }
    a = Math.sqrt(1.0/s);
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        for (k=0; k<NNz; k++) {
          ph[i][j][k] = a*ph[i][j][k];
        }
      }
    }
  }

/* --------------------- steepest descent method -------------- */

  void SD(int stateMax, double v[][][], int iterMax, double dump) {
    int i,ist;

    for (i=0; i<iterMax; i++) {
	  for (ist=0; ist<stateMax; ist++) {
        sdEnergy[ist] = steepestDescent(sdState[ist], v, dump);
      }
	  GramSchmidt(stateMax);
    }
  }

  void initSD(int stateMax) {
    int i,j,k,ist;
    double pai = 3.14159265;

    for (ist=0; ist<stateMax; ist++) {
      for (i=4; i<NNx-4; i++) {
        for (j=0; j<NNy; j++) {
          for (k=0; k<NNz; k++) {
            sdState[ist][i][j][k] = Math.random()-0.5;
          }
        }
      }
	  normalize(sdState[ist]);
    }
  }

  double steepestDescent(double psi[][][], double v[][][], double dump) {
    int i,j,k,ip,im,jp,jm,kp,km;
    double h2,ei,wk;

    h2 = 2.0*dx*dx;
    ei = energyOf(psi,vv);
    for (i=1; i<NNx-1; i++) {
	/*
      ip = i+1; if (ip>=NNx) ip = 0;
      im = i-1; if (im<0) im = NNx-1;
	*/
      for (j=0; j<NNy; j++) {
        jp = j+1; if (jp>=NNy) jp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;
        for (k=0; k<NNz; k++) {
          kp = k+1; if (kp>=NNz) kp = 0;
          km = k-1; if (km<0) km = NNz-1;
          psi[i][j][k] -= dump*((6.0*psi[i][j][k]-psi[i+1][j][k]-psi[i-1][j][k]-psi[i][jp][k]-psi[i][jm][k]-psi[i][j][kp]-psi[i][j][km])/h2 + (v[i][j][k]-ei)*psi[i][j][k]);
        }
      }
    }
    normalize(psi);
    return (ei);
  }

  double energyOf(double psi[][][], double v[][][]) {
    int i,j,k,ip,im,jp,jm,kp,km;
    double h2,s,sn,ps;

    h2 = 2.0*dx*dx;
    s = 0.0; sn=0.0;
    for (i=1; i<NNx-1; i++) {
	/*
      ip = i+1; if (ip>=NNx) ip = 0;
      im = i-1; if (im<0) im = NNx-1;
	  */
      for (j=0; j<NNy; j++) {
        jp = j+1; if (jp>=NNy) jp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;
        for (k=0; k<NNz; k++) {
          kp = k+1; if (kp>=NNz) kp = 0;
          km = k-1; if (km<0) km = NNz-1;
		  ps = psi[i][j][k];
          s += ps*((6.0*ps-psi[i+1][j][k]-psi[i-1][j][k]-psi[i][jp][k]-psi[i][jm][k]-psi[i][j][kp]-psi[i][j][km])/h2 + v[i][j][k]*ps);
          sn += ps*ps;
        }
      }
    }
    return ( s/sn );
  }

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

}