ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/madProps/GofR.c
(Generate patch)

Comparing trunk/madProps/GofR.c (file contents):
Revision 45 by mmeineke, Tue Jul 23 16:08:35 2002 UTC vs.
Revision 103 by mmeineke, Thu Aug 29 17:21:54 2002 UTC

# Line 12 | Line 12 | void GofR( char* out_prefix, char* atom1, char* atom2,
12   // constant boxSize
13  
14   void GofR( char* out_prefix, char* atom1, char* atom2,
15 <           struct xyz_frame* frames, int nFrames ){
15 >           struct xyz_frame* frames, int nFrames,
16 >           int startFrame, int endFrame ){
17  
18 <  int i,j,k,l;
18 >  int i,j,k;
19  
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;
25 <  double allDens;
24 >  double pairDens;
25 >  double pairConstant;
26  
27 –  double atom1Constant;
28 –  double atom2Constant;
29 –  double allConstant;
30 –
27    double nAtom1;
28    double nAtom2;
29    double nAtoms;
# Line 36 | Line 32 | void GofR( char* out_prefix, char* atom1, char* atom2,
32    double boxVol;
33    double shortBox;
34  
35 +  double rxj, ryj, rzj;
36    double dx, dy, dz;
37    double distSqr;
38    double dist;
# Line 102 | Line 99 | void GofR( char* out_prefix, char* atom1, char* atom2,
99      for( i=0; i<frames[0].nAtoms; i++ ){
100    
101        if( !strcmp( frames[0].names[i], atom1 ) ) nAtom1++;
102 <      else if( !strcmp( frames[0].names[i], atom2 ) ) nAtom2++;
102 >      if( !strcmp( frames[0].names[i], atom2 ) ) nAtom2++;
103      }
104      
105      if( !nAtom1 ){
# Line 125 | Line 122 | void GofR( char* out_prefix, char* atom1, char* atom2,
122  
123      // calculate some of the constants;
124      
125 <    atom1Dens = nAtom1 / boxVol;
126 <    atom2Dens = nAtom2 / boxVol;
127 <    
128 <    atom1Constant = ( 4.0 * M_PI * atom1Dens ) / 3.0;
129 <    atom2Constant = ( 4.0 * M_PI * atom2Dens ) / 3.0;
133 <    
125 >    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      // calculate the histogram
131  
132 <    for( i=0; i<nFrames; i++){
132 >    for( i=startFrame; i<endFrame; i++){
133        for( j=0; j<(frames[i].nAtoms-1); j++ ){
134 <        for( k=j+1; k< frames[i].nAtoms; k++ ){
134 >        
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            
141 <          if( !strcmp( frames[0].names[j], atom1 ) ){
141 >          for( k=j+1; k< frames[i].nAtoms; k++ ){
142 >            
143              if( !strcmp( frames[0].names[k], atom2 ) ){
144            
145 <              dx = frames[i].r[j].x - frames[i].r[k].x;
146 <              dy = frames[i].r[j].y - frames[i].r[k].y;
147 <              dz = frames[i].r[j].z - frames[i].r[k].z;
145 >              dx = rxj - frames[i].r[k].x;
146 >              dy = ryj - frames[i].r[k].y;
147 >              dz = rzj - frames[i].r[k].z;
148                
149                map( &dx, &dy, &dz,
150                     frames[i].boxX, frames[i].boxY, frames[i].boxZ );
# Line 155 | Line 157 | void GofR( char* out_prefix, char* atom1, char* atom2,
157                if( bin < GofRBins ) histogram[bin] += 2;
158              }
159            }
160 <
161 <          else if( !strcmp( frames[0].names[j], atom2 ) ){
160 >        }
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              if( !strcmp( frames[0].names[k], atom1 ) ){
171                
172 <              dx = frames[i].r[j].x - frames[i].r[k].x;
173 <              dy = frames[i].r[j].y - frames[i].r[k].y;
174 <              dz = frames[i].r[j].z - frames[i].r[k].z;
172 >              dx = rxj - frames[i].r[k].x;
173 >              dy = ryj - frames[i].r[k].y;
174 >              dz = rzj - frames[i].r[k].z;
175                
176                map( &dx, &dy, &dz,
177                     frames[i].boxX, frames[i].boxY, frames[i].boxZ );
# Line 186 | Line 196 | void GofR( char* out_prefix, char* atom1, char* atom2,
196        rUpper = rLower + delR;
197        
198        volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
199 <      nIdeal = volSlice * ( atom1Constant + atom2Constant );
199 >      nIdeal = volSlice * pairConstant;
200        
201 <      the_GofR[i] = histogram[i] / ( nFrames * ( nAtom1 + nAtom2 ) * nIdeal );
201 >      the_GofR[i] = histogram[i] / ( nFrames * nIdeal );
202        rValue[i] = rLower + ( delR / 2.0 );
203      }
204  
# Line 229 | Line 239 | void GofR( char* out_prefix, char* atom1, char* atom2,
239  
240      // calculate some of the constants;
241      
242 <    atom1Dens = nAtom1 / boxVol;
243 <    allDens  = frames[0].nAtoms / boxVol;
244 <    
235 <    atom1Constant = ( 4.0 * M_PI * atom1Dens ) / 3.0;
236 <    allConstant  = ( 4.0 * M_PI * allDens ) / 3.0;
237 <    
242 >    pairDens = frames[0].nAtoms * (nAtom1 - 1) / boxVol;
243 >    pairConstant =  ( 4.0 * M_PI * pairDens ) / 3.0;
244 >
245      // calculate the histogram
246  
247 <    for( i=0; i<nFrames; i++){
247 >    for( i=startFrame; i<endFrame; i++){
248        for( j=0; j<(frames[i].nAtoms-1); j++ ){
249 <        for( k=j+1; k< frames[i].nAtoms; k++ ){
249 >        
250 >        if( !strcmp( frames[0].names[j], allAtom ) ){
251            
252 <          if( !strcmp( frames[0].names[j], allAtom ) ){
253 <                  
254 <            dx = frames[i].r[j].x - frames[i].r[k].x;
255 <            dy = frames[i].r[j].y - frames[i].r[k].y;
256 <            dz = frames[i].r[j].z - frames[i].r[k].z;
252 >          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              
262              map( &dx, &dy, &dz,
263                   frames[i].boxX, frames[i].boxY, frames[i].boxZ );
# Line 257 | Line 269 | void GofR( char* out_prefix, char* atom1, char* atom2,
269              bin = (int)( dist / delR );
270              if( bin < GofRBins ) histogram[bin] += 2;
271            }
272 +        }
273            
274 <          else if( !strcmp( frames[0].names[k], allAtom ) ){
274 >        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              
282 <            dx = frames[i].r[j].x - frames[i].r[k].x;
283 <            dy = frames[i].r[j].y - frames[i].r[k].y;
284 <            dz = frames[i].r[j].z - frames[i].r[k].z;
285 <            
286 <            map( &dx, &dy, &dz,
287 <                 frames[i].boxX, frames[i].boxY, frames[i].boxZ );
288 <            
289 <            distSqr = (dx * dx) + (dy * dy) + (dz * dz);
290 <            dist = sqrt( distSqr );
291 <            
292 <            // add to the appropriate bin
293 <            bin = (int)( dist / delR );
294 <            if( bin < GofRBins ) histogram[bin] += 2;
282 >            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            }
299          }        
300        }
# Line 286 | Line 308 | void GofR( char* out_prefix, char* atom1, char* atom2,
308        rUpper = rLower + delR;
309        
310        volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
311 <      nIdeal = volSlice * ( allConstant );
311 >      nIdeal = volSlice * pairConstant;
312        
313 <      the_GofR[i] = histogram[i] / ( nFrames * frames[0].nAtoms * nIdeal );
313 >      the_GofR[i] = histogram[i] / ( nFrames * nIdeal );
314        rValue[i] = rLower + ( delR / 2.0 );
315      }
316  
# Line 312 | Line 334 | void GofR( char* out_prefix, char* atom1, char* atom2,
334      
335      // calculate some of the constants;
336      
337 <    allDens  = frames[0].nAtoms / boxVol;
337 >    pairDens  = frames[0].nAtoms * (frames[0].nAtoms - 1) / boxVol;
338 >    pairConstant  = ( 4.0 * M_PI * pairDens ) / 3.0;
339      
317 –    allConstant  = ( 4.0 * M_PI * allDens ) / 3.0;
318 –    
340      // calculate the histogram
341  
342 <    for( i=0; i<nFrames; i++){
342 >    for( i=startFrame; i<endFrame; i++){
343        for( j=0; j<(frames[i].nAtoms-1); j++ ){
344 +        
345 +        rxj = frames[i].r[j].x;
346 +        ryj = frames[i].r[j].y;
347 +        rzj = frames[i].r[j].z;
348 +        
349          for( k=j+1; k< frames[i].nAtoms; k++ ){
350 <          
351 <          dx = frames[i].r[j].x - frames[i].r[k].x;
352 <          dy = frames[i].r[j].y - frames[i].r[k].y;
353 <          dz = frames[i].r[j].z - frames[i].r[k].z;
354 <          
350 >
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            map( &dx, &dy, &dz,
356                 frames[i].boxX, frames[i].boxY, frames[i].boxZ );
357            
# Line 347 | Line 373 | void GofR( char* out_prefix, char* atom1, char* atom2,
373        rUpper = rLower + delR;
374        
375        volSlice = pow( rUpper, 3.0 ) - pow( rLower, 3.0 );
376 <      nIdeal = volSlice * ( allConstant );
376 >      nIdeal = volSlice * pairConstant;
377        
378 <      the_GofR[i] = histogram[i] / ( nFrames * frames[0].nAtoms * nIdeal );
378 >      the_GofR[i] = histogram[i] / ( nFrames * nIdeal );
379        rValue[i] = rLower + ( delR / 2.0 );
380      }
381  
# Line 379 | Line 405 | void GofR( char* out_prefix, char* atom1, char* atom2,
405    fclose( out_file );
406   }
407  
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 – }

Diff Legend

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