# Optimal control OC-M2S pulse # by Christoffer Laustsen, Aarhus University 2012 # Spin system 2 spin = 1/2 spinsys { channels 13C nuclei 13C 13C shift 1 0 0 0 0 0 0 shift 2 0.13p 0 0 0 0 0 jcoupling 1 2 178.8 0 0 0 0 0 } # Parameters par { spin_rate 0 crystal_file alpha0beta0 gamma_angles 1 # Magnetic field proton_frequency 400e6 # define start operators in transfer start_operator Inz # define detect operators in transfer detect_operator I1z*I1z-I1x*I2x-I1y*I2y-I1z*I2z # number of elements in the shape variable NOC 1024 # duration of the shape (us) variable duration 80000 # J width for excitation in Hz variable J 20 variable Joff 178.8 # Chemical shift difference in ppm variable CS 0.07 variable CSoff 0.13 # width for excitation in ppm variable RF 1.0 variable RFoff 0.0 # Nutation scaling factor % variable B1 1.0 variable B1off 1.0 # Plotting variable cpI1 64 variable cpI2 64 # for creation of empty fid we need this sw 120000 # Parameters for optimization variable niter 20 oc_tol_cg 1e-5 oc_tol_ls 1e-3 oc_mnbrak_step 10 oc_max_iter 500 oc_var_save_proc rfstore # this MUST be used when working with optimal control!!! conjugate_fid false } proc pulseq {} { global par rfsh reset pulse_shaped $par(duration) $rfsh oc_acq_hermit } proc pulseqmon {} { global par rfsh reset pulse_shaped $par(duration) $rfsh acq } ### generate list of cp values distributed over +/-SW/2 range proc get_lims {SWH cp offset} { global par if { $cp <= 1} { set Res $offset } else { set step [expr double($SWH/double($cp-1))] set Res {} for {set i 0} {$i < $cp} {incr i} { set shft [expr $offset + double($SWH)/2.0-double($i)*double($step)] lappend Res $shft } } return $Res } ### optimal control procedures proc gradient {} { global par limsI1 limsI2 limsI3 set fs [fcreate -np $par(NOC) -sw $par(sw) ] # looping over J foreach shftI1 $limsI1 { # looping over cs foreach shftI2 $limsI2 { foreach shftI3 $limsI3 { set par(np) $par(NOC) #oc_gradmode on set f [fsimpson [list [list jcoupling_1_2_iso $shftI1] \ [list shift_1_iso [expr $shftI3]\p] \ [list shift_2_iso [expr $shftI2+$shftI3]\p]]] # sum up gradient over shifts and quadrupoles fadd $fs $f funload $f } } } return $fs } proc target_function {} { global par limsI1 limsI2 limsI3 ares set par(np) 1 set Res 0.0 set ares 0.0 set bres 0.0 # looping over J foreach shftI1 $limsI1 { # looping over cs foreach shftI2 $limsI2 { foreach shftI3 $limsI3 { set f [fsimpson [list [list jcoupling_1_2_iso $shftI1] \ [list shift_1_iso [expr $shftI3]\p] \ [list shift_2_iso [expr $shftI2+$shftI3]\p]]] set a [findex $f 1 -re] funload $f set Res [expr $Res + $a] set ares [expr $ares + $a] } } } return [format "%.20f" $Res] } ### MAIN to control the optimization and plotting proc main {} { global par rfsh limsI1 limsI2 limsI3 # optimization parameters (points in profile) set cpJ 0 set cpCS 0 set cpRF 11 # RF inhomogeniety #set par(rfprof_file) lorentzRF.rf #maximum RF power set maxRF 800 puts " " puts "****RUNNING OPTIMIZATION with : " puts " " # prepare list of offset control points set limsI1 [get_lims $par(J) $cpJ $par(Joff)] set limsI2 [get_lims $par(CS) $cpCS $par(CSoff)] set limsI3 [get_lims $par(RF) $cpRF $par(RFoff)] puts "J-couplings (Hz): $limsI1" puts "Chemical shifts differences (ppm): $limsI2" puts "RFoffset (ppm) : $limsI3" puts "****creating random sequence: " set rfsh [rand_shape 500 $par(NOC) [expr round($par(NOC)*0.2)] ] puts "****running optimization: " set tfopt [oc_optimize $rfsh -max $maxRF] puts "****saving shape : $par(name).dat " save_shape $rfsh $par(name).dat set rmspwr [shape_ampl $rfsh -rms] set avgpwr [shape_ampl $rfsh -avg] set maxpwr [shape_ampl $rfsh -max] puts "RMS power: $rmspwr" puts "AVG power: $avgpwr" puts "Max power: $maxpwr" puts "****cleaning shapes : " free_shape $rfsh free_all_shapes puts "****optimization done : " }