import mmf_setup
mmf_setup.nbinit()
This cell adds /home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:
- Choose "Trust Notebook" from the "File" menu.
- Re-execute this cell.
- Reload the notebook.
Travelling Waves#
A traveling wave with velocity \(v\) has the following form:
\[\begin{split}
\psi(x, t) = \sqrt{\rho}(x - vt) e^{\I\phi(x-vt) - \I\omega t}, \\
\I\hbar \dot{\psi} = -\frac{\hbar^2\psi''}{2m} + (g\rho - \mu) \psi.\\
-\I\hbar v\psi' + \hbar\omega\psi = -\frac{\hbar^2\psi''}{2m} + (g\rho - \mu) \psi.
\end{split}\]
By appropriately adjusting \(\mu\) we may set \(\omega = 0\), hence we may drop the factor \(e^{-\I\omega t}\) from the wavefunction:
\[
\psi(x - vt) = \sqrt{\rho}(x - vt) e^{\I\phi(x-vt)}.
\]
The last equation above now represents a stationary solution in a moving frame with coordinate \(y = x-vt\):
\[
0 = -\frac{\hbar^2\psi''(y)}{2m} + \I\hbar v\psi'(y) + [g\rho(y)- \mu]\psi(y)
= \left[\frac{\op{p}^2}{2m} -\op{p} v + g\rho(y) - \mu\right]\psi(y)
= \left[\frac{(\op{p} - p_v)^2}{2m} - \frac{p_v^2}{2m} + g\rho(y) - \mu\right]\psi(y)
\]
where \(p_v = mv\). Finally, we may recast this in terms of the original GPE by transforming \(\psi(y)\) as follows:
\[\begin{split}
\psi(y) = e^{\I p_v x/\hbar}\tilde{\psi}(y)\\
(\op{p}-p_v)\psi(y) = e^{-\I p_v x/\hbar}\op{p}\tilde{\psi}(y)\\
\left[\frac{\op{p}^2}{2m} + g\rho(y) - \left(\mu + \frac{p_v^2}{2m}\right)\right]\tilde{\psi}(y)
\end{split}\]
If the original function \(\psi(y+L) = \psi(y)\) was periodic, then \(\tilde{\psi}(y+L) = e^{\I p_v L/\hbar}\tilde{\psi}(y)\). In other words, twisted boundary conditions should be applied with a twist \(\theta = p_vL/\hbar\).
GPE#
\[\begin{split}
\psi(x-vt) = e^{\I\phi(x-vt)}f(x-vt), \qquad f(x) = \sqrt{\rho(x)}\\
\psi' = \I\phi'\psi + e^{\I\phi}f', \qquad
\psi'' = \I\phi''\psi - (\phi')^2\psi +2\I\phi'\psi'
+ e^{\I\phi}\I\phi'f' + e^{\I\phi}f'',\\
\I\dot{\psi} = -v\I\psi' = v\phi'\psi - v\I e^{\I\phi}f'
\end{split}\]
\[\begin{split}
2v\phi'\psi - 2v\I e^{\I\phi}f' =
\I\phi''\psi - (\phi')^2\psi +2\I\phi'\psi' + e^{\I\phi}\I\phi'f' + e^{\I\phi}f'' + 2e^{\I\phi}f^3\\
[3\phi'+2v]f' = -\phi''f\\
[3\phi' + 2v]\phi'f = f'' + 2f^3\\
-\left(\frac{(\phi')^2}{2}\right)'
= \frac{f'f''}{f^2} + (f^2)'
= \frac{f'f''}{f^2} + \rho'\\
\end{split}\]
\[
f^2 = \rho, \qquad
2ff' = \rho', \qquad
2f'f' + 2ff'' = \rho'', \qquad
\]
Here we consider the problem of finding traveling wave solutions in a BEC. From a numerical perspective, it is highly beneficial if the problem can be stated in terms of a well-defined minimization problem. To start, we consider the available analytic solution for the conventional GPE. These solutions are presented in [El:2016]:
\[\begin{split}
\DeclareMathOperator{\sn}{sn}
\I\dot{\psi} = \frac{-\psi''}{2} + \abs{\psi}^2\psi\\
\rho = \abs{\psi}^2 = \frac{(r_4-r_3-r_2+r_1)^2}{4}
+ (r_4-r_3)(r_2-r_1)\sn^2(\sqrt{(r_4-r_2)(r_3-r_1)}\xi; m)\\
C = \frac{(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)(r_1-r_2-r_3+r_4)}{8}\\
u = v - \frac{C}{\rho}, \qquad
\xi = x-vt - \xi_0, \qquad
v = \frac{1}{2}\sum r_i, \\
a = (r_4-r_3)(r_2-r_1), \qquad
r_1 \leq r_2 \leq r_3 \leq r_4, \qquad
m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)},\\
L = \frac{1}{2}\oint \frac{\d\lambda}
{\sqrt{(\lambda - r_1)(\lambda - r_2)(\lambda - r_3)(\lambda - r_4)}}
= \frac{2K(m)}{\sqrt{(r_4-r_2)(r_3-r_1)}}.
\end{split}\]
The physical interpretations are:
\(a\): Amplitude
\(v\): Phase velocity
\(u\):
\[\begin{split}
\DeclareMathOperator{\sn}{sn}
\I\dot{\psi} = \frac{-\psi''}{2} + \abs{\psi}^2\psi\\
\rho = \abs{\psi}^2 = \alpha + a\sn^2(z; m)\\
\alpha = \frac{(r_4-r_3-r_2+r_1)^2}{4}, \\
D = \sqrt{(r_4-r_2)(r_3-r_1)}, \qquad
z = D\xi\\
C = \frac{(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)(r_1-r_2-r_3+r_4)}{8}\\
u = v - \frac{C}{\rho}, \qquad
\xi = x-vt - \xi_0, \qquad
v = \frac{1}{2}\sum r_i, \\
a = (r_4-r_3)(r_2-r_1), \qquad
r_1 \leq r_2 \leq r_3 \leq r_4, \qquad
m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)},\\
L = \frac{1}{2}\oint \frac{\d\lambda}
{\sqrt{(\lambda - r_1)(\lambda - r_2)(\lambda - r_3)(\lambda - r_4)}}
= \frac{2K(m)}{\sqrt{(r_4-r_2)(r_3-r_1)}}.
\end{split}\]
\[\begin{split}
\DeclareMathOperator{\cn}{cn}
\DeclareMathOperator{\dn}{dn}
\sn'(z) = \cn(z)\dn(z), \qquad
\cn'(z) = -\sn(z)\dn(z), \qquad
\dn'(z) = -k^2\sn(z)\cn(z)\\
\sn''(z) = -\sn(z)[\dn^2(z) +k^2\cn^2(z)], \qquad
\cn'(z) = -\sn(z)\dn(z), \qquad
\dn'(z) = -k^2\sn(z)\cn(z)\\
\cn^2 + \sn^2 = \dn^2 + k^2\sn^2 = 1.
\end{split}\]
A special limit is when \(m=1\):
\[
\sn(z;m=1) = \tanh(z), \qquad
\cn(z;m=1) = \dn(z;m=1) = \sech(z).
\]
\[\begin{split}
\psi = \sqrt{1+\alpha\sn^2(z)}\\
\psi' = \frac{\alpha\sn(z)\sn'(z)}{\psi}\\
\psi'' = \frac{\alpha}{\psi}\left(
\sn'(z)\sn'(z) + \sn(z)\sn''(z)
-\frac{\alpha\sn^2(z)[\sn'(z)]^2}{\psi^2}\right)\\
\end{split}\]
\[\begin{split}
\psi' = \frac{\alpha\sn(z)\cn(z)\dn(z)}{\psi}\\
\psi'' = \frac{\alpha}{\psi}\left(
\cn^2(z)\dn^2(z) - \sn^2(z)\dn^2(z) - k^2\sn^2(z)\cn^2(z)
-\frac{\alpha\sn^2(z)\cn^2(z)\dn^2(z)}{\psi^2}\right)\\
\end{split}\]
First we test these. To match units we set \(\hbar = m = g = 1\).
K(0.9)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[2], line 1
----> 1 K(0.9)
NameError: name 'K' is not defined
from gpe.imports import *
from scipy.special import ellipj, ellipk
import gpe.bec
def sn(u, m):
return ellipj(u, m)[0]
def K(m):
return ellipk(m)
class State(gpe.bec.State):
def __init__(self, Nx=32, rs=[0.1, 0.2, 0.3, 0.4], xi0=0):
r1, r2, r3, r4 = rs = sorted(rs)
k = np.sqrt((r4 - r2) * (r3 - r1))
m = (r2 - r1) * (r4 - r3) / (r4 - r2) / (r3 - r1)
v = sum(rs) / 2.0
L = 2.0 * K(m) / k
gpe.bec.State.__init__(self, Nxyz=(Nx,), Lxyz=(L,), m=1.0, g=1.0, hbar=1.0)
x = self.xyz[0]
self.T = abs(L / v)
t = 0
xi = x - v * t - xi0
a = (r4 - r3) * (r2 - r1)
# a = m**2*k**2
# rho_0 = (r4-r3-r2+r1)**2/4.0
rho_0 = (r4 - r3 + r2 - r1) ** 2 / 4.0
rho = rho_0 + a * sn(k * xi, m) ** 2
self.a_ = self.a = a
self.m_ = m
self.k_ = self.k = k
self.rho_0_ = self.rho_0 = rho_0
k_ = self.m * v / self.hbar
self[...] = np.exp(1j * k_ * x) * np.sqrt(rho)
def get_Vext(self):
return 0.0
s = State(Nx=256, rs=[-2.0, -1.0, 1.0, 2.0])
[I 03:38:46 numexpr.utils] NumExpr defaulting to 2 threads.
[I 03:38:46 root] Patching zope.interface.document.asReStructuredText to format code
/home/docs/checkouts/readthedocs.org/user_builds/gpe/checkouts/latest/src/gpe/utils.py:2755: PerformanceWarning: Default State.get_Vext_GPU() delegates to non-GPU version.
This is a potential performance issue, even if not using the GPU. To ensure optimal
performance, you should make sure that you overload get_Vext_GPU() instead.
warnings.warn(
/tmp/ipykernel_4984/1645578402.py:23: RuntimeWarning: divide by zero encountered in scalar divide
self.T = abs(L / v)
n = s.get_density()
x = s.xyz[0]
x_ = x[1:-1]
psi_ = s[1:-1]
ddpsi_ = np.diff(np.diff(s[...])) / np.diff(x)[1:] ** 2
plt.plot(x_, ddpsi_ / psi_ / 2 + abs(psi_) ** 2)
/home/docs/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/cbook.py:1719: ComplexWarning: Casting complex values to real discards the imaginary part
return math.isfinite(val)
/home/docs/checkouts/readthedocs.org/user_builds/gpe/conda/latest/lib/python3.14/site-packages/matplotlib/cbook.py:1355: ComplexWarning: Casting complex values to real discards the imaginary part
return np.asarray(x, float)
[<matplotlib.lines.Line2D at 0x76bb13d0d940>]
Check the soliton limit (2.122)
s = State(Nx=128, rs=[-4.000001, 1.0, 1.000001, 2.0])
x = s.xyz[0]
s.plot()
rho_min = abs(s[...]).min() ** 2
rho_max = abs(s[...]).max() ** 2
plt.axhspan(rho_min, rho_min + s.a, fc="y", alpha=0.5)
# plt.plot(x, rho_max-s.a/np.cosh(np.sqrt(s.a)*x)**2, '+:')
<matplotlib.patches.Rectangle at 0x76bb13c020d0>
rho_0 = s.rho_0
s[...] = np.sqrt(rho_0) * np.tanh(np.sqrt(rho_0) * x)
s.plot()
[<Axes: >]
s.cooling_phase = 1
e = EvolverABM(s, dt=0.2 * s.t_scale)
s = State(Nx=256, rs=[-2.0, -1.0, 1.0, 2.0])
with NoInterrupt(ignore=True) as interrupted:
while e.t < e.y.T and not interrupted:
e.evolve(1000)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
%debug
Solitons#
The GPE admits grey solitons with infinite period. These solitons arise from the previous solutions in the limit where \(m\rightarrow 1\) so that \(\sn\rightarrow\tanh\) and \(\cn\rightarrow\dn\rightarrow \sech\).
The soliton solution is:
\[\begin{split}
\psi(x, t) = \sqrt{\bar{\rho}}e^{-\I\mu t/\hbar}\left[
\I\frac{v}{c}
+ \frac{u}{c}\tanh\left(\frac{x-vt}{l}\right)\right], \qquad
\rho = \bar{\rho}\left[\frac{v^2}{c^2}
+ \frac{u^2}{c^2}\tanh^2\left(\frac{x-vt}{l}\right)\right]\\
u^2+v^2 = c^2, \qquad \mu=g\bar{\rho} = mc^2,\qquad l=\hbar/m u.
\end{split}\]
This can be expressed as
\[
\psi(x, t) = \sqrt{\rho}e^{-\I\mu t/\hbar + \I\phi(x)}, \qquad
\tan\phi(x) = \frac{v}{u\tanh\left(\frac{x-vt}{l}\right)}, \qquad
\phi_{\mathrm{twist}} = \phi(\infty) - \phi(-\infty) = 2\tan^{-1}\frac{v}{u}
\]
\[
\rho(0) = \bar{\rho}\frac{v^2}{c^2}
\]
or, in the same units \(m=\hbar = g = 1\) as El and Hoefer:
\[\begin{split}
\rho = \bar{\rho}\left[\frac{v^2}{c^2} + \left(1-\frac{v^2}{c^2}\right)\tanh^2[u(x-vt]\right]\\
u= \sqrt{c^2-v^2}, \qquad \mu=\bar{\rho} = c^2
\end{split}\]
To compare with (2.117):
\[
\rho = \frac{(r_4-r_3-r_2+r_1)^2}{4} +(r_4-r_3)(r_2-r_1)\sn^2[\sqrt{(r_4-r_2)(r_3-r_1)}(x-Vt), m]
\]
we have \(m=1\), \(\sn = \tanh\) and
\[\begin{split}
\frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)} = m = 1,\\
\frac{(r_4-r_3-r_2+r_1)^2}{4} = \bar{\rho}\frac{v^2}{c^2} = v^2\\
a = (r_4-r_3)(r_2-r_1) = \bar{\rho}\left(1-\frac{v^2}{c^2}\right) = \bar{\rho} - v^2\qquad
\sqrt{(r_4-r_2)(r_3-r_1)} = u, \\
V = \frac{r_1+r_2+r_3+r_4}{2}
\end{split}\]
Everything here is reasonable except the definition of the phase velocity \(V\) which differs from the soliton velocity \(v\).
The background density is
\[
\bar{\rho} = \frac{(r_4-r_3-r_2+r_1)^2}{4} + (r_4-r_3)(r_2-r_1)
= \frac{(r_4-r_3+r_2-r_1)^2}{4}.
\]
To compare with (2.122):
, the stationary solution must have \(\bar{\rho} = a_s\) or:
\[
(r_4-r_1)^2 = 4(r_4-r_2)(r_2-r_1)
\]
If \(r_1 = -r_4-2r_2\), then
\[
4(r_4+r_2)^2 = 4(r_4-r_2)(r_4+3r_2)
\]
Issues#
I was having some issues finding solutions so I considered stationary solutions:
\[
\psi = \sqrt{\rho_0 + a\sn^2(kx,m)}.
\]
Symbolically solving for solutions to the GPE (in Maple) gives the following solutions:
(28)#\[\begin{align}
a&=k^2m^2, &&& \mu &= \rho_0 - \frac{k^4m^2}{2\rho_0}, \\
a&=0, &&& \mu &= \rho_0, & \psi &= \sqrt{\rho_0}.\tag{Constant}\\
a&=\rho_0 m^2 , & k&=\sqrt{\rho_0}, & \mu &= \rho_0- \frac{\rho_0m^2}{2}, & \psi &= \sqrt{\rho_0}\sqrt{1+m^2\sn^2(\sqrt{\rho_0}x,m)}.\tag{Constant}\\
a&=\rho_0 , & k&=\frac{\sqrt{\rho_0}}{m}, & \mu &= \rho_0 - \frac{\rho_0}{2m^2}, & \psi &= \sqrt{\rho_0}\sqrt{1+\sn^2\left(\frac{\sqrt{\rho_0}}{m}x,m\right)}.\tag{Constant}\\
\end{align}\]
(assuming \(a,k\neq 0\) and real \(k\)) gives:
\[
m^2 = \frac{a}{k^2}, \qquad
\mu = \rho_0 - \frac{ak^2}{2\rho_0}, \qquad
\rho = \rho_0[1 - \sn^2(kx,m)].
\]
restart;
X_:=JacobiSN(k*x,sqrt(m2));
psi:=sqrt(rho[0] + a*X_^2);
Rho:=psi^2:
res:=collect(numer(simplify(-diff(psi, x, x)/2/psi + Rho - mu)), X_):
mu:=expand(solve(coeff(res, X_, 0), mu));
m2:=solve(coeff(res, X_^6), m2);
solve([coeffs(res, X_)], [a,k^2]);
The dark soliton solution has \(m=1\), \(a=k^2\)
s.m
\[\begin{split}
a = (r_4-r_3)(r_2-r_1)\\
a = mk^2\\
k^2 = (r_4-r_2)(r_3-r_1)\\
\end{split}\]
s = State(rs=[-4.1, 1.0, 1.1, 2.0])
x = s.xyz[0]
m = s.m_
a = s.a_
k = s.k_
rho_0 = s.rho_0_
print (m, a, m ** 2 * k ** 2, rho_0, k)
k = np.sqrt(rho_0)
plt.plot(x, 1 + sn(k * x, m) ** 2)
# s[...] = np.sqrt(rho_0*(1-sn(x,m)))
\[
m^2 = \frac{a\alpha^2}{D^2}, \qquad
(3m^2 - a) = 6\frac{a\alpha^2}{D^2}, \qquad
-2a(m^2 + 1) = 6\frac{a\alpha^2}{D^2}, \qquad
\mu = \frac{-D^2a}{2} + \alpha^2
\]
\[\begin{split}
-2a(m^2 + 1) = (3m^2 - a) = 6m^2\\
a = -3m^2\\
m^2 + 1 = 1
\end{split}\]
To do this, we first consider boosting to a moving frame with velocity \(v\) so that in this frame, the traveling wave solution is stationary and periodic. We second Here we are looking for solutions that satisfy:
\[
\psi(x + L) = e^{\I m v L / \hbar}\psi(x)
\]
Minimization#
Here we consider the following hypothesis:
The traveling waves can be found as minimum energy solutions in a box of period \(L\) with twisted boundary conditions holding both the total particle number fixed and the value of the wavefunction at one point.
We start with the soliton solution:
\[\begin{split}
\rho = \bar{\rho}\left[\frac{v^2}{c^2} + \left(1-\frac{v^2}{c^2}\right)\tanh^2[u(x-vt]\right]\\
u= \sqrt{c^2-v^2}, \qquad \mu=\bar{\rho} = c^2.
\end{split}\]
At the core, the density is:
\[
\bar{\rho}\frac{v^2}{c^2}
\]
from gpe.imports import *
from scipy.special import ellipj, ellipk
import gpe.bec, gpe.minimize
reload(gpe.bec)
class State(gpe.bec.StateTwist):
def __init__(
self,
Nx=32,
L=10.0,
mu=1.0,
psi_0=1.0,
ind=None,
v=0.0,
twist=None,
m=1.0,
hbar=1.0,
g=1.0,
):
if ind is None:
ind = Nx // 2
self.ind = ind
self.psi_0 = psi_0
self.p_v = p_v = m * v
if twist is None:
twist = np.pi + p_v * L / hbar
gpe.bec.StateTwist.__init__(
self, Nxyz=(Nx,), Lxyz=(L,), m=m, hbar=hbar, g=g, mu=mu, twist=(twist,)
)
self[self.ind] = psi_0
def exact_psi(self):
"""Return the analytic soliton"""
v = abs(self.psi_0)
u = np.sqrt(self.mu - v ** 2)
return 1j * v + u * np.tanh(u * self.xyz[0])
def get_Vext(self):
return -(self.mu + self.p_v ** 2 / 2.0 / self.m)
def _compute_dy_dt(self, dy, subtract_mu=False):
dy = gpe.bec.StateTwist.compute_dy_dt(self, dy=dy, subtract_mu=subtract_mu)
# dy[self.ind] = 0
return dy
class MinimizeState(gpe.minimize.MinimizeState):
def __init__(self, state, **kw):
self.psi_0 = state.psi_0
self.ind = state.Nxyz[0] // 2
gpe.minimize.MinimizeState.__init__(self, state, **kw)
def pack(self, psi):
psi[self.ind] = self.psi_0
fact = 1 if self.real else 2
x = gpe.minimize.MinimizeState.pack(self, psi)
return np.concatenate([x[: self.ind * fact], x[(self.ind + 1) * fact :]])
def unpack(self, x, state=None):
fact = 1 if self.real else 2
x = np.concatenate([x[: self.ind * fact], [0, 0], x[self.ind * fact :]])
state = gpe.minimize.MinimizeState.unpack(self, x, state=state)
assert np.allclose(state[self.ind], 0)
state[self.ind] = self.psi_0
# plt.clf()
# state.plot()
# plt.twinx()
# plt.plot(state.xyz[0], np.angle(state[...]/state.twist_phase), 'r--')
# display(plt.gcf())
# clear_output(wait=True)
return state
s = State(Nx=128 * 8, L=120.0, mu=1.0, psi_0=0.5, v=v)
psi_s = s.exact_psi()
v = (
(np.angle(psi_s[-1]) - np.angle(psi_s[0]) + np.pi)
* s.hbar
/ s.m
/ (s.Lxyz[0] - s.Lxyz[0] / s.Nxyz[0])
)
s = State(Nx=128 * 8, L=120.0, mu=1.0, psi_0=0.5, v=v)
s[...] = s.exact_psi()
print (v)
plt.plot(s[...] / s.twist_phase)
v
s = s1
s.cooling_phase = 1.0
e = EvolverABM(s, dt=0.5 * s.t_scale)
with NoInterrupt(ignore=True) as interrupted:
while not interrupted:
e.evolve(100)
plt.clf()
e.y.plot()
display(plt.gcf())
clear_output(wait=True)
s = State(Nx=128 * 4, L=20.0, mu=1.0, psi_0=0.1, v=0.1)
s[...] = 1.0
m = MinimizeState(s, fix_N=False)
s1 = m.minimize(use_scipy=True)
plt.plot(s1.xyz[0], s1.get_density() - abs(s1.exact_psi()) ** 2)
# s1.plot()
abs(s1.get_density() - abs(s1.exact_psi()) ** 2).max()
plt.plot(s1.xyz[0], np.angle(s1[...]))
plt.plot(s1.xyz[0], np.angle(s1.exact_psi()))
vs = np.linspace(0, 0.4, 10)
errs = []
for v in vs:
s = State(Nx=128 * 4, L=20.0, mu=1.0, psi_0=0.2, v=v)
m = MinimizeState(s, fix_N=False)
s1 = m.minimize(use_scipy=True)
errs.append(abs(s1.get_density() - abs(s1.exact_psi()) ** 2).max())
plt.plot(vs, errs, "-+")
s = gpe.bec.State(Nxyz=(64,), Lxyz=(10.0,))
s[...] = 1.0
m = gpe.minimize.MinimizeState(
fix_N=False,
)
s1 = m.minimize(use_scipy=True)
s1.plot()
Es = []
m = MinimizeState(s)
m.check()
s1 = m.minimize(use_scipy=True, fix_N=False, callback=callback)
# plt.plot(s1.xyz[0], np.log10(abs(abs(s1[...])**2-s.mu)))
s1.plot()
plt.plot(s1.xyz[0], abs(s1.exact_psi()) ** 2)
plt.plot(Es)
s1.get_density().max(), 1 + (s1.k_B[0] ** 2) / 2
-s1.mu, s1.get_Vext()
Es = []
twists = np.linspace(0, 2 * np.pi, 10)
for twist in twists:
s = State(Nx=128, L=50.0, mu=1.0, psi_0=0.05, twist=twist)
m = MinimizeState(s)
s1 = m.minimize(use_scipy=True, fix_N=False)
Es.append(s1.get_energy())
print twist, Es[-1]
s1.plot()
plt.plot(twists, Es)
m = MinimizeState(s)
s1 = m.minimize(use_scipy=True)
plt.clf()
s1.plot()
s1.get_energy()