/** applet No. 1164
 *
 * jellium Slab
 * - y,z-periodic - external charge density - real space DFT-LDA
 * jellium(Z=36) - electron 36 system
 * LDA:  J. P. Perdew and A. Zunger; Phys. Rev., B23, 5048 (1981)
 * - multi thread  ECD :: display - asynchronous
 *
 * Created by Ikeuchi Mitsuru on April 28 2007.
 * Copyright (c) 2007 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2007.04.28  created
 * ver 0.0.2   2007.06.23  improved code
 *
 */

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

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

  // ------------------------------------ preset field -----------

  Thread th = null; // for run()-paint() thread
  RSDFTPyzECD rs = null;
  DrawGraph3d dg3d = new DrawGraph3d();

  // for event
  Choice ch_noe, ch_view;
  Button bt_reset, bt_startStop, bt_step;
  Scrollbar sc_mixing, sc_broadening, sc_orbit, sc_kPos;

  // for off-paint buffer
  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int dgX, dgY, dgXb,dgYb; // mouse

  int sleepTime = 50;
  int dispMode = 9;
  int thCount = 0;


  // ------------------------------ applet main thread -----------

  public void init() {

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

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

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

    ch_view = new Choice();
    ch_view.add("orbit dens");ch_view.add("orbit phase");
    ch_view.add("ph-grid");ch_view.add("Vext/4");
    ch_view.add("VH/4"); ch_view.add("Vx");
    ch_view.add("Vc"); ch_view.add("Veff");
    ch_view.add("rho");ch_view.add("rho 3D");
    ch_view.add("rho /zPos");
    ch_view.add("rho+charge");
    ch_view.addItemListener(this);
    ch_view.select("rho 3D");

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

    bt_startStop = new Button("start/stop");
    bt_startStop.addActionListener(this);

    bt_step = new Button("step");
    bt_step.addActionListener(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,17,10,0,rs.Nstate+10-1);
    sc_orbit.addAdjustmentListener(this);

    sc_kPos = new Scrollbar(Scrollbar.HORIZONTAL,rs.NNz/2,10,0,rs.NNz+10-1);
    sc_kPos.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(2,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(bt_step);
    pnl.add(new Label("  "));
    pnl.add(new Label("  view"));
    pnl.add(ch_view);

    pnl.add(new Label("  "));
    pnl.add(ch_noe);
    pnl.add(sc_mixing);
    pnl.add(sc_broadening);
    pnl.add(sc_orbit);
    pnl.add(sc_kPos);
    add(pnl,"North");

  }

  public void start() {
    if (rs == null) {
      rs = new RSDFTPyzECD();
      rs.start();
    }
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

  public void stop() {
    if (th != null) th = null;
    if (rs != null) rs = null;
  }


  // ---------------------------------- event listener -----------

  public void itemStateChanged(ItemEvent ev) {
    double zz,ne;

    if (ev.getSource() == ch_noe){
	  ne = rs.getNuclearCarge() - 1.0 + 0.5*ch_noe.getSelectedIndex();
	  rs.setNumberOfElectron(ne);
    } else if (ev.getSource() == ch_view){
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      rs.reset();
    } else if (ev.getSource() == bt_startStop){
      if (rs.getStartSW()==0) {
        rs.setStartSW(1);
      } else {
        rs.setStartSW(0);
      }
    } else if (ev.getSource() == bt_step){
      if (rs.getStartSW()==0) {
        rs.setStepSW();
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_broadening) {
	  rs.setBroadening( 0.001*sc_broadening.getValue() );
    } else if (ev.getSource() == sc_mixing) { 
      rs.setMixing( 0.01*(double)(sc_mixing.getValue()) );
    } else if (ev.getSource() == sc_orbit) { 
      dg3d.setDispObj( sc_orbit.getValue() );
    } else if (ev.getSource() == sc_kPos) { 
      dg3d.setkPos( sc_kPos.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;
    }
    dg3d.setViewDelta( (dgX-dgXb), (dgY-dgYb) );
  }



  // ========================= run() - paint() loop ==============

  public void run() {
    while (th != null) {
      thCount += 1;
      offPaint();
      repaint();
      try {
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }

  public void offPaint() {
    if (rs != null) {
      dg3d.plotGraph(gOff, rs, dispMode, thCount);
    }
  }

}



// =============================== draw graph 3d class ===========

class DrawGraph3d {

  private Graphics gOff;
  private RSDFTPyzECD rs;

  private int NNx = rs.NNx;
  private int NNy = rs.NNy;
  private int NNz = rs.NNz;
  private double dx = rs.dx;
  private double dy = rs.dy;
  private double dz = rs.dz;
  private int Nstate = rs.Nstate;

  private double xMax = NNx*dx;
  private double yMax = NNy*dy;
  private double zMax = NNz*dz;
  private double cxx = xMax/2.0;
  private double cyy = yMax/2.0;
  private double czz = zMax/2.0;

  private double dispScale = 300.0/xMax;
  private double viewScale = 1.0;
  private double scale = dispScale*viewScale;

  private int xImageLoc = 50;
  private int yImageLoc = 120;

  private int kPos = NNz/2;  
  private int dispObj = 17;

  private double viewTheta = -15.0*3.14/180.0;
  private double viewFai = -72.0*3.14/180.0;
  private double cosTh = Math.cos(viewTheta); 
  private double sinTh = Math.sin(viewTheta);
  private double cosFi = Math.cos(viewFai);
  private double sinFi = Math.sin(viewFai);

  private double wkx[] = new double[8]; 
  private double wky[] = new double[8]; 
  private double wkz[] = new double[8]; 

  private double pwkx[] = new double[8]; 
  private double pwky[] = new double[8]; 
  private double pwkz[] = new double[8]; 

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

  private int xpts[] = new int[500];
  private int ypts[] = new int[500];

  private int drawnCount = 0;


  // ------------------------------- class constructor -----------

  DrawGraph3d() {
  }


  // ----------------------------------- class methods -----------

  public void plotGraph(Graphics gg, RSDFTPyzECD rsdft, int dispMode, int thCount) {
    int rsCount,gb,i;

    gOff = gg;
    rs = rsdft;
    if (rs==null) return;
    rsCount = rs.getCount();
    //if (rsCount==drawnCount) return;
    //drawnCount = rsCount;

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,630,460);

    if (dispMode==0) {
      densityInterPlot( 2, dispObj, 0);

    } else if (dispMode==1) {
      densityInterPlot( 3, dispObj, 0);

    } else if (dispMode==2) {
      gridFnPlot(rs.getOrbit(dispObj), 20.0);

    } else if (dispMode==3) {
      gridFnPlot(rs.getVext(), 0.25);

    } else if (dispMode==4) {
      gridFnPlot(rs.getVh(), 0.25);

    } else if (dispMode==5) {
      gridFnPlot(rs.getVx(), 1.0);

    } else if (dispMode==6) {
      gridFnPlot(rs.getVc(), 1.0);

    } else if (dispMode==7) {
      gridFnPlot(rs.getV(), 1.0);

    } else if (dispMode==8) {
      gridFnPlot(rs.getDensity(), 1.0);

    } else if (dispMode==9) {
      densityInterPlot(0, dispObj, 0);

    } else if (dispMode==10) {
      densityInterPlot(1, dispObj, 0);

    } else if (dispMode==11) {
      densityInterPlot(0, dispObj ,1);

    }

    gOff.setColor(Color.black);
    gOff.drawString("Z="+rs.getNuclearCarge()+" ",630/6*0+10,70);
    gOff.drawString("ne="+rs.getNumberOfElectron()+" ",630/6*1+20,70);
    gOff.drawString("mix="+(int)(rs.getMixing()*100+0.5)/100.0+" ",630/6*2+10,70);
    gOff.drawString("kT0="+(int)(rs.getkT0()*1000+0.5)/1000.0+" ",630/6*3+10,70);
    gOff.drawString("orbit="+dispObj+" ",630/6*4+10,70);
    gOff.drawString("zPos="+kPos+" ",630/6*5+20,70);
    gOff.drawString("W slab",10,90);

    gb = 400;
    gOff.setColor(Color.black);
    gOff.drawString("kT="+(int)(rs.getkT()*1000.0+0.5)/1000.0+" ",gb,90);
    gOff.drawString("iteration="+rs.getCount()+" ",gb+120,90);
    gOff.drawString("box="+(int)(rs.NNx*rs.dx)+"x"+ +(int)(rs.NNy*rs.dy)+"x"+(int)(rs.NNz*rs.dz)+" au",gb,110);
    gOff.drawString("dx="+(int)(dx*100.0)/100.0+" au",gb+120,110);
    gOff.drawString("total energy="+(float)rs.getTotalEnergy()+" au",gb,130);

    dispTable(gb, 150);
  }

  private void dispTable(int gx, int gy) {
    int i,i0,gyi;

    gOff.setColor(Color.black);
    gOff.drawString("No.",gx,gy);
    gOff.drawString("occupation",gx+30,gy);
    gOff.drawString("energy",gx+150,gy);

    i0 = dispObj-10; if (i0<0) i0 = 0;
    if (i0+20>Nstate) i0 = Nstate - 20;

    for (i=i0; i<i0+20; i++) {
      gyi = gy + 15 + 12*(i-i0);
      gOff.setColor(Color.getHSBColor((float)(0.66-0.33*rs.getOccupation(i)), 0.9f,0.8f));
      gOff.drawString(""+i+"",gx,gyi);
      gOff.drawString(" "+(float)(rs.getOccupation(i))+"",gx+30,gyi);
      gOff.drawString(" "+(float)(rs.getLevel(i))+"",gx+120,gyi);
    }
  }


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

  //
  // densityPlot(int mode, int orbit)
  //
  // mode: 0-density 1-half  2-orbit density 3-orbit phase
  //

  private void densityPlot(int mode, int orbit) {
    int i,j,k,ix,iy,ir,fp;
    int iBegin, iEnd, iInc, jBegin, jEnd, jInc, kBegin, kEnd, kInc;
    float orbden;
    double x,y,z,ppx,ppy,ppz,dens;
    double ptr[][][];

    ptr = rs.getOrbit(orbit);
    if (mode<=1) ptr = rs.getDensity();
    orbden = (float)(0.3+0.3*rs.getOccupation(orbit));

    fp = drawBox(0);
    iBegin = 0; iEnd = NNx; iInc = 1;
    if (wkx[fp]>10.0) {
      iBegin = NNx-1; iEnd = 0; iInc = -1;
    }
    jBegin = 0; jEnd = NNy; jInc = 1;
    if (wky[fp]>10.0) {
      jBegin = NNy-1; jEnd = 0; jInc = -1;
    }
    kBegin = 0; kEnd = NNz; kInc = 1;
    if (wkz[fp]>10.0) {
      kBegin = NNz-1; kEnd = 0; kInc = -1;
    }

    for (i=iBegin; ((i<iEnd && iInc>0) || (i>=iEnd && iInc<0)) ; i += iInc) {
      x = dx*i;
      for (j=jBegin; ((j<jEnd && jInc>0) || (j>=jEnd && jInc<0)) ; j += jInc) {
        y = dy*j;
        for (k=kBegin; ((k<kEnd && kInc>0) || (k>=kEnd && kInc<0)) ; k += kInc) {
          z = dz*k;

          dens = 0.0;
          if (mode<=1) {
            dens = 1.0*ptr[i][j][k]; if (dens>0.999) dens = 0.999;
            gOff.setColor(Color.getHSBColor((float)(0.66-0.66*dens), 0.9f,0.9f));
            if (mode==1 && k>kPos) dens = 0.0;
          } else if (mode==2) {
            dens = 5.0*ptr[i][j][k]*ptr[i][j][k]; if (dens>0.999) dens = 0.999;
            gOff.setColor(Color.getHSBColor((float)(0.66-0.66*dens), orbden,0.9f));
          } else if (mode==3) {
            dens = 5.0*ptr[i][j][k]*ptr[i][j][k]; if (dens>0.999) dens = 0.999;
            if (ptr[i][j][k]>0.0) {
              gOff.setColor(Color.getHSBColor(0.33f, orbden,0.9f));
            } else {
              gOff.setColor(Color.getHSBColor(0.01f, orbden,0.9f));
            }
          }
          ir = (int)(10.0*dens+0.98);

          if (ir>0) {
            ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
            ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
            ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
            ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

            ix = (int)(scale*ppx)+xImageLoc;
            iy = (int)(scale*ppy)+yImageLoc;
			
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
        }
      }
    }

    drawBox(1);
  }

  private void densityInterPlot(int mode, int orbit, int nuc) {
    int i,j,k,ix,iy,ir,fp, ii,im, jj,jm, kk,km;
    int iBegin, iEnd, iInc, jBegin, jEnd, jInc, kBegin, kEnd, kInc;
    float orbden;
    double x,y,z,ppx,ppy,ppz,dens,d;
    double[][][] ptr, rhoExt;

    ptr = rs.getOrbit(orbit);
	rhoExt = rs.getrhoExt();
    if (mode<=1) ptr = rs.getDensity();
    orbden = (float)(0.3+0.3*rs.getOccupation(orbit));

    fp = drawBox(0);
    iBegin = 0; iEnd = 2*NNx-1; iInc = 1;
    if (wkx[fp]>10.0) {
      iBegin = 2*NNx-2; iEnd = 0; iInc = -1;
    }
    jBegin = 0; jEnd = 2*NNy-1; jInc = 1;
    if (wky[fp]>10.0) {
      jBegin = 2*NNy-2; jEnd = 0; jInc = -1;
    }
    kBegin = 0; kEnd = 2*NNz-1; kInc = 1;
    if (wkz[fp]>10.0) {
      kBegin = 2*NNz-2; kEnd = 0; kInc = -1;
    }

    for (i=iBegin; ((i<iEnd && iInc>0) || (i>=iEnd && iInc<0)) ; i += iInc) {
      x = 0.5*dx*i; ii = i/2; im = i/2+i%2;
      for (j=jBegin; ((j<jEnd && jInc>0) || (j>=jEnd && jInc<0)) ; j += jInc) {
        y = 0.5*dy*j; jj = j/2; jm = j/2+j%2;
        for (k=kBegin; ((k<kEnd && kInc>0) || (k>=kEnd && kInc<0)) ; k += kInc) {
          z = 0.5*dz*k; kk = k/2; km = k/2+k%2;
          d = ptr[ii][jj][kk]+ptr[im][jj][kk]+ptr[ii][jm][kk]+ptr[ii][jj][km]
		     +ptr[im][jm][kk]+ptr[im][jj][km]+ptr[ii][jm][km]+ptr[im][jm][km];
          dens = 0.0;
          if (mode<=1) {
            dens = 0.125*d; if (dens>0.999) dens = 0.999;
            gOff.setColor(Color.getHSBColor((float)(0.66-0.66*dens), 0.9f,0.9f));
            if (mode==1 && k>2*kPos) dens = 0.0;
          } else if (mode==2) {
            dens = 5.0*0.125*d*d; if (dens>0.999) dens = 0.999;
            gOff.setColor(Color.getHSBColor((float)(0.66-0.66*dens), orbden,0.9f));
          } else if (mode==3) {
            dens = 5.0*0.125*d*d; if (dens>0.999) dens = 0.999;
            if (d>0.0) {
              gOff.setColor(Color.getHSBColor(0.33f, orbden,0.9f));
            } else {
              gOff.setColor(Color.getHSBColor(0.01f, orbden,0.9f));
            }
          }
          ir = (int)(10.0*dens+0.98);

          if (ir>0) {
            ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
            ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
            ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
            ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

            ix = (int)(scale*ppx)+xImageLoc;
            iy = (int)(scale*ppy)+yImageLoc;
			
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
          if (nuc==1) {
	  if (i%2==0 && j%2==0 && k%2==0 && (rhoExt[ii][jj][kk]<0.0)) {
            ppy = cosFi*(y-cyy)+sinFi*(z-czz) + cyy; 
            ppz = -sinFi*(y-cyy)+cosFi*(z-czz) + czz;
            ppx = cosTh*(x-cxx)+sinTh*(ppz-czz) + cxx; 
            ppz = -sinTh*(x-cxx)+cosTh*(ppz-czz) + czz;

            ix = (int)(scale*ppx)+xImageLoc;
            iy = (int)(scale*ppy)+yImageLoc;
		    ir = 3;
            gOff.setColor(Color.getHSBColor(0.85f, 0.9f,0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
	  }
        }
      }
    }

    drawBox(1);
  }

  // draw box

  private int drawBox(int imode) {
    int i,gx,gy,g2x,g2y,farPoint;
    double sc;

    if (imode==0) {
      setBox();
      rotateBox();
    }
    sc = dispScale*viewScale;
    farPoint = findFarPoint();

    for (i=0; i<12; i++) {
      gx = (int)(sc*pwkx[boxp[i][0]])+xImageLoc;
      gy = (int)(sc*pwky[boxp[i][0]])+yImageLoc;
      g2x = (int)(sc*pwkx[boxp[i][1]])+xImageLoc;
      g2y = (int)(sc*pwky[boxp[i][1]])+yImageLoc;
      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.darkGray);
          gOff.drawLine(gx,gy, g2x, g2y);
        }
      }
    }
    return farPoint;
  }

  private int findFarPoint() {
    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);
  }

  private void setBox() {
    int i,ix,iy,iz;
    double ballSize;

    i=0;
    ballSize = 0.0*xMax;

    for (ix=0; ix<2; ix++) {
      for (iy=0; iy<2; iy++) {
        for (iz=0; iz<2; iz++) {
          wkx[i] = -0.5*ballSize+(xMax+ballSize)*(double)(ix);
          wky[i] = -0.5*ballSize+(yMax+ballSize)*(double)(iy);
          wkz[i] = -0.5*ballSize+(zMax+ballSize)*(double)(iz);
          i += 1;
        }
      }
    }
  }

  private void rotateBox() {
    int i;
 
    for (i=0; i<8; i++) {
      pwky[i] = cosFi*(wky[i]-cyy)+sinFi*(wkz[i]-czz) + cyy; 
      pwkz[i] = -sinFi*(wky[i]-cyy)+cosFi*(wkz[i]-czz) + czz;
      pwkx[i] = cosTh*(wkx[i]-cxx)+sinTh*(pwkz[i]-czz) + cxx; 
      pwkz[i] = -sinTh*(wkx[i]-cxx)+cosTh*(pwkz[i]-czz) + czz;
    }  
  }

  //

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

    nx = fn.length; ny = fn[0].length; nz = fn[0][0].length;

    cx = nx/2;
    cy = ny/2;
    cz = 0;
    sc = 300.0/NNx;

    gbx = 50+(int)(cx/sc);
    gby = 120+(int)(cy/sc);
    k = kPos;

    for(j=0; j<ny; j += 1) {
      for (i=0; i<nx; i++) { 
        z = weight*fn[i][j][k] + (kPos-nz/2) + 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<nx-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<nx; i += 1) {
      for (j=0; j<ny; j++) {
        z = weight*fn[i][j][k] + (kPos-nz/2) + 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<ny-1; j++) { 
        gOff.setColor(Color.getHSBColor(0.66f,0.6f,0.9f));
        gOff.drawLine(xpts[j],ypts[j],xpts[j+1],ypts[j+1]);
      }
    }
  }


  // ------------------------------------- I/O methods -----------

  public void setViewDelta(double xDisplacement, double yDisplacement) {
    viewTheta += 0.5*3.14159/180.0*xDisplacement;
    viewFai += 0.5*3.14159/180.0*yDisplacement;
    cosTh = Math.cos(viewTheta); 
    sinTh = Math.sin(viewTheta);
    cosFi = Math.cos(viewFai);
    sinFi = Math.sin(viewFai);
  }

  public void setDispObj(int obj) { dispObj = obj; }
  public int getDispObj() { return dispObj; }
  public void setkPos(int kp) { kPos = kp; }

}



// ======================== RSDFT - periodic ECD class ===========

class RSDFTPyzECD extends Thread {

  public static double dx = 1.0;
  public static double dy = 1.0;
  public static double dz = 1.0;

  public static int NNx = 28;
  public static int NNy = 12;
  public static int NNz = 12;

  public static int Nstate = 32;

  private double SDdump = 0.2;
  private double mixing = 0.2;
  private double broadening = 0.001;

  private double NOE = 36.0;
  private double Znuc = 36.0;
  private int setVhQ = 0;


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

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

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

  private double sdEnergy[] = new double[Nstate];
  private double sdState[][][][] = new double[Nstate][NNx][NNy][NNz];
  private double occ[] = new double[Nstate];

  private double nuc[][] = {
    { 6.0,  6.0,  3.0,  3.0  },
    { 6.0,  6.0,  3.0,  9.0  },
    { 6.0,  6.0,  9.0,  3.0  },
    { 6.0,  6.0,  9.0,  9.0  },

    { 6.0,  9.0,  0.0,  0.0  },
    { 6.0,  9.0,  0.0,  6.0  },
    { 6.0,  9.0,  6.0,  0.0  },
    { 6.0,  9.0,  6.0,  6.0  },

    { 6.0, 12.0,  3.0,  3.0  },
    { 6.0, 12.0,  3.0,  9.0  },
    { 6.0, 12.0,  9.0,  3.0  },
    { 6.0, 12.0,  9.0,  9.0  },

    { 0.0,  0.0 , 0.0 , 0.0  }
  };

  private int resetSW = 0;
  private int startSW = 1;
  private int stepSW = 0;
  private int sleepTime = 1;
  private int stopSleepTime = 100;
  private int iterCount = 0;


  // ------------------------------- class constructor -----------

  RSDFTPyzECD() {
    setInitialCondition();
  }


  // -------------------------------------- thread run -----------

  public void run() {
    while (true) {
      iteration();
      if (sleepTime>0) {
        try {
          Thread.sleep(sleepTime);
        } catch (InterruptedException e) { }
      }
    }
  }


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

  private void setInitialCondition() {

    iterCount = 0;
    setPotential();
    setOcc(Nstate,NOE);
    setState(Nstate);
  }

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

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

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

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

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

  private void setJelliumSlab(double x0, double x1) {
    int i,j,k;
    double s,x;

    s = 0.0;
    for (i=0;i<NNx;i++) {
      x = i*dx;
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          if (x>=x0 && x<=x1) {
            rhoExt[i][j][k] = 1.0;
            s += 1.0;
          }
        }
      }
    }

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          rhoExt[i][j][k] *= -Znuc/s/(dx*dy*dz);
        }
      }
    }
  }

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

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


  // --------------------------------------- iteration -----------

  private void iteration() {

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

    if (startSW==1) {
      iterStep();
    } else {
      if (stepSW==1) {
        iterStep();
        stepSW = 0;
      }
      try {
        Thread.sleep(stopSleepTime);
      } catch (InterruptedException e) { }
    }
  }

  private void iterStep() {
    iterCount += 1;
    setVeff();
    SD(Nstate, vv, 2, SDdump);
    sortState(Nstate);
    setOcc(Nstate,NOE);
  }

  // (1) set Veff[][][]

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

  private 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<Nstate; ie++) {
            if (occ[ie]<=0.0) break;
            rho[i][j][k] += mx*occ[ie]*(sdState[ie][i][j][k]*sdState[ie][i][j][k]);
          }
        }
      }
    }
  }

  private 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)

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

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

  // (2) steepest descent - sdEnergy[ist], sdState[ist][][][]

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

  // 2-1 steepest descent

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

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

  // 2-2 Gram-Schmidt

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

  private 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];
        }
      }
    }
    return ( s*dv );
  }

  private 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];
        }
      }
    }
    s *= 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];
        }
      }
    }
  }

  // (3) sort state

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

  // (4) set occ[ist]

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

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

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

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


  // ----------------------------------------- utility -----------

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

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

    enn = 0.0;
	/*
    for (i=0;i<100;i++) {
      if (nuc[i][0]<=0.0) break;
      for (j=i+1;j<100;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 );
  }


  // ------------------------------------- I/O methods -----------

  // control
  public void reset() { resetSW = 1; }
  public void setStartSW(int sw) { startSW = sw; }
  public int getStartSW() { return startSW; }
  public void setStepSW() { stepSW = 1; }
  public int getCount() { return iterCount; }

  public void setNumberOfElectron(double ne) {
    NOE = ne;
    iterCount = 0;
  }
  public double getNumberOfElectron() { return NOE; }
  public double getNuclearCarge() { return Znuc; }

  public void setBroadening(double br) { broadening = br; }
  public double getkT0() { return broadening; }
  public double getkT() { return levelWidth(); }
  public void setMixing(double mx) { mixing = mx; }
  public double getMixing() { return mixing; }

  // RSDFT get data
  public double getTotalEnergy() { return totalEnergy(); }
  public double getLevel(int i) { return sdEnergy[i]; }
  public double getOccupation(int i) { return occ[i]; }
  public double[][][] getOrbit(int i) { return sdState[i]; }

  public double[][][] getDensity() { return rho; }
  public double[][][] getrhoExt() { return rhoExt; }
  public double[][][] getV() { return vv; }
  public double[][][] getVext() { return vvext; }
  public double[][][] getVh() { return vvh; }
  public double[][][] getVx() { return vvx; }
  public double[][][] getVc() { return vvc; }
}