/*   Written by Ralph C. Merkle, May 13 1991.
 *   Xerox PARC
 *   3333 Coyote Hill Road
 *   Palo Alto, CA 94304
 *   merkle@xerox.com
 *
 *   Copyright (c) Xerox Corporation 1991. All rights reserved.
 *
 *   License to copy and use this software is granted provided that
 *   appropriate credit is given to both Xerox and the author.
 *
 *   License is also granted to make and use derivative works provided
 *   that appropriate credit is given to both Xerox and the author.
 *
 *   Xerox Corporation makes no representations concerning either the
 *   merchantability of this software or the suitability of this software
 *   for any particular purpose.  It is provided "as is" without express or
 *   implied warranty of any kind.
 *
 *   These notices must be retained in any copies of any part of this software.
 *
 *   This program generates the description of a diamond block in
 *   PolyGraf or Brookhaven format.  The block can be wrapped into a
 *   tube, and the tube can then be wrapped into a toroid (ring).
 *   The tube is basically a sheet of diamond crystal wrapped into a
 *   tubular shape.  A few lines of additional code wraps the tube into
 *   a toroid (ring).
 *
 *   The program is invoked by saying:
 *    tube  ["brookhaven"] ["addH"] ["scale" <value>]  "block" | "tube" | "ring"
 *        OX OY OZ   SX SY SZ   TX TY TZ   AX AY AZ
 *        [ "grid" V1X V1Y V1Z   V2X V2Y V2Z
 *           [ "delete" P1X P1Y P1Z [box P2X P2Y P3Z]] |
 *           [ "add"    P1X P1Y P1Z [box P2X P2Y P3Z]] |
 *           [ "change" <from element name> "to" <to element name>
 *              [bondswanted <int>] P1X P1Y P1Z [box P2X P2Y P2Z]]
 *        ]*
 *   where all parameters are floating point numbers with the
 *   following meanings:
 *      "brookhaven":  The output format is the Brookhaven file format.
 *         (This format was deduced by examining a few Brookhaven
 *         formatted files.  It looks like it works, but....)
 *      "addH":        Undercoordinated atoms have hydrogens added
 *         in semi-reasonable places.
 *      "scale":  multiply the output coordinates by the scale factor,
 *         to increase and decrease the size.  Very useful when two
 *         structures are supposed to interlock.  The inner structure can
 *         be slightly shrunk before joint minimization is started.
 *      "block", "tube", or "ring":  specifies the general shape of
 *         the structure being generated.
 *      OX OY OZ -- (OX, OY, OZ) the origin of the structure.
 *      SX SY SZ -- (SX, SY, SZ) the plane of the tube surface.
 *         To generate a tube with a 111 surface, specify 1 1 1
 *         for these three parameters.
 *      TX TY TZ -- (TX, TY, TZ) the direction tangent to the tube
 *         surface.  S and T define a plane that bisects the tube
 *         longitudinally, dividing it into two tubular sections
 *         of half the length.
 *           The vectors S and T should be orthogonal.  If they
 *              are not orthogonal, an error message is generated
 *      AX AY AZ -- (AX, AY, AZ) the axis of the tube.  It must be
 *         orthogonal to the surface vector.  It defines the direction
 *         of the axis of the tube.  T and A need not be orthogonal
 *         If they are not orthogonal, the resulting toroid will have
 *         additional strain.
 *      "grid" is the literal keyword "grid", and indicates that two vectors
 *         defining the two-dimensional grid follow the keyword.
 *      V[12][XYZ] -- define two vectors that specify the grid of
 *         atoms for special processing.  This can be used to delete
 *         specific atoms in a plane, while not deleting other atoms in the same
 *         plane.
 *      "delete" is the literal keyword "delete", and indicates that a
 *         point is to be deleted.  The point [specified by three
 *         numbers that specify the coordinates of the point] follows the
 *         keyword.  Deletion of the point is repeated at all grid spacings,
 *         as indicated by the two grid vectors.
 *      "add" is the literal keyword "add", and behaves similarly but in the
 *         opposite sense as "delete"
 *      "change" indicates that <from element name> at the selected grid point
 *         is to be altered to <to element name>.  <from element name> can
 *         be "any", in which case any element at the grid point is changed.
 *      "bondswanted"  specify the integer number of bonds that the modified
 *         atom wants, e.g., carbon wants 4, nitrogen wants 3, etc.  If left
 *         unspecified, the default value of the previous atom will be used.
 *         This produces undesirable results when changing carbon to hydrogen,
 *         for example, where "bondswanted" should be specified as "1"
 *      PI[XYZ]  -- define the points to which the action is to be applied.
 *         All points P + K*V1 + J*V2 will be operated on for all integer
 *         values of K and J.  If "box" is specified, then the range of
 *         points inside the box defined by the first and second points
 *         will be operated on.
 *
 *    The output of this program is sometimes physically and chemically
 *    plausible, and sometimes not.  It is expected that the output will
 *    at least be entered into a molecular mechanics program for energy
 *    minimization. Whether or not the bonds in the resulting structure
 *    are too strained to be stable depends on the circumference and thickness.
 *    As the thickness increases and the circumference decreases, the result
 *    will be increasing strain.  The point at which this strain is too great
 *    for the structure to remain stable is not always obvious.  The program,
 *    though, will happily generate a structure that is quite unstable.
 *    Furthermore, surface atoms might require modification.  For example,
 *    a carbon atom might need to be replaced by a nitrogen atom at a site where
 *    there are only three other adjacent atoms (instead of the usual four).
 *    An alternative modification would be to add a hydrogen to the under-
 *    coordinated carbon.  Inspection of the output is required.
 *    Note that modifying the routine "AtomName" will allow the type of
 *    atom at each node to be modified during print out.  Further, by
 *    changing the "hydrogenation loop" that attempts to add hydrogens in
 *    reasonable places after the main atoms have been positioned, the
 *    types of atoms being produced can be further modified.
 *
 *    A simple example of a correct command line (we assume the executable
 *    is in the file "tube"):
tube block   0 0 0    2 0 0    0 8 0    0 0 4
       |       |        |        |        |
       |       |        |        |        Specify the axis
       |       |        |        Specify the tangent
       |       |        Specify the surface
       |       The origin.  Often 0 0 0
       The type of block being produced ("block", "tube", or "ring")

This program can be used to generate "buckytubes" in the following
(slightly indirect) fashion:  First, generate a "diamond-like" tube
which has a very thin 111 surface.  Then edit the output file,
replacing the sp3 coordinated atoms with sp2 coordinated
atoms.  Finally, do an energy minimization to produce the appropriate
tubular graphite sheet structure.

A simple example of an input that will produce an output which can
be easily edited into a "buckytube" is as follows:

tube tube  0.25 0.25 0.25  0.5 0.5 0.5  0 10 -10  -8 4 4


Another example (whose chemical stability is not known, but
which looks pretty):
tube ring 0 0 0   1 0 0   0 8 -6   0 24 22 \
    grid 0 0.5 -0.5 0 1 1  delete 0.75 0.75 0.25

The following two examples are the inner and outer races of
a simple molecular bearing.  Note that the inner race has been
reduced in size ("scale" is set to 0.9).  This allows both the
inner and outer race to be jointly minimized, starting from a
configuration in which the inner race fits within the outer race.

Specification of the inner race:
tube addH scale 0.9 tube  0 0.75 0.75  1.75 0 0  0 17 -17   0 5.25 5.25 \
  grid  0 0.5 -0.5  0 1 1 \
    delete  1.25 0.0 0.0 \
  grid  0 0.25 0  0 0 0.25 \
    change O_3 to S_3  1.5 0 0 \
    change N_3 to C_3 bondswanted 4   0 0 0 box 0.75 0 0

Specification of the outer race:
tube addH tube 1.25 1.5 1.5   1.75 0 0   0 23 -23  0 2.75 2.75 \
  grid 0 0.25 -0.25  0 1 1 \
    delete 0 0.5 0.5 \
    delete 0.25 0.5 0.5 \
  grid 0 0.25 -0.25  0 0.25 0.25 \
    change N_3 to C_3 bondswanted 4  0.75 0 0  box  4 0 0 \
  grid 0.25 0 0  0 0.25 -0.25 \
    delete 0 0 0 box 0 0.25 0.25

An example of a small toroid (of dubious physical stability):
tube  ring 0 0 0   1.50 0 0   0 5 3   0 53 -45 \
  grid  0 0.5 0.5   0 1.0 -1.0 \
      delete 1.25 0.75 0.75 \
  grid 0 0.25 0.25  0 0.25 -0.25 \
      change O_3 to S_3 1.25 0.25 0.25 \
  grid 0 5 3  0 7.25 -5.25 \
      delete 1.25 0 0  box 1.25 6.90625  1.78125


A larger toroid with several fiddly little changes:
-------------------------------------------------
tube brookhaven addH \
tube 0 0 0   1.75 0 0   0 8 -6   0  96 94 \
grid 0 0.5 -0.5   0 1 1  delete 1.5 0.5 0.00 \
grid 0 8 -6  0 3.25 2.75 \
   delete 1 0 0 box 2 9.5 -5.0 \
grid 0 8 -6   0 3.25 2.75 \
   delete 1.5 8.5 1.0 \
   delete 1.5 8.0 1.5 \
   delete 1.5 7.5 2.0 \
   delete 1.5 7.0 2.5 \
 \
   delete 1.5 13 0.5 \
   delete 1.5 13.5 0.0 \
   delete 1.5 6.0 5.5 \
   delete 1.5 6.5 5.0 \
  \
   delete 1.5 3.0 2.5 \
   delete 1.5 10.5 -3.0 \
   delete 1.5 10 -2.5 \
   delete 1.5 9.5 -2 \
  \
   delete 1.5 10.0 5.5 \
   delete 1.5 10.5 5.0 \
   delete 1.5 11.0 4.5 \
   delete 1.5 11.5 4.0 \
grid 0 0.5 -0.5  0 0.5 0.5 \
   change N_3 to C_3 bondswanted 4 0.75 0.75 0.25 \
   change N_3 to C_3 bondswanted 4 1.00 0.00 0.00 \
   change N_3 to C_3 bondswanted 4 1.25 0.25 0.25 \
grid 0 4 -3  0 3.25 2.75 \
   change C_3 to N_3 bondswanted 3  0.75  10.75  -0.75 \
   change C_3 to N_3 bondswanted 3  0.75   7.25  -3.25 \
   change C_3 to N_3 bondswanted 3  0.75   9.75   8.25 \
   change C_3 to N_3 bondswanted 3  0.75   6.75   5.25
-------------------------------------------------
 */
#define MAXATOMS     20000
#define MAXBONDS     8
#define VALENCE      4
#define DIMENSION    3
#define X 0
#define Y 1
#define Z 2
#define SURFACE 0
#define TANGENT 1
#define AXIS    2
#define XHASHSIZE 113
#define YHASHSIZE 23
#define ZHASHSIZE 11
#define WARNINGBUCKETSIZE 100
#define TRUE  1
#define EXTERIOR1 1
#define EXTERIOR2 2
#define INTERIOR1 -1
#define FALSE 0
#define NILATOM (-1)
/*  lattice consant for carbon.  */
#define LATTICECONSTANT 3.56683
#define TWOPI (2*3.14159265358979323846)
#define EPSILON 0.001
#define ODD 1
#define EVEN 0
#define BLOCKSHAPE  0
#define TUBESHAPE   1
#define RINGSHAPE   2
#define MAXSHAPE    3
#define MAXGRIDPOINTS 50
#define DELETEACTION   0
#define CHANGEACTION   1
#define ADDACTION      2
#define MAXACTION      3
#define MAXNAMELENGTH  23
#define POLYGRAFNAMELENGTH 5
#define BROOKHAVENNAMELENGTH 4
#define POLYGRAF       0
#define BROOKHAVEN     1

#include <stdio.h>
#include <math.h>
#include <string.h>


struct gridRecord {
    int    action;
    char   toName[MAXNAMELENGTH];
    int    toBondsWanted;
    char   fromName[MAXNAMELENGTH];
    double edge[DIMENSION][DIMENSION];
    double point1[DIMENSION];
    double point2[DIMENSION];
  };
/*  The following global definition is of a struct that holds the
 *  parameters that define the shape of the object being constructed.
 */
struct shapeDef {
int     shape;
double  scale;
double  origin[DIMENSION];
double  basis [DIMENSION][DIMENSION];
int     gridCount;
struct  gridRecord grid[MAXGRIDPOINTS];
  };

/*  The following global variables describe a collection of bonded atoms.
 *  The "coord" array holds the coordinates of each atom.
 *  The "atomType" array specifies whether the atom is "odd" or "even".
 *  The "bond" array holds the indices of the other atoms to which
 *    this atom is bonded.  bondCount[3][2] == 7 means atom 3 is
 *    bonded to atom 7 by the second bond coming from atom 3.  The
 *    bond numbers (0, 1, 2, 3, etc.) are semantically irrelevant.  They
 *    mean nothing.  In particular, just because bondCount[3][2]==7
 *    does not imply bondCount[7][2]==3.  It might instead be
 *    bondCount[7][1]==3 or bondCount[7][0]==3.
 *  The bondCount array holds the count of the number of bonds from
 *    this atom to other atoms.  bondCount[3] == 4 means atom 3 has
 *    4 bonds coming from it.
 *  Note that bonds are symmetric.  If i is bonded to j, then j is
 *    bonded to i.  This property is insured by the routine "BondAtoms"
 *  "lastAtomPlusOne" gives the index of the last atom in the array plus one.
 *  "firstAtom" indicates the first atom that has actual data.
 */
double	coord[MAXATOMS][DIMENSION];
int     atomType[MAXATOMS];
int	bond[MAXATOMS][MAXBONDS];
int	bondCount[MAXATOMS];
int     bondsWanted[MAXATOMS];
char    name[MAXATOMS][MAXNAMELENGTH];
int	lastAtomPlusOne = 1;
int     firstAtom = 1;

/*  The hash table global variables */
int     bucket[XHASHSIZE][YHASHSIZE][ZHASHSIZE];
int     nextLink[MAXATOMS];

/*  The following array gives the direction of the four neighboring
 *  atoms from a given atom in a diamond lattice.  By subtracting
 *  instead of adding, the odd/even atom distinction can be
 *  correctly maintained.  These four directions are mapped onto
 *  four new directions to take into account the rotation of the
 *  crystal structure specified by the user.  This modified set
 *  of vectors is then put into the array "delta".
 */
double delta[VALENCE][DIMENSION] =
      {
         { 0.25,  0.25,  0.25},
         {-0.25, -0.25,  0.25},
         { 0.25, -0.25, -0.25},
         {-0.25,  0.25, -0.25}
      };

char *shapeList[MAXSHAPE] = {
      "block",
      "tube",
      "ring"
      };
         

/*  A simple error routine.  Print an error message and die. */
ErrAbort(s)
char *s;
{
      fprintf(stderr,"%s\n", s);
      exit(0);
}

/*  Add two vectors: a = b + c */
VectorAdd(a,b,c)
double a[DIMENSION];
double b[DIMENSION];
double c[DIMENSION];
{int i;
   for(i=0; i<DIMENSION; i++)
     a[i] = b[i] + c[i];
}

/*  Substract two vectors: a = b - c */
VectorSub(a,b,c)
double a[DIMENSION];
double b[DIMENSION];
double c[DIMENSION];
{int i;
   for(i=0; i<DIMENSION; i++)
     a[i] = b[i] - c[i];
}

/*  copy a vector: a = b */
VectorCopy(a,b)
double a[DIMENSION];
double b[DIMENSION];
{int i;
    for(i=0; i<DIMENSION; i++) a[i]=b[i];
}

/*  print a vector */
VectorPrint(f,s,a)
FILE *f;
char *s;
double a[DIMENSION];
{
 /*  fprintf(f, " %s %15f %15f %15f\n", s, a[X], a[Y], a[Z]); */
}

/*  Computes the dot product of two three-dimensional vectors. */
double Dot(a,b)
double a[DIMENSION], b[DIMENSION];
{int i;
 double sum;
   sum = 0.0;
   for (i=0; i<DIMENSION; i++) sum += a[i]*b[i];
   return(sum);
}

/*  Multiply a vector times a scalar:  v = v * s  */
ScalarMul(v, s)
double v[DIMENSION];
double s;
{int i;
  for(i=0;i<DIMENSION;i++) v[i] *= s;
}

/* clear the "bottom" bit.  Note that little-endian/
 * big-endian machines behave differently.  This
 * routine should be portable, but hasn't been
 * tested.
 */
int ClearBottomBit(i)
int i;
{int bit;
  if(i >= 0) {
   bit = i % 2;
   return(i-bit);
  }
  else
  {
   bit = (-i) % 2;
   return(i-bit);
  };
}

/* Mod computes a correct mod function.  Note that %
 * does not compute a correct mod function.  This
 * mod function should be portable.
 */
int Mod(i,m)
int i,m;
{
   if (m <= 0) ErrAbort("Mod called with negative modulus");
   if (i >= 0)
     return(i%m);
   else
     return( (m-((-i)%m))%m );
}

/*  set the bottom bit.  */
int SetBottomBit(i, bit)
int i,bit;
{
   return(ClearBottomBit(i)+bit);
}

/*  computes the length of a vector. */
double Length(a)
double a[DIMENSION];
{
   return(sqrt(Dot(a,a)));
}

/*  computes the distance between two points */
double Distance(a,b)
double a[DIMENSION];
double b[DIMENSION];
{double temp[DIMENSION];
    temp[X] = a[X]-b[X];
    temp[Y] = a[Y]-b[Y];
    temp[Z] = a[Z]-b[Z];
    return(Length(temp));
}

/*  Initialize the hash table.
 *  The hash table preserves spatial relationships. If two
 *  points A and B are adjacent to each other in three-space,
 *  they will be in adjacent buckets in the integerized and
 *  folded version of three-space held in the array "bucket".
 *  This is very convenient when looking for "nearby" coordinates.
 */
InitializeHashTable()
{int i,j,k;
   for(i=0; i<XHASHSIZE; i++)
   for(j=0; j<YHASHSIZE; j++)
   for(k=0; k<ZHASHSIZE; k++)
     bucket[i][j][k] = NILATOM;
   for(i=0; i< MAXATOMS; i++) nextLink[i] = NILATOM;
}

/*  Determine if two coordinates for atoms are too close to each other */
TooClose(coord1, coord2)
double coord1[DIMENSION];
double coord2[DIMENSION];
{int i;
 int veryClose;
   veryClose = TRUE;
   for(i=0; i<DIMENSION; i++) 
     if(fabs(coord1[i]-coord2[i]) > EPSILON) veryClose = FALSE;
   return(veryClose);
}

/*  Compute an integer index value suitable for use in the hash table.
 *  Basically, simply integerize one of the three spatial
 *  coordinates and fold it into the proper range (given by
 *  the parameter "size").
 */
HashIndex(r,size)
double r;
int size;
{int temp;
   temp = (r*2.3458565130673);
   temp %= size;
   if (temp < 0) temp += size;
   if (temp < 0) ErrAbort("temp < 0");
   if (temp >= size) ErrAbort("temp >= size");
   return(temp);
}

/*  look in the hash bucket associated with the integer coordinates */
LookInHashBucket(i,j,k,atomCoord)
int i,j,k;
double atomCoord[DIMENSION];
{int linkPointer;
 int bucketSize;
 static int notWarned = TRUE;
    bucketSize = 0;
    for(linkPointer  = bucket[i][j][k];
        linkPointer != NILATOM;
        linkPointer  = nextLink[linkPointer]) {
           /*  We've found an atom in this bucket which is "too close"
            *  to the target coordinates we're searching for.
            *  Return the pointer to the atom and exit.
            */
           if (TooClose(coord[linkPointer], atomCoord)) return(linkPointer);
           /*  Diagnostics code.  If a bucket gets too big, performance
            *  goes through the floor.  Warn the user.  If this is a problem,
            *  increase the size of the array "bucket" and re-compile.
            */
           if (notWarned)
             if (bucketSize++ > WARNINGBUCKETSIZE) {
               fprintf(stderr,"Over %d entries in one bucket\n",WARNINGBUCKETSIZE);
               notWarned = FALSE;
             };
    };
    return(NILATOM);
}

/*  add an atom to the hash table */
AddToHashTable(atomNumber)
int atomNumber;
{int i,j,k;
    i=HashIndex(coord[atomNumber][X],XHASHSIZE);
    j=HashIndex(coord[atomNumber][Y],YHASHSIZE);
    k=HashIndex(coord[atomNumber][Z],ZHASHSIZE);
    if (nextLink[atomNumber] != NILATOM) {
       fprintf(stderr,"Attempt to add atom %d to hash table twice\n", atomNumber);
       ErrAbort("");
       };
    nextLink[atomNumber] = bucket[i][j][k];
    bucket[i][j][k] = atomNumber;
}

/*  Given the index of an atom in the data structure, add that atom
 *  to the collection of atoms in the global data structure.
 */
AddAtom(atomCoord,localAtomType)
double atomCoord[DIMENSION];
int    localAtomType;
{int i;
     if(lastAtomPlusOne >= MAXATOMS) ErrAbort("Too many atoms");
     for(i=0;i<DIMENSION;i++) coord[lastAtomPlusOne][i] = atomCoord[i];
     atomType[lastAtomPlusOne] = localAtomType;
     bondCount[lastAtomPlusOne]=0;
     bondsWanted[lastAtomPlusOne] = 0;
     strncpy(name[lastAtomPlusOne],"bad",MAXNAMELENGTH);
     AddToHashTable(lastAtomPlusOne);
     return(lastAtomPlusOne++);
}

/* Given the index of an atom previously added by AddAtom,
 * delete the atom from the global data structure.
 */
Delete(atom)
int  atom;
{
/*  Not implemented.  Any volunteers? */
}

/*  Returns the number of bonds that this atom "wants" to have.
 *  This is used when hydrogen atoms are automatically added --
 *  the number of additional hydrogen atoms bonded to the atom
 *  is simply the number of bonds that atom wants to have minus
 *  the number of bonds the atom already has.
 */
AtomWantsThisManyBonds(atomNumber)
int atomNumber;
{
   return(bondsWanted[atomNumber]);
}

/*  The following routine returns a 5 character string that
 *  specifies the "label" of the atom, given the atom number in
 *  the global data structure.
 */
char *AtomLabel(atomNumber,s,length)
int  atomNumber;
struct shapeDef s;
int  length;
{
   static char tempS[MAXNAMELENGTH];
   int i;
   strcpy(tempS,name[atomNumber]);
   for (i=strlen(name[atomNumber]);i<MAXNAMELENGTH;i++) tempS[i] = ' ';
   tempS[length]=0;
   return(tempS);
}

/*  The following routine returns a string that
 *  specifies the name of the atom, given the atom number in
 *  the global data structure.
 */
char *AtomName(atomNumber, s, length)
int  atomNumber;
struct shapeDef s;
int  length;
{
   static char tempS[MAXNAMELENGTH];
   int i;
   strcpy(tempS,name[atomNumber]);
   for (i=strlen(name[atomNumber]);i<MAXNAMELENGTH;i++) tempS[i] = ' ';
   tempS[length]=0;
   return(tempS);
}

BondAtoms(atom1, atom2)
/*  Create a bond between atoms 1 and 2.  If they are
 *  already bonded, do nothing.
 *  This routine works by adding the indices of the atoms
 *  to the "bond" array and bumping the bondCount value
 *  by 1 (if the atoms aren't already bonded).
 */
int atom1, atom2;
{int hit1, hit2;
 int i;
 int atomTypeOK;
     if(atom1 >= lastAtomPlusOne) ErrAbort("Atom1 number out of range");
     if(atom2 >= lastAtomPlusOne) ErrAbort("Atom2 number out of range");
     if(atom1 < firstAtom) ErrAbort("Atom1 number less than firstAtom");
     if(atom2 < firstAtom) ErrAbort("Atom2 number less than firstAtom");
     /*  If we want to bond two atoms, then one has to be "odd" and
      *  the other "even."  If this isn't true, something's wrong.
      */
     atomTypeOK = FALSE;
     if((atomType[atom1] == EVEN) && (atomType[atom2]==ODD)) atomTypeOK = TRUE;
     if((atomType[atom2] == EVEN) && (atomType[atom1]==ODD)) atomTypeOK = TRUE;
     if(atomTypeOK != TRUE) {
          printf(" types 1 and 2: %d %d\n",atomType[atom1], atomType[atom2]);
          printf(" atoms 1 and 2: %d %d\n", atom1, atom2);
          ErrAbort("Atom types being bonded are wrong");
     };
     /*  See if we've already got pointers from atom1 to atom2
      *  and vice versa.
      */
     hit1 = FALSE;
     hit2 = FALSE;
     for(i=0;i<bondCount[atom1];i++)
         if (bond[atom1][i] == atom2) hit1 = TRUE;
     for(i=0;i<bondCount[atom2];i++)
         if (bond[atom2][i] == atom1) hit2 = TRUE;
     if (hit1 != hit2) ErrAbort("Logic error, asymmetric bond");
     if(hit1 == FALSE) {
        /*  There's no bond from atom1 to atom2.  Add one.  */
        bond[atom1][bondCount[atom1]++] = atom2;
        if (bondCount[atom1] > MAXBONDS)
             ErrAbort("Bond overflow, logic error");
        if (bondCount[atom1] > VALENCE)
           fprintf(stderr,"Atom %d has more than %d bonds\n", atom1, VALENCE);
     };
     if(hit2 == FALSE) {
        /*  There's no bond from atom2 to atom2.  Add one.  */
        bond[atom2][bondCount[atom2]++] = atom1;
        if (bondCount[atom2] > MAXBONDS)
             ErrAbort("Bond overflow, logic error");
        if (bondCount[atom2] > VALENCE)
           fprintf(stderr,"Atom %d has more than %d bonds\n", atom2, VALENCE);
     };
}

/*  Get a parameter from the command line */
double GetParam(argc, argv, argNum, s)
int argc;
char *argv[];
int argNum;
char *s;
{int status;
 double temp;
 int i;
  if(argc <= argNum) {
      fprintf(stderr," argc   = %d\n argNum = %d\n",argc,argNum);
      ErrAbort("Invalid parameter number in GetParam");
  };
  /*  See if we're looking for "block", "tube", or "shape" */
  if (strcmp(s,"shape") == 0) {
     for (i=0;i<MAXSHAPE;i++)
     if (strcmp(argv[argNum], shapeList[i]) == 0)  return(i);
     fprintf(stderr,"Shape must be: block, tube, or ring\n");
  };
  status = sscanf(argv[argNum],"%lf", &temp);
  if (status != 1) {
    fprintf(stderr," Error reading parameter %d, status = %d.\n", argNum, status);
    ErrAbort("Ending run.");
  };
 fprintf(stderr,"%s:", s);
 for(i=0;i<14-strlen(s);i++) fprintf(stderr," ");
 fprintf(stderr,"%11f\n", temp);
 return(temp);
}

ExpectStringParam(argc, argv, argNum, s)
int argc;
char *argv[];
int argNum;
char *s;
{
  if(argc <= argNum) {
      fprintf(stderr," argc   = %d\n argNum = %d\n",argc,argNum);
      ErrAbort("Invalid parameter number in ExpectStringParam");
  };
  if(strcmp(argv[argNum],s)!=0){
    fprintf(stderr,"%s expected as parameter %d\n",s,argNum);
    ErrAbort("");
  };
}

/*
 *  If the values a, b, and c are linear, (e.g., are in
 *  the order a, b, c or c, b, a) then return TRUE.
 *  If they are not linear, return FALSE.
 */
InOrder(a, b, c)
double a,b,c;
{
  if ((a <= b) && (b <= c)) return(TRUE);
  if ((c <= b) && (b <= a)) return(TRUE);
  return(FALSE);
}


GetStringParam(argc, argv, argNum, target, length, s)
int argc;
char *argv[];
int argNum;
char *target;
int  length;
char *s;
{
  int i;
  if(argc <= argNum) {
      fprintf(stderr," argc   = %d\n argNum = %d\n",argc,argNum);
      ErrAbort("Invalid parameter number in GetStringParam");
  };
  strncpy(target, argv[argNum], length);
  target[length-1] = 0;
  fprintf(stderr,"%s:", s);
  for(i=0;i<14-strlen(s);i++) fprintf(stderr," ");
  fprintf(stderr,"%s\n", target);
}

/* Compute the greatest common divisor (GCD) of a and b */
long Gcd(a,b)
long a,b;
{
   long c;
   if (a<0) a = -a;
   if (b<0) b = -b;
   if (a==0) return(b);
   if (b==0) return(a);
   c = a%b;
   while (c != 0) {
     a=b;
     b=c;
     c = a%b;
   };
   return(b);
}

/*  This program generates a PolyGraf file describing a diamond tube.
 *  The tube can have grooves in it.
 */
main(argc, argv)
int     argc;
char    *argv[];
{
static  double seed[DIMENSION] = {0.0, 0.0, 0.0};
double  tempCoord[DIMENSION];
int	processNext; /*  the atom to process next */
int     tempAtom;
int     i,j,k;
int     tempType;
int     limitAtom;
int     currentAtom;
int     currentBond;
int     seedNum;
int     seedType;
int     argi;
int     lastV1V2;
double  tempVector[DIMENSION];
static  double zeroVector[DIMENSION] = {0.0, 0.0, 0.0};
double  cross[DIMENSION][DIMENSION];
long    longCross[DIMENSION][DIMENSION];
FILE    *streamOut;
int     addH;
int     outputFormat;

struct shapeDef s;



if (argc < 13){
     fprintf(stderr,"The command line should have the following format:\n");
     fprintf(stderr,"tube  [\"brookhaven\"] [\"addH\"]   \"block\" | \"tube\" | \"ring\"\n");
     fprintf(stderr,"   OX OY OZ   SX SY SZ   TX TY TZ   AX AY AZ\n");
     fprintf(stderr,"   [ \"grid\" V1X V1Y V1Z   V2X V2Y V2Z\n");
     fprintf(stderr,"      [ \"delete\" P1X P1Y P1Z [box P2X P2Y P3Z]] |\n");
     fprintf(stderr,"      [ \"add\"    P1X P1Y P1Z [box P2X P2Y P3Z]] |\n");
     fprintf(stderr,"      [ \"change\" <from element name> \"to\" <to element name>\n");
     fprintf(stderr,"         [bondswanted <int>] P1X P1Y P1Z [box P2X P2Y P2Z]]\n");
     fprintf(stderr,"   ]*\n");
     fprintf(stderr,"A vector O describing the origin :  0  0  0\n");
     fprintf(stderr,"A vector S describing the surface:  1  1  1\n");
     fprintf(stderr,"A vector T describing the tangent:  0  7 -7\n");
     fprintf(stderr,"A vector A describing the axis : -4  2  2\n");
     fprintf(stderr,"The surface and tangent should be orthogonal.\n");
     fprintf(stderr,"Example:\n");
     fprintf(stderr,"makeTube  tube 0 0 0   1 1 1   0 17 -17   -4 2 2\n");
     fprintf(stderr,"See the program for further examples\n");
     };
  argi = 1;

  outputFormat = POLYGRAF;
  if ( 
         strcmp(argv[argi],"brook") == 0
      || strcmp(argv[argi],"brookhaven") == 0
      || strcmp(argv[argi],"bkv") == 0
     ) {
      outputFormat = BROOKHAVEN;
      argi++;
     };
  addH = FALSE;
  if (
         strcmp(argv[argi],"addH") == 0
      || strcmp(argv[argi],"addHydrogen") == 0
     ) {
    addH = TRUE;
    argi++;
  };

  s.scale = 1.0;
  if (
         strcmp(argv[argi],"scale") == 0
     ) {
    argi++;
    s.scale = GetParam(argc, argv, argi++, "scale");
  };


  s.shape                = GetParam(argc, argv, argi++, "shape");

  s.origin[X]            = GetParam(argc, argv, argi++, "origin[X]");
  s.origin[Y]            = GetParam(argc, argv, argi++, "origin[Y]");
  s.origin[Z]            = GetParam(argc, argv, argi++, "origin[Z]");
  fprintf(stderr,"\n");

  s.basis[SURFACE][X]    = GetParam(argc, argv, argi++, "surface[X]");
  s.basis[SURFACE][Y]    = GetParam(argc, argv, argi++, "surface[Y]");
  s.basis[SURFACE][Z]    = GetParam(argc, argv, argi++, "surface[Z]");
  fprintf(stderr,"\n");
  s.basis[TANGENT][X]    = GetParam(argc, argv, argi++, "tangent[X]");
  s.basis[TANGENT][Y]    = GetParam(argc, argv, argi++, "tangent[Y]");
  s.basis[TANGENT][Z]    = GetParam(argc, argv, argi++, "tangent[Z]");
  fprintf(stderr,"\n");
  s.basis[AXIS][X]       = GetParam(argc, argv, argi++, "axis[X]");
  s.basis[AXIS][Y]       = GetParam(argc, argv, argi++, "axis[Y]");
  s.basis[AXIS][Z]       = GetParam(argc, argv, argi++, "axis[Z]");
  fprintf(stderr,"\n");

/*  Compute "cross" just for the interest of the user.  Don't actually
 *  use it in internal computations.
 */
  for (i=0;i<DIMENSION;i++)
         CrossProduct(cross[i],
                      s.basis[(i+1)%DIMENSION],
                      s.basis[(i+2)%DIMENSION]);
  fprintf(stderr," The surface, tangent, and axis vectors\n");
  fprintf(stderr," along with their reduced integerized cross product matrix\n");
  for (i=0;i<DIMENSION;i++) {
    for (j=0;j<DIMENSION;j++) longCross[i][j] = 16.0*cross[i][j];
    for (j=0;j<DIMENSION;j++)
       if (longCross[i][j] != (16.0*cross[i][j]) )
         fprintf(stderr, "GCD won't work, values aren't exact integers:  %ld %10f.\n",
           longCross[i][j], cross[i][j] );
    k = Gcd(Gcd(longCross[i][X],longCross[i][Y]),longCross[i][Z]);
    for (j=0;j<DIMENSION;j++) longCross[i][j] /= k;
  };

  for (i=0;i<DIMENSION;i++) {
    for (j=0;j<DIMENSION;j++) {
       fprintf(stderr," %8.2f",s.basis[i][j]);
    };
    fprintf(stderr,"\n");
    for (j=0;j<DIMENSION;j++) fprintf(stderr," %5ld   ",  longCross[i][j]);
    fprintf(stderr,"\n\n");
  };

  /*  Now read in the optional "grid point parameters" */
  s.gridCount = 0;
  while (  (argi < argc) && ((strcmp(argv[argi],"grid") == 0))) {
    argi++;
    s.grid[s.gridCount].edge[1][X] = GetParam(argc, argv, argi++, "edge[1][X]");
    s.grid[s.gridCount].edge[1][Y] = GetParam(argc, argv, argi++, "edge[1][Y]");
    s.grid[s.gridCount].edge[1][Z] = GetParam(argc, argv, argi++, "edge[1][Z]");
    fprintf(stderr,"\n");
    s.grid[s.gridCount].edge[2][X] = GetParam(argc, argv, argi++, "edge[2][X]");
    s.grid[s.gridCount].edge[2][Y] = GetParam(argc, argv, argi++, "edge[2][Y]");
    s.grid[s.gridCount].edge[2][Z] = GetParam(argc, argv, argi++, "edge[2][Z]");
    fprintf(stderr,"\n");
    /*  Compute edge[0], which is always at right angles
     *  to edge[1] and edge[2].  Usually, edge[0] is parallel
     *  to the surface vector (though this is not required)
     *  and edges 1 and 2 are parallel to the plain of the surface.
     */
    CrossProduct(s.grid[s.gridCount].edge[0],
                 s.grid[s.gridCount].edge[1],
                 s.grid[s.gridCount].edge[2]);

/*  Try to give some useful information about whether or not this
 *  series of grid points will or will not "match up" around the
 *  tube or ring.
 */
    VectorCopy(s.grid[s.gridCount].point1, zeroVector);
    VectorCopy(s.grid[s.gridCount].point2, zeroVector);
    VectorCopy(tempVector,s.basis[TANGENT]);
    if ( (s.shape == RINGSHAPE) || (s.shape == TUBESHAPE) ) {
      if (InBox(tempVector,s.grid[s.gridCount]) != TRUE)
           fprintf(stderr," Grid not compatible with tangent of ring/tube\n");
      if (s.shape == RINGSHAPE) {
        VectorCopy(tempVector,s.basis[AXIS]);
        if (InBox(tempVector,s.grid[s.gridCount]) != TRUE)
           fprintf(stderr," Grid not compatible with axis of ring\n");
      };
    };

/*
 *  Now read in the "boxes"
 */

    lastV1V2 = s.gridCount;
    while  (argi < argc) {
           if (strcmp(argv[argi],"delete") == 0) {
           fprintf(stderr,"Deleting atoms:\n");
           argi++;
           s.grid[s.gridCount].action    = DELETEACTION;
      }
      else if (strcmp(argv[argi],"change") == 0) {
           fprintf(stderr,"Changing atoms:\n");
           argi++;
           s.grid[s.gridCount].action    = CHANGEACTION;
           GetStringParam(argc, argv, argi++,
                       s.grid[s.gridCount].fromName, MAXNAMELENGTH, "Change from");
           ExpectStringParam(argc, argv, argi++, "to");
           GetStringParam(argc, argv, argi++,
                       s.grid[s.gridCount].toName,   MAXNAMELENGTH, "Change to");
           s.grid[s.gridCount].toBondsWanted = 0;
           if (strcmp(argv[argi],"bondswanted") == 0) {
              argi++;
              s.grid[s.gridCount].toBondsWanted = GetParam(argc, argv, argi++,
                       "bonds wanted");
           };
      }
      else if (strcmp(argv[argi],"add") == 0) {
           fprintf(stderr,"Adding atoms:\n");
           argi++;
           s.grid[s.gridCount].action    = ADDACTION;
      }
      else break;
      s.grid[s.gridCount].point1[X] = GetParam(argc, argv, argi++, "point1[X]");
      s.grid[s.gridCount].point1[Y] = GetParam(argc, argv, argi++, "point1[Y]");
      s.grid[s.gridCount].point1[Z] = GetParam(argc, argv, argi++, "point1[Z]");
      fprintf(stderr,"\n");
      if ((argi <argc) && (strcmp(argv[argi],"box")) == 0) {
         argi++;
         s.grid[s.gridCount].point2[X] = GetParam(argc, argv, argi++, "point2[X]");
         s.grid[s.gridCount].point2[Y] = GetParam(argc, argv, argi++, "point2[Y]");
         s.grid[s.gridCount].point2[Z] = GetParam(argc, argv, argi++, "point2[Z]");
         fprintf(stderr,"\n");
      }
      else
        VectorCopy(s.grid[s.gridCount].point2, s.grid[s.gridCount].point1);
      for (i=0;i<DIMENSION;i++)
          VectorCopy(s.grid[s.gridCount].edge[i], s.grid[lastV1V2].edge[i]);
      fprintf(stderr,"\n");
      if (++s.gridCount >= MAXGRIDPOINTS)
          ErrAbort("Too many grid points, increase MAXGRIDPOINTS and recompile\n");
    };
  };

  if (argi != argc) {
     fprintf(stderr," Only %d command line parameters were read.\n",argi-1);
     ErrAbort("Change parameters and try again");
  };

  if(Dot(s.basis[TANGENT],s.basis[SURFACE]) != 0.0) {
       fprintf(stderr,"The tangent and surface vectors are not orthogonal\n");
       fprintf(stderr,"The dot product is: %15f\n",
             Dot(s.basis[TANGENT], s.basis[SURFACE]));
       ErrAbort("Change vectors and try again.");
  };
  if(
     Dot(s.basis[AXIS],s.basis[SURFACE]) != 0.0
    ) {
        CrossProduct(tempCoord, s.basis[SURFACE], s.basis[TANGENT]);
        fprintf(stderr, "The axis is not orthogonal to the surface vector\n");
        ErrAbort("Change vectors and try again.");
  };

/*  Yes, this is a hack.  Print out the internal coordinates of each atom.
 *  This is useful in specifying the coordinates of atoms that are to be
 *  deleted or changed.  A graphical user interface would obviously be\
 *  much superior.  In the meantime, from the atom number it's possible
 *  to determine the atom coordinates, which can be used to modify
 *  the atom.
 */
  streamOut = fopen("coords","w");


/*  print some headers to the output file */
  if (outputFormat == POLYGRAF) {
     printf("BIOGRF 160\n");
     printf("DESCRP %s\n", shapeList[s.shape]);
  }
  else if (outputFormat == BROOKHAVEN) {
     printf("HEADER Brookhaven                                             BGRF             \n");
     printf("REMARK %s\n", shapeList[s.shape]);
  }
  else ErrAbort("variable \"outputFormat\" has illegal value");

  printf("REMARK Version 1.02 of tube generation program\n");
  printf("REMARK  SX= %15f SY= %15f SZ= %15f\n",
             s.basis[SURFACE][X], s.basis[SURFACE][Y], s.basis[SURFACE][Z]);
  printf("REMARK  TX= %15f TY= %15f TZ= %15f\n",
             s.basis[TANGENT][X], s.basis[TANGENT][Y], s.basis[TANGENT][Z]);
  printf("REMARK  AX= %15f AY= %15f AZ= %15f\n",
             s.basis[AXIS][X], s.basis[AXIS][Y], s.basis[AXIS][Z]);

  printf("REMARK  surface*tangent = %15f\n",
             Dot(s.basis[SURFACE],s.basis[TANGENT]));
  printf("REMARK  surface*axis    = %15f\n",
             Dot(s.basis[SURFACE],s.basis[AXIS]));
  printf("REMARK  axis*tangent    = %15f\n",
             Dot(s.basis[AXIS],s.basis[TANGENT]));

/* processNext is the index of the atom to be processed next.
 * Notice that firstAtom <= processNext <= lastAtomPlusOne.
 * If processNext = lastAtomPlusOne, there are no more atoms
 * to process.
 * To start, processNext and lastAtomPlusOne are all set
 * equal to first atom, because there are no atoms in the
 * list and no atoms to process.
 */
   firstAtom = 1;
   processNext = firstAtom;
   lastAtomPlusOne = firstAtom;

/* Initialize the hash table */
   InitializeHashTable();

/* Add the seed coordinate */
   seed[X] = s.origin[X];
   seed[Y] = s.origin[Y];
   seed[Z] = s.origin[Z];
   seedNum = AddAtom(seed, EVEN);
   seedType = atomType[seedNum];

   printf("REMARK  seed   = %15f, %15f, %15f\n",seed[X], seed[Y], seed[Z]);
   printf("REMARK  seedType = %d\n", seedType);
   if (TestCoordinates(seed, s)>0)
      ErrAbort("Could not create seed atom within the region.");

/* Loop until the atoms processed equals the atoms in the array */
  while (processNext < lastAtomPlusOne) {
    for(currentBond=0; currentBond<VALENCE; currentBond++) {
      /*  loop over each neighbor of atom "processNext" */
      MakeNeighbor(processNext, currentBond, tempCoord, &tempType);
      MapCoordinates(tempCoord, s);
      /*  See if the neighbor atom is in the region we're creating */
      if (TestCoordinates(tempCoord, s)<0) {
      /*  See if the neighbor atom is already there. */
        tempAtom = AlreadyThere(tempCoord);
        if (tempAtom == NILATOM) {
           /*  The neighbor is NOT already there.  Add it */
           tempAtom = AddAtom(tempCoord, tempType);
        };
        /*  Bond the neighbor to the current atom. */
        BondAtoms(processNext, tempAtom);
      };
    };
    processNext++;
  };

/*  Loop over all the atoms and determine the element type for each atom.
 */
  limitAtom = processNext;
  for(currentAtom=firstAtom; currentAtom<limitAtom; currentAtom++) {
         if (bondCount[currentAtom] == 4) {
           strcpy(name[currentAtom],"C_3");
           bondsWanted[currentAtom] = 4;
         }
    else if (bondCount[currentAtom] == 3) {
           strcpy(name[currentAtom],"N_3");
           bondsWanted[currentAtom] = 3;
         }
    else if (bondCount[currentAtom] == 2) {
           strcpy(name[currentAtom],"O_3");
           bondsWanted[currentAtom] = 2;
         }
    else if (bondCount[currentAtom] == 1) {
           strcpy(name[currentAtom],"H_");
           bondsWanted[currentAtom] = 1;
         }
    else fprintf(stderr," ERROR: atom %d is bonded to %d other atoms\n",
           currentAtom,bondCount[currentAtom]);

/*  See if it's a grid point with a "change" command
 *  and if it is, change from the "fromName" to the "toName"
 */
    for(i=0; i<s.gridCount; i++) {
       if (  InBox(coord[currentAtom],s.grid[i])==TRUE  )
         if (s.grid[i].action == CHANGEACTION) {
            if (  (strcmp(s.grid[i].fromName,name[currentAtom]) == 0) ||
                  (strcmp(s.grid[i].fromName,"any"              ) == 0)
               ) {
              strcpy(name[currentAtom],s.grid[i].toName);
              if (s.grid[i].toBondsWanted != 0)
                bondsWanted[currentAtom] = s.grid[i].toBondsWanted;
            }
         };
    };

/*  Add hydrogens when there are dangling bonds.
 */
    if (bondCount[currentAtom] < bondsWanted[currentAtom] && addH)
       /*  Under coordinated atom.  Add hydrogen  */
       for(currentBond=0;currentBond<VALENCE;currentBond++) {
         MakeNeighbor(currentAtom, currentBond, tempCoord, &tempType);
         MapCoordinates(tempCoord, s);
      
         /*  See if the neighbor atom is already there. */
         tempAtom = AlreadyThere(tempCoord);
         if (tempAtom == NILATOM) {
              /*  No neighbor.  Add a hydrogen.  */
              /*  don't want mapped coordinate, so redo MakeNeighbor   */
              MakeNeighbor(currentAtom, currentBond, tempCoord, &tempType);
              /*  move the hydrogen atom closer */
              for(i=0;i<DIMENSION;i++)
                 tempCoord[i] = 0.7*tempCoord[i]+0.3*coord[currentAtom][i];
              tempAtom = AddAtom(tempCoord, tempType);
              BondAtoms(currentAtom, tempAtom);
              strcpy(name[tempAtom],"H_");
              bondsWanted[tempAtom] = 1;
           };
       };
  };  /*  end of loop over all atoms */

/*  We've now created the data structure describing the region of diamond
 *  that we're interested in.  Print out the data in PolyGraf format.
 */

/*  The PolyGraf file format was deduced by simply taking an existing PolyGraf
 *  file and making the obvious assumptions about which fields meant what.  Some
 *  of the fields which have no obvious immediate applicability to this program
 *  (the charge on the atom, various mysterious integer values, residue labels,
 *  etc) have been left in the format string.
 */
  printf("REMARK This structure was generated by a computer program.\n");
  printf("REMARK It may or may not be chemically reasonable.\n");
  if (outputFormat == POLYGRAF) {
    printf("FORCEFIELD DREIDING\n");
    printf("FORMAT ATOM   (a6,1x,i5,1x,a5,1x,a3,1x,a1,1x,a5,3f10.5,1x,a5,i3,i2,1x,f8.5)\n");
  };
  for (j=firstAtom; j<lastAtomPlusOne; j++) {
    fprintf(streamOut, "%8d, %8.2f, %8.2f, %8.2f\n", j,
        coord[j][X], coord[j][Y], coord[j][Z]);
    Bend(tempCoord, coord[j], s);
    if (outputFormat == POLYGRAF)
       printf("HETATM %5d %5s ORG B    1 %10f%10f%10f %5s%3d 0  0.00000\n",
          j,
          AtomLabel(j,s,POLYGRAFNAMELENGTH),
          tempCoord[X],
          tempCoord[Y],
          tempCoord[Z],
          AtomName(j,s,POLYGRAFNAMELENGTH),
          AtomWantsThisManyBonds(j)
        );
    else if (outputFormat == BROOKHAVEN)
       printf("HETATM%5d %4s ORG B   1    %8.3f%8.3f%8.3f\n",
          j,
          AtomName(j,s,BROOKHAVENNAMELENGTH),
          tempCoord[X],
          tempCoord[Y],
          tempCoord[Z]
        );
    else ErrAbort("Illegal value for outputFormat");
  };
  if (outputFormat == POLYGRAF)
     printf("%s\n","FORMAT CONECT (a6,12i6)");
  for (j=firstAtom; j<lastAtomPlusOne; j++) {
    if (outputFormat == POLYGRAF) {
       printf("CONECT%6d", j);
       for(k=0; k<bondCount[j]; k++) printf("%6d", bond[j][k]);
       printf("\n");
    }
    else if (outputFormat == BROOKHAVEN) {
       printf("CONECT%5d", j);
       for(k=0; k<bondCount[j]; k++) {
           if ( (k%4 == 0) && (k>0) ) printf("\nCONECT%5d", j);
         printf("%5d", bond[j][k]);
       };
       printf("\n");
    }
    else ErrAbort("Illegal outputFormat value");
  };
  printf("%s\n","END");
  return(0);

}

/*
 *  Project p1, p2, p3, and p4 along the vector v.  If
 *  the resulting values are "linear", e.g., fall in
 *  sequential order or reverse sequential order,
 *  return TRUE.
 *  If the values can be MADE to fall in sequential order
 *  by adding or subtracting an integer multiple of the
 *  projection of p4 against 4, then also return TRUE.
 *  If this doesn't work return FALSE.
 */
ModInOrder(v,p1,p2,p3,p4)
double v[DIMENSION];
double p1[DIMENSION];
double p2[DIMENSION];
double p3[DIMENSION];
double p4[DIMENSION];
{
   double a1,a2,a3,a4;
   int    n;
/*  project all four points onto the vector v */
   a1 = Dot(v,p1);
   a2 = Dot(v,p2);
   a3 = Dot(v,p3);
   a4 = Dot(v,p4);
/*  make sure the "modulus" value is positive, a negative
 *  value would foul things up later
 */
   if (a4 < 0.0) a4 = -a4;
/*  Compute the magic integer "n" */
   n = 0;
   if (a4 != 0.0) {
    if (a1 <= a3)
       n = floor((a2-a1)/a4);
    else
       n = floor((a2-a3)/a4);
   };

/*  put a2 between a1 and a3, if possible, by adding a multiple of a4.  */
   a2 -= n*a4;
   if (InOrder(a1,a2,a3)) return(TRUE);
/*  A little paranoid checking, just in case...*/
   a2 -= a4;
   if (InOrder(a1,a2,a3)) ErrAbort("Logic Error, n off by one");
   a2 += 2*a4;
   if (InOrder(a1,a2,a3)) ErrAbort("Logic Error, n off by one");
   return(FALSE);
}

/*
 *  Check to see if the point p lies within the "bounding box"
 *  defined by grid.point1 and grid.point2.  The box edges are
 *  defined by edge[1], edge[2], and edge[1] X edge[2].  ("X" is the cross
 *  product).  We also check to see if p lies within boxes
 *  that are periodically offset from the "base box" by
 *  increments of the edge vectors.  Of course, edge[0] was
 *  a fake, generated by taking the cross product of edge[1]
 *  and edge[2], and we don't really want periodicity in
 *  that direction at all.  This leads to some special-case
 *  handling of edge[0] to prevent the undesired periodicity
 *  from occuring along that axis.
 */
InBox(p, grid)
double p[DIMENSION];
struct gridRecord grid;
{
   double cross[DIMENSION][DIMENSION];
   double temp [DIMENSION][DIMENSION];
   static double zeroVector[DIMENSION] = {0.0, 0.0, 0.0};
   int i;

   for (i=0;i<DIMENSION;i++)
      CrossProduct(cross[i],grid.edge[(i+1)%DIMENSION],
                            grid.edge[(i+2)%DIMENSION]);
/*  We almost don't need the "temp" matrix, but require it
 *  because temp[0] has to be the 0 vector, rather than
 *  edge[0].  In some sense, we want the edge[0] vector
 *  to be infinitely long, so that the "period" in that
 *  direction is also infinitely long.  We can't do that
 *  in C, so we have to "special case" it.
 */
   VectorCopy(temp[0],zeroVector);
   VectorCopy(temp[1],grid.edge[1]);
   VectorCopy(temp[2],grid.edge[2]);
/*  For each dimension, see if the given point falls outside the
 *  bounding box in that dimension
 */
   for (i=0;i<DIMENSION;i++)
      if (ModInOrder(cross[i],grid.point1, p, grid.point2, temp[i])   ==FALSE)
         return(FALSE);
   return(TRUE);
}

MakeNeighbor(atom, bondNumber, newCoord, newType)
int atom;
int bondNumber;
double newCoord[DIMENSION];
int    *newType;
{
/*  This algorithm generates the 4 neighbors of an atom in a diamond
 *  lattice.  There are two types of atoms in a diamond lattice,
 *  arbitrarily called "odd" and "even."  The prototypical even
 *  site is at 0,0,0.  The prototypical odd site is 1,1,1
 *  The four neighbors for an even site are at:
 *  coords + [ 1/4,  1/4,  1/4]
 *  coords + [-1/4, -1/4,  1/4]
 *  coords + [ 1/4, -1/4, -1/4]
 *  coords + [-1/4,  1/4, -1/4]
 *  The four neighbors for an odd site are at the coords - the
 *  deltas given above, instead of +.
 *
 *  WARNING:  The floating point representation of the lattice points
 *  is exact.  This is deliberate.  The code will produce subtle
 *  errors if the floating point representation of the lattice
 *  points cannot be represented exactly in a floating point number.
 *  There might also be problems when using floating point operations
 *  that are not as neat and elegant as the IEEE floating point
 *  operations.
 *  newCoord and newType are output variables, and are modified by
 *  this program.
 *
 */
   int dim;
   if(bondNumber < 0) ErrAbort("Logic Error bondNumber < 0");
   if(bondNumber >= VALENCE) ErrAbort("Logic Error bondNumber >= VALENCE");
   for (dim=0; dim<DIMENSION; dim++) {
     if(atomType[atom]==EVEN)
         newCoord[dim] = coord[atom][dim] + delta[bondNumber][dim];
     else
         newCoord[dim] = coord[atom][dim] - delta[bondNumber][dim];
   };
   /*  The new neighbor has an atom type opposite to the current
    *  atom type.
    */
        if (atomType[atom]==EVEN) *newType = ODD;
   else if (atomType[atom]==ODD)  *newType = EVEN;
   else ErrAbort("Bad atom type.");
}

/*
 *  This subroutine determines whether an atom with the given coordinates
 *  is already present.  If it is, it returns the atom number.  If it
 *  isn't, it returns NILATOM.
 *  It uses the hash table set up by AddAtom.  Note that the hashing
 *  method used here preserves spatial relationships.  If coordinates
 *  C1 and C2 are close to each other in physical space, then their
 *  integer hash indices will differ by at most 1.
 *
 */
AlreadyThere(atomCoord)
double   atomCoord[DIMENSION];
{ int xLower, xUpper, yLower, yUpper, zLower, zUpper;
  int i,j,k;
  int tempAtom;
  int atomToReturn;
  xLower = HashIndex(atomCoord[X]-EPSILON,XHASHSIZE);
  xUpper = HashIndex(atomCoord[X]+EPSILON,XHASHSIZE)+1;
  yLower = HashIndex(atomCoord[Y]-EPSILON,YHASHSIZE);
  yUpper = HashIndex(atomCoord[Y]+EPSILON,YHASHSIZE)+1;
  zLower = HashIndex(atomCoord[Z]-EPSILON,ZHASHSIZE);
  zUpper = HashIndex(atomCoord[Z]+EPSILON,ZHASHSIZE)+1;
  xUpper %= XHASHSIZE; /*  just in case xUpper = XHASHSIZE */
  yUpper %= YHASHSIZE;
  zUpper %= ZHASHSIZE;
  atomToReturn = NILATOM;
  for(i=xLower; i != xUpper; i = (i+1)%XHASHSIZE)
  for(j=yLower; j != yUpper; j = (j+1)%YHASHSIZE)
  for(k=zLower; k != zUpper; k = (k+1)%ZHASHSIZE) {
    tempAtom = LookInHashBucket(i,j,k,atomCoord);
    if (atomToReturn == NILATOM)
       atomToReturn = tempAtom;
    else if (tempAtom != NILATOM)
       fprintf(stderr,"Atoms %d and %d too close near: (%f, %f, %f)\n",
              atomToReturn, tempAtom, atomCoord[X], atomCoord[Y], atomCoord[Z]);
  };
  return(atomToReturn);
}

/*  The following routine returns a simple integer that is either
 *  positive or negative.  (0 is never returned).  By modifying
 *  the code in this routine you can easily modify the shape of
 *  the generated block of diamond.  If you want a block of diamond
 *  shaped like an elephant, just make the appropriate modifications
 *  here.  Note that the "seed" in the main program must fall inside
 *  the region being created.  If you define a region that does
 *  not include 0,0,0, then you have to modify the seed appropriately.
 *  More obscurely, at least one neighbor of the seed should also lie
 *  within the region.  Sometimes the seed lies on the boundary of
 *  the region, but none of its neighbors fall inside the region.
 *  The value EXTERIOR2 can be used to indicate that some boundary
 *  is be treated in a special fashion.  By changing which test
 *  returns EXTERIOR2 (and by altering the sequence of testing)
 *  it is possible to indicate in a very flexible manner which
 *  atoms should receive special handling.  In general, it is easy
 *  to have this return return information about the region in which
 *  that atom falls by adding new variables EXTERIOR3, EXTERIOR4,...
 *  and INTERIOR2, INTERIOR3, ...
 */
TestCoordinates(atomCoord, s)
double atomCoord[DIMENSION];
struct shapeDef s;
{
 int    i;
    VectorPrint(stdout,"atom coord", atomCoord);

/*  See if it's a grid point */
    for(i=0; i<s.gridCount; i++)
       if (InBox(atomCoord,s.grid[i])==TRUE) {
         if (s.grid[i].action == DELETEACTION) return(EXTERIOR1);
         if (s.grid[i].action == ADDACTION)    return(INTERIOR1);
         if (s.grid[i].action != CHANGEACTION)
              ErrAbort("Action must be delete, add, or change");
       };

    if (Dot(atomCoord,s.basis[SURFACE]) < 0) return(EXTERIOR1);
    if (Dot(atomCoord,s.basis[SURFACE]) >= Dot(s.basis[SURFACE], s.basis[SURFACE]))
                    return(EXTERIOR1);
    if (s.shape == RINGSHAPE) return(INTERIOR1);
    if (Dot(atomCoord,s.basis[AXIS]) < 0)    return(EXTERIOR1);
    if (Dot(atomCoord,s.basis[AXIS]) >= Dot(s.basis[AXIS],s.basis[AXIS]))
                    return(EXTERIOR1);
    if (s.shape == TUBESHAPE) return(INTERIOR1);
    if (Dot(atomCoord,s.basis[TANGENT]) < 0) return(EXTERIOR1);
    if (Dot(atomCoord,s.basis[TANGENT]) >= Dot(s.basis[TANGENT],s.basis[TANGENT]))
                    return(EXTERIOR1);
    if (s.shape != BLOCKSHAPE) ErrAbort("Bad shape in TestCoordinates");

    return(INTERIOR1);
}

/*  a is set to the cross product of b and c.
 *  By definition, this means that a is orthogonal to
 *  both b and c, and that
 *  Length(a) = Length(b) * Length(c) * sin(theta)
 *  where theta is the angle from b to c.
 */
CrossProduct(a,b,c)
double a[DIMENSION], b[DIMENSION], c[DIMENSION];
{
/*  use a temporary just in case a and b or c overlap in memory */
double t[DIMENSION];
   /*  Cross products only work in three dimensions anyway,
    *  so this routine is locked into DIMENSION = 3.
    */
   t[X] =  b[Y]*c[Z]-c[Y]*b[Z];
   t[Y] =  c[X]*b[Z]-b[X]*c[Z];
   t[Z] =  b[X]*c[Y]-c[X]*b[Y];
   a[X] = t[X];
   a[Y] = t[Y];
   a[Z] = t[Z];
}

/*  This routine maps atoms around into a tube.
 *  That is, the mapping of a rectangular object
 *  into a tube involves connecting two opposing
 *  faces of the object.  This routine performs
 *  that connection by mapping atoms that "fall
 *  off" the rectangle along the X axis back into
 *  the region occupied by the rectangle, but
 *  starting over again near X = 0.  This routine
 *  is essential for getting the bonds around the tube
 *  edge correct.
 */
MapCoordinates(atomCoord, s)
double atomCoord[DIMENSION];
struct shapeDef s;
{
double orthoTangent[DIMENSION];
double orthoAxis[DIMENSION];
double tempVector[DIMENSION];
double dot1, dot2;
double k;
 if (s.shape < 0) ErrAbort("Bad shape parameter: negative.");
 if (s.shape >= MAXSHAPE) ErrAbort("Bad shape parameter: too large.");
 if (s.shape == BLOCKSHAPE) return;

/*  turn it into a toroid */
 if (s.shape == RINGSHAPE) {
   CrossProduct(orthoTangent, s.basis[SURFACE], s.basis[TANGENT]);
   dot1 = Dot(atomCoord,orthoTangent);
   dot2 = Dot(s.basis[AXIS],   orthoTangent);
   k = floor(dot1/dot2);
   VectorCopy(tempVector, s.basis[AXIS]);
   ScalarMul (tempVector, k);
   VectorSub (atomCoord,atomCoord,tempVector);
 };

 if (s.shape == RINGSHAPE || s.shape == TUBESHAPE) {
   CrossProduct(orthoAxis, s.basis[AXIS],s.basis[SURFACE]);
   dot1 = Dot(atomCoord,orthoAxis);
   dot2 = Dot(s.basis[TANGENT],   orthoAxis);
   k = floor(dot1/dot2);
   VectorCopy(tempVector, s.basis[TANGENT]);
   ScalarMul (tempVector, k);
   VectorSub (atomCoord,atomCoord,tempVector);
 };
}

/*  This routine bends the nice neat crystal
 *  coordinates that describe a nice neat diamond
 *  crystal in three space into the bent or strained
 *  structure of a tube or toroid.
 */
Bend(atomCoordOut, atomCoordIn, s)
double atomCoordOut[DIMENSION];
double atomCoordIn [DIMENSION];
struct shapeDef s;
{ double radius1, radius2, angle1, angle2;
  double orthoTangent[DIMENSION];
  double orthoAxis[DIMENSION];
  int i;

  VectorCopy(atomCoordOut,atomCoordIn);
  if (s.shape == BLOCKSHAPE) goto adjustLatticeConstant;
  CrossProduct(orthoTangent, s.basis[TANGENT], s.basis[SURFACE]);
  CrossProduct(orthoAxis,    s.basis[AXIS],    s.basis[SURFACE]);
  angle1  = (Dot(atomCoordIn,orthoAxis)/Dot(s.basis[TANGENT],orthoAxis))*TWOPI;
  radius1 = Length(s.basis[TANGENT])/TWOPI +
             Dot(atomCoordIn,s.basis[SURFACE])/Length(s.basis[SURFACE]);
  atomCoordOut[X]  = sin(angle1)*radius1;
  atomCoordOut[Y]  = Dot(atomCoordIn,orthoTangent)/Dot(s.basis[AXIS],orthoTangent)
                         * Length(s.basis[AXIS]);
  atomCoordOut[Z]  = cos(angle1)*radius1;

  if(s.shape == TUBESHAPE) goto adjustLatticeConstant;
/*  We've finished bending it into a tube.  Now bend it into a toroid. */
  angle2  = (Dot(atomCoordIn,orthoTangent)/
                     Dot(s.basis[AXIS],orthoTangent))*TWOPI;
  radius2 = Length(s.basis[AXIS])/TWOPI + cos(angle1)*radius1;
  atomCoordOut[X]  = sin(angle2)*radius2;
  atomCoordOut[Y]  = sin(angle1)*radius1;
  atomCoordOut[Z]  = cos(angle2)*radius2;

  if (s.shape != RINGSHAPE) ErrAbort("Bad shape parameter in Bend");

adjustLatticeConstant:
  for(i=0;i<DIMENSION;i++) atomCoordOut[i] *=
         LATTICECONSTANT*s.scale;
}


/* This routine orthogonalizes vector a using vector b */
orthogonalize(a,b)
double a[DIMENSION];
double b[DIMENSION];
{int i;
double ftemp;
    ftemp = Dot(a,b)/(Length(b)*Length(b));
    for(i=0;i<DIMENSION;i++) a[i] -= ftemp*b[i];
}