gpe.bec2 ======== .. py:module:: gpe.bec2 .. autoapi-nested-parse:: Two component BEC in 1D. Here we demonstrate the dynamics a two component BEC such as 87Rb. Simplifications: * g_aa = g_bb = g_ab = g * m_a = m_b = m * Assume that only total particle number is conserved. Attributes ---------- .. autoapisummary:: gpe.bec2.u Classes ------- .. autoapisummary:: gpe.bec2.StateDFTBase gpe.bec2._StateBase gpe.bec2.GPEMixin gpe.bec2.StateGPEBase gpe.bec2.HOMixin gpe.bec2.StateScaleBase gpe.bec2.StateScaleHO Module Contents --------------- .. py:data:: u .. py:class:: StateDFTBase(**kw) Bases: :py:obj:`gpe.bec.StateDFTBase` Underlying implementation of features needed for IStateDFT. .. py:method:: axpy1(x, a) axpy where a is a scalar. .. py:method:: axpy2(x, a) Version of axpy where a has 2 components. Uses `self._axpy` for each component. .. py:method:: axpy(x, a=1) Version of axpy where a has 2 components. Uses `self._axpy` for each component. .. py:method:: _axpy_numpy(y, x, a=1) Default numpy version of axpy. .. py:method:: apply_V(Va, Vb, Vab, exp=False) .. py:method:: fill(value) .. py:method:: braket_GPU(y) Return the dot product of `self.conjugate()` with `y`. :param y: The states for which the braket will be computed. :type y: State .. py:method:: get_Vext_mu_GPU() Return `(Va, Vb, Vab)` with the chemical potential subtracted if initializing. The chemical potential should be subtracted if initializing or minimizing to get the initial state. Minimization might be done with imaginary time evolution, which should be done with negative times. Thus, we check for `self.initializing`, `self.t < 0`, and make sure that `self.mu` is valid. .. py:method:: get_Vext_GPU() Return (Va, Vb, Vab), the external potentials. .. py:class:: _StateBase(basis=None, Nxyz=(2**5, 2**5, 2**5), Lxyz=(30 * u.micron, 50 * u.micron, 50 * u.micron), symmetric_grid=False, twist=None, v_x=0, kwz2=None, t=0, hbar=u.hbar, ms=(u.m, u.m), mus=None, x_TF=None, cooling_phase=1.0, constraint='N', PGPE=False, mu_Bs=(0 * u.mu_B, 0 * u.mu_B), k_r=0.0, rotating_phases=False, basis_args={}, t_final=None) Bases: :py:obj:`StateDFTBase` Two component state. This state is designed for a spin-orbit (SO) coupled BEC in a harmonic confining potential with the SO field along the x-axis, and a possible magnetic field gradient also along the x-axis. :param basis: If provided, then this is used to specify the basis. If not, then a periodic basis is constructed from `Nxyz` and `Lxyz`. :type basis: IBasis, None :param constraint: If 'N', then constrain the total density allowing particles to convert from one species to the other. If 'Nab' Then independently constrain `Na` and `Nb`. :type constraint: 'N' | 'Nab' :param twist: Twisted boundary conditions for species a and b. This is the twist angle. :type twist: (float, float) :param x_TF: Thomas Fermi "radius" for setting the initial state. The initial state will be set so that the density will fall to zero at this point `x=x_TF` (with y and z in the middle of the trap). If `None`, then we default to 80% of the length. :type x_TF: float :param PGPE: If `True`, then the code implements the PGPE, imposing a cutoff at `kc` :type PGPE: bool .. py:attribute:: t0 :value: 63.50779925891489 .. py:attribute:: basis :value: None .. py:attribute:: data .. py:attribute:: t :value: 0 .. py:attribute:: hbar :value: 1.0 .. py:attribute:: ms .. py:attribute:: mus :value: None .. py:attribute:: x_TF :value: None .. py:attribute:: cooling_phase :value: 1.0 .. py:attribute:: constraint :value: 'N' .. py:attribute:: PGPE :value: False .. py:attribute:: k_r :value: 0.0 .. py:attribute:: twist .. py:attribute:: v_x :value: 0 .. py:attribute:: kwz2 :value: None .. py:attribute:: mu_Bs .. py:property:: rotating_phases Get rotating_phase. .. py:attribute:: _rotating_phases :value: False .. py:method:: init() Initialize the state. This method defines the basis positions, momenta, etc. for use later on. We define these here rather than in the constructor `__init__()` so that the user can change them later and the reinitialize the state. We also call this function from the `pre_evolve_hook()` so that it is called before any evolution takes place. For this reason, we should not modify the state here. .. py:property:: x Flat x abscissa as a numpy array. .. py:method:: set_initial_state(mus=None, x_TF=None) Set the state using the Thomas Fermi (TF) approximation from either `mus` or `x_TF` (pick only one or the other). :param mus: Fixed chemical potentials. :type mus: (float, float) :param x_TF: Position defining the Thomas Fermi "radius". (The external potential is evaluated at this position and this is used to set `mu`.) :type x_TF: float .. py:method:: get_mus_from_Vs_TF(Vs_TF) Return the corrected chemical potential from Vs_TF. In some cases, the chemical potential may differ from the value of the external potentials at V(x_TF) due to kinetic energy shifts (in the SOC case for example) or due to the energy of radial excitations (see the tube codes). This function adds the appropriate correction. :param Vs_TF: External potentials at the Thomas Fermi "radius". (The external potential is evaluated at this position and this is used to get `mus`.) :type Vs_TF: (float, float) .. py:method:: get_Vs_TF_from_mus(mus) Return Vs_TF from the chemical potentials mus. :param mus: Physical chemical potentials (i.e. what you would pass to the minimizer). :type mus: (float, float) .. py:method:: get_Vs_TF(x_TF, V_ext=None) Return the Thomas Fermi chemical potential at x_TF. :param x_TF: Position defining the Thomas Fermi "radius". (The external potential is evaluated at this position and this is used to get `mu`.) :type x_TF: float .. py:method:: _get_kx2(kx=None, Lx=None) Return the effective `kx**2` and `twist_phase_x` for the kinetic energy matrix. This will be passed as the `kx2` and `twist_phase_x` arguments to the basis.lagrangian function. .. py:method:: get_xyz_GPU() .. py:property:: xyz .. py:method:: get_metric_GPU() .. py:property:: metric .. py:method:: get_psi_GPU() Return the physical wavefunction (applying any twist). .. py:method:: set_psi(psi) Set the state from a physical wavefunction (removing any twist). .. py:method:: get_density_GPU() .. py:method:: get_V_GPU() Return the complete potentials `(Va, Vb, Vab)` - internal and external. .. py:property:: dim Dimension of the state. .. py:property:: shape Shape of the state Nxyz (not including components). .. py:property:: bcast Return a set of indices suitable for broadcasting masses etc. .. py:property:: E_max Return the maximum kinetic energy in the basis. This is useful when using evolvers as convergence should be obtained when the time-step is roughly:: dt = 0.1 * state.hbar / state.E_max See `t_scale`. .. py:property:: t_scale Return the smallest time-scale for the problem. Evolvers - especially the ABM evolvers - should use a `dt=0.1*t_scale` or so. If much smaller `dt` values are required for convergence, then it usually indicates that your lattice spacing is too larger. Likewise, if you can get away with much larger `dt` values, then your lattice spacing might be unnecessarily small. Required by `Simulation`. .. py:method:: get_Ns_GPU() .. py:method:: get_N_GPU() .. py:method:: integrate_GPU(a) Integrate over both components, but do not sum. .. py:method:: compute_dy_dt(dy, subtract_mu=True) Return `dy_dt` storing the results in `dy`. :param subtract_mu: If `True`, then subtract the chemical potential such that `dy_dt` is orthogonal to the original state `y`. This will minimize the evolution of the overall phase during real-time evolution (which can reduce numerical errors) and will ensure that evolution under imaginary or complex time will preserve particle number. This should not be set if computing physical energy of the state, however, which is why it is a parameter. :type subtract_mu: bool .. py:attribute:: linear :value: False .. py:method:: apply_exp_K(dt) Apply $e^{-i K dt}$ in place. .. py:method:: apply_exp_V(dt, state) Apply $e^{-i V dt}$ in place using `state` for any nonlinear dependence in V. (Linear problems should ignore `state`.) .. py:method:: normalize(s=None) Normalize the state, return the scale factors and number `(s, N)`. .. py:method:: apply_laplacian(factor, exp=False, **_kw) Apply the laplacian multiplied by `factor` to the state. :param factor: The result will be multiplied by this factor. :type factor: array-like :param exp: If `True` then `exp(factor*laplacian)(y)` will be computed instead. :type exp: bool .. py:method:: get_energy(**kw) .. py:method:: get_Hy(subtract_mu=False) .. py:method:: get_energy_density(a=True, b=True, ab=True) Return the energy density. :param a: Include the energies associated with this species. This includes the kinetic energy and the self-interaction, but no inter-species interactions. :type a: bool :param b: Include the energies associated with this species. This includes the kinetic energy and the self-interaction, but no inter-species interactions. :type b: bool :param ab: Include only the interaction energy. :type ab: bool .. py:method:: get_mu() Compute the chemical potential for convenience only. .. py:method:: plot(log=False) .. py:method:: plot_k() Plot the momentum distribution. .. py:method:: evolve_to(t_end, dt_t_scale=0.2, Evolver=evolvers.EvolverABM, evolve_steps=200, callback=None) Evolve state to `self.t = t_end`. :param t_end: Evolves the state for `t_end` time units. Remember to convert times by `self.t_unit`. :type t_end: float :param dt_t_scale: Scales dt. :type dt_t_scale: float :param Evolver: Pick the relevent evolver, either `evolvers.EvolverABM` or `evolvers.EvolverSplit`. :type Evolver: IEvolver :param evolve_steps: Number of evolution steps between callbacks if defined. :type evolve_steps: float :param callback: Any function that takes the state as an argument, like plotting or visualization. :type callback: function .. py:class:: GPEMixin Mixing providing the 2-component GPE equation of state. .. py:method:: init() .. py:method:: get_Vint_GPU(state=None) Return (Va, Vb), the "internal" mean-field potential. .. py:method:: get_Eint(state=None) Return the "internal" mean-field energy-densities (Ea, Eb, Eab). This version implements the standard GPE where the energy-density has $gn^2/2$. The method get_Vint() should return the appropriate derivative of this. .. py:method:: get_ns_TF(Vs_TF, V_ext=None, gs=None, state=None) Return the Thomas Fermi density profile n from V_TF. :param V_TF: Value of V(x_TF) where the density should vanish in the TF limit. :type V_TF: float :param Note: :type Note: Ignores `gab`. .. py:method:: _get_n_TF(V_TF, V_ext, g, state=None) Helper the total TF density for a single component. :param V_TF: Value of V(x_TF) where the density should vanish in the TF limit. :type V_TF: float .. py:class:: StateGPEBase(gs=None, **kw) Bases: :py:obj:`GPEMixin`, :py:obj:`_StateBase` Two component state with GPE equation of state. .. py:attribute:: gs :value: (1.0, 1.0, 1.0) .. py:class:: HOMixin(ws=None, Omega=None, delta=None, **kw) Bases: :py:obj:`gpe.utils.GPUHelper` Helper mixin class for backwards compatibility. This provides a default get_Vext_GPU(). .. py:attribute:: _Vext :value: None .. py:method:: init() .. py:method:: get_Vext_GPU(state=None) Return (V_a, V_b, V_ab), the external potentials. .. py:class:: StateScaleBase(basis=None, Nxyz=(2**5, 2**5, 2**5), Lxyz=(30 * u.micron, 50 * u.micron, 50 * u.micron), symmetric_grid=False, twist=None, v_x=0, kwz2=None, t=0, hbar=u.hbar, ms=(u.m, u.m), mus=None, x_TF=None, cooling_phase=1.0, constraint='N', PGPE=False, mu_Bs=(0 * u.mu_B, 0 * u.mu_B), k_r=0.0, rotating_phases=False, basis_args={}, t_final=None) Bases: :py:obj:`_StateBase` This state implements the scaling from Castin and Dum. To use this class, you must provide `get_lambdas()` which defines the scaling factors at the current time. .. py:method:: init() Initialize the state. This method defines the basis positions, momenta, etc. for use later on. We define these here rather than in the constructor `__init__()` so that the user can change them later and the reinitialize the state. We also call this function from the `pre_evolve_hook()` so that it is called before any evolution takes place. For this reason, we should not modify the state here. .. py:method:: _get_scale_factor_and_phase_GPU() .. py:method:: get_lambdas(t=None) Return `(lams, dlams, ddlams)`: the scale factors and derivatives. These should be computed at time `t` which should be `self.t` if `t` is None. There should be exactly `self.dim` scale factors - one for each dimension. If a dimension is not to be scaled, then the scale factor should be kept constant at 1. .. py:method:: get_psi_GPU() Return the physical wavefunction (applying any twist). .. py:method:: set_psi(psi) Set the state from a physical wavefunction (removing any twist). .. py:method:: get_xyz_GPU() .. py:method:: get_metric_GPU() .. py:method:: apply_laplacian(factor, exp=False, **_kw) Apply the laplacian multiplied by `factor` to the state. :param factor: The result will be multiplied by this factor. :type factor: array-like :param exp: If `True` then `exp(factor*laplacian)(y)` will be computed instead. :type exp: bool .. py:method:: get_V_GPU() Return the complete potential `V` - internal and external. This version includes the correction from the scaled coordinates. .. py:class:: StateScaleHO(basis=None, Nxyz=(2**5, 2**5, 2**5), Lxyz=(30 * u.micron, 50 * u.micron, 50 * u.micron), symmetric_grid=False, twist=None, v_x=0, kwz2=None, t=0, hbar=u.hbar, ms=(u.m, u.m), mus=None, x_TF=None, cooling_phase=1.0, constraint='N', PGPE=False, mu_Bs=(0 * u.mu_B, 0 * u.mu_B), k_r=0.0, rotating_phases=False, basis_args={}, t_final=None) Bases: :py:obj:`StateScaleBase` This state implements the scaling from Castin and Dum for harmonic oscillators. This provides the required get_lambdas() function but now requires the user specify the time-dependence of the trapping frequencies in get_ws(). .. py:method:: get_ws(t=None) Return the trapping frequencies at time t. .. py:method:: init() Initialize the state. This method defines the basis positions, momenta, etc. for use later on. We define these here rather than in the constructor `__init__()` so that the user can change them later and the reinitialize the state. We also call this function from the `pre_evolve_hook()` so that it is called before any evolution takes place. For this reason, we should not modify the state here. .. py:method:: get_lambdas(t=None) Return `(lams, dlams, ddlams)`: the scale factors and derivatives. These should be computed at time `t` which should be `self.t` if `t` is None. There should be exactly `self.dim` scale factors - one for each dimension. If a dimension is not to be scaled, then the scale factor should be kept constant at 1. .. py:method:: _rhs(q, t) RHS for lambda(t) ODE.