program pfadd
      implicit none
      integer*4 nat,n,iargc,i,j,ig,icav,iz,ico
      real*8 hartreetoJ,bohr,pau,p,d,cc
      real*8,allocatable::r(:),r0(:),g(:),f(:,:),f0(:,:)
      integer*4,allocatable::z(:),ib(:),io(:)
      character*80 s80
      hartreetoJ=4.3597482d-18
      bohr=0.529177d0

      write(6,600)
600   format(' Makes gradients in FILE.OUT from pressure',/,
     1' and PCM cavity, makes FF and adds it to FILE.FC')
      if(iargc().lt.1)call report('Usage: pfad <p/Gp>')
      call getarg(1,s80)
      read(s80,*)p
      pau=p/hartreetoJ*1.0d-30*bohr**3*1.0d9
      cc=0.0d0
      if(iargc().gt.1)then
       call getarg(2,s80)
       read(s80,*)cc
      endif

      allocate(z(1),r0(1))
      call rg('FILE.OUT',0,nat,z,r0,.true.)
      deallocate(z,r0)
      n=3*nat
      allocate(z(nat),r0(n),g(n),r(n),f(n,n),ib(n),io(n),
     1f0(n,n))
      call rg('FILE.OUT',1,nat,z,r0,.true.)
      f=0.0d0

      ig=0
      icav=0
      iz=0
      d=0.0d0
      ib=0
      io=0
      open(2,file='FILE.OUT')
1     call rgf(nat,z,r,.true.,ig)
      write(6,*)
      write(6,*)
      if(nat.ne.0)write(6,601)'   Geometry',ig
601   format(a11,i6)
      i=ico(n,r,r0,d)
      call mkg(n,g,pau,r,r0,icav,cc)
      if(nat.ne.0)write(6,601)'     Cavity',icav
      if(nat.eq.0)goto 2
      write(6,607)i,d
607   format(' Coordinate',i6,', shift',f12.4)
      if(i.ne.0)then
       write(6,602)i,d
       d=d/bohr
602    format(i6,' coordinate, step ',f10.4,' A')
       if(d.lt.0.0d0)then
        ib(i)=ib(i)+1
        do 3 j=1,n
3       f(i,j)=f(i,j)-g(j)/2.0d0/dabs(d)
       else
        io(i)=io(i)+1
        do 31 j=1,n
31      f(i,j)=f(i,j)+g(j)/2.0d0/dabs(d)
       endif
      else
       write(6,601)' Zero geom.'
       iz=iz+1
      endif
      goto 1
2     close(2)
      write(6,608)ig,icav
608   format(i6,' gradients,',i6,' cavities')
      if(iz.ne.1)then
       write(6,604)iz
604    format(i6,' zero geometries ... check FILE.OUT')
       stop
      endif
      do 4 i=1,n
      if(ib(i).ne.1)then
       write(6,605)i,ib(i)
605    format(' Coord.',i6,':',i6,' back geoms, check FILE.OUT')
       stop
      endif
      if(io(i).ne.1)then
       write(6,606)i,io(i)
606    format(' Coord.',i6,':',i6,' forward geoms, check FILE.OUT')
       stop
      endif
4     continue
      do 5 i=1,n
      do 5 j=1,i
c     write(6,609)i,j,f(i,j),f(j,i)
c09   format(2i5,2f12.6)
      f(i,j)=(f(i,j)+f(j,i))/2.0d0
5     f(j,i)=f(i,j)
      call project(F,r,N)
      call READFF(N,f0)
      do 6 i=1,n
      do 6 j=1,n
c     minus here because the gradient is actually the force:
6     f0(i,j)=f0(i,j)-f(i,j)
      call WRITEFF(N,f0)
      end

      SUBROUTINE PROJECT(F,r,N)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION F(N,N),r(N)
      real*8,allocatable::TEM(:,:),A(:,:)
C
      allocate(TEM(N,N),A(N,N))
      WRITE(*,*)'Projecting transl/rotations from the force field...'
      call DOMA(A,n,r)
      DO 11 I=1,N
      DO 11 J=1,N
11    TEM(I,J)=0.0d0

      DO 1 I=1,N
      DO 1 II=1,N
      if(dabs(A(I,II)).gt.1.0d-10)then
       DO 12 J=1,N
12     TEM(I,J)=TEM(I,J)+A(I,II)*F(II,J)
      endif
1     continue

      DO 22 I=1,N
      DO 22 J=1,N
22    F(I,J)=0.0d0

      DO 2 J=1,N
      DO 2 II=1,N
      if(dabs(A(J,II)).gt.1.0d-10)then
       DO 23 I=1,N
23     F(I,J)=F(I,J)+TEM(I,II)*A(J,II)
      endif
2     continue
      WRITE(*,*)'  ... done.'
      RETURN
      END


      function ico(n,r,r0,d)
      implicit none
      integer*4 ico,n,i
      real*8 tol,r(n),r0(n),d
      tol=1.0d-5
      do 1 i=1,n
      d=r(i)-r0(i)
      if(dabs(d).gt.tol)then
       ico=i
       return
      endif
1     continue
      ico=0
      d=0.0d0
      return
      end
      
      subroutine rg(fo,ic,nat,z,r,lzmat)
c     read first geometry
      implicit none
      integer*4 ic,nat,z(*),I,l,ix
      real*8 r(*),q
      character*(*) fo
      character*80 s80
      logical lzmat
      open(2,file=fo)
1     read(2,2000,end=99,err=99)s80
2000  format(a80)
      IF((lzmat.and.(s80(19:39).EQ.'Z-Matrix orientation:'.OR.
     1               s80(26:46).EQ.'Z-Matrix orientation:'.OR.
     1               s80(20:37).EQ.'Input orientation:'.OR.
     1               s80(27:44).EQ.'Input orientation:'))
     1    .OR.
     1    ((.not.lzmat).and.
     2              (s80(20:40).EQ.'Standard orientation:'.OR.
     2               s80(26:46).EQ.'Standard orientation:')))THEN
       DO 2004 I=1,4
2004   READ(2,*)
       if(ic.eq.0)then
        l=0
2005    READ(2,2000)s80
        IF(s80(2:4).NE.'---')THEN
         l=l+1
         BACKSPACE 2
         READ(2,*)q,q
         IF(q.EQ.-1)l=l-1
         GOTO 2005
        ENDIF
        nat=l
       else
        do 2 l=1,nat
   2    read(2,*)z(l),z(l),r(1+3*(l-1)),(r(ix+3*(l-1)),ix=1,3)
        write(6,*)nat,' atoms read'
       endif
       close(2)
       return
      endif
      goto 1
99    nat=0
      close(2)
      call report('geometry not found')
      end

      subroutine rgf(nat,z,r,lzmat,ig)
c     read the next geometry
      implicit none
      integer*4 nat,z(*),I,l,ix,ig
      real*8 r(*)
      character*80 s80
      logical lzmat
1     read(2,2000,end=99,err=99)s80
2000  format(a80)
      IF((lzmat.and.(s80(19:39).EQ.'Z-Matrix orientation:'.OR.
     1               s80(26:46).EQ.'Z-Matrix orientation:'.OR.
     1               s80(20:37).EQ.'Input orientation:'.OR.
     1               s80(27:44).EQ.'Input orientation:'))
     1    .OR.
     1    ((.not.lzmat).and.
     2              (s80(20:40).EQ.'Standard orientation:'.OR.
     2               s80(26:46).EQ.'Standard orientation:')))THEN
       DO 2004 I=1,4
2004   READ(2,*)
       do 2 l=1,nat
   2   read(2,*)z(l),z(l),r(1+3*(l-1)),(r(ix+3*(l-1)),ix=1,3)
       write(6,*)nat,' atoms read'
       ig=ig+1
       return
      endif
      goto 1
99    nat=0
      return
      end

      subroutine mkg(n,g,p,r,r0,icav,cc)
c     read tesserea nad make p-gradient
      implicit none
      integer*4 n,np,i,j,isp,nat,icav
      real*8 g(n),x,y,z,s,st,bohr,xij,yij,zij,d,r(n),p,rc,cc,r0(n),dd
      real*8,allocatable::rs(:)
      character*80 s80
      nat=n/3
      allocate(rs(nat))
      bohr=0.529177d0
      np=0
      isp=0
      st=0.0d0
      g=0.0d0

1     read(2,80,end=88,err=88)s80
 80   format(a80)
      if(s80(4:14).eq.'ISph     Re'.or.
     1   s80(4:13).eq.'ISph    Re')then
       isp=nat
       do 3 i=1,isp
3      read(2,*)rs(i),rs(i)
      endif
      if(s80(2:26).eq.'GePol: Number of points  ')read(s80(56:62),*)np
      if(s80(4:28).eq.'ITs           Coordinates'.and.np.ne.0.
     1and.isp.ne.0)then
       read(2,*)
       read(2,*)
       write(6,603)np
603    format(' Number of points:',i6)
       do 2 i=1,np
       read(2,*)x,x,y,z,s
c      which atom:
       do 4 j=1,nat
       rc=rs(j)
       xij=x-r(1+3*(j-1))
       yij=y-r(2+3*(j-1))
       zij=z-r(3+3*(j-1))
       d=dsqrt(xij**2+yij**2+zij**2)
4      if(dabs(d-rc).lt.1.0d-3)goto 41
       write(6,604)i,x,y,z
604    format(i6,3f12.6)
       call report('unassigned tile')
41     st=st+s
c      pressure force:
       g(1+3*(j-1))=g(1+3*(j-1))-xij*s*p/bohr**2/rc
       g(2+3*(j-1))=g(2+3*(j-1))-yij*s*p/bohr**2/rc
       g(3+3*(j-1))=g(3+3*(j-1))-zij*s*p/bohr**2/rc
c      distance-dependent, ad hoc term
       dd=d-dsqrt((x-r0(1+3*(j-1)))**2+
     1            (y-r0(2+3*(j-1)))**2+
     1            (z-r0(3+3*(j-1)))**2)
       g(1+3*(j-1))=g(1+3*(j-1))-xij*s*cc/bohr**2/rc*dd
       g(2+3*(j-1))=g(2+3*(j-1))-yij*s*cc/bohr**2/rc*dd
2      g(3+3*(j-1))=g(3+3*(j-1))-zij*s*cc/bohr**2/rc*dd
       write(6,600)st
600    format(' Cavity surface area (A^2):',f10.2)
c      call wg(' Pressure force:',n,g)
       call pg(n,g,r)
c      call wg(' Projected:',n,g)
       icav=icav+1
       return
      endif
      goto 1
 88   nat=0
      end

      subroutine pg(n,g,r)
c     project translations/rotations from gradient
      implicit none
      integer*4 n,i,j
      real*8 g(n),r(n)
      real*8,allocatable::A(:,:),t(:)
      allocate(A(n,n),t(n))
      call DOMA(A,n,r)
      t=0.0d0
      do 1 i=1,N
      do 1 j=1,N
1     t(i)=t(i)+A(i,j)*g(j)
      g=t
      return
      end

      subroutine wg(s,n,g)
      implicit none
      character*(*)s
      integer*4 n,i
      real*8 g(n)
      write(6,'(a)')s
      do 5 i=1,n/3
5     write(6,602)i,g(1+3*(i-1)),g(2+3*(i-1)),g(3+3*(i-1))
602   format(i6,3f12.6)
      return
      end

      SUBROUTINE DOMA(A,N,r)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(N,N),r(N)
      real*8,allocatable::C(:,:),TEM(:,:)
C
      allocate(C(3,N/3),TEM(N,N))
      AMACH=0.00000000000001d0
      NAT=N/3
      do 11 I=1,NAT
      do 11 J=1,3
11    C(J,I)=r(J+3*(I-1))
      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,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

      SUBROUTINE ORT(A,N6,NAT)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      DIMENSION A(3*NAT,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 WRITEFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      open(20,file='FILE.NEW.FC')
      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)
      close(20)
      RETURN
      END

      SUBROUTINE READFF(N,FCAR)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(N,N)
      OPEN(20,FILE='FILE.FC',STATUS='OLD')
      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)
      DO 31 I=1,N
      DO 31 J=I+1,N
   31 FCAR(I,J)=FCAR(J,I)
      CLOSE(20)
      WRITE(*,*)' Cartesian FF read in ... '
      RETURN
      END

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