[ViewVC] Diff of: group/trunk/madProps/GofR.c

Comparing trunk/madProps/GofR.c (file contents):
Revision 44 by mmeineke, Fri Jul 19 19:05:59 2002 UTC vs.
Revision 45 by mmeineke, Tue Jul 23 16:08:35 2002 UTC

+#include <string.h>
+#include "GofR.h"
-+
-+
+void map( double *x, double *y, double *z,
-+
+          double boxX, double boxY, double boxZ );
-+
-+
+// this algoritm assumes constant number of atoms between frames, and
-+
+// constant boxSize
-+
-+
+void GofR( char* out_prefix, char* atom1, char* atom2,
-+
+           struct xyz_frame* frames, int nFrames ){
-+
-+
+  int i,j,k,l;
-+
-+
+  enum g_types { all_all, atom_all, atom_atom };
-+
+  enum g_types the_type;
-+
+  char* allAtom;
-+
-+
+  double atom1Dens;
-+
+  double atom2Dens;
-+
+  double allDens;
-+
-+
+  double atom1Constant;
-+
+  double atom2Constant;
-+
+  double allConstant;
-+
-+
+  double nAtom1;
-+
+  double nAtom2;
-+
+  double nAtoms;
-+
-+
+  double delR;
-+
+  double boxVol;
-+
+  double shortBox;
-+
-+
+  double dx, dy, dz;
-+
+  double distSqr;
-+
+  double dist;
-+
-+
+  double rLower, rUpper;
-+
+  double nIdeal;
-+
+  double volSlice;
-+
-+
+  int bin;
-+
+  int histogram[GofRBins];
-+
+  double the_GofR[GofRBins];
-+
+  double rValue[GofRBins];
-+
-+
+  char out_name[500];
-+
+  char tempString[100];
-+
+  FILE *out_file;
-+
-+
-+
+   // figure out the type of GofR we are calculating
-+
-+
+  if( !strcmp( atom1, "all" ) ){
-+
-+
+    if( !strcmp( atom2, "all" ) ) the_type = all_all;
-+
+    else{
-+
+      the_type = atom_all;
-+
+      allAtom = atom2;
-+
+    }
-+
+  }
-+
+  else{
-+
-+
+    if( !strcmp( atom2, "all" ) ){
-+
+      the_type = atom_all;
-+
+      allAtom = atom1;
-+
+    }
-+
+    else the_type = atom_atom;
-+
+  }
-+
-+
+  // find the box size and delR;
-+
-+
+  shortBox = frames[0].boxX;
-+
+  if( shortBox > frames[0].boxY ) shortBox = frames[0].boxY;
-+
+  if( shortBox > frames[0].boxZ ) shortBox = frames[0].boxZ;
-+
-+
+  delR = ( shortBox / 2.0 ) / GofRBins;
-+
+  boxVol = frames[0].boxX * frames[0].boxY * frames[0].boxZ;
-+
-+
+  // zero the histograms;
-+
-+
+  for(i=0; i<GofRBins; i++ ){
-+
-+
+    rValue[i]    = 0.0;
-+
+    the_GofR[i]  = 0.0;
-+
+    histogram[i] = 0;
-+
+  }
-+
-+
+  switch( the_type ){
-+
-+
+  case atom_atom:
-+
-+
+    // find the number of each type;
-+
-+
+    nAtom1 = 0;
-+
+    nAtom2 = 0;
-+
+    for( i=0; i<frames[0].nAtoms; i++ ){
-+
-+
+      if( !strcmp( frames[0].names[i], atom1 ) ) nAtom1++;
-+
+      else if( !strcmp( frames[0].names[i], atom2 ) ) nAtom2++;
-+
+    }
-+
-+
+    if( !nAtom1 ){
-+
-+
+      fprintf( stderr,
-+
+               "\n"
-+
+               "GofR error, \"%s\" was not found in the trajectory.\n",
-+
+               atom1 );
-+
+      exit(8);
-+
+    }
-+
-+
+    if( !nAtom2 ){
-+
-+
+      fprintf( stderr,
-+
+               "\n"
-+
+               "GofR error, \"%s\" was not found in the trajectory.\n",
-+
+               atom2 );
-+
+      exit(8);
-+
+    }
-+
-+
+    // calculate some of the constants;
-+
-+
+    atom1Dens = nAtom1 / boxVol;
-+
+    atom2Dens = nAtom2 / boxVol;
-+
-+
+    atom1Constant = ( 4.0 * M_PI * atom1Dens ) / 3.0;
-+
+    atom2Constant = ( 4.0 * M_PI * atom2Dens ) / 3.0;
-+
-+
+    // calculate the histogram
-+
-+
+    for( i=0; i<nFrames; i++){
-+
+      for( j=0; j<(frames[i].nAtoms-1); j++ ){
-+
+        for( k=j+1; k< frames[i].nAtoms; k++ ){
-+
-+
+          if( !strcmp( frames[0].names[j], atom1 ) ){
-+
+            if( !strcmp( frames[0].names[k], atom2 ) ){
-+
-+
+              dx = frames[i].r[j].x - frames[i].r[k].x;
-+
+              dy = frames[i].r[j].y - frames[i].r[k].y;
-+
+              dz = frames[i].r[j].z - frames[i].r[k].z;
-+
-+
+              map( &dx, &dy, &dz,
-+
+                   frames[i].boxX, frames[i].boxY, frames[i].boxZ );
-+
-+
+              distSqr = (dx * dx) + (dy * dy) + (dz * dz);
-+
+              dist = sqrt( distSqr );
-+
-+
+              // add to the appropriate bin
-+
+              bin = (int)( dist / delR );
-+
+              if( bin < GofRBins ) histogram[bin] += 2;
-+
+            }
-+
+          }
-+
-+
+          else if( !strcmp( frames[0].names[j], atom2 ) ){
-+
+            if( !strcmp( frames[0].names[k], atom1 ) ){
-+
-+
+              dx = frames[i].r[j].x - frames[i].r[k].x;
-+
+              dy = frames[i].r[j].y - frames[i].r[k].y;
-+
+              dz = frames[i].r[j].z - frames[i].r[k].z;
-+
-+
+              map( &dx, &dy, &dz,
-+
+                   frames[i].boxX, frames[i].boxY, frames[i].boxZ );
-+
-+
+              distSqr = (dx * dx) + (dy * dy) + (dz * dz);
-+
+              dist = sqrt( distSqr );
-+
-+
+              // add to the appropriate bin
-+
+              bin = (int)( dist / delR );
-+
+              if( bin < GofRBins ) histogram[bin] += 2;
-+
+            }
-+
+          }
-+
+        }
-+
+      }
-+
+    }
-+
-+
+    // calculate the GofR
-+
-+
+    for( i=0; i<GofRBins; i++ ){
-+
-+
+      rLower = i * delR;
-+
+      rUpper = rLower + delR;
-+
-+
+      volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
-+
+      nIdeal = volSlice * ( atom1Constant + atom2Constant );
-+
-+
+      the_GofR[i] = histogram[i] / ( nFrames * ( nAtom1 + nAtom2 ) * nIdeal );
-+
+      rValue[i] = rLower + ( delR / 2.0 );
-+
+    }
-+
-+
+    // make the out_name;
-+
-+
+    strcpy( out_name, out_prefix );
-+
+    sprintf( tempString, "-%s-%s.GofR", atom1, atom2 );
-+
+    strcat( out_name, tempString );
-+
-+
+    out_file = fopen( out_name, "w" );
-+
+    if( out_file == NULL ){
-+
+      fprintf( stderr,
-+
+               "\n"
-+
+               "GofR error, unable to open \"%s\" for writing.\n",
-+
+               out_name );
-+
+      exit(8);
-+
+    }
-+
+    break;
-+
-+
+  case atom_all:
-+
-+
+    // find the number of AllAtoms;
-+
-+
+    nAtom1 = 0;
-+
+    for( i=0; i<frames[0].nAtoms; i++ ){
-+
-+
+      if( !strcmp( frames[0].names[i], allAtom ) ) nAtom1++;
-+
+    }
-+
-+
+    if( !nAtom1 ){
-+
-+
+      fprintf( stderr,
-+
+               "\n"
-+
+               "GofR error, \"%s\" was not found in the trajectory.\n",
-+
+               atom1 );
-+
+      exit(8);
-+
+    }
-+
-+
+    // calculate some of the constants;
-+
-+
+    atom1Dens = nAtom1 / boxVol;
-+
+    allDens  = frames[0].nAtoms / boxVol;
-+
-+
+    atom1Constant = ( 4.0 * M_PI * atom1Dens ) / 3.0;
-+
+    allConstant  = ( 4.0 * M_PI * allDens ) / 3.0;
-+
-+
+    // calculate the histogram
-+
-+
+    for( i=0; i<nFrames; i++){
-+
+      for( j=0; j<(frames[i].nAtoms-1); j++ ){
-+
+        for( k=j+1; k< frames[i].nAtoms; k++ ){
-+
-+
+          if( !strcmp( frames[0].names[j], allAtom ) ){
-+
-+
+            dx = frames[i].r[j].x - frames[i].r[k].x;
-+
+            dy = frames[i].r[j].y - frames[i].r[k].y;
-+
+            dz = frames[i].r[j].z - frames[i].r[k].z;
-+
-+
+            map( &dx, &dy, &dz,
-+
+                 frames[i].boxX, frames[i].boxY, frames[i].boxZ );
-+
-+
+            distSqr = (dx * dx) + (dy * dy) + (dz * dz);
-+
+            dist = sqrt( distSqr );
-+
-+
+            // add to the appropriate bin
-+
+            bin = (int)( dist / delR );
-+
+            if( bin < GofRBins ) histogram[bin] += 2;
-+
+          }
-+
-+
+          else if( !strcmp( frames[0].names[k], allAtom ) ){
-+
-+
+            dx = frames[i].r[j].x - frames[i].r[k].x;
-+
+            dy = frames[i].r[j].y - frames[i].r[k].y;
-+
+            dz = frames[i].r[j].z - frames[i].r[k].z;
-+
-+
+            map( &dx, &dy, &dz,
-+
+                 frames[i].boxX, frames[i].boxY, frames[i].boxZ );
-+
-+
+            distSqr = (dx * dx) + (dy * dy) + (dz * dz);
-+
+            dist = sqrt( distSqr );
-+
-+
+            // add to the appropriate bin
-+
+            bin = (int)( dist / delR );
-+
+            if( bin < GofRBins ) histogram[bin] += 2;
-+
+          }
-+
+        }
-+
+      }
-+
+    }
-+
-+
+    // calculate the GofR
-+
-+
+    for( i=0; i<GofRBins; i++ ){
-+
-+
+      rLower = i * delR;
-+
+      rUpper = rLower + delR;
-+
-+
+      volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
-+
+      nIdeal = volSlice * ( allConstant );
-+
-+
+      the_GofR[i] = histogram[i] / ( nFrames * frames[0].nAtoms * nIdeal );
-+
+      rValue[i] = rLower + ( delR / 2.0 );
-+
+    }
-+
-+
+    // make the out_name;
-+
-+
+    strcpy( out_name, out_prefix );
-+
+    sprintf( tempString, "-%s-%s.GofR", allAtom, "all" );
-+
+    strcat( out_name, tempString );
-+
-+
+    out_file = fopen( out_name, "w" );
-+
+    if( out_file == NULL ){
-+
+      fprintf( stderr,
-+
+               "\n"
-+
+               "GofR error, unable to open \"%s\" for writing.\n",
-+
+               out_name );
-+
+      exit(8);
-+
+    }
-+
+    break;
-+
-+
+  case all_all:
-+
-+
+    // calculate some of the constants;
-+
-+
+    allDens  = frames[0].nAtoms / boxVol;
-+
-+
+    allConstant  = ( 4.0 * M_PI * allDens ) / 3.0;
-+
-+
+    // calculate the histogram
-+
-+
+    for( i=0; i<nFrames; i++){
-+
+      for( j=0; j<(frames[i].nAtoms-1); j++ ){
-+
+        for( k=j+1; k< frames[i].nAtoms; k++ ){
-+
-+
+          dx = frames[i].r[j].x - frames[i].r[k].x;
-+
+          dy = frames[i].r[j].y - frames[i].r[k].y;
-+
+          dz = frames[i].r[j].z - frames[i].r[k].z;
-+
-+
+          map( &dx, &dy, &dz,
-+
+               frames[i].boxX, frames[i].boxY, frames[i].boxZ );
-+
-+
+          distSqr = (dx * dx) + (dy * dy) + (dz * dz);
-+
+          dist = sqrt( distSqr );
-+
-+
+          // add to the appropriate bin
-+
+          bin = (int)( dist / delR );
-+
+          if( bin < GofRBins ) histogram[bin] += 2;
-+
+        }
-+
+      }
-+
+    }
-+
-+
+    // calculate the GofR
-+
-+
+    for( i=0; i<GofRBins; i++ ){
-+
-+
+      rLower = i * delR;
-+
+      rUpper = rLower + delR;
-+
-+
+      volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
-+
+      nIdeal = volSlice * ( allConstant );
-+
-+
+      the_GofR[i] = histogram[i] / ( nFrames * frames[0].nAtoms * nIdeal );
-+
+      rValue[i] = rLower + ( delR / 2.0 );
-+
+    }
-+
-+
+    // make the out_name;
-+
-+
+    strcpy( out_name, out_prefix );
-+
+    sprintf( tempString, "-%s-%s.GofR", "all", "all" );
-+
+    strcat( out_name, tempString );
-+
-+
+    out_file = fopen( out_name, "w" );
-+
+    if( out_file == NULL ){
-+
+      fprintf( stderr,
-+
+               "\n"
-+
+               "GofR error, unable to open \"%s\" for writing.\n",
-+
+               out_name );
-+
+      exit(8);
-+
+    }
-+
+    break;
-+
+  }
-+
-+
+  for( i=0; i<GofRBins; i++ ){
-+
+    fprintf( out_file,
-+
+             "%lf\t%lf\n",
-+
+             rValue[i], the_GofR[i] );
-+
+  }
-+
-+
+  fclose( out_file );
-+
+}
-+
-+
-+
+void map( double *x, double *y, double *z,
-+
+          double boxX, double boxY, double boxZ ){
-+
-+
+  *x -= boxX * copysign(1.0,*x) * floor( fabs( *x/boxX ) + 0.5  );
-+
+  *y -= boxY * copysign(1.0,*y) * floor( fabs( *y/boxY ) + 0.5  );
-+
+  *z -= boxZ * copysign(1.0,*z) * floor( fabs( *z/boxZ ) + 0.5  );
-+
-+
+}

Diff Legend

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

Comparing trunk/madProps/GofR.c (file contents): Revision 44 by mmeineke, Fri Jul 19 19:05:59 2002 UTC vs. Revision 45 by mmeineke, Tue Jul 23 16:08:35 2002 UTC

Diff Legend

Comparing trunk/madProps/GofR.c (file contents):
Revision 44 by mmeineke, Fri Jul 19 19:05:59 2002 UTC vs.
Revision 45 by mmeineke, Tue Jul 23 16:08:35 2002 UTC