/
generic.py
396 lines (330 loc) · 16.1 KB
/
generic.py
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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
# coding: utf-8
# Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
# Distributed under the terms of "New BSD License", see the LICENSE file.
import numpy as np
from pyiron.atomistics.structure.atoms import Atoms
from pyiron.vasp.structure import write_poscar
__author__ = "Sudarsan Surendralal, Su-Hyun Yoo"
__copyright__ = (
"Copyright 2020, Max-Planck-Institut für Eisenforschung GmbH "
"- Computational Materials Design (CM) Department"
)
__version__ = "1.0"
__maintainer__ = "Sudarsan Surendralal"
__email__ = "surendralal@mpie.de"
__status__ = "development"
__date__ = "Sep 1, 2017"
class VolumetricData(object):
"""
A new class to handle 3-dimensional volumetric data elegantly (charge densities, electrostatic potentials etc) based
on the numpy.ndarray instance. This module is adapted from the pymatgen vasp VolumtricData class
http://pymatgen.org/_modules/pymatgen/io/vasp/outputs.html#VolumetricData
Attributes:
total_data (numpy.ndarray): A 3D array containing the data
"""
def __init__(self):
self._total_data = None
self._atoms = None
@property
def atoms(self):
"""
The structure related to the volumeric data
Returns:
pyiron.atomistics.structure.Atoms: The structure associated with the data
"""
return self._atoms
@atoms.setter
def atoms(self, val):
self._atoms = val
@property
def total_data(self):
"""
numpy.ndarray: The Nx x Ny x Nz sized array for the total data
"""
return self._total_data
@total_data.setter
def total_data(self, val):
if not (isinstance(val, (np.ndarray, list))):
raise TypeError(
"Attribute total_data should be a numpy.ndarray instance or a list and "
"not {}".format(type(val))
)
val = np.array(val)
shape = np.array(np.shape(val))
if not (len(shape) == 3):
raise ValueError("Attribute total_data should be a 3D array")
self._total_data = val
@staticmethod
def gauss_f(d, fwhm=0.529177):
"""
Generates a Gaussian distribution for a given distance and full width half maximum value
Args:
d (float): distance between target point and reference point
fwhm (float): Full width half maximum in angstrom
Returns:
float: Gaussian reduction constant
"""
sigma = fwhm / (2 * np.sqrt(2 * np.log(2)))
d2 = d * d
return np.exp(-1 / (2 * sigma ** 2) * d2)
@staticmethod
def dist_between_two_grid_points(target_grid_point, n_grid_at_center, lattice, grid_shape):
"""
Calculates the distance between a target grid point and another grid point
Args:
target_grid_point (numpy.ndarray/list): Target grid point
n_grid_at_center (numpy.ndarray/list): coordinate of center of sphere
lattice (numpy.ndarray/list): lattice vector
grid_shape (tuple/list/numpy.ndarray): size of grid
Returns:
float: Distance between target grid and center of sphere in angstrom
"""
unit_dist_in_grid = [np.sqrt(np.dot(lattice[0], lattice[0])) / grid_shape[0],
np.sqrt(np.dot(lattice[1], lattice[1])) / grid_shape[1],
np.sqrt(np.dot(lattice[2], lattice[2])) / grid_shape[2]]
dn = np.multiply(np.subtract(target_grid_point, n_grid_at_center), unit_dist_in_grid)
dist = np.linalg.norm(dn)
return dist
def spherical_average_potential(self, structure, spherical_center, rad=2, fwhm=0.529177):
"""
Calculates the spherical average about a given point in space
Args:
structure (pyiron.atomistics.structure.Atoms): Input structure
spherical_center (list/numpy.ndarray): position of spherical_center in direct coordinate
rad (float): radius of sphere to be considered in Angstrom (recommended value: 2)
fwhm (float): Full width half maximum of gaussian function in Angstrom (recommended value: 0.529177)
Returns:
float: Spherical average at the target center
"""
grid_shape = self._total_data.shape
# Position of center of sphere at grid coordinates
n_grid_at_center = [int(np.ceil(spherical_center[0] * grid_shape[0])),
int(np.ceil(spherical_center[1] * grid_shape[1])),
int(np.ceil(spherical_center[2] * grid_shape[2]))]
# Unit distance between grids
dist_in_grid = [np.linalg.norm(structure.cell[0]) / grid_shape[0],
np.linalg.norm(structure.cell[1]) / grid_shape[1],
np.linalg.norm(structure.cell[2]) / grid_shape[2]]
# Range of grids to be considered within the provided radius w.r.t. center of sphere
num_grid_in_sph = [[], []]
for i, dist in enumerate(dist_in_grid):
num_grid_in_sph[0].append(n_grid_at_center[i] - int(np.ceil(rad / dist)))
num_grid_in_sph[1].append(n_grid_at_center[i] + int(np.ceil(rad / dist)))
sph_avg_tmp = []
weight = 0
for k in range(num_grid_in_sph[0][0], num_grid_in_sph[1][0]):
for l in range(num_grid_in_sph[0][1], num_grid_in_sph[1][1]):
for m in range(num_grid_in_sph[0][2], num_grid_in_sph[1][2]):
target_grid_point = [k, l, m]
dist = self.dist_between_two_grid_points(target_grid_point,
n_grid_at_center, structure.cell, grid_shape)
if dist <= rad:
sph_avg_tmp.append(
self._total_data[k % grid_shape[0], l % grid_shape[1], m % grid_shape[2]]
* self.gauss_f(dist, fwhm))
weight += self.gauss_f(dist, fwhm)
else:
pass
sum_list = np.sum(sph_avg_tmp)
sph_avg = sum_list / weight
return sph_avg
@staticmethod
def dist_between_two_grid_points_cyl(target_grid_point, n_grid_at_center, lattice, grid_shape, direction_of_cyl):
"""
Distance between a target grid point and the center of a cylinder
Args:
target_grid_point (numpy.ndarray/list): Target grid point
n_grid_at_center (numpy.ndarray/list): coordinate of center of sphere
lattice (numpy.ndarray/list): lattice vector
grid_shape (tuple/list/numpy.ndarray): size of grid
direction_of_cyl (int): Axis of cylinder (0 (x) or 1 (y) or 2 (z))
Returns:
float: Distance between target grid and in-plane center of cylinder
"""
unit_dist_in_grid = [np.sqrt(np.dot(lattice[0], lattice[0])) / grid_shape[0],
np.sqrt(np.dot(lattice[1], lattice[1])) / grid_shape[1],
np.sqrt(np.dot(lattice[2], lattice[2])) / grid_shape[2]]
dn = np.multiply(np.subtract(target_grid_point, n_grid_at_center), unit_dist_in_grid)
if direction_of_cyl == 0:
dn[0] = 0
elif direction_of_cyl == 1:
dn[1] = 0
elif direction_of_cyl == 2:
dn[2] = 0
else:
print("check the direction of cylindrical axis")
dist = np.linalg.norm(dn)
return dist
def cylindrical_average_potential(self, structure, spherical_center, axis_of_cyl, rad=2, fwhm=0.529177):
"""
Calculates the cylindrical average about a given point in space
Args:
structure (pyiron.atomistics.structure.Atoms): Input structure
spherical_center (list/numpy.ndarray): position of spherical_center in direct coordinate
rad (float): radius of sphere to be considered in Angstrom (recommended value: 2)
fwhm (float): Full width half maximum of gaussian function in Angstrom (recommended value: 0.529177)
axis_of_cyl (int): Axis of cylinder (0 (x) or 1 (y) or 2 (z))
Returns:
float: Cylindrical average at the target center
"""
grid_shape = self._total_data.shape
# Position of center of sphere at grid coordinates
n_grid_at_center = [int(np.ceil(spherical_center[0] * grid_shape[0])),
int(np.ceil(spherical_center[1] * grid_shape[1])),
int(np.ceil(spherical_center[2] * grid_shape[2]))]
# Unit distance between grids
dist_in_grid = [np.linalg.norm(structure.cell[0]) / grid_shape[0],
np.linalg.norm(structure.cell[1]) / grid_shape[1],
np.linalg.norm(structure.cell[2]) / grid_shape[2]]
# Range of grids to be considered within the provided radius w.r.t. center of sphere
num_grid_in_cyl = [[], []]
for i, dist in enumerate(dist_in_grid):
if i == axis_of_cyl:
num_grid_in_cyl[0].append(0)
num_grid_in_cyl[1].append(grid_shape[i])
else:
num_grid_in_cyl[0].append(n_grid_at_center[i] - int(np.ceil(rad / dist)))
num_grid_in_cyl[1].append(n_grid_at_center[i] + int(np.ceil(rad / dist)))
cyl_avg_tmp = []
weight = 0
for k in range(num_grid_in_cyl[0][0], num_grid_in_cyl[1][0]):
for l in range(num_grid_in_cyl[0][1], num_grid_in_cyl[1][1]):
for m in range(num_grid_in_cyl[0][2], num_grid_in_cyl[1][2]):
target_grid_point = [k, l, m]
dist = self.dist_between_two_grid_points_cyl(target_grid_point, n_grid_at_center, structure.cell,
grid_shape, axis_of_cyl)
if dist <= rad:
cyl_avg_tmp.append(
self._total_data[k % grid_shape[0], l % grid_shape[1], m % grid_shape[2]]
* self.gauss_f(dist, fwhm))
weight += self.gauss_f(dist, fwhm)
else:
pass
sum_list = np.sum(cyl_avg_tmp)
cyl_avg = sum_list / weight
return cyl_avg
def get_average_along_axis(self, ind=2):
"""
Get the lateral average along a certain axis direction. This function is adapted from the pymatgen vasp
VolumetricData class
http://pymatgen.org/_modules/pymatgen/io/vasp/outputs.html#VolumetricData.get_average_along_axis
Args:
ind (int): Index of axis (0, 1 and 2 for the x, y, and z axis respectively)
Returns:
numpy.ndarray: A 1D vector with the laterally averaged values of the volumetric data
"""
if ind == 0:
return np.average(np.average(self._total_data, axis=1), 1)
elif ind == 1:
return np.average(np.average(self._total_data, axis=0), 1)
else:
return np.average(np.average(self._total_data, axis=0), 0)
def to_hdf(self, hdf5, group_name="volumetric_data"):
"""
Writes the data as a group to a HDF5 file
Args:
hdf5 (pyiron_base.generic.hdfio.ProjectHDFio): The HDF file/path to write the data to
group_name (str): The name of the group under which the data must be stored as
"""
with hdf5.open(group_name) as hdf_vd:
hdf_vd["TYPE"] = str(type(self))
hdf_vd["total"] = self.total_data
def from_hdf(self, hdf5, group_name="volumetric_data"):
"""
Recreating the VolumetricData instance by reading data from the HDF5 files
Args:
hdf5 (pyiron_base.generic.hdfio.ProjectHDFio): The HDF file/path to write the data to
group_name (str): The name of the group under which the data must be stored as
Returns:
pyiron.atomistics.volumetric.generic.VolumetricData: The VolumetricData instance
"""
with hdf5.open(group_name) as hdf_vd:
self._total_data = hdf_vd["total"]
def write_cube_file(self, filename="cube_file.cube", cell_scaling=1.0):
"""
Write the volumetric data into the CUBE file format
Args:
filename (str): Filename
cell_scaling (float): Scale the cell by this fraction
"""
if self._atoms is None:
raise ValueError(
"The volumetric data object must have a valid structure assigned to it before writing "
"to the cube format"
)
data = self.total_data
n_x, n_y, _ = data.shape
origin = np.zeros(3)
flattened_data = np.hstack(
[data[i, j, :] for i in range(n_x) for j in range(n_y)]
)
n_atoms = len(self.atoms)
total_lines = int(len(flattened_data) / 6) * 6
reshaped_data = np.reshape(flattened_data[0:total_lines], (-1, 6))
last_line = [flattened_data[total_lines:]]
head_array = np.zeros((4, 4))
head_array[0] = np.append([n_atoms], origin)
head_array[1:, 0] = data.shape
head_array[1:, 1:] = self.atoms.cell / data.shape * cell_scaling
position_array = np.zeros((len(self.atoms.positions), 5))
position_array[:, 0] = self.atoms.get_atomic_numbers()
position_array[:, 2:] = self.atoms.positions
with open(filename, "w") as f:
f.write("Cube file generated by pyiron (http://pyiron.org) \n")
f.write("z is the fastest index \n")
np.savetxt(f, head_array, fmt="%4d %.6f %.6f %.6f")
np.savetxt(f, position_array, fmt="%4d %.6f %.6f %.6f %.6f")
np.savetxt(f, reshaped_data, fmt="%.5e")
np.savetxt(f, last_line, fmt="%.5e")
def read_cube_file(self, filename="cube_file.cube"):
"""
Generate data from a CUBE file
Args:
filename (str): Filename to parse
"""
with open(filename, "r") as f:
lines = f.readlines()
n_atoms = int(lines[2].strip().split()[0])
cell_data = np.genfromtxt(lines[3:6])
cell_grid = cell_data[:, 1:]
grid_shape = np.array(cell_data[:, 0], dtype=int)
# total_data = np.zeros(grid_shape)
cell = np.array([val * grid_shape[i] for i, val in enumerate(cell_grid)])
pos_data = np.genfromtxt(lines[6: n_atoms + 6])
if n_atoms == 1:
pos_data = np.array([pos_data])
atomic_numbers = np.array(pos_data[:, 0], dtype=int)
positions = pos_data[:, 2:]
self._atoms = Atoms(numbers=atomic_numbers, positions=positions, cell=cell)
end_int = n_atoms + 6 + int(np.prod(grid_shape) / 6)
data = np.genfromtxt(lines[n_atoms+6: end_int])
data_flatten = np.hstack(data)
if np.prod(grid_shape) % 6 > 0:
data_flatten = np.append(
data_flatten, [float(val) for val in lines[end_int].split()]
)
n_x, n_y, n_z = grid_shape
self._total_data = data_flatten.reshape((n_x, n_y, n_z))
def write_vasp_volumetric(self, filename="CHGCAR", normalize=False):
"""
Writes volumetric data into a VASP CHGCAR format
Args:
filename (str): Filename of the new file
normalize (bool): True if the data is to be normalized by the volume
"""
write_poscar(structure=self.atoms, filename=filename)
with open(filename, "a") as f:
f.write("\n")
f.write(" ".join(list(np.array(self.total_data.shape, dtype=str))))
f.write("\n")
_, n_y, n_z = self.total_data.shape
flattened_data = np.hstack(
[self.total_data[:, i, j] for j in range(n_z) for i in range(n_y)]
)
if normalize:
flattened_data /= self.atoms.get_volume()
num_lines = int(len(flattened_data) / 5) * 5
reshaped_data = np.reshape(flattened_data[0:num_lines], (-1, 5))
np.savetxt(f, reshaped_data, fmt="%.12f")
if len(flattened_data) % 5 > 0:
np.savetxt(f, [flattened_data[num_lines:]], fmt="%.12f")