Comparison of magn. susceptibility results from ES and Phi
Posted: Mon May 30, 2016 11:31 am
I have now tried several models to fit magn. susceptibility data. To make sure my code is correct, I have treated the same data with the Phi package.
In particular, I fit an all-ferric triangle (S1 = S2 = S3 = 5/2) under an isosceles magnetic model, including antisymmetric exchange: H = -2J (S1S2 + S2S3) - 2J'S1S3 - 2d(S1xS2 + S2xS3 + S3xS1) + HZeeman. For this type of antiferromagnetic triangles, there are normally two best-fit solutions, one with|J| > |J'| and one with |J| < |J'|. I always fix g = 2.
On the upside the results from both programs are comparable in the isotropic version of the Hamiltonian (d = 0). However, when I liberate d, while Phi gives reasonable values of 1.6-1.7 cm-1, ES gives 2.86 cm-1 for the one solution and 6.5 cm-1 for the other:
J J' d
A -22.9 -19.5 2.86
B -20.9 -25.6 6.49
While 2.86 cm-1 could be considered borderline reasonable, 6.5 cm-1 is way too high.
So:
-somewhere in my code there is a small error
-there are some details in ES's implementation of calculating magn. susceptibility that affects systems with antisymmetric exchange.
The ES code is:
And my custom fitting function is:
The data for ES is here:
In the following post I will upload the the Phi files.
Thanks a lot!
UPDATE: On the above, I must also note that I have plotted the Zeeman diagrams for the same parameter sets and they are identical.
In particular, I fit an all-ferric triangle (S1 = S2 = S3 = 5/2) under an isosceles magnetic model, including antisymmetric exchange: H = -2J (S1S2 + S2S3) - 2J'S1S3 - 2d(S1xS2 + S2xS3 + S3xS1) + HZeeman. For this type of antiferromagnetic triangles, there are normally two best-fit solutions, one with|J| > |J'| and one with |J| < |J'|. I always fix g = 2.
On the upside the results from both programs are comparable in the isotropic version of the Hamiltonian (d = 0). However, when I liberate d, while Phi gives reasonable values of 1.6-1.7 cm-1, ES gives 2.86 cm-1 for the one solution and 6.5 cm-1 for the other:
J J' d
A -22.9 -19.5 2.86
B -20.9 -25.6 6.49
While 2.86 cm-1 could be considered borderline reasonable, 6.5 cm-1 is way too high.
So:
-somewhere in my code there is a small error
-there are some details in ES's implementation of calculating magn. susceptibility that affects systems with antisymmetric exchange.
The ES code is:
Code: Select all
clear all;
close all;
cm=100*clight/1e6; % Conversion constant from cm-1 to MHz
% [T,chiSI] = textread('data_easyspin_chiSI.dat','%f %f'); % Data is x vs T in SI units
[T,chi] = textread('data_easyspin_chi.dat','%f %f'); chiSI = chi*4*pi*1e-6; % Data is x vs T in cgs units
chiSIT=chiSI.*T;
[m,n]=size(chiSI); % size of experimental data array
Exp.Field = 1000; Exp.Temperature = T;
format long;
Ja = -20.2;
Jb = -25.3;
gxy = 2;
gz = 2;
Gx = 0; Gy = 0; Gz = 3.0;
Sys1.S=[5/2 5/2 5/2];
Sys1.g=[gxy gz; gxy gz; gxy gz];
Sys1.Ja = -2*cm*Ja;
Sys1.Jb = -2*cm*Jb;
Sys1.Gz = 3*cm*Gz;
Vary1.g = [0 0];
Vary1.Ja = 2*cm;
Vary1.Jb = 2*cm;
Vary1.Gz = 2*cm;
% Fitting arguments
FitOpt.OutArg = [2 2]; % to fit the magnetic susceptibility chizz (second output)
FitOpt.Method = 'simplex fcn';
FitOpt.Scaling = 'none';
% Opt.nKnots = 501; % I have tried with no Opt.nKnots, with Opt.nKnots = 101, 301 and 501
% % Get graph
% esfit('constrainJgd_isosceles_xT',chiSIT,Sys1,Vary1,Exp,[],FitOpt);
% Get ascii
tic
[BestSys,BestSpc] = esfit('constrainJgd_isosceles_xT',chiSIT,Sys1,Vary1,Exp,[],FitOpt);
toc
calcdata = [T(:) BestSpc(:)/(4*pi*1e-6)]; % Arrange calculated data as an ascii file with chi in emu
datafname = ['Ja=',num2str(BestSys.Ja/(-2*cm),'%.3f'),'_Jb=',num2str(BestSys.Jb/(-2*cm),'%.3f'),'_d=',num2str(BestSys.Gz/(2*cm),'%.3f'),'_bestfit.txt'];
save(datafname,'calcdata','-ascii'); % Save the calculated curve as an ascii file
% Calculate error
chicalc = BestSpc.'; % Transpose BestSpc matrix to subtract from chiSI
residuals = (chiSIT-chicalc).^2./chiSIT.^2;
ressum = sum(residuals,1); % sum along the column
rfactor = ressum/m %divide by number of experimental points
Code: Select all
function [x,y] = constrainJgd_isosceles_xT(Sys1,Exp,Opt);
Gx = 0; Gy = 0;
fullSys = Sys1;
fullSys.g = [Sys1.g];
fullSys.ee = [
Sys1.Ja Sys1.Gz -Gy; -Sys1.Gz Sys1.Ja Gx; Gy -Gx Sys1.Ja;
Sys1.Ja Sys1.Gz -Gy; -Sys1.Gz Sys1.Ja Gx; Gy -Gx Sys1.Ja;
Sys1.Jb -Sys1.Gz Gy; Sys1.Gz Sys1.Jb -Gx; -Gy Gx Sys1.Jb]; % Because S1xS3 = -S3xS1
[x,y] = curry(fullSys,Exp,Opt);
T = Exp.Temperature;
[x,y] = curry(fullSys,Exp,Opt);
y = y.';
y = y.*T;
y = y.';
returnThanks a lot!
UPDATE: On the above, I must also note that I have plotted the Zeeman diagrams for the same parameter sets and they are identical.