PROGRAM GAR_CASTEP 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,MENDELEV PARAMETER (NT0=500,MX3=3*NT0,MENDELEV=89) integer*4 nat,i,kind(NT0),j,k,nat3,nwv,ii,it,iargc, 1ix,jx,kx,igeo real*8 x(NT0,3),e(MX3),s(MX3,MX3),CONSTS,zmass(NT0),d(MX3), 1uv(3,3),a,b,c,ab,ac,bc,alpha,beta,gamma,ev,bohr, 1o11,o12,o22,o23,o33,o13,P1,down,doc, 1f(MX3,MX3),w2,c1,ram(MX3),alphaq(MX3,3,3),alphac(MX3,3,3), 1AMU,YDY,YDX,SPAL,SpA3,beta2,ram1,SAL0,SAL1,z0 character*80 st,fn CHARACTER*2 atsy(MENDELEV),s2 data atsy/' H','He','Li','Be',' B',' C',' N',' O',' F','Ne', 3'Na','Mg','Al','Si',' P',' S','Cl','Ar', 4' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn', 4 'Ga','Ge','As','Se','Br','Kr', 5'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd', 5 'In','Sn','Sb','Te',' I','Xe', 6'Cs','Ba','La', 6 'Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho', 6 'Er','Tm','Yb','Lu', 6'Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg', 6 'Tl','Pb','Bi','Po','At','Rn', 7'Fr','Ra','Ac'/ logical lram,lrd,lsm,l1,l2,lff,lmm lram=.false. lmm=.false. lrd=.false. lsm=.false. lff=.false. c1=1302.828d0 ev=27.2107d0 bohr=0.5291772d0 CONSTS=4.359828d0/bohr**2 nat=0 nat3=0 igeo=0 AMU=1822.0d0 z0=0.0d0 c if(iargc().gt.0)then call getarg(1,fn) else fn='out.phonon' endif open(89,file='FILE.X') open(88,file=fn) 1 read(88,880,end=99,err=99)st 880 format(a80) if(st(30:52).eq.'Raman Intensity Tensors'.or. 1 st(25:52).eq.'Raman Susceptibility Tensors')then lrd=.true. write(6,*)'trying to read raman derivatives' read(88,*) i=0 2 read(88,880,end=99,err=99)st if(st(21:32).eq.'Raman tensor')then i=i+1 if(MX3.lt.i)call report('too many modes') read(88,880,end=99,err=99)st read(st(4:33),*)(alphaq(i,1,j),j=1,3) read(88,880,end=99,err=99)st read(st(4:33),*)(alphaq(i,2,j),j=1,3) read(88,880,end=99,err=99)st read(st(4:33),*)(alphaq(i,3,j),j=1,3) read(88,*) goto 2 endif write(6,*)i,' modes found' endif if(st(1:23).eq.' + N Frequency'.and..not.lram)then lram=.true. read(88,*) read(88,*) i=0 8821 read(88,880)st i=i+1 if(st(6:8).ne.'???'.and.st(6:8).ne.'==='.and.st(6:8).ne.'...') 1 then read(st(10:23),*,err=8822)e(i) read(st(30:46),*,err=8822,end=8822)d(i) read(st(52:70),*,err=8822)ram(i) goto 8821 endif 8822 nat3=i-1 endif if(st(1:31).eq.' + N Frequency irrep. ')then lram=.true. read(88,*) read(88,*) read(88,*) i=0 882 read(88,880)st i=i+1 if(st(6:8).ne.'...')then read(st(10:23),*)e(i) read(st(30:46),*)d(i) read(st(52:70),*)ram(i) goto 882 endif nat3=i endif if(st(1:16).eq.' Number of ions ')read(st(23:27),*)nat if(st(1:16).eq.' Number of branc')read(st(23:27),*)nat3 if(st(1:16).eq.' Number of wavev')then read(st(23:27),*)nwv if(nwv.ne.1)call report('nwv<>1 not implemented') endif if(st(1:16).eq.' Unit cell vecto'.or. 1 st(9:23).eq.'Real Lattice(A)')then do 201 i=1,3 201 read(88,*)(uv(i,j),j=1,3) open(14,file='CELL.TXT') do 202 i=1,3 202 write(14,1414)(uv(i,j),j=1,3) 1414 format(3f12.6) write(6,*)'unit cell read, CELL.TXT written' close(14) endif if(st(26:55).eq.'Total number of ions in cell =') 1read(st(56:len(st)),*)nat if(st(1:11).eq.' Fractional')then if(nat.gt.NT0)call report('too many atoms') call setcp(uv,a,b,c,ab,ac,bc,alpha,beta,gamma, 1 o11,o12,o13,o22,o23,o33) do 203 i=1,nat read(88,881)j,(x(i,j),j=1,3),s2,zmass(i) 881 format(i6,1x,3f12.6,3x,a2,f17.6) if(s2(2:2).eq.' ')then s2(2:2)=s2(1:1) s2(1:1)=' ' endif kind(i)=0 do 204 j=1,MENDELEV 204 if(atsy(j).eq.s2)kind(i)=j x(i,1)=x(i,1)*o11+x(i,2)*o12+x(i,3)*o13 x(i,2)= x(i,2)*o22+x(i,3)*o23 203 x(i,3)= x(i,3)*o33 write(89,*)'CASTEP' write(89,*)nat do 91 i=1,nat 91 write(89,890)kind(i),(x(i,k),k=1,3) write(6,*)'coordinates read' igeo=igeo+1 lmm=.true. endif l1=st(39:69).eq.'Fractional coordinates of atoms' l2=st(44:74).eq.'Fractional coordinates of atoms' if(l1.or.l2)then if(nat.gt.NT0)call report('too many atoms') call setcp(uv,a,b,c,ab,ac,bc,alpha,beta,gamma, 1 o11,o12,o13,o22,o23,o33) read(88,*) read(88,*) do 205 i=1,nat if(l1)read(88,883)s2,(x(i,j),j=1,3) 883 format(14x,a2,19x,3f11.6) if(l2)read(88,884)s2,(x(i,j),j=1,3) 884 format(14x,a2,24x,3f11.6) if(igeo.eq.0)then c read atom types only in the beginning: if(s2(2:2).eq.' ')then s2(2:2)=s2(1:1) s2(1:1)=' ' endif kind(i)=0 do 209 j=1,MENDELEV 209 if(atsy(j).eq.s2)kind(i)=j endif x(i,1)=x(i,1)*o11+x(i,2)*o12+x(i,3)*o13 x(i,2)= x(i,2)*o22+x(i,3)*o23 205 x(i,3)= x(i,3)*o33 write(6,*)'coordinates read' write(89,*)'CASTEP' 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') igeo=igeo+1 endif if(st(6:10).eq.'q-pt='.and..not.lram)then do 206 i=1,nat3 206 read(88,*)j,e(j),d(j),ram(j) lram=.true. write(6,*)nat3,' Raman intensities read' endif if(st(25:43).eq.'Phonon Eigenvectors'.and.nat3.gt.0)then read(88,*) do 215 j=1,nat3 do 215 i=1,nat 215 read(88,*)(s(3*(i-1)+1,j),ix=1,3),(s(3*(i-1)+2,j),ix=1,2), 1 (s(3*(i-1)+3,j),ix=1,2) write(6,*)'s-matrix read',nat,nat3 lsm=.true. endif if(st(13:33).eq.'Force Constant matrix')then nat3=3*nat read(88,*) read(88,*) do 213 i=1,nat do 213 ix=1,3 213 read(88,*)(f(ix+3*(i-1),j),j=1,2),(f(ix+3*(i-1),j),j=1,nat3) write(6,*)' Force field read' lff=.true. write(6,*)nat endif goto 1 99 close(88) if(lsm.and.lmm)then open(99,file='F.INP') write(99,991)nat3,nat3,nat 991 format(3i5) do 208 i=1,nat 208 write(99,990)kind(i),(x(i,j),j=1,3) 990 format(i4,3f10.5) write(99,*)' Atom Mode X-disp y-disp z-disp' do 210 i=1,nat do 210 j=1,nat3 210 write(99,993)i,j,(s(3*(i-1)+k,j)/dsqrt(zmass(i)),k=1,3) 993 format(2i5,3f12.3) write(99,*)nat3 write(99,992)(e(j),j=1,nat3) 992 format(6f11.2) write(6,*)'F.INP written' close(99) endif if(lsm.and.lrd)then write(6,*)'Transforming normal Raman derivatives to Cartesian' do 3 ix=1,3 do 3 jx=1,3 do 3 j=1,nat3 alphac(j,ix,jx)=0.0d0 do 3 i=1,nat3 c ia=(i+2)/3 3 alphac(j,ix,jx)=alphac(j,ix,jx)+s(j,i)*alphaq(i,ix,jx)*14.1d0 c 1 zmass(ia) open(99,file='FILE.TTT') write(99,*) write(99,*)nat write(99,*) write(99,*) do 4 ix=1,3 write(99,*)ix do 4 i=1,nat3 4 write(99,994)(i+2)/3,i-3*((i+2)/3-1),(alphac(i,ix,jx),jx=1,3) 994 FORMAT(I5,1H ,I1,3F15.7) write(99,*) write(99,*) do 5 ix=1,3 write(99,*)ix do 5 i=1,nat3 5 write(99,994)(i+2)/3,i-3*((i+2)/3-1),(0.0d0,jx=1,3) write(99,*) write(99,*) do 6 ix=1,3 do 6 kx=1,3 write(99,*)ix,kx do 6 i=1,nat3 6 write(99,994)(i+2)/3,i-3*((i+2)/3-1),(0.0d0,jx=1,3) write(99,*) write(99,*) do 7 ix=1,3 write(99,*)ix do 7 i=1,nat3 7 write(99,994)(i+2)/3,i-3*((i+2)/3-1),(alphac(i,ix,jx),jx=1,3) close(99) endif if(lff)then do 214 i=1,nat3 do 214 j=1,nat3 214 f(i,j)=f(i,j)/ev*bohr**2 else do 211 i=1,nat3 do 211 j=1,nat3 f(i,j)=0.0d0 do 211 ii=1,nat3 w2=(e(ii)/c1)**2 if(e(ii).lt.0.0d0)w2=-w2 211 f(i,j)=f(i,j)+s(i,ii)*w2*s(j,ii) do 212 i=1,nat3 do 212 j=1,nat3 212 f(i,j)=f(i,j)*dsqrt(zmass((i+2)/3)*zmass((j+2)/3))/CONSTS endif call WRITEFF(MX3,nat3,f) write(6,*)'FILE.FC written',nat3 close(89) write(6,*)'FILE.X written',igeo,' geometries' open(85,file='DOG.TAB') write(85,891) 891 format('Spectrum from CP2K',/) write(85,892) 892 format(60(1h-)) do 10 i=1,nat3 10 write(85,893)i,e(i),d(i) 893 format(i8,5f15.6) write(85,892) close(85) write(6,*)'DOG.TAB written' if(lram)then open(85,file='ROA.TAB') write(85,891) write(85,892) do 11 i=1,nat3 11 if(e(i).gt.1.0d0)write(85,893)i,e(i),ram(i),ram(i),ram(i) write(85,892) close(85) write(6,*)'ROA.TAB written' if(lrd)then open(85,file='ROAi.TAB') write(85,891) write(85,892) do 12 i=1,nat3 SPAL=alphaq(i,1,1)+alphaq(i,2,2)+alphaq(i,3,3) SAL0=SPAL**2 SAL1=alphaq(i,1,1)*alphaq(i,1,1)+alphaq(i,2,2)*alphaq(i,2,2)+ 1 alphaq(i,3,3)*alphaq(i,3,3)+2.0d0*(alphaq(i,1,2)*alphaq(i,1,2)+ 1 alphaq(i,1,3)*alphaq(i,1,3)+ alphaq(i,2,3)*alphaq(i,2,3)) c c polarized backscattering components (YDY+YDX=S180): YDY=3.0d0*SAL1-SAL0 YDX=4.0d0*SAL1+2.0d0*SAL0 if(dabs(YDX).gt.1.0d-9)then P1=YDY/YDX else P1=0.0d0 endif c calculate degree of circularity down= 7.0d0*SAL1+ 1.0d0*SAL0 doc=5.0d0*( 1.0d0*SAL1- 1.0d0*SAL0) if(dabs(down).gt.1.0d-10)doc=doc/down YDY=6.0d0*YDY*(AMU*BOHR**4) YDX=6.0d0*YDX*(AMU*BOHR**4) SpA3=SPAL/3.0d0*dsqrt((AMU*BOHR**4)) beta2=0.5d0*(3*SAL1-SAL0)*(AMU*BOHR**4) ram1=4.0d0*(45.0d0*SpA3**2+7.0d0*beta2) 12 if(e(i).gt.1.0d0)WRITE(85,3000)i,e(i),YDX,YDY,ram1,z0,z0,z0,P1, 1 z0,doc 3000 FORMAT(I5,f9.2,6e12.4,f6.3,e12.4,f7.3,5e12.4) write(85,892) close(85) write(6,*)'ROAi.TAB written' endif endif STOP END SUBROUTINE WRITEFF(MX3,N,FCAR) IMPLICIT INTEGER*4 (I-N) IMPLICIT REAL*8(A-H,O-Z) DIMENSION FCAR(MX3,N) open(20,file='FILE.FC') 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) close(20) RETURN END SUBROUTINE REPORT(RS) CHARACTER*(*) RS WRITE(6,3000) 3000 FORMAT(/,80(1H*)) WRITE(6,*)RS WRITE(6,3001) 3001 FORMAT(80(1H*),/,/,'PROGRAM STOPPED') STOP END subroutine setcp(uv,a,b,c,ab,ac,bc,alpha,beta,gamma, 1o11,o12,o13,o22,o23,o33) implicit none real*8 uv(3,3),a,b,c,ab,ac,bc,alpha,beta,gamma,pi, 1ca,cb,sb,cg,sg,car,sar,o11,o12,o13,o22,o23,o33 c 2,sa,o11i,o12i,v,o13i,o22i,o23i,o33i pi=4.0d0*atan(1.0d0) a=dsqrt(uv(1,1)**2+uv(1,2)**2+uv(1,3)**2) b=dsqrt(uv(2,1)**2+uv(2,2)**2+uv(2,3)**2) c=dsqrt(uv(3,1)**2+uv(3,2)**2+uv(3,3)**2) ab=uv(1,1)*uv(2,1)+uv(1,2)*uv(2,2)+uv(1,3)*uv(2,3) ac=uv(1,1)*uv(3,1)+uv(1,2)*uv(3,2)+uv(1,3)*uv(3,3) bc=uv(2,1)*uv(3,1)+uv(2,2)*uv(3,2)+uv(2,3)*uv(3,3) alpha=180.0d0/pi*acos(bc/b/c) beta=180.0d0/pi*acos(ac/a/c) gamma=180.0d0/pi*acos(ab/b/a) write(6,600)alpha,a,beta,b,gamma,c 600 format(/,' Crystall cell parameters:',/,/, 1 ' alpha: ',f10.2,' a: ',f10.4,/, 1 ' beta: ',f10.2,' b: ',f10.4,/, 1 ' gamma: ',f10.2,' c: ',f10.4,/,/) ca=cos(pi*alpha/180.0d0) c sa=sin(pi*alpha/180.0d0) cb=cos(pi*beta /180.0d0) sb=sin(pi*beta /180.0d0) cg=cos(pi*gamma/180.0d0) sg=sin(pi*gamma/180.0d0) car=(cb*cg-ca)/(sb*sg) sar=dsqrt(1.0d0-car*car) o11=a o12=b*cg o13=c*cb o22=b*sg o23=-c*sb*car o33=c*sb*sar c inverse trafo, not used here: c o11i=1.0d0/a c o12i=-cg/a/sg c v=dsqrt(1.0d0-ca**2-cb**2-cg**2+2.0d0*ca*cb*cg) c o13i=(ca*cg-cb)/a/v/sg c o22i=1.0d0/b/sg c o23i=(cb*cg-ca)/b/v/sg c o33i=sg/c/v return end