[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
        epsA = 1 + -----------------
                   3 kB <T> <V> eps0
   
        epsA: 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).

Also, optionally applies a correction factor as follows:

                ((Q + 2) (epsA - 1) + 3)
          eps = ------------------------
                ((Q - 1) (epsA - 1) + 3)

Where Q depends on the method used to compute the electrostatic interaction.
See M. Neumann and O. Steinhauser, Chem. Phys. Lett. 95, 417 (1983))

Usage: stat2dielectric 

Options:
  -h, --help              show this help
  -f, --stat-file=...     use specified stat file
  -o, --output-file=...   use specified output (.dielectric) file
  -q, --Q-value=...       use the specified Q value to correct the dielectric

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)

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


        corrected1 = ((Q+2.0)*(d1 - 1.0) + 3.0)/((Q-1.0) * (d1 - 1.0) + 3.0)
        corrected2 = ((Q+2.0)*(d2 - 1.0) + 3.0)/((Q-1.0) * (d2 - 1.0) + 3.0)
        
        outFile.write("%lf\t%lf\t%lf\t%lf\t%lf\n" % (time[i], d1, d2, corrected1, corrected2))

    outFile.close()


def main(argv):
    global haveStatFileName
    global haveOutputFileName
    global haveQValue
    
    haveStatFileName = False
    haveOutputFileName = False
    haveQValue = False
 
    try:                                
        opts, args = getopt.getopt(argv, "hq:f:o:", ["help", "Q-value=", "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 ("-q", "--Q-value="): 
            Q = float(arg)
            haveQValue = 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)
    if (not haveQValue):
        computeAverages(outputFileName, 1.0)
    else:
        computeAverages(outputFileName, Q)


if __name__ == "__main__":
    if len(sys.argv) == 1:
        usage()
        sys.exit()
    main(sys.argv[1:])
Revision:	2008
Committed:	Tue Jul 8 14:54:55 2014 UTC (10 years, 10 months ago) by gezelter
File size:	5801 byte(s)
Log Message:	Added correction factor to dielectric script
#	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	gezelter	2008	<M^2> - <M>^2
11			epsA = 1 + -----------------
12			3 kB <T> <V> eps0
13	gezelter	2007
14	gezelter	2008	epsA: dielectric constant
15	gezelter	2007	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	2008	Also, optionally applies a correction factor as follows:
24
25			((Q + 2) (epsA - 1) + 3)
26			eps = ------------------------
27			((Q - 1) (epsA - 1) + 3)
28
29			Where Q depends on the method used to compute the electrostatic interaction.
30			See M. Neumann and O. Steinhauser, Chem. Phys. Lett. 95, 417 (1983))
31
32	gezelter	2006	Usage: stat2dielectric
33
34			Options:
35			-h, --help show this help
36			-f, --stat-file=... use specified stat file
37			-o, --output-file=... use specified output (.dielectric) file
38	gezelter	2008	-q, --Q-value=... use the specified Q value to correct the dielectric
39	gezelter	2006
40			To use this, the OpenMD file for the run should specify:
41
42			statFileFormat = "TIME\|TOTAL_ENERGY\|POTENTIAL_ENERGY\|KINETIC_ENERGY\|TEMPERATURE\|PRESSURE\|VOLUME\|CONSERVED_QUANTITY\|SYSTEM_DIPOLE";
43
44	gezelter	2007	This line will get OpenMD to compute the total system dipole and place
45			it in the final three columns of the stat file. The stat2dielectric
46			script looks for the temperature in column 5, the volume in column 7,
47			and the dipole vector in the last 3 columns of the stat file.
48	gezelter	2006
49			Example:
50			stat2dielectric -f stockmayer.stat -o stockmayer.dielectric
51
52			"""
53
54			__author__ = "Dan Gezelter (gezelter@nd.edu)"
55			__version__ = "$Revision: $"
56			__date__ = "$Date: $"
57
58			__copyright__ = "Copyright (c) 2014 by the University of Notre Dame"
59			__license__ = "OpenMD"
60
61			import sys
62			import getopt
63			import string
64			import math
65
66			def usage():
67			print __doc__
68
69			def readStatFile(statFileName):
70			global time
71			global temperature
72			global volume
73			global boxDipole
74			time = []
75			temperature = []
76			volume = []
77			boxDipole = []
78
79			statFile = open(statFileName, 'r')
80			line = statFile.readline()
81
82			print "reading File"
83			line = statFile.readline()
84			while 1:
85			L = line.split()
86
87			time.append(float(L[0]))
88			temperature.append(float(L[4]))
89			volume.append(float(L[6]))
90			dipX = float(L[-3])
91			dipY = float(L[-2])
92			dipZ = float(L[-1])
93			boxDipole.append([dipX, dipY, dipZ])
94
95			line = statFile.readline()
96			if not line: break
97
98			statFile.close()
99
100			def computeAverages(outFileName):
101
102	gezelter	2007	# permittivity in C^2 N^-1 m^-2
103			eps0 = 8.854187817620E-12
104			# Boltzmann constant in J / K
105			kB = 1.380648E-23
106			# Convert angstroms^3 to m^3
107			a3tom3 = 1.0E-30
108
109	gezelter	2006	M2 = 0.0
110			M = 0.0
111			Dx = 0.0
112			Dy = 0.0
113			Dz = 0.0
114			Temp = 0.0
115			Vol = 0.0
116
117	gezelter	2007	print "computing Dielectrics"
118	gezelter	2006	outFile = open(outFileName, 'w')
119
120			for i in range(len(time)):
121
122			dipX = boxDipole[i][0]
123			dipY = boxDipole[i][1]
124			dipZ = boxDipole[i][2]
125
126			length2 = dipXdipX + dipYdipY + dipZ*dipZ
127	gezelter	2007	M2 += length2
128			M2avg = M2 / float(1+i)
129	gezelter	2006
130			Dx += dipX
131			Dy += dipY
132			Dz += dipZ
133	gezelter	2007
134			# the average of each of the three components of the box dipole
135	gezelter	2006	Mx = Dx / float(1+i)
136			My = Dy / float(1+i)
137			Mz = Dz / float(1+i)
138	gezelter	2007	Mavg = MxMx + MyMy + Mz*Mz
139	gezelter	2006
140	gezelter	2007	# the average of the box dipole vector length
141	gezelter	2006	M += math.sqrt(length2)
142			Mavg2 = M / float(1+i)
143
144			Temp += temperature[i]
145	gezelter	2007	Tavg = Temp / float(1+i)
146
147			Vol += volume[i] * a3tom3
148			Vavg = Vol / float(1+i)
149	gezelter	2006
150	gezelter	2008	d1 = 1.0 + (M2avg - Mavg) / (3.0eps0kBTavgVavg)
151			d2 = 1.0 + (M2avg - Mavg2Mavg2) / (3.0eps0kBTavg*Vavg)
152	gezelter	2007
153	gezelter	2006
154	gezelter	2008	corrected1 = ((Q+2.0)(d1 - 1.0) + 3.0)/((Q-1.0) (d1 - 1.0) + 3.0)
155			corrected2 = ((Q+2.0)(d2 - 1.0) + 3.0)/((Q-1.0) (d2 - 1.0) + 3.0)
156
157			outFile.write("%lf\t%lf\t%lf\t%lf\t%lf\n" % (time[i], d1, d2, corrected1, corrected2))
158
159	gezelter	2006	outFile.close()
160
161
162			def main(argv):
163			global haveStatFileName
164			global haveOutputFileName
165	gezelter	2008	global haveQValue
166	gezelter	2006
167			haveStatFileName = False
168			haveOutputFileName = False
169	gezelter	2008	haveQValue = False
170	gezelter	2006
171			try:
172	gezelter	2008	opts, args = getopt.getopt(argv, "hq:f:o:", ["help", "Q-value=", "stat-file=", "output-file="])
173	gezelter	2006	except getopt.GetoptError:
174			usage()
175			sys.exit(2)
176			for opt, arg in opts:
177			if opt in ("-h", "--help"):
178			usage()
179			sys.exit()
180	gezelter	2008	elif opt in ("-q", "--Q-value="):
181			Q = float(arg)
182			haveQValue = True
183	gezelter	2006	elif opt in ("-f", "--stat-file"):
184			statFileName = arg
185			haveStatFileName = True
186			elif opt in ("-o", "--output-file"):
187			outputFileName = arg
188			haveOutputFileName = True
189			if (not haveStatFileName):
190			usage()
191			print "No stat file was specified"
192			sys.exit()
193			if (not haveOutputFileName):
194			usage()
195			print "No output file was specified"
196			sys.exit()
197	gezelter	2008
198
199	gezelter	2006	readStatFile(statFileName)
200	gezelter	2008	if (not haveQValue):
201			computeAverages(outputFileName, 1.0)
202			else:
203			computeAverages(outputFileName, Q)
204	gezelter	2006
205
206			if __name__ == "__main__":
207			if len(sys.argv) == 1:
208			usage()
209			sys.exit()
210			main(sys.argv[1:])