program bderq implicit none integer*4 nat,i,n,ii,idiff,ic,nm real*8 step,bohr,al(9),al0(9),a(27),a0(27),g(9),g0(9), 1ff,f1,SFACR real*8,allocatable::aq(:),alq(:),gq(:),r(:,:),r0(:,:),s(:,:), 1e(:),mass(:),stepq(:) character*4 key integer*4,allocatable::iqlist(:) write(6,600) 600 format(/, 1' Polarizability derivatives from polarizabilities',/, 2' at displaced geometries',/,/, 3' Input : PMZ.PAR',/, 4' F.INP',/, 4' */POL.TTT etc.',/, 4' Output: FILE.Q.TTT',/) step=0.05d0 SFACR=0.02342179d0 bohr=0.529177d0 idiff=2 open(7,FILE='PMZ.PAR',status='old') 72 read(7,7000,end=700,err=700)key 7000 format(a4) if(key.eq.'STEP')read(7,*)step if(key.eq.'IDIF')read(7,*)idiff if(key(1:2).eq.'IC')read(7,*)ic if(key(1:3).eq.'NMD')then read(7,*)nm allocate(iqlist(nm)) do 4 i=1,nm 4 read(7,*)iqlist(i) endif goto 72 700 close(7) if(idiff.ne.2)call report('idiff<>2') if(ic.ne.1)call report('ic<>1') write(6,601)nm,nm*2+1 601 format(i4,' modes,',i4,' points') allocate(alq(9*nm),gq(9*nm),aq(27*nm),stepq(nm)) call vz(alq, 9*nm) call vz(gq , 9*nm) call vz(aq ,27*nm) c find out number of atoms open(7,file='1/GEO.X',status='old') read(7,*) read(7,*)nat close(7) n=3*nat allocate(r0(nat,3),r(nat,3),s(n,n),mass(nat),E(n)) c read S and e, transform into a.u.: call reads(S,n,E,mass) c read zero point call rd(al0,g0,a0,r0,nat,0) write(6,*)'zero point read' c start to read the shifted geometries: c DIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIF do 1 IC=1,2*nm I=IC if(IC.GT.nm)I=IC-nm if(IC.le.nm)then ff=1.0d0 else ff=-1.0d0 endif call rd(al,g,a,r,nat,IC) c check, that we are at the right point, eventually extract step: step=0.0d0 call chkgeos(r0,nat,step,S,r,iqlist(I),mass,E) stepq(I)=step do 1 ii=1,9 alq(9*(I-1)+ii ) =alq( 9*(I-1)+ii )+ff*al(ii ) gq( 9*(I-1)+ii ) =gq( 9*(I-1)+ii )+ff*g( ii ) aq(27*(I-1)+ii ) =aq( 27*(I-1)+ii )+ff*a( ii ) aq(27*(I-1)+ii+ 9) =aq( 27*(I-1)+ii+ 9)+ff*a( ii+ 9) 1 aq(27*(I-1)+ii+18) =aq( 27*(I-1)+ii+18)+ff*a( ii+18) c DIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIF do 2 I=1,nm f1=0.5d0* SFACR*bohr/stepq(I) do 2 ii=1,9 alq(9*(I-1)+ii ) =alq( 9*(I-1)+ii )*f1 gq( 9*(I-1)+ii ) =gq( 9*(I-1)+ii )*f1 aq(27*(I-1)+ii ) =aq( 27*(I-1)+ii )*f1 aq(27*(I-1)+ii+ 9) =aq( 27*(I-1)+ii+ 9)*f1 2 aq(27*(I-1)+ii+18) =aq( 27*(I-1)+ii+18)*f1 call WRITETTTQ(nm,ALQ,AQ,GQ,E,iqlist) end c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE vz(v,n) real*8 v(*) integer*4 n,i do 1 i=1,n 1 v(i)=0.0d0 return end c ============================================================ subroutine report(s) character*(*) s write(6,*)s stop end c ============================================================ subroutine rd(al,g,a,r,nat,io) implicit none integer*4 nat,io,ia,ix,iy,iz,is real*8 al(*),g(*),a(*),r(nat,3) character*10 s10 write(s10,10)io+1 10 format(i10) do 1 is=1,len(s10) 1 if(s10(is:is).ne.' ')goto 2 write(6,*)s10(is:len(s10))//'/POL.TTT' 2 open(9,file=s10(is:len(s10))//'/POL.TTT') read(9,*) read(9,*) read(9,*)((al(IY+(IX-1)*3),IY=1,IX),IX=1,3) AL(2+(1-1)*3)=AL(1+(2-1)*3) AL(3+(1-1)*3)=AL(1+(3-1)*3) AL(3+(2-1)*3)=AL(2+(3-1)*3) close(9) open(9,file=s10(is:len(s10))//'/GP.TTT') read(9,*) read(9,*) read(9,*)((g(IY+(IX-1)*3),IY=1,3),IX=1,3) close(9) open(9,file=s10(is:len(s10))//'/A.TTT') read(9,*) read(9,*) do 3 iZ=1,3 read(9,*)((A(IZ+3*(IY-1)+9*(IX-1)),IY=1,IX),IX=1,3) A(IZ+3*(2-1)+9*(1-1))=A(IZ+3*(1-1)+9*(2-1)) A(IZ+3*(3-1)+9*(1-1))=A(IZ+3*(1-1)+9*(3-1)) 3 A(IZ+3*(3-1)+9*(2-1))=A(IZ+3*(2-1)+9*(3-1)) close(9) open(9,file=s10(is:len(s10))//'/GEO.X') read(9,*) read(9,*) do 4 ia=1,nat 4 read(9,*)r(ia,1),(r(ia,ix),ix=1,3) close(9) return end c ============================================================ SUBROUTINE READS(S,N,E,mass) implicit none integer*4 N,NM,nat,I,IA,K,iq,iz real*8 S(N,N),E(*),cm,SFACR,mass(*),amave cm=219470.0d0 SFACR=0.02342179d0 call vz(E,N) call vz(s,N*N) OPEN(7,FILE='F.INP') READ(7,*)NM,N,nat DO 1 I=1,nat READ(7,*)iz 1 mass(I)=amave(iz) READ(7,*) DO 3 IA=1,NAT DO 3 I=1,NM 3 READ(7,*)(S(3*(IA-1)+K,I),K=1,2),(S(3*(IA-1)+K,I),K=1,3) READ(7,*) READ(7,*)(E(I),I=1,NM) CLOSE(7) do 9 iq=1,nm E(iq)=E(iq)/cm do 9 i=1,nat do 9 k=1,3 9 S(3*(i-1)+k,iq)=S(3*(i-1)+k,iq)*SFACR RETURN END c ============================================================ subroutine chkgeos(r0,nat,step,S,r,ireq,mass,E) implicit none integer*4 nat,l,ic,ix,iq,ii,ireq real*8 r0(nat,3),step,r(nat,3),dq,S(3*nat,3*nat),CM,AM, 1SFACR,bohr,geotol,mass(*),E(*) real*8,allocatable:: rc(:,:) allocate(rc(nat,3)) c CM=219470.0d0 SFACR=0.02342179d0 bohr=0.529177d0 geotol=1.0d-4 ic=0 c c Cartesian displacements from equilibrium DO 2 l=1,nat DO 2 ix=1,3 2 rc(l,ix)=r(l,ix)-r0(l,ix) c c normal mode displacements do 3 iq=1,3*nat dq=0.0d0 do 4 l=1,nat AM=mass(l) do 4 ix=1,3 ii=3*(l-1)+ix 4 dq=dq+AM*S(ii,iq)*rc(l,ix) c c adapt to the strange pmz units dq=dq/SFACR: dq=dq/SFACR if(dabs(dq).gt.geotol)then write(6,6000)iq,e(iq)*CM,dq,dq/(SFACR*bohr) 6000 format('mode',i4,f10.2,' cm-1',' dQ =',f10.4,' au:',f10.4) if(ireq.ne.iq)call report('errorneous mode numbering') ic=ic+1 step=dabs(dq) endif 3 continue c if(ic.ne.1) call report('too many modes alike') return end c ============================================================ function amave(iz) c atomic mass of element with atomiz number iz c approximate - to be compatible with new2 implicit none integer*4 iz,MENDELEV PARAMETER (MENDELEV=89) real*8 amave,amas(MENDELEV) data amas/1.007825d0,4.003, 2 6.941, 9.012, 10.810,12.0d0,14.003074d0,15.999d0, 2 18.9984d0, 20.179, 3 22.990,24.305, 26.981,28.086,30.9737d0, 3 31.972d0,34.968852d0,39.948, 4 39.098,40.080,44.956,47.900,50.941,51.996,54.938,55.847, 4 58.933,58.700,63.546,65.380, 4 69.720,72.590,74.922,78.960,79.904,83.800, 5 85.468,87.620,88.906,91.220,92.906,95.940,98.906,101.070, 5 102.906,106.400,107.868,112.410, 5 114.82,118.69,121.75,127.600,126.905,131.300, 6 132.905,137.330,138.906, 6 140.120,140.908,144.240,145.000,150.400, 6 151.960,157.250,158.925,162.500,164.930,167.260,168.934, 6 173.040,174.970, 6 178.490,180.948,183.850,186.207,190.200,192.220,195.090, 6 196.967,207.590,204.370,207.200,208.980,210.000,210.001, 6 222.02, 7 223.000,226.025,227.028/ c if(iz.lt.0.or.iz.gt.MENDELEV)then write(6,*)iz call report('Unknown atomic number') else amave=amas(iz) return endif end C ========================================== SUBROUTINE WRITETTTQ(N,ALQ,AQ,GQ,E,iqlist) implicit none integer*4 N,im,I,J,K,iqlist(*) real*8 ALQ(*),AQ(*),GQ(*),E(*),CM real*8,allocatable::ALPHA(:,:,:),G(:,:,:),A(:,:,:,:) allocate(ALPHA(N,3,3),G(N,3,3),A(N,3,3,3)) do 101 im=1,N do 101 I=1,3 do 101 J=1,3 ALPHA(im,I,J)=ALQ( 9*(im-1)+I+3*(J-1)) G(im,I,J)=GQ( 9*(im-1)+I+3*(J-1)) do 101 K=1,3 101 A(im,I,J,K)=AQ(27*(im-1)+I+3*(J-1)+9*(K-1)) c CM=219470.0d0 open(22,file='FILE.Q.TTT') WRITE(22,2003)N 2003 FORMAT(' ROA tensors, normal modes derivatives',/,I4,' modes',/, 1' The electric-dipolar electric-dipolar polarizability:',/, 2 ' mode e(cm-1) jx jy jz') DO 1 im=1,N write(22,221)im,E(iqlist(im))*CM 221 format(i5,f12.2) DO 1 I=1,3 1 WRITE(22,2001)I,(ALPHA(im,I,J),J=1,3) 2001 FORMAT(I3,3g14.6) WRITE(22,2004) 2004 FORMAT(' The electric dipole magnetic dipole polarizability:',/, 1 ' mode e(cm-1) jx(Bx) jy(By) jz(Bz)') DO 2 im=1,N write(22,221)im,E(iqlist(im))*CM DO 2 I=1,3 2 WRITE(22,2001)I,(G(im,J,I),J=1,3) WRITE(22,2005) 2005 FORMAT(' The electric dipole electric quadrupole polarizability:', 2/, ' mode e(cm-1) kx ky kz') DO 3 im=1,N write(22,221)im,E(iqlist(im))*CM DO 3 I=1,3 DO 3 J=1,3 3 WRITE(22,2002)I,J,(A(im,I,J,K),K=1,3) 2002 FORMAT(2I3,3g14.6) WRITE(22,2009) 2009 FORMAT('The velocity form of alpha',/, 1' mode e(cm-1) vx vy vz') DO 5 im=1,N write(22,221)im,E(iqlist(im))*CM DO 5 I=1,3 5 WRITE(22,2001)I,(ALPHA(im,I,J),J=1,3) close(22) RETURN END