[ViewVC] Annotation of: OpenMD/branches/heatflux/src/applications/utilities/stat2thcond

#!/usr/bin/env python
"""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
  -g, --green-kubo        use Green-Kubo formulae (noisy!)
  -e, --einstein          use Einstein relation (based on Hess 2002 paper)

The Green-Kubo formulae option will compute: V*<(S(t)*(S(0)>/kT^2,
which may be integrated to give a slowly-converging value for the viscosity.

The Einstein relation option will compute: V*<(\int_0^t S(t')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 -e -o ring5.pcorr

"""

__author__ = "Gianluca Puliti (gpuliti@nd.edu) and Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1667 $"
__date__ = "$Date: 2012-02-23 11:25:26 -0400 (Thu, 23 February 2012) $"

__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

    # Boltzmann's constant kcal/K
    kB = 3.29762e-27

    # converts (amu/fs^3)^2*fs^2 --> kcal^2/angstrom^4  ----> FOR HESS-EINSTEIN
    intSSdtdtConvert = 1.57288e-41

    # converts (amu/fs^3)^2*fs --> kcal^2/(fs *angstrom^4)   ----> FOR GREEN-KUBO
    intSSdtConvert = 1.57288e-41

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

    # Without unit conversions as follows it should be in Ang^3 / (kcal K)
    preV = volAve / (kB * tempAve * tempAve)
    
    if doGreenKubo:
        gkXcorr = []
        gkYcorr = []
        gkZcorr = []
        print "computing Green-Kubo-style Correlation Function"        
        # i corresponds to dt
        for i in range(len(time)):
            # j is the starting time for the correlation
            pp = 0.0
            
            ggX = 0.0
            ggY = 0.0
            ggZ = 0.0
            for j in range( len(time) - i ):
                
                ggX = ggX + Sx[j+i]*Sx[j]
                ggY = ggY + Sy[j+i]*Sy[j] 
                ggZ = ggZ + Sz[j+i]*Sz[j] 
           
            gkXcorr.append(ggX / float(len(time)-i))
            gkYcorr.append(ggY / float(len(time)-i))
            gkZcorr.append(ggZ / float(len(time)-i))

    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 doGreenKubo:
            outputFile.write("%f\t%13e\n" % (time[i], preV * intSSdtConvert * (gkXcorr[i]+gkYcorr[i]+gkZcorr[i])/3))
            
        if doEinstein:
            # outputFile.write("%f\t%13e\n" % (time[i], 0.5 * preV * heatfluxConvert * heatfluxConvert * dtConvert * dtConvert * (eXcorr[i] + eYcorr[i] + eZcorr[i])))
            outputFile.write("%f\t%13e\n" % (time[i], 0.5 * preV * intSSdtdtConvert * (eXcorr[i] + eYcorr[i] + eZcorr[i])/3))
                                              # Ang^3 / (kcal K)   *       kcal^2/angstrom^4
    outputFile.close()

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

    haveStatFileName = False
    haveOutputFileName = False
    doGreenKubo = 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 ("-g", "--green-kubo"):
            doGreenKubo = True
        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:	1678
Committed:	Thu Feb 23 18:20:53 2012 UTC (13 years, 9 months ago) by chuckv
File size:	9021 byte(s)
Log Message:	Added revision to svn for Gianluca.
#	User	Rev	Content
1	chuckv	1678	#!/usr/bin/env python
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			-g, --green-kubo use Green-Kubo formulae (noisy!)
15			-e, --einstein use Einstein relation (based on Hess 2002 paper)
16
17			The Green-Kubo formulae option will compute: V<(S(t)(S(0)>/kT^2,
18			which may be integrated to give a slowly-converging value for the viscosity.
19
20			The Einstein relation option will compute: V*<(\int_0^t S(t')dt')^2>/2kT^2,
21			which will grow approximately linearly in time. The long-time slope of this
22			function will be the viscosity.
23
24			Example:
25			stat2thcond -f ring5.stat -e -o ring5.pcorr
26
27			"""
28
29			__author__ = "Gianluca Puliti (gpuliti@nd.edu) and Dan Gezelter (gezelter@nd.edu)"
30			__version__ = "$Revision: 1667 $"
31			__date__ = "$Date: 2012-02-23 11:25:26 -0400 (Thu, 23 February 2012) $"
32
33			__copyright__ = "Copyright (c) 2007 by the University of Notre Dame"
34			__license__ = "OpenMD"
35
36			import sys
37			import getopt
38			import string
39			import math
40
41			def usage():
42			print __doc__
43
44			def readStatFile(statFileName):
45
46			global time
47			global temperature
48			global pressure
49			global volume
50			global Sx
51			global Sy
52			global Sz
53			time = []
54			temperature = []
55			pressure = []
56			volume = []
57			Sx = []
58			Sy = []
59			Sz = []
60
61			statFile = open(statFileName, 'r')
62			line = statFile.readline()
63
64			print "reading File"
65			pressSum = 0.0
66			volSum = 0.0
67			tempSum = 0.0
68			line = statFile.readline()
69			while 1:
70			L = line.split()
71			time.append(float(L[0]))
72			temperature.append(float(L[4]))
73			#
74			# OpenMD prints out pressure in units of atm.
75			#
76			pressure.append(float(L[5]))
77			volume.append(float(L[6]))
78			#
79			# OpenMD prints out heatflux in units of kcal / (mol s Ang^2).
80			#
81			Sx.append(float(L[8]))
82			Sy.append(float(L[9]))
83			Sz.append(float(L[10]))
84
85			line = statFile.readline()
86			if not line: break
87
88			statFile.close()
89
90			def computeAverages():
91
92			global tempAve
93			global pressAve
94			global volAve
95			global pvAve
96
97			print "computing Averages"
98
99			tempSum = 0.0
100			pressSum = 0.0
101			volSum = 0.0
102			pvSum = 0.0
103
104			temp2Sum = 0.0
105			press2Sum = 0.0
106			vol2Sum = 0.0
107			pv2Sum = 0.0
108
109			# converts amufs^-2Ang^-1 -> atm
110			pressureConvert = 1.63882576e8
111
112			for i in range(len(time)):
113			tempSum = tempSum + temperature[i]
114			pressSum = pressSum + pressure[i]
115			volSum = volSum + volume[i]
116			# in units of amu Ang^2 fs^-1
117			pvTerm = pressure[i]*volume[i] / pressureConvert
118			pvSum = pvSum + pvTerm
119			temp2Sum = temp2Sum + math.pow(temperature[i],2)
120			press2Sum = press2Sum + math.pow(pressure[i],2)
121			vol2Sum = vol2Sum + math.pow(volume[i],2)
122			pv2Sum = pv2Sum + math.pow(pvTerm,2)
123
124			tempAve = tempSum / float(len(time))
125			pressAve = pressSum / float(len(time))
126			volAve = volSum / float(len(time))
127			pvAve = pvSum / float(len(time))
128
129			tempSdev = math.sqrt(temp2Sum / float(len(time)) - math.pow(tempAve,2))
130			pressSdev = math.sqrt(press2Sum / float(len(time)) - math.pow(pressAve,2))
131			if (vol2Sum / float(len(time)) < math.pow(volAve,2)):
132			volSdev = 0.0
133			else:
134			volSdev = math.sqrt(vol2Sum / float(len(time)) - math.pow(volAve,2))
135			pvSdev = math.sqrt(pv2Sum / float(len(time)) - math.pow(pvAve,2))
136
137			print " Average pressure = %f +/- %f (atm)" % (pressAve, pressSdev)
138			print " Average volume = %f +/- %f (Angst^3)" % (volAve, volSdev)
139			print "Average temperature = %f +/- %f (K)" % (tempAve, tempSdev)
140			print " Average PV product = %f +/- %f (amu Angst^2 fs^-1)" % (pvAve, pvSdev)
141
142			def computeCorrelations(outputFileName):
143
144			# converts amufs^-2Ang^-1 -> atm
145			pressureConvert = 1.63882576e8
146
147			# converts Ang^-3 * kcal/mol * Ang / fs to m^-3 * J/mol * m /s (= W / mol m^2)
148			heatfluxConvert = 4.187e38
149
150			# converts fs to s
151			dtConvert = 1e-15
152
153			# Boltzmann's constant amuAng^2fs^-2/K
154			# kB = 8.31451e-7
155
156			# Boltzmann's constant kcal/K
157			kB = 3.29762e-27
158
159			# converts (amu/fs^3)^2*fs^2 --> kcal^2/angstrom^4 ----> FOR HESS-EINSTEIN
160			intSSdtdtConvert = 1.57288e-41
161
162			# converts (amu/fs^3)^2fs --> kcal^2/(fs angstrom^4) ----> FOR GREEN-KUBO
163			intSSdtConvert = 1.57288e-41
164
165			# preV = thcondConvert * volAve / (kB * tempAve * tempAve)
166
167			# Without unit conversions as follows it should be in Ang^3 / (kcal K)
168			preV = volAve / (kB * tempAve * tempAve)
169
170			if doGreenKubo:
171			gkXcorr = []
172			gkYcorr = []
173			gkZcorr = []
174			print "computing Green-Kubo-style Correlation Function"
175			# i corresponds to dt
176			for i in range(len(time)):
177			# j is the starting time for the correlation
178			pp = 0.0
179
180			ggX = 0.0
181			ggY = 0.0
182			ggZ = 0.0
183			for j in range( len(time) - i ):
184
185			ggX = ggX + Sx[j+i]*Sx[j]
186			ggY = ggY + Sy[j+i]*Sy[j]
187			ggZ = ggZ + Sz[j+i]*Sz[j]
188
189			gkXcorr.append(ggX / float(len(time)-i))
190			gkYcorr.append(ggY / float(len(time)-i))
191			gkZcorr.append(ggZ / float(len(time)-i))
192
193			if doEinstein:
194			print "computing Einstein-style Correlation Function"
195
196			# Precompute sum variables to aid integration.
197			# The integral from t0 -> t0 + t can be easily obtained
198			# from the precomputed sum variables: sum[t0+t] - sum[t0-1]
199			#for i in range(1, len(time)):
200			xSum = []
201			xSum.append(Sx[0])
202			ySum = []
203			ySum.append(Sy[0])
204			zSum = []
205			zSum.append(Sz[0])
206			for i in range(1, len(time)):
207			xSum.append(xSum[i-1] + Sx[i])
208			ySum.append(ySum[i-1] + Sy[i])
209			zSum.append(zSum[i-1] + Sz[i])
210
211			dt = time[1] - time[0]
212
213			eXcorr = []
214			eYcorr = []
215			eZcorr = []
216
217			# i corresponds to the total duration of the integral
218			for i in range(len(time)):
219
220			xIntSum = 0.0
221			yIntSum = 0.0
222			zIntSum = 0.0
223			# j corresponds to the starting point of the integral
224			for j in range(len(time) - i):
225			if (j == 0):
226
227			xInt = dt*xSum[j+i]
228			yInt = dt*ySum[j+i]
229			zInt = dt*zSum[j+i]
230			else:
231			xInt = dt*(xSum[j+i] - xSum[j-1])
232			yInt = dt*(ySum[j+i] - ySum[j-1])
233			zInt = dt*(zSum[j+i] - zSum[j-1])
234
235			xIntSum = xIntSum + xInt*xInt
236			yIntSum = yIntSum + yInt*yInt
237			zIntSum = zIntSum + zInt*zInt
238
239			eXcorr.append(xIntSum / float(len(time)-i))
240			eYcorr.append(yIntSum / float(len(time)-i))
241			eZcorr.append(zIntSum / float(len(time)-i))
242
243
244			outputFile = open(outputFileName, 'w')
245			for i in range(len(time)):
246			if doGreenKubo:
247			outputFile.write("%f\t%13e\n" % (time[i], preV * intSSdtConvert * (gkXcorr[i]+gkYcorr[i]+gkZcorr[i])/3))
248
249			if doEinstein:
250			# outputFile.write("%f\t%13e\n" % (time[i], 0.5 * preV * heatfluxConvert * heatfluxConvert * dtConvert * dtConvert * (eXcorr[i] + eYcorr[i] + eZcorr[i])))
251			outputFile.write("%f\t%13e\n" % (time[i], 0.5 * preV * intSSdtdtConvert * (eXcorr[i] + eYcorr[i] + eZcorr[i])/3))
252			# Ang^3 / (kcal K) * kcal^2/angstrom^4
253			outputFile.close()
254
255			def main(argv):
256			global doGreenKubo
257			global doEinstein
258			global haveStatFileName
259			global haveOutputFileName
260
261			haveStatFileName = False
262			haveOutputFileName = False
263			doGreenKubo = False
264			doEinstein = False
265
266			try:
267			opts, args = getopt.getopt(argv, "hgesf:o:", ["help", "einstein", "stat-file=", "output-file="])
268			except getopt.GetoptError:
269			usage()
270			sys.exit(2)
271			for opt, arg in opts:
272			if opt in ("-h", "--help"):
273			usage()
274			sys.exit()
275			elif opt in ("-g", "--green-kubo"):
276			doGreenKubo = True
277			elif opt in ("-e", "--einstein"):
278			doEinstein = True
279			elif opt in ("-f", "--stat-file"):
280			statFileName = arg
281			haveStatFileName = True
282			elif opt in ("-o", "--output-file"):
283			outputFileName = arg
284			haveOutputFileName = True
285			if (not haveStatFileName):
286			usage()
287			print "No stat file was specified"
288			sys.exit()
289			if (not haveOutputFileName):
290			usage()
291			print "No output file was specified"
292			sys.exit()
293
294			readStatFile(statFileName);
295			computeAverages();
296			computeCorrelations(outputFileName);
297
298			if __name__ == "__main__":
299			if len(sys.argv) == 1:
300			usage()
301			sys.exit()
302			main(sys.argv[1:])