/* Gibbs-ensemble simulation of a Lennard-Jones fluid Copyright 2011-2014 Cheng Zhang This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. A copy of the GNU General Public License can be found in . */ import java.applet.*; import java.awt.*; import java.awt.event.*; import java.awt.geom.*; import java.awt.image.*; import static java.lang.Math.*; import java.text.*; import javax.swing.*; import java.util.Hashtable; import java.util.Random; public class LJ3GibbsApp extends JApplet implements ActionListener { double rhoInit = 0.35; int NInit = 108; LJ3Gibbs mc1 = new LJ3Gibbs(rhoInit, NInit, NInit*2); LJ3Gibbs mc2 = new LJ3Gibbs(rhoInit, NInit, NInit*2); int delay = 100; int speed = 1000; // mc steps per frame int nstVMove = 10; double volAmp = 0.15; int nstNMove = 1; // frequency of particle transfers Timer timer; XYZCanvasMC canvas1, canvas2; // 3D drawing here JPanel cpnl, cpnl1, cpnl2, cpnl3, cpnl4, tpnl; JTextField tNum = new JTextField(" " + (mc1.N + mc2.N) + " "); JTextField tNum1 = new JTextField(" " + mc1.N + " "); JTextField tNum2 = new JTextField(" " + mc2.N + " "); JTextField tVol1 = new JTextField(" " + mc1.Vol + " "); JTextField tVol2 = new JTextField(" " + mc2.Vol + " "); JTextField tTemp = new JTextField(" 1.2 "); JTextField tRho = new JTextField(" " + rhoInit + " "); JTextField tSpeed = new JTextField(" " + speed + " "); JTextField tMCAmp = new JTextField(" " + mc1.mcAmp + " "); JTextField tVolAmp = new JTextField(" " + volAmp + " "); JTextField tVolFreq = new JTextField(" " + nstVMove + " "); JTextField tAccVol = new JTextField(" 0 "); JTextField tNumFreq = new JTextField(" " + nstNMove + " "); JTextField tAccNum = new JTextField(" 0 "); JButton bReset = new JButton("Reset"); JButton bRetime = new JButton("Clear data"); JToggleButton bStart = new JToggleButton("Start", false); JLabel lStatus = new JLabel("Status"); JTextField tRho1 = new JTextField(" 0"); JTextField tAvU1 = new JTextField(" 0"); JTextField tAvp1 = new JTextField(" 0"); JTextField tAcc1 = new JTextField(" 0"); JTextField tRho2 = new JTextField(" 0"); JTextField tAvU2 = new JTextField(" 0"); JTextField tAvp2 = new JTextField(" 0"); JTextField tAcc2 = new JTextField(" 0"); public static final long serialVersionUID = 1L; public void init() { Container box = getContentPane(); box.setLayout(new BorderLayout()); cpnl = new JPanel(); // create a panel for controls cpnl.setLayout(new GridLayout(5, 1)); box.add(cpnl, BorderLayout.SOUTH); cpnl1 = new JPanel(); cpnl.add(cpnl1); cpnl2 = new JPanel(); cpnl.add(cpnl2); cpnl3 = new JPanel(); cpnl3.setLayout(new GridLayout(1, 10)); cpnl.add(cpnl3); cpnl4 = new JPanel(); cpnl4.setLayout(new GridLayout(1, 10)); cpnl.add(cpnl4); // add controls cpnl1.add(bStart); bStart.addActionListener(this); cpnl1.add(bReset); bReset.addActionListener(this); cpnl1.add(bRetime); bRetime.addActionListener(this); cpnl1.add(new JLabel(" N:")); tNum.addActionListener(this); cpnl1.add(tNum); cpnl1.add(new JLabel(" T:")); tTemp.addActionListener(this); cpnl1.add(tTemp); cpnl1.add(new JLabel(" \u03c1 (overall):")); tRho.addActionListener(this); cpnl1.add(tRho); cpnl1.add(new JLabel(" Speed:")); tSpeed.addActionListener(this); cpnl1.add(tSpeed); cpnl1.add(new JLabel(" \u0394X (trial move size):")); tMCAmp.addActionListener(this); cpnl1.add(tMCAmp); cpnl2.add(new JLabel(" Volume exchange:")); cpnl2.add(new JLabel(" \u0394V:")); tVolAmp.addActionListener(this); cpnl2.add(tVolAmp); cpnl2.add(new JLabel(", ")); cpnl2.add(new JLabel(" every:")); tVolFreq.addActionListener(this); cpnl2.add(tVolFreq); cpnl2.add(new JLabel("steps, ")); cpnl2.add(new JLabel(" AR:")); tAccVol.setEditable(false); cpnl2.add(tAccVol); cpnl2.add(new JLabel(". ")); cpnl2.add(new JLabel(" Particle exchange:")); cpnl2.add(new JLabel(" every:")); tNumFreq.addActionListener(this); cpnl2.add(tNumFreq); cpnl2.add(new JLabel("steps, ")); cpnl2.add(new JLabel(" AR:")); tAccNum.setEditable(false); cpnl2.add(tAccNum); cpnl2.add(new JLabel(". ")); cpnl3.add(new JLabel(" < N\u2081 >:")); tNum1.setEditable(false); cpnl3.add(tNum1); cpnl3.add(new JLabel(" < V\u2081 >:")); tVol1.setEditable(false); cpnl3.add(tVol1); cpnl3.add(new JLabel(" < \u03c1\u2081 >:")); tRho1.setEditable(false); cpnl3.add(tRho1); cpnl3.add(new JLabel(" < U\u2081/N\u2081 >:")); tAvU1.setEditable(false); cpnl3.add(tAvU1); cpnl3.add(new JLabel(" < p\u2081 >:")); tAvp1.setEditable(false); cpnl3.add(tAvp1); cpnl3.add(new JLabel(" AR\u2081:")); tAcc1.setEditable(false); cpnl3.add(tAcc1); cpnl4.add(new JLabel(" < N\u2082 >:")); tNum2.setEditable(false); cpnl4.add(tNum2); cpnl4.add(new JLabel(" < V\u2082 >:")); tVol2.setEditable(false); cpnl4.add(tVol2); cpnl4.add(new JLabel(" < \u03c1\u2082 >:")); tRho2.setEditable(false); cpnl4.add(tRho2); cpnl4.add(new JLabel(" < U\u2082/N\u2082 >:")); tAvU2.setEditable(false); cpnl4.add(tAvU2); cpnl4.add(new JLabel(" < p\u2082 >:")); tAvp2.setEditable(false); cpnl4.add(tAvp2); cpnl4.add(new JLabel(" AR\u2082:")); tAcc2.setEditable(false); cpnl4.add(tAcc2); lStatus.setFont(new Font("Courier", 0, 12)); cpnl.add(lStatus); tpnl = new JPanel(); // two canvas panel tpnl.setLayout(new GridLayout(1, 2)); box.add(tpnl, BorderLayout.CENTER); canvas1 = new XYZCanvasMC(); tpnl.add(canvas1); canvas2 = new XYZCanvasMC(); tpnl.add(canvas2); timer = new Timer(delay, this); // start and stop the timer, so we can call restart() later timer.start(); timer.stop(); } public void start() { } public void update(Graphics g) { paint(g); } DecimalFormat df = new DecimalFormat("###0.000"); public void paint(Graphics g) { lStatus.setText(" step: " + mc1.step + ", " + "N\u2081: " + mc1.N + ", " + "N\u2082: " + mc2.N + ", " + "V\u2081: " + df.format(mc1.Vol) + ", " + "V\u2082: " + df.format(mc2.Vol) + ", " + "\u03c1\u2081: " + df.format(mc1.rho) + ", " + "\u03c1\u2082: " + df.format(mc2.rho) + ", " + "U\u2081/N\u2081: " + df.format(mc1.U/mc1.N) + ", " + "U\u2082/N\u2082: " + df.format(mc2.U/mc2.N) + ", " + "p\u2081: " + df.format(mc1.p) + ", " + "p\u2082: " + df.format(mc2.p) + ";"); double rate; tAccVol.setText(" " + df.format(mc1.vMvTot > 0 ? 100.0*mc1.vMvAcc/mc1.vMvTot : 0.0) + "%"); rate = (mc1.nMvTot+mc2.nMvTot) > 0 ? 1.0*(mc1.nMvAcc + mc2.nMvAcc)/(mc1.nMvTot + mc2.nMvTot) : 0.0; tAccNum.setText(" " + df.format(rate) + "%"); tNum1.setText(" " + df.format(mc1.avNum.getAve())); tVol1.setText(" " + df.format(mc1.avVol.getAve())); tRho1.setText(" " + df.format(mc1.avRho.getAve())); tAvU1.setText(" " + df.format(mc1.avU.getAve()/mc1.N)); tAvp1.setText(" " + df.format(mc1.avp.getAve())); tAcc1.setText(" " + df.format(mc1.mcTot > 0 ? 100.0*mc1.mcAcc/mc1.mcTot : 0.0) + "%"); tNum2.setText(" " + df.format(mc2.avNum.getAve())); tVol2.setText(" " + df.format(mc2.avVol.getAve())); tRho2.setText(" " + df.format(mc2.avRho.getAve())); tAvU2.setText(" " + df.format(mc2.avU.getAve()/mc2.N)); tAvp2.setText(" " + df.format(mc2.avp.getAve())); tAcc2.setText(" " + df.format(mc2.mcTot > 0 ? 100.0*mc2.mcAcc/mc2.mcTot : 0.0) + "%"); refreshCanvases(false); cpnl.repaint(); } Random rng = new Random(); public void actionPerformed(ActionEvent e) { Object src = e.getSource(); boolean adjScale = false; if (src == timer) { for (int i = 1; i <= speed; i++) { // sample a few steps mc1.metropolis(); mc2.metropolis(); if (nstVMove > 0 && i % nstVMove == 0) mc1.gibbsVolMove(mc2, volAmp); if (nstNMove > 0 && i % nstNMove == 0) { if ((int) (rng.nextDouble() * 2.0) == 1) { mc2.gibbsNumMove(mc1); } else { mc1.gibbsNumMove(mc2); } } } repaint(); return; } if (src == tTemp || src == bReset) { double kT = Double.parseDouble(tTemp.getText().trim()); if (kT < 1e-8) { kT = 1e-8; tTemp.setText(" " + kT); } mc1.kT = kT; mc1.clearData(); mc2.kT = kT; mc2.clearData(); } if (src == tRho) { double rho = Double.parseDouble(tRho.getText().trim()); if (rho < 1e-3) { rho = 1e-3; tRho.setText(" " + rho); } if (rho > 1.2) { rho = 1.2; tRho.setText(" " + rho); } double vn = (mc1.N + mc2.N)/rho; double vo1 = mc1.L*mc1.L*mc1.L; double vo2 = mc2.L*mc2.L*mc2.L; double vn1 = vn*vo1/(vo1 + vo2); double vn2 = vn*vo2/(vo1 + vo2); if (src == bReset) { mc1.setDensity(rho); mc2.setDensity(rho); } else { mc1.setDensity(mc1.N/vn1); mc2.setDensity(mc2.N/vn2); } mc1.clearData(); mc2.clearData(); adjScale = true; } if (src == tSpeed || src == bReset) { speed = Integer.parseInt(tSpeed.getText().trim()); if (speed < 1) { speed = 1; tSpeed.setText(" " +speed); } } if (src == tMCAmp) { double amp = Double.parseDouble(tMCAmp.getText().trim()); if (amp < 0.) { amp = 0; tMCAmp.setText(" " +amp); } if (amp > mc1.L) { amp = mc1.L; tMCAmp.setText(" " + amp); } if (amp > mc2.L) { amp = mc2.L; tMCAmp.setText(" " + amp); } mc1.mcAmp = amp; mc1.clearData(); mc2.mcAmp = amp; mc2.clearData(); } if (src == tVolAmp) { double amp = Double.parseDouble(tVolAmp.getText().trim()); if (amp < 0.) { amp = 0; tVolAmp.setText(" " +amp); } volAmp = amp; } if (src == tVolFreq) { int freq = Integer.parseInt(tVolFreq.getText().trim()); if (freq <= 0) { freq = 0; tVolFreq.setText(" " +freq); } nstVMove = freq; } if (src == tNumFreq) { int freq = Integer.parseInt(tNumFreq.getText().trim()); if (freq <= 0) { freq = 0; tNumFreq.setText(" " +freq); } nstNMove = freq; } if (src == bRetime) { mc1.clearData(); mc2.clearData(); } if (src == bStart) { boolean on = bStart.getModel().isSelected(); if (on) { timer.restart(); bStart.setText("Pause"); } else { timer.stop(); bStart.setText("Resume"); } } if (src == tNum || src == bReset) { // change total # of particles if (timer.isRunning()) timer.stop(); bStart.setSelected(false); bStart.setText("Start"); int n = Integer.parseInt(tNum.getText().trim()); if (n < 2) { n = 2; tNum.setText(" " + n); } if (n % 2 == 1) { n = n/2*2; tNum.setText(" " + n); } n /= 2; double rho = Double.parseDouble(tRho.getText().trim()); if (rho < 1e-3) { rho = 1e-3; tRho.setText(" " + rho); } if (rho > 1.2) { rho = 1.2; tRho.setText(" " + rho); } mc1.init(rho, n, n*2); mc2.init(rho, n, n*2); adjScale = false; } refreshCanvases(adjScale); repaint(); } private void refreshCanvases(boolean adjScale) { refreshCanvas(canvas1, mc1, adjScale); refreshCanvas(canvas2, mc2, adjScale); } private void refreshCanvas(XYZCanvasMC canvas, LJ3Gibbs mc, boolean adjScale) { // make sure the canvas has a model if (canvas.model == null) { canvas.model = new XYZModelMC(); adjScale = true; } // change the color of the default atom double den = mc.N / (mc.L * mc.L * mc.L); double denMax = 0.9, kTMax = 10.0; // the blue-to-red transition shows the density double blue = min(den/denMax, 1.); double red = max(1 - den/denMax, 0); // the green component shows the temperature double green = min(mc.kT/kTMax, 1); canvas.model.atomDefault = new Atom(red, green, blue, 1.0, 0.7, 0.5, 1.0); canvas.refresh(mc.getXWrap(), mc.N, mc.xo, mc.xi, true, mc.moveAtom, mc.moveAcc, adjScale); } } /** Always rc = .5*L */ class LJ3Gibbs { static final int D = 3; public int N = 108; // number of particles 6^3/2 public int NMAX = 216; // number of particles in both boxes double rho = 0.2; double Vol = 540; // volume; double dt = 0.002; // time step for integrating Newton's equation public double L = 10.0; // side of the box double rc = L/2; // it's necessary to use half-box cutoff to use fast volume move public double kT = 1.2; // temperature times Boltzmann constant double x[][]; // position, from 0 to 1 double xi[] = new double [D], xo[] = new double [D]; double xwrap[][]; // wrapped coordinates double U, U12, U6; // potential energy double Vir, p; // virial and pressure double Utail; // potential energy tail correction double ptail; // pressure tail correction int step; // simulation steps double mcAmp = 0.1; // amplitude of randomly moving a particle long mcAcc = 0, mcTot = 0; Ave avU = new Ave(), avp = new Ave(), avRho = new Ave(), avNum = new Ave(), avVol = new Ave(); /** constructor */ LJ3Gibbs(double den, int n, int nmax) { init(den, n, nmax); } /** initialize system */ public void init(double den, int n, int nmax) { int i, j, k, id; step = 0; N = n; NMAX = nmax; x = new double[NMAX][D]; xwrap = new double[NMAX][D]; int N1 = (int) (pow(2*N, 1.0/D) + .999999); // # of particles per side (cubic lattic) double a = 1./N1; for (id = 0, i = 0; i < N1 && id < N; i++) // place particles on a cubic lattice for (j = 0; j < N1 && id < N; j++) for (k = 0; k < N1 && id < N; k++) { if ((i + j + k) % 2 != 0) continue; x[id][0] = a * (i + .5); x[id][1] = a * (j + .5); x[id][2] = a * (k + .5); id++; } setDensity(den); } /** set density and recalculate tail corrections */ void setDensity(double den) { rho = den; Vol = N/rho; L = pow(Vol, 1.0/D); /* compute shifts and tail corrections */ rc = .5*L; double irc = 1./rc, irc3 = irc*irc*irc, irc6 = irc3*irc3; Utail = 8.0*PI/9*rho*N*(irc6 - 3)*irc3; ptail = 16.0*PI/9*rho*rho*(irc6 - 1.5)*irc3; clearData(); energy(L, rc); /* re-compute energy */ } void clearData() { step = 0; mcAcc = mcTot = 0; vMvAcc = vMvTot = 0; nMvAcc = nMvTot = 0; avU.clear(); avp.clear(); avRho.clear(); avVol.clear(); avNum.clear(); } public int moveAtom = -1; public boolean moveAcc = false; Random rng = new Random(); /** Metropolis algorithm */ double xij[] = new double[D]; double dxi[] = new double[D]; void metropolis() { double dU = 0, dU12 = 0, dU6 = 0, dVir = 0, r2, invr2, invr6, fscal, amp; int i, j, d; /* randomly move a particle */ i = (int) (rng.nextDouble() * N); amp = mcAmp / L / rho; /* rho adjusted */ rv3Copy(xo, x[i]); rv3Rand(dxi, .5, -2*amp); rv3Diff(xi, x[i], dxi); energy_i(i, x[i], L, rc); dU12 -= dU12_i; dU6 -= dU6_i; dU -= dU_i; dVir -= dVir_i; energy_i(i, xi, L, rc); dU12 += dU12_i; dU6 += dU6_i; dU += dU_i; dVir += dVir_i; mcTot++; if (dU <= 0 || rng.nextDouble() < exp(-dU/kT) ) { mcAcc++; rv3Copy(x[i], xi); U += dU; U12 += dU12; U6 += dU6; Vir += dVir; p = calcPressure(); moveAcc = true; //System.out.println("a U " + U + " dU " + dU); //energy(L, rc); //System.out.println("b U " + U); } else { moveAcc = false; } step++; avU.add(U); avp.add(p); avRho.add(rho); avVol.add(Vol); avNum.add(1.0*N); } /** calculate potential energy and force * half-box cutoff, minimal distance image */ void energy(double l, double r_c) { int i, j, d; double invr2, invr6, r2; // clear the energy and pressure U = U12 = U6 = Vir = 0.0; // loop over pairs of particles and compute energy for (i = 0; i < N - 1; i++) { for (j = i + 1; j < N; j++) { rv3Diff(xij, x[i], x[j]); vpbc(xij, l); r2 = vsqr(xij); if (r2 > r_c*r_c) continue; invr2 = 1.0/r2; invr6 = invr2 * invr2 * invr2; U6 += invr6; U12 += invr6*invr6; } } Vir = 48.0 * U12 - 24.0 * U6; U = 4.0 * U12 - 4.0 * U6 + Utail; p = calcPressure(); } /* partial energy invovled in particle i * no tail correction */ double dU_i, dU12_i, dU6_i, dVir_i; void energy_i(int i, double ri[], double newl, double newrc) { int j; double r2, invr6, invr2; dU12_i = dU6_i = 0.0; for (j = 0; j < N; j++) { if (i == j) continue; rv3Diff(xij, ri, x[j]); vpbc(xij, newl); r2 = vsqr(xij); if (r2 < newrc*newrc) { invr2 = 1.0/r2; invr6 = invr2 * invr2 * invr2; dU12_i += invr6*invr6; dU6_i += invr6; } } dVir_i = 48.0 * dU12_i - 24.0 * dU6_i; dU_i = 4.0 * dU12_i - 4.0 * dU6_i; } /** compute the energy under the new volume * save temporary quantities to Voln, Utailn, etc. * return the potential energy change */ double Voln, Utailn, ptailn, Un, U12n, U6n, Virn, Ln, rcn; int Nn; private double setDen2(int nn, double voln, double dU12, double dU6) { Nn = nn; Voln = voln; Ln = pow(Voln, 1.0/D); double denn = Nn/Voln; /* compute shifts and tail corrections */ rcn = .5*Ln; double irc = 1./rcn, irc3 = irc*irc*irc, irc6 = irc3*irc3; Utailn = 8.0*PI/9*denn*Nn*(irc6 - 3)*irc3; ptailn = 16.0*PI/9*denn*denn*(irc6 - 1.5)*irc3; if (Nn != N) { U6n = U6 + dU6; U12n = U12 + dU12; } else { double sc6 = (Vol*Vol)/(Voln*Voln); U6n = U6*sc6; U12n = U12*sc6*sc6; } Un = 4.0*(U12n - U6n) + Utailn; Virn = 48.0*U12n - 24.0*U6n; return Un - U; } /** commit to new number of particles/volume */ private void commitNew(int newN, double newV, boolean refresh) { boolean nmove = (newN != N); boolean verbose = false; if (verbose) { System.out.println("N: " + N + " vs " + newN + " Vol: " + Vol + " vs " + newV + " U: " + U + " vs " + Un); } N = newN; Vol = newV; L = pow(Vol, 1.0/3); rho = N/Vol; U12 = U12n; U6 = U6n; U = Un; Vir = Virn; rc = rcn; Utail = Utailn; ptail = ptailn; if (refresh) energy(L, rc); if (verbose) { System.out.println(" U(correct): " + U + " U: " + Un + " Utail: " + Utail + ", " + Utailn + " Vir: " + Vir + " Virn: " + Virn); System.out.println(" U12: " + U12 + " U12n: " + U12n + " U6: " + U6 + " U6n: " + U6n); } } /** change volume */ long vMvTot = 0, vMvAcc = 0; void gibbsVolMove(LJ3Gibbs lj2, double volstep) { double lnV, Vtot, Vo1, Vo2, Vn1, Vn2, sc6a, sc6b, dU1, dU2; Vo1 = Vol; Vo2 = lj2.Vol; Vtot = Vo1 + Vo2; lnV = log(Vo1/Vo2) + (2*rng.nextDouble() - 1) * volstep; Vn1 = Vtot/(exp(-lnV) + 1); Vn2 = Vtot - Vn1; dU1 = setDen2(N, Vn1, 0, 0); dU2 = lj2.setDen2(lj2.N, Vn2, 0, 0); double xp = -(dU1 + dU2)/kT + (N + 1)*log(Vn1/Vo1) + (lj2.N + 1)*log(Vn2/Vo2); vMvTot++; if (xp > 0.0 || rng.nextDouble() < exp(xp)) { vMvAcc++; boolean refresh = (vMvAcc % 100 == 0); commitNew(N, Vn1, refresh); lj2.commitNew(lj2.N, Vn2, refresh); } } long nMvTot = 0, nMvAcc = 0; /** try to move a particle from here to lj2 */ void gibbsNumMove(LJ3Gibbs lj2) { int i; double dU1, dU2; nMvTot++; if (N == 0) return; i = (int) (rng.nextDouble() * N); /* energy of removing i */ energy_i(i, x[i], L, rc); setDen2(N - 1, Vol, -dU12_i, -dU6_i); dU1 = Un - U; /* energy of inserting i */ double rn[] = new double[3]; rv3Rand(rn, 0, 1); /* randomly in the box */ lj2.energy_i(lj2.N+1000, rn, L, rc); lj2.setDen2(lj2.N + 1, lj2.Vol, lj2.dU12_i, lj2.dU6_i); dU2 = lj2.Un - lj2.U; double xp = -(dU1 + dU2)/kT + log(N/(lj2.N + 1.0)*Vol/(lj2.Vol)); if (rng.nextDouble() < exp(xp)) { if (lj2.N >= lj2.NMAX) { System.out.print("Fatal: lj2 over flow, N " + lj2.N + " NMAX" + lj2.NMAX); return; } nMvAcc++; boolean refresh = (vMvAcc % 100 == 0); /* add particle to box 2*/ rv3Copy(lj2.x[lj2.N], rn); lj2.commitNew(lj2.N+1, lj2.Vol, refresh); /* remove this particle */ if (i != N - 1) rv3Copy(x[i], x[N - 1]); commitNew(N - 1, Vol, refresh); } } double calcPressure() { return rho*kT + Vir/(3*Vol) + ptail; } private void rv3Diff(double c[], double a[], double b[]) { c[0] = a[0] - b[0]; c[1] = a[1] - b[1]; c[2] = a[2] - b[2]; } private void rv3Copy(double a[], double b[]) { a[0] = b[0]; a[1] = b[1]; a[2] = b[2]; } private void rv3Rand(double a[], double offset, double amp) { a[0] = amp*(rng.nextDouble() - offset); a[1] = amp*(rng.nextDouble() - offset); a[2] = amp*(rng.nextDouble() - offset); } /** image with the minimal distance */ double pbc(double a, double l) { return l * (a - ((int)(a + 1000.5) - 1000)); } void vpbc(double v[], double l) { v[0] = pbc(v[0], l); v[1] = pbc(v[1], l); v[2] = pbc(v[2], l); } double vsqr(double v[]) { return v[0]*v[0] + v[1]*v[1] + v[2]*v[2]; } /** wrap coordinates back into box */ public double [][] getXWrap() { double xx; for (int i = 0; i < N; i++) for (int d = 0; d < D; d++) { xx = x[i][d] + 1000.0; xwrap[i][d] = (xx - (int) xx) * L; } for (int i = N; i < NMAX; i++) for (int d = 0; d < D; d++) { xwrap[i][d] = 0; } for (int d = 0; d < D; d++) { xx = xi[d] + 1000.0; xi[d] = (xx - (int) xx) * L; xx = xo[d] + 1000.0; xo[d] = (xx - (int) xx) * L; } return xwrap; } } /** Average and standard deviation */ class Ave { double cnt, xsm, x2sm; public void clear() { cnt = xsm = x2sm = 0; } public void add(double x) { cnt += 1; xsm += x; x2sm += x*x; } public double getAve() { return (cnt > 0) ? xsm/cnt : 0; } public double getVar() { if (cnt <= 1) return 0; double av = xsm/cnt; return x2sm/cnt - av*av; } public double getStd() { return Math.sqrt( getVar() ); } } /** Panel for drawing particles */ class XYZCanvas extends JPanel implements MouseListener, MouseMotionListener, MouseWheelListener { Image img; Graphics imgG; Dimension imgSize; double realSpan; // size of a real span (saved as a copy) double zoomScale = 1.0; public Matrix3D viewMatrix = new Matrix3D(); // view matrix private Matrix3D tmpMatrix = new Matrix3D(); // temporary matrix int mouseX, mouseY; // mouse position XYZModel model; public static final long serialVersionUID = 2L; public XYZCanvas() { super(); addMouseListener(this); addMouseMotionListener(this); addMouseWheelListener(this); } /** Prepare a buffer for the image, return if the canvas is ready * Only need to call this when the size of the canvas is changed, * Since we automatically detect the size change in paint() * only call this on start */ public boolean newImgBuf() { Dimension sz = getSize(); if (sz.width == 0 || sz.height == 0) return false; // quit if the current image already has the right size if (img != null && imgG != null && sz.equals(imgSize)) return true; img = createImage(sz.width, sz.height); if (imgG != null) imgG.dispose(); imgG = img.getGraphics(); imgSize = sz; return true; } public void update(Graphics g) { if (img == null) g.clearRect(0, 0, getSize().width, getSize().height); paintComponent(g); } protected void paintComponent(Graphics g) { super.paintComponent(g); if (model != null) { newImgBuf(); // refresh the image buffer if necessary // compute the real-to-screen ratio, this variable differs // from model.real2Screen by zoomScale Dimension sz = getSize(); double real2Screen0 = model.getScaleFromSpan(realSpan, sz.width, sz.height); model.setMatrix(viewMatrix, real2Screen0 * zoomScale, sz.width/2, sz.height/2); imgG.setColor(Color.BLACK); imgG.fillRect(0, 0, sz.width, sz.height); model.paint(imgG); g.drawImage(img, 0, 0, this); } } /** Refresh the coordinates * x[][] is the wrapped coordinates * `n' may be less than x.length */ public void refresh(double x[][], int n, boolean center, boolean adjScale) { if (model == null) { model = new XYZModel(); adjScale = true; } model.updateXYZ(x, n, center); if ( adjScale ) realSpan = model.getSpan(x, n); repaint(); } /** Event handling */ public void mouseClicked(MouseEvent e) { } public void mousePressed(MouseEvent e) { mouseX = e.getX(); mouseY = e.getY(); // consume this event so that it will not be processed // in the default manner e.consume(); } public void mouseReleased(MouseEvent e) { } public void mouseEntered(MouseEvent e) { } public void mouseExited(MouseEvent e) { } public void mouseDragged(MouseEvent e) { int x = e.getX(); int y = e.getY(); tmpMatrix.unit(); tmpMatrix.xrot(360.0 * (mouseY - y) / getSize().height); tmpMatrix.yrot(360.0 * (x - mouseX) / getSize().width); viewMatrix.mult(tmpMatrix); repaint(); mouseX = x; mouseY = y; e.consume(); } public void mouseMoved(MouseEvent e) { } public void mouseWheelMoved(MouseWheelEvent e) { int notches = e.getWheelRotation(); if ((zoomScale -= 0.05f*notches) < 0.09999f) zoomScale = 0.1f; repaint(); } } /** Panel for drawing particles */ class XYZCanvasMC extends XYZCanvas { public static final long serialVersionUID = 3L; /** Get 3D position * `n' may be less than x.length */ private double [][] appendXMove(double x[][], int n, double xMove[]) { int D = x[0].length; double r[][] = new double[n + 1][3]; // copy the direct coordinates for (int i = 0; i < n; i++) { r[i][0] = x[i][0]; r[i][1] = x[i][1]; r[i][2] = (D == 2) ? 0.0 : x[i][2]; } if (xMove == null) xMove = x[0]; // the last is the trial position r[n][0] = xMove[0]; r[n][1] = xMove[1]; r[n][2] = (D == 2) ? 0 : xMove[2]; return r; } /** Refresh the coordinates * x[][] is the wrapped coordinates * `n' may be less than x.length */ public void refresh(double x[][], int n, double xo[], double xi[], boolean center, int mvAtom, boolean mvAcc, boolean adjScale) { if (model == null) { model = new XYZModelMC(); adjScale = true; } XYZModelMC modelMC = (XYZModelMC) model; double xyz[][] = appendXMove(x, n, mvAcc ? xo : xi); modelMC.updateXYZ(xyz, n + 1, center); modelMC.setMoveAtom(mvAtom, mvAcc); if ( adjScale ) realSpan = model.getSpan(x, n); repaint(); } } /** A set of atoms with 3D coordinates */ class XYZModel { double realXYZ[][]; // 3D real coordinates [np][3] int screenXYZ[][]; // 3D screen coordinates in pixels [np][3] // only the first two dimensions are used in drawing int zOrder[]; // z-order, larger is closer to the viewer int np = -1; boolean transformed; // rotation/scaling/translation matrix for the conversion from realXYZ[] to screenXYZ[] Matrix3D mat = new Matrix3D(); double real2Screen = 50.0; // real size --> screen size double ballSize = 0.5; // ball size (radius) in terms of the real coordinates // 0.5 for hard spheres Atom atomDefault = new Atom(0.1, 0.7, 0.1, 1.0, 0.5, 0.5, 1.0); Atom atoms[]; // for colors of atoms XYZModel() {} /** Set the color of a particular atom */ void setAtom(int i, Atom atom) { if ( i >= 0 && i < atoms.length ) atoms[i] = atom; } /** Refresh coordinates * x[0..n-1][3] is a three-dimensional vector * n can be less than x.length */ void updateXYZ(double x[][], int n, boolean center) { if (n != np) { //System.out.printf("XYZModel.updateXYZ n %d --> %d, (%g, %g, %g) (%g, %g, %g)\n", np, n, x[0][0], x[0][1], x[0][2], x[n-1][0], x[n-1][1], x[n-1][2]); np = n; realXYZ = new double [np][3]; screenXYZ = new int [np][3]; zOrder = new int [np]; atoms = new Atom [np]; for (int i = 0; i < np; i++) atoms[i] = atomDefault; } for (int d = 0; d < 3; d++) { double xc = 0; if ( center ) { for (int i = 0; i < np; i++) xc += x[i][d]; xc /= np; } for (int i = 0; i < np; i++) realXYZ[i][d] = x[i][d] - xc; } transformed = false; } /** Set the view matrix * `s' is the scaling factor of translating real coordinates * to the screen coordinates * (x0, y0) the screen coordinates of the center */ void setMatrix(Matrix3D viewMat, double s, double x0, double y0) { mat.unit(); mat.mult(viewMat); mat.scale(s, s, s); real2Screen = s; mat.translate(x0, y0, 0); transformed = false; } /** Get the span of the model * `n' may be less than x.length */ double getSpan(double x[][], int n) { int dim = x[0].length; double realSpan = 0, del, fw, fh; for (int d = 0; d < dim; d++) { double xmin = 1e30, xmax = -1e30; for (int i = 0; i < n; i++) if ( x[i][d] < xmin ) xmin = x[i][d]; else if ( x[i][d] > xmax ) xmax = x[i][d]; if ( (del = xmax - xmin) > realSpan ) realSpan = del; } return realSpan; } /** Translate a given span, return the real-to-screen ratio * `w' and `h' are the width and height of the screen in pixels */ double getScaleFromSpan(double realSpan, int w, int h) { realSpan += ballSize * 2; // add two radii double fw = w / realSpan; double fh = h / realSpan; double facShrink = 0.9; // shrink a bit for the margin return (fw < fh ? fw : fh) * facShrink; } /** Compute the Z-order */ void getZOrder() { // transform the coordinates if ( !transformed ) { mat.transform(realXYZ, screenXYZ, np); transformed = true; } // bubble sort z-order // zOrder[0] is the fartherest from the viewer // zOrder[np - 1] is the nearest to the viewer for (int i = 0; i < np; i++) zOrder[i] = i; for (int i = 0; i < np; i++) { // find the particle with the smallest z int jm = i, k; for (int j = i + 1; j < np; j++) if (screenXYZ[zOrder[j]][2] < screenXYZ[zOrder[jm]][2]) jm = j; if (jm != i) { k = zOrder[i]; zOrder[i] = zOrder[jm]; zOrder[jm] = k; } } } /** Draw atom */ void drawAtom(Graphics g, int iz) { int i = zOrder[iz]; Atom atom = atoms[i]; if (atom == null) return; double zMin = screenXYZ[zOrder[0]][2]; // fartherest from the viewer double zMax = screenXYZ[zOrder[np - 1]][2]; // closest to the viewer int greyScale = Atom.nBalls - 1; if (zMin != zMax) // zMin == zMax means the two-dimensional case greyScale = (int) (Atom.nBalls * (screenXYZ[i][2] - zMin) / (zMax - zMin) - 1e-6); // the atom closest to the viewer has a greyScale of Atom.nBalls - 1 // the atom fartherest from the viewer has a greyScale of 0 double radius = ballSize * atom.relRadius * real2Screen; atom.paint(g, screenXYZ[i][0], screenXYZ[i][1], greyScale, radius); } /** Paint this model to the graphics context `g' */ void paint(Graphics g) { if (realXYZ == null || np <= 0) return; getZOrder(); for (int iz = 0; iz < np; iz++) drawAtom(g, iz); } } /** XYZModel for a typical MC simulation */ class XYZModelMC extends XYZModel { int moveAtom = -1; boolean moveAcc = false; // failed move Atom atomTrial = new Atom(1.0, 0.1, 0.1, 1.0, 0.5, 0.3, 2.0); // failed trial postion Atom atomFailed = new Atom(0.0, 0.0, 0.4); // new and old position // successful move Atom atomMoved = new Atom(0.1, 1.0, 0.5); // trial and new position Atom atomOldPos = new Atom(0.3, 0.3, 0.3, 1.0, 0.5, 0.5, 2.0); // old position /** Constructor from the real coordinates * x[][3] is a three-dimensional vector */ XYZModelMC() { atomDefault = new Atom(0.1, 0.2, 1.0, 1.0, 0.6, 0.4, 1.0); } /** Set the atom */ public void setMoveAtom(int mvAtom, boolean mvAcc) { moveAtom = mvAtom; moveAcc = mvAcc; if (mvAtom < 0) return; for (int i = 0; i < np; i++) atoms[i] = atomDefault; if ( moveAtom < 0 || moveAtom >= np - 1 ) { atoms[np - 1] = null; } else { atoms[moveAtom] = moveAcc ? atomMoved : atomFailed; atoms[np - 1] = moveAcc ? atomOldPos : atomTrial; } } } /** Atom: a ball */ class Atom { private final static int R = 120; private final static int hx = 45; // (hx, hy) is the offset of the spot light from the center private final static int hy = 45; private static int maxr; // maximal distance from the spotlight public final static int nBalls = 16; // shades of grey // spotlight brightness (0.0, 1.0) // 1.0 means the spotlight is pure white // 0.0 means the spotlight is the same as the solid color double spotlightAmp = .4; // color damping along the distance from the spotlight double rContrast = 0.7; // z-depth contrast (0.0, inf) // inf means the fartherest atom is completely dark // 0.0 means the fartherest atom is the same as the // nearest atom double zContrast = 2.0; private static byte data[]; static { // data[] is a bitmap image of the ball of radius R data = new byte[R * 2 * R * 2]; for (int Y = -R; Y < R; Y++) { int x0 = (int) (Math.sqrt(R * R - Y * Y) + 0.5); for (int X = -x0; X < x0; X++) { // sqrt(x^2 + y^2) gives distance from the spot light int x = X + hx, y = Y + hy; int r = (int) (Math.sqrt(x * x + y * y) + 0.5); // set the maximal intensity to the maximal distance // (in pixels) from the spot light if (r > maxr) maxr = r; data[(Y + R) * (R * 2) + (X + R)] = (r <= 0) ? 1 : (byte) r; } } } // the following variables are atom dependent private int Rl = 100, Gl = 100, Bl = 100; private Image balls[]; // 0..nBalls-1, at different z distances /** Constructor */ Atom(int r, int g, int b) { setRGB(r, g, b); } /** colors that range from 0 to 1 */ Atom(double r, double g, double b) { setRGB(r, g, b); } double relRadius = 1; // only used to save the radius Atom(double r, double g, double b, double rad) { this(r, g, b); relRadius = rad; } Atom(double r, double g, double b, double rad, double rcontrast, double spotlight, double zcontrast) { relRadius = rad; rContrast = rcontrast; spotlightAmp = spotlight; zContrast = zcontrast; setRGB(r, g, b); } /** Set color */ void setRGB(int r, int g, int b) { Rl = r; Gl = g; Bl = b; makeBalls(); } void setRGB(double r, double g, double b) { Rl = (int) (256*r - 1e-6); Gl = (int) (256*g - 1e-6); Bl = (int) (256*b - 1e-6); makeBalls(); } /** Linearly interpolate colors */ private int blend(double fg, double bg, double fgfactor) { return (int) (bg + (fg - bg) * fgfactor); } // need a component instance to call createImage() private static Component component = new Applet(); /** Prepare ball images with different sizes */ private void makeBalls() { balls = new Image[nBalls]; byte red[] = new byte[256]; byte green[] = new byte[256]; byte blue[] = new byte[256]; for (int id = 0; id < nBalls; id++) { // smaller `b' means closer to black // if id == 0 (fartherest from the viewer) // b = 1/(1 + zContrast) // the outer blend() gives a color close to bgGrey // if id == nBalls - 1 (closest to the viewer), // b = 1, the outer blend() gives the color of // the inner blend() double b = (zContrast*id/(nBalls - 1) + 1) / (zContrast + 1); for (int i = maxr; i >= 0; --i) { // closeness to the spotlight double q = 1 - 1. * i / maxr; // dampness of the color along the radius double p = 1 - rContrast * i / maxr; // contrast of the spotlight // if i == 0 (closest to the spotlight), // d = 1.0 - spotlightAmp, the inner // blend() gives a color close to 255 // (if spotlightAmp == 1). // if i == maxr (fartherest from the spotlight), // d = 1.0, the inner blend() gives // the foreground color, i.e., Rl, Gl, Bl // Thus, the inner blend() depends on the distance // from the spotlight, i == 0 means to be closest // to the spotlight double d = 1 - q * spotlightAmp; red[i] = (byte) blend(blend(Rl*p, 255, d), 0, b); green[i] = (byte) blend(blend(Gl*p, 255, d), 0, b); blue[i] = (byte) blend(blend(Bl*p, 255, d), 0, b); } // 256 color model IndexColorModel model = new IndexColorModel( 8, maxr + 1, red, green, blue, 0); balls[id] = component.createImage( new MemoryImageSource(R * 2, R * 2, model, data, 0, R * 2) ); } } /** Draw a ball at screen coordinate (x, y) with a ball index `id' * (0, 0) represents the top-left corner * `x', `y', `radius' are given in pixels * the ball index (gray code) `id' can be 0 to 15 */ void paint(Graphics gc, int x, int y, int id, double radius) { if (balls == null) makeBalls(); Image img = balls[id]; // id = [0..15] int size = (int) (radius * 2 + .5); gc.drawImage(img, x - size/2, y - size/2, size, size, null); //System.out.println("" + x + " " + y + " " + id + " " + radius); } } class Matrix3D { double xx, xy, xz, xo; double yx, yy, yz, yo; double zx, zy, zz, zo; static final double pi = 3.14159265; /** Create a new unit matrix */ Matrix3D() { xx = 1.0; yy = 1.0; zz = 1.0; } /** Scale along each axis independently */ void scale(double xf, double yf, double zf) { xx *= xf; xy *= xf; xz *= xf; xo *= xf; yx *= yf; yy *= yf; yz *= yf; yo *= yf; zx *= zf; zy *= zf; zz *= zf; zo *= zf; } /** Translate the origin */ void translate(double x, double y, double z) { xo += x; yo += y; zo += z; } /** Rotate theta degrees around the y axis */ void yrot(double theta) { theta *= (pi / 180); double ct = Math.cos(theta); double st = Math.sin(theta); double Nxx = (xx * ct + zx * st); double Nxy = (xy * ct + zy * st); double Nxz = (xz * ct + zz * st); double Nxo = (xo * ct + zo * st); double Nzx = (zx * ct - xx * st); double Nzy = (zy * ct - xy * st); double Nzz = (zz * ct - xz * st); double Nzo = (zo * ct - xo * st); xo = Nxo; xx = Nxx; xy = Nxy; xz = Nxz; zo = Nzo; zx = Nzx; zy = Nzy; zz = Nzz; } /** Rotate theta degrees about the x axis */ void xrot(double theta) { theta *= (pi / 180); double ct = Math.cos(theta); double st = Math.sin(theta); double Nyx = (yx * ct + zx * st); double Nyy = (yy * ct + zy * st); double Nyz = (yz * ct + zz * st); double Nyo = (yo * ct + zo * st); double Nzx = (zx * ct - yx * st); double Nzy = (zy * ct - yy * st); double Nzz = (zz * ct - yz * st); double Nzo = (zo * ct - yo * st); yo = Nyo; yx = Nyx; yy = Nyy; yz = Nyz; zo = Nzo; zx = Nzx; zy = Nzy; zz = Nzz; } /** Rotate theta degrees about the z axis */ void zrot(double theta) { theta *= pi / 180; double ct = Math.cos(theta); double st = Math.sin(theta); double Nyx = (yx * ct + xx * st); double Nyy = (yy * ct + xy * st); double Nyz = (yz * ct + xz * st); double Nyo = (yo * ct + xo * st); double Nxx = (xx * ct - yx * st); double Nxy = (xy * ct - yy * st); double Nxz = (xz * ct - yz * st); double Nxo = (xo * ct - yo * st); yo = Nyo; yx = Nyx; yy = Nyy; yz = Nyz; xo = Nxo; xx = Nxx; xy = Nxy; xz = Nxz; } /** Multiply this matrix by a second: M = M*R */ void mult(Matrix3D rhs) { double lxx = xx * rhs.xx + yx * rhs.xy + zx * rhs.xz; double lxy = xy * rhs.xx + yy * rhs.xy + zy * rhs.xz; double lxz = xz * rhs.xx + yz * rhs.xy + zz * rhs.xz; double lxo = xo * rhs.xx + yo * rhs.xy + zo * rhs.xz + rhs.xo; double lyx = xx * rhs.yx + yx * rhs.yy + zx * rhs.yz; double lyy = xy * rhs.yx + yy * rhs.yy + zy * rhs.yz; double lyz = xz * rhs.yx + yz * rhs.yy + zz * rhs.yz; double lyo = xo * rhs.yx + yo * rhs.yy + zo * rhs.yz + rhs.yo; double lzx = xx * rhs.zx + yx * rhs.zy + zx * rhs.zz; double lzy = xy * rhs.zx + yy * rhs.zy + zy * rhs.zz; double lzz = xz * rhs.zx + yz * rhs.zy + zz * rhs.zz; double lzo = xo * rhs.zx + yo * rhs.zy + zo * rhs.zz + rhs.zo; xx = lxx; xy = lxy; xz = lxz; xo = lxo; yx = lyx; yy = lyy; yz = lyz; yo = lyo; zx = lzx; zy = lzy; zz = lzz; zo = lzo; } /** Recover the unit matrix */ void unit() { xo = 0; xx = 1; xy = 0; xz = 0; yo = 0; yx = 0; yy = 1; yz = 0; zo = 0; zx = 0; zy = 0; zz = 1; } /** Transform np points from v into tv. * v contains the input coordinates in floating point. * Three successive entries in the array constitute a point. * tv ends up holding the transformed points as integers; * three successive entries per point */ void transform(double v[][], int tv[][], int np) { // np can be different from v.length for (int i = 0; i < np; i++) { double x = v[i][0], y = v[i][1], z = v[i][2]; tv[i][0] = (int) (xx * x + xy * y + xz * z + xo); tv[i][1] = (int) (yx * x + yy * y + yz * z + yo); tv[i][2] = (int) (zx * x + zy * y + zz * z + zo); } } }