| 1 |
|
#!@PYTHON_EXECUTABLE@ |
| 2 |
< |
"""program that accumulates the static dielectric constant from a stat file. |
| 2 |
> |
"""A script that computes the static dielectric constant |
| 3 |
|
|
| 4 |
|
Accumulates the static dielectric constant from an OpenMD stat file |
| 5 |
< |
that has the SYSTEM_DIPOLE added to the statFileFormat. |
| 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 |
+ |
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 |
|
Usage: stat2dielectric |
| 24 |
|
|
| 25 |
|
Options: |
| 31 |
|
|
| 32 |
|
statFileFormat = "TIME|TOTAL_ENERGY|POTENTIAL_ENERGY|KINETIC_ENERGY|TEMPERATURE|PRESSURE|VOLUME|CONSERVED_QUANTITY|SYSTEM_DIPOLE"; |
| 33 |
|
|
| 34 |
< |
This will compute the total system dipole and place it in the final |
| 35 |
< |
three columns of the stat file. This program looks for the time in |
| 36 |
< |
column 1, the temperature in column 5, the volume in column 7, and the |
| 37 |
< |
dipole vector in the last 3 columns of the stat file. |
| 34 |
> |
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 |
|
|
| 39 |
|
Example: |
| 40 |
|
stat2dielectric -f stockmayer.stat -o stockmayer.dielectric |
| 89 |
|
|
| 90 |
|
def computeAverages(outFileName): |
| 91 |
|
|
| 92 |
+ |
# 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 |
|
M2 = 0.0 |
| 100 |
|
M = 0.0 |
| 101 |
|
Dx = 0.0 |
| 104 |
|
Temp = 0.0 |
| 105 |
|
Vol = 0.0 |
| 106 |
|
|
| 107 |
< |
print "computing Averages" |
| 107 |
> |
print "computing Dielectrics" |
| 108 |
|
outFile = open(outFileName, 'w') |
| 109 |
|
|
| 110 |
|
for i in range(len(time)): |
| 114 |
|
dipZ = boxDipole[i][2] |
| 115 |
|
|
| 116 |
|
length2 = dipX*dipX + dipY*dipY + dipZ*dipZ |
| 117 |
+ |
M2 += length2 |
| 118 |
+ |
M2avg = M2 / float(1+i) |
| 119 |
|
|
| 120 |
|
Dx += dipX |
| 121 |
|
Dy += dipY |
| 122 |
|
Dz += dipZ |
| 123 |
+ |
|
| 124 |
+ |
# the average of each of the three components of the box dipole |
| 125 |
|
Mx = Dx / float(1+i) |
| 126 |
|
My = Dy / float(1+i) |
| 127 |
|
Mz = Dz / float(1+i) |
| 128 |
+ |
Mavg = Mx*Mx + My*My + Mz*Mz |
| 129 |
|
|
| 130 |
< |
M2 += length2 |
| 130 |
> |
# the average of the box dipole vector length |
| 131 |
|
M += math.sqrt(length2) |
| 104 |
– |
Mavg = Mx*Mx + My*My + Mz*Mz; |
| 132 |
|
Mavg2 = M / float(1+i) |
| 106 |
– |
M2avg = M2 / float(1+i) |
| 133 |
|
|
| 134 |
|
Temp += temperature[i] |
| 135 |
< |
Vol += volume[i] * 1E-30 |
| 136 |
< |
n2 = float(1+i)*float(1+i) |
| 135 |
> |
Tavg = Temp / float(1+i) |
| 136 |
> |
|
| 137 |
> |
Vol += volume[i] * a3tom3 |
| 138 |
> |
Vavg = Vol / float(1+i) |
| 139 |
|
|
| 140 |
< |
dielectric1 = 1.0 + 2.7267411E33*(M2avg - Mavg)*n2/(Temp*Vol) |
| 141 |
< |
dielectric2 = 1.0 + 2.7267411E33*(M2avg - Mavg2*Mavg2)*n2/(Temp*Vol) |
| 140 |
> |
dielectric1 = 1.0 + (M2avg - Mavg) / (3.0*eps0*kB*Tavg*Vavg) |
| 141 |
> |
dielectric2 = 1.0 + (M2avg - Mavg2*Mavg2) / (3.0*eps0*kB*Tavg*Vavg) |
| 142 |
> |
|
| 143 |
|
outFile.write("%lf\t%lf\t%lf\n" % (time[i], dielectric1, dielectric2)) |
| 144 |
|
|
| 145 |
|
outFile.close() |
