--- trunk/OOPSE/utils/quickLate.c 2003/04/03 13:43:02 444 +++ trunk/OOPSE/utils/quickLate.c 2003/07/16 17:35:18 624 @@ -6,8 +6,12 @@ struct coords{ #include #include +inline double roundMe( double x ){ + return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); +} + struct coords{ - double x,y,z; // cartesian coords + double pos[3]; // cartesian coords double q[4]; // the quanternions char name[30]; }; @@ -15,7 +19,8 @@ struct linked_xyz{ struct linked_xyz{ struct coords *r; double time; - double boxX, boxY, boxZ; + double Hmat[3][3]; + double HmatI[3][3]; struct linked_xyz *next; }; @@ -30,9 +35,10 @@ void rotMe( double A[3][3], double r[3] ); 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 ); +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; @@ -42,7 +48,7 @@ int main(argc, argv) struct coords *out_coords; - int i, j, k, l; /* loop counters */ + int i, j, k, l, m; /* loop counters */ mode_t dir_mode = S_IRWXU; int lineCount = 0; //the line number @@ -259,22 +265,70 @@ int main(argc, argv) printf("error in reading file at line: %d\n", lineCount); exit(8); } - (void)sscanf(foo, "%lf",¤t_frame->boxX); + (void)sscanf(foo, "%lf",¤t_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",¤t_frame->boxY); + (void)sscanf(foo, "%lf",¤t_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",¤t_frame->boxZ); + (void)sscanf(foo, "%lf",¤t_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",¤t_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",¤t_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",¤t_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",¤t_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",¤t_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",¤t_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); @@ -297,7 +351,7 @@ int main(argc, argv) in_name, nAtoms, lineCount ); exit(8); } - (void)sscanf( foo, "%lf", ¤t_frame->r[i].x ); + (void)sscanf( foo, "%lf", ¤t_frame->r[i].pos[0] ); foo = strtok(NULL, " ,;\t"); if(foo == NULL){ @@ -306,7 +360,7 @@ int main(argc, argv) in_name, nAtoms, lineCount ); exit(8); } - (void)sscanf( foo, "%lf", ¤t_frame->r[i].y ); + (void)sscanf( foo, "%lf", ¤t_frame->r[i].pos[1] ); foo = strtok(NULL, " ,;\t"); if(foo == NULL){ @@ -315,7 +369,7 @@ int main(argc, argv) in_name, nAtoms, lineCount ); exit(8); } - (void)sscanf( foo, "%lf", ¤t_frame->r[i].z ); + (void)sscanf( foo, "%lf", ¤t_frame->r[i].pos[2] ); // get the velocities @@ -471,20 +525,14 @@ int main(argc, argv) } */ // else - map( ¤t_frame->r[i].x, - ¤t_frame->r[i].y, - ¤t_frame->r[i].z, - cX, cY, cZ, - current_frame->boxX, - current_frame->boxY, - current_frame->boxZ ); + + wrapVector(current_frame->r[i].pos, + current_frame->Hmat, + current_frame->HmatI); } - } - + } - - nframes++; // write out here @@ -499,10 +547,14 @@ int main(argc, argv) if( !(nframes%skipFrames) ){ fprintf( out_file, "%d\n" - "%lf\t%lf\t%lf\t%lf\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->boxX, current_frame->boxY, - current_frame->boxZ ); + 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 )); @@ -521,10 +573,14 @@ int main(argc, argv) for(l=0; l<(nRepeatZ+1); l++){ for(i=0; ir[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); - + + 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] ); @@ -539,31 +595,31 @@ int main(argc, argv) fprintf( out_file, "%s\t%lf\t%lf\t%lf\n", rCopy[i].name, - rCopy[i].x, - rCopy[i].y, - rCopy[i].z ); + 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; @@ -573,9 +629,9 @@ int main(argc, argv) eof_test = fgets(read_buffer, sizeof(read_buffer), in_file); lineCount++; } - + (void)fclose(in_file); - + return 0; } @@ -591,8 +647,7 @@ void rotWrite( FILE* out_file, struct coords* theDipol 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; @@ -617,37 +672,55 @@ void rotWrite( FILE* out_file, struct coords* theDipol fprintf( out_file, "X\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", - theDipole->x, - theDipole->y, - theDipole->z, + 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->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] ); + 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] ); + } } } @@ -657,18 +730,18 @@ void rotWrite( FILE* out_file, struct coords* theDipol fprintf( out_file, "HEAD\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", - theDipole->x, - theDipole->y, - theDipole->z, + 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->x, - theDipole->y, - theDipole->z ); + theDipole->pos[0], + theDipole->pos[1], + theDipole->pos[2] ); } } } @@ -745,22 +818,69 @@ void map( x, y, z, centerX, centerY, centerZ, boxX, bo 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; -{ +void wrapVector( double thePos[3], double Hmat[3][3], double HmatInv[3][3]){ - *x -= centerX; - *y -= centerY; - *z -= centerZ; + int i; + double scaled[3]; - 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)); + // calc the scaled coordinates. - if(*z < 0) *z -= boxZ * (double)( (int)( (*z / boxZ) - 0.5 ) ); - else *z -= boxZ * (double)( (int)( (*z / boxZ ) + 0.5)); + 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; +}