trunk/madProps/cosCorr.c

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "cosCorr.h"

void cosCorr( char* out_prefix, char* atom1, char* atom2,
              struct xyz_frame* frames, int nFrames ){
  
  
  int i,j,k;
    
  double atom1Dens;
  double atom2Dens;
  
  double atom1Constant;
  double atom2Constant;
  
  double nAtom1;
  double nAtom2;
  
  double delR;
  double boxVol;
  double shortBox;

  double dx, dy, dz;
  double rxj, ryj, rzj;
  double uxj, uyj, uzj;
  double uxk, uyk, uzk;
  double uj, uk;
  double distSqr;
  double dist;
  double dotProd;
  
  double rLower, rUpper;
  double nIdeal;
  double volSlice;

  int bin;
  int histogram[histBins];
  double the_cosCorr[histBins];
  double rValue[histBins];

  char out_name[500];
  char tempString[100];
  FILE *out_file;

  // 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 ) / histBins;
  boxVol = frames[0].boxX * frames[0].boxY * frames[0].boxZ;
  
  // zero the histograms;
  
  for(i=0; i<histBins; i++ ){
    
    rValue[i]      = 0.0;
    the_cosCorr[i] = 0.0;
    histogram[i]   = 0;
  }
  
  // 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"
             "cosCorr error, \"%s\" was not found in the trajectory.\n",
             atom1 );
    exit(8);
  }

  if( !nAtom2 ){
    
    fprintf( stderr, 
             "\n"
             "cosCorr 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++ ){
      
      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;
        
        // normalize the unit vector

        uxj = frames[i].v[j].x;
        uyj = frames[i].v[j].y;
        uzj = frames[i].v[j].z;

        uj = sqrt( (uxj * uxj) + (uyj * uyj) + (uzj * uzj) );
        
        uxj /= uj;
        uyj /= uj;
        uzj /= uj;

        for( k=j+1; k< frames[i].nAtoms; k++ ){
          
          if( !strcmp( frames[0].names[k], atom2 ) ){
            
            // normalize the unit vector

            uxk = frames[i].v[k].x;
            uyk = frames[i].v[k].y;
            uzk = frames[i].v[k].z;
            
            uk = sqrt( (uxk * uxk) + (uyk * uyk) + (uzk * uzk) );
            
            uxk /= uk;
            uyk /= uk;
            uzk /= uk;

            // calculate the distance
            
            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 );
            
            // calculate the cos correllation
            
            dotProd = (uxj * uxk) + (uyj * uyk) + (uzj * uzk);

            // add to the appropriate bin
            bin = (int)( dist / delR );
            if( bin < histBins ){
              histogram[bin] += 2;
              the_cosCorr[bin]   += 2.0 * dotProd;
            }
          }
        }
      }
      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;
        
        // normalize the unit vector

        uxj = frames[i].v[j].x;
        uyj = frames[i].v[j].y;
        uzj = frames[i].v[j].z;

        uj = sqrt( (uxj * uxj) + (uyj * uyj) + (uzj * uzj) );
        
        uxj /= uj;
        uyj /= uj;
        uzj /= uj;

        for( k=j+1; k< frames[i].nAtoms; k++ ){
          
          if( !strcmp( frames[0].names[k], atom1 ) ){
            
            // normalize the unit vector

            uxk = frames[i].v[k].x;
            uyk = frames[i].v[k].y;
            uzk = frames[i].v[k].z;
            
            uk = sqrt( (uxk * uxk) + (uyk * uyk) + (uzk * uzk) );
            
            uxk /= uk;
            uyk /= uk;
            uzk /= uk;

            // calculate the distance
            
            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 );
            
            // calculate the cos correllation
            
            dotProd = (uxj * uxk) + (uyj * uyk) + (uzj * uzk);

            // add to the appropriate bin
            bin = (int)( dist / delR );
            if( bin < histBins ){
              histogram[bin] += 2;
              the_cosCorr[bin]   += 2.0 * dotProd;
            }
          }
        }
      }
    }
  }
  
  // calculate the cosCorr
  
  for( i=0; i<histBins; i++ ){
    
    rLower = i * delR;
    rUpper = rLower + delR;
    
    if( histogram[i] != 0 ) the_cosCorr[i] /= histogram[i];
    
    rValue[i] = rLower + ( delR / 2.0 );
  }
  
  // make the out_name;
  
  strcpy( out_name, out_prefix );
  sprintf( tempString, "-%s-%s.cosCorr", atom1, atom2 );
  strcat( out_name, tempString );
  
  out_file = fopen( out_name, "w" );
  if( out_file == NULL ){
    fprintf( stderr,
             "\n"
             "cosCorr error, unable to open \"%s\" for writing.\n",
             out_name );
    exit(8);
  }
  
  // write out the findings

  for( i=0; i<histBins; i++ ){
    fprintf( out_file,
             "%lf\t%lf\n",
             rValue[i], the_cosCorr[i] );
  }
  
  fclose( out_file );

}

Revision:	46
Committed:	Tue Jul 23 20:10:49 2002 UTC (21 years, 11 months ago) by mmeineke
Content type:	text/plain
File size:	5427 byte(s)
Log Message:	added the cos correlation function
#	User	Rev	Content
1	mmeineke	46	#include <stdio.h>
2			#include <stdlib.h>
3			#include <string.h>
4			#include <math.h>
5
6			#include "cosCorr.h"
7
8			void cosCorr( char* out_prefix, char* atom1, char* atom2,
9			struct xyz_frame* frames, int nFrames ){
10
11
12			int i,j,k;
13
14			double atom1Dens;
15			double atom2Dens;
16
17			double atom1Constant;
18			double atom2Constant;
19
20			double nAtom1;
21			double nAtom2;
22
23			double delR;
24			double boxVol;
25			double shortBox;
26
27			double dx, dy, dz;
28			double rxj, ryj, rzj;
29			double uxj, uyj, uzj;
30			double uxk, uyk, uzk;
31			double uj, uk;
32			double distSqr;
33			double dist;
34			double dotProd;
35
36			double rLower, rUpper;
37			double nIdeal;
38			double volSlice;
39
40			int bin;
41			int histogram[histBins];
42			double the_cosCorr[histBins];
43			double rValue[histBins];
44
45			char out_name[500];
46			char tempString[100];
47			FILE *out_file;
48
49			// find the box size and delR;
50
51			shortBox = frames[0].boxX;
52			if( shortBox > frames[0].boxY ) shortBox = frames[0].boxY;
53			if( shortBox > frames[0].boxZ ) shortBox = frames[0].boxZ;
54
55			delR = ( shortBox / 2.0 ) / histBins;
56			boxVol = frames[0].boxX * frames[0].boxY * frames[0].boxZ;
57
58			// zero the histograms;
59
60			for(i=0; i<histBins; i++ ){
61
62			rValue[i] = 0.0;
63			the_cosCorr[i] = 0.0;
64			histogram[i] = 0;
65			}
66
67			// find the number of each type;
68
69			nAtom1 = 0;
70			nAtom2 = 0;
71			for( i=0; i<frames[0].nAtoms; i++ ){
72
73			if( !strcmp( frames[0].names[i], atom1 ) ) nAtom1++;
74			if( !strcmp( frames[0].names[i], atom2 ) ) nAtom2++;
75			}
76
77			if( !nAtom1 ){
78
79			fprintf( stderr,
80			"\n"
81			"cosCorr error, \"%s\" was not found in the trajectory.\n",
82			atom1 );
83			exit(8);
84			}
85
86			if( !nAtom2 ){
87
88			fprintf( stderr,
89			"\n"
90			"cosCorr error, \"%s\" was not found in the trajectory.\n",
91			atom2 );
92			exit(8);
93			}
94
95			// calculate some of the constants;
96
97			atom1Dens = nAtom1 / boxVol;
98			atom2Dens = nAtom2 / boxVol;
99
100			atom1Constant = ( 4.0 * M_PI * atom1Dens ) / 3.0;
101			atom2Constant = ( 4.0 * M_PI * atom2Dens ) / 3.0;
102
103			// calculate the histogram
104
105			for( i=0; i<nFrames; i++){
106			for( j=0; j<(frames[i].nAtoms-1); j++ ){
107
108			if( !strcmp( frames[0].names[j], atom1 ) ){
109
110			rxj = frames[i].r[j].x;
111			ryj = frames[i].r[j].y;
112			rzj = frames[i].r[j].z;
113
114			// normalize the unit vector
115
116			uxj = frames[i].v[j].x;
117			uyj = frames[i].v[j].y;
118			uzj = frames[i].v[j].z;
119
120			uj = sqrt( (uxj * uxj) + (uyj * uyj) + (uzj * uzj) );
121
122			uxj /= uj;
123			uyj /= uj;
124			uzj /= uj;
125
126			for( k=j+1; k< frames[i].nAtoms; k++ ){
127
128			if( !strcmp( frames[0].names[k], atom2 ) ){
129
130			// normalize the unit vector
131
132			uxk = frames[i].v[k].x;
133			uyk = frames[i].v[k].y;
134			uzk = frames[i].v[k].z;
135
136			uk = sqrt( (uxk * uxk) + (uyk * uyk) + (uzk * uzk) );
137
138			uxk /= uk;
139			uyk /= uk;
140			uzk /= uk;
141
142			// calculate the distance
143
144			dx = rxj - frames[i].r[k].x;
145			dy = ryj - frames[i].r[k].y;
146			dz = rzj - frames[i].r[k].z;
147
148			map( &dx, &dy, &dz,
149			frames[i].boxX, frames[i].boxY, frames[i].boxZ );
150
151			distSqr = (dx * dx) + (dy * dy) + (dz * dz);
152			dist = sqrt( distSqr );
153
154			// calculate the cos correllation
155
156			dotProd = (uxj * uxk) + (uyj * uyk) + (uzj * uzk);
157
158			// add to the appropriate bin
159			bin = (int)( dist / delR );
160			if( bin < histBins ){
161			histogram[bin] += 2;
162			the_cosCorr[bin] += 2.0 * dotProd;
163			}
164			}
165			}
166			}
167			else if( !strcmp( frames[0].names[j], atom2 ) ){
168
169			rxj = frames[i].r[j].x;
170			ryj = frames[i].r[j].y;
171			rzj = frames[i].r[j].z;
172
173			// normalize the unit vector
174
175			uxj = frames[i].v[j].x;
176			uyj = frames[i].v[j].y;
177			uzj = frames[i].v[j].z;
178
179			uj = sqrt( (uxj * uxj) + (uyj * uyj) + (uzj * uzj) );
180
181			uxj /= uj;
182			uyj /= uj;
183			uzj /= uj;
184
185			for( k=j+1; k< frames[i].nAtoms; k++ ){
186
187			if( !strcmp( frames[0].names[k], atom1 ) ){
188
189			// normalize the unit vector
190
191			uxk = frames[i].v[k].x;
192			uyk = frames[i].v[k].y;
193			uzk = frames[i].v[k].z;
194
195			uk = sqrt( (uxk * uxk) + (uyk * uyk) + (uzk * uzk) );
196
197			uxk /= uk;
198			uyk /= uk;
199			uzk /= uk;
200
201			// calculate the distance
202
203			dx = rxj - frames[i].r[k].x;
204			dy = ryj - frames[i].r[k].y;
205			dz = rzj - frames[i].r[k].z;
206
207			map( &dx, &dy, &dz,
208			frames[i].boxX, frames[i].boxY, frames[i].boxZ );
209
210			distSqr = (dx * dx) + (dy * dy) + (dz * dz);
211			dist = sqrt( distSqr );
212
213			// calculate the cos correllation
214
215			dotProd = (uxj * uxk) + (uyj * uyk) + (uzj * uzk);
216
217			// add to the appropriate bin
218			bin = (int)( dist / delR );
219			if( bin < histBins ){
220			histogram[bin] += 2;
221			the_cosCorr[bin] += 2.0 * dotProd;
222			}
223			}
224			}
225			}
226			}
227			}
228
229			// calculate the cosCorr
230
231			for( i=0; i<histBins; i++ ){
232
233			rLower = i * delR;
234			rUpper = rLower + delR;
235
236			if( histogram[i] != 0 ) the_cosCorr[i] /= histogram[i];
237
238			rValue[i] = rLower + ( delR / 2.0 );
239			}
240
241			// make the out_name;
242
243			strcpy( out_name, out_prefix );
244			sprintf( tempString, "-%s-%s.cosCorr", atom1, atom2 );
245			strcat( out_name, tempString );
246
247			out_file = fopen( out_name, "w" );
248			if( out_file == NULL ){
249			fprintf( stderr,
250			"\n"
251			"cosCorr error, unable to open \"%s\" for writing.\n",
252			out_name );
253			exit(8);
254			}
255
256			// write out the findings
257
258			for( i=0; i<histBins; i++ ){
259			fprintf( out_file,
260			"%lf\t%lf\n",
261			rValue[i], the_cosCorr[i] );
262			}
263
264			fclose( out_file );
265
266			}
267