PROGRAM GAR_CP2K
C     this program reads vibrational output of CP2K
c     and saves Hessian to FILE.FC and geometry to FILE.X
      IMPLICIT none
      integer*4 NT0,MX3
      PARAMETER (NT0=500,MX3=3*NT0)
      integer*4 nat,n,i,kind(NT0),n1,n2,j,k
      real*8 x(NT0,3),e(MX3),sm(MX3,MX3),CONSTS,zmass(NT0),d(MX3)
      character*80 st
      character*22 s22
      character*10 s10

c     c1=1302.828
      CONSTS=4.359828d0/0.5291772d0/0.5291772d0
c
      open(88,file='vibanal.out') 
1     read(88,880,end=99,err=99)st
880   format(a80)

      if(st(29:38).eq.' - Atoms: ')then
       read(st(61:80),*)nat
       n=3*nat
       write(6,*)nat,'atoms'
       if(nat.gt.NT0)then
        write(6,*)'too many atoms'
        stop
       endif
      endif

      if(st(1:38).eq.' VIB| Hessian in cartesian coordinates')then
       write(6,*)'Hessian found'
       n1=1
102    n2=min(n1+4,n)
       read(88,*)
       read(88,*)
       do 7 i=1,n
7      read(88,2222)(sm(i,j),j=n1,n2)
2222   format(12x,5f13.6)
       if(n2.lt.n)then
        n1=min(n1+5,n)
        goto 102
       endif
       do 8 i=1,n 
       do 8 j=1,n
8      sm(i,j)=sm(i,j)*dsqrt(zmass((i+2)/3)*zmass((j+2)/3))/548.58d0
       open(20,file='FILE.FC')
       call WRITEFF(MX3,N,sm)
       close(20)
       write(6,*)'FILE.FC written'
      endif

      if(st(1:21).eq.'  Atom  Kind  Element')then
       read(88,880)st
       do 2 i=1,nat
       read(88,880)st
       read(st(19:21),*)kind(i)
       read(st(68:80),*)zmass(i)
       read(st(22:33),*)x(i,1)
       read(st(34:45),*)x(i,2)
2      read(st(46:57),*)x(i,3)
       open(89,file='FILE.X')
       write(89,*)'CP2K'
       write(89,*)nat
       do 90 i=1,nat
90     write(89,890)kind(i),(x(i,k),k=1,3)
890    format(i3,3f12.6,' 0 0 0 0 0 0 0 0.0')
       close(89)
       write(6,*)'FILE.X written'
      endif

      if(st(1:67).eq.
     1' VIB|                        NORMAL MODES - CARTESIAN DISPLACEMEN
     2TS')then

       n1=1
22     n2=min(n1+2,n)
       read(88,*)
       read(88,*)
       read(88,2211,end=33,err=33)(e(i),i=n1,n2)
2211   format(22x,f13.6,2f21.6)
       read(88,2211,end=33,err=33)(d(i),i=n1,n2)
       read(88,*)
       read(88,*)
       read(88,*)
       do 4 i=1,nat
4      read(88,2212)((sm(3*(i-1)+k,j),k=1,3),j=n1,n2)
2212   format(17x,3(3x,3f6.2))
       read(88,*)
       if(n2.lt.n)then 
        n1=min(n1+3,n)
        goto 22
       endif
33     n=n2
      endif

      goto 1

99    close(88)

      open(85,file='D.TAB')
      write(85,891)
891   format('Spectrum from CP2K',/)
      write(85,892)
892   format(60(1h-))
      do 10 i=1,n
10    write(85,893)i,e(i),d(i)
893   format(i8,2f12.6)
      write(85,892)
      close(85)
      write(6,*)'D.TAB written'

      open(88,file='F.INP')
      write(88,881)n,n,nat
881   format(3I10)
      do 5 i=1,nat
5     write(88,882)kind(i),(x(i,k),k=1,3)
882   format(I3,3f10.5)
      write(88,*)' Atom Mode X-disp y-disp z-disp'
      do 6 i=1,nat
      do 6 j=1,n
6     write(88,883)i,j,(sm(3*(i-1)+k,j),k=1,3)
883   format(2i5,3f12.6)
      write(88,*)n
      write(88,884)(e(j),j=1,n)
884   format(6f12.3)
      close(88)
      write(6,*)'F.inp written'
      write(6,*)nat,' atoms'
      
      STOP
      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
      CONST=1.0d0
      DO 6 I=1,N
      DO 6 J=1,N
6     FCAR(I,J)=FCAR(I,J)/CONST
      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