/* applet No. 981
 *
 * Lennard-Jones potential --> Morse potential 
 * fitting - Metropolis sampling - Monte-Carlo simulation
 *
 * Created by Ikeuchi Mitsuru on July 22 2006.
 * Copyright (c) 2006 Ikeuchi Mitsuru. All rights reserved.
 *
 * ver 0.0.1   2006.07.22  created
 *
 */

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

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

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

  Choice ch_mat, ch_walker, ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_beta, sc_ss0, sc_ss1;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int dispWidth = 630;
  int dispHeight = 440;

  int sleepTime = 50;

  double t = 0.0;
  double dt = 1.0;

  int resetSW = 0;
  int started = 1;
  int dispMode = 2;
  int maxWalker = 10;
  int maxWalkerPow = 1;

  double ddMorse = 0.01;
  double r0Morse = 4.0;
  double aaMorse = 1.0;


  int nnWidth = 400;
  int nnHight = 300;

  double ss0 = 0.1;
  double ss1 = 0.1;
  double betaInv = 0.5;

  int varMax = 3;
  double xx[][] = new double[1000][varMax];
  double xxTrial[] = new double[varMax];
  double xxWork[] = new double[varMax];
  double xxMin[] = new double[varMax];
  double xxMax[] = new double[varMax];
  double sss[] = new double[varMax];
  double ppp[] = new double[1000];


  double AU = 1.661e-27;
  double EE = 1.602e-19;
  double KB = 1.38e-23;
  double AA = 1.0e-10;
  double FS = 1.0e-15;
  double LJdata[][] = {
  /*             mass       charge   E(in kB)  sigma(A)   dt(s)  */
  /* 0 He*/   {   4.003*AU,  0.0*EE,  10.2*KB, 2.576*AA,  5.0*FS },
  /* 1 Ne*/   {  20.183*AU,  0.0*EE,  36.2*KB, 2.976*AA, 10.0*FS },
  /* 2 Ar*/   {  39.948*AU,  0.0*EE, 124.0*KB, 3.418*AA, 20.0*FS },
  /* 3 Kr*/   {  83.50 *AU,  0.0*EE, 190.0*KB, 3.610*AA, 20.0*FS },
  /* 4 Xe*/   { 131.30 *AU,  0.0*EE, 229.0*KB, 4.055*AA, 20.0*FS },
  /* 5 Hg*/   { 200.59 *AU,  0.0*EE, 851.0*KB, 2.898*AA, 20.0*FS },
  /* 6 H2*/   {   2.016*AU,  0.0*EE,  33.3*KB, 2.968*AA,  5.0*FS },
  /* 7 N2*/   {  28.013*AU,  0.0*EE,  91.5*KB, 3.681*AA, 20.0*FS },
  /* 8 O2*/   {  31.999*AU,  0.0*EE, 113.0*KB, 3.433*AA, 20.0*FS },
  /* 9 H2O*/  {  18.015*AU,  0.0*EE, 809.1*KB, 2.641*AA, 10.0*FS },
  /*10 CH4*/  {  16.043*AU,  0.0*EE, 137.0*KB, 3.822*AA, 10.0*FS },
  /*11 CO2*/  {  44.010*AU,  0.0*EE, 190.0*KB, 3.996*AA, 20.0*FS },
  /*12 CO */  {  28.011*AU,  0.0*EE, 110.0*KB, 3.590*AA, 20.0*FS }
  };
  int kindMax = 13;
  int MAT1 = 2;

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

  public void init() {

    resize(dispWidth,dispHeight);
    setBackground(Color.white);

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

    ch_mat = new Choice();
    ch_mat.addItem("He");
	ch_mat.addItem("Ne");
	ch_mat.addItem("Ar");
    ch_mat.addItem("Kr");
	ch_mat.addItem("Xe");
	ch_mat.addItem("Hg");
    ch_mat.addItem("H2");
	ch_mat.addItem("N2");
	ch_mat.addItem("O2"); 
    ch_mat.addItem("H2O");
	ch_mat.addItem("CH4");
	ch_mat.addItem("CO2");
    ch_mat.addItem("CO");
    ch_mat.addItemListener(this);
    ch_mat.select("Ar");

    ch_walker = new Choice();
    ch_walker.addItem("1");
    ch_walker.addItem("10");
    ch_walker.addItem("100");
    //ch_walker.addItem("1000");
    ch_walker.addItemListener(this);
    ch_walker.select("10");

    ch_view = new Choice();
    ch_view.addItem("x1");
    ch_view.addItem("x10");
    ch_view.addItem("x100");
    ch_view.addItem("x1000");
    ch_view.addItemListener(this);
    ch_view.select("x100");

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

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 
    
    sc_beta = new Scrollbar(Scrollbar.HORIZONTAL,500,10,0,1010);
    sc_beta.addAdjustmentListener(this);

    sc_ss0 = new Scrollbar(Scrollbar.HORIZONTAL,1000,10,0,2010);
    sc_ss0.addAdjustmentListener(this);

    sc_ss1 = new Scrollbar(Scrollbar.HORIZONTAL,1000,10,0,2010);
    sc_ss1.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
    //pnl.add(new Label("  "));
    pnl.add(ch_mat);
    pnl.add(new Label("  "));
    pnl.add(ch_walker);
    pnl.add(sc_beta);
    pnl.add(sc_ss0);
    pnl.add(ch_view);
    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_mat) {
      MAT1 = ch_mat.getSelectedIndex();
      resetSW = 1;
    } else if (ev.getSource() == ch_walker) {
      maxWalkerPow = ch_walker.getSelectedIndex();
      if (maxWalkerPow==0) {
        maxWalker = 1;
      } else if (maxWalkerPow==1) {
        maxWalker = 10;
      } else if (maxWalkerPow==2) {
        maxWalker = 100;
      } else if (maxWalkerPow==3) {
        maxWalker = 1000;
      }
      resetSW = 1;
    } else if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_startStop){
      if (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_beta) { 
      betaInv = 0.001*(double)(sc_beta.getValue());
    } else if (ev.getSource() == sc_ss0) { 
      ss0 = 0.0001*(double)(sc_ss0.getValue());
      setsss(ss0);
    } else if (ev.getSource() == sc_ss1) { 
      ss1 = 0.0001*(double)(sc_ss1.getValue());
    }
  }

  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 gx,gy,gDiv;
    double col;

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

    gx = 40;
    gy = 100;

    drawGrid(gx,gy);
    if (dispMode==0) {
      drawTrial(gx,gy, 1.0, 0.33);
      drawLJ(gx,gy, 1.0, 0.01);

    } else if (dispMode==1) {
      drawTrial(gx,gy, 10.0, 0.33);
      drawLJ(gx,gy, 10.0, 0.01);

    } else if (dispMode==2) {
      drawTrial(gx,gy, 100.0, 0.33);
      drawLJ(gx,gy, 100.0, 0.01);

    } else if (dispMode==3) {
      drawTrial(gx,gy, 1000.0, 0.33);
      drawLJ(gx,gy, 1000.0, 0.01);

    } else if (dispMode==4) {
      ;

    } else if (dispMode==5) {
      ;

    }


    gOff.setColor(Color.black);
    gOff.drawString("t="+(int)(t+0.5)+" ",10,40);
    gOff.drawString("No. of walkers",dispWidth/6*2+0,40);
    gOff.drawString("kT="+(int)(betaInv*1000.0+0.5)/1000.0+"",dispWidth/6*3+10,40);
    gOff.drawString("s="+(int)(ss0*10000.0+0.5)/100.0+"%",dispWidth/6*4+10,40);
    gOff.drawString("view",dispWidth/6*5+10,40);

    gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.4f));
    gOff.drawString("Lennard-Jones V(r) = 4 eps { (sgm/r)^12-(sgm/r)^6 }",20,60);
	gOff.drawString("eps="+(int)((LJdata[MAT1][2]/1.6e-19)*10000.0+0.5)/10.0+"meV",400,60);
	gOff.drawString("sgm="+(int)((LJdata[MAT1][3]/1.0e-10)*100.0+0.5)/100.0+"A",500,60);

    gOff.setColor(Color.getHSBColor(0.33f,0.8f,0.4f));
    gOff.drawString("Morse V(r) = D exp( -2A(r-r0) ) - 2 D exp( -A(r-r0) )",20,80);
    //gOff.drawString("TBP V(r) = 2U -f/G + f^2/(4G^3), U=A*exp(-p(r/r0-1)), f=xx^2*exp(-2q(r/r0-1))",20,80);

    gOff.setColor(Color.getHSBColor(0.33f,0.8f,0.4f));
    printPotMorse(gx,gy);
  }

  void printPotMorse(int gbx, int gby) {
    int i,tb1,tb2,tb3,tb4;

    tb1=gbx+250; tb2=tb1+60; tb3=tb2+60; tb4=tb3+60;
    for (i=0; i<18; i++) {
      if (i>=maxWalker) break;
      gOff.drawString("D(meV)",tb1,gby+20);
      gOff.drawString("r0(A)",tb2,gby+20);
      gOff.drawString("A",tb3,gby+20);
      gOff.drawString("err",tb4,gby+20);
      gOff.drawString(""+(int)(xx[i][0]*1000000.0+0.5)/1000.0+"",tb1,gby+40+15*i);
      gOff.drawString(""+(int)(xx[i][1]*1000.0+0.5)/1000.0+"",tb2,gby+40+15*i);
      gOff.drawString(""+(int)(xx[i][2]*1000.0+0.5)/1000.0+"",tb3,gby+40+15*i);
      gOff.drawString(""+(int)(ppp[i]*1000.0+0.5)/1000.0+"",tb4,gby+40+15*i);
    }
  }

  void printPotTBP(int gbx, int gby) {
    int i,tb1,tb2,tb3,tb4,tb5,tb6;

    tb1=gbx+200; tb2=tb1+60; tb3=tb2+60; tb4=tb3+40; tb5=tb4+40; tb6=tb5+40;
    for (i=0; i<18; i++) {
      if (i>=maxWalker) break;
      gOff.drawString("A(meV)",tb1,gby+20);
      gOff.drawString("x(meV)",tb2,gby+20);
      gOff.drawString("p",tb3,gby+20);
      gOff.drawString("q",tb4,gby+20);
      gOff.drawString("r0(A)",tb5,gby+20);
      gOff.drawString("err",tb6,gby+20);
      gOff.drawString(""+(int)(xx[i][0]*10000.0+0.5)/10.0+"",tb1,gby+40+15*i);
      gOff.drawString(""+(int)(xx[i][1]*10000.0+0.5)/10.0+"",tb2,gby+40+15*i);
      gOff.drawString(""+(int)(xx[i][2]*100.0+0.5)/100.0+"",tb3,gby+40+15*i);
      gOff.drawString(""+(int)(xx[i][3]*100.0+0.5)/100.0+"",tb4,gby+40+15*i);
      gOff.drawString(""+(int)(xx[i][4]*100.0+0.5)/100.0+"",tb5,gby+40+15*i);
      gOff.drawString(""+(int)(ppp[i]*1000.0+0.5)/1000.0+"",tb6,gby+40+15*i);
    }
  }

/*---------------------- potential plot ---------------*/

  void drawGrid(int gbx, int gby) {
    int i,j,nnx,nny,nnDiv,nny0;

    nnx = 500; nny = 300; nnDiv = 50; 
    nny0 = 100;

    gOff.setColor(Color.lightGray);
    for (i=0; i<nnx; i+=nnDiv) {
      gOff.drawLine(gbx+i,gby, gbx+i, gby+nny);
    }
    for (j=0; j<nny; j+=nnDiv) {
      gOff.drawLine(gbx,gby+j, gbx+nnx, gby+j);
    }
    gOff.setColor(Color.black);
    gOff.drawRect(gbx,gby,nnx,nny);

    for (i=0;i<=10;i++) {
      gOff.drawString(""+i+"",gbx-5+nnDiv*i,gby+nny+20);
    }
	gOff.drawString("A",gbx-5+nnDiv*10+20,gby+nny+20);

    gOff.drawString("-2",gbx-25,gby+300+5);
    gOff.drawString("0",gbx-15,gby+200+5);
    gOff.drawString("2",gbx-15,gby+100+5);
	gOff.drawString("4eV",gbx-30,gby+5);
  }

  void drawLJ(int gbx, int gby, double mag,double col) {
    int i,j,nnx,nny,gy0;
    int f0,f1;
    double ra;

    nnx = 500; nny = 300;
    gy0 = gby+200;


    gOff.setColor(Color.getHSBColor((float)col,0.8f,0.9f));
    for (i=50; i<nnx-1; i++) {
	  ra = i*0.02;
	  f0 = (int)(50.0*mag*potentialLJ(MAT1,ra));
      if (f0>250) continue;
	  if (f0<-150) continue;
	  f1 = (int)(50.0*mag*potentialLJ(MAT1,ra+0.02));
      if (f1>250) continue;
	  if (f1<-150) continue;
      gOff.drawLine(gbx+i,gy0-f0, gbx+i+1, gy0-f1);
    }
  }

  void drawTrial(int gbx, int gby,double mag, double col) {
    int i;
    double dd,r0,aa;

    for (i=0; i<maxWalker; i++) {
      drawTrial(gbx,gby,xx[i],mag,col+(0.1*i)/maxWalker);
    }
  }

  void drawTrial(int gbx, int gby, double xi[], double mag, double col) {
    int i,j,nnx,nny,gy0;
    int f0,f1;
    double ra;

    nnx = 500; nny = 300;
    gy0 = gby+200;

    gOff.setColor(Color.getHSBColor((float)col,0.8f,0.9f));
    for (i=50; i<nnx-1; i++) {
      ra = i*0.02;
      f0 = (int)(50.0*mag*trialMorse(xi,ra));
      if (f0>250) continue;
      if (f0<-150) continue;
      f1 = (int)(50.0*mag*trialMorse(xi,ra+0.02));
      if (f1>250) continue;
      if (f1<-150) continue;
      gOff.drawLine(gbx+i,gy0-f0, gbx+i+1, gy0-f1);
    }
  }

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

  void setInitialCondition() {
    int i,j;

    t = 0.0;

/* Morse */
    xxMin[0] = 0.0001; xxMax[0] = 0.1; sss[0] = ss0; // dd
    xxMin[1] = 2.0; xxMax[1] = 6.0;	sss[1] = ss0; // r0
    xxMin[2] = 0.5; xxMax[2] = 3.0;	sss[2] = ss0; // aa
/* */

/* TBP 
    xxMin[0] = 0.0001; xxMax[0] = 0.1; sss[0] = ss0; // aa
    xxMin[1] = 0.0001; xxMax[1] = 4.0;	sss[1] = ss0; // xxx
    xxMin[2] = 6.0; xxMax[2] =15.0;	sss[2] = ss0; // pp
    xxMin[3] = 0.8; xxMax[3] = 4.0;	sss[3] = ss0; // qq
    xxMin[4] = 2.0; xxMax[4] = 6.0;	sss[4] = ss0; // r0
  */

    for (i=0; i<maxWalker; i++) {
      for (j=0; j<varMax; j++) {
        xx[i][j] = xxMin[j]+(xxMax[j]-xxMin[j])*Math.random();
      }
    }
  }

  void setsss(double ss) {
    int i;

    for (i=0; i<varMax; i++) {
      sss[i] = ss;
    }
  }

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

  void timeEvolution() {
    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
      started = 1;
    }

    if (started==1) {
      timeStep();
    }
  }

  void timeStep() {
    int i;

    t = t + dt;

    for (i=0; i<maxWalker; i++) {
      ppp[i] = MetropolisSampling(i, betaInv);
    }
  }

  double MetropolisSampling(int ii, double kT) {
    int j;
    double xij,e0,e1,de;

    e0 = potential(xx[ii]);
    for (j=0; j<varMax; j++) {
      xij = xx[ii][j] + sss[j]*(xxMax[j]-xxMin[j])*(Math.random()-0.5);
      if (xij<xxMin[j]) xij = xxMin[j];
      if (xij>xxMax[j]) xij = xxMax[j];
      xxTrial[j] = xij;
    }
    e1 = potential(xxTrial);
    de = e1 - e0;
    if (de<=0.0) {
      set(xx[ii],xxTrial);
      e0 = e1;
    } else {
      if (Math.exp(-de/kT)>Math.random() ) {
        set(xx[ii],xxTrial);
        e0 = e1;
      }
    }
    return ( e0 );
  }

  void set(double x0[], double x1[] ) {
    int j;

    for (j=0; j<varMax; j++) {
      x0[j] = x1[j];
    }
  }

  double potential(double xi[]) {
    int i;
    double sgmInA,pMin,pLJ,s,ra,d,w;

    sgmInA = LJdata[MAT1][3]/1.0e-10;
    pMin = potentialLJ(MAT1, 1.12246*sgmInA);

    s = 0.0;
    for (i=10; i<100; i++) {
      ra = i*0.1;
      pLJ = potentialLJ(MAT1, ra);
      d = pLJ - trialMorse(xi,ra);
      w = 0.01*i/(pLJ-2.0*pMin);
      s += d*d*w*w*w;
    }
    return ( s );
  }


  double potentialLJ(int knd, double ra) {
    double epsIneV,sgmInA,svr,svr3,svr6;

    epsIneV = LJdata[knd][2]/1.602e-19;
    sgmInA = LJdata[knd][3]/1.0e-10;
    svr = sgmInA/ra;
    svr3 = svr*svr*svr;
    svr6 = svr3*svr3;
    return( 4.0*epsIneV*(svr6*(svr6-1.0)) );
  }

  double trialMorse(double xi[], double ra) {
    double dd,r0,aa,p;
	
    dd = xi[0]; r0 = xi[1]; aa = xi[2];
    p = Math.exp(-aa*(ra-r0));

    return( dd*p*(p-2.0) );
  }

/*
  double trialTBP(double xi[], double ra) {
    double aa,xxx,pp,qq,r0, uu,ff,gg,b;

    gg = 3.5;
    aa = xi[0]; xxx = xi[1]; pp = xi[2]; qq = xi[3]; r0 = xi[4];

    uu = aa*Math.exp(-pp*(ra/r0-1.0));
    ff = xxx*xxx*Math.exp(-2.0*qq*(ra/r0-1.0));
    b = ff/gg;

    return( 2.0*uu - b + b*b/(4.0*gg) );
  }
*/

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

}