[ViewVC] Annotation of: OpenMD/trunk/src/applications/utilities/stat2visco

#!/usr/bin/env python
"""Pressure Correlation function

Computes various correlation functions of the pressure and pressure tensor
that have been stored in a stat file.   These can be used to compute
shear and bulk viscosities. 

Usage: stat2visco 

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 (best)
  -s, --shear             compute the shear viscosity (the off-diagonal 
                          pressure tensor values must be present in the .stat
                          file)

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

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

Example:
   stat2visco -f ring5.stat -o ring5.pcorr -e -s

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision: 1.1 $"
__date__ = "$Date: 2008-01-23 21:22:18 $"

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

import sys
import getopt
import string
import math

def usage():
    print __doc__

def readStatFile(statFileName):

    global time
    global temperature
    global pressure
    global volume
    time = []
    temperature = []
    pressure = []
    volume = []

    if (doShear):
        global Pxx
        global Pyy
        global Pzz
        global Pxy
        global Pxz
        global Pyz
        
        Pxx = []
        Pyy = []
        Pzz = []
        Pxy = []
        Pxz = []
        Pyz = []

    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]))
        #
        # OOPSE prints out pressure in units of atm.
        #
        pressure.append(float(L[5]))
        volume.append(float(L[6]))

        if doShear:
            if (len(L) > 16):
                #
                # OOPSE prints out the pressure tensor in units of amu*fs^-2*Ang^-1
                #
                Pxx.append(float(L[8]))
                Pyy.append(float(L[12]))
                Pzz.append(float(L[16]))
                #
                # symmetrize the off-diagonal terms in the pressure tensor
                #
                Pxy.append(0.5*(float(L[9])  + float(L[11])))
                Pxz.append(0.5*(float(L[10]) + float(L[14])))
                Pyz.append(0.5*(float(L[13]) + float(L[15])))
            else:
                print "Not enough columns are present in the .stat file"
                print "to calculate the shear viscosity..."
                print 
                print "stat2visco expects to find all 9 elements of the"
                print "pressure tensor in columns 9-17 of the .stat file"
                print
                print "You may need to set the statFileFormat string"
                print "explicitly in your .md file when running OOPSE."
                print "Consult the OOPSE documentation for more details."
                sys.exit()
                
        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
    
    # Boltzmann's constant amu*Ang^2*fs^-2/K
    kB = 8.31451e-7

    # converts amu Ang^-1 fs^-1  ->  g cm^-1 s^-1
    viscoConvert = 0.16605387
    
    preV = viscoConvert * volAve / (kB * tempAve)
    preVi = viscoConvert / (volAve * kB * tempAve)
    
    if doGreenKubo:
        gkPcorr = []
        if doShear:
            gkXYcorr = []
            gkXZcorr = []
            gkYZcorr = []
        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
            if doShear:
                ppXY = 0.0
                ppXZ = 0.0
                ppYZ = 0.0
            for j in range( len(time) - i ):
                pv1 = pressure[j]*volume[j]/pressureConvert - pvAve
                pv2 = pressure[j+i]*volume[j+i]/pressureConvert - pvAve
                pp = pp + pv1*pv2
                if doShear:
                    ppXY = ppXY + Pxy[j+i]*Pxy[j]
                    ppXZ = ppXZ + Pxz[j+i]*Pxz[j] 
                    ppYZ = ppYZ + Pyz[j+i]*Pyz[j] 

            gkPcorr.append(pp / float(len(time) - i))
            if doShear:
                gkXYcorr.append(ppXY / float(len(time)-i))
                gkXZcorr.append(ppXZ / float(len(time)-i))
                gkYZcorr.append(ppYZ / 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]
        pSum = []
        pSum.append( (pressure[0] - pressAve) / pressureConvert)
        for i in range(1, len(time)):
            pSum.append(pSum[i-1] + (pressure[i]-pressAve)/pressureConvert )

        if doShear:
            xySum = []
            xySum.append(Pxy[0])
            xzSum = []
            xzSum.append(Pxz[0])
            yzSum = []
            yzSum.append(Pyz[0])
            for i in range(1, len(time)):
                xySum.append(xySum[i-1] + Pxy[i])
                xzSum.append(xzSum[i-1] + Pxz[i])
                yzSum.append(yzSum[i-1] + Pyz[i])


        ePcorr = []
        dt = time[1] - time[0]

        if doShear:
            eXYcorr = []
            eXZcorr = []
            eYZcorr = []

        # i corresponds to the total duration of the integral
        for i in range(len(time)):
            pIntSum = 0.0
            if doShear:
                xyIntSum = 0.0
                xzIntSum = 0.0
                yzIntSum = 0.0
            # j corresponds to the starting point of the integral
            for j in range(len(time) - i):
                if (j == 0):
                    pInt = dt*pSum[j+i]
                    if doShear:
                        xyInt = dt*xySum[j+i]
                        xzInt = dt*xzSum[j+i]
                        yzInt = dt*yzSum[j+i]
                else:
                    pInt = dt*(pSum[j+i] - pSum[j-1])
                    if doShear:
                        xyInt = dt*(xySum[j+i] - xySum[j-1])
                        xzInt = dt*(xzSum[j+i] - xzSum[j-1])
                        yzInt = dt*(yzSum[j+i] - yzSum[j-1])
                    
                pIntSum = pIntSum + pInt*pInt
                if doShear:
                    xyIntSum = xyIntSum + xyInt*xyInt
                    xzIntSum = xzIntSum + xzInt*xzInt
                    yzIntSum = yzIntSum + yzInt*yzInt
            ePcorr.append(pIntSum / float(len(time)-i))
            if doShear:
                eXYcorr.append(xyIntSum / float(len(time)-i))
                eXZcorr.append(xzIntSum / float(len(time)-i))
                eYZcorr.append(yzIntSum / float(len(time)-i))


    outputFile = open(outputFileName, 'w')
    for i in range(len(time)):
        if doGreenKubo:
            if doShear:
                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]))
            else:
                outputFile.write("%f\t%13e\n" % (time[i], preVi*gkPcorr[i]))
            
        if doEinstein:
            if doShear:
                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]))
            else:
                outputFile.write("%f\t%13e\n" % (time[i], 0.5*preV*ePcorr[i]))
    outputFile.close()

def main(argv):
    global doGreenKubo
    global doEinstein
    global doShear
    global haveStatFileName
    global haveOutputFileName
    
    haveStatFileName = False
    haveOutputFileName = False
    doShear = False
    doGreenKubo = False
    doEinstein = False
 
    try:                                
        opts, args = getopt.getopt(argv, "hgesf:o:", ["help", "green-kubo", "einstein", "shear", "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 ("-s", "--shear"):
            doShear = 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:	1214
Committed:	Wed Jan 23 21:22:18 2008 UTC (17 years, 9 months ago) by xsun
File size:	11456 byte(s)
Log Message:	added stat2visco
#	User	Rev	Content
1	xsun	1214	#!/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 amufs^-2Ang^-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 amufs^-2Ang^-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 amufs^-2Ang^-1 -> atm
177			pressureConvert = 1.63882576e8
178
179			# Boltzmann's constant amuAng^2fs^-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], preVigkPcorr[i], preVgkXYcorr[i], preVgkXZcorr[i], preVgkYZcorr[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.5preVePcorr[i], 0.5preVeXYcorr[i], 0.5preVeXZcorr[i], 0.5preVeYZcorr[i]))
296			else:
297			outputFile.write("%f\t%13e\n" % (time[i], 0.5preVePcorr[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:])