Add crystal structure renders using ASE and povray

scripts/cm_crystal_structures.py

@@ -0,0 +1,51 @@
+from formulary import *
+from util.aseutil import set_atom_color, get_pov_settings
+
+"""
+Create crystal structures using ase and render them with povray
+
+Rotation angle:
+To get the rotation angle, open the structure in the ase.visualize.view
+and use "View->Rotation" to get the desired angles
+"""
+
+set_atom_color("Na", COLORSCHEME["fg-red"])
+set_atom_color("Cl", COLORSCHEME["fg-blue"])
+
+set_atom_color("Zn", COLORSCHEME["fg-blue"])
+set_atom_color("S", COLORSCHEME["fg-yellow"])
+
+from ase.lattice import compounds
+from ase.build import cut, bulk
+from ase import Atom, Atoms
+
+
+def zincblende():
+    zns = compounds.Zincblende(("Zn", "S"), latticeconstant=5.0, size=(1,1,1))
+    zns_cell = cut(zns, b=(0,0,1), origo=(0,0,0), extend=1.1)
+    return zns_cell
+
+# NaCl cut
+def nacl():
+    nacl = compounds.NaCl(("Na", "Cl"), latticeconstant=5.0, size=(1,1,1))
+    nacl_cell = cut(nacl, b=(0,0,1), origo=(0,0,0), extend=1.1)
+    return nacl_cell
+
+def wurtzite():
+    compounds.L1_2
+    wurtzite = bulk('SZn', 'wurtzite', a=3.129, c=5.017)
+    wurtzite_cell = cut(wurtzite,
+        a=[1, 0, 0],
+        b=[-1, -1, 0],
+        c=[0, 0, 1], extend=1.1)
+    return wurtzite_cell
+
+
+if __name__ == "__main__":
+    export_atoms(nacl(), "cm_crystal_NaCl", size_formula_half_quadratic)
+    export_atoms(wurtzite(), "cm_crystal_wurtzite", size_formula_half_quadratic, rotation="70x,20y,174z")
+    export_atoms(zincblende(), "cm_crystal_zincblende", size_formula_half_quadratic, rotation="-155x,70y,24z")
+    w = wurtzite()
+    from ase.visualize import view
+    view(w)
+

scripts/formulary.py

@@ -31,9 +31,11 @@ COLORSCHEME = cs.gruvbox_light()
 # cs.p_tum["fg0"] = cs.p_tum["alt-blue"]
 # COLORSCHEME = cs.tum()
 # COLORSCHEME = cs.legacy()
+# COLORSCHEME = cs.stupid()
 
+tex_aux_path = "../.aux/"
 tex_src_path = "../src/"
-img_out_dir = os.path.join(tex_src_path, "img")
+img_out_dir = os.path.abspath(os.path.join(tex_src_path, "img"))
 filetype = ".pdf"
 skipasserts = False
 
@@ -74,6 +76,34 @@ def export(fig, name, tight_layout=True):
         fig.tight_layout()
     fig.savefig(filename, bbox_inches="tight", pad_inches=0.0)
 
+
+def export_atoms(atoms, name, size, rotation="-30y,20x", get_bonds=True):
+    """Export a render of ase atoms object"""
+    assert_directory()
+    wd = os.getcwd()
+    from util.aseutil import get_bondatoms, get_pov_settings
+    from ase import io
+
+    tmp_dir = os.path.join(os.path.abspath(tex_aux_path), "scripts_aux")
+    os.makedirs(tmp_dir, exist_ok=True)
+    os.chdir(tmp_dir)
+
+    out_filename = f"{name}.png"
+
+    bondatoms = None
+    if get_bonds:
+        bondatoms = get_bondatoms(atoms)
+    renderer = io.write(f'{name}.pov', atoms,
+            rotation=rotation,# text string with rotation (default='' )
+            radii=0.4,          # float, or a list with one float per atom
+            show_unit_cell=2,   # 0, 1, or 2 to not show, show, and show all of cell
+            colors=None,        # List: one (r, g, b, t) tuple per atom  
+            povray_settings=get_pov_settings(size, COLORSCHEME, bondatoms),
+    )
+    renderer.render()
+    os.chdir(wd)
+    os.rename(os.path.join(tmp_dir, out_filename), os.path.join(img_out_dir, out_filename))
+
 @np.vectorize
 def smooth_step(x: float, left_edge: float, right_edge: float):
     x = (x - left_edge) / (right_edge - left_edge)

scripts/util/aseutil.py

@@ -0,0 +1,47 @@
+from util.colorschemes import hex_to_rgb_float
+
+def set_atom_color(symbol, hexcolor):
+    from ase.data import atomic_numbers
+    from ase.data.colors import jmol_colors, cpk_colors
+    float_color = hex_to_rgb_float(hexcolor)
+    n = atomic_numbers[symbol]
+    jmol_colors[n] = float_color
+    cpk_colors[n] = float_color
+
+
+from scipy.spatial.distance import pdist, squareform
+import numpy as np
+def get_bondatoms(atoms):
+    site_positions = [site.position for site in atoms]
+    pair_distances = squareform(pdist(np.stack(site_positions)))
+    vs = pair_distances
+    bondatoms = []
+    for i in range(vs.shape[0]):
+        for j in range(i):
+            if vs[i, j] < 3: # up to 3 angstrom distance show a bond TODO
+                bondatoms.append((i, j))
+    return bondatoms
+# returns to many
+# from ase.io.pov import get_bondpairs
+# bondatoms=get_bondpairs(lat, 5)
+
+
+TARGET_DPI = 300
+# doc: https://github.com/WMD-group/ASE-Tutorials/blob/master/povray-tools/ase_povray.py
+def get_pov_settings(size, COLORSCHEME, bondatoms=None):
+    white = hex_to_rgb_float(COLORSCHEME["bg0"])
+    other = hex_to_rgb_float(COLORSCHEME["fg-yellow"])
+    pixels = TARGET_DPI * size[0]
+    pov_settings=dict(
+        transparent=True,
+        display=False,
+        # camera_type='orthographic',
+        camera_type='perspective',
+        canvas_width=pixels,
+        # point_lights : [], #[(18,20,40), 'White'],[(60,20,40),'White'],             # [[loc1, color1], [loc2, color2],...]
+        point_lights=[[(18,20,40), white],[(60,20,40),other]],             # [[loc1, color1], [loc2, color2],...]
+        background=(0, 0, 0, 1.,),
+        bondlinewidth=0.07,
+        bondatoms=bondatoms
+    )
+    return pov_settings

scripts/util/colorschemes.py

@@ -9,6 +9,26 @@ from math import floor
 
 colors = ["red", "orange", "yellow", "green", "aqua", "blue", "purple", "gray"]
 
+def duplicate_letters(color: str):
+    return ''.join([c+c for c in color])
+
+def hex_to_rgb_int(color: str) -> list[int]:
+    color = color.strip("#")
+    ctuple = []
+    # turn RGBA to RRGGBBAA
+    if len(color) == 3 or len(color) == 4:
+        color = duplicate_letters(color)
+    for i in range(len(color)//2):
+        ctuple.append(int(color[i*2:i*2+2], 16))
+    return ctuple
+
+def hex_to_rgb_float(color: str) -> list[float]:
+    clist = hex_to_rgb_int(color)
+    fclist = [float(c) / 255 for c in clist]
+    return fclist
+
+
+
 def brightness(color:str, percent:float):
     if color.startswith("#"):
         color = color.strip("#")