S-Wave Scattering

The Diagrammatic Monte Carlo method (DiagMC) is a relatively new approach of solving differential equations, which are based on integral expressions. Since my PhD thesis studies the application of a more advanced version (the BMC, or Bold Diagrammatic Monte Carlo, method) to examine observables within the theory of Quantum Chromodynamics, I have to become familiar with the DiagMC method.

To apply the algorithm we need to prepare our analytical expression. We need to identify various parts and calculate some values. In this case we pick a spherical symmetrical potential U(r) with the exact form

Now we have to calculate the Fourier transform of this potential, i.e.

which results in

Here we also need to calculate an expression for small q, i.e. q→0,

Now that we prepared a few functions we need to work out our differential equation, which needs to be evaluated by using the Diagrammatic Monte Carlo method.

We want to solve a differential equation for f(q),

Here we can use that

So we identify:

• D₁(q) = -u(q)
• D₂(q,q',χ) = - 1/π u((q²+q'²-2qq'χ)^1/2)f(q')

Calculating the normalization I can be a tough problem. Introducing another (non-physical) diagram we can simplify this:

• D₀ = const.

Now the number of order 0 diagrams Z₀ is helpful. Using Z₀ / Z_steps approaches D₀ / I (here I is the normalization, which takes D₀ into account!). Using this we can estimate

where H_s is the value (sum of signs) in the s-th bin.

Finally we have to plugin everything into a working code, thinking of some (statistical) important steps:

• Everything has to be ergodic, i.e. we need to generate Markov chains that follow the equation of detailed balance.
• We need to normalize everything properly.
• Sufficient statistics are required before using the DiagMC sampling.

Considering these steps and following the DiagMC method we obtain the following algorithm:

1. Initialize all variables
2. Perform an infinite (or long enough) outer loop
3. Measure in each step
4. Perform an accept / reject step in every iteration
5. Maybe print or save the statistics
6. Maybe end the loop

The measuring is quite straight forward. Here we have to do the following:

1. If we are in the 0th-order we exclude the measurement (this diagram was brought in artificially and must be excluded from measurements) and increment the value of the number of 0th-order measurements
2. Otherwise we look at the current value of q and select the corresponding index for binning. The index can be calculated by (int)(q / Δq), where Δq = q_max / N_bins. The value of the selected bin has to add the value of the current sign.

The accept / reject step is just an applied Metropolis algorithm, with respect to the current and the old diagram. Here we perform:

1. If we are in the lowest order we compute the value of the next higher order
2. If we are in the highest order we compute the value of the next lower order
3. Otherwise we decide between all possible n changes with a chance of 1/n for each change and compute the value of the picked order
4. Now we calculate the ratio of the computed value over the current value
5. If a random number (between 0 and 1) is smaller than this ratio we accept the step: the new order is the picked one and the current value is the computed one
6. Finally we have to save the sign of the current order

The full C code can be downloaded below. Using the GNU C compiler (gcc) we can built the code (swave.c) by calling:

gcc swave.c -lm

So we must include the math library for the code to run. If you use a different compiler (like the one from Intel, icc), then remember to specify the correct parameters.

Download swave.c (5 kB)