org.openscience.cdk.geometry.CrystalGeometryTools Maven / Gradle / Ivy
/* Copyright (C) 2002-2007 The Chemistry Development Kit (CDK) project
*
* Contact: [email protected]
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public License
* as published by the Free Software Foundation; either version 2.1
* of the License, or (at your option) any later version.
* All we ask is that proper credit is given for our work, which includes
* - but is not limited to - adding the above copyright notice to the beginning
* of your source code files, and to any copyright notice that you may distribute
* with programs based on this work.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
* */
package org.openscience.cdk.geometry;
import java.util.Iterator;
import org.openscience.cdk.interfaces.IAtom;
import org.openscience.cdk.interfaces.IAtomContainer;
import org.openscience.cdk.interfaces.ICrystal;
import javax.vecmath.Point3d;
import javax.vecmath.Vector3d;
/**
* A set of static methods for working with crystal coordinates.
*
* @cdk.module standard
* @cdk.githash
*
* @author Egon Willighagen <[email protected]>
*
* @cdk.keyword fractional coordinates, crystal
*/
public class CrystalGeometryTools {
/**
* Inverts three cell axes.
*
* @return a 3x3 matrix with the three Cartesian vectors representing
* the unit cell axes. The a axis is the first row.
*/
public static Vector3d[] calcInvertedAxes(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis) {
double det = aAxis.x * bAxis.y * cAxis.z - aAxis.x * bAxis.z * cAxis.y - aAxis.y * bAxis.x * cAxis.z + aAxis.y
* bAxis.z * cAxis.x + aAxis.z * bAxis.x * cAxis.y - aAxis.z * bAxis.y * cAxis.x;
Vector3d[] invaxes = new Vector3d[3];
invaxes[0] = new Vector3d();
invaxes[0].x = (bAxis.y * cAxis.z - bAxis.z * cAxis.y) / det;
invaxes[0].y = (bAxis.z * cAxis.x - bAxis.x * cAxis.z) / det;
invaxes[0].z = (bAxis.x * cAxis.y - bAxis.y * cAxis.x) / det;
invaxes[1] = new Vector3d();
invaxes[1].x = (aAxis.z * cAxis.y - aAxis.y * cAxis.z) / det;
invaxes[1].y = (aAxis.x * cAxis.z - aAxis.z * cAxis.x) / det;
invaxes[1].z = (aAxis.y * cAxis.x - aAxis.x * cAxis.y) / det;
invaxes[2] = new Vector3d();
invaxes[2].x = (aAxis.y * bAxis.z - aAxis.z * bAxis.y) / det;
invaxes[2].y = (aAxis.z * bAxis.x - aAxis.x * bAxis.z) / det;
invaxes[2].z = (aAxis.x * bAxis.y - aAxis.y * bAxis.x) / det;
return invaxes;
}
/**
* @cdk.dictref blue-obelisk:convertCartesianIntoFractionalCoordinates
*/
public static Point3d cartesianToFractional(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis, Point3d cartPoint) {
Vector3d[] invaxis = calcInvertedAxes(aAxis, bAxis, cAxis);
Point3d frac = new Point3d();
frac.x = invaxis[0].x * cartPoint.x + invaxis[0].y * cartPoint.y + invaxis[0].z * cartPoint.z;
frac.y = invaxis[1].x * cartPoint.x + invaxis[1].y * cartPoint.y + invaxis[1].z * cartPoint.z;
frac.z = invaxis[2].x * cartPoint.x + invaxis[2].y * cartPoint.y + invaxis[2].z * cartPoint.z;
return frac;
}
/**
* @cdk.dictref blue-obelisk:convertFractionIntoCartesianCoordinates
*/
public static Point3d fractionalToCartesian(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis, Point3d frac) {
Point3d cart = new Point3d();
cart.x = frac.x * aAxis.x + frac.y * bAxis.x + frac.z * cAxis.x;
cart.y = frac.x * aAxis.y + frac.y * bAxis.y + frac.z * cAxis.y;
cart.z = frac.x * aAxis.z + frac.y * bAxis.z + frac.z * cAxis.z;
return cart;
}
/**
* Calculates Cartesian vectors for unit cell axes from axes lengths and angles
* between axes.
*
* To calculate Cartesian coordinates, it places the a axis on the x axes,
* the b axis in the xy plane, making an angle gamma with the a axis, and places
* the c axis to fulfill the remaining constraints. (See also
* the
* CCL archive.)
*
* @param alength length of the a axis
* @param blength length of the b axis
* @param clength length of the c axis
* @param alpha angle between b and c axes in degrees
* @param beta angle between a and c axes in degrees
* @param gamma angle between a and b axes in degrees
* @return an array of Vector3d objects with the three Cartesian vectors representing
* the unit cell axes.
*
* @cdk.keyword notional coordinates
* @cdk.dictref blue-obelisk:convertNotionalIntoCartesianCoordinates
*/
public static Vector3d[] notionalToCartesian(double alength, double blength, double clength, double alpha,
double beta, double gamma) {
Vector3d[] axes = new Vector3d[3];
/* 1. align the a axis with x axis */
axes[0] = new Vector3d();
axes[0].x = alength;
axes[0].y = 0.0;
axes[0].z = 0.0;
double toRadians = Math.PI / 180.0;
/* some intermediate variables */
double cosalpha = Math.cos(toRadians * alpha);
double cosbeta = Math.cos(toRadians * beta);
double cosgamma = Math.cos(toRadians * gamma);
double singamma = Math.sin(toRadians * gamma);
/* 2. place the b is in xy plane making a angle gamma with a */
axes[1] = new Vector3d();
axes[1].x = blength * cosgamma;
axes[1].y = blength * singamma;
axes[1].z = 0.0;
/* 3. now the c axis, with more complex maths */
axes[2] = new Vector3d();
double volume = alength
* blength
* clength
* Math.sqrt(1.0 - cosalpha * cosalpha - cosbeta * cosbeta - cosgamma * cosgamma + 2.0 * cosalpha
* cosbeta * cosgamma);
axes[2].x = clength * cosbeta;
axes[2].y = clength * (cosalpha - cosbeta * cosgamma) / singamma;
axes[2].z = volume / (alength * blength * singamma);
return axes;
}
/**
* @cdk.dictref blue-obelisk:convertCartesianIntoNotionalCoordinates
*/
public static double[] cartesianToNotional(Vector3d aAxis, Vector3d bAxis, Vector3d cAxis) {
double[] notionalCoords = new double[6];
notionalCoords[0] = aAxis.length();
notionalCoords[1] = bAxis.length();
notionalCoords[2] = cAxis.length();
notionalCoords[3] = bAxis.angle(cAxis) * 180.0 / Math.PI;
notionalCoords[4] = aAxis.angle(cAxis) * 180.0 / Math.PI;
notionalCoords[5] = aAxis.angle(bAxis) * 180.0 / Math.PI;
return notionalCoords;
}
/**
* Determines if this model contains fractional (crystal) coordinates.
*
* @return boolean indication that 3D coordinates are available
*/
public static boolean hasCrystalCoordinates(IAtomContainer container) {
Iterator atoms = container.atoms().iterator();
while (atoms.hasNext()) {
if (atoms.next().getFractionalPoint3d() == null) {
return false;
}
}
return true;
}
/**
* Creates Cartesian coordinates for all Atoms in the Crystal.
*/
public static void fractionalToCartesian(ICrystal crystal) {
Iterator atoms = crystal.atoms().iterator();
Vector3d aAxis = crystal.getA();
Vector3d bAxis = crystal.getB();
Vector3d cAxis = crystal.getC();
while (atoms.hasNext()) {
IAtom atom = atoms.next();
Point3d fracPoint = atom.getFractionalPoint3d();
if (fracPoint != null) {
atom.setPoint3d(fractionalToCartesian(aAxis, bAxis, cAxis, fracPoint));
}
}
}
}