1""" 

2Module for povray file format support. 

3 

4See http://www.povray.org/ for details on the format. 

5""" 

6from collections.abc import Sequence 

7from subprocess import check_call, DEVNULL 

8from os import unlink 

9from pathlib import Path 

10 

11import numpy as np 

12 

13from ase.io.utils import PlottingVariables 

14from ase.constraints import FixAtoms 

15from ase import Atoms 

16 

17 

18def pa(array): 

19 """Povray array syntax""" 

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

21 

22 

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

33 

34 

35def get_bondpairs(atoms, radius=1.1): 

36 """Get all pairs of bonding atoms 

37 

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

41 

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

43 

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

45 

46 _ _ _ 

47 i c + i c + i c , 

48 1 1 2 2 3 3 

49 

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

51 integers.""" 

52 

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 

64 

65 

66def set_high_bondorder_pairs(bondpairs, high_bondorder_pairs=None): 

67 """Set high bondorder pairs 

68 

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

70 bondorder pairs. 

71 

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

85 

86 if high_bondorder_pairs is None: 

87 high_bondorder_pairs = dict() 

88 bondpairs_ = [] 

89 for pair in bondpairs: 

90 (a, b) = (pair[0], pair[1]) 

91 if (a, b) in high_bondorder_pairs.keys(): 

92 bondpair = [a, b] + [item for item in high_bondorder_pairs[(a, b)]] 

93 bondpairs_.append(bondpair) 

94 elif (b, a) in high_bondorder_pairs.keys(): 

95 bondpair = [a, b] + [item for item in high_bondorder_pairs[(b, a)]] 

96 bondpairs_.append(bondpair) 

97 else: 

98 bondpairs_.append(pair) 

99 return bondpairs_ 

100 

101 

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 ) 

128 

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 ) 

150 

151 def __init__(self, cell, cell_vertices, positions, diameters, colors, 

152 image_width, image_height, constraints=tuple(), 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?""" 

230 

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 

249 

250 self.image_width = image_width 

251 self.image_height = image_height 

252 self.colors = colors 

253 self.cell = cell 

254 self.diameters = diameters 

255 

256 # calculations based on passed inputs 

257 

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

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

260 self.positions = positions - self.offset 

261 

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 

267 

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

281 

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 

289 

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] 

296 

297 @classmethod 

298 def from_PlottingVariables(cls, pvars, **kwargs): 

299 cell = pvars.cell 

300 cell_vertices = pvars.cell_vertices 

301 if 'colors' in kwargs.keys(): 

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) 

314 

315 @classmethod 

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

317 return cls.from_plotting_variables( 

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

319 

320 def write_ini(self, path): 

321 """Write ini file.""" 

322 

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 

342 

343 def write_pov(self, path): 

344 """Write pov file.""" 

345 

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

347 for loc, rgb in self.point_lights) 

348 

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

355 

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

366 

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) 

370 

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

378 

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

380 if distance < 1e-12: 

381 continue 

382 

383 cell_vertices += f'cylinder {{{pa(p1)}, {pa(p2)}, '\ 

384 f'Rcell pigment {{Black}}}}\n' 

385 # all strings are f-strings for consistency 

386 cell_vertices = cell_vertices.strip('\n') 

387 

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

402 

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) 

433 

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) 

443 

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

445 # bond_offset for all bonds. 

446 

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 

460 

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' 

469 

470 if self.transmittances is not None: 

471 transa = self.transmittances[a] 

472 transb = self.transmittances[b] 

473 else: 

474 transa = transb = 0. 

475 

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. 

486 

487 posa = self.positions[a] 

488 posb = self.positions[b] 

489 cola = self.colors[a] 

490 colb = self.colors[b] 

491 

492 if bond_order == 1: 

493 draw_tuples = (posa, mida, cola, transa, texa),\ 

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

495 

496 elif bond_order == 2: 

497 bs = [x / 2 for x in bond_offset] 

498 draw_tuples = (posa - bs, mida - bs, cola, transa, texa),\ 

499 (posb - bs, midb - bs, colb, transb, texb),\ 

500 (posa + bs, mida + bs, cola, transa, texa),\ 

501 (posb + bs, midb + bs, colb, transb, texb) 

502 

503 elif bond_order == 3: 

504 bs = bond_offset 

505 draw_tuples = (posa, mida, cola, transa, texa),\ 

506 (posb, midb, colb, transb, texb),\ 

507 (posa + bs, mida + bs, cola, transa, texa),\ 

508 (posb + bs, midb + bs, colb, transb, texb),\ 

509 (posa - bs, mida - bs, cola, transa, texa),\ 

510 (posb - bs, midb - bs, colb, transb, texb) 

511 

512 bondatoms += ''.join(f'cylinder {{{pa(p)}, ' 

513 f'{pa(m)}, Rbond texture{{pigment ' 

514 f'{{color {pc(c)} ' 

515 f'transmit {tr}}} finish{{{tx}}}}}}}\n' 

516 for p, m, c, tr, tx in 

517 draw_tuples) 

518 

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

520 

521 # Draw constraints if requested 

522 constraints = '' 

523 if self.exportconstraints: 

524 for a in self.constrainatoms: 

525 dia = self.diameters[a] 

526 loc = self.positions[a] 

527 trans = 0.0 

528 if self.transmittances is not None: 

529 trans = self.transmittances[a] 

530 constraints += f'constrain({pa(loc)}, {dia/2.:.2f}, Black, '\ 

531 f'{trans}, {tex}) // #{a:n} \n' 

532 constraints = constraints.strip('\n') 

533 

534 pov = f"""#version 3.6; 

535#include "colors.inc" 

536#include "finish.inc" 

537 

538global_settings {{assumed_gamma 2.2 max_trace_level 6}} 

539background {{{pc(self.background)}{' transmit 1.0' if self.transparent else ''}}} 

540camera {{{self.camera_type} 

541 right -{self.image_width:.2f}*x up {self.image_height:.2f}*y 

542 direction {self.image_plane:.2f}*z 

543 location <0,0,{self.camera_dist:.2f}> look_at <0,0,0>}} 

544{point_lights} 

545{area_light if area_light != '' else '// no area light'} 

546{fog if fog != '' else '// no fog'} 

547{mat_style_keys} 

548#declare Rcell = {self.celllinewidth:.3f}; 

549#declare Rbond = {self.bondlinewidth:.3f}; 

550 

551#macro atom(LOC, R, COL, TRANS, FIN) 

552 sphere{{LOC, R texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}} 

553#end 

554#macro constrain(LOC, R, COL, TRANS FIN) 

555union{{torus{{R, Rcell rotate 45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}} 

556 torus{{R, Rcell rotate -45*z texture{{pigment{{color COL transmit TRANS}} finish{{FIN}}}}}} 

557 translate LOC}} 

558#end 

559 

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

561{atoms} 

562{bondatoms} 

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

564""" # noqa: E501 

565 

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

567 fd.write(pov) 

568 

569 return path 

570 

571 def write(self, pov_path): 

572 pov_path = require_pov(pov_path) 

573 ini_path = pov_path.with_suffix('.ini') 

574 self.write_ini(ini_path) 

575 self.write_pov(pov_path) 

576 if self.isosurfaces is not None: 

577 with open(pov_path, 'a') as fd: 

578 for iso in self.isosurfaces: 

579 fd.write(iso.format_mesh()) 

580 return POVRAYInputs(ini_path) 

581 

582 

583def require_pov(path): 

584 path = Path(path) 

585 if path.suffix != '.pov': 

586 raise ValueError(f'Expected .pov path, got {path}') 

587 return path 

588 

589 

590class POVRAYInputs: 

591 def __init__(self, path): 

592 self.path = path 

593 

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

595 clean_up=False): 

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

597 

598 check_call(cmd, stderr=stderr) 

599 png_path = self.path.with_suffix('.png').absolute() 

600 if not png_path.is_file(): 

601 raise RuntimeError(f'Povray left no output PNG file "{png_path}"') 

602 

603 if clean_up: 

604 unlink(self.path) 

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

606 

607 return png_path 

608 

609 

610class POVRAYIsosurface: 

611 def __init__(self, density_grid, cut_off, cell, cell_origin, 

612 closed_edges=False, gradient_ascending=False, 

613 color=(0.85, 0.80, 0.25, 0.2), material='ase3'): 

614 """ 

615 density_grid: 3D float ndarray 

616 A regular grid on that spans the cell. The first dimension 

617 corresponds to the first cell vector and so on. 

618 cut_off: float 

619 The density value of the isosurface. 

620 cell: 2D float ndarray or ASE cell object 

621 The 3 vectors which give the cell's repetition 

622 cell_origin: 4 float tuple 

623 The cell origin as used by POVRAY object 

624 closed_edges: bool 

625 Setting this will fill in isosurface edges at the cell boundaries. 

626 Filling in the edges can help with visualizing 

627 highly porous structures. 

628 gradient_ascending: bool 

629 Lets you pick the area you want to enclose, i.e., should the denser 

630 or less dense area be filled in. 

631 color: povray color string, float, or float tuple 

632 1 float is interpreted as grey scale, a 3 float tuple is rgb, 

633 4 float tuple is rgbt, and 5 float tuple is rgbft, where 

634 t is transmission fraction and f is filter fraction. 

635 Named Povray colors are set in colors.inc 

636 (http://wiki.povray.org/content/Reference:Colors.inc) 

637 material: string 

638 Can be a finish macro defined by POVRAY.material_styles 

639 or a full Povray material {...} specification. Using a 

640 full material specification willoverride the color parameter. 

641 """ 

642 

643 self.gradient_direction = 'ascent' if gradient_ascending else 'descent' 

644 self.color = color 

645 self.material = material 

646 self.closed_edges = closed_edges 

647 self._cut_off = cut_off 

648 

649 if self.gradient_direction == 'ascent': 

650 cv = 2 * cut_off 

651 else: 

652 cv = 0 

653 

654 if closed_edges: 

655 shape_old = density_grid.shape 

656 # since well be padding, we need to keep the data at origin 

657 cell_origin += -(1.0 / np.array(shape_old)) @ cell 

658 density_grid = np.pad( 

659 density_grid, pad_width=( 

660 1,), mode='constant', constant_values=cv) 

661 shape_new = density_grid.shape 

662 s = np.array(shape_new) / np.array(shape_old) 

663 cell = cell @ np.diag(s) 

664 

665 self.cell = cell 

666 self.cell_origin = cell_origin 

667 self.density_grid = density_grid 

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

669 

670 scaled_verts, faces, normals, values = self.compute_mesh( 

671 self.density_grid, 

672 self.cut_off, 

673 self.spacing, 

674 self.gradient_direction) 

675 

676 # The verts are scaled by default, this is the super easy way of 

677 # distributing them in real space but it's easier to do affine 

678 # transformations/rotations on a unit cube so I leave it like that 

679 # verts = scaled_verts.dot(self.cell) 

680 self.verts = scaled_verts 

681 self.faces = faces 

682 

683 @property 

684 def cut_off(self): 

685 return self._cut_off 

686 

687 @cut_off.setter 

688 def cut_off(self, value): 

689 raise Exception("Use the set_cut_off method") 

690 

691 def set_cut_off(self, value): 

692 self._cut_off = value 

693 

694 if self.gradient_direction == 'ascent': 

695 cv = 2 * self.cut_off 

696 else: 

697 cv = 0 

698 

699 if self.closed_edges: 

700 shape_old = self.density_grid.shape 

701 # since well be padding, we need to keep the data at origin 

702 self.cell_origin += -(1.0 / np.array(shape_old)) @ self.cell 

703 self.density_grid = np.pad( 

704 self.density_grid, pad_width=( 

705 1,), mode='constant', constant_values=cv) 

706 shape_new = self.density_grid.shape 

707 s = np.array(shape_new) / np.array(shape_old) 

708 self.cell = self.cell @ np.diag(s) 

709 

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

711 

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

713 self.density_grid, 

714 self.cut_off, 

715 self.spacing, 

716 self.gradient_direction) 

717 

718 self.verts = scaled_verts 

719 self.faces = faces 

720 

721 @classmethod 

722 def from_POVRAY(cls, povray, density_grid, cut_off, **kwargs): 

723 return cls(cell=povray.cell, 

724 cell_origin=povray.cell_vertices[0, 0, 0], 

725 density_grid=density_grid, 

726 cut_off=cut_off, **kwargs) 

727 

728 @staticmethod 

729 def wrapped_triples_section(triple_list, 

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

731 triples_per_line=4): 

732 

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

734 n = len(triples) 

735 s = '' 

736 tpl = triples_per_line 

737 c = 0 

738 

739 while c < n - tpl: 

740 c += tpl 

741 s += '\n ' 

742 s += ', '.join(triples[c - tpl:c]) 

743 s += '\n ' 

744 s += ', '.join(triples[c:]) 

745 return s 

746 

747 @staticmethod 

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

749 """ 

750 

751 Import statement is in this method and not file header 

752 since few users will use isosurface rendering. 

753 

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

755 

756 """ 

757 

758 # marching_cubes name was changed in skimage v0.19 

759 try: 

760 # New skimage 

761 from skimage.measure import marching_cubes 

762 except ImportError: 

763 # Old skimage (remove at some point) 

764 from skimage.measure import ( 

765 marching_cubes_lewiner as marching_cubes) 

766 

767 return marching_cubes( 

768 density_grid, 

769 level=cut_off, 

770 spacing=spacing, 

771 gradient_direction=gradient_direction, 

772 allow_degenerate=False) 

773 

774 def format_mesh(self): 

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

776 

777 Example: 

778 material = ''' 

779 material { // This material looks like pink jelly 

780 texture { 

781 pigment { rgbt <0.8, 0.25, 0.25, 0.5> } 

782 finish{ diffuse 0.85 ambient 0.99 brilliance 3 specular 0.5 roughness 0.001 

783 reflection { 0.05, 0.98 fresnel on exponent 1.5 } 

784 conserve_energy 

785 } 

786 } 

787 interior { ior 1.3 } 

788 } 

789 photons { 

790 target 

791 refraction on 

792 reflection on 

793 collect on 

794 }''' 

795 """ # noqa: E501 

796 

797 if self.material in POVRAY.material_styles_dict: 

798 material = f"""material {{ 

799 texture {{ 

800 pigment {{ {pc(self.color)} }} 

801 finish {{ {self.material} }} 

802 }} 

803 }}""" 

804 else: 

805 material = self.material 

806 

807 # Start writing the mesh2 

808 vertex_vectors = self.wrapped_triples_section( 

809 triple_list=self.verts, 

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

811 triples_per_line=4) 

812 

813 face_indices = self.wrapped_triples_section( 

814 triple_list=self.faces, 

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

816 triples_per_line=5) 

817 

818 cell = self.cell 

819 cell_or = self.cell_origin 

820 mesh2 = f"""\n\nmesh2 {{ 

821 vertex_vectors {{ {len(self.verts):n}, 

822 {vertex_vectors} 

823 }} 

824 face_indices {{ {len(self.faces):n}, 

825 {face_indices} 

826 }} 

827{material if material != '' else '// no material'} 

828 matrix < {cell[0][0]:f}, {cell[0][1]:f}, {cell[0][2]:f}, 

829 {cell[1][0]:f}, {cell[1][1]:f}, {cell[1][2]:f}, 

830 {cell[2][0]:f}, {cell[2][1]:f}, {cell[2][2]:f}, 

831 {cell_or[0]:f}, {cell_or[1]:f}, {cell_or[2]:f}> 

832 }} 

833 """ 

834 return mesh2 

835 

836 

837def pop_deprecated(dct, name): 

838 import warnings 

839 if name in dct: 

840 del dct[name] 

841 warnings.warn(f'The "{name}" keyword of write_pov() is deprecated ' 

842 'and has no effect; this will raise an error in the ' 

843 'future.', FutureWarning) 

844 

845 

846def write_pov(filename, atoms, *, 

847 povray_settings=None, isosurface_data=None, 

848 **generic_projection_settings): 

849 

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

851 pop_deprecated(generic_projection_settings, name) 

852 

853 if povray_settings is None: 

854 povray_settings = {} 

855 

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

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

858 

859 if isosurface_data is None: 

860 isosurface_data = [] 

861 elif not isinstance(isosurface_data, Sequence): 

862 isosurface_data = [isosurface_data] 

863 

864 isosurfaces = [] 

865 for isodata in isosurface_data: 

866 if isinstance(isodata, POVRAYIsosurface): 

867 iso = isodata 

868 else: 

869 iso = POVRAYIsosurface.from_POVRAY(pov_obj, **isodata) 

870 isosurfaces.append(iso) 

871 pov_obj.isosurfaces = isosurfaces 

872 

873 return pov_obj.write(filename)