accuracy problem

[Johanna]

Dear all,

I'm trying to track down a problem with my simulations. Therefore I created a very simple and stupid input file which gives strange results. The file you can find below.

This pulse sequence shouldn't do anything and I would expect the result to be exactly 1.000000. Here a few results I obtained with different timings (everything else I did not change):

Input                                             Result                              Comment
a)  timing(usec)        0                    0.9999999                         Shouldn't that be exactly 1.0000000?
b)  timing(usec)      100                   0.9999999
c)  timing(usec)     (100)2                0.9999913   0.9999913       Why is this result different from a) and b) ?
d)  timing(usec)     1000000             0.9999719                         Why is this result different from 1.00000 ?

The results depend on the powder averaging used. However I couldn't find a powder that works as expected.
I would be very grateful if somebody has an explanation for the behaviour I observe.

Thank you very much
Best wishes
Johanna


****** The System *******
spectrometer(MHz)  850
spinning_freq(kHz) 10
channels           N15
nuclei             N15
atomic_coords      *
cs_isotropic       0
csa_parameters     1 -109 0.2 0 0 0 ppm
j_coupling         *
quadrupole         *
dip_switchboard    *
csa_switchboard    *
exchange_nuclei    *
bond_len_nuclei    *
bond_ang_nuclei    *
tors_ang_nuclei    *
groups_nuclei      *
******* Pulse Sequence ******************************
CHN 1
timing(usec)        0
power(kHz)          0
phase(deg)          0
freq_offs(kHz)      0
******* Options *************************************
rho0               F1z
observables        F1z
EulerAngles        asgind30
n_gamma            *
line_broaden(Hz)   *
zerofill           *
FFT_dimensions     *
options            -t
******************************************************

[Mikhail Veshtort]

Hi Johanna,

The issue you are observing has to do with the fact that these computations are numerical, not analytical. They are done in finite floating point precision. If you run your file with the option -s, you will see a line that says "Machine epsilon is 1.192093e-07", which is a way to characterize the round-off error in these computations. In your simulation, the program first computes the propagator for one rotor cycle, uses this propagator to compute the propagator for the whole evolution time, applies the latter to the density matrix to propagate it, and then computes the signal with this propagated density matrix. This computation is performed for each crystallite and the results are summed up with the weights specified by the powder averaging set.

Now, 0.9999999 (in the floating point world) is different from 1 by exactly 1.192093e-07, which is actually very good for this kind of computation.

In your last example, the evolution time is 1 s, and so the rotor cycle propagator is raised to the 10000th power to get the propagator for the whole evolution time. Therefore, if the matrix elements of the propagator have even a single round-off error (1.192093e-07), it may grow by a factor of 10000 during this computation. Actually, if you run your simulation with the option -ve, the program will report you its estimates of these errors. Powder averaging improves these errors quite a bit because they are randomly distributed. So the final error observed in that computation (4e-4) looks quite reasonable.

Overall, one should realize that in the last example, you are asking the program to integrate the motion with a time-dependent (periodic) Hamiltonian to a time that is 10,000-fold of the period. This is quite a fit to perform in the single-precision floating point arithmetic used by the program. Also, note that if you were measuring the FID for 1 s and then plotting a spectrum, that would be a completely different computation, and this kind of errors would not be seen in the results.

Hope this helps.

Best regards,
Mikhail

> --
>
>

[Johanna]

Dear Mikhail,

thank you very much for your quick and detailed answer.

Do you also have an explanation why in the version c of my example the result is different from what I obtain in version a and b (which calculate the same time points but individually)?

Best wishes
Johanna

[Mikhail Veshtort]

Hi Johanna,

In your example c) a different method is used in the calculation. Here the evolution is computed in the spectral domain. In this setup, it is a bit less accurate than the other method.

Mikhail

--