UiB : MatNat : Informatikk : Undervisning : Kurs : BER200 : BER200 Thierry Matthey

BER200 - Problem 3

Ion Crystalization and Dynamics


The purpose of this exercise is to investigate the structure and dynamics of so-called ions trapped in a electric and magnetic field (Paul Trap). Ions in free space repel each other due to Coulomb force repulsion and expand to infinity. In a trap (Fig. below) an additional attractive potential keeps the system bounded. For certain parameters of the potential the most favorable energetic configuration (the state with lowest energy) of the ions will be to organize them self in a linear string. Such a string is at present one of the most promising candidates for implementing a quantum processor, which have the potential of solving classical exponential problems in linear time. In this exercise we will apply a classical model and simulation package to explore how many atoms we are able to align before our string system changes structure for a given set of experimental realistic trap parameters.

Paul Trap

Physical Problem Description:

On the atomic scale a set of particles is most accurately described with Schrödingers equation,


where H_hat is the Hamiltonian of the ions interacting with the time dependent electro magnetic fields and Psi - Wave funtion is the wave function. This problem is extremely difficult for more than 2 particles but fortunately for sufficient massive particles (like ions) at sufficient high temperatures, we may skip all quantum effects and attack it classically. In classical mechanics the dynamics of the particles is determined by Newtons 2. law,

Newton's equations of motion

The forces can be calculated from the potential energy function


The potential is a sum of electrostatic repulsions and trap attractions,

Potential (electrostatic repulsions and trap attractions)

where Positions is the coordinate of particle i.


The mass Mi is defined in [amu], qi, qj are defined in [e] and Coulomb factor.

From experimentalists we learn that a typical trap parameter:

Experimental Omega value

In the present exercise we apply a 10 times weaker . To obtain a stable configuration the experimentalists apply lasers in the z-direction which can cause the ions to slow down (Nobel Prize in Physics 1997). We will model this effect by introducing a friction force:

Friction force

With this extra force each particle will be slowed down and the system will move to equilibrium where all particles are stationary. Note that too strong friction will freeze the system before reaching equilibrium, and with too weak friction the system will oscillate without converging to equilibrium. An additional check of whether the final structure has reached equilibrium is provided by the fact that for spherical system the Coulomb energy of the system is equal to twice the trapping energy.
In order to cool down the system by an extra force, one could control the temperature with a Nosé-Hoover thermostat (1).

ProtoMol Quickstart:


ProtoMol is designed to be an important tool for Molecular Dynamics (MD) simulations. With a simple configuration setup, compatibility with other popular MD file formats, and the ability to run on parallel systems, ProtoMol combines high performance and ease of use.

What ProtoMol provides is a generic, object-oriented component framework for MD simulations. To meet these high performance expectations, ProtoMol uses cell algorithms, grid techniques and well-optimized libraries to challenge the most computationally expensive forces. The design of ProtoMol also includes parallelization, based on components to distribute the work and data. The approach follows an incremental and partial parallelization scheme, which allows the developer to start with a sequential implementation and then do step by step parallelization.

The overall framework of ProtoMol is designed for non-bonded, bonded, short-range and long-range forces for systems with tens of thousands of atoms representing water and several large molecules. It is designed for high flexibility, ease in extension and maintainence, and high performance demands.

MD describes a molecular system as a function of time based on the integration of Newton's equations of motion and interacting forces. The integrators are the part of the program that solve the differential equations that describe the system. Specifically, the integrators provide a set of forces, that describe the system at each time. Thus, it is easy to see that the central part of the entire system is the integrators with their force definitions.

ProtoMol provides several integrators (including multiple-time-stepping) and a variety of forces (fast approximations and exact algorithms) to describe and solve Newton's equations of motion or more advanced MD equations. The definition of integrators and forces are part of the configuration file, such that there is no need of re-compilation to change the integrator scheme or the set of forces. All this together makes ProtoMol a favorable MD framework to experiment with different integrator schemes and/ or different force(-algorithms)/ potentials in a sequential or parallel environment.

Download - Binaries and Source: Release 1.8.3

Simple Test Case I: 2-Type/Shell 3Mg+24 7Ca+40 (Leapfrog and friction)

3mg7ca.conf definition of I/O files, number of simulation steps, integrators, forces.
3mg7ca.pos.xyz initial positions of the particles, XYZ-format.
3mg7ca.psf charge and mass definitions (also the topology of the molecules).
3mg7ca.par non-bonded (CHARMm) parameter definitions.
2K archive
protomol 3mg7ca.conf --XYZPOSFILE 3mg7ca.out.trajectory.xyz

Simple Test Case II: 1-Type/Shell 20,288Ca+40 (Nosé-Hoover and Multigrid)

20288ca.config definition of I/O files, number of simulation steps, integrators, forces.
20288ca.cmp.config configuration file to compare multigrid and direct method.
20288ca.pos.xyz initial positions of the particles, XYZ-format.
20288ca.psf charge and mass definitions (also the topology of the molecules).
20288ca.par non-bonded (CHARMm) parameter definitions.
56K archive

Running ProtoMol

./protomol 3mg7ca.conf
In order to get help:
./protomol -h
To figure out the supported forces:
./protomol -f
To figure out all possible keywords:
./protomol -k

ProtoMol writes DCD files, which can be read by VMD or XYZ trajectory files visualized by the OpenGL/GLUT viewer (xyzviz.cpp).

ProtoMol Units

ProtoMol Units

1 Velocities in PDB files are scaled to get a more accurate representation since the PDB format has a limited representation of floating numbers.

More Information

  • TTP5: Application - 2 Coulomb Crystals

  • User Guide 1.08: PROTOMOL - An Object-Oriented Molecular Dynamics Framework (PS) (PDF) (Jesus A. Izaguirre with Thierry Matthey).

  • ProtoMol Tutorial (PDF) (T. Matthey), Nov. 2002.

  • Consult the ProtoMol site at sourceforge, you may also find more information at UiB or ND, USA.

Problem 1:

  • For a two particles, calculate the equilibrium distance between the two Ca+40 ions, m=40.08, q=1.0.
  • Verify your calculation by running ProtoMol with the two-ion example. Download VMD (DCD trajectories) to visualize the behavior or use the OpenGL/GLUT viewer (xyzviz.cpp) (XYZ trajectories).
  • Set up a series of simulations with increasing number of ions and run ProtoMol until the initial string forms a bulk.

Problem 2:

  • Find the right Paul trap parameters such that up to 41 Ca+40 ions a convergence to a string is seen, where 42 or more ions form a zig-zag line or a bulk. Down load and use the system files (PSF, XYZ-positions and PAR) for 41 and 42 ions. For the configuration files set = 1e-10[fs-1], temperature = 1e-6K and use the Leapfrog integrator and a friction force.
  • Movie of 41 Ca+40 ions converging to a string (MPEG, 322KB).
  • Movie of 42 Ca+40 ions converging to a zig-zag line (MPEG, 320KB).
  • Analyze and visualize the structure using MatLab and a XYZ-reader (xyzread.m):

    x = xyzread('42ca.out.fin.pos.xyz');

    or compile the OpenGL/GLUT viewer (xyzviz.cpp):

    xyzviz -movie 50ca50a.out.all.trajectory.pos.xyz

    xyzviz reads and displays/animates ordinary XYZ position or XYZ trajectory files, and has several output possibilities and features (PPM, EPS), e.g., MPEG:50 Ca+40 + 50 A2+80 ions convering to a sphere, 0.1ms, 13MB.

Problem 3:

  • The real scientific interest is in large number of particles. Run some test for different values of n and develop a model for the time complexity for the code. Assume you got the new IBM regatta system at your disposal for a week, how large system can you simulate?
  • If we increase the number of ions beyond the breakover from string to bulk other forms will appear. Set the number of ions to 10000 and compute until a stable form is achieved. How does this form look like?
  • Now the problem has become computational expensive. Thus we need fast force calculation (cut-off, multigrid) and efficient time integration. Play with different force calculation and time integrators (impulse) and find one that balance efficiency with accuracy well.

Related Links:


The report should be published on the web and must include:
  • Description of the project.
  • Description of numerical methods used and answers to the given questions.
  • Conclusion.
  • References.