1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
from .vgrid import (
startVio,
Vgrid_ctor,
Vgrid_readDX,
Vgrid_value,
)
import sys
import math
from sys import stdout, stderr
"""
average.py - An example script for interfacing Python with APBS
Vgrid routines
"""
header = "\n\n\
----------------------------------------------------------------------\n\
Adaptive Poisson-Boltzmann Solver (APBS)\n\
----------------------------------------------------------------------\n\
\n\n"
usage = "python[2] average.py file.dx\n"
def main():
# *************** CHECK INVOCATION *******************
startVio()
if len(sys.argv) != 2:
stderr.write(
"\n*** Syntax error: got %d arguments, expected 2.\n\n"
% len(sys.argv)
)
stderr.write("%s\n" % usage)
sys.exit(2)
else:
inpath = sys.argv[1]
# *************** APBS INITIALIZATION *******************
stdout.write(header)
data = []
stdout.write("main: Reading data from %s... \n" % inpath)
grid = Vgrid_ctor(0, 0, 0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, data)
Vgrid_readDX(grid, "FILE", "ASC", "", inpath)
nx = grid.nx
ny = grid.ny
nz = grid.nz
hx = grid.hx
hy = grid.hy
hzed = grid.hzed
xmin = grid.xmin
ymin = grid.ymin
zmin = grid.zmin
stdout.write("# nx = %d, ny = %d, nz = %d\n" % (nx, ny, nz))
stdout.write("# hx = %g, hy = %g, hz = %g\n" % (hx, hy, hzed))
stdout.write(
"# xmin = %g, ymin = %g, zmin = %g\n" % (xmin, ymin, zmin)
)
# *************** SETTINGS ********************
xcentAVG = 112.160
ycentAVG = 63.5
zcentAVG = 137.245
xlenAVG = 70.0
zlenAVG = 70.0
ylenAVG = hy * (ny - 1)
# *************** AVERAGE **********************
xminAVG = xcentAVG - 0.5 * xlenAVG
xmaxAVG = xcentAVG + 0.5 * xlenAVG
yminAVG = ycentAVG - 0.5 * ylenAVG
ymaxAVG = ycentAVG + 0.5 * ylenAVG
zminAVG = zcentAVG - 0.5 * zlenAVG
zmaxAVG = zcentAVG + 0.5 * zlenAVG
imin = int(math.floor((xminAVG - xmin) / hx))
imin = max(imin, 0)
jmin = int(math.floor((yminAVG - ymin) / hy))
jmin = max(jmin, 0)
kmin = int(math.floor((zminAVG - zmin) / hzed))
kmin = max(kmin, 0)
imax = int(math.ceil((xmaxAVG - xmin) / hx))
imax = min(imax, nx - 1)
jmax = int(math.ceil((ymaxAVG - ymin) / hy))
jmax = min(jmax, ny - 1)
kmax = int(math.ceil((zmaxAVG - zmin) / hzed))
kmax = min(kmax, nz - 1)
stdout.write("# \tY POS\t\tAVERAGE\n")
for j in range(jmin, jmax):
avg = 0.0
navg = 0
for k in range(kmin, kmax):
for i in range(imin, imax):
pt = [i, j, k]
val = 0.0
ret, value = Vgrid_value(grid, pt, val)
if ret:
avg = avg + value
navg = navg + 1
if navg != 0:
avg = avg / navg
stdout.write(" \t%e\t\t%e\n" % ((hy * j + ymin), avg))
else:
stdout.write(" \t%e\t\t%s\n" % ((hy * j + ymin), "nan"))
if __name__ == "__main__":
main()
|
