Application Center - Maplesoft

App Preview:

Vibrational and rotational spectrometry of diatomic molecules

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application


 

dimolspc.mws

Symbolic Programmes for

Analysis of Vibrational and Rotational Spectra of Diatomic Molecules

F. M. Fernandez and J. F. Ogilvie @ 2001 May 12

This Maple work sheet contains procedures in four suites described in journal MapleTech [volume 5 , 1998, No. 1, pages 42 - 46] for symbolic computation in Maple 6 to form the principal expressions required to undertake a comprehensive analysis of vibration-rotational spectra of diatomic molecules in an electronic state ``^`1` 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, coefficients C[0] ``^`v'` and D[0] ``^`v'` in 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 indicated monograph.

term coefficients Y[k,l] for vibration-rotational energies

Procedures tp1 and tp2 in this suite generate algebraic expressions for Dunham's term coefficients Y[kl] for vibration-rotational energies of diatomic molecules in electronic state ``^`1`*Sigma^`+` with a method hypervirial perturbation theory (i.e. perturbation theory based on hypervirial and Hellmann-Feynman theorems); quantity JB is considered to be independent of perturbation parameter lambda ; because gamma has a dedicated meaning in Maple , we use simply g to represent this conventional quantity, according to these definitions:

lambda^2 = 2*B[e]/W[e] = gamma = g , vh = v+1/2 , JB = JJ*g/2 , JJ = J*(J+1) .

To use procedure tp1 to generate components Y[k,l,m] , in which k is a power of ( v+1/2 ), l is a power of [ J*(J+1) ] and m is an order of contribution ( m = 1 for the leading contribution to Y[k,l] , m = 2 for the first correction etc.), of term coefficient Y[k,l] according to a sum

Y[k,l] = sum(Y[k,l,m],m = 1 .. n) ,

invoke tp1(p1): in which argument p1 is the required order of perturbation theory. This procedure works with values of p1 up to at least p1 = 14 even with Maple V 5; with an extended work space available with Maple 6 a larger order of perturbation theory is practicable. All expressions Y[k,l,m] according to a selected order of perturbation theory are automatically generated, and are displayed 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 containing coefficients a[1] ``^14 ... a[14] to k = 0, l = 8 containing with a[1] ``^6 ... a[6] .
To use procedure
tp2 to generate components Y[k,l,m] of term coefficients Y[k,l] with symbols having the same meaning as for tp1 , invoke tp2(p2): in which argument p2 is the required order of perturbation theory. This procedure should work with values of p2 up to at least p2 = 24 , but the duration of execution of procedure 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] contain primary rotational parameter B[e] , W1 that is reciprocal of harmonic vibrational parameter omega[e] , and coefficients a[j] in Dunham's function V(x) for potential energy in terms of x = (R-R[e])/R[e] . If a user prefers these expression to contain coefficients b[j] of V(y) with y = (R-R[e])/R or coefficients c[j] of V(z) with z = 2*(R-R[e])/(R+R[e]) instead of coefficients a[j] , then, after forming expressions Y[k,l,m] with procedure either tp1 or tp2 , invoke 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, according to examples below.

> 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;
end if:
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)):
end if:
end if:
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))):
end if:
end do:
end do:
En[0] := En[0] + JJ;
for i from 0 to p by 2 do
xtrct(i)
end do;
gc();
end proc:

> trunca := proc(n) # to eliminate unwanted terms for negative parameters
if n<0 then 0;
else 1;
end if:
end proc:

> xtrct := proc(n) # to extract Y[k,l,m] from eigenvalues En[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)));
end do;
end do;
En[n] := 0
end proc:

> converta := proc(l) # to convert from coefficients of V(x) to those of V(y) or V(z)
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 + add((-1)^j*(k+1)!*
c||(k-j)/(2^j*j!*(k+1-j)!), j=0..k-1);
end do;
elif l=b then
for k from 1 to 22 do
a||k := (-1)^k*(k+1) + add((-1)^j*(k+1)!*
b||(k-j)/(j!*(k+1-j)!), j=0..k-1);
end do;
end if;
end proc:

> 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
end if:
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:
end do:
end do:
end do:
end do:
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)):
end if:
end if:
end do:
end do:
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))):
end if;
end do;
end do;
end do;
end do;
e[0,1,0] := 1:
e[1,0,0] := 1:
extrct(p);
gc();
end proc:

> extrct := proc(n) # to extract Y[k,l,m] from eigenvalues En[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]));
end do;
end do;
end do;
end proc:

> tp1(12): # for instance for perturbation theory of order 12
Y[4,1,1] := Y[4,1,1]; # to display leading contribution to Y[4,1] in terms of a[j]

Y[4,1,1] := Be^5*W1^4*(2430855/256*a1^7+2430855/256...
Y[4,1,1] := Be^5*W1^4*(2430855/256*a1^7+2430855/256...
Y[4,1,1] := Be^5*W1^4*(2430855/256*a1^7+2430855/256...
Y[4,1,1] := Be^5*W1^4*(2430855/256*a1^7+2430855/256...
Y[4,1,1] := Be^5*W1^4*(2430855/256*a1^7+2430855/256...

> converta(c): # to convert from coefficients a[j] to coefficients c[j]

> Y[4,1,1] := sort(simplify(Y[4,1,1]));
# to display leading contribution to Y[4,1] in terms of c[j]

Y[4,1,1] := 15/256*(162057*c1^7+143010*c1^6-743568*...
Y[4,1,1] := 15/256*(162057*c1^7+143010*c1^6-743568*...
Y[4,1,1] := 15/256*(162057*c1^7+143010*c1^6-743568*...
Y[4,1,1] := 15/256*(162057*c1^7+143010*c1^6-743568*...
Y[4,1,1] := 15/256*(162057*c1^7+143010*c1^6-743568*...

>

term coefficients Z[k,l] for extra-mechanical effects

Procedures in this suite 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 ``^`1`*Sigma^`+` .
To use procedure
tp3 to generate term coefficients Z[k,l] , in which k is a power of ( v+1/2 ) and l is a power of [ J*(J+1) ], invoke tp3(p3): , in which argument p3 is the required order of perturbation theory. This procedure works with values of p3 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 p3 = 12 , the range of Z[k,l] is k+l <= 6 .
To convert coefficients
a[j] , K[j], q[j] and d[j] as coefficients of x in V(x) etc. to c[j], s[j], t[j] and u[j] as coefficients of z in V(z) etc., invoke convertacstu(): after running tp3 ; then display an expression with collect( ... ); , according to examples below.

> 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
end if;
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)):
end if:
end if;
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))):
end if:
end do;
end do;
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;
end do;
# sum of contributions adiabatic and nonadiabatic in the energy
DE := simplify(We*(s1 + s2) + s3);
xtrct(P);
gc();
end proc:

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

> 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))
end do;
end do;
end do;
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];
end if;
end do;
end do;
end proc:

> 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;
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);
end do:
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);
end do;
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);
end do;
q||0 := t||0;
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)
end do;
d||0 := s||0;
end proc:

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

Z[2,1] := (-3/2*a1*d2+3/2*d3+(3/8*a1^2-3/2*a2+9/4*a...
Z[2,1] := (-3/2*a1*d2+3/2*d3+(3/8*a1^2-3/2*a2+9/4*a...
Z[2,1] := (-3/2*a1*d2+3/2*d3+(3/8*a1^2-3/2*a2+9/4*a...
Z[2,1] := (-3/2*a1*d2+3/2*d3+(3/8*a1^2-3/2*a2+9/4*a...
Z[2,1] := (-3/2*a1*d2+3/2*d3+(3/8*a1^2-3/2*a2+9/4*a...

> convertacstu():

> Z[2,1] := collect(expand(Z[2,1]),
[g,Be,s0,s1,s2,s3,s4,s5,s6,t0,t1,t2,t3,t4,t5,t6,u1,u2,u3,u4,u5,u6]);

Z[2,1] := ((423/8*c1^3-45/2*c4+285/4*c3*c1+33/8*c1+...
Z[2,1] := ((423/8*c1^3-45/2*c4+285/4*c3*c1+33/8*c1+...
Z[2,1] := ((423/8*c1^3-45/2*c4+285/4*c3*c1+33/8*c1+...
Z[2,1] := ((423/8*c1^3-45/2*c4+285/4*c3*c1+33/8*c1+...

>

vibrational matrix elements

Procedures in this suite generate vibrational matrix elements <`v' |x^l | v`> of the reduced displacement coordinate x to various non-negative powers between vibrational states v and v' .
To use this suite, specify 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] in function V(x) for potential energy, type P:=6: state(0): state(1): . Particular matrix elements, in an 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, according to examples below.

> 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));
end if:
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)));
end if:
end do:
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));
end if:
end do:
end proc:

> Q := proc(i,n,j) # matrix elements <i|q^n|j> in basis set of harmonic oscillator
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))
end if:
end proc:

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

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

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

> mep := proc(m,k,n,p) # contributions to matrix elements <m|q^k|n>
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 proc:

> convertac := proc() # to convert from a_j to c_j of V(z)
local j,k;
global c;
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);
end do;
end proc:

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

me(0,1,1) := 1/2*sqrt(g)*sqrt(2)+g^(3/2)*(11/32*sqr...
me(0,1,1) := 1/2*sqrt(g)*sqrt(2)+g^(3/2)*(11/32*sqr...
me(0,1,1) := 1/2*sqrt(g)*sqrt(2)+g^(3/2)*(11/32*sqr...

>

coefficients in Herman-Wallis factor

Procedures in this suite generate quantities used to form a Herman-Wallis factor F[0,v*`'`](iota) for electric-dipolar transitions with J ' - J = +1 or -1, with

iota = 1/2 [ J ' ( J ' + 1) - J ( J+ 1)]

> F[0,`v'`](iota) = 1 + 2 *alpha||` v' `* iota / `<0| M(x) |v'>` +
(alpha||`v' `^2 / `<0| M(x) |v'>`^2 + 2* beta||` v'` / `<0| M(x) |v'>`) *iota^2 ;

F[0,`v'`](iota) = 1+2*`alpha v' `*iota/`<0| M(x) |v...

Conventional Herman-Wallis coefficients C[0] ``^`v'` and D[0] ``^`v'` according to a formula

F[0,`v'`](iota) = (1 + C[0] ``^`v'` iota + D[0] ``^`v'` iota^2 )

hence have these expressions:

C[0] ``^`v'` = 2*alpha || `v'`/<`0 | M(x) | v'`>

D[0] ``^`v'` = alpha || `v'`^2/(<`0 | M(x) | v'`>^2) + 2*beta || `v'`/<`0 | M(x) | v'`>

To use this suite to form a Herman-Wallis factor for a band in absorption connecting vibrational ground state v =0 to vibrationally excited state v ', invoke hw0(n): in which n is the numerical value of v '; results of this calculation are alpha (0, v ') and beta (0, v '), displayed on typing alpha||n; and beta||n; . For instance, for the band v =2 <-- v =0, input hw0(2): , then alpha||2; and beta||2; . Numerators of coefficients C[0] ``^`v'` and D[0] ``^`v'` are formed in a succeeding stage; to maintain consistency in powers of gamma in this calculation of D[0] ``^`v'` , only the first term of C[0] ``^`v'` must be included.

> restart:

> hw0 := proc(v) # main procedure to form Herman-Wallis coefficients
local i,lista,dip:
global g,l,JJ,m,P:
P := v + 4: # greatest order of perturbation
# m = iota = 1/2 [J'(J'+1) - J(J+1)], with 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).
alpha||v := sort(collect(combine(coeff(dip,m,1)/l^v,radical), lista)):
alpha||0 := 0:
beta||v := sort(collect(combine(coeff(dip,m,2)/l^v,radical), lista)):
alpha||v := collect(subs(l=g^(1/2), g^(v/2)*alpha||v), g):
beta||v := collect(subs(l=g^(1/2), g^v*beta||v), g):
NULL:
end proc:

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

> Q := proc(i,n,j) # matrix elements <i|q^n|j> in basis set for harmonic oscillator
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))
end if:
end proc:

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

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

> state := proc(n) # corrections to energy e[n,p] and wavefunction c[j,n,p]
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));
end if:
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)));
end if:
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):
end do:
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):
end if:
end do:
end proc:

> mep := proc(mm,k,n,p) # corrections to matrix elements <m|q^k|n>
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 proc:

> convertac:=proc() # to convert from a||j of V(x) to c||j of V(z)
local j,k;
global c;
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);
end do;
end proc:

Here follow some instances of calling these procedures to form numerators of Herman-Wallis coefficients.

> for j from 0 to 4 do # for instance
hw0(j):
alpha||j := alpha||j;
beta||j := beta||j;
end do;

alpha0 := 0

beta0 := g^2*M1

alpha1 := (-1/2*sqrt(2)*M2+(-41/16*sqrt(2)*a1^2+9/4...

beta1 := g^3*(sqrt(2)*M2+(-3/8*sqrt(2)*a1-3/8*sqrt(...

alpha2 := (-3/2*sqrt(2)*M3+(33/8*sqrt(2)*a1+9/8*sqr...
alpha2 := (-3/2*sqrt(2)*M3+(33/8*sqrt(2)*a1+9/8*sqr...

beta2 := g^4*(3/2*sqrt(2)*M3+(-1/4*sqrt(2)*a1-3/4*s...

alpha3 := ((27/4*sqrt(3)*a1+9/4*sqrt(3))*M3+(-7/2*s...
alpha3 := ((27/4*sqrt(3)*a1+9/4*sqrt(3))*M3+(-7/2*s...
alpha3 := ((27/4*sqrt(3)*a1+9/4*sqrt(3))*M3+(-7/2*s...
alpha3 := ((27/4*sqrt(3)*a1+9/4*sqrt(3))*M3+(-7/2*s...

beta3 := g^5*((3/8*sqrt(3)*a1-9/8*sqrt(3))*M3+(-39/...
beta3 := g^5*((3/8*sqrt(3)*a1-9/8*sqrt(3))*M3+(-39/...

alpha4 := ((5/8*sqrt(6)*a1^2+10*sqrt(6)*a2-105/8*sq...
alpha4 := ((5/8*sqrt(6)*a1^2+10*sqrt(6)*a2-105/8*sq...
alpha4 := ((5/8*sqrt(6)*a1^2+10*sqrt(6)*a2-105/8*sq...
alpha4 := ((5/8*sqrt(6)*a1^2+10*sqrt(6)*a2-105/8*sq...
alpha4 := ((5/8*sqrt(6)*a1^2+10*sqrt(6)*a2-105/8*sq...
alpha4 := ((5/8*sqrt(6)*a1^2+10*sqrt(6)*a2-105/8*sq...
alpha4 := ((5/8*sqrt(6)*a1^2+10*sqrt(6)*a2-105/8*sq...
alpha4 := ((5/8*sqrt(6)*a1^2+10*sqrt(6)*a2-105/8*sq...

beta4 := g^6*((-57/16*sqrt(6)*a1^2+15/4*sqrt(6)*a2-...
beta4 := g^6*((-57/16*sqrt(6)*a1^2+15/4*sqrt(6)*a2-...
beta4 := g^6*((-57/16*sqrt(6)*a1^2+15/4*sqrt(6)*a2-...
beta4 := g^6*((-57/16*sqrt(6)*a1^2+15/4*sqrt(6)*a2-...

> convertac(): # to convert from a||j to c||j

> for j from 0 to 4 do
alpha||j := collect(expand(alpha||j), [g,M0,M1,M2,M3,M4]):
beta||j := collect(expand(beta||j), [g,M0,M1,M2,M3,M4]):
end do;

alpha0 := 0

beta0 := g^2*M1

alpha1 := ((-1/8*sqrt(2)-2*sqrt(2)*c1-41/16*sqrt(2)...

beta1 := (sqrt(2)*M2-3/8*M1*sqrt(2)*c1)*g^3

alpha2 := ((-11/16*sqrt(2)+267/64*sqrt(2)*c1+63/16*...
alpha2 := ((-11/16*sqrt(2)+267/64*sqrt(2)*c1+63/16*...

beta2 := (sqrt(2)*M0+(-1/4*sqrt(2)-3/8*sqrt(2)*c1-9...

alpha3 := ((3/16*sqrt(3)*c3*c1+37/32*sqrt(3)*c2*c1^...
alpha3 := ((3/16*sqrt(3)*c3*c1+37/32*sqrt(3)*c2*c1^...
alpha3 := ((3/16*sqrt(3)*c3*c1+37/32*sqrt(3)*c2*c1^...
alpha3 := ((3/16*sqrt(3)*c3*c1+37/32*sqrt(3)*c2*c1^...

beta3 := ((-2*sqrt(3)+1/2*sqrt(3)*c1)*M0+(3/2*sqrt(...
beta3 := ((-2*sqrt(3)+1/2*sqrt(3)*c1)*M0+(3/2*sqrt(...

alpha4 := ((-331/48*sqrt(6)-73/3*sqrt(6)*c2+1425/12...
alpha4 := ((-331/48*sqrt(6)-73/3*sqrt(6)*c2+1425/12...
alpha4 := ((-331/48*sqrt(6)-73/3*sqrt(6)*c2+1425/12...
alpha4 := ((-331/48*sqrt(6)-73/3*sqrt(6)*c2+1425/12...
alpha4 := ((-331/48*sqrt(6)-73/3*sqrt(6)*c2+1425/12...
alpha4 := ((-331/48*sqrt(6)-73/3*sqrt(6)*c2+1425/12...
alpha4 := ((-331/48*sqrt(6)-73/3*sqrt(6)*c2+1425/12...

beta4 := ((43/12*sqrt(6)-2*sqrt(6)*c1+3/16*sqrt(6)*...
beta4 := ((43/12*sqrt(6)-2*sqrt(6)*c1+3/16*sqrt(6)*...
beta4 := ((43/12*sqrt(6)-2*sqrt(6)*c1+3/16*sqrt(6)*...
beta4 := ((43/12*sqrt(6)-2*sqrt(6)*c1+3/16*sqrt(6)*...

>

>