/* 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; /** Normal mode analysis */ public class NormalModeApp extends JApplet implements ActionListener { boolean isWin = false; NMode nmode; // normal mode parameters double rCutoff = 13.0; // cutoff distance for contacts int nearBy = 3; // residues apart for using nearFac double nearFac = 10000.0; // factor for scaling neighboring atoms boolean modeInit = false; // normal mode has been initialized int modeID = 7; // current mode id Timer timer; int delay = 100; int step; String pdbIDLocal; XYZCanvasProtein canvas; JPanel cpnl = new JPanel(); JTextField tURL = new JTextField("1MBN"); JTextField tRCutoff = new JTextField("" + rCutoff); JTextField tNearBy = new JTextField("" + nearBy); JTextField tNearFac = new JTextField("" + nearFac); JButton bLoad = new JButton("Compute"); JTextField tNumRes = new JTextField(""); JTextField tModeID = new JTextField("" + modeID); JButton bPlay = new JButton("Play"); JTextField tEigVal = new JTextField(""); JTextField tBallSize = new JTextField("0.5"); JTextField tStickWidth = new JTextField("4.0"); private static final long serialVersionUID = 1L; public void init() { tURL.setHorizontalAlignment(JTextField.CENTER); tRCutoff.setHorizontalAlignment(JTextField.CENTER); tNearBy.setHorizontalAlignment(JTextField.CENTER); tNearFac.setHorizontalAlignment(JTextField.CENTER); tBallSize.setHorizontalAlignment(JTextField.CENTER); tStickWidth.setHorizontalAlignment(JTextField.CENTER); tModeID.setHorizontalAlignment(JTextField.CENTER); tNumRes.setHorizontalAlignment(JTextField.RIGHT); tEigVal.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(21, 2)); cpnl.add(new JLabel(" PDB:")); tURL.addActionListener(this); cpnl.add(tURL); cpnl.add(lb = new JLabel(" Cutoff (\u212b): ")); lb.setToolTipText("Cutoff distance within which two particles are connected by a spring"); tRCutoff.addActionListener(this); cpnl.add(tRCutoff); cpnl.add(lb = new JLabel(" Nres. (nb.): ")); lb.setToolTipText("Number of successive residues to be treated as neighbors"); tNearBy.addActionListener(this); cpnl.add(tNearBy); cpnl.add(lb = new JLabel(" K (nb.) ")); lb.setToolTipText("Spring constant for neighboring atoms"); tNearFac.addActionListener(this); cpnl.add(tNearFac); bLoad.addActionListener(this); bLoad.setToolTipText("Compute all normal modes"); cpnl.add(bLoad); cpnl.add(new JLabel(" # of res.: ")); tNumRes.setEditable(false); cpnl.add(tNumRes); cpnl.add(new JLabel(" Mode:")); cpnl.add(tModeID); tModeID.addActionListener(this); cpnl.add(new JLabel(" Eigval.:")); tEigVal.setEditable(false); cpnl.add(tEigVal); bPlay.setEnabled(false); cpnl.add(bPlay); bPlay.addActionListener(this); 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(); nmode = new NMode( XYZModelProtein.getPDBInputStream(pdbID, this) ); pdbIDLocal = pdbID; } catch (Exception e) { System.out.printf("Cannot load %s, %s\n", tURL.getText(), e); } } public void start() { } public void update(Graphics g) { paint(g); } /** temperory position array */ private int np = -1; private double rp[][]; /** show current position */ public void paint(Graphics g) { if (nmode != null && nmode.n > 0) { if (modeInit && timer != null) { // start animation if (np != nmode.n) { // reallocate space np = nmode.n; rp = new double [np][3]; } // use the radius of gyration to estimate the amplitude of motion double ampMax = sqrt(nmode.n) * 3.0, amp = ampMax, freq = 1.0; // search for the first nonzero mode int i0 = 6; while (nmode.val[i0] <= 0) i0++; double eval0 = nmode.val[i0]; double eval = nmode.val[modeID - 1]; // adjust the vibrational amplitude and frequency if (modeID - 1 > i0) { amp = ampMax * sqrt(eval0/eval); freq = sqrt(eval/eval0); } amp *= sin(step * 2 * PI * freq / 100); for (int i = 0; i < np; i++) for (int j = 0; j < 3; j++) rp[i][j] = nmode.r[i][j] + nmode.vec[3*i+j][modeID-1]*amp; canvas.refresh(rp, np, nmode.res, false); step++; } else { // draw the reference structure canvas.refresh(nmode.r, nmode.n, nmode.res, false); } } cpnl.repaint(); } /** initialize normal mode from PDBID */ private void initPDB(String pdbID) { try { // open local pdb if (nmode == null || pdbID != pdbIDLocal) // initialize normal modes nmode = new NMode( XYZModelProtein.getPDBInputStream(pdbID, this) ); canvas.refresh(nmode.r, nmode.n, nmode.res, true); nmode.eig(rCutoff, nearBy, nearFac); pdbIDLocal = pdbID; modeInit = true; } catch (Exception e) { System.out.println("cannot open " + pdbID + ", error: " + e); } } DecimalFormat df = 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 == tRCutoff) { rCutoff = Double.parseDouble(tRCutoff.getText().trim()); if (rCutoff < 2.0) { rCutoff = 2.0; tRCutoff.setText(" " + rCutoff); } } if (src == tNearBy) { nearBy = Integer.parseInt(tNearBy.getText().trim()); if (nearBy < 0) { nearBy = 0; tNearBy.setText(" 0 "); } } if (src == tNearFac) { nearFac = Double.parseDouble(tNearFac.getText().trim()); if (nearFac < 0.0) { nearFac = 0.0; tNearFac.setText(" 0 "); } } if (src == bLoad || src == tURL || src == tRCutoff || src == tNearBy || src == tNearFac) { if (timer != null) { timer.stop(); timer = null; } String pdbID = tURL.getText().trim(); bPlay.setEnabled(false); initPDB(pdbID); tNumRes.setText(nmode.n + " "); bPlay.setEnabled(true); bPlay.setText("Play"); } if (src == tModeID) { modeID = Integer.parseInt(tModeID.getText().trim()); if (modeID < 1) { modeID = 1; tModeID.setText("" + modeID); } if (modeInit && modeID > nmode.n * 3) { modeID = nmode.n * 3; tModeID.setText("" + modeID); } } boolean timerStarted = (src == bPlay && timer == null) || src == tModeID; if (src == bPlay) { if (timer == null) { // currently paused bPlay.setText("Stop"); } else { // currently playing timer.stop(); timer = null; bPlay.setText("Play"); repaint(); } } if (timerStarted) { if ( !modeInit ) return; tEigVal.setText(df.format(nmode.val[modeID - 1]) + " "); if (timer != null) { timer.stop(); timer = null; } step = 0; repaint(); timer = new Timer(delay, this); timer.start(); } if (src == timer) { if ( !modeInit ) return; repaint(); } } } /** normal mode anaylsis */ class NMode { int n; // number of particles int nmax; // maximal number String res[]; // residue names double r[][]; // position 3xn double mat[][]; // matrix double val[]; // eigenvalues double vec[][]; // eigenvectors NMode() { n = 0; nmax = 1024; r = new double [nmax][3]; res = new String [nmax]; } /** constructor with PDB input stream */ NMode(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(); r[n][0] = Double.parseDouble( s.substring(31,39).trim() ); r[n][1] = Double.parseDouble( s.substring(39,47).trim() ); r[n][2] = Double.parseDouble( s.substring(47,55).trim() ); n++; if (n >= nmax) { nmax += 1024; double rn[][] = new double [nmax][3]; String resn[] = new String [nmax]; for (int i = 0; i < n; i++) { for (int j = 0; j < 3; j++) rn[i][j] = r[i][j]; resn[i] = res[i]; } r = rn; res = resn; } } System.out.println("" + n + " residues"); } double rcutoff = 13.0; double nearfac = 10000.0; int nearby = 3; /** build the Hessian matrix */ private void buildMatrix() { int i, j, d1, d2; double rji[] = new double [3], nm, fac; mat = new double[3*n][3*n]; for (i = 0; i < n - 1; i++) { for (j = i + 1; j < n; j++) { // for particle pair i and j rv3Diff(rji, r[j], r[i]); nm = rv3Norm(rji); if (nm > rcutoff) continue; rv3SMul(rji, 1.0/nm); // normalize the vector fac = (j <= i + nearby) ? nearfac : 1.0; for (d1 = 0; d1 < 3; d1++) for (d2 = 0; d2 < 3; d2++) { double rr = rji[d1] * rji[d2] * fac; mat[3*i + d1][3*i + d2] += rr; mat[3*i + d1][3*j + d2] -= rr; mat[3*j + d1][3*i + d2] -= rr; mat[3*j + d1][3*j + d2] += rr; } } } } /** Solve the eigenvalue problem * rc: cutoff distance */ public void eig(double rc, int nby, double nfac) { rcutoff = rc; nearby = nby; nearfac = nfac; buildMatrix(); Eig e = new Eig(mat, 3*n); val = e.val; vec = e.vec; for (int i = 0; i < 3*n; i++) { // print the first few eigenvalues if (i < 20 || i >= 3*n - 3) System.out.println("" + (i+1) + ": " + val[i]); if (i == 19) System.out.println("..."); } } private void rv3Diff(double z[], double x[], double y[]) { z[0] = x[0] - y[0]; z[1] = x[1] - y[1]; z[2] = x[2] - y[2]; } private double rv3Norm(double x[]) { return sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]); } private void rv3SMul(double x[], double s) { x[0] *= s; x[1] *= s; x[2] *= s; } } /** Eigenvalue and eigenvectors */ class Eig { double val[]; double vec[][]; double e[]; int n; /** constructor */ Eig(double mat[][], int n1) { n = n1; val = new double [n]; vec = new double [n][n]; e = new double [n]; for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) vec[i][j] = mat[i][j]; tridiag(vec, val, e); triqr(val, e, vec); sort(true); } /** To reduce a real symmetric matrix 'm' to a tridiagonal one by Householder transformations * The diagonal elements are saved in vector 'd' and off-diagonal elements 'e'. */ private void tridiag(double m[][], double d[], double e[]) { int i, j, k; double H, sigma, p, K, x[]; // use d[i] to indicate if the i'th Householder transformation is performed for (i = 0; i < n; i++) d[i] = 0; // n-2 Householder transformations for (i = 0; i < n-2; i++) { x = m[i]; // alias x[k] == m[i][k] for (H = 0, k = i+1; k < n; k++) H += x[k]*x[k]; sigma = x[i+1] > 0 ? sqrt(H) : -sqrt(H); // sigma = sgn(x1) |x| e[i] = -sigma; // P x = - sigma e1 H += sigma*x[i+1]; // H= (1/2) |u|^2 = |x|^2 + sigma x1 // To avoid singularity due to (partially) diagonal matrix as input if (sigma + m[i][i] == m[i][i]) { e[i] = m[i][i+1]; continue; } x[i+1] += sigma; // u = x + sigma e1, we now switch to 'u' for (j = i+1; j < n; j++) m[j][i] = x[j]/H; // save u/H in column i // CALCULATE P A P K = 0; for (j = i+1; j < n; j++) { // calculate p=A u /H, we only use the up triangle for (p = 0, k = i+1; k <= j; k++) p += m[k][j]*x[k]; for (k = j+1; k < n; k++) p += m[j][k]*x[k]; e[j] = (p /= H); // save p temporarily to e[j], notice e[i+1..n-1] are not used yet. K += x[j]*p; // K = u' p / (2H) } K /= (2*H); // K = u' p / (2H) for (j = i+1; j < n; j++) e[j] -= K*x[j]; // form q = p - K u for (j = i+1; j < n; j++) // calculate A' = A - q u' - u q' (only right-top triangle) for (k = j; k < n; k++) m[j][k] -= e[j]*x[k]+x[j]*e[k]; d[i] = 1; // indicate that the transformation is performed } e[n-2] = m[n-2][n-1]; /* for i == n-2 */ e[n-1] = 0; // if only eigenvalues are required, enable the above line and ignore the rest // To form Q = P1 ... Pn-2 d[n-2] = m[n-2][n-2]; d[n-1] = m[n-1][n-1]; // copy last two eigenvalues m[n-2][n-2] = 1; m[n-2][n-1] = 0; // initialize the right-bottom corner m[n-1][n-2] = 0; m[n-1][n-1] = 1; // P Q = (1 - u u'/H) Q = Q - (u/H) (u' Q) for (i = n-3; i >= 0; i--) { // for each P x = m[i]; // alias x[k] == m[i][k] // Form eigenvector, ONLY if i'th transformation is performed if (d[i] != 0) { for (j = i+1; j < n; j++) { // form K = u'Q for (K = 0, k = i+1; k < n; k++) K += x[k]*m[k][j]; // Q = Q - K (u/H) for (k = i+1; k < n; k++) m[k][j] -= K*m[k][i]; } } // copy the diagonal element and proceed d[i] = m[i][i]; m[i][i] = 1; for (j = i+1; j < n; j++) m[i][j] = m[j][i] = 0.; } } static final int itermax = 1000; /** diagonize a tridiagonal matrix by QR algorithm, * whose diagonal is d[0..n-1], off-diagonal is e[0..n-2]; * reduce from the left-top to right-left */ private void triqr(double d[], double e[], double mat[][]) { int i, j, k, m, iter, sgn; double ks = 0, r, c, s, delta, f = 0.0, t1, t2, tol; e[n-1] = 0; tol = 1e-12f; for (i = 0; i < n; i++) { // for each eigenvalue for (iter = 0; iter < itermax; iter++) { // Look for a single small subdiagonal element to split the matrix for (m = i; m < n-1; m++) { if (abs(e[m]) < (abs(d[m+1])+abs(d[m]))*tol) break; } // I have isolated d[i] from other matrix elements // so that d[i] is the eigenvalue. // stop iteration and look for next(i+1) eigenvalue if (m == i) break; // form shift ks delta = d[m]-d[m-1]; sgn = ((delta > 0) ? 1: -1); delta /= e[m-1]; r = hypot(delta, 1); ks = d[m] + sgn*e[m-1]/(r + abs(delta)); // Rotations for (j = i; j <= m-1; j++) { // calculate c and s if (j == i) { // First rotation r = hypot(d[i]-ks, e[i]); c = (d[i]-ks)/r; s = e[i]/r; } else { // Givens rotations r = hypot(e[j-1], f); c = e[j-1]/r; s = f/r; e[j-1] = r; } // update the diagonal and off-diagonal elements r = s*(d[j+1]-d[j]) + 2*c*e[j]; d[j] += s*r; d[j+1] -= s*r; e[j] = c*r - e[j]; f = s*e[j+1]; e[j+1] *= c; // update eigenvectors for (k = 0; k < n; k++) { t1 = mat[k][j]; t2 = mat[k][j+1]; mat[k][j] = c*t1+s*t2; mat[k][j+1] = -s*t1+c*t2; } } // end of rotations } // end for iteration }// end for each eigenvalue } /** sort eigen values and vectors into ascending order */ public void sort(boolean ac) { int i, j, im; double max, tmp; for (i = 0; i < n - 1; i++) { for (max = val[i], im = i, j = i+1; j < n; j++) { if (ac) { if (val[j] < max) max = val[im = j]; } else { if (val[j] > max) max = val[im = j]; } } if (im != i) { /* change column im and i */ tmp = val[i]; val[i] = val[im]; val[im] = tmp; for (j = 0; j < n; j++) { tmp = vec[j][i]; vec[j][i] = vec[j][im]; vec[j][im] = tmp; } } } } /** returns sqrt(x*x + y*y) */ private double hypot(double x, double y) { x = abs(x); y = abs(y); if (x > y) { y = y/x; return x * sqrt(1 + y * y); } else if (y > 0.0) { x = x/y; return y * sqrt(1 + x * x); } else return 0.0; } /** test routine */ public static final void main() { int N = 3; double [][] mat = new double [N][N]; int i, j; for (i = 0; i < N; i++) { for (j = 0; j < N; j++) mat[i][j] = 1.; mat[i][i] = i + 1; } Eig eig = new Eig(mat, N); for (i = 0; i < N; i++) { System.out.println("" + i + ": " + eig.val[i] + ";"); for (j = 0; j < N; j++) System.out.println("" + eig.vec[j][i] + " "); } } } /** 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 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); } } }