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 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 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 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 and, in the fourth suite, coefficients and in the Herman-Wallis factor 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 for vibration-rotational energies
Procedures tp1 and tp2 in this suite generate algebraic expressions for Dunham's term coefficients for vibration-rotational energies of diatomic molecules in electronic state 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 ; because has a dedicated meaning in Maple , we use simply to represent this conventional quantity, according to these definitions:
= , , , .
To use procedure tp1 to generate components , in which k is a power of ( ), l is a power of [ ] and m is an order of contribution ( for the leading contribution to , for the first correction etc.), of term coefficient according to a sum
= ,
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 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 is from k=8 , l=0 containing coefficients ... to containing with ... . To use procedure tp2 to generate components of term coefficients with symbols having the same meaning as for tp1 , invoke tp2(p2): in which argument is the required order of perturbation theory. This procedure should work with values of p2 up to at least , but the duration of execution of procedure tp2 for a given order is about three times that of tp1 . For , the range of is from with ... to with ... . By default, expressions contain primary rotational parameter , that is reciprocal of harmonic vibrational parameter , and coefficients in Dunham's function for potential energy in terms of . If a user prefers these expression to contain coefficients of with or coefficients of with instead of coefficients , then, after forming expressions with procedure either tp1 or tp2 , invoke procedure converta(l): with argument as either b or c as required. To ensure that resulting expressions in terms of or 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]
> 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]
>
term coefficients for extra-mechanical effects
Procedures in this suite generate additional term coefficients to take account of adiabatic and nonadiabatic contributions to vibration-rotational energies of diatomic molecules in electronic state . To use procedure tp3 to generate term coefficients , in which k is a power of ( ) and l is a power of [ ], invoke tp3(p3): , in which argument p3 is the required order of perturbation theory. This procedure works with values of up to at least 12. All expressions are automatically generated, and are displayed simply on typing Z[k,l]; with appropriate values of k and l . For , the range of is . To convert coefficients , and as coefficients of in etc. to and as coefficients of in etc., invoke convertacstu(): after running tp3 ; then display an expression with collect( ... ); , according to examples below.
> 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];
> 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]);
vibrational matrix elements
Procedures in this suite generate vibrational matrix elements 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 and up to in function for potential energy, type P:=6: state(0): state(1): . Particular matrix elements, in an array , are formed on typing me(v,l,v'): . If matrix elements are required to be expressed in terms of coefficients of function for potential energy, after matrix elements 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);
coefficients in Herman-Wallis factor
Procedures in this suite generate quantities used to form a Herman-Wallis factor for electric-dipolar transitions with J ' - J = +1 or -1, with
= [ 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 ;
Conventional Herman-Wallis coefficients and according to a formula
= (1 + + )
hence have these expressions:
=
= +
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 (0, v ') and (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 and are formed in a succeeding stage; to maintain consistency in powers of in this calculation of , only the first term of must be included.
> 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:
> 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;
> 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;
