small changes
This commit is contained in:
parent
eade0f6f2e
commit
ac4ff1d16b
39
scripts/ase_spacegroup.py
Normal file
39
scripts/ase_spacegroup.py
Normal file
@ -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()
|
||||
50
scripts/bz.py
Normal file
50
scripts/bz.py
Normal file
@ -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()
|
||||
@ -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)
|
||||
|
||||
@ -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")
|
||||
|
||||
@ -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`.
|
||||
|
||||
@ -3,4 +3,4 @@ scipy
|
||||
matplotlib
|
||||
scqubits
|
||||
qutip
|
||||
|
||||
ase
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user