SIMPSON coding examples

SIMPSON Tutorial #7: Fitting an experimental spectrum

by Thomas Vosegaard

In this snippet, we will first create an "experimental" spectrum and then fit it using SIMPSON.

Add comment

SIMPSON Tutorial #7: Creating an "experimental" spectrum

by Thomas Vosegaard

We will create our own “experimental” spectrum by adding Gaussian noise to a simulated spectrum. Open the file 7_Exp.in and run the simulation. You may play around with the amount of noise by changing the line variable noise 0.003 in the par section of the input file.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
spinsys {
  channels 13C
  nuclei 13C
  shift 1 0 4000 0.4 0 0 0
}

par {
  variable index  1
  variable noise  0.003

  method          gcompute
  gamma_angles    16
  spin_rate       1000
  crystal_file    rep66
  start_operator  I1x
  detect_operator I1p
  np              512
  sw              gamma_angles*spin_rate
  verbose         1101
  variable tsw    1e6/sw
}

proc pulseq {} {
  delay 9999
}

proc grand {} {
  # Generate gaussian noise
  set s 2
  while {$s >= 1} {
    set v1 [expr 2.0*rand()-1.0]
    set v2 [expr 2.0*rand()-1.0]
    set s [expr $v1*$v1 + $v2*$v2]
  }
  if {$s != 0} {
    set s [expr $v1*sqrt(-2.0*log($s)/$s)]
  }
  return $s
}

proc faddnoise {f {level 0.01}} {
  for {set i 1} {$i <= [fget $f -np]} {incr i} {
    set re [findex $f $i -re]
    set im [findex $f $i -im]
    fsetindex $f $i [expr $re + $level*[grand]] [expr $im + $level*[grand]]
  }
}

proc main {} {
  global par
  set f [fsimpson]
  faddlb $f 100 0
  faddnoise $f $par(noise)
  fphase $f -scale 2
  fsave $f $par(name),$par(index).fid
  fzerofill $f 1024
  fft $f 
  fsave $f $par(name),$par(index).spe
}
Files:7_Exp,1.spe  preview
7_Exp.in
Add comment

SIMPSON Tutorial #7: Fitting the experimental spectrum

by Thomas Vosegaard

Note that this part requires the opt-1.0 library.
Open the file 7_Fit.in and inspect the lines containing the instructions opt::newpar, opt::minimize, and opt::scan. The input file loads the spectrum 7_Exp,1.spe and optimizes the chemical shielding anisotropy magnitude and asymmetry parameter. Run the simulation and compare the spectra 7_Exp,1.spe and 7_Fit,1.spe. You may try to run the simulation for experimental spectra with different noise levels.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
# Remember to download the opt-1.0 package
lappend ::auto_path ./opt
package require opt

spinsys {
  channels 13C
  nuclei 13C
  shift 1 0 4000 0.2 0 0 0
}

par {
  variable index  1

  method          gcompute
  gamma_angles    16
  spin_rate       1000
  crystal_file    rep66
  start_operator  I1x
  detect_operator I1p
  np              512
  sw              gamma_angles*spin_rate
  verbose         0
  variable tsw    1e6/sw
}

proc pulseq {} {
  delay 10000
}

proc progress {} {
  global par
  if ![info exists par(progress)] { set par(progress) -1 }
  incr par(progress)
  set str {"*   " \
           " *  " \
           "  * " \
           "   *" \
           "  * " \
	   " *  "}
  return [lindex $str [expr $par(progress)%6]]
}

proc rms {{save 0}} {
  global par
  set f [fsimpson [list \
    [list shift_1_iso $opt::iso] \
    [list shift_1_aniso $opt::csa] \
    [list shift_1_eta $opt::eta] \
  ]]
  faddlb $f $opt::lb 0
  fzerofill $f 1024
  fft $f
  fautoscale $f $par(exp) -re
  set rms [frms $f $par(exp) -re]
  if {$save == 1} {
    puts [format " \[%s\]  %10.3f %10.3f %10.3f %10.3f" \
      FINAL $opt::csa $opt::eta $opt::lb $rms]
    fsave $f $par(name),$par(index).spe
  } else {
    puts -nonewline [format " \[%s\]  %10.3f %10.3f %10.3f %10.3f\015" \
       [progress] $opt::csa $opt::eta $opt::lb $rms]
  }
  flush stdout
  funload $f
  return $rms
}

proc main {} {
  global par mn
  set par(exp) [fload 7_Exp,1.spe]

  opt::function rms
  puts " Progress   csa       eta       lb          rms"
  opt::newpar iso 0 1 -10 10
  opt::newpar csa 3500 10 0 5000
  opt::newpar eta 0.4 0.1 0 1
  opt::newpar lb 150 10 0 300

  opt::scan iso
  opt::minimize 1.0e-6

  rms 1
}
Files:7_Exp_Sim.html
7_Fit,1.spe  preview
7_Fit.in
Add comment

Snippets version 1.0.1