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

1# fmt: off 

2 

3""" 

4Module for povray file format support. 

5 

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 

12 

13import numpy as np 

14 

15from ase import Atoms 

16from ase.constraints import FixAtoms 

17from ase.io.utils import PlottingVariables 

18 

19 

20def pa(array): 

21 """Povray array syntax""" 

22 return '<' + ', '.join(f"{x:>6.2f}" for x in tuple(array)) + '>' 

23 

24 

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)) + '>' 

35 

36 

37def get_bondpairs(atoms, radius=1.1): 

38 """Get all pairs of bonding atoms 

39 

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:: 

43 

44 (a, b, (i1, i2, i3)) 

45 

46 so that atoms a bonds to atom b displaced by the vector:: 

47 

48 _ _ _ 

49 i c + i c + i c , 

50 1 1 2 2 3 3 

51 

52 where c1, c2 and c3 are the unit cell vectors and i1, i2, i3 are 

53 integers.""" 

54 

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 

66 

67 

68def set_high_bondorder_pairs(bondpairs, high_bondorder_pairs=None): 

69 """Set high bondorder pairs 

70 

71 Modify bondpairs list (from get_bondpairs((atoms)) to include high 

72 bondorder pairs. 

73 

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 """ 

87 

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_ 

102 

103 

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 ) 

130 

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 ) 

152 

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?""" 

232 

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 

251 

252 self.image_width = image_width 

253 self.image_height = image_height 

254 self.colors = colors 

255 self.cell = cell 

256 self.diameters = diameters 

257 

258 # calculations based on passed inputs 

259 

260 z0 = positions[:, 2].max() 

261 self.offset = (image_width / 2, image_height / 2, z0) 

262 self.positions = positions - self.offset 

263 

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 

269 

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!") 

283 

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 

291 

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] 

298 

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) 

316 

317 @classmethod 

318 def from_atoms(cls, atoms, **kwargs): 

319 return cls.from_plotting_variables( 

320 PlottingVariables(atoms, scale=1.0), **kwargs) 

321 

322 def write_ini(self, path): 

323 """Write ini file.""" 

324 

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 

344 

345 def write_pov(self, path): 

346 """Write pov file.""" 

347 

348 point_lights = '\n'.join(f"light_source {{{pa(loc)} {pc(rgb)}}}" 

349 for loc, rgb in self.point_lights) 

350 

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}}""" 

357 

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)}}}' 

368 

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) 

372 

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:])] 

380 

381 distance = np.linalg.norm(p2 - p1) 

382 if distance < 1e-12: 

383 continue 

384 

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') 

389 

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') 

404 

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) 

435 

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) 

445 

446 # Up to here, we should have all a, b, offset, bond_order, 

447 # bond_offset for all bonds. 

448 

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 

462 

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' 

471 

472 if self.transmittances is not None: 

473 transa = self.transmittances[a] 

474 transb = self.transmittances[b] 

475 else: 

476 transa = transb = 0. 

477 

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. 

488 

489 posa = self.positions[a] 

490 posb = self.positions[b] 

491 cola = self.colors[a] 

492 colb = self.colors[b] 

493 

494 if bond_order == 1: 

495 draw_tuples = ( 

496 (posa, mida, cola, transa, texa), 

497 (posb, midb, colb, transb, texb)) 

498 

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)) 

506 

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)) 

516 

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) 

523 

524 bondatoms = bondatoms.strip('\n') 

525 

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') 

538 

539 pov = f"""#version 3.6; 

540#include "colors.inc" 

541#include "finish.inc" 

542 

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}; 

555 

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 

564 

565{cell_vertices if cell_vertices != '' else '// no cell vertices'} 

566{atoms} 

567{bondatoms} 

568{constraints if constraints != '' else '// no constraints'} 

569""" # noqa: E501 

570 

571 with open(path, 'w') as fd: 

572 fd.write(pov) 

573 

574 return path 

575 

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) 

586 

587 

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 

593 

594 

595class POVRAYInputs: 

596 def __init__(self, path): 

597 self.path = path 

598 

599 def render(self, povray_executable='povray', stderr=DEVNULL, 

600 clean_up=False): 

601 cmd = [povray_executable, str(self.path)] 

602 

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}"') 

607 

608 if clean_up: 

609 unlink(self.path) 

610 unlink(self.path.with_suffix('.pov')) 

611 

612 return png_path 

613 

614 

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 """ 

647 

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 

653 

654 if self.gradient_direction == 'ascent': 

655 cv = 2 * cut_off 

656 else: 

657 cv = 0 

658 

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) 

669 

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)) 

674 

675 scaled_verts, faces, _normals, _values = self.compute_mesh( 

676 self.density_grid, 

677 self.cut_off, 

678 self.spacing, 

679 self.gradient_direction) 

680 

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 

687 

688 @property 

689 def cut_off(self): 

690 return self._cut_off 

691 

692 @cut_off.setter 

693 def cut_off(self, value): 

694 raise Exception("Use the set_cut_off method") 

695 

696 def set_cut_off(self, value): 

697 self._cut_off = value 

698 

699 if self.gradient_direction == 'ascent': 

700 cv = 2 * self.cut_off 

701 else: 

702 cv = 0 

703 

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) 

714 

715 self.spacing = tuple(1.0 / np.array(self.density_grid.shape)) 

716 

717 scaled_verts, faces, _, _ = self.compute_mesh( 

718 self.density_grid, 

719 self.cut_off, 

720 self.spacing, 

721 self.gradient_direction) 

722 

723 self.verts = scaled_verts 

724 self.faces = faces 

725 

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) 

732 

733 @staticmethod 

734 def wrapped_triples_section(triple_list, 

735 triple_format="<{:f}, {:f}, {:f}>".format, 

736 triples_per_line=4): 

737 

738 triples = [triple_format(*x) for x in triple_list] 

739 n = len(triples) 

740 s = '' 

741 tpl = triples_per_line 

742 c = 0 

743 

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 

751 

752 @staticmethod 

753 def compute_mesh(density_grid, cut_off, spacing, gradient_direction): 

754 """ 

755 

756 Import statement is in this method and not file header 

757 since few users will use isosurface rendering. 

758 

759 Returns scaled_verts, faces, normals, values. See skimage docs. 

760 

761 """ 

762 

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 

770 

771 return marching_cubes( 

772 density_grid, 

773 level=cut_off, 

774 spacing=spacing, 

775 gradient_direction=gradient_direction, 

776 allow_degenerate=False) 

777 

778 def format_mesh(self): 

779 """Returns a formatted data output for POVRAY files 

780 

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 

800 

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 

810 

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) 

816 

817 face_indices = self.wrapped_triples_section( 

818 triple_list=self.faces, 

819 triple_format="<{:n}, {:n}, {:n}>".format, 

820 triples_per_line=5) 

821 

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 

839 

840 

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) 

848 

849 

850def write_pov(filename, atoms, *, 

851 povray_settings=None, isosurface_data=None, 

852 **generic_projection_settings): 

853 

854 for name in ['run_povray', 'povray_path', 'stderr', 'extras']: 

855 pop_deprecated(generic_projection_settings, name) 

856 

857 if povray_settings is None: 

858 povray_settings = {} 

859 

860 pvars = PlottingVariables(atoms, scale=1.0, **generic_projection_settings) 

861 pov_obj = POVRAY.from_PlottingVariables(pvars, **povray_settings) 

862 

863 if isosurface_data is None: 

864 isosurface_data = [] 

865 elif not isinstance(isosurface_data, Sequence): 

866 isosurface_data = [isosurface_data] 

867 

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 

876 

877 return pov_obj.write(filename)