Improper casting of magnetic force in OBE force()

[SteveEckel]

For some reason, it can sometimes be left as complex. See the following report from the Google Group:

For the example online for real atoms in a MOT, the cell for plotting 87Rb result is in fact for 7Li. It'd be nice to correct it.

Now, I add the following segment to solve 87Rb MOT using obe after the trap_D2.generate_force_profile statement:

obe_D2 = pylcp.obe(laserBeams_D2, magField, hamiltonian_D2,
                     transform_into_re_im=True,
                     include_mag_forces=True)
obe_D2.generate_force_profile(
    [np.zeros(X.shape), np.zeros(X.shape), X],
    [np.zeros(V.shape), np.zeros(V.shape), V],
    name='MOT_1', deltat_tmax=2*np.pi*100, deltat_r=4/alpha,
    itermax=1000, progress_bar=True
)

I get an error at f+=f_mag in force() in obe.py: UFuncTypeError: Cannot cast ufunc 'add' output from dtype('complex128') to dtype('float64') with casting rule 'same_kind'

[SteveEckel]

I cannot reproduce this issue, and I do not want to force a fix without being able to reproduce the bug.

[SteveEckel]

Actually, I can reproduce this bug by exactly replicating exactly what was suggested on the Google Group. With debug enabled on the generate_force_profile statement, I get the following:

Finding equilbrium force at r=(0.00, 0.00, -5.00) v=(0.00, 0.00, -5.00) with deltat = 5.03, itermax = 1000, Npts = 5001, rel = 1.0e-05 and abs = 1.0e-09
0 [ 2.82544548e-04 -8.55377144e-05  1.88325669e-01] 0.03546664462912284
1 [-0.00030088 -0.00059065  0.1792337 ] 0.032125160028288906
2 [-0.00135071 -0.00421136 -0.04742452] 0.002268644868615372
3 [-0.02640673  0.00398102 -0.43183555] 0.18719510205354298
4 [-0.11755991 -0.05992344 -1.26343037] 1.6136674483604287
5 [-0.80805931 -0.11903923 -2.26899962] 5.815489439318571
6 [ 3.88685817 -0.08798319 -5.9180601 ] 50.13884279294571
7 [-15.22338883   4.36235643  15.58482576] 493.6685151862424
8 [-58.76806442  30.1547223   26.40886731] 5060.420944808086
9 [   9.32313569 -244.31696433   16.58299033] 60052.69548822518
10 [976.71597017  17.11726695 -20.12871358] 954672.2523142454

So clearly it is not really a bug, because it was working through multiple interactions, until it stopped working. It's really unclear what is going on, except that the force is growing and growing and growing. That instability is probably due to the fact that we are using an inconsistent coordinate system, or at least a coordinate system where the strength of the magnetic field is way too big.

[dsbarker]

I've now seen forces growing without bound in simulations of frequency-swept MOTs. For a small number of states, the issue could be fixed with transform_into_re_im=False, but forces for atoms with more complicated hyperfine structure still blew up.

I suspect that the issue is that the evolve methods use the RK45 integrator with default precision: rtol=1e-3, atol=1e-6. As with all numeric integrators, the RK45 integrator will not strictly conserve physical quantities and this issue appears to get worse with more complicated systems of differential equations (see this stack exchange post).

I've been able to solve this issue in some systems by setting more stringent limits on rtol and atol, which can be passed as keyword arguments to any of the pylcp functions that rely on solve_ivp. Switching to the DOP853 integrator should save some computation time after doing so (see the solve_ivp docs).

To really solve this issue, we may need to investigate using a symplectic integrator (see this stack exchange post). However, scipy does not offer any as methods for solve_ivp and I'm not certain that they can be applied to density matrix evolution in a straightforward way.

It looks like the only python package with symplectic integrators is diffeqpy, but using that would mean learning (and requiring) Julia, because diffeqpy is a python wrapper for DifferentialEquations.jl.

[dsbarker]

See also #71 for discussion of relative performance of DOP853 and RK45 integrators.