(* Calculating of the critical charge for helium isoelectronic series by finding a stationary point of the variational functional *) (* With advanced differentiation over a, b, c *) (* With optimized performance *) (********** INPUT STARTS **********) {nmin,nmax}={7,8}; choice=0; (* 0 for i+j+k<=N, 1 for i+j^2+k^2<=N *) ndig=96; (* Precision of calculations *) {astart,bstart,cstart}={0.35,1.03,0}; (* guess *) eps=0.03; (* limit of guessing *) (********** INPUT ENDS **********) ndig1=Min[ndig-10,22]; (* AccuracyGoal in FindRoot *) ndim=3; (* dimensionality - 3 or 5 or 7 ...*) miter=100; (* max n. of iterations *) x=" "; x3=" "; {$MinPrecision,$MaxPrecision}={0,5000}; (* To avoid errors at running the following line *) $MinPrecision=$MaxPrecision=ndig; Off[Precision::"precsm"]; Off[Solve::"precsm"]; Clear[t]; chop[t_]=If[t==0,0,t]; {astart,bstart,cstart,eps}=SetPrecision[{astart,bstart,cstart,eps},ndig]; Print["{Nrange,choice,precision} = ",{{nmin,nmax},If[choice==0," i+j+k<=N "," i+j^2+k^2<=N "],ndig}]; p=ToFileName["..","packages"]; If[!MemberQ[$Path,p],$Path=Append[$Path,p]]; <int1[i,j,k][a,b,c], hylleraas`private`int2[i_,j_,k_]->int2[i,j,k][a,b,c], hylleraas`private`ints[i_,j_,k_]->ints[i,j,k][a,b,c] }; {Ta,V0a,V1a,Sa}=D[{T,V0,V1,S},a]/.{ Derivative[1,0,0][int1[i_,j_,k_]][a,b,c]->int1a[i,j,k], Derivative[1,0,0][int2[i_,j_,k_]][a,b,c]->int2a[i,j,k], Derivative[1,0,0][ints[i_,j_,k_]][a,b,c]->intsa[i,j,k], int1[i_,j_,k_][a,b,c]->int1[i,j,k], int2[i_,j_,k_][a,b,c]->int2[i,j,k], ints[i_,j_,k_][a,b,c]->ints[i,j,k] }; {Tb,V0b,V1b,Sb}=D[{T,V0,V1,S},b]/.{ Derivative[0,1,0][int1[i_,j_,k_]][a,b,c]->int1b[i,j,k], Derivative[0,1,0][int2[i_,j_,k_]][a,b,c]->int2b[i,j,k], Derivative[0,1,0][ints[i_,j_,k_]][a,b,c]->intsb[i,j,k], int1[i_,j_,k_][a,b,c]->int1[i,j,k], int2[i_,j_,k_][a,b,c]->int2[i,j,k], ints[i_,j_,k_][a,b,c]->ints[i,j,k] }; {Tc,V0c,V1c,Sc}=D[{T,V0,V1,S},c]/.{ Derivative[0,0,1][int1[i_,j_,k_]][a,b,c]->int1c[i,j,k], Derivative[0,0,1][int2[i_,j_,k_]][a,b,c]->int2c[i,j,k], Derivative[0,0,1][ints[i_,j_,k_]][a,b,c]->intsc[i,j,k], int1[i_,j_,k_][a,b,c]->int1[i,j,k], int2[i_,j_,k_][a,b,c]->int2[i,j,k], ints[i_,j_,k_][a,b,c]->ints[i,j,k] }; {T,V0,V1,S}={T,V0,V1,S}/.{ int1[i_,j_,k_][a,b,c]->int1[i,j,k], int2[i_,j_,k_][a,b,c]->int2[i,j,k], ints[i_,j_,k_][a,b,c]->ints[i,j,k] }; Do[ (* nm = nmin, ..., nmax *) Print[nm//Definition]; outf="lambopt.dat"; Clear[ind,mmax]; mm=hmapping[nm,choice,ind,mmax]; Print[Definition[mm]]; L=M=Table[Null,{mm},{mm}]; (* *) Clear[a,b,c]; time0=Timing[ hintegrals[a,b,c,ndim,nm]; (* Diffirentiating integrals *) nmax=2 nm; nmadd=2 (ndim-3); nmax3=nmax+nmadd+3; Do[ Do[namax=n-nc;namax1=Floor[namax/2]; Do[nb=namax-na; integ1[na-1,nb-1,nc-1]=integ2[nb-1,na-1,nc-1]=hylleraas`private`int1[na-1,nb-1,nc-1]; integ2[na-1,nb-1,nc-1]=integ1[nb-1,na-1,nc-1]=hylleraas`private`int2[na-1,nb-1,nc-1]; integs[na-1,nb-1,nc-1]=integs[nb-1,na-1,nc-1]=hylleraas`private`ints[na-1,nb-1,nc-1]; integ1a[na-1,nb-1,nc-1]=integ2a[nb-1,na-1,nc-1]= D[hylleraas`private`int1[na-1,nb-1,nc-1],a]//Together; integ2a[na-1,nb-1,nc-1]=integ1a[nb-1,na-1,nc-1]= D[hylleraas`private`int2[na-1,nb-1,nc-1],a]//Together; integsa[na-1,nb-1,nc-1]=integsa[nb-1,na-1,nc-1]= D[hylleraas`private`ints[na-1,nb-1,nc-1],a]//Together; integ1b[na-1,nb-1,nc-1]=integ2b[nb-1,na-1,nc-1]= D[hylleraas`private`int1[na-1,nb-1,nc-1],b]//Together; integ2b[na-1,nb-1,nc-1]=integ1b[nb-1,na-1,nc-1]= D[hylleraas`private`int2[na-1,nb-1,nc-1],b]//Together; integsb[na-1,nb-1,nc-1]=integsb[nb-1,na-1,nc-1]= D[hylleraas`private`ints[na-1,nb-1,nc-1],b]//Together; integ1c[na-1,nb-1,nc-1]=integ2c[nb-1,na-1,nc-1]= D[hylleraas`private`int1[na-1,nb-1,nc-1],c]//Together; integ2c[na-1,nb-1,nc-1]=integ1c[nb-1,na-1,nc-1]= D[hylleraas`private`int2[na-1,nb-1,nc-1],c]//Together; integsc[na-1,nb-1,nc-1]=integsc[nb-1,na-1,nc-1]= D[hylleraas`private`ints[na-1,nb-1,nc-1],c]//Together; ,{na,0,namax1}],{nc,0,n}],{n,0,nmax3}]; ]//First;Print["time-matrix-general form = ",time0]; (* *) guess={{atrial,{astart,astart+eps}}, {btrial,{bstart,bstart+eps}}, {ctrial,{cstart,cstart+eps}}}; (******************************* START find-root ******************************) dampf=1; iter=0; abciter={}; rt=FindRoot[ iter++; {a,b,c}= SetPrecision[{atrial,btrial,ctrial},ndig]//N; (* increase of precision *) (* finding the eigenvalue *) time1=Timing[ Clear[int1,int2,ints,int1a,int2a,intsa,int1b,int2b,intsb,int1c,int2c,intsc,na,nb,nc]; int1[na_,nb_,nc_]:=int1[na,nb,nc]=int2[nb,na,nc]=integ1[na,nb,nc]; int2[na_,nb_,nc_]:=int2[na,nb,nc]=int1[nb,na,nc]=integ2[na,nb,nc]; ints[na_,nb_,nc_]:=ints[na,nb,nc]=ints[nb,na,nc]=integs[na,nb,nc]; int1a[na_,nb_,nc_]:=int1a[na,nb,nc]=int2a[nb,na,nc]=integ1a[na,nb,nc]; int2a[na_,nb_,nc_]:=int2a[na,nb,nc]=int1a[nb,na,nc]=integ2a[na,nb,nc]; intsa[na_,nb_,nc_]:=intsa[na,nb,nc]=intsa[nb,na,nc]=integsa[na,nb,nc]; int1b[na_,nb_,nc_]:=int1b[na,nb,nc]=int2b[nb,na,nc]=integ1b[na,nb,nc]; int2b[na_,nb_,nc_]:=int2b[na,nb,nc]=int1b[nb,na,nc]=integ2b[na,nb,nc]; intsb[na_,nb_,nc_]:=intsb[na,nb,nc]=intsb[nb,na,nc]=integsb[na,nb,nc]; int1c[na_,nb_,nc_]:=int1c[na,nb,nc]=int2c[nb,na,nc]=integ1c[na,nb,nc]; int2c[na_,nb_,nc_]:=int2c[na,nb,nc]=int1c[nb,na,nc]=integ2c[na,nb,nc]; intsc[na_,nb_,nc_]:=intsc[na,nb,nc]=intsc[nb,na,nc]=integsc[na,nb,nc]; Do[{i2,j2,k2}=ind[m2]; Do[{i1,j1,k1}=ind[m1]; L[[m1,m2]]=L[[m2,m1]]=T-V0+1/2 S; M[[m1,m2]]=M[[m2,m1]]=V1; ,{m1,m2,mm}],{m2,mm}]; A=Inverse[M].L; A=Chop[A,10.^(-ndig/2)]; (* necessary to avoid overflow *) ]//First;(*Print["time-matrix = ",time1];*) time2=Timing[ EA=Re[Transpose[Eigensystem[A]]]; {val0,vec0}=Sort[EA,#1[[1]]<#2[[1]]&][[1]]; ]//First;(*Print["time-eigenvalues = ",time2];*) (* finding derivative of the eigenvalue d/da *) time3=Timing[ Do[{i2,j2,k2}=ind[m2]; Do[{i1,j1,k1}=ind[m1]; L[[m1,m2]]=L[[m2,m1]]=Ta-V0a+1/2 Sa; M[[m1,m2]]=M[[m2,m1]]=V1a; ,{m1,m2,mm}],{m2,mm}]; vala=(vec0.L.vec0)/(vec0.M.vec0)-val0; (* finding derivative of the eigenvalue d/db *) Do[{i2,j2,k2}=ind[m2]; Do[{i1,j1,k1}=ind[m1]; L[[m1,m2]]=L[[m2,m1]]=Tb-V0b+1/2 Sb; M[[m1,m2]]=M[[m2,m1]]=V1b; ,{m1,m2,mm}],{m2,mm}]; valb=(vec0.L.vec0)/(vec0.M.vec0)-val0; (* finding derivative of the eigenvalue d/dc *) Do[{i2,j2,k2}=ind[m2]; Do[{i1,j1,k1}=ind[m1]; L[[m1,m2]]=L[[m2,m1]]=Tc-V0c+1/2 Sc; M[[m1,m2]]=M[[m2,m1]]=V1c; ,{m1,m2,mm}],{m2,mm}]; valc=(vec0.L.vec0)/(vec0.M.vec0)-val0; ]//First; delt={vala,valb,valc}; deltabs=Plus@@Abs/@delt; deltpr=PaddedForm[deltabs//FortranForm,{3,1},ExponentFunction->Identity]; abciter=Append[abciter,{deltabs,{a,b,c},val0}]; If[iter==1,Print["Time (one step): matrix, eigenv., derivatives: ",time1,", ",time2,", ",time3]]; Print[iter," --- ",deltpr," ",{printF[a,9,6],printF[b,9,6],printF[c,9,6]} ];delt=={0,0,0}, Sequence@@guess//Release, WorkingPrecision->ndig,AccuracyGoal->ndig1,MaxIterations->miter,DampingFactor->dampf]; (******************************* END find-root ******************************) {delt,abc0,eig0}=Sort[abciter,#1[[1]]<#2[[1]]&][[1]]; {apr,bpr,cpr}={printF[abc0[[1]],ndig1+2,ndig1-2], printF[abc0[[2]],ndig1+2,ndig1-2], printF[abc0[[3]],ndig1+2,ndig1-2]}; lambpr=printF[-eig0,ndig1+4,ndig1]; zpr=printF[-1/eig0,ndig1+4,ndig1]; deltpr=PaddedForm[delt//FortranForm,{3,1},ExponentFunction->Identity]; Print["Non-linear parameters: ",{apr,bpr,cpr}]; Print["Lambda: ",lambpr]; Print["L.h.s. of the equations (should be zero), n. of iter.: ",deltpr,", ",iter]; If[iter>=miter, Print["Number of iterations exceeded the limit of ",miter," iterations."]; Break[] ]; os=OpenAppend[outf,PageWidth->Infinity,FormatType->OutputForm]; Write[os,nm,x3,apr,x,bpr,x,cpr,x3,lambpr,x,zpr,x3,deltpr]; Close[os]; ,{nm,nmin,nmax}]; If[$MachineType!="PC",Exit[]];