c mo.f c generate density dot diagrams for molecular 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 molecular orbital do you want to generate?' write(6,*)'1 sigma1s = { Ca*1s(a) + Cb*1s(b) }' write(6,*)'2 sigma*1s = { Ca*1s(a) - Cb*1s(b) }' write(6,*)'3 sigma2p = { Ca*2px(a) + Cb*2px(b) }' write(6,*)'4 sigma*2p = { Ca*2px(a) - Cb*2px(b) }' write(6,*)'5 pi2p = { Ca*2py(a) + Cb*2py(b) }' write(6,*)'6 pi2p* = { Ca*2py(a) - Cb*2py(b) }' write(6,*)'7 sigma1s2p= { Ca*1s(a) + Cb*2px(b) }' read(5,*)icon write(6,*) &'What are the effective atomic numbers Za and Zb,resp?' read(5,*)Za,Zb write(6,*) &'What are the weighting coefficient Ca and Cb for LCAO?' read(5,*)Ca,Cb write(6,*)' What is the internuclear distance in Angstroms?' read(5,*)Re WRITE(6,*)'Give a name for the output file:' READ(5,111)NFILE 111 format(a) call parse (NFILE,'hdr',HEAD) OPEN(UNIT=1,FILE=HEAD) y=0. z=0. x=0.-Re/2. xtra1=0. write(1,101)x,y,z,xtra1,xtra1 101 format('ATOM 1 N GLY 1 ',3f8.3,2f6.2) x=0.+Re/2. write(1,102)x,y,z,xtra1,xtra1 102 format('ATOM 2 C 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 Zee=Za if(Zb.lt.Za)Zee=Zb if(icon.eq.1.or.icon.eq.2)width=(Re/2.+3.2/Zee) !sigma-1s if(icon.ge.3.and.icon.le.7)width=(Re/2.+6.4/Zee) !sigma/pi-2p write(6,*)'Cube width = ',width raddeg=180./3.14159 xtra1=0. xtra2=0. Weight=1 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 Find the value of the wave function squared at this point if(icon.eq.1) &call mo_sig1s1s(1,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) if(icon.eq.2) &call mo_sig1s1s(-1,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) if(icon.eq.3) &call mo_sig2p2p(1,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) if(icon.eq.4) &call mo_sig2p2p(-1,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) if(icon.eq.5) &call mo_pi2p2p(-1,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) if(icon.eq.6) &call mo_pi2p2p(+1,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) if(icon.eq.7) &call mo_sig1s2p(+1,Za,Zb,Ca,Cb,Re,x,y,z,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 mo_sig1s1s &(isn,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) static ifrst data ifrst/0/ c write(6,*)'Za,Zb,Ca,Cb,Re,x,y,z' c write(6,*)Za,Zb,Ca,Cb,Re,x,y,z c Convert this point to spherical polar coordinates in each atomic ref frame x=x+Re/2. ra=sqrt(x*x+y*y+z*z) if(ra.eq.0.0)ra=0.0001 c thetaa=acos(z/ra) c if(x.eq.0.00)x=0.0001 c phia=atan(y/x) call orb_1s(Za,ra,psia,psi2a,psimaxa,psi2maxa) psia=Ca*psia c write(6,*)'ra,psia=',ra,psia x=x-Re rb=sqrt(x*x+y*y+z*z) if(rb.eq.0.0)rb=0.0001 c thetab=acos(z/rb) c if(x.eq.0.00)x=0.0001 c phib=atan(y/x) call orb_1s(Zb,rb,psib,psi2b,psimaxb,psi2maxb) psib=Cb*psib c write(6,*)'rb,psib=',rb,psib psi=Ca*psia+isn*Cb*psib psi2=psi*psi x=x+Re/2. if(ifrst.eq.0)then ifrst=1 psimax=Ca*psimaxa+Cb*psimaxa psi2max=psimax*psimax c write(6,*)'psimaxa,psimaxb,psimax,psi2max' c write(6,*)psimaxa,psimaxb,psimax,psi2max end if return end subroutine mo_sig2p2p &(isn,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) static ifrst data ifrst/0/ c write(6,*)'Za,Zb,Ca,Cb,Re,x,y,z' c write(6,*)Za,Zb,Ca,Cb,Re,x,y,z c Convert this point to spherical polar coordinates in each atomic ref frame x=x+Re/2. ra=sqrt(x*x+y*y+z*z) if(ra.eq.0.0)ra=0.0001 thetaa=acos(z/ra) if(x.eq.0.00)x=0.0001 phia=atan(y/x) call orb_2px(Za,ra,thetaa,phia,psia,psi2a,psimaxa,psi2maxa) psia=Ca*psia c write(6,*)'ra,psia=',ra,psia x=x-Re rb=sqrt(x*x+y*y+z*z) if(rb.eq.0.0)rb=0.0001 thetab=acos(z/rb) if(x.eq.0.00)x=0.0001 phib=atan(y/x) call orb_2px(Zb,rb,thetab,phib,psib,psi2b,psimaxb,psi2maxb) psib=Cb*psib c write(6,*)'rb,psib=',rb,psib psi=Ca*psia+isn*Cb*psib psi2=psi*psi x=x+Re/2. if(ifrst.eq.0)then ifrst=1 psimax=Ca*psimaxa+Cb*psimaxa psi2max=psimax*psimax c write(6,*)'psimaxa,psimaxb,psimax,psi2max' c write(6,*)psimaxa,psimaxb,psimax,psi2max end if return end subroutine mo_pi2p2p &(isn,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) static ifrst data ifrst/0/ c write(6,*)'Za,Zb,Ca,Cb,Re,x,y,z' c write(6,*)Za,Zb,Ca,Cb,Re,x,y,z c Convert this point to spherical polar coordinates in each atomic ref frame x=x+Re/2. ra=sqrt(x*x+y*y+z*z) if(ra.eq.0.0)ra=0.0001 thetaa=acos(z/ra) if(x.eq.0.00)x=0.0001 phia=atan(y/x) call orb_2py(Za,ra,thetaa,phia,psia,psi2a,psimaxa,psi2maxa) psia=Ca*psia c write(6,*)'ra,psia=',ra,psia x=x-Re rb=sqrt(x*x+y*y+z*z) if(rb.eq.0.0)rb=0.0001 thetab=acos(z/rb) if(x.eq.0.00)x=0.0001 phib=atan(y/x) call orb_2py(Zb,rb,thetab,phib,psib,psi2b,psimaxb,psi2maxb) psib=Cb*psib c write(6,*)'rb,psib=',rb,psib psi=Ca*psia+isn*Cb*psib psi2=psi*psi x=x+Re/2. if(ifrst.eq.0)then ifrst=1 psimax=Ca*psimaxa+Cb*psimaxa psi2max=psimax*psimax c write(6,*)'psimaxa,psimaxb,psimax,psi2max' c write(6,*)psimaxa,psimaxb,psimax,psi2max end if return end subroutine mo_sig1s2p &(isn,Za,Zb,Ca,Cb,Re,x,y,z,psi,psi2,psimax,psi2max) static ifrst data ifrst/0/ c write(6,*)'Za,Zb,Ca,Cb,Re,x,y,z' c write(6,*)Za,Zb,Ca,Cb,Re,x,y,z c Convert this point to spherical polar coordinates in each atomic ref frame x=x+Re/2. ra=sqrt(x*x+y*y+z*z) if(ra.eq.0.0)ra=0.0001 c thetaa=acos(z/ra) c if(x.eq.0.00)x=0.0001 c phia=atan(y/x) call orb_1s(Za,ra,psia,psi2a,psimaxa,psi2maxa) psia=Ca*psia c write(6,*)'ra,psia=',ra,psia x=x-Re rb=sqrt(x*x+y*y+z*z) if(rb.eq.0.0)rb=0.0001 thetab=acos(z/rb) if(x.eq.0.00)x=0.0001 phib=atan(y/x) call orb_2px(Zb,rb,thetab,phib,psib,psi2b,psimaxb,psi2maxb) psib=Cb*psib c write(6,*)'rb,psib=',rb,psib psi=Ca*psia+isn*Cb*psib psi2=psi*psi x=x+Re/2. if(ifrst.eq.0)then ifrst=1 psimax=Ca*psimaxa+Cb*psimaxb psi2max=psimax*psimax c write(6,*)'psimaxa,psimaxb,psimax,psi2max' c write(6,*)psimaxa,psimaxb,psimax,psi2max end if return 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