PROGRAM FIELDFIX
      implicit none
      INTEGER*4 NOAT,NQ,N,i,j,ix,ia,iq,ISYM,NOSC
      REAL*8 :: df,CM,AMU
      REAL*8,allocatable::SV(:,:,:),R(:,:),FRQ(:),s(:,:),FF(:,:),
     1ZIM(:),ZM(:)
      character*4 JOB
      character*20 s20
      df=100.0d0
      CM=219474.63d0
      AMU=1822.0d0

      write(6,601)
 601  format(' Projection of negative frequencies from a froce field',/,
     1 /,    '      Input: F.INP',/,
     1       '             FTRY.INP',/,  
     1       '             FILE.FC  ... will be overwritten!',/,/)  
      if(iargc().gt.0)then
       call getarg(1,s20)
       read(s20,*)df
      endif
      write(6,600)df
 600  format(' Small frequencies will be changed to',f8.2,' cm-1')

c     dimensions from S-matrix:
      call readsd(NQ,NOAT)
      N=3*NOAT
      allocate(SV(NOAT,NQ,3),R(3,NOAT),FRQ(NQ),ZIM(NOAT),ZM(N),
     1s(N,N),FF(N,N))
      call reads(NQ,NOAT,SV,FRQ,R)
C     Read-in scalefactors etc from FTRY.INP:
      call readinp(JOB,ISYM,NOSC,NOAT,ZIM,ZM)
      if(NQ.ne.N)then
       write(6,*)'NQ <> N, go to hell'
       stop
      endif

c     reduced s-matrix, s^t.s = 1 (s^t=s^-1):
      s=0.0d0
      do 1 ia=1,NOAT
      do 1 ix=1,3
      i=ix+3*(ia-1)
      do 1 iq=1,NQ
  1   s(i,iq)=SV(ia,iq,ix)*ZM(i)
      write(6,*)'s done'

c     replace low/negative frequencies,then replace them by w^2 in au:
      j=0
      do 2 i=1,NQ
      if(FRQ(i).lt.df)then
       j=j+1
       FRQ(i)=df
      endif
  2   FRQ(i)=(FRQ(i)/CM)**2
      write(6,602)j
602   format(i6,' frequencies replaced')
      
c     new mass-weighted force field  st w^2 s:
      do 3 i=1,N
      do 3 j=1,N
      FF(i,j)=0.0d0
      do 3 iq=1,N
3     FF(i,j)=FF(i,j)+s(i,iq)*FRQ(iq)*s(j,iq)
      write(6,*)'FF done'

c     force field in au:
      do 4 i=1,N
      do 4 j=1,N
4     FF(i,j)=FF(i,j)*ZM(i)*ZM(j)*AMU
      call WRITEFF(NOAT,FF,'FILE.FC')
      END 

      subroutine reads(NQ,NOAT,SV,FRQ,R)
      implicit none
      integer*4 NQ,NOAT,K,I,J
      real*8 SV(NOAT,NQ,3),R(3,NOAT),FRQ(*),BOHR
      BOHR=0.52917715d0
      OPEN(17,FILE='F.INP',STATUS='OLD')
      READ(17,*)NQ,NOAT,NOAT
      DO 3 K=1,NOAT
      READ(17,*)R(1,K),(R(I,K),I=1,3)
      DO 3 I=1,3
3     R(I,K)=R(I,K)/BOHR
      WRITE(*,*)' Geometry read in'
      READ(17,*)
      DO 4 K=1,NOAT
      DO 4 I=1,NQ
4     READ (17,*)(SV(K,I,J),J=1,2),(SV(K,I,J),J=1,3)
      READ(17,*)
      READ(17,*)(FRQ(I),I=1,NQ)
c1717 FORMAT(3F15.7)
      WRITE(*,*)' S-matrix read in'
      return
      end
     
      subroutine readsd(NQ,NOAT)
      implicit none
      integer*4 NQ,NOAT
      OPEN(17,FILE='F.INP',STATUS='OLD')
      READ(17,*)NQ,NOAT,NOAT
      close(17)
      return
      end


      subroutine readinp(JOB,ISYM,NOSC,NOAT,ZIM,ZM)
      implicit none
      character*4 JOB
      integer*4 ISYM,NOSC,NOAT,I,NMOL,J
      real*8 CHATOT,ZIM(*),BM,ZM(*)
      OPEN(15,FILE='FTRY.INP')
      READ(15,1500)JOB,ISYM
1500  FORMAT(A4,I4)
      READ(15,*)NOSC
      DO 222 I=1,NOSC
222   READ(15,*) 
1040  FORMAT(6G12.6)
      WRITE(*,*)' Scale factors read in...'
C
C     Read-in atomic charges for FPC:
c     READ(15,1040) (CHR(I),I=1,NOAT)
      READ(15,*) (ZIM(I),I=1,NOAT)
      CHATOT=0.0d0
      DO 77 I=1,NOAT
77    CHATOT=CHATOT+ZIM(I)
      WRITE(*,777)CHATOT,NOAT
777   FORMAT(' TOTAL CHARGE',F20.6,/I6,' atoms')
C
C     Read-in atomic masses:
      READ(15,*)NMOL
      IF(NMOL.NE.1)THEN
       CLOSE(15)
c       call report(' Only one molecule allowed here !')
      ENDIF
      READ(15,1040)(ZIM(I),I=1,NOAT)
      BM=0.0d0
      DO 102 I=1,NOAT
102   BM=BM+ZIM(I)
      WRITE(*,778)BM
778   FORMAT(' TOTAL MASS  ',F20.6)
      DO 220 I=1,NOAT
      DO 220 J=1,3
220   ZM(3*I-3+J)=DSQRT(ZIM(I))
C
      CLOSE(15)
      WRITE(*,*)' FTRY.INP READ IN'
      END
      
      SUBROUTINE WRITEFF(NA,FCAR,s)
      IMPLICIT INTEGER*4 (I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION FCAR(3*NA,3*NA)
      character*(*) s
      open(2,FILE=s)
      N=NA*3
      N1=1
1     N3=N1+4
      IF(N3.GT.N)N3=N
      DO 130 LN=N1,N
130   WRITE(2,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(2)
      RETURN
      END