1 |
#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 |
} |