OOPSE/utils/quickLate.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include <sys/types.h>
#include <sys/stat.h>

inline double roundMe( double x ){
  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
}

struct coords{
  double pos[3]; // cartesian coords
  double q[4];  // the quanternions
  char name[30];
};

struct linked_xyz{
  struct coords *r;
  double time;
  double Hmat[3][3];
  double HmatI[3][3];
  struct linked_xyz *next;
};

char *program_name; /*the name of the program */
int printWater = 0;
int printDipole = 0;

void usage(void);

void rotWrite( FILE* out_file, struct coords* theSSD );

void setRot( double A[3][3], double q[4] );

void rotMe( double A[3][3], double r[3] );
void wrapVector( double thePos[3], double Hmat[3][3], double HmatI[3][3]);
double matDet3(double a[3][3]);
void invertMat3(double a[3][3], double b[3][3]);
void matVecMul3(double m[3][3], double inVec[3], double outVec[3]);

int main(argc, argv)
     int argc;
     char *argv[];
{


  struct coords *out_coords;

  int i, j, k, l, m; /* loop counters */
  mode_t dir_mode = S_IRWXU;

  int lineCount = 0;  //the line number
  int newN; // the new number of atoms
  int nSSD=0; // the number of SSD per frame

  unsigned int nAtoms; /*the number of atoms in each time step */
  char read_buffer[1000]; /*the line buffer for reading */
  char *eof_test; /*ptr to see when we reach the end of the file */
  char *foo; /*the pointer to the current string token */
  FILE *in_file; /* the input file */
  FILE *out_file; /*the output file */
  char *out_prefix = NULL; /*the prefix of the output file */
  char out_name[500]; /*the output name */
  int have_outName =0;
  char in_name[500]; // the input file name
  char temp_name[500];
  char *root_name = NULL; /*the root name of the input file */
  unsigned int n_out = 0; /*keeps track of which output file is being written*/


  int skipFrames = 1;

  struct linked_xyz *current_frame;
  struct linked_xyz *temp_frame;

  int nframes =0;

  int wrap = 0;

  double cX = 0.0;
  double cY = 0.0;
  double cZ = 0.0;

  double mcX, mcY, mcZ;
  double newCX, newCY, newCZ;
  double dx, dy, dz;
  int mcount;

  int repeatX = 0;
  int repeatY = 0;
  int repeatZ = 0;

  int nRepeatX = 0;
  int nRepeatY = 0;
  int nRepeatZ = 0;

  struct coords* rCopy;
  
  int done;
  char current_flag;

  program_name = argv[0]; /*save the program name in case we need it*/
  
  for( i = 1; i < argc; i++){
    
    if(argv[i][0] =='-'){
      
      if(argv[i][1] == '-' ){
        
        // parse long word options
        
        if( !strcmp( argv[i], "--repeatX" ) ){
          repeatX = 1;
          i++;
          nRepeatX = atoi( argv[i] );
        }

        else if( !strcmp( argv[i], "--repeatY" ) ){
          repeatY = 1;
          i++;
          nRepeatY = atoi( argv[i] );
        }

        else if( !strcmp( argv[i], "--repeatZ" ) ){
          repeatZ = 1;
          i++;
          nRepeatZ = atoi( argv[i] );
        }

        else{
          fprintf( stderr, 
                   "Invalid option \"%s\"\n", argv[i] );
          usage();
        }
      }
      
      else{

        // parse single character options
        
        done =0;
        j = 1;
        current_flag = argv[i][j];
        while( (current_flag != '\0') && (!done) ){
          
          switch(current_flag){
            
          case 'o':
            // -o <prefix> => the output
            
            i++;
            strcpy( out_name, argv[i] );
            have_outName = 1;
            done = 1;
            break;
            
          case 'n':
            // -n <#> => print every <#> frames
            
            i++;
            skipFrames = atoi(argv[i]);
            done = 1;
            break;
            
          case 'h':
            // -h => give the usage
            
            usage();
            break;
            
          case 'd':
            // print the dipoles
            
            printDipole = 1;
            break;
            
          case 'w':
            // print the waters
            
            printWater = 1;
            break;
            
          case 'm':
            // map to a periodic box
            
            wrap = 1;
            break;
            
          default:
            
            (void)fprintf(stderr, "Bad option \"-%s\"\n", current_flag);
            usage();
          }
          j++;
          current_flag = argv[i][j];
        }
      }
    }
    
    else{
      
      if( root_name != NULL ){
        fprintf( stderr, 
                 "Error at \"%s\", program does not currently support\n"
                 "more than one input file.\n"
                 "\n",
                 argv[i]);
        usage();
      }
      
      root_name = argv[i];
    }
  }
  
  if(root_name == NULL){
    usage();
  }

  sprintf( in_name, "%s", root_name );
  if( !have_outName ) sprintf( out_name, "%s.xyz", root_name );

  in_file = fopen(in_name, "r");
  if(in_file == NULL){
    printf("Cannot open file \"%s\" for reading.\n", in_name);
    exit(8);
  }
  
  out_file = fopen( out_name, "w" );
  if( out_file == NULL ){
    printf("Cannot open file \"%s\" for writing.\n", out_name);
    exit(8);
  }

  // start reading the first frame

  eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
  lineCount++;
  
  current_frame = (struct linked_xyz *)malloc(sizeof(struct linked_xyz));
  current_frame->next = NULL;
  

  while(eof_test != NULL){
    
    (void)sscanf(read_buffer, "%d", &nAtoms);
    current_frame->r = 
      (struct coords *)calloc(nAtoms, sizeof(struct coords));

    // read and the comment line and grab the time and box dimensions

    eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
    lineCount++;
    if(eof_test == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }

    foo = strtok( read_buffer, " ,;\t" );
    sscanf( read_buffer, "%lf", &current_frame->time );
    
    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[0][0]);

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[1][0]);

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[2][0]);

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[0][1]);

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[1][1]);

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[2][1]);

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[0][2]);

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[1][2]);

    foo = strtok(NULL, " ,;\t");
    if(foo == NULL){
      printf("error in reading file at line: %d\n", lineCount);
      exit(8);
    }
    (void)sscanf(foo, "%lf",&current_frame->Hmat[2][2]);


    // Find HmatI:

    invertMat3(current_frame->Hmat, current_frame->HmatI);


    for( i=0; i < nAtoms; i++){
      
      eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
      lineCount++;
      if(eof_test == NULL){
        printf("error in reading file at line: %d\n", lineCount);
        exit(8);
      }

      foo = strtok(read_buffer, " ,;\t");
      if ( !strcmp( "SSD", foo ) ) nSSD++;
      (void)strcpy(current_frame->r[i].name, foo); //copy the atom name 
                          
      // next we grab the positions 
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading postition x from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &current_frame->r[i].pos[0] );
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading postition y from %s\n"
             "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &current_frame->r[i].pos[1] );
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading postition z from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &current_frame->r[i].pos[2] );
    
      // get the velocities
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading velocity x from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      
    
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading velocity y from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading velocity z from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      
      // get the quaternions
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading quaternion 0 from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &current_frame->r[i].q[0] );
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading quaternion 1 from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &current_frame->r[i].q[1] );
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading quaternion 2 from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &current_frame->r[i].q[2] );

      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading quaternion 3 from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      (void)sscanf( foo, "%lf", &current_frame->r[i].q[3] );
      
      // get the angular velocities
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading angular momentum jx from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading angular momentum jy from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
      
      foo = strtok(NULL, " ,;\t");
      if(foo == NULL){
        printf("error in reading angular momentum jz from %s\n"
               "natoms  = %d, line = %d\n",
               in_name, nAtoms, lineCount );
        exit(8);
      }
    }

    // if necassary, wrap into the periodic image
    
    if(wrap){
      
      for( i=0; i<nAtoms; i++ ){
        
        /*
        if( !strcmp( current_frame->r[i].name, "HEAD" ) ){
          
          // find the center of the lipid
          
          mcX = 0.0;
          mcY = 0.0;
          mcZ = 0.0;
          mcount = 0;

          mcX += current_frame->r[i].x;
          mcY += current_frame->r[i].y;
          mcZ += current_frame->r[i].z;
          mcount++;
          i++;

          while( strcmp( current_frame->r[i].name, "HEAD" ) &&
                 strcmp( current_frame->r[i].name, "SSD" ) ){
            
            mcX += current_frame->r[i].x;
            mcY += current_frame->r[i].y;
            mcZ += current_frame->r[i].z;
            mcount++;
            i++;
          }
          
          mcX /= mcount;
          mcY /= mcount;
          mcZ /= mcount;

          newCX = mcX;
          newCY = mcY;
          newCZ = mcZ;
          
          map( &newCX, &newCY, &newCZ,
               cX, cY, cZ,
               current_frame->boxX,
               current_frame->boxY,
               current_frame->boxZ );
          
          dx = mcX - newCX;
          dy = mcY - newCY;
          dz = mcZ - newCZ;

          for(j=(i-mcount); j<mcount; j++ ){
            
            current_frame->r[j].x -= dx;
            current_frame->r[j].y -= dy;
            current_frame->r[j].z -= dz;
          }
          
          i--; // so we don't skip the next atom
        }
        */
        // else 
        
        wrapVector(current_frame->r[i].pos, 
                   current_frame->Hmat, 
                   current_frame->HmatI);
        
      }
    }   

    nframes++; 
    
    // write out here    
    
    if(printWater) newN = nAtoms + ( nSSD * 3 ); 
    else newN = nAtoms - nSSD;
    
    newN *= (nRepeatX+1);
    newN *= (nRepeatY+1);
    newN *= (nRepeatZ+1);

    if( !(nframes%skipFrames) ){
      fprintf( out_file,        
               "%d\n"
               "%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\n",
               newN, current_frame->time,
               current_frame->Hmat[0][0], current_frame->Hmat[1][0],
               current_frame->Hmat[2][0],
               current_frame->Hmat[0][1], current_frame->Hmat[1][1],
               current_frame->Hmat[2][1],
               current_frame->Hmat[0][2], current_frame->Hmat[1][2],
               current_frame->Hmat[2][2] );
    
      rCopy = (struct coords*)
        calloc( nAtoms, sizeof( struct coords ));

      for(i=0; i<nAtoms; i++){
        rCopy[i].q[0] = current_frame->r[i].q[0];
        rCopy[i].q[1] = current_frame->r[i].q[1];
        rCopy[i].q[2] = current_frame->r[i].q[2];
        rCopy[i].q[3] = current_frame->r[i].q[3];

        strcpy(rCopy[i].name, current_frame->r[i].name);
      }

      for(j=0; j<(nRepeatX+1); j++){
        for(k=0; k<(nRepeatY+1); k++){
          for(l=0; l<(nRepeatZ+1); l++){

            for(i=0; i<nAtoms; i++){

              for(m = 0; m < 3; m++) {
                rCopy[i].pos[m] = current_frame->r[i].pos[m] +
                  j * current_frame->Hmat[m][0] + 
                  k * current_frame->Hmat[m][1] + 
                  l * current_frame->Hmat[m][2];
              }
              
              if( !strcmp( "SSD", rCopy[i].name ) ){
                
                rotWrite( out_file, &rCopy[i] );
              }
              
              else if( !strcmp( "HEAD", rCopy[i].name ) ){
                
                rotWrite( out_file, &rCopy[i] );
              }
              else{
                //modify
                fprintf( out_file,
                         "%s\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
                         rCopy[i].name,
                         rCopy[i].pos[0],
                         rCopy[i].pos[1],
                         rCopy[i].pos[2] );
              }
            }
          }
        }
      }
      
      free( rCopy );
    }
    
    /*free up memory */
    
    temp_frame = current_frame->next;
    current_frame->next = NULL;
    
    if(temp_frame != NULL){
      
      free(temp_frame->r);
      free(temp_frame);
    }
    
    /* make a new frame */
    
    temp_frame = (struct linked_xyz *)malloc(sizeof(struct linked_xyz));
    temp_frame->next = current_frame;
    current_frame = temp_frame;
    
    nSSD = 0;

    eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
    lineCount++;
  }
  
  (void)fclose(in_file);
  
  return 0;
  
}


void rotWrite( FILE* out_file, struct coords* theDipole ){

  double u[3];
  double h1[3];
  double h2[3];
  double ox[3];
  
  double A[3][3];
  
 
  u[0] = 0.0; u[1] = 0.0; u[2] = 1.0;

  h1[0] = 0.0; h1[1] = -0.75695; h1[2] = 0.5206;

  h2[0] = 0.0; h2[1] = 0.75695;  h2[2] = 0.5206;

  ox[0] = 0.0; ox[1] = 0.0;      ox[2] = -0.0654;

  
  setRot( A, theDipole->q );
  rotMe( A, u );

  if( !strcmp( "SSD", theDipole->name )){
    
    rotMe( A, h1 );
    rotMe( A, h2 );
    rotMe( A, ox );

    if( printWater ){
      
      if(printDipole) {
        
        fprintf( out_file,
                 "X\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
                 theDipole->pos[0],
                 theDipole->pos[1],
                 theDipole->pos[2],
                 u[0], u[1], u[2] );

        fprintf( out_file,
                 "O\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
                 theDipole->pos[0] + ox[0],
                 theDipole->pos[1] + ox[1],
                 theDipole->pos[2] + ox[2] );
        
        fprintf( out_file,
                 "H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
                 theDipole->pos[0] + h1[0],
                 theDipole->pos[1] + h1[1],
                 theDipole->pos[2] + h1[2] );
        
        fprintf( out_file,
                 "H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
                 theDipole->pos[0] + h2[0],
                 theDipole->pos[1] + h2[1],
                 theDipole->pos[2] + h2[2] );
        }
      else{
        
        fprintf( out_file,
                 "X\t%lf\t%lf\t%lf\n",
                 theDipole->pos[0],
                 theDipole->pos[1],
                 theDipole->pos[2] );

        fprintf( out_file,
                 "O\t%lf\t%lf\t%lf\n",
                 theDipole->pos[0] + ox[0],
                 theDipole->pos[1] + ox[1],
                 theDipole->pos[2] + ox[2] );
        
        fprintf( out_file,
                 "H\t%lf\t%lf\t%lf\n",
                 theDipole->pos[0] + h1[0],
                 theDipole->pos[1] + h1[1],
                 theDipole->pos[2] + h1[2] );
        
        fprintf( out_file,
                 "H\t%lf\t%lf\t%lf\n",
                 theDipole->pos[0] + h2[0],
                 theDipole->pos[1] + h2[1],
                 theDipole->pos[2] + h2[2] );
      }     
    }
  }
  
  else if( !strcmp( "HEAD", theDipole->name )) {
    
    if( printDipole ){
      
      fprintf( out_file,
               "HEAD\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
               theDipole->pos[0],
               theDipole->pos[1],
               theDipole->pos[2],
               u[0], u[1], u[2] );
    }
    else{
     
      fprintf( out_file,
               "HEAD\t%lf\t%lf\t%lf\n",
               theDipole->pos[0],
               theDipole->pos[1],
               theDipole->pos[2] );
    }
  }
}


void setRot( double A[3][3], double the_q[4] ){

  double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
  
  q0Sqr = the_q[0] * the_q[0];
  q1Sqr = the_q[1] * the_q[1];
  q2Sqr = the_q[2] * the_q[2];
  q3Sqr = the_q[3] * the_q[3];
  
  A[0][0] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
  A[0][1] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
  A[0][2] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
  
  A[1][0]  = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
  A[1][1] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
  A[1][2] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
  
  A[2][0] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
  A[2][1] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
  A[2][2] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
}


void rotMe( double A[3][3], double r[3] ){
  
  int i,j;
  double rb[3]; // the body frame vector 
 
  for(i=0; i<3; i++ ){
    rb[i] = r[i];
  }
  
  for( i=0; i<3; i++ ){
    r[i] = 0.0;
    
    for( j=0; j<3; j++ ){
      r[i] += A[j][i] * rb[j];
    }
  }
}


/***************************************************************************
 * prints out the usage for the command line arguments, then exits.
 ***************************************************************************/

void usage(){
  (void)fprintf(stderr, 
                "The proper usage is: %s [options] <dumpfile>\n"
                "\n"
                "Options:\n"
                "\n"
                "   -h              Display this message\n"
                "   -o <out_name>   The output file (Defaults to <dumpfile>.xyz)\n"
                "   -d              print the dipole moments\n"
                "   -w              print the waters\n"
                "   -m              map to the periodic box\n"
                "   -n <#>          print every <#> frames\n" 
                "\n"
                "  --repeatX <#>    The number of images to repeat in the x direction\n"
                "  --repeatY <#>    The number of images to repeat in the y direction\n"
                "  --repeatZ <#>    The number of images to repeat in the z direction\n"
                
                "\n",
                
                program_name);
  exit(8);
}

void wrapVector( double thePos[3], double Hmat[3][3], double HmatInv[3][3]){

  int i;
  double scaled[3];

  // calc the scaled coordinates.
  
  matVecMul3(HmatInv, thePos, scaled);
  
  for(i=0; i<3; i++)
    scaled[i] -= roundMe(scaled[i]);
  
  // calc the wrapped real coordinates from the wrapped scaled coordinates
  
  matVecMul3(Hmat, scaled, thePos);
    
}

double matDet3(double a[3][3]) {
  int i, j, k;
  double determinant;

  determinant = 0.0;

  for(i = 0; i < 3; i++) {
    j = (i+1)%3;
    k = (i+2)%3;

    determinant += a[0][i] * (a[1][j]*a[2][k] - a[1][k]*a[2][j]);
  }
  return determinant;
}

void invertMat3(double a[3][3], double b[3][3]) {

  int  i, j, k, l, m, n;
  double determinant;

  determinant = matDet3( a );

  if (determinant == 0.0) {
    printf("Can't invert a matrix with a zero determinant!\n");
  }

  for (i=0; i < 3; i++) {
    j = (i+1)%3;
    k = (i+2)%3;
    for(l = 0; l < 3; l++) {
      m = (l+1)%3;
      n = (l+2)%3;

      b[l][i] = (a[j][m]*a[k][n] - a[j][n]*a[k][m]) / determinant;
    }
  }
}


void matVecMul3(double m[3][3], double inVec[3], double outVec[3]) {
  double a0, a1, a2;

  a0 = inVec[0];  a1 = inVec[1];  a2 = inVec[2];

  outVec[0] = m[0][0]*a0 + m[0][1]*a1 + m[0][2]*a2;
  outVec[1] = m[1][0]*a0 + m[1][1]*a1 + m[1][2]*a2;
  outVec[2] = m[2][0]*a0 + m[2][1]*a1 + m[2][2]*a2;
}
Revision:	1108
Committed:	Wed Apr 14 15:37:41 2004 UTC (20 years, 2 months ago) by tim
Content type:	text/plain
File size:	20554 byte(s)
Log Message:	Change DumpWriter and InitFromFile
#	User	Rev	Content
1	gezelter	444	#include <stdio.h>
2			#include <stdlib.h>
3			#include <string.h>
4			#include <math.h>
5			#include <unistd.h>
6			#include <sys/types.h>
7			#include <sys/stat.h>
8
9	gezelter	624	inline double roundMe( double x ){
10			return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
11			}
12
13	gezelter	444	struct coords{
14	gezelter	624	double pos[3]; // cartesian coords
15	gezelter	444	double q[4]; // the quanternions
16			char name[30];
17			};
18
19			struct linked_xyz{
20			struct coords *r;
21			double time;
22	gezelter	624	double Hmat[3][3];
23			double HmatI[3][3];
24	gezelter	444	struct linked_xyz *next;
25			};
26
27			char program_name; /the name of the program */
28			int printWater = 0;
29			int printDipole = 0;
30
31			void usage(void);
32
33			void rotWrite( FILE* out_file, struct coords* theSSD );
34
35			void setRot( double A[3][3], double q[4] );
36
37			void rotMe( double A[3][3], double r[3] );
38	gezelter	624	void wrapVector( double thePos[3], double Hmat[3][3], double HmatI[3][3]);
39			double matDet3(double a[3][3]);
40			void invertMat3(double a[3][3], double b[3][3]);
41			void matVecMul3(double m[3][3], double inVec[3], double outVec[3]);
42	gezelter	444
43			int main(argc, argv)
44			int argc;
45			char *argv[];
46			{
47
48
49			struct coords *out_coords;
50
51	gezelter	624	int i, j, k, l, m; /* loop counters */
52	gezelter	444	mode_t dir_mode = S_IRWXU;
53
54			int lineCount = 0; //the line number
55			int newN; // the new number of atoms
56			int nSSD=0; // the number of SSD per frame
57
58			unsigned int nAtoms; /the number of atoms in each time step /
59			char read_buffer[1000]; /the line buffer for reading /
60			char eof_test; /ptr to see when we reach the end of the file */
61			char foo; /the pointer to the current string token */
62			FILE in_file; / the input file */
63			FILE out_file; /the output file */
64			char out_prefix = NULL; /the prefix of the output file */
65			char out_name[500]; /the output name /
66			int have_outName =0;
67			char in_name[500]; // the input file name
68			char temp_name[500];
69			char root_name = NULL; /the root name of the input file */
70			unsigned int n_out = 0; /keeps track of which output file is being written/
71
72
73			int skipFrames = 1;
74
75			struct linked_xyz *current_frame;
76			struct linked_xyz *temp_frame;
77
78			int nframes =0;
79
80			int wrap = 0;
81
82			double cX = 0.0;
83			double cY = 0.0;
84			double cZ = 0.0;
85
86			double mcX, mcY, mcZ;
87			double newCX, newCY, newCZ;
88			double dx, dy, dz;
89			int mcount;
90
91			int repeatX = 0;
92			int repeatY = 0;
93			int repeatZ = 0;
94
95			int nRepeatX = 0;
96			int nRepeatY = 0;
97			int nRepeatZ = 0;
98
99			struct coords* rCopy;
100
101			int done;
102			char current_flag;
103
104			program_name = argv[0]; /save the program name in case we need it/
105
106			for( i = 1; i < argc; i++){
107
108			if(argv[i][0] =='-'){
109
110			if(argv[i][1] == '-' ){
111
112			// parse long word options
113
114			if( !strcmp( argv[i], "--repeatX" ) ){
115			repeatX = 1;
116			i++;
117			nRepeatX = atoi( argv[i] );
118			}
119
120			else if( !strcmp( argv[i], "--repeatY" ) ){
121			repeatY = 1;
122			i++;
123			nRepeatY = atoi( argv[i] );
124			}
125
126			else if( !strcmp( argv[i], "--repeatZ" ) ){
127			repeatZ = 1;
128			i++;
129			nRepeatZ = atoi( argv[i] );
130			}
131
132			else{
133			fprintf( stderr,
134			"Invalid option \"%s\"\n", argv[i] );
135			usage();
136			}
137			}
138
139			else{
140
141			// parse single character options
142
143			done =0;
144			j = 1;
145			current_flag = argv[i][j];
146			while( (current_flag != '\0') && (!done) ){
147
148			switch(current_flag){
149
150			case 'o':
151			// -o <prefix> => the output
152
153			i++;
154			strcpy( out_name, argv[i] );
155			have_outName = 1;
156			done = 1;
157			break;
158
159			case 'n':
160			// -n <#> => print every <#> frames
161
162			i++;
163			skipFrames = atoi(argv[i]);
164			done = 1;
165			break;
166
167			case 'h':
168			// -h => give the usage
169
170			usage();
171			break;
172
173			case 'd':
174			// print the dipoles
175
176			printDipole = 1;
177			break;
178
179			case 'w':
180			// print the waters
181
182			printWater = 1;
183			break;
184
185			case 'm':
186			// map to a periodic box
187
188			wrap = 1;
189			break;
190
191			default:
192
193			(void)fprintf(stderr, "Bad option \"-%s\"\n", current_flag);
194			usage();
195			}
196			j++;
197			current_flag = argv[i][j];
198			}
199			}
200			}
201
202			else{
203
204			if( root_name != NULL ){
205			fprintf( stderr,
206			"Error at \"%s\", program does not currently support\n"
207			"more than one input file.\n"
208			"\n",
209			argv[i]);
210			usage();
211			}
212
213			root_name = argv[i];
214			}
215			}
216
217			if(root_name == NULL){
218			usage();
219			}
220
221			sprintf( in_name, "%s", root_name );
222			if( !have_outName ) sprintf( out_name, "%s.xyz", root_name );
223
224			in_file = fopen(in_name, "r");
225			if(in_file == NULL){
226			printf("Cannot open file \"%s\" for reading.\n", in_name);
227			exit(8);
228			}
229
230			out_file = fopen( out_name, "w" );
231			if( out_file == NULL ){
232			printf("Cannot open file \"%s\" for writing.\n", out_name);
233			exit(8);
234			}
235
236			// start reading the first frame
237
238			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
239			lineCount++;
240
241			current_frame = (struct linked_xyz *)malloc(sizeof(struct linked_xyz));
242			current_frame->next = NULL;
243
244
245			while(eof_test != NULL){
246
247			(void)sscanf(read_buffer, "%d", &nAtoms);
248			current_frame->r =
249			(struct coords *)calloc(nAtoms, sizeof(struct coords));
250
251			// read and the comment line and grab the time and box dimensions
252
253			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
254			lineCount++;
255			if(eof_test == NULL){
256			printf("error in reading file at line: %d\n", lineCount);
257			exit(8);
258			}
259
260			foo = strtok( read_buffer, " ,;\t" );
261			sscanf( read_buffer, "%lf", &current_frame->time );
262
263			foo = strtok(NULL, " ,;\t");
264			if(foo == NULL){
265			printf("error in reading file at line: %d\n", lineCount);
266			exit(8);
267			}
268	gezelter	624	(void)sscanf(foo, "%lf",&current_frame->Hmat[0][0]);
269	gezelter	444
270			foo = strtok(NULL, " ,;\t");
271			if(foo == NULL){
272			printf("error in reading file at line: %d\n", lineCount);
273			exit(8);
274			}
275	gezelter	624	(void)sscanf(foo, "%lf",&current_frame->Hmat[1][0]);
276	gezelter	444
277			foo = strtok(NULL, " ,;\t");
278			if(foo == NULL){
279			printf("error in reading file at line: %d\n", lineCount);
280			exit(8);
281			}
282	gezelter	624	(void)sscanf(foo, "%lf",&current_frame->Hmat[2][0]);
283	gezelter	444
284	gezelter	624	foo = strtok(NULL, " ,;\t");
285			if(foo == NULL){
286			printf("error in reading file at line: %d\n", lineCount);
287			exit(8);
288			}
289			(void)sscanf(foo, "%lf",&current_frame->Hmat[0][1]);
290
291			foo = strtok(NULL, " ,;\t");
292			if(foo == NULL){
293			printf("error in reading file at line: %d\n", lineCount);
294			exit(8);
295			}
296			(void)sscanf(foo, "%lf",&current_frame->Hmat[1][1]);
297
298			foo = strtok(NULL, " ,;\t");
299			if(foo == NULL){
300			printf("error in reading file at line: %d\n", lineCount);
301			exit(8);
302			}
303			(void)sscanf(foo, "%lf",&current_frame->Hmat[2][1]);
304
305			foo = strtok(NULL, " ,;\t");
306			if(foo == NULL){
307			printf("error in reading file at line: %d\n", lineCount);
308			exit(8);
309			}
310			(void)sscanf(foo, "%lf",&current_frame->Hmat[0][2]);
311
312			foo = strtok(NULL, " ,;\t");
313			if(foo == NULL){
314			printf("error in reading file at line: %d\n", lineCount);
315			exit(8);
316			}
317			(void)sscanf(foo, "%lf",&current_frame->Hmat[1][2]);
318
319			foo = strtok(NULL, " ,;\t");
320			if(foo == NULL){
321			printf("error in reading file at line: %d\n", lineCount);
322			exit(8);
323			}
324			(void)sscanf(foo, "%lf",&current_frame->Hmat[2][2]);
325
326
327			// Find HmatI:
328
329			invertMat3(current_frame->Hmat, current_frame->HmatI);
330
331
332	gezelter	444	for( i=0; i < nAtoms; i++){
333
334			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
335			lineCount++;
336			if(eof_test == NULL){
337			printf("error in reading file at line: %d\n", lineCount);
338			exit(8);
339			}
340
341			foo = strtok(read_buffer, " ,;\t");
342			if ( !strcmp( "SSD", foo ) ) nSSD++;
343			(void)strcpy(current_frame->r[i].name, foo); //copy the atom name
344
345			// next we grab the positions
346
347			foo = strtok(NULL, " ,;\t");
348			if(foo == NULL){
349			printf("error in reading postition x from %s\n"
350			"natoms = %d, line = %d\n",
351			in_name, nAtoms, lineCount );
352			exit(8);
353			}
354	gezelter	624	(void)sscanf( foo, "%lf", &current_frame->r[i].pos[0] );
355	gezelter	444
356			foo = strtok(NULL, " ,;\t");
357			if(foo == NULL){
358			printf("error in reading postition y from %s\n"
359			"natoms = %d, line = %d\n",
360			in_name, nAtoms, lineCount );
361			exit(8);
362			}
363	gezelter	624	(void)sscanf( foo, "%lf", &current_frame->r[i].pos[1] );
364	gezelter	444
365			foo = strtok(NULL, " ,;\t");
366			if(foo == NULL){
367			printf("error in reading postition z from %s\n"
368			"natoms = %d, line = %d\n",
369			in_name, nAtoms, lineCount );
370			exit(8);
371			}
372	gezelter	624	(void)sscanf( foo, "%lf", &current_frame->r[i].pos[2] );
373	gezelter	444
374			// get the velocities
375
376			foo = strtok(NULL, " ,;\t");
377			if(foo == NULL){
378			printf("error in reading velocity x from %s\n"
379			"natoms = %d, line = %d\n",
380			in_name, nAtoms, lineCount );
381			exit(8);
382			}
383
384
385			foo = strtok(NULL, " ,;\t");
386			if(foo == NULL){
387			printf("error in reading velocity y from %s\n"
388			"natoms = %d, line = %d\n",
389			in_name, nAtoms, lineCount );
390			exit(8);
391			}
392
393
394			foo = strtok(NULL, " ,;\t");
395			if(foo == NULL){
396			printf("error in reading velocity z from %s\n"
397			"natoms = %d, line = %d\n",
398			in_name, nAtoms, lineCount );
399			exit(8);
400			}
401
402			// get the quaternions
403
404			foo = strtok(NULL, " ,;\t");
405			if(foo == NULL){
406			printf("error in reading quaternion 0 from %s\n"
407			"natoms = %d, line = %d\n",
408			in_name, nAtoms, lineCount );
409			exit(8);
410			}
411			(void)sscanf( foo, "%lf", &current_frame->r[i].q[0] );
412
413			foo = strtok(NULL, " ,;\t");
414			if(foo == NULL){
415			printf("error in reading quaternion 1 from %s\n"
416			"natoms = %d, line = %d\n",
417			in_name, nAtoms, lineCount );
418			exit(8);
419			}
420			(void)sscanf( foo, "%lf", &current_frame->r[i].q[1] );
421
422			foo = strtok(NULL, " ,;\t");
423			if(foo == NULL){
424			printf("error in reading quaternion 2 from %s\n"
425			"natoms = %d, line = %d\n",
426			in_name, nAtoms, lineCount );
427			exit(8);
428			}
429			(void)sscanf( foo, "%lf", &current_frame->r[i].q[2] );
430
431			foo = strtok(NULL, " ,;\t");
432			if(foo == NULL){
433			printf("error in reading quaternion 3 from %s\n"
434			"natoms = %d, line = %d\n",
435			in_name, nAtoms, lineCount );
436			exit(8);
437			}
438			(void)sscanf( foo, "%lf", &current_frame->r[i].q[3] );
439
440			// get the angular velocities
441
442			foo = strtok(NULL, " ,;\t");
443			if(foo == NULL){
444			printf("error in reading angular momentum jx from %s\n"
445			"natoms = %d, line = %d\n",
446			in_name, nAtoms, lineCount );
447			exit(8);
448			}
449
450			foo = strtok(NULL, " ,;\t");
451			if(foo == NULL){
452			printf("error in reading angular momentum jy from %s\n"
453			"natoms = %d, line = %d\n",
454			in_name, nAtoms, lineCount );
455			exit(8);
456			}
457
458			foo = strtok(NULL, " ,;\t");
459			if(foo == NULL){
460			printf("error in reading angular momentum jz from %s\n"
461			"natoms = %d, line = %d\n",
462			in_name, nAtoms, lineCount );
463			exit(8);
464			}
465			}
466
467			// if necassary, wrap into the periodic image
468
469			if(wrap){
470
471			for( i=0; i<nAtoms; i++ ){
472
473			/*
474			if( !strcmp( current_frame->r[i].name, "HEAD" ) ){
475
476			// find the center of the lipid
477
478			mcX = 0.0;
479			mcY = 0.0;
480			mcZ = 0.0;
481			mcount = 0;
482
483			mcX += current_frame->r[i].x;
484			mcY += current_frame->r[i].y;
485			mcZ += current_frame->r[i].z;
486			mcount++;
487			i++;
488
489			while( strcmp( current_frame->r[i].name, "HEAD" ) &&
490			strcmp( current_frame->r[i].name, "SSD" ) ){
491
492			mcX += current_frame->r[i].x;
493			mcY += current_frame->r[i].y;
494			mcZ += current_frame->r[i].z;
495			mcount++;
496			i++;
497			}
498
499			mcX /= mcount;
500			mcY /= mcount;
501			mcZ /= mcount;
502
503			newCX = mcX;
504			newCY = mcY;
505			newCZ = mcZ;
506
507			map( &newCX, &newCY, &newCZ,
508			cX, cY, cZ,
509			current_frame->boxX,
510			current_frame->boxY,
511			current_frame->boxZ );
512
513			dx = mcX - newCX;
514			dy = mcY - newCY;
515			dz = mcZ - newCZ;
516
517			for(j=(i-mcount); j<mcount; j++ ){
518
519			current_frame->r[j].x -= dx;
520			current_frame->r[j].y -= dy;
521			current_frame->r[j].z -= dz;
522			}
523
524			i--; // so we don't skip the next atom
525			}
526			*/
527			// else
528	gezelter	624
529			wrapVector(current_frame->r[i].pos,
530			current_frame->Hmat,
531			current_frame->HmatI);
532	gezelter	444
533			}
534	gezelter	624	}
535	gezelter	444
536	tim	1074	nframes++;
537	gezelter	444
538			// write out here
539
540			if(printWater) newN = nAtoms + ( nSSD * 3 );
541			else newN = nAtoms - nSSD;
542
543			newN *= (nRepeatX+1);
544			newN *= (nRepeatY+1);
545			newN *= (nRepeatZ+1);
546
547			if( !(nframes%skipFrames) ){
548			fprintf( out_file,
549			"%d\n"
550	gezelter	624	"%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\n",
551	gezelter	444	newN, current_frame->time,
552	gezelter	624	current_frame->Hmat[0][0], current_frame->Hmat[1][0],
553			current_frame->Hmat[2][0],
554			current_frame->Hmat[0][1], current_frame->Hmat[1][1],
555			current_frame->Hmat[2][1],
556			current_frame->Hmat[0][2], current_frame->Hmat[1][2],
557			current_frame->Hmat[2][2] );
558	gezelter	444
559			rCopy = (struct coords*)
560			calloc( nAtoms, sizeof( struct coords ));
561
562			for(i=0; i<nAtoms; i++){
563			rCopy[i].q[0] = current_frame->r[i].q[0];
564			rCopy[i].q[1] = current_frame->r[i].q[1];
565			rCopy[i].q[2] = current_frame->r[i].q[2];
566			rCopy[i].q[3] = current_frame->r[i].q[3];
567
568			strcpy(rCopy[i].name, current_frame->r[i].name);
569			}
570
571			for(j=0; j<(nRepeatX+1); j++){
572			for(k=0; k<(nRepeatY+1); k++){
573			for(l=0; l<(nRepeatZ+1); l++){
574
575			for(i=0; i<nAtoms; i++){
576	gezelter	624
577			for(m = 0; m < 3; m++) {
578			rCopy[i].pos[m] = current_frame->r[i].pos[m] +
579			j * current_frame->Hmat[m][0] +
580			k * current_frame->Hmat[m][1] +
581			l * current_frame->Hmat[m][2];
582			}
583
584	gezelter	444	if( !strcmp( "SSD", rCopy[i].name ) ){
585
586			rotWrite( out_file, &rCopy[i] );
587			}
588
589			else if( !strcmp( "HEAD", rCopy[i].name ) ){
590
591			rotWrite( out_file, &rCopy[i] );
592			}
593			else{
594	tim	1108	//modify
595	gezelter	444	fprintf( out_file,
596	tim	1108	"%s\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
597	gezelter	444	rCopy[i].name,
598	gezelter	624	rCopy[i].pos[0],
599			rCopy[i].pos[1],
600			rCopy[i].pos[2] );
601	gezelter	444	}
602			}
603			}
604			}
605			}
606	gezelter	624
607	gezelter	444	free( rCopy );
608			}
609	gezelter	624
610	gezelter	444	/free up memory /
611	gezelter	624
612	gezelter	444	temp_frame = current_frame->next;
613			current_frame->next = NULL;
614	gezelter	624
615	gezelter	444	if(temp_frame != NULL){
616	gezelter	624
617	gezelter	444	free(temp_frame->r);
618			free(temp_frame);
619			}
620	gezelter	624
621	gezelter	444	/* make a new frame */
622	gezelter	624
623	gezelter	444	temp_frame = (struct linked_xyz *)malloc(sizeof(struct linked_xyz));
624			temp_frame->next = current_frame;
625			current_frame = temp_frame;
626
627			nSSD = 0;
628
629			eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
630			lineCount++;
631			}
632	gezelter	624
633	gezelter	444	(void)fclose(in_file);
634	gezelter	624
635	gezelter	444	return 0;
636
637			}
638
639
640
641			void rotWrite( FILE* out_file, struct coords* theDipole ){
642
643			double u[3];
644			double h1[3];
645			double h2[3];
646			double ox[3];
647
648			double A[3][3];
649
650	gezelter	624
651	gezelter	444	u[0] = 0.0; u[1] = 0.0; u[2] = 1.0;
652
653			h1[0] = 0.0; h1[1] = -0.75695; h1[2] = 0.5206;
654
655			h2[0] = 0.0; h2[1] = 0.75695; h2[2] = 0.5206;
656
657			ox[0] = 0.0; ox[1] = 0.0; ox[2] = -0.0654;
658
659
660			setRot( A, theDipole->q );
661			rotMe( A, u );
662
663			if( !strcmp( "SSD", theDipole->name )){
664
665			rotMe( A, h1 );
666			rotMe( A, h2 );
667			rotMe( A, ox );
668
669			if( printWater ){
670
671			if(printDipole) {
672
673			fprintf( out_file,
674			"X\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
675	gezelter	624	theDipole->pos[0],
676			theDipole->pos[1],
677			theDipole->pos[2],
678	gezelter	444	u[0], u[1], u[2] );
679	gezelter	624
680			fprintf( out_file,
681			"O\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
682			theDipole->pos[0] + ox[0],
683			theDipole->pos[1] + ox[1],
684			theDipole->pos[2] + ox[2] );
685
686			fprintf( out_file,
687			"H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
688			theDipole->pos[0] + h1[0],
689			theDipole->pos[1] + h1[1],
690			theDipole->pos[2] + h1[2] );
691
692			fprintf( out_file,
693			"H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
694			theDipole->pos[0] + h2[0],
695			theDipole->pos[1] + h2[1],
696			theDipole->pos[2] + h2[2] );
697	gezelter	444	}
698			else{
699
700			fprintf( out_file,
701			"X\t%lf\t%lf\t%lf\n",
702	gezelter	624	theDipole->pos[0],
703			theDipole->pos[1],
704			theDipole->pos[2] );
705
706			fprintf( out_file,
707			"O\t%lf\t%lf\t%lf\n",
708			theDipole->pos[0] + ox[0],
709			theDipole->pos[1] + ox[1],
710			theDipole->pos[2] + ox[2] );
711
712			fprintf( out_file,
713			"H\t%lf\t%lf\t%lf\n",
714			theDipole->pos[0] + h1[0],
715			theDipole->pos[1] + h1[1],
716			theDipole->pos[2] + h1[2] );
717
718			fprintf( out_file,
719			"H\t%lf\t%lf\t%lf\n",
720			theDipole->pos[0] + h2[0],
721			theDipole->pos[1] + h2[1],
722			theDipole->pos[2] + h2[2] );
723			}
724	gezelter	444	}
725			}
726
727			else if( !strcmp( "HEAD", theDipole->name )) {
728
729			if( printDipole ){
730
731			fprintf( out_file,
732			"HEAD\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
733	gezelter	624	theDipole->pos[0],
734			theDipole->pos[1],
735			theDipole->pos[2],
736	gezelter	444	u[0], u[1], u[2] );
737			}
738			else{
739
740			fprintf( out_file,
741			"HEAD\t%lf\t%lf\t%lf\n",
742	gezelter	624	theDipole->pos[0],
743			theDipole->pos[1],
744			theDipole->pos[2] );
745	gezelter	444	}
746			}
747			}
748
749
750			void setRot( double A[3][3], double the_q[4] ){
751
752			double q0Sqr, q1Sqr, q2Sqr, q3Sqr;
753
754			q0Sqr = the_q[0] * the_q[0];
755			q1Sqr = the_q[1] * the_q[1];
756			q2Sqr = the_q[2] * the_q[2];
757			q3Sqr = the_q[3] * the_q[3];
758
759			A[0][0] = q0Sqr + q1Sqr - q2Sqr - q3Sqr;
760			A[0][1] = 2.0 * ( the_q[1] * the_q[2] + the_q[0] * the_q[3] );
761			A[0][2] = 2.0 * ( the_q[1] * the_q[3] - the_q[0] * the_q[2] );
762
763			A[1][0] = 2.0 * ( the_q[1] * the_q[2] - the_q[0] * the_q[3] );
764			A[1][1] = q0Sqr - q1Sqr + q2Sqr - q3Sqr;
765			A[1][2] = 2.0 * ( the_q[2] * the_q[3] + the_q[0] * the_q[1] );
766
767			A[2][0] = 2.0 * ( the_q[1] * the_q[3] + the_q[0] * the_q[2] );
768			A[2][1] = 2.0 * ( the_q[2] * the_q[3] - the_q[0] * the_q[1] );
769			A[2][2] = q0Sqr - q1Sqr -q2Sqr +q3Sqr;
770			}
771
772
773
774			void rotMe( double A[3][3], double r[3] ){
775
776			int i,j;
777			double rb[3]; // the body frame vector
778
779			for(i=0; i<3; i++ ){
780			rb[i] = r[i];
781			}
782
783			for( i=0; i<3; i++ ){
784			r[i] = 0.0;
785
786			for( j=0; j<3; j++ ){
787			r[i] += A[j][i] * rb[j];
788			}
789			}
790			}
791
792
793
794			/***************************************************************************
795			* prints out the usage for the command line arguments, then exits.
796			***************************************************************************/
797
798			void usage(){
799			(void)fprintf(stderr,
800			"The proper usage is: %s [options] <dumpfile>\n"
801			"\n"
802			"Options:\n"
803			"\n"
804			" -h Display this message\n"
805			" -o <out_name> The output file (Defaults to <dumpfile>.xyz)\n"
806			" -d print the dipole moments\n"
807			" -w print the waters\n"
808			" -m map to the periodic box\n"
809			" -n <#> print every <#> frames\n"
810			"\n"
811			" --repeatX <#> The number of images to repeat in the x direction\n"
812			" --repeatY <#> The number of images to repeat in the y direction\n"
813			" --repeatZ <#> The number of images to repeat in the z direction\n"
814
815			"\n",
816
817			program_name);
818			exit(8);
819			}
820
821	gezelter	624	void wrapVector( double thePos[3], double Hmat[3][3], double HmatInv[3][3]){
822	gezelter	444
823	gezelter	624	int i;
824			double scaled[3];
825	gezelter	444
826	gezelter	624	// calc the scaled coordinates.
827	gezelter	444
828	gezelter	624	matVecMul3(HmatInv, thePos, scaled);
829
830			for(i=0; i<3; i++)
831			scaled[i] -= roundMe(scaled[i]);
832
833			// calc the wrapped real coordinates from the wrapped scaled coordinates
834
835			matVecMul3(Hmat, scaled, thePos);
836
837	gezelter	444	}
838	gezelter	624
839			double matDet3(double a[3][3]) {
840			int i, j, k;
841			double determinant;
842
843			determinant = 0.0;
844
845			for(i = 0; i < 3; i++) {
846			j = (i+1)%3;
847			k = (i+2)%3;
848
849			determinant += a[0][i] * (a[1][j]a[2][k] - a[1][k]a[2][j]);
850			}
851			return determinant;
852			}
853
854			void invertMat3(double a[3][3], double b[3][3]) {
855
856			int i, j, k, l, m, n;
857			double determinant;
858
859			determinant = matDet3( a );
860
861			if (determinant == 0.0) {
862			printf("Can't invert a matrix with a zero determinant!\n");
863			}
864
865			for (i=0; i < 3; i++) {
866			j = (i+1)%3;
867			k = (i+2)%3;
868			for(l = 0; l < 3; l++) {
869			m = (l+1)%3;
870			n = (l+2)%3;
871
872			b[l][i] = (a[j][m]a[k][n] - a[j][n]a[k][m]) / determinant;
873			}
874			}
875			}
876
877
878			void matVecMul3(double m[3][3], double inVec[3], double outVec[3]) {
879			double a0, a1, a2;
880
881			a0 = inVec[0]; a1 = inVec[1]; a2 = inVec[2];
882
883			outVec[0] = m[0][0]a0 + m[0][1]a1 + m[0][2]*a2;
884			outVec[1] = m[1][0]a0 + m[1][1]a1 + m[1][2]*a2;
885			outVec[2] = m[2][0]a0 + m[2][1]a1 + m[2][2]*a2;
886			}