/* oscillators time dependent Khon-Sham 1D */
/* coded by Ikeuchi Mitsuru     */
/* ver 0.0.1   2004.11.18       */

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

public class oscillatorsTDKS1D extends Applet
  implements ActionListener, AdjustmentListener, Runnable {

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

  Button bt_reset, bt_start, bt_stop, bt_loss;
  Scrollbar sc_x0, sc_x1;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int started = 1;
  int lossSW = 0;
  int setkFlag = 0;
  double xPos0 = 8.0;
  double xPos1 = 3.0;
  double xPos2 = 5.5;
  double xPos3 =10.0;
  double xPos4 =12.0;

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

  int NNx = 256+1;
  int NNe = 5;
  double energy[] = new double[NNe];  
  double psi[][][] = new double[NNe][NNx][2];
  double vv[] = new double[NNx];

  double wrk[][] = new double[NNx][2];

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

  int NN = NNx+1;
  double aaRe[] = new double[NN];
  double aaIm[] = new double[NN];
  double bbRe[] = new double[NN];
  double bbIm[] = new double[NN];
  double bRe[] = new double[NN];
  double bIm[] = new double[NN];
  double uRe[] = new double[NN];
  double uIm[] = new double[NN];

  int xpts[] = new int[NNx+1];
  int ypts[] = new int[NNx+1];

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

  public void init() {

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

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

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

    bt_loss = new Button("loss");
    bt_loss.addActionListener(this);

    sc_x0= new Scrollbar(Scrollbar.HORIZONTAL,80,10,30,140);
    sc_x0.addAdjustmentListener(this);

    sc_x1= new Scrollbar(Scrollbar.HORIZONTAL,30,10,30,140);
    sc_x1.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    pnl.add(bt_reset);
    pnl.add(bt_start);
    pnl.add(bt_stop);
    pnl.add(new Label("  "));
    pnl.add(new Label("  "));
    pnl.add(bt_loss);
    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 actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      started = -1;
	  lossSW = 0;
    } else if(ev.getSource() == bt_start){ 
      started = 1;
    } else if(ev.getSource() == bt_stop){
      started = 0;
    } else if(ev.getSource() == bt_loss){
      if (lossSW==0) {
        lossSW = 1;
      } else {
        lossSW = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_x0) { 
      xPos0 = 0.1*(double)(sc_x0.getValue());
      setInitialCondition();
    } else if (ev.getSource() == sc_x1) { 
      xPos1 = 0.1*(double)(sc_x1.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,ie,gb,gbx,gby,gby2;
    double v1,v2,ph2,y,z;
    double ek0,ev0,ek1,ev1;

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

    gbx = 20;
    gby = 220;

    gOff.setColor(Color.getHSBColor(0.10f,0.5f,0.9f));
    for (i=1; i<NNx-1; i++) {
      ph2 = 100.0*lylz*rho[i];
      gby2 = gby-(int)(vv[i]);
      if (ph2>=1.0) {
        gOff.drawLine(gbx+i,gby2-(int)(ph2), gbx+i, gby2);
      }
    }

    for (ie=0; ie<NNe; ie++) {
      drawPsi(ie, 100, gbx, gby-(int)(20.0*energy[ie]) );
    }

    gOff.setColor(Color.getHSBColor(0.33f,0.9f,0.9f));
    drawPot(vv, 1.0, gbx, gby);
    gOff.drawString("Veff", 280, 160);

    gOff.setColor(Color.getHSBColor(0.66f,0.9f,0.9f));
    drawPot(vvh, 100.0, gbx, gby);
    gOff.drawString("x100 VH", 280, 180);

    gOff.setColor(Color.getHSBColor(0.50f,0.9f,0.9f));
    drawPot(vvx, 100.0, gbx, gby);
    gOff.drawString("x100 Vx", 280, 200);

    gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.9f));
    drawPot(vvc, 100.0, gbx, gby);
    gOff.drawString("x100 Vc", 280, 220);
	
    for (ie=0; ie<NNe; ie++) {
      drawPhase(ie, 100, gbx, gby+80-(int)(10.0*ie));
    }
    gOff.drawString("phase", 280, 280);

    gb = 360;
    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*2418.9+0.5)/100.0+" as",10,40);

    gOff.setColor(Color.black);
    gOff.drawString("lossSW="+lossSW+" ",630/6*5+10,40);

    gOff.drawString("Vext(x) = 1.0*x^2 ",gb,60);

    gOff.setColor(Color.blue);
    for (ie=0; ie<NNe; ie++) {
      gOff.drawString("|<"+ie+"|0>|="+(float)(innerNorm(psi[ie],psi[0]))+" ",gb,100+15*ie);
	}
    for (ie=0; ie<NNe; ie++) {
	  gOff.setColor(Color.getHSBColor((float)(0.15+0.15*ie),0.9f,(float)(0.5+0.3*psi[ie][i][1])));
      gOff.drawString("E"+ie+"="+(float)(energy[ie])+" ",gb,280-15*ie);
    }

  }


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

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

  void drawPsi(int ie, double mag, int gx, int gy) {
    int i;
    double ph2;

    for (i=0; i<NNx; i++) { 
      ph2 = mag*(psi[ie][i][0]*psi[ie][i][0]+psi[ie][i][1]*psi[ie][i][1]);
      xpts[i] = gx+i;
      ypts[i] = gy-(int)(ph2);
    }
    for (i=0; i<NNx-2; i++) {
      gOff.setColor(Color.getHSBColor((float)(0.15+0.15*ie),0.9f,(float)(0.5+0.3*psi[ie][i][1])));
      gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
    }
  }

  void drawPhase(int ie, double mag, int gx, int gy) {
    int i;
    double y,z;

    for (i=0; i<NNx; i++) { 
      y = 20.0*psi[ie][i][0];
      z = 20.0*psi[ie][i][1];
      xpts[i] = gx+i-(int)(0.5*z);
      ypts[i] = gy+(int)(0.5*z)-(int)(y);
    }
    for (i=0; i<NNx-2; i++) {
      gOff.setColor(Color.getHSBColor((float)(0.15+0.15*ie),0.9f,(float)(0.5+0.3*psi[ie][i][1])));
      gOff.drawLine(xpts[i],ypts[i],xpts[i+1],ypts[i+1]);
    }
  }

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

  double norm(double ph[][]) {
    int i;
    double p;

    p = 0.0;
    for (i=0; i<NNx; i++) {
      p += (ph[i][0]*ph[i][0]+ph[i][1]*ph[i][1])*dx;
    }
    return ( p );
  }

  double kineticEnergy(double ph[][]) {
    int i;
    double s,hphRe,hphIm;

    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      hphRe = (2.0*ph[i][0]-ph[i+1][0]-ph[i-1][0])/(2.0*dx*dx);
      hphIm = (2.0*ph[i][1]-ph[i+1][1]-ph[i-1][1])/(2.0*dx*dx);
      s += (ph[i][0]*hphRe + ph[i][1]*hphIm)*dx;
    }
    return ( s );
  }

  double momentum(double ph[][]) {
    int i;
    double s,pphRe,pphIm;

    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      pphRe = (ph[i+1][1]-ph[i-1][1])/(2.0*dx);
      pphIm = (-ph[i+1][0]+ph[i-1][0])/(2.0*dx);
      s += (ph[i][0]*pphRe + ph[i][1]*pphIm)*dx;
    }
    return ( s );
  }

  double potentialEnergy(double ph[][]) {
    int i;
    double s;

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

  double meanPosX(double ph[][]) {
    int i;
    double s;

    s = 0.0;
    for (i=1; i<NNx-1; i++) {
      s += i*dx*(ph[i][0]*ph[i][0] + ph[i][1]*ph[i][1])*dx;
    }
    return ( s - 8.0 );
  }

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

  void steepestDescent(double ph[][], double v[], double dump) {
    int i;
    double h2,ee;

    h2 = 2.0*dx*dx;
    ee = kineticEnergy(ph)+potentialEnergy(ph);

    for (i=1; i<NNx-1; i++) {
      wrk[i][0] = (2.0*ph[i][0]-ph[i+1][0]-ph[i-1][0])/h2+(v[i]-ee)*ph[i][0];
      wrk[i][1] = (2.0*ph[i][1]-ph[i+1][1]-ph[i-1][1])/h2+(v[i]-ee)*ph[i][0];

    }

    for (i=1; i<NNx-1; i++) {
      ph[i][0] -= dump*wrk[i][0];
      ph[i][1] -= dump*wrk[i][1];
    }
  }

/*-------------------------- set loss energy ----------*/

  void setLoss(double ph[][]) {
    steepestDescent(ph, vv, 0.003);
  }

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

  void GramSchmidt(int stateMax) {
    int i,istate,ist;
    double pRe,pIm;

    normalize(psi[0]);

    for (istate=1; istate<stateMax; istate++) {
      for (ist=0; ist<istate; ist++) {
        pRe = 0.0; pIm = 0.0;
        for (i=1;i<NNx-1;i++) {
          pRe += (psi[ist][i][0]*psi[istate][i][0] + psi[ist][i][1]*psi[istate][i][1])*dx;
          pIm += (psi[ist][i][0]*psi[istate][i][1] - psi[ist][i][1]*psi[istate][i][0])*dx;
        }
        for (i=1;i<NNx-1;i++) {
          psi[istate][i][0] -= (pRe*psi[ist][i][0]-pIm*psi[ist][i][1]);
          psi[istate][i][1] -= (pRe*psi[ist][i][1]+pIm*psi[ist][i][0]);
        }
      }
      normalize(psi[istate]);
    }
  }

  void normalize(double f[][]) {
    int i;
    double a;

    a = Math.sqrt(norm(f));
    for (i=1;i<NNx-1;i++) {
      f[i][0] = f[i][0]/a;
      f[i][1] = f[i][1]/a;
    }
  }

  double innerNorm(double f[][], double g[][]) {
    int i;
    double pRe,pIm;

    pRe = 0.0;
    pIm = 0.0;
    for (i=1;i<NNx-1;i++) {
      pRe += (f[i][0]*g[i][0] + f[i][1]*g[i][1])*dx;
      pIm += (f[i][0]*g[i][1] - f[i][1]*g[i][0])*dx;
    }
  
    return ( Math.sqrt(pRe*pRe+pIm*pIm) );
  }

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

  void setInitialCondition() {
    t = 0.0;
    setWave(psi[0],xPos0, Math.sqrt(1/Math.sqrt(8.0)), 0.0);
    setWave(psi[1],xPos1, Math.sqrt(1/Math.sqrt(8.0)), 0.0);
    setWave(psi[2],xPos2, Math.sqrt(1/Math.sqrt(8.0)), 0.0);
    setWave(psi[3],xPos3, Math.sqrt(1/Math.sqrt(8.0)), 0.0);
    setWave(psi[4],xPos4, Math.sqrt(1/Math.sqrt(8.0)), 0.0);
    GramSchmidt(5);
    setPotential();
    setRho(NNe);
    poisson(10000);
  }


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

  void setVext() {
    int i;
    double x,x0;

    x0 = 8.0;
    for (i=1;i<NNx-1;i++) {
      x = i*dx;
      vvext[i] = 1.0*(x-x0)*(x-x0);
    }
  }


  void setWave(double ph[][],double wavePos, double waveWidth, double momentum) {
    int i;
    double a,x,phAb,phPh;

    a = Math.pow(2.0*Math.PI*waveWidth*waveWidth,-0.25);
    for (i=1;i<NNx-1;i++) {
      x = i*dx;
      phAb = a*Math.exp(-((x-wavePos)/(2.0*waveWidth))*((x-wavePos)/(2.0*waveWidth)));
      phPh = momentum*x;
      ph[i][0] = phAb*Math.cos(phPh);
      ph[i][1] = phAb*Math.sin(phPh);
    }
    ph[0][0] = 0.0;
    ph[0][1] = 0.0;
    ph[NNx-1][0] = 0.0;
    ph[NNx-1][1] = 0.0;

    a = Math.sqrt(norm(ph));
    for (i=1;i<NNx-1;i++) {
      ph[i][0] = ph[i][0]/a;
      ph[i][1] = ph[i][1]/a;
    }
  }


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

  void timeEvolution() {
    int i,ie;

    if (started==1) {
      for (i=0; i<20; i++) {

        t = t + 0.5*dt;
        for (ie=0; ie<NNe; ie++) {
          kxStep(psi[ie],0.5*dt);
        }

        setVeff();

        for (ie=0; ie<NNe; ie++) {
          phaseStep(psi[ie],dt);
        }

        t = t + 0.5*dt;
        for (ie=0; ie<NNe; ie++) {
          kxStep(psi[ie],0.5*dt);
        }
	
        for (ie=0; ie<NNe; ie++) {
          energy[ie]=kineticEnergy(psi[ie])+potentialEnergy(psi[ie]);
        }


      }

    } else if (started==-1) {
      t = 0.0;
      setInitialCondition();
      started = 0;
    }
    if (lossSW==1) {
      setLoss(psi[0]);
      setLoss(psi[1]);
      setLoss(psi[2]);
      setLoss(psi[3]);
      setLoss(psi[4]);
      GramSchmidt(5);
    }
  }


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

  void setVeff() {
    int i;

    setRho(NNe);
    poisson(100);
    setVcx();
    for (i=1;i<NNx-1;i++) {
      vv[i] = vvext[i]+vvh[i]+vvx[i]+vvc[i];
    }
  }

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

    g = 1.0;

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

  void poisson(int iterMax) {
    int i,iter;
    double h2;

    h2 = dx*dx;
    for (iter=1; iter<iterMax; iter++) {
      for (i=1; i<NNx-1; i++) {
        vvh[i] = vvh[i]+0.66*(vvh[i+1]+vvh[i-1]-2.0*vvh[i]+h2*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 );
  }


/* ----------------------------- time step -------------- */


  void kxStep(double ph[][],double deltat) {
    int i;
    double a,aaAb,auAb;

    a = 4.0*dx*dx/deltat;
    for (i=1; i<NNx-1; i++) {
      bbRe[i] = 2.0;
      bbIm[i] = a;
      aaRe[i] = -2.0;
      aaIm[i] = a;
      bRe[i] = bbRe[i]*ph[i][0]-bbIm[i]*ph[i][1] - ph[i+1][0] - ph[i-1][0];
      bIm[i] = bbRe[i]*ph[i][1]+bbIm[i]*ph[i][0] - ph[i+1][1] - ph[i-1][1];
    }

    aaAb = aaRe[1]*aaRe[1]+aaIm[1]*aaIm[1];
    uRe[1] = aaRe[1]/aaAb;
    uIm[1] = -aaIm[1]/aaAb;
    ph[1][0] = bRe[1]*uRe[1] - bIm[1]*uIm[1];
    ph[1][1] = bIm[1]*uRe[1] + bRe[1]*uIm[1];

    for (i=2; i<NNx-1; i++) {
      auAb = (aaRe[i]-uRe[i-1])*(aaRe[i]-uRe[i-1])+(aaIm[i]-uIm[i-1])*(aaIm[i]-uIm[i-1]);
      uRe[i] = (aaRe[i]-uRe[i-1])/auAb;
      uIm[i] = -(aaIm[i]-uIm[i-1])/auAb;
      ph[i][0] = (bRe[i]-ph[i-1][0])*uRe[i] - (bIm[i]-ph[i-1][1])*uIm[i];
      ph[i][1] = (bRe[i]-ph[i-1][0])*uIm[i] + (bIm[i]-ph[i-1][1])*uRe[i];
    }

    for (i=NNx-3; i>=1; i--) {
      ph[i][0] -= ph[i+1][0]*uRe[i] - ph[i+1][1]*uIm[i];
      ph[i][1] -= ph[i+1][0]*uIm[i] + ph[i+1][1]*uRe[i];
    }

  }

  void phaseStep(double ph[][],double deltat) {
    int i;
    double th,cs,sn,phr,phi;

    for (i=1; i<NNx-1; i++) {
      th = deltat*vv[i];
      cs = Math.cos(th); sn = Math.sin(th);
      phr = ph[i][0];
      phi = ph[i][1];
      ph[i][0] = cs*phr+sn*phi;
      ph[i][1] = cs*phi-sn*phr;
    }
  }

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

}