Aarhus University Seal


SIMPSON 4.1.1 released

SIMPSON v 4.1.1 is a bug fix release. The major changes concern the acq_block command. Nov 25, 2014


SIMPSON v 4.1 is a minor release. The major changes are

- Compiled with Tcl 8.6

- Fixed bug in storage of unit propagator

- Fixed bugs in fint

SIMPSON version 4

Released September 9, 2013

We are pleased to release a major update of simpson. The simulations are significantly faster than previously, and many new features are added. Use the links on the right to download your flavour or download the source if you want to compile SIMPSON yourself. Please note that the new version of SIMPSON uses two new libraries, fftw3 and nfft3.

Section spinsys

dipole_ave N1 N2 bIS eta alpha beta gamma - for motionally averaged dipole-dipole interactions this allows to define asymmetry parametereta

Section par

method value1 value2 value3 ... - defines calculation methods for fsimpson.

value feature
direct, idirect both values equivalent, uses direct evaluation of pulseq, possibly with COMPUTE method if acq_block is used*)
gcompute, igcompute both values equivalent, uses gamma-COMPUTE method, either in original pulseq syntax or in new syntax with acq_block – allows for manipulation before acquisition)*)
time forces calculations to be done in time domain
frequency forces calculations to be done in frequency domain – can be used only in combination with acq_block (both direct or gcompute)
block_diag tries to use block-diagonal Hamiltonians whenever it is possible; considerable speedup for free evolution periods and when rf not applied on all channels
FWTinterpolation triggers FWT interpolation method to be used for both amplitudes and frequencies – works with gcompute only
FWT2interpolation triggers FWT interpolation method to be used just for frequencies, amplitudes are taken the same as in the source orientational set (nearest) – works with gcompute only
ASGingterpolation triggers ASG interpolation method used after spectrum generation - works with gamma-COMPUTE in frequency domain only
FWTASGinterpolation, FWT2ASGinterpolation combinations of the interpolation methods
sparse uses sparse matix algorithms, needs simpson linked with Intel MKL libraries
diag, pade, chebyshev, taylor, lanczos experimenting with different methods to calculate propagators (matrix exponentials)
*) For a description of possible combinations of method instructions see below.

num_cores N – defines number of threads to be used for parallel calculations of parameter averages (mainly crystallite orientations). If not defined then all processor cores are used. Number of threads can also be limited by environment variable SIMPSON_NUM_CORES

crystal_file filename N1 N2 – optional arguments N1 and N2 define index of the first and the last orientation to be read from the filename

crystal_file filename1 filename2 – used for FWT interpolation, defines source set of orientations (filename1) and larger target set of orientations (filename2). At the moment, only hemisphere Lebedev sets are allowed.

rfprof_file filename

rfprof_file filename N1 N2 – optional arguments N1 and N2 define index of the first and the last parameter to be read from the filename containing rf inhomogeneity profiles

averaging_file filename

averaging_file filename N1 N2 – optional arguments N1 and N2 define index of the first and the last parameter to be read from the filename containing set of parameter values over which the calculation should average. The filename is formatted as follows: first line defines parameters in the same way as accepted for the fsimpson command, other lines define their values. Example:

shift_1_iso dipole_1_2_aniso weight
-10p -10000 0.2
+10p -10000 0.2
0 -10000 0.2
0 -8000 0.2
0 -12000 0.2

point_per_cycle N – used for direct method in frequency domain, defines number of points to which the one cycle of acquisition Hamiltonian should be split

ED_symmetry val – controls how to deal with possible excitation-detection symmetry: if val is -1 then it is not considered, if val is 0 then it is internally checked and used, if val is 1 then it is not checked and it is directly considered (user responsibility). Default is -1.

oc_grad_level N – defines accuracy level of optimal control gradient, default is 1 but L-BFGS methods requires at least 2 or higher

oc_method – defines optimal control method, can be CG for conjugated gradients (default), or L-BFGS

Section pulseq

acq_block {…} – allows easy definition of pulse manipulations during acquisition. The pulse sequence defined in parenthesis should contain one period of the total Hamiltonian. For gamma-COMPUTE the total length has to be the rotor period divided by the number of gamma-angles. For COMPUTE the pulse sequence is repeated such that it fills up multiple of rotor periods.


gcompute can be used whenever the pulse sequence is synchronized over a period equal to the rotor period divided by the number of gamma angles. Otherwise, the following methods are possible

static spinning gcompute
direct direct gcompute
direct freq direct freq gcompute time FWTinterpolation
direct freq ASGinterpolation gcompute freq
direct freq FWTinterpolation gcompute freq ASGinterpolation
direct freq FWTASGinterpolation gcompute freq FWTinterpolation
gcompute freq FWTASGinterpolation
gcompute freq FWT2interpolation
gcompute freq FWT2ASGinterpolation

all of these methods can be combined with your choice of block-diagonalization, sparse matrix handling, acq_block, etc.

Version 3.1 released March 25, 2011

Version 3.1 is a bug-fix release to version 3.0. The bugs were in the OC routines and for the second-order quadrupolar interaction during delays.

New features in version 3.0:

This section describes add-ons to the functionality of SIMPSON. All previous features are kept without any change except everything executes much faster for larger spin systems (more than 2 spins). Such speed up was gained thanks to replacing original matrix manipulations hard coded in SIMPSON by full use of BLAS and LAPACK pre-built (and performance optimized) libraries. Extra speed is achieved by implementing fast acquisition algorithms and parallelization.

Parallel calculations

There are implemented two ways of parallel calculations. The first one makes use of fork mechanism supported on Linux and Mac (not MS windows). The program automatically detects number of processor cores and starts appropriate number of instances. This number can be modified by setting the environment variable SIMPSON_NUM_CORES to a given number.

The second way is based on MPI and needs SIMPSON properly built from sources.  Before installing form sources, make sure the following prerequisites are installed:

  • TCL
  • OpenMPI

According to the specific versions and flavours of these prerequisites, it may be necessary to adjust the makefile. The makefile directives below show how to include and link to these packages. To enable MPI, add '-I/path/to/mpi/include' to the include list, '-lmpi' to the library list, and '-DMPI' to the flags list, i.e.:

INCLUDES = -I/usr/include/tcl8.5 -I/usr/lib/openmpi/include
LIBRARIES = -lm -lgslcblas -llapack -ltcl8.5 -lmpi
FLAGS = -DMPI -c -O3

Now, type make to build SIMPSON. Once built, SIMPSON can be run with MPI using mpirun. Make sure to start as many MPI processes as there are processor cores. For instance, in a cluster setup with 8 machines, each equipped with 4 processor cores, one should start 8*4 processes (the fork mechanism from above is not used when using MPI): 

mpirun -np 32 simpson <inputfile>

New direct acquisition method

Works for static as well as MAS simulations. Pulse sequence includes new command
acq_block {...}
which defines the whole acquisition period of the simulated experiment. It calculates the whole fid at once. In the curly braces one can define any pulse sequence applied during the acquisition. The sequence is coded using standard simpson commands but must obey certain rules.

  • Static spectra: durations of pulses and delays must sum up to one dwell-time. (this may be released in coming versions, was not priority this time). Propagator of the dwell-time is diagonalized and re-used in the diagonal form.
  • MAS: pulses and delays must be separable into durations of the dwell-time. Simpson then uses variant of COMPUTE algorithm for fid calculations (making use of Hamiltonian periodicity).

In the input file, use coding like this:

par {

  method idirect



proc pulseq {} {

  global par


  pulse ...


  acq_block {delay $dw}

for decoupling on some channel, or

  acq_block {

    pulse $tp1 ...

    pulse $tp2 ...

    delay $td3 ...



for some periodic rf pulse sequence repeated every two dwell-time periods (and that would work just for MAS at the moment). Here one should be careful with synchronization, when the pulse sequence has to be repeated many times to match MAS, efficiency is lost.

New version of g-COMPUTE

The original simpson did not allow for preparation of density matrix before acquisition. The spin state had to be independent of the g-angle when acquisition starts. This is now released, again using acq_block {} command. Pulse sequence in parenthesis has to have duration of one rotor period divided by the number of g-angles (this originates from ideas behind g-COMPUTE, only then efficient averaging can be used). However, spectral width does not need to be spin_rate*gamma_angles - it just needs to stay simply synchronized like in the example. The following code is possible

par {


  method igcompute

  gamma_angles 10

  spin_rate 12000

  sw spin_rate*gamma_angles/3



proc pulseq {} {

  global par


  pulse ...


  acq_block {

    delay [expr 1.0e6/$par(spin_rate)/$par(gamma_angles)]



Some tricks cutting unnecessary operations has been used and this implementation is not only more general but also slightly faster (20% for 5 spins) than the original gcompute code.

General notes

  • For safe calculations, use only one acq_block inside the pulseq. Think of applications where you would like to use it repeatedly and tell me. One obvious candidate is calculation of 2D spectra which is not implemented yet.
  • Do not use acq_block inside any loop.
  • At the moment, there are no error checks for pulseq code. Combinations of acq_block with original direct, gcompute and gammarep as well as combining idirect  and igcompute with original acq leads to undefined behavior.
  • Inside the parenthesis of acq_block one can use also any variables defined in previous lines of pulseq code (I hope). This example should work:

      proc pulseq {} {
        global par
        set dw [expr 1.0e6/$par(sw)]
        acq_block {delay $dw}

    as well as this (but there is no advantage of doing it like this, actually it is recommended not to do so due to propagator synchronization with MAS)
      pulseq {} {
       global par
       delay [expr 1.0e6/$par(sw)]
       store 1
       acq_block { prop 1}
  • The new algorithms were developed with method names direct_new and gcompute2_new. It may be scrambled in error messages or even in input files. Please be patient, the SIMPSON code is not yet ready for stable 3.0 release.