/*------------------------------------------------------------- Copyright (C) 2000 Peter Clote. All Rights Reserved. Permission to use, copy, modify, and distribute this software and its documentation for NON-COMMERCIAL purposes and without fee is hereby granted provided that this copyright notice appears in all copies. THE AUTHOR AND PUBLISHER MAKE NO REPRESENTATIONS OR WARRANTIES ABOUT THE SUITABILITY OF THE SOFTWARE, EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, OR NON-INFRINGEMENT. THE AUTHORS AND PUBLISHER SHALL NOT BE LIABLE FOR ANY DAMAGES SUFFERED BY LICENSEE AS A RESULT OF USING, MODIFYING OR DISTRIBUTING THIS SOFTWARE OR ITS DERIVATIVES. -------------------------------------------------------------*/ Here are 2 references for 3 dimensional lattices, pertinent to our project. Molecular Solid State Physics G.G. Hall Springer-Verlag 1991 An Introduction to Mathematical Crystallography M.A. Jaswon American Elsevier Publishing Co. 1965 Neither of these references has anything similar to what I've worked out below, but some drawings of crystals and the description of packing with interstices in the layer above were helpful. The mathematical theory of crystallography involves some group theory (symmetries form groups), and appears to be very elegant. Three dimensional lattices have a classification scheme to to Bravais, who conjectured that macroscopic symmetry of crystals is closely related to the microscopic symmetry. Cubic closed packed (ccp) and hexagonal closed packed (hcp) crystals are the closest packings of balls in 3 dimensions; for each the number of nearest neighbors of a ball is 12. CCP is formed by packing balls in a plane, each ball having 6 nearest neighbors (thus the neighbors of a center ball form a hexagon); in the second layer, there are 6 ``holes'' or interstices where balls may tightly fit -- only 3 may be used in packing the second layer; in the third layer, one may either pack balls using the interstices where the third layer balls are directly above the first layer balls (hexagonal closed packed) or not (cubic closed packed). Let's begin by relegating ouselves to the xy-plane (i.e. one layer). We assume that the vector (1,0) lies along the north compass direction (really due east, but we take this as north). Other directions are NE, SE, S, SW, NW, where S is not allowed, since we would backtrack on the vector from (0,0) to (1,0), thus not producing a self-avoiding walk. Direction Angle THETA ============================ N 0 NE PI/3 SE 2PI/3 S PI SW PI + PI/3 NW -PI/3 Recall that angles are positive when the rotation is counterclockwise, and recall that the coordinate transformation matrix in 2 dimensions is cos THETA -sin THETA sin THETA cos THETA Now to extend this to 3 dimensions, we now have the following directions: west up WU, north northeast NNEU, south southeast SSEU, east down ED, north northwest down NNWD, south southwest down SSWD. When transforming coordinate systems, by moving north N to one of the 11 directions (only S is disallowed): N,NE,SE,SW,NW,WU,NNEU,SSEU,ED,NNWD,SSWD we first perform a coordinate transformation by rotation about the z-axis, so the first angle THETA is in the xy-plane, then we perform a rotation about the y-axis; i.e. the positive x-axis is NOW identified with the new north direction, and we must incline this upwards or downwards. The angular inclination is PSI. We have the following table. Angles: THETA in xy-plane with axis of rotation z-axis, PSI in xz-plane with axis of rotation y-axis. We first perform the coordinate transformation for THETA, then PSI. CCP lattice Direction THETA PSI ==================================== N 0 0 NE PI/3 0 SE 2PI/3 0 SW PI + PI/3 0 NW -PI/3 0 WU PI/2 PI/3 NNEU -PI/6 PI/3 SSEU PI+PI/6 PI/3 ED -PI/2 -PI/3 NNWD PI/6 -PI/3 SSWD PI-PI/6 -PI/3 Let A be the transformation matrix cos THETA -sin THETA 0 sin THETA cos THETA 0 0 0 1 and B the transformation matrix cos PSI 0 -sin PSI 0 1 0 sin PSI 0 cos PSI Then given the current matrix M, perform the update M = B x A x M recalling that matrix multiplication is NOT commutative. When starting, first initialize M to be the identity matrix 1 0 0 0 1 0 0 0 1 If possible, I think it would be best to incorporate the matrix approach with angles THETA and PSI into the entire program, RATHER than having a switch statement as presently the case. One could introduce variables cosTHETA, sinTHETA, cosPSI, sinPSI whose values are cos(THETA), etc. In the cubic lattice case, these values are 0,1 or -1. Nevertheless, by introducing variables as indicated, one could lift the current program to allow the user to indicate by a flag the type of lattice; eg a.out -ccp 3 would mean cubic close packed in 3 dimensions, a.out -ccp 2 would mean cubic close packed in 2 dimensions (i.e. triangular, where each point has 6 nearest neighbors), a.out -c 3 a.out -c 2 would mean cubic lattice in 3 dimensions, respectively cubic lattice in 2 dimensions (i.e. rectangular 2 dimensional lattice). We could incorporate the hexagonal close packed (hcp), where the directions would be N,NE,SE,SW,NW,WU,NNEU,SSEU,WD,NNED,SSED and we would have to work out the 3 new angles. This should yield the following. HCP lattice Direction THETA PSI ==================================== N 0 0 NE PI/3 0 SE 2PI/3 0 SW PI + PI/3 0 NW -PI/3 0 WU PI/2 PI/3 NNEU -PI/6 PI/3 SSEU PI+PI/6 PI/3 WD PI/2 -PI/3 NNED -PI/6 -PI/3 SSED PI+PI/6 -PI/3 The transformation matrices A,B will be exactly as above, but with these different angles. However, we should first incorporate the ccp rather than the hcp lattice, since the former is used in a paper on approximate protein folding on the HP model (see Journal of Computational Biology 1997 in the library, where Martin Farach is one of the co-authors). Now, everything above concerned only the angles for the transformation matrices. In going from a sequence of relative directions in the ccp lattice, we need distances in x,y,z direction as well. This is easy from trigonometry. Notice that y = -sin(THETA), rather than +sin(THETA), since we are taking positive x axis to be north pole, and positive THETA rotations begin by sweeping into the -y direction. CCP lattice Dir THETA PSI x y z ================================================================== - - - cos(THETA) -sin(THETA) sin(PSI) N 0 0 1 0 0 NE PI/3 0 0.5 -0.866025404 0 SE 2PI/3 0 -0.5 -0.866025404 0 SW PI+PI/3 0 -0.5 0.8660254504 0 NW -PI/3 0 0.5 0.8660254504 0 WU PI/2 PI/3 0 -1 0.8660254504 NNEU -PI/6 PI/3 0.8660254504 0.5 0.8660254504 SSEU PI+PI/6 PI/3 -0.8660254504 0.5 0.8660254504 ED -PI/2 -PI/3 0 1 -0.8660254504 NNWD PI/6 -PI/3 0.8660254504 -0.5 -0.8660254504 SSWD PI-PI/6 -PI/3 -0.8660254504 -0.5 -0.8660254504 I did this just now, rather quickly, using Excel. In place of 0.8660254504 we should let the program compute sqrt(3)/2 to have as much precision as the C compiler allows. Also would some of you please check these values independently of me. We can go over this next class, but need to have an independent check in case of error.