program bderq2 implicit none integer*4 n7 parameter (n7=18) integer*4 nat,i,n,ii,idiff,ic,nm,ns real*8 step,bohr,ff,f1,SFACR,g,w real*8,allocatable::aq(:),alq(:),gq(:),r(:,:),r0(:,:),s(:,:), 1e(:),mass(:),stepq(:),es(:),es0(:),esq(:) character*4 key integer*4,allocatable::iqlist(:) write(6,600) 600 format(/, 1' Polarizability derivatives from moments and freqs',/, 2' at displaced geometries',/,/, 3' Input : PMZ.PAR',/, 4' F.INP',/, 4' */JN.SCR.TXT 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 number of electronic states open(55,file='1/JN.SCR.TXT') read(55,*)ns close(55) allocate(es0(n7*ns),es(n7*ns),esq(nm*n7*ns)) call vz(esq,n7*ns*nm) c read zero point call rd(n7,es0,r0,nat,0,g,w) 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(n7,es,r,nat,IC,g,w) 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,n7*ns 1 esq(n7*ns*(I-1)+ii) =esq(n7*ns*(I-1)+ii)+ff*es(ii) c DIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIFDIF do 2 I=1,nm f1=0.5d0* SFACR*bohr/stepq(I) do 2 ii=1,n7*ns 2 esq(n7*ns*(I-1)+ii) =esq(n7*ns*(I-1)+ii)*f1 call cder(n7,esq,ns,nm,alq,aq,gq,es0,g,w) call WRITETTTQ(nm,ALQ,aq,GQ,E,iqlist) end c ============================================================ subroutine tqq(qf,u) implicit none real*8 qf(*),u(3,3),xx,xy,xz,yy,yz,zz xx=qf(1) xy=qf(2) xz=qf(3) yy=qf(4) yz=qf(5) zz=qf(6) c traceless quadrupole: u(1,1)=0.5d0*(xx+xx-yy-zz) u(2,2)=0.5d0*(yy+yy-zz-xx) u(3,3)=0.5d0*(zz+zz-xx-yy) u(1,2)=1.5d0*xy u(1,3)=1.5d0*xz u(2,1)=1.5d0*xy u(2,3)=1.5d0*yz u(3,1)=1.5d0*xz u(3,2)=1.5d0*yz return end c ============================================================ subroutine qbq(q,t) implicit none real*8 q(*),t(3,3) q(1)=t(1,1) q(2)=t(1,2) q(3)=t(1,3) q(4)=t(2,2) q(5)=t(2,3) q(6)=t(3,3) return end c =========================================================== 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(n7,es,r,nat,io,g,w) implicit none integer*4 n7,nat,io,ia,ix,is,ns,ii,kk real*8 es(*),r(nat,3),g,w 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 2 write(6,*)s10(is:len(s10))//'/JN.SCR.TXT' open(9,file=s10(is:len(s10))//'/JN.SCR.TXT') read(9,*)ns,g,w do 5 ii=1,ns 5 read(9,*)(es(kk+n7*(ii-1)),kk=1,n7) 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 C ========================================== subroutine cder(n7,esq,ns,nm,alq,aq,gq,es0,g,w) implicit none integer*4 n7,ns,nm,j,I,ix,jj,jx,kx real*8 esq(*),alq(*),aq(*),gq(*),es0(*),f1,f2,dl(3),r(3),qv(3,3), 1dlq(3),rq(3),qvq(3,3),qf(6),tq(3,3),f1q,f2q,g,w,dd,d,wq,w0, 1f1w,f2w,f1qa,f2qa do 307 j=1,ns f1 =es0( 1+n7*(j-1)) f2 =es0( 2+n7*(j-1)) do 11 ix=1,3 dl(ix) =es0(2+ix+n7*(j-1)) 11 r( ix) =es0(5+ix+n7*(j-1)) jj=0 do 12 jx=1,3 do 12 ix=1,3 jj=jj+1 12 qv(ix,jx)=es0(8+jj+n7*(j-1)) w0 =es0(18+n7*(j-1)) dd=w0**2-w**2 d=dd**2+g**2*w**2 c df1/dwjn,df2/dwjn: f1w=2.0d0/d*(dd+2.0d0*w0**2-4.0d0*(w0*dd)**2/d) f2w=4.0d0/d*( w0 -2.0d0* w0* dd**2/d) do 307 I=1,nm f1q=esq(1+n7*(j-1)+n7*ns*(I-1)) f2q=esq(2+n7*(j-1)+n7*ns*(I-1)) do 1 ix=1,3 dlq(ix)=esq(2+ix+n7*(j-1)+n7*ns*(I-1)) 1 rq( ix)=esq(5+ix+n7*(j-1)+n7*ns*(I-1)) jj=0 do 2 jx=1,3 do 2 ix=1,3 jj=jj+1 2 qvq(ix,jx)=esq(8+jj+n7*(j-1)+n7*ns*(I-1)) c dwjn/dq: wq=esq(18+n7*(j-1)+n7*ns*(I-1)) c take analytical f1q and f2q derivatives: f1qa=f1w*wq f2qa=f2w*wq c write(6,*)f1q,f1qa,f2q,f2qa do 3 ix=1,3 do 3 jx=1,3 jj=ix+3*(jx-1)+9*(I-1) alq(jj)=alq(jj)+f1qa*dl( ix)*dl( jx) 1 +f1 *(dlq(ix)*dl( jx)+dl( ix)*dlq(jx)) gq( jj)=gq (jj)+f2qa*dl( ix)*r( jx) 1 +f2 *(dlq(ix)*r( jx)+dl( ix)*rq (jx)) do 3 kx=1,3 jj=ix+3*(jx-1)+9*(kx-1)+27*(I-1) 3 aq(jj)=aq(jj)+f1qa*dl( ix)*qv( jx,kx) 1 +f1 *(dlq(ix)*qv( jx,kx)+dl( ix)*qvq(jx,kx)) 307 continue do 8 I=1,nm do 8 ix=1,3 c Quadrupole order: 'XX','XY','XZ','YY','YZ','ZZ' jj=0 do 81 jx=1,3 do 81 kx=jx,3 jj=jj+1 81 qf(jj)=aq(ix+3*(jx-1)+9*(kx-1)+27*(I-1)) c get traceless quadrupole, multiply by 3/2: call tqq(qf,tq) c get back, same order: call qbq(qf,tq) jj=0 do 8 jx=1,3 do 8 kx=jx,3 jj=jj+1 aq(ix+3*(jx-1)+9*(kx-1)+27*(I-1))=qf(jj) 8 aq(ix+3*(kx-1)+9*(jx-1)+27*(I-1))=qf(jj) return end