teamatomecs / atomecs Goto Github PK
View Code? Open in Web Editor NEWCold atom simulation code
License: GNU General Public License v3.0
Cold atom simulation code
License: GNU General Public License v3.0
I just realized that the current implementation of the ramps does not allow for ramping two components of the same entity at the same time. This is difficult for me since I would like to replace my dipole::transition_switcher::RampMOTBeamsSystem
with them to ramp detuning and power of MOT beams. I believe there is no simple way to allow for two Ramps
to be associated with an entity, right?
One solution would be to make the wavelength an attribute of GaussianBeam
, so only one component needs to be ramped or we need a different approach to ramping in general. Or I leave it roughly like it is currently. (Custom systems to ramp them)
Also, I find it a little unhandy that I have to re-initiate the entire Component for every keyframe, even if I only want to change one attribute of it. This happens, if you, for example, want to ramp the power of a GaussianBeam
component. Just takes a lot of code, I think.
Any ideas?
Would anyone object if we added the num = "0.4.0"
to the list of project dependencies? It enables the use of complex numbers with:
num::complex
See: https://docs.rs/num-complex/0.4.0/num_complex/index.html
Makes much more sense, this is already half done for the Matlab scripts (some of which are being used to analyse examples for the accompanying paper).
Long-ranged forces are important for 3D MOTs with lots of atoms, as the repulsive forces reduce the phase space density.
Some thoughts:
MacroAtom
approach, where relevant quantities like repulsive forces are scaled so that a small number of atoms can simulate the PSD of an (intractable) larger number of atoms.Possible implementation:
leaf.force
+ force calculated pair-wise between other atoms in the leaf.Building tree will be O(L) O(N).
Following a good discussion with Tiffany:
Especially for dipole trapping, it can actually become problematic to pretend the Atoms have only one transition that is relevant. For example, what if you have a dipole beam for Strontium, close to the red transition? You'd have to respect the blue because it has a huge linewidth but the red one also, because you are close.
A far more realistic treatment would be to use polarizabilities instead and also E-fields instead of intensity gradients. I'd propose this as an enhancement for the foreseeable future. What do you think?
The GaussianReferenceFrame
already exists (and has an "ellipticity" attribute) in the new-dipole-force branch but elliptical Gaussian beams are yet to be implemented. Could be useful not only for dipole traps.
Hi all,
sorry for being quiet - I have been desperately trying to find the cause for a heating effect which seems to be still in the code but is most likely a dipole trap specific feature. I apologize in advance, this is super nasty and I'd be glad for new ideas how this could happen.
First observation: If you simply put atoms in a dipole trap, but also have some CoolingLight Entities in the game, however, very far detuned, so they shouldn't do anything - you'll notice that there appears to be a preferred direction for their escape - the x-direction. It can best be observed in the (branch: new-dipole-force) cross_beam_dipole_trap.rs example - where there should be no difference between x and y since we have two identical, horizontal beams in exactly these two directions (and the MOT beams are on the diagonals). However, if you plot the absolute coordinate values of the ensemble (or just have a very sharp eye when plotting the data in 3D in Blender) you'll see that they are more likely to be further away from the centre in (plus or minus) x-direction.
In addition to that, there is still a good amount of numerical heating happening on top of that, this effect, however can be reduced by increasing the detuning of the MOT beams even further:
Before you ask: I checked in a pure symmetrical 3D MOT if there was a difference between x and y and I don't think it is present there. Of course, I also checked the way the dipole force is calculated for any symmetry breaking but couldn't find one. At that point, however, I am running a bit out of ideas where to look for it.
For comparison, here are the same diagrams as above just without any CoolingLight entities in the code:
How did I find this behaviour in the first place? Well, I was running MOT -> DipoleTrap transitions with three stages:
(see file: example/red_mod_xodt_transition.rs)
I'd be extremely grateful for any, even far-fetched, ideas what could cause this behaviour.
Hi,
@tiffanyharte and myself would like to know what other people are using AtomECS for. For planning future features it would be very helpful to know what kinds of experiments people are interested in simulating - so please, anyone, feel welcome to comment on this issue :)
I'll add information from the Oxford setups below and tag others who have been using it.
I was wondering if the far detuned dipole beams of a FORT / ODT lead to an A.C. Stark shift of the optical transitions used for a MOT. For the example of a 1064nm dipole trap that is briefly superimposed with a red Sr-MOT, I tried to calculate the light shift according to 7.7 Foot (see PDF attached). As you can see, the result I get is pretty much zero but I'm confused by the sentence "Normally light shifts are most important at large frequency detuning [...]" in 7.7 Foot.
What do you think - is my rough estimation of the magnitude of the light shift correct in principle or have I missed something crucial in my approach? Does someone have practical experience with MOT-->ODT transitions and knows if that effect becomes relevant in real life?
Why I am needing this: I want to simulate MOT --> ODT transitions with AtomECS and I want to know if an a.c. Stark shift is something we'd have to implement for this to work correctly.
Stark-shift_ODT.pdf
The QuadrupoleField3D
and QuadrupoleField2D
have their position associated with an external Position
component, while PrecalculatedMagneticFieldGrid
has its position defined internally. I was wondering if there is any reason for such inconsistency?
I guess this could be extended to laser beam as well because GaussianBeam
also has its "position" (intersection) defined internally.
The code seems to assume a Gaussian beam profile, I would like to simulate a MOT geometry using laser beams with a tophat profile or more generally be able to modify the beam profile to be an abritrary function. Adding a tophat module in lasers is not sufficient as the code continues to assume a gaussian beam in other instances.
We are unsure if the total scattering rate in photons_scattered.rs
should be calculated with the linewidth in Hertz or angular units. The same question has been raised for the two-level excited state population in twolevel.rs
where the Einstein A coefficient is used.
The corresponding section in C. Foot: 7.5.2
Starting to profile the latest changes to red-mot
to get a head start on optimisations.
Interestingly, it seems the vast majority of program time was being taken up by rayon implementations - specifically SwitchToThread
and SleepConditionVariableSRW
. There was also poor spin time.
This seems likely related to a known issue in rayon, where sleeping threads would burn a lot of cycles/power. For example, see
rayon-rs/rayon#642
rayon-rs/rayon#795
rayon-rs/rfcs#5
It looks like some of these issues were fixed in rayon 1.4.0
, although we are currently on rayon 1.3.0
. I'll try updating rayon (new version is 1.5.0
) and see if it improves things.
There are two problems though the converted power in the current implementation is correct.
Line 71 in f2a5804
Shouldn't the right hand side be multiplying e_radius
by sqrt(2)
?
Line 72 in f2a5804
Shouldn't the right hand side be like divide the product of peak_intensity
, pi
e^2_radius
by 2 instead of multiplying it by 2?
Given that
(Edmund's optics, RP Photonics)
Rearranging would give
According to the documentation
Line 37 in f2a5804
making the changes suggested in Q1 would give the same result because following the suggestion above would give
And the current implementation is
which gives the same result.
However, this leads to the problem of a small optimization where the square root and the multiply by two can be simply neglected.
1d MOT is currently broken - it's probably related to the fact that ScatteringFluctuationsOption::Off
can give NaNs
I have discovered that the length of the integration time-step has a very significant effect on the "numerical" heating in a dipole trap. I attached a picture, showing the simulation results of 3 (otherwise identical) simulations with:
as the delta - timestep. They started with 100 atoms and ran for a total real time of 0.2 seconds each.
As you can see from the picture, for ODT one should choose a smaller timestep than in a MOT setup. I just felt this is something to be aware of. I'll try to find the largest time-step with negligible heating over a couple of seconds. Any comments?
There are quite a few branches left over, it needs a tidy up.
We have to redocument the newer classes, as sometimes doc coverage has been lost due to changes (all in the name of progress, which is good :) )
Make sure the final code is consistent with cargo doc --open
.
I would like to use the uom
crate to define quantities in atomecs.
https://github.com/iliekturtles/uom
Advantages:
I had tried previously but was not able to because of the incompatibility of Vector3<>
from nalgebra
with uom
. There have been a few updates to the package since then and it is worth revisiting.
Hey Maurice,
We now have ApplyEmissionForceOption
which enables to use of random emission kicks required for recoil limit and part of the Doppler limit.
I say part of Doppler T, because even with this option disabled the simulation will still draw random numbers from the mean dists, and so there will be fluctuations in the absorption forces. That contributes part of the Doppler T.
Suggested course of action, we change ApplyEmissionForceOption
to something else (eg, EnableScatteringFluctuations
). When this resource is present, the systems which transfer 'mean numbers' to 'actual numbers' use Poisson dist; when it isnt, they just transfer the numbers exactly (so the execution is faster).
To simulate evaporation in optical/magnetic traps, AtomECS requires a way to calculate the effects of short-ranged (eg, s-wave) collisions.
Given that the collisions are short-ranged, all methods used should a spatial binning method for speed.
Beyond that, there are subtle variations in the implementation details, and we could support a number of approaches:
Monte-Carlo method: determine a number of collisions to occur within a box, eg based on density and mean velocity. Perform collisions between velocities of random atomic pairs.
'Pair-wise' Monte-Carlo method; Iterate over each pair within a box. Calculate the chance for a collision to occur, then draw a random number to determine if it did.
All models should also have the ability to specify some kind of 'macro atom', so that the behaviour of a large number of atoms can be represented by a numerically tractable ensemble. There's a few ways we could do this:
MacroAtom
components.I like the global scaling parameter because it is straight forward. Also, I can see how we might be able to simulate multiple orders of magnitude of atom loss by rescaling the parameter as atoms are lost. For instance, this could allow one to start with 1e4 atoms representing 1e10 atoms and end with 1e4 atoms representing 1e4 atoms.
bevy
has a really nice implementation of 'Plugins', where you can define individual parts of a program to add/remove. A similar plugin model could be used in AtomECS to add different models of laser cooling, different forces, collisions, etc. This would improve the usability and increase program flexibility (possibly also performance by only running required/declared modules).
Related to #3.
Add support for magnetic traps.
For the purposes of this issue, this refers to traps where the potential energy can be written in the form E = g_F \mu_B | B(r) |
, and where B is a real magnetic field. This includes static magnetic field traps, such as quadrupole traps, and traps where the magnetic field is slowly varying, e.g. Time-Orbiting Potential (TOP) traps.
This issue does not include support for RF-dressed traps (where the magnetic field is oscillating at a rate comparable to the Zeeman splitting).
Hi,
I've been using AtomECS a bit lately and I just had a question on how it is possible to handle linear polarization. In the doc it's mentionned that it should be possible to have a linear polarization by combining sigma_+/- beams, but I don't understand how that can be.
For example I'd like to consider the following situation: a magnetic field along e_z
and a cooling light with wavevector along e_x
linearly polarized along e_y
. This light can induce sigma_+/-
transition, but not pi.
However, when I look at the calculation for CalculateRateCoefficientsSystem
, the term scatter3
for the pi component will always be non zero whatever the configuration I set for the polarization of the cooling lights.
What would be the correct way to handle this situation?
Thanks
Sorry for spamming so many issues here - let me know when to stop. ^^
Can you debug AtomECS in VS Code? It seems it never stops at my breakpoints in systems... am I getting something wrong? Hast it something to do with the strong parallelism? It would be great to stop at points in the simulation have a look at the entries in all components, is that possible?
Thanks in advance!
There are currently a number of ways of configuring simulation options which should be standardised in format. For instance, we currently have EnableScatteringFluctuations
and ApplyEmissionForceOption
, which are dissimilar naming conventions and should be made consistent.
The options also currently favor speed over accuracy. I propose we change it so that accuracy is favored over speed; users disable functionality to improve simulation speed when required, but they only compromise physics when they choose to.
It looks like the scattering force abruptly cuts off at weaker scattering. Running plot_molasses.m
with more points illustrates this well:
You can see that for velocities above ~22 the force abruptly cuts to zero.
I expect this is related to the Poisson implementation when the mean is low, but will investigate more.
I would like to see a high level API for creating entities in the simulation, using something like factory structs with a fluent interface. e.g.:
BeamBuilder()
.in_direction(Vector3::x)
.with_gaussian_profile(waist, power)
.with_sigma_plus_polarisation()
.for_cooling(Strontium)
.with_detuning(-32.0)
.build(&mut world);
The advantages of this are:
Disadvantages:
atomecs::prelude::*
mod.Originally posted by @ElliotB256 in #29 (comment)
Performance is slightly better with Legion than specs (Amethyst has now dropped specs and is moving to Legion). Consider moving to Legion.
https://github.com/amethyst/legion
Hi @MauriceZeuner and @d-garrick ,
I added a WIP implementation of ramps using a generic system here: https://github.com/TeamAtomECS/AtomECS/blob/collisions/src/ramp.rs
I intend to define a macro, rampable
, so you can auto derive the required trait impl.
Euler integration is currently used - it's straightforward but has terrible integration error. I'll change it so that AtomECS uses Velocity-Verlet by default. VV has much smaller integration error. This is particularly important for future thermodynamic simulations of atoms in traps, as integration error leads to heating similar to parametric excitation.
I'll admit I thought we were already using the Velocity-Verlet integrator - I confused atomecs with the other ECS code (Unity/C#). My other code uses the follwing implementation:
/// <summary>
/// Update velocity according to velocity verlet.
/// </summary>
[UpdateAfter(typeof(ForceCalculationSystems))]
[UpdateInGroup(typeof(FixedUpdateGroup))]
public class UpdateVelocitySystem : JobComponentSystem
{
protected override JobHandle OnUpdate(JobHandle inputDependencies)
{
float DeltaTime = FixedUpdateGroup.FIXED_TIME_DELTA;
return Entities.ForEach(
(
ref Velocity velocity,
in Mass mass,
in Force force,
in PrevForce oldForce
) =>
velocity.Value = velocity.Value + 0.5f * (force.Value + oldForce.Value) * DeltaTime / mass.Value
).Schedule(inputDependencies);
}
}
/// <summary>
/// Update position according to velocity verlet.
/// </summary>
[UpdateBefore(typeof(ForceCalculationSystems))]
[UpdateInGroup(typeof(FixedUpdateGroup))]
public class UpdatePositionSystem : JobComponentSystem
{
protected override void OnCreateManager()
{
Enabled = false;
}
protected override JobHandle OnUpdate(JobHandle inputDependencies)
{
float DeltaTime = FixedUpdateGroup.FIXED_TIME_DELTA;
return Entities.ForEach(
(
ref Translation translation,
in Mass mass,
in PrevForce force,
in Velocity velocity
) =>
translation.Value = translation.Value + velocity.Value * DeltaTime + 0.5f * (force.Value / mass.Value) * DeltaTime * DeltaTime
).Schedule(inputDependencies);
}
}
I'll add this to AtomECS when I get time, it is a straightforward addition and worth it.
On most recent toolchain, running AtomECS gives a long error which ends in:
thread panicked while panicking. aborting.
As a temporary fix, use the stable toolchain from 2020-06-17. Type the following into the console:
rustup install stable-2020-07-16
rustup default stable-2020-07-16
There is a known issue related to specs
, and it is caused by one of its dependencies, hibitset
. This restricts the maximum number of entities to usize^4
(ie. 32^4
for 32-bit systems or 64^4
for 64-bit systems). I don't think there really is an easy solution other than writing a PR for hibitset
. However, I am just wondering if you happened to have some workaround in the occasions where there is a large number of atoms? Plus, is there some kind of plan to overcome this in the future?
I was trying to increase the number of cooling lasers to replicate sidebands.
Every time I add more than 16 lasers to the world, a thread panicks.
Any thoughts on what could be causing the problem ?
Hi here,
I am new to here, and I want to fix compiling and clippy first, if there is any other issues, I may take one.
Hey @MauriceZeuner ,
Could you check and fix this unit test? It's been broken since I fixed the scattering rates, so I'm confident it's the test is outdated rather than the implementation (for instance, the integration test for scattering rate is working fine and agrees with the analytic formulae).
Cheers!
What is the correct way to generate atoms from an oven with multiple transitions ?
I assumed that I could just add an AtomicTransition component to the oven creation, but it doesn't seem to work.
For example, replacing in 2d_plus_mot_from_oven the line
l156: .with(AtomicTransition::strontium())
either by
l156: .with(AtomicTransition::strontium())
l157: .with(AtomicTransition::strontium_red())
or
l157: .with(AtomicTransition::strontium_red())
l156: .with(AtomicTransition::strontium())
doesn't give the same simulation results in the two cases.
It seems that in the first case, adding the red transition overwrites the blue.
I actually just wanted to see if an error that I encountered on the dipole branch also happens to be in the master branch (good news: it's not) but I noticed that when I ran cargo test
, two tests fail, namely:
magnetic::tests::test_field_samplers_are_added
magnetic::tests::test_magnetics_systems
I am happy to look for a solution if no one else feels responsible or has time but I thought it would be nice to make you aware of these fails first.
Full terminal output is attached and also a screenshot of the relevant section.
magnetic_test_fail_full_output.txt
I think it would be great to have a very simple way to write not only Position and Velocity to a textfile and analyse it in Python/Matlab but maybe a general WriteComponentToFileSystem
that users can configure to write the values they are interested in.
For example: If I wanted to print the contents of the TotalPhotonsScattered
component for all atoms over time, I might want to plot that against the power
of all GaussianBeam
entities. Currently (please correct me if I'm wrong) I write my own writer-system for that every time I need it and there exists no nice and most importantly run-time optimized way of doing that. (See the super-slow xyz-writer, ahem, ahem... ).
But especially if one wants to monitor all sorts of interesting parameters over complicated runs (red Sr-MOT to dipole transition, for example), it would be so cool to get that data printed to a text file in one line of code.
What do you think?
Hi,
is it a known problem that Anti-Virus software (such as Kaspersky Internet Security for instance) think that the
AtomECS-master\target\rls\debug\build\crossbeam-epoch-9c14771e7fefb252\build-script-build.exe
is some sort of malware? It is a bit annoying since it always forces me to "desinfect" (aka delete) files from the project directory. Moreover, because the Anti-Virus software "keeps" the file occupied, cargo run --release
fails due to
error: failed to run custom build command for crossbeam-epoch v0.8.2
Caused by:
could not execute process D:...\AtomECS-master\target\release\build\crossbeam-epoch-576fdc92c8f92a03\build-script-build (never executed)Caused by:
Access is denied. (os error 5)
warning: build failed, waiting for other jobs to finish...
error: build failed
The software gives me the following information on what it thinks about the "malware":
https://threats.kaspersky.com/en/threat/Trojan.Win32.Agent/
Best wishes and thank you!
Maurice
Are you considering publishing this on crates.io? There have been some "bad actors" who claims the crate's name by just publishing an empty crate on crates.io, and you should probably consider claiming this name to avoid any potential issues in the future?
Is there a package or algorithm you could recommend for calculating numerical gradients of a function that takes the form:
f(pos: Vector3<f64>) -> f64
Strictly speaking, I need it to get the intensity gradient (E-field-squared-gradient) for cases of many-beam interference, where an analytical solution is hard to obtain (Although I will think about it some more... for some special cases (like linear polarized Gaussian beams only).
I've seen https://docs.rs/GSL/4.0.0/rgsl/numerical_differentiation/index.html but I don't know how to "install" it.
Obviously, I could just write a 2,3 or 5 -point algorithm myself, and implement it in maths
- but I was wondering if there existed an optimized algorithm for that already.
I need to update to the latest specs+rayon to fix the panic issue, so we can get back to nightly toolchain.
The specs API has some breaking changes, so I need to update all examples aswell.
Something Tiffany and I would be very interested to see some results of in the mid-future.
Hi!
I finally started properly fixing the first of the numerical heating problems. For this, I want to replace the faulty maths::random_direction
with rand_distr::UnitSphere
since Elliot recommended that
Hmm, I didn't realise we were defining our own function for it! I assumed we were using one from another crate. Can you replace ours with https://docs.rs/rand_distr/0.4.0/rand_distr/struct.UnitSphere.html ?
see #27
So, I was trying that, however, I encountered a problem where the compiler complains about the argument of UnitSphere.sample
in force.rs
use rand_distr::{Distribution, UnitSphere};
// explicit random walk implementation
for _i in 0..total {
let v: [f64; 3] = UnitSphere.sample(&mut rng);
force.force = force.force
+ force_one_kick * Vector3::new(v[0], v[1], v[2]);
}
saying:
error[E0277]: the trait bound `atom_sources::rand::prelude::ThreadRng: rand_core::RngCore` is not satisfied
--> src\laser_cooling\force.rs:174:77
|
174 | ... let v: [f64; 3] = UnitSphere.sample(&mut rng);
| ^^^^^^^^ the trait `rand_core::RngCore` is not implemented for `atom_sources::rand::prelude::ThreadRng`
|
= note: required because of the requirements on the impl of `rand::rng::Rng` for `atom_sources::rand::prelude::ThreadRng`
which makes no sense to me since I don't see any rand
being imported from atom_sources
nor is it existing there at all. I don't quite know where to start, to be honest, so I'd be glad about suggestions. ^^
Thank you!
@MauriceZeuner is working on a dipole force trap implementation
The current ApplyEmissionForceSystem
uses a normal distribution when the number of scattered photons N > 5. I don't think this respect recoil limit for reasons we have discussed beforehand.
The code runs faster when you approximate N kicks with a normal distribution, and the recoil limit isn't applicable for some systems. Let's add an optional flag to either always perform individual kicks, or use the distribution. This could be an enum in the ApplyEmissionForceOption
struct.
Just to let people know that I am working on extending the Gaussian beam description to correct treatment of (at least fully) polarized light with arbitrary ellipticity of the Jones vector.
If you wonder why: Well, I think we do already have 50% of the stuff we need to model Sisyphus-cooling and I suppose that this would definitely be a nice feature in the mid-future. Any thoughts?
The branch will be "alternative-dipole-force" where I also switch to a full complex electric field treatment for dipole traps.
At the moment, I imagine it working like this: Your basic component always is the LinearGaussianEBeam
and then you just add a (name provisional) JonesVector
component to the beam entity. Since you can separate the rotation of the electric field vector and the other gaussian beam characteristics, you can then in a SampleElectricFieldSystem
just (complex) multiply the jones vector (which is basically a rotation-operator) to the other terms of the gaussian E-field and you should get a (for example circularly) polarized gaussian beam.
For Sysiphus cooling I think, you should be able to get a standing wave in polarization (not intensity!) by just adding two counterpropagating beams with orthogonal polarization to the simulation. At positions where the counter-propagating beams have a phase difference of π/2 the polarization is circular, and where there is no phase difference, the polarization is linear. In the intermediate regions, there is a gradient elipticity of the superposed fields.
What, of course would also be required for this to work (will be a separate issue, probably) is a slightly more complicated atom-component which can respect the important M_J substrates of the atom.
Following the nasty heating error in the dipole trap, I'd propose we should think about how to generally clean component data (in each iteration step). Two things come to mind:
Force
. Pro: clean, Con: additional code that runs every stepI believe this topic becomes more and more important as we move towards simulation files that combine various side-arm stages.
What do you think?
Hey @d-garrick,
I noticed your recent merge incorrectly added some matlab code to the repository under the +util
folder. Can you please delete it? That code should be added to https://github.com/TeamAtomECS/matlab instead.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.