small changes

scripts/ase_spacegroup.py

@@ -0,0 +1,39 @@
+import ase.io as io
+from ase.build import cut
+from ase.spacegroup import crystal
+
+
+
+a = 9.04
+skutterudite = crystal(('Co', 'Sb'),
+                       basis=[(0.25, 0.25, 0.25), (0.0, 0.335, 0.158)],
+                       spacegroup=204,
+                       cellpar=[a, a, a, 90, 90, 90])
+
+# Create a new atoms instance with Co at origo including all atoms on the
+# surface of the unit cell
+cosb3 = cut(skutterudite, origo=(0.25, 0.25, 0.25), extend=1.01)
+
+# Define the atomic bonds to show
+bondatoms = []
+symbols = cosb3.get_chemical_symbols()
+for i in range(len(cosb3)):
+    for j in range(i):
+        if (symbols[i] == symbols[j] == 'Co' and
+                cosb3.get_distance(i, j) < 4.53):
+            bondatoms.append((i, j))
+        elif (symbols[i] == symbols[j] == 'Sb' and
+              cosb3.get_distance(i, j) < 2.99):
+            bondatoms.append((i, j))
+
+# Create nice-looking image using povray
+renderer = io.write('spacegroup-cosb3.pov', cosb3,
+                    rotation='90y',
+                    radii=0.4,
+                    povray_settings=dict(transparent=False,
+                                         camera_type='perspective',
+                                         canvas_width=320,
+                                         bondlinewidth=0.07,
+                                         bondatoms=bondatoms))
+
+renderer.render()

scripts/bz.py

@@ -0,0 +1,50 @@
+# creates: bztable.rst
+# creates: 00.CUB.svg 01.FCC.svg 02.BCC.svg 03.TET.svg 04.BCT1.svg
+# creates: 05.BCT2.svg 06.ORC.svg 07.ORCF1.svg 08.ORCF2.svg 09.ORCF3.svg
+# creates: 10.ORCI.svg 11.ORCC.svg 12.HEX.svg 13.RHL1.svg 14.RHL2.svg
+# creates: 15.MCL.svg 16.MCLC1.svg 17.MCLC3.svg 18.MCLC5.svg 19.TRI1a.svg
+# creates: 20.TRI1b.svg 21.TRI2a.svg 22.TRI2b.svg
+# creates: 23.OBL.svg 24.RECT.svg 25.CRECT.svg 26.HEX2D.svg 27.SQR.svg
+
+# taken from https://wiki.fysik.dtu.dk/ase/gallery/gallery.html
+
+
+from formulary import *
+
+from ase.lattice import all_variants
+
+from ase.data import colors
+
+
+header = """\
+
+Brillouin zone data
+-------------------
+
+.. list-table::
+    :widths: 10 15 45
+"""
+
+
+entry = """\
+    * - {name} ({longname})
+      - {bandpath}
+      - .. image:: {fname}
+            :width: 40 %
+"""
+
+with open('bztable.rst', 'w') as fd:
+    print(header, file=fd)
+
+    for i, lat in enumerate(all_variants()):
+        id = f'{i:02d}.{lat.variant}'
+        imagefname = f'out/{id}.svg'
+        txt = entry.format(name=lat.variant,
+                           longname=lat.longname,
+                           bandpath=lat.bandpath().path,
+                           fname=imagefname)
+        print(txt, file=fd)
+        ax = lat.plot_bz()
+        fig = ax.get_figure()
+        fig.savefig(imagefname, bbox_inches='tight')
+        fig.clear()

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/cm_phonons.py

@@ -1,5 +1,9 @@
 #!/usr/bin env python3
 from formulary import *
+from scipy.constants import Boltzmann as kB, hbar
+
+hbar = 1
+kB = 1
 
 def fone_atom_basis(q, a, M, C1, C2):
     return np.sqrt(4*C1/M * (np.sin(q*a/2)**2 + C2/C1 * np.sin(q*a)**2))
@@ -57,8 +61,56 @@ def two_atom_basis():
     ax.set_xticklabels([f"${i}\\pi/a$" if i != 0 else "0" for i in range(-2, 3)])
     ax.legend()
     ax.grid()
-
     return fig
+
+
+def fcv_einstein(T, N, omegaE):
+    ThetaT = hbar * omegaE / (kB * T)
+    return 3 * N * kB * ThetaT**2 * np.exp(ThetaT) / (np.exp(ThetaT) - 1)**2
+
+def fcv_debye_integral(x):
+    print(np.exp(x), (np.exp(x) - 1)**2)
+    return x**4 * np.exp(x) / ((np.exp(x) - 1)**2)
+
+def heat_capacity_einstein_debye():
+    Ts = np.linspace(0, 10, 500)
+    omegaD = 1e1
+    omegaE = 1
+    # N = 10**23
+    N = 1
+    cvs_einstein = fcv_einstein(Ts, N, omegaE)
+    cvs_debye = np.zeros(Ts.shape, dtype=float)
+    integral = np.zeros(Ts.shape, dtype=float)
+    # cvs_debye = [0.0 for _ in range(Ts.shape[0])] # np.zeros(Ts.shape, dtype=float)
+    # integral = [0.0 for _ in range(Ts.shape[0])] # np.zeros(Ts.shape, dtype=float)
+    dT = Ts[1] - Ts[0]
+    dThetaT = kB*dT/(hbar*omegaD)
+    for i, T in enumerate(Ts):
+        if i == 0: continue
+        ThetaT = kB*T/(hbar*omegaD)
+        dIntegral = fcv_debye_integral(ThetaT) * dThetaT
+        integral[i] = dIntegral
+        # print(integral)
+        integral[i] += integral[i-1]
+        C_debye = 9 * N * kB * ThetaT**3 * integral[i]
+        cvs_debye[i] = C_debye
+        print(i, T, ThetaT, dIntegral, C_debye, integral[i])
+    fig, ax = plt.subplots(1, 1, figsize=size_formula_normal_default)
+    ax.set_xlabel("$T$")
+    ax.set_ylabel("$c_V$")
+    ax.plot(Ts, cvs_einstein, label="Einstein")
+    ax.plot(Ts, cvs_debye, label="Debye")
+    ax.plot(Ts, integral, label="integral")
+    ax.hlines([3*N*kB], xmin=0, xmax=Ts[-1], colors=COLORSCHEME["fg1"], linestyles="dashed")
+    # print(cvs_debye)
+    ax.legend()
+    return fig
+
+
+
+
 if __name__ == '__main__':
-    export(one_atom_basis(), "cm_phonon_dispersion_one_atom_basis")
+    export(one_atom_basis(), "cm_vib_dispersion_one_atom_basis")
-    export(two_atom_basis(), "cm_phonon_dispersion_two_atom_basis")
+    export(two_atom_basis(), "cm_vib_dispersion_two_atom_basis")
+    export(heat_capacity_einstein_debye(), "cm_vib_heat_capacity_einstein_debye")
+    print(kB)

scripts/cm_superconductivity.py

@@ -45,6 +45,74 @@ def n_s_boundary():
     ax.legend(loc='center right', handles=lines)
     return fig
 
+from mpl_toolkits.mplot3d import Axes3D
+from scipy.interpolate import griddata
+
+def critical_type2():
+    Jc0 = 100
+    Bc2_0 = 30
+    Tc = 90
+
+    T = np.linspace(0, Tc, 100)
+    Jc_T = Jc0 * (1 - (T / Tc)**2)
+    Bc2_T = Bc2_0 * (1 - (T / Tc)**2)
+    B = np.linspace(0, Bc2_0, 100)
+    Jc_B = Jc0 * (1 - B / Bc2_0)
+
+    fig = plt.figure(figsize=size_formula_normal_default)
+    ax = fig.add_subplot(111, projection='3d')
+
+    ax.plot(T, np.zeros_like(Jc_T), Jc_T, label='$J_c(T)$', color='r')
+    ax.plot(T, Bc2_T, np.zeros_like(Bc2_T), label='$B_{c2}(T)$', color='g')
+    ax.plot(np.zeros_like(Jc_B), B, Jc_B, label='$J_c(B)$', color='b')
+
+    ax.set_xlim(0, Tc)
+    ax.set_ylim(0, Bc2_0)
+    ax.set_zlim(0, Jc0)
+
+    # surface
+    # T_grid, B_grid = np.meshgrid(T, B)
+    # Jc_grid = Jc0 * (1 - (T_grid / Tc)**2) * (1 - B_grid / Bc2_0)
+    # surf = ax.plot_surface(T_grid, B_grid, Jc_grid, color='cyan', alpha=0.5)
+    ax.set_xlabel('$T$')
+    ax.set_ylabel('$B_{c2}$')
+    ax.set_zlabel('$J_c$')
+    # ax.legend()
+    ax.grid(True)
+
+    ax.view_init(elev=30., azim=45)
+    ax.set_box_aspect(None, zoom=0.85)
+    return fig
+
+
+
+def heat_capacity():
+    fig, ax = plt.subplots(1, 1, figsize=size_formula_small_quadratic)
+
+    T_max  = 1.7
+    Cn_max = 3
+    f_Cn = lambda T: T * Cn_max/T_max
+    Delta_C = f_Cn(1.0) * 1.43  # BCS prediction
+    CsTc = f_Cn(1.0) * (1+1.43)  # BCS prediction
+    # exp decay from there
+    f_Cs = lambda T:  np.exp(-1 / T + 1) * CsTc
+
+    Tns = np.linspace(0.0, T_max, 100)
+    Tss = np.linspace(0.0, 1.0, 100)
+    Cns = f_Cn(Tns)
+    Css = f_Cs(Tss)
+    ax.plot(Tns, Cns, label=r"$c_\text{n}$")
+    ax.plot(Tss, Css, label=r"$c_\text{s}$")
+    ax.vlines([1.0], ymin=f_Cn(1.0), ymax=(CsTc), color=COLORSCHEME["fg1"], linestyles="dashed")
+    ax.text(1.05, CsTc - Delta_C/2, "$\\Delta c$",  color=COLORSCHEME["fg1"])
+    ax.set_xlabel(r"$T/T_\text{c}$")
+    ax.set_ylabel(r"$c$ [a.u.]")
+    ax.legend()
+    return fig
+
+
 
 if __name__ == "__main__":
-    export(n_s_boundary(), "cm_sc_n_s_boundary")
+    export(n_s_boundary(), "cm_super_n_s_boundary")
+    export(critical_type2(), "cm_super_critical_type2")
+    export(heat_capacity(), "cm_super_heat_capacity")

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/readme.md

@@ -3,11 +3,19 @@ Put all scripts that generate plots or tex files here.
 You can run all files at once using `make scripts`
 
 ## Plots
+### `matplotlib`
 For plots with `matplotlib`:
 1. import `formulary.py` 
 2. use one of the preset figsizes
 3. save the image using the `export` function in the `if __name__ == '__main__'` part
 
+### `ase` - Atomic Simulation Environment
+For plots with `ase`:
+1. import `formulary.py` and `util.aseutil`
+2. Use `util.aseutil.set_atom_color` to change the color of all used atoms to one in the colorscheme
+3. export the render using the `export_atoms` function in the `if __name__ == '__main__'` part.
+   Pass one of the preset figsizes as size.
+
 ## Colorscheme
 To ensure a uniform look of the tex source and the python plots, 
 the tex and matplotlib colorschemes are both handled in `formulary.py`.

scripts/requirements.txt

@@ -3,4 +3,4 @@ scipy
 matplotlib
 scqubits
 qutip
-
+ase

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("#")

src/img_static_svgs/convertToPdf.sh

@@ -1,7 +1,7 @@
 IFS=$'\n'
 for d in $(find . -type d); do
-    mkdir -p "../img/$d"
+    mkdir -p "../img_static/$d"
 done
 for file in $(find . -type f -name '*.svg'); do
-    inkscape -o "../img/${file%.*}.pdf" "$file"
+    inkscape -o "../img_static/${file%.*}.pdf" "$file"
 done