c ao.f c generate density dot diagrams for atomic orbitals. Prepare output that c Midas can read and plot. CHARACTER*3 ANAME(1),RESNAME(1) CHARACTER*6 SPECI(1) CHARACTER*4 DUMMY CHARACTER*40 UFILE,NFILE,HEAD INTEGER COUNT,INUMBER(1),IRESNUM(1) REAL*4 X,Y,Z double precision dseed ANAME(1)='O ' RESNAME(1)='HOH' SPECI(1)='HETATM' inumber(1)=1 iresnum(1)=1 psi2max=0.0 call system('clear') write(6,*)'Which atomic orbital do you want to generate?' write(6,*)'1 1s' write(6,*)'2 2s' write(6,*)'3 2px' write(6,*)'4 2py' write(6,*)'5 3s' write(6,*)'6 3px' write(6,*)'7 3pz' write(6,*)'8 3dz2' write(6,*)'9 3dxz' write(6,*)'10 3dyz' write(6,*)'11 3dx2-y2' write(6,*)'12 3dxy' read(5,*)icon write(6,*)'What is the atomic number Z?' read(5,*)Za write(6,*)'What weighting coefficient (0 to 1); for LCAO?' read(5,*)Weight WRITE(6,*)'Give a name for the output file:' READ(5,111)NFILE 111 format(a) call parse(NFILE,'hdr',HEAD) OPEN(1,FILE=HEAD) y=0. z=0. x=0. xtra1=0. write(1,101)x,y,z,xtra1,xtra1 101 format('ATOM 1 N GLY 1 ',3f8.3,2f6.2) close(1) OPEN(UNIT=2,FILE=NFILE) write(2,110) 110 format('GLY 1 N 0.000 0.000 0.000 A') dseed=1.0 if(icon.eq.1)width=3.2/Za !1s if(icon.ge.2.and.icon.le.4)width=6.4/Za !n=2 if(icon.ge.5.and.icon.le.12)width=12./Za !n=3 raddeg=180./3.14159 xtra1=0. xtra2=0. ndots=Weight*3000 DO 5 I=1,ndots c Randomly generate a point in cartesian space in a cube of large width. 2 call ggubs(dseed,1,x) call ggubs(dseed,1,y) call ggubs(dseed,1,z) x=(2.*x-1.)*width y=(2.*y-1.)*width z=(2.*z-1.)*width c Convert this point to spherical polar coordinates. r=sqrt(x*x+y*y+z*z) if(r.eq.0.0)r=0.0001 theta=acos(z/r) if(x.eq.0.00)x=0.0001 phi=atan(y/x) cdebug write(6,*)r,theta,phi,theta*raddeg,phi*raddeg c Find the value of the wave function squared at this point if(icon.eq.1)call orb_1s(Za,r,psi,psi2,psimax,psi2max) if(icon.eq.2)call orb_2s(Za,r,psi,psi2,psimax,psi2max) if(icon.eq.3) & call orb_2px(Za,r,theta,phi,psi,psi2,psimax,psi2max) if(icon.eq.4) & call orb_2py(Za,r,theta,phi,psi,psi2,psimax,psi2max) if(icon.eq.5)call orb_3s(Za,r,psi,psi2,psimax,psi2max) if(icon.eq.6) & call orb_3px(Za,r,theta,phi,psi,psi2,psimax,psi2max) if(icon.eq.7) & call orb_3pz(Za,r,theta,phi,psi,psi2,psimax,psi2max) if(icon.eq.8) & call orb_3dz2(Za,r,theta,phi,psi,psi2,psimax,psi2max) if(icon.eq.9) & call orb_3dxz(Za,r,theta,phi,psi,psi2,psimax,psi2max) if(icon.eq.10) & call orb_3dyz(Za,r,theta,phi,psi,psi2,psimax,psi2max) if(icon.eq.11) & call orb_3dx2y2(Za,r,theta,phi,psi,psi2,psimax,psi2max) if(icon.eq.12) & call orb_3dxy(Za,r,theta,phi,psi,psi2,psimax,psi2max) cdebug if(I.le.10)write(6,*)r,psi,psi2,psimax,psi2max cdebug if(I.le.10)write(6,*)r,psi,psi2,psimax,psi2max c Test to see if a density dot will be placed at this point. If so, write c this point out in PDB format for Midas to read. f=psi2/psi2max call ggubs(dseed,1,random) if(random.gt.f)goto 2 write (2,130) x,y,z 130 format('GLY 1 O',F8.3,1X,F8.3,1X,F8.3,1X,'SR0 0.185') 5 CONTINUE STOP END subroutine orb_1s(Za,r,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms factor=(1./sqrt(3.14159))*(za/a)**1.5 psimax=factor psi2max=psimax*psimax alpha=za/a ifirst=1 end if psi=factor*exp(-alpha*r) psi2=psi*psi return end subroutine orb_2s(Za,r,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(2.*a) factor=(1./(4.*sqrt(2.*3.14159)))*(za/a)**1.5 c search for maximum: rmax=0. psimax=factor*2. psi2max=psimax*psimax write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(2.-Za*r/a)*exp(-alpha*r) psi2=psi*psi return end subroutine orb_2px(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/a/2. factor=(1./4.0/sqrt(2.*3.14159))*(za/a)**2.5 rmax=1./alpha psimax=factor*rmax*exp(-alpha*rmax) psi2max=psimax*psimax ifirst=1 end if psi=factor*r*exp(-alpha*r)*sin(theta)*cos(phi) psi2=psi*psi return end subroutine orb_2py(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/a/2. factor=(1./4.0/sqrt(2.*3.14159))*(za/a)**2.5 rmax=1./alpha psimax=factor*rmax*exp(-alpha*rmax) psi2max=psimax*psimax ifirst=1 end if psi=factor*r*exp(-alpha*r)*sin(theta)*sin(phi) psi2=psi*psi return end subroutine orb_3s(Za,r,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(3.*a) factor=(1./(18.*sqrt(3.*3.14159)))*(za/a)**1.5 c search for maximum: rmax=0. psimax=factor*6. psi2max=psimax*psimax write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(6.-4.*Za*r/a+(2.*Za*r/(3.*a))**2.)*exp(-alpha*r) psi2=psi*psi return end subroutine orb_3px(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(3.*a) factor=(sqrt(2.)/(81.*sqrt(3.14159)))*(za/a)**1.5 c search for maximum: do rtry=0.01,40.,0.01 psi=factor*(6.*Za*rtry/a-(za*rtry/a)**2) & *exp(-alpha*rtry) psi2=psi*psi if(psi2.gt.psi2max)then psi2max=psi2 psimax=psi rmax=rtry end if end do write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(6.*Za*r/a-(za*r/a)**2)*exp(-alpha*r)* & sin(theta)*cos(phi) psi2=psi*psi return end subroutine orb_3pz(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(3.*a) factor=(sqrt(2.)/(81.*sqrt(3.14159)))*(za/a)**1.5 c search for maximum: do rtry=0.01,40.,0.01 psi=factor*(6.*Za*rtry/a-(za*rtry/a)**2) & *exp(-alpha*rtry) psi2=psi*psi if(psi2.gt.psi2max)then psi2max=psi2 psimax=psi rmax=rtry end if end do write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(6.*Za*r/a-(za*r/a)**2)*exp(-alpha*r)* & cos(theta) psi2=psi*psi return end subroutine orb_3dz2(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(3.*a) factor=(1./(81.*sqrt(6.*3.14159)))*(za/a)**1.5 c search for maximum: do rtry=0.01,40.,0.01 psi=factor*(Za*rtry/a)**2 & *exp(-alpha*rtry)*3.0 psi2=psi*psi if(psi2.gt.psi2max)then psi2max=psi2 psimax=psi rmax=rtry end if end do write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(Za*r/a)**2*exp(-alpha*r)* & (3.*cos(theta)**2-1.) psi2=psi*psi return end subroutine orb_3dxz(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(3.*a) factor=(sqrt(2.)/(81.*sqrt(3.14159)))*(za/a)**1.5 c search for maximum: do rtry=0.01,40.,0.01 psi=factor*(Za*rtry/a)**2 & *exp(-alpha*rtry) psi2=psi*psi if(psi2.gt.psi2max)then psi2max=psi2 psimax=psi rmax=rtry end if end do write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(Za*r/a)**2*exp(-alpha*r)* & sin(theta)*cos(theta)*cos(phi) psi2=psi*psi return end subroutine orb_3dyz(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(3.*a) factor=(sqrt(2.)/(81.*sqrt(3.14159)))*(za/a)**1.5 c search for maximum: do rtry=0.01,40.,0.01 psi=factor*(Za*rtry/a)**2 & *exp(-alpha*rtry) psi2=psi*psi if(psi2.gt.psi2max)then psi2max=psi2 psimax=psi rmax=rtry end if end do write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(Za*r/a)**2*exp(-alpha*r)* & sin(theta)*cos(theta)*sin(phi) psi2=psi*psi return end subroutine orb_3dx2y2(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(3.*a) factor=(1./(81.*sqrt(2.*3.14159)))*(za/a)**1.5 c search for maximum: do rtry=0.01,40.,0.01 psi=factor*(Za*rtry/a)**2 & *exp(-alpha*rtry) psi2=psi*psi if(psi2.gt.psi2max)then psi2max=psi2 psimax=psi rmax=rtry end if end do write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(Za*r/a)**2*exp(-alpha*r)* & sin(theta)**2*cos(2.*phi) psi2=psi*psi return end subroutine orb_3dxy(Za,r,theta,phi,psi,psi2,psimax,psi2max) static ifirst,factor,a,alpha data ifirst/0/ if(ifirst.eq.0)then a=0.529 !Bohr radius in Angstroms alpha=za/(3.*a) factor=(1./(81.*sqrt(2.*3.14159)))*(za/a)**1.5 c search for maximum: do rtry=0.01,40.,0.01 psi=factor*(Za*rtry/a)**2 & *exp(-alpha*rtry) psi2=psi*psi if(psi2.gt.psi2max)then psi2max=psi2 psimax=psi rmax=rtry end if end do write(6,*)'Psi and Psi2 = ',psimax,psi2max,' at r = ',rmax ifirst=1 end if psi=factor*(Za*r/a)**2*exp(-alpha*r)* & sin(theta)**2*sin(2.*phi) psi2=psi*psi return end subroutine ggubs(dseed,n,unif) c This is written by Scott H. Northrup and seems to work for the IRIS 4D/70. c When the IMSL routine of the same name is available, this routine should be c taken out. c gererates uniformly distributed random numbers. real*4 unif double precision u1,rand,dseed u1=rand() unif=u1 return end subroutine parse (name,newext,newname) character*40 name,newname character*3 newext integer letter,position do letter=40,1,-1 if (name(letter:letter).ne.' ') then if (name(letter:letter).eq.'.') then position=letter goto 200 endif endif end do 200 continue newname=name(1:position) // newext return end