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