/* H2 molecules - first principle molecular dynamics */
/* DFT-LDA, steepest descent - Gram-Schumidt, SCF    */
/* coded by Ikeuchi Mitsuru    */
/* ver 0.0.1   2004.09.04      */
/* ver 0.0.2   2004.09.25  setVcx() bug fixed */

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

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

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

  Choice cvw, ch_md, ch_noe;
  Button bt_reset, bt_start, bt_stop;
  Scrollbar spx,spy,spz;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 30;

  int dispMode = 1;
  int dispObj = 0;
  double px0 = 0.0;
  double py0 = 0.0;
  double pz0 = 0.0;
  int started = 1;

  double NOE = 2.0;
  double ppDistance = 2.0;
  double sdE0mem = 0.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/4.0;
  double dy = 1.0/4.0;
  double dz = 1.0/4.0;

  double dt = 2.0;

  int NNx = 64+1;
  int NNy = 64+1;
  int NNz = 64+1;
  int NN = NNx+1; /* max(NNx,NNy,NNz)+1 */
  int NNr = 500;

  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 rho[][][] = new double[NNx][NNy][NNz];

  double sdEnergy[] = new double[20];
  double sdState[][][][] = new double[20][NNx][NNy][NNz];


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

/* ----------------------- md global------ */

  double AUx = 5.291772e-11; /* length (m) */
  double AUt = 2.418884e-17; /* time (s) */
  double AUe = 4.359748e-18; /* energy (J) */

  int Nmt = 2;
  double qq[] = new double[Nmt]; /* charge */
  double mm[] = new double[Nmt];
  double xx[] = new double[Nmt];
  double yy[] = new double[Nmt];
  double zz[] = new double[Nmt];   
  double vx[] = new double[Nmt];
  double vy[] = new double[Nmt];
  double vz[] = new double[Nmt];
  double ffx[] = new double[Nmt];
  double ffy[] = new double[Nmt];
  double ffz[] = new double[Nmt];

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

  public void init() {

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

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

    cvw = new Choice();
    cvw.addItem("ph-0");cvw.addItem("ph-1");
    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"); ch_md.addItem("VH x2"); ch_md.addItem("Vx x2");
    ch_md.addItem("Vc x2"); ch_md.addItem("Veff");
    ch_md.addItemListener(this);
    ch_md.select("phase");

    ch_noe = new Choice();
    ch_noe.addItem("n=1.5");ch_noe.addItem("n=2");
    ch_noe.addItemListener(this);
    ch_noe.select("n=2");

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

    
    spx= new Scrollbar(Scrollbar.HORIZONTAL,200,10,100,410);
    spx.addAdjustmentListener(this);

    spy= new Scrollbar(Scrollbar.HORIZONTAL,0,10,0,110);
    spy.addAdjustmentListener(this);

    spz= new Scrollbar(Scrollbar.HORIZONTAL,0,10,0,110);
    spz.addAdjustmentListener(this);

    addMouseListener(this);  
    addMouseMotionListener(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(cvw);
    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() == cvw){
      dispObj = cvw.getSelectedIndex();
    } else if (ev.getSource() == ch_md){
      dispMode = ch_md.getSelectedIndex();
    } else if (ev.getSource() == ch_noe){
      NOE = 0.5*ch_noe.getSelectedIndex()+1.5;
    }
  }

  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) { 
      ppDistance = 0.01*(double)(spx.getValue());
	  setNuclearPosition();
	  setVext();
    } else if (ev.getSource() == spy) { 
      py0 = 0.1*(double)(spy.getValue());
    } else if (ev.getSource() == spz) { 
      pz0 = 0.1*(double)(spz.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;
    double ev,ek;
    double dens;

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

    dens = 500.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, 5.0);

    } else if (dispMode==4) {

      gridFnPlot(vvh, NNz/2, 10.0);

    } else if (dispMode==5) {

      gridFnPlot(vvx, NNz/2, 10.0);

    } else if (dispMode==6) {

      gridFnPlot(vvc, NNz/2, 10.0);

    } else if (dispMode==7) {

      gridFnPlot(vv, NNz/2, 5.0);

    }

    gOff.setColor(Color.black);
    gOff.drawString("  t="+t+" au   "+(float)(1.0e15*t*AUt)+" fs",30,40);

    gOff.drawString("orbit",630/6*4+10,40);
    gOff.drawString("view",630/6*5+10,40);

    gb = 320;
    gOff.setColor(Color.black);
    gOff.drawString("box="+(int)(NNx*dx)+"x"+ +(int)(NNy*dy)+"x"+(int)(NNz*dz)+" au",gb,60);
    gOff.drawString("dx="+(int)(dx*100.0)/100.0+" au",gb+150,60);
    gOff.drawString("Vext(r)= -1/|R0-r| -1/|R1-r|,",gb,80);
    gOff.drawString("            (|Ri-r|>=dx/2)",gb,100);

    gOff.setColor(Color.blue);
    gOff.drawString("x[0]="+(float)(xx[0]/AUx)+"",gb,120);
    gOff.drawString("vx[0]="+(float)vx[0]+"",gb,140);
    gOff.drawString("fx[0]="+(float)ffx[0]+"",gb,160);

    gOff.drawString("x[1]="+(float)(xx[1]/AUx)+" au",gb+150,120);
    gOff.drawString("vx[1]="+(float)vx[1]+" m/s",gb+150,140);
    gOff.drawString("fx[1]="+(float)ffx[1]+" N",gb+150,160);

    gOff.drawString("distance="+(float)(distanceAU())+" au",gb,220);
    gOff.drawString("total Energy="+(float)totalEnergy()+" au",gb,240);

    for (i=0; i<2; i++) {
	  if (i<1) {
	    gOff.setColor(Color.red);
	  } else {
	    gOff.setColor(Color.blue);
	  }
      gOff.drawString("E["+i+"]="+(float)(sdEnergy[i])+"",gb,260+20*i);
	}
  }

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

  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;

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

    for (i=1; i<NNx-1; i+=2) {
      for(j=1; j<NNy-1; j +=2) {
        for (k=1; k<NNz-1; k+=2) {
          if (vv[i][j][k]>0.0) {
            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;
            ir = 2;
            ix = (int)(ppx)+gbx;
            iy = (int)(ppy)+gby;
            gOff.setColor(Color.getHSBColor(0.33f, 0.9f,0.9f));
            gOff.fillOval(ix-ir/2,iy-ir/2,ir,ir);
          }
          dens = 1.0*sdState[istate][i][j][k]*sdState[istate][i][j][k];
          ir = (int)(weight*dens+0.5);
          if (ir>0) {
            if (ir>10) {
              ir = 10;
            }
            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;

            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 (sdState[istate][i][j][k]>=0) {
                col = 0.33;
              } else {
                col = 0.01;
              }
            }
			
            gOff.setColor(Color.getHSBColor((float)col, 0.9f,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 = 4.0;

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

    for(j=0; j<NNy; j += 2) {
      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 += 2) {
      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 = 60;
    gby = 80;
    mag = 3.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;
    double p,d2ph,h2;

    h2 = dx*dx;
    p = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;k++) {
          d2ph = (ph[i+1][j][k]+ph[i-1][j][k]+ph[i][j+1][k]+ph[i][j-1][k]+ph[i][j][k+1]+ph[i][j][k-1]-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=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;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 s,sei,enn;

    sei = 2.0*sdEnergy[0];
    enn = 1.0/distanceAU();
    s = 0.0;
    for (i=1;i<NNx-1;i++) {
      for (j=1;j<NNy-1;j++) {
        for (k=1;k<NNz-1;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 );
  }

  double distanceAU() {
    double d;

    d = Math.sqrt( (xx[0]-xx[1])*(xx[0]-xx[1])+(yy[0]-yy[1])*(yy[0]-yy[1])+(zz[0]-zz[1])*(zz[0]-zz[1]) );
  
    return ( d/AUx );

  }

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

  void setInitialCondition() {

    t = 0.0;
    setNuclearPosition();
    setVext();

    initSD(2);
  }


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

  void timeEvolution() {
    int i;
    double edE0;

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

      setVeff();
      SD(2, vv, 2, 0.03);
      sortState(2);

      if (Math.abs(sdEnergy[0]-sdE0mem)<0.0001) {
        t += dt;
        moveNuclear();
        setVext();
      }
      sdE0mem = sdEnergy[0];

    }
  }

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

    for (ist=0; ist<maxState-1; ist++) {
      if (sdEnergy[ist]>sdEnergy[ist+1]+0.001) {
        for (i=1;i<NNx-1;i++) {
          for (j=1;j<NNy-1;j++) {
            for (k=1;k<NNz-1;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;
      }
    }

  }


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

  void setVext() {
    int i,j,k,inuc;
    double x,y,z,x0,y0,z0,r;

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        for (k=0;k<NNz;k++) {
          vvext[i][j][k] = 0.0;
          x = i*dx;
          y = j*dy;
          z = k*dz;
          for (inuc=0; inuc<Nmt; inuc++) {
            x0 = xx[inuc]/AUx; y0 = yy[inuc]/AUx; z0 = zz[inuc]/AUx;
            r = Math.sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)+(z-z0)*(z-z0));
            if (r<0.5*dx) r = 0.5*dx;
            vvext[i][j][k] += -qq[inuc]/r;
          }
        }
      }
    }

    for (i=0;i<NNx;i++) {
      for (j=0;j<NNy;j++) {
        k = 0;
        vvh[i][j][k] = -vvext[i][j][k]; if (vvh[i][j][k]>1.0) vvh[i][j][k]=1.0;
        k = NNz-1;
        vvh[i][j][k] = -vvext[i][j][k]; if (vvh[i][j][k]>1.0) vvh[i][j][k]=1.0;
      }
    }
    for (i=0;i<NNx;i++) {
      for (k=0;k<NNz;k++) {
        j = 0;
        vvh[i][j][k] = -vvext[i][j][k]; if (vvh[i][j][k]>1.0) vvh[i][j][k]=1.0;
        j = NNy-1;
        vvh[i][j][k] = -vvext[i][j][k]; if (vvh[i][j][k]>1.0) vvh[i][j][k]=1.0;
      }
    }
    for (j=0;j<NNy;j++) {
      for (k=0;k<NNz;k++) {
        i = 0;
        vvh[i][j][k] = -vvext[i][j][k]; if (vvh[i][j][k]>1.0) vvh[i][j][k]=1.0;
        i = NNx-1;
        vvh[i][j][k] = -vvext[i][j][k]; if (vvh[i][j][k]>1.0) vvh[i][j][k]=1.0;
      }
    }
  }

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

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

    setRho(1);
    poisson(10);
    setVcx();

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

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

    g = NOE;

    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          rho[i][j][k] = 0.0;
          for (ie=0; ie<istateMax; ie++) {
            rho[i][j][k] += g*(sdState[ie][i][j][k]*sdState[ie][i][j][k]);
          }
        }
      }
    }
  }

  void poisson(int iterMax) {
    int i,j,k,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++) {
        for (j=1; j<NNy-1; j++) {
          for (k=1; k<NNz-1; k++) {
            vvh[i][j][k] = vvh[i][j][k]+w*(vvh[i+1][j][k]+vvh[i-1][j][k]+vvh[i][j+1][k]+vvh[i][j-1][k]+vvh[i][j][k+1]+vvh[i][j][k-1]-6.0*vvh[i][j][k] +h2*rho[i][j][k]);
          }
        }
      }
    }
  }

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

  void setVcx() {
    int i,j,k;
    double c1,rh,drho;

    c1 = -0.984745022;
    drho = 0.01;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          rh = rho[i][j][k];
          vvx[i][j][k] = c1*Math.pow(rh,0.33333333);
          vvc[i][j][k] = eeCorrelation(rh)+rh*(eeCorrelation(rh+drho)-eeCorrelation(rh))/drho;
        }
      }
    }
  }

  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=1; i<NNx-1; i++) {
          for (j=1; j<NNy-1; j++) {
            for (k=1; k<NNz-1; 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=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; 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=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          s += ph[i][j][k]*ph[i][j][k]*dv;
        }
      }
    }
    a = Math.sqrt(1.0/s);
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; 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=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]);
    }
  }

  double steepestDescent(double psi[][][], double v[][][], double dump) {
    int i,j,k;
    double h2,ei,wk;

    h2 = 2.0*dx*dx;
    ei = energyOf(psi,vv);
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          psi[i][j][k] -= dump*((6.0*psi[i][j][k]-psi[i+1][j][k]-psi[i-1][j][k]-psi[i][j+1][k]-psi[i][j-1][k]-psi[i][j][k+1]-psi[i][j][k-1])/h2 + (v[i][j][k]-ei)*psi[i][j][k]);
        }
      }
    }
    normalize(psi);
    return (ei);
  }

  double energyOf(double psi[][][], double v[][][]) {
    int i,j,k;
    double h2,dv,s,sn;

    h2 = 2.0*dx*dx;
    dv = dx*dy*dz;
    s = 0.0; sn=0.0;
    for (i=1; i<NNx-1; i++) {
      for (j=1; j<NNy-1; j++) {
        for (k=1; k<NNz-1; k++) {
          s += psi[i][j][k]*((6.0*psi[i][j][k]-psi[i+1][j][k]-psi[i-1][j][k]-psi[i][j+1][k]-psi[i][j-1][k]-psi[i][j][k+1]-psi[i][j][k-1])/h2 + v[i][j][k]*psi[i][j][k])*dv;
          sn += psi[i][j][k]*psi[i][j][k]*dv;
        }
      }
    }
    return ( s/sn );
  }


/*----------------------- nuclear motion (classical MD) -------------*/

  void setNuclearPosition(){
    int i; 
    double s;

    xx[0] = AUx*(8.0-0.5*ppDistance);
    yy[0] = AUx*(8.0);
    zz[0] = AUx*(8.0);
    mm[0] = 1.6726485e-27;
    qq[0] = 1.0;

    xx[1] = AUx*(8.0+0.5*ppDistance);
    yy[1] = AUx*(8.0);
    zz[1] = AUx*(8.0);
    mm[1] = 1.6726485e-27;
    qq[1] = 1.0;

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

  void moveNuclear() {
    int i; 
    double dtSI, dtSIv2, rr, tmp;

    dtSI= AUt*dt;
    dtSIv2 = dtSI/2.0;

    for (i=0; i<Nmt; i++) {
      vx[i] += dtSIv2*ffx[i]/mm[i];
      vy[i] += dtSIv2*ffy[i]/mm[i];
      vz[i] += dtSIv2*ffz[i]/mm[i];
      xx[i] += vx[i]*dtSI;
      yy[i] += vy[i]*dtSI;
      zz[i] += vz[i]*dtSI;
    }
    forceCalc();
    for (i=0; i<Nmt; i++) {
      vx[i] += dtSIv2*ffx[i]/mm[i];
      vy[i] += dtSIv2*ffy[i]/mm[i];
      vz[i] += dtSIv2*ffz[i]/mm[i];
    }
  }

  void forceCalc() {
    int i,j;
    double rij,f,fxij,fyij,fzij;

    for(i=0;i<Nmt;i++) {
      ffx[i] = electronForceX(xx[i], yy[i], zz[i]);
      ffy[i] = electronForceY(xx[i], yy[i], zz[i]);
      ffz[i] = electronForceZ(xx[i], yy[i], zz[i]);
    }

    for(i=0;i<Nmt;i++) {
      for(j=i+1;j<Nmt;j++) {
        rij = Math.sqrt((xx[i]-xx[j])*(xx[i]-xx[j])+(yy[i]-yy[j])*(yy[i]-yy[j])+(zz[i]-zz[j])*(zz[i]-zz[j]));
        f = force(rij,i,j);
        fxij = f*(xx[i]-xx[j])/rij;
        fyij = f*(yy[i]-yy[j])/rij;
        fzij = f*(zz[i]-zz[j])/rij;

        ffx[i] += fxij;
        ffy[i] += fyij;
        ffz[i] += fzij;

        ffx[j] += -fxij;
        ffy[j] += -fyij;
        ffz[j] += -fzij;
      }
    }
  }

  double force(double r, int ii, int jj) {
    double qi,qj;

    qi = qq[ii]*1.602e-19;
    qj = qq[jj]*1.602e-19;

    return ( qi*qj/(4.0*3.1416*8.854e-12*r*r) );
  }

  double electronForceX(double x, double y, double z) {
    int j,k;
    double r;

    j = (int)((y/AUx)/dy+0.5);
    k = (int)((z/AUx)/dz+0.5);

	r = -(vvelec(x/AUx+1.0*dx,j,k)-vvelec(x/AUx-1.0*dx,j,k))/(2.0*dx);
    return( -r*AUe/AUx );
  }


  double vvelec(double x,int j, int k) {
    int i;
    double a;

    i = (int)(x/dx);
    a = x - i*dx;

    return ( ((dx-a)*vvh[i][j][k] + a*vvh[i+1][j][k])/dx );
  }

  double electronForceY(double x, double y, double z) {
    int i,k;
    double r;

    i = (int)((x/AUx)/dy+0.5);
    k = (int)((z/AUx)/dz+0.5);

	r = -(vvelecY(i,y/AUx+1.0*dy,k)-vvelecY(i,y/AUx-1.0*dy,k))/(2.0*dx);
    return( -r*AUe/AUx );
  }


  double vvelecY(int i,double y, int k) {
    int j;
    double a;

    j = (int)(y/dy);
    a = y - j*dy;

    return ( ((dy-a)*vvh[i][j][k] + a*vvh[i][j+1][k])/dy );
  }


  double electronForceZ(double x, double y, double z) {
    int i,j;
    double r;

    i = (int)((x/AUx)/dx+0.5);
    j = (int)((y/AUx)/dy+0.5);

	r = -(vvelecZ(i,j,z/AUx+1.0*dz)-vvelecZ(i,j,z/AUx-1.0*dz))/(2.0*dz);
    return( -r*AUe/AUx );
  }


  double vvelecZ(int i, int j, double z) {
    int k;
    double a;

    k = (int)(z/dz);
    a = z - k*dz;

    return ( ((dz-a)*vvh[i][j][k] + a*vvh[i][j][k+1])/dz );
  }

/*----- end of nuclear motion -----*/

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

}