| 1 |
#!@PYTHON_EXECUTABLE@ |
| 2 |
"""A script that computes the static dielectric constant |
| 3 |
|
| 4 |
Accumulates the static dielectric constant from an OpenMD stat file |
| 5 |
that has been run with the SYSTEM_DIPOLE added to the statFileFormat. |
| 6 |
|
| 7 |
This assumes the fluctuation formula appropriate for conducting |
| 8 |
boundaries: |
| 9 |
|
| 10 |
<M^2> - <M>^2 |
| 11 |
epsA = 1 + ----------------- |
| 12 |
3 kB <T> <V> eps0 |
| 13 |
|
| 14 |
epsA: 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 |
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 |
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 |
-q, --Q-value=... use the specified Q value to correct the dielectric |
| 39 |
|
| 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 |
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 |
|
| 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, Q): |
| 101 |
|
| 102 |
# 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 |
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 |
print "computing Dielectrics" |
| 118 |
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 = dipX*dipX + dipY*dipY + dipZ*dipZ |
| 127 |
M2 += length2 |
| 128 |
M2avg = M2 / float(1+i) |
| 129 |
|
| 130 |
Dx += dipX |
| 131 |
Dy += dipY |
| 132 |
Dz += dipZ |
| 133 |
|
| 134 |
# the average of each of the three components of the box dipole |
| 135 |
Mx = Dx / float(1+i) |
| 136 |
My = Dy / float(1+i) |
| 137 |
Mz = Dz / float(1+i) |
| 138 |
Mavg = Mx*Mx + My*My + Mz*Mz |
| 139 |
|
| 140 |
# the average of the box dipole vector length |
| 141 |
M += math.sqrt(length2) |
| 142 |
Mavg2 = M / float(1+i) |
| 143 |
|
| 144 |
Temp += temperature[i] |
| 145 |
Tavg = Temp / float(1+i) |
| 146 |
|
| 147 |
Vol += volume[i] * a3tom3 |
| 148 |
Vavg = Vol / float(1+i) |
| 149 |
|
| 150 |
x1 = (M2avg - Mavg) / (9.0*eps0*kB*Tavg*Vavg) |
| 151 |
x2 = (M2avg - Mavg2*Mavg2) / (9.0*eps0*kB*Tavg*Vavg) |
| 152 |
|
| 153 |
# Clausius-Mossotti |
| 154 |
d1 = (2.0 * x1 + 1.0) / (1.0 - x1) |
| 155 |
d2 = (2.0 * x2 + 1.0) / (1.0 - x2) |
| 156 |
|
| 157 |
# Conducting boundary conditions |
| 158 |
d3 = 1.0 + 3.0*x1 |
| 159 |
d4 = 1.0 + 3.0*x2 |
| 160 |
|
| 161 |
#d1 = 1.0 + (M2avg - Mavg) / (3.0*eps0*kB*Tavg*Vavg) |
| 162 |
#d2 = 1.0 + (M2avg - Mavg2*Mavg2) / (3.0*eps0*kB*Tavg*Vavg) |
| 163 |
|
| 164 |
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.0*x1 + Q*x1)/(1.0 - x1 - Q*x1) |
| 168 |
#corrected2 = (1.0 + 2.0*x2 + Q*x2)/(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 |
outFile.close() |
| 173 |
|
| 174 |
|
| 175 |
def main(argv): |
| 176 |
global haveStatFileName |
| 177 |
global haveOutputFileName |
| 178 |
global haveQValue |
| 179 |
|
| 180 |
haveStatFileName = False |
| 181 |
haveOutputFileName = False |
| 182 |
haveQValue = False |
| 183 |
|
| 184 |
try: |
| 185 |
opts, args = getopt.getopt(argv, "hq:f:o:", ["help", "Q-value=", "stat-file=", "output-file="]) |
| 186 |
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 |
elif opt in ("-q", "--Q-value="): |
| 194 |
Q = float(arg) |
| 195 |
haveQValue = True |
| 196 |
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 |
|
| 211 |
|
| 212 |
readStatFile(statFileName) |
| 213 |
if (not haveQValue): |
| 214 |
computeAverages(outputFileName, 1.0) |
| 215 |
else: |
| 216 |
computeAverages(outputFileName, Q) |
| 217 |
|
| 218 |
|
| 219 |
if __name__ == "__main__": |
| 220 |
if len(sys.argv) == 1: |
| 221 |
usage() |
| 222 |
sys.exit() |
| 223 |
main(sys.argv[1:]) |
