| 1 |
#!@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 = dipX*dipX + dipY*dipY + 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 = Mx*Mx + My*My + 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 - Mavg2*Mavg2)*n2/(Temp*Vol) |
| 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:]) |
