// ComPhys File: Eigen2.java // Chapter 18: Bound state solutions of Schroedinger's Equation // Improved program uses matching to find bound states of // the Lennard-Jones potential import java.applet.*; import java.awt.*; import java.awt.event.*; import java.util.*; import comphys.*; public class Eigen2 extends Applet implements ActionListener, AdjustmentListener, ItemListener { double gamma = 0.005; // quantum parameter double dx = 0.01; // integration step size double xmin = 0.9; // minimum x to be plotted double xmax = 2.5; // maximum x to be plotted double V (double x) { // Scaled Lennard-Jones potential return 4 * (Math.pow(1/x, 12) - Math.pow(1/x, 6)); } Graphics gr; int leftPixel; int rightPixel; int topPixel; int bottomPixel; void plot_potential () { gr.setColor(Color.red); int ixOld = 0; int iyOld = 0; boolean first = true; for (int ix = leftPixel; ix <= rightPixel; ix++) { double x = xmin + ix * (xmax - xmin) / (double) (rightPixel - leftPixel); double y = V(x); y = 0.3 - y; y *= (bottomPixel - topPixel) / 1.5; int iy = (int) y; if (first) gr.drawLine(ix, iy, ix, iy); else gr.drawLine(ixOld, iyOld, ix, iy); first = false; ixOld = ix; iyOld = iy; } int yZero = (int) (0.3 * (bottomPixel - topPixel) / 1.5); gr.drawLine(leftPixel, yZero, rightPixel, yZero); } double E = -0.9; int yE; // y coordinate of E on plot void plot_energy () { double y = 0.3 - E; y *= (bottomPixel - topPixel) / 1.5; yE = (int) y; gr.setColor(Color.green); gr.drawLine(leftPixel, yE, rightPixel, yE); } double dE = 0.01; // step size in E int parity = 1; double xLeft = 0.95; double xMatch = 1.2; double xRight = 2.0; int nLeft, nRight; double[] psiLeft, psiRight, dpsiLeft, dpsiRight; void plot_psi () { if (psiLeft == null) return; gr.setColor(Color.blue); int ixOld = 0; int iyOld = 0; boolean first = true; for (int n = 0; n <= nLeft; n++) { double x = xLeft + n * (xMatch - xLeft) / (double) nLeft; x = leftPixel + (x - xmin) / (xmax - xmin) * (rightPixel - leftPixel); int ix = (int) x; double y = E + psiLeft[n]; y = 0.3 - y; y *= (bottomPixel - topPixel) / 1.5; int iy = (int) y; if (first) gr.drawLine(ix, iy, ix, iy); else gr.drawLine(ixOld, iyOld, ix, iy); first = false; ixOld = ix; iyOld = iy; } if (!drawSpectrum) gr.setColor(Color.magenta); ixOld = 0; iyOld = 0; first = true; for (int n = nRight; n >= 0; n--) { double x = xMatch + n * (xRight - xMatch) / (double) nRight; x = leftPixel + (x - xmin) / (xmax - xmin) * (rightPixel - leftPixel); int ix = (int) x; double y = E + psiRight[n]; y = 0.3 - y; y *= (bottomPixel - topPixel) / 1.5; int iy = (int) y; if (first) gr.drawLine(ix, iy, ix, iy); else gr.drawLine(ixOld, iyOld, ix, iy); first = false; ixOld = ix; iyOld = iy; } } double misMatch; int signOfPsiLeft; int signOfPsiRight; void integrate () { setLeftRightMatch(); nLeft = (int) Math.round((xMatch - xLeft) / dx); if (psiLeft == null || psiLeft.length < nLeft + 1) { psiLeft = new double[nLeft + 1]; dpsiLeft = new double[nLeft + 1]; } nRight = (int) Math.round((xRight - xMatch) / dx); if (psiRight == null || psiRight.length < nRight + 1) { psiRight = new double[nRight + 1]; dpsiRight = new double[nRight + 1]; } // integrate to find solution on left of matching point psiLeft[0] = 0; dpsiLeft[0] = 1e-10; for (int n = 0; n < nLeft; n++) { double dx = (xMatch - xLeft) / nLeft; double x = xLeft + n * dx; double d2psi = (V(x) - E) / gamma * psiLeft[n]; dpsiLeft[n + 1] = dpsiLeft[n] + d2psi * dx; psiLeft[n + 1] = psiLeft[n] + dpsiLeft[n + 1] * dx; } signOfPsiLeft = psiLeft[nLeft] > 0 ? 1 : -1; // integrate to find solution on right of matching point psiRight[nRight] = 0; dpsiRight[nRight] = 1e-10; for (int n = nRight; n > 0; n--) { double dx = (xRight - xMatch) / nRight; double x = xMatch + n * dx; double d2psi = (V(x) - E) / gamma * psiRight[n]; dpsiRight[n - 1] = dpsiRight[n] - d2psi * dx; psiRight[n - 1] = psiRight[n] - dpsiRight[n - 1] * dx; } signOfPsiRight = psiRight[0] > 0 ? 1 : -1; // get solutions to match in magnitude double scale = psiLeft[nLeft] / psiRight[0]; for (int n = nRight; n >= 0; n--) { psiRight[n] *= scale; dpsiRight[n] *= scale; } // scale for drawing on screen scale = 0; for (int n = 0; n <= nLeft; n++) if (Math.abs(psiLeft[n]) > scale) scale = Math.abs(psiLeft[n]); for (int n = 0; n <= nRight; n++) if (Math.abs(psiRight[n]) > scale) scale = Math.abs(psiRight[n]); scale *= 4; for (int n = 0; n <= nLeft; n++) { psiLeft[n] /= scale; dpsiLeft[n] /= scale; } for (int n = 0; n <= nRight; n++) { psiRight[n] /= scale; dpsiRight[n] /= scale; } misMatch = dpsiLeft[nLeft] / psiLeft[nLeft] - dpsiRight[0] / psiRight[0]; } // search for energy eigenstate class MatchingCondition implements DoubleFunctionOfDouble { // Need to make f a continuous function of energy int sign = 1; int previousSignOfPsiLeft = 1; int previousSignOfPsiRight = 1; public double f (double energy) { E = energy; integrate(); if (signOfPsiLeft != previousSignOfPsiLeft) { previousSignOfPsiLeft = -previousSignOfPsiLeft; sign = -sign; } if (signOfPsiRight != previousSignOfPsiRight) { previousSignOfPsiRight = -previousSignOfPsiRight; sign = -sign; } return misMatch * sign; } } MatchingCondition matchingCondition = new MatchingCondition(); SimpleSearch energySimpleSearch = new SimpleSearch(); void find_eigenstate () { energySimpleSearch.guessRoot(E); energySimpleSearch.setStepSize(dE); energySimpleSearch.setAccuracy(0.1 * dE); E = energySimpleSearch.findRoot(matchingCondition); EtextField.setText(CP.format(E, precision)); calibrateSlider(); } // Find energy spectrum Vector spectrum; boolean drawSpectrum; void find_spectrum () { if (spectrum == null) spectrum = new Vector(); else spectrum.removeAllElements(); // first estimate ground state energy as the zero point // energy using a harmonic oscillator approximation double x0 = Math.pow(2, 1.0/6.0); // minimum of potential double dx = 0.001; double VdoublePrime = V(x0 + dx) + V(x0 - dx) - 2 * V(x0); VdoublePrime *= gamma / (dx * dx); double m = 1; double omega = Math.sqrt(VdoublePrime / m); double hbar = 1; double E0 = 0.5 * hbar * omega; double saveE = E; E = V(x0) + 0.5 * E0; System.out.println("Spectrum of Lennard-Jones Potential"); System.out.println("-----------------------------------"); System.out.println("Quantum parameter gamma: " + gamma); System.out.println("Zero point energy estimate: " + CP.format(E0, precision)); int nLevel = 0; while (true) { find_eigenstate(); if (E < 0) { spectrum.addElement(new Double(E)); System.out.print("Level No: " + nLevel); System.out.println("\tE = " + CP.format(E, precision)); } else break; E += dE; ++nLevel; } E = saveE; } // Locate turning points of classical motion given E double leftTurningPoint; double rightTurningPoint; class TurningPoint implements DoubleFunctionOfDouble { public double f (double x) { if (x > xmax) return 0; // simulate infinite wall at this end else return E - V(x); } } SimpleSearch ss = new SimpleSearch(); TurningPoint tp = new TurningPoint(); void findTurningPoints () { double xBottom = Math.pow(2, 1.0 / 6.0); ss.guessRoot(xBottom); ss.setStepSize(-dx); ss.setAccuracy(0.1 * dx); leftTurningPoint = ss.findRoot(tp); ss.guessRoot(xBottom); ss.setStepSize(+dx); rightTurningPoint = ss.findRoot(tp); } // set left and right boundaries and match point void setLeftRightMatch () { findTurningPoints(); xLeft = leftTurningPoint - 0.1; xRight = rightTurningPoint; if (xRight < xmax) xRight += 0.5 * (xmax - rightTurningPoint); if (E < V(xmax)) xMatch = rightTurningPoint; else xMatch = 0.5 * (leftTurningPoint + rightTurningPoint); } // matching function as function of E void matchingFunction () { double de = 0.001; dx = 0.001; for (double e = -1 + de; e < de; e += de) { E = e; integrate(); System.out.println(E + "\t" + misMatch); } } class Picture extends Canvas { int pixels = 400; Picture () { setSize(pixels, pixels); setBackground(new Color(255, 250, 240)); leftPixel = topPixel = 0; rightPixel = bottomPixel = pixels; } public void paint (Graphics g) { gr = g; if (drawSpectrum) { plot_potential(); Enumeration en = spectrum.elements(); while (en.hasMoreElements()) { E = ((Double) en.nextElement()).doubleValue(); plot_energy(); integrate(); plot_psi(); } } else { plot_potential(); plot_energy(); plot_psi(); } } } Picture picture; TextField gammaTextField; TextField dxTextField; TextField EtextField; Scrollbar Eslider; TextField dEtextField; Button findButton; Button spectrumButton; int sliderValues = 100; public void init () { add(picture = new Picture()); Panel p = new Panel(); p.setLayout(new GridLayout(6, 2)); p.add(new Label("Gamma =")); p.add(gammaTextField = new TextField("" + gamma)); gammaTextField.addActionListener(this); p.add(new Label("Step Size dx =")); p.add(dxTextField = new TextField("" + dx)); dxTextField.addActionListener(this); p.add(new Label("Energy E =")); p.add(EtextField = new TextField("" + E)); EtextField.addActionListener(this); p.add(new Label("Adjust E using:")); p.add(Eslider = new Scrollbar(Scrollbar.HORIZONTAL, 10, 0, 1, sliderValues)); Eslider.addAdjustmentListener(this); p.add(new Label("in steps of dE =")); p.add(dEtextField = new TextField("" + dE)); dEtextField.addActionListener(this); p.add(findButton = new Button("Find Eigenstate")); findButton.addActionListener(this); p.add(spectrumButton = new Button("Find Spectrum")); spectrumButton.addActionListener(this); add(p); calibrateSlider(); picture.repaint(); } double Emin; double Emax; int precision = 6; void setPrecision () { precision = 1; precision += (int) (Math.abs(Math.log(dE) / Math.log(10)) + 0.5); precision += (int) (Math.abs(Math.log(Math.abs(E)) / Math.log(10)) + 0.5); } void calibrateSlider () { Emin = E - 0.5 * sliderValues * dE; Emax = E + 0.5 * sliderValues * dE; Eslider.setValue(sliderValues / 2); setPrecision(); } public void actionPerformed (ActionEvent event) { if (event.getSource() == gammaTextField) { try { gamma = new Double(event.getActionCommand()).doubleValue(); if (gamma < 0) gamma = -gamma; } catch (NumberFormatException nfe) { } gammaTextField.setText("" + gamma); } else if (event.getSource() == dxTextField) { try { dx = new Double(event.getActionCommand()).doubleValue(); } catch (NumberFormatException nfe) { dxTextField.setText("" + dx); } } else if (event.getSource() == EtextField) { try { E = new Double(event.getActionCommand()).doubleValue(); calibrateSlider(); } catch (NumberFormatException nfe) { EtextField.setText(CP.format(E, precision)); } } else if (event.getSource() == dEtextField) { try { dE = new Double(event.getActionCommand()).doubleValue(); if (dE < 0) dE = -dE; calibrateSlider(); } catch (NumberFormatException nfe) { dEtextField.setText("" + dE); } } else if (event.getSource() == findButton) { find_eigenstate(); drawSpectrum = false; } else if (event.getSource() == spectrumButton) { find_spectrum(); drawSpectrum = true; picture.repaint(); } if (!drawSpectrum) picture.repaint(); } public void adjustmentValueChanged (AdjustmentEvent event) { if (event.getSource() == Eslider) { int es = Eslider.getValue(); E = Emin + es * dE; setPrecision(); EtextField.setText("" + CP.format(E, precision)); integrate(); picture.repaint(); } } public void itemStateChanged (ItemEvent item) { picture.repaint(); } public static void main (String[] args) { Eigen2 eigen2 = new Eigen2(); CPFrame aFrame = new CPFrame("Bound State Solutions"); aFrame.add(eigen2); eigen2.init(); aFrame.setSize(700, 450); aFrame.setLocation(50, 50); aFrame.setVisible(true); if (args.length > 0) eigen2.matchingFunction(); } }