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