1 |
xsun |
3323 |
#!/usr/bin/env python |
2 |
|
|
"""Pressure Correlation function |
3 |
|
|
|
4 |
|
|
Computes various correlation functions of the pressure and pressure tensor |
5 |
|
|
that have been stored in a stat file. These can be used to compute |
6 |
|
|
shear and bulk viscosities. |
7 |
|
|
|
8 |
|
|
Usage: stat2visco |
9 |
|
|
|
10 |
|
|
Options: |
11 |
|
|
-h, --help show this help |
12 |
|
|
-f, --stat-file=... use specified stat file |
13 |
|
|
-o, --output-file=... use specified output (.pcorr) file |
14 |
|
|
-g, --green-kubo use Green-Kubo formulae (noisy!) |
15 |
|
|
-e, --einstein use Einstein relation (best) |
16 |
|
|
-s, --shear compute the shear viscosity (the off-diagonal |
17 |
|
|
pressure tensor values must be present in the .stat |
18 |
|
|
file) |
19 |
|
|
|
20 |
|
|
The Green-Kubo formulae option will compute: V*<(P(t)-<P>)*(P(0)-<P>)>/kT , |
21 |
|
|
which may be integrated to give a slowly-converging value for the viscosity. |
22 |
|
|
|
23 |
|
|
The Einstein relation option will compute: V*<(\int_0^t (P(t')-<P>)dt')^2>/2kT, |
24 |
|
|
which will grow approximately linearly in time. The long-time slope of this |
25 |
|
|
function will be the viscosity. |
26 |
|
|
|
27 |
|
|
Example: |
28 |
|
|
stat2visco -f ring5.stat -o ring5.pcorr -e -s |
29 |
|
|
|
30 |
|
|
""" |
31 |
|
|
|
32 |
|
|
__author__ = "Dan Gezelter (gezelter@nd.edu)" |
33 |
|
|
__version__ = "$Revision: 1.1 $" |
34 |
|
|
__date__ = "$Date: 2008-01-23 21:22:18 $" |
35 |
|
|
|
36 |
|
|
__copyright__ = "Copyright (c) 2007 by the University of Notre Dame" |
37 |
|
|
__license__ = "OOPSE" |
38 |
|
|
|
39 |
|
|
import sys |
40 |
|
|
import getopt |
41 |
|
|
import string |
42 |
|
|
import math |
43 |
|
|
|
44 |
|
|
def usage(): |
45 |
|
|
print __doc__ |
46 |
|
|
|
47 |
|
|
def readStatFile(statFileName): |
48 |
|
|
|
49 |
|
|
global time |
50 |
|
|
global temperature |
51 |
|
|
global pressure |
52 |
|
|
global volume |
53 |
|
|
time = [] |
54 |
|
|
temperature = [] |
55 |
|
|
pressure = [] |
56 |
|
|
volume = [] |
57 |
|
|
|
58 |
|
|
if (doShear): |
59 |
|
|
global Pxx |
60 |
|
|
global Pyy |
61 |
|
|
global Pzz |
62 |
|
|
global Pxy |
63 |
|
|
global Pxz |
64 |
|
|
global Pyz |
65 |
|
|
|
66 |
|
|
Pxx = [] |
67 |
|
|
Pyy = [] |
68 |
|
|
Pzz = [] |
69 |
|
|
Pxy = [] |
70 |
|
|
Pxz = [] |
71 |
|
|
Pyz = [] |
72 |
|
|
|
73 |
|
|
statFile = open(statFileName, 'r') |
74 |
|
|
line = statFile.readline() |
75 |
|
|
|
76 |
|
|
print "reading File" |
77 |
|
|
pressSum = 0.0 |
78 |
|
|
volSum = 0.0 |
79 |
|
|
tempSum = 0.0 |
80 |
|
|
line = statFile.readline() |
81 |
|
|
while 1: |
82 |
|
|
L = line.split() |
83 |
|
|
time.append(float(L[0])) |
84 |
|
|
temperature.append(float(L[4])) |
85 |
|
|
# |
86 |
|
|
# OOPSE prints out pressure in units of atm. |
87 |
|
|
# |
88 |
|
|
pressure.append(float(L[5])) |
89 |
|
|
volume.append(float(L[6])) |
90 |
|
|
|
91 |
|
|
if doShear: |
92 |
|
|
if (len(L) > 16): |
93 |
|
|
# |
94 |
|
|
# OOPSE prints out the pressure tensor in units of amu*fs^-2*Ang^-1 |
95 |
|
|
# |
96 |
|
|
Pxx.append(float(L[8])) |
97 |
|
|
Pyy.append(float(L[12])) |
98 |
|
|
Pzz.append(float(L[16])) |
99 |
|
|
# |
100 |
|
|
# symmetrize the off-diagonal terms in the pressure tensor |
101 |
|
|
# |
102 |
|
|
Pxy.append(0.5*(float(L[9]) + float(L[11]))) |
103 |
|
|
Pxz.append(0.5*(float(L[10]) + float(L[14]))) |
104 |
|
|
Pyz.append(0.5*(float(L[13]) + float(L[15]))) |
105 |
|
|
else: |
106 |
|
|
print "Not enough columns are present in the .stat file" |
107 |
|
|
print "to calculate the shear viscosity..." |
108 |
|
|
print |
109 |
|
|
print "stat2visco expects to find all 9 elements of the" |
110 |
|
|
print "pressure tensor in columns 9-17 of the .stat file" |
111 |
|
|
print |
112 |
|
|
print "You may need to set the statFileFormat string" |
113 |
|
|
print "explicitly in your .md file when running OOPSE." |
114 |
|
|
print "Consult the OOPSE documentation for more details." |
115 |
|
|
sys.exit() |
116 |
|
|
|
117 |
|
|
line = statFile.readline() |
118 |
|
|
if not line: break |
119 |
|
|
|
120 |
|
|
statFile.close() |
121 |
|
|
|
122 |
|
|
def computeAverages(): |
123 |
|
|
|
124 |
|
|
global tempAve |
125 |
|
|
global pressAve |
126 |
|
|
global volAve |
127 |
|
|
global pvAve |
128 |
|
|
|
129 |
|
|
print "computing Averages" |
130 |
|
|
|
131 |
|
|
tempSum = 0.0 |
132 |
|
|
pressSum = 0.0 |
133 |
|
|
volSum = 0.0 |
134 |
|
|
pvSum = 0.0 |
135 |
|
|
|
136 |
|
|
temp2Sum = 0.0 |
137 |
|
|
press2Sum = 0.0 |
138 |
|
|
vol2Sum = 0.0 |
139 |
|
|
pv2Sum = 0.0 |
140 |
|
|
|
141 |
|
|
# converts amu*fs^-2*Ang^-1 -> atm |
142 |
|
|
pressureConvert = 1.63882576e8 |
143 |
|
|
|
144 |
|
|
for i in range(len(time)): |
145 |
|
|
tempSum = tempSum + temperature[i] |
146 |
|
|
pressSum = pressSum + pressure[i] |
147 |
|
|
volSum = volSum + volume[i] |
148 |
|
|
# in units of amu Ang^2 fs^-1 |
149 |
|
|
pvTerm = pressure[i]*volume[i] / pressureConvert |
150 |
|
|
pvSum = pvSum + pvTerm |
151 |
|
|
temp2Sum = temp2Sum + math.pow(temperature[i],2) |
152 |
|
|
press2Sum = press2Sum + math.pow(pressure[i],2) |
153 |
|
|
vol2Sum = vol2Sum + math.pow(volume[i],2) |
154 |
|
|
pv2Sum = pv2Sum + math.pow(pvTerm,2) |
155 |
|
|
|
156 |
|
|
tempAve = tempSum / float(len(time)) |
157 |
|
|
pressAve = pressSum / float(len(time)) |
158 |
|
|
volAve = volSum / float(len(time)) |
159 |
|
|
pvAve = pvSum / float(len(time)) |
160 |
|
|
|
161 |
|
|
tempSdev = math.sqrt(temp2Sum / float(len(time)) - math.pow(tempAve,2)) |
162 |
|
|
pressSdev = math.sqrt(press2Sum / float(len(time)) - math.pow(pressAve,2)) |
163 |
|
|
if (vol2Sum / float(len(time)) < math.pow(volAve,2)): |
164 |
|
|
volSdev = 0.0 |
165 |
|
|
else: |
166 |
|
|
volSdev = math.sqrt(vol2Sum / float(len(time)) - math.pow(volAve,2)) |
167 |
|
|
pvSdev = math.sqrt(pv2Sum / float(len(time)) - math.pow(pvAve,2)) |
168 |
|
|
|
169 |
|
|
print " Average pressure = %f +/- %f (atm)" % (pressAve, pressSdev) |
170 |
|
|
print " Average volume = %f +/- %f (Angst^3)" % (volAve, volSdev) |
171 |
|
|
print "Average temperature = %f +/- %f (K)" % (tempAve, tempSdev) |
172 |
|
|
print " Average PV product = %f +/- %f (amu Angst^2 fs^-1)" % (pvAve, pvSdev) |
173 |
|
|
|
174 |
|
|
def computeCorrelations(outputFileName): |
175 |
|
|
|
176 |
|
|
# converts amu*fs^-2*Ang^-1 -> atm |
177 |
|
|
pressureConvert = 1.63882576e8 |
178 |
|
|
|
179 |
|
|
# Boltzmann's constant amu*Ang^2*fs^-2/K |
180 |
|
|
kB = 8.31451e-7 |
181 |
|
|
|
182 |
|
|
# converts amu Ang^-1 fs^-1 -> g cm^-1 s^-1 |
183 |
|
|
viscoConvert = 0.16605387 |
184 |
|
|
|
185 |
|
|
preV = viscoConvert * volAve / (kB * tempAve) |
186 |
|
|
preVi = viscoConvert / (volAve * kB * tempAve) |
187 |
|
|
|
188 |
|
|
if doGreenKubo: |
189 |
|
|
gkPcorr = [] |
190 |
|
|
if doShear: |
191 |
|
|
gkXYcorr = [] |
192 |
|
|
gkXZcorr = [] |
193 |
|
|
gkYZcorr = [] |
194 |
|
|
print "computing Green-Kubo-style Correlation Function" |
195 |
|
|
# i corresponds to dt |
196 |
|
|
for i in range(len(time)): |
197 |
|
|
# j is the starting time for the correlation |
198 |
|
|
pp = 0.0 |
199 |
|
|
if doShear: |
200 |
|
|
ppXY = 0.0 |
201 |
|
|
ppXZ = 0.0 |
202 |
|
|
ppYZ = 0.0 |
203 |
|
|
for j in range( len(time) - i ): |
204 |
|
|
pv1 = pressure[j]*volume[j]/pressureConvert - pvAve |
205 |
|
|
pv2 = pressure[j+i]*volume[j+i]/pressureConvert - pvAve |
206 |
|
|
pp = pp + pv1*pv2 |
207 |
|
|
if doShear: |
208 |
|
|
ppXY = ppXY + Pxy[j+i]*Pxy[j] |
209 |
|
|
ppXZ = ppXZ + Pxz[j+i]*Pxz[j] |
210 |
|
|
ppYZ = ppYZ + Pyz[j+i]*Pyz[j] |
211 |
|
|
|
212 |
|
|
gkPcorr.append(pp / float(len(time) - i)) |
213 |
|
|
if doShear: |
214 |
|
|
gkXYcorr.append(ppXY / float(len(time)-i)) |
215 |
|
|
gkXZcorr.append(ppXZ / float(len(time)-i)) |
216 |
|
|
gkYZcorr.append(ppYZ / float(len(time)-i)) |
217 |
|
|
|
218 |
|
|
|
219 |
|
|
if doEinstein: |
220 |
|
|
print "computing Einstein-style Correlation Function" |
221 |
|
|
|
222 |
|
|
# Precompute sum variables to aid integration. |
223 |
|
|
# The integral from t0 -> t0 + t can be easily obtained |
224 |
|
|
# from the precomputed sum variables: sum[t0+t] - sum[t0-1] |
225 |
|
|
pSum = [] |
226 |
|
|
pSum.append( (pressure[0] - pressAve) / pressureConvert) |
227 |
|
|
for i in range(1, len(time)): |
228 |
|
|
pSum.append(pSum[i-1] + (pressure[i]-pressAve)/pressureConvert ) |
229 |
|
|
|
230 |
|
|
if doShear: |
231 |
|
|
xySum = [] |
232 |
|
|
xySum.append(Pxy[0]) |
233 |
|
|
xzSum = [] |
234 |
|
|
xzSum.append(Pxz[0]) |
235 |
|
|
yzSum = [] |
236 |
|
|
yzSum.append(Pyz[0]) |
237 |
|
|
for i in range(1, len(time)): |
238 |
|
|
xySum.append(xySum[i-1] + Pxy[i]) |
239 |
|
|
xzSum.append(xzSum[i-1] + Pxz[i]) |
240 |
|
|
yzSum.append(yzSum[i-1] + Pyz[i]) |
241 |
|
|
|
242 |
|
|
|
243 |
|
|
ePcorr = [] |
244 |
|
|
dt = time[1] - time[0] |
245 |
|
|
|
246 |
|
|
if doShear: |
247 |
|
|
eXYcorr = [] |
248 |
|
|
eXZcorr = [] |
249 |
|
|
eYZcorr = [] |
250 |
|
|
|
251 |
|
|
# i corresponds to the total duration of the integral |
252 |
|
|
for i in range(len(time)): |
253 |
|
|
pIntSum = 0.0 |
254 |
|
|
if doShear: |
255 |
|
|
xyIntSum = 0.0 |
256 |
|
|
xzIntSum = 0.0 |
257 |
|
|
yzIntSum = 0.0 |
258 |
|
|
# j corresponds to the starting point of the integral |
259 |
|
|
for j in range(len(time) - i): |
260 |
|
|
if (j == 0): |
261 |
|
|
pInt = dt*pSum[j+i] |
262 |
|
|
if doShear: |
263 |
|
|
xyInt = dt*xySum[j+i] |
264 |
|
|
xzInt = dt*xzSum[j+i] |
265 |
|
|
yzInt = dt*yzSum[j+i] |
266 |
|
|
else: |
267 |
|
|
pInt = dt*(pSum[j+i] - pSum[j-1]) |
268 |
|
|
if doShear: |
269 |
|
|
xyInt = dt*(xySum[j+i] - xySum[j-1]) |
270 |
|
|
xzInt = dt*(xzSum[j+i] - xzSum[j-1]) |
271 |
|
|
yzInt = dt*(yzSum[j+i] - yzSum[j-1]) |
272 |
|
|
|
273 |
|
|
pIntSum = pIntSum + pInt*pInt |
274 |
|
|
if doShear: |
275 |
|
|
xyIntSum = xyIntSum + xyInt*xyInt |
276 |
|
|
xzIntSum = xzIntSum + xzInt*xzInt |
277 |
|
|
yzIntSum = yzIntSum + yzInt*yzInt |
278 |
|
|
ePcorr.append(pIntSum / float(len(time)-i)) |
279 |
|
|
if doShear: |
280 |
|
|
eXYcorr.append(xyIntSum / float(len(time)-i)) |
281 |
|
|
eXZcorr.append(xzIntSum / float(len(time)-i)) |
282 |
|
|
eYZcorr.append(yzIntSum / float(len(time)-i)) |
283 |
|
|
|
284 |
|
|
|
285 |
|
|
outputFile = open(outputFileName, 'w') |
286 |
|
|
for i in range(len(time)): |
287 |
|
|
if doGreenKubo: |
288 |
|
|
if doShear: |
289 |
|
|
outputFile.write("%f\t%13e\t%13e\t%13e\t%13e\n" % (time[i], preVi*gkPcorr[i], preV*gkXYcorr[i], preV*gkXZcorr[i], preV*gkYZcorr[i])) |
290 |
|
|
else: |
291 |
|
|
outputFile.write("%f\t%13e\n" % (time[i], preVi*gkPcorr[i])) |
292 |
|
|
|
293 |
|
|
if doEinstein: |
294 |
|
|
if doShear: |
295 |
|
|
outputFile.write("%f\t%13e\t%13e\t%13e\t%13e\n" % (time[i], 0.5*preV*ePcorr[i], 0.5*preV*eXYcorr[i], 0.5*preV*eXZcorr[i], 0.5*preV*eYZcorr[i])) |
296 |
|
|
else: |
297 |
|
|
outputFile.write("%f\t%13e\n" % (time[i], 0.5*preV*ePcorr[i])) |
298 |
|
|
outputFile.close() |
299 |
|
|
|
300 |
|
|
def main(argv): |
301 |
|
|
global doGreenKubo |
302 |
|
|
global doEinstein |
303 |
|
|
global doShear |
304 |
|
|
global haveStatFileName |
305 |
|
|
global haveOutputFileName |
306 |
|
|
|
307 |
|
|
haveStatFileName = False |
308 |
|
|
haveOutputFileName = False |
309 |
|
|
doShear = False |
310 |
|
|
doGreenKubo = False |
311 |
|
|
doEinstein = False |
312 |
|
|
|
313 |
|
|
try: |
314 |
|
|
opts, args = getopt.getopt(argv, "hgesf:o:", ["help", "green-kubo", "einstein", "shear", "stat-file=", "output-file="]) |
315 |
|
|
except getopt.GetoptError: |
316 |
|
|
usage() |
317 |
|
|
sys.exit(2) |
318 |
|
|
for opt, arg in opts: |
319 |
|
|
if opt in ("-h", "--help"): |
320 |
|
|
usage() |
321 |
|
|
sys.exit() |
322 |
|
|
elif opt in ("-g", "--green-kubo"): |
323 |
|
|
doGreenKubo = True |
324 |
|
|
elif opt in ("-e", "--einstein"): |
325 |
|
|
doEinstein = True |
326 |
|
|
elif opt in ("-s", "--shear"): |
327 |
|
|
doShear = True |
328 |
|
|
elif opt in ("-f", "--stat-file"): |
329 |
|
|
statFileName = arg |
330 |
|
|
haveStatFileName = True |
331 |
|
|
elif opt in ("-o", "--output-file"): |
332 |
|
|
outputFileName = arg |
333 |
|
|
haveOutputFileName = True |
334 |
|
|
if (not haveStatFileName): |
335 |
|
|
usage() |
336 |
|
|
print "No stat file was specified" |
337 |
|
|
sys.exit() |
338 |
|
|
if (not haveOutputFileName): |
339 |
|
|
usage() |
340 |
|
|
print "No output file was specified" |
341 |
|
|
sys.exit() |
342 |
|
|
|
343 |
|
|
readStatFile(statFileName); |
344 |
|
|
computeAverages(); |
345 |
|
|
computeCorrelations(outputFileName); |
346 |
|
|
|
347 |
|
|
if __name__ == "__main__": |
348 |
|
|
if len(sys.argv) == 1: |
349 |
|
|
usage() |
350 |
|
|
sys.exit() |
351 |
|
|
main(sys.argv[1:]) |