[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, Q):

    # 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)

        x1 = (M2avg - Mavg) / (9.0*eps0*kB*Tavg*Vavg)
        x2 = (M2avg - Mavg2*Mavg2) / (9.0*eps0*kB*Tavg*Vavg)       

        # Clausius-Mossotti
        d1 = (2.0 * x1 + 1.0) / (1.0 - x1)
        d2 = (2.0 * x2 + 1.0) / (1.0 - x2)

        # Conducting boundary conditions
        d3 = 1.0 + 3.0*x1 
        d4 = 1.0 + 3.0*x2 
        
        #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)*(d3 - 1.0) + 3.0)/((Q-1.0) * (d3 - 1.0) + 3.0)
        corrected2 = ((Q+2.0)*(d4 - 1.0) + 3.0)/((Q-1.0) * (d4 - 1.0) + 3.0)

        #corrected1 = (1.0 + 2.0*x1 + Q*x1)/(1.0 - x1 - Q*x1)
        #corrected2 = (1.0 + 2.0*x2 + Q*x2)/(1.0 - x2 - Q*x2)
        
        outFile.write("%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n" % (time[i], x1, x2, d1, d2, d3, d4, 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:	2010
Committed:	Wed Jul 16 14:42:50 2014 UTC (10 years, 10 months ago) by gezelter
File size:	6308 byte(s)
Log Message:	Fixed a selection bug for integer tokens Dielectric script now uses both clausius-mossotti and conducting boundaries
#	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	gezelter	2009	def computeAverages(outFileName, Q):
101	gezelter	2006
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	2010	x1 = (M2avg - Mavg) / (9.0eps0kBTavgVavg)
151			x2 = (M2avg - Mavg2Mavg2) / (9.0eps0kBTavg*Vavg)
152	gezelter	2007
153	gezelter	2010	# Clausius-Mossotti
154			d1 = (2.0 * x1 + 1.0) / (1.0 - x1)
155			d2 = (2.0 * x2 + 1.0) / (1.0 - x2)
156	gezelter	2006
157	gezelter	2010	# Conducting boundary conditions
158			d3 = 1.0 + 3.0*x1
159			d4 = 1.0 + 3.0*x2
160	gezelter	2008
161	gezelter	2010	#d1 = 1.0 + (M2avg - Mavg) / (3.0eps0kBTavgVavg)
162			#d2 = 1.0 + (M2avg - Mavg2Mavg2) / (3.0eps0kBTavg*Vavg)
163	gezelter	2008
164	gezelter	2010	corrected1 = ((Q+2.0)(d3 - 1.0) + 3.0)/((Q-1.0) (d3 - 1.0) + 3.0)
165			corrected2 = ((Q+2.0)(d4 - 1.0) + 3.0)/((Q-1.0) (d4 - 1.0) + 3.0)
166
167			#corrected1 = (1.0 + 2.0x1 + Qx1)/(1.0 - x1 - Q*x1)
168			#corrected2 = (1.0 + 2.0x2 + Qx2)/(1.0 - x2 - Q*x2)
169
170			outFile.write("%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n" % (time[i], x1, x2, d1, d2, d3, d4, corrected1, corrected2))
171
172	gezelter	2006	outFile.close()
173
174
175			def main(argv):
176			global haveStatFileName
177			global haveOutputFileName
178	gezelter	2008	global haveQValue
179	gezelter	2006
180			haveStatFileName = False
181			haveOutputFileName = False
182	gezelter	2008	haveQValue = False
183	gezelter	2006
184			try:
185	gezelter	2008	opts, args = getopt.getopt(argv, "hq:f:o:", ["help", "Q-value=", "stat-file=", "output-file="])
186	gezelter	2006	except getopt.GetoptError:
187			usage()
188			sys.exit(2)
189			for opt, arg in opts:
190			if opt in ("-h", "--help"):
191			usage()
192			sys.exit()
193	gezelter	2008	elif opt in ("-q", "--Q-value="):
194			Q = float(arg)
195			haveQValue = True
196	gezelter	2006	elif opt in ("-f", "--stat-file"):
197			statFileName = arg
198			haveStatFileName = True
199			elif opt in ("-o", "--output-file"):
200			outputFileName = arg
201			haveOutputFileName = True
202			if (not haveStatFileName):
203			usage()
204			print "No stat file was specified"
205			sys.exit()
206			if (not haveOutputFileName):
207			usage()
208			print "No output file was specified"
209			sys.exit()
210	gezelter	2008
211
212	gezelter	2006	readStatFile(statFileName)
213	gezelter	2008	if (not haveQValue):
214			computeAverages(outputFileName, 1.0)
215			else:
216			computeAverages(outputFileName, Q)
217	gezelter	2006
218
219			if __name__ == "__main__":
220			if len(sys.argv) == 1:
221			usage()
222			sys.exit()
223			main(sys.argv[1:])