From f71c4f73dcb0766c0ad94d3288a90bf83d7f4688 Mon Sep 17 00:00:00 2001 From: Scheiermann Date: Fri, 8 Jan 2021 15:23:17 +0100 Subject: [PATCH] Adds V_interaction for 3D. Updates builds to new release. --- PKGBUILD | 2 +- setup.py | 2 +- supersolids/MayaviAnimation.py | 2 +- supersolids/Schroedinger.py | 31 ++++++++++++++++++++++++++----- supersolids/functions.py | 12 ++++++++++++ supersolids/main.py | 15 ++++++++++----- 6 files changed, 51 insertions(+), 13 deletions(-) diff --git a/PKGBUILD b/PKGBUILD index d2fa8de..c1cfcad 100644 --- a/PKGBUILD +++ b/PKGBUILD @@ -1,7 +1,7 @@ # Maintainer: Daniel Scheiermann _name=supersolids pkgname=python-${_name} -pkgver=0.1.16 +pkgver=0.1.17 pkgrel=1 pkgdesc="Notes and script to supersolids" url="https://github.com/Scheiermann/${_name}" diff --git a/setup.py b/setup.py index 7e39220..b340cae 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ setup( name="supersolids", - version="0.1.16", + version="0.1.17", packages=["", "supersolids"], package_data={"supersolids": ["results/anim.mp4", ]}, diff --git a/supersolids/MayaviAnimation.py b/supersolids/MayaviAnimation.py index ea3fab9..abd75c2 100644 --- a/supersolids/MayaviAnimation.py +++ b/supersolids/MayaviAnimation.py @@ -93,7 +93,7 @@ def animate(System: Schroedinger.Schroedinger, accuracy: float = 10 ** -6, extent=[*x_lim, *y_lim, *z_lim]) plot_V = mlab.contour3d(System.x_mesh, System.y_mesh, System.z_mesh, System.V_val, - colormap="hot", opacity=System.alpha_V, transparent=True) + colormap="hot", opacity=System.alpha_V, transparent=True) if System.psi_sol_val is not None: plot_psi_sol = mlab.contour3d(System.x_mesh, System.y_mesh, System.z_mesh, System.psi_sol_val, diff --git a/supersolids/Schroedinger.py b/supersolids/Schroedinger.py index 09bd5c1..8f4c987 100644 --- a/supersolids/Schroedinger.py +++ b/supersolids/Schroedinger.py @@ -37,6 +37,7 @@ def __init__(self, resolution: int, timesteps: int, L: float, dt: float, g: floa dim: int = 3, psi_0: Callable = functions.psi_gauss_3d, V: Callable = functions.v_harmonic_3d, + V_interaction: Callable = None, psi_sol: Callable = functions.thomas_fermi_3d, mu_sol: Callable = functions.mu_3d, alpha_psi: float = 0.8, @@ -131,7 +132,18 @@ def __init__(self, resolution: int, timesteps: int, L: float, dt: float, g: floa # here a number (U) is multiplied elementwise with an (1D, 2D or 3D) array (k_squared) self.H_kin = np.exp(self.U * (0.5 * self.k_squared) * self.dt) - # attributes for animation + if V_interaction is None: + # For no interaction the identity is needed with respect to 2D * 2D (array with 1.0 everywhere) + self.V_k_val = np.full(self.psi_val.shape, 1.0) + else: + self.V_k_val = V_interaction(kx_mesh, ky_mesh, kz_mesh, g=self.g) + + density = np.abs(self.psi_val) ** 2.0 + density_k = np.fft.fftn(density) + density_k_interact = self.V_k_val * density_k + U_interaction = np.fft.ifftn(density_k_interact) + + # attributes for animation self.t = 0.0 self.alpha_psi = alpha_psi @@ -153,8 +165,12 @@ def get_norm(self, p: float = 2.0) -> float: def time_step(self): # Here we use half steps in real space, but will use it before and after H_kin with normal steps + + # Calculate the interaction by appling it to the density in k-space (transform back and forth) + density = np.abs(self.psi_val) ** 2.0 + U_interaction = np.fft.ifftn(self.V_k_val * np.fft.fftn(density)) # update H_pot before use - H_pot = np.exp(self.U * (self.V_val + self.g * np.abs(self.psi_val) ** 2.0) * (0.5 * self.dt)) + H_pot = np.exp(self.U * (self.V_val + self.g * density + U_interaction) * (0.5 * self.dt)) # multiply element-wise the (1D, 2D or 3D) arrays with each other self.psi_val = H_pot * self.psi_val @@ -164,8 +180,10 @@ def time_step(self): self.psi_val = self.H_kin * self.psi_val self.psi_val = np.fft.ifftn(self.psi_val) - # update H_pot before use - H_pot = np.exp(self.U * (self.V_val + self.g * np.abs(self.psi_val) ** 2.0) * (0.5 * self.dt)) + # update H_pot, density, U_interaction before use + density = np.abs(self.psi_val) ** 2.0 + U_interaction = np.fft.ifftn(self.V_k_val * np.fft.fftn(density)) + H_pot = np.exp(self.U * (self.V_val + self.g * density + U_interaction) * (0.5 * self.dt)) # multiply element-wise the (1D, 2D or 3D) arrays with each other self.psi_val = H_pot * self.psi_val @@ -182,4 +200,7 @@ def time_step(self): self.E = self.mu - 0.5 * self.g * psi_quadratic_integral print(f"mu: {self.mu}") - print(f"E: {self.E}, E_sol: {self.mu_sol - 0.5 * self.g * psi_quadratic_integral}") + if self.g != 0: + print(f"E: {self.E}, E_sol: {self.mu_sol - 0.5 * self.g * psi_quadratic_integral}") + else: + print(f"E: {self.E}") diff --git a/supersolids/functions.py b/supersolids/functions.py index 293ed93..1c11a09 100644 --- a/supersolids/functions.py +++ b/supersolids/functions.py @@ -289,6 +289,7 @@ def thomas_fermi_3d(x, y, z, g: float = 0.0): else: return None + def mu_1d(g: float = 0.0): # mu is the chemical potential mu = ((3.0 * g) / (4.0 * np.sqrt(2.0))) ** (2.0 / 3.0) @@ -329,6 +330,17 @@ def v_harmonic_3d(x, y, z, alpha_y: float = 1.0, alpha_z: float = 1.0): return 0.5 * (x ** 2 + (alpha_y * y) ** 2 + (alpha_z * z) ** 2) +def dipol_dipol_interaction(kx_mesh: float, ky_mesh: float, kz_mesh: float, + g: float = 1.0, d: float = 1.0, epsilon_dd: float = 1.0): + k_squared = kx_mesh ** 2.0 + ky_mesh ** 2.0 + kz_mesh ** 2.0 + factor = 3.0 * (kz_mesh ** 2.0) + V_k_val = epsilon_dd * g * (4.0 * np.pi / 3.0) * d ** 2.0 * ((factor / k_squared) - 1.0) + # Remove singularity + V_k_val[np.isnan(V_k_val)] = 0.0 + + return V_k_val + + def camera_func_r(frame: int, r_0: float = 10.0, phi_0: float = 45.0, diff --git a/supersolids/main.py b/supersolids/main.py index 015ef08..c6980fd 100644 --- a/supersolids/main.py +++ b/supersolids/main.py @@ -28,6 +28,7 @@ def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, s=1.1, E=1.0 dim=1, psi_0=functions.psi_gauss_3d, V=functions.v_harmonic_3d, + V_interaction=None, psi_sol=functions.thomas_fermi_3d, mu_sol=functions.mu_3d, alpha_psi=0.8, @@ -52,6 +53,7 @@ def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, s=1.1, E=1.0 dim=dim, psi_0=psi_0, V=V, + V_interaction=V_interaction, psi_sol=psi_sol, mu_sol=mu_sol, alpha_psi=alpha_psi, @@ -103,9 +105,9 @@ def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, s=1.1, E=1.0 resolution: int = 2 ** datapoints_exponent # constants needed for the Schroedinger equation - g = 0.0 + g = 100.0 g_step = 10 - dt = 0.4 + dt = 0.001 # box length in 1D: [-L,L], in 2D: [-L,L, -L,L], , in 3D: [-L,L, -L,L, -L,L] # generators for L, g, dt to compute for different parameters @@ -120,6 +122,8 @@ def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, s=1.1, E=1.0 V_2d = functools.partial(functions.v_harmonic_2d, alpha_y=1.0) V_3d = functools.partial(functions.v_harmonic_3d, alpha_y=1.0, alpha_z=1.0) + V_3d_ddi = functools.partial(functions.dipol_dipol_interaction, d=1.0, epsilon_dd=1.0) + # functools.partial sets all arguments except x, y, z, as multiple arguments for Schroedinger aren't implement yet # psi_0 = functools.partial(functions.psi_0_rect, x_min=-0.25, x_max=-0.25, a=2.0) psi_0_1d = functools.partial(functions.psi_gauss_1d, a=3.0, x_0=2.0, k_0=0.0) @@ -133,20 +137,21 @@ def simulate_case(resolution, timesteps, L, g, dt, imag_time=False, s=1.1, E=1.0 # TODO: get mayavi lim to work # 3D works in single core mode - simulate_case(resolution, timesteps=200, L=L_generator[0], g=g, dt=dt, imag_time=True, + simulate_case(resolution, timesteps=800, L=L_generator[0], g=g, dt=dt, imag_time=True, s=1.1, E=1.0, dim=3, psi_0=psi_0_3d, V=V_3d, + V_interaction=V_3d_ddi, psi_sol=psi_sol_3d, mu_sol=functions.mu_3d, - accuracy=10**-6, + accuracy=10**-8, alpha_psi=0.8, alpha_psi_sol=0.5, alpha_V=0.3, file_name="anim.mp4", x_lim=(-2.0, 2.0), y_lim=(-2.0, 2.0), z_lim=(0, 0.5), - slice_x_index=resolution//2, slice_y_index=resolution//2, # just for mayavi (3D) + slice_x_index=resolution//3, slice_y_index=resolution//3, # just for mayavi (3D) r_func=functools.partial(functions.camera_func_r, r_0=10.0, r_per_frame=0.0), # from here just 2D phi_func=functools.partial(functions.camera_func_phi, phi_0=45.0, phi_per_frame=10.0), z_func=functools.partial(functions.camera_func_r, r_0=20.0, r_per_frame=0.0),