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