Menu ▾ ▴

Home

simonvg Pierre Dupré radbal

SIMBUCA

Formerly known as the Simonion project. SIMBUCA is a simulation package that stands for: Simulation of Ion Motion in a Penning trap with realistic BUffer gas collisions and Coulomb interaction using A Graphics Card. However it can and has been applied for other applications that require charged particle tracking under the influence of Electric and Magnetic fields with/without Coulomb interaction included.

Table of contents:

Installation and first usage

The installation is relatively easy and does not require much linux knowledge. The files can be found on http://sourceforge.net/projects/simbuca/

  1. Simbuca first usage
    • install Simbuca
    • select cpu/gpu and icpc/g++/mpic++ compiler in the Makefile
    • type make all in the bin directory. The executable simbuca will be constructed in the map simulations.
    • cd to the map simulations
    • configure the parameters of the simulation in a .sim script file
    • excecute the executable by typing ./simbuca myscript.sim (without mpi)
  2. Compiling Simbuca
    • cd bin
    • make all
  3. Working with the GPU or CPU
    • pico Makefile (or use another editor)
    • change the pu variable to cpu or gpu
    • type make clean to remove object files
    • remake all the object files and the executable by typing make all in the bin directory
  4. Choose your N-body algorithm for the Coulomb interaction calculation
    • pico Makefile (or use another editor)
    • change the NBODY ALGO variable to NBODY for the NVIDIA NBODY ALGORITHM
    • change the NBODY ALGO variable to CUNBODY for the CUNBODY ALGORITHM
    • type make clean to remove object files
    • remake all the object files and the executable by typing make all in the bin directory
  5. Compiling with icpc or g++ or mpic++
    • pico Makefile (or use another editor)
    • change the compiler variable to icpc or g++ or mpic++
    • remove the old object files by typing make clean in the bin directory
    • remake all the object files and the executable by typing make all in the bin directory
  6. Doing a trap simulation
    • All the information a user needs is provided in the .sim file for Penning trap simulation. Few advanced options must be hard coded (such as Image Charge calculation).
    • The functions a user can use, p.e. for trap excitations, are described in IonFly.h .
    • Or a user can use a input file ”.sim” with all the parameters of the simulations (see the User Interface section for more details).

FAQ

  1. How to get started?
    A Linux preferentially. Windows should work aswell but was not tried with a graphics card and linking can be different in this case. If you want to enable Coulomb interaction calculation on a GPU, a Linux 64- bit is needed and a CUDA ready GPU. In this case also pay attention that there is enough cooling for your GPU. A fast CPU is also recommended since the time stepper is calculated on the CPU but do not waste too much money on this. concerning RAM & HD there are no constraints.
  2. How do I start a simulation?
    See section Simbuca first usage.
  3. How to use a graphics processor?
    This is definately the hardest part. Nevertheless what one has to do is installing the NVIDIA driver then the CUDA Toolkit and finally CUDA SDK.
    Use google to see what extra classes and libraries are needed in order to install everything correctly.. Edit the Makefile (in the bin directory) so that pu=gpu, remove old object files by typing make clean and make them again for the GPU by typing makeall
  4. Why and how to use the NBODY algorithm instead of CUNBODY?
    The NBODY algorithm from NVIDIA SDK is an alternative to the CUN- BODY algorithm. It’s faster since it calculates the Coulomb interaction on the GPU meanwhile the CPU computes the external forces exerted on ions. However NBODY is more restrictive since the number of ions has to be a power a 2.
    • Test the NBODY algorithm with the test nbody.cpp code on tools directory.
    • Set the NBODY ALGO to NBODY
  5. How to use MPI (Message Passing Interface) for the CPU parallelization of Simbuca?
    CPU parallelization can be done using MPI, on a cluster or on a multi-core CPU.
    • install a MPI package such as MPICH2
    • configure MPI (deamons) by typing: mpd −−ncpus=Nn &. Nn is the number of CPU or cores used. For example Nn =4 with a quad- core CPU. Then type: mpdboot &.
    • Use the compiler option mpic++ in the Makefile
    • Set the good path of the MPI include directory named MPICHDIR in Makefile
    • Launch Simbuca by typing: mpirun -l -n Nn ./MyFirstSimulation
      **The number of ions has to be divisible by the Nn. **
  6. How does the Makefile work?
    Simbuca consists of different files that are first (as always with a c++ project) compiled separatly into object (.o) files and after this linked to- gether by the linker. Normally this would mean that the end user hast to type the compile and link command in the linux terminal. However since this is rather cumbersome and not easy if one wants to change between running the simulation on a GPU or on a CPU a Makefile was made. The makefile hides compilation and linking for the user. Some commands for this makefile are:
    • make clean removes the object files and the executable(do this if you go from GPU to CPU or vica versa)
    • make all makes all the object files and the executable
    • make export makes a backup of the files and saves this in bin/backup.
    • make parameters prints out the parameters that can be set in the Makefile: Processing Unit (GPU or CPU), name executable, compiler (g++ or icpc or mpic++)

  7. What is subversion?
    Subversion (svn) allows synchronisation with the latest changes. After in- stalling svn on a linux machine one can get the latest version of Simbuca by typing: svn checkout svn://svn.code.sf.net/p/simbuca/code/ simbuca. A directory "simbuca", will be made where all the files are in. To checkout if some files are updated and possibly download them. Type svn up in the directory Simbuca to update to the latest file changes. This way the version you‘re using is always up date and the least bugfree. More information about subversion can be found by typing svn help or online http://wiki.apertium.org/wiki/Using_SVN.

  8. What is the following error: "Problem to open the file of particle 1020"?
    This is related to the linux permissions. On most linux distros (for example Ubuntu) you can open by default only 1024 files simultaneously, when simulating > 1020 particles and when selecting the option to print the information of each particle in a separate file this error will hence pop up.
    The administrator of the machine can change this number with the "ulimit -n x" command were x is the max number of files to be open simultaneously. Alternatively it can be coded in the "
    /etc/security/limits.conf" file.
    You can check with the command "ulimit -a" how many files you can open.

Simulation questions

  1. How do I include my Penning trap dimensions?
    Trap specific parameters are in trapparameters.h. Here one can change the trapdepth and the magnetic field. Note that this will be changed in a future version such that Simbuca doesn`t need to be recompiled if one changes the magnetic field.
  2. Can I use the program also for RFQ’s or Paul traps?
    YES. In fact it is rather easy. This only means that the user has to change the equations of motion in force gpu.cpp and force cpu.cpp
  3. How is the output of the program formatted?
    It can be chosen to print each particle in one file or all particles together in one file. Information is printed out after a certain timestep in the simu- lation (this can be chosen aswell). Information printed out is: index(ID), mass, x,y,z (in mm), vx,vy,vz (in m/s), Energy (eV), magnetron radius, cylotron radius, total radius (mm).
  4. Are the phases of ions used in the simulation?
    NO. However if there is a need for it can be inserted in the program (this means just a simple extention is the Ion or Particle class). Here it was at first not implemented since in large ionclouds the phases are cancelling out.

User Interface

Description

The parameters of the simulation can be directly stored in a ”.sim” input file read by Simbuca. This input file must be the first argument of the executable and its extension must be ”.sim” :

no MPI : ./simbuca foo.som
MPI : mpirun -n 4 ./simbuca foo.sim

.sim file parameters

The sim file contains 4 kinds of information:

  • the parameters of the cloud: number, isotopes, dimensions, energy...
  • the parameters of the ODE: order, time step, Coulomb interaction
  • the configuration of the trap
  • the different operations: No excitation, dipole/quadrupole/octupole excitation, rotating wall, transfer

The general layout of a .sim file stores the following information in this order:

1: Particle number N
2: initial cloud parameters
3: buffergas parameters
4: integrator parameters
5: Coulomb interaction parameters
6: output parameters
7: Trap configuration parameters
8: rf-excitations parameters

  • 1: Particle number N: an integer

  • 2: initial cloud parameters: The cloud can be defined with one of the following commands

    • CREATECLOUD: creates a speroidal cloud of one or more specieswith a Maxwell Boltzmann energy distribution
      For example: 

      CREATECLOUD 100 3
CLOUDPARTS 7Li 5. 6Li 80. 39K 15. 
CLOUDCOORD 0.002 0.001 0.007 0.0001 0. 0.0005
TEMP 1.5
      • first CREATECLOUD 100 3 means to create a cloud of 100 particles with 3 species.
      • the isotopes used are specified with CLOUDPARTS, that is 7Li (5%) 6Li (80%) and 39K (15%). The total amount has to be 100 %. The masses are loaded from the AME file: mass.mas12 and the Particle data group (PDG) file.
      • the cloud has a spheroidal shape with semi axis lx=0.002, ly=0.001 and lz=0.007 m, the center of the cloud is at the position (0.0001, 0., 0.0005). This is configured with CLOUDCOORD.
      • A Maxwell-Boltzman energy distribution is used with a mean energy of 1.5 eV. This is configured with TEMP.

    • CREATEPARTICLES: each line defining particles along with position and velocity
      For example:

      CREATEPARTICLES 3
PARTICLES proton 0.0 0.01 -0.15 0. 0. 0. 1H 0.0 -0.01 -0.15 0. 0. 0. 35Ar6+ 0.0 -0.01 -0.15 500. 0. 0.
      • 3 particles are defined, their number has to be put after CREATEPARTICLES.
      • With PARTICLES key a proton, a hydrogen atom (which is the same in principle), and a 35Ar atom with a charge of 6+ are defined. Their names and parameters must be in one line, following each other. After the name the positions x,y,z are the following first 3 values (m), and finally the velocities vx,vy,vz are the next 3 values (m/s), for each particle. Note: The masses are loaded from the AME file: mass.mas12 and the Particle data group.

    • IMPORTDATA: Import data from a previous simulation. Follow IMPORTDATA followed by the prefix of the importfile

  • 3: buffer gas parameters, BUFFER : 1 integer + 1 double

    • First parameter is 1/0 to switch on/off the usage of buffergas. The second parameters is the buffergas pressure in mbar.
      For example: 

      BUFFER 1 1e-4
      • With BUFFER key a buffergas pressure of 1e-4 mbar is defined

  • 4: integrator parameters, ODE: 2 integers + 1 double

    • the first integer is the order of the ODE: 0 for the Velocity Verlet integrator, 1 for Gear method, 4 for Runge Kutta 4th order, 5 for Dormand Price 5th order.
    • the second integer indicates if the adaptive stepsize is enable (1) or disable (0)
    • the third number is the time step 10−9 s. In case of adaptive stepsize it is your initial timestep guess afterwards it will be adapted automatically.
    • The Velocity Verlet (VV) integrator calculates wc=qB/m instead of calculating it (see: http://dx.doi.org/10.1006/jcph.1999.6237 ). Since wc is normally the fastest motion in such a simulation usage of VV is around 7-15 times faster than using another integrator. However one should be aware that the energy can run away as well. As a rule of thumb; for a 39K 1+ ion in a 6T magnetic field 1e-8 s is a good timestep. Adaptive stepsize is not implemented.
    • The Gear integrator is a fifth order leapfrog algorithm. This means the force is calculated just once every timestep making it a fast algorithm. The integration if fifth order though. Adaptive stepsize is not possible so be aware of energy non-conserving, i.e. check how low you can make the timestep

    • The Dormand Price and Runge Kutta algorithm do include adaptive stepsize and are therefore also energy conserving. Using the DP should be the standard choise if you do not care about speed.
      For example: 

      ODE 4 1 1e-9
      • RK4 with adaptative step size is used with an initial time step of 1 ns

  • 5: Coulomb interaction parameters, COULOMB : 2 integer

    • the first integer indicates if the Coulomb interaction is calculated (1) or not (0).
    • If Coulomb interaction is calculated, a second number, which is the Coulomb scale factor, has to be added

  • 6: output parameters, OUTPUT :

    • first parameter: The prefix name of the output files
    • second parameter: the interval time between each output (printstep)
    • third parameter: printing all output in one large file at each printstep (0) or print create a separate outputfile for each particle (1) or print all output in one large file at the end of each operation (2)

    (trunk version only) Data can be recorded at a Z crossing plane using (3) as third parameter. The second parameter is then the Z position of this plane. The data are recorded in one file.

  • 7: Trap configuration parameters: either one of the following parameters:

    • IDEALTRAP: followed by 3 doubles. The electromagnetic field in the trap is ideal:
      • the first double is the radius of electrodes (m)
      • the second double is Ud2=U0/d^2. U0=Voltage difference between ring electrode and endcap electrode, d is the characteristic trap depth. Ud2 hence defines the longitudinal frequency ω_z^2=Ud2*e/m
      • the third double is the magnetic field strength (Tesla)
    • NONIDEALTRAP: followed by the name of the trap config file. The electromagnectic field takes into account the correction factors c4, c6 and b2 (see Trap configuration file for the used equations)
      • a trap config file example is given in the simulation folder.
    • REALTRAP: followed by the number of files and the names of 2 or 3 map files:
      • 2 files: first file: Electric field map, second file: Magnetic field map
      • 3 files: first file: radial Electric field map, longitudinal Electric field map and third file: Magnetic field map
      • NOTE: The nr. rows and columns of the fieldmap is hard-wired in force.cpp (line ~20-30). One needs to change these values and recompile if these values change.
    • Example of the output of a COMSOL magnetic field file:
      % r z mf.Br (T) mf.Bz (T)
      5e-4 -0.7485 8.491617e-5 -5.095829e-4
      5e-4 -0.7455 8.499792e-5 -0.001529
      5e-4 -0.7425 8.509695e-5 -0.00255
      ...

  • 8: rf-excitations parameters: multiple line input. Excitations are possible via the following parameters:

    • NE No excitation:
      NE time_operation
    • DEW/QEW/OEW Dipole/Quadrupole/Octupole Excitation without buffergas
    • DEB/QEB/OEB Dipole/Quadrupole/Octupole Excitation with buffergas
      DEW time_operation eigen_freq_letter Element frequency_bias amplitude
      or
      DEW time_operation f frequency amplitude
    • RWB Rotating wall excitation with buffer gas:
      RWB time_operation order f frequency amplitude
    • RWW Rotating wall excitation without buffer gas:
      RWW time_operation order f frequency amplitude
    • AWB Anti-Rotating wall excitation with buffer gas:
      AWB time_operation order f frequency amplitude
    • AWW Anti-Rotating wall excitation without buffer gas:
      AWW time_operation order f frequency amplitude

    • AR auto resonance excitation:
      AR time_operation amplitude frequency_start frequency_stop:

    • For the excitations, the parameters are:

      • time operation : time of operation (in s)
      • eigen_freq_letter : c for cyclotron frequency, m for magnetron frequency, p for reduced cyclotron frequency. With Pc, Pm and Pp, the time of operation is set at a period of the chosen eigen frequency multiplied by the given time (trunk version only).
      • Element: the eigenfrequency of this isotope is used (ex: 133Cs, 7Li etc...)
      • frequency bias: with this the frequency shift (rad/s)
      • frequency : frequency of the excitation
      • amplitude: with this amplitude (V)

    • For the rotating wall and anti-RW, the parameters are:

      • time operation : time of operation (in s)
      • order: 1 for dipolar, 2 for quadrupolar, 3 for octupolar
      • frequency (rad/s)
      • amplitude (V)

    • note: Simbuca calculates the frequencies for a excitation (for example the wc frequency). However when using the NONIDEALTRAP option and a non-number excitation frequency Simbuca doesn`t calculate the magnetic/Electric field to get the excitation frequency (too time consuming). So one should type in the frequency as a number. Another possibility is changing the default trapparameters in the default constructor in trapparameters.cpp.

Trap configuration file

The NONIDEALTRAP flag of the .sim script file needs a Trap configuration file composed as follows:
Trap_config.dat:

#
# COMMENTS can be inserted if they start with the # character
#
#
MAGNETIC FIELD
B0 (B2)
ELECTRIC POTENTIAL
U0 c2 (c4) (c6)
ELECTRODE RADIUS
Relec
CORRECTION FACTOR
a
PUMPING DIAPHRAGM
0.02
  • Variables in brackets are optionnal
  • B0, magnetic strengh (T)
  • B2, magnetic bottle factor (T/m^2)
  • U0, potential (V)
  • c2, harmonic factor (1/m^2)
  • c4, c6 anharmonic correction factor (1/m^4, 1/m^6)

  • Relec, electrode radius (m)

  • inhomogeneous magnetic field expression:
    B_z(r,z) = B0 (1 - B2 (z^2 - r^2/2) )

  • Real trap potential expression:
    V(r,z) = U0 c0 + U0 c2 (z^2 - r^2/2)+ U0 c4 (z^4 - 3r^2 z^2 + 3/8 r^4) + U0 c6 (z^6 - 15/2 r^2 z^4 + 45/8 r^4 z^2 - 5/16 r^6)

.sim file example

You can find below 2 example of .sim files. The first one which is steps-by-step explained with the parameters as described in the previous section.

First example (explained)

foo.sim:

# this is a comment
CREATECLOUD 2048 3
CLOUDPARTS 7Li 5. 6Li 80. 39K 15. 
CLOUDCOORD 0.002 0.001 0.007 0.0001 0. 0.0005
TEMP 1.5
BUFFER 1 1e-8
ODE 4 1 1e-9
COULOMB 1 10
OUTPUT output_file 1e-5 1
IDEALTRAP 0.02 5e4 4.7
NE 0.001
DEW 0.01 m 7Li 10 0.5 
QEB 0.05 c 7Li -5 1.5
  • CREATECLOUD: 2048 is the number of simulated particles, there are 3 species
  • CLOUDPARTS: cloud is created with the following parameters:
    • the isotopes used are 7Li (5%) 6Li (80%) and 39K (15%). The total amount has to be 100 %. The masses are loaded from the AME file: mass.mas03.
  • CLOUDCOORD: specifying the cloud coordinates
    • the cloud has a spheroidal shape with semi axis lx=0.002, ly=0.001 and lz=0.007 m.
    • the center of the cloud is at the position (0.0001, 0., 0.0005)
  • TEMP: specifying the mean energy
    • A Maxwell-Boltzman energy distribution is used with a mean energy of 1.5 eV.
  • BUFFER: specifying the buffer gas
    • 1 means that buffergas is enabled. The pressure of the of the buffer gas is 10−8 mbar, buffergas is included (1) with a pressure of 1e-8 mbar.
  • ODE: The 4th order Runge Kutta integrator is used (4) with adaptive stepsize (1)and a first step of 1e-9 s.
  • COULOMB: Coulomb interaction is enable(1) with a scaled coulomb factor of 10.
  • OUTPUT: The outputfiles are created with prefix "output_file". Particles are printed out every 10 us and the information of each particle is printed out in a separate file (1),
  • IDEALTRAP: The electromagnetic field in the trap is ideal, its parameters are : electrode radius = 2 cm, Ud2 =5e4, B=4.7 T.
  • At first no excitation (NE) is applied during 1 ms.
  • Then a dipolar excitation (DE) without buffer gas (W) is done during 10 ms at the magetron frequency of 7Li with a shift of 10 rad/s and an amplitude of 0.5 V.
  • Afterwards a quadrupolar excitation (QE) with buffer gas (B) is done during 50 ms at the cyclotron frequency of 7Li with a shift of -5 rad/s and an amplitude of 1.5 V.

Second example

foo2.sim:

IMPORTDATA prev_simulation_prefix
ODE 4 1 1e-9
COULOMB 1 1
OUTPUT output_file4 1e-7 0
IDEALTRAP 0.02 18000 6
OEW 0.01 c 7Li 10 0.5

Simbuca, a simulation program for Penning traps

Simbuca overview

Simbuca is a modular Penning (ion) trap simulation program written in C++. Due to the modularity different parts of the program can be used for other purposes. For example one can use the collision methods implemented in coll.cpp and compile it with another program to check how helium gas affects ion beam optics in beam tuning. Another example would be to use the fourth order Runga Kutta method in ode.cpp as an integrator for a totally different application where particle trajectories are calculated.
The different files that make up Simbuca will be discussed only briefly here. For more detailed information thre reader is referred to their implementation, described in the References section below.

  • mtrand.h is a function to create random numbers between 0 and 1.
  • globals.h contains all the physics constants used in Simbuca.
  • trapparameters.h contains the characteristic parameters, e.g. trap depth and magnetic field.
  • matrix.h is a C++ container of the matrix class.
  • particle.cpp is a class that contains positions and speeds in the three directions x,y and z.
  • ion.cpp is a class that contains the mass of the ion. This class uses the constants from Trapparameters.h to calculate the trap specific frequencies ω+,ω− and ωc.
  • coll.cpp handles the collisions. The speed of the ion is given as input and, if a collision has occurred, the updated speed after collision is returned. See the NIM paper for more information about the buffer gas collisions.
  • ioncloud.cpp is a class that contains an array (vector) of ions and an array of (corresponding) ions. The ioncloud is a static global variable, i.e. there is one ioncloud that all the seperate files are using.
  • fieldmap.cpp is used to read the fieldmaps of the magnetic and electrical fields in both radial and axial direction. The magnetic fieldmap can obtained from the manufacturer of the magnet. Electrical fieldmaps can be constructed in programs like Comsol.
  • force.cpp consists of two files. One for a simulation on the CPU and one for a simulation on the GPU. Depending on the variable pu in the Makefile either force gpu.cpp or force cpu.cpp will be copied to force.cpp and then compiled. force.cpp calculates the Coulomb force, the force due to the electrical and magnetic field in the Penning trap and, if applied, additional forces from excitations.
  • ode.cpp is the heart of the program. The ordinary differential integrators are implemented here. See section about integrators in the NIM paper for more information about the implementation of these integrators.
  • ionFly.cpp initialises all the seperate files when the program starts, e.g. the order of the integrator, or what buffer gas model to use. ionFly.h gives an overview of all the functions that can be used in main.cpp.
  • parser.cpp is used to import data from previous simulations and use these as the starting point for new simulations.
  • logger.cpp constructs a log file function used by all the other files in the program.
  • main.cpp is the only file that is changed if different simulations are performed. The user can set different environment parameters here like Coulomb interaction, boundary inclusion, order of the integrator, ...
  • the Makefile compiles Simbuca and can create backups of the program. One can change the compiler (g++ or icpc or mpic++) and the processing unit (GPU or GPU) here.

Contact

  • Simon van Gorp: simon.van.gorp@cern.ch // simonvgorp@gmail.com
  • Pierre Dupre: pierre.dupre@csnsm.in2p3.fr

Contributors

  • Logger classes: Attila Krasznahorkay
  • SimParser classes: Balint Radics

References

  • First publication on Simbuca:
    S. Van Gorp, et all. "Simbuca, using a graphics card to simulate coulomb interactions in a penning trap." Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, 638(1):192 – 200, 2011: http://dx.doi.org/10.1016/j.nima.2010.11.032
  • Follow up publication about developments between 2010 and 2012
    S.Van Gorp and P. Dupr ́e. Improvements of the simbuca trapped charged particle simulation program. In 10th workshop on non neutral plasmas, 2013: http://dx.doi.org/10.1063/1.4796087

Screenshot thumbnail
Simbuca is also used in othe applications were the Coulomb interaction and EM field tracking is required like e.g. an MR-TOF: Simulated ToF spectra of an initial mixture of CO + (red) and N 2 + (blue) ions (ToF compared to the non-interacting N 2 + reference ion). The trapping time increases from top to bottom of every part of the picture (details in picture). The Coulomb scaling factors are a: 0 (no Coulomb interaction), b: 5, c: 10, d: 20. Picture taken from M.Rosenbush: AIP Conf. Proc. 1521 , 53 (2013); doi: 10.1063/1.4796061
Screenshot thumbnail
Test of Simbuca by comparing it with the equations of motion of a dipolar excitation at the eigenfrequency n. Perfect correspondence between simulation and theory is found
Screenshot thumbnail
Movement of the ion cloud in the X,Y-plane of the trap (in mm). The left column represents 1000 ions without Coulomb interaction. Coulomb interactions between ions are taken into account in the right column. Here 1000 ions are each given a charge of 10^4+ to represent 10^7 ions. Ions move 50 μs between the frames. Due to the Coulomb interaction in the cloud, the ion rotates around the center of the cloud (right column). Without Coulomb interaction taken into account the ion stays at the same position inside the cloud (left column).

Project Admins: