program contour dimension ia(51,84),hol(20),ktor(10),xaa(51,84) data hol(1)/4h 1s/,hol(2)/4h 2s/,hol(3)/4h 2pz/,hol(4)/4h 2px/, 1hol(5)/4h 2py/,hol(6)/4h 3s/,hol(7)/4h 3pz/,hol(8)/4h 3px/, 2hol(9)/4h 3py/,hol(10)/4h3dz2/,hol(11)/4h3dxz/,hol(12)/4h3dyz/, 3hol(13)/4hdxmy/ ,hol(14)/4h3dxy/ data ktor(1)/1h /,ktor(2)/1h./,ktor(3)/1h,/,ktor(4)/1h-/,ktor(5) 1/1h+/,ktor(6)/1h=/,ktor(7)/1h*/,ktor(8)/1ho/,ktor(9)/1hz/,ktor(10) 2/1hx/ integer outop a0=0.0529 rtpi=1.7725 aais=sqrt(1/(a0*a0*a0))/rtpi open(23,file='contour.out') 990 continue print*,'Give iop, outop, nz,wdth,z' read(5,*)iop,outop,nz,wdth,z print*,iop,outop,nz,z if(iop.eq.0)stop if(outop.eq.0)write(2,1073) if(outop.eq.1)write(2,1074) if(outop.eq.2)write(2,2075) 2075 format(1h1,1x, 77(1h*),/,30x,'electron density map') 1073 format(1h1,1x, 77(1h*),/,30x,'electron density maps relative to', 1' psi is at r=0') 1074 format(1h1,1x, 77(1h*),/,30x,'map of wave function relative to', 1' psi is at r=0') 1020 format(1x,'option=',i2,2x,'nuclear charge=+',i1,2x,'plot width' 1,'/nm = ',1f10.4) write(2,1020)iop,nz,wdth 1075 format(1x, 77(1h*)) xx=sqrt(float(nz*nz*nz)/(a0*a0*a0)) nup=25 if(outop.eq.2)nup=75 do 2 j=1,46 y=(23.5-j)*wdth/45.0 do 2 k=1,nup if(outop.eq.2)x=(k-38.0)*wdth/74.0 if(outop.ne.2)x=(k-13.0)*wdth/24.0 r=sqrt(x*x+y*y+z*z) xxz=nz*r/a0 aa=0.0 if(r.lt.0.00001)go to 701 cost=z/r sint=sqrt(1.-cost**2) if(iop.eq.3)aa=0.17677*xx/rtpi*cost*exp(-xxz/2)*xxz if(iop.eq.7)aa=0.017459*xx/rtpi*(6.0-xxz)*xxz*cost*exp(-xxz/3) if(iop.eq.10)aa=0.005040*xx/rtpi*xxz*xxz*(3*cost*cost-1)*exp(-xxz/ 13) if(sint.lt.0.00001)go to 701 d=r*sint cosp=x/d sinp=y/d if(iop.eq.4)aa=0.176776*xx/rtpi*cosp*sint*xxz*exp(-xxz/2) if(iop.eq.5)aa=0.176776*xx/rtpi*sinp*sint*xxz*exp(-xxz/2) if(iop.eq.8)aa=0.017459*xx/rtpi*(6.0-xxz)*xxz*sint*cosp*exp(-xxz/ 13) if(iop.eq.9)aa=0.017459*xx/rtpi*(6.0-xxz)*xxz*sint*sinp*exp(-xxz/ 13) if(iop.ge.11)afac=0.017459*xx/rtpi*xxz*xxz*exp(-xxz/3) if(iop.eq.11)aa=afac*sint*cost*cosp if(iop.eq.12)aa=afac*sint*cost*sinp if(iop.eq.13)aa=afac*sint*sint*(cosp*cosp-sinp*sinp)/2 if(iop.eq.14)aa=afac*sint*sint*sinp*cosp 701 if(iop.eq.1)aa=xx*exp(-xxz)/rtpi if(iop.eq.2)aa=xx/rtpi*(2.0-xxz)*exp(-xxz/2)*0.176776 if(iop.eq.6)aa=0.0071127*xx/rtpi*(27.0-18*xxz+2*xxz*xxz)*exp(-xxz/ 1 3) if(outop.eq.0.or.outop.eq.2)xaa(j,k)=aa**2 if(outop.eq.1)xaa(j,k)=aa write(23,122) x,y,xaa(j,k) 2 continue c 1 continue 122 format(1x,(3f10.5,2x)) zmax=0.0 do 60 j=1,46 do 60 k=1,nup 60 if(abs(xaa(j,k)).gt.zmax)zmax=abs(xaa(j,k)) print*,zmax,aais xfac=zmax/aais if(outop.eq.0)xfac=zmax/aais**2 xfac=1.0/xfac write(2,1003)z 1003 format(1x,'in the plane z=',1f10.4,'nm',/) if(outop.eq.2)go to 497 if(outop.eq.0)go to 631 write(2,1004)hol(iop),xfac 1004 format(1x,'psi map for the ',1a4,'-orbital times ',f10.3,' re' 1'lative to psi is at r=0',/,1x,77(1h*)) go to 632 631 continue write(2,1008)hol(iop),xfac 1008 format(1x,'psi squared map for the ',1a4,'-orbital times',f10.3, 1' relative to psi**2 is at r=0',/,1x,77(1h*)) 632 continue do 62 j=1,46 do 62 k=1,25 ia(j,k)=nint(99*xaa(j,k)/zmax) 62 continue do 10 j=1,46 10 write(2,1006)(ia(j,k),k=1,25) 1006 format(1x,1h*,25i3,1h*) 1010 format(1x,77(1h*),/,1h1) write(2,1010) go to 61 497 continue write(2,3008)hol(iop) write(2,1075) 3008 format(1x,'nap for the ',1a4,'orbital') do 82 j=1,46 do 82 k=1,75 xjjk=(9.9*xaa(j,k)/zmax) jk=xjjk+1 82 ia(j,k)=ktor(jk) do 83 j=1,46 write(2,1740)(ia(j,k),k=1,75) 83 continue 1740 format(1x,1h*,75a1,1h*) write(2,1010) 61 continue go to 990 end