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 pairDens;
  double pairConstant;

  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;
    
    if( !strcmp( atom1, atom2 ) ) pairDens = nAtom1 * (nAtom1 - 1) / boxVol;
    else pairDens = nAtom1 * nAtom2 / boxVol;

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