Coverage for /builds/kinetik161/ase/ase/io/pov.py: 79.78%
356 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-12-10 11:04 +0000
1"""
2Module for povray file format support.
4See http://www.povray.org/ for details on the format.
5"""
6from collections.abc import Sequence
7from os import unlink
8from pathlib import Path
9from subprocess import DEVNULL, check_call
11import numpy as np
13from ase import Atoms
14from ase.constraints import FixAtoms
15from ase.io.utils import PlottingVariables
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 = {}
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 # These new styles were an attempt to port the old styles o the correct
104 # gamma, many had or still have unphysical light properties inorder to
105 # acheive a certain look.
106 material_styles_dict = dict(
107 simple='finish {phong 0.7 ambient 0.4 diffuse 0.55}',
108 # In general, 'pale' doesn't conserve energy and can look
109 # strange in many cases.
110 pale=('finish {ambient 0.9 diffuse 0.30 roughness 0.001 '
111 'specular 0.2 }'),
112 intermediate=('finish {ambient 0.4 diffuse 0.6 specular 0.1 '
113 'roughness 0.04}'),
114 vmd=(
115 'finish {ambient 0.2 diffuse 0.80 phong 0.25 phong_size 10.0 '
116 'specular 0.2 roughness 0.1}'),
117 jmol=('finish {ambient 0.4 diffuse 0.6 specular 1 roughness 0.001 '
118 'metallic}'),
119 ase2=('finish {ambient 0.2 brilliance 3 diffuse 0.6 metallic '
120 'specular 0.7 roughness 0.04 reflection 0.15}'),
121 ase3=('finish {ambient 0.4 brilliance 2 diffuse 0.6 metallic '
122 'specular 1.0 roughness 0.001 reflection 0.0}'),
123 glass=('finish {ambient 0.4 diffuse 0.35 specular 1.0 '
124 'roughness 0.001}'),
125 glass2=('finish {ambient 0.3 diffuse 0.3 specular 1.0 '
126 'reflection 0.25 roughness 0.001}'),
127 )
129 # These styles were made when assumed_gamma was 1.0 which gives poor color
130 # reproduction, the correct gamma is 2.2 for the sRGB standard.
131 material_styles_dict_old = dict(
132 simple='finish {phong 0.7}',
133 pale=('finish {ambient 0.5 diffuse 0.85 roughness 0.001 '
134 'specular 0.200 }'),
135 intermediate=('finish {ambient 0.3 diffuse 0.6 specular 0.1 '
136 'roughness 0.04}'),
137 vmd=('finish {ambient 0.0 diffuse 0.65 phong 0.1 phong_size 40.0 '
138 'specular 0.5 }'),
139 jmol=('finish {ambient 0.2 diffuse 0.6 specular 1 roughness 0.001 '
140 'metallic}'),
141 ase2=('finish {ambient 0.05 brilliance 3 diffuse 0.6 metallic '
142 'specular 0.7 roughness 0.04 reflection 0.15}'),
143 ase3=('finish {ambient 0.15 brilliance 2 diffuse 0.6 metallic '
144 'specular 1.0 roughness 0.001 reflection 0.0}'),
145 glass=('finish {ambient 0.05 diffuse 0.3 specular 1.0 '
146 'roughness 0.001}'),
147 glass2=('finish {ambient 0.01 diffuse 0.3 specular 1.0 '
148 'reflection 0.25 roughness 0.001}'),
149 )
151 def __init__(self, cell, cell_vertices, positions, diameters, colors,
152 image_width, image_height, constraints=(), isosurfaces=[],
153 display=False, pause=True, transparent=True, canvas_width=None,
154 canvas_height=None, camera_dist=50., image_plane=None,
155 camera_type='orthographic', point_lights=[],
156 area_light=[(2., 3., 40.), 'White', .7, .7, 3, 3],
157 background='White', textures=None, transmittances=None,
158 depth_cueing=False, cue_density=5e-3,
159 celllinewidth=0.05, bondlinewidth=0.10, bondatoms=[],
160 exportconstraints=False):
161 """
162 # x, y is the image plane, z is *out* of the screen
163 cell: ase.cell
164 cell object
165 cell_vertices: 2-d numpy array
166 contains the 8 vertices of the cell, each with three coordinates
167 positions: 2-d numpy array
168 number of atoms length array with three coordinates for positions
169 diameters: 1-d numpy array
170 diameter of atoms (in order with positions)
171 colors: list of str
172 colors of atoms (in order with positions)
173 image_width: float
174 image width in pixels
175 image_height: float
176 image height in pixels
177 constraints: Atoms.constraints
178 constraints to be visualized
179 isosurfaces: list of POVRAYIsosurface
180 composite object to write/render POVRAY isosurfaces
181 display: bool
182 display while rendering
183 pause: bool
184 pause when done rendering (only if display)
185 transparent: bool
186 make background transparent
187 canvas_width: int
188 width of canvas in pixels
189 canvas_height: int
190 height of canvas in pixels
191 camera_dist: float
192 distance from camera to front atom
193 image_plane: float
194 distance from front atom to image plane
195 camera_type: str
196 if 'orthographic' perspective, ultra_wide_angle
197 point_lights: list of 2-element sequences
198 like [[loc1, color1], [loc2, color2],...]
199 area_light: 3-element sequence of location (3-tuple), color (str),
200 width (float), height (float),
201 Nlamps_x (int), Nlamps_y (int)
202 example [(2., 3., 40.), 'White', .7, .7, 3, 3]
203 background: str
204 color specification, e.g., 'White'
205 textures: list of str
206 length of atoms list of texture names
207 transmittances: list of floats
208 length of atoms list of transmittances of the atoms
209 depth_cueing: bool
210 whether or not to use depth cueing a.k.a. fog
211 use with care - adjust the camera_distance to be closer
212 cue_density: float
213 if there is depth_cueing, how dense is it (how dense is the fog)
214 celllinewidth: float
215 radius of the cylinders representing the cell (Ang.)
216 bondlinewidth: float
217 radius of the cylinders representing bonds (Ang.)
218 bondatoms: list of lists (polymorphic)
219 [[atom1, atom2], ... ] pairs of bonding atoms
220 For bond order > 1 = [[atom1, atom2, offset,
221 bond_order, bond_offset],
222 ... ]
223 bond_order: 1, 2, 3 for single, double,
224 and triple bond
225 bond_offset: vector for shifting bonds from
226 original position. Coordinates are
227 in Angstrom unit.
228 exportconstraints: bool
229 honour FixAtoms and mark?"""
231 # attributes from initialization
232 self.area_light = area_light
233 self.background = background
234 self.bondatoms = bondatoms
235 self.bondlinewidth = bondlinewidth
236 self.camera_dist = camera_dist
237 self.camera_type = camera_type
238 self.celllinewidth = celllinewidth
239 self.cue_density = cue_density
240 self.depth_cueing = depth_cueing
241 self.display = display
242 self.exportconstraints = exportconstraints
243 self.isosurfaces = isosurfaces
244 self.pause = pause
245 self.point_lights = point_lights
246 self.textures = textures
247 self.transmittances = transmittances
248 self.transparent = transparent
250 self.image_width = image_width
251 self.image_height = image_height
252 self.colors = colors
253 self.cell = cell
254 self.diameters = diameters
256 # calculations based on passed inputs
258 z0 = positions[:, 2].max()
259 self.offset = (image_width / 2, image_height / 2, z0)
260 self.positions = positions - self.offset
262 if cell_vertices is not None:
263 self.cell_vertices = cell_vertices - self.offset
264 self.cell_vertices.shape = (2, 2, 2, 3)
265 else:
266 self.cell_vertices = None
268 ratio = float(self.image_width) / self.image_height
269 if canvas_width is None:
270 if canvas_height is None:
271 self.canvas_width = min(self.image_width * 15, 640)
272 self.canvas_height = min(self.image_height * 15, 640)
273 else:
274 self.canvas_width = canvas_height * ratio
275 self.canvas_height = canvas_height
276 elif canvas_height is None:
277 self.canvas_width = canvas_width
278 self.canvas_height = self.canvas_width / ratio
279 else:
280 raise RuntimeError("Can't set *both* width and height!")
282 # Distance to image plane from camera
283 if image_plane is None:
284 if self.camera_type == 'orthographic':
285 self.image_plane = 1 - self.camera_dist
286 else:
287 self.image_plane = 0
288 self.image_plane += self.camera_dist
290 self.constrainatoms = []
291 for c in constraints:
292 if isinstance(c, FixAtoms):
293 # self.constrainatoms.extend(c.index) # is this list-like?
294 for n, i in enumerate(c.index):
295 self.constrainatoms += [i]
297 @classmethod
298 def from_PlottingVariables(cls, pvars, **kwargs):
299 cell = pvars.cell
300 cell_vertices = pvars.cell_vertices
301 if 'colors' in kwargs:
302 colors = kwargs.pop('colors')
303 else:
304 colors = pvars.colors
305 diameters = pvars.d
306 image_height = pvars.h
307 image_width = pvars.w
308 positions = pvars.positions
309 constraints = pvars.constraints
310 return cls(cell=cell, cell_vertices=cell_vertices, colors=colors,
311 constraints=constraints, diameters=diameters,
312 image_height=image_height, image_width=image_width,
313 positions=positions, **kwargs)
315 @classmethod
316 def from_atoms(cls, atoms, **kwargs):
317 return cls.from_plotting_variables(
318 PlottingVariables(atoms, scale=1.0), **kwargs)
320 def write_ini(self, path):
321 """Write ini file."""
323 ini_str = f"""\
324Input_File_Name={path.with_suffix('.pov').name}
325Output_to_File=True
326Output_File_Type=N
327Output_Alpha={'on' if self.transparent else 'off'}
328; if you adjust Height, and width, you must preserve the ratio
329; Width / Height = {self.canvas_width/self.canvas_height:f}
330Width={self.canvas_width}
331Height={self.canvas_height}
332Antialias=True
333Antialias_Threshold=0.1
334Display={self.display}
335Display_Gamma=2.2
336Pause_When_Done={self.pause}
337Verbose=False
338"""
339 with open(path, 'w') as fd:
340 fd.write(ini_str)
341 return path
343 def write_pov(self, path):
344 """Write pov file."""
346 point_lights = '\n'.join(f"light_source {{{pa(loc)} {pc(rgb)}}}"
347 for loc, rgb in self.point_lights)
349 area_light = ''
350 if self.area_light is not None:
351 loc, color, width, height, nx, ny = self.area_light
352 area_light += f"""\nlight_source {{{pa(loc)} {pc(color)}
353 area_light <{width:.2f}, 0, 0>, <0, {height:.2f}, 0>, {nx:n}, {ny:n}
354 adaptive 1 jitter}}"""
356 fog = ''
357 if self.depth_cueing and (self.cue_density >= 1e-4):
358 # same way vmd does it
359 if self.cue_density > 1e4:
360 # larger does not make any sense
361 dist = 1e-4
362 else:
363 dist = 1. / self.cue_density
364 fog += f'fog {{fog_type 1 distance {dist:.4f} '\
365 f'color {pc(self.background)}}}'
367 mat_style_keys = (f'#declare {k} = {v}'
368 for k, v in self.material_styles_dict.items())
369 mat_style_keys = '\n'.join(mat_style_keys)
371 # Draw unit cell
372 cell_vertices = ''
373 if self.cell_vertices is not None:
374 for c in range(3):
375 for j in ([0, 0], [1, 0], [1, 1], [0, 1]):
376 p1 = self.cell_vertices[tuple(j[:c]) + (0,) + tuple(j[c:])]
377 p2 = self.cell_vertices[tuple(j[:c]) + (1,) + tuple(j[c:])]
379 distance = np.linalg.norm(p2 - p1)
380 if distance < 1e-12:
381 continue
383 cell_vertices += f'cylinder {{{pa(p1)}, {pa(p2)}, '\
384 'Rcell pigment {Black}}\n'
385 # all strings are f-strings for consistency
386 cell_vertices = cell_vertices.strip('\n')
388 # Draw atoms
389 a = 0
390 atoms = ''
391 for loc, dia, col in zip(self.positions, self.diameters, self.colors):
392 tex = 'ase3'
393 trans = 0.
394 if self.textures is not None:
395 tex = self.textures[a]
396 if self.transmittances is not None:
397 trans = self.transmittances[a]
398 atoms += f'atom({pa(loc)}, {dia/2.:.2f}, {pc(col)}, '\
399 f'{trans}, {tex}) // #{a:n}\n'
400 a += 1
401 atoms = atoms.strip('\n')
403 # Draw atom bonds
404 bondatoms = ''
405 for pair in self.bondatoms:
406 # Make sure that each pair has 4 componets: a, b, offset,
407 # bond_order, bond_offset
408 # a, b: atom index to draw bond
409 # offset: original meaning to make offset for mid-point.
410 # bond_oder: if not supplied, set it to 1 (single bond).
411 # It can be 1, 2, 3, corresponding to single,
412 # double, triple bond
413 # bond_offset: displacement from original bond position.
414 # Default is (bondlinewidth, bondlinewidth, 0)
415 # for bond_order > 1.
416 if len(pair) == 2:
417 a, b = pair
418 offset = (0, 0, 0)
419 bond_order = 1
420 bond_offset = (0, 0, 0)
421 elif len(pair) == 3:
422 a, b, offset = pair
423 bond_order = 1
424 bond_offset = (0, 0, 0)
425 elif len(pair) == 4:
426 a, b, offset, bond_order = pair
427 bond_offset = (self.bondlinewidth, self.bondlinewidth, 0)
428 elif len(pair) > 4:
429 a, b, offset, bond_order, bond_offset = pair
430 else:
431 raise RuntimeError('Each list in bondatom must have at least '
432 '2 entries. Error at %s' % pair)
434 if len(offset) != 3:
435 raise ValueError('offset must have 3 elements. '
436 'Error at %s' % pair)
437 if len(bond_offset) != 3:
438 raise ValueError('bond_offset must have 3 elements. '
439 'Error at %s' % pair)
440 if bond_order not in [0, 1, 2, 3]:
441 raise ValueError('bond_order must be either 0, 1, 2, or 3. '
442 'Error at %s' % pair)
444 # Up to here, we should have all a, b, offset, bond_order,
445 # bond_offset for all bonds.
447 # Rotate bond_offset so that its direction is 90 deg. off the bond
448 # Utilize Atoms object to rotate
449 if bond_order > 1 and np.linalg.norm(bond_offset) > 1.e-9:
450 tmp_atoms = Atoms('H3')
451 tmp_atoms.set_cell(self.cell)
452 tmp_atoms.set_positions([
453 self.positions[a],
454 self.positions[b],
455 self.positions[b] + np.array(bond_offset),
456 ])
457 tmp_atoms.center()
458 tmp_atoms.set_angle(0, 1, 2, 90)
459 bond_offset = tmp_atoms[2].position - tmp_atoms[1].position
461 R = np.dot(offset, self.cell)
462 mida = 0.5 * (self.positions[a] + self.positions[b] + R)
463 midb = 0.5 * (self.positions[a] + self.positions[b] - R)
464 if self.textures is not None:
465 texa = self.textures[a]
466 texb = self.textures[b]
467 else:
468 texa = texb = 'ase3'
470 if self.transmittances is not None:
471 transa = self.transmittances[a]
472 transb = self.transmittances[b]
473 else:
474 transa = transb = 0.
476 # draw bond, according to its bond_order.
477 # bond_order == 0: No bond is plotted
478 # bond_order == 1: use original code
479 # bond_order == 2: draw two bonds, one is shifted by bond_offset/2,
480 # and another is shifted by -bond_offset/2.
481 # bond_order == 3: draw two bonds, one is shifted by bond_offset,
482 # and one is shifted by -bond_offset, and the
483 # other has no shift.
484 # To shift the bond, add the shift to the first two coordinate in
485 # write statement.
487 posa = self.positions[a]
488 posb = self.positions[b]
489 cola = self.colors[a]
490 colb = self.colors[b]
492 if bond_order == 1:
493 draw_tuples = (
494 (posa, mida, cola, transa, texa),
495 (posb, midb, colb, transb, texb))
497 elif bond_order == 2:
498 bs = [x / 2 for x in bond_offset]
499 draw_tuples = (
500 (posa - bs, mida - bs, cola, transa, texa),
501 (posb - bs, midb - bs, colb, transb, texb),
502 (posa + bs, mida + bs, cola, transa, texa),
503 (posb + bs, midb + bs, colb, transb, texb))
505 elif bond_order == 3:
506 bs = bond_offset
507 draw_tuples = (
508 (posa, mida, cola, transa, texa),
509 (posb, midb, colb, transb, texb),
510 (posa + bs, mida + bs, cola, transa, texa),
511 (posb + bs, midb + bs, colb, transb, texb),
512 (posa - bs, mida - bs, cola, transa, texa),
513 (posb - bs, midb - bs, colb, transb, texb))
515 bondatoms += ''.join(f'cylinder {{{pa(p)}, '
516 f'{pa(m)}, Rbond texture{{pigment '
517 f'{{color {pc(c)} '
518 f'transmit {tr}}} finish{{{tx}}}}}}}\n'
519 for p, m, c, tr, tx in
520 draw_tuples)
522 bondatoms = bondatoms.strip('\n')
524 # Draw constraints if requested
525 constraints = ''
526 if self.exportconstraints:
527 for a in self.constrainatoms:
528 dia = self.diameters[a]
529 loc = self.positions[a]
530 trans = 0.0
531 if self.transmittances is not None:
532 trans = self.transmittances[a]
533 constraints += f'constrain({pa(loc)}, {dia/2.:.2f}, Black, '\
534 f'{trans}, {tex}) // #{a:n} \n'
535 constraints = constraints.strip('\n')
537 pov = f"""#version 3.6;
538#include "colors.inc"
539#include "finish.inc"
541global_settings {{assumed_gamma 2.2 max_trace_level 6}}
542background {{{pc(self.background)}{' transmit 1.0' if self.transparent else ''}}}
543camera {{{self.camera_type}
544 right -{self.image_width:.2f}*x up {self.image_height:.2f}*y
545 direction {self.image_plane:.2f}*z
546 location <0,0,{self.camera_dist:.2f}> look_at <0,0,0>}}
547{point_lights}
548{area_light if area_light != '' else '// no area light'}
549{fog if fog != '' else '// no fog'}
550{mat_style_keys}
551#declare Rcell = {self.celllinewidth:.3f};
552#declare Rbond = {self.bondlinewidth:.3f};
554#macro atom(LOC, R, COL, TRANS, FIN)
555 sphere{{LOC, R texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
556#end
557#macro constrain(LOC, R, COL, TRANS FIN)
558union{{torus{{R, Rcell rotate 45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
559 torus{{R, Rcell rotate -45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
560 translate LOC}}
561#end
563{cell_vertices if cell_vertices != '' else '// no cell vertices'}
564{atoms}
565{bondatoms}
566{constraints if constraints != '' else '// no constraints'}
567""" # noqa: E501
569 with open(path, 'w') as fd:
570 fd.write(pov)
572 return path
574 def write(self, pov_path):
575 pov_path = require_pov(pov_path)
576 ini_path = pov_path.with_suffix('.ini')
577 self.write_ini(ini_path)
578 self.write_pov(pov_path)
579 if self.isosurfaces is not None:
580 with open(pov_path, 'a') as fd:
581 for iso in self.isosurfaces:
582 fd.write(iso.format_mesh())
583 return POVRAYInputs(ini_path)
586def require_pov(path):
587 path = Path(path)
588 if path.suffix != '.pov':
589 raise ValueError(f'Expected .pov path, got {path}')
590 return path
593class POVRAYInputs:
594 def __init__(self, path):
595 self.path = path
597 def render(self, povray_executable='povray', stderr=DEVNULL,
598 clean_up=False):
599 cmd = [povray_executable, str(self.path)]
601 check_call(cmd, stderr=stderr)
602 png_path = self.path.with_suffix('.png').absolute()
603 if not png_path.is_file():
604 raise RuntimeError(f'Povray left no output PNG file "{png_path}"')
606 if clean_up:
607 unlink(self.path)
608 unlink(self.path.with_suffix('.pov'))
610 return png_path
613class POVRAYIsosurface:
614 def __init__(self, density_grid, cut_off, cell, cell_origin,
615 closed_edges=False, gradient_ascending=False,
616 color=(0.85, 0.80, 0.25, 0.2), material='ase3'):
617 """
618 density_grid: 3D float ndarray
619 A regular grid on that spans the cell. The first dimension
620 corresponds to the first cell vector and so on.
621 cut_off: float
622 The density value of the isosurface.
623 cell: 2D float ndarray or ASE cell object
624 The 3 vectors which give the cell's repetition
625 cell_origin: 4 float tuple
626 The cell origin as used by POVRAY object
627 closed_edges: bool
628 Setting this will fill in isosurface edges at the cell boundaries.
629 Filling in the edges can help with visualizing
630 highly porous structures.
631 gradient_ascending: bool
632 Lets you pick the area you want to enclose, i.e., should the denser
633 or less dense area be filled in.
634 color: povray color string, float, or float tuple
635 1 float is interpreted as grey scale, a 3 float tuple is rgb,
636 4 float tuple is rgbt, and 5 float tuple is rgbft, where
637 t is transmission fraction and f is filter fraction.
638 Named Povray colors are set in colors.inc
639 (http://wiki.povray.org/content/Reference:Colors.inc)
640 material: string
641 Can be a finish macro defined by POVRAY.material_styles
642 or a full Povray material {...} specification. Using a
643 full material specification willoverride the color parameter.
644 """
646 self.gradient_direction = 'ascent' if gradient_ascending else 'descent'
647 self.color = color
648 self.material = material
649 self.closed_edges = closed_edges
650 self._cut_off = cut_off
652 if self.gradient_direction == 'ascent':
653 cv = 2 * cut_off
654 else:
655 cv = 0
657 if closed_edges:
658 shape_old = density_grid.shape
659 # since well be padding, we need to keep the data at origin
660 cell_origin += -(1.0 / np.array(shape_old)) @ cell
661 density_grid = np.pad(
662 density_grid, pad_width=(
663 1,), mode='constant', constant_values=cv)
664 shape_new = density_grid.shape
665 s = np.array(shape_new) / np.array(shape_old)
666 cell = cell @ np.diag(s)
668 self.cell = cell
669 self.cell_origin = cell_origin
670 self.density_grid = density_grid
671 self.spacing = tuple(1.0 / np.array(self.density_grid.shape))
673 scaled_verts, faces, normals, values = self.compute_mesh(
674 self.density_grid,
675 self.cut_off,
676 self.spacing,
677 self.gradient_direction)
679 # The verts are scaled by default, this is the super easy way of
680 # distributing them in real space but it's easier to do affine
681 # transformations/rotations on a unit cube so I leave it like that
682 # verts = scaled_verts.dot(self.cell)
683 self.verts = scaled_verts
684 self.faces = faces
686 @property
687 def cut_off(self):
688 return self._cut_off
690 @cut_off.setter
691 def cut_off(self, value):
692 raise Exception("Use the set_cut_off method")
694 def set_cut_off(self, value):
695 self._cut_off = value
697 if self.gradient_direction == 'ascent':
698 cv = 2 * self.cut_off
699 else:
700 cv = 0
702 if self.closed_edges:
703 shape_old = self.density_grid.shape
704 # since well be padding, we need to keep the data at origin
705 self.cell_origin += -(1.0 / np.array(shape_old)) @ self.cell
706 self.density_grid = np.pad(
707 self.density_grid, pad_width=(
708 1,), mode='constant', constant_values=cv)
709 shape_new = self.density_grid.shape
710 s = np.array(shape_new) / np.array(shape_old)
711 self.cell = self.cell @ np.diag(s)
713 self.spacing = tuple(1.0 / np.array(self.density_grid.shape))
715 scaled_verts, faces, _, _ = self.compute_mesh(
716 self.density_grid,
717 self.cut_off,
718 self.spacing,
719 self.gradient_direction)
721 self.verts = scaled_verts
722 self.faces = faces
724 @classmethod
725 def from_POVRAY(cls, povray, density_grid, cut_off, **kwargs):
726 return cls(cell=povray.cell,
727 cell_origin=povray.cell_vertices[0, 0, 0],
728 density_grid=density_grid,
729 cut_off=cut_off, **kwargs)
731 @staticmethod
732 def wrapped_triples_section(triple_list,
733 triple_format="<{:f}, {:f}, {:f}>".format,
734 triples_per_line=4):
736 triples = [triple_format(*x) for x in triple_list]
737 n = len(triples)
738 s = ''
739 tpl = triples_per_line
740 c = 0
742 while c < n - tpl:
743 c += tpl
744 s += '\n '
745 s += ', '.join(triples[c - tpl:c])
746 s += '\n '
747 s += ', '.join(triples[c:])
748 return s
750 @staticmethod
751 def compute_mesh(density_grid, cut_off, spacing, gradient_direction):
752 """
754 Import statement is in this method and not file header
755 since few users will use isosurface rendering.
757 Returns scaled_verts, faces, normals, values. See skimage docs.
759 """
761 # marching_cubes name was changed in skimage v0.19
762 try:
763 # New skimage
764 from skimage.measure import marching_cubes
765 except ImportError:
766 # Old skimage (remove at some point)
767 from skimage.measure import \
768 marching_cubes_lewiner as marching_cubes
770 return marching_cubes(
771 density_grid,
772 level=cut_off,
773 spacing=spacing,
774 gradient_direction=gradient_direction,
775 allow_degenerate=False)
777 def format_mesh(self):
778 """Returns a formatted data output for POVRAY files
780 Example:
781 material = '''
782 material { // This material looks like pink jelly
783 texture {
784 pigment { rgbt <0.8, 0.25, 0.25, 0.5> }
785 finish{ diffuse 0.85 ambient 0.99 brilliance 3 specular 0.5 roughness 0.001
786 reflection { 0.05, 0.98 fresnel on exponent 1.5 }
787 conserve_energy
788 }
789 }
790 interior { ior 1.3 }
791 }
792 photons {
793 target
794 refraction on
795 reflection on
796 collect on
797 }'''
798 """ # noqa: E501
800 if self.material in POVRAY.material_styles_dict:
801 material = f"""material {{
802 texture {{
803 pigment {{ {pc(self.color)} }}
804 finish {{ {self.material} }}
805 }}
806 }}"""
807 else:
808 material = self.material
810 # Start writing the mesh2
811 vertex_vectors = self.wrapped_triples_section(
812 triple_list=self.verts,
813 triple_format="<{:f}, {:f}, {:f}>".format,
814 triples_per_line=4)
816 face_indices = self.wrapped_triples_section(
817 triple_list=self.faces,
818 triple_format="<{:n}, {:n}, {:n}>".format,
819 triples_per_line=5)
821 cell = self.cell
822 cell_or = self.cell_origin
823 mesh2 = f"""\n\nmesh2 {{
824 vertex_vectors {{ {len(self.verts):n},
825 {vertex_vectors}
826 }}
827 face_indices {{ {len(self.faces):n},
828 {face_indices}
829 }}
830{material if material != '' else '// no material'}
831 matrix < {cell[0][0]:f}, {cell[0][1]:f}, {cell[0][2]:f},
832 {cell[1][0]:f}, {cell[1][1]:f}, {cell[1][2]:f},
833 {cell[2][0]:f}, {cell[2][1]:f}, {cell[2][2]:f},
834 {cell_or[0]:f}, {cell_or[1]:f}, {cell_or[2]:f}>
835 }}
836 """
837 return mesh2
840def pop_deprecated(dct, name):
841 import warnings
842 if name in dct:
843 del dct[name]
844 warnings.warn(f'The "{name}" keyword of write_pov() is deprecated '
845 'and has no effect; this will raise an error in the '
846 'future.', FutureWarning)
849def write_pov(filename, atoms, *,
850 povray_settings=None, isosurface_data=None,
851 **generic_projection_settings):
853 for name in ['run_povray', 'povray_path', 'stderr', 'extras']:
854 pop_deprecated(generic_projection_settings, name)
856 if povray_settings is None:
857 povray_settings = {}
859 pvars = PlottingVariables(atoms, scale=1.0, **generic_projection_settings)
860 pov_obj = POVRAY.from_PlottingVariables(pvars, **povray_settings)
862 if isosurface_data is None:
863 isosurface_data = []
864 elif not isinstance(isosurface_data, Sequence):
865 isosurface_data = [isosurface_data]
867 isosurfaces = []
868 for isodata in isosurface_data:
869 if isinstance(isodata, POVRAYIsosurface):
870 iso = isodata
871 else:
872 iso = POVRAYIsosurface.from_POVRAY(pov_obj, **isodata)
873 isosurfaces.append(iso)
874 pov_obj.isosurfaces = isosurfaces
876 return pov_obj.write(filename)