/* Elastic normal mode demonstration * Copyright (c) 2011-2014 Cheng Zhang */ 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; import java.io.*; import java.net.URL; /** Go model analysis */ public class GoApp extends JApplet implements ActionListener { boolean isWin = false; Go go; boolean sysInit = false; // model has been initialized Timer timer; int delay = 100; int speed = 100; int step = 0; XYZCanvasProtein canvas; JPanel cpnl = new JPanel(); JTextField tURL = new JTextField("1VII"); JTextField tRCutoff = new JTextField("" + Go.rc); JTextField tKb = new JTextField("" + Go.kb); JTextField tKa = new JTextField("" + Go.ka); JTextField tKd1 = new JTextField("" + Go.kd1); JTextField tKd3 = new JTextField("" + Go.kd3); JTextField tNbe = new JTextField("" + Go.nbe); JTextField tNbc = new JTextField("" + Go.nbc); JButton bLoad = new JButton("Load"); JTextField tNumRes = new JTextField(""); JButton bPlay = new JButton("Start"); JCheckBox bUseWL = new JCheckBox("Wang-Landau"); JTextField tLnf0 = new JTextField("" + WL.lnf0); JTextField tPerc = new JTextField("" + (WL.perc * 100)); JTextField tLnf = new JTextField(""); JTextField tFlatness = new JTextField(""); JTextField tSpeed = new JTextField("" + speed); JTextField tTemp = new JTextField("" + Go.getT()); JTextField tBallSize = new JTextField(" 0.5"); JTextField tStickWidth = new JTextField(" 4.0"); JTextField tEpot = new JTextField(""); JTextField tEkin = new JTextField(""); private static final long serialVersionUID = 1L; public void init() { tURL.setHorizontalAlignment(JTextField.CENTER); tRCutoff.setHorizontalAlignment(JTextField.CENTER); tKb.setHorizontalAlignment(JTextField.CENTER); tKa.setHorizontalAlignment(JTextField.CENTER); tKd1.setHorizontalAlignment(JTextField.CENTER); tKd3.setHorizontalAlignment(JTextField.CENTER); tNbe.setHorizontalAlignment(JTextField.CENTER); tNbc.setHorizontalAlignment(JTextField.CENTER); tNumRes.setHorizontalAlignment(JTextField.CENTER); tLnf0.setHorizontalAlignment(JTextField.CENTER); tPerc.setHorizontalAlignment(JTextField.CENTER); tLnf.setHorizontalAlignment(JTextField.CENTER); tFlatness.setHorizontalAlignment(JTextField.CENTER); tSpeed.setHorizontalAlignment(JTextField.CENTER); tTemp.setHorizontalAlignment(JTextField.CENTER); tBallSize.setHorizontalAlignment(JTextField.CENTER); tStickWidth.setHorizontalAlignment(JTextField.CENTER); tEpot.setHorizontalAlignment(JTextField.RIGHT); tEkin.setHorizontalAlignment(JTextField.RIGHT); JLabel lb; Container box = getContentPane(); box.setLayout(new BorderLayout()); canvas = new XYZCanvasProtein(); box.add(canvas, BorderLayout.CENTER); box.add(cpnl, BorderLayout.EAST); cpnl.setLayout(new GridLayout(22, 2)); cpnl.add(new JLabel(" PDB:")); tURL.addActionListener(this); cpnl.add(tURL); bLoad.addActionListener(this); bLoad.setToolTipText("Load PDB structure"); cpnl.add(bLoad); bPlay.addActionListener(this); bPlay.setToolTipText("Start MD simulation"); cpnl.add(bPlay); cpnl.add(new JLabel(" # of res.: ")); tNumRes.setEditable(false); cpnl.add(tNumRes); cpnl.add(lb = new JLabel(" Cutoff (\u212b): ")); lb.setToolTipText("Cutoff distance within which two particles are connected"); tRCutoff.addActionListener(this); cpnl.add(tRCutoff); cpnl.add(lb = new JLabel(" kb ")); lb.setToolTipText("Spring constant for harmonic bonds"); tKb.addActionListener(this); cpnl.add(tKb); cpnl.add(lb = new JLabel(" ka ")); lb.setToolTipText("Spring constant for harmonic bond angles"); tKa.addActionListener(this); cpnl.add(tKa); cpnl.add(lb = new JLabel(" kd1 ")); lb.setToolTipText("Spring constant for dihedral angles"); tKd1.addActionListener(this); cpnl.add(tKd1); cpnl.add(lb = new JLabel(" kd3 ")); lb.setToolTipText("Spring constant for dihedral angles"); tKd3.addActionListener(this); cpnl.add(tKd3); cpnl.add(lb = new JLabel(" \u03b5 ")); lb.setToolTipText("eplison for nonboned interaction"); tNbe.addActionListener(this); cpnl.add(tNbe); cpnl.add(lb = new JLabel(" rp (\u212b) ")); lb.setToolTipText("radius of repulsion"); tNbc.addActionListener(this); cpnl.add(tNbc); cpnl.add(lb = new JLabel(" Speed: ")); tSpeed.addActionListener(this); cpnl.add(tSpeed); cpnl.add(lb = new JLabel(" T (K): ")); tTemp.addActionListener(this); cpnl.add(tTemp); cpnl.add(lb = new JLabel(" Epot. ")); tEpot.setEditable(false); cpnl.add(tEpot); cpnl.add(lb = new JLabel(" Tkin (K) ")); tEkin.setEditable(false); cpnl.add(tEkin); // Wang-Landau stuff bUseWL.addActionListener(this); cpnl.add(bUseWL); cpnl.add(new JLabel("")); cpnl.add(lb = new JLabel(" lnf0:")); tLnf0.addActionListener(this); cpnl.add(tLnf0); cpnl.add(lb = new JLabel(" H. cutoff (%):")); tPerc.addActionListener(this); cpnl.add(tPerc); cpnl.add(lb = new JLabel(" lnf:")); tLnf.setEditable(false); cpnl.add(tLnf); cpnl.add(lb = new JLabel(" H. flatness:")); tFlatness.setEditable(false); cpnl.add(tFlatness); cpnl.add(new JLabel(" Ball radius: ")); tBallSize.addActionListener(this); cpnl.add(tBallSize); cpnl.add(new JLabel(" Stick width: ")); tStickWidth.addActionListener(this); cpnl.add(tStickWidth); try { // try to load the initial structure, but don't compute modes String pdbID = tURL.getText(); go = new Go( XYZModelProtein.getPDBInputStream(pdbID, this) ); } catch (Exception e) { System.out.printf("Cannot load %s, %s\n", tURL, e); } } /** Note: for application, this function is not called */ public void start() { } public void update(Graphics g) { paint(g); } /** show current position */ public void paint(Graphics g) { if (go != null && go.n > 0) { if (sysInit) { step++; canvas.refresh(go.x, go.n, go.res, false); } else { canvas.refresh(go.xref, go.n, go.res, false); } } cpnl.repaint(); } /** initialize normal mode from PDBID */ private void initPDB(String pdbID) { try { // open local pdb go = new Go( XYZModelProtein.getPDBInputStream(pdbID, this) ); go.initSys(); canvas.refresh(go.x, go.n, go.res, false); sysInit = true; } catch (Exception e) { System.out.println("cannot open " + pdbID + ", error: " + e); } } DecimalFormat df = new DecimalFormat("####.###"), dfsci = new DecimalFormat("0.0E0"), dfper = new DecimalFormat("##.##%"); /** handle actions */ public void actionPerformed(ActionEvent e) { Object src = e.getSource(); if (src == tBallSize) { float ballSize = Float.parseFloat(tBallSize.getText().trim()); if (ballSize < 0.0f) ballSize = 0.0f; XYZModelProtein m = (XYZModelProtein) canvas.model; m.ballSize = ballSize; repaint(); } if (src == tStickWidth) { float stickWidth = Float.parseFloat(tStickWidth.getText().trim()); if (stickWidth < 0.0) stickWidth = 0.0f; XYZModelProtein m = (XYZModelProtein) canvas.model; m.stickWidth = stickWidth; repaint(); } if (src == tTemp) { double temp = Double.parseDouble(tTemp.getText().trim()); if (temp < 0.0) { temp = 0.0; tTemp.setText(" " + temp); } Go.setT(temp); } if (src == tSpeed) { speed = Integer.parseInt(tSpeed.getText().trim()); } boolean reset = false; // reset the model if (src == tRCutoff) { double rCutoff = Double.parseDouble(tRCutoff.getText().trim()); if (rCutoff < 6.0) { rCutoff = 6.0; tRCutoff.setText("" + rCutoff); } Go.rc = rCutoff; reset = true; } if (src == tKb) { double kb = Double.parseDouble(tKb.getText().trim()); if (kb < 0.0) { kb = 0.0; tKb.setText("0"); } Go.kb = kb; reset = true; } if (src == tKa) { double ka = Double.parseDouble(tKa.getText().trim()); if (ka < 0.0) { ka = 0.0; tKa.setText("0"); } Go.ka = ka; reset = true; } if (src == tKd1) { double kd1 = Double.parseDouble(tKd1.getText().trim()); if (kd1 < 0.0) { kd1 = 0.0; tKd1.setText(" 0 "); } Go.kd1 = kd1; reset = true; } if (src == tKd3) { double kd3 = Double.parseDouble(tKd3.getText().trim()); if (kd3 < 0.0) { kd3 = 0.0; tKd3.setText(" 0 "); } Go.kd3 = kd3; reset = true; } if (src == tNbe) { double nbe = Double.parseDouble(tNbe.getText().trim()); if (nbe < 0.0) { nbe = 0.0; tNbe.setText(" 0 "); } Go.nbe = nbe; reset = true; } if (src == tNbc) { double nbc = Double.parseDouble(tNbc.getText().trim()); if (nbc < 0.0) { nbc = 0.0; tNbc.setText(" 0 "); } Go.nbc = nbc; reset = true; } if (src == bLoad || src == tURL || reset) { if (timer != null) { timer.stop(); timer = null; } bPlay.setText("Start"); bPlay.setEnabled(false); sysInit = false; String pdbID = tURL.getText().trim(); initPDB(pdbID); tNumRes.setText("" + go.n); bPlay.setEnabled(true); } // Wang Landau stuff boolean resetWL = false; if (src == bUseWL) { go.useWL = !go.useWL; if (go.useWL) resetWL = true; } if (src == tPerc) { double p = Double.parseDouble(tPerc.getText().trim()); p *= 0.01; if (p < 1e-4) { p = 1e-4; tPerc.setText(" " + dfper.format(p)); } WL.perc = p; } if (src == tLnf0) { double lnf0 = Double.parseDouble(tLnf0.getText().trim()); WL.lnf0 = lnf0; if (go.useWL) resetWL = true; } if (resetWL) { go.initWL(); } boolean startTimer = (src == bPlay && timer == null); if (src == bPlay) { if (timer == null) { // currently paused bPlay.setText("Stop"); } else { // currently playing timer.stop(); timer = null; bPlay.setText("Start"); repaint(); } } if (startTimer) { if (!sysInit) return; if (timer != null) { timer.stop(); timer = null; } step = 0; repaint(); timer = new Timer(delay, this); timer.start(); } // scheduled MD step if (src == timer) { if (!sysInit) return; for (int i = 0; i < speed; i++) go.vv(); tEpot.setText(df.format(go.epot) + " "); tEkin.setText(df.format(2*go.ekin/go.dof/Go.kB) + " "); if (go.useWL) { tLnf.setText(dfsci.format(go.wl.lnf) + " "); tFlatness.setText(dfper.format(go.wl.flatness) + " "); } repaint(); } } } /** Go model */ class Go { // Go model parameters static double kb = 200.0; // .5 kb (b - b0)^2 static double ka = 40; // .5 ka (a - a0)^2 static double kd1 = 1, kd3 = .5; static double nbe = 1.0, nbc = 4.0; static double rc = 10.0; // cutoff distance for contacts // other parameters double mddt = 2e-3; // MD step size double thermdt = 1e-2; // thermostat step int n; // number of particles int nmax; // maximal number int dof; // degrees of freedom String res[]; // residue names double xref[][]; // reference position double x[][]; // position 3xn double v[][]; // velocity double f[][]; // force double epot, ekin; double epotRef; static final double mass = 1.0; static final double kB = 0.0083145 / 4.184; // Boltzmann constant static double kT = kB * 300; // kB T boolean useWL = false; // use Wang-Landau to compute free energy WL wl; Random rng = new Random(); Go() { n = dof = 0; nmax = 1024; res = new String [nmax]; xref = new double [nmax][3]; } /** constructor with PDB input stream */ Go(InputStream is) throws Exception { this(); BufferedReader br = new BufferedReader(new InputStreamReader(is)); String s; while ((s = br.readLine()) != null) { if ( s.startsWith("ENDMDL") || s.startsWith("TER") ) break; // single chain if ( !s.startsWith("ATOM") ) continue; String atnm = s.substring(12,16).trim(); if (!atnm.equals("CA")) continue; res[n] = s.substring(17, 20).trim(); xref[n][0] = Double.parseDouble( s.substring(31,39).trim() ); xref[n][1] = Double.parseDouble( s.substring(39,47).trim() ); xref[n][2] = Double.parseDouble( s.substring(47,55).trim() ); n++; if (n >= nmax) { nmax += 1024; double xn[][] = new double [nmax][3]; String resn[] = new String [nmax]; for (int i = 0; i < n; i++) { for (int j = 0; j < 3; j++) xref[i][j] = xn[i][j]; resn[i] = res[i]; } xref = xn; res = resn; } } System.out.println("" + n + " residues"); dof = 3 * n - 6; } double bref[]; // distance between i, i + 1 double aref[]; // distance between i, i + 2 double dref[]; // distance between i, i + 3 double r2ref[][]; boolean iscont[][]; /** initialize potential energy */ void initPot() { int i, j; bref = new double [n - 1]; for (i = 0; i < n - 1; i++) bref[i] = rv3.dist(xref[i], xref[i+1]); aref = new double [n - 2]; for (i = 0; i < n - 2; i++) aref[i] = rv3.ang(xref[i], xref[i+1], xref[i+2], null, null, null); dref = new double [n - 3]; for (i = 0; i < n - 3; i++) dref[i] = rv3.dih(xref[i], xref[i+1], xref[i+2], xref[i+3], null, null, null, null); r2ref = new double [n][n]; iscont = new boolean [n][n]; for (i = 0; i < n - 1; i++) { for (j = i + 1; j < n; j++) { double dr = rv3.dist(xref[i], xref[j]); r2ref[i][j] = r2ref[j][i] = dr * dr; if (j >= i + 4 && dr < rc) { iscont[i][j] = iscont[j][i] = true; System.out.println("contact " + i + " and " + j); } } r2ref[i][i] = 0.; } } /** initialize an MD system */ void initMD() { x = new double [n][3]; v = new double [n][3]; f = new double [n][3]; for (int i = 0; i < n; i++) { for (int d = 0; d < 3; d++) { x[i][d] = xref[i][d] + (rng.nextDouble() - 0.5) * .2; // xref + noise v[i][d] = 0; f[i][d] = 0; } } epotRef = force(f, xref); epot = force(f, x); ekin = ekinet(v); } void initWL() { double emin, emax; if (epotRef < 0.0) { emin = epotRef * 0.9; emax = -emin; } else { // a very bad set up emin = epotRef; emax = epotRef * 3; } wl = new WL(emin, emax); } void initSys() { initPot(); initMD(); initWL(); } double dx[] = new double [3]; /** compute go force */ private double force(double f[][], double x[][]) { double dr2, dr4, dr6, dr10, invr2, amp, nbc2 = nbc * nbc, ene = 0.0; int i, j; for (i = 0; i < n; i++) rv3.zero(f[i]); // bonded for (i = 0; i < n - 1; i++) { ene += harmonic(x[i], x[i+1], bref[i], kb, f[i], f[i+1]); } // angles for (i = 0; i < n - 2; i++) { ene += angle(x[i], x[i+1], x[i+2], aref[i], kb, f[i], f[i+1], f[i+2]); } // dihedrals for (i = 0; i < n - 3; i++) { ene += dihedral(x[i], x[i+1], x[i+2], x[i+3], dref[i], kd1, kd3, f[i], f[i+1], f[i+2], f[i+3]); } // nonbonded for (i = 0; i < n - 4; i++) { for (j = i + 4; j < n; j++) { rv3.diff(dx, x[i], x[j]); dr2 = rv3.sqr(dx); invr2 = 1.0/dr2; if (iscont[i][j]) { dr2 = r2ref[i][j] * invr2; dr4 = dr2 * dr2; dr6 = dr4 * dr2; dr10 = dr4 * dr6; amp = nbe * 60 * (dr2 - 1) * dr10 * invr2; ene += nbe * (5 * dr2 - 6) * dr10; } else { dr2 = nbc2 / dr2; dr6 = dr2 * dr2 * dr2; dr6 *= dr6; amp = nbe * 12 * dr6 * invr2; ene += nbe * dr6; } rv3.sinc(f[i], dx, amp); rv3.sinc(f[j], dx, -amp); } } return ene; } double gi[] = new double [3]; double gj[] = new double [3]; double gk[] = new double [3]; double gl[] = new double [3]; /** harmonic force */ private double harmonic(double xi[], double xj[], double ref, double k, double fi[], double fj[]) { rv3.diff(dx, xi, xj); double dr = rv3.norm(dx); double del = dr - ref; double amp = k * del/dr; rv3.sinc(fi, dx, -amp); rv3.sinc(fj, dx, +amp); return .5 * k * del * del; } /** bond-angle force */ private double angle(double xi[], double xj[], double xk[], double ref, double k, double fi[], double fj[], double fk[]) { double ang = rv3.ang(xi, xj, xk, gi, gj, gk); double del = ang - ref; double amp = - k * del; rv3.sinc(fi, gi, amp); rv3.sinc(fj, gj, amp); rv3.sinc(fk, gk, amp); return .5 * k * del * del; } /** dihedral-angle force */ private double dihedral(double xi[], double xj[], double xk[], double xl[], double ref, double k1, double k3, double fi[], double fj[], double fk[], double fl[]) { double dih = rv3.dih(xi, xj, xk, xl, gi, gj, gk, gl); double del = dih - ref; if (del > PI) del -= 2*PI; else if (del < -PI) del += 2*PI; double amp = -k1 * sin(del) - k3 * 3 * sin(3*del); rv3.sinc(fi, gi, amp); rv3.sinc(fj, gj, amp); rv3.sinc(fk, gk, amp); rv3.sinc(fl, gl, amp); return k1 * (1 - cos(del)) + k3 * (1 - cos(3*del)); } double fscal = 1.0; /** velocity scaling */ void vv() { for (int i = 0; i < n; i++) { rv3.sinc(v[i], f[i], fscal * .5 * mddt/mass); rv3.sinc(x[i], v[i], mddt); } epot = force(f, x); fscal = (wl != null && useWL) ? wl.add(epot) * kT : 1.0; for (int i = 0; i < n; i++) { rv3.sinc(v[i], f[i], fscal *.5 * mddt/mass); } rmcom(v); ekin = vrescale(v); } /** return the kinetic energy */ private double ekinet(double v[][]) { double ene = 0; for (int i = 0; i < n; i++) ene += .5 * mass * rv3.sqr(v[i]); return ene; } /** remove motion of the center of mass */ private void rmcom(double v[][]) { for (int d = 0; d < 3; d++) { // remove the center of mass motion double vc = 0.0; for (int id = 0; id < n; id++) vc += v[id][d]; vc /= n; for (int id = 0; id < n; id++) v[id][d] -= vc; } } /** velocity rescaling thermostat */ private double vrescale(double v[][]) { double ek1 = ekinet(v); double ekav = .5f*kT*dof; double amp = 2*sqrt(ek1*ekav*thermdt/dof); double ek2 = ek1 + (ekav - ek1)*thermdt + amp*rng.nextGaussian(); if (ek2 < 0) ek2 = ek1; double s = sqrt(ek2/ek1); for (int i = 0; i < n; i++) for (int d = 0; d < 3; d++) v[i][d] *= s; return ek2; } /** set the current temperature */ static void setT(double T) { kT = kB * T; } /** get the current temperature */ static double getT() { return kT / kB; } } /** 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 XYZCanvasProtein extends XYZCanvas { public static final long serialVersionUID = 4L; /** Refresh coordinates * x[][] is the wrapped coordinates * `n' may be less than x.length */ public void refresh(double x[][], int n, String res[], boolean adjScale) { if (model == null) { model = new XYZModelProtein(); adjScale = true; } XYZModelProtein modelPro = (XYZModelProtein) model; modelPro.updateXYZ(x, n, true); modelPro.updateRes(res, n); if ( adjScale ) realSpan = model.getSpan(x, n); repaint(); } void loadPDB(String pdb, Applet applet) { XYZModelProtein modelPro = new XYZModelProtein(); if ((realSpan = modelPro.loadPDB(pdb, applet)) > 0) model = modelPro; } } /** 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 proteins */ class XYZModelProtein extends XYZModel { boolean useStick = true; double stickWidth = 4.0f; // line width static Hashtable resTable = new Hashtable(); static { resTable.put("GLY", new Atom(1.0, 1.0, 1.0, 0.9*2)); resTable.put("ALA", new Atom(0.8, 1.0, 1.0, 1.0*2)); resTable.put("VAL", new Atom(0.5, 0.9, 0.9, 1.3*2)); resTable.put("LEU", new Atom(0.6, 0.8, 0.8, 1.4*2)); resTable.put("ILE", new Atom(0.6, 0.8, 0.8, 1.4*2, 0.7, 0.4, 2.0)); resTable.put("PRO", new Atom(0.6, 0.8, 0.8, 1.3*2, 0.7, 0.4, 2.0)); resTable.put("SER", new Atom(0.9, 0.7, 1.0, 1.2*2, 0.7, 0.4, 2.0)); resTable.put("THR", new Atom(0.7, 0.6, 0.8, 1.3*2, 0.7, 0.4, 2.0)); resTable.put("CYS", new Atom(1.0, 1.0, 0.0, 1.2*2, 0.7, 0.4, 2.0)); resTable.put("MET", new Atom(0.6, 0.6, 0.0, 1.4*2, 0.7, 0.4, 2.0)); resTable.put("ASN", new Atom(0.2, 1.0, 0.4, 1.4*2, 0.7, 0.4, 2.0)); resTable.put("GLN", new Atom(0.0, 0.8, 0.2, 1.5*2, 0.7, 0.4, 2.0)); resTable.put("ASP", new Atom(1.0, 0.0, 0.0, 1.4*2, 0.7, 0.4, 2.0)); resTable.put("GLU", new Atom(1.0, 0.2, 0.0, 1.5*2, 0.7, 0.4, 2.0)); resTable.put("LYS", new Atom(0.0, 0.4, 1.0, 1.5*2, 0.7, 0.4, 2.0)); resTable.put("ARG", new Atom(0.0, 0.0, 1.0, 1.7*2, 0.7, 0.4, 2.0)); resTable.put("HIS", new Atom(0.0, 0.4, 0.8, 1.5*2, 0.7, 0.4, 2.0)); resTable.put("PHE", new Atom(0.5, 0.7, 0.7, 1.5*2, 0.7, 0.4, 2.0)); resTable.put("TYR", new Atom(0.4, 0.6, 0.6, 1.8*2, 0.7, 0.4, 2.0)); resTable.put("TRP", new Atom(0.3, 0.5, 0.5, 2.0*2, 0.7, 0.4, 2.0)); } static Hashtable atomTable = new Hashtable(); { // Atom(red, green, blue, relRadius, rContrast, spotlightMag, zContrast atomTable.put("H", new Atom(0.7, 0.7, 0.7, 1.00, 0.7, 0.4, 2.0)); atomTable.put("C", new Atom(0.0, 0.7, 0.7, 1.70, 0.7, 0.4, 2.0)); atomTable.put("N", new Atom(0.0, 0.0, 0.9, 1.55, 0.7, 0.4, 2.0)); atomTable.put("O", new Atom(1.0, 0, 0, 1.52, 0.7, 0.4, 2.0)); atomTable.put("S", new Atom(0.7, 0.7, 0, 1.80, 0.7, 0.4, 2.0)); atomTable.put("F", new Atom(0.7, 0.1, 0.7, 1.47, 0.7, 0.4, 2.0)); atomTable.put("P", new Atom(0.7, 0, 0.7, 1.80, 0.7, 0.4, 2.0)); atomTable.put("CL", new Atom(0.7, 0.7, 0.4, 1.75, 0.7, 0.4, 2.0)); atomTable.put("CU", new Atom(0.6, 0.3, 0, 1.40, 0.7, 0.4, 2.0)); } Atom atomGrey = new Atom(0.5, 0.5, 0.5); XYZModelProtein() { atomDefault = atomGrey; } static InputStream getPDBInputStream(String pdb, Applet applet) { pdb = pdb.trim().toUpperCase() + ".pdb"; InputStream is = null; URL base = null; try { // web mode base = applet.getCodeBase(); } catch (NullPointerException e1) { // window application // "user.dir" means the current working directory File f = new File(System.getProperty("user.dir")); try { base = f.toURI().toURL(); } catch (Exception e) {} } // try to load the local file try { is = new URL(base, pdb).openStream(); if (is != null) return is; } catch (Exception e) { } // try to download PDB from the internet try { is = new URL("http://www.rcsb.org/pdb/files/" + pdb).openStream(); } catch (Exception e) { System.out.println("Failed to download " + pdb); } return is; } /** Parse PDB from the internet * load all atoms * return the span */ double parsePDB(InputStream is) { BufferedReader br = new BufferedReader(new InputStreamReader(is)); String s; maxVert = np = 0; try { while ((s = br.readLine()) != null) { if ( s.startsWith("TER") || s.startsWith("ENDMDL") ) continue; if ( !s.startsWith("ATOM") ) continue; String atnm = s.substring(12, 16).trim(); String atomnm = Character.isDigit(atnm.charAt(0)) ? atnm.substring(1, 2) : atnm.substring(0, 1); double x = Double.parseDouble( s.substring(31, 39).trim() ); double y = Double.parseDouble( s.substring(39, 47).trim() ); double z = Double.parseDouble( s.substring(47, 55).trim() ); addVert(atomnm, x, y, z); } } catch (IOException e) { e.printStackTrace(); } updateXYZ(realXYZ, np, true); useStick = false; // turn off sticks return getSpan(realXYZ, np); } double loadPDB(String pdb, Applet applet) { InputStream is = getPDBInputStream(pdb, applet); double realSpan = 0; if (is != null) { realSpan = parsePDB(is); try { is.close(); } catch (Exception e){} } return realSpan; } int maxVert = 0; /** Add a vertex to this model */ int addVert(String name, double x, double y, double z) { int i = np; if (i >= maxVert) { if (realXYZ == null) { maxVert = 100; realXYZ = new double[maxVert][3]; atoms = new Atom[maxVert]; screenXYZ = new int [maxVert][3]; zOrder = new int [maxVert]; } else { maxVert *= 2; double arr[][] = new double[maxVert][3]; System.arraycopy(realXYZ, 0, arr, 0, realXYZ.length); realXYZ = arr; Atom na[] = new Atom[maxVert]; System.arraycopy(atoms, 0, na, 0, atoms.length); atoms = na; screenXYZ = new int [maxVert][3]; zOrder = new int [maxVert]; } } Atom a = atomTable.get(name.toUpperCase()); atoms[i] = (a == null) ? atomDefault : a; realXYZ[i][0] = x; realXYZ[i][1] = y; realXYZ[i][2] = z; return np++; } void updateRes(String res[], int n) { for (int i = 0; i < n; i++) { Atom a = resTable.get(res[i].toUpperCase()); atoms[i] = (a == null) ? atomDefault : a; } } /** Return i, such that zOrder[i] == j */ private int zOrderWho(int j) { int i; for (i = 0; i < np; i++) if (zOrder[i] == j) return i; return -1; } /** Draw a stick */ void drawStick(Graphics g, int iz) { Graphics2D g2 = (Graphics2D) g; // let farther sticks darker int grey = 50 + 150 * iz/np; g.setColor(new Color(grey, grey, grey)); int i = zOrder[iz]; if ( i < np - 1 && zOrderWho(i + 1) > iz ) // i -- i + 1 g2.drawLine(screenXYZ[i][0], screenXYZ[i][1], screenXYZ[i + 1][0], screenXYZ[i + 1][1]); if ( i > 0 && zOrderWho(i - 1) > iz ) // i -- i - 1 g2.drawLine(screenXYZ[i][0], screenXYZ[i][1], screenXYZ[i - 1][0], screenXYZ[i - 1][1]); } /** Paint this model to a graphics context. * It uses the matrix associated with this model * to map from model space to screen space. */ void paint(Graphics g) { if (realXYZ == null || np <= 0) return; getZOrder(); if ( useStick ) { // preparing a stroke for drawing sticks BasicStroke stroke = new BasicStroke( (float) (.1f * real2Screen * stickWidth), BasicStroke.CAP_ROUND, BasicStroke.JOIN_ROUND); Graphics2D g2 = (Graphics2D) g; g2.setStroke(stroke); } for (int iz = 0; iz < np; iz++) { drawAtom(g, iz); // draw the atom first if ( useStick ) drawStick(g, iz); } } } /** 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 WL { double emin, emax; double de; int ecnt = 40; // number of energy bins double hist[]; // histogram double lng[]; // estimated entropy function double lnf = lnf0; int step = 0; double flatness; double epot; // last epot; double beta; // last beta; static double lnf0 = 0.1; static double lnffac = 0.5; static double perc = 0.2; static int nstcheck = 10000; /** constructor */ WL(double e0, double e1) { emin = e0; emax = e1; de = (emax - emin)/ecnt; hist = new double [ecnt + 1]; lng = new double [ecnt + 1]; lnf = lnf0; step = 0; flatness = 1.0; } /** compute histogram flatness */ double histogramFlatness() { double hmin = 1e9, hmax = -1e9; for (int i = 0; i < ecnt; i++) { if (hist[i] > hmax) hmax = hist[i]; else if (hist[i] < hmin) hmin = hist[i]; } if (hmax > 0.0) return (hmax - hmin) / (hmax + hmin); else return 1.0; } DecimalFormat df = new DecimalFormat("####.00"); DecimalFormat dfh = new DecimalFormat("######.0"); /** dump histogram information etc */ void dump() { for (int i = 0; i <= ecnt; i++) System.out.println("E " + df.format(emin + (i+.5)*de) + " hist " + dfh.format(hist[i]) + " S " + df.format(lng[i])); System.out.println("flatness " + flatness + " lnf " + lnf + " epot " + epot + " beta " + beta); } /** check histogram flatness, switch stage if possible */ int check() { for (int i = ecnt; i >= 0; i--) lng[i] -= lng[0]; // let lng start from 0 flatness = histogramFlatness(); dump(); if (flatness < perc) { System.out.println("switch stages step " + step); lnf *= lnffac; for (int i = 0; i < ecnt; i++) hist[i] = 0.0; // clear histogram return 1; } return 0; } // beta for 300K is 1.6774 static double betmin = -1.0; static double betmax = 10.0; private double limit(double b) { if (b < betmin) return betmin; else if (b > betmax) return betmax; else return b; } /** add energy to histogram, update lng, * return mean force */ double add(double e) { epot = e; // register e int ies, ie = -1, ieb; // index of updating lng, histogram, and computing beta if (e < emin) { ies = 0; ieb = 0; } else if (e >= emax) { ies = ecnt; ieb = ecnt - 1; ie = ecnt; } else { ies = (int)((e - emin)/de + .5); ieb = ie = (int)((e - emin)/de); } lng[ies] += lnf; if (ie > 0) hist[ie] += 1.0; step++; if (step % nstcheck == 0) check(); beta = limit( (lng[ieb+1] - lng[ieb])/de ); return beta; } } class rv3 { static void zero(double x[]) { x[0] = x[1] = x[2] = 0.0; } static double[] diff(double z[], double x[], double y[]) { z[0] = x[0] - y[0]; z[1] = x[1] - y[1]; z[2] = x[2] - y[2]; return z; } static void inc(double x[], double dx[]) { x[0] += dx[0]; x[1] += dx[1]; x[2] += dx[2]; } static void sinc(double x[], double dx[], double s) { x[0] += s * dx[0]; x[1] += s * dx[1]; x[2] += s * dx[2]; } static double dot(double x[], double y[]) { return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; } static double dx[] = new double [3]; static double dist2(double x[], double y[]) { diff(dx, x, y); return sqr(dx); } static double dist(double x[], double y[]) { return sqrt(dist2(x, y)); } static double sqr(double x[]) { return x[0]*x[0] + x[1]*x[1] + x[2]*x[2]; } static double norm(double x[]) { return sqrt(sqr(x)); } static void smul(double x[], double s) { x[0] *= s; x[1] *= s; x[2] *= s; } static void smul2(double x[], double y[], double s) { x[0] = y[0]*s; x[1] = y[1]*s; x[2] = y[2]*s; } static double[] lincomb2(double z[], double x[], double y[], double a, double b) { z[0] = x[0]*a + y[0]*b; z[1] = x[1]*a + y[1]*b; z[2] = x[2]*a + y[2]*b; return z; } static double [] cross(double z[], double x[], double y[]) { z[2] = x[0]*y[1] - x[1]*y[0]; z[0] = x[1]*y[2] - x[2]*y[1]; z[1] = x[2]*y[0] - x[0]*y[2]; return z; } /* temporary vectors */ static double x12[] = new double [3]; static double x32[] = new double [3]; static double x34[] = new double [3]; static double m[] = new double [3]; static double n[] = new double [3]; static double uvec[] = new double [3]; static double vvec[] = new double [3]; static double svec[] = new double [3]; /** angle and gradients of cos(x1-x2-x3) */ static double cosang(double x1[], double x2[], double x3[], double g1[], double g2[], double g3[]) { double ra, rb, dot; ra = norm(diff(x12, x1, x2)); smul(x12, 1.f/ra); rb = norm(diff(x32, x3, x2)); smul(x32, 1.f/rb); dot = rv3.dot(x12, x32); if (dot > 1) dot = 1; else if (dot < -1) dot = -1; if (g1 != null && g2 != null && g3 != null) { lincomb2(g1, x32, x12, 1.f/ra, -dot/ra); lincomb2(g3, x12, x32, 1.f/rb, -dot/rb); lincomb2(g2, g1, g3, -1.0, -1.0); } return dot; } /** angle and gradients of x1-x2-x3 */ static double ang(double x1[], double x2[], double x3[], double g1[], double g2[], double g3[]) { double dot, sn; dot = cosang(x1, x2, x3, g1, g2, g3); sn = sqrt(1 - dot*dot); if (sn < 1e-7) sn = 1; else sn = -1.f/sn; if (g1 != null && g2 != null && g3 != null) { smul(g1, sn); smul(g2, sn); smul(g3, sn); } return acos(dot); } /** dihedral */ static double dih(double x1[], double x2[], double x3[], double x4[], double g1[], double g2[], double g3[], double g4[]) { diff(x12, x1, x2); diff(x32, x3, x2); diff(x34, x3, x4); double r32sqr = sqr(x32); double r32 = sqrt(r32sqr); double tol = r32sqr * 1e-16f; cross(m, x12, x32); double m2 = sqr(m); cross(n, x32, x34); double n2 = sqr(n); double cosphi = 1; if (m2 > tol && n2 > tol) { cosphi = dot(m, n) / sqrt(m2 * n2); if (cosphi > 1.0) cosphi = 1.0; else if (cosphi < -1.0) cosphi = -1.0; } double phi = acos(cosphi); double vol = dot(n, x12); double sgn = (vol > 0.0) ? 1.0 : -1.0; phi *= sgn; if (g1 != null && g2 != null && g3 != null && g4 != null) { smul2(g1, m, r32/m2); smul2(g4, n, -r32/n2); double p = dot(x12, x32)/r32sqr; smul2(uvec, g1, p); double q = dot(x34, x32)/r32sqr; smul2(vvec, g4, q); diff(svec, uvec, vvec); lincomb2(g2, svec, g1, 1.0, -1.0); lincomb2(g3, svec, g4, -1.0, -1.0); } return phi; } } 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); } } }