// ComPhys File: MDSimulation.java
// Chapter 8: A Molecular Dynamics Program

import java.awt.*;
import java.awt.event.*;
import java.applet.*;
import java.io.*;
 
public class MDSimulation extends Applet implements ActionListener, ItemListener, Runnable {

    double[] x = new double[36];
    double[] y = new double[36];
    double[] vx = new double[36];
    double[] vy = new double[36];
    double[] ax = new double[36];
    double[] ay = new double[36];

    int N = 16;
    double Lx = 6.0;
    double Ly = 6.0;
    double dt;
    double dt2;
    double E;
    int ncum;

    double[] x_initial = {1.09,3.12,0.08,0.54,2.52,3.03,4.25,0.89,
                          2.76,3.14,0.23,1.91,4.77,5.10,4.97,3.90};
    double[] y_initial = {0.98,5.25,2.38,4.08,4.39,2.94,3.01,3.11,
                          0.31,1.91,5.71,2.46,0.96,4.63,5.88,0.20};
    double[] vx_initial = {-0.33,0.12,-0.08,-1.94,0.75,1.70,0.84,-1.04,
                           1.64,0.38,-1.58,-1.55,-0.23,-0.31,1.18,0.46};
    double[] vy_initial = {-0.78,-1.19,-0.10,-0.56,0.34,-1.08,0.47,0.06,
                           1.36,-1.24,0.55,-0.16,-0.83,0.65,1.48,-0.51};

    // variables set by initial()
    double t;
    double ke;
    double kecum;
    double pecum;
    double vcum;
    double area;

    void initial () {

        dt = 0.01;
        dt2 = dt * dt;

        for (int i = 0; i < N; ++i) {
            x[i] = x_initial[i];
            y[i] = y_initial[i];
            vx[i] = vx_initial[i];
            vy[i] = vy_initial[i];
        }

        ke = 0;
        for (int i = 0; i < N; ++i)
            ke += vx[i] * vx[i] + vy[i] * vy[i];
        ke *= 0.5;
        area = Lx * Ly;
        t = 0;
        kecum = 0;
        pecum = 0;
        vcum = 0;

        accel();
        E = ke + pe;
        ncum = 0;

        show_output();
        show_positions();
        
    }

    // variables set by accel()
    double pe;
    double virial;

    // variables used by accel() and set by force()
    double fx;
    double fy;
    double pot;

    void accel () {

        for (int i = 0; i < N; ++i) {
            ax[i] = 0;
            ay[i] = 0;
        }

        pe = 0;
        virial = 0;

        for (int i = 0; i < N - 1; ++i) {
            for (int j = i + 1; j < N; ++j) {

                double dx = separation(x[i] - x[j], Lx);
                double dy = separation(y[i] - y[j], Ly);

                force(dx, dy);

                ax[i] += fx;
                ay[i] += fy;
                ax[j] -= fx;
                ay[j] -= fy;
                pe += pot;
                virial += dx * ax[i] * dy * ay[i];
                
            }
        }

    }

    double separation (double ds, double L) {

        if (ds > 0.5 * L)
            ds -= L;
        if (ds < - 0.5 * L)
            ds += L;

        return ds;

    }

    void force (double dx, double dy) {

        double r2 = dx * dx + dy * dy;
        double rm2 = 1 / r2;
        double rm6 = rm2 * rm2 * rm2;
        double f_over_r = 24 * rm6 * (2 * rm6 - 1) * rm2;
        fx = f_over_r * dx;
        fy = f_over_r * dy;
        pot = 4 * (rm6 * rm6 - rm6);

    }

    void Verlet () {

        for (int i = 0; i < N; ++i) {

            double xnew = x[i] + vx[i] * dt + 0.5 * ax[i] * dt2;
            double ynew = y[i] + vy[i] * dt + 0.5 * ay[i] * dt2;

            vx[i] += 0.5 * ax[i] * dt;
            vy[i] += 0.5 * ay[i] * dt;

            x[i] = pbc(xnew, Lx);
            y[i] = pbc(ynew, Ly);

        }

        accel();

        ke = 0;
        for (int i = 0; i < N; ++i) {

            vx[i] += 0.5 * ax[i] * dt;
            vy[i] += 0.5 * ay[i] * dt;
            ke += vx[i] * vx[i] + vy[i] * vy[i];

        }

        ke *= 0.5;
        t += dt;

    }

    double pbc (double pos, double L) {

        if (pos < 0)
            pos = pos + L;
        if (pos > L)
            pos = pos - L;

        return pos;

    }

    Label numberLabel, densityLabel, ncumLabel, tLabel, ELabel, TLabel, pLabel;

    String format (double d, int decimalPlaces) {

        for (int i = 0; i < decimalPlaces; ++i)
            d *= 10;
        d = Math.round(d);
        for (int i = 0; i < decimalPlaces; ++i)
            d /= 10;
        String str = new String("" + d);
        int len = Math.min(str.length(),
                           str.indexOf('.') + decimalPlaces + 1); 
        return str.substring(0, len);
                                         
    }

    double T;

    void show_output () {

        ++ncum;
        ncumLabel.setText("" + ncum);
        tLabel.setText(format(t, 2));
        E = ke + pe;
        ELabel.setText(format(E, 3));
        kecum += ke;
        vcum += virial;
        double mean_ke = kecum / ncum;
        T = mean_ke / N;
        TLabel.setText(format(T, 3));
        double p = mean_ke + 0.5 * vcum / ncum;
        p /= area;
        pLabel.setText(format(p, 3));
        
    }

    void show_positions () {

        mdPicture.clearPicture = true;
        mdPicture.repaint();

    }

    void save_config () {

        Frame frame = new Frame("Save Configuration");
        FileDialog fileDialog = new FileDialog(frame);
        fileDialog.show();
        String fileName = fileDialog.getFile();
        frame.dispose();

        if (fileName == null)
            return;

        try {

            PrintWriter pw = new PrintWriter(new FileWriter(fileName));
            pw.println(N);
            pw.println(Lx);
            pw.println(Ly);

            for (int n = 0; n < N; ++n)
                pw.println(x[n] + "\t" + y[n] + "\t" + vx[n] + "\t" + vy[n]);

            pw.close();

        } catch (IOException ioe) {

            System.err.println("Error writing file: " + fileName);

        }

    }

    String str;     // set by load_config() and reset by next_double()

    double next_double () throws NumberFormatException {

        // extract next double from str and convert to double

        str = str.trim();
        int len = str.length();
        char c;
        int ind = 0;
        while ( (c = str.charAt(ind)) != ' '
                && c != '\t' && ind < len - 1 )
            ++ind;
        String s = str.substring(0, ind);
        str = str.substring(ind);
        double d = new Double(s).doubleValue();
        return d;

    }

    void load_config () {

        Frame frame = new Frame("Load Configuration");
        FileDialog fileDialog = new FileDialog(frame);
        fileDialog.show();
        String fileName = fileDialog.getFile();
        frame.dispose();

        if (fileName == null)
            return;

        try {

            BufferedReader br = new BufferedReader(new FileReader(fileName));
            str = br.readLine().trim();
            N = Integer.parseInt(str);
            str = br.readLine().trim();
            Lx = new Double(str).doubleValue();
            str = br.readLine().trim();
            Ly = new Double(str).doubleValue();

            if (N > 0) {

                x_initial = new double[N];
                y_initial = new double[N];
                vx_initial = new double[N];
                vy_initial = new double[N];

            }

            for (int n = 0; n < N; ++n) {

                str = br.readLine().trim();

                x_initial[n] = next_double();
                y_initial[n] = next_double();
                vx_initial[n] = next_double();
                vy_initial[n] = next_double();

            }

            br.close();

        } catch (IOException ioe) {

            System.err.println("Error reading file: " + fileName);

        }

        numberLabel.setText("N Particles = " + N);
        densityLabel.setText("Density = " + format(N / Lx / Ly, 3));
        show_positions();

    }

    MDPicture mdPicture;        // to draw molecule paths
                                // class MDPicture defined at end of file
    Button runButton;           // to start and stop motion
    Button resetButton;         // to reset initial conditions
    Button saveButton;          // to save configuration in file
    Button loadButton;          // to load configuration from file
    Checkbox checkbox;          // to show or hide molecules
    int stepsToSkip;            // time steps to skip when drawing
    TextField skipTextField;    // to enter stepsToSkip
    TextField newTTextField;    // to enter new T
    Thread runThread;           // animation in a separate thread
    boolean running;            // animation running or not

    // Applet initialization
    public void init () {

        add(mdPicture = new MDPicture(this, 400, 400));

        Panel p = new Panel();
        p.setLayout(new GridLayout(11, 2));
                    
        p.add(numberLabel = new Label("N Particles = " + N));
        p.add(densityLabel = new Label("Density = " + format(N / Lx / Ly, 3)));

        p.add(new Label("Time Steps = "));
        p.add(ncumLabel = new Label("      "));

        p.add(new Label("Time t = "));
        p.add(tLabel = new Label("      "));

        p.add(new Label("Energy E = "));
        p.add(ELabel = new Label("      "));

        p.add(new Label("Temperature <T> = "));
        p.add(TLabel = new Label("      "));

        p.add(new Label("Enter new T:"));
        p.add(newTTextField = new TextField(""));
        newTTextField.addActionListener(this);

        p.add(new Label("Pressure <P> = "));
        p.add(pLabel = new Label("      "));

        p.add(new Label("Click in Pic to Clear"));
              
        p.add(checkbox = new Checkbox("Show Molecules"));
        checkbox.addItemListener(this);

        p.add(new Label("Enter Steps to skip:"));
        
        p.add(skipTextField = new TextField("0"));
        skipTextField.addActionListener(this);
        
        runButton = new Button("Start");
        p.add(runButton);
        runButton.addActionListener(this);

        p.add(resetButton = new Button("Reset"));
        resetButton.addActionListener(this);

        p.add(saveButton = new Button("Save Config"));
        saveButton.addActionListener(this);

        p.add(loadButton = new Button("Load Config"));
        loadButton.addActionListener(this);

        add(p);

        initial();              // simulation initialization

    }

    // must define this to implement ActionListener interface
    public void actionPerformed (ActionEvent event) {

        String arg = event.getActionCommand();

        if (arg.equals("Reset")) {

            running = false;
            try { Thread.sleep(100); } catch (InterruptedException e) { }
            initial();

        } else if (arg.equals("Save Config")) {

            running = false;
            try { Thread.sleep(100); } catch (InterruptedException e) { }
            save_config();

        } else if (arg.equals("Load Config")) {

            running = false;
            try { Thread.sleep(100); } catch (InterruptedException e) { }
            load_config();
            initial();

        } else if (arg.equals("Start") || arg.equals("Stop")) {

            if (!running) {
                
                running = true;
                runButton.setLabel("Stop");
                runThread = new Thread(this);
                runThread.start();

            } else {

                running = false;

            }

        } else {

            // check text fields for changed values
            
            if ( event.getSource() == skipTextField ) {

                int skip = stepsToSkip;
                try {
                    
                    skip = Integer.parseInt(arg);
                    stepsToSkip = skip;
            
                } catch (NumberFormatException nfe) {
                    
                    skipTextField.setText("" + stepsToSkip);

                }
            
            } else if ( event.getSource() == newTTextField ) {

                double newT = T;
                try {

                    newT = new Double(arg).doubleValue();

                } catch (NumberFormatException nfe) {

                    newTTextField.setText(format(T, 3));
                    
                }

                double scale = Math.sqrt(newT / T);
                for (int n = 0; n < N; ++n) {

                    vx[n] *= scale;
                    vy[n] *= scale;

                }

            }
                
        }

        if (!running)
            runButton.setLabel("Start");

    }

    // must define this to implement ItemListener interface
    public void itemStateChanged (ItemEvent evt) {

        mdPicture.showMolecules = checkbox.getState();
        show_positions();

    }

    // must define this to implement Runnable interface
    public void run () {

        while (running) {

            runThread.setPriority(Thread.MIN_PRIORITY);

            // skip output for stepsToSkip to make program run faster
            for (int step = 0; step < stepsToSkip; ++step) {
                Verlet();
                ++ncum;
                E += ke + pe;
                kecum += ke;
                vcum += virial;
            }
            Verlet();
            mdPicture.repaint();
            show_output();

            try {
                Thread.sleep(30);
            } catch (InterruptedException e) { }

        }

    }

    public void paint (Graphics g) {

        mdPicture.repaint();

    }

    public static void main (String[] args) {

        MDSimulation md = new MDSimulation();
 
        Frame aFrame = new Frame("Molecular Dynamics Simulation");
 
        // anonymous inner class to destroy window
        aFrame.addWindowListener(new WindowAdapter () {
            public void windowClosing(WindowEvent e) {
                System.exit(0);
            }
        });
        
        aFrame.add(md);
        md.init();
        aFrame.setSize(700, 450);
        aFrame.setLocation(100, 100);
        aFrame.setVisible(true);
 
    }

}

class MDPicture extends Canvas implements MouseListener {

    int xSize, ySize;
    Color bgColor = new Color(255, 250, 240);
    MDSimulation md;            // simulation to which this picture belongs

    MDPicture (MDSimulation mds, int dx, int dy) {

        md = mds;
        xSize = dx;
        ySize = dy;
        setSize(xSize, ySize);
        setBackground(bgColor);
        addMouseListener(this);

    }

    boolean clearPicture = true;

    // must define the following fucntions to implement MouseListener
    public void mouseClicked (MouseEvent me) {

        clearPicture = true;
        repaint();

    }

    public void mouseEntered (MouseEvent me) { }
    public void mouseExited (MouseEvent me) { }
    public void mousePressed (MouseEvent me) { }
    public void mouseReleased (MouseEvent me) { }

    public void paint (Graphics g) {

        update(g);

    }

    boolean showMolecules;      // draw circles at molecule positions
    Image offScreen;            // for double buffering to show molecules

    public void update (Graphics g) {

        double xScale = xSize / md.Lx;
        double yScale = ySize / md.Ly;
        int scale = (int) Math.round(Math.min(xScale, yScale));
        int xMax = (int) Math.round(md.Lx * scale);
        int yMax = (int) Math.round(md.Ly * scale);

        if (offScreen == null)
            offScreen = createImage(xSize, ySize);
        
        Graphics gr = offScreen.getGraphics();

        if (showMolecules) {

            gr.clearRect(0, 0, xSize, ySize);
            int r = scale / 2;
            gr.setColor(Color.red);
            for (int n = 0; n < md.N; ++n) {

                int ix = (int) (md.x[n] * scale);
                int iy = (int) (md.y[n] * scale);
                gr.fillOval(ix - r, iy - r, 2 * r, 2 * r);

                // draw periodic images if partly visible on picture
                if (ix < r)
                    gr.fillOval(ix + xMax - r, iy - r, 2 * r, 2 * r);
                if (iy < r)
                    gr.fillOval(ix - r, iy + yMax - r, 2 * r, 2 * r);
                if (ix > xMax - r)
                    gr.fillOval(ix - xMax - r, iy - r, 2 * r, 2 * r);
                if (iy > yMax - r)
                    gr.fillOval(ix - r, iy - yMax - r, 2 * r, 2 * r);
                if (ix < r && iy < r)
                    gr.fillOval(ix + xMax - r, iy + yMax - r, 2 * r, 2 * r);
                if (ix < r && iy > yMax - r)
                    gr.fillOval(ix + xMax - r, iy - yMax - r, 2 * r, 2 * r);
                if (ix > xMax - r && iy < r)
                    gr.fillOval(ix - xMax - r, iy + yMax - r, 2 * r, 2 * r);
                if (ix > xMax - r && iy > yMax - r)
                    gr.fillOval(ix - xMax - r, iy - yMax - r, 2 * r, 2 * r);

            }

        } else {

            if (clearPicture) {

                gr.clearRect(0, 0, xSize, ySize);

                gr.setColor(Color.lightGray);
                for (int ix = scale; ix < xMax; ix += scale)
                    gr.drawLine(ix, 0, ix, yMax);
                for (int iy = scale; iy < yMax; iy += scale)
                    gr.drawLine(0, iy, xMax, iy);

                clearPicture = false;

            }

            gr.setColor(Color.red);
            for (int n = 0; n < md.N; ++n) {

                int ix = (int) (md.x[n] * scale);
                int iy = (int) (md.y[n] * scale);
                gr.drawLine(ix, iy, ix, iy);

            }

        }

        if (xMax < xSize)
            gr.clearRect(xMax, 0, xSize, ySize);
        if (yMax < ySize)
            gr.clearRect(0, yMax, xSize, ySize);
        gr.setColor(Color.black);
        gr.drawRect(0, 0, xMax - 1, yMax - 1);

        g.drawImage(offScreen, 0, 0, null);

    }

}