From b5450708c10b680a533e42148e27d3e7aa0c2c34 Mon Sep 17 00:00:00 2001 From: "matthias@quintern.xyz" Date: Sun, 9 Mar 2025 20:24:15 +0100 Subject: [PATCH] Add crystal structure renders using ASE and povray --- scripts/cm_crystal_structures.py | 51 ++++++++++++++++++++++++++++++++ scripts/formulary.py | 32 +++++++++++++++++++- scripts/util/aseutil.py | 47 +++++++++++++++++++++++++++++ scripts/util/colorschemes.py | 20 +++++++++++++ 4 files changed, 149 insertions(+), 1 deletion(-) create mode 100644 scripts/cm_crystal_structures.py create mode 100644 scripts/util/aseutil.py diff --git a/scripts/cm_crystal_structures.py b/scripts/cm_crystal_structures.py new file mode 100644 index 0000000..698bcea --- /dev/null +++ b/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) + diff --git a/scripts/formulary.py b/scripts/formulary.py index 8fda69f..14da93e 100644 --- a/scripts/formulary.py +++ b/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) diff --git a/scripts/util/aseutil.py b/scripts/util/aseutil.py new file mode 100644 index 0000000..6fa832a --- /dev/null +++ b/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 diff --git a/scripts/util/colorschemes.py b/scripts/util/colorschemes.py index 373b2e8..f9be34e 100644 --- a/scripts/util/colorschemes.py +++ b/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("#")