CHAPTER 12

Unit Cell and Coordinate Systems

 

Although we mentioned slightly about the volume of a unit cell and polar orientation, we have not seriously involved crystallographic topics yet.  With the tools introduced in the previous chapters, this chapter provides an ¡®orientation¡¯ for the rest of the crystallographic concepts.  Several coordinate systems are used in crystallography, three of which are directly related to the concept of unit cell.  Fractional coordinate system is often used for symmetry operation (See Chapter 13).  Cartesian or orthogonal coordinate system services as an interface for crystallography to communicate with other scientific fields.  The coordinate system in the reciprocal space is heavily used without a name.  Miller indices are the coordinates in the discrete reciprocal space.  See Chapter 13 for details.

 

12.1 ccl::cell Class

 

The class ccl::cell implements the concept of unit cell.  It can be either a direct or a reciprocal cell.  This class itself contains the information about the type, and provides functions for the inter-conversion between direct and reciprocal cells.  This class also provides matrices for conversions between various coordinate systems.

 

12.1.1 Constructors, assignment and parameter access

 

cell(const double a     = 0,

     const double b     = 0,

     const double c     = 0,

     const double alpha = 0,

     const double beta  = 0,

     const double gamma = 0, const bool d = 1);

cell(const mac::vector3D<double> l,

     const mac::vector3D<double> a, const bool d = 1);

The general constructors take either two mac::vector3D<double>s and a Boolean, or six doubles and a Boolean as the cell lengths a, b, c and angles a, b, g.  The Boolean indicates whether the cell is direct.  The default constructor constructs a cell with 0-lengths and 0-angles, which will cause assertion error if further calculation is carried out.  This will prevent unintended result.

 

cell(const cell& c);

A copy constructor is available.

 

const cell& operator=(const cell& c);

The assignment operator = is defined.

 

void setA(const double a);

void setB(const double b);

void setC(const double c);

void setAlpha(const double a);

void setBeta(const double b);

void setGamma(const double g);

void setLengths(const mac::vector3D<double> l);

void setAngles(const  mac::vector3D<double> a);

void setCell(const double a,

             const double b,

             const double c,

             const double alpha,

             const double beta,

             const double gamma);

void setCell(const mac::vector3D<double> l,

             const mac::vector3D<double> a);

Each individual cell parameter, the cell lengths, angles, or all six parameters can be set by access functions.

 

double getA()     const;

double getB()     const;

double getC()     const;

double getAlpha() const;

double getBeta()  const;

double getGamma() const;

mac::vector3D<double> getLengths() const;

mac::vector3D<double> getAngles()  const;

Return individual cell parameter, the cell lengths, and angles.

 

isDirect()

isReciprocal()

Return Boolean that tests the cell type.

 

12.1.2 Cell volume and comparison

 

static double volume(const cell&);

double        volume() const;

The cell volume is calculated by a static and a member function volume().

 

bool operator==(const cell&) const;

bool operator!=(const cell&) const;

bool operator&=(const cell&) const;

The operators ==, &= and != test absolutely equal, approximately equal, and not equal between two cells, respectively.

 

bool operator< (const cell&) const;

bool operator<=(const cell&) const;

bool operator> (const cell&) const;

bool operator>=(const cell&) const;

The comparison operator < is defined by comparing the cell volumes.  If two cell volumes are equal, use the comparison of cell lengths.  If the cell lengths are all equal, use the comparison of cell angles.  The other comparison operators >, <=, and >= are defined accordingly.

 

12.1.3 Conversion between direct and reciprocal cells

 

Every cell has its reciprocal cell.  The reciprocal cell of a direct cell is, of course, called reciprocal cell.  The reciprocal cell of a reciprocal cell becomes the direct cell again.

 

static cell reciprocal(const cell&,

                       mac::matrix3D<double>* = NULL);

static cell reciprocalScalar(const cell&);

Two static functions reciprocal() and reciprocalScalar() are defined to convert a cell to its reciprocal cell.  The former is implemented in a vector form.  The latter is a scalar form and obsolete, but it serves as an example to show a scalar implementation.  The vector form reciprocal() may also modify a mac::matrix3D<double>, if a valid pointer points to it when the vector form is called.  The matrix will store three vectors by row, which can be accessed by using the operator [] of mac::matrix3D<double>.  These vectors are a*, b*, and c* in the orthogonal coordinate system of the cell.  See below for details of coordinate systems.

 

cell reciprocal(mac::matrix3D<double>* = NULL);

The member function reciprocal() converts the cell itself to its reciprocal cell, and optionally modifies a matrix if a valid pointer is given.

 

cell operator*() const;

cell operator~() const;

Two unary operators * and ~ have the exactly same definition that returns the reciprocal cell, and leaves the original cell unchanged.

 

12.1.4 Coordinate systems and orthogonalization

 

The three cell edges a, b, c naturally form a coordinate system, often called crystallographic coordinate system.  This system could be orthogonal for some crystal systems, but it is not in general.  Either the absolute length in Å or the three cell lengths a, b, c can be the units.  If the cell lengths are used as units for their axes, this crystallographic coordinate system is a fractional coordinate system.  Symmetry operation from one atomic position to its symmetry-related positions is very convenient in the fractional system.

 

However, arbitrary rotation is difficult using fractional coordinate system.  An orthogonal system is therefore required.  It is entirely a convention in the field to choose one unique orthogonal coordinate system respect to the crystallographic coordinate systems.  The unique orthogonal system also serves as an interface when structural results are communicated to other scientific fields.  We therefore use the definition of orthogonal system from the Protein Data Bank (PDB; http://www.rcsb.org/pdb).  The PDB convention is to select a-axis as the x-axis, a ´ b as the z-axis, and of course z ´ x as y-axis.  Considering the relationship between direct and reciprocal cells, this is also to say that c* is selected as z-axis, b* ¡Á c* as x-axis, and z ¡Á x as y-axis.  Sometimes other convention may also be used.  For example, Rossmann convention selects a as x-axis, c ¡Á a as y-axis, and x ¡Á y as z-axis, that is equivalent to b* as y-axis, b* ¡Á c* as x-axis, and x ¡Á y as z-axis.

 

static mac::matrix3D<double>   orthogonal(

   const cell&, const bool = 1);

static mac::matrix3D<double> deorthogonal(

   const cell&, const bool = 1);

Two static functions orthogonal() and deorthogonal() take a cell and an optional Boolean, and return the orthogonalization and deorthogonalization matrices, respectively.  The Boolean specifies whether PDB convention is used.  The default is certainly the PDB convention.  If the Boolean is unset, Rossmann convention will be used.  The orthogonalization matrix pre-multiplies a vector in crystallographic coordinate system, and results in a vector in the orthogonal system.  The deorthogonalization matrix is the inverse of orthogonalization matrix.  The required argument of these functions is ccl::cell, which can be either a direct or reciprocal cell.  If it is a direct cell, the returned matrices can ONLY pre-multiply a vector in the direct coordinate system.  It is an error to pre-multiply vectors in the reciprocal system.  USE CAUTION.

 

mac::matrix3D<double>   orthogonal(

   const bool = 1) const;

mac::matrix3D<double> deorthogonal(

   const bool = 1) const;

Two member functions calculate the orthogonalization and deorthogonalization matrices, respectively.  If the optional Boolean is set, PDB convention will used; otherwise, use Rossmann convention.

 

mac::matrix3D<double> operator+() const;

mac::matrix3D<double> operator-() const;

The unary operators + and ¨C also return these matrices, but under the PDB convention only.

 

12.1.5 Transformation between direct and reciprocal coordinate systems

 

The previous section is about the transformation between a crystallographic coordinate system, direct or reciprocal, and an orthogonal system.  This section describes the transformation between two crystallographic coordinate systems.  Transformation of one location in a system to the other does not make much sense, since these two systems are reciprocal to each other and do not share the same unit.  However, the relative orientation between these two systems are fixed.  Transformation of one orientation in one system to the other is useful in the process of indexing diffraction patterns, for example.

 

static mac::matrix3D<double> directToReciprocal(cell);

static mac::matrix3D<double> reciprocalToDirect(cell);

mac::matrix3D<double> directToReciprocal() const;

mac::matrix3D<double> reciprocalToDirect() const;

Static and member functions directToReciprocal() and reciprocalToDirect() are available.  They return matrices that can be pre-multiplied to vectors in direct and reciprocal systems.  The static functions require an argument of direct or reciprocal cell.  The type of the cell does not affect the result.

 

12.1.6 Output

 

The operator << is overloaded to print the cell parameters to an ostream.

 

The following program tests the class ccl::cell.

 

#include <assert.h>

#include <CCL/cell.h>

 

main()

{

   ccl::cell c(40, 30, 80, 90, 90, 120); // Construct a direct cell

   ccl::cell d = c; // Assignment

   assert (d == c); // Absolutely equal

   assert (d &= c); // Approximately equal

   c.setA(50);      // Set a

   assert (d != c); // Not equal

 

   assert (c.isDirect());     // c is a direct cell.

   cout << c << endl << endl;

   c.reciprocal();

   assert (c.isReciprocal()); // c becomes a reciprocal cell now.

   cout << c << endl << endl;

 

   assert (ccl::cell::reciprocal(d) &= ccl::cell::reciprocalScalar(d));

      // The vector and scalar versions.  The scalar version is obsolete.

   assert (*d == ~d); // The unary operators * and ~ are the same.

 

   c = ~d; // c is the reciprocal cell of d.

   mac::vector3D<double> v(1, 1);

 

   // Orthogonalization

   assert (+d * v != +c * v);

      // Orthogonalization matrix +d can only pre-multiply a vector in direct

      // space.  Orthogonalization matrix +c pre-multiplies a numerically same

      // vector v, but now v is considered in reciprocal space, since c is the

      // reciprocal cell.

 

   assert (d.reciprocalToDirect() * (d.directToReciprocal() * v) &= v);

   assert (c.reciprocalToDirect() * (d.directToReciprocal() * v) &= v);

   assert (d.reciprocalToDirect() * (c.directToReciprocal() * v) &= v);

   assert (c.reciprocalToDirect() * (c.directToReciprocal() * v) &= v);

   assert (d.reciprocalToDirect() * (c.directToReciprocal() * v) &= v);

   assert (c.reciprocalToDirect() * (d.directToReciprocal() * v) &= v);

      // Transform vector v from direct space to reciprocal and back to

      // direct.  It does not make a difference whether a direct or reciprocal

      // cell is used.

 

   return 0;

}

 

Listing 12.1.0.0.1 Testing ccl:cell.

 

12.2 The Shadow Class cpl.ccl.symm.cell

 

The shadow class cpl.ccl.symm.cell has almost all functions of the original class ccl::cell.

 

12.2.1 Constructors and parameter access

 

cell(a = 0, b = 0, c = 0,

     alpha = 0, beta = 0, gamma = 0, direct = 1)

Constructs a cell object.  direct indicates whether the cell is direct cell.  If no arguments are given, the default constructor returns a cell with 0-lengths and 0-angles.  Further calculation based on the default cell may cause assertion error.

 

cell(len, ang, direct = 1)

Constructs a cell object with two cpl.mac.vm.vector3DDBLs being cell lengths and angles.  If the Boolean direct is set, returns a direct cell; otherwise a reciprocal cell.

 

If cl is a cell object, the following functions are defined:

 

cell(cl)

Constructs a new cell object identical to cl.

 

cl.setA(a)

cl.setB(b)

cl.setC(c)

cl.setAlpha(alpha)

cl.setBeta(beta)

cl.setGamma(gamma)

cl.setLengths(len)

cl.setAngles(ang)

cl.setCell(a, b, c, alpha, beta, gamma)

cl.setCell(len, ang)

Set or change the cell parameters of cl.  len and ang are cpl.mac.vm.vector3DDBL, the others are all floats.

 

cl.getA()

cl.getB()

cl.getC()

cl.getAlpha()

cl.getBeta()

cl.getGamma()

cl.getLengths()

cl.getAngles()

Return the cell parameters of cl.  getLengths and getAngles return cpl.mac.vm.vector3DDBL, the others return floats.

 

cl.isDirect()

cl.isReciprocal()

Return Booleans that test whether cl is direct or reciprocal cell.

 

12.2.2 Cell volume and comparison

 

cell.volume(cl)

cl.volume()

The static and member functions compute the cell volume.

 

cl == c

cl.ap_eq(c)

cl != c

Return true if cl is absolutely, approximately, or not equal to another cell c, respectively.

 

cl < c

Returns true if the volume of cl is less than that of c.  If these cells are equal in volume, returns true if the vector of cl cell lengths is less than that of c.  If the cell length vectors are also equal, returns true if the vector of cl cell angles is less than that of c.  Returns false otherwise.

 

cl <= c

Returns true if cl < c or cl == c is true.

 

cl > c

Returns true if cl <= c is false.

 

cl >= c

Returns true if cl > c or cl == c is true.

 

12.2.3 Conversion between direct and reciprocal cells

 

cell_reciprocalScalar(cl)

cell.reciprocalScalar(cl)

Computes the reciprocal cell of cl by a scalar implementation.  This static function is obsolete.

 

cell_recip(cl, m = Null)

cell.recip(cl, m = Null)

Computes the reciprocal cell of cl by a vector implementation.  If m is a valid cpl.mac.vm.matrix3DDBL, it is modified to store vectors a*, b*, c* by rows in the orthogonal coordinate system of the cell cl.

 

cl.reciprocal(m = Null)

Converts cl to its reciprocal cell and returns the reciprocal cell.  See above for m.

 

~cl

Returns the reciprocal cell of cl.  cl remains unchanged.

 

12.2.4 Orthogonalization

 

cell.orthogonal(cl, pdb = 1)

cell.deorthogonal(cl, pdb = 1)

Return the orthogonalization and deorthogonalization matrices, respectively.  cl can be either a direct or a reciprocal cell.  If cl is a direct cell, the returned matrices can be used to pre-multiply vectors in the direct space only.  It will result in unintended solution if vectors in reciprocal space are pre-multiplied by these matrices.  Use a reciprocal cell cl to generate matrices for orthogonalization or deorthogonalization of vectors in reciprocal space.  pdb indicates whether the PDB convention is used.  If it is unset, Rossmann convention will be used.

 

cl.orthogonal(pdb = 1)

cl.deorthogonal(pdb = 1)

Return the orthogonalization and deorthogonalization matrices, respectively.  See above.

 

+cl

-cl

Return the orthogonalization and deorthogonalization matrices, respectively.  Always use the PDB convention.  See above.

 

12.2.5 Transformation between direct and reciprocal coordinate systems

 

cell_directToReciprocal(cl)

cell_reciprocalToDirect(cl)

Return matrices that can transform between direct and reciprocal coordinate systems.  cl can be either direct or reciprocal cell, and does not affect the results.

 

cl.directToReciprocal()

cl.reciprocalToDirect()

Same as above.

 

12.2.6 Output

 

cl.__str__()

The special function __str__ is added, so that print cl in Python will print out the cell parameters.

 

This example tests the shadow class cpl.ccl.symm.cell.  This program shall print a direct cell and its reciprocal cell, as well as the vectors a*, b*, c* in the PDB orthogonal coordinate system.

 

#!/usr/bin/python

 

import cpl

 

C = cpl.ccl.symm.cell

 

c = C(30, 50, 40, 90, 110, 90)  # Constructs a direct cell

d = cpl.ccl.symm.cell(c)        # Copy constructor

assert d == c                   # Absolutely equal

assert d.ap_eq(c)               # Approximately equal

c.setA(60)                      # Set a

assert d != c                   # Not equal

assert c.isDirect()             # c is a direct cell.

print c                         # Output

c.reciprocal()                  # Reciprocalize

print c

assert c.isReciprocal()         # c becomes a reciprocal cell now.

d = ~c # The unary operator ~.  Reciprocal of reciprocal.

m = cpl.mac.vm.matrix3DDBL()

assert cpl.ccl.symm.cell_recip(d, m).ap_eq(\

       cpl.ccl.symm.cell_reciprocalScalar(d))

assert cpl.ccl.symm.cell.recip(d, m).ap_eq(\

       cpl.ccl.symm.cell.reciprocalScalar(d))

   # The vector and scalar versions.  The scalar version is obsolete.

print 'a* =', m(1)

print 'b* =', m(2)

print 'c* =', m(3)

 

# Orthogonalization

v = cpl.mac.vm.vector3DDBL(1, 1)

assert +d * v != +c * v

   # Orthogonalization matrix +d can only pre-multiply a vector in direct

   # space.  Orthogonalization matrix +c pre-multiplies a numerically same

   # vector v, but now v is considered in reciprocal space, since c is the

   # reciprocal cell.

assert d.orthogonal() * v != c.orthogonal() * v

assert d.deorthogonal() * v != c.deorthogonal() * v

assert cpl.ccl.symm.cell.orthogonal(d) * v !=\

       cpl.ccl.symm.cell.orthogonal(c) * v

assert cpl.ccl.symm.cell.deorthogonal(d) * v !=\

       cpl.ccl.symm.cell.deorthogonal(c) * v

 

assert v.ap_eq(d.reciprocalToDirect() * (d.directToReciprocal() * v))

assert v.ap_eq(c.reciprocalToDirect() * (d.directToReciprocal() * v))

assert v.ap_eq(d.reciprocalToDirect() * (c.directToReciprocal() * v))

assert v.ap_eq(c.reciprocalToDirect() * (c.directToReciprocal() * v))

assert v.ap_eq(d.reciprocalToDirect() * (c.directToReciprocal() * v))

assert v.ap_eq(c.reciprocalToDirect() * (d.directToReciprocal() * v))

   # Transform vector v from direct space to reciprocal and back to

   # direct.  It does not make a difference whether a direct or reciprocal

   # cell is used.

 

Listing 12.2.0.0.1 Testing cpl.ccl.symm.cell.

 

12.3 Coordinates

 

A location in 3-dimensional space is described by three coordinates.  The numerical values of the coordinates depend on coordinate system and unit, so that one location may have different sets of coordinates.  Conversion between these sets is therefore required.  CCL implements three types of the most frequently-used coordinates and the inter-conversion between them: grid coordinates, fractional coordinates, and Cartesian coordinates.

 

12.3.1 Grid coordinates

 

If a unit cell length is divided into a number of grids, any location in the unit cell can be represented by the closest grid point.  The grid numbers along three cell edges are the grid coordinates.  Electron density map usually has the density values only on the grids, since it is often calculated by discrete Fourier transform.  Grid coordinates can only be integer.  A negative grid can be translated one or several unit cells to be converted into a non-negative grid.

 

12.3.1.1 ccl::gridCoordinates

The class ccl::gridCoordinates is derived from mac::vector3D<short>.

 

Constructors and destructor

gridCoordinates(const int x = 0,

                const int y = 0,

                const int z = 0,

                const unsigned short u = 1,

                const unsigned short v = 1,

                const unsigned short w = 1,

                const bool unity = 1);

The general constructor takes six integers.  The first three are the grid coordinates, the last three are the total grid numbers that all three cell edges are divided into.  The total grid numbers are also the cell dimensions in grid instead of in Angstrom.  By default, the constructor returns the origin.  However, the cell dimensions are set to 1, which is usually not intended.  It is not a good idea to use the default constructor to generate the origin.  If the Boolean unity is set, this constructor returns grid coordinates in the first unit cell.

 

gridCoordinates(

   const mac::vector3D<int>&,

   const mac::vector3D<unsigned short>&,

   const bool unity = 1);

Obviously, the first mac::vector3D<int> is the coordinates, and the mac::vector3D<unsigned short> is the cell dimensions.

 

gridCoordinates(

   const gridCoordinates&,

   const bool = 0);

gridCoordinates(

   const fractionalCoordinates&,

   const mac::vector3D<unsigned short>&,

   const bool = 1);

gridCoordinates(

   const CartesianCoordinates&,

   const mac::vector3D<unsigned short>&,

   const cell&,

   const bool = 1);

These three are the copy constructors.  The first is the ordinary one that copies one gridCoordinates to another.  The other two copy fractional coordinates and Cartesian coordinates to grid coordinates, respectively.  Both of these require the total grid dimensions, and the last copy constructor requires the unit cell.  See 12.3.2 and 12.3.3 for more detail.

 

All constructors including the copy constructor may accept an optional Boolean unity.  All defaults of this argument are true, except for the copy constructor.  By default, all constructors return gridCoordinates in the first unit cell, while the copy constructor does not translate the location of the coordinates unless explicitly indicated.

 

const gridCoordinates& operator=(

   const gridCoordinates&);

Assignment operator = is defined.

 

operator mac::vector3D<int>() const;

Converts to mac::vector3D<int>.

 

~gridCoordinates();

There is nothing to do in destructor.

 

Equivalent positions

gridCoordinates equivalent(const spaceGroup&,

                           const int i) const;

Returns the i-th equivalent position, where i is a 1-offset subscript.  The function requires the space group information.  i must be from 1 to the number of symmetry operators of the space group.  See the next chapter for space group.

 

vector<gridCoordinates> equivalent(

   const spaceGroup&,

   int* multiplicity = NULL) const;

Returns all equivalent positions in a vector<gridCoordinates>.  If the last pointer is valid, the multiplicity of this grid coordinates will also be returned.  This function obviously requires the symmetry information of the space group.  See the next chapter for space group.

 

vector<gridCoordinates> equivalentAU(

   const spaceGroup&,

   int* multiplicity = NULL);

This function has the exact same interface as the previous one, except this is a non-const version.  This function returns the same as the previous one, but after return, the original grid coordinates will be replaced by the equivalent position in the default asymmetric unit.  The default asymmetric unit is defined by the volume in which a gridCoordinates is less than any of its equivalent ones.  See 3.1.1.19 for vector comparison.

 

int multiplicity(const spaceGroup&) const;

Returns the multiplicity of this grid coordinates, given a space group.

 

Cell dimensions and grid distance

mac::vector3D<unsigned short>

   getDimensions() const;

Returns the cell dimensions in grid.

 

double gridDistance(

   const gridCoordinates& gc = gridCoordinates())

   const;

Returns grid distance from this gridCoordinates to the given one gc.  If no one is given, returns the grid distance to the origin.  Please note that grid distance is not the normal sense distance.  Its unit is grid.  Greater value may not mean longer distance.  Avoid using this function if possible or use caution.  This function may be removed or put into protected members in future release.

 

Output

friend ostream &operator<<(

   ostream&, const gridCoordinates&);

Operator << is defined to return an ostream.

 

#include <assert.h>

#include <ccl/coordinates.h>

 

main()

{

   typedef ccl::gridCoordinates       GC;

   typedef ccl::fractionalCoordinates FC;

   typedef ccl::CartesianCoordinates  CC;

   typedef mac::vector3D<unsigned short> V;

 

   GC gc; // Construct a default grid coordinates

   assert (!gc);

   assert (gc.getDimensions() == V(1, 1, 1));

   gc = GC(9, -28, 20, 30, 40, 50);

   assert (gc.gridDistance() == 25);

      // Grid distance is not a distance in the normal sense.

 

   const ccl::spaceGroup c2("c2");

   int m;

   const vector<GC> eq = gc.equivalent(c2, &m);

   assert (GC(V(9, 12, 20), V(30, 40, 50)) == gc); // Didn't change gc.

   assert (m == 1);                                // Multiplicity is 1.

   for (int i = 0; i < eq.size(); i++) cout << eq[i] << endl;

 

   gc.equivalentAU(c2, &m);

   assert (GC(V(9, 12, 20), V(30, 40, 50)) != gc); // gc is changed in AU.

 

   GC(0, 12, 25, 30, 40, 50).equivalent(c2, &m);

   assert (m == 2);