| 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 |
| 10 |
> |
<M^2> - <M>^2 |
| 11 |
> |
epsA = 1 + ----------------- |
| 12 |
> |
3 kB <T> <V> eps0 |
| 13 |
|
|
| 14 |
< |
eps: dielectric constant |
| 14 |
> |
epsA: dielectric constant |
| 15 |
|
M: total dipole moment of the box |
| 16 |
|
<V>: average volume of the box |
| 17 |
|
<T>: average temperature |
| 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 |
|
|
| 147 |
|
Vol += volume[i] * a3tom3 |
| 148 |
|
Vavg = Vol / float(1+i) |
| 149 |
|
|
| 150 |
< |
dielectric1 = 1.0 + (M2avg - Mavg) / (3.0*eps0*kB*Tavg*Vavg) |
| 151 |
< |
dielectric2 = 1.0 + (M2avg - Mavg2*Mavg2) / (3.0*eps0*kB*Tavg*Vavg) |
| 150 |
> |
d1 = 1.0 + (M2avg - Mavg) / (3.0*eps0*kB*Tavg*Vavg) |
| 151 |
> |
d2 = 1.0 + (M2avg - Mavg2*Mavg2) / (3.0*eps0*kB*Tavg*Vavg) |
| 152 |
|
|
| 143 |
– |
outFile.write("%lf\t%lf\t%lf\n" % (time[i], dielectric1, dielectric2)) |
| 153 |
|
|
| 154 |
+ |
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 |
|
outFile.close() |
| 160 |
|
|
| 161 |
|
|
| 162 |
|
def main(argv): |
| 163 |
|
global haveStatFileName |
| 164 |
|
global haveOutputFileName |
| 165 |
+ |
global haveQValue |
| 166 |
|
|
| 167 |
|
haveStatFileName = False |
| 168 |
|
haveOutputFileName = False |
| 169 |
+ |
haveQValue = False |
| 170 |
|
|
| 171 |
|
try: |
| 172 |
< |
opts, args = getopt.getopt(argv, "hf:o:", ["help", "stat-file=", "output-file="]) |
| 172 |
> |
opts, args = getopt.getopt(argv, "hq:f:o:", ["help", "Q-value=", "stat-file=", "output-file="]) |
| 173 |
|
except getopt.GetoptError: |
| 174 |
|
usage() |
| 175 |
|
sys.exit(2) |
| 177 |
|
if opt in ("-h", "--help"): |
| 178 |
|
usage() |
| 179 |
|
sys.exit() |
| 180 |
+ |
elif opt in ("-q", "--Q-value="): |
| 181 |
+ |
Q = float(arg) |
| 182 |
+ |
haveQValue = True |
| 183 |
|
elif opt in ("-f", "--stat-file"): |
| 184 |
|
statFileName = arg |
| 185 |
|
haveStatFileName = True |
| 194 |
|
usage() |
| 195 |
|
print "No output file was specified" |
| 196 |
|
sys.exit() |
| 197 |
< |
|
| 197 |
> |
|
| 198 |
> |
|
| 199 |
|
readStatFile(statFileName) |
| 200 |
< |
computeAverages(outputFileName) |
| 200 |
> |
if (not haveQValue): |
| 201 |
> |
computeAverages(outputFileName, 1.0) |
| 202 |
> |
else: |
| 203 |
> |
computeAverages(outputFileName, Q) |
| 204 |
|
|
| 205 |
|
|
| 206 |
|
if __name__ == "__main__": |
