******************************************************************************** * * * 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 ### prepg94 ### c c This routine reads the output from the GAUSSSIAN94 program c (gauss94.out and gauss94.pun files) to determine parameters c of the basis functions (n, l, m, exponents) and coefficients of c molecular orbitals. c The basis can contain functions defined on centres A and B and c also at the bond centre. c subroutine prepg94 implicit integer*4 (i-n) implicit real*8 (a-h,o-z) character*46 matchstr include 'commons8.inc' parameter(maxbasis=_MAXBASIS_) dimension excoeff(0:60,0:maxbasis),ccoeff(maxbasis) dimension ifbord(60),ifdord(60) c nforb - 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) 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 opening a file with the corresponding GAUSSSIAN94 output open(7,file='gauss94.out', status='old',form='formatted') do ilines=1,1000 read(7,1001,end=990,err=910) matchstr if (matchstr.eq. & ' Basis set in the form of general basis input:') goto 900 enddo write(*,*) 'prepg94: ' write(*,*) ' "Basis set in the form of general basis input:"', & ' not found in gaussian94.output file' stop 'prepg94' 00900 continue ipbasis=1 npbasis=0 c start extracting data from the output 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 ibc=0 ibp=0 n1prim=0 n2prim=0 n3prim=0 call r_exponents(ibc,ibp,istop) n1prim=ibp if (istop.eq.0) then call r_exponents(ibc,ibp,istop) n3prim=ibp-n1prim if (istop.eq.0) then icent=2 call r_exponents(ibc,ibp,istop) n2prim=ibp-n3prim-n1prim endif endif 00910 continue nexpon=ibc npbasis=ibp write(*,1140) nexpon,npbasis,n1prim,n2prim,n3prim close(7) open(7,file='gauss94.pun', status='old',form='formatted') read(7,1001,end=991,err=990) c start extracting basis set expansion coefficients from c gaussian94.pun file c determine number of molecular orbitals as defined in the GAUSSIAN94 c program (one nonsigma finite difference orbital corresponds c to two GAUSSIAN94 ones 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(*,1142) norb,nfborb c retrieve expansion coefficients of the contracted gaussians if (iprint1.ne.0) then print *,'no. of basis function no. of exp. coeff.' do ibp=1,npbasis ibc=ixref(ibp) enddo endif do ifbo=1,nfborb read(7,1011,end=991,err=992) oe do ibc=1,nexpon,5 read(7,1012,end=991,err=992) (ccoeff(ibc+i),i=0,4) enddo c transform expansion coefficients of the contracted gaussians c into coefficients of primitive gaussians do ibp=1,npbasis ibc=ixref(ibp) excoeff(ifbo,ibp)=ccoeff(ibc)*coeff(ibp) c write(*,'(3i5,3e16.6,i5)') ifbo,ibp,ibc,ccoeff(ibc), c & coeff(ibp),excoeff(ifbo,ibp),mprim(ibp) enddo enddo c prepare expansion coefficients of a molecular orbital in the basis c set of primitive gaussian functions do ib=1,npbasis 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 icgau(ib)=icent c write(*,1050) ib,ixref(ib),lprim(ib),mprim(ib),icgau(ib), c & primexp(ib) enddo c determine symmetry of GAUSSIAN94 molecular orbitals symthresh=0.00001d0 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),ibc ', & ifbo,ib,mprim(ib),ixref(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 ib=1,npbasis primcoef(iorb,ib)=0.d0 enddo enddo do iorb=1,norb do ifbo=nfborb,1,-1 icount=0 if (ifdord(iorb).eq.ifbord(ifbo)) then do ib=1,npbasis primcoef(iorb,ib)=excoeff(ifbo,ib) enddo 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 return 00990 print *,'PREPG94: end of file encountered when reading ', & 'gaussian94.output' stop 00991 print *,'PREPG94: end of file encountered when reading ', & 'gaussian94.pun' stop 00992 stop 'PREPG94: error encountered when reading gaussian94.pun' c 01001 format(a46) 01011 format(15x,d15.8) 01012 format(5d15.8) 01050 format(5i5,d15.8) 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'/) 01142 format(15x,i4,' finite difference orbitals', & /15x,i4,' finite basis set orbitals'/) end