[ViewVC] Contents 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
#	Content
1	#!/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:])