program splitq
c     divides INP.NEW into more inputs
      parameter (nat0=600,nf0=10,nh0=40)
      character*80 s80,hsi(nh0),hnsi(nh0),cmsi,hosi(nh0),
     1hst(nh0,nf0),hnst(nh0,nf0),cmst(nf0),host(nh0,nf0)
      dimension na(nf0),xi(3*nat0),xt(3*nat0,nf0),nq(nf0),
     1xqi(3*nat0),xqt(3*nat0,nf0),cqi(nat0),cqt(nat0,nf0),
     2ityi(nat0),ityt(nat0,nf0),
     4gxi(3*nat0),nht(nf0),nst(nf0),not(nf0)

c
      open(3,file='SPLITQ.OUT')
      open(7,file='GLUEQ.OPT')
      call ru7 
      read(7,*)nf
      if(nf.gt.nf)call bye('Too many fragments.')
      write(3,*)nf,' fragments'
      nat=0
      nqt=0
      es=0.0
      do 1 i=1,nf
c
      call ru7 
      read(7,*)nai
      if(nai.gt.nat0)call bye('Too many atoms')
      write(3,*)i,nai,' atoms'
c
c     gradient
      call ru7 
      read(7,7000)s80
      write(3,7000)s80
7000  format(a80)
      do 6 istart=1,len(s80)
6     if(s80(istart:istart).ne.' ')goto 7
7     do 8 iend=len(s80),1,-1
8     if(s80(iend:iend).ne.' ')goto 9
9     call loadx1(nat0,xi,nai,ityi,nqi,xqi,cqi,s80(istart:iend))
      call loadg1(nat0,gxi,nai,ei,s80(istart:iend))
      na(i)=nai
      nat=nat+nai
      es=es+ei
      do 4 ii=1,nai
4     ityt(ii,i)=ityi(ii)
      nq(i)=nqi
      nqt=nqt+nqi
c
c     gaussian input
      call ru7 
      read(7,7000)s80
      write(3,7000)s80
      do 92 istart=1,len(s80)
92    if(s80(istart:istart).ne.' ')goto 93
93    do 94 iend=len(s80),1,-1
94    if(s80(iend:iend).ne.' ')goto 95
95    call ginp1(nqi,s80(istart:iend),nhi,hsi,nsi,hnsi,cmsi,noi,hosi)
      nht(i)=nhi
      do 96 ii=1,nhi
96    hst(ii,i)=hsi(ii)
      nst(i)=nsi
      do 97 ii=1,nsi
97    hnst(ii,i)=hnsi(ii)
      cmst(i)=cmsi
      not(i)=noi
      do 98 ii=1,noi
98    host(ii,i)=hosi(ii)
c
c     force field
      call ru7 
      read(7,7000)s80
1     continue
      close(7)
c
c     coordinates from INP.NEW
      call ginp2(nqt,xt,xqt,cqt,nf,na,nq)
c
c     write them into sepparate inputs
      open(7,file='GLUEQ.OPT')
      call ru7 
      read(7,*)nf
      do 100 i=1,nf
      call ru7 
      read(7,*)nai
      call ru7 
      read(7,7000)s80
c     gaussian input
      call ru7 
      read(7,7000)s80
      write(3,7000)s80
      do 192 istart=1,len(s80)
192   if(s80(istart:istart).ne.' ')goto 193
193   do 194 iend=len(s80),1,-1
194   if(s80(iend:iend).ne.' ')goto 195
195   open(71,file=s80(istart:iend))
      lchk=0
      loniom=0
      do 50 ii=1,nht(i)
      call chkchek(hst(ii,i),lchk)
      call chkoniom(hst(ii,i),loniom)
50    write(71,7000)hst(ii,i)
      if(lchk.eq.0)write(71,*)'guess=checkpoint'
      if(loniom.ne.0)call bye('oniom not allowed yet')
      write(71,*)
      do 51 ii=1,nst(i)
51    write(71,7000)hnst(ii,i)
      write(71,*)
      write(71,7000)cmst(i)
      do 52 ia=1,na(i)
52    write(71,7171)ityt(ia,i),(xt(3*(ia-1)+ix,i),ix=1,3)
7171  format(i4,3f15.8)
      write(71,*)
      if(nq(i).ne.0)then
       do 53 iq=1,nq(i)
53     write(71,7172)(xqt(3*(iq-1)+ix,i),ix=1,3),cqt(iq,i)
7172   format(4f15.7)
       write(71,*)
      endif
      do 54 ii=1,not(i)
54    write(71,7000)host(ii,i)
      write(71,*)
      close(71)
c
c     force field
      call ru7 
100   read(7,7000)s80
      close(7)

      call msg('Sepparate inputs written created')
   
      close(3)
      stop
      end
c ======================================================================
      subroutine chkchek(ts,lchk)
      CHARACTER*(*) ts
      CHARACTER*16 s
c
      l=len(ts)
      do 1 i=1,l-15
      do 2 ii=1,16
      ia=ichar(ts(i+ii-1:i+ii-1))
      if(ia.ge.65.and.ia.le.90)ia=ia+32
2     s(ii:ii)=char(ia)
1     if(s.eq.'guess=checkpoint')lchk=1
      return
      end
c ======================================================================
      subroutine chkoniom(ts,loniom)
      CHARACTER*(*) ts
      CHARACTER*16 s
c
      l=len(ts)
      do 1 i=1,l-4
      do 2 ii=1,5
      ia=ichar(ts(i+ii-1:i+ii-1))
      if(ia.ge.65.and.ia.le.90)ia=ia+32
2     s(ii:ii)=char(ia)
1     if(s.eq.'oniom')loniom=1
      return
      end
c ======================================================================
      subroutine loadg1(NAT0,gx,nat,e,filename)
      DIMENSION gx(*)
      character*80 filename,ts
c
      load=0
      open(10,file=filename,status='old')
2     read(10,1000,end=899,err=899)ts
1000  format(a80)
      if(ts(2:29).eq.'ONIOM: extrapolated energy =')then
       write(3,1000)ts
       open(12,file='scr')
       write(12,'(a)')ts(30:80)
       rewind 12
       read(12,*)e
       close(12)
      endif
      if(ts(1:8).eq.' Energy='.and.ts(28:32).eq.'NIter')then
       write(3,1000)ts
       open(12,file='scr')
       write(12,'(a)')ts(9:27)
       rewind 12
       read(12,*)e
       close(12)
      endif
      if(ts(1:10).eq.' SCF Done:')then
       write(3,1000)ts
       e=exe(ts)
      endif
      if(ts(38:59).eq.'Forces (Hartrees/Bohr)')then
       read(10,*)
       read(10,*)
       load=1
       ia=0
3      read(10,1000)ts
       if(ts(2:3).ne.'--')then
        ia=ia+1
        if(ia.gt.nat0)then
         call bye('Too many atoms')
         stop
        endif
        backspace 10
        read(10,*)idum,idum,(gx(3*(ia-1)+ix),ix=1,3)
        goto 3
       endif
      endif
      goto 2
899   close(10)
      nat=ia
      if(load.eq.0)then
       call bye('Gradient not found')
      else
       call msg('Cartesian gradient loaded.')
       cnorm=0.0
       cmax=gx(1)
       do 4 i=1,3*nat
       cnorm=cnorm+gx(i)**2
4      if(abs(gx(i)).gt.abs(cmax))cmax=gx(i)
       cnorm=sqrt(cnorm)
       rms=cnorm/sqrt(real(3*nat))
       write(3,3003)cnorm,cmax,rms
3003   format(' Norm ',f15.8,/,' Max  ',f15.8,/,' RMS  ',f15.8,/)
      endif
      return
      end
c ======================================================================
C     ============================================================ 
      subroutine loadx1(nat0,x,nat,ity,nq,xq,cq,filename)
      DIMENSION x(*),ity(*),xq(*),cq(*)
      character*80 ts
      character*(*) filename
      character*4 s4
      logical*4 lex,lcharges
c
      nq=0
      lcharges=.false.
      open(10,file='Q.OPT')
1     read(10,1000,end=999,err=999)ts
1000  format(a80)
      if(ts(1:6).eq.'CHARGE')read(10,*)lcharges
      goto 1
999   close(10)
c
      inquire(file=filename,exist=lex)
      if(.not.lex)call bye('Output does not exist')
      open(10,file=filename,status='old')
c
c     gaussian output
2     read(10,1000,end=899,err=899)ts
      if(ts(27:44).eq.'Input orientation:'.or.
     1   ts(26:46).eq.'Z-Matrix orientation:')then
       do 4 i=1,4
4      read(10,*)
       do 5 i=1,nat
       read(10,*)idum,ity(i),idum,(x(3*(i-1)+ix),ix=1,3)
5      if(ity(i).lt.0)read(10,*)idum,ity(i),idum,(x(3*(i-1)+ix),ix=1,3)
      endif
      if(lcharges.and.ts(2:15).eq.'Point Charges:')then
       i=0
88     read(10,112)s4
112    format(1x,a4)
       if(s4.eq.'XYZ=')then
        backspace 10
        i=i+1
        if(i.gt.NAT0)call bye('Too many charges.')
        read(10,10012)(xq(3*(i-1)+ix),ix=1,3),cq(i)
10012   format(5x,3f10.4,3x,f10.4)
        goto 88
       endif 
       nq=i
       goto 899
      endif
      goto 2
899   close(10)
      call msg('Coordinates loaded')
      if(lcharges.and.nq.eq.0)then
       call bye('No charges found')
      else
       call msg('Charges:')
       do 7 i=1,nq
7      write(3,301)i,(xq(3*(i-1)+ix),ix=1,3),cq(i)
301    format(i5,3f10.4,' Q = ',f10.4)
      endif
      return
      end
C     ============================================================ 
      subroutine msg(ts)
      character*(*) ts
c
      do 1 i=min(80,len(ts)),1,-1
1     if(ts(i:i).ne.' ')goto 2
2     write(3,*)ts(1:i)
      return
      end
C     ============================================================ 
      subroutine bye(ts)
      character*(*) ts
      write(3,6001)
      write(3,'(a)')ts
      write(3,6001)
6001  format(80(1h*))
      write(3,6000)
6000  format('program stopped')
      close(3)
      stop
      end
c ======================================================================
      function exe(ts)
      character*(*) ts
      ns=len(ts)
      do 1 i1=1,ns
1     if(ts(i1:i1).eq.'=')goto 2
2     do 3 i2=i1,ns
3     if(ts(i2:i2).eq.'A')goto 4
4     open(12,file='scr')
      write(12,'(a)')ts(i1+1:i2-1)
      rewind 12
      read(12,*)e
      exe=e
      close(12)
      return
      end
c ======================================================================
      subroutine ru7 
      character*1 s1
c
c     read file 7 until line does not begin with #
1     read(7,7000,err=99,end=99)s1
7000  format(a1)
      if(s1.ne.'#')then
       backspace 7
       return
      endif
      goto 1
99    call bye('End of option file')
      end
c ======================================================================
      function notemptyline(ts)
      LOGICAL*4 notemptyline
      CHARACTER*(*) ts
c
      l=len(ts)
      notemptyline=.false.
      do 1 i=1,l
1     if(ts(i:i).ne.' ')notemptyline=.true.
      return
      end
c ======================================================================
      subroutine ginp1(nq,filename,nh,hs,ns,hns,cms,no,hos)
      PARAMETER (nh0=40)
      LOGICAL*4 notemptyline
      CHARACTER*80 filename,ts,hs(nh0),hns(nh0),cms,hos(nh0)
c
      open(10,file=filename)
c
c     header
      nh=0
      lchk=0
2     read(10,1000,end=899,err=899)ts
1000  format(a80)
      call chkchek(ts,lchk)
      if(notemptyline(ts))then
       nh=nh+1
       if(nh.gt.nh0)call bye('too many header lines')
       hs(nh)=ts
       goto 2
      else
       if(lchk.eq.0)then
        nh=nh+1
        if(nh.gt.nh0)call bye('too many header lines')
        hs(nh)='guess=checkpoint'
       endif
      endif 
c
c     name
      ns=0
21    read(10,1000,end=899,err=899)ts
      if(notemptyline(ts))then
       ns=ns+1
       hns(ns)=ts 
       goto 21
      endif
c
c     charge and multiplicity
      read(10,1000,end=899,err=899)cms
c
c     Cartesian coordinates skip
23    read(10,1000,end=899,err=899)ts
      if(notemptyline(ts)) goto 23
c
c     fix point charges-skip
      if(nq.gt.0)then
24     read(10,1000,end=899,err=899)ts
       if(notemptyline(ts)) goto 24
      endif
c
c     optional parameters until end of file:
      no=0
31    read(10,1000,end=899,err=899)ts
      if(notemptyline(ts))then
       no=no+1
       hos(no)=ts
       goto 31
      endif
c      
899   close(10)
      call msg('Input read')
      return
      end
c ======================================================================
      subroutine ginp2(nqt,xt,xqt,cqt,nf,na,nq)
      PARAMETER (nat0=600,nf0=10)
      LOGICAL*4 notemptyline
      CHARACTER*80 ts
      dimension xt(3*nat0,nf0),xqt(3*nat0,nf0),cqt(nat0,nf0),
     1na(nf0),nq(nf0)
c
      open(10,file='INP.NEW')
c     header
2     read(10,1000,end=899,err=899)ts
1000  format(a80)
      if(notemptyline(ts))goto 2
c     name
21    read(10,1000,end=899,err=899)ts
      if(notemptyline(ts))goto 21
c     charge and multiplicity
      read(10,1000,end=899,err=899)ts
c     Cartesian coordinates 
      do 23 i=1,nf
      do 23 ia=1,na(i)
      read(10,1000,end=899,err=899)ts
c
c     get rid of the character:
      do 3 is=1,len(ts)
3     if(ts(is:is).ne.' ')goto 4
4     do 5 ie=is,len(ts)
5     if(ts(ie:ie).eq.' ')goto 6
6     open(45,file='scr')
      write(45,*)ts(ie+1:len(ts))
      rewind 45
      read(45,*,end=899,err=899)(xt(3*(ia-1)+ix,i),ix=1,3)
23    close(45)
      read(10,*)
c     fix point charges
      if(nqt.gt.0)then
       do 24 i=1,nf
       do 24 iq=1,nq(i)
24     read(10,*,end=899,err=899)(xqt(3*(iq-1)+ix,i),ix=1,3),cqt(iq,i)
       read(10,*)
      endif
c     optional parameters skipped
      close(10)
      call msg('Input read')
      return
899   close(10)
      call bye('INP.NEW could not be read')
      end
c ======================================================================