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

#!@PYTHON_EXECUTABLE@
"""A script that computes the static dielectric constant

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

This assumes the fluctuation formula appropriate for conducting
boundaries:

                    <M^2> - <M>^2
        eps = 1 + -----------------
                  3 kB <T> <V> eps0
   
        eps:  dielectric constant
        M:    total dipole moment of the box
        <V>:  average volume of the box
        <T>:  average temperature
        eps0: vacuum permittivity
        kB:   Boltzmann constant

See M. Neumann, O. Steinhauser, and G. S. Pawley, Mol. Phys. 52, 97 (1984).

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 line will get OpenMD to compute the total system dipole and place
it in the final three columns of the stat file.  The stat2dielectric
script looks for 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):

    # permittivity in C^2 N^-1 m^-2
    eps0 = 8.854187817620E-12
    # Boltzmann constant in J / K
    kB = 1.380648E-23
    # Convert angstroms^3 to m^3
    a3tom3 = 1.0E-30

    M2 = 0.0
    M = 0.0
    Dx = 0.0
    Dy = 0.0
    Dz = 0.0
    Temp = 0.0
    Vol = 0.0
    
    print "computing Dielectrics"
    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
        M2 += length2
        M2avg = M2 / float(1+i)

        Dx += dipX
        Dy += dipY
        Dz += dipZ

        # the average of each of the three components of the box dipole
        Mx = Dx / float(1+i)
        My = Dy / float(1+i)
        Mz = Dz / float(1+i)
        Mavg = Mx*Mx + My*My + Mz*Mz

        # the average of the box dipole vector length
        M += math.sqrt(length2)
        Mavg2 = M / float(1+i)

        Temp += temperature[i]
        Tavg = Temp / float(1+i)
        
        Vol += volume[i] * a3tom3
        Vavg = Vol / float(1+i)

        dielectric1 = 1.0 + (M2avg - Mavg) / (3.0*eps0*kB*Tavg*Vavg)
        dielectric2 = 1.0 + (M2avg - Mavg2*Mavg2) / (3.0*eps0*kB*Tavg*Vavg)

        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:	2007
Committed:	Wed Jul 2 19:03:53 2014 UTC (11 years, 4 months ago) by gezelter
File size:	4982 byte(s)
Log Message:	Most recent version
#	User	Rev	Content
1	gezelter	2006	#!@PYTHON_EXECUTABLE@
2	gezelter	2007	"""A script that computes the static dielectric constant
3	gezelter	2006
4			Accumulates the static dielectric constant from an OpenMD stat file
5	gezelter	2007	that has been run with the SYSTEM_DIPOLE added to the statFileFormat.
6	gezelter	2006
7	gezelter	2007	This assumes the fluctuation formula appropriate for conducting
8			boundaries:
9
10			<M^2> - <M>^2
11			eps = 1 + -----------------
12			3 kB <T> <V> eps0
13
14			eps: dielectric constant
15			M: total dipole moment of the box
16			<V>: average volume of the box
17			<T>: average temperature
18			eps0: vacuum permittivity
19			kB: Boltzmann constant
20
21			See M. Neumann, O. Steinhauser, and G. S. Pawley, Mol. Phys. 52, 97 (1984).
22
23	gezelter	2006	Usage: stat2dielectric
24
25			Options:
26			-h, --help show this help
27			-f, --stat-file=... use specified stat file
28			-o, --output-file=... use specified output (.dielectric) file
29
30			To use this, the OpenMD file for the run should specify:
31
32			statFileFormat = "TIME\|TOTAL_ENERGY\|POTENTIAL_ENERGY\|KINETIC_ENERGY\|TEMPERATURE\|PRESSURE\|VOLUME\|CONSERVED_QUANTITY\|SYSTEM_DIPOLE";
33
34	gezelter	2007	This line will get OpenMD to compute the total system dipole and place
35			it in the final three columns of the stat file. The stat2dielectric
36			script looks for the temperature in column 5, the volume in column 7,
37			and the dipole vector in the last 3 columns of the stat file.
38	gezelter	2006
39			Example:
40			stat2dielectric -f stockmayer.stat -o stockmayer.dielectric
41
42			"""
43
44			__author__ = "Dan Gezelter (gezelter@nd.edu)"
45			__version__ = "$Revision: $"
46			__date__ = "$Date: $"
47
48			__copyright__ = "Copyright (c) 2014 by the University of Notre Dame"
49			__license__ = "OpenMD"
50
51			import sys
52			import getopt
53			import string
54			import math
55
56			def usage():
57			print __doc__
58
59			def readStatFile(statFileName):
60			global time
61			global temperature
62			global volume
63			global boxDipole
64			time = []
65			temperature = []
66			volume = []
67			boxDipole = []
68
69			statFile = open(statFileName, 'r')
70			line = statFile.readline()
71
72			print "reading File"
73			line = statFile.readline()
74			while 1:
75			L = line.split()
76
77			time.append(float(L[0]))
78			temperature.append(float(L[4]))
79			volume.append(float(L[6]))
80			dipX = float(L[-3])
81			dipY = float(L[-2])
82			dipZ = float(L[-1])
83			boxDipole.append([dipX, dipY, dipZ])
84
85			line = statFile.readline()
86			if not line: break
87
88			statFile.close()
89
90			def computeAverages(outFileName):
91
92	gezelter	2007	# permittivity in C^2 N^-1 m^-2
93			eps0 = 8.854187817620E-12
94			# Boltzmann constant in J / K
95			kB = 1.380648E-23
96			# Convert angstroms^3 to m^3
97			a3tom3 = 1.0E-30
98
99	gezelter	2006	M2 = 0.0
100			M = 0.0
101			Dx = 0.0
102			Dy = 0.0
103			Dz = 0.0
104			Temp = 0.0
105			Vol = 0.0
106
107	gezelter	2007	print "computing Dielectrics"
108	gezelter	2006	outFile = open(outFileName, 'w')
109
110			for i in range(len(time)):
111
112			dipX = boxDipole[i][0]
113			dipY = boxDipole[i][1]
114			dipZ = boxDipole[i][2]
115
116			length2 = dipXdipX + dipYdipY + dipZ*dipZ
117	gezelter	2007	M2 += length2
118			M2avg = M2 / float(1+i)
119	gezelter	2006
120			Dx += dipX
121			Dy += dipY
122			Dz += dipZ
123	gezelter	2007
124			# the average of each of the three components of the box dipole
125	gezelter	2006	Mx = Dx / float(1+i)
126			My = Dy / float(1+i)
127			Mz = Dz / float(1+i)
128	gezelter	2007	Mavg = MxMx + MyMy + Mz*Mz
129	gezelter	2006
130	gezelter	2007	# the average of the box dipole vector length
131	gezelter	2006	M += math.sqrt(length2)
132			Mavg2 = M / float(1+i)
133
134			Temp += temperature[i]
135	gezelter	2007	Tavg = Temp / float(1+i)
136
137			Vol += volume[i] * a3tom3
138			Vavg = Vol / float(1+i)
139	gezelter	2006
140	gezelter	2007	dielectric1 = 1.0 + (M2avg - Mavg) / (3.0eps0kBTavgVavg)
141			dielectric2 = 1.0 + (M2avg - Mavg2Mavg2) / (3.0eps0kBTavg*Vavg)
142
143	gezelter	2006	outFile.write("%lf\t%lf\t%lf\n" % (time[i], dielectric1, dielectric2))
144
145			outFile.close()
146
147
148			def main(argv):
149			global haveStatFileName
150			global haveOutputFileName
151
152			haveStatFileName = False
153			haveOutputFileName = False
154
155			try:
156			opts, args = getopt.getopt(argv, "hf:o:", ["help", "stat-file=", "output-file="])
157			except getopt.GetoptError:
158			usage()
159			sys.exit(2)
160			for opt, arg in opts:
161			if opt in ("-h", "--help"):
162			usage()
163			sys.exit()
164			elif opt in ("-f", "--stat-file"):
165			statFileName = arg
166			haveStatFileName = True
167			elif opt in ("-o", "--output-file"):
168			outputFileName = arg
169			haveOutputFileName = True
170			if (not haveStatFileName):
171			usage()
172			print "No stat file was specified"
173			sys.exit()
174			if (not haveOutputFileName):
175			usage()
176			print "No output file was specified"
177			sys.exit()
178
179			readStatFile(statFileName)
180			computeAverages(outputFileName)
181
182
183			if __name__ == "__main__":
184			if len(sys.argv) == 1:
185			usage()
186			sys.exit()
187			main(sys.argv[1:])