[ViewVC] Contents of: OpenMD/trunk/src/applications/utilities/stat2thcond

#!@PYTHON_EXECUTABLE@
"""Heat Flux Correlation function

Computes the correlation function of the heat flux vector
that has been stored in a stat file.   These can be used to compute
the thermal conductivity.

Usage: stat2thcond

Options:
  -h, --help              show this help
  -f, --stat-file=...     use specified stat file
  -o, --output-file=...   use specified output (.pcorr) file
  -e, --einstein          use Einstein relation (based on Hess 2002 paper)

The Einstein relation option will compute: V*<(\int_0^t (S(t')-<S>)dt')^2>/2kT^2,
which will grow approximately linearly in time.  The long-time slope of this
function will be the viscosity.

Example:
   stat2thcond -f ring5.stat -o ring5.pcorr

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1665 $"
__date__ = "$Date: 2011-12-08 15:25:26 -0400 (Thu, 9 December 2011) $"

__copyright__ = "Copyright (c) 2007 by the University of Notre Dame"
__license__ = "OpenMD"

import sys
import getopt
import string
import math

def usage():
    print __doc__

def readStatFile(statFileName):

    global time
    global temperature
    global pressure
    global volume
    global Sx
    global Sy
    global Sz
    time = []
    temperature = []
    pressure = []
    volume = []
    Sx = []
    Sy = []
    Sz = []

    statFile = open(statFileName, 'r')
    line = statFile.readline()

    print "reading File"
    pressSum = 0.0
    volSum = 0.0
    tempSum = 0.0
    line = statFile.readline()
    while 1:
        L = line.split()
        time.append(float(L[0]))
        temperature.append(float(L[4]))
        #
        # OpenMD prints out pressure in units of atm.
        #
        pressure.append(float(L[5]))
        volume.append(float(L[6]))
        #
        # OpenMD prints out heatflux in units of kcal / (mol s Ang^2).
        #
        Sx.append(float(L[8]))
        Sy.append(float(L[9]))
        Sz.append(float(L[10]))

        line = statFile.readline()
        if not line: break

    statFile.close()

def computeAverages():

    global tempAve
    global pressAve
    global volAve
    global pvAve

    print "computing Averages"

    tempSum = 0.0
    pressSum = 0.0
    volSum = 0.0
    pvSum = 0.0

    temp2Sum = 0.0
    press2Sum = 0.0
    vol2Sum = 0.0
    pv2Sum = 0.0

    # converts amu*fs^-2*Ang^-1 -> atm
    pressureConvert = 1.63882576e8

    for i in range(len(time)):
        tempSum = tempSum + temperature[i]
        pressSum = pressSum + pressure[i]
        volSum = volSum + volume[i]
        # in units of amu Ang^2 fs^-1
        pvTerm = pressure[i]*volume[i] / pressureConvert
        pvSum = pvSum + pvTerm
        temp2Sum = temp2Sum + math.pow(temperature[i],2)
        press2Sum = press2Sum + math.pow(pressure[i],2)
        vol2Sum = vol2Sum + math.pow(volume[i],2)
        pv2Sum = pv2Sum + math.pow(pvTerm,2)

    tempAve = tempSum / float(len(time))
    pressAve = pressSum / float(len(time))
    volAve = volSum / float(len(time))
    pvAve = pvSum / float(len(time))

    tempSdev = math.sqrt(temp2Sum / float(len(time)) - math.pow(tempAve,2))
    pressSdev = math.sqrt(press2Sum / float(len(time)) - math.pow(pressAve,2))
    if (vol2Sum / float(len(time)) < math.pow(volAve,2)):
        volSdev = 0.0
    else:
        volSdev = math.sqrt(vol2Sum / float(len(time)) - math.pow(volAve,2))
    pvSdev = math.sqrt(pv2Sum / float(len(time)) - math.pow(pvAve,2))

    print "   Average pressure = %f +/- %f (atm)" % (pressAve, pressSdev)
    print "     Average volume = %f +/- %f (Angst^3)" % (volAve, volSdev)
    print "Average temperature = %f +/- %f (K)" % (tempAve, tempSdev)
    print " Average PV product = %f +/- %f (amu Angst^2 fs^-1)" % (pvAve, pvSdev)

def computeCorrelations(outputFileName):

    # converts amu*fs^-2*Ang^-1 -> atm
    pressureConvert = 1.63882576e8

    # converts Ang^-3 * kcal/mol * Ang / fs to m^-3 * J/mol * m /s (= W / mol m^2)
    heatfluxConvert = 4.187e38

    # converts fs to s
    dtConvert = 1e-15

    # Boltzmann's constant amu*Ang^2*fs^-2/K
    kB = 8.31451e-7

    # converts Ang^3 / (amu*Ang^2*fs^-2/K * K^2)  ->  m s^2 / (kg K^2)
    thcondConvert = 6.0224e-14

    preV = thcondConvert * volAve / (kB * tempAve * tempAve)


    if doEinstein:
        print "computing Einstein-style Correlation Function"

        # Precompute sum variables to aid integration.
        # The integral from t0 -> t0 + t  can be easily obtained
        # from the precomputed sum variables:  sum[t0+t] - sum[t0-1]
        for i in range(1, len(time)):

        xSum = []
        xSum.append(Sx[0])
        ySum = []
        ySum.append(Sy[0])
        zSum = []
        zSum.append(Sz[0])
        for i in range(1, len(time)):
            xSum.append(xSum[i-1] + Sx[i])
            ySum.append(ySum[i-1] + Sy[i])
            zSum.append(zSum[i-1] + Sz[i])

        dt = time[1] - time[0]

        eXcorr = []
        eYcorr = []
        eZcorr = []

        # i corresponds to the total duration of the integral
        for i in range(len(time)):

            xIntSum = 0.0
            yIntSum = 0.0
            zIntSum = 0.0
            # j corresponds to the starting point of the integral
            for j in range(len(time) - i):
                if (j == 0):

                    xInt = dt*xSum[j+i]
                    yInt = dt*ySum[j+i]
                    zInt = dt*zSum[j+i]
                else:
                    xInt = dt*(xSum[j+i] - xSum[j-1])
                    yInt = dt*(ySum[j+i] - ySum[j-1])
                    zInt = dt*(zSum[j+i] - zSum[j-1])

                xIntSum = xIntSum + xInt*xInt
                yIntSum = yIntSum + yInt*yInt
                zIntSum = zIntSum + zInt*zInt

            eXcorr.append(xIntSum / float(len(time)-i))
            eYcorr.append(yIntSum / float(len(time)-i))
            eZcorr.append(zIntSum / float(len(time)-i))


    outputFile = open(outputFileName, 'w')
    for i in range(len(time)):

        if doEinstein:
            outputFile.write("%f\t%13e\n" % (time[i], 0.5 * preV * heatfluxConvert * heatfluxConvert * dtConvert * dtConvert * (eXcorr[i] + eYcorr[i] + eZcorr[i]))
    outputFile.close()

def main(argv):
    global doEinstein
    global haveStatFileName
    global haveOutputFileName

    haveStatFileName = False
    haveOutputFileName = False
    doEinstein = False

    try:
        opts, args = getopt.getopt(argv, "hgesf:o:", ["help", "einstein", "stat-file=", "output-file="])
    except getopt.GetoptError:
        usage()
        sys.exit(2)
    for opt, arg in opts:
        if opt in ("-h", "--help"):
            usage()
            sys.exit()
        elif opt in ("-e", "--einstein"):
            doEinstein = True
        elif opt in ("-f", "--stat-file"):
            statFileName = arg
            haveStatFileName = True
        elif opt in ("-o", "--output-file"):
            outputFileName = arg
            haveOutputFileName = True
    if (not haveStatFileName):
        usage()
        print "No stat file was specified"
        sys.exit()
    if (not haveOutputFileName):
        usage()
        print "No output file was specified"
        sys.exit()

    readStatFile(statFileName);
    computeAverages();
    computeCorrelations(outputFileName);

if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	1795
Committed:	Fri Sep 7 18:13:55 2012 UTC (12 years, 7 months ago) by gezelter
File size:	7336 byte(s)
Log Message:	Windows fixes
#	Content
1	#!@PYTHON_EXECUTABLE@
2	"""Heat Flux Correlation function
3
4	Computes the correlation function of the heat flux vector
5	that has been stored in a stat file. These can be used to compute
6	the thermal conductivity.
7
8	Usage: stat2thcond
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	-e, --einstein use Einstein relation (based on Hess 2002 paper)
15
16	The Einstein relation option will compute: V*<(\int_0^t (S(t')-<S>)dt')^2>/2kT^2,
17	which will grow approximately linearly in time. The long-time slope of this
18	function will be the viscosity.
19
20	Example:
21	stat2thcond -f ring5.stat -o ring5.pcorr
22
23	"""
24
25	__author__ = "Dan Gezelter (gezelter@nd.edu)"
26	__version__ = "$Revision: 1665 $"
27	__date__ = "$Date: 2011-12-08 15:25:26 -0400 (Thu, 9 December 2011) $"
28
29	__copyright__ = "Copyright (c) 2007 by the University of Notre Dame"
30	__license__ = "OpenMD"
31
32	import sys
33	import getopt
34	import string
35	import math
36
37	def usage():
38	print __doc__
39
40	def readStatFile(statFileName):
41
42	global time
43	global temperature
44	global pressure
45	global volume
46	global Sx
47	global Sy
48	global Sz
49	time = []
50	temperature = []
51	pressure = []
52	volume = []
53	Sx = []
54	Sy = []
55	Sz = []
56
57	statFile = open(statFileName, 'r')
58	line = statFile.readline()
59
60	print "reading File"
61	pressSum = 0.0
62	volSum = 0.0
63	tempSum = 0.0
64	line = statFile.readline()
65	while 1:
66	L = line.split()
67	time.append(float(L[0]))
68	temperature.append(float(L[4]))
69	#
70	# OpenMD prints out pressure in units of atm.
71	#
72	pressure.append(float(L[5]))
73	volume.append(float(L[6]))
74	#
75	# OpenMD prints out heatflux in units of kcal / (mol s Ang^2).
76	#
77	Sx.append(float(L[8]))
78	Sy.append(float(L[9]))
79	Sz.append(float(L[10]))
80
81	line = statFile.readline()
82	if not line: break
83
84	statFile.close()
85
86	def computeAverages():
87
88	global tempAve
89	global pressAve
90	global volAve
91	global pvAve
92
93	print "computing Averages"
94
95	tempSum = 0.0
96	pressSum = 0.0
97	volSum = 0.0
98	pvSum = 0.0
99
100	temp2Sum = 0.0
101	press2Sum = 0.0
102	vol2Sum = 0.0
103	pv2Sum = 0.0
104
105	# converts amufs^-2Ang^-1 -> atm
106	pressureConvert = 1.63882576e8
107
108	for i in range(len(time)):
109	tempSum = tempSum + temperature[i]
110	pressSum = pressSum + pressure[i]
111	volSum = volSum + volume[i]
112	# in units of amu Ang^2 fs^-1
113	pvTerm = pressure[i]*volume[i] / pressureConvert
114	pvSum = pvSum + pvTerm
115	temp2Sum = temp2Sum + math.pow(temperature[i],2)
116	press2Sum = press2Sum + math.pow(pressure[i],2)
117	vol2Sum = vol2Sum + math.pow(volume[i],2)
118	pv2Sum = pv2Sum + math.pow(pvTerm,2)
119
120	tempAve = tempSum / float(len(time))
121	pressAve = pressSum / float(len(time))
122	volAve = volSum / float(len(time))
123	pvAve = pvSum / float(len(time))
124
125	tempSdev = math.sqrt(temp2Sum / float(len(time)) - math.pow(tempAve,2))
126	pressSdev = math.sqrt(press2Sum / float(len(time)) - math.pow(pressAve,2))
127	if (vol2Sum / float(len(time)) < math.pow(volAve,2)):
128	volSdev = 0.0
129	else:
130	volSdev = math.sqrt(vol2Sum / float(len(time)) - math.pow(volAve,2))
131	pvSdev = math.sqrt(pv2Sum / float(len(time)) - math.pow(pvAve,2))
132
133	print " Average pressure = %f +/- %f (atm)" % (pressAve, pressSdev)
134	print " Average volume = %f +/- %f (Angst^3)" % (volAve, volSdev)
135	print "Average temperature = %f +/- %f (K)" % (tempAve, tempSdev)
136	print " Average PV product = %f +/- %f (amu Angst^2 fs^-1)" % (pvAve, pvSdev)
137
138	def computeCorrelations(outputFileName):
139
140	# converts amufs^-2Ang^-1 -> atm
141	pressureConvert = 1.63882576e8
142
143	# converts Ang^-3 * kcal/mol * Ang / fs to m^-3 * J/mol * m /s (= W / mol m^2)
144	heatfluxConvert = 4.187e38
145
146	# converts fs to s
147	dtConvert = 1e-15
148
149	# Boltzmann's constant amuAng^2fs^-2/K
150	kB = 8.31451e-7
151
152	# converts Ang^3 / (amuAng^2fs^-2/K * K^2) -> m s^2 / (kg K^2)
153	thcondConvert = 6.0224e-14
154
155	preV = thcondConvert * volAve / (kB * tempAve * tempAve)
156
157
158	if doEinstein:
159	print "computing Einstein-style Correlation Function"
160
161	# Precompute sum variables to aid integration.
162	# The integral from t0 -> t0 + t can be easily obtained
163	# from the precomputed sum variables: sum[t0+t] - sum[t0-1]
164	for i in range(1, len(time)):
165
166	xSum = []
167	xSum.append(Sx[0])
168	ySum = []
169	ySum.append(Sy[0])
170	zSum = []
171	zSum.append(Sz[0])
172	for i in range(1, len(time)):
173	xSum.append(xSum[i-1] + Sx[i])
174	ySum.append(ySum[i-1] + Sy[i])
175	zSum.append(zSum[i-1] + Sz[i])
176
177	dt = time[1] - time[0]
178
179	eXcorr = []
180	eYcorr = []
181	eZcorr = []
182
183	# i corresponds to the total duration of the integral
184	for i in range(len(time)):
185
186	xIntSum = 0.0
187	yIntSum = 0.0
188	zIntSum = 0.0
189	# j corresponds to the starting point of the integral
190	for j in range(len(time) - i):
191	if (j == 0):
192
193	xInt = dt*xSum[j+i]
194	yInt = dt*ySum[j+i]
195	zInt = dt*zSum[j+i]
196	else:
197	xInt = dt*(xSum[j+i] - xSum[j-1])
198	yInt = dt*(ySum[j+i] - ySum[j-1])
199	zInt = dt*(zSum[j+i] - zSum[j-1])
200
201	xIntSum = xIntSum + xInt*xInt
202	yIntSum = yIntSum + yInt*yInt
203	zIntSum = zIntSum + zInt*zInt
204
205	eXcorr.append(xIntSum / float(len(time)-i))
206	eYcorr.append(yIntSum / float(len(time)-i))
207	eZcorr.append(zIntSum / float(len(time)-i))
208
209
210	outputFile = open(outputFileName, 'w')
211	for i in range(len(time)):
212
213	if doEinstein:
214	outputFile.write("%f\t%13e\n" % (time[i], 0.5 * preV * heatfluxConvert * heatfluxConvert * dtConvert * dtConvert * (eXcorr[i] + eYcorr[i] + eZcorr[i]))
215	outputFile.close()
216
217	def main(argv):
218	global doEinstein
219	global haveStatFileName
220	global haveOutputFileName
221
222	haveStatFileName = False
223	haveOutputFileName = False
224	doEinstein = False
225
226	try:
227	opts, args = getopt.getopt(argv, "hgesf:o:", ["help", "einstein", "stat-file=", "output-file="])
228	except getopt.GetoptError:
229	usage()
230	sys.exit(2)
231	for opt, arg in opts:
232	if opt in ("-h", "--help"):
233	usage()
234	sys.exit()
235	elif opt in ("-e", "--einstein"):
236	doEinstein = True
237	elif opt in ("-f", "--stat-file"):
238	statFileName = arg
239	haveStatFileName = True
240	elif opt in ("-o", "--output-file"):
241	outputFileName = arg
242	haveOutputFileName = True
243	if (not haveStatFileName):
244	usage()
245	print "No stat file was specified"
246	sys.exit()
247	if (not haveOutputFileName):
248	usage()
249	print "No output file was specified"
250	sys.exit()
251
252	readStatFile(statFileName);
253	computeAverages();
254	computeCorrelations(outputFileName);
255
256	if __name__ == "__main__":
257	if len(sys.argv) == 1:
258	usage()
259	sys.exit()
260	main(sys.argv[1:])