Coverage for /builds/ase/ase/ase/io/pov.py: 79.78%
356 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +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 self.canvas_height = min(self.image_height * 15, 640)
275 else:
276 self.canvas_width = canvas_height * ratio
277 self.canvas_height = canvas_height
278 elif canvas_height is None:
279 self.canvas_width = canvas_width
280 self.canvas_height = self.canvas_width / ratio
281 else:
282 raise RuntimeError("Can't set *both* width and height!")
284 # Distance to image plane from camera
285 if image_plane is None:
286 if self.camera_type == 'orthographic':
287 self.image_plane = 1 - self.camera_dist
288 else:
289 self.image_plane = 0
290 self.image_plane += self.camera_dist
292 self.constrainatoms = []
293 for c in constraints:
294 if isinstance(c, FixAtoms):
295 # self.constrainatoms.extend(c.index) # is this list-like?
296 for n, i in enumerate(c.index):
297 self.constrainatoms += [i]
299 @classmethod
300 def from_PlottingVariables(cls, pvars, **kwargs):
301 cell = pvars.cell
302 cell_vertices = pvars.cell_vertices
303 if 'colors' in kwargs:
304 colors = kwargs.pop('colors')
305 else:
306 colors = pvars.colors
307 diameters = pvars.d
308 image_height = pvars.h
309 image_width = pvars.w
310 positions = pvars.positions
311 constraints = pvars.constraints
312 return cls(cell=cell, cell_vertices=cell_vertices, colors=colors,
313 constraints=constraints, diameters=diameters,
314 image_height=image_height, image_width=image_width,
315 positions=positions, **kwargs)
317 @classmethod
318 def from_atoms(cls, atoms, **kwargs):
319 return cls.from_plotting_variables(
320 PlottingVariables(atoms, scale=1.0), **kwargs)
322 def write_ini(self, path):
323 """Write ini file."""
325 ini_str = f"""\
326Input_File_Name={path.with_suffix('.pov').name}
327Output_to_File=True
328Output_File_Type=N
329Output_Alpha={'on' if self.transparent else 'off'}
330; if you adjust Height, and width, you must preserve the ratio
331; Width / Height = {self.canvas_width / self.canvas_height:f}
332Width={self.canvas_width}
333Height={self.canvas_height}
334Antialias=True
335Antialias_Threshold=0.1
336Display={self.display}
337Display_Gamma=2.2
338Pause_When_Done={self.pause}
339Verbose=False
340"""
341 with open(path, 'w') as fd:
342 fd.write(ini_str)
343 return path
345 def write_pov(self, path):
346 """Write pov file."""
348 point_lights = '\n'.join(f"light_source {{{pa(loc)} {pc(rgb)}}}"
349 for loc, rgb in self.point_lights)
351 area_light = ''
352 if self.area_light is not None:
353 loc, color, width, height, nx, ny = self.area_light
354 area_light += f"""\nlight_source {{{pa(loc)} {pc(color)}
355 area_light <{width:.2f}, 0, 0>, <0, {height:.2f}, 0>, {nx:n}, {ny:n}
356 adaptive 1 jitter}}"""
358 fog = ''
359 if self.depth_cueing and (self.cue_density >= 1e-4):
360 # same way vmd does it
361 if self.cue_density > 1e4:
362 # larger does not make any sense
363 dist = 1e-4
364 else:
365 dist = 1. / self.cue_density
366 fog += f'fog {{fog_type 1 distance {dist:.4f} '\
367 f'color {pc(self.background)}}}'
369 mat_style_keys = (f'#declare {k} = {v}'
370 for k, v in self.material_styles_dict.items())
371 mat_style_keys = '\n'.join(mat_style_keys)
373 # Draw unit cell
374 cell_vertices = ''
375 if self.cell_vertices is not None:
376 for c in range(3):
377 for j in ([0, 0], [1, 0], [1, 1], [0, 1]):
378 p1 = self.cell_vertices[tuple(j[:c]) + (0,) + tuple(j[c:])]
379 p2 = self.cell_vertices[tuple(j[:c]) + (1,) + tuple(j[c:])]
381 distance = np.linalg.norm(p2 - p1)
382 if distance < 1e-12:
383 continue
385 cell_vertices += f'cylinder {{{pa(p1)}, {pa(p2)}, '\
386 'Rcell pigment {Black}}\n'
387 # all strings are f-strings for consistency
388 cell_vertices = cell_vertices.strip('\n')
390 # Draw atoms
391 a = 0
392 atoms = ''
393 for loc, dia, col in zip(self.positions, self.diameters, self.colors):
394 tex = 'ase3'
395 trans = 0.
396 if self.textures is not None:
397 tex = self.textures[a]
398 if self.transmittances is not None:
399 trans = self.transmittances[a]
400 atoms += f'atom({pa(loc)}, {dia / 2.:.2f}, {pc(col)}, '\
401 f'{trans}, {tex}) // #{a:n}\n'
402 a += 1
403 atoms = atoms.strip('\n')
405 # Draw atom bonds
406 bondatoms = ''
407 for pair in self.bondatoms:
408 # Make sure that each pair has 4 componets: a, b, offset,
409 # bond_order, bond_offset
410 # a, b: atom index to draw bond
411 # offset: original meaning to make offset for mid-point.
412 # bond_oder: if not supplied, set it to 1 (single bond).
413 # It can be 1, 2, 3, corresponding to single,
414 # double, triple bond
415 # bond_offset: displacement from original bond position.
416 # Default is (bondlinewidth, bondlinewidth, 0)
417 # for bond_order > 1.
418 if len(pair) == 2:
419 a, b = pair
420 offset = (0, 0, 0)
421 bond_order = 1
422 bond_offset = (0, 0, 0)
423 elif len(pair) == 3:
424 a, b, offset = pair
425 bond_order = 1
426 bond_offset = (0, 0, 0)
427 elif len(pair) == 4:
428 a, b, offset, bond_order = pair
429 bond_offset = (self.bondlinewidth, self.bondlinewidth, 0)
430 elif len(pair) > 4:
431 a, b, offset, bond_order, bond_offset = pair
432 else:
433 raise RuntimeError('Each list in bondatom must have at least '
434 '2 entries. Error at %s' % pair)
436 if len(offset) != 3:
437 raise ValueError('offset must have 3 elements. '
438 'Error at %s' % pair)
439 if len(bond_offset) != 3:
440 raise ValueError('bond_offset must have 3 elements. '
441 'Error at %s' % pair)
442 if bond_order not in [0, 1, 2, 3]:
443 raise ValueError('bond_order must be either 0, 1, 2, or 3. '
444 'Error at %s' % pair)
446 # Up to here, we should have all a, b, offset, bond_order,
447 # bond_offset for all bonds.
449 # Rotate bond_offset so that its direction is 90 deg. off the bond
450 # Utilize Atoms object to rotate
451 if bond_order > 1 and np.linalg.norm(bond_offset) > 1.e-9:
452 tmp_atoms = Atoms('H3')
453 tmp_atoms.set_cell(self.cell)
454 tmp_atoms.set_positions([
455 self.positions[a],
456 self.positions[b],
457 self.positions[b] + np.array(bond_offset),
458 ])
459 tmp_atoms.center()
460 tmp_atoms.set_angle(0, 1, 2, 90)
461 bond_offset = tmp_atoms[2].position - tmp_atoms[1].position
463 R = np.dot(offset, self.cell)
464 mida = 0.5 * (self.positions[a] + self.positions[b] + R)
465 midb = 0.5 * (self.positions[a] + self.positions[b] - R)
466 if self.textures is not None:
467 texa = self.textures[a]
468 texb = self.textures[b]
469 else:
470 texa = texb = 'ase3'
472 if self.transmittances is not None:
473 transa = self.transmittances[a]
474 transb = self.transmittances[b]
475 else:
476 transa = transb = 0.
478 # draw bond, according to its bond_order.
479 # bond_order == 0: No bond is plotted
480 # bond_order == 1: use original code
481 # bond_order == 2: draw two bonds, one is shifted by bond_offset/2,
482 # and another is shifted by -bond_offset/2.
483 # bond_order == 3: draw two bonds, one is shifted by bond_offset,
484 # and one is shifted by -bond_offset, and the
485 # other has no shift.
486 # To shift the bond, add the shift to the first two coordinate in
487 # write statement.
489 posa = self.positions[a]
490 posb = self.positions[b]
491 cola = self.colors[a]
492 colb = self.colors[b]
494 if bond_order == 1:
495 draw_tuples = (
496 (posa, mida, cola, transa, texa),
497 (posb, midb, colb, transb, texb))
499 elif bond_order == 2:
500 bs = [x / 2 for x in bond_offset]
501 draw_tuples = (
502 (posa - bs, mida - bs, cola, transa, texa),
503 (posb - bs, midb - bs, colb, transb, texb),
504 (posa + bs, mida + bs, cola, transa, texa),
505 (posb + bs, midb + bs, colb, transb, texb))
507 elif bond_order == 3:
508 bs = bond_offset
509 draw_tuples = (
510 (posa, mida, cola, transa, texa),
511 (posb, midb, colb, transb, texb),
512 (posa + bs, mida + bs, cola, transa, texa),
513 (posb + bs, midb + bs, colb, transb, texb),
514 (posa - bs, mida - bs, cola, transa, texa),
515 (posb - bs, midb - bs, colb, transb, texb))
517 bondatoms += ''.join(f'cylinder {{{pa(p)}, '
518 f'{pa(m)}, Rbond texture{{pigment '
519 f'{{color {pc(c)} '
520 f'transmit {tr}}} finish{{{tx}}}}}}}\n'
521 for p, m, c, tr, tx in
522 draw_tuples)
524 bondatoms = bondatoms.strip('\n')
526 # Draw constraints if requested
527 constraints = ''
528 if self.exportconstraints:
529 for a in self.constrainatoms:
530 dia = self.diameters[a]
531 loc = self.positions[a]
532 trans = 0.0
533 if self.transmittances is not None:
534 trans = self.transmittances[a]
535 constraints += f'constrain({pa(loc)}, {dia / 2.:.2f}, Black, '\
536 f'{trans}, {tex}) // #{a:n} \n'
537 constraints = constraints.strip('\n')
539 pov = f"""#version 3.6;
540#include "colors.inc"
541#include "finish.inc"
543global_settings {{assumed_gamma 2.2 max_trace_level 6}}
544background {{{pc(self.background)}{' transmit 1.0' if self.transparent else ''}}}
545camera {{{self.camera_type}
546 right -{self.image_width:.2f}*x up {self.image_height:.2f}*y
547 direction {self.image_plane:.2f}*z
548 location <0,0,{self.camera_dist:.2f}> look_at <0,0,0>}}
549{point_lights}
550{area_light if area_light != '' else '// no area light'}
551{fog if fog != '' else '// no fog'}
552{mat_style_keys}
553#declare Rcell = {self.celllinewidth:.3f};
554#declare Rbond = {self.bondlinewidth:.3f};
556#macro atom(LOC, R, COL, TRANS, FIN)
557 sphere{{LOC, R texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
558#end
559#macro constrain(LOC, R, COL, TRANS FIN)
560union{{torus{{R, Rcell rotate 45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
561 torus{{R, Rcell rotate -45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}}
562 translate LOC}}
563#end
565{cell_vertices if cell_vertices != '' else '// no cell vertices'}
566{atoms}
567{bondatoms}
568{constraints if constraints != '' else '// no constraints'}
569""" # noqa: E501
571 with open(path, 'w') as fd:
572 fd.write(pov)
574 return path
576 def write(self, pov_path):
577 pov_path = require_pov(pov_path)
578 ini_path = pov_path.with_suffix('.ini')
579 self.write_ini(ini_path)
580 self.write_pov(pov_path)
581 if self.isosurfaces is not None:
582 with open(pov_path, 'a') as fd:
583 for iso in self.isosurfaces:
584 fd.write(iso.format_mesh())
585 return POVRAYInputs(ini_path)
588def require_pov(path):
589 path = Path(path)
590 if path.suffix != '.pov':
591 raise ValueError(f'Expected .pov path, got {path}')
592 return path
595class POVRAYInputs:
596 def __init__(self, path):
597 self.path = path
599 def render(self, povray_executable='povray', stderr=DEVNULL,
600 clean_up=False):
601 cmd = [povray_executable, str(self.path)]
603 check_call(cmd, stderr=stderr)
604 png_path = self.path.with_suffix('.png').absolute()
605 if not png_path.is_file():
606 raise RuntimeError(f'Povray left no output PNG file "{png_path}"')
608 if clean_up:
609 unlink(self.path)
610 unlink(self.path.with_suffix('.pov'))
612 return png_path
615class POVRAYIsosurface:
616 def __init__(self, density_grid, cut_off, cell, cell_origin,
617 closed_edges=False, gradient_ascending=False,
618 color=(0.85, 0.80, 0.25, 0.2), material='ase3'):
619 """
620 density_grid: 3D float ndarray
621 A regular grid on that spans the cell. The first dimension
622 corresponds to the first cell vector and so on.
623 cut_off: float
624 The density value of the isosurface.
625 cell: 2D float ndarray or ASE cell object
626 The 3 vectors which give the cell's repetition
627 cell_origin: 4 float tuple
628 The cell origin as used by POVRAY object
629 closed_edges: bool
630 Setting this will fill in isosurface edges at the cell boundaries.
631 Filling in the edges can help with visualizing
632 highly porous structures.
633 gradient_ascending: bool
634 Lets you pick the area you want to enclose, i.e., should the denser
635 or less dense area be filled in.
636 color: povray color string, float, or float tuple
637 1 float is interpreted as grey scale, a 3 float tuple is rgb,
638 4 float tuple is rgbt, and 5 float tuple is rgbft, where
639 t is transmission fraction and f is filter fraction.
640 Named Povray colors are set in colors.inc
641 (http://wiki.povray.org/content/Reference:Colors.inc)
642 material: string
643 Can be a finish macro defined by POVRAY.material_styles
644 or a full Povray material {...} specification. Using a
645 full material specification willoverride the color parameter.
646 """
648 self.gradient_direction = 'ascent' if gradient_ascending else 'descent'
649 self.color = color
650 self.material = material
651 self.closed_edges = closed_edges
652 self._cut_off = cut_off
654 if self.gradient_direction == 'ascent':
655 cv = 2 * cut_off
656 else:
657 cv = 0
659 if closed_edges:
660 shape_old = density_grid.shape
661 # since well be padding, we need to keep the data at origin
662 cell_origin += -(1.0 / np.array(shape_old)) @ cell
663 density_grid = np.pad(
664 density_grid, pad_width=(
665 1,), mode='constant', constant_values=cv)
666 shape_new = density_grid.shape
667 s = np.array(shape_new) / np.array(shape_old)
668 cell = cell @ np.diag(s)
670 self.cell = cell
671 self.cell_origin = cell_origin
672 self.density_grid = density_grid
673 self.spacing = tuple(1.0 / np.array(self.density_grid.shape))
675 scaled_verts, faces, _normals, _values = self.compute_mesh(
676 self.density_grid,
677 self.cut_off,
678 self.spacing,
679 self.gradient_direction)
681 # The verts are scaled by default, this is the super easy way of
682 # distributing them in real space but it's easier to do affine
683 # transformations/rotations on a unit cube so I leave it like that
684 # verts = scaled_verts.dot(self.cell)
685 self.verts = scaled_verts
686 self.faces = faces
688 @property
689 def cut_off(self):
690 return self._cut_off
692 @cut_off.setter
693 def cut_off(self, value):
694 raise Exception("Use the set_cut_off method")
696 def set_cut_off(self, value):
697 self._cut_off = value
699 if self.gradient_direction == 'ascent':
700 cv = 2 * self.cut_off
701 else:
702 cv = 0
704 if self.closed_edges:
705 shape_old = self.density_grid.shape
706 # since well be padding, we need to keep the data at origin
707 self.cell_origin += -(1.0 / np.array(shape_old)) @ self.cell
708 self.density_grid = np.pad(
709 self.density_grid, pad_width=(
710 1,), mode='constant', constant_values=cv)
711 shape_new = self.density_grid.shape
712 s = np.array(shape_new) / np.array(shape_old)
713 self.cell = self.cell @ np.diag(s)
715 self.spacing = tuple(1.0 / np.array(self.density_grid.shape))
717 scaled_verts, faces, _, _ = self.compute_mesh(
718 self.density_grid,
719 self.cut_off,
720 self.spacing,
721 self.gradient_direction)
723 self.verts = scaled_verts
724 self.faces = faces
726 @classmethod
727 def from_POVRAY(cls, povray, density_grid, cut_off, **kwargs):
728 return cls(cell=povray.cell,
729 cell_origin=povray.cell_vertices[0, 0, 0],
730 density_grid=density_grid,
731 cut_off=cut_off, **kwargs)
733 @staticmethod
734 def wrapped_triples_section(triple_list,
735 triple_format="<{:f}, {:f}, {:f}>".format,
736 triples_per_line=4):
738 triples = [triple_format(*x) for x in triple_list]
739 n = len(triples)
740 s = ''
741 tpl = triples_per_line
742 c = 0
744 while c < n - tpl:
745 c += tpl
746 s += '\n '
747 s += ', '.join(triples[c - tpl:c])
748 s += '\n '
749 s += ', '.join(triples[c:])
750 return s
752 @staticmethod
753 def compute_mesh(density_grid, cut_off, spacing, gradient_direction):
754 """
756 Import statement is in this method and not file header
757 since few users will use isosurface rendering.
759 Returns scaled_verts, faces, normals, values. See skimage docs.
761 """
763 # marching_cubes name was changed in skimage v0.19
764 try:
765 # New skimage
766 from skimage.measure import marching_cubes
767 except ImportError:
768 # Old skimage (remove at some point)
769 from skimage.measure import marching_cubes_lewiner as marching_cubes
771 return marching_cubes(
772 density_grid,
773 level=cut_off,
774 spacing=spacing,
775 gradient_direction=gradient_direction,
776 allow_degenerate=False)
778 def format_mesh(self):
779 """Returns a formatted data output for POVRAY files
781 Example:
782 material = '''
783 material { // This material looks like pink jelly
784 texture {
785 pigment { rgbt <0.8, 0.25, 0.25, 0.5> }
786 finish{ diffuse 0.85 ambient 0.99 brilliance 3 specular 0.5 roughness 0.001
787 reflection { 0.05, 0.98 fresnel on exponent 1.5 }
788 conserve_energy
789 }
790 }
791 interior { ior 1.3 }
792 }
793 photons {
794 target
795 refraction on
796 reflection on
797 collect on
798 }'''
799 """ # noqa: E501
801 if self.material in POVRAY.material_styles_dict:
802 material = f"""material {{
803 texture {{
804 pigment {{ {pc(self.color)} }}
805 finish {{ {self.material} }}
806 }}
807 }}"""
808 else:
809 material = self.material
811 # Start writing the mesh2
812 vertex_vectors = self.wrapped_triples_section(
813 triple_list=self.verts,
814 triple_format="<{:f}, {:f}, {:f}>".format,
815 triples_per_line=4)
817 face_indices = self.wrapped_triples_section(
818 triple_list=self.faces,
819 triple_format="<{:n}, {:n}, {:n}>".format,
820 triples_per_line=5)
822 cell = self.cell
823 cell_or = self.cell_origin
824 mesh2 = f"""\n\nmesh2 {{
825 vertex_vectors {{ {len(self.verts):n},
826 {vertex_vectors}
827 }}
828 face_indices {{ {len(self.faces):n},
829 {face_indices}
830 }}
831{material if material != '' else '// no material'}
832 matrix < {cell[0][0]:f}, {cell[0][1]:f}, {cell[0][2]:f},
833 {cell[1][0]:f}, {cell[1][1]:f}, {cell[1][2]:f},
834 {cell[2][0]:f}, {cell[2][1]:f}, {cell[2][2]:f},
835 {cell_or[0]:f}, {cell_or[1]:f}, {cell_or[2]:f}>
836 }}
837 """
838 return mesh2
841def pop_deprecated(dct, name):
842 import warnings
843 if name in dct:
844 del dct[name]
845 warnings.warn(f'The "{name}" keyword of write_pov() is deprecated '
846 'and has no effect; this will raise an error in the '
847 'future.', FutureWarning)
850def write_pov(filename, atoms, *,
851 povray_settings=None, isosurface_data=None,
852 **generic_projection_settings):
854 for name in ['run_povray', 'povray_path', 'stderr', 'extras']:
855 pop_deprecated(generic_projection_settings, name)
857 if povray_settings is None:
858 povray_settings = {}
860 pvars = PlottingVariables(atoms, scale=1.0, **generic_projection_settings)
861 pov_obj = POVRAY.from_PlottingVariables(pvars, **povray_settings)
863 if isosurface_data is None:
864 isosurface_data = []
865 elif not isinstance(isosurface_data, Sequence):
866 isosurface_data = [isosurface_data]
868 isosurfaces = []
869 for isodata in isosurface_data:
870 if isinstance(isodata, POVRAYIsosurface):
871 iso = isodata
872 else:
873 iso = POVRAYIsosurface.from_POVRAY(pov_obj, **isodata)
874 isosurfaces.append(iso)
875 pov_obj.isosurfaces = isosurfaces
877 return pov_obj.write(filename)