******************************************************************************** * * * 2DHF version 1-2003 * * Copyright (C) 1996 Jacek Kobus, Leif Laaksonen, Dage Sundholm * * * * This software may be used and distributed according to the terms * * of the GNU General Public License, see README and COPYING. * * * ******************************************************************************** c ### prepg94l ### c c Reads the output from the GAUSSSIAN94 program (gauss94l.out file) c to determine parameters of the basis functions (n, l, m, exponents) c and coefficients of molecular orbitals. c The basis can contain functions defined on centres A and B and c also at the bond centre. c c Last modified 7/09/02 c This routine extracts data from the customized (long) output c (see ISKIP below) c subroutine prepg94l implicit integer*4 (i-n) implicit real*8 (a-h,o-z) character*2 st character*3 conbt character*9 matchstr include 'commons8.inc' parameter (symthresh=0.001d0) parameter(maxbasis=200) dimension sexpon(3,maxbasis),excoeff(0:60,maxbasis), & conbt(maxbasis),ibfres(maxbasis) dimension ifbord(60),ifdord(60) c nfborb - number of finite basis (fb) set orbitals c ifbord - ordering of fb orbitals (0=sigma, 1=pi, 2=delta, etc) c the order is determined from the gaussian94.output file c ifdord - ordering of finite difference (fd) orbitals c the order of fd orbitals is the same as in the input data c (phi orbitals first, then delta, etc) c c BF cc data nfborb /7/ cc data ifbord /0,0,0,0,1,1,0,53*100/ cc data ifdord /1,0,0,0,0,0,54*100/ c c InF cc data nfborb /29/ cc data ifbord /0,0,1,1,0, cc & 0,0,1,1,0, cc & 2,2,1,1,0, cc & 0,1,1,0,0, cc & 2,2,1,1,0, cc & 0,1,1,0,100, cc & 30*100/ cc & cc data ifdord /2*2,6*1,13*0,39*100/ c c opening a file with the corresponding (customized) GAUSSIAN94 output open(7,file='gauss94l.out', status='old',form='formatted') c ibc counts contracted basis functions c ibp counts primitive basis functions c n1prim - number of primitive gaussian function at the centre A c n3prim - number of primitive gaussian function at the centre B c n2prim - number of primitive gaussian function at the bond centre iprint0=0 iprint1=0 iprint2=0 if (idbg(553).ne.0) iprint0=1 if (idbg(554).ne.0) iprint1=1 if (idbg(555).ne.0) iprint2=1 c format of data should be detected automatically c if Gaussian94 output contains eigenvectors in 2-column format iskip=2 c if Gaussian94 output contains eigenvectors in 5-column format iskip=5 iskip=2 c if (idbg(560).ne.0) iskip=5 c determine number of molecular orbitals as defined in the GAUSSIAN94 c program (one nonsigma finite difference orbital corresponds c to two GAUSSIAN94 ones do i=0,60 do j=1,maxbasis excoeff(i,j)=0.d0 enddo enddo nfborb=0 do iorb=1,norb ifdord(iorb)=mm(iorb) if (mm(iorb).eq.0) then nfborb=nfborb+1 else nfborb=nfborb+2 endif enddo write(*,*) icentre=0 ilabel=1 ib=0 ibpt=0 icent=0 do ilines=1,100000 c read(7,1001,end=991,err=991) matchstr,icentt read(7,1001,end=991,err=991) matchstr c write(*,1001) matchstr 01000 format(a9) 01001 format(a9,i5) if (matchstr.eq.' Centers:') then icent=icent+1 icentre=icentre+1 ibc=0 ipb=0 c start extracting data from the output 00100 ib=ib+1 ibc=ibc+1 c read(7,*,end=910,err=910) st,expon read(7,1010,end=910,err=910) st if (st.eq.'**') goto 910 01010 format(1x,a2,6x,f18.9) 01011 format(14x,f18.9) read(7,1011,end=910,err=910) expon sexpon(icent,ibc) =expon if (st.eq.'s '.or.st.eq.'S ') then ipb=ipb+1 ibpt=ibpt+1 primexp(ibpt)=expon lprim(ibpt)=0 mprim(ibpt)=0 icgau(ibpt)=icent elseif (st.eq.'p '.or.st.eq.'P ') then ipb=ipb+3 do k=1,3 ibpt=ibpt+1 primexp(ibpt)=expon lprim(ibpt)=1 icgau(ibpt)=icent if (k.eq.1) mprim(ibpt)=+1 if (k.eq.2) mprim(ibpt)=-1 if (k.eq.3) mprim(ibpt)= 0 enddo elseif (st.eq.'d '.or.st.eq.'D ') then ipb=ipb+5 do k=1,5 ibpt=ibpt+1 primexp(ibpt)=expon lprim(ibpt)=2 icgau(ibpt)=icent if (k.eq.1) mprim(ibpt)= 0 if (k.eq.2) mprim(ibpt)=+1 if (k.eq.3) mprim(ibpt)=-1 if (k.eq.4) mprim(ibpt)=+2 if (k.eq.5) mprim(ibpt)=-2 enddo elseif (st.eq.'f '.or.st.eq.'F ') then ipb=ipb+7 do k=1,7 ibpt=ibpt+1 primexp(ibpt)=expon lprim(ibpt)=3 icgau(ibpt)=icent if (k.eq.1) mprim(ibpt)= 0 if (k.eq.2) mprim(ibpt)=+1 if (k.eq.3) mprim(ibpt)=-1 if (k.eq.4) mprim(ibpt)=+2 if (k.eq.5) mprim(ibpt)=-2 if (k.eq.6) mprim(ibpt)=+3 if (k.eq.7) mprim(ibpt)=-3 enddo elseif (st.eq.'g '.or.st.eq.'G ') then ipb=ipb+9 do k=1,9 ibpt=ibpt+1 primexp(ibpt)=expon lprim(ibpt)=4 icgau(ibpt)=icent if (k.eq.1) mprim(ibpt)= 0 if (k.eq.2) mprim(ibpt)=+1 if (k.eq.3) mprim(ibpt)=-1 if (k.eq.4) mprim(ibpt)=+2 if (k.eq.5) mprim(ibpt)=-2 if (k.eq.6) mprim(ibpt)=+3 if (k.eq.7) mprim(ibpt)=-3 if (k.eq.8) mprim(ibpt)=+4 if (k.eq.9) mprim(ibpt)=-4 enddo elseif (st.eq.'h '.or.st.eq.'H ') then ipb=ipb+11 do k=1,11 ibpt=ibpt+1 primexp(ibpt)=expon lprim(ibpt)=5 icgau(ibpt)=icent if (k.eq. 1) mprim(ibpt)= 0 if (k.eq. 2) mprim(ibpt)=+1 if (k.eq. 3) mprim(ibpt)=-1 if (k.eq. 4) mprim(ibpt)=+2 if (k.eq. 5) mprim(ibpt)=-2 if (k.eq. 6) mprim(ibpt)=+3 if (k.eq. 7) mprim(ibpt)=-3 if (k.eq. 8) mprim(ibpt)=+4 if (k.eq. 9) mprim(ibpt)=-4 if (k.eq.10) mprim(ibpt)=+5 if (k.eq.11) mprim(ibpt)=-5 enddo elseif (st.eq.'sp'.or.st.eq.'SP') then ipb=ipb+4 ibpt=ibpt+1 primexp(ibpt)=expon lprim(ibpt)=0 mprim(ibpt)=0 icgau(ibpt)=icent do k=1,3 ibpt=ibpt+1 primexp(ibpt)=expon lprim(ibpt)=1 icgau(ibpt)=icent if (k.eq.1) mprim(ibpt)=+1 if (k.eq.2) mprim(ibpt)=-1 if (k.eq.3) mprim(ibpt)= 0 enddo endif if (iprint0.ne.0) then print *,ibpt,icgau(ibpt), & lprim(ibpt),mprim(ibpt),primexp(ibpt) endif goto 100 910 continue ib=ib-1 if (icentre.eq.1) then ibret=ib ibrett=ibret n1prim=ipb elseif (icentre.eq.2) then ibret=ib-ibrett ibrett=ibrett+ibret n2prim=ipb elseif (icentre.eq.3) then ibret=ib-ibrett ibrett=ibrett+ibret n3prim=ipb endif write(*,1141) icent,ibret endif 991 continue if (icentre.eq.3) goto 992 enddo 992 continue rewind (7) nexpon=ib npbasis=n1prim+n2prim+n3prim write(*,1140) nexpon,npbasis,n1prim,n2prim,n3prim idavid=1 c start extracting basis set expansion coefficientsdata from the output c since David Moncrieff supplied data in various formats they are all c tried one by one 900 continue rewind(7) do ilines=1,1000000 read(7,1000,end=990,err=990) matchstr if (iprint1.ne.0) write(*,*) '1',matchstr,'2' if (matchstr.eq.' EIGE') then do iorb=1,nfborb,iskip do ib=1,npbasis if (idavid.eq.1) then c iskip.eq.2 read(7,1015,end=940,err=800) ibt,ibr,conbt(ib), & excoeff(iorb,ib),excoeff(iorb+1,ib) 1015 format(i4,6x,i3,a3,11x,f19.16,5x,f19.16) elseif (idavid.eq.2) then c iskip.eq.2 read(7,1115,end=940,err=801) ibt,ibr,conbt(ib), & excoeff(iorb,ib),excoeff(iorb+1,ib) c 1115 format(a4,6x,i3,a3,5x,2f25.16) 1115 format(i5,6x,i3,a3,5x,2f20.16) elseif (idavid.eq.10) then c iskip.eq.5 read(7,1016,end=940,err=810) ibt,ibr,conbt(ib), & excoeff(iorb ,ib),excoeff(iorb+1,ib), & excoeff(iorb+2,ib),excoeff(iorb+3,ib), & excoeff(iorb+4,ib) 1016 format(a4,6x,i3,a3,5x,5f10.5) elseif (idavid.eq.11) then c iskip.eq.5 read(7,1016,end=940,err=811) ibt,ibr,conbt(ib), & excoeff(iorb ,ib),excoeff(iorb+1,ib), & excoeff(iorb+2,ib),excoeff(iorb+3,ib), & excoeff(iorb+4,ib) 1116 format(a4,6x,i3,a3,5x,5f10.5) endif if (abs(excoeff(iorb,ib)).gt.1.d0.and. & iorb.le.nfborb) then write(*,1130) iorb,ib,excoeff(iorb,ib) endif if (abs(excoeff(iorb+1,ib)).gt.1.d0.and. & iorb+1.le.nfborb) then write(*,1130) iorb+1,ib,excoeff(iorb+1,ib) endif ibfres(ib)=ibr if (ib.le.n1prim) icent=1 if (ib.gt.n1prim.and.ib.le.(n1prim+n2prim)) icent=2 if (ib.gt.(n1prim+n2prim) & .and.ib.le.(n1prim+n2prim+n3prim)) icent=3 if (iprint1.ne.0) then if (iskip.eq.2) then write(*,1013,err=990) iorb,ib,ibr,icent, & conbt(ib), & excoeff(iorb,ib),excoeff(iorb+1,ib), & sexpon(icent,ibr) elseif(iskip.eq.5) then write(*,1014,err=990) iorb,ib,ibr,icent, & conbt(ib), & excoeff(iorb, ib),excoeff(iorb+1,ib), & excoeff(iorb+2,ib),excoeff(iorb+3,ib), & excoeff(iorb+4,ib), & sexpon(icent,ibr) endif endif enddo read(7,1000,end=950,err=950) matchstr read(7,1000,end=950,err=950) matchstr read(7,1000,end=950,err=950) matchstr enddo endif enddo 800 continue idavid=2 iskip=2 goto 900 801 continue idavid=10 iskip=5 goto 900 810 continue idavid=11 iskip=5 goto 900 811 continue WRITE(*,*) 'prepg94l: check format of gauss94l.out file' stop 'gauss94l' 940 continue write(*,*) 'prepg94l: uncomplete gauss94l file', ib-1,ibt stop 'prepg94l' 950 continue if (iprint1.ne.0) write(*,*) 'warning: 950' 990 continue if (iprint1.ne.0) write(*,*) 'warning: 990' write(*,1142) norb,nfborb c determine symmetry of GAUSSIAN94 molecular orbitals c icount=0 do ifbo=nfborb,1,-1 do ib=npbasis,1,-1 if (abs(excoeff(ifbo,ib)).gt.symthresh) then if(icount.ne.ifbo) then icount=ifbo if (iprint2.ne.0) print *,'ifbo,ib,mprim(ib) ', & ifbo,ib,mprim(ib) c ifbord(ifbo)=lprim(ib) ifbord(ifbo)=abs(mprim(ib)) endif endif enddo enddo c fb orbitals are being associated with the corresponding fd ones do iorb=1,norb do ifbo=nfborb,1,-1 icount=0 if (ifdord(iorb).eq.ifbord(ifbo)) then if (iprint2.ne.0) then print *,'iorb,ifbo', & iorb,ifdord(iorb),ifbo,ifbord(ifbo) endif ib=0 500 continue ib=ib+1 if (ib.le.npbasis) then if (mprim(ib).eq.0) then primcoef(iorb,ib)=excoeff(ifbo,ib) else excoeffp=excoeff(ifbo,ib) excoeffm=excoeff(ifbo-1,ib) primcoef(iorb,ib)= & sqrt(excoeffp*excoeffp+excoeffm*excoeffm) if (excoeffp.lt.0.d0.or.excoeffm.lt.0.d0) & primcoef(iorb,ib)=-primcoef(iorb,ib) endif else goto 510 endif goto 500 00510 continue if (ifbord(ifbo).gt.0) ifbord(ifbo-1)=-1 ifbord(ifbo)=-1 ifdord(iorb)=-2 endif enddo enddo c normalization factor for a spherical harmonic Gaussian-type functions do ib=1,npbasis d1=primexp(ib) l1=lprim (ib) m1=mprim (ib) m1abs=iabs(m1) fngau2(ib)=( d1**(2*l1+3) * & 2**(4*l1+7)/pii/(factor2(2*l1+1))**2)**0.25d0 c normalization factor for spherical harmonics shngau(ib)=(-1.d0)**((m1+m1abs)/2) /sqrt(4.d0*pii)* & sqrt((2*l1+1)*factor(l1-m1abs)/factor(l1+m1abs)) enddo c 01011 format(14x,f18.9) 01012 format(13x,a3,5x,5f10.5) 01013 format(4i5,3x,a3,2f19.9,d15.6) 01014 format(4i5,3x,a3,5f12.5,d15.6) 01042 format(i4,i3,1x,5x,a3,5x,5f10.5) 01052 format(2i4,5x,a3,5x,5f10.5) 01020 format(i3,2i4,5x,f20.9,f15.5) 01032 format(2i4,i3,a2,2x,f18.9) 01035 format(4i4,5x,2d25.10) 01130 format(1x,'... expansion coefficient greater than 1.0', & ' detected ...',/,1x,' orbital:',i3, & ', basis function:',i4,', coefficient:',f6.3,/) 01140 format(/15x,i4,' exponents ', & /15x,i4,' primitive basis functions ', & /15x,i4,' primitives on centre 1' & /15x,i4,' primitives on centre 2' & /15x,i4,' primitives on centre 3'/) 01141 format(6x,'centre',i2,':',i4,' exponents ') 01142 format(15x,i4,' finite difference orbitals', & /15x,i4,' finite basis set orbitals'/) return end