Coverage for gpaw/lrtddft2/tools.py: 4%
157 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1import sys
2import re
3import numpy as np
4import datetime
7def write_parallel_cube(basename, data, gd, atoms, rank):
8 hx = gd.h_cv[0, 0]
9 hy = gd.h_cv[1, 1]
10 hz = gd.h_cv[2, 2]
12 f = open('%s_%06d.cube' % (basename, rank), 'w', 256 * 1024)
14 f.write(
15 ('GPAW Global (natom nx ny nz hx hy hz): '
16 '%6d %6d %6d %6d %12.6lf %12.6lf %12.6lf\n')
17 % (len(atoms), gd.N_c[0], gd.N_c[1], gd.N_c[2], hx, hy, hz))
18 f.write(
19 ('GPAW Local (xbeg xend ybeg yend zbeg zend): '
20 '%6d %6d %6d %6d %6d %6d\n')
21 % (gd.beg_c[0], gd.end_c[0], gd.beg_c[1], gd.end_c[1], gd.beg_c[2],
22 gd.end_c[2]))
24 f.write('%6d %12.6lf %12.6lf %12.6lf\n' %
25 (len(atoms), hx * gd.beg_c[0], hy * gd.beg_c[1], hz * gd.beg_c[2]))
26 f.write('%6d %12.6lf %12.6lf %12.6lf\n' % (gd.n_c[0], hx, 0.0, 0.0))
27 f.write('%6d %12.6lf %12.6lf %12.6lf\n' % (gd.n_c[1], 0.0, hy, 0.0))
28 f.write('%6d %12.6lf %12.6lf %12.6lf\n' % (gd.n_c[2], 0.0, 0.0, hz))
30 for (i, atom) in enumerate(atoms):
31 f.write('%6d %12.6lf %12.6lf %12.6lf %12.6lf\n' %
32 (atom.number, 0.0, atom.position[0] / 0.529177,
33 atom.position[1] / 0.529177, atom.position[2] / 0.529177))
34 for i in range(data.shape[0]):
35 for j in range(data.shape[1]):
36 for k in range(data.shape[2]):
37 f.write('%18.8lf\n' % data[i, j, k])
38 f.close()
41def cubify(out_filename, filenames):
42 sys.stderr.write('Reading partial cube files...%s\n' %
43 datetime.datetime.now())
44 data = None
45 for fname in filenames:
46 f = open(fname, 'r', 256 * 1024)
47 line0 = f.readline()
48 elems = re.compile(r'[\d.e+-]+').findall(line0)
49 natom = int(elems[0])
50 (nx, ny, nz) = (int(n) for n in elems[1:4])
51 (hx, hy, hz) = (float(h) for h in elems[4:7])
52 line = f.readline()
53 (xb, xe, yb, ye, zb,
54 ze) = (int(x) for x in re.compile(r'\d+').findall(line))
55 sys.stderr.write('%d %d %d %d %d %d\n' % (xb, xe, yb, ye, zb, ze))
57 line = f.readline()
58 line = f.readline()
59 line = f.readline()
60 line = f.readline()
62 atom_lines = []
63 for n in range(natom):
64 atom_lines.append(f.readline())
66 if data is None:
67 data = np.zeros((nx, ny, nz))
69 for i in range(xb, xe):
70 for j in range(yb, ye):
71 for k in range(zb, ze):
72 data[i, j, k] = float(f.readline())
74 f.close()
76 sys.stderr.write('Reading done. %s\n' % datetime.datetime.now())
77 sys.stderr.write('Writing the cube file... %d %d %d\n' % (nx, ny, nz))
79 out = open(out_filename, 'w', 256 * 1024)
80 out.write(line0)
81 out.write('OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n')
82 out.write('%8d %12.6lf %12.6lf %12.6lf\n' % (natom, 0.0, 0.0, 0.0))
83 out.write('%8d %12.6lf %12.6lf %12.6lf\n' % (nx, hx, 0.0, 0.0))
84 out.write('%8d %12.6lf %12.6lf %12.6lf\n' % (ny, 0.0, hy, 0.0))
85 out.write('%8d %12.6lf %12.6lf %12.6lf\n' % (nz, 0.0, 0.0, hz))
86 for n in range(natom):
87 out.write(atom_lines[n])
88 for i in range(nx):
89 sys.stderr.write('.')
90 for j in range(ny):
91 for k in range(nz):
92 out.write('%18.8lf\n' % data[i, j, k])
93 out.close()
94 sys.stderr.write('\n')
95 sys.stderr.write('Writing done. %s\n' % datetime.datetime.now())
98def isocubes(filename, isoval, scale):
99 """Requires PIL"""
100 import Image
102 f = open(filename, 'r', 256 * 1024)
103 f.readline()
104 f.readline()
106 natom = int(f.readline().split()[0])
108 elems = f.readline().split()
109 nx = int(elems[0])
110 # x0 = float(elems[1])
112 elems = f.readline().split()
113 ny = int(elems[0])
114 # y0 = float(elems[2])
116 elems = f.readline().split()
117 nz = int(elems[0])
118 # z0 = float(elems[3])
120 for i in range(natom):
121 f.readline()
123 # data = np.zeros((nx,ny,nz))
124 dataxy = [[0.0, 0.0] for i in range(nx * ny)]
125 dataxz = [[0.0, 0.0] for i in range(nx * nz)]
126 datayz = [[0.0, 0.0] for i in range(ny * nz)]
127 for i in range(nx):
128 for j in range(ny):
129 for k in range(nz):
130 val = float(f.readline())
131 if val > isoval:
132 dataxy[i * ny + j][0] += val * scale
133 dataxy[i * ny + j][1] += val * scale
134 dataxz[i * nz + k][0] += val * scale
135 dataxz[i * nz + k][1] += val * scale
136 datayz[j * nz + k][0] += val * scale
137 datayz[j * nz + k][1] += val * scale
138 elif -val > isoval:
139 dataxy[i * ny + j][0] += -val * scale
140 dataxy[i * ny + j][1] -= -val * scale
141 dataxz[i * nz + k][0] += -val * scale
142 dataxz[i * nz + k][1] -= -val * scale
143 datayz[j * nz + k][0] += -val * scale
144 datayz[j * nz + k][1] -= -val * scale
146 image_xy = Image.new('RGB', (nx, ny))
147 image_xy_data = image_xy.load()
149 for i in range(nx):
150 for j in range(ny):
151 n = dataxy[i * ny + j][0] * 100
152 d = dataxy[i * ny + j][1] * 100
153 if n > 255:
154 n = 255
155 d = 255 * d / n
156 n = int(n + .5)
157 d = int(d)
158 if d > 0:
159 dataxy[i * ny + j] = (n, n - d, n - d)
160 elif d < 0:
161 dataxy[i * ny + j] = (n + d, n + d, n)
162 else:
163 dataxy[i * ny + j] = (0, 0, 0)
165 image_xy_data[i, j] = dataxy[i * ny + j]
167 image_xy.save(filename + '_iso_xy.png')
169 image_xz = Image.new('RGB', (nx, nz))
170 image_xz_data = image_xz.load()
172 for i in range(nx):
173 for j in range(nz):
174 n = dataxz[i * nz + j][0] * 100
175 d = dataxz[i * nz + j][1] * 100
176 if n > 255:
177 n = 255
178 d = 255 * d / n
179 n = int(n + .5)
180 d = int(d)
181 if d > 0:
182 dataxz[i * nz + j] = (n, n - d, n - d)
183 elif d < 0:
184 dataxz[i * nz + j] = (n + d, n + d, n)
185 else:
186 dataxz[i * nz + j] = (0, 0, 0)
188 image_xz_data[i, j] = dataxz[i * nz + j]
190 image_xz.save(filename + '_iso_xz.png')
192 image_yz = Image.new('RGB', (ny, nz))
193 image_yz_data = image_yz.load()
195 for i in range(ny):
196 for j in range(nz):
197 n = datayz[i * nz + j][0] * 100
198 d = datayz[i * nz + j][1] * 100
199 if n > 255:
200 n = 255
201 d = 255 * d / n
202 n = int(n + .5)
203 d = int(d)
204 if d > 0:
205 datayz[i * nz + j] = (n, n - d, n - d)
206 elif d < 0:
207 datayz[i * nz + j] = (n + d, n + d, n)
208 else:
209 datayz[i * nz + j] = (0, 0, 0)
211 image_yz_data[i, j] = datayz[i * nz + j]
213 image_yz.save(filename + '_iso_yz.png')