Symbolic Programmes for Analysis of Vibrational and Rotational Spectra of Diatomic Molecules

This work sheet in Maple contains procedures in four suites described

by F. M. Fernandez and J. F. Ogilvie

in journal MapleTech , volume 5, 1998, No. 1, pages 42 - 46. These procedures for symbolic computation in Maple (version V, release 5.1) allow one to form the principal expressions required to undertake a comprehensive analysis of vibration-rotational spectra of diatomic molecules in an electronic state singlet-SIGMA-+ according to conventional analytic treatment in the spirit of Dunham (cf. J. F. Ogilvie, The Vibrational and Rotational Spectrometry of Diatomic Molecules , Academic Press, London, 1998). To relate to measurements of frequencies or wavenumbers of spectral transitions, two procedures form, in the first suite, expressions Y_kl for spectral term coefficients associated with mechanical effects, i.e. vibrational and rotational motions of atomic centres about the centre of molecular mass, and, in the second suite, expressions Z_kl for further spectral term coefficients associated with extra-mechanical effects, i.e. in relation to failure of electrons to follow perfectly one or other nucleus. To relate to measurements of intensities of spectral lines, two additional procedures form, in the third suite, vibrational matrix elements <v|x^l|v'> and, in the fourth suite, the Herman-Wallis factor F_0^v' that takes into account vibration-rotational interactions that affect relative intensities within bands. Use of these suites is described specifically in the paper in MapleTech mentioned above; their applications are described at length in the book specified above.

1 Term coefficients Y_kl for mechanical effects

The first suite DUNHAM of procedures below generates Dunham's term coefficients Y[k,l] for vibration-rotational energies of diatomic molecules in electronic state singlet-SIGMA+ with the method HPM (TP-hypervirial-Hellmann-Feynman perturbation theory); the quantity JB is considered to be independent of perturbation parameter lambda;

definitions: lambda^2=2*Be/We=gamma=g, vh=v+1/2, JB=JJ*g/2, JJ=J(J+1).
To use procedure tp1 to generate term coefficients Y[k,l,m], in which k is the power of (v+1/2), l is the power of [J(J+1)] and m is the order of the contribution (m=1 for the leading contribution to Y[k,l], m=2 for the first correction etc.), invoke "tp1(p1):" in which the argument "p1" is the required order of perturbation theory. This procedure should work with values of p1 up to at least p1=14. All expressions Y[k,l,m] according to the selected order of perturbation theory are automatically generated, and are invoked simply on typing "Y[k,l,m];" with appropriate values of k, l and m. For p1=14, the range of Y[k,l,1] is from k=8,l=0 with a_1^14 ... a_14 to k=0,l=8 with a_1^6 ... a_6.
To use procedure tp2 to generate term coefficients Y[k,l,m], with the same meaning of symbols as for tp1, invoke "tp2(p2):" in which the argument "p2" is the required order of perturbation theory. This procedure should work with values of p2 up to at least p2=20, but the duration of execution of tp2 for a given order is about three times that of tp1. For p2=20, the range of Y[k,l,1] is from k=11,l=0 with a_1^20 ... a_20 to k=0,l=11 with a_1^9 ... a_9.
By default expressions Y[k,l,m] appear in terms of primary rotational parameter Be, reciprocal W1 of primary vibrational parameter (harmonic vibrational wavenumber) We, and coefficients a_j in Dunham's function V(x) for potential energy with x=(R-Re)/Re. If these expressions are preferred in terms of coefficients b_j of V(y) with y=(R-Re)/R or c_j of V(z) with z=2(R-Re)/(R+Re) instead of a_j, then after the expressions Y[k,l,m] are formed with tp1 or tp2 call procedure "converta(l):" with argument l as either b or c as required. To ensure that resulting expressions Y[k,l,m] in terms of b_j or c_j appear in an attractive form, type an instruction "sort(simplify(Y[k,l,m]));" to display the desired expression.

> restart: tp1:=proc(p1) # perturbation theory of order p1
local i,s,N,Q,p: global En,JJ,vh:
if type(p1,odd) then p:=p1-1 else p:=p1 fi:
Q[0,0]:=1: En[0]:=vh:
for s from 0 to p do
if s> 0 then
Q[0,s]:=0:
if type(s,odd) then En[s]:=0 else
En[s]:=expand(1/s*add(i*(a.i/2*Q[i+2,s-i]
+JJ*(-1)^i*(i+1)*Q[i,s-i]),i=1..s)):
fi:
fi:
for N from 0 to p-s+1 do
if type(N+1+s,odd) then Q[N+1,s]:=0 else
Q[N+1,s]:=expand(1/(N+1)*(trunca(N-3)*N*(N-1)*(N-2)*Q[N-3,s]/4
+2*N*En[0]*trunca(N-1)*Q[N-1,s]
+add(2*N*En[i]*trunca(N-1)*Q[N-1,s-i]
-(N+i/2+1)*a.i*Q[N+i+1,s-i]
-(-1)^i*(i+1)*(2*N+i)*JJ*Q[N+i-1,s-i],i=1..s))):
fi:
od:
od:
En[0]:=En[0]+JJ;
for i from 0 to p by 2 do xtrct(i) od; gc() end:

> trunca:=proc(n) if n<0 then 0 else 1 fi: end:

> xtrct:=proc(n) local k,l,p; global En,Y,JJ,vh,Be,W1;
for k from 0 to n/2+1 do p:=subs(vh=0,En[n]); En[n]:=expand((En[n]-p)/vh);
for l from 0 to n/2+1 do Y[k,l,3/2+n/4-k/2-l/2]:=combine(
(Be*W1)^(l+n/2)/W1*sort(2^(n/2)*coeff(p,JJ,l))) od
od; En[n]:=0
end:

> converta:=proc(l) local j,k; global b,c;
if l=c then
for k from 1 to 22 do a.k:=(-1)^k*(k+1)/2^k+mysum('(-1)^j*(k+1)!*
c.(k-j)/(2^j*j!*(k+1-j)!)', 'j'=0..k-1) od
elif l=b then
for k from 1 to 22 do a.k:=(-1)^k*(k+1)+mysum('(-1)^j*(k+1)!*
b.(k-j)/(j!*(k+1-j)!)','j'=0..k-1) od
fi end:

> tp2:=proc(p2) # perturbation theory of order p2
local i,k,l,k1,l1,s,N,p: global e, Q:
if type(p2,odd) then p:=p2-1 else p:=p2 fi:
for s from 0 to p do
for N from 0 to p-s+1 do
for k from 0 to p do
for l from 0 to p do
Q[N,k,l,s]:=0: e[k,l,s]:=0:
od:
od:
od:
od:
Q[0,0,0,0]:=1: e[0,0,0]:=1:
for s from 0 to p do
for k from 0 to s/2+1 do
for l from 0 to s/2+1 do
if s> 0 then
if type(s,odd) then e[k,l,s]:=0 else
e[k,l,s]:=simplify(1/s*add(i*(a.i/2*trunca(s/2+1-k-l)
*Q[i+2,k,l,s-i]
+(-1)^i*(i+1)*trunca(l-1)*
trunca(s/2-k-l+1)*Q[i,k,l-1,s-i]),i=1..s)):
fi:
fi:
od:
od:
for N from 0 to p-s+1 do
for k from 0 to (N+1+s)/2 do
for l from 0 to (N+1+s)/2 do
if type(N+1+s,odd) then Q[N+1,k,l,s]:=0 else
Q[N+1,k,l,s]:=simplify(1/(N+1)*(trunca(N-3)*N*(N-1)*(N-2)
*Q[N-3,k,l,s]/4
+2*N*trunca(N-1)*trunca(k-1)*Q[N-1,k-1,l,s]
+add(2*N*add(add(e[k1,l1,i]*
trunca(N-1)*Q[N-1,k-k1,l-l1,s-i],
l1=0..l),k1=0..k)
-(N+i/2+1)*a.i*Q[N+i+1,k,l,s-i]
-(-1)^i*(i+1)*(2*N+i)*trunca(l-1)*Q[N+i-1,k,l-1,s-i],i=1..s))):
fi
od
od
od
od: e[0,1,0]:=1: e[1,0,0]:=1: extrct(p); gc()
end:

> extrct:=proc(n) local k,l,m; global e,Y,Be,W1;
for k from 0 to n/2+1 do
for l from 0 to n/2+1 do
for m from 0 to n do
Y[k,l,3/2+m/4-k/2-l/2]:=combine((Be*W1)^(l+m/2)/W1*
sort(2^(m/2)*e[k,l,m]))
od
od
od
end:

> tp1(12): Y[2,1,1]; # for instance

2 Term coefficients Z_kl for extra-mechanical effects

Procedures in this suite ADNONAD allow one to generate additional term coefficients Z[k,l] to take account of adiabatic and nonadiabatic contributions to vibration-rotational energies of diatomic molecules in electronic state singlet-SIGMA+
To use procedure tp3 to generate term coefficients Z[k,l], in which k is the power of (v+1/2) and l is the power of [J(J+1)], invoke "tp3(P3):", in which the argument P3 is the required order of perturbation theory. This procedure should work with values of P1 up to at least 12. All expressions Z[k,l] are automatically generated, and are displayed simply on typing "Z[k,l];" with appropriate values of k and l. For P1=12, the range of Z[k,l] is k+l<=6.
To convert coefficients a_j, K_j, q_j and d_j in V(x) etc. to c_j, u_j, t_j and s_j in V(z) etc. invoke "convertacstu():" after running tp3; then display an expression with "sort(simplify(Z[k,l]));".

> restart: tp3:=proc(P3) # perturbation theory of order P3
local l,n,p,i,k,p1,la,P,s1,s2,s3,JB,ad1,ad2,ad3;
global lista,Be,DE,En,JJ,We,en,g,Q,vh;
if type(P3,odd) then P:=P3-1 else P:=P3 fi;
lista:=[g,Be,K.(1..12),q.(0..8),d.(0..8),a.(1..14)];
Q[0,0]:=1; en[0]:=vh;
for p from 0 to P do
if p>0 then
Q[0,p]:=0:
if type(p,odd) then en[p]:=0 else
en[p]:=expand(1/p*add(i*(a.i/2*Q[i+2,p-i]
+JB*(-1)^i*(i+1)*Q[i,p-i]),i=1..p)):
fi:
fi;
for n from 0 to P-p+1 do
if type(n+1+p,odd) then Q[n+1,p]:=0 else
Q[n+1,p]:=expand(1/(n+1)*(trunca(n-3)*n*(n-1)*(n-2)*Q[n-3,p]/4
+2*n*en[0]*trunca(n-1)*Q[n-1,p]
+add(2*n*en[i]*trunca(n-1)*Q[n-1,p-i]
-(n+i/2+1)*a.i*Q[n+i+1,p-i]
-(-1)^i*(i+1)*(2*n+i)*JB*Q[n+i-1,p-i],i=1..p))):
fi: od; od; JB:=JJ/2*g:
# Form corrections adiabatic and nonadiabatic.
la:=g^(1/2); s1:=0; s2:=0; s3:=0;
for k from 0 to P do
# correction nonadiabatic vibrational
ad1[k]:=expand(k*(k-1)/4*add(Q[k-2,p1]*la^(p1+k),p1=0..P-k)
+1/2/(k+1)*(1/2*add((i+2)*a.i*add(Q[i+k+2,p1]*la^(p1+k+i),
p1=0..P-k-i),i=1..P-k)+JB*add((-1)^i*i*(i+1)
*add(Q[i+k,p1]*la^(p1+k+i),p1=0..P-k-i),
i=1..P-k)+add(Q[k+2,p1]*la^(k+p1),p1=0..P-k)));
# correction nonadiabatic rotational
ad2[k]:=expand(JB*add((-1)^i*(i+1)*la^(i+k)
*add(Q[i+k,l]*la^l,l=0..P-k-i),i=0..P-k));
# correction adiabatic multiplied by We
ad3[k]:=add(Q[k,i]*la^(i+k),i=0..P-k);
s1:=s1+ad1[k]*d.k; s2:=s2+ad2[k]*q.k; s3:=s3+ad3[k]*K.k;
od;
# sum of contributions adiabatic and nonadiabatic in the energy
DE:=simplify(We*(s1+s2)+s3); xtrct(P) end:

> trunca:=proc(n) if n<0 then 0 else 1 fi: end: # utility to truncate terms

> xtrct:=proc(P) # extraction of Z[k,l] from energy
local k,l,p,s,s1,z; global lista,Z,DE,vh,Be,We,JJ,g;
for k from 0 to P/2 do s:=subs(vh=0,DE); DE:=expand((DE-s)/vh);
for l from 0 to P/2 do s1:=subs(JJ=0,s); s:=expand((s-s1)/JJ);
for p from 0 to P do z[k,l,p]:=
sort(collect(expand(subs(We=2*Be/g,coeff(s1,g,p)*g^p)),lista))
od od od;
for k from 0 to P/2 do
for l from 0 to P/2 do if k=0 and l=0 then Z[0,0]:=z[0,0,0]
elif z[k,l,k+2*l-1]<>0 and z[k,l,k+2*l]<>0
then Z[k,l]:=z[k,l,k+2*l-1]+z[k,l,k+2*l] fi od od end:

> convertacstu:=proc() # to convert from a_j to c_j, K_j to u_j etc.
local j,l1,l2,l3,l4; global a,K,q,d,c,s,t,u,q0,t0,b0,s0;
for l1 from 1 to 12 do a.l1:=(-1)^l1*(l1+1)/2^l1+add((-1)^j*(l1+1)!
*c.(l1-j)/(2^j*j!*(l1+1-j)!), j=0..l1-1) od:
for l2 from 1 to 8 do K.l2:=add((-1)^(l2-j)*2^(j-l2)*(l2-1)!
*u.j/((l2-j)!*(j-1)!),j=1..l2) od;
for l3 from 1 to 8 do q.l3:=add((-1)^(l3-j)*2^(j-l3)*(l3-1)!
*t.j/((l3-j)!*(j-1)!),j=1..l3) od; q0:=t0;
for l4 from 1 to 8 do d.l4:=add((-1)^(l4-j)*2^(j-l4)*(l4-1)!
*s.j/((l4-j)!*(j-1)!),j=1..l4) od; b0:=s0; end:

> tp3(6): Z[2,1]; # for instance

3 Vibrational matrix elements

Procedures in this suite VIBMATEL allow one to generate vibrational matrix elements <v'|x^l|v> of the coordinate x to various non-negative powers between vibrational states v and v'.
To use this suite, specify the required order P of perturbation theory; then specify values of vibrational quantum number for vibrational states v and v' as arguments of "state". For instance, to generate matrix elements of x to various powers between vibrational states v=0 and v'=1 up to a_6, type "P:=6: state(0): state(1):". Particular matrix elements, in the array mel[v,l,v'], are formed on typing "me(v,l,v'):". If matrix elements <v'|x^l|v> are required to
be expressed in terms of coefficients c_j of function V(z) for potential energy, after matrix elements <v'|x^l|v> are generated they can be so converted by invoking "convertac():". To ensure that resulting expressions appear in attractive form, type "series(sort(simplify(mel[v',l,v])),g);" to display a specific expression.

> restart: P:=6: # The default value for greatest order of perturbation is 6.

> state:=proc(n) # corrections to energy e[n,p] and wavefunction c[j,n,p]
local j,k,p,s,e: global c,P,a:
c[n,n,0]:=1:
for p from 0 to P do
if p>0 then
e[n,p]:=expand(add(a.s/2*add(Q(n,s+2,k)*c[k,n,p-s]*cut(k,n,p-s),
k=max(n-s-2,0)..n+s+2),s=1..p)- add(e[n,s]*c[n,n,p-s],s=1..p-1))
fi:
for j from 0 to n+3*p do
if p>0 and j<>n then
c[j,n,p]:=expand(1/(j-n)*(add(e[n,s]*c[j,n,p-s]*cut(j,n,p-s),s=1..p)
- add(a.s/2*add(c[k,n,p-s]*Q(j,s+2,k)*cut(k,n,p-s)
,k=max(j-s-2,0)..j+s+2),s=1..p))) fi: od:
if p>0 then c[n,n,p]:=-expand(1/2*add(add(c[k,n,s]*c[k,n,p-s],
k=0..n+3*min(s,p-s)),s=1..p-1)) fi: od: end:

> # matrix elements <i|q^n|j> in a basis set of an harmonic oscillator
Q:=proc(i,n,j) option remember:
if i<0 or j<0 then 0
elif n=0 then delta(i,j) else
expand(sqrt(i/2)*Q(i-1,n-1,j)+sqrt((i+1)/2)*Q(i+1,n-1,j))
fi: end:

> # cutoff function to achieve c[j,n,p]=0 if |n-j|>3p
cut:=proc(n,j,p) if abs(n-j)<=3*p then 1 else 0 fi: end:

> delta:=proc(i,j) if i=j then 1 else 0 fi end: # Kronecker delta

> me:=proc(j,k,n) global g,l,p,P,mel; # matrix elements <m|q^k|n>
l:=g^(1/2); mel[j,k,n]:=sum('mep(j,k,n,p)*l^p','p'=0..P):
mel[n,k,j]:=mel[j,k,n]: end:

> # contributions to matrix elements <m|q^k|n>
mep:=proc(m,k,n,p) local s,i,j: global c,l:
l^k*sort(combine(expand(add(add(add(c[i,m,s]*c[j,n,p-s]*cut(i,m,s)
*cut(j,n,p-s)*Q(i,k,j),j=0..n+3*(p-s)),i=0..m+3*s),s=0..p)),
radical)): end:

> convertac:=proc() local j,k; global c; # to convert from a_j to c_j of V(z)
for k from 1 to 12 do a.k:=(-1)^k*(k+1)/2^k+add((-1)^j*(k+1)!*
c.(k-j)/(2^j*j!*(k+1-j)!), j=0..k-1) od end:

> P:=6: state(0): state(1): me(0,1,1); # for instance

4 Coefficients in the Herman-Wallis factor

Procedures in this suite HERMWALL allow one to generate quantities used to form the Herman-Wallis factor F[0,v'](iota) for electric-dipolar transitions with J'-J= +1 or -1.
To use this suite to form the Herman-Wallis factor for the band from v=0 to v=v', invoke "hw0(n):" in which n is the numerical value of v'; then the results of the calculation are alpha(0,v') and beta(0,v'), displayed on typing "aln;" and "ben;". For instance, for the band v=2 -- v=0, input "hw0(2):", then "al2;" and "be2;". F[0,v'](iota) = 1 + 2 alv' iota / <0|M(x)|v'> + (alv'^2 + 2 bev') iota^2 / <0|M(x)|v'>^2

> restart: hw0:=proc(v)
local i,lista,dip: global g,l,JJ,m,P:
P:=v+4: # Greatest perturbation order
# m=iota=1/2 [J'(J'+1)-J(J+1)], J'=J+1,J-1 for electric-dipolar transitions
# Matrix elements, energies etc. are reduced to terms in m^0, m^1, and m^2.
lista:=[l,M.(0..v)]:
JJ:=m*(m-1): state(0):
JJ:=m*(m+1): state(v):
# first two terms of the matrix element <0J|M(x)|v'J'>
dip:=l^(v+2)*add(M.i*mep(0,i,v,v+2-i),i=0..v+2)
+l^(v+4)*add(M.i*mep(0,i,v,v+4-i),i=0..v+4):
dip:=collect(dip,m):
dip:=dip-coeff(dip,m,0):
# Form alpha(0,v) and beta(0,v).
al.v:=sort(collect(combine(coeff(dip,m,1)/l^v,radical),lista)): al.0:=0:
be.v:=sort(collect(combine(coeff(dip,m,2)/l^v,radical),lista)):
al.v:=subs(l=g^(1/2),al.v):
be.v:=subs(l=g^(1/2),be.v): end:

> delta:=proc(i,j) if i=j then 1 else 0 fi end: #Kronecker delta

> # matrix elements <i|q^n|j> in a basis set for an harmonic oscillator
Q:=proc(i,n,j) option remember:
if i<0 or j<0 or n<0 then 0
elif n=0 then delta(i,j) else
expand(sqrt(i/2)*Q(i-1,n-1,j)+sqrt((i+1)/2)*Q(i+1,n-1,j)) fi: end:

> # truncation function to achieve c[j,n,p]=0 if |n-j|>3p
cut:=proc(n,j,p) if abs(n-j)<=3*p then 1 else 0 fi: end:

> # truncation to shift the the rotational part two orders
cut2:=proc(s) if s>2 then 1 else 0 fi end:

> # corrections to energy e[n,p] and wavefunction c[j,n,p]
state:=proc(n)
local j,k,p,s,e: global JJ,P,m,c:
c[n,n,0]:=1:
for p from 0 to P do
if p>0 then
e[n,p]:=expand(add(a.s/2*add(Q(n,s+2,k)*c[k,n,p-s]*cut(k,n,p-s),
k=max(n-s-2,0)..n+s+2) + (-1)^s*(s-1)/2*JJ*cut2(s)*add(Q(n,s-2,k)
*c[k,n,p-s]*cut(k,n,p-s),k=max(n-s,0)..n+s),s=1..p)
-add(e[n,s]*c[n,n,p-s],s=1..p-1)) fi:
e[n,p]:=collect(e[n,p],m):
e[n,p]:=coeff(e[n,p],m,0)+m*coeff(e[n,p],m,1)+m^2*coeff(e[n,p],m,2):
for j from 0 to n+3*p do
if p>0 and j<>n then
c[j,n,p]:=expand(1/(j-n)*(add(e[n,s]*c[j,n,p-s]
*cut(j,n,p-s),s=1..p) - add(a.s/2*add(c[k,n,p-s]
*Q(j,s+2,k)*cut(k,n,p-s),k=max(j-s-2,0)..j+s+2)
+(-1)^s*(s-1)/2*JJ*cut2(s)*add(c[k,n,p-s]*Q(j,s-2,k)
*cut(k,n,p-s),k=max(j-s,0)..j+s),s=1..p))) fi:
c[j,n,p]:=collect(c[j,n,p],m):
c[j,n,p]:=coeff(c[j,n,p],m,0)+m*coeff(c[j,n,p],m,1)
+m^2*coeff(c[j,n,p],m,2): od:
if p>0 then
c[n,n,p]:=-expand(1/2*add(add(c[k,n,s]*c[k,n,p-s],
k=0..n+3*min(s,p-s)),s=1..p-1)):
c[n,n,p]:=collect(c[n,n,p],m):
c[n,n,p]:=coeff(c[n,n,p],m,0)+m*coeff(c[n,n,p],m,1)
+m^2*coeff(c[n,n,p],m,2):
fi: od: end:

> # corrections to matrix elements <m|q^k|n>
mep:=proc(mm,k,n,p)
local s,i,j,s1: global c,m:
s1:=expand(add(add(add(c[i,mm,s]*c[j,n,p-s]*cut(i,mm,s)
*cut(j,n,p-s)*Q(i,k,j),j=0..n+3*(p-s)),i=0..mm+3*s),s=0..p)):
s1:=collect(s1,m):
s1:=coeff(s1,m,0)+m*coeff(s1,m,1)+m^2*coeff(s1,m,2):
end:

> convertac:=proc() local j,k; global c; # to convert from a_j to c_j of V(z)
for k from 1 to 12 do a.k:=(-1)^k*(k+1)/2^k+add((-1)^j*(k+1)!*
c.(k-j)/(2^j*j!*(k+1-j)!), j=0..k-1) od end:

> hw0(2): al2; be2; # for instance

Here endeth this worksheet.