57Fe Mossbauer spectrum model of bcc FexCo1-x alloy

Fit function libraries that can be linked to and utilized in the MossWinn program.
Z. Klencsar
Posts: 10
Joined: Fri May 20, 2016 11:17 pm

57Fe Mossbauer spectrum model of bcc FexCo1-x alloy

Postby Z. Klencsar » Tue May 24, 2016 5:13 pm

library sub_feco86a;

{

Written by Zoltán Klencsár, z.klencsar@mosswinn.hu

The routine calculates the 57Fe Mossbauer spectrum model of bcc FexCo1-x alloy by
taking into account additive effects of the 8 nearest neighbor (n.n.) and 6 next-nearest-neighbor (n.n.n.) iron atoms
on the hyperfine magnetic field, isomer shift and quadrupole shift experienced by the central Fe atom.
The relative occurrence of a given number of iron atoms in the first and in the second coordination sphere of
a central iron atom is assumed to follow the binomial probability distribution function.

Various degrees of order/disorder can be modeled by setting the proportion of alpha sites taken by Fe ( P(alpha) )
independently from the total x atomic concentration of Fe (CFe). P(alpha)=CFe means full disorder (alpha and beta sites being equivalent), CFe=1/2 and P(alpha)=1 means perfect B2-type order (alpha sites all taken by Fe, beta sites all taken by Co). Though in the program CFe and P(alpha) can be set independently, not all combinations of CFe and P(alpha) are physically possible. For impossible combinations the code will calculate nothing so that to inform you about the problem.

The model is based on the corresponding theory put forward in the work

J.E. Frackowiak: Phys. Stat. Sol. (a) 87 (1985) 109.

The code can be compiled with 32 bit Delphi compilers to a DLL library that can be used together
with the MossWinn program (version 4.0Pre or later) to fit corresponding Mossbauer spectra.

The source code of this library is also available for download from
http://www.mosswinn.hu/dlls/sub_feco86a.dpr

A compiled version of this code is available for download from
http://www.mosswinn.hu/dlls/sub_feco86a.dll

Copy the compiled DLL into the folder ..\MossWinn 4.0\DLLs\ in order to be able to use it in MossWinn.

If you use this code to fit spectra and publish the results, please,
reference in your work the corresponding paper of

Z. Klencsár, P. Németh, Z. Sándor, T. Horváth, I.E. Sajó, S. Mészáros,
J. Mantilla, J.A.H. Coaquira, V.K. Garg, E. Kuzmann, Gy. Tolnai:
Structure and magnetism of Fe-Co alloy nanoparticles
Journal of Alloys and Compounds 674 (2016) 153-161.

http://dx.doi.org/10.1016/j.jallcom.2016.03.068

as well as the above paper of J.E. Frackowiak.

If you improve upon this library code, please, post the improvement here as answer.

}

uses
Math;

procedure DLL1_SET_INTERACTION_NAME(var s:ShortString;source:ShortString);export;
begin
if (source='57FE') then s:='Fe-Co BINOM 86a' else s:='Undefined';
end;

procedure DLL1_SET_DESCRIPTION(var s:ShortString;source:ShortString);export;
begin
DLL1_SET_INTERACTION_NAME(s,source);{For simplicity description is set equal to the interaction name.}
end;

procedure DLL1_SET_SUBSPECTRUM_INTERNAL_VALUES(var AN,PN,WN:byte;
var Name:ShortString;Source:ShortString;
p:pointer);export;
begin
AN:=2;{Number of amplitude type parameters}
PN:=11;{Number of position type parameters}
WN:=1;{Number of width type parameters}
if (source='57FE') then Name:='Fe-Co BINOM 86a (1)' else Name:='Undefined (1)';
end;


procedure DLL1_SET_PARAMETERS_INTERNAL_VALUES(n:longint;
Cmin,Cmax,Vmin,Vmax:double;
var Name:ShortString;
var Minimum_Value,
Maximum_Value,
Initial_Value:double;
Source:ShortString;
p:pointer);export;
begin

Case n of
1:begin
Name:='TOTAL AREA'; {Total spectral area of the component associated with bcc FexCo1-x.}
Minimum_Value:=0;
Maximum_Value:=(Vmax-Vmin)*(Cmax-Cmin)/4;
Initial_Value:=Maximum_Value/5;
end;

2:begin
Name:='A25 / A34';{Relative area of the 2-5 peaks of all the sextets with respect to that of the 3-4 peaks.}
Minimum_Value:=0;
Maximum_Value:=4;
Initial_Value:=2;
end;

3:begin
Name:='IS (14 Co)'; {Isomer shift of 57Fe with only Co as n.n. and n.n.n.}
Minimum_Value:=-0.1;
Maximum_Value:=0.1;
Initial_Value:=0.0;
end;

4:begin
Name:='DELTA IS / Fe (8nn)'; {Additive effect of one n.n. Fe to the isomer shift of the central Fe atom.}
Minimum_Value:=0.00;
Maximum_Value:=0.01;
Initial_Value:=0.0025;
end;

5:begin
Name:='DELTA IS / Fe (6nnn)'; {Additive effect of one n.n.n. Fe to the isomer shift of the central Fe atom.}
Minimum_Value:=0.00;
Maximum_Value:=0.01;
Initial_Value:=0.0025;
end;

6:begin
Name:='HMF (14 Co)'; {Hyperfine magnetic field of 57Fe with only Co as n.n. and n.n.n.}
Minimum_Value:=27;
Maximum_Value:=32;
Initial_Value:=29;
end;

7:begin
Name:='DELTA HMF / Fe (8nn)'; {Additive effect of one n.n. Fe to the hyperfine magnetic field of the central Fe atom.}
Minimum_Value:=0.0;
Maximum_Value:=1.0;
Initial_Value:=0.77;
end;

8:begin
Name:='DELTA HMF / Fe (6nnn)'; {Additive effect of one n.n.n. Fe to the hyperfine magnetic field of the central Fe atom.}
Minimum_Value:=0.0;
Maximum_Value:=1.0;
Initial_Value:=0.77;
end;

9:begin
Name:='QS (14 Co)'; {2*quadrupole shift of 57Fe with only Co as n.n. and n.n.n.}
Minimum_Value:=-0.1;
Maximum_Value:=0.1;
Initial_Value:=0.0;
end;


10:begin
Name:='DELTA QS / Fe (8nn)'; {Additive effect of one n.n. Fe to the (2*quadrupole shift) of the central Fe atom.}
Minimum_Value:=-0.02;
Maximum_Value:=0.02;
Initial_Value:=0.0;
end;

11:begin
Name:='DELTA QS / Fe (6nnn)'; {Additive effect of one n.n.n. Fe to the (2*quadrupole shift) of the central Fe atom.}
Minimum_Value:=-0.02;
Maximum_Value:=0.02;
Initial_Value:=0.0;
end;

12:begin
Name:='Fe concentration'; {x in FexCo1-x}
Minimum_Value:=0.01;
Maximum_Value:=1.00;
Initial_Value:=0.5;
end;

13:begin
Name:='P(alpha)'; {Proportion of alpha sites taken by Fe. Constrain it to be equal to "Fe concentration" in order to model full disorder.}
Minimum_Value:=0;
Maximum_Value:=1;
Initial_Value:=0.5;
end;

14:begin
Name:='LINE WIDTH';
Minimum_Value:=0.23;
Maximum_Value:=0.35;
Initial_Value:=0.26;
end;

end;

end;

procedure DLL1_SET_SUBSPECTRUM(n:longint;
pPARAMETERS,pDATA,pVELOCITY:pointer;
AreaMode:boolean;Source:ShortString;p:pointer);export;
CONST
twopi = pi+pi;

{Physical constants}
speed_of_light = 2.99792458E8;
elementary_charge = 1.60217653E-19;

gFe32 = -0.103267; {I=3/2 giromagnetic factor of 57Fe}
gFe12 = +0.18088; {I=1/2 giromagnetic factor of 57Fe}

nuclear_magneton = 5.05078343E-27; {Joule/Tesla}
hreduced = 1.05457168E-34; {Js/rad} {Planck constant divided by 2*pi}

EgammaFe = 14.4129; {keV}
Fe_Egamma = EgammaFe*1000.0*elementary_charge; {Joule}
mmpstoJoule57Fe = 0.001*Fe_Egamma/speed_of_light;

HMe = -gFe32*nuclear_magneton/mmpstoJoule57Fe;
HMg = -gFe12*nuclear_magneton/mmpstoJoule57Fe;

{Sextet line positions for 1 T hyperfine magnetic field, in mm/s}
PB1 = HMe*(-3/2)-HMg*(-1/2);
PB2 = HMe*(-1/2)-HMg*(-1/2);
PB3 = HMe*(1/2)-HMg*(-1/2);
PB4 = HMe*(-1/2)-HMg*(1/2);
PB5 = HMe*(1/2)-HMg*(1/2);
PB6 = HMe*(3/2)-HMg*(1/2);

TYPE
TARRAY = array[1..8192] of double;

TSiteRecord = Record
Weight:double;
ISH:double;
QS:double;
HMF:double;
LW:double;
End;

VAR
i,j:longint;

PARAMETERS,DATA,VELOCITY:^TARRAY;
Area,RELA:double;
FeSites:Array[0..8,0..6] of TSiteRecord; {0..8 Fe neighbors in the 1st and 0..6 Fe neighbors in the 2nd coordination sphere of Fe}

IS0,DELTAISnn8,DELTAISnnn6,
HMF0,DELTAHMFnn8,DELTAHMFnnn6,
QS0,DELTAQSnn,DELTAQSnnn,
LWcommon,cFe,cCo,palpha,pbeta:double;

Procedure AddSinglet(const A,V0,LW:double); {Add Lorentzian singlet with area A, center V0 and FWHM of LW.}
var
i:longint;
xa,x,y,v,LW2p4:double;
begin
xa:=A*LW/twopi;
LW2p4:=LW*LW/4;
For i:=1 to n do
begin
v:=VELOCITY^[i];
x:=(v-v0);
y:=x*x+LW2p4;
DATA^[i]:=DATA^[i]+xa/y;
end;
end;

{A is the full area of the sextet.}
Procedure AddSextet(const A,ISH,HMF,QS,LW:double);
var
QSp2:double;
LP1,LP2,LP3,LP4,LP5,LP6:double;
LA16,LA25,LA34:double;
begin
QSp2:=QS/2;
LP1:=ISH+QSp2+PB1*HMF;
LP2:=ISH-QSp2+PB2*HMF;
LP3:=ISH-QSp2+PB3*HMF;
LP4:=ISH-QSp2+PB4*HMF;
LP5:=ISH-QSp2+PB5*HMF;
LP6:=ISH+QSp2+PB6*HMF;

LA34:=A/(8+RELA+RELA);
LA16:=3*LA34;
LA25:=RELA*LA34;

AddSinglet(LA16,LP1,LW);AddSinglet(LA16,LP6,LW);
AddSinglet(LA25,LP2,LW);AddSinglet(LA25,LP5,LW);
AddSinglet(LA34,LP3,LW);AddSinglet(LA34,LP4,LW);
end;

Function WeightOfSite(const nn,nnn:longint):double;
var
ii,nCo:longint;
bfnn8,bfnnn6,xx,yy,q:double;
begin
{First we consider iron atoms at the alpha sites.}

{First coordination sphere.}

{Binomial factors (8|nn)}
case nn of
0:bfnn8:=1;
1:bfnn8:=8;
2:bfnn8:=28;
3:bfnn8:=56;
4:bfnn8:=70;
5:bfnn8:=56;
6:bfnn8:=28;
7:bfnn8:=8;
8:bfnn8:=1;
end;

nCo:=8-nn; {Number of Co as n.n.}

xx:=1.0;
if nn>0 then for ii:=1 to nn do xx:=xx*pbeta;

yy:=1.0;
if nCo>0 then for ii:=1 to nCo do yy:=yy*(1-pbeta);

Result:=bfnn8*xx*yy; {Intermediate result}

{Second coordination sphere.}

{Binomial factors (6|nnn)}
case nnn of
0:bfnnn6:=1;
1:bfnnn6:=6;
2:bfnnn6:=15;
3:bfnnn6:=20;
4:bfnnn6:=15;
5:bfnnn6:=6;
6:bfnnn6:=1;
end;

nCo:=6-nnn; {Number of Co as n.n.n.}

xx:=1.0;
if nnn>0 then for ii:=1 to nnn do xx:=xx*palpha;

yy:=1.0;
if nCo>0 then for ii:=1 to nCo do yy:=yy*(1-palpha);

Result:=(palpha/(cFe+cFe))*bfnnn6*xx*yy*Result; {First term of Eq.(3)'s w0(i,j) in J.E. Frackowiak: Phys. Stat. Sol. (a) 87 (1985) 109.}

{Now the second term: iron atoms at the beta sites.}

{First coordination sphere.}

nCo:=8-nn;
xx:=1.0;
if nn>0 then for ii:=1 to nn do xx:=xx*palpha;

yy:=1.0;
if nCo>0 then for ii:=1 to nCo do yy:=yy*(1-palpha);

q:=bfnn8*xx*yy; {Intermediate result}

{Second coordination sphere.}

nCo:=6-nnn;
xx:=1.0;
if nnn>0 then for ii:=1 to nnn do xx:=xx*pbeta;

yy:=1.0;
if nCo>0 then for ii:=1 to nCo do yy:=yy*(1-pbeta);

{Final result:}
Result:=Result+(pbeta/(cFe+cFe))*bfnnn6*xx*yy*q; {Final weight, Eq.(3)'s w0(i,j) in J.E. Frackowiak: Phys. Stat. Sol. (a) 87 (1985) 109.}
end;

begin
if NOT(source='57FE') then exit;

PARAMETERS:=pPARAMETERS;DATA:=pDATA;VELOCITY:=pVELOCITY;

cFe:=MAX(MIN(PARAMETERS^[12],1),0); {x in FexCo1-x}
palpha:=MAX(MIN(PARAMETERS^[13],1),0);{Proportion of alpha sites taken by Fe.}
pbeta:=cFe+cFe-palpha;{Proportion of beta sites taken by Fe.}
if ((pbeta>1) or (pbeta<0)) then exit;{Not valid combination of cFe and Palpha. Nothing is calculated.}

{Note that would fixing/setting the iron concentration not be our preference, it would be better to use palpha and pbeta as fit parameters,
and calculate cFe = (palpha+pbeta)/2 , which would be valid for all valid values of palpha and pbeta.}

cCo:=1.0-cFe;

Area:=MAX(PARAMETERS^[1],0);{Total spectrum area of the FexCo1-x component.}
RELA:=MIN(MAX(PARAMETERS^[2],0),4);{Relative area of the 2-5 peaks with respect to that of the 3-4 peaks for all sextets.}

IS0:=PARAMETERS^[3];
HMF0:=PARAMETERS^[6];
QS0:=PARAMETERS^[9];

DELTAISnn8:=PARAMETERS^[4];
DELTAISnnn6:=PARAMETERS^[5];

DELTAHMFnn8:=PARAMETERS^[7];
DELTAHMFnnn6:=PARAMETERS^[8];

DELTAQSnn:=PARAMETERS^[10];
DELTAQSnnn:=PARAMETERS^[11];

LWCommon:=PARAMETERS^[14];

for j:=0 to 6 do {nnn}
begin
for i:=0 to 8 do {nn}
begin
With FeSites[i,j] do {i nn and j nnn Fe neighbors}
begin
ISH:=IS0+i*DELTAISnn8+j*DELTAISnnn6;
HMF:=HMF0+i*DELTAHMFnn8+j*DELTAHMFnnn6;
QS:=QS0+i*DELTAQSnn+j*DELTAQSnnn;
LW:=LWCommon;
Weight:=WeightOfSite(i,j);
end;
end;
end;

{The sum of weights is equal to 1}
for j:=0 to 6 do {nnn}
begin
for i:=0 to 8 do {nn}
begin
With FeSites[i,j] do AddSextet(Weight*Area,ISH,HMF,QS,LW);{i nn and j nnn Fe neighbors}
end;
end;

end;

procedure DLL1_SET_AREA_OF_SUBSPECTRUM(var AREA:double;
var CanDO:boolean;
pPARAMETERS:pointer;
AreaMode:boolean;Source:ShortString;
p:pointer);export;
TYPE
TARRAY = array[1..8192] of double;
var
PARAMETERS:^TARRAY;
begin
CanDO:=True;

if (source='57FE') then
begin
PARAMETERS:=pPARAMETERS;
Area:=MAX(PARAMETERS^[1],0);
end
else Area:=0.0;
end;

exports

DLL1_SET_INTERACTION_NAME,
DLL1_SET_DESCRIPTION,
DLL1_SET_SUBSPECTRUM_INTERNAL_VALUES,
DLL1_SET_PARAMETERS_INTERNAL_VALUES,
DLL1_SET_SUBSPECTRUM,
DLL1_SET_AREA_OF_SUBSPECTRUM;

begin

end.