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 nat,n,i,n1,n2,j,k,ix,iy,iz,iat,ia,iargc real*8 step,ss,aa(3,3),z0,rx,ry,s1,s0,doc,dcpI,dcpII real*8,allocatable::r(:),p(:),u(:),x(:,:),e(:),sm(:,:), 1d(:),zmass(:),alpha(:),Gt(:),At(:) integer*4,allocatable::kind(:) character*80 st c z0=0.0d0 iz=0 step=z0 iat=0 if(iargc().gt.0)then call getarg(1,st) open(88,file=st,status='old') else open(88,file='vibanal.out',status='old') endif 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 allocate(r(n),u(n),p(n),x(nat,3),e(n),sm(n,n),zmass(nat), 1 kind(nat),d(n),alpha(9*n),Gt(9*n),At(27*n)) zmass=z0 alpha=z0 Gt=z0 At=z0 d=z0 r=z0 p=z0 u=z0 e=z0 write(6,*)nat,'atoms' endif if(st(25:33).eq.'step_size')read(st(34:len(st)),*)step if(st(7:17).eq.'REPLICA Nr.'.and.step.ne.z0)then read(st(56:60),*)ia if(st(75:75).eq.'X')iz=1 if(st(75:75).eq.'Y')iz=2 if(st(75:75).eq.'Z')iz=3 if(iz.eq.0)call report('iz not defined') iat=iz+3*(ia-1) ss=1.0d0 if(st(77:77).eq.'-')ss=-1.0d0 endif if(st(9:36).eq.'Polarizability tensor [a.u.]')then if(iz.eq.0)call report('iz not defined') if(iat.eq.0)call report('iat not defined') read(88,880,end=99,err=99)st read(st(18:len(st)),*)aa(1,1),aa(2,2),aa(3,3) read(88,880,end=99,err=99)st read(st(18:len(st)),*)aa(1,2),aa(1,3),aa(2,3) read(88,880,end=99,err=99)st read(st(18:len(st)),*)aa(2,1),aa(3,1),aa(3,2) aa(1,2)=(aa(1,2)+aa(2,1))/2.0d0 aa(2,1)=aa(1,2) aa(1,3)=(aa(1,3)+aa(3,1))/2.0d0 aa(3,1)=aa(1,3) aa(2,3)=(aa(2,3)+aa(3,2))/2.0d0 aa(3,2)=aa(2,3) c record only the second value - already the difference: if(ss.lt.z0)then do 11 ix=1,3 do 11 iy=1,3 i=iy+3*(ix-1)+9*(iat-1) 11 alpha(i)=aa(ix,iy) endif if(ia.eq.nat.and.iz.eq.3.and.ss.lt.z0)then write(6,*)'End of numerical differentiation for alpha reached' alpha=alpha/2.0d0/step call WRITETTT(n,alpha,Gt,At) endif endif if(st(8:15).eq.'Vector a')then open(9,file='CELL.TXT') write(9,*)st(28:59) read(88,880,end=99,err=99)st write(9,*)st(28:59) read(88,880,end=99,err=99)st write(9,*)st(28:59) close(9) write(6,*)'CELL.TXT' 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(N,sm) close(20) write(6,*)'FILE.FC written' endif if(st(4:20).eq.'Atom Kind Element')then do 2 i=1,nat read(88,880)st 2 read(st(16:len(st)),*)kind(i),x(i,1),x(i,2),x(i,3),zmass(i), 1 zmass(i) 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) 101 read(88,880,end=99,err=99)st if(st(6:14).eq.'Frequency')read(st(24:len(st)),*)(e(i),i=n1,n2) c if(st(6:14).eq.'Frequency')read(st(24:len(st)),*)(d(i),i=n1,n2) if(st(6:10).eq.'Raman' )read(st(24:len(st)),*)(r(i),i=n1,n2) if(st(15:20).eq.'io (P)' )read(st(24:len(st)),*)(p(i),i=n1,n2) if(st(15:20).eq.'io (U)' )read(st(24:len(st)),*)(u(i),i=n1,n2) if(st(3:6).ne.'ATOM')goto 101 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 n=n2 endif goto 1 99 close(88) open(85,file='DOG.TAB') write(85,891) 891 format('Spectrum from CP2K',/) write(85,892) 892 format(60(1h-)) do 10 i=1,n 10 if(e(i).gt.5.0d0)write(85,893)i,e(i),d(i),d(i) 893 format(i8,f12.6,12G12.4) write(85,892) close(85) write(6,*)'DOG.TAB written' open(85,file='ROA.TAB') write(85,3304) 3304 FORMAT( 1'MODE FREQ Ramx (180) Ramy',6x,'Ramtot',11x,'CID90', 27x,'CID-X',6x,'CID180',3x,'DEP',5x,'ROA180',5x,'DOC', 15x,'RamDCPI',5x,'RoaDCPI',3x,'RamnDCPII',4x,'IRR(180)', 14x,'IRL(180)',/,10x,'cm-1') write(85,892) do 110 i=1,n rx=r(i)/(1.0d0+p(i)) ry=r(i)*p(i)/(1.0d0+p(i)) s1=(2.0d0*ry+rx)/60.0d0 s0=(3.0d0*rx-4.0d0*ry)/60.0d0 doc=(3.0d0*p(i)-1.0d0)/(p(i)+1.0d0) dcpI =6.0d0*(6.0d0*s1-2.0d0*s0) dcpII=6.0d0*( s1+3.0d0*s0) 110 if(e(i).gt.5.0d0)write(85,3000)i,e(i),rx,ry,r(i),z0,z0,z0,p(i),z0, 1doc,dcpI,z0,dcpII,z0,z0 3000 FORMAT(I5,f9.2,6e12.4,f6.3,e12.4,f7.3,7e12.4) write(85,892) close(85) write(6,*)'ROA.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(N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION FCAR(N,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 subroutine report(s) character*(*) s write(6,*)s stop end SUBROUTINE WRITETTT(N,ALPHA,G,A) IMPLICIT none INTEGER*4 N,I,J,K,L,IX,IIND,NAT,ii real*8 ALPHA(9*N),G(9*N),A(27*N),p p=3.0d0/2.0d0 NAT=N/3 OPEN(2,FILE='FILE.TTT') WRITE(2,2000)NAT,1,0.0d0 2000 FORMAT(' ROA tensors, cartesian derivatives',/, 1I4,' atoms, freq. ',i2,f11.6/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2' Atom/x jx jy jz') DO 1 I=1,3 WRITE(2,2002)I 2002 FORMAT(' Alpha(',I1,',J):') DO 1 L=1,NAT DO 1 IX=1,3 IIND=3*(L-1)+IX 1 WRITE(2 ,2001)L,IX,(ALPHA(I+3*(J-1)+9*(IIND-1)),J=1,3) 2001 FORMAT(I5,1H ,I1,6g15.7) WRITE(2,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' Atom/x jx(Bx) jy(By) jz(Bz)') DO 2 I=1,3 WRITE(2,2003)I 2003 FORMAT(' G(',I1,',J):') DO 2 L=1,NAT DO 2 IX=1,3 IIND=3*(L-1)+IX 2 WRITE(2 ,2001)L,IX,(G(I+3*(J-1)+9*(IIND-1)),J=1,3) WRITE(2,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' Atom/x kx ky kz') DO 3 I=1,3 DO 3 J=1,3 WRITE(2,2006)I,J 2006 FORMAT(' A(',I1,',',I1,',K):') DO 3 L=1,NAT DO 3 IX=1,3 ii=3*(J-1)+9*(I-1)+27*(3*(L-1)+IX-1) 3 WRITE(2,2007)L,IX,(A(K+ii)*p,K=1,3),L,IX,I,J 2007 FORMAT(I5,1H ,I1,3g15.7,' ',4i3) write(2,*) write(2,*)'dummy alpha v:' DO 4 I=1,3 WRITE(2,2002)I DO 4 L=1,NAT DO 4 IX=1,3 IIND=3*(L-1)+IX 4 WRITE(2 ,2001)L,IX,(ALPHA(I+3*(J-1)+9*(IIND-1)),J=1,3) CLOSE(2) write(6,*)'FILE.TTT written' RETURN END