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>

struct coords{
  double x,y,z; // cartesian coords
  double q[4];  // the quanternions
  char name[30];
};

struct linked_xyz{
  struct coords *r;
  double time;
  double boxX, boxY, boxZ;
  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 map( double *x, double *y, double *z, double centerX, double centerY, 
          double centerZ, double boxX, double boxY, double boxZ );

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


  struct coords *out_coords;

  int i, j, k, l; /* 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->boxX);

    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->boxY);

    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->boxZ);

    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].x );
      
      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].y );
      
      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].z );
    
      // 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 
        map( &current_frame->r[i].x,
                  &current_frame->r[i].y,
                  &current_frame->r[i].z,
                  cX, cY, cZ,
                  current_frame->boxX,
                  current_frame->boxY,
                  current_frame->boxZ );
        
      }
    }
    

    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\n",
               newN, current_frame->time,
               current_frame->boxX, current_frame->boxY,
               current_frame->boxZ );
    
      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++){
              rCopy[i].x = current_frame->r[i].x + (j * current_frame->boxX);
              rCopy[i].y = current_frame->r[i].y + (k * current_frame->boxY);
              rCopy[i].z = current_frame->r[i].z + (l * current_frame->boxZ);
            
              if( !strcmp( "SSD", rCopy[i].name ) ){
                
                rotWrite( out_file, &rCopy[i] );
              }
              
              else if( !strcmp( "HEAD", rCopy[i].name ) ){
                
                rotWrite( out_file, &rCopy[i] );
              }
              else{
                
                fprintf( out_file,
                         "%s\t%lf\t%lf\t%lf\n",
                         rCopy[i].name,
                         rCopy[i].x,
                         rCopy[i].y,
                         rCopy[i].z );
              }
            }
          }
        }
      }

      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->x,
                 theDipole->y,
                 theDipole->z,
                 u[0], u[1], u[2] );
        }
      else{
        
        fprintf( out_file,
                 "X\t%lf\t%lf\t%lf\n",
                 theDipole->x,
                 theDipole->y,
                 theDipole->z );
      }
      
      fprintf( out_file,
               "O\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
               theDipole->x + ox[0],
               theDipole->y + ox[1],
               theDipole->z + ox[2] );
      
      fprintf( out_file,
               "H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
               theDipole->x + h1[0],
               theDipole->y + h1[1],
               theDipole->z + h1[2] );
      
      fprintf( out_file,
               "H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
               theDipole->x + h2[0],
               theDipole->y + h2[1],
               theDipole->z + 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->x,
               theDipole->y,
               theDipole->z,
               u[0], u[1], u[2] );
    }
    else{
     
      fprintf( out_file,
               "HEAD\t%lf\t%lf\t%lf\n",
               theDipole->x,
               theDipole->y,
               theDipole->z );
    }
  }
}


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 map( x, y, z, centerX, centerY, centerZ, boxX, boxY, boxZ )
     double *x, *y, *z;
     double centerX, centerY, centerZ;
     double boxX, boxY, boxZ;
{ 

  *x -= centerX;
  *y -= centerY;
  *z -= centerZ;

  if(*x < 0) *x -= boxX * (double)( (int)( (*x / boxX)  - 0.5 ) );
  else *x -= boxX * (double)( (int)( (*x / boxX ) + 0.5));

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