• comment on precision of H and L
• comment on precision of HH and LL
• comment on efficiency of numerical integration of wave functions of harmonic oscillator
• comment on efficiency of numerical integration with wave functions of Morse oscillator

Estimation of Franck-Condon factors with model wave functions

by G. J. Fee , J. W. Nibler and J. F. Ogilvie @ 2001 June 6

Centre for Experimental and Constructive Mathematics, Department of Mathematics, 8888 University Drive,
Simon Fraser University, Burnaby, BC V5A 1S6 Canada, gjfee@cecm.sfu.ca, ogilvie@cecm.sfu.ca

Department of Chemistry, Oregon State University, 101 Gilbert Hall, Corvallis, OR 97331-4003 USA, joseph.nibler@orst.edu

The integrated intensity of each band in a system corresponding to transitions between two electronic states of a diatomic molecule is proportional to the square of a matrix element that one can form according to wave mechanics by integration of electronic and vibrational wave functions over electronic and nuclear coordinates. This matrix element can be taken to be a product of a factor involving electronic wave functions and a factor that is the square of an integral of wave functions of vibrational states; this squared integral is called a Franck-Condon factor, and this worksheet enables estimation of this quantity using vibrational wave functions of either an harmonic oscillator or a Morse oscillator.

> restart:

Here are three procedures useful in estimating Franck-Condon factors of diatomic molecules. Of these HH to produce Hermite's polynomials is a replacement for procedure H in Maple 's library orthopoly , and LL to produce Laguerre's polynomials is an analogous replacement for procedure L therein. For large values of vibrational quantum number these replacements are necessary because of errors of roundoff with use of hardware precision at 16 digits and because type checking in H is not supported when running inside evalhf .

>

comment on precision of H and L

For example, with Hermite's polynomial of order 70,

> with(orthopoly):

> H70 := sort(H(70,x));

of which the largest coefficient is,

> maxnorm(H70);

containing this many digits,

> length(%);

evaluation of that polynomial at as a decimal number yields a value

> eval(H70, x=9/2);

> evalf(%,16);

in contrast with the same value with approximate arithmetic at the same value of x .

> Digits := 16:

> eval(H70, x=4.5);

> (%%%-%)/%%%;

The error at this value of x for Hermite's polynomial of order 70 is thus of order 10 per cent.
Likewise for Laguerre's polynomial of order 70,

> L70 := sort(L(70,x));

the largest coefficient is only

> maxnorm(L70);

> evalf(%);

Evaluation of this polynomial at yields, with exact arithmetic,

> eval(L70,x=11/2);

> evalf(%,16);

and with approximate arithmetic,

> Digits := 16:

> eval(L70, x=5.5);

> (%%%-%)/%%%;

In this case the error is of order 10 per cent.

>

Another procedure numintn performs numerical integration more rapidly but less accurately than procedures provided directly in Maple .

> HH := proc(n, X) # This function is intended for only numeric values of X.
local x,i,s0,s1,s2;
x := evalf(X);
if 1 < n then
s0 := 1.;
s1 := 2.*x;
for i from 2 to n do
s2 := 2.*x*s1 - 2*(i-1)*s0;
s0 := s1;
s1 := s2;
end do;
s2;
elif n=1 then 2*x;
elif n=0 then 1.;
else ERROR(`n must be a non-negative integer`);
end if;
end:
LL := proc(n, A, R) # generalised Laguerre polynomials
local a,x,i,s0,s1,s2; # This procedure is intended for only numeric values of R.
a := evalf(A);
x := evalf(R);
if 1<n then
s0 := 1.;
s1 := 1 + a - x;
for i from 2 to n do
s2 := ((2.*i+a-1-x)*s1 - (i+a-1)*s0)/i;
s0 := s1;
s1 := s2;
end do;
s2;
elif n=1 then 1 + a - x;
elif n=0 then 1.;
else ERROR(`n must be a non-negative integer`);
end if;
end:
numintn := proc(f,a,b,n) # procedure for quick and accurate numerical integration
# fitting with polynomial of degree 7 at each subinterval
# f = function, a = lower limit, b = upper limit, n = number of subintervals
# evaluation at 8*n points
local ba,w,s,aw,i;
ba := evalf(b-a);
w := ba/n;
s := 0;
aw := a-w;
for i to n do
aw := aw + w;
s := s + 295627*evalf(f(aw+0.0625*w))
+ 71329*evalf(f(aw+0.1875*w))
+ 471771*evalf(f(aw+0.3125*w))
+ 128953*evalf(f(aw+0.4375*w))
+ 128953*evalf(f(aw+0.5625*w))
+ 471771*evalf(f(aw+0.6875*w))
+ 71329*evalf(f(aw+0.8125*w))
+ 295627*evalf(f(aw+0.9375*w));
end do;
s*w/1935360;
end proc:

>

comment on precision of HH and LL

> Digits := 16;

> ex1 := H(70, 9/2);

> ap1 := HH(70, 4.5);

> er1 := (ap1 - ex1)/ex1;

> evalf(ex1);

> ex2 := L(70, 93, 11/2);

> ap2 := LL(70, 93, 5.5);

> evalf(ex2);

> er2 := (ap2 - ex2)/ex2;

With precision 16 digits, the relative error is of order rather than with H and L.

>

These fundamental physical constants enter the calculations. SI units are used consistently throughout, with wave number hence as m , and distances are taken in units m/ .

> h := 6.62606876e-34: # fundamental physical constants in SI units
c := 299792458:
N[A] := 6.02214199e23:
k[B] := 1.3806503e-23:
cB := evalf(2*Pi/10^10*sqrt(2*c*mu/h)): # conversion factor

To illustrate this approach we employ spectral data of I , namely harmonic vibrational parameter , primary rotational parameter and electronic binding energy , of each electronic state. With these data for I in states X and B, these procedures are useful up to " = 70 and ' = 69.

> mu := (0.12690446842/2)/N[A]; # properties of diiodine I_2
we1 := 21450.2; # in electronic ground state X
Be1 := 3.7363;
De1 := 1254700;
Req1 := 1/cB/sqrt(Be1);
alpha1 := evalf((2*Pi)^2*mu*we1*c/h/10^20):
d1 := 2*De1/we1:
beta1 := cB*Req1*sqrt(we1/2):
alphaM1 := beta1/sqrt(d1):

> mu := (0.12690446842/2)/N[A]; # properties of diiodine I_2
we2 := 12569.7; # in electronically excited state B
Be2 := 2.9039;
De2 := 439100;
Req2 := 1/cB/sqrt(Be2);
alpha2 := evalf((2*Pi)^2*mu*we2*c/h/10^20):
d2 := 2*De2/we2:
beta2 := cB*Req2*sqrt(we2/2):
alphaM2 := beta2/sqrt(d2):

This procedure FCFho calculates a Franck-Condon factor with wave functions of an harmonic oscillator; in the last executable instruction of the procedure, in which appears the name of procedure numintn , the range of integration as lower and upper limits in a sequence, and the number of steps of integration -- here 2, 4, 8*24 -- must be set aptly for each molecular system; the number of steps of integration is actually eight times the number specified in that statement.

> FCFho := proc(v1,v2) # harmonic oscillator
global HO1,HO2,HO1HO2;
HO1 := subs(v1_=v1, R -> 1/sqrt(2^v1_*v1_!)*(alpha1/Pi)^(1/4)*exp(-alpha1*((R-Req1)^2)/2)
*HH(v1_,sqrt(alpha1)*(R-Req1)));
HO2 := subs(v2_=v2, R -> 1/sqrt(2^v2_*v2_!)*(alpha2/Pi)^(1/4)*exp(-alpha2*((R-Req2)^2)/2)
*HH(v2_,sqrt(alpha2)*(R-Req2)));
HO1HO2:= x -> HO1(x)*HO2(x):
evalhf(numintn(HO1HO2,2,4,24)^2);
end proc:

The Franck-Condon factor for a transition between two states of I , X v " = 0 and B ' = 15, calculated with wave functions of an harmonic oscillator is

> FCFho(0,15);

>

comment on efficiency of numerical integration of wave functions of harmonic oscillator

> st := time(): ans1 := FCFho(2,70): ans1time := time() - st; ans1;

If evalhf be replaced with evalf for the numerical integration, the duration of calculation might increase a factor 50 or more. For example,

> FCFhoslow := proc(v1,v2)
global HO1,HO2,HO1HO2;
HO1 := subs(v1_=v1, R -> 1/sqrt(2^v1_*v1_!)*(alpha1/Pi)^(1/4)*exp(-alpha1*((R-Req1)^2)/2)
*HH(v1_,sqrt(alpha1)*(R-Req1)));
HO2 := subs(v2_=v2, R -> 1/sqrt(2^v2_*v2_!)*(alpha2/Pi)^(1/4)*exp(-alpha2*((R-Req2)^2)/2)
*HH(v2_,sqrt(alpha2)*(R-Req2)));
HO1HO2:= x -> HO1(x)*HO2(x):
evalf(numintn(HO1HO2,2,4,24)^2);
end proc:

> Digits := 16;
st := time(): ans2 := FCFhoslow(2,70): ans2time := time() - st; ans2;

> timeratio := ans2time/ans1time;

>

These procedures calculate a Franck-Condon factor with wave functions of a Morse oscillator; in the last executable instruction of the procedure, in which appears the name of procedure numintn , the range of numerical integration -- here 2 .. 4, and number of steps of integration -- here , must be set aptly for each molecular system.

> NMO1 := proc(v1) # Morse oscillator
local j,s,t,u;
u := 2*d1 - 2*v1 - 1;
t := GAMMA(u);
s := t:
for j from 1 to v1 do
t := t*u/j;
u := u + 1;
s := s + t;
end do;
1/sqrt(s*Req1/alphaM1);
end proc:
NMO2 := proc(v2)
local j,s,t,u;
u := 2*d2 - 2*v2 - 1;
t := GAMMA(u);
s := t:
for j from 1 to v2 do
t := t*u/j;
u := u + 1;
s := s + t;
end do;
1/sqrt(s*Req2/alphaM2);
end proc:
a1 := v1 -> 2*d1 - 2*v1 - 1;
a2 := v2 -> 2*d2 - 2*v2 - 1:
y1 := R -> 2*d1*exp(-alphaM1*(R/Req1-1)):
y2 := R -> 2*d2*exp(-alphaM2*(R/Req2-1)):
lnNMO1_ := hfarray(0..10):
lnNMO2_ := hfarray(0..70):
FCFMo := proc(v1,v2)
global MO1,MO2,MO1MO2,lnNMO1_,lnNMO2_;
lnNMO1_[v1] := evalf(ln(NMO1(v1))):
lnNMO2_[v2] := evalf(ln(NMO2(v2))):
MO1 := subs(v1_=v1, R -> LL(v1_,a1(v1_),y1(R))*
exp(lnNMO1_[v1_]+a1(v1_)/2*ln(y1(R))-y1(R)/2)):
MO2 := subs(v2_=v2, R -> LL(v2_,a2(v2_),y2(R))*
exp(lnNMO2_[v2_]+a2(v2_)/2*ln(y2(R))-y2(R)/2)):
MO1MO2:= x -> MO1(x)*MO2(x):
evalhf(numintn(MO1MO2,2,5,iquo(5*v2+100,6))^2);
end proc:

The Franck-Condon factor for a transition between two states of I , X v " = 0 and B ' = 35, calculated with wave functions of a Morse oscillator is

> FCFMo(0,35);

>

comment on efficiency of numerical integration with wave functions of Morse oscillator

> st := time(): ans1 := FCFMo(2,65): ans1time := time() - st; ans1;

If evalhf be replaced with evalf for the numerical integration, the duration of calculation might increase a factor 40 or more. For example,

> FCFMoslow := proc(v1,v2)
global MO1,MO2,MO1MO2,lnNMO1_,lnNMO2_;
lnNMO1_[v1] := evalf(ln(NMO1(v1))):
lnNMO2_[v2] := evalf(ln(NMO2(v2))):
MO1 := subs(v1_=v1, R -> LL(v1_,a1(v1_),y1(R))*
exp(lnNMO1_[v1_]+a1(v1_)/2*ln(y1(R))-y1(R)/2)):
MO2 := subs(v2_=v2, R -> LL(v2_,a2(v2_),y2(R))*
exp(lnNMO2_[v2_]+a2(v2_)/2*ln(y2(R))-y2(R)/2)):
MO1MO2:= x -> MO1(x)*MO2(x):
evalf(numintn(MO1MO2,2,5,iquo(5*v2+100,6))^2);
end proc:

> st := time(): ans2 := FCFMoslow(2,65): ans2time := time() - st; ans2;

> timeratio := ans2time/ans1time;

>

>