/** * lterm.java - MM2-style length energy term * Copyright (c) 1997 Will Ware, all rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and other materials provided with the distribution. * 3. All advertising materials mentioning features or use of this software * or its derived works must display the following acknowledgement: * This product includes software developed by Will Ware. * * This software is provided "as is" and any express or implied warranties, * including, but not limited to, the implied warranties of merchantability * or fitness for any particular purpose are disclaimed. In no event shall * Will Ware be liable for any direct, indirect, incidental, special, * exemplary, or consequential damages (including, but not limited to, * procurement of substitute goods or services; loss of use, data, or * profits; or business interruption) however caused and on any theory of * liability, whether in contract, strict liability, or tort (including * negligence or otherwise) arising in any way out of the use of this * software, even if advised of the possibility of such damage. */ import java.lang.Math; import java.util.Vector; import atom; public class lterm extends term { public static final String rcsid = "$Id: lterm.java,v 1.8 1997/08/13 02:51:11 wware Exp $"; private double ks, r0; public lterm (atom a1, atom a2) { int i; boolean found; myAtoms = new atom[2]; myAtoms[0] = a1; myAtoms[1] = a2; for (i = 0, found = false; i < lengthCoeffs.length && !found; i++) if ((a1.atomicNumber () == lengthCoeffs[i][0] && a1.hybridization == lengthCoeffs[i][1] && a2.atomicNumber () == lengthCoeffs[i][2] && a2.hybridization == lengthCoeffs[i][3]) || (a1.atomicNumber () == lengthCoeffs[i][2] && a1.hybridization == lengthCoeffs[i][3] && a2.atomicNumber () == lengthCoeffs[i][0] && a2.hybridization == lengthCoeffs[i][1])) { found = true; ks = lengthCoeffs[i][4]; r0 = lengthCoeffs[i][5]; } if (!found) { // something innocuous ks = 2.0; r0 = 1.2; } } public boolean redundant (lterm t) { if ((t.myAtoms[0] == myAtoms[0] && t.myAtoms[1] == myAtoms[1]) || (t.myAtoms[1] == myAtoms[0] && t.myAtoms[0] == myAtoms[1])) return true; return false; } public void enumerate (Vector atomList, Vector termList) { int i, j; for (i = 0; i < atomList.size (); i++) { atom a1 = (atom) atomList.elementAt (i); for (j = 0; j < a1.bonds.size (); j++) { atom a2 = (atom) a1.bonds.elementAt (j); lterm t = new lterm (a1, a2); if (!t.redundant (termList)) termList.addElement (t); } } } /* I've been hacking a bit with acceptable forms for the length-energy term. * _Nanosystems_ lists one simple expression with a square and cubic term * for close distances, and various kinds of exponential expressions for * farther away. But these depend on numbers called De and beta and I don't * have tables for those (just the few examples in _Nanosystems_). So after * some puttering and graphing and algebra, I came up with the following. * Energies are in attojoules (10^-18 joules), distances are in angstroms * (10^-10 m), forces are in attojoules per angstrom (10^-8 N), and * spring constants are in attojoules per angstrom-squared (100 N/m). This * is consistent with the units in _Nanosystems_, and also the units used * in the 'mm2.prm' parameter file I found on Jay Ponder's ftp site. * Expressions for energy and force, where ks and r0 are a function of * the two atoms involved: * * dr = r - r0; * * Energy (r) := * if (dr < rthresh) * { return 0.5 * ks * dr * dr * (1 - kc * dr) - z * ks; } * else * { return (-a * ks / beta) * exp (-beta * dr); } * * Force (r) := * if (dr < rthresh) * { return -(ks * dr * (1 - 1.5 * kc * dr)); } * else * { return -a * ks * exp (-beta * dr); } * * * A, beta, z, kcubic, and rthresh are defined below. Notice my use of * 'beta' differs from the use in _Nanosystems_, but plays a similar role. * * These follow Drexler's quadratic-cubic form for dr < rthresh, and for * dr > rthresh, they follow an exponential drop-off which is pretty * continuous for both energy and force. The only point that I think would * be worth quibbling over would be the value of beta, which I picked so * that it would look "reasonable" on a graph. */ private static final double kcubic = 2.0; /* (0.5 angstrom)^-1 */ private static final double rthresh = 1 / (3 * kcubic); private static final double beta = 3.0; private static final double a = rthresh * (1 - 1.5 * kcubic * rthresh) * Math.exp (beta * rthresh); private static final double z = (a / beta) * Math.exp(-beta * rthresh) + 0.5 * rthresh * rthresh * (1 - kcubic * rthresh); public void computeForces () { int i; // compute forces on each atom, add it to the atom's force vector double[] diff = new double[3]; double r, m; for (i = 0, r = 0.0; i < 3; i++) { diff[i] = myAtoms[0].x[i] - myAtoms[1].x[i]; r += diff[i] * diff[i]; } r = Math.sqrt (r); double rdiff = (r - r0), expr; if (rdiff < rthresh) m = ks * rdiff * (1 - 1.5 * kcubic * rdiff); else m = a * ks * Math.exp (-beta * rdiff); // at this point, m is du/dr m /= r; // m > 0 attract, m < 0 repel for (i = 0; i < 3; i++) { myAtoms[0].f[i] -= m * diff[i]; myAtoms[1].f[i] += m * diff[i]; } } // Coefficient data // Forces are in aJ per square angstrom (spring constants) private final static double[][] lengthCoeffs = { { C, atom.SP3, C, atom.SP3, 4.400, 1.523 }, { C, atom.SP3, C, atom.SP2, 4.400, 1.497 }, { C, atom.SP3, C, atom.SP, 5.200, 1.470 }, { C, atom.SP3, H, atom.NONE, 4.600, 1.113 }, { C, atom.SP3, O, atom.SP3, 5.360, 1.402 }, { C, atom.SP3, N, atom.SP3, 5.100, 1.438 }, { C, atom.SP3, N, atom.SP2, 3.520, 1.437 }, { C, atom.SP3, N, atom.SP3, 5.100, 1.499 }, { C, atom.SP3, O, atom.SP2, 5.360, 1.414 }, { C, atom.SP2, C, atom.SP2, 9.600, 1.337 }, { C, atom.SP2, C, atom.SP, 9.900, 1.313 }, { C, atom.SP2, H, atom.NONE, 4.600, 1.101 }, { C, atom.SP2, O, atom.SP3, 6.000, 1.355 }, { C, atom.SP2, N, atom.SP3, 6.320, 1.377 }, { C, atom.SP2, N, atom.SP2, 5.000, 1.410 }, { C, atom.SP2, O, atom.SP2, 10.000, 1.225 }, { C, atom.SP, C, atom.SP, 15.600, 1.212 }, { C, atom.SP, H, atom.NONE, 5.900, 1.090 }, { C, atom.SP, N, atom.SP, 17.730, 1.158 }, { O, atom.SP3, O, atom.SP3, 3.950, 1.428 }, { N, atom.SP3, N, atom.SP3, 5.600, 1.381 }, { N, atom.SP3, H, atom.NONE, 6.100, 1.045 } }; }