program checkNormalization c The purpose is to decide which orbital has the smallest energy, 3d or 4s, c in a simple model of the atom where filled orbitals are taken from hydrogen implicit none integer i1,i2,j,Nr,Nx,lpass parameter(Nr=50) parameter(Nx=16) double precision a0,R10,R20,R21,R30,R31,R32,R40,Rcut,Z, & auxJacr1,auxJacr2,r1,r2,x,intr1,intr2,intx,chargeDens,pi, & Q10,Q20,Q21,Q30,Q31,Q32,Q40,SquaredSphHarmonic,auxJacx,hx double precision Vscreened,hr1,hr2 external SquaredSphHarmonic pi=3.14159265d0 c Bohr radius in nanometers a0=0.053d0 c Maximum extent of the integration grid Rcut= /a0 c integration steps in r and x hr1=Rcut/dble(Nr+1) hx=2d0/dble(Nx+1) c Z is 19 but it can be set to 0 to test the code without the nucleus Z=19d0 c a file to save the screened potential open(unit=15,file='DATA.DIR/Vscreened.dat') c here is the choice between 3d or 4s (pass 2 or 0 argument) lpass=0 **************************************************************** c This block to compute the triple integral **************************************************************** c initialize all integrals to zero intr1=0d0 c Start integral in r1 do i1=1,Nr r1=dble(i1)*hr1 call HydLikeFunctions(r1,R10,R20,R21,R30,R31,R32,R40) if(i1.eq.1.or.i1.eq.Nr) then auxJacr1=r1**2*hr1*1.5d0 else auxJacr1=r1**2*hr1 endif c Start integral in r2, here choose interval 0 to r1 or infty. hr2=r1/dble(Nr+1) c hr2=hr1 intr2=0d0 do i2=1,Nr r2=dble(i2)*hr2 call HydLikeFunctions(r2,Q10,Q20,Q21,Q30,Q31,Q32,Q40) chargeDens=2d0*(Q10**2+Q20**2+3d0*Q21**2+Q30**2+3d0*Q31**2) if(i2.eq.1.or.i2.eq.Nr) then auxJacr2=r2**2*hr2*1.5d0 else auxJacr2=r2**2*hr2 endif c Start integral in cos(theta) intx=0d0 c We separate the first and last point of the loop due to the 1.5 c to avoid having to evaluate an "if" condition in the inner loop auxJacx=hx*chargeDens do j=2,Nx-1 x=-1d0+dble(j)*hx intx=intx+auxJacx/dsqrt(r1**2+r2**2-2d0*r1*r2*x) & *SquaredSphHarmonic(lpass,x) enddo c here are the border points to complete the trapezoidal rule auxJacx=auxJacx*1.5d0 x=-1d0+hx intx=intx+auxJacx/dsqrt(r1**2+r2**2-2d0*r1*r2*x) & *SquaredSphHarmonic(lpass,x) x=-1d0+dble(Nx)*hx intx=intx+auxJacx/dsqrt(r1**2+r2**2-2d0*r1*r2*x) & *SquaredSphHarmonic(lpass,x) c Can use this line to set the integral to zero and leave only the nuclear c charge for studying the pieces separately c intx=0d0 c Ended integral in cos(theta), continue with r' integral intr2=intr2+2d0*pi*intx*auxJacr2 enddo print*, i1 c Ended integral in r', continue with r integral c Choose between the two lines depending on 3d or 4s Vscreened=(intr2-Z/r1)*auxJacr1 if(lpass.eq.0) then intr1=intr1+Vscreened*R40**2 elseif(lpass.eq.2) then intr1=intr1+Vscreened*R32**2 else print*, 'only set for 3d and 4s' stop endif c this line saves to disk the effective potential write(15,*) r1, Vscreened/auxJacr1 c can also use the next line to save the effective charge c write(15,*) r1, Z-intr2*r1 enddo c Ended integral in r1. Print energy in Hartree and in eV to screen print*, intr1, intr1*27.2d0 close(15) ************************************************************* end double precision function SquaredSphHarmonic(l,costheta) implicit none integer l double precision costheta,pi pi=3.14159265d0 if(l.eq.0) then SquaredSphHarmonic=1d0/(4d0*pi) return endif if(l.eq.2) then SquaredSphHarmonic=5d0/(16d0*pi)*(3d0*costheta**2-1d0)**2 return endif print*, 'only set for 0,2' stop end