Adds V_interaction for 3D. Updates builds to new release.

PKGBUILD

@@ -1,7 +1,7 @@
 # Maintainer: Daniel Scheiermann  <[email protected]>
 _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}"

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",
                                   ]},

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,

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}")

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,

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),