(**************************************************************************) (* *) (* RSperturb v. 1.5 *) (* Written by John M. Herbert *) (* Last modified 1/26/98 *) (* *) (**************************************************************************) (* This is a general code for application of arbitrary-order Rayleigh Schrödinger perturbation theory *) (* Two variables must be specified prior to using this program: HighestState is the numerical value of the highest quantum state for which formulae are to be derived, while HighestOrder is the highest order of correction for which formulae will be derived. Specifying larger values of these variables than are necessary will not adversely affect the output, although the derivations may require more CPU time than would otherwise be necessary. *) Energy::usage = "Energy[n,v] derives the nth-order energy correction formula for quantum state v using the most general form of the perturbation energy formula (see J. M. Herbert, Technical Memorandum No. 222, Argonne National Laboratory, Mathematics and Computer Science Division, 1997)."; Egy::usage = "Egy[i,v] is the ith-order energy correction for quantum state v as referenced by the function Energy[n,v]. This notation allows a recursive formula for Energy[n,v] to be generated without explicit evaluation of the lower-order energies on which it depends."; KDelta::usage = "The Kronecker delta function. Takes two arguments and returns 1 if they are equal, 0 if they are not."; SumStates::usage = "SumStates[expr, {v,v’,vmax}] sums expr over the quantum numbers from v to vmax, skipping v = v’."; Qdelta::usage = "Qdelta[n] derives expressions for the elements of the matrices Q^1, Q^2, ..., Q^n. (see M. M. Dudas, H. C. Hsieh, and W. C. Ermler, Mathematica Journal, V. 2, p. 66, 1992)"; GenerateCoefficients::usage = "GenerateCoefficients[n,v] returns a matrix called 'coefficients' containing all of the Rayleigh-Schrödinger expansion coefficients necessary to expand the nth perturbed wave function for quantum state v. This function is useful for obtaining the explicit form of the wave function as a linear combination of the zeroth-order wave functions. Note that coefficients[[1]] corresponds to the quantum state v = 0."; RSCoefficients1::usage = "RSCoefficients1[n,j,{v}] generates the nth-order Rayleigh-Schrödinger expansion coefficient associated with quantum state j such that the summation excludes j = v. This function is intended to be called by other Mathematica functions."; RSCoefficients2::usage = "Usage is the same as the function RSCoefficients1; both functions are necessary for matrix elements that require two different expansion coefficients."; Psi::usage = "Psi[n,v] represents the nth-order perturbed wave function for quantum state v."; H::usage = "H[n] represents the nth-order perturbed Hamiltonian operator."; k::usage = "k[i] represents the ith molecular force constant."; Int::usage = "Int[Psi[n1,v1], H[x], Psi[n2,v2]] is a Hamiltonian matrix element involving a perturbed Hamiltonian operator. Int[Psi[n1,v1], Psi[n2,v2]] is an overlap integral."; (* Other symbols appearing in the symbolic output are: ke, the equilibrium molecular force constant; ve, the harmonic frequency; h, Planck’s constant; Be, the equilibrium rotational constant; Re, the equilibrium internuclear separation; a, which represents a constant appearing in the harmonic oscillator wave functions. a is equal to 4 Pi^2 mu ve/h, where mu is the molecular reduced mass (including electrons). *) (* Note that SI units are assumed in these constants *) Off[General::spell]; Off[Array::ilsmn]; Highest = HighestState + 3 + Sum[index, {index, 3, HighestOrder + 1}]; AA = Array[A, Highest]; BB = Array[B, Highest]; Unprotect[Part]; Attributes[KDelta] = {Orderless}; KDelta[a_, b_] := Which[a == b, 1, a != b, 0] /; a == b || a != b SumStates[expr__, {var_, hole_, max_}] := Sum[expr, {var, -max, hole - 1}] + Sum[expr, {var, hole + 1, max}] (* Expansion of perturbed wave functions in Hamiltonian matrix elements *) Int /: Int[Psi[n1_Integer, v1_], h:_H, Psi[n2_Integer, v2_]] := Block[{pdt}, pdt = SumStates[ RSCoefficients1[n2, i, {v2}] Int[Psi[0, j], h, Psi[0, i]], {i, v2, Highest}]; SumStates[ RSCoefficients2[n1, j, {v1}] pdt, {j, v1, Highest}] ] /; (n1 > 0) && (n2 > 0) Int /: Int[Psi[n1_Integer, v1_], h:_H, Psi[n2_Integer, v2_]] := SumStates[ RSCoefficients1[n1, j, {v1}] Int[Psi[0, j], h, Psi[n2, v2]], {j, v1, Highest}] /; (n2 == 0) && (n1 > 0) Int /: Int[Psi[n1_Integer, v1_], h:_H, Psi[n2_Integer, v2_]] := SumStates[ RSCoefficients2[n2, j, {v2}] Int[Psi[n1, v1], h, Psi[0, j]], {j, v2, Highest}] /; (n1 == 0) && (n2 > 0) (* Expansion of perturbed wave functions in overlap integrals *) Int /: Int[Psi[n1_Integer, v_], Psi[n2_Integer, v_]] := SumStates[ RSCoefficients1[n1, j, {v}] RSCoefficients2[n2, j, {v}], {j, v, Highest}] /; (n1 > 0) || (n2 > 0) (* Hermitian property of quantum-mechanical integrals *) Int /: Int[Psi[n1_, v1_], h:_H, Psi[n2_, v2_]] := Int[Psi[n2, v2], h, Psi[n1, v1]] /; Order[Psi[n1, v1], Psi[n2, v2]] == -1 Int /: Int[Psi[n1_, v1_], Psi[n2_, v2_]] := Int[Psi[n2, v2], Psi[n1, v1]] /; Order[Psi[n1, v1], Psi[n2, v2]] == -1 (* Orthogonality of zeroth-order wave functions *) Int /: Int[Psi[0, v1_], Psi[0, v2_]] := KDelta[v1, v2] (* Generalized energy formula *) Energy[n_, V_] := Block[{term1, term2, term3, term4, term5, span, ii, jj}, If[EvenQ[n] == True, span = n/2, span = (n - 1)/2]; term1 = Sum[ (2 - KDelta[n, 2 ii]) Int[Psi[jj - 1, V], H[n - ii - jj + 1], Psi[ii, V]], {jj, 1, span}, {ii, jj, span}]; term2 = Sum[ Int[Psi[ii, V], H[n - 2 ii], Psi[ii, V]], {ii, 0, span - 1}]; term3 = Sum[ (2 - KDelta[n, 2 ii]) Egy[n - ii - jj + 1, V] Int[Psi[jj - 1, V], Psi[ii, V]], {jj, 2, span}, {ii, jj, span}]; term4 = Sum[ Egy[n - 2 ii, V] Int[Psi[ii, V], Psi[ii, V]], {ii, 1, span - 1}]; term5 = KDelta[n, 2 span + 1] (Int[Psi[span, V], H[1], Psi[span, V]] - Egy[1, V] Int[Psi[span, V], Psi[span, V]] (1 - KDelta[1, n])); term1 + term2 - term3 - term4 + term5 ] /; n > 0 (* Normal coordinate integrals *) Qdelta[n_Integer] := Block[{zz, u, yy, gg, dr, y, zed, ii, jj, ym}, ss = Array[s, {16, Highest}]; ss[[1,1]] = ((qnum + 1)/(2 a))^(1/2); ss[[1,2]] = (qnum/(2 a))^(1/2); u = 2; While[u < n + 1, yy = u - 1; ss[[u,1]] = (ss[[yy,1]] /. qnum -> qnum + 1) ss[[1,1]]; zed = 2; gg = u + 1; While[zed < gg, zz = zed - 1; ss[[u,zed]] = (ss[[yy,zz]] /. qnum -> qnum - 1) ss[[1,2]] + (ss[[yy,zed]] /. qnum -> qnum + 1) ss[[1,1]]; zed++]; ym = u + 1; ss[[u,ym]] = (ss[[yy,u]] /. qnum -> qnum - 1) ss[[1,2]]; u++]; dr = Highest - 1 - n; Drop[ss[[n]], -dr]; qdel = Array[qd, n]; Do[ qdel[[ii]] = Sum[ ss[[ii, jj]] KDelta[qp, qnum + ii - 2 (jj - 1)], {jj, ii + 1}], {ii, n}]; qdel[[n]] ] (* Calculation of wavefunction expansion coefficients *) GenerateCoefficients[ord_Integer, state_Integer] := Block[{counter, index1, index2, temp, count2}, temp = Array[t, ord]; coefficients = Array[cc, {ord, state + ord + 3}]; Do[ term1 = - Int[Psi[0, m], H[counter], Psi[0, state]]/ DeltaE[m, state]; term2 = 1/DeltaE[m, state]* Sum[ Egy[counter - index1, state] temp[[index1]], {index1, counter - 1}]; term3 = -1/DeltaE[m, state]* Sum[ Sum[ (temp[[index1]] /. {m -> index2}) Int[Psi[0, m], H[counter - index1], Psi[0, index2]], {index2, state - counter + index1 - 2, state - 1, 2}] + Sum[ (temp[[index1]] /. {m -> index2}) Int[Psi[0, m], H[counter - index1], Psi[0, index2]], {index2, state + counter - index1 + 2, state + 1, -2}], {index1, counter - 1}]; temp[[counter]] = term1 + term2 + term3; Do[ coefficients[[counter, count2]] = temp[[counter]] /. m -> (count2 - 1), {count2, state}]; Do[ coefficients[[counter, count2]] = temp[[counter]] /. m -> (count2 - 1), {count2, state + 2, state + counter + 3}]; coefficients[[counter, state + 1]] = 0, {counter, ord}]; ] RSCoefficients1[n_Integer, index_, {not_}] := Block[{m, vv, coeffs, cc, counter, term1, term2, term3}, coeffs = Array[cc, Highest]; coeffs[[1]] = -(Int[Psi[0, m], H[1], Psi[0, vv]]/DeltaE[m, vv]); counter = 2; While[counter <= n, term1 = -(Int[Psi[0, m], H[counter], Psi[0, vv]]/DeltaE[m, vv]); term2 = Sum[ Egy[counter - AA[[2*counter - 1]], vv] coeffs[[AA[[2*counter - 1]]]], {AA[[2*counter - 1]], 1, counter - 1}]/DeltaE[m, vv]; term3 = -(Sum[ SumStates[ Int[Psi[0, m], H[counter - AA[[2*counter - 1]]], Psi[0, AA[[2*counter]]]] (coeffs[[AA[[2*counter - 1]]]] /. m -> AA[[2*counter]]), {AA[[2*counter]], vv, Highest}], {AA[[2*counter - 1]], 1, counter-1}]/DeltaE[m, vv]); coeffs[[counter]] = term1 + term2 + term3; counter++]; coeffs[[n]] /. {m -> index, vv -> not} ] RSCoefficients2[n_Integer, index_, {not_}] := Block[{m, vv, coeffs, cc, counter, term1, term2, term3}, coeffs = Array[cc, Highest]; coeffs[[1]] = -(Int[Psi[0, m], H[1], Psi[0, vv]]/DeltaE[m, vv]); counter = 2; While[counter <= n, term1 = -(Int[Psi[0, m], H[counter], Psi[0, vv]]/DeltaE[m, vv]); term2 = Sum[ Egy[counter - BB[[2*counter - 1]], vv] coeffs[[BB[[2*counter - 1]]]], {BB[[2*counter - 1]], 1, counter - 1}]/DeltaE[m, vv]; term3 = -(Sum[ SumStates[ Int[Psi[0, m], H[counter - BB[[2*counter - 1]]], Psi[0, BB[[2*counter]]]] (coeffs[[BB[[2*counter - 1]]]] /. m -> BB[[2*counter]]), {BB[[2*counter]], vv, Highest}], {BB[[2*counter - 1]], 1, counter-1}]/DeltaE[m, vv]); coeffs[[counter]] = term1 + term2 + term3; counter++]; coeffs[[n]] /. {m -> index, vv -> not} ]