## File: Gaunt.c

package info (click to toggle)
openmx 3.7.6-1
 `123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222` ``````#include #include #include #include "openmx_common.h" #define S3J_0 1e-10 #define S3J_MAX_FACT 40 #define S3J_EQUAL(a,b) (fabs((a)-(b))(b)?(ris=(a)):(ris=(b)))>(c)?ris:(ris=(c))) #define S3J_MIN(a,b,c,ris) (((a)<(b)?(ris=(a)):(ris=(b)))<(c)?ris:(ris=(c))) static double Clebsch_Gordan(int j1, int m1, int j2, int m2, int j, int m); static double s3j(double j1, double j2, double j3, double m1, double m2, double m3); double Gaunt(int l, int m, int l1, int m1, int l2, int m2) { /************************************************************ Ref. Eq.(3.7.73) in Modern Quantum Mechanics by J.J.Sakurai ************************************************************/ int Ls; double tmp0,tmp1,tmp2,tmp3; double result,cleb1,cleb2; cleb1 = Clebsch_Gordan(l1,0,l2,0,l,0); cleb2 = Clebsch_Gordan(l1,m1,l2,m2,l,m); tmp0 = 2.0*(double)l1 + 1.0; tmp1 = 2.0*(double)l2 + 1.0; tmp2 = 4.0*PI*(2.0*(double)l + 1.0); tmp3 = sqrt(tmp0*tmp1/tmp2); result = tmp3*cleb1*cleb2; return result; } double Clebsch_Gordan(int j1, int m1, int j2, int m2, int j, int m) { int esp; double cgris; esp=(int)(j1-j2+m); if (!S3J_EQUAL(esp,j1-j2+m)) return 0; if (esp%2==0) cgris=1.0; else cgris=-1.0; cgris*=sqrt(2*j+1)*s3j((double)j1,(double)j2,(double)j, (double)m1,(double)m2,(double)(-m)); return cgris; } double s3j(double j1, double j2, double j3, double m1, double m2, double m3) { /****************************************************************************** ************************************************************** This program was adopted from the following webpage by T.Ozaki at Sep. 13 2002. T.Ozaki greatly thanks the original author (Dr. Gusmeroli). ************************************************************** http://www.ph.surrey.ac.uk/~phs3ps/cleb.html The below is the original message. Name: s3j Evaluates 3j symbol Author: Riccardo Gusmeroli (rgusmero@elet.polimi.it) Notes: - defining S3J_TEST enables the compilation of a very small test suite. - the maximum allowed factorial is S3J_MAX_FACT (currently 25!). This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. ( j1 j2 j3 ) ( ) = delta(m1+m2+m3,0) * (-1)^(j1-j2-m3) * ( m1 m2 m3 ) +- | (j1+j2-j3)! (j1-j2+j3)! (-j1+j2+j3)! * | -------------------------------------- ... | +- -+ 1/2 (j1-m1)! (j1+m1)! (j2-m2)! (j2+m2)! (j3-m3)! (j3+m3)! | ... ------------------------------------------------------- | * (j1+j2+j3+1)! | -+ +--- \ (-1)^k * | --------------------------------------------------------------------- / k! (j1+j2-j3-k)! (j1-m1-k)! (j2+m2-k)! (j3-j2+m1+k)! (j3-j1-m2+k)! +--- k Where factorials must have non-negative integral values: j1+j2-j3 >= 0 j1-j2+j3 >= 0 -j1+j2+j3 >= 0 j1+j2+j3+1 >= 0 k >= 0 j1+j2-j3-k >= 0 j1-m1-k >= 0 j2+m2-k >= 0 j3-j2+m1+k >= 0 j3-j1-m2+k >= 0 The 3j symbol is therefore non-null if j1+j2 >= j3 (1) j1+j3 >= j2 (2) j2+j3 >= j1 (3) and k values in the sum must be such that k <= j1+j2-j3 (4) k >= 0 (7) k <= j1-m1 (5) k >= -j3+j2-m1 (8) k <= j2+m2 (6) k >= -j3+j1+m2 (9) If no values of k satisfy the (4) to (9), the result is null because the sum is null, otherwise one can find kmin < kmax such that kmin <= k <= kmax (4) to (6) => kmin=MAX(j1+j2-j3, j1-m1, j2+m2 ) (7) to (9) => kmax=MIN(0, -j3+j2-m1, -j3+j1+m2 ) The condition kmin < kmax includes (1) to (3) because (4) and (7) => (1) (5) and (8) => (2) (6) and (9) => (3) Once the values of kmin and kmax are found, the only "selection rule" is kminkmax) return 0.0; ris=0.0; if (kmin%2==0) mult=1.0; else mult=-1.0; for (k=kmin; k<=kmax; ++k) { ris+=mult/(f[k]*f[j1pj2mj3-k]*f[jmm1-k]*f[jpm2-k]*f[j3mj2pm1+k]*f[j3mj1mm2+k]); mult=-mult; } /* (-1)^(j1-j2-m3)=(-1)^(j1-j2-m3+m1+m2+m3)=(-1)^(jpm1-jmm2) */ if ((jpm1-jmm2)%2!=0) ris=-ris; ris*=sqrt(f[j1pj2mj3]*f[jpm1-jmm2+jpm3]*f[-jmm1+jpm2+jpm3]* f[jpm1]*f[jpm2]*f[jpm3]*f[jmm1]*f[jmm2]*f[jmm3]/ f[jpm1+jpm2+jpm3+1]); return ris; } ``````