trunk/tcProps/gofz.c

#define _FILE_OFFSET_BITS 64

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


#include "params.h"
#include "tcProps.h"
#include "readWrite.h"
#include "gofz.h"


struct gofzStr{
  double gofz[GOFZNBINS];
};

struct gofzStr gofzCorrs[5];
double gofzDl;
int halfGzBins;

void accumGofZ( struct atomCoord *atoms, double Hmat[3][3], double maxL );

double wrapMyZ( double z, double box );

void calcGofz(double startTime, struct atomCoord *atoms,
              char* outPrefix ){

    // list of 'a priori' constants

  const int nLipAtoms = NL_ATOMS;
  const int nBonds    = NBONDS;
  const int nLipids   = NLIPIDS;
  const int nSSD      = NSSD;
  const int nAtoms    = nLipAtoms * nLipids + nSSD;

  // variables

  char outName[500];
  FILE* outFile;
  int i, j, k;
  int startFrame, corrFrames, framesFinished;
  int startFound, percentComplete;
  
  double Hmat[3][3];
  double maxL, halfMax;

  framesFinished = 0;
  
  startFound = 0;
  startFrame = -1;
  while( !startFound ){

    startFrame++;

    if(startFrame >= nFrames){
      
      fprintf( stderr,
               "Start Time, %G, was not found in the dump file.\n",
               startTime );
      exit(0);
    }
    
    if(startTime <= frameTimes[startFrame])
      startFound = 1;


  }

  corrFrames = nFrames - startFrame;
  
  maxL = scanSmallestZ( startFrame );
  halfMax = maxL * 0.5;
  gofzDl = maxL / GOFZNBINS;
  halfGzBins = (GOFZNBINS-1)*0.5;
  
  for(i=0;i<GOFZNBINS;i++)
    for(j=0;j<5;j++)
      gofzCorrs[j].gofz[i] = 0.0;
  

  for(i=startFrame; i<nFrames; i++){

    percentComplete = 
      (int)( 100.0 * (double)framesFinished / (double) corrFrames );
    
    fprintf( stdout,
             "\rg(z) corr %3d%% complete.",
             percentComplete );
    fflush( stdout );

    readFrame( i, atoms, Hmat );
    
    accumGofZ( atoms, Hmat, maxL );

    framesFinished++;
  }

  for(i=0;i<GOFZNBINS;i++)
    for(j=0;j<5;j++)
      gofzCorrs[j].gofz[i] /= corrFrames;

  sprintf( outName, "%s.gofz", outPrefix );
  outFile = fopen( outName, "w" );

  fprintf( outFile,
           "#z\tHEAD\tCH\tCH2\tCH3\tSSD\n");
  for(i=0;i<GOFZNBINS;i++)
    fprintf( outFile,
             "%6G\t%6G\t%6G\t%6G\t%6G\t%6G\n",
             -halfMax + i*gofzDl,
             gofzCorrs[HEAD].gofz[i],
             gofzCorrs[CH].gofz[i],
             gofzCorrs[CH2].gofz[i],
             gofzCorrs[CH3].gofz[i],
             gofzCorrs[SSD].gofz[i] );
  
  fflush(outFile);
  fclose(outFile);


  percentComplete = 
    (int)( 100.0 * (double)framesFinished / (double) corrFrames );
  
  fprintf( stdout,
           "\rg(z) corr %3d%% complete.\n"
           "done.\n",
           percentComplete );
  fflush( stdout );
}

void accumGofZ( struct atomCoord *atoms, double Hmat[3][3], double maxL ){
  
  // list of 'a priori' constants

  const int nLipAtoms = NL_ATOMS;
  const int nBonds    = NBONDS;
  const int nLipids   = NLIPIDS;
  const int nSSD      = NSSD;
  const int nAtoms    = nLipAtoms * nLipids + nSSD;

  int i,j,k,l;
  int bin;

  double totMass;
  double comZ;
  double halfMax;
  double wrapZ;
  double scale;
  double volume, density, densCons;
  
  struct gofzStr tempCorrs[5];
  
  for(i=0;i<GOFZNBINS;i++)
    for(j=0;j<5;j++)
      tempCorrs[j].gofz[i] = 0.0;

  // calculate comZ;
  
  totMass = 0.0;
  comZ = 0.0;
  for(i=0;i<nAtoms;i++){
    comZ += atoms[i].pos[2]*atoms[i].mass;
    totMass += atoms[i].mass;
  }

  comZ /= totMass;

  halfMax = maxL * 0.5;

  // calc the gofZ;

  for(i=0;i<nAtoms;i++){
    
    wrapZ = wrapMyZ( (atoms[i].pos[2] - comZ), Hmat[2][2] );
    
    if(fabs(wrapZ) < halfMax ){
      
      scale = wrapZ / halfMax;
      bin = (int)( (scale*halfGzBins) + halfGzBins );
      
      tempCorrs[atoms[i].type].gofz[bin] += 1.0;
    }
  }

  volume = Hmat[0][0] * Hmat[1][1] * Hmat[2][2];
  density = (double)nAtoms / volume;
  
  densCons = Hmat[0][0] * Hmat[1][1] * gofzDl * density;

  for(i=0;i<GOFZNBINS;i++)
    for(j=0;j<5;j++)
      tempCorrs[j].gofz[i] /= densCons;
  

  for(i=0;i<GOFZNBINS;i++)
    for(j=0;j<5;j++)
      gofzCorrs[j].gofz[i] += tempCorrs[j].gofz[i];
}


double wrapMyZ( double z, double box ){

  double scale, round;

  scale = z / box;
  
  round = ( scale >= 0 ) ? floor( scale + 0.5 ) : ceil( scale - 0.5 );
  
  scale -= round;
  
  return scale*box;
}
Revision:	1080
Committed:	Wed Mar 3 15:16:15 2004 UTC (21 years, 8 months ago) by mmeineke
Content type:	text/plain
File size:	4231 byte(s)
Log Message:	added gofz code
#	User	Rev	Content
1	mmeineke	1080	#define _FILE_OFFSET_BITS 64
2
3			#include <stdio.h>
4			#include <stdlib.h>
5			#include <string.h>
6			#include <math.h>
7
8
9			#include "params.h"
10			#include "tcProps.h"
11			#include "readWrite.h"
12			#include "gofz.h"
13
14
15
16			struct gofzStr{
17			double gofz[GOFZNBINS];
18			};
19
20			struct gofzStr gofzCorrs[5];
21			double gofzDl;
22			int halfGzBins;
23
24			void accumGofZ( struct atomCoord *atoms, double Hmat[3][3], double maxL );
25
26			double wrapMyZ( double z, double box );
27
28			void calcGofz(double startTime, struct atomCoord *atoms,
29			char* outPrefix ){
30
31			// list of 'a priori' constants
32
33			const int nLipAtoms = NL_ATOMS;
34			const int nBonds = NBONDS;
35			const int nLipids = NLIPIDS;
36			const int nSSD = NSSD;
37			const int nAtoms = nLipAtoms * nLipids + nSSD;
38
39			// variables
40
41			char outName[500];
42			FILE* outFile;
43			int i, j, k;
44			int startFrame, corrFrames, framesFinished;
45			int startFound, percentComplete;
46
47			double Hmat[3][3];
48			double maxL, halfMax;
49
50			framesFinished = 0;
51
52			startFound = 0;
53			startFrame = -1;
54			while( !startFound ){
55
56			startFrame++;
57
58			if(startFrame >= nFrames){
59
60			fprintf( stderr,
61			"Start Time, %G, was not found in the dump file.\n",
62			startTime );
63			exit(0);
64			}
65
66			if(startTime <= frameTimes[startFrame])
67			startFound = 1;
68
69
70			}
71
72			corrFrames = nFrames - startFrame;
73
74			maxL = scanSmallestZ( startFrame );
75			halfMax = maxL * 0.5;
76			gofzDl = maxL / GOFZNBINS;
77			halfGzBins = (GOFZNBINS-1)*0.5;
78
79			for(i=0;i<GOFZNBINS;i++)
80			for(j=0;j<5;j++)
81			gofzCorrs[j].gofz[i] = 0.0;
82
83
84			for(i=startFrame; i<nFrames; i++){
85
86			percentComplete =
87			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
88
89			fprintf( stdout,
90			"\rg(z) corr %3d%% complete.",
91			percentComplete );
92			fflush( stdout );
93
94			readFrame( i, atoms, Hmat );
95
96			accumGofZ( atoms, Hmat, maxL );
97
98			framesFinished++;
99			}
100
101			for(i=0;i<GOFZNBINS;i++)
102			for(j=0;j<5;j++)
103			gofzCorrs[j].gofz[i] /= corrFrames;
104
105			sprintf( outName, "%s.gofz", outPrefix );
106			outFile = fopen( outName, "w" );
107
108			fprintf( outFile,
109			"#z\tHEAD\tCH\tCH2\tCH3\tSSD\n");
110			for(i=0;i<GOFZNBINS;i++)
111			fprintf( outFile,
112			"%6G\t%6G\t%6G\t%6G\t%6G\t%6G\n",
113			-halfMax + i*gofzDl,
114			gofzCorrs[HEAD].gofz[i],
115			gofzCorrs[CH].gofz[i],
116			gofzCorrs[CH2].gofz[i],
117			gofzCorrs[CH3].gofz[i],
118			gofzCorrs[SSD].gofz[i] );
119
120			fflush(outFile);
121			fclose(outFile);
122
123
124			percentComplete =
125			(int)( 100.0 * (double)framesFinished / (double) corrFrames );
126
127			fprintf( stdout,
128			"\rg(z) corr %3d%% complete.\n"
129			"done.\n",
130			percentComplete );
131			fflush( stdout );
132			}
133
134			void accumGofZ( struct atomCoord *atoms, double Hmat[3][3], double maxL ){
135
136			// list of 'a priori' constants
137
138			const int nLipAtoms = NL_ATOMS;
139			const int nBonds = NBONDS;
140			const int nLipids = NLIPIDS;
141			const int nSSD = NSSD;
142			const int nAtoms = nLipAtoms * nLipids + nSSD;
143
144			int i,j,k,l;
145			int bin;
146
147			double totMass;
148			double comZ;
149			double halfMax;
150			double wrapZ;
151			double scale;
152			double volume, density, densCons;
153
154			struct gofzStr tempCorrs[5];
155
156			for(i=0;i<GOFZNBINS;i++)
157			for(j=0;j<5;j++)
158			tempCorrs[j].gofz[i] = 0.0;
159
160			// calculate comZ;
161
162			totMass = 0.0;
163			comZ = 0.0;
164			for(i=0;i<nAtoms;i++){
165			comZ += atoms[i].pos[2]*atoms[i].mass;
166			totMass += atoms[i].mass;
167			}
168
169			comZ /= totMass;
170
171			halfMax = maxL * 0.5;
172
173			// calc the gofZ;
174
175			for(i=0;i<nAtoms;i++){
176
177			wrapZ = wrapMyZ( (atoms[i].pos[2] - comZ), Hmat[2][2] );
178
179			if(fabs(wrapZ) < halfMax ){
180
181			scale = wrapZ / halfMax;
182			bin = (int)( (scale*halfGzBins) + halfGzBins );
183
184			tempCorrs[atoms[i].type].gofz[bin] += 1.0;
185			}
186			}
187
188			volume = Hmat[0][0] * Hmat[1][1] * Hmat[2][2];
189			density = (double)nAtoms / volume;
190
191			densCons = Hmat[0][0] * Hmat[1][1] * gofzDl * density;
192
193			for(i=0;i<GOFZNBINS;i++)
194			for(j=0;j<5;j++)
195			tempCorrs[j].gofz[i] /= densCons;
196
197
198			for(i=0;i<GOFZNBINS;i++)
199			for(j=0;j<5;j++)
200			gofzCorrs[j].gofz[i] += tempCorrs[j].gofz[i];
201			}
202
203
204			double wrapMyZ( double z, double box ){
205
206			double scale, round;
207
208			scale = z / box;
209
210			round = ( scale >= 0 ) ? floor( scale + 0.5 ) : ceil( scale - 0.5 );
211
212			scale -= round;
213
214			return scale*box;
215			}