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

#!@PYTHON_EXECUTABLE@
"""program that accumulates the static dielectric constant from a stat file.

Accumulates the static dielectric constant from an OpenMD stat file
that has the SYSTEM_DIPOLE added to the statFileFormat.

Usage: stat2dielectric 

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

To use this, the OpenMD file for the run should specify:

statFileFormat = "TIME|TOTAL_ENERGY|POTENTIAL_ENERGY|KINETIC_ENERGY|TEMPERATURE|PRESSURE|VOLUME|CONSERVED_QUANTITY|SYSTEM_DIPOLE";

This will compute the total system dipole and place it in the final
three columns of the stat file.  This program looks for the time in
column 1, the temperature in column 5, the volume in column 7, and the
dipole vector in the last 3 columns of the stat file.

Example:
   stat2dielectric -f stockmayer.stat -o stockmayer.dielectric

"""

__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

def usage():
    print __doc__

def readStatFile(statFileName):
    global time
    global temperature
    global volume
    global boxDipole
    time = []
    temperature = []
    volume = []
    boxDipole = []

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

    print "reading File"
    line = statFile.readline()
    while 1:
        L = line.split()

        time.append(float(L[0]))
        temperature.append(float(L[4]))
        volume.append(float(L[6]))
        dipX = float(L[-3])
        dipY = float(L[-2])
        dipZ = float(L[-1])
        boxDipole.append([dipX, dipY, dipZ])
                       
        line = statFile.readline()
        if not line: break
        
    statFile.close()
    
def computeAverages(outFileName):

    M2 = 0.0
    M = 0.0
    Dx = 0.0
    Dy = 0.0
    Dz = 0.0
    Temp = 0.0
    Vol = 0.0
    
    print "computing Averages"
    outFile = open(outFileName, 'w')

    for i in range(len(time)):

        dipX = boxDipole[i][0]
        dipY = boxDipole[i][1]
        dipZ = boxDipole[i][2]

        length2 = dipX*dipX + dipY*dipY + dipZ*dipZ

        Dx += dipX
        Dy += dipY
        Dz += dipZ
        Mx = Dx / float(1+i)
        My = Dy / float(1+i)
        Mz = Dz / float(1+i)

        M2 += length2
        M += math.sqrt(length2)
        Mavg = Mx*Mx + My*My + Mz*Mz;
        Mavg2 = M / float(1+i)
        M2avg = M2 / float(1+i)

        Temp += temperature[i]
        Vol += volume[i] * 1E-30
        n2 = float(1+i)*float(1+i)

        dielectric1 = 1.0 + 2.7267411E33*(M2avg - Mavg)*n2/(Temp*Vol)
        dielectric2 = 1.0 + 2.7267411E33*(M2avg - Mavg2*Mavg2)*n2/(Temp*Vol)
        outFile.write("%lf\t%lf\t%lf\n" % (time[i], dielectric1, dielectric2))

    outFile.close()


def main(argv):
    global haveStatFileName
    global haveOutputFileName
    
    haveStatFileName = False
    haveOutputFileName = False
 
    try:                                
        opts, args = getopt.getopt(argv, "hf:o:", ["help", "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 ("-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(outputFileName)


if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	2006
Committed:	Mon Jun 30 20:55:21 2014 UTC (10 years, 11 months ago) by gezelter
File size:	4147 byte(s)
Log Message:	Added a script to compute dielectric constants
#	User	Rev	Content
1	gezelter	2006	#!@PYTHON_EXECUTABLE@
2			"""program that accumulates the static dielectric constant from a stat file.
3
4			Accumulates the static dielectric constant from an OpenMD stat file
5			that has the SYSTEM_DIPOLE added to the statFileFormat.
6
7			Usage: stat2dielectric
8
9			Options:
10			-h, --help show this help
11			-f, --stat-file=... use specified stat file
12			-o, --output-file=... use specified output (.dielectric) file
13
14			To use this, the OpenMD file for the run should specify:
15
16			statFileFormat = "TIME\|TOTAL_ENERGY\|POTENTIAL_ENERGY\|KINETIC_ENERGY\|TEMPERATURE\|PRESSURE\|VOLUME\|CONSERVED_QUANTITY\|SYSTEM_DIPOLE";
17
18			This will compute the total system dipole and place it in the final
19			three columns of the stat file. This program looks for the time in
20			column 1, the temperature in column 5, the volume in column 7, and the
21			dipole vector in the last 3 columns of the stat file.
22
23			Example:
24			stat2dielectric -f stockmayer.stat -o stockmayer.dielectric
25
26			"""
27
28			__author__ = "Dan Gezelter (gezelter@nd.edu)"
29			__version__ = "$Revision: $"
30			__date__ = "$Date: $"
31
32			__copyright__ = "Copyright (c) 2014 by the University of Notre Dame"
33			__license__ = "OpenMD"
34
35			import sys
36			import getopt
37			import string
38			import math
39
40			def usage():
41			print __doc__
42
43			def readStatFile(statFileName):
44			global time
45			global temperature
46			global volume
47			global boxDipole
48			time = []
49			temperature = []
50			volume = []
51			boxDipole = []
52
53			statFile = open(statFileName, 'r')
54			line = statFile.readline()
55
56			print "reading File"
57			line = statFile.readline()
58			while 1:
59			L = line.split()
60
61			time.append(float(L[0]))
62			temperature.append(float(L[4]))
63			volume.append(float(L[6]))
64			dipX = float(L[-3])
65			dipY = float(L[-2])
66			dipZ = float(L[-1])
67			boxDipole.append([dipX, dipY, dipZ])
68
69			line = statFile.readline()
70			if not line: break
71
72			statFile.close()
73
74			def computeAverages(outFileName):
75
76			M2 = 0.0
77			M = 0.0
78			Dx = 0.0
79			Dy = 0.0
80			Dz = 0.0
81			Temp = 0.0
82			Vol = 0.0
83
84			print "computing Averages"
85			outFile = open(outFileName, 'w')
86
87			for i in range(len(time)):
88
89			dipX = boxDipole[i][0]
90			dipY = boxDipole[i][1]
91			dipZ = boxDipole[i][2]
92
93			length2 = dipXdipX + dipYdipY + dipZ*dipZ
94
95			Dx += dipX
96			Dy += dipY
97			Dz += dipZ
98			Mx = Dx / float(1+i)
99			My = Dy / float(1+i)
100			Mz = Dz / float(1+i)
101
102			M2 += length2
103			M += math.sqrt(length2)
104			Mavg = MxMx + MyMy + Mz*Mz;
105			Mavg2 = M / float(1+i)
106			M2avg = M2 / float(1+i)
107
108			Temp += temperature[i]
109			Vol += volume[i] * 1E-30
110			n2 = float(1+i)*float(1+i)
111
112			dielectric1 = 1.0 + 2.7267411E33(M2avg - Mavg)n2/(Temp*Vol)
113			dielectric2 = 1.0 + 2.7267411E33(M2avg - Mavg2Mavg2)n2/(TempVol)
114			outFile.write("%lf\t%lf\t%lf\n" % (time[i], dielectric1, dielectric2))
115
116			outFile.close()
117
118
119			def main(argv):
120			global haveStatFileName
121			global haveOutputFileName
122
123			haveStatFileName = False
124			haveOutputFileName = False
125
126			try:
127			opts, args = getopt.getopt(argv, "hf:o:", ["help", "stat-file=", "output-file="])
128			except getopt.GetoptError:
129			usage()
130			sys.exit(2)
131			for opt, arg in opts:
132			if opt in ("-h", "--help"):
133			usage()
134			sys.exit()
135			elif opt in ("-f", "--stat-file"):
136			statFileName = arg
137			haveStatFileName = True
138			elif opt in ("-o", "--output-file"):
139			outputFileName = arg
140			haveOutputFileName = True
141			if (not haveStatFileName):
142			usage()
143			print "No stat file was specified"
144			sys.exit()
145			if (not haveOutputFileName):
146			usage()
147			print "No output file was specified"
148			sys.exit()
149
150			readStatFile(statFileName)
151			computeAverages(outputFileName)
152
153
154			if __name__ == "__main__":
155			if len(sys.argv) == 1:
156			usage()
157			sys.exit()
158			main(sys.argv[1:])