ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/utils/quickLate.c
(Generate patch)

Comparing trunk/OOPSE/utils/quickLate.c (file contents):
Revision 623 by gezelter, Thu Apr 3 13:43:02 2003 UTC vs.
Revision 624 by gezelter, Wed Jul 16 17:35:18 2003 UTC

# Line 6 | Line 6
6   #include <sys/types.h>
7   #include <sys/stat.h>
8  
9 + inline double roundMe( double x ){
10 +  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
11 + }
12 +
13   struct coords{
14 <  double x,y,z; // cartesian coords
14 >  double pos[3]; // cartesian coords
15    double q[4];  // the quanternions
16    char name[30];
17   };
# Line 15 | Line 19 | struct linked_xyz{
19   struct linked_xyz{
20    struct coords *r;
21    double time;
22 <  double boxX, boxY, boxZ;
22 >  double Hmat[3][3];
23 >  double HmatI[3][3];
24    struct linked_xyz *next;
25   };
26  
# Line 30 | Line 35 | void rotMe( double A[3][3], double r[3] );
35   void setRot( double A[3][3], double q[4] );
36  
37   void rotMe( double A[3][3], double r[3] );
38 <
39 < void map( double *x, double *y, double *z, double centerX, double centerY,
40 <          double centerZ, double boxX, double boxY, double boxZ );
38 > 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  
43   int main(argc, argv)
44       int argc;
# Line 42 | Line 48 | int main(argc, argv)
48  
49    struct coords *out_coords;
50  
51 <  int i, j, k, l; /* loop counters */
51 >  int i, j, k, l, m; /* loop counters */
52    mode_t dir_mode = S_IRWXU;
53  
54    int lineCount = 0;  //the line number
# Line 259 | Line 265 | int main(argc, argv)
265        printf("error in reading file at line: %d\n", lineCount);
266        exit(8);
267      }
268 <    (void)sscanf(foo, "%lf",&current_frame->boxX);
268 >    (void)sscanf(foo, "%lf",&current_frame->Hmat[0][0]);
269  
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 <    (void)sscanf(foo, "%lf",&current_frame->boxY);
275 >    (void)sscanf(foo, "%lf",&current_frame->Hmat[1][0]);
276  
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 <    (void)sscanf(foo, "%lf",&current_frame->boxZ);
282 >    (void)sscanf(foo, "%lf",&current_frame->Hmat[2][0]);
283 >
284 >    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      for( i=0; i < nAtoms; i++){
333        
334        eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
# Line 297 | Line 351 | int main(argc, argv)
351                 in_name, nAtoms, lineCount );
352          exit(8);
353        }
354 <      (void)sscanf( foo, "%lf", &current_frame->r[i].x );
354 >      (void)sscanf( foo, "%lf", &current_frame->r[i].pos[0] );
355        
356        foo = strtok(NULL, " ,;\t");
357        if(foo == NULL){
# Line 306 | Line 360 | int main(argc, argv)
360                 in_name, nAtoms, lineCount );
361          exit(8);
362        }
363 <      (void)sscanf( foo, "%lf", &current_frame->r[i].y );
363 >      (void)sscanf( foo, "%lf", &current_frame->r[i].pos[1] );
364        
365        foo = strtok(NULL, " ,;\t");
366        if(foo == NULL){
# Line 315 | Line 369 | int main(argc, argv)
369                 in_name, nAtoms, lineCount );
370          exit(8);
371        }
372 <      (void)sscanf( foo, "%lf", &current_frame->r[i].z );
372 >      (void)sscanf( foo, "%lf", &current_frame->r[i].pos[2] );
373      
374        // get the velocities
375        
# Line 471 | Line 525 | int main(argc, argv)
525          }
526          */
527          // else
528 <        map( &current_frame->r[i].x,
529 <                  &current_frame->r[i].y,
530 <                  &current_frame->r[i].z,
531 <                  cX, cY, cZ,
478 <                  current_frame->boxX,
479 <                  current_frame->boxY,
480 <                  current_frame->boxZ );
528 >        
529 >        wrapVector(current_frame->r[i].pos,
530 >                   current_frame->Hmat,
531 >                   current_frame->HmatI);
532          
533        }
534 <    }
484 <    
534 >    }  
535  
486 –
487 –
536      nframes++;
537      
538      // write out here    
# Line 499 | Line 547 | int main(argc, argv)
547      if( !(nframes%skipFrames) ){
548        fprintf( out_file,        
549                 "%d\n"
550 <               "%lf\t%lf\t%lf\t%lf\n",
550 >               "%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\t%lf\t%lf\t%lf;\n",
551                 newN, current_frame->time,
552 <               current_frame->boxX, current_frame->boxY,
553 <               current_frame->boxZ );
552 >               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      
559        rCopy = (struct coords*)
560          calloc( nAtoms, sizeof( struct coords ));
# Line 521 | Line 573 | int main(argc, argv)
573            for(l=0; l<(nRepeatZ+1); l++){
574  
575              for(i=0; i<nAtoms; i++){
576 <              rCopy[i].x = current_frame->r[i].x + (j * current_frame->boxX);
577 <              rCopy[i].y = current_frame->r[i].y + (k * current_frame->boxY);
578 <              rCopy[i].z = current_frame->r[i].z + (l * current_frame->boxZ);
579 <            
576 >
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                if( !strcmp( "SSD", rCopy[i].name ) ){
585                  
586                  rotWrite( out_file, &rCopy[i] );
# Line 539 | Line 595 | int main(argc, argv)
595                  fprintf( out_file,
596                           "%s\t%lf\t%lf\t%lf\n",
597                           rCopy[i].name,
598 <                         rCopy[i].x,
599 <                         rCopy[i].y,
600 <                         rCopy[i].z );
598 >                         rCopy[i].pos[0],
599 >                         rCopy[i].pos[1],
600 >                         rCopy[i].pos[2] );
601                }
602              }
603            }
604          }
605        }
606 <
606 >      
607        free( rCopy );
608      }
609 <          
609 >    
610      /*free up memory */
611 <
611 >    
612      temp_frame = current_frame->next;
613      current_frame->next = NULL;
614 <
614 >    
615      if(temp_frame != NULL){
616 <
616 >      
617        free(temp_frame->r);
618        free(temp_frame);
619      }
620 <
620 >    
621      /* make a new frame */
622 <
622 >    
623      temp_frame = (struct linked_xyz *)malloc(sizeof(struct linked_xyz));
624      temp_frame->next = current_frame;
625      current_frame = temp_frame;
# Line 573 | Line 629 | int main(argc, argv)
629      eof_test = fgets(read_buffer, sizeof(read_buffer), in_file);
630      lineCount++;
631    }
632 <
632 >  
633    (void)fclose(in_file);
634 <
634 >  
635    return 0;
636    
637   }
# Line 591 | Line 647 | void rotWrite( FILE* out_file, struct coords* theDipol
647    
648    double A[3][3];
649    
650 <  
595 <
650 >
651    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;
# Line 617 | Line 672 | void rotWrite( FILE* out_file, struct coords* theDipol
672          
673          fprintf( out_file,
674                   "X\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
675 <                 theDipole->x,
676 <                 theDipole->y,
677 <                 theDipole->z,
675 >                 theDipole->pos[0],
676 >                 theDipole->pos[1],
677 >                 theDipole->pos[2],
678                   u[0], u[1], u[2] );
679 +
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          }
698        else{
699          
700          fprintf( out_file,
701                   "X\t%lf\t%lf\t%lf\n",
702 <                 theDipole->x,
703 <                 theDipole->y,
704 <                 theDipole->z );
705 <      }
706 <      
707 <      fprintf( out_file,
708 <               "O\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
709 <               theDipole->x + ox[0],
710 <               theDipole->y + ox[1],
711 <               theDipole->z + ox[2] );
712 <      
713 <      fprintf( out_file,
714 <               "H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
715 <               theDipole->x + h1[0],
716 <               theDipole->y + h1[1],
717 <               theDipole->z + h1[2] );
718 <      
719 <      fprintf( out_file,
720 <               "H\t%lf\t%lf\t%lf\t0.0\t0.0\t0.0\n",
721 <               theDipole->x + h2[0],
722 <               theDipole->y + h2[1],
723 <               theDipole->z + h2[2] );
702 >                 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      }
725    }
726    
# Line 657 | Line 730 | void rotWrite( FILE* out_file, struct coords* theDipol
730        
731        fprintf( out_file,
732                 "HEAD\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n",
733 <               theDipole->x,
734 <               theDipole->y,
735 <               theDipole->z,
733 >               theDipole->pos[0],
734 >               theDipole->pos[1],
735 >               theDipole->pos[2],
736                 u[0], u[1], u[2] );
737      }
738      else{
739      
740        fprintf( out_file,
741                 "HEAD\t%lf\t%lf\t%lf\n",
742 <               theDipole->x,
743 <               theDipole->y,
744 <               theDipole->z );
742 >               theDipole->pos[0],
743 >               theDipole->pos[1],
744 >               theDipole->pos[2] );
745      }
746    }
747   }
# Line 745 | Line 818 | void usage(){
818    exit(8);
819   }
820  
821 < 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 < {
821 > void wrapVector( double thePos[3], double Hmat[3][3], double HmatInv[3][3]){
822  
823 <  *x -= centerX;
824 <  *y -= centerY;
756 <  *z -= centerZ;
823 >  int i;
824 >  double scaled[3];
825  
826 <  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));
826 >  // calc the scaled coordinates.
827    
828 <  if(*z < 0) *z -= boxZ * (double)( (int)( (*z / boxZ)  - 0.5 ) );
829 <  else *z -= boxZ * (double)( (int)( (*z / boxZ ) + 0.5));
828 >  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   }
838 +
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 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines