Coverage for /builds/debichem-team/python-ase/ase/io/pov.py : 68.41%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Module for povray file format support.
4See http://www.povray.org/ for details on the format.
5"""
6from collections.abc import Mapping, Sequence
7from subprocess import check_call, DEVNULL
8from os import unlink
9from pathlib import Path
11import numpy as np
13from ase.io.utils import PlottingVariables
14from ase.constraints import FixAtoms
15from ase import Atoms
18def pa(array):
19 """Povray array syntax"""
20 return '<' + ', '.join(f"{x:>6.2f}" for x in tuple(array)) + '>'
23def pc(array):
24 """Povray color syntax"""
25 if isinstance(array, str):
26 return 'color ' + array
27 if isinstance(array, float):
28 return f'rgb <{array:.2f}>*3'.format(array)
29 l = len(array)
30 if l > 2 and l < 6:
31 return f"rgb{'' if l == 3 else 't' if l == 4 else 'ft'} <" +\
32 ', '.join(f"{x:.2f}" for x in tuple(array)) + '>'
35def get_bondpairs(atoms, radius=1.1):
36 """Get all pairs of bonding atoms
38 Return all pairs of atoms which are closer than radius times the
39 sum of their respective covalent radii. The pairs are returned as
40 tuples::
42 (a, b, (i1, i2, i3))
44 so that atoms a bonds to atom b displaced by the vector::
46 _ _ _
47 i c + i c + i c ,
48 1 1 2 2 3 3
50 where c1, c2 and c3 are the unit cell vectors and i1, i2, i3 are
51 integers."""
53 from ase.data import covalent_radii
54 from ase.neighborlist import NeighborList
55 cutoffs = radius * covalent_radii[atoms.numbers]
56 nl = NeighborList(cutoffs=cutoffs, self_interaction=False)
57 nl.update(atoms)
58 bondpairs = []
59 for a in range(len(atoms)):
60 indices, offsets = nl.get_neighbors(a)
61 bondpairs.extend([(a, a2, offset)
62 for a2, offset in zip(indices, offsets)])
63 return bondpairs
66def set_high_bondorder_pairs(bondpairs, high_bondorder_pairs=None):
67 """Set high bondorder pairs
69 Modify bondpairs list (from get_bondpairs((atoms)) to include high
70 bondorder pairs.
72 Parameters:
73 -----------
74 bondpairs: List of pairs, generated from get_bondpairs(atoms)
75 high_bondorder_pairs: Dictionary of pairs with high bond orders
76 using the following format:
77 { ( a1, b1 ): ( offset1, bond_order1, bond_offset1),
78 ( a2, b2 ): ( offset2, bond_order2, bond_offset2),
79 ...
80 }
81 offset, bond_order, bond_offset are optional.
82 However, if they are provided, the 1st value is
83 offset, 2nd value is bond_order,
84 3rd value is bond_offset """
86 if high_bondorder_pairs is None:
87 high_bondorder_pairs = dict()
88 bondpairs_ = []
89 for pair in bondpairs:
90 (a, b) = (pair[0], pair[1])
91 if (a, b) in high_bondorder_pairs.keys():
92 bondpair = [a, b] + [item for item in high_bondorder_pairs[(a, b)]]
93 bondpairs_.append(bondpair)
94 elif (b, a) in high_bondorder_pairs.keys():
95 bondpair = [a, b] + [item for item in high_bondorder_pairs[(b, a)]]
96 bondpairs_.append(bondpair)
97 else:
98 bondpairs_.append(pair)
99 return bondpairs_
102class POVRAY:
103 material_styles_dict = dict(
104 simple='finish {phong 0.7}',
105 pale=('finish {ambient 0.5 diffuse 0.85 roughness 0.001 '
106 'specular 0.200 }'),
107 intermediate=('finish {ambient 0.3 diffuse 0.6 specular 0.1 '
108 'roughness 0.04}'),
109 vmd=('finish {ambient 0.0 diffuse 0.65 phong 0.1 phong_size 40.0 '
110 'specular 0.5 }'),
111 jmol=('finish {ambient 0.2 diffuse 0.6 specular 1 roughness 0.001 '
112 'metallic}'),
113 ase2=('finish {ambient 0.05 brilliance 3 diffuse 0.6 metallic '
114 'specular 0.7 roughness 0.04 reflection 0.15}'),
115 ase3=('finish {ambient 0.15 brilliance 2 diffuse 0.6 metallic '
116 'specular 1.0 roughness 0.001 reflection 0.0}'),
117 glass=('finish {ambient 0.05 diffuse 0.3 specular 1.0 '
118 'roughness 0.001}'),
119 glass2=('finish {ambient 0.01 diffuse 0.3 specular 1.0 '
120 'reflection 0.25 roughness 0.001}'),
121 )
123 def __init__(self, cell, cell_vertices, positions, diameters, colors,
124 image_width, image_height, constraints=tuple(), isosurfaces=[],
125 display=False, pause=True, transparent=True, canvas_width=None,
126 canvas_height=None, camera_dist=50., image_plane=None,
127 camera_type='orthographic', point_lights=[],
128 area_light=[(2., 3., 40.), 'White', .7, .7, 3, 3],
129 background='White', textures=None, transmittances=None,
130 depth_cueing=False, cue_density=5e-3,
131 celllinewidth=0.05, bondlinewidth=0.10, bondatoms=[],
132 exportconstraints=False):
133 """
134 # x, y is the image plane, z is *out* of the screen
135 cell: ase.cell
136 cell object
137 cell_vertices: 2-d numpy array
138 contains the 8 vertices of the cell, each with three coordinates
139 positions: 2-d numpy array
140 number of atoms length array with three coordinates for positions
141 diameters: 1-d numpy array
142 diameter of atoms (in order with positions)
143 colors: list of str
144 colors of atoms (in order with positions)
145 image_width: float
146 image width in pixels
147 image_height: float
148 image height in pixels
149 constraints: Atoms.constraints
150 constraints to be visualized
151 isosurfaces: list of POVRAYIsosurface
152 composite object to write/render POVRAY isosurfaces
153 display: bool
154 display while rendering
155 pause: bool
156 pause when done rendering (only if display)
157 transparent: bool
158 make background transparent
159 canvas_width: int
160 width of canvas in pixels
161 canvas_height: int
162 height of canvas in pixels
163 camera_dist: float
164 distance from camera to front atom
165 image_plane: float
166 distance from front atom to image plane
167 camera_type: str
168 if 'orthographic' perspective, ultra_wide_angle
169 point_lights: list of 2-element sequences
170 like [[loc1, color1], [loc2, color2],...]
171 area_light: 3-element sequence of location (3-tuple), color (str),
172 width (float), height (float),
173 Nlamps_x (int), Nlamps_y (int)
174 example [(2., 3., 40.), 'White', .7, .7, 3, 3]
175 background: str
176 color specification, e.g., 'White'
177 textures: list of str
178 length of atoms list of texture names
179 transmittances: list of floats
180 length of atoms list of transmittances of the atoms
181 depth_cueing: bool
182 whether or not to use depth cueing a.k.a. fog
183 use with care - adjust the camera_distance to be closer
184 cue_density: float
185 if there is depth_cueing, how dense is it (how dense is the fog)
186 celllinewidth: float
187 radius of the cylinders representing the cell (Ang.)
188 bondlinewidth: float
189 radius of the cylinders representing bonds (Ang.)
190 bondatoms: list of lists (polymorphic)
191 [[atom1, atom2], ... ] pairs of bonding atoms
192 For bond order > 1 = [[atom1, atom2, offset,
193 bond_order, bond_offset],
194 ... ]
195 bond_order: 1, 2, 3 for single, double,
196 and triple bond
197 bond_offset: vector for shifting bonds from
198 original position. Coordinates are
199 in Angstrom unit.
200 exportconstraints: bool
201 honour FixAtoms and mark?"""
203 # attributes from initialization
204 self.area_light = area_light
205 self.background = background
206 self.bondatoms = bondatoms
207 self.bondlinewidth = bondlinewidth
208 self.camera_dist = camera_dist
209 self.camera_type = camera_type
210 self.celllinewidth = celllinewidth
211 self.cue_density = cue_density
212 self.depth_cueing = depth_cueing
213 self.display = display
214 self.exportconstraints = exportconstraints
215 self.isosurfaces = isosurfaces
216 self.pause = pause
217 self.point_lights = point_lights
218 self.textures = textures
219 self.transmittances = transmittances
220 self.transparent = transparent
222 self.image_width = image_width
223 self.image_height = image_height
224 self.colors = colors
225 self.cell = cell
226 self.diameters = diameters
228 # calculations based on passed inputs
230 z0 = positions[:, 2].max()
231 self.offset = (image_width / 2, image_height / 2, z0)
232 self.positions = positions - self.offset
234 if cell_vertices is not None:
235 self.cell_vertices = cell_vertices - self.offset
236 self.cell_vertices.shape = (2, 2, 2, 3)
237 else:
238 self.cell_vertices = None
240 ratio = float(self.image_width) / self.image_height
241 if canvas_width is None:
242 if canvas_height is None:
243 self.canvas_width = min(self.image_width * 15, 640)
244 self.canvas_height = min(self.image_height * 15, 640)
245 else:
246 self.canvas_width = canvas_height * ratio
247 self.canvas_height = canvas_height
248 elif canvas_height is None:
249 self.canvas_width = canvas_width
250 self.canvas_height = self.canvas_width / ratio
251 else:
252 raise RuntimeError("Can't set *both* width and height!")
254 # Distance to image plane from camera
255 if image_plane is None:
256 if self.camera_type == 'orthographic':
257 self.image_plane = 1 - self.camera_dist
258 else:
259 self.image_plane = 0
260 self.image_plane += self.camera_dist
262 self.constrainatoms = []
263 for c in constraints:
264 if isinstance(c, FixAtoms):
265 # self.constrainatoms.extend(c.index) # is this list-like?
266 for n, i in enumerate(c.index):
267 self.constrainatoms += [i]
269 @classmethod
270 def from_PlottingVariables(cls, pvars, **kwargs):
271 cell = pvars.cell
272 cell_vertices = pvars.cell_vertices
273 if 'colors' in kwargs.keys():
274 colors = kwargs.pop('colors')
275 else:
276 colors = pvars.colors
277 diameters = pvars.d
278 image_height = pvars.h
279 image_width = pvars.w
280 positions = pvars.positions
281 constraints = pvars.constraints
282 return cls(cell=cell, cell_vertices=cell_vertices, colors=colors,
283 constraints=constraints, diameters=diameters,
284 image_height=image_height, image_width=image_width,
285 positions=positions, **kwargs)
287 @classmethod
288 def from_atoms(cls, atoms, **kwargs):
289 return cls.from_plotting_variables(
290 PlottingVariables(atoms, scale=1.0), **kwargs)
292 def write_ini(self, path):
293 """Write ini file."""
295 ini_str = f"""\
296Input_File_Name={path.with_suffix('.pov').name}
297Output_to_File=True
298Output_File_Type=N
299Output_Alpha={'on' if self.transparent else 'off'}
300; if you adjust Height, and width, you must preserve the ratio
301; Width / Height = {self.canvas_width/self.canvas_height:f}
302Width={self.canvas_width}
303Height={self.canvas_height}
304Antialias=True
305Antialias_Threshold=0.1
306Display={self.display}
307Pause_When_Done={self.pause}
308Verbose=False
309"""
310 with open(path, 'w') as fd:
311 fd.write(ini_str)
312 return path
314 def write_pov(self, path):
315 """Write pov file."""
317 point_lights = '\n'.join(f"light_source {{{pa(loc)} {pc(rgb)}}}"
318 for loc, rgb in self.point_lights)
320 area_light = ''
321 if self.area_light is not None:
322 loc, color, width, height, nx, ny = self.area_light
323 area_light += f"""\nlight_source {{{pa(loc)} {pc(color)}
324 area_light <{width:.2f}, 0, 0>, <0, {height:.2f}, 0>, {nx:n}, {ny:n}
325 adaptive 1 jitter}}"""
327 fog = ''
328 if self.depth_cueing and (self.cue_density >= 1e-4):
329 # same way vmd does it
330 if self.cue_density > 1e4:
331 # larger does not make any sense
332 dist = 1e-4
333 else:
334 dist = 1. / self.cue_density
335 fog += f'fog {{fog_type 1 distance {dist:.4f} '\
336 f'color {pc(self.background)}}}'
338 mat_style_keys = (f'#declare {k} = {v}'
339 for k, v in self.material_styles_dict.items())
340 mat_style_keys = '\n'.join(mat_style_keys)
342 # Draw unit cell
343 cell_vertices = ''
344 if self.cell_vertices is not None:
345 for c in range(3):
346 for j in ([0, 0], [1, 0], [1, 1], [0, 1]):
347 p1 = self.cell_vertices[tuple(j[:c]) + (0,) + tuple(j[c:])]
348 p2 = self.cell_vertices[tuple(j[:c]) + (1,) + tuple(j[c:])]
350 distance = np.linalg.norm(p2 - p1)
351 if distance < 1e-12:
352 continue
354 cell_vertices += f'cylinder {{{pa(p1)}, {pa(p2)}, '\
355 f'Rcell pigment {{Black}}}}\n'
356 # all strings are f-strings for consistency
357 cell_vertices = cell_vertices.strip('\n')
359 # Draw atoms
360 a = 0
361 atoms = ''
362 for loc, dia, col in zip(self.positions, self.diameters, self.colors):
363 tex = 'ase3'
364 trans = 0.
365 if self.textures is not None:
366 tex = self.textures[a]
367 if self.transmittances is not None:
368 trans = self.transmittances[a]
369 atoms += f'atom({pa(loc)}, {dia/2.:.2f}, {pc(col)}, '\
370 f'{trans}, {tex}) // #{a:n}\n'
371 a += 1
372 atoms = atoms.strip('\n')
374 # Draw atom bonds
375 bondatoms = ''
376 for pair in self.bondatoms:
377 # Make sure that each pair has 4 componets: a, b, offset,
378 # bond_order, bond_offset
379 # a, b: atom index to draw bond
380 # offset: original meaning to make offset for mid-point.
381 # bond_oder: if not supplied, set it to 1 (single bond).
382 # It can be 1, 2, 3, corresponding to single,
383 # double, triple bond
384 # bond_offset: displacement from original bond position.
385 # Default is (bondlinewidth, bondlinewidth, 0)
386 # for bond_order > 1.
387 if len(pair) == 2:
388 a, b = pair
389 offset = (0, 0, 0)
390 bond_order = 1
391 bond_offset = (0, 0, 0)
392 elif len(pair) == 3:
393 a, b, offset = pair
394 bond_order = 1
395 bond_offset = (0, 0, 0)
396 elif len(pair) == 4:
397 a, b, offset, bond_order = pair
398 bond_offset = (self.bondlinewidth, self.bondlinewidth, 0)
399 elif len(pair) > 4:
400 a, b, offset, bond_order, bond_offset = pair
401 else:
402 raise RuntimeError('Each list in bondatom must have at least '
403 '2 entries. Error at %s' % pair)
405 if len(offset) != 3:
406 raise ValueError('offset must have 3 elements. '
407 'Error at %s' % pair)
408 if len(bond_offset) != 3:
409 raise ValueError('bond_offset must have 3 elements. '
410 'Error at %s' % pair)
411 if bond_order not in [0, 1, 2, 3]:
412 raise ValueError('bond_order must be either 0, 1, 2, or 3. '
413 'Error at %s' % pair)
415 # Up to here, we should have all a, b, offset, bond_order,
416 # bond_offset for all bonds.
418 # Rotate bond_offset so that its direction is 90 deg. off the bond
419 # Utilize Atoms object to rotate
420 if bond_order > 1 and np.linalg.norm(bond_offset) > 1.e-9:
421 tmp_atoms = Atoms('H3')
422 tmp_atoms.set_cell(self.cell)
423 tmp_atoms.set_positions([
424 self.positions[a],
425 self.positions[b],
426 self.positions[b] + np.array(bond_offset),
427 ])
428 tmp_atoms.center()
429 tmp_atoms.set_angle(0, 1, 2, 90)
430 bond_offset = tmp_atoms[2].position - tmp_atoms[1].position
432 R = np.dot(offset, self.cell)
433 mida = 0.5 * (self.positions[a] + self.positions[b] + R)
434 midb = 0.5 * (self.positions[a] + self.positions[b] - R)
435 if self.textures is not None:
436 texa = self.textures[a]
437 texb = self.textures[b]
438 else:
439 texa = texb = 'ase3'
441 if self.transmittances is not None:
442 transa = self.transmittances[a]
443 transb = self.transmittances[b]
444 else:
445 transa = transb = 0.
447 # draw bond, according to its bond_order.
448 # bond_order == 0: No bond is plotted
449 # bond_order == 1: use original code
450 # bond_order == 2: draw two bonds, one is shifted by bond_offset/2,
451 # and another is shifted by -bond_offset/2.
452 # bond_order == 3: draw two bonds, one is shifted by bond_offset,
453 # and one is shifted by -bond_offset, and the
454 # other has no shift.
455 # To shift the bond, add the shift to the first two coordinate in
456 # write statement.
458 posa = self.positions[a]
459 posb = self.positions[b]
460 cola = self.colors[a]
461 colb = self.colors[b]
463 if bond_order == 1:
464 draw_tuples = (posa, mida, cola, transa, texa),\
465 (posb, midb, colb, transb, texb)
467 elif bond_order == 2:
468 bs = [x / 2 for x in bond_offset]
469 draw_tuples = (posa - bs, mida - bs, cola, transa, texa),\
470 (posb - bs, midb - bs, colb, transb, texb),\
471 (posa + bs, mida + bs, cola, transa, texa),\
472 (posb + bs, midb + bs, colb, transb, texb)
474 elif bond_order == 3:
475 bs = bond_offset
476 draw_tuples = (posa, mida, cola, transa, texa),\
477 (posb, midb, colb, transb, texb),\
478 (posa + bs, mida + bs, cola, transa, texa),\
479 (posb + bs, midb + bs, colb, transb, texb),\
480 (posa - bs, mida - bs, cola, transa, texa),\
481 (posb - bs, midb - bs, colb, transb, texb)
483 bondatoms += ''.join(f'cylinder {{{pa(p)}, '
484 f'{pa(m)}, Rbond texture{{pigment '
485 f'{{color {pc(c)} '
486 f'transmit {tr}}} finish{{{tx}}}}}}}\n'
487 for p, m, c, tr, tx in
488 draw_tuples)
490 bondatoms = bondatoms.strip('\n')
492 # Draw constraints if requested
493 constraints = ''
494 if self.exportconstraints:
495 for a in self.constrainatoms:
496 dia = self.diameters[a]
497 loc = self.positions[a]
498 trans = 0.0
499 if self.transmittances is not None:
500 trans = self.transmittances[a]
501 constraints += f'constrain({pa(loc)}, {dia/2.:.2f}, Black, '\
502 f'{trans}, {tex}) // #{a:n} \n'
503 constraints = constraints.strip('\n')
505 pov = f"""#include "colors.inc"
506#include "finish.inc"
508global_settings {{assumed_gamma 1 max_trace_level 6}}
509background {{{pc(self.background)}{' transmit 1.0' if self.transparent else ''}}}
510camera {{{self.camera_type}
511 right -{self.image_width:.2f}*x up {self.image_height:.2f}*y
512 direction {self.image_plane:.2f}*z
513 location <0,0,{self.camera_dist:.2f}> look_at <0,0,0>}}
514{point_lights}
515{area_light if area_light != '' else '// no area light'}
516{fog if fog != '' else '// no fog'}
517{mat_style_keys}
518#declare Rcell = {self.celllinewidth:.3f};
519#declare Rbond = {self.bondlinewidth:.3f};
521#macro atom(LOC, R, COL, TRANS, FIN)
522 sphere{{LOC, R texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
523#end
524#macro constrain(LOC, R, COL, TRANS FIN)
525union{{torus{{R, Rcell rotate 45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
526 torus{{R, Rcell rotate -45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
527 translate LOC}}
528#end
530{cell_vertices if cell_vertices != '' else '// no cell vertices'}
531{atoms}
532{bondatoms}
533{constraints if constraints != '' else '// no constraints'}
534""" # noqa: E501
536 with open(path, 'w') as fd:
537 fd.write(pov)
539 return path
541 def write(self, pov_path):
542 pov_path = require_pov(pov_path)
543 ini_path = pov_path.with_suffix('.ini')
544 self.write_ini(ini_path)
545 self.write_pov(pov_path)
546 if self.isosurfaces is not None:
547 with open(pov_path, 'a') as fd:
548 for iso in self.isosurfaces:
549 fd.write(iso.format_mesh())
550 return POVRAYInputs(ini_path)
553def require_pov(path):
554 path = Path(path)
555 if path.suffix != '.pov':
556 raise ValueError(f'Expected .pov path, got {path}')
557 return path
560class POVRAYInputs:
561 def __init__(self, path):
562 self.path = path
564 def render(self, povray_executable='povray', stderr=DEVNULL,
565 clean_up=False):
566 cmd = [povray_executable, str(self.path)]
568 check_call(cmd, stderr=stderr)
569 png_path = self.path.with_suffix('.png').absolute()
570 if not png_path.is_file():
571 raise RuntimeError(f'Povray left no output PNG file "{png_path}"')
573 if clean_up:
574 unlink(self.path)
575 unlink(self.path.with_suffix('.pov'))
577 return png_path
580class POVRAYIsosurface:
581 def __init__(self, density_grid, cut_off, cell, cell_origin,
582 closed_edges=False, gradient_ascending=False,
583 color=(0.85, 0.80, 0.25, 0.2), material='ase3'):
584 """
585 density_grid: 3D float ndarray
586 A regular grid on that spans the cell. The first dimension
587 corresponds to the first cell vector and so on.
588 cut_off: float
589 The density value of the isosurface.
590 cell: 2D float ndarray or ASE cell object
591 The 3 vectors which give the cell's repetition
592 cell_origin: 4 float tuple
593 The cell origin as used by POVRAY object
594 closed_edges: bool
595 Setting this will fill in isosurface edges at the cell boundaries.
596 Filling in the edges can help with visualizing
597 highly porous structures.
598 gradient_ascending: bool
599 Lets you pick the area you want to enclose, i.e., should the denser
600 or less dense area be filled in.
601 color: povray color string, float, or float tuple
602 1 float is interpreted as grey scale, a 3 float tuple is rgb,
603 4 float tuple is rgbt, and 5 float tuple is rgbft, where
604 t is transmission fraction and f is filter fraction.
605 Named Povray colors are set in colors.inc
606 (http://wiki.povray.org/content/Reference:Colors.inc)
607 material: string
608 Can be a finish macro defined by POVRAY.material_styles
609 or a full Povray material {...} specification. Using a
610 full material specification willoverride the color parameter.
611 """
613 self.gradient_direction = 'ascent' if gradient_ascending else 'descent'
614 self.color = color
615 self.material = material
616 self.closed_edges = closed_edges
617 self._cut_off = cut_off
619 if self.gradient_direction == 'ascent':
620 cv = 2 * cut_off
621 else:
622 cv = 0
624 if closed_edges:
625 shape_old = density_grid.shape
626 # since well be padding, we need to keep the data at origin
627 cell_origin += -(1.0 / np.array(shape_old)) @ cell
628 density_grid = np.pad(
629 density_grid, pad_width=(
630 1,), mode='constant', constant_values=cv)
631 shape_new = density_grid.shape
632 s = np.array(shape_new) / np.array(shape_old)
633 cell = cell @ np.diag(s)
635 self.cell = cell
636 self.cell_origin = cell_origin
637 self.density_grid = density_grid
638 self.spacing = tuple(1.0 / np.array(self.density_grid.shape))
640 scaled_verts, faces, normals, values = self.compute_mesh(
641 self.density_grid,
642 self.cut_off,
643 self.spacing,
644 self.gradient_direction)
646 # The verts are scaled by default, this is the super easy way of
647 # distributing them in real space but it's easier to do affine
648 # transformations/rotations on a unit cube so I leave it like that
649 # verts = scaled_verts.dot(self.cell)
650 self.verts = scaled_verts
651 self.faces = faces
653 @property
654 def cut_off(self):
655 return self._cut_off
657 @cut_off.setter
658 def cut_off(self, value):
659 raise Exception("Use the set_cut_off method")
661 def set_cut_off(self, value):
662 self._cut_off = value
664 if self.gradient_direction == 'ascent':
665 cv = 2 * self.cut_off
666 else:
667 cv = 0
669 if self.closed_edges:
670 shape_old = self.density_grid.shape
671 # since well be padding, we need to keep the data at origin
672 self.cell_origin += -(1.0 / np.array(shape_old)) @ self.cell
673 self.density_grid = np.pad(
674 self.density_grid, pad_width=(
675 1,), mode='constant', constant_values=cv)
676 shape_new = self.density_grid.shape
677 s = np.array(shape_new) / np.array(shape_old)
678 self.cell = self.cell @ np.diag(s)
680 self.spacing = tuple(1.0 / np.array(self.density_grid.shape))
682 scaled_verts, faces, _, _ = self.compute_mesh(
683 self.density_grid,
684 self.cut_off,
685 self.spacing,
686 self.gradient_direction)
688 self.verts = scaled_verts
689 self.faces = faces
691 @classmethod
692 def from_POVRAY(cls, povray, density_grid, cut_off, **kwargs):
693 return cls(cell=povray.cell,
694 cell_origin=povray.cell_vertices[0, 0, 0],
695 density_grid=density_grid,
696 cut_off=cut_off, **kwargs)
698 @staticmethod
699 def wrapped_triples_section(triple_list,
700 triple_format="<{:f}, {:f}, {:f}>".format,
701 triples_per_line=4):
703 triples = [triple_format(*x) for x in triple_list]
704 n = len(triples)
705 s = ''
706 tpl = triples_per_line
707 c = 0
709 while c < n - tpl:
710 c += tpl
711 s += '\n '
712 s += ', '.join(triples[c - tpl:c])
713 s += '\n '
714 s += ', '.join(triples[c:])
715 return s
717 @staticmethod
718 def compute_mesh(density_grid, cut_off, spacing, gradient_direction):
719 """
721 Import statement is in this method and not file header
722 since few users will use isosurface rendering.
724 Returns scaled_verts, faces, normals, values. See skimage docs.
726 """
728 from skimage import measure
729 return measure.marching_cubes_lewiner(
730 density_grid,
731 level=cut_off,
732 spacing=spacing,
733 gradient_direction=gradient_direction,
734 allow_degenerate=False)
736 def format_mesh(self):
737 """Returns a formatted data output for POVRAY files
739 Example:
740 material = '''
741 material { // This material looks like pink jelly
742 texture {
743 pigment { rgbt <0.8, 0.25, 0.25, 0.5> }
744 finish{ diffuse 0.85 ambient 0.99 brilliance 3 specular 0.5 roughness 0.001
745 reflection { 0.05, 0.98 fresnel on exponent 1.5 }
746 conserve_energy
747 }
748 }
749 interior { ior 1.3 }
750 }
751 photons {
752 target
753 refraction on
754 reflection on
755 collect on
756 }'''
757 """ # noqa: E501
759 if self.material in POVRAY.material_styles_dict:
760 material = f"""material {{
761 texture {{
762 pigment {{ {pc(self.color)} }}
763 finish {{ {self.material} }}
764 }}
765 }}"""
766 else:
767 material = self.material
769 # Start writing the mesh2
770 vertex_vectors = self.wrapped_triples_section(
771 triple_list=self.verts,
772 triple_format="<{:f}, {:f}, {:f}>".format,
773 triples_per_line=4)
775 face_indices = self.wrapped_triples_section(
776 triple_list=self.faces,
777 triple_format="<{:n}, {:n}, {:n}>".format,
778 triples_per_line=5)
780 cell = self.cell
781 cell_or = self.cell_origin
782 mesh2 = f"""\n\nmesh2 {{
783 vertex_vectors {{ {len(self.verts):n},
784 {vertex_vectors}
785 }}
786 face_indices {{ {len(self.faces):n},
787 {face_indices}
788 }}
789{material if material != '' else '// no material'}
790 matrix < {cell[0][0]:f}, {cell[0][1]:f}, {cell[0][2]:f},
791 {cell[1][0]:f}, {cell[1][1]:f}, {cell[1][2]:f},
792 {cell[2][0]:f}, {cell[2][1]:f}, {cell[2][2]:f},
793 {cell_or[0]:f}, {cell_or[1]:f}, {cell_or[2]:f}>
794 }}
795 """
796 return mesh2
799def pop_deprecated(dct, name):
800 import warnings
801 if name in dct:
802 del dct[name]
803 warnings.warn(f'The "{name}" keyword of write_pov() is deprecated '
804 'and has no effect; this will raise an error in the '
805 'future.', FutureWarning)
808def write_pov(filename, atoms, *,
809 povray_settings=None, isosurface_data=None,
810 **generic_projection_settings):
812 for name in ['run_povray', 'povray_path', 'stderr', 'extras']:
813 pop_deprecated(generic_projection_settings, name)
815 if povray_settings is None:
816 povray_settings = {}
818 pvars = PlottingVariables(atoms, scale=1.0, **generic_projection_settings)
819 pov_obj = POVRAY.from_PlottingVariables(pvars, **povray_settings)
821 if isinstance(isosurface_data, Mapping):
822 pov_obj.isosurfaces = [POVRAYIsosurface.from_POVRAY(
823 pov_obj, **isosurface_data)]
824 elif isinstance(isosurface_data, Sequence):
825 pov_obj.isosurfaces = [POVRAYIsosurface.from_POVRAY(
826 pov_obj, **isodata) for isodata in isosurface_data]
828 return pov_obj.write(filename)