r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1"""

2Module for povray file format support.

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

5"""

6from collections.abc import Sequence

7from subprocess import check_call, DEVNULL

9from pathlib import Path

11import numpy as np

13from ase.io.utils import PlottingVariables

14from ase.constraints import FixAtoms

15from ase import Atoms

18def pa(array):

19 """Povray array syntax"""

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

23def pc(array):

24 """Povray color syntax"""

25 if isinstance(array, str):

26 return 'color ' + array

27 if isinstance(array, float):

28 return f'rgb <{array:.2f}>*3'.format(array)

29 L = len(array)

30 if L > 2 and L < 6:

31 return f"rgb{'' if L == 3 else 't' if L == 4 else 'ft'} <" +\

32 ', '.join(f"{x:.2f}" for x in tuple(array)) + '>'

36 """Get all pairs of bonding atoms

38 Return all pairs of atoms which are closer than radius times the

39 sum of their respective covalent radii. The pairs are returned as

40 tuples::

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

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

46 _ _ _

47 i c + i c + i c ,

48 1 1 2 2 3 3

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

51 integers."""

54 from ase.neighborlist import NeighborList

56 nl = NeighborList(cutoffs=cutoffs, self_interaction=False)

57 nl.update(atoms)

58 bondpairs = []

59 for a in range(len(atoms)):

60 indices, offsets = nl.get_neighbors(a)

61 bondpairs.extend([(a, a2, offset)

62 for a2, offset in zip(indices, offsets)])

63 return bondpairs

66def set_high_bondorder_pairs(bondpairs, high_bondorder_pairs=None):

67 """Set high bondorder pairs

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

70 bondorder pairs.

72 Parameters:

73 -----------

74 bondpairs: List of pairs, generated from get_bondpairs(atoms)

75 high_bondorder_pairs: Dictionary of pairs with high bond orders

76 using the following format:

77 { ( a1, b1 ): ( offset1, bond_order1, bond_offset1),

78 ( a2, b2 ): ( offset2, bond_order2, bond_offset2),

79 ...

80 }

81 offset, bond_order, bond_offset are optional.

82 However, if they are provided, the 1st value is

83 offset, 2nd value is bond_order,

84 3rd value is bond_offset """

86 if high_bondorder_pairs is None:

87 high_bondorder_pairs = 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_

102class POVRAY:

103 # These new styles were an attempt to port the old styles o the correct

104 # gamma, many had or still have unphysical light properties inorder to

105 # acheive a certain look.

106 material_styles_dict = dict(

107 simple='finish {phong 0.7 ambient 0.4 diffuse 0.55}',

108 # In general, 'pale' doesn't conserve energy and can look

109 # strange in many cases.

110 pale=('finish {ambient 0.9 diffuse 0.30 roughness 0.001 '

111 'specular 0.2 }'),

112 intermediate=('finish {ambient 0.4 diffuse 0.6 specular 0.1 '

113 'roughness 0.04}'),

114 vmd=(

115 'finish {ambient 0.2 diffuse 0.80 phong 0.25 phong_size 10.0 '

116 'specular 0.2 roughness 0.1}'),

117 jmol=('finish {ambient 0.4 diffuse 0.6 specular 1 roughness 0.001 '

118 'metallic}'),

119 ase2=('finish {ambient 0.2 brilliance 3 diffuse 0.6 metallic '

120 'specular 0.7 roughness 0.04 reflection 0.15}'),

121 ase3=('finish {ambient 0.4 brilliance 2 diffuse 0.6 metallic '

122 'specular 1.0 roughness 0.001 reflection 0.0}'),

123 glass=('finish {ambient 0.4 diffuse 0.35 specular 1.0 '

124 'roughness 0.001}'),

125 glass2=('finish {ambient 0.3 diffuse 0.3 specular 1.0 '

126 'reflection 0.25 roughness 0.001}'),

127 )

129 # These styles were made when assumed_gamma was 1.0 which gives poor color

130 # reproduction, the correct gamma is 2.2 for the sRGB standard.

131 material_styles_dict_old = dict(

132 simple='finish {phong 0.7}',

133 pale=('finish {ambient 0.5 diffuse 0.85 roughness 0.001 '

134 'specular 0.200 }'),

135 intermediate=('finish {ambient 0.3 diffuse 0.6 specular 0.1 '

136 'roughness 0.04}'),

137 vmd=('finish {ambient 0.0 diffuse 0.65 phong 0.1 phong_size 40.0 '

138 'specular 0.5 }'),

139 jmol=('finish {ambient 0.2 diffuse 0.6 specular 1 roughness 0.001 '

140 'metallic}'),

141 ase2=('finish {ambient 0.05 brilliance 3 diffuse 0.6 metallic '

142 'specular 0.7 roughness 0.04 reflection 0.15}'),

143 ase3=('finish {ambient 0.15 brilliance 2 diffuse 0.6 metallic '

144 'specular 1.0 roughness 0.001 reflection 0.0}'),

145 glass=('finish {ambient 0.05 diffuse 0.3 specular 1.0 '

146 'roughness 0.001}'),

147 glass2=('finish {ambient 0.01 diffuse 0.3 specular 1.0 '

148 'reflection 0.25 roughness 0.001}'),

149 )

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

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

231 # attributes from initialization

232 self.area_light = area_light

233 self.background = background

234 self.bondatoms = bondatoms

235 self.bondlinewidth = bondlinewidth

236 self.camera_dist = camera_dist

237 self.camera_type = camera_type

238 self.celllinewidth = celllinewidth

239 self.cue_density = cue_density

240 self.depth_cueing = depth_cueing

241 self.display = display

242 self.exportconstraints = exportconstraints

243 self.isosurfaces = isosurfaces

244 self.pause = pause

245 self.point_lights = point_lights

246 self.textures = textures

247 self.transmittances = transmittances

248 self.transparent = transparent

250 self.image_width = image_width

251 self.image_height = image_height

252 self.colors = colors

253 self.cell = cell

254 self.diameters = diameters

256 # calculations based on passed inputs

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

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

260 self.positions = positions - self.offset

262 if cell_vertices is not None:

263 self.cell_vertices = cell_vertices - self.offset

264 self.cell_vertices.shape = (2, 2, 2, 3)

265 else:

266 self.cell_vertices = None

268 ratio = float(self.image_width) / self.image_height

269 if canvas_width is None:

270 if canvas_height is None:

271 self.canvas_width = min(self.image_width * 15, 640)

272 self.canvas_height = min(self.image_height * 15, 640)

273 else:

274 self.canvas_width = canvas_height * ratio

275 self.canvas_height = canvas_height

276 elif canvas_height is None:

277 self.canvas_width = canvas_width

278 self.canvas_height = self.canvas_width / ratio

279 else:

280 raise RuntimeError("Can't set *both* width and height!")

282 # Distance to image plane from camera

283 if image_plane is None:

284 if self.camera_type == 'orthographic':

285 self.image_plane = 1 - self.camera_dist

286 else:

287 self.image_plane = 0

288 self.image_plane += self.camera_dist

290 self.constrainatoms = []

291 for c in constraints:

292 if isinstance(c, FixAtoms):

293 # self.constrainatoms.extend(c.index) # is this list-like?

294 for n, i in enumerate(c.index):

295 self.constrainatoms += [i]

297 @classmethod

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

299 cell = pvars.cell

300 cell_vertices = pvars.cell_vertices

301 if 'colors' in kwargs.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)

315 @classmethod

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

317 return cls.from_plotting_variables(

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

320 def write_ini(self, path):

321 """Write ini file."""

323 ini_str = f"""\

324Input_File_Name={path.with_suffix('.pov').name}

325Output_to_File=True

326Output_File_Type=N

327Output_Alpha={'on' if self.transparent else 'off'}

328; if you adjust Height, and width, you must preserve the ratio

329; Width / Height = {self.canvas_width/self.canvas_height:f}

330Width={self.canvas_width}

331Height={self.canvas_height}

332Antialias=True

333Antialias_Threshold=0.1

334Display={self.display}

335Display_Gamma=2.2

336Pause_When_Done={self.pause}

337Verbose=False

338"""

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

340 fd.write(ini_str)

341 return path

343 def write_pov(self, path):

344 """Write pov file."""

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

347 for loc, rgb in self.point_lights)

349 area_light = ''

350 if self.area_light is not None:

351 loc, color, width, height, nx, ny = self.area_light

352 area_light += f"""\nlight_source {{{pa(loc)} {pc(color)}

353 area_light <{width:.2f}, 0, 0>, <0, {height:.2f}, 0>, {nx:n}, {ny:n}

356 fog = ''

357 if self.depth_cueing and (self.cue_density >= 1e-4):

358 # same way vmd does it

359 if self.cue_density > 1e4:

360 # larger does not make any sense

361 dist = 1e-4

362 else:

363 dist = 1. / self.cue_density

364 fog += f'fog {{fog_type 1 distance {dist:.4f} '\

365 f'color {pc(self.background)}}}'

367 mat_style_keys = (f'#declare {k} = {v}'

368 for k, v in self.material_styles_dict.items())

369 mat_style_keys = '\n'.join(mat_style_keys)

371 # Draw unit cell

372 cell_vertices = ''

373 if self.cell_vertices is not None:

374 for c in range(3):

375 for j in ([0, 0], [1, 0], [1, 1], [0, 1]):

376 p1 = self.cell_vertices[tuple(j[:c]) + (0,) + tuple(j[c:])]

377 p2 = self.cell_vertices[tuple(j[:c]) + (1,) + tuple(j[c:])]

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

380 if distance < 1e-12:

381 continue

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

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

385 # all strings are f-strings for consistency

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

388 # Draw atoms

389 a = 0

390 atoms = ''

391 for loc, dia, col in zip(self.positions, self.diameters, self.colors):

392 tex = 'ase3'

393 trans = 0.

394 if self.textures is not None:

395 tex = self.textures[a]

396 if self.transmittances is not None:

397 trans = self.transmittances[a]

398 atoms += f'atom({pa(loc)}, {dia/2.:.2f}, {pc(col)}, '\

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

400 a += 1

401 atoms = atoms.strip('\n')

403 # Draw atom bonds

404 bondatoms = ''

405 for pair in self.bondatoms:

406 # Make sure that each pair has 4 componets: a, b, offset,

407 # bond_order, bond_offset

408 # a, b: atom index to draw bond

409 # offset: original meaning to make offset for mid-point.

410 # bond_oder: if not supplied, set it to 1 (single bond).

411 # It can be 1, 2, 3, corresponding to single,

412 # double, triple bond

413 # bond_offset: displacement from original bond position.

414 # Default is (bondlinewidth, bondlinewidth, 0)

415 # for bond_order > 1.

416 if len(pair) == 2:

417 a, b = pair

418 offset = (0, 0, 0)

419 bond_order = 1

420 bond_offset = (0, 0, 0)

421 elif len(pair) == 3:

422 a, b, offset = pair

423 bond_order = 1

424 bond_offset = (0, 0, 0)

425 elif len(pair) == 4:

426 a, b, offset, bond_order = pair

427 bond_offset = (self.bondlinewidth, self.bondlinewidth, 0)

428 elif len(pair) > 4:

429 a, b, offset, bond_order, bond_offset = pair

430 else:

431 raise RuntimeError('Each list in bondatom must have at least '

432 '2 entries. Error at %s' % pair)

434 if len(offset) != 3:

435 raise ValueError('offset must have 3 elements. '

436 'Error at %s' % pair)

437 if len(bond_offset) != 3:

438 raise ValueError('bond_offset must have 3 elements. '

439 'Error at %s' % pair)

440 if bond_order not in [0, 1, 2, 3]:

441 raise ValueError('bond_order must be either 0, 1, 2, or 3. '

442 'Error at %s' % pair)

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

445 # bond_offset for all bonds.

447 # Rotate bond_offset so that its direction is 90 deg. off the bond

448 # Utilize Atoms object to rotate

449 if bond_order > 1 and np.linalg.norm(bond_offset) > 1.e-9:

450 tmp_atoms = Atoms('H3')

451 tmp_atoms.set_cell(self.cell)

452 tmp_atoms.set_positions([

453 self.positions[a],

454 self.positions[b],

455 self.positions[b] + np.array(bond_offset),

456 ])

457 tmp_atoms.center()

458 tmp_atoms.set_angle(0, 1, 2, 90)

459 bond_offset = tmp_atoms[2].position - tmp_atoms[1].position

461 R = np.dot(offset, self.cell)

462 mida = 0.5 * (self.positions[a] + self.positions[b] + R)

463 midb = 0.5 * (self.positions[a] + self.positions[b] - R)

464 if self.textures is not None:

465 texa = self.textures[a]

466 texb = self.textures[b]

467 else:

468 texa = texb = 'ase3'

470 if self.transmittances is not None:

471 transa = self.transmittances[a]

472 transb = self.transmittances[b]

473 else:

474 transa = transb = 0.

476 # draw bond, according to its bond_order.

477 # bond_order == 0: No bond is plotted

478 # bond_order == 1: use original code

479 # bond_order == 2: draw two bonds, one is shifted by bond_offset/2,

480 # and another is shifted by -bond_offset/2.

481 # bond_order == 3: draw two bonds, one is shifted by bond_offset,

482 # and one is shifted by -bond_offset, and the

483 # other has no shift.

484 # To shift the bond, add the shift to the first two coordinate in

485 # write statement.

487 posa = self.positions[a]

488 posb = self.positions[b]

489 cola = self.colors[a]

490 colb = self.colors[b]

492 if bond_order == 1:

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

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

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)

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)

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)

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

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

534 pov = f"""#version 3.6;

535#include "colors.inc"

536#include "finish.inc"

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

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

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

561{atoms}

562{bondatoms}

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

564""" # noqa: E501

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

567 fd.write(pov)

569 return path

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)

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

590class POVRAYInputs:

591 def __init__(self, path):

592 self.path = path

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

595 clean_up=False):

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

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

603 if clean_up:

607 return png_path

610class POVRAYIsosurface:

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

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.

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

644 self.color = color

645 self.material = material

646 self.closed_edges = closed_edges

647 self._cut_off = cut_off

650 cv = 2 * cut_off

651 else:

652 cv = 0

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

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)

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

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

671 self.density_grid,

672 self.cut_off,

673 self.spacing,

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

683 @property

684 def cut_off(self):

685 return self._cut_off

687 @cut_off.setter

688 def cut_off(self, value):

689 raise Exception("Use the set_cut_off method")

691 def set_cut_off(self, value):

692 self._cut_off = value

695 cv = 2 * self.cut_off

696 else:

697 cv = 0

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

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)

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

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

713 self.density_grid,

714 self.cut_off,

715 self.spacing,

718 self.verts = scaled_verts

719 self.faces = faces

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)

728 @staticmethod

729 def wrapped_triples_section(triple_list,

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

731 triples_per_line=4):

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

734 n = len(triples)

735 s = ''

736 tpl = triples_per_line

737 c = 0

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

747 @staticmethod

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

749 """

751 Import statement is in this method and not file header

752 since few users will use isosurface rendering.

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

756 """

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)

767 return marching_cubes(

768 density_grid,

769 level=cut_off,

770 spacing=spacing,

772 allow_degenerate=False)

774 def format_mesh(self):

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

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

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

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)

813 face_indices = self.wrapped_triples_section(

814 triple_list=self.faces,

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

816 triples_per_line=5)

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

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)

846def write_pov(filename, atoms, *,

847 povray_settings=None, isosurface_data=None,

848 **generic_projection_settings):

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

851 pop_deprecated(generic_projection_settings, name)

853 if povray_settings is None:

854 povray_settings = {}

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

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

859 if isosurface_data is None:

860 isosurface_data = []

861 elif not isinstance(isosurface_data, Sequence):

862 isosurface_data = [isosurface_data]

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

873 return pov_obj.write(filename)