import java.awt.*;
import java.awt.event.*;
import java.applet.Applet;
import java.util.*;
import java.lang.Math;

class Atom {
  public double x;
  public double y;
  public double vx;
  public double vy;
  public double fx;
  public double fy;
  public Atom(double x, double y, double vx, double vy) {
    this.x = x;
    this.y = y;
    this.vx = vx;
    this.vy = vy;
    this.fx = 0.0;
    this.fy = 0.0;
  }
}

public class MDDemo extends Applet {
  MDPanel panel;
  Simulation s;

  public void init() {
    setLayout(new GridLayout(1,0));
    s = new Simulation(150.0,0.069);
    panel = new MDPanel(this,s);
    s.setPanel(panel);
    add(panel);
  }

  //The following functions should be called by events from the panel.
  // Before they are called, the simulation thread should be suspended.

  public void resetSim() {
    Enumeration e;
    Atom a;
    double x,y,vx,vy;
    int atomd = (int)Math.sqrt((double)s.N_ATOMS);
    double cellEdge = s.regionEdge/atomd;

    e = s.atoms.elements();

    for(int i = 0; i < atomd; i++) {
      for(int j = 0 ; j < atomd ; j++) {
	x = (((double)i + 0.5) * cellEdge) - s.hregionEdge;
	y = (((double)j + 0.5) * cellEdge) - s.hregionEdge;
	vx = Math.sqrt( (s.KB * s.applTemp) / s.mass ) * s.rand.nextGaussian();
	vy = Math.sqrt( (s.KB * s.applTemp) / s.mass ) * s.rand.nextGaussian();
	a = (Atom)e.nextElement();
	a.x = x;
	a.y = y;
	a.vx = vx;
	a.vy = vy;
	a.fx = 0.0;
	a.fy = 0.0;
      }
    }
  }

  public void changeDensity(double new_dens) {
    Enumeration e;
    Atom a;
    double factor;
    double oldedge;

    oldedge = s.regionEdge;
    e = s.atoms.elements();
    s.regionEdge = Math.sqrt((double)s.N_ATOMS / new_dens);
    s.hregionEdge = 0.5 * s.regionEdge;
    s.density = new_dens;

    factor = s.regionEdge / oldedge;
    while(e.hasMoreElements()) {
      a = (Atom)e.nextElement();
      a.x = a.x * factor;
      a.y = a.y * factor;
    }
  }
}

class Simulation extends Thread {
  static final int N_ATOMS = 100;  //must be a perfect square
  MDPanel panel = null;
  LinkedList atoms;
  Random rand;
  double applTemp,temp;
  double density;
  double regionEdge;
  double hregionEdge;
  double pe;
  double ke;
  double te;
  int ts = 0;
  static final double mass = 39.948;
  static final double epsilon = 398.4;
  static final double sigma = 3.41;
  static final double cutoff = 6.0;
  static final double AVOGAD = 6.0221367e23;
  static final double MASSUNIT = 1.6605402e-27;
  static final double LENGTHUNIT = 1.0e-10;
  static final double TIMEUNIT = 1.0e-12;
  static final double ENERGYUNIT = (MASSUNIT * (LENGTHUNIT / TIMEUNIT) *
				    (LENGTHUNIT / TIMEUNIT));
  static final double KB = 1.380658e-23 / ENERGYUNIT;
  static final double TIMESTEP = 0.01;

  public Simulation(double t, double d) {
    Atom a;
    double x,y,cellEdge;
    double vx,vy;
    int atomd;
    
    atomd = (int)Math.sqrt((double)N_ATOMS);
    this.density = d;
    this.applTemp = t;
    regionEdge = Math.sqrt((double)N_ATOMS/this.density);
    hregionEdge = 0.5 * regionEdge;
    cellEdge = regionEdge/atomd;
    rand = new Random();
    atoms = new LinkedList();
    for(int i = 0; i < atomd; i++) {
      for(int j = 0 ; j < atomd ; j++) {
	x = (((double)i + 0.5) * cellEdge) - hregionEdge;
	y = (((double)j + 0.5) * cellEdge) - hregionEdge;
	vx = Math.sqrt( (KB * applTemp) / mass ) * rand.nextGaussian();
	vy = Math.sqrt( (KB * applTemp) / mass ) * rand.nextGaussian();
	a = new Atom(x,y,vx,vy);
	atoms.append(a);
      }
    }
  }

  public void run() {
    calc_forces();
    measurements();
    updatePanel();
    while(true) {
      integrate_step1();
      zero_forces();
      calc_forces();
      integrate_step2();
      boundary_cond();
      measurements();
      updatePanel();
      panel.redraw();
      try {sleep(2);} catch (InterruptedException e) {}
    } 
  }

  public void updatePanel() {
    String str;
    int n;

    str = String.valueOf(temp);
    if(str.length() < 5) n = str.length();
    else n = 5;
    panel.tempField.setText(str.substring(0,n));
    str = String.valueOf(te);
    if(str.length() < 8) n = str.length();
    else n = 8;
    panel.energyField.setText(str.substring(0,n));
    str = String.valueOf(ts);
    panel.timeField.setText(str);
    ts++;
  }
  
  public void calc_forces() {
    Enumeration e1,e2;
    Atom a1,a2;
    double rij2,rij;
    double rijx,rijy;
    double ri2,ri12,ri6;
    double cutoff2 = cutoff * cutoff;
    double force;
    
    e1 = atoms.elements();
    e2 = atoms.elements();
    pe = 0.0;

    while(e1.hasMoreElements()) {
      a1 = (Atom)e1.nextElement();
      e2 = atoms.elements();
      while(e2.hasMoreElements()) {
	a2 = (Atom)e2.nextElement();
	if(a2 != a1) {
	  rijx = a1.x - a2.x;
	  if(rijx > hregionEdge) rijx = rijx - regionEdge;
	  if(rijx < -hregionEdge) rijx = rijx + regionEdge;
	  rijy = a1.y - a2.y;
	  if(rijy > hregionEdge) rijy = rijy - regionEdge;
	  if(rijy < -hregionEdge) rijy = rijy + regionEdge;
	  rij2 = (rijx * rijx) + (rijy * rijy);
	  if(rij2 < cutoff2) {
	    ri2 = 1.0 / rij2;
	    ri6 = sigma * sigma * ri2;
	    ri6 = ri6 * ri6 * ri6;
	    ri12 = ri6 * ri6;
	    pe += epsilon * (ri12 - ri6);
	    force = ri2 * 6.0 * epsilon * ((2.0 * ri12) - ri6);
	    a1.fx += force * rijx;
	    a1.fy += force * rijy;
	  }
	}
      }
    }
    pe = 0.5 * pe;
  }

  public void zero_forces() {
    Enumeration e;
    Atom a;

    e = atoms.elements();
    while(e.hasMoreElements()) {
      a = (Atom)e.nextElement();
      a.fx = 0.0;
      a.fy = 0.0;
    }
  }

  public void integrate_step1() {
    Enumeration e;
    Atom a;
    double dt2 = TIMESTEP * TIMESTEP;

    e = atoms.elements();
    while(e.hasMoreElements()) {
      a = (Atom)e.nextElement();
      a.x = a.x + (TIMESTEP * a.vx) + (0.5 * dt2 * (a.fx / mass));
      a.y = a.y + (TIMESTEP * a.vy) + (0.5 * dt2 * (a.fy / mass));
      a.vx = a.vx + (0.5 * TIMESTEP * (a.fx / mass));
      a.vy = a.vy + (0.5 * TIMESTEP * (a.fy / mass));
    }
  }

  public void integrate_step2() {
    Enumeration e;
    Atom a;

    e = atoms.elements();
    while(e.hasMoreElements()) {
      a = (Atom)e.nextElement();
      a.vx = a.vx + (0.5 * TIMESTEP * (a.fx / mass));
      a.vy = a.vy + (0.5 * TIMESTEP * (a.fy / mass));
    }
  }

  public void boundary_cond() {
    Enumeration e;
    Atom a;

    e = atoms.elements();
    while(e.hasMoreElements()) {
      a = (Atom)e.nextElement();
      if(a.x > hregionEdge) a.x = a.x - regionEdge;
      if(a.x <= -hregionEdge) a.x = a.x + regionEdge;
      if(a.y > hregionEdge) a.y = a.y - regionEdge;
      if(a.y <= -hregionEdge) a.y = a.y + regionEdge;
    }
  }

  public void measurements() {
    Enumeration e;
    Atom a;
    double vsqr = 0.0;

    ke = 0.0;
    te = 0.0;
    e = atoms.elements();
    while(e.hasMoreElements()) {
      a = (Atom)e.nextElement();
      vsqr += (a.vx * a.vx) + (a.vy * a.vy);
    }
    ke = 0.5 * mass * vsqr;
    te = ke + pe;
    temp = (mass * vsqr) / (2.0 * N_ATOMS * KB);
  }

  public void setPanel(MDPanel p) {
    this.panel = p;
  }
}

class MDPanel extends Panel {
  Label edgeLabel;
  TextField timeField;
  TextField tempField;
  TextField applTempField;
  TextField energyField;
  TextField densField;
  FramedArea regionCanvas;
  Scrollbar tempSlider;
  Scrollbar densSlider;
  MDDemo controller;
  Button startButton;
  Button stopButton;
  Button resetButton;
  Simulation s;
  private boolean threadRunning = false;
  private boolean threadSuspended = false;

  public MDPanel(MDDemo mycontroller, Simulation s) {
    this.s = s;

    GridBagLayout gridBag = new GridBagLayout();
    GridBagConstraints c = new GridBagConstraints();
    setLayout(gridBag);

    //Save arguments in instance variables
    controller = mycontroller;

    //Default constraints
    c.fill = GridBagConstraints.HORIZONTAL;

    Label label1 = new Label("Java Molecular Dynamics",Label.CENTER);
    c.gridwidth = GridBagConstraints.REMAINDER;
    gridBag.setConstraints(label1,c);
    add(label1);
    
    regionCanvas = new FramedArea(this,s);
    c.fill = GridBagConstraints.BOTH;
    c.gridwidth = GridBagConstraints.REMAINDER;
    c.weighty = 1.0;
    regionCanvas.setBackground(Color.white);
    gridBag.setConstraints(regionCanvas, c);
    add(regionCanvas);

    String str = String.valueOf(s.regionEdge);
    int n;
    if(str.length() > 6) n = 6;
    else n = str.length();
    edgeLabel = new Label(str.substring(0,n) + " Angstroms",Label.CENTER);
    c.fill = GridBagConstraints.HORIZONTAL;
    c.weightx = 1.0;
    c.weighty = 0.0;
    gridBag.setConstraints(edgeLabel,c);
    add(edgeLabel);

    Label label2 = new Label("Timestep:",Label.LEFT);
    c.fill = GridBagConstraints.HORIZONTAL;
    c.weightx = 1.0;
    c.weighty = 0.0;
    c.gridwidth = 1;
    c.gridheight = 1;
    gridBag.setConstraints(label2,c);
    add(label2);

    timeField = new TextField("0",5);
    c.weightx = 1.0;
    gridBag.setConstraints(timeField,c);
    add(timeField);
    timeField.setEditable(false);

    Label label3 = new Label("Temperature:",Label.LEFT);
    gridBag.setConstraints(label3,c);
    add(label3);
    
    tempField = new TextField("0",6);
    c.gridwidth = 1;
    gridBag.setConstraints(tempField,c);
    add(tempField);
    tempField.setEditable(false);

    Label label4 = new Label("Energy:",Label.LEFT);
    gridBag.setConstraints(label4,c);
    add(label4);
    
    energyField = new TextField("0",8);
    c.gridwidth = GridBagConstraints.REMAINDER;
    gridBag.setConstraints(energyField,c);
    add(energyField);
    energyField.setEditable(false);

    Label label5 = new Label("Applied Temp:",Label.LEFT);
    c.gridwidth = 1;
    gridBag.setConstraints(label5,c);
    add(label5);

    tempSlider = new Scrollbar(Scrollbar.HORIZONTAL,150,50,1,150);
    c.gridwidth = 4;
    gridBag.setConstraints(tempSlider,c);
    add(tempSlider);

    applTempField = new TextField("150",6);
    c.gridwidth = GridBagConstraints.REMAINDER;
    c.fill = GridBagConstraints.NONE;
    gridBag.setConstraints(applTempField,c);
    add(applTempField);
    applTempField.setEditable(false);

    Label label6 = new Label("Density (atoms/A^2):",Label.LEFT);
    c.gridwidth = 1;
    gridBag.setConstraints(label6,c);
    add(label6);

    densSlider = new Scrollbar(Scrollbar.HORIZONTAL,69,150,10,150);
    c.fill = GridBagConstraints.HORIZONTAL;
    c.gridwidth = 4;
    gridBag.setConstraints(densSlider,c);
    add(densSlider);

    densField = new TextField("0.069",6);
    c.gridwidth = GridBagConstraints.REMAINDER;
    c.fill = GridBagConstraints.NONE;
    gridBag.setConstraints(densField,c);
    add(densField);
    densField.setEditable(false);

    startButton = new Button("Start");
    c.gridwidth = 2;
    gridBag.setConstraints(startButton,c);
    add(startButton);

    stopButton = new Button("Stop");
    gridBag.setConstraints(stopButton,c);
    add(stopButton);

    resetButton = new Button("Reset");
    c.gridwidth = GridBagConstraints.REMAINDER;
    gridBag.setConstraints(resetButton,c);
    add(resetButton);

  }
  
  public void redraw() {
    regionCanvas.redraw();
  }

  public boolean action(Event e, Object arg) {
    if(arg.equals("Stop")) {
      if(threadRunning) {
	controller.s.suspend();
	threadSuspended = true;
	threadRunning = false;
	return true;
      }
    } else if(arg.equals("Start")) {
      if(!threadRunning && !threadSuspended) {
	controller.s.start();
	threadRunning = true;
	threadSuspended = false;
	return true;
      }
      if(threadSuspended) {
	controller.s.resume();
	threadSuspended = false;
	threadRunning = true;
	return true;
      }
    } else if(arg.equals("Reset")) {
      if(threadRunning) {
	controller.s.suspend();
	threadSuspended = true;
	threadRunning = false;
      }
      controller.resetSim();
      redraw();
    }
    return true;
  }

  public boolean handleEvent(Event e) {
    if(e.id == Event.SCROLL_ABSOLUTE ||
       e.id == Event.SCROLL_LINE_DOWN ||
       e.id == Event.SCROLL_LINE_UP ||
       e.id == Event.SCROLL_PAGE_DOWN ||
       e.id == Event.SCROLL_PAGE_UP) {
      int n;
      if (e.target == (Object)densSlider) {
	double new_dens = getDensity();
	if(threadRunning) {
	  controller.s.suspend();
	}
	controller.changeDensity(new_dens);
	String s = String.valueOf(new_dens);
	if(s.length() > 6) n = 6;
	else n = s.length();
	densField.setText(s.substring(0,n));
	String str = String.valueOf(controller.s.regionEdge);
	if(str.length() > 6) n = 6;
	else n = str.length();
	this.edgeLabel.setText(str.substring(0,n) + " Angstroms");
	if(threadRunning) {
	  controller.s.resume();
	}
	return true;
      }
      if (e.target == (Object)tempSlider) {
	if(threadRunning) {
	  controller.s.suspend();
	}
	applTempField.setText(String.valueOf(tempSlider.getValue()));
	if(threadRunning) {
	  controller.s.resume();
	}
	return true;
      }
    }
    return super.handleEvent(e);
  }
  
  public double getDensity() {
    double dens;
    int d;

    d = densSlider.getValue();
    dens = (double)d * 0.001;
    return dens;
  }
}

class FramedArea extends Panel {
  MDPanel controller;
  SimulationRegion sr;

  public FramedArea(MDPanel controller, Simulation s) {
    super();
    this.controller = controller;

    setLayout(new GridLayout(1,0));
    sr = new SimulationRegion(controller,s);
    add(sr);
  }

  public Insets getInsets() {
    return new Insets(4,4,5,5);
  }

  public void redraw() {
    repaint();
    sr.repaint();
  }
    
  public void paint(Graphics g) {
    Dimension d = size();

    g.draw3DRect(0,0,d.width - 1, d.height - 1, true);
    g.draw3DRect(3,3,d.width - 7, d.height - 7, false);
  }
}

class SimulationRegion extends Canvas {
  Point point = null;
  MDPanel controller;
  Simulation s;

  public SimulationRegion(MDPanel controller,Simulation s) {
    super();
    this.controller = controller;
    this.s = s;
  }

  public void paint(Graphics g) {
    Atom a;
    Dimension d = size();
    double xfactor = d.width / s.regionEdge ;
    double yfactor = d.height / s.regionEdge ;
    int x,y;
    int radius = 3;

    g.setColor(Color.red);
    Enumeration e = s.atoms.elements();
    while(e.hasMoreElements()) {
      a = (Atom)e.nextElement();
      x = (int)((a.x + s.hregionEdge) * xfactor);
      y = (int)((a.y + s.hregionEdge) * yfactor);
      g.fillOval(x,y,2 * radius, 2 * radius);
    }
  }
}