
(*     Mathematica Package for Evaluation 
       of Matrix Elements of Coulomb Interaction 
        between Electrons

        Package for Clebsch Gordans; 6-j package
        Matrix elements evaluation  
  
                  Prepared by: L. Kocbach
                  based on work of 
    James M. Feagin, Ladislav  Kocbach, Richard Liska   and  Alonso Nunez
    There is also a Notebook and HTML document
 *)
BeginPackage["`CoulMEpack`"]
(*  
 PART I:  Package of Clebsch Gordans : Origin James M. Feagin
 *)

(**** Uses Racah s formula to compute Clebsch-Gordan coefficient 
      and its related 3j symbol. Definitions of CG coefficient and 
      3j symbol are as in 
	  BRINK & SATCHLER, RAngular Momentum,S 2nd ed., p.34, eq. 2.34,
		and App. I, p.136.   
****)

(**** Entered 1/26/91 constraints:
		C < Abs[A-B] || C > A+B.
		Abs[a] > A || Abs[b] > B || Abs[c] > C.  
****)


CoulME::usage =
        "Clebsch Gordans, 6-j and Coulomb 1 and 2 electron integrals
         Clebsch  Threej   Sixj       Redmat   \n
         RadHyd   MatEl    DoubMatEl  Corremat"
Clebsch::usage =
	"Clebsch[{A,a},{B,b},{C,c}] gives the Clebsch-Gordan coefficient 
	for the addition of angular momenta C = A + B using RacahUs formula 
	(see Brink&Satchler, RAngular Momentum,S 2nd ed.(Oxford), p.34)."

Threej::usage =	
	"Threej[{A,a},{B,b},{C,c}] gives the 3-j symbol for 
	combining angular momenta, which is proportional to the CG coefficient."
Sixj::usage =
	" Sixj[{A ,B , C },{a ,b ,c }]  gives the 6-j symbol"
Redmat::usage =
	"Redmat[A ,K ,B ]  gives the reduced M.E. of spherical harmonics
	The usual notation  (A|| K ||B)"
RadHyd::usage = 
	"RadHyd[n  ,  l  ,  r  , Z   ]   
	Radial Functions of Hydrogen-like ion with Z atomic number"	
MatEl::usage =
      "MatEl[ Z , n1  ,  l1  , n2  ,  l2  , LL  , rr   ] 
      Matrix Element (Radial part) of a particle - electron interaction	" 
DoubMatEl::usage =  
       "DoubMatEl[Z ,n1 ,l1 ,n2 ,l2 ,n3 ,l3 ,n4 ,l4 ,L ] 
             Reduced  Matrix Element (Radial part) 
             of  electron - electron  repulsion" 
Corremat::usage =
  "Corremat[Z ,na ,la ,nb ,lb ,nc ,lc ,nd ,ld ,L ] 
  Reduced  Matrix Element (Total) of electron - electron  repulsion"                           

Begin["`private`"]

(* see Brink&Satchler, p.140 *)
Clebsch[{A_,a_},{A_,b_},{0,0}] := (-1)^(A-a)/Sqrt[2A+1] /; b == -a 
Clebsch[{A_,a_},{B_,b_},{0,0}] := 0 

f1[C_,A_,B_,a_,b_] := 
	Sum[ (-1)^k (k! (A+B-C-k)! (A-a-k)! (B+b-k)! (C-B+a+k)! 
	     (C-A-b+k)!)^-1, 
		{ k, Max[0,B-C-a,A-C+b], Min[A+B-C,A-a,B+b] } ]  
		(* These Max and Min prevent the arguments of 
	      Factorial (!) from becoming negative when the sum over k is 
	      performed. *)
	      
f2[C_,A_,B_] := (-C+A+B)! (C-A+B)! (C+A-B)! /(1+C+A+B)!

f3[C_,A_,B_,c_,a_,b_] := (1+2C) (C-c)! (C+c)! (A-a)! (A+a)! (B-b)! (B+b)!

(* see Brink&Satchler, p.34, eq. 2.34 *)
Clebsch[ {A_, a_}, {B_, b_}, {C_, c_} ] :=
		If[ c!=a+b || C < Abs[A-B] || C > A+B || 
			Abs[a] > A || Abs[b] > B || Abs[c] > C, 0 , 
			Block[ {}, f1[C,A,B,a,b] Sqrt[ f2[C,A,B] 
			f3[C,A,B,c,a,b] ]
		   ]
		]

(* see Brink&Satchler, App. I, p.136 *)
Threej[ {A_,a_},{B_,b_},{C_,c_} ] :=
	(-1)^(B-A+c)/Sqrt[2C+1] Clebsch[{A,a},{B,b},{C,-c}]
	
	

(*This is a modification of the matrix elements evaluation
 to include  Z- dependence     :  
 Origin    L. Kocbach and R. Liska;  
 Modified by Alonso Nunez
 *)
Rad[n_ ,  l_ ,  r_ , Z_  ]:=
LaguerreL[n-l -1 , 2 l+1,2 Z r/n] Exp[- Z r/ n] ( Z r )^l

RadHyd[n_ ,  l_ ,  r_ , Z_  ] :=  
    Rad[n ,  l  ,  r  , Z ] / Sqrt[Norm2 [ n , l , Z    ] ]

MatEl[ Z_, n1_ ,  l1_ , n2_ ,  l2_ , LL_ , rr_  ] :=
Block[
 { fun1, fun2, v, res1, res2, int1, int2 , res},
  fun1=(RadHyd[n1,l1, v , Z]/. E->enat);
   fun2=(RadHyd[n2,l2, v, Z ]/. E->enat); 
   res2=Integrate[Expand[v^2 fun1 fun2 / v^(LL+1)],v];
   res1=Integrate[Expand[v^2 v^LL   fun1  fun2],v];
   int1=( res1/. v->rr) ;
    int1=  int1   -  ( res1/.v->0);
   int2=  -  res2/. v->rr  ;
   res = (   (  int1/rr^(LL+1)  +  rr^(LL)  int2   )/. Log[enat]->1);  
   res/. enat->E 
 ]

Norm2[n1_ ,  l1_ , Z_   ] := 
Block[
 { fun,  v, res,  int },
 fun=(Rad[n1,l1, v, Z]/. E->enat);
   res=Integrate[Expand[v^2 fun fun],v];
   int =  -  res/. v->0  ;
   res =  int /. Log[enat]->1 ;  
   res/. enat->E 
  ]

DoubMatEl[Z_,n1_,l1_,n2_,l2_,n3_,l3_,n4_,l4_,L_] :=
Block[
  {fun1,fun2,  v, res,  int },
  fun1=(RadHyd[n1,l1, v ,Z]/. E->enat);
   fun2=(RadHyd[n4,l4, v ,Z]/. E->enat);
   fun3=(MatEl[Z,n2,l2,n3,l3,L,v]/. E->enat); 
   res=Integrate[Simplify[Expand[v^2 fun1 fun2 fun3 ]],v];
   int =  -  res/. v->0  ;
   res =  int /. Log[enat]->1 ;  
   res/. enat->E 
 ]

(* 
First version of 6-j package :Using Cowan, page 147   
Origin:  Alonso Nunez 
*)

(* Special value for one J zero *)
Clebshing1[C_,A_,B_,a_,b_,0] := (-1)^(A+B+C)/Sqrt[(2A+1)(2B+1)]
(* Used in general expression *)
ClebDelta2[C_,A_,B_] := (-C+A+B)! (C-A+B)! (C+A-B)! /(1+C+A+B)!
(* Used in general expression, Big Sum *)
ClebSum3[C_,A_,B_,c_,a_,b_] := 
Sum[ (-1)^k (k+1)! ((A+B+a+b-k)! (B+C+b+c-k)! (C+A+c+a-k)! 
(k-A-B-C)! (k-A-b-c)! (k-a-B-c)! (k-a-b-C)!)^-1, 
		{ k, Max[0,A+B+C,A+b+c,a+B+c,a+b+C], 
		Min[A+B+b+a,B+C+b+c,C+A+c+a] } ]
		
Sixj[{A_,B_, C_},{a_,b_,c_}] := 

If[ C < Abs[A-B] || C > A+B , 0  (* if not satisfied *),
	Block[  (* if satisfied *)
		{},
		If [ c==0, Clebshing1[C,A,B,a,b,0] , 
			 Sqrt[ ClebDelta2[C,A,B]ClebDelta2[c,b,A]
			       ClebDelta2[c,a,B]  ClebDelta2[C,a,b]] 
             ClebSum3[C,A,B,c,a,b] ]
       	] (* end block *)
  ] (* end if *)
		


Redmat[A_,K_,B_] := 
(-1)^A Sqrt[(2A+1)(2B+1)] Threej[{A,0},{K,0},{B,0}]

Smrk[la_,lb_,lc_,ld_,L_,K_] :=
 (-1)^(lc+lb+L) Sixj[{la,lb,L},{ld,lc,K}] Redmat[la,K,lc]  Redmat[lb,K,ld]
                        
Corremat[Z_,na_,la_,nb_,lb_,nc_,lc_,nd_,ld_,L_] :=
           Sum[ Smrk[la,lb,lc,ld,L,K]
          DoubMatEl[Z,na,la,nb,lb,nc,lc,nd,ld,K]
         ,{K, Max[0,Abs[la-lc],Abs[lb-ld]],
              Min[la+lc,lb+ld] } ]
              
End[ ]
EndPackage[]              



