1 /* Copyright (C) 2002-2007 The Chemistry Development Kit (CDK) project 2 * 3 * Contact: cdk-devel@lists.sourceforge.net 4 * 5 * This program is free software; you can redistribute it and/or 6 * modify it under the terms of the GNU Lesser General Public License 7 * as published by the Free Software Foundation; either version 2.1 8 * of the License, or (at your option) any later version. 9 * All we ask is that proper credit is given for our work, which includes 10 * - but is not limited to - adding the above copyright notice to the beginning 11 * of your source code files, and to any copyright notice that you may distribute 12 * with programs based on this work. 13 * 14 * This program is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17 * GNU Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public License 20 * along with this program; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. 22 * */ 23 package org.openscience.cdk.geometry; 24 25 import java.util.Iterator; 26 27 import org.openscience.cdk.interfaces.IAtom; 28 import org.openscience.cdk.interfaces.IAtomContainer; 29 import org.openscience.cdk.interfaces.ICrystal; 30 31 import javax.vecmath.Point3d; 32 import javax.vecmath.Vector3d; 33 34 /** 35 * A set of static methods for working with crystal coordinates. 36 * 37 * @cdk.module standard 38 * @cdk.githash 39 * 40 * @author Egon Willighagen <egonw@sci.kun.nl> 41 * 42 * @cdk.keyword fractional coordinates, crystal 43 */ 44 public class CrystalGeometryTools { 45 46 /** 47 * Inverts three cell axes. 48 * 49 * @return a 3x3 matrix with the three Cartesian vectors representing 50 * the unit cell axes. The a axis is the first row. 51 */ calcInvertedAxes(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis)52 public static Vector3d[] calcInvertedAxes(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis) { 53 double det = aAxis.x * bAxis.y * cAxis.z - aAxis.x * bAxis.z * cAxis.y - aAxis.y * bAxis.x * cAxis.z + aAxis.y 54 * bAxis.z * cAxis.x + aAxis.z * bAxis.x * cAxis.y - aAxis.z * bAxis.y * cAxis.x; 55 Vector3d[] invaxes = new Vector3d[3]; 56 invaxes[0] = new Vector3d(); 57 invaxes[0].x = (bAxis.y * cAxis.z - bAxis.z * cAxis.y) / det; 58 invaxes[0].y = (bAxis.z * cAxis.x - bAxis.x * cAxis.z) / det; 59 invaxes[0].z = (bAxis.x * cAxis.y - bAxis.y * cAxis.x) / det; 60 61 invaxes[1] = new Vector3d(); 62 invaxes[1].x = (aAxis.z * cAxis.y - aAxis.y * cAxis.z) / det; 63 invaxes[1].y = (aAxis.x * cAxis.z - aAxis.z * cAxis.x) / det; 64 invaxes[1].z = (aAxis.y * cAxis.x - aAxis.x * cAxis.y) / det; 65 66 invaxes[2] = new Vector3d(); 67 invaxes[2].x = (aAxis.y * bAxis.z - aAxis.z * bAxis.y) / det; 68 invaxes[2].y = (aAxis.z * bAxis.x - aAxis.x * bAxis.z) / det; 69 invaxes[2].z = (aAxis.x * bAxis.y - aAxis.y * bAxis.x) / det; 70 return invaxes; 71 } 72 73 /** 74 * @cdk.dictref blue-obelisk:convertCartesianIntoFractionalCoordinates 75 */ cartesianToFractional(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis, Point3d cartPoint)76 public static Point3d cartesianToFractional(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis, Point3d cartPoint) { 77 Vector3d[] invaxis = calcInvertedAxes(aAxis, bAxis, cAxis); 78 Point3d frac = new Point3d(); 79 frac.x = invaxis[0].x * cartPoint.x + invaxis[0].y * cartPoint.y + invaxis[0].z * cartPoint.z; 80 frac.y = invaxis[1].x * cartPoint.x + invaxis[1].y * cartPoint.y + invaxis[1].z * cartPoint.z; 81 frac.z = invaxis[2].x * cartPoint.x + invaxis[2].y * cartPoint.y + invaxis[2].z * cartPoint.z; 82 return frac; 83 } 84 85 /** 86 * @cdk.dictref blue-obelisk:convertFractionIntoCartesianCoordinates 87 */ fractionalToCartesian(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis, Point3d frac)88 public static Point3d fractionalToCartesian(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis, Point3d frac) { 89 Point3d cart = new Point3d(); 90 cart.x = frac.x * aAxis.x + frac.y * bAxis.x + frac.z * cAxis.x; 91 cart.y = frac.x * aAxis.y + frac.y * bAxis.y + frac.z * cAxis.y; 92 cart.z = frac.x * aAxis.z + frac.y * bAxis.z + frac.z * cAxis.z; 93 return cart; 94 } 95 96 /** 97 * Calculates Cartesian vectors for unit cell axes from axes lengths and angles 98 * between axes. 99 * 100 * <p>To calculate Cartesian coordinates, it places the a axis on the x axes, 101 * the b axis in the xy plane, making an angle gamma with the a axis, and places 102 * the c axis to fulfill the remaining constraints. (See also 103 * <a href="http://server.ccl.net/cca/documents/molecular-modeling/node4.html">the 104 * CCL archive</a>.) 105 * 106 * @param alength length of the a axis 107 * @param blength length of the b axis 108 * @param clength length of the c axis 109 * @param alpha angle between b and c axes in degrees 110 * @param beta angle between a and c axes in degrees 111 * @param gamma angle between a and b axes in degrees 112 * @return an array of Vector3d objects with the three Cartesian vectors representing 113 * the unit cell axes. 114 * 115 * @cdk.keyword notional coordinates 116 * @cdk.dictref blue-obelisk:convertNotionalIntoCartesianCoordinates 117 */ notionalToCartesian(double alength, double blength, double clength, double alpha, double beta, double gamma)118 public static Vector3d[] notionalToCartesian(double alength, double blength, double clength, double alpha, 119 double beta, double gamma) { 120 Vector3d[] axes = new Vector3d[3]; 121 122 /* 1. align the a axis with x axis */ 123 axes[0] = new Vector3d(); 124 axes[0].x = alength; 125 axes[0].y = 0.0; 126 axes[0].z = 0.0; 127 128 double toRadians = Math.PI / 180.0; 129 130 /* some intermediate variables */ 131 double cosalpha = Math.cos(toRadians * alpha); 132 double cosbeta = Math.cos(toRadians * beta); 133 double cosgamma = Math.cos(toRadians * gamma); 134 double singamma = Math.sin(toRadians * gamma); 135 136 /* 2. place the b is in xy plane making a angle gamma with a */ 137 axes[1] = new Vector3d(); 138 axes[1].x = blength * cosgamma; 139 axes[1].y = blength * singamma; 140 axes[1].z = 0.0; 141 142 /* 3. now the c axis, with more complex maths */ 143 axes[2] = new Vector3d(); 144 double volume = alength 145 * blength 146 * clength 147 * Math.sqrt(1.0 - cosalpha * cosalpha - cosbeta * cosbeta - cosgamma * cosgamma + 2.0 * cosalpha 148 * cosbeta * cosgamma); 149 axes[2].x = clength * cosbeta; 150 axes[2].y = clength * (cosalpha - cosbeta * cosgamma) / singamma; 151 axes[2].z = volume / (alength * blength * singamma); 152 153 return axes; 154 } 155 156 /** 157 * @cdk.dictref blue-obelisk:convertCartesianIntoNotionalCoordinates 158 */ cartesianToNotional(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis)159 public static double[] cartesianToNotional(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis) { 160 double[] notionalCoords = new double[6]; 161 notionalCoords[0] = aAxis.length(); 162 notionalCoords[1] = bAxis.length(); 163 notionalCoords[2] = cAxis.length(); 164 notionalCoords[3] = bAxis.angle(cAxis) * 180.0 / Math.PI; 165 notionalCoords[4] = aAxis.angle(cAxis) * 180.0 / Math.PI; 166 notionalCoords[5] = aAxis.angle(bAxis) * 180.0 / Math.PI; 167 return notionalCoords; 168 } 169 170 /** 171 * Determines if this model contains fractional (crystal) coordinates. 172 * 173 * @return boolean indication that 3D coordinates are available 174 */ hasCrystalCoordinates(IAtomContainer container)175 public static boolean hasCrystalCoordinates(IAtomContainer container) { 176 Iterator<IAtom> atoms = container.atoms().iterator(); 177 while (atoms.hasNext()) { 178 if (atoms.next().getFractionalPoint3d() == null) { 179 return false; 180 } 181 } 182 return true; 183 } 184 185 /** 186 * Creates Cartesian coordinates for all Atoms in the Crystal. 187 */ fractionalToCartesian(ICrystal crystal)188 public static void fractionalToCartesian(ICrystal crystal) { 189 Iterator<IAtom> atoms = crystal.atoms().iterator(); 190 Vector3d aAxis = crystal.getA(); 191 Vector3d bAxis = crystal.getB(); 192 Vector3d cAxis = crystal.getC(); 193 while (atoms.hasNext()) { 194 IAtom atom = atoms.next(); 195 Point3d fracPoint = atom.getFractionalPoint3d(); 196 if (fracPoint != null) { 197 atom.setPoint3d(fractionalToCartesian(aAxis, bAxis, cAxis, fracPoint)); 198 } 199 } 200 } 201 } 202