Add support for simulate metal plate

Author: allegro0132

pyqhe/core/structure.py

@@ -336,6 +336,7 @@ def __init__(self,
         # dimension as 'universal grid'.
         self._universal_grid = None
         self.dim = None
+        self.bound_dirichlet = None
         # Structure's properties
         self.temp = temp
         self.fi = None
@@ -410,6 +411,10 @@ def _prepare_structure_stroage(self, delta=1, spline_storage=False):
         self.cb_meff = np.broadcast_to(cb_meff, self.dim)
         self.doping = np.broadcast_to(doping, self.dim)
 
+    def add_dirichlet_boundary(self, bound: np.ndarray):
+        if bound.shape == tuple(self.dim):
+            self.bound_dirichlet = bound
+
 
 class Structure3D:
     """Class for modeling 3d material structure.

pyqhe/equation/poisson.py

@@ -85,7 +85,12 @@ def __init__(self,
             self.eps = eps
         else:
             raise ValueError('The dimension of eps is not match.')
-        self.bound_dirichlet = bound_dirichlet
+        if bound_dirichlet is None:
+            self.bound_dirichlet = None
+        elif bound_dirichlet.shape == tuple(self.dim) or len(self.dim) == 1:
+            self.bound_dirichlet = bound_dirichlet
+        else:
+            raise ValueError('The dimension of eps is not match.')
 
     def build_d_matrix(self, loc):
         """Build 1D time independent Schrodinger equation kinetic operator.
@@ -123,10 +128,11 @@ def calc_poisson(self, **kwargs):
 
         if self.bound_dirichlet is not None:
             delta = self.grid[0][1] - self.grid[0][0]
+            # adjust coefficient let matrix looks good
             bound_b = self.bound_dirichlet * self.eps / delta**2
             bound_a = np.zeros_like(b_vec)
-            bound_loc = np.flatnonzero(self.bound_dirichlet)
-            bound_a[bound_loc] = self.eps.flatten()[bound_loc] / delta**2
+            bound_loc = np.flatnonzero(~np.isnan(self.bound_dirichlet))
+            bound_a[bound_loc] = 1 * self.eps.flatten()[bound_loc] / delta**2
             # tensor contraction
             bound_mat = np.diag(bound_a)
             a_mat[bound_loc] = bound_mat[bound_loc]
@@ -162,17 +168,18 @@ def calc_poisson(self, **kwargs):
     top_plate = (yv <= 0.1 + delta / 2) * (yv >= 0.1 - delta / 2)
     bottom_plate = (yv <= -0.1 + delta /2) * (yv >= -0.1 -delta / 2)
     length = (xv <= 0.5) * (xv >= -0.5)
-    bound = np.zeros_like(xv)
+    bound = np.empty_like(xv)
+    bound[:] = np.nan
     bound[top_plate * length] = 1
     bound[bottom_plate * length] = -1
-    sol = PoissonFDM([x, y], np.zeros_like(bound),
-                     np.ones_like(bound) * const.eps0, bound)
+    sol = PoissonFDM([x, y], np.zeros_like(xv),
+                     np.ones_like(xv) * const.eps0, bound)
     # sol = PoissonFDM([x, y], bound * 10,
     #                   np.ones_like(bound) * const.eps0)
     v_p = sol.calc_poisson()
     # v potential
     plt.pcolormesh(xv, yv, v_p)
     plt.show()
     # e field
-    plt.pcolormesh(xv, yv, np.sqrt(sol.e_field[0]**2 +  sol.e_field[1]**2))
+    plt.pcolormesh(xv, yv, np.sqrt(sol.e_field[0]**2 + sol.e_field[1]**2))
 # %%

pyqhe/pytest/test_quantum_well_2d.py

@@ -51,6 +51,17 @@ def calc_omega(thickness=10, tol=5e-5):
     layer_list.append(Layer(20, 0.24, 0.0, name='barrier'))
 
     model = Structure2D(layer_list, width=50, temp=10, delta=1)
+    # add boundary condition
+    grid = model.universal_grid
+    delta = grid[0][1] - grid[0][0]
+    xv, yv = np.meshgrid(*grid, indexing='ij')
+    top_plate = (yv <= 15) * (yv >= 10)
+    bottom_plate = (yv <= 65) * (yv >= 60)
+    bound = np.empty_like(xv)
+    bound[:] = np.nan
+    bound[top_plate] = 0.05
+    # bound[bottom_plate] = 0
+    model.add_dirichlet_boundary(bound)
     # instance of class SchrodingerPoisson
     schpois = SchrodingerPoisson(
         model,
@@ -74,50 +85,6 @@ def calc_omega(thickness=10, tol=5e-5):
                            antialiased=False)
     plt.show()
 
-    # # fit a normal distribution to the data
-    # def gaussian(x, amplitude, mean, stddev):
-    #     return amplitude * np.exp(-(x - mean)**2 / 2 / stddev**2)
-
-    # popt, _ = optimize.curve_fit(gaussian,
-    #                              res.grid,
-    #                              res.electron_density,
-    #                              p0=[1, np.mean(res.grid), 10])
-    # wf2 = res.wave_function[0] * np.conj(
-    #     res.wave_function[0])  # only ground state
-    # symmetry_axis = popt[1]
-    # # calculate standard deviation
-    # # the standard deviation of the charge distribution from its center, from PRL, 127, 056801 (2021)
-    # charge_distribution = res.electron_density / np.trapz(
-    #     res.electron_density, res.grid)
-    # sigma = np.sqrt(
-    #     np.trapz(charge_distribution * (res.grid - symmetry_axis)**2, res.grid))
-    # plt.plot(res.grid,
-    #          wf2,
-    #          label=r'$|\Psi(z)|^2$')
-    # plt.plot(res.grid,
-    #          charge_distribution,
-    #          label='Charge distribution',
-    #          color='r')
-    # plt.axvline(symmetry_axis - sigma, ls='--', color='y')
-    # plt.axvline(symmetry_axis + sigma, ls='--', color='y')
-    # plt.xlabel('Position (nm)')
-    # plt.ylabel('Distribution')
-    # plt.legend()
-    # plt.show()
-    # # plot factor_q verses q
-    # q_list = np.linspace(0, 1, 20)
-    # f_fh = []
-    # f_self = []
-    # for q in q_list:
-    #     f_fh.append(factor_q_fh(thickness, q))
-    #     f_self.append(factor_q(res.grid, res.wave_function[0], q))
-    # plt.plot(q_list, f_fh, label='the Fang-Howard wave function')
-    # plt.plot(q_list, f_self, label='self-consistent wave function', color='r')
-    # plt.legend()
-    # plt.xlabel('wave vector q')
-    # plt.ylabel(r'$F(q)$')
-    # plt.show()
-
     return res
 
 
@@ -140,29 +107,14 @@ def calc_omega(thickness=10, tol=5e-5):
 plt.show()
 plt.plot(res.grid[0], res.sigma[:, shape[1]] * 20 * 1e14)
 # %%
-thickness_list = np.linspace(10, 80, 30)
-res_list = []
-omega_list = []
-for thick in thickness_list:
-    omega, res = calc_omega(thick)
-    res_list.append(res)
-    omega_list.append(omega)
-
-
-# %%
-def line(x, a, b):
-    return a * np.asarray(x) + b
-
-
-popt1, _ = optimize.curve_fit(line, omega_list[:3], thickness_list[:3])
-# %%
-plt.plot(np.asarray(omega_list), thickness_list, label='PyQHE')
-plt.plot(np.asarray(stick_list) * 7.1, stick_nm, label='Shayegan')
-plt.xlabel(r'Layer thickness $\bar{\omega}$ (nm)')
-plt.ylabel(r'Geometry thickness $\omega$ (nm)')
-plt.legend()
-# %%
-res_list[-1].plot_quantum_well()
-# %%
-plt.plot(res.grid, res.params)
+xv, yv = np.meshgrid(*res.grid, indexing='ij')
+fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
+surf = ax.plot_surface(xv,
+                       yv,
+                       res.v_potential,
+                       cmap=cm.coolwarm,
+                       linewidth=0,
+                       antialiased=False)
+plt.show()
+plt.plot(res.grid[1], res.v_potential[shape[0]])
 # %%

pyqhe/schrodinger_poisson.py

@@ -5,7 +5,7 @@
 from pyqhe.equation.schrodinger import SchrodingerSolver, SchrodingerShooting
 from pyqhe.equation.poisson import PoissonSolver, PoissonODE, PoissonFDM
 from pyqhe.utility.fermi import FermiStatistic
-from pyqhe.core.structure import Structure1D
+from pyqhe.core.structure import Structure2D
 
 
 class OptimizeResult:
@@ -72,7 +72,7 @@ class SchrodingerPoisson:
     """
 
     def __init__(self,
-                 model: Structure1D,
+                 model: Structure2D,
                  schsolver: SchrodingerSolver = SchrodingerShooting,
                  poisolver: PoissonSolver = PoissonFDM,
                  learning_rate=0.5,
@@ -87,6 +87,8 @@ def __init__(self,
         self.doping = model.doping  # doping profile
         # load grid configure
         self.grid = model.universal_grid
+        # load boundary condition
+        self.bound_dirichlet = model.bound_dirichlet
         # Setup Quantum region
         # if quantum_region is not None and len(quantum_region) == 2:
         #     self.quantum_mask = (self.grid > quantum_region[0]) * (
@@ -102,7 +104,8 @@ def __init__(self,
         self.fermi_util = FermiStatistic(self.grid,
                                          self.cb_meff,
                                          self.doping)
-        self.poi_solver = poisolver(self.grid, self.doping, self.eps)
+        self.poi_solver = poisolver(self.grid, self.doping, self.eps,
+                                    self.bound_dirichlet)
         # accumulate charge density
         self.accumulate_q = self.doping
         for grid in self.grid[::-1]: