program numforce
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      parameter (nat0=1000)
      character*80 filename,st
      real*8 r0(3,nat0),g0(3,nat0),r(3),g(3,nat0),f(3*nat0,3*nat0),
     1pp(nat0,3,3),JJ(nat0,3,3),I0(nat0,3,3),pv(nat0,3,3),p0(3),
     2p(3),QQ(nat0*3,6),q0(6),q(6),qqcheck(3,6)
      dimension iz(nat0)
      logical lex
c
      bohr=0.529177d0
      write(6,*)' FF from gradients'
      write(6,*)' Z-Matrix orientation required !!!'
      write(6,*)
      write(6,*)' Zero-point output filename: '
      open(11,file='TEM.OUT')
      read(5,5000)filename
5000  format(a80)
      open(15,file=filename)
1     read(15,5000,err=900,end=900)st
      if(st(26:46).eq.'Z-Matrix orientation:'.or.
     1   st(27:44).eq.'Input orientation:')then
       write(11,5000)st
       do 14 i=1,4
       write(11,*)
14     read(15,*)
       i=0
2      read(15,5000)st
       write(11,5000)st
       if(st(2:3).ne.'--')then
        backspace 15
        i=i+1
        if(i.gt.nat0)then
         write(6,*)'too many atoms'
         close(15)
         stop
        endif
        read(15,*)idum,iz(i),idum,(r0(j,i),j=1,3)
        goto 2
       endif
       nat=i
       write(6,*)nat,' atoms read'
      endif
      if(st(38:59).eq.'Forces (Hartrees/Bohr)')then
       write(11,5000)st
       do 24 i=1,2
       write(11,*)
24     read(15,*)
       do 25 i=1,nat
25     read(15,*)idum,idum,(g0(j,i),j=1,3)
       do 251 i=1,nat
251    write(11,10000)i,iz(i),(g0(j,i),j=1,3)
10000  format(2i8,3f17.9)
       write(6,*)'Gradient read'
       goto 900
      endif
c     
c     PM3 dipole moment:
      if(st(1:15).eq.' Dipole moment=')then
       write(11,*)'point 0'
       write(11,5000)st
       backspace 15
       read(15,156)(p0(i),i=1,3)
156    format(15x,3f10.6)
      endif
c     
c     
c     SCF dipole moment:
      if(st(1:23).eq.' Dipole moment (Debye):'.or.
     1 st(1:23).eq.' Dipole moment (field-i')then
       write(11,*)'point 0'
       write(11,5000)st
       read(15,158)(p0(i),i=1,3)
158    format(3(6x,f11.4))
      endif
c     
c     SCF quadrupole moment:
      if(st(1:18).eq.' Quadrupole moment')then
       write(11,*)'point 0'
       write(11,5000)st
       read(15,158)(q0(i),i=1,6)
       write(11,158)(q0(i),i=1,6)
c      XX YY ZZ XY XZ YZ
      endif
      goto 1
900   close(15)
c
      inquire(file='FILE.FC',exist=lex)
      if(lex)then
       OPEN(20,file='FILE.FC')
       CALL READFF(3*NAT0,3*NAT,F)
       CLOSE(20)
       write(6,*)'Initial FC read from FILE.FC'
      else
       do 3 i=1,3*nat
       do 3 j=1,3*nat
3      f(i,j)=0.0d0
       write(6,*)'Initial FC zeroed-out'
      endif
C
      inquire(file='FILE.TEN',exist=lex)
      if(lex)then
       OPEN(15,FILE='FILE.TEN',STATUS='OLD')
       READ (15,*)NOAT
       if(NOAT.ne.nat)then
        write(6,*)'nat <> noat'
	close(15)
	stop
       endif
       DO 10 L=1,NAT
       DO 10 J=1,3
10     READ (15,*) (PP(L,J,I),I=1,3)
       DO 220 L=1,NAT
       DO 220 J=1,3
 220   READ (15,*) (I0(L,J,I),I=1,3)
       DO 230 L=1,NAT
       DO 230 J=1,3
230    READ (15,*) (JJ(L,J,I),I=1,3)
       DO 100 L=1,NAT
       DO 100 J=1,3
100    READ (15,*) (PV(L,J,I),I=1,3)
       CLOSE(15)
       WRITE(*,*)' Dipole derivatives read in'
      else 
       DO 101 L=1,NAT
       DO 101 J=1,3
       DO 101 I=1,3
       JJ(L,J,I)=0.0d0
       PV(L,J,I)=0.0d0
       I0(L,J,I)=0.0d0
101    PP(L,J,I)=0.0d0
       WRITE(*,*)' Dipole derivatives zeroed-out'
      endif

      DO 102 L=1,NAT
      DO 102 J=1,3
      DO 102 I=1,6
102   QQ(3*(L-1)+J,I)=0.0d0
      WRITE(*,*)' Quadrupole derivatives zeroed-out'
c
      write(6,*)' Num. diff. output filename: '
      read(5,5000)filename
      ixdif=0
      iadif=0
      diff=0.0d0
      open(15,file=filename)
4     read(15,5000,err=901,end=901)st
      if(st(26:46).eq.'Z-Matrix orientation:'.or.
     1   st(27:44).eq.'Input orientation:')then
       write(11,5000)st
       do 34 i=1,4
       write(11,*)
34     read(15,*)
       do 35 i=1,nat
       read(15,*)idum,izi,idum,(r(j),j=1,3)
       write(11,2000)i,iz(i),0,(r(j),j=1,3)
2000   format(3i6,3f15.6)
       if(izi.ne.iz(i))then
        write(6,*)'Zs different'
        close(15)
        stop
       endif
       do 36 j=1,3
       del=r(j)-r0(j,i)
       if(abs(del).gt.1.0d-5)then
        iadif=i
        ixdif=j
        diff=del
       endif
36     continue
35     continue
       write(11,1111)
1111   format(' -----')
       dd=diff/bohr
       write(6,600)iadif,ixdif,diff,dd
       write(11,600)iadif,ixdif,diff,dd
600    format('atom ',i6,' coord ',i2,f15.5,' A,',f15.5,' Bohr')
      endif
      if(st(38:59).eq.'Forces (Hartrees/Bohr)')then
       write(11,5000)st
       if(diff.lt.1.0d-5)then
        write(6,*)'skipped zero point'
        goto 4
       endif
       do 44 i=1,2
       write(11,*)
44     read(15,*)
       do 45 i=1,nat
       read(15,*)idum,idum,(g(j,i),j=1,3)
45     write(11,10000)i,iz(i),(g(j,i),j=1,3)
       write(6,*)'Gradient read'
       jr=3*(iadif-1)+ixdif
       do 38 i=1,nat
       do 38 ix=1,3
       ii=3*(i-1)+ix
       fnew=-(g(ix,i)-g0(ix,i))/dd
       f1old=f(jr,ii)
       f2old=f(ii,jr)
       oo=0.0d0
       if(abs(f1old).gt.0.0d0)oo=oo+1.0d0
       if(abs(f2old).gt.0.0d0)oo=oo+1.0d0
       if(oo.gt.0.0d0)then
        fold=(f1old+f2old)/oo
c       if(abs((fold-fnew)/fold).gt.1.2d0)then
c        write(6,6666)iadif,ixdif,i,ix,fold,fnew
c6666     format(i4,i2,i4,i2,2f15.10,' big difference - set to 0')
c        fold=0.0d0
c        fnew=0.0d0
c       endif
        fnew=(fold+fnew)*0.05d0
       endif
       f(jr,ii)=fnew
38     f(ii,jr)=fnew
       write(6,*)'Forces updated'
      endif
c     
c     SCF quadrupole moment:
      if(st(1:18).eq.' Quadrupole moment')then
       write(11,5000)st
       if(abs(diff).lt.1.0d-5)then
        write(6,*)'skipped zero point q'
        goto 4
       endif
       read(15,158)(Q(i),i=1,6)
       write(11,158)(Q(i),i=1,6)
c      debye-angstrom to au:
c      quad(ours)=(1/2)quad(from gaussian)
       const=1.0d0/bohr/2.541767d0*0.5d0
       if(dd.gt.0)then
        do 392 i=1,6
392     QQ((iadif-1)*3+ixdif,i)=(Q(i)-q0(i))/dd*const
       else
        write(6,6392)(QQ((iadif-1)*3+ixdif,i),i=1,6)
6392    format(6f11.4)
        do 393 i=1,6
393     QQ((iadif-1)*3+ixdif,i)=
     1  0.5d0*(QQ((iadif-1)*3+ixdif,i)+(Q(i)-q0(i))/dd*const)
       endif
       write(6,*)'Quadrupole derivatives updated'
      endif

      if(st(1:15).eq.' Dipole moment=')then
       write(11,5000)st
       if(diff.lt.1.0d-5)then
        write(6,*)'skipped zero point dipole'
        goto 4
       endif
       backspace 15
       read(15,156)(p(i),i=1,3)
       do 39 i=1,3
39     pp(iadif,ixdif,i)=(p(i)-p0(i))/dd
       write(6,*)'Dipole derivatives updated'
      endif

      if(st(1:23).eq.' Dipole moment (Debye):'.or.
     1 st(1:23).eq.' Dipole moment (field-i')then
       write(11,5000)st
       if(diff.lt.1.0d-5)then
        write(6,*)'skipped zero point dipole'
        goto 4
       endif
       read(15,158)(p(i),i=1,3)
       if(dd.gt.0)then
        do 391 i=1,3
391     pp(iadif,ixdif,i)=(p(i)-p0(i))/dd/2.541767d0
       else
        do 3911 i=1,3
3911    pp(iadif,ixdif,i)=0.5d0*(pp(iadif,ixdif,i)
     1  +(p(i)-p0(i))/dd/2.541767d0 )
       endif
       write(6,*)'Dipole derivatives updated'
      endif
      goto 4
901   close(15)
      close(11)
      OPEN(20,file='FILENEW.FC')
      CALL WRITEFF(3*NAT0,3*NAT,F)
      CLOSE(20)
      write(6,*)'Forces written into FILENEW.FC'
c      
      OPEN(15,FILE='FILEQQ.TEN')
      write(15,1515)NAT
1515  format(i5,' atoms, quadrupole derivatives (au)',/,
     1'at. ix        XX        YY        ZZ        XY        XZ ',
     2'     YZ')

      write(6,*)'Checksums:'
      do 53 ic=1,2
      if(ic.eq.2) CALL PROJECTQ(QQ,r0,nat,6)
      DO 52 J=1,3
      do 52 I=1,6
      qqcheck(J,I)=0.0d0
      DO 52 L=2,NAT
52    qqcheck(J,I)=qqcheck(J,I)-QQ(3*(L-1)+J,I)       
      do 53 k=1,3
      write(6,1571)k,k,(QQ(k,I),I=1,6)
      write(6,1571)k,k,(qqcheck(k,I),I=1,6)
53    write(6,*)

      DO 50 L=1,NAT
      DO 50 J=1,3
50    write(15,1571)L,J, (QQ(3*(L-1)+J,I),I=1,6)
1571  format(i4,i2,6f10.4)
      CLOSE(15)
      WRITE(*,*)' Quadrupole derivatives written into FILEQQ.TEN'
c
      OPEN(15,FILE='FILENEW.TEN')
      write(15,*)NAT
      DO 61 L=1,NAT
      DO 61 J=1,3
61    write(15,157) (PP(L,J,I),I=1,3),L
157   format(3f14.8,i5)
      DO 420 L=1,NAT
      DO 420 J=1,3
 420  write(15,157) (I0(L,J,I),I=1,3),l
      DO 430 L=1,NAT
      DO 430 J=1,3
430   write(15,157) (JJ(L,J,I),I=1,3),l
      DO 400 L=1,NAT
      DO 400 J=1,3
400   write(15,157) (PV(L,J,I),I=1,3),l
      CLOSE(15)
      WRITE(*,*)' Dipole derivatives written into FILENEW.TEN'
      stop
      end
C
      SUBROUTINE READFF(MX3,N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(MX3,N)
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   READ(20,17)(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(4X,5D14.6)
C
c     CONST=4.359828/0.5291772/0.5291772
      DO 3 I=1,N
      DO 3 J=I,N
3     FCAR(J,I)=FCAR(J,I)
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      WRITE(*,*)' Cartesian FF read in ... '
      RETURN
      END

      SUBROUTINE WRITEFF(MX3,N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(MX3,N)
c     CONST=4.359828/0.5291772**2
      DO 6 I=1,N
      DO 6 J=1,N
6     FCAR(I,J)=FCAR(I,J)
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   WRITE(20,17)LN,(FCAR(LN,J),J=N1,MIN(LN,N3))
      N1=N1+5
      IF(N3.LT.N)GOTO 1
17    FORMAT(I4,5D14.6)
      RETURN
      END


C
      SUBROUTINE PROJECTQ(Q,r0,nat,n1)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      parameter (nat0=1000,MX3=3*nat0)
      DIMENSION Q(nat0*3,n1),T(nat0*3,6),r0(3,nat0),A(MX3,MX3),
     1TEM(MX3,MX3)
      N=3*nat
      WRITE(*,*)'Projecting transl/rotations'
      if(n1.gt.6)call report('string too large')
      CALL DOMA(r0,nat,TEM,A)
      DO 2 I=1,n1
      DO 1 J=1,N
      T(J,I)=0.0d0
      DO 1 II=1,N
1     T(J,I)=T(J,I)+A(J,II)*Q(II,I)
      do 2 J=1,N
2     Q(J,I)=T(J,I)
      RETURN
      END

      SUBROUTINE DOMA(C,nat,TEM,A)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      PARAMETER (nat0=1000,MX3=3*nat0,MX=MX3/3)
      DIMENSION C(3,MX),A(MX3,MX3),TEM(MX3,MX3)
      N=3*nat
      AMACH=0.00000000000001d0
      DO 12 I=1,N
      DO 12 J=1,6
12    A(I,J)=0.0d0
      DO 1 I=1,NAT
      A(3*I-2,1)= 1.00d0
      A(3*I-1,2)= 1.00d0
      A(3*I  ,3)= 1.00d0
      A(3*I-2,4)=-C(2,I)
      A(3*I-1,4)= C(1,I)
      A(3*I-1,5)=-C(3,I)
      A(3*I  ,5)= C(2,I)
      A(3*I-2,6)= C(3,I)
1     A(3*I  ,6)=-C(1,I)
      N6=6
      CALL ORT(A,MX3,N6,NAT)
      ICOL=1
5     IF(ICOL.GT.N6)GOTO 8
      S=0.0d0
      DO 6 J=1,N
6     S=S+A(J,ICOL)**2
      IF(S.LT.1000.0d0*AMACH)THEN
       DO 7 J=1,N
7      A(J,ICOL)=A(J,N6)
       N6=N6-1
       GOTO 5
      ENDIF
      ICOL=ICOL+1
      GOTO 5
8     CONTINUE
      WRITE(*,*)N6,' coordinates projected'
      DO 2 I=1,N
      DO 2 J=1,N
      TEM(I,J)=0.0d0
      DO 2 II=1,N6
2     TEM(I,J)=TEM(I,J)+A(I,II)*A(J,II)
      DO 3 I=1,N
      DO 4 J=1,N
4     A(I,J)=-TEM(I,J)
3     A(I,I)=A(I,I)+1.0d0
      RETURN
      END
C     ============================================================ 
      SUBROUTINE ORT(A,NDIM,N6,NAT)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(NDIM,N6)
      AMACH=0.00000000000001d0
      NAT3=3*NAT
      C=0.0d0
      DO 1 I=1,NAT3
1     C=C+A(I,1)**2
      IF(C.GT.AMACH)C=1.0d0/SQRT(C)
      DO 2 I=1,NAT3
2     A(I,1)=A(I,1)*C
      DO 3 IC=2,N6
C     Column IC, orthonormalize to previous columns:
      DO 3 ICP=1,IC-1
      SP=0.0d0
      DO 5 J=1,NAT3
5     SP=SP+A(J,IC)*A(J,ICP)
C     Subtract projection to ICP:
      DO 6 J=1,NAT3
6     A(J,IC)=A(J,IC)-SP*A(J,ICP)
C     Normalize :
      SP=0.0d0
      DO 7 J=1,NAT3
7     SP=SP+A(J,IC)**2
      IF(SP.GT.AMACH)SP=1.0d0/SQRT(SP)
      DO 3 J=1,NAT3
3     A(J,IC)=A(J,IC)*SP
      RETURN
      END

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