program glueq
c     makes G98.OUT and G98.INP
c     from smaller files
      parameter (nat0=600,nf0=10,nh0=40)
      character*80 s80,hsi(nh0),hnsi(nh0),cmsi,hosi(nh0)
      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),fft(3*nat0,3*nat0,nf0),
     3ffi(3*nat0,3*nat0),nfs(nf0),nfe(nf0),fxa(3*nat0*nf0),
     4gxi(3*nat0),gxt(3*nat0,nf0)
      logical lex
c
      open(3,file='GLUEQ.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
      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
      nfs(i)=nat+1
      nat=nat+nai
      nfe(i)=nat
      es=es+ei
      do 2 ii=1,3*nai
      gxt(ii,i)=gxi(ii)
2     xt(ii,i)=xi(ii)
      do 4 ii=1,nai
4     ityt(ii,i)=ityi(ii)
      nq(i)=nqi
      do 5 ii=1,nqi
5     cqt(ii,i)=cqi(ii)
      do 3 ii=1,3*nqi
3     xqt(ii,i)=xqi(ii)
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)
c
c     force field
      call ru7 
      read(7,7000)s80
      write(3,7000)s80
      do 91 istart=1,len(s80)
91    if(s80(istart:istart).ne.' ')goto 10
10    do 11 iend=len(s80),1,-1
11    if(s80(iend:iend).ne.' ')goto 12
12    call loadff1(NAT0,ffi,nai,s80(istart:iend))
      do 13 ii=1,3*nai
      do 13 jj=1,3*nai
13    fft(ii,jj,i)=ffi(ii,jj)

1     continue
      close(7)
c
      call wrg98(nat0,nf0,xt,ityt,nq,xqt,cqt,es,gxt,nf,na)

      inquire(file='FILE.X',exist=lex)
      if(lex)then
       call msg('FILE.X exists and is not changed')
      else
       open(7,file='FILE.X')
       write(7,*)nat
       write(7,*)nat
       do 14 i=1,nf
       do 14 ia=1,na(i)
14     write(7,7001)ityt(ia,i),(xt(3*(ia-1)+ix,i),ix=1,3)
7001   format(i4,3f15.8,' 0 0 0 0 0 0 0 0.0')
       close(7)
       call msg('Phony FILE.X written') 
      endif

      inquire(file='FILE.FC',exist=lex)
      if(lex)then
       call msg('FILE.FC exists and is not changed')
      else
       open(10,file='FILE.FC')
       N=3*nat
       N1=1
111    N3=N1+4
       IF(N3.GT.N)N3=N

       DO 130 LN=N1,N

       IA=(LN-1)/3+1
       do 132 i=1,nf
132    if(IA.ge.nfs(i).and.IA.le.nfe(i))igi=i
       noff=3*(nfs(igi)-1)
c
       do 131 J=N1,MIN(LN,N3)
       JA=(J-1)/3+1
       do 133 i=1,nf
133    if(JA.ge.nfs(i).and.JA.le.nfe(i))igj=i
       if(igi.ne.igj)then
        fxa(J)=0.0
       else
        fxa(J)=fft(LN-noff,J-noff,igi)
       endif
131    continue
c      
130    WRITE(10,17)LN,(fxa(J),J=N1,MIN(LN,N3))

       N1=N1+5
       IF(N3.LT.N)GOTO 111
17     FORMAT(I4,5D14.6)
       close(10)
       call msg('Force field written.')
      endif

c
c     G98.INP from the last file read
      open(10,file='G98.INP')
      do 18 i=1,nhi
18    write(10,7000)hsi(i)
      write(10,*)
      do 19 i=1,nsi
19    write(10,7000)hnsi(i)
      write(10,*)
      write(10,7000)cmsi
c
      do 20 i=1,nf
      do 20 ia=1,na(i)
20    write(10,1100)ityt(ia,i),(xt(3*(ia-1)+ix,i),ix=1,3)
1100  format(i4,2x,3f15.8)
      write(10,*)
c
      if(nqi.gt.0)then 
       do 22 i=1,nf
       do 22 ii=1,nq(i)
22     write(10,119)(xqt(3*(ii-1)+ix,i),ix=1,3),cqt(ii,i)
119    format(4f12.6)
      endif
      write(10,*)
      do 21 i=1,noi
21    write(10,7000)hosi(i) 
      write(10,*)
      close(10)
c
      call msg('Phony gaussian input written in G98.INP')
   
      close(3)
      stop
      end
c ======================================================================
      subroutine chkchek(ts,lchk)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      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)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      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 wrg98(nat0,nf0,xt,ityt,nq,xqt,cqt,es,gxt,nf,na)
      DIMENSION xt(3*NAT0,nf0),ityt(NAT0,nf0),xqt(3*nat0,nf0),
     1gxt(3*NAT0,nf0),nq(nf0),na(nf0),cqt(nat0,nf0)
      character*80 filename,ts
      logical*4 lcharges
c
      lcharges=.false.
      filename='G98.OUT'
      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
c     gaussian output
      open(10,file=filename)
      do 2 i=1,80
2     ts(i:i)=' '
      ts(27:44)='Input orientation:'
      write(10,1000)ts
      do 4 i=1,4
4     write(10,*)
      ii=0
      do 5 i=1,nf
      do 5 ia=1,na(i)
      ii=ii+1
5     write(10,101)ii,ityt(ia,i),i,(xt(3*(ia-1)+ix,i),ix=1,3)
101   format(3i5,3f15.8)
      write(10,103)
103   format(1x,71(1h-))
c
      if(lcharges)then
       do 6 i=1,80
6      ts(i:i)=' '
       ts(2:15)='Point Charges:'
       write(10,1000)ts
       do 7 i=1,nf
       do 7 iq=1,nq(i)
7      write(10,10012)(xqt(3*(iq-1)+ix,i),ix=1,3),cqt(iq,i)
10012  format(' XYZ=',3f10.4,3x,f10.4)
      endif

      write(10,10001)es
10001 format(' SCF Done:   SUM  = ',f20.8,' A')

      do 61 i=1,80
61    ts(i:i)=' '
      ts(38:59)='Forces (Hartrees/Bohr)'
      write(10,1000)ts
      write(10,*)
      write(10,*)
      ii=0
      do 62 i=1,nf
      do 62 ia=1,na(i)
      ii=ii+1
62    write(10,1004)ii,ityt(ia,i),(gxt(3*(ia-1)+ix,i),ix=1,3)
1004  format(2i5,3f20.8)
      write(10,1005)
1005  format(' -------------')

      close(10)
      call msg('Phony G98.OUT made')
      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 ======================================================================
      subroutine loadff1(NAT0,fx,nat,filename)
      DIMENSION fx(3*NAT0,3*NAT0)
      character*(*) filename
c
      open(20,file=filename)
      N=3*nat
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(20,17)(fx(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)

c     to comply with new4 transfer ff in H/bohr^2 to mdyn/A^2 or something
      DO 3 I=1,N
      DO 3 J=I,N
3     fx(J,I)=fx(J,I)
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 fx(I,J)=fx(J,I)
      close(20)
      call msg('Cartesian FF read')
      RETURN
      END
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)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      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)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      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 ======================================================================