From ac4ff1d16b6fbc8a14bb0a4e428ad927998f36c2 Mon Sep 17 00:00:00 2001 From: "matthias@quintern.xyz" Date: Sun, 9 Mar 2025 20:25:09 +0100 Subject: [PATCH] small changes --- scripts/ase_spacegroup.py | 39 ++++++++++++++++++ scripts/bz.py | 50 +++++++++++++++++++++++ scripts/cm_phonons.py | 58 +++++++++++++++++++++++++-- scripts/cm_superconductivity.py | 70 ++++++++++++++++++++++++++++++++- scripts/readme.md | 8 ++++ scripts/requirements.txt | 2 +- 6 files changed, 222 insertions(+), 5 deletions(-) create mode 100644 scripts/ase_spacegroup.py create mode 100644 scripts/bz.py diff --git a/scripts/ase_spacegroup.py b/scripts/ase_spacegroup.py new file mode 100644 index 0000000..0cc0266 --- /dev/null +++ b/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() diff --git a/scripts/bz.py b/scripts/bz.py new file mode 100644 index 0000000..8ef45b2 --- /dev/null +++ b/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() diff --git a/scripts/cm_phonons.py b/scripts/cm_phonons.py index 10db99a..d200a85 100644 --- a/scripts/cm_phonons.py +++ b/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(two_atom_basis(), "cm_phonon_dispersion_two_atom_basis") + export(one_atom_basis(), "cm_vib_dispersion_one_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) diff --git a/scripts/cm_superconductivity.py b/scripts/cm_superconductivity.py index 586c546..318294a 100644 --- a/scripts/cm_superconductivity.py +++ b/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") diff --git a/scripts/readme.md b/scripts/readme.md index bfc5912..18fc6c6 100644 --- a/scripts/readme.md +++ b/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`. diff --git a/scripts/requirements.txt b/scripts/requirements.txt index 4761a58..7e35c47 100644 --- a/scripts/requirements.txt +++ b/scripts/requirements.txt @@ -3,4 +3,4 @@ scipy matplotlib scqubits qutip - +ase