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

#!@PYTHON_EXECUTABLE@
"""A script that processes a velocity autocorrelation function into a
amplitude spectrum

Post-processes a vcorr file and creates a normalized spectrum.

Usage: vcorr2spectrum 

Options:
  -h, --help              show this help
  -f, --vcorr-file=...    use specified vcorr file
  -o, --output-file=...   use specified output (.pspect) file

Example:
   vcorr2spectrum -f stockmayer.vcorr -o stockmayer.pspect

"""

__author__ = "Dan Gezelter (gezelter@nd.edu)"
__version__ = "$Revision:  $"
__date__ = "$Date:  $"

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

import sys
import getopt
import string
import math
import numpy as np

def usage():
    print __doc__

def readVcorrFile(vcorrFileName):
    global time
    global vcorr
    time = []
    vcorr = []

    vcorrFile = open(vcorrFileName, 'r')
    line = vcorrFile.readline()

    print "reading File"
    line = vcorrFile.readline()
    while 1:

        if not line.startswith("#"):
            L = line.split()
        
            time.append(float(L[0]))
            vcorr.append(float(L[1]))
            
        line = vcorrFile.readline()
        if not line: break
        
    vcorrFile.close()
    
def computeSpectrum(outFileName):
    global tnew
    global vnew
    tnew = []
    vnew = []

    # the speed of light in cm / fs
    c = 2.99792458e-5
    
    # symmetrize arrays

    tlen = len(time)
    for i in range(1,len(time)):
        tnew.append( -time[tlen-i] )
        vnew.append( vcorr[tlen-i] )
    for i in range(len(time)):
        tnew.append( time[i] )
        vnew.append( vcorr[i] )
   
    spect = np.fft.fft( vnew )  
    n = len(vnew)
    timestep = tnew[1] - tnew[0]
    freq = np.fft.fftfreq(n, d=timestep)  

    # freq = np.fft.fftshift(freq)
    # y = np.fft.fftshift(spect)

    outFile = open(outFileName, 'w')

    freqScale = 1.0 / c
    dfreq = (freq[1]-freq[0]) * freqScale
    s = 0.0
    
    for i in range(n/2):
        s = s + np.abs(spect[i])*dfreq
    
    for i in range(n/2):
         outFile.write("%lf\t%lf\n" % (freq[i] * freqScale, np.abs(spect[i])/s))

    outFile.close()
    
def main(argv):
    global haveVcorrFileName
    global haveOutputFileName
    
    haveVcorrFileName = False
    haveOutputFileName = False
    haveQValue = False
 
    try:                                
        opts, args = getopt.getopt(argv, "hf:o:", ["help", "Q-value=", "vcorr-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 ("-f", "--vcorr-file"): 
            vcorrFileName = arg
            haveVcorrFileName = True
        elif opt in ("-o", "--output-file"): 
            outputFileName = arg
            haveOutputFileName = True
    if (not haveVcorrFileName):
        usage() 
        print "No vcorr file was specified"
        sys.exit()
    if (not haveOutputFileName):
        usage()
        print "No output file was specified"
        sys.exit()
    
    
    readVcorrFile(vcorrFileName)
    computeSpectrum(outputFileName)


if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	2051
Committed:	Thu Jan 8 16:48:14 2015 UTC (11 years ago) by gezelter
File size:	3399 byte(s)
Log Message:	Fixed a frequency computation bug in a utility script
#	User	Rev	Content
1	gezelter	2045	#!@PYTHON_EXECUTABLE@
2			"""A script that processes a velocity autocorrelation function into a
3			amplitude spectrum
4
5			Post-processes a vcorr file and creates a normalized spectrum.
6
7			Usage: vcorr2spectrum
8
9			Options:
10			-h, --help show this help
11			-f, --vcorr-file=... use specified vcorr file
12			-o, --output-file=... use specified output (.pspect) file
13
14			Example:
15			vcorr2spectrum -f stockmayer.vcorr -o stockmayer.pspect
16
17			"""
18
19			__author__ = "Dan Gezelter (gezelter@nd.edu)"
20			__version__ = "$Revision: $"
21			__date__ = "$Date: $"
22
23			__copyright__ = "Copyright (c) 2014 by the University of Notre Dame"
24			__license__ = "OpenMD"
25
26			import sys
27			import getopt
28			import string
29			import math
30			import numpy as np
31
32			def usage():
33			print __doc__
34
35			def readVcorrFile(vcorrFileName):
36			global time
37			global vcorr
38			time = []
39			vcorr = []
40
41			vcorrFile = open(vcorrFileName, 'r')
42			line = vcorrFile.readline()
43
44			print "reading File"
45			line = vcorrFile.readline()
46			while 1:
47
48			if not line.startswith("#"):
49			L = line.split()
50
51			time.append(float(L[0]))
52			vcorr.append(float(L[1]))
53
54			line = vcorrFile.readline()
55			if not line: break
56
57			vcorrFile.close()
58
59			def computeSpectrum(outFileName):
60			global tnew
61			global vnew
62			tnew = []
63			vnew = []
64
65			# the speed of light in cm / fs
66	gezelter	2051	c = 2.99792458e-5
67	gezelter	2045
68			# symmetrize arrays
69
70			tlen = len(time)
71			for i in range(1,len(time)):
72			tnew.append( -time[tlen-i] )
73			vnew.append( vcorr[tlen-i] )
74			for i in range(len(time)):
75			tnew.append( time[i] )
76			vnew.append( vcorr[i] )
77
78			spect = np.fft.fft( vnew )
79			n = len(vnew)
80			timestep = tnew[1] - tnew[0]
81			freq = np.fft.fftfreq(n, d=timestep)
82
83			# freq = np.fft.fftshift(freq)
84			# y = np.fft.fftshift(spect)
85
86			outFile = open(outFileName, 'w')
87
88	gezelter	2051	freqScale = 1.0 / c
89	gezelter	2045	dfreq = (freq[1]-freq[0]) * freqScale
90			s = 0.0
91
92			for i in range(n/2):
93			s = s + np.abs(spect[i])*dfreq
94
95			for i in range(n/2):
96			outFile.write("%lf\t%lf\n" % (freq[i] * freqScale, np.abs(spect[i])/s))
97
98			outFile.close()
99
100			def main(argv):
101			global haveVcorrFileName
102			global haveOutputFileName
103
104			haveVcorrFileName = False
105			haveOutputFileName = False
106			haveQValue = False
107
108			try:
109			opts, args = getopt.getopt(argv, "hf:o:", ["help", "Q-value=", "vcorr-file=", "output-file="])
110			except getopt.GetoptError:
111			usage()
112			sys.exit(2)
113			for opt, arg in opts:
114			if opt in ("-h", "--help"):
115			usage()
116			sys.exit()
117			elif opt in ("-f", "--vcorr-file"):
118			vcorrFileName = arg
119			haveVcorrFileName = True
120			elif opt in ("-o", "--output-file"):
121			outputFileName = arg
122			haveOutputFileName = True
123			if (not haveVcorrFileName):
124			usage()
125			print "No vcorr file was specified"
126			sys.exit()
127			if (not haveOutputFileName):
128			usage()
129			print "No output file was specified"
130			sys.exit()
131
132
133			readVcorrFile(vcorrFileName)
134			computeSpectrum(outputFileName)
135
136
137			if __name__ == "__main__":
138			if len(sys.argv) == 1:
139			usage()
140			sys.exit()
141			main(sys.argv[1:])