program gnnj
      implicit none
      integer*4 nat,i,ia,ig,is,ix,iy,iz,nao,nmo,n2
      real*8 step,x,g0(9),gt(9),si,stepau,d0(3),dt(3),mt(3),m0(3),p,
     1q0(6),qt(6),v0(3),vt(3)
      character*180 s80
      character*1 xyz(3)
      data xyz/'x','y','z'/
      real*8,allocatable::g(:),d(:),m(:),v(:),q(:)
      real*8,allocatable::c0(:),c1(:),c2(:),s(:),bc0(:),bc1(:),bc2(:)
      logical lphase
c     debug option:
      lphase=.false.
c
      write(6,600)
600   format(' Derivatives of MCD Gnj tensor')
c
c     zero geometry:
      nat=0
      step=0.001d0
      nao=0
      nmo=0
      open(9,file='0/GV.OUT',status='old')
11    read(9,940,end=88,err=88)s80
940   format(a180)
      if(s80(2:14).eq.'Nuclear step=')read(s80(15:23),*)step
      if(s80(27:44).eq.'Input orientation:')then
       write(6,900)s80
       read(9,*)
       read(9,*)
       read(9,*)
       read(9,*)
       nat=0
200    read(9,*,err=87)x,x,x,x,x,x
       nat=nat+1
       goto 200
87     continue
      endif
      if(s80(1:8).eq.' NBasis=')read(s80(9:14),*)nao
      if(s80(1:8).eq.' NBsUse=')then
       read(s80(9:14),*)nmo
       goto 88
      endif
      goto 11
88    close(9)
      write(6,*)nat,' atoms'
      write(6,*)nao,' atomic orbitals'
      write(6,*)nmo,' molecular orbitals'

      if(lphase)then
       n2=nao*nao
       allocate(c0(n2),c1(n2),c2(n2),s(n2),bc0(n2),bc1(n2),bc2(n2))
c      read S-matrix of AOs, suppose all the same for other geometries
       call rds(s,nao)
c      read C and c
       call rdc(0,c0,bc0,nao,nmo)
       call gp(p,c0,bc0,c0,bc0,s,nao,nmo)
       write(6,609)0,0,p
609    format(' phase ',2i4,f10.3)
      endif

      call rg(0,g0,d0,m0,q0,v0)
      allocate(g(9*3*nat),d(9*nat),m(9*nat),v(9*nat),q(18*nat))
      call vz(g,9*3*nat)
      call vz(d,9*nat)
      call vz(m,9*nat)
      call vz(v,9*nat)
      call vz(q,18*nat)
      ig=0
      do 1 ia=1,nat
      do 1 ix=1,3
      si=-1.0d0
      do 1 is=1,2
      write(6,*)ia,ix,is
      si=-si
      ig=ig+1
      call rg(ig,gt,dt,mt,qt,vt)

c     read C c coef for phase determination
      if(lphase)then
       if(is.eq.1)then
        call rdc(ig,c1,bc1,nao,nmo)
        call gp(p,c1,bc1,c0,bc0,s,nao,nmo)
        write(6,609)0,ig,p
       else
        call rdc(ig,c2,bc2,nao,nmo)
        call gp(p,c2,bc2,c0,bc0,s,nao,nmo)
        write(6,609)0,ig,p
        call gp(p,c1,bc1,c2,bc2,s,nao,nmo)
        write(6,609)ig-1,ig,p
       endif
      endif

      do 1 iy=1,9
      if(iy.lt.4)then
       d(iy+3*(ix-1)+ 9*(ia-1))=d(iy+3*(ix-1)+ 9*(ia-1))+si*dt(iy)
       m(iy+3*(ix-1)+ 9*(ia-1))=m(iy+3*(ix-1)+ 9*(ia-1))+si*mt(iy)
       v(iy+3*(ix-1)+ 9*(ia-1))=v(iy+3*(ix-1)+ 9*(ia-1))+si*vt(iy)
      endif
      if(iy.lt.7)then
       q(iy+6*(ix-1)+18*(ia-1))=q(iy+6*(ix-1)+18*(ia-1))+si*qt(iy)
      endif
1     g(iy+9*(ix-1)+27*(ia-1))=g(iy+9*(ix-1)+27*(ia-1))+si*gt(iy)
      stepau=step/0.529177d0
      do 3 i=1,27*nat
3     g(i)=g(i)/(2.0d0*stepau)
      do 13 i=1,18*nat
13    q(i)=q(i)/(2.0d0*stepau)
      do 5 i=1,9*nat
      v(i)=v(i)/(2.0d0*stepau)
      m(i)=m(i)/(2.0d0*stepau)
5     d(i)=d(i)/(2.0d0*stepau)

      open(90,file='GNJD.TEN')
      do 2 ia=1,nat
      do 2 iz=1,3
      write(90,890)ia,xyz(iz)
890   format(' Atom ',i6,1x,a1)
      do 2 ix=1,3
  2   write(90,900)(g(ix+3*(iy-1)+9*(iz-1)+27*(ia-1)),iy=1,3)
c      ix=electric
c      iy=magnetic
c      iz=atomic
 900  format(3g14.6)
      write(90,*)'g0:'
      do 6 ix=1,3
6     write(90,401)(g0(ix+3*(iy-1)),iy=1,3)
401   format(3g12.4)
      close(90)

      open(40,file='DE.TEN')
      call wd('  dipx         dipy         dipz:',nat,d,3,'u0:',d0)
      call wd(' dipvx        dipvy        dipvz:',nat,v,3,'v0:',v0)
      call wd(' mdipx        mdipy        mdipz:',nat,m,3,'m0:',m0)
      call wd('  Q       xx  xy  xz  yy  yz  zz:',nat,q,6,'q0:',q0)
      close(40)

      end

      subroutine wd(s1,n,v,m,s2,u)
      implicit none
      character*(*)s1,s2
      integer*4 n,ia,ix,iy,m
      real*8 v(*),u(*)
      write(40,*)'atom coord'//s1
      do 1 ia=1,n
      do 1 ix=1,3
c     change sign, i.e. record <r>, not <er>, etc:
1     write(40,400)ia,ix,(-v(iy+m*(ix-1)+3*m*(ia-1)),iy=1,m)
400   format(i5,i2,6g12.4)
      write(40,*)s2
      write(40,400)0,0,(-u(iy),iy=1,m)
      return
      end

      subroutine rg(io,g,d,m,q,v)
      implicit none
      integer*4 io,is,iy,ix
      real*8 g(*),d(*),m(*),q(*),v(*)
      character*10 s10
      write(s10,10)io
10    format(i10)
      do 3 is=1,len(s10)
3     if(s10(is:is).ne.' ')goto 4
4     open(90,file=s10(is:len(s10))//'/GNJ.TEN')
      read(90,*)
      do 1 ix=1,3
1     read(90,*)(g(ix+3*(iy-1)),iy=1,3)
c     ix index-electric
      read(90,*)
      read(90,*)
      read(90,*)(d(ix),ix=1,3)
      read(90,*)
      read(90,*)(v(ix),ix=1,3)
      read(90,*)
      read(90,*)(m(ix),ix=1,3)
      read(90,*)
      read(90,*)(q(ix),ix=1,6)
      close(90)
      return
      end

      subroutine vz(g,n)
      integer*4 n,i
      real*8 g(*)
      do 1 i=1,n
1     g(i)=0.0d0
      return
      end

      subroutine rds(s,nao)
      implicit none
      real*8 s(*)
      integer*4 nao,i,j
c     reads s-matrix (AO|AO)
      open(33,file='0/SAO.SCR',form='unformatted')
      do 1 i=1,nao
      read(33)(s(j+nao*(i-1)),j=1,i)
      do 2 j=1,i-1
2     s(i+nao*(j-1))=s(j+nao*(i-1))
1     continue
      close(33)
      return
      end
            
      subroutine rdc(io,c,bc,nao,nmo)
      implicit none
      integer*4 io,is,ieigen,i1,i2,iex,ir,nao,nmo,ia,ib,i,ie,nd,ni
      real*8 c(*),bc(*)
      real*8,allocatable::e(:),be(:),b(:),d(:)
      character*10 s10
      character*35 key
      character*80 s80

      open(9,file='0/GNJ.TEN')
      read(9,*)iex
      close(9)

      do 200 i1=1,nao*nao
      c(i1)=0.0d0
200   bc(i1)=0.0d0

      allocate(e(nmo),be(nmo),b(nao*nao),d(nao*nao))
      write(s10,10)io
10    format(i10)
      ieigen=0
      ie=0
      ir=0
      do 3 is=1,len(s10)
3     if(s10(is:is).ne.' ')goto 4
4     open(2,file=s10(is:len(s10))//'/GV.OUT')
1     read(2,900,end=99,err=99)key
900   format(a35)
      if(key(6:35).eq.'Molecular Orbital Coefficients'.or.
     1   key(2:11).eq.'Alpha MOs:') then
       ieigen=ieigen+1
c      eigenvectors listed either as a result of IOP(6/7=1) or IOP(5/33=1)
       call readorb(key,e,c,nao,nmo)
       do 377 i1=1,nmo
       be(i1)=e(i1)
       do 377 i2=1,nao
377    b(i2+nao*(i1-1))=c(i2+nao*(i1-1))
       read(2,900)key
       if(key(2:10).eq.'Beta MOs:')then
        write(6,900)key
        call readorb(key,be,b,nao,nmo)
       endif
      endif

      if(key(2:15).eq.'Excited State ')then
       do 203 i=1,len(key)
203    if(key(i:i).eq.':')read(key(15:i-1),*)nd
       ie=ie+1
       if(iex.eq.ie.and.ie.eq.nd)then
        ir=ir+1
        ni=0
111     read(2,2000,end=299,err=299)s80
2000    format(a80)
        if(s80(10:11).eq.'->')then
         if(s80(8:8).eq.'A')then
c         open shell, alpha:
          ni=ni+1
          read(s80( 1: 7),*)ia
          if(s80(15:15).eq.'A')then
           read(s80(12:14),*)ib
          else
           if(s80(16:16).eq.'A')then
            read(s80(12:15),*)ib
           else
            write(6,2000)s80
            call report('Unknown format of excitations')
           endif
          endif
          read(s80(18:30),*)bc(ia+nao*(ib-1))
         else
          if(s80(8:8).eq.'B')then
c          open shell, beta:
           ni=ni+1
           read(s80( 1: 7),*)ia
           if(s80(15:15).eq.'B')then
            read(s80(12:14),*)ib
           else
            if(s80(16:16).eq.'B')then
             read(s80(12:15),*)ib
            else
             write(6,2000)s80
             call report('Unknown format of excitations')
            endif
           endif
           read(s80(18:30),*) d(ia+ib*(nao-1))
          else
c          closed shell:
           ni=ni+1
           read(s80( 1: 8),*)ia
           read(s80(12:15),*)ib
           read(s80(18:30),*)bc(ia+nao*(ib-1))
          endif
         endif
         goto 111
        endif
299     continue
        write(6,*)ni,' CI coefficients'
       endif
      endif

      if(ir.eq.0.or.ieigen.eq.0) goto 1

99    close(2)
      if(ieigen.eq.0)call report('eigenvectors not found')
      if(ir.eq.0)call report('CI vectors not found')
      return
      end

      subroutine gp(phase,c1,bc1,c2,bc2,s,nao,nmo)
      implicit none
      real*8 phase,c1(*),bc1(*),c2(*),bc2(*),s(*),f
      integer*4 i,ii,jj,j,nao,nmo
      real*8,allocatable::t(:),k(:)
      allocate(t(nao*nao),k(nao*nao))
      do 1 i=1,nmo
      do 1 ii=1,nao
      t(i+nmo*(ii-1))=0.0d0
      do 1 jj=1,nao
1     t(i+nmo*(ii-1))=t(i+nmo*(ii-1))+c1(jj+nao*(i-1))*s(jj+nao*(ii-1))
      do 2 i=1,nmo
      do 2 j=1,nmo
      k(i+nmo*(j-1))=0.0d0
      do 2 ii=1,nao
2     k(i+nmo*(j-1))=k(i+nmo*(j-1))+c2(ii+nao*(j-1))*t(i+nao*(ii-1))
      phase=0.0d0
      do 3 i=1,nmo
      do 3 j=i+1,nmo
      f=bc1(i+nao*(j-1))
      if(f.ne.0.0d0)then
       do 31 ii=1,nmo
       do 31 jj=ii+1,nmo
31     phase=phase+bc1(i+nao*(j-1))*bc2(ii+nao*(jj-1))*
     1(k(i+nmo*(jj-1))*k(j+nmo*(ii-1))+k(i+nmo*(ii-1))*k(j+nmo*(jj-1)))
      endif
3     continue
      phase=phase*2.0d0
      return
      end

      subroutine readorb(key,e,c,nao,nmo)
      implicit none
      integer*4 i,j,i1,j1,i2,nao,nmo,ic
      real*8 e(*),c(*)
      logical number
      character*(*) key
      character*160 s160
      i=1
      read(2,*)
377   j=i+4
      if(j.gt.nmo)j=nmo
      if(key(8:11).ne.'MOs:'.and.key(7:10).ne.'MOs:')read(2,*)
c
c     to enable reading of various formats, use s160:
      read(2,160)s160
160   format(a160)
c
c     add space before negative numbers:
      i1=22
3775  if(s160(i1:i1).eq.'-'.and.s160(i1-1:i1-1).ne.'E')then
       do 1 j1=160,i1+1,-1
1      s160(j1:j1)=s160(j1-1:j1-1)
       s160(i1:i1)=' '
       i1=i1+1
      endif
      i1=i1+1
      if(i1.lt.160)goto 3775

      read(s160(22:160),*)(e(i1),i1=i,j)
      do 39 i2=1,nao
      read(2,160)s160

c     add space before negative numbers:
      ic=0
      i1=22
3776  if(s160(i1:i1).eq.'-'.
     1and.number(s160(i1-1:i1-1)).and.number(s160(i1+1:i1+1)))then
       ic=ic+1
       write(6,*)  s160(j1-1:j1-1)
       do 2 j1=160,i1+1,-1
2      s160(j1:j1)=s160(j1-1:j1-1)
       s160(i1:i1)=' '
       i1=i1+1
      endif
      i1=i1+1
      if(i1.lt.160)goto 3776
      if(ic.gt.0)then
       write(6,*)'Corrected:',ic
       write(6,160)s160
      endif

39    read(s160(22:160),*,err=99)(c(i2+nao*(i1-1)),i1=i,j)

      if(j.lt.nmo)then
       i=i+5
       read(2,*)
       goto 377
      endif
      write(6,*)'MOs read'
      return
99    write(6,160)s160
      write(6,*)i2,i,j
      call report('input error')
      end

      subroutine report(s)
      character*(*) s
      write(6,*)s
      stop
      end

      function number(s)
      character*(*) s
      logical number
      number=s.eq.'1'.or.s.eq.'2'.or.s.eq.'3'.or.s.eq.'4'.or.s.eq.'5'
     1   .or.s.eq.'6'.or.s.eq.'7'.or.s.eq.'8'.or.s.eq.'9'.or.s.eq.'0'
      return
      end