trunk/madProps/GofR.c

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#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  );

}
Revision:	45
Committed:	Tue Jul 23 16:08:35 2002 UTC (21 years, 11 months ago) by mmeineke
Content type:	text/plain
File size:	9168 byte(s)
Log Message:	got the GofR code working, and slimmed down the memory usage.
#	User	Rev	Content
1	mmeineke	43	#include <math.h>
2			#include <stdlib.h>
3			#include <stdio.h>
4			#include <string.h>
5
6			#include "GofR.h"
7	mmeineke	45
8			void map( double x, double y, double *z,
9			double boxX, double boxY, double boxZ );
10
11			// this algoritm assumes constant number of atoms between frames, and
12			// constant boxSize
13
14			void GofR( char* out_prefix, char* atom1, char* atom2,
15			struct xyz_frame* frames, int nFrames ){
16
17			int i,j,k,l;
18
19			enum g_types { all_all, atom_all, atom_atom };
20			enum g_types the_type;
21			char* allAtom;
22
23			double atom1Dens;
24			double atom2Dens;
25			double allDens;
26
27			double atom1Constant;
28			double atom2Constant;
29			double allConstant;
30
31			double nAtom1;
32			double nAtom2;
33			double nAtoms;
34
35			double delR;
36			double boxVol;
37			double shortBox;
38
39			double dx, dy, dz;
40			double distSqr;
41			double dist;
42
43			double rLower, rUpper;
44			double nIdeal;
45			double volSlice;
46
47			int bin;
48			int histogram[GofRBins];
49			double the_GofR[GofRBins];
50			double rValue[GofRBins];
51
52			char out_name[500];
53			char tempString[100];
54			FILE *out_file;
55
56
57			// figure out the type of GofR we are calculating
58
59			if( !strcmp( atom1, "all" ) ){
60
61			if( !strcmp( atom2, "all" ) ) the_type = all_all;
62			else{
63			the_type = atom_all;
64			allAtom = atom2;
65			}
66			}
67			else{
68
69			if( !strcmp( atom2, "all" ) ){
70			the_type = atom_all;
71			allAtom = atom1;
72			}
73			else the_type = atom_atom;
74			}
75
76			// find the box size and delR;
77
78			shortBox = frames[0].boxX;
79			if( shortBox > frames[0].boxY ) shortBox = frames[0].boxY;
80			if( shortBox > frames[0].boxZ ) shortBox = frames[0].boxZ;
81
82			delR = ( shortBox / 2.0 ) / GofRBins;
83			boxVol = frames[0].boxX * frames[0].boxY * frames[0].boxZ;
84
85			// zero the histograms;
86
87			for(i=0; i<GofRBins; i++ ){
88
89			rValue[i] = 0.0;
90			the_GofR[i] = 0.0;
91			histogram[i] = 0;
92			}
93
94			switch( the_type ){
95
96			case atom_atom:
97
98			// find the number of each type;
99
100			nAtom1 = 0;
101			nAtom2 = 0;
102			for( i=0; i<frames[0].nAtoms; i++ ){
103
104			if( !strcmp( frames[0].names[i], atom1 ) ) nAtom1++;
105			else if( !strcmp( frames[0].names[i], atom2 ) ) nAtom2++;
106			}
107
108			if( !nAtom1 ){
109
110			fprintf( stderr,
111			"\n"
112			"GofR error, \"%s\" was not found in the trajectory.\n",
113			atom1 );
114			exit(8);
115			}
116
117			if( !nAtom2 ){
118
119			fprintf( stderr,
120			"\n"
121			"GofR error, \"%s\" was not found in the trajectory.\n",
122			atom2 );
123			exit(8);
124			}
125
126			// calculate some of the constants;
127
128			atom1Dens = nAtom1 / boxVol;
129			atom2Dens = nAtom2 / boxVol;
130
131			atom1Constant = ( 4.0 * M_PI * atom1Dens ) / 3.0;
132			atom2Constant = ( 4.0 * M_PI * atom2Dens ) / 3.0;
133
134			// calculate the histogram
135
136			for( i=0; i<nFrames; i++){
137			for( j=0; j<(frames[i].nAtoms-1); j++ ){
138			for( k=j+1; k< frames[i].nAtoms; k++ ){
139
140			if( !strcmp( frames[0].names[j], atom1 ) ){
141			if( !strcmp( frames[0].names[k], atom2 ) ){
142
143			dx = frames[i].r[j].x - frames[i].r[k].x;
144			dy = frames[i].r[j].y - frames[i].r[k].y;
145			dz = frames[i].r[j].z - frames[i].r[k].z;
146
147			map( &dx, &dy, &dz,
148			frames[i].boxX, frames[i].boxY, frames[i].boxZ );
149
150			distSqr = (dx * dx) + (dy * dy) + (dz * dz);
151			dist = sqrt( distSqr );
152
153			// add to the appropriate bin
154			bin = (int)( dist / delR );
155			if( bin < GofRBins ) histogram[bin] += 2;
156			}
157			}
158
159			else if( !strcmp( frames[0].names[j], atom2 ) ){
160			if( !strcmp( frames[0].names[k], atom1 ) ){
161
162			dx = frames[i].r[j].x - frames[i].r[k].x;
163			dy = frames[i].r[j].y - frames[i].r[k].y;
164			dz = frames[i].r[j].z - frames[i].r[k].z;
165
166			map( &dx, &dy, &dz,
167			frames[i].boxX, frames[i].boxY, frames[i].boxZ );
168
169			distSqr = (dx * dx) + (dy * dy) + (dz * dz);
170			dist = sqrt( distSqr );
171
172			// add to the appropriate bin
173			bin = (int)( dist / delR );
174			if( bin < GofRBins ) histogram[bin] += 2;
175			}
176			}
177			}
178			}
179			}
180
181			// calculate the GofR
182
183			for( i=0; i<GofRBins; i++ ){
184
185			rLower = i * delR;
186			rUpper = rLower + delR;
187
188			volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
189			nIdeal = volSlice * ( atom1Constant + atom2Constant );
190
191			the_GofR[i] = histogram[i] / ( nFrames * ( nAtom1 + nAtom2 ) * nIdeal );
192			rValue[i] = rLower + ( delR / 2.0 );
193			}
194
195			// make the out_name;
196
197			strcpy( out_name, out_prefix );
198			sprintf( tempString, "-%s-%s.GofR", atom1, atom2 );
199			strcat( out_name, tempString );
200
201			out_file = fopen( out_name, "w" );
202			if( out_file == NULL ){
203			fprintf( stderr,
204			"\n"
205			"GofR error, unable to open \"%s\" for writing.\n",
206			out_name );
207			exit(8);
208			}
209			break;
210
211			case atom_all:
212
213			// find the number of AllAtoms;
214
215			nAtom1 = 0;
216			for( i=0; i<frames[0].nAtoms; i++ ){
217
218			if( !strcmp( frames[0].names[i], allAtom ) ) nAtom1++;
219			}
220
221			if( !nAtom1 ){
222
223			fprintf( stderr,
224			"\n"
225			"GofR error, \"%s\" was not found in the trajectory.\n",
226			atom1 );
227			exit(8);
228			}
229
230			// calculate some of the constants;
231
232			atom1Dens = nAtom1 / boxVol;
233			allDens = frames[0].nAtoms / boxVol;
234
235			atom1Constant = ( 4.0 * M_PI * atom1Dens ) / 3.0;
236			allConstant = ( 4.0 * M_PI * allDens ) / 3.0;
237
238			// calculate the histogram
239
240			for( i=0; i<nFrames; i++){
241			for( j=0; j<(frames[i].nAtoms-1); j++ ){
242			for( k=j+1; k< frames[i].nAtoms; k++ ){
243
244			if( !strcmp( frames[0].names[j], allAtom ) ){
245
246			dx = frames[i].r[j].x - frames[i].r[k].x;
247			dy = frames[i].r[j].y - frames[i].r[k].y;
248			dz = frames[i].r[j].z - frames[i].r[k].z;
249
250			map( &dx, &dy, &dz,
251			frames[i].boxX, frames[i].boxY, frames[i].boxZ );
252
253			distSqr = (dx * dx) + (dy * dy) + (dz * dz);
254			dist = sqrt( distSqr );
255
256			// add to the appropriate bin
257			bin = (int)( dist / delR );
258			if( bin < GofRBins ) histogram[bin] += 2;
259			}
260
261			else if( !strcmp( frames[0].names[k], allAtom ) ){
262
263			dx = frames[i].r[j].x - frames[i].r[k].x;
264			dy = frames[i].r[j].y - frames[i].r[k].y;
265			dz = frames[i].r[j].z - frames[i].r[k].z;
266
267			map( &dx, &dy, &dz,
268			frames[i].boxX, frames[i].boxY, frames[i].boxZ );
269
270			distSqr = (dx * dx) + (dy * dy) + (dz * dz);
271			dist = sqrt( distSqr );
272
273			// add to the appropriate bin
274			bin = (int)( dist / delR );
275			if( bin < GofRBins ) histogram[bin] += 2;
276			}
277			}
278			}
279			}
280
281			// calculate the GofR
282
283			for( i=0; i<GofRBins; i++ ){
284
285			rLower = i * delR;
286			rUpper = rLower + delR;
287
288			volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
289			nIdeal = volSlice * ( allConstant );
290
291			the_GofR[i] = histogram[i] / ( nFrames * frames[0].nAtoms * nIdeal );
292			rValue[i] = rLower + ( delR / 2.0 );
293			}
294
295			// make the out_name;
296
297			strcpy( out_name, out_prefix );
298			sprintf( tempString, "-%s-%s.GofR", allAtom, "all" );
299			strcat( out_name, tempString );
300
301			out_file = fopen( out_name, "w" );
302			if( out_file == NULL ){
303			fprintf( stderr,
304			"\n"
305			"GofR error, unable to open \"%s\" for writing.\n",
306			out_name );
307			exit(8);
308			}
309			break;
310
311			case all_all:
312
313			// calculate some of the constants;
314
315			allDens = frames[0].nAtoms / boxVol;
316
317			allConstant = ( 4.0 * M_PI * allDens ) / 3.0;
318
319			// calculate the histogram
320
321			for( i=0; i<nFrames; i++){
322			for( j=0; j<(frames[i].nAtoms-1); j++ ){
323			for( k=j+1; k< frames[i].nAtoms; k++ ){
324
325			dx = frames[i].r[j].x - frames[i].r[k].x;
326			dy = frames[i].r[j].y - frames[i].r[k].y;
327			dz = frames[i].r[j].z - frames[i].r[k].z;
328
329			map( &dx, &dy, &dz,
330			frames[i].boxX, frames[i].boxY, frames[i].boxZ );
331
332			distSqr = (dx * dx) + (dy * dy) + (dz * dz);
333			dist = sqrt( distSqr );
334
335			// add to the appropriate bin
336			bin = (int)( dist / delR );
337			if( bin < GofRBins ) histogram[bin] += 2;
338			}
339			}
340			}
341
342			// calculate the GofR
343
344			for( i=0; i<GofRBins; i++ ){
345
346			rLower = i * delR;
347			rUpper = rLower + delR;
348
349			volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
350			nIdeal = volSlice * ( allConstant );
351
352			the_GofR[i] = histogram[i] / ( nFrames * frames[0].nAtoms * nIdeal );
353			rValue[i] = rLower + ( delR / 2.0 );
354			}
355
356			// make the out_name;
357
358			strcpy( out_name, out_prefix );
359			sprintf( tempString, "-%s-%s.GofR", "all", "all" );
360			strcat( out_name, tempString );
361
362			out_file = fopen( out_name, "w" );
363			if( out_file == NULL ){
364			fprintf( stderr,
365			"\n"
366			"GofR error, unable to open \"%s\" for writing.\n",
367			out_name );
368			exit(8);
369			}
370			break;
371			}
372
373			for( i=0; i<GofRBins; i++ ){
374			fprintf( out_file,
375			"%lf\t%lf\n",
376			rValue[i], the_GofR[i] );
377			}
378
379			fclose( out_file );
380			}
381
382
383			void map( double x, double y, double *z,
384			double boxX, double boxY, double boxZ ){
385
386			x -= boxX copysign(1.0,x) floor( fabs( *x/boxX ) + 0.5 );
387			y -= boxY copysign(1.0,y) floor( fabs( *y/boxY ) + 0.5 );
388			z -= boxZ copysign(1.0,z) floor( fabs( *z/boxZ ) + 0.5 );
389
390			}