(* Calculations for arbitrary atom, including p-orbitals *) Clear[f,l,m,n,r,theta,phi,fs,fp,x,la,ma,na,lb,mb,nb,int,A,a,b,intne,int0, intt,int0ee,intee,l1a,m1a,n1a,l2a,m2a,n2a,l1b,m1b,n1b,l2b,m2b,n2b,a1,a2,b1,b2]; npart=2; (* int0[n,a]=Integral[Exp[-a r] r^n dr] *) int0[n_,a_]=(n+2)!/a^(n+3); (* int[la,ma,na,lb,mb,nb]= *) int[la_,ma_,na_,lb_,mb_,nb_]:=0; int[l_,m_,na_,l_,m_,nb_]=int0[2 l,A[na]+A[nb]]; (* intne[la,ma,na,lb,mb,nb]= *) intne[la_,ma_,na_,lb_,mb_,nb_]:=0; intne[l_,m_,na_,l_,m_,nb_]=int0[2 l-1,A[na]+A[nb]]; (* intt[la,ma,na,lb,mb,nb]= where T=L^2/(2 r^2) -1/(2 r^2) d/dr(r^2 d/dr) *) (* psia=r^l Exp[-A[na] r];psib=r^l Exp[-A[nb] r]; Unprotect[Re];Re[x_]:=1; i=Integrate[r^2 psia (-1/(2 r^2) D[r^2 D[psib,r],r]),{r,0,Infinity}]; Simplify[(l(l+1)/2 int0[2 l-2,A[na]+A[nb]]+ Integrate[r^2 psia (-1/(2 r^2) D[r^2 D[psib,r],r]),{r,0,Infinity}])/.{ Gamma[2 l]->(2 l)!/(2 l),Gamma[2(l+1)]->(2 l)! (2 l+1)}] *) intt[la_,ma_,na_,lb_,mb_,nb_]:=0; intt[l_,m_,na_,l_,m_,nb_]= (1 + 3*l + 2*l^2)*A[na]*A[nb]*(A[na] + A[nb])^(-3 - 2*l)*(2*l)!; (* int0ee[l1a,m1a,l2a,m2a,l1b,m1b,l2b,m2b,a1,a2]= where A[n1a]+A[n1b]=a1, A[n2a]+A[n2b]=a2 *) int0ee[l1a_,m1a_,l2a_,m2a_,l1b_,m1b_,l2b_,m2b_,a1_,a2_]:=0; (* (2,2) *) int0ee[0,0,0,0, 0,0,0,0, a1_,a2_]=2*(a1^2 + 3*a1*a2 + a2^2)/(a1^2*a2^2*(a1 + a2)^3); (* 3*(10,10) *) int0ee[1,0,0,0, 1,0,0,0, a1_,a2_]=(12*(2*a1^4 + 10*a1^3*a2 + 10*a1^2*a2^2 + 5*a1*a2^3 + a2^4))/(a1^4*a2^2*(a1 + a2)^5); int0ee[0,0,1,0, 0,0,1,0, a1_,a2_]=int0ee[1,0,0,0, 1,0,0,0, a1,a2]; (* 3*(10,17) *) int0ee[1,0,0,0, 0,0,1,0, a1_,a2_]=(8*(a1^2 + 5*a1*a2 + a2^2))/(a1^2*a2^2*(a1 + a2)^5); int0ee[0,0,1,0, 1,0,0,0, a1_,a2_]=int0ee[1,0,0,0, 0,0,1,0, a1,a2]; (* 3/2*((8,8)-(10,10)) *) int0ee[1,1,0,0, 1,1,0,0, a1_,a2_]=int0ee[1,0,0,0, 1,0,0,0, a1,a2]; int0ee[0,0,1,1, 0,0,1,1, a1_,a2_]=int0ee[1,0,0,0, 1,0,0,0, a1,a2]; (* 3/2*(20,3) *) int0ee[1,1,0,0, 0,0,1,1, a1_,a2_]=int0ee[1,0,0,0, 0,0,1,0, a1,a2]; int0ee[0,0,1,1, 1,1,0,0, a1_,a2_]=int0ee[1,0,0,0, 0,0,1,0, a1,a2]; (* 9/2*((22,22)-(24,24)) *) int0ee[1,0,1,1, 1,0,1,1, a1_,a2_]=1/(5*a1^4*a2^4*(a1 + a2)^7)* 144*(5*a1^6 + 35*a1^5*a2 + 103*a1^4*a2^2 + 161*a1^3*a2^3 + 103*a1^2*a2^4 + 35*a1*a2^5 + 5*a2^6); int0ee[1,1,1,0, 1,1,1,0, a1_,a2_]=int0ee[1,0,1,1, 1,0,1,1, a1,a2]; (* 9/2*((22,22)-(24,24)) *) int0ee[1,0,1,1, 1,1,1,0, a1_,a2_]=(432*(a1^2 + 7*a1*a2 + a2^2))/(5*a1^2*a2^2*(a1 + a2)^7); int0ee[1,1,1,0, 1,0,1,1, a1_,a2_]=int0ee[1,0,1,1, 1,1,1,0, a1,a2]; (* (*int0ee[l1a,l2a,l1b,l2b,m,a1,a2]= where m1a+m1b=m2a+m2b=m, A[n1a]+A[n1b]=a1, A[n2a]+A[n2b]=a2 *) int0ee[l1a_,l2a_,l1b_,l2b_,m_,a1_,a2_]:=0; (* (2,2) *) int0ee[0,0,0,0,0,a1_,a2_]=2*(a1^2 + 3*a1*a2 + a2^2)/(a1^2*a2^2*(a1 + a2)^3); (* 3*(10,10) *) int0ee[1,0,1,0,0,a1_,a2_]=(12*(2*a1^4 + 10*a1^3*a2 + 10*a1^2*a2^2 + 5*a1*a2^3 + a2^4))/(a1^4*a2^2*(a1 + a2)^5); int0ee[0,1,0,1,0,a1_,a2_]=int0ee[1,0,1,0,0,a1,a2]; (* 3*(10,17) *) int0ee[1,0,0,1,0,a1_,a2_]=(8*(a1^2 + 5*a1*a2 + a2^2))/(a1^2*a2^2*(a1 + a2)^5); int0ee[0,1,1,0,0,a1_,a2_]=int0ee[1,0,0,1,0,a1,a2]; (* 3/2*((8,8)-(10,10)) *) int0ee[1,0,1,0,1,a1_,a2_]=int0ee[1,0,1,0,0,a1,a2]; int0ee[0,1,0,1,1,a1_,a2_]=int0ee[1,0,1,0,1,a1,a2]; (* 3/2*(20,3) *) int0ee[1,0,0,1,1,a1_,a2_]=int0ee[1,0,0,1,0,a1,a2]; int0ee[0,1,1,0,1,a1_,a2_]=int0ee[1,0,0,1,1,a1,a2]; (* 9/2*((22,22)-(24,24)) *) int0ee[1,1,1,1,0,a1_,a2_]=1/(5*a1^4*a2^4*(a1 + a2)^7)* 144*(5*a1^6 + 35*a1^5*a2 + 103*a1^4*a2^2 + 161*a1^3*a2^3 + 103*a1^2*a2^4 + 35*a1*a2^5 + 5*a2^6); (* 9/2*((22,22)-(24,24)) *) int0ee[1,1,1,1,1,a1_,a2_]=;*) i1a[a_]=2/a^3; (* ic1a[a]=Integral[Exp[-a r] r^-1 dr]/(4 Pi) *) ic1a[a_]=1/a^2; (* it1ab[a,b]=Integral[Exp[-a r] (-1/(2 r^2) D[r^2 D[Exp[-b r],r],r]) dr]/(4 Pi) *) (*psi=Exp[-b r]; Integrate[r^2 Exp[-a r] (-1/(2 r^2) D[r^2 D[psi,r],r]),{r,0,Infinity}]*) it1ab[a_,b_]=a b/(a+b)^3; (* ic2a[a1,a2]=Integral[Exp[-a1 r1] Exp[-a2 r2] r12^(-1) dr]/(16 Pi^2) *) ic2a[a1_,a2_]=2*(a1^2 + 3*a1*a2 + a2^2)/(a1^2*a2^2*(a1 + a2)^3); iNa[a_]=Product[i1a[a[n]],{n,npart}]; icNa[a_]=(iNa[a] Sum[ic1a[a[n]]/i1a[a[n]],{n,npart}])//Simplify; itNab[a_,b_]=(iNa[(a[#]+b[#])&] Sum[it1ab[a[n],b[n]]/i1a[a[n]+b[n]],{n,npart}])//Simplify; ieeNa[a_]=(iNa[a] Sum[ic2a[a[n1],a[n2]]/i1a[a[n1]]/i1a[a[n2]], {n1,npart-1},{n2,n1+1,npart}])//Simplify; (*f[l_,m_,n_,theta_,phi_,r_]=SphericalHarmonicY[l,m,theta,phi] r^l Exp[-A[n] r]; f[l_,m_,n_]=Array[f[l,m,n,theta[#],phi[#],r[#]]&,npart]; fs[n_]=Array[f[0,0,n,theta[#],phi[#],r[#]]&,npart]; fp[n_]=Array[f[1,1,n,theta[#],phi[#],r[#]]&,npart]; *) f[l_,m_,n_]=Array[f[l,m,n,x[#]]&,npart]; fs[n_]=Array[f[0,0,n,x[#]]&,npart]; fp[n_]=Array[f[1,1,n,x[#]]&,npart]; alp=Array[alph,npart]; bet=Array[beta,npart]; psi=Det[{alp fs[1],bet fp[2]}]- Det[{bet fs[1],alp fp[2]}]; spins=Flatten[Table[{alph[n],beta[n]},{n,npart}]]; psi=Flatten[CoefficientList[psi,spins]]; psi1=Complement[psi,{0}]; (* Test that elements of psi1 enter psi equal number of times *) If[Length[Union[Map[Count[psi,#]&,psi1]]]==1,psi=psi1, Print["Warning: Something strange with wavefunction!"]]; If[Length[Complement[Map[Count[psi1,-#]&,psi1],{0}]]>0, Print["Further removing of elements with an opposite sign is possible (but not done)."]]; mpsi=Length[psi]; rs=Array[r,npart]; {mean,cmean,eemean}= (* , , *) Sum[meansum=psi[[m]]^2//Expand; lm=Length[meansum]; Sum[meanl=meansum[[l]]; {cl,el}=meanl/.c_. E^r_->{c,r}; el=Collect[el,rs]; alist=-Operate[List&,el]/.r[n_]->1; Do[al[n]=alist[[n]],{n,npart}]; cl {iNa[al],icNa[al],ieeNa[al]},{l,lm}], {m,mpsi}]; (* *) tmean= Sum[meansum=psi[[m]]//Expand; lm=Length[meansum]; Sum[meanl=meansum[[l1]]; {cl1,el}=meanl/.c_. E^r_->{c,r}; alist=-Operate[List&,el]/.r[n_]->1; Do[al[n]=alist[[n]],{n,npart}]; cl1 Sum[meanl=meansum[[l2]]; {cl2,el}=meanl/.c_. E^r_->{c,r}; alist=-Operate[List&,el]/.r[n_]->1; Do[bl[n]=alist[[n]],{n,npart}]; cl2 itNab[al,bl] ,{l2,lm}] ,{l1,lm}] ,{m,mpsi}]; W=(tmean-cmean+lambda eemean)/mean; (* Table of energy as a function of Z *) Clear[en,den]; <ToString[npart]<>".dat"; os=OpenWrite[file]; SetOptions[os,PageWidth->Infinity]; Write[os,"Variational energy of ions with npart, dependence on Z"]; Save[os,npart]; Do[lambda=1/z; z0=Round[100 z]/100; Do[A[n]=a[n][npart,z0],{n,npart}]; w=W; {en[npart,z],den[npart,z]}={z^2 w,2 z w-eemean/mean} ,{z,1/100,5,1/100}]; Save[os,en,den]; Close[os];