Something is not right in my script...

[Svist]

Dear all,

What am I doing wrong?

I am trying to optimise a free radical X-band CW spectrum simulation with six paramagnetic nuclei,
Sys.Nucs = '1H,1H,1H,1H,14N,1H'

My problem is that when I run this m-file (see below) to optimise just the three g-values, the script performs exactly 11 steps (this number is independent of the starting gx, gy, gz values) and reports the initial parameters as the best ones. Which is clearly wrong.

I suspect I am making some syntax errors. Any suggestions where I should be looking into? Many thanks for your help.

My script:

clear;
[B,spc] = textread('myfile3.txt','%f %f'); % reads my exper data file saved as two columns, B and spc
field = B'; % makes the field a row
a4 = spc'; % makes the experimental spectrum values a row
ScaleExp = 0.3; % Scaling factor for experimental spectrum when plotting the two together
% Experimental parameters for simulation xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
Exp.mwFreq = 9.394 ; % <-- this is in GHz
Exp.Range = [328 339.98]; % <-- alternative, left and right fields, in mT
Exp.nPoints = 1024;
% Experimental parameters for simulation end xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Am = 3; % modulation amplitude in Gauss
cAm = 0; % positive or negative correction value to Am, in Gauss
h = planck; % Planck constant in m² kg / s (or J s)
Bohr = 9.27401E-24; % Bohr magneton in J/T

Sys.gx = 2.004;
Sys.gy = 2.004;
Sys.gz = 2.004;

giso=(Sys.gx+Sys.gy+Sys.gz)/3; % isotropic component

GaussToMHz= giso*Bohr/h/10000/1000000; % Gauss to MHz conversion factor, in MHz/G, 10⁴ converts from T to G, 10⁶ - from Hz to MHz

deltaH_Wx= 3.8384; % Individual width values, in Gausses
deltaH_Wy= 6.3277;
deltaH_Wz= 3.7467;

deltaHx=deltaH_WxGaussToMHz; % Individual values of delta H, in the W version, in MHz
deltaHy=deltaH_WyGaussToMHz;%
deltaHz=deltaH_Wz*GaussToMHz;%

% All A-values are in MHz (the values are results of calculations, therefore so many decimals…):

A11 = 76.207424038009280; % A-values for 1st proton (beta 1)
A12 = 68.100361913024200;
A13 = 68.100361913024200;

A21 = 36.394251806515356; % A-values for 2nd proton (beta 1)
A22 = 31.941868915428213;
A23 = 31.941868915428213;

A31 = -15.1343430849; % A-values for 3rd proton (on C5)
A32 = -12.4432445592;
A33 = -4.93779546;

A41 = -26.985989745; % A-values for 4th proton (on C7)
A42 = -19.15739631;
A43 = -6.109740585;

AN1 = 1.1781203175; % A-values for nitrogen (N)
AN2 = 1.3277228975;
AN3 = 33.249173405;

A51 = 0.648; % A-values for 5th proton (on C6)
A52 = 1.836;
A53 = 3.924;

euler11 = -5.904431025945578; % Euler angles for proton beta1
euler12 = 30.609359340719950;
euler13 = 83.530674424679120;

euler21 = -23.887785712648018; % Euler angles for proton beta2
euler22 = -9.929009124573753;
euler23 = 39.434871913301170;

euler31 = -6.56; % Euler angles for third proton (on C5)
euler32 = 0;
euler33 = 90;

euler41 = -77; % Euler angles for fourth proton (on C7)
euler42 = 0;
euler43 = 90;

eulerN1 = -61.26; % Euler angles for Nitrogen
eulerN2 = 84.97;
eulerN3 = 16.27;

euler51 = -82; % Euler angles for fifth proton (on C6)
euler52 = 0;
euler53 = 0;

Sys.S = 1/2;
Sys.HStrain = [deltaHx deltaHy deltaHz];
Sys.Nucs = '1H,1H,1H,1H,14N,1H';
Sys.A = [A11 A12 A13;
A21 A22 A23;
A31 A32 A33;
A41 A42 A43;
AN1 AN2 AN3;
A51 A52 A53;
];

Sys.AFrame = [euler11 euler12 euler13;
euler21 euler22 euler23;
euler31 euler32 euler33;
euler41 euler42 euler43;
eulerN1 eulerN2 eulerN3;
euler51 euler52 euler53;
]*pi/180;

Vary.gx = 0.002;
Vary.gy = 0.002;
Vary.gz = 0.002;

% Vary.A11 = 0;
% Vary.A12 = 0;
% Vary.A13 = 0;

% Vary.A21 = 0;
% Vary.A22 = 0;
% Vary.A23 = 0;

% Vary.A31 = 0;
% Vary.A32 = 0;
% Vary.A33 = 0;

% Vary.A41 = 0;
% Vary.A42 = 0;
% Vary.A43 = 0;

% Vary.AN1 = 0;
% Vary.AN2 = 0;
% Vary.AN3 = 0;

% Vary.A51 = 0;
% Vary.A52 = 0;
% Vary.A53 = 0;

% Vary.euler11 = 0;
% Vary.euler12 = 0;
% Vary.euler13 = 0;

% Vary.euler21 = 0;
% Vary.euler22 = 0;
% Vary.euler23 = 0;

% Vary.euler31 = 0;
% Vary.euler32 = 0;
% Vary.euler33 = 0;

% Vary.euler41 = 0;
% Vary.euler42 = 0;
% Vary.euler43 = 0;

% Vary.eulerN1 = 0;
% Vary.eulerN2 = 0;
% Vary.eulerN3 = 0;

% Vary.euler51 = 0;
% Vary.euler52 = 0;
% Vary.euler53 = 0;

Opt.Method='perturb';
[BestSys,BestSpc] = esfit('pepper',spc,Sys,Vary,Exp,Opt);

Script ends

(I don’t think you would need to see the myfile3.txt file, the experimental spectrum.

Your comments will be appreciated
Dima Svistunenko

[ykliao]

Hello,
I think instead of

Code: Select all

Sys.gx = 2.004;
Sys.gy = 2.004;
Sys.gz = 2.004;

you should use

Code: Select all

gx = 2.004; gy = 2.004; gz = 2.004;
Sys.g = [gx, gy, gz];

and by the way, there are typos where you calculate the inhomogeneous linewidth

Code: Select all

deltaHx=deltaH_WxGaussToMHz; % Individual values of delta H, in the W version, in MHz
deltaHy=deltaH_WyGaussToMHz;%
deltaHz=deltaH_Wz*GaussToMHz;%

where multipications between the variables are missing for Hx and Hy.

I hope this helps you solve the problem.

[Svist]

Dear ykliao,

Thank you for your reply and suggestions, Certainly thank you for spotting a typo. I have advanced a bit during the weekend in understanding of what I am doing, so I hope I will get there soon. I will report when I am happy (or when I stuck again...)

Dima