Sean, Here are some notes I sent to you earlier in the summer. The points made here are quite important..
So there are some things you should know about the Wigner D matrix: * While I implemented it, I have 0 confidence that the current implementation is free of bugs. Because of this, one of the first tasks should be to test this part of the code extensively, both for numerical and symbolic arguments. I would test the D and d matrix code thoroughly *before* moving on to anything else. It would be great to see a pull request for that. * The current method of computing the D matrix is not optimal. This is because we probably want to use a different form for the integer and half integer spin values. This is to get the canonical forms where, for 1/2 integer spins there are things like sin(theta/2) and cos(theta/2), but for integer spins, sin(theta) and cos(theta) appear. This should be implemented and then tested extensively. This will also help in testing and comparing to the tables in Varsh. * Let's use the convention of Varshalovich for the D-matrices. He follows the convention of most of the other main people in the field (see Table 4.2 for a list). * We should use his convention for Euler angles as well. * Because of the intimate relationship between CG, Wigner 3j6j, etc and the rotation matrices, we probably want to have a symbolic DMatrix class that we can use when we are representing various sum and integrals symbolically. This should be separate from the rotation matrix operator and should be in the same category object as the GC and Wigner symbols. This will allow us to implement all of the symmetry relationships that D satisfies as well. * Currently I have implemented one of the differential representations of d(J,m,m') from section 4.3.2. I think that is probably OK for integer spin, but for half integer spins we should use one of the formulas from section 4.3.1 that are a sum over k. This very important I think... * When testing d(j,m,m') we probably want to simply go through the specific cases from Tables 4.3, 4.4, 4.5 (j=1/2,1,3/2). It looks like Ondrej has helped to track down the actual bug in the current implementation, BUT that does solve the issue of the tests and implementing a different approach for 1/2 in integer spins. Cheers, Brian On Thu, Jun 16, 2011 at 8:34 PM, Ondrej Certik <[email protected]> wrote: > On Thu, Jun 16, 2011 at 6:38 PM, Ondrej Certik <[email protected]> > wrote: >> On Thu, Jun 16, 2011 at 12:04 AM, Sean Vig <[email protected]> wrote: >>> Hi all, >>> In working on some stuff with spin states, I ran into some problems with the >>> current implementation of the Wigner small-d matrix, Rotation.d in >>> sympy.physics.quantum.spin. I had written methods to change bases using the >>> Wigner D-function [0] and in testing decided to try >>>>>>qapply(JzBra(1,1)*JzKet(1,1).rewrite('Jx')) >>> which should be 1, as the spin ket is rewritten in the Jx basis and then >>> back to the Jz basis to apply the innerproduct, but I have that it gives 0. >>> I traced this back to a bug in the Rotation.d function, which currently has >>> an open issue [1]. For the Wigner D-function and the small d-matrix, the >>> conventions laid out in Varshalovich "Quantum Theory of Angular Momentum". >>> It seems the d-matrix is fine for positive values of angle argument, but >>> does not obey the symmetry d(j,m,mp,-beta)=(-1)**(m-mp)*d(j,m,mp,beta) and >>> does not agree with the tables in Varshalovich. Those terms that fail for >>> the j=1 case are in an XFAIL test in my branch. What is odd is that when I >> >> >> I think there is a bug in the code. See my comment here: >> >> http://code.google.com/p/sympy/issues/detail?id=2423#c3 >> >> The correct results for general beta are: >> >> d(1, 0, 1, beta) = sin(beta)/sqrt(2) >> d(1, 1, 0, beta) = -sin(beta)/sqrt(2) >> >> However, sympy gives: >> >>>>> Rotation.d(1,0,1,beta) >> ⎽⎽⎽ >> ╲╱ 2 ⋅(2⋅cos(β) + 2) >> ──────────────────── >> 4 >>>>> Rotation.d(1,1,0,beta) >> ⎽⎽⎽ >> -╲╱ 2 ⋅(cos(β) + 1) >> ─────────────────── >> 2 >> >> Which is wrong (it looks quite ok, that it changes sign, as it should, >> but something is wrong with the cos(beta) thing). The sympy code >> implements the "Eq. 7 in Section 4.3.2 of Varshalovich." >> >> >>> ran the equation used to define the d-matrix through Mathematica, I got >>> results that agreed with the sympy output, so the problem may be in the >>> equation and not a bug in the code. If anyone could take a look at that, I'd >>> appreciate it. >> >> Which *exact* equation did you run through Mathematica? Eq. 7 in >> section 4.3.2? What *exactly* did you get? Did you get an expression >> involving cos(beta), just like sympy above? Can you paste here the >> Mathematica code? I'll run it with my Mathematica, to verify, that we >> didn't make a mistake. >> >> >> Now we just need to systematically look at it, and nail it down. We >> need to get the correct expressions: >> >> d(1, 0, 1, beta) = sin(beta)/sqrt(2) >> d(1, 1, 0, beta) = -sin(beta)/sqrt(2) >> >> one way or the other, for general beta. Then things will start to work. >> >> Sean, let me know if you have any questions to the above. >> >> Once we fix this, we'll move on to the other problems you raised. >> >> I am CCing Brian, who implemented that code in SymPy. However, we >> should be able to fix this ourselves Sean ---- all we need to do is to >> take the eq. (7), and see what expression we get for J=1, M=1, M'=0, >> just put it there by hand (don't use mathematica) and see what you >> get. >> >> Post here your results, and I'll verify them with my independent >> calculation and we'll nail it down. > > Ok, I think that I have nailed it down. There are actually several problems: > > 1) First of all, this is the range for beta: > > 0 <= beta <= pi > > see Varshalovich, page 74. > > 2) See the attached screenshot of my calculation, which calculates > d(j, 0, 1, beta) from the equation (7) on page 75, and shows, that it > is equal to > sin(beta)/sqrt(2), consistent with Varshalovich result in the table 4.4. > > > As such, from 2) it follows, that the sympy result (see my previous > email) is *wrong*, as can be checked by substituting beta = pi/3 and > checking against sin(beta)/sqrt(2). For beta=pi/2, it happens to be > equal, but that is pure accident. > > From 1) it follows, that you can't call it for beta=-pi/2, you need to > adjust alpha and gamma instead. > > So all problems with the wigner function seems to be resolved. Sean, > would you be able to take it from here and find the bug in the code > and fix it? > > Ondrej > > -- > You received this message because you are subscribed to the Google Groups > "sympy" group. > To post to this group, send email to [email protected]. > To unsubscribe from this group, send email to > [email protected]. > For more options, visit this group at > http://groups.google.com/group/sympy?hl=en. > > -- Brian E. Granger Cal Poly State University, San Luis Obispo [email protected] and [email protected] -- You received this message because you are subscribed to the Google Groups "sympy" group. To post to this group, send email to [email protected]. To unsubscribe from this group, send email to [email protected]. For more options, visit this group at http://groups.google.com/group/sympy?hl=en.
