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 startFrame, int endFrame ){

  int i,j,k;

  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 rxj, ryj, rzj;
  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++;
      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=startFrame; i<endFrame; i++){
      for( j=0; j<(frames[i].nAtoms-1); j++ ){
        
        if( !strcmp( frames[0].names[j], atom1 ) ){

          rxj = frames[i].r[j].x;
          ryj = frames[i].r[j].y;
          rzj = frames[i].r[j].z;
          
          for( k=j+1; k< frames[i].nAtoms; k++ ){
            
            if( !strcmp( frames[0].names[k], atom2 ) ){
          
              dx = rxj - frames[i].r[k].x;
              dy = ryj - frames[i].r[k].y;
              dz = rzj - 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 ) ){
          
          rxj = frames[i].r[j].x;
          ryj = frames[i].r[j].y;
          rzj = frames[i].r[j].z;
          
          for( k=j+1; k< frames[i].nAtoms; k++ ){
            
            if( !strcmp( frames[0].names[k], atom1 ) ){
              
              dx = rxj - frames[i].r[k].x;
              dy = ryj - frames[i].r[k].y;
              dz = rzj - 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=startFrame; i<endFrame; i++){
      for( j=0; j<(frames[i].nAtoms-1); j++ ){
        
        if( !strcmp( frames[0].names[j], allAtom ) ){
          
          rxj = frames[i].r[j].x;
          ryj = frames[i].r[j].y;
          rzj = frames[i].r[j].z;
          
          for( k=j+1; k< frames[i].nAtoms; k++ ){
                          
            dx = rxj - frames[i].r[k].x;
            dy = ryj - frames[i].r[k].y;
            dz = rzj - 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{
          
          rxj = frames[i].r[j].x;
          ryj = frames[i].r[j].y;
          rzj = frames[i].r[j].z;
          
          for( k=j+1; k< frames[i].nAtoms; k++ ){
            
            if( !strcmp( frames[0].names[k], allAtom ) ){
              
              dx = rxj - frames[i].r[k].x;
              dy = ryj - frames[i].r[k].y;
              dz = rzj - 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=startFrame; i<endFrame; i++){
      for( j=0; j<(frames[i].nAtoms-1); j++ ){
        
        rxj = frames[i].r[j].x;
        ryj = frames[i].r[j].y;
        rzj = frames[i].r[j].z;
        
        for( k=j+1; k< frames[i].nAtoms; k++ ){

          dx = rxj - frames[i].r[k].x;
          dy = ryj - frames[i].r[k].y;
          dz = rzj - 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 );
}

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