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