Add function

- add electorn_density
- deprecate areal electron density
- use trapz for samples integrating
- add factor_q

pyqhe/equation/poisson.py

@@ -42,24 +42,16 @@ def __init__(self, grid: np.ndarray, charge_density: np.ndarray,
  def calc_poisson(self, **kwargs):
  """Calculate electric field."""
 
- def righthand(z, y):
- # interpolate
- sigma = np.interp(z, self.grid, self.charge_density)
- return sigma * -const.q
-
- sol = solve_ivp(righthand, (self.grid[0], self.grid[-1]), [0],
- t_eval=self.grid,# method='DOP853',
- **kwargs)
- # divide dielectric `eps`
- # sol = odeint(righthand, [0], self.grid, tfirst=False)
- # self.e_field = sol.flatten() / self.eps
- self.e_field = sol.y.flatten() / self.eps
+ d_z = cumulative_trapezoid(const.q * self.charge_density,
+ self.grid,
+ initial=0)
+ self.e_field = d_z / self.eps
  # integral the potential
  self.v_potential = cumulative_trapezoid(self.e_field,
  self.grid,
  initial=0)
 
- return self.v_potential
+ return self.v_potential * const.q
 
 
 class PoissonFDM(PoissonSolver):
@@ -87,9 +79,7 @@ def calc_poisson(self, **kwargs):
  tmp = (np.hstack(
  ([0.0], self.charge_density[:-1])) + self.charge_density
  ) # using trapezium rule for integration (?).
- tmp *= (
- const.q / 2.0
- )
+ tmp *= (const.q / 2.0)
  # Note: sigma is a number density per unit area, needs to be converted
  # to Couloumb per unit area
  tmp[0] = F0

pyqhe/equation/schrodinger.py

@@ -129,7 +129,8 @@ def calc_esys(self, **kwargs):
  wave_function = []
  for energy in eig_val:
  self._psi_iteration(energy)
- norm = np.linalg.norm(self.psi)
+ norm = np.sqrt(np.trapz(self.psi * np.conj(self.psi),
+ x=self.grid)) # l2-norm
  wave_function.append(self.psi / norm)
 
  return eig_val, np.asarray(wave_function)

pyqhe/pytest/test_quantum_well.py

@@ -12,7 +12,30 @@
 stick_nm = [20, 30, 40, 45, 50, 60, 70, 80]
 stick_list = [
  0.66575952, 0.97971301, 1.36046512, 1.64324592, 1.88372093, 2.59401286,
- 3.25977239, 3.98119743] # unit in l_B
+ 3.25977239, 3.98119743
+] # unit in l_B
+
+
+def factor_q_fh(thickness, q):
+ """Using the Fang-Howard variational wave function to describe the electron
+ wave function in the perpendicular direction.
+
+ Args:
+ thickness: physical meaning of electron layer thickness.
+ """
+ b = 1 / thickness
+ return (1 + 9 / 8 * q / b + 3 / 8 * (q / b)**2) * (1 + q / b)**(-3)
+
+
+def factor_q(grid, wave_func, q):
+ """Calculate `F(q)` in reduced Coulomb interaction."""
+ # make 2-d coordinate matrix
+ z_1, z_2 = np.meshgrid(grid, grid)
+ exp_term = np.exp(-q * np.abs(z_1 - z_2))
+ wf2_z1, wf2_z2 = np.meshgrid(wave_func**2, wave_func**2)
+ factor_matrix = wf2_z1 * wf2_z2 * exp_term
+ # integrate using the composite trapezoidal rule
+ return np.trapz(np.trapz(factor_matrix, grid), grid)
 
 
 def calc_omega(thickness=10, tol=5e-5):
@@ -27,32 +50,69 @@ def calc_omega(thickness=10, tol=5e-5):
 
  model = Structure1D(layer_list, temp=10, dz=0.2)
  # instance of class SchrodingerPoisson
- schpois = SchrodingerPoisson(model,
- poisolver=PoissonFDM,
- # quantum_region=(255 - 20, 255 + thickness + 30),
- )
+ schpois = SchrodingerPoisson(
+ model,
+ poisolver=PoissonODE,
+ # quantum_region=(255 - 20, 255 + thickness + 30),
+ )
  # perform self consistent optimization
  res, loss = schpois.self_consistent_minimize(tol=tol)
  if loss > tol:
  res, loss = schpois.self_consistent_minimize(tol=tol)
+ # plot 2DES areal electron density
+ plt.plot(res.grid, res.sigma * thickness * 1e14)
+ 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[schpois.quantum_mask],
- res.wave_function[0][schpois.quantum_mask],
+ res.electron_density[schpois.quantum_mask],
  p0=[1, np.mean(res.grid), 10])
- plt.plot(res.grid[schpois.quantum_mask], res.wave_function[0][schpois.quantum_mask], label='Self consistent')
+ 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[schpois.quantum_mask],
- gaussian(res.grid[schpois.quantum_mask], *popt),
- label='Norm fit')
+ wf2[schpois.quantum_mask],
+ label=r'$|\Psi(z)|^2$')
+ plt.plot(res.grid[schpois.quantum_mask],
+ charge_distribution[schpois.quantum_mask],
+ 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('Charge distribution')
+ 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 popt[2], res
+ return sigma, res
+
+
+# %%
+_, res = calc_omega(50)
+res.plot_quantum_well()
 # %%
 thickness_list = np.linspace(10, 65, 20)
 res_list = []
@@ -61,19 +121,20 @@ def gaussian(x, amplitude, mean, stddev):
  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(omega_list, line(omega_list, *popt1))
 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[15].plot_quantum_well()
 # %%
-plt.plot(res_list[15].grid, res_list[15].wave_function[0])
-# %%

pyqhe/schrodinger_poisson.py

@@ -26,6 +26,8 @@ def __init__(self) -> None:
  # electric field properties
  self.v_potential = None
  self.e_field = None
+ # Accumulate electron density
+ self.electron_density = None
 
  def plot_quantum_well(self):
  """Plot dressed conduction band of quantum well, and electrons'
@@ -108,17 +110,16 @@ def __init__(self,
  def _calc_net_density(self, n_states, wave_func):
  """Calculate the net charge density."""
 
- # firstly, convert to areal charge density
- grid_dist = np.diff(self.grid)
- grid_dist = np.append(grid_dist, grid_dist[-1]) # adjust shape
- doping_2d = self.doping * grid_dist
  # Accumulate electron areal density in the subbands
  elec_density = np.zeros_like(self.grid[self.quantum_mask])
  for i, distri in enumerate(n_states):
  elec_density += distri * wave_func[i] * np.conj(wave_func[i])
+ # normalize by electric neutrality
+ norm = np.trapz(self.doping, self.grid) / np.trapz(
+ elec_density, self.grid)
+ elec_density *= norm
  # Let dopants density minus electron density
- net_density = doping_2d
- net_density[self.quantum_mask] = doping_2d[self.quantum_mask] - elec_density
+ net_density = self.doping[self.quantum_mask] - elec_density
 
  return net_density
 
@@ -172,6 +173,8 @@ def self_consistent_minimize(self,
  if i and loss < tol:
  break
  # save optimize result
+ # optimal_index = np.argmin(loss_list[1:])
+ # self.params = param_list[optimal_index]
  res = OptimizeResult()
  res.params = self.params
  res.grid = self.grid
@@ -192,5 +195,10 @@ def self_consistent_minimize(self,
  self.poi_solver.charge_density = res.sigma
  self.poi_solver.calc_poisson()
  res.e_field = self.poi_solver.e_field
+ # Accumulate electron areal density in the subbands
+ res.electron_density = np.zeros_like(self.grid[self.quantum_mask])
+ for i, distri in enumerate(res.n_states):
+ res.electron_density += distri * res.wave_function[i] * np.conj(
+ res.wave_function[i])
 
  return res, loss

pyqhe/utility/fermi.py

@@ -4,10 +4,10 @@
 import pyqhe.utility.constant as const
 
 
-def calc_meff_state(wave_function, cb_meff):
+def calc_meff_state(grid, wave_function, cb_meff):
  """Calculate subband effective mass."""
  wave_function = np.asarray(wave_function)
- return 1.0 / np.sum(wave_function**2 / cb_meff, axis=1)
+ return 1.0 / np.trapz(wave_function**2 / cb_meff, x=grid, axis=1)
 
 
 class FermiStatistic:
@@ -30,7 +30,7 @@ def __init__(
  # Calculate 2d doping density
  # integrate 3-d density along axis z
  # in discrete method, just sum over r'n_3d[i] * (grid[i + 1] - grid[i])'
- self.doping_2d = np.trapz(doping, grid)
+ self.max_occupation = np.trapz(doping, grid)
  # Cache parameters
  self.fermi_energy = None
  self.n_states = None
@@ -46,14 +46,15 @@ def integral_fermi_dirac(self, energy, fermi_energy, temp):
  def fermilevel_0k(self, eig_val, wave_function):
  """Calculate Fermi level at 0 K."""
 
- self.meff_state = calc_meff_state(wave_function, self.cb_meff)
+ self.meff_state = calc_meff_state(self.grid, wave_function,
+ self.cb_meff)
  # list all fermi energy candidate
  candidate_fermi_energy = []
  for i, _ in enumerate(eig_val):
  accumulate_energy = np.sum(eig_val[:i + 1])
  candidate_fermi_energy.append(
- (self.doping_2d * const.hbar**2 * np.pi / self.meff_state[i] +
- accumulate_energy) / (i + 1))
+ (self.max_occupation * const.hbar**2 * np.pi /
+ self.meff_state[i] + accumulate_energy) / (i + 1))
  # check true Fermi energy
  fermi_idx = np.argwhere(
  (np.asarray(candidate_fermi_energy) - eig_val) < 0)
@@ -82,7 +83,7 @@ def func(f_energy):
  csb_meff * self.integral_fermi_dirac(eig_v, f_energy, temp)
  for eig_v, csb_meff in zip(eig_val, self.meff_state)
  ]
- return self.doping_2d - np.sum(dist) / (const.hbar**2 * np.pi)
+ return self.max_occupation - np.sum(dist) / (const.hbar**2 * np.pi)
 
  f_energy_0k, _ = self.fermilevel_0k(eig_val, wave_function)
  # sol = optimize.root_scalar(func, x0=f_energy_0k, method='halley', **kwargs)