gpe.bec_basic#
Code for simulating a single component BEC.
This is the simplest code in the project, implementing all of the basic features without excess complications. It is somewhat excessively commented to provide a pedagogical introduction for users unfamiliar with the techniques and using python for scientific computing.
A more complete implementation can be found in bec.py, but this relies on external packages to implement, for example, different bases, making it less suitable for teaching.
Classes#
Module Contents#
- class Units[source]#
Units and physical constants.
This class is simply a container for physical constants and units. We have chosen values here that are relevant for 87Rb which is studied at WSU in Peter Engels’ lab.
- class State(Nxyz=(2**5, 2**5, 2**5), Lxyz=(30 * u.micron, 50 * u.micron, 50 * u.micron), ws=np.array([np.sqrt(8) * 126.0, 126.0, 126.0]) * u.Hz, g=4 * np.pi * u.hbar**2 * u.a / u.m, m=u.m, mu=None, hbar=u.hbar, x_TF=None, cooling_phase=1.0, symmetric_grid=False, t=0)[source]#
Bases:
gpe.mixins.StateMixin,pytimeode.mixins.ArrayStateMixinState class implementing the GPE for a single-component BEC.
To provide time-evolution functionality, the pytimeode package requires that you provide a State class that implements the required behaviour of at least one of the IStateFor*Evolvers interfaces. Most of the required behaviour can be provided by the mixins.ArrayStateMixin class which we inherit from as long as we provide an array called data which represents the wavefunction.
On top of this, we need to provide the following three methods used by the pytimeode.evolvers:
compute_dy_dt(dy): Required by interface IStateForABMEvolvers
This computes the rhs of the time-dependent differential equation and is all that is needed by the pytimeode.evolvers.EvolverABM evolver. This is a high-order generic evolver. It requires a rather small step size to ensure that it behaves well, but once the step size is small enough, it is converges to high accuracy. For example: there is often a critical step size above which the evolution will blow up. If you reduce this by a factor of 2, then you will typically obtain accuracy to 5 digits or so, at least for modest evolution lengths.
The ABM evolver requires about 10 copies of the state to start up, and 8 copies of the state to be maintained for evolution. This can be a memory issue if the states are very large.
apply_exp_K(dt), apply_exp_V(dt, state): Required by interface IStateForSplitEvolvers
These two methods allow you to use pytimeode.evolvers.EvolverSplit based on the Trotter decomposition. This is rather specialized to systems like the GPE where the exponential of the kinetic and potential terms can be computed exactly. The ODE solver is not as high order as the ABM method, but is manifestly unitary. Thus, you can often get away with very large step sizes and the system will still behave “reasonably” meaning that the evolution will not be accurate quantitatively, but will be qualitatively reasonable so you can get an idea of what is happening. Unfortunately, to achieve quantitative accuracy, you must generally go to very small step sizes.
Another advantage of the split operator methods is that they require use only 2 or 3 copies of the state, and so are better for large states that might be a memory issue
normalize(): Required by interface IStateWithNormalize
This method allows the evolvers to continually normalize the state, which can improve numerical performance, or be of use when finding initial states.
The reset of the methods in the code play a supporting role to these methods required by the interfaces and for our convenience (for example, plotting the state).
- Parameters:
x_TF (float) – 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.
- init()[source]#
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.
- property E_max[source]#
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.
- property t_scale[source]#
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.
- pre_evolve_hook()[source]#
This method is called by the evolvers at the start of evolution to ensure that the state is properly initialized.
- set_initial_state(mu=None, x_TF=None)[source]#
Set the state using the Thomas Fermi (TF) approximation from either mu or x_TF (pick only one or the other).
- get_mu_TF(x_TF, V_ext=None)[source]#
Return the Thomas Fermi chemical potential at x_TF.
- Parameters:
x_TF (float) – Position defining the Thomas Fermi “radius”. (The external potential is evaluated at this position and this is used to get mu.)
- get_mu_HO(N)[source]#
Return the chemical potential required to set the particle number in the Thomas Fermi (TF) approximation for a 3D harmonic oscillator.
- get_Vext()[source]#
Return the external potential.
Overload this method if you want to change the external potential. If the potential should be time dependent, use self.t which will be updated by the evolvers.
If a chemical potential self.mu is defined, then this is subtracted from Vext. This allows gpe.minimize to find states at a constant chemical potential. Note: for general evolution, it is better not to set the chemical potential as this is automatically set by compute_dy_dt and will then cause get_energy to return the thermodynamic potential instead.
- braket(a, b)[source]#
Return the inner product of the functions a and b.
Note: a is conjugated and the metric is applied.
- get_V()[source]#
Return the complete potential V - internal and external.
This function defines the non-linear interaction of the equations. For the GPE, the energy-density has \(gn^2/2\), so we have the derivative gn here.
- _fftn(y)[source]#
Return the FFT of y over the spatial indices. These versions with an underscore may be overloaded for performance.
- _ifftn(y)[source]#
Return the IFFT of y over the spatial indices. These versions with an underscore may be overloaded for performance.
- compute_dy_dt(dy, subtract_mu=True)[source]#
Return dy_dt storing the results in dy.
- Parameters:
subtract_mu (bool) –
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.
- apply_exp_V(dt, state)[source]#
Apply \(e^{-i V dt}\) in place using state for any nonlinear dependence in V. (Linear problems should ignore state.)
- get_H(subtract_mu=False)[source]#
Return the single-particle Hamiltonian for convenience only.
Note
May be huge!
This is a square matrix whose sides are of dimension prod(Nxyz) which can be huge. It is not needed for general computations, but might be useful for analysis such as exploring the BdG in 1D systems.