/* Fermi sea DFT LDA 1D        */
/* coded by Ikeuchi Mitsuru    */
/* ver 0.0.1   2004.07.04      */
/* ver 0.0.2   2004.07.09  rename bandDFT1D to FermiSeaDFT1D */
/* ver 0.0.3   2004.09.21  setVcx() bug fixed */

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

public class FermiSeaDFT1D extends Applet
  implements ItemListener, ActionListener, AdjustmentListener, Runnable {

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

  Choice cpt,cne,cvw;
  Button bt_reset, bt_start, bt_stop;
  Scrollbar sc_ds, sc_ww, sc_wh, spx;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int started = 1;
  double px0 = 12.0;
  double vDist = 2.5;
  double vWidth = 2.22;
  double vHight = -1.5;
  double energy = 0.5;
  int potentialMode = 5;
  int dispMode = 0;
  int numberOfElectrons = 4;

  double t = 0.0;
  double dt = 1.0/256.0;
  double dx = 1.0/16.0;
  double lylz = 4.0*4.0;

  int NNx = 256;
  double ph[] = new double[NNx];
  double vv[] = new double[NNx];

  double vvext[] = new double[NNx];
  double vvh[] = new double[NNx];
  double vvx[] = new double[NNx];
  double vvc[] = new double[NNx];
  
  double rho[] = new double[NNx];

  double eigenEnergy[] = new double[20];
  double eigenState[][] = new double[20][NNx];

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

  public void init() {

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

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


    cpt = new Choice();
    cpt.addItem("well"); cpt.addItem("parabolic"); cpt.addItem("triangle");
    cpt.addItem("band");
    cpt.addItemListener(this);
    cpt.select("band");

    cne = new Choice();
    cne.addItem("2"); cne.addItem("4"); cne.addItem("6"); cne.addItem("8");
    cne.addItemListener(this);
    cne.select("4");

    cvw = new Choice();
    cvw.addItem("wave func"); cvw.addItem("density"); cvw.addItem("both");
    cvw.addItemListener(this); 

    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_ds= new Scrollbar(Scrollbar.HORIZONTAL,250,10,50,1210);
    sc_ds.addAdjustmentListener(this);

    sc_ww= new Scrollbar(Scrollbar.HORIZONTAL,222,10,100,260);
    sc_ww.addAdjustmentListener(this);

    sc_wh= new Scrollbar(Scrollbar.HORIZONTAL,150,10,0,310);
    sc_wh.addAdjustmentListener(this);

    spx= new Scrollbar(Scrollbar.HORIZONTAL,500,10,0,2010);
    spx.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(new Label("  "));
    pnl.add(sc_ww);
    pnl.add(sc_wh);
    pnl.add(cne);
    pnl.add(new Label("  "));
    pnl.add(cvw);
    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() == cpt){
      potentialMode = 1+cpt.getSelectedIndex();
      setInitialCondition();
    } else if (ev.getSource() == cne){
      numberOfElectrons = 2*(cne.getSelectedIndex()+1);
    } else if (ev.getSource() == cvw){
      dispMode = cvw.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() == spx) { 
      energy = 0.001*(double)(spx.getValue())-1.0;
    } else if (ev.getSource() == sc_ds) { 
      vDist = 0.01*(double)(sc_ds.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_ww) { 
      vWidth = 0.01*(double)(sc_ww.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_wh) { 
      vHight = -0.01*(double)(sc_wh.getValue());
      setInitialCondition();
    }
  }

  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,istate,gbx,gby,gby2,gb;
    double emag, dmag;

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

    gbx = 30;
    gby = 160;
    emag = 100.0;
    dmag = 50.0;

    gOff.setColor(Color.getHSBColor( 0.85f, 0.4f, 0.9f ));
    fillGraph(rho, gbx, gby-10, lylz*dmag);

    for (istate=1;istate<7; istate++) {
      gby2 = gby-(int)(emag*eigenEnergy[istate]);
      setClr(istate, 0);
      gOff.drawString("E"+istate+" = "+(float)eigenEnergy[istate]+" au",400,160-15*istate);

      if ((dispMode==1 || dispMode==2) && istate*2<=numberOfElectrons) {
        setClr(istate, 1);
        drawDensity(istate, gbx, gby2, dmag);
      }

      if (dispMode==0 || dispMode==2) {
        setClr(istate, 0);
        drawGraph(eigenState[istate], gbx, gby2, dmag);
      }
    }

    gb = 300;
    gOff.setColor(Color.getHSBColor( 0.6f, 0.6f, 0.7f ));
    drawGraph(vvext, gbx, gby, emag);
    gOff.drawString("Vext",gb,100);

    gOff.setColor(Color.getHSBColor( 0.7f, 0.6f, 0.7f ));
    drawGraph(vvh, gbx, gby, emag);
    gOff.drawString("VH",gb,120);

    gOff.setColor(Color.getHSBColor( 0.8f, 0.6f, 0.7f ));
    drawGraph(vvx, gbx, gby, emag);
    gOff.drawString("Vx",gb,140);

    gOff.setColor(Color.getHSBColor( 0.01f, 0.6f, 0.7f ));
    drawGraph(vvc, gbx, gby, emag);
    gOff.drawString("Vc",gb,160);

    gOff.setColor(Color.getHSBColor( 0.33f, 0.9f, 0.9f ));
    drawGraph(vv, gbx, gby, emag);
    gOff.drawString("Veff",gb,180);

    gOff.setColor(Color.black);
    gOff.drawString("dist.="+(int)(vDist*1000.0)/1000.0+" au",630/6*0+10,40);
    gOff.drawString("width="+(int)(vWidth*1000.0)/1000.0+" au",630/6*1+10,40);
    gOff.drawString("depth="+(int)(vHight*1000.0)/1000.0+" au",630/6*2+10,40);
    gOff.drawString("N of electrons",630/6*3+10,40);
    gOff.drawString("view",630/6*5+10,40);

    for (istate=1;istate<=3; istate++) {
      for (i=1;i<=3; i++) {
        gOff.drawString("<"+istate+"|"+i+">="+(float)intersection(istate,i)+" ",400,180+(istate*3+i-4)*15);
      }
    }

  }

  void setClr(int istate, int mode) {
    if (mode==0) {
      if (istate*2<=numberOfElectrons) {
        gOff.setColor(Color.getHSBColor( 0.01f, 0.9f, (float)(0.4+0.1*istate) ));
      } else {
        gOff.setColor(Color.lightGray);
      }
    } else if (mode==1) {
      if (istate*2<=numberOfElectrons) {
        gOff.setColor(Color.getHSBColor( 0.6f, 0.6f, 0.9f));
      } else {
        gOff.setColor(Color.getHSBColor( 0.6f, 0.1f, 0.9f));
      }
    }
  }

  void drawDensity(int istate, int gx, int gy, double mag) {
    int i;
    double ph2;

    for (i=1; i<NNx-2; i++) {
      ph2 = mag*eigenState[istate][i]*eigenState[istate][i];
      if (ph2>=1.0) {
        gOff.drawLine(gx+i,gy-(int)(ph2), gx+i, gy);
      }
    }
  }

  void drawGraph(double fn[],int gx, int gy, double mag) {
    int i;
    double f1,f2;

    for (i=1; i<NNx-2; i++) {
      f1 = mag*fn[i]; f2 = mag*fn[i+1];
	  gOff.drawLine(gx+i,gy-(int)(f1), gx+i, gy-(int)(f2));
    }
  }

  void fillGraph(double fn[],int gx, int gy, double mag) {
    int i;
    double f;

    for (i=1; i<NNx-2; i++) {
      f = mag*fn[i];
      if (f>=1.0) {
        gOff.drawLine(gx+i,gy-(int)(f), gx+i, gy);
      }
    }
  }


/*------------------------- <i|j> integral ----------*/

  double intersection(int st1,int st2) {
    int i;
    double s;

    s = 0.0;
    for (i=0; i<NNx; i++) {
      s += eigenState[st1][i]*eigenState[st2][i]*dx;
    }
    return( s );
  }
  
/*-------------------------- set initial condition ----------*/

  void setInitialCondition() {
    int i;

    for (i=0; i<NNx; i++) {
      vv[i] = 0.0;
      vvext[i] = 0.0;
      vvh[i] = 0.0;
      vvx[i] = 0.0;
      vvc[i] = 0.0;
    }

    setVext(potentialMode, vWidth, vHight);
  }

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

  void setVext(int mode, double vw, double vh) {
    int i,j;
    double x,x0,a;

    x0 = 8.0;
    if (mode==1) {

      for (i=1;i<NNx-1;i++) {
        x = i*dx;
        if (x>x0-0.5*vw && x<(x0+0.5*vw)) {
          vvext[i] = vh;
        } else {
          vvext[i] = 0.0;
        }
      }

    } else if (mode==2) {

      a = -vh*4.0/(vw*vw);
      for (i=1;i<NNx-1;i++) {
        x = i*dx;
        if (x>x0-0.5*vw && x<(x0+0.5*vw)) {
          vvext[i] = a*(x-x0)*(x-x0)+vh;
        } else {
          vvext[i] = 0.0;
        }
      }

    } else if (mode==3) {

      for (i=1;i<NNx-1;i++) {
        x = i*dx;
        if (x>x0-0.5*vw && x<=x0) {
          vvext[i] = 2.0*((0.5*vw-(x0-x))/vw)*vh;
        } else if (x>x0 && x<=x0+0.5*vw) {
          vvext[i] = 2.0*((0.5*vw-(x-x0))/vw)*vh;
        } else {
          vvext[i] = 0.0;
        }
      }

    } else if (mode==4) {

      for (i=1;i<NNx-1;i++) {
        x = i*dx;
        if (x>x0-0.5*vDist-0.5*vw && x<=x0-0.5*vDist+0.5*vw) {
          vvext[i] = vh;
        } else if (x>x0+0.5*vDist-0.5*vw && x<=x0+0.5*vDist+0.5*vw) {
          vvext[i] = vh;
        } else {
          vvext[i] = 0.0;
        }
      }

    } else if (mode==5) {

      for (i=1;i<NNx-1;i++) {
        x = i*dx;
        vvext[i] = 0.0;
        for (j=0;j<4;j++) {
          if (x>4.0+j*vDist-0.5*vw && x<=4.0+j*vDist+0.5*vw) {
            vvext[i] = vh;
          }
        }
      }

    }

  }

/* ----------------------------- time evolution -------------- */

  void timeEvolution() {
    int i,istate,nstate;

    if (started==1) {
      for (i=1;i<NNx-1;i++) {
        vv[i] = vvext[i]+vvh[i]+vvx[i]+vvc[i];
      }

      getEigenValue(vHight);
      setRho(numberOfElectrons/2);
      poisson();
      setVcx();

    } else if (started==-1) {
      t = 0.0;
      setInitialCondition();
      started = 0;
    }
  }

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

  void setRho(int istateMax) {
    int i,istate;
    double g;

    g = 2.0;

    rho[0] = 0.0; rho[NNx-1] = 0.0;
    for (i=1; i<NNx-1; i++) {
      rho[i] = 0.0;
      for (istate=1; istate<=istateMax; istate++) {
        rho[i] += g*eigenState[istate][i]*eigenState[istate][i]/lylz;
      }
    }
  }

  void poisson() {
    int i,iter;

    for (iter=1; iter<1000; iter++) {
      for (i=1; i<NNx-1; i++) {
        vvh[i] = vvh[i]+0.66*(vvh[i+1]+vvh[i-1]-2.0*vvh[i]+dx*dx*rho[i]);
      }
    }
  }

  void setVcx() {
    int i;
    double drho;

    drho = 0.01;
    for (i=1; i<NNx-1; i++) {
      vvx[i] = eeExchange(rho[i])+rho[i]*(eeExchange(rho[i]+drho)-eeExchange(rho[i]))/drho;
      vvc[i] = eeCorrelation(rho[i])+rho[i]*(eeCorrelation(rho[i]+drho)-eeCorrelation(rho[i]))/drho;
    }
  }

/* exchange term :  ex(rho) = -(3/4)*(3/pai)^(1/3)*rho^(1/3)  */

  double eeExchange(double rh) {
    return( -0.7368*Math.pow(rh,0.33333333) );
  }

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

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

/* ----------------------- get eigen value and eigen state -------- */

  void getEigenValue(double energyLow) {
    int i;
    double err;

    err = 1.0e-14;
    eigenEnergy[0] = energyLow-0.001;
    for (i=1; i<10; i++) {
      eigenEnergy[i] = seekEigenValue(eigenEnergy[i-1]+0.001,err);
      setEigenState(i,eigenEnergy[i]);
    }
  }

  void setEigenState(int istate,double energy) {
    int i;

    trialPh(energy);
    normalize();
    for (i=1; i<NNx-1; i++) {
      eigenState[istate][i] = ph[i];
    }
  }

  double seekEigenValue(double energyLow,double error) {
    double eh,el,em, de;
    double phLow,phHigh,phMid;

    de = 0.001;
    el = energyLow;
    phLow = trialPh(el);
    phHigh = phLow;
    while (phLow*phHigh>0.0) {
      el = el + de;
      phHigh = trialPh(el);
    }

    eh = el;
    el = el - de;
    em = el;
    while (eh-el>error) {
      em = (el+eh)/2.0;
      phMid = trialPh(em);
      if (phMid*phLow>0.0) {
        el = em;
      } else {
        eh = em;
      }
    }
    return( em );
  }

  double trialPh(double energy) {
    int i;
    double h2;

    ph[0] = 0.0;
    ph[1] = dx;
    h2 = dx*dx;
    for (i=1; i<NNx-1; i++) {
      ph[i+1] = 2.0*ph[i]-ph[i-1]+2.0*(vv[i]-energy)*ph[i]*h2;
    }
    return( ph[NNx-1] );
  }

  void normalize() {
    int i;
    double s,a;

    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      s += ph[i]*ph[i]*dx;
    }
    a = Math.sqrt(1.0/s);
    for (i=1; i<NNx-1; i++) {
      ph[i] = a*ph[i];
    }
  }


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

}